Ignore:
Files:
11 added
12 deleted
24 edited

Legend:

Unmodified
Added
Removed

  • branches/FLEXPART_9.1.3/src/FLEXPART.f90

    r14 r14  
    4949  integer :: i,j,ix,jy,inest
    5050  integer :: idummy = -320
     51  character(len=256) :: pathfile
    5152
    5253  ! Generate a large number of random numbers
    … …  
    6061  ! Print the GPL License statement
    6162  !*******************************************************
    62   print*,'Welcome to FLEXPART Version 9.0'
     63  ! print*,'Welcome to FLEXPART Version 9.1 (Build 20121029)'
     64  print*,'Welcome to FLEXPART Version 9.1.3 (Build 2013158)'
    6365  print*,'FLEXPART is free software released under the GNU Genera'// &
    6466       'l Public License.'
    … …  
    6769  !*******************************************************
    6870
    69   call readpaths
     71  select case (iargc())
     72  case (1)
     73    call getarg(1,pathfile)
     74    if (trim(pathfile).eq.'-v') then
     75        !write(*,*) 'FLEXPART Version 9.1.3 (Build 2013158)'
     76        stop
     77    endif 
     78  case (0)
     79    write(pathfile,'(a11)') './pathnames'
     80  end select
     81
     82  call readpaths(pathfile)
     83  print*,length(4)
    7084
    7185  ! Read the user specifications for the current model run

  • branches/FLEXPART_9.1.3/src/com_mod.f90

    • Property svn:executable deleted
    r14 r14  
    77!        June 1996                                                             *
    88!                                                                              *
    9 !        Last update: 9 August 2000                                            *
     9!        Last update:15 August 2013 IP                                         *
    1010!                                                                              *
    1111!*******************************************************************************
    … …  
    139139  real :: decay(maxspec)
    140140  real :: weta(maxspec),wetb(maxspec)
     141! NIK: 31.01.2013- parameters for in-cloud scavening
     142  real :: weta_in(maxspec), wetb_in(maxspec), wetc_in(maxspec), wetd_in(maxspec)
    141143  real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
    142144  real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
    … …  
    172174
    173175  ! WET DEPOSITION
    174   ! weta, wetb              parameters for determining wet scavenging coefficients
     176  ! weta, wetb              parameters for determining below-cloud wet scavenging coefficients
     177  ! weta_in, wetb_in       parameters for determining in-cloud wet scavenging coefficients
     178  ! wetc_in, wetd_in       parameters for determining in-cloud wet scavenging coefficients
    175179
    176180  ! GAS DEPOSITION
    … …  
    331335  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    332336  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
     337  !scavenging NIK, PS
    333338  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    334339  integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
     340      integer icloudbot(0:nxmax-1,0:nymax-1,2)
     341      integer icloudthck(0:nxmax-1,0:nymax-1,2)
     342
    335343
    336344  ! uu,vv,ww [m/2]       wind components in x,y and z direction
    … …  
    346354  !      rainout  conv/lsp dominated  2/3
    347355  !      washout  conv/lsp dominated  4/5
     356! PS 2013
     357!c icloudbot (m)        cloud bottom height
     358!c icloudthck (m)       cloud thickness     
     359
    348360  ! pplev for the GFS version
    349361

  • branches/FLEXPART_9.1.3/src/concoutput.f90

    r14 r14  
    158158        yl=outlat0+real(jy)*dyout
    159159        xl=(xl-xlon0)/dx
    160         yl=(yl-ylat0)/dx
     160        yl=(yl-ylat0)/dy !v9.1.1
    161161        iix=max(min(nint(xl),nxmin1),0)
    162162        jjy=max(min(nint(yl),nymin1),0)

  • branches/FLEXPART_9.1.3/src/erf.f90

    • Property svn:executable deleted
    r14 r14  
    4242
    4343  integer :: j
    44   real :: x,tmp,ser,xx,gammln
    45   real :: cof(6) = (/ &
     44  real(kind=dp) :: x,tmp,ser,xx,gammln
     45  real(KIND=dp) :: cof(6) = (/ &
    4646       76.18009173_dp, -86.50532033_dp, 24.01409822_dp, &
    4747       -1.231739516_dp, .120858003e-2_dp, -.536382e-5_dp    /)
    48   real :: stp = 2.50662827465_dp
    49   real :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
     48  real(KIND=dp) :: stp = 2.50662827465_dp
     49  real(KIND=dp) :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
    5050
    5151  x=xx-one
    … …  
    6262function gammp(a,x)
    6363
    64   implicit none
    65 
    66   real :: a, x, gln, gamser, gammp, gammcf
     64  use par_mod, only: dp
     65
     66  implicit none
     67
     68  real(KIND=dp) :: a, x, gln, gamser, gammp, gammcf
    6769
    6870  if(x .lt. 0. .or. a .le. 0.) then
    … …  
    8183function gammq(a,x)
    8284
    83   implicit none
    84 
    85   real :: a, x, gln, gamser, gammq, gammcf
     85  use par_mod, only: dp
     86
     87  implicit none
     88
     89  real(KIND=dp) :: a, x, gln, gamser, gammq, gammcf
    8690
    8791  if(x.lt.0..or.a.le.0.) then
    … …  
    100104subroutine gser(gamser,a,x,gln)
    101105
     106  use par_mod, only: dp
     107
    102108  implicit none
    103109
    104110  integer :: n
    105   real :: gamser, a, x, gln, ap, summ, del
    106   real, external :: gammln
     111  real(KIND=dp) :: gamser, a, x, gln, ap, summ, del
     112  real(KIND=dp), external :: gammln
    107113
    108114  integer,parameter :: itmax=100
    … …  
    134140subroutine gcf(gammcf,a,x,gln)
    135141
     142  use par_mod, only: dp
     143
    136144  implicit none
    137145
    138146  integer :: n
    139   real :: gammcf, a, x, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
    140   real, external :: gammln
     147  real(KIND=4) :: gammcf, x, a, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
     148  real(KIND=dp), external :: gammln
    141149
    142150  integer,parameter :: itmax=100
    … …  
    172180function erf(x)
    173181
    174   implicit none
    175 
    176   real :: x, erf
    177   real, external :: gammp
     182  use par_mod, only: dp
     183
     184  implicit none
     185
     186  real(KIND=dp) :: x, erf
     187  real(KIND=dp), external :: gammp
    178188
    179189  if(x.lt.0.)then
    … …  
    186196function erfc(x)
    187197
    188   implicit none
    189 
    190   real :: x, erfc
    191   real, external :: gammp, gammq
     198  use par_mod, only: dp
     199
     200  implicit none
     201
     202  real(KIND=dp) :: x, erfc
     203  real(KIND=dp), external :: gammp, gammq
    192204
    193205  if(x.lt.0.)then
    … …  
    200212function erfcc(x)
    201213
    202   implicit none
    203 
    204   real :: x, z, t, erfcc
     214  use par_mod, only: dp
     215
     216  implicit none
     217
     218  real(KIND=dp) :: x, z, t, erfcc
    205219
    206220  z=abs(x)

  • branches/FLEXPART_9.1.3/src/getvdep.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/gridcheck.f90

    • Property svn:executable deleted
    r14 r14  
    9898
    9999  integer :: isec1(56),isec2(22+nxmax+nymax)
    100   real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
     100  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
     101  real(kind=4) :: zsec2(184),zsec4(jpunp)
    101102  character(len=1) :: opt
    102103
    … …  
    466467    k=nlev_ec+1+numskip+i
    467468    akm(nwz-i+1)=zsec2(j)
    468   !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
    469469    bkm(nwz-i+1)=zsec2(k)
    470470  end do

  • branches/FLEXPART_9.1.3/src/gridcheck_gfs.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/gridcheck_nests.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/interpol_rain.f90

    r14 r14  
    2020!**********************************************************************
    2121
    22 subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
    23        ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
    24   !                          i   i   i    i    i     i   i
    25   !i    i    i  i    i     i      i      i     o     o     o
     22!subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
     23!       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
     24!  !                          i   i   i    i    i     i   i
     25!  !i    i    i  i    i     i      i      i     o     o     o
     26
     27      subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
     28     ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
     29     intiy1,intiy2,icmv)
     30!                               i   i   i    i    i     i   i
     31!     i    i    i  i    i     i      i      i     o     o     o
     32
    2633  !****************************************************************************
    2734  !                                                                           *
    … …  
    4047  !                                                                           *
    4148  !     30 August 1996                                                        *
    42   !                                                                           *
     49  !
     50  !*  Petra Seibert, 2011/2012:
     51  !*  Add interpolation of cloud bottom and cloud thickness
     52  !*  for fix to SE's new wet scavenging scheme  *
    4353  !****************************************************************************
    4454  !                                                                           *
    … …  
    7080
    7181  integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
    72   integer :: itime,itime1,itime2,level,indexh
     82  !integer :: itime,itime1,itime2,level,indexh
     83  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
     84  integer :: intiy1,intiy2,ipsum,icmv 
    7385  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    7486  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    7587  real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
    76   real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
    77   real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
     88  integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
     89  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
     90  !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
     91  real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4 
     92  !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
    7893
    7994
    … …  
    127142         + p3*yy3(ix ,jyp,level,indexh) &
    128143         + p4*yy3(ixp,jyp,level,indexh)
     144
     145!CPS clouds:
     146        ip1=1
     147        ip2=1
     148        ip3=1
     149        ip4=1
     150        if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
     151        if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
     152        if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
     153        if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
     154        ipsum= ip1+ip2+ip3+ip4
     155        if (ipsum .eq. 0) then
     156          yi1(m)=icmv
     157        else
     158          yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh)    &
     159             + ip2*p2*iy1(ixp,jy ,indexh)       &
     160             + ip3*p3*iy1(ix ,jyp,indexh)       &
     161             + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
     162        endif
     163       
     164        ip1=1
     165        ip2=1
     166        ip3=1
     167        ip4=1
     168        if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
     169        if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
     170        if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
     171        if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
     172        ipsum= ip1+ip2+ip3+ip4
     173        if (ipsum .eq. 0) then
     174          yi2(m)=icmv
     175        else
     176          yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh)  &
     177             + ip2*p2*iy2(ixp,jy ,indexh)     &
     178             + ip3*p3*iy2(ix ,jyp,indexh)     &
     179             + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
     180        endif
     181!CPS end clouds
     182
    129183  end do
    130184
    … …  
    134188  !************************************
    135189
     190      if (itime .lt. itime1) then
     191        print*,'interpol_rain.f90'
     192        print*,itime,itime1,itime2
     193        stop 'ITIME PROBLEM'
     194      endif
     195
     196
    136197  dt1=real(itime-itime1)
    137198  dt2=real(itime2-itime)
    … …  
    143204
    144205
     206!PS clouds:
     207      intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
     208      if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
     209      if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
     210
     211      intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
     212      if (yi2(1) .eq. float(icmv)) intiy2=yi2(2) 
     213      if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
     214     
     215      if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
     216        intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
     217      else
     218        intiy1=icmv
     219        intiy2=icmv
     220      endif
     221!PS end clouds
     222
    145223end subroutine interpol_rain

  • branches/FLEXPART_9.1.3/src/makefile.ecmwf_gfortran

    r14 r14  
    99LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
    1010#LIBPATH2 =   /flex_wrk/flexpart/lib64/gfortran/lib/
    11 #FFLAGS   =   -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    1211FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    1312#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)

  • branches/FLEXPART_9.1.3/src/par_mod.f90

    • Property svn:executable deleted
    r14 r14  
    2828!        1997                                                                  *
    2929!                                                                              *
    30 !        Last update 10 August 2000                                            *
     30!        Last update 15 August 2013 IP                                         *
    3131!                                                                              *
    3232!*******************************************************************************
    … …  
    122122
    123123  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
    124   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
     124  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    125125  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    126   integer,parameter :: nxshift=359  ! for ECMWF
     126  !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
     127  integer,parameter :: nxshift=359 ! for ECMWF
    127128  !integer,parameter :: nxshift=0     ! for GFS
    128129
    … …  
    150151  !*********************************************
    151152
    152   integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    153   !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
     153  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
     154  integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
    154155
    155156  ! maxnests                maximum number of nested grids
    … …  
    194195  !**************************************************
    195196
    196   integer,parameter :: maxpart=2000000
    197   integer,parameter :: maxspec=4
     197  integer,parameter :: maxpart=150000
     198  integer,parameter :: maxspec=1
    198199
    199200
    … …  
    258259  integer,parameter :: unitboundcond=89
    259260
     261!******************************************************
     262! integer code for missing values, used in wet scavenging (PS, 2012)
     263!******************************************************
     264
     265      integer icmv
     266      parameter(icmv=-9999)
     267
    260268end module par_mod

  • branches/FLEXPART_9.1.3/src/partdep.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/readageclasses.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/readavailable.f90

    • Property svn:executable deleted
    r14 r14  
    123123
    124124  do k=1,numbnests
     125  print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
     126  print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
    125127    open(unitavailab,file=path(numpath+2*(k-1)+2) &
    126128         (1:length(numpath+2*(k-1)+2)),status='old',err=998)
    … …  
    275277  return
    276278
    277 998   write(*,*) ' #### FLEXPART MODEL ERROR! FILE   #### '
     279998   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE   #### '
    278280  write(*,'(a)') '     '//path(numpath+2*(k-1)+2) &
    279281       (1:length(numpath+2*(k-1)+2))
    … …  
    281283  stop
    282284
    283 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
     285999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
    284286  write(*,'(a)') '     '//path(4)(1:length(4))
    285287  write(*,*) ' #### CANNOT BE OPENED           #### '

  • branches/FLEXPART_9.1.3/src/readcommand.f90

    r14 r14  
    8080  character(len=50) :: line
    8181  logical :: old
    82 
     82  logical :: nmlout=.false.
     83  integer :: readerror
     84
     85  namelist /command/ &
     86    ldirect, &
     87    ibdate,ibtime, &
     88    iedate,ietime, &
     89    loutstep, &
     90    loutaver, &
     91    loutsample, &
     92    itsplit, &
     93    lsynctime, &
     94    ctl, &
     95    ifine, &
     96    iout, &
     97    ipout, &
     98    lsubgrid, &
     99    lconvection, &
     100    lagespectra, &
     101    ipin, &
     102    ioutputforeachrelease, &
     103    iflux, &
     104    mdomainfill, &
     105    ind_source, &
     106    ind_receptor, &
     107    mquasilag, &
     108    nested_output, &
     109    linit_cond
     110
     111  ! Presetting namelist command
     112  ldirect=1
     113  ibdate=20000101
     114  ibtime=0
     115  iedate=20000102
     116  ietime=0
     117  loutstep=10800
     118  loutaver=10800
     119  loutsample=900
     120  itsplit=999999999
     121  lsynctime=900
     122  ctl=-5.0
     123  ifine=4
     124  iout=3
     125  ipout=0
     126  lsubgrid=1
     127  lconvection=1
     128  lagespectra=0
     129  ipin=1
     130  ioutputforeachrelease=0
     131  iflux=1
     132  mdomainfill=0
     133  ind_source=1
     134  ind_receptor=1
     135  mquasilag=0
     136  nested_output=0
     137  linit_cond=0
    83138
    84139  ! Open the command file and read user options
    85   !********************************************
    86 
    87 
     140  ! Namelist input first: try to read as namelist file
     141  !**************************************************************************
    88142  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    89        err=999)
    90 
    91   ! Check the format of the COMMAND file (either in free format,
    92   ! or using formatted mask)
    93   ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
    94   !**************************************************************************
    95 
    96   call skplin(9,unitcommand)
    97   read (unitcommand,901) line
    98 901   format (a)
    99   if (index(line,'LDIRECT') .eq. 0) then
    100     old = .false.
    101   else
    102     old = .true.
    103   endif
    104   rewind(unitcommand)
    105 
    106   ! Read parameters
    107   !****************
    108 
    109   call skplin(7,unitcommand)
    110   if (old) call skplin(1,unitcommand)
    111 
    112   read(unitcommand,*) ldirect
    113   if (old) call skplin(3,unitcommand)
    114   read(unitcommand,*) ibdate,ibtime
    115   if (old) call skplin(3,unitcommand)
    116   read(unitcommand,*) iedate,ietime
    117   if (old) call skplin(3,unitcommand)
    118   read(unitcommand,*) loutstep
    119   if (old) call skplin(3,unitcommand)
    120   read(unitcommand,*) loutaver
    121   if (old) call skplin(3,unitcommand)
    122   read(unitcommand,*) loutsample
    123   if (old) call skplin(3,unitcommand)
    124   read(unitcommand,*) itsplit
    125   if (old) call skplin(3,unitcommand)
    126   read(unitcommand,*) lsynctime
    127   if (old) call skplin(3,unitcommand)
    128   read(unitcommand,*) ctl
    129   if (old) call skplin(3,unitcommand)
    130   read(unitcommand,*) ifine
    131   if (old) call skplin(3,unitcommand)
    132   read(unitcommand,*) iout
    133   if (old) call skplin(3,unitcommand)
    134   read(unitcommand,*) ipout
    135   if (old) call skplin(3,unitcommand)
    136   read(unitcommand,*) lsubgrid
    137   if (old) call skplin(3,unitcommand)
    138   read(unitcommand,*) lconvection
    139   if (old) call skplin(3,unitcommand)
    140   read(unitcommand,*) lagespectra
    141   if (old) call skplin(3,unitcommand)
    142   read(unitcommand,*) ipin
    143   if (old) call skplin(3,unitcommand)
    144   read(unitcommand,*) ioutputforeachrelease
    145   if (old) call skplin(3,unitcommand)
    146   read(unitcommand,*) iflux
    147   if (old) call skplin(3,unitcommand)
    148   read(unitcommand,*) mdomainfill
    149   if (old) call skplin(3,unitcommand)
    150   read(unitcommand,*) ind_source
    151   if (old) call skplin(3,unitcommand)
    152   read(unitcommand,*) ind_receptor
    153   if (old) call skplin(3,unitcommand)
    154   read(unitcommand,*) mquasilag
    155   if (old) call skplin(3,unitcommand)
    156   read(unitcommand,*) nested_output
    157   if (old) call skplin(3,unitcommand)
    158   read(unitcommand,*) linit_cond
     143         form='formatted',iostat=readerror)
     144  ! If fail, check if file does not exist
     145  if (readerror.ne.0) then
     146
     147    print*,'***ERROR: file COMMAND not found. Check your pathnames file.'
     148    stop
     149
     150  endif
     151
     152  read(unitcommand,command,iostat=readerror)
    159153  close(unitcommand)
    160154
     155  ! If error in namelist format, try to open with old input code
     156  if (readerror.ne.0) then
     157
     158    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
     159         err=999)
     160
     161    ! Check the format of the COMMAND file (either in free format,
     162    ! or using formatted mask)
     163    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
     164    !**************************************************************************
     165
     166    call skplin(9,unitcommand)
     167    read (unitcommand,901) line
     168  901   format (a)
     169    if (index(line,'LDIRECT') .eq. 0) then
     170      old = .false.
     171    else
     172      old = .true.
     173    endif
     174    rewind(unitcommand)
     175
     176    ! Read parameters
     177    !****************
     178
     179    call skplin(7,unitcommand)
     180    if (old) call skplin(1,unitcommand)
     181
     182    read(unitcommand,*) ldirect
     183    if (old) call skplin(3,unitcommand)
     184    read(unitcommand,*) ibdate,ibtime
     185    if (old) call skplin(3,unitcommand)
     186    read(unitcommand,*) iedate,ietime
     187    if (old) call skplin(3,unitcommand)
     188    read(unitcommand,*) loutstep
     189    if (old) call skplin(3,unitcommand)
     190    read(unitcommand,*) loutaver
     191    if (old) call skplin(3,unitcommand)
     192    read(unitcommand,*) loutsample
     193    if (old) call skplin(3,unitcommand)
     194    read(unitcommand,*) itsplit
     195    if (old) call skplin(3,unitcommand)
     196    read(unitcommand,*) lsynctime
     197    if (old) call skplin(3,unitcommand)
     198    read(unitcommand,*) ctl
     199    if (old) call skplin(3,unitcommand)
     200    read(unitcommand,*) ifine
     201    if (old) call skplin(3,unitcommand)
     202    read(unitcommand,*) iout
     203    if (old) call skplin(3,unitcommand)
     204    read(unitcommand,*) ipout
     205    if (old) call skplin(3,unitcommand)
     206    read(unitcommand,*) lsubgrid
     207    if (old) call skplin(3,unitcommand)
     208    read(unitcommand,*) lconvection
     209    if (old) call skplin(3,unitcommand)
     210    read(unitcommand,*) lagespectra
     211    if (old) call skplin(3,unitcommand)
     212    read(unitcommand,*) ipin
     213    if (old) call skplin(3,unitcommand)
     214    read(unitcommand,*) ioutputforeachrelease
     215    if (old) call skplin(3,unitcommand)
     216    read(unitcommand,*) iflux
     217    if (old) call skplin(3,unitcommand)
     218    read(unitcommand,*) mdomainfill
     219    if (old) call skplin(3,unitcommand)
     220    read(unitcommand,*) ind_source
     221    if (old) call skplin(3,unitcommand)
     222    read(unitcommand,*) ind_receptor
     223    if (old) call skplin(3,unitcommand)
     224    read(unitcommand,*) mquasilag
     225    if (old) call skplin(3,unitcommand)
     226    read(unitcommand,*) nested_output
     227    if (old) call skplin(3,unitcommand)
     228    read(unitcommand,*) linit_cond
     229    close(unitcommand)
     230
     231  endif ! input format
     232
     233  ! write command file in namelist format to output directory if requested
     234  if (nmlout.eqv..true.) then
     235    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',status='new',err=999)
     236    write(unitcommand,nml=command)
     237    close(unitcommand)
     238  endif
     239
    161240  ifine=max(ifine,1)
    162 
    163241
    164242  ! Determine how Markov chain is formulated (for w or for w/sigw)

  • branches/FLEXPART_9.1.3/src/readdepo.f90

    • Property svn:executable deleted

  • branches/FLEXPART_9.1.3/src/readpartpositions.f90

    r14 r14  
    7777
    7878  read(unitpartin) numpointin
    79   if (numpointin.ne.numpoint) then
    80      write(*,*) ' ####              WARNING                     #### '
    81      write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
    82      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
    83      write(*,*) ' #### THIS IS JUST A CHECK TO SEE IF YOU READ  #### '
    84      write(*,*) ' #### THE CORRECT PARTICLE DUMP FILE.          #### '
    85   endif
     79  if (numpointin.ne.numpoint) goto 995
    8680  do i=1,numpointin
    8781    read(unitpartin)
    … …  
    152146  stop
    153147
     148995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
     149  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
     150  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
     151  stop
    154152
    155153996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '

  • branches/FLEXPART_9.1.3/src/readpaths.f90

    • Property svn:executable deleted
    r14 r14  
    2020!**********************************************************************
    2121
    22 subroutine readpaths
     22subroutine readpaths(pathfile)
    2323
    2424  !*****************************************************************************
    … …  
    3030  !                                                                            *
    3131  !     1 February 1994                                                        *
     32  !     last modified                                                          *
     33  !     HS, 7.9.2012                                                           *
     34  !     option to give pathnames file as command line option                   *
    3235  !                                                                            *
    3336  !*****************************************************************************
    … …  
    4750  implicit none
    4851
    49   integer :: i
     52  integer   :: i
     53  character(256) :: pathfile
    5054
    5155  ! Read the pathname information stored in unitpath
    5256  !*************************************************
    5357
    54 
    55   open(unitpath,file='pathnames',status='old',err=999)
     58  open(unitpath,file=trim(pathfile),status='old',err=999)
    5659
    5760  do i=1,numpath
    … …  
    7073    length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
    7174  end do
     75  print*,length(5),length(6)
    7276
    7377

  • branches/FLEXPART_9.1.3/src/readreleases.f90

    • Property svn:executable deleted
    r14 r14  
    3333  !     Update: 29 January 2001                                                *
    3434  !     Release altitude can be either in magl or masl                         *
    35   !                                                                            *
     35  !     Update: Aug 2013 IP                                                                        *
    3636  !*****************************************************************************
    3737  !                                                                            *
    … …  
    5555  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
    5656  !                     area                                                   *
    57   ! weta, wetb          parameters to determine the wet scavenging coefficient *
     57  ! before: weta, wetb          parameters to determine the wet scavenging coefficient *
     58  ! weta, wetb          parameters to determine the below-cloud scavenging (NIK)     *
     59  ! weta_in, wetb_in    parameters to determine the in-cloud scavenging    (NIK)    *
     60  ! wetc_in, wetd_in    parameters to determine the in-cloud scavenging    (NIK)    * 
    5861  ! zpoint1,zpoint2     height range, over which release takes place           *
     62  ! num_min_discrete    if less, release cannot be randomized and happens at   *
     63  !                     time mid-point of release interval                     *
    5964  !                                                                            *
    6065  !*****************************************************************************
    … …  
    6873
    6974  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
     75  integer,parameter :: num_min_discrete=100
    7076  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    71   real(kind=dp) :: jul1,jul2,juldate
     77  real(kind=dp) :: jul1,jul2,julm,juldate
    7278  character(len=50) :: line
    7379  logical :: old
    … …  
    277283    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
    278284
     285
     286
    279287  ! Check if wet deposition or OH reaction shall be calculated
    280288  !***********************************************************
    281289    if (weta(i).gt.0.)  then
    282290      WETDEP=.true.
    283       write (*,*) 'Wetdeposition switched on: ',weta(i),i
    284     endif
     291      !write (*,*) 'Wetdeposition switched on: ',weta(i),i
     292       write (*,*) 'Wet deposition switched on'
     293       write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
     294   endif
     295
     296! NIK 31.01.2013
     297    if (weta_in(i).gt.0.)  then
     298      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
     299    endif
     300
    285301    if (ohreact(i).gt.0) then
    286302      OHREA=.true.
    … …  
    376392  jul1=juldate(id1,it1)
    377393  jul2=juldate(id2,it2)
     394  julm=(jul1+jul2)/2.
    378395  if (jul1.gt.jul2) then
    379396    write(*,*) 'FLEXPART MODEL ERROR'
    … …  
    391408        stop
    392409      endif
    393       ireleasestart(numpoint)=int((jul1-bdate)*86400.)
    394       ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     410      if (npart(numpoint).gt.num_min_discrete) then
     411        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     412        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     413      else
     414        ireleasestart(numpoint)=int((julm-bdate)*86400.)
     415        ireleaseend(numpoint)=int((julm-bdate)*86400.)
     416      endif
    395417    else if (ldirect.eq.-1) then
    396418      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
    … …  
    401423        stop
    402424      endif
    403       ireleasestart(numpoint)=int((jul1-bdate)*86400.)
    404       ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     425      if (npart(numpoint).gt.num_min_discrete) then
     426        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     427        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     428      else
     429        ireleasestart(numpoint)=int((julm-bdate)*86400.)
     430        ireleaseend(numpoint)=int((julm-bdate)*86400.)
     431      endif
    405432    endif
    406433  endif

  • branches/FLEXPART_9.1.3/src/readspecies.f90

    r14 r14  
    3131  !   11 July 1996                                                             *
    3232  !                                                                            *
     33  !   Changes:                                                                 *
     34  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
     35  !                                                                            *
    3336  !*****************************************************************************
    3437  !                                                                            *
    … …  
    3639  ! decaytime(maxtable)  half time for radiological decay                      *
    3740  ! specname(maxtable)   names of chemical species, radionuclides              *
    38   ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
     41  ! weta, wetb           Parameters for determining below-cloud scavenging     *
     42  ! weta_in              Parameter for determining in-cloud scavenging         *
     43  ! wetb_in              Parameter for determining in-cloud scavenging         *
     44  ! wetc_in              Parameter for determining in-cloud scavenging         *
     45  ! wetd_in              Parameter for determining in-cloud scavenging         *
    3946  ! ohreact              OH reaction rate                                      *
    4047  ! id_spec              SPECIES number as referenced in RELEASE file          *
    … …  
    7885    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    7986  !  write(*,*) wetb(pos_spec)
     87
     88!*** NIK 31.01.2013: including in-cloud scavening parameters
     89   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
     90  !  write(*,*) weta_in(pos_spec)
     91   read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
     92  !  write(*,*) wetb_in(pos_spec)
     93   read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
     94  !  write(*,*) wetc_in(pos_spec)
     95   read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
     96  !  write(*,*) wetd_in(pos_spec)
     97
    8098    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    8199  !  write(*,*) reldiff(pos_spec)

  • branches/FLEXPART_9.1.3/src/readwind.f90

    r14 r14  
    101101  character(len=24) :: gribErrorMsg = 'Error reading grib file'
    102102  character(len=20) :: gribFunction = 'readwind'
     103
     104  !HSO conversion of ECMWF etadot to etadot*dp/deta
     105  logical :: etacon=.false.
     106  real,parameter :: p00=101325.
     107  real :: dak,dbk
    103108
    104109  hflswitch=.false.
    … …  
    365370  endif
    366371
     372  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
     373  if(etacon.eqv..true.) then
     374    do k=1,nwzmax
     375      dak=akm(k+1)-akm(k)
     376      dbk=bkm(k+1)-bkm(k)
     377      do i=0,nxmin1
     378        do j=0,nymin1
     379          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
     380          if (k.gt.1) then
     381            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
     382          endif
     383        end do
     384      end do
     385    end do
     386  endif
     387
    367388  ! For global fields, assign the leftmost data column also to the rightmost
    368389  ! data column; if required, shift whole grid by nxshift grid points

  • branches/FLEXPART_9.1.3/src/readwind_gfs.f90

    r14 r14  
    703703  endif
    704704
    705 
    706705  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
    707706  if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'

  • branches/FLEXPART_9.1.3/src/verttransform.f90

    r14 r14  
    4949  ! Sabine Eckhardt, March 2007
    5050  ! added the variable cloud for use with scavenging - descr. in com_mod
     51  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
     52  ! note that also other subroutines are affected by the fix
    5153  !*****************************************************************************
    5254  !                                                                            *
    … …  
    7072
    7173  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    72   integer :: rain_cloud_above,kz_inv
     74 integer :: rain_cloud_above,kz_inv !SE
     75  integer icloudtop !PS
    7376  real :: f_qvsat,pressure
    74   real :: rh,lsp,convp
     77 !real :: rh,lsp,convp
     78  real :: rh,lsp,convp,prec,rhmin 
    7579  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
    7680  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    7781  real :: xlon,ylat,xlonr,dzdx,dzdy
    78   real :: dzdx1,dzdx2,dzdy1,dzdy2
     82  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
    7983  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8084  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8387  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8488  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
     89  logical lconvectprec
    8590  real,parameter :: const=r_air/ga
     91  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    8692
    8793  logical :: init = .true.
    … …  
    545551  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    546552  !   total cloudheight is stored at level 0
     553
     554
     555
    547556  do jy=0,nymin1
    548557    do ix=0,nxmin1
    549       rain_cloud_above=0
    550       lsp=lsprec(ix,jy,1,n)
    551       convp=convprec(ix,jy,1,n)
    552       cloudsh(ix,jy,n)=0
    553       do kz_inv=1,nz-1
    554          kz=nz-kz_inv+1
    555          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
    556          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
    557          clouds(ix,jy,kz,n)=0
    558          if (rh.gt.0.8) then ! in cloud
    559             if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    560                rain_cloud_above=1
    561                cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
    562                     height(kz)-height(kz-1)
    563                if (lsp.ge.convp) then
    564                   clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
    565                else
    566                   clouds(ix,jy,kz,n)=2 ! convp dominated rainout
    567                endif
    568             else ! no precipitation
    569                   clouds(ix,jy,kz,n)=1 ! cloud
     558
     559
     560
     561  !    rain_cloud_above=0
     562  !    lsp=lsprec(ix,jy,1,n)
     563  !    convp=convprec(ix,jy,1,n)
     564  !    cloudsh(ix,jy,n)=0
     565  !    do kz_inv=1,nz-1
     566  !       kz=nz-kz_inv+1
     567  !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     568  !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     569  !       clouds(ix,jy,kz,n)=0
     570  !       if (rh.gt.0.8) then ! in cloud
     571  !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
     572  !             rain_cloud_above=1
     573  !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
     574  !                  height(kz)-height(kz-1)
     575  !             if (lsp.ge.convp) then
     576  !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
     577  !             else
     578  !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
     579  !             endif
     580  !          else ! no precipitation
     581  !                clouds(ix,jy,kz,n)=1 ! cloud
     582  !          endif
     583  !       else ! no cloud
     584  !          if (rain_cloud_above.eq.1) then ! scavenging
     585  !             if (lsp.ge.convp) then
     586  !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
     587  !             else
     588  !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
     589  !             endif
     590  !          endif
     591  !       endif
     592  !    end do
     593
     594
     595   ! PS 3012
     596
     597             lsp=lsprec(ix,jy,1,n)
     598          convp=convprec(ix,jy,1,n)
     599          prec=lsp+convp
     600          if (lsp.gt.convp) then !  prectype='lsp'
     601            lconvectprec = .false.
     602          else ! prectype='cp '
     603            lconvectprec = .true.
     604          endif
     605          rhmin = 0.90 ! standard condition for presence of clouds
     606!PS       note that original by Sabine Eckhart was 80%
     607!PS       however, for T<-20 C we consider saturation over ice
     608!PS       so I think 90% should be enough         
     609          icloudbot(ix,jy,n)=icmv
     610          icloudtop=icmv ! this is just a local variable
     61198        do kz=1,nz
     612            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
     613            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
     614!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
     615            if (rh .gt. rhmin) then
     616              if (icloudbot(ix,jy,n) .eq. icmv) then
     617                icloudbot(ix,jy,n)=nint(height(kz))
     618              endif
     619              icloudtop=nint(height(kz)) ! use int to save memory
    570620            endif
    571          else ! no cloud
    572             if (rain_cloud_above.eq.1) then ! scavenging
    573                if (lsp.ge.convp) then
    574                   clouds(ix,jy,kz,n)=5 ! lsp dominated washout
    575                else
    576                   clouds(ix,jy,kz,n)=4 ! convp dominated washout
    577                endif
     621          enddo
     622
     623!PS try to get a cloud thicker than 50 m
     624!PS if there is at least .01 mm/h  - changed to 0.002 and put into
     625!PS parameter precpmin       
     626          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
     627              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
     628              prec .gt. precmin) then
     629            rhmin = rhmin - 0.05
     630            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
     631          endif
     632!PS implement a rough fix for badly represented convection
     633!PS is based on looking at a limited set of comparison data
     634          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
     635             prec .gt. precmin) then
     636
     637            if (convp .lt. 0.1) then
     638              icloudbot(ix,jy,n) = 500
     639              icloudtop =         8000
     640            else
     641              icloudbot(ix,jy,n) = 0
     642              icloudtop =      10000
    578643            endif
    579          endif
    580       end do
     644          endif
     645          if (icloudtop .ne. icmv) then
     646            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
     647          else
     648            icloudthck(ix,jy,n) = icmv
     649          endif
     650!PS  get rid of too thin clouds     
     651          if (icloudthck(ix,jy,n) .lt. 50) then
     652            icloudbot(ix,jy,n)=icmv
     653            icloudthck(ix,jy,n)=icmv
     654          endif
     655
     656
    581657    end do
    582658  end do
    … …  
    605681  !104   continue
    606682  ! close(4)
     683
     684
    607685end subroutine verttransform

  • branches/FLEXPART_9.1.3/src/wetdepo.f90

    r14 r14  
    3939  ! Code may not be correct for decay of deposition!                           *
    4040  !                                                                            *
     41  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
     42  ! scheme, not dated
     43  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
    4144  !*****************************************************************************
    4245  !                                                                            *
    … …  
    7174
    7275  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    73   integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
    74   integer :: ks, kp
     76  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
     77  integer :: kz !PS scheme
     78  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
     79  integer :: scheme_number ! NIK==1, PS ==2
    7580  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    7681  real :: clouds_h ! cloud height for the specific grid point
    77   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
     82  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
    7883  real :: wetdeposit(maxspec),restmass
    7984  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    146151    interp_time=nint(itime-0.5*ltsample)
    147152
    148     if (ngrid.eq.0) then
    149       call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    150            1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
    151            memtime(1),memtime(2),interp_time,lsp,convp,cc)
    152     else
    153       call interpol_rain_nests(lsprecn,convprecn,tccn, &
    154            nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
    155            memtime(1),memtime(2),interp_time,lsp,convp,cc)
    156     endif
    157 
    158     if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
     153!    PS nest case still needs to be implemented!!   
     154!    if (ngrid.eq.0) then
     155!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
     156!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
     157!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     158           call interpol_rain(lsprec,convprec,tcc,     &
     159            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
     160            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
     161            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv) 
     162!    else
     163!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
     164!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
     165!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
     166!    endif
     167
     168
     169 !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
     170 !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
     171        prec = lsp+convp
     172        precsub = 0.01
     173        if (prec .lt. precsub) then
     174          goto 20
     175        else
     176          f = (prec-precsub)/prec
     177          lsp = f*lsp
     178          convp = f*convp
     179        endif
     180
    159181
    160182  ! get the level were the actual particle is in
    161183      do il=2,nz
    162184        if (height(il).gt.ztra1(jpart)) then
    163           hz=il-1
     185          !hz=il-1
     186          kz=il-1
    164187          goto 26
    165188        endif
    … …  
    174197  ! scavenging is done
    175198
    176   if (ngrid.eq.0) then
    177      clouds_v=clouds(ix,jy,hz,n)
    178      clouds_h=cloudsh(ix,jy,n)
    179   else
    180      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
    181      clouds_h=cloudsnh(ix,jy,n,ngrid)
    182   endif
    183   !write(*,*) 'there is
    184   !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
    185   if (clouds_v.le.1) goto 20
    186   !write (*,*) 'there is scavenging'
     199!old scheme
     200!  if (ngrid.eq.0) then
     201!     clouds_v=clouds(ix,jy,hz,n)
     202!     clouds_h=cloudsh(ix,jy,n)
     203!  else
     204!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
     205!     clouds_h=cloudsnh(ix,jy,n,ngrid)
     206!  endif
     207!  !write(*,*) 'there is
     208!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
     209!  if (clouds_v.le.1) goto 20
     210!  !write (*,*) 'there is scavenging'
     211
     212  ! PS: part of 2011/2012 fix
     213
     214        if (ztra1(jpart) .le. float(ictop)) then
     215          if (ztra1(jpart) .gt. float(icbot)) then
     216            indcloud = 2 ! in-cloud
     217          else
     218            indcloud = 1 ! below-cloud
     219          endif
     220        elseif (ictop .eq. icmv) then
     221          indcloud = 0 ! no cloud found, use old scheme
     222        else
     223          goto 20 ! above cloud
     224        endif
     225
    187226
    188227  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    228267  !**********************************************************
    229268
    230     do ks=1,nspec                                  ! loop over species
     269    do ks=1,nspec            ! loop over species
    231270      wetdeposit(ks)=0.
     271
     272   
     273    !conflicting changes to the same routine: 1=NIK 2 =PS
     274    scheme_number=2   
     275    if (scheme_number.eq.1) then !NIK
     276
    232277      if (weta(ks).gt.0.) then
    233278        if (clouds_v.ge.4) then
    … …  
    236281          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    237282  !        write(*,*) 'bel. wetscav: ',wetscav
     283
    238284        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    239285  !        IN CLOUD SCAVENGING
    … …  
    245291             act_temp=tt(ix,jy,hz,n)
    246292          endif
    247           cl=2E-7*prec**0.36
     293
     294! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
     295! weta_in=2.0E-07 (default)
     296! wetb_in=0.36 (default)
     297! wetc_in=0.9 (default)
     298! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
     299          cl=weta_in(ks)*prec**wetb_in(ks)
    248300          if (dquer(ks).gt.0) then ! is particle
    249             S_i=0.9/cl
     301            S_i=wetc_in(ks)/cl
    250302           else ! is gas
    251303            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    252304            S_i=1/cle
    253305           endif
    254            wetscav=S_i*prec/3.6E6/clouds_h
     306           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
    255307  !         write(*,*) 'in. wetscav:'
    256308  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
    … …  
    289341         wetdeposit(ks)=0.
    290342      endif ! weta(k)
     343
     344     elseif (scheme_number.eq.2) then ! PS
     345
     346!PS          indcloud=0 ! Use this for FOR TESTING,
     347!PS          will skip the new in/below cloud method !!!             
     348
     349          if (weta(ks).gt.0.) then
     350            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
     351!C               for aerosols and not highliy soluble substances weta=5E-6
     352              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
     353!c             write(*,*) 'bel. wetscav: ',wetscav
     354            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
     355              if (ngrid.gt.0) then
     356                 act_temp=ttn(ix,jy,kz,n,ngrid)
     357              else
     358                 act_temp=tt(ix,jy,kz,n)
     359              endif
     360
     361! from NIK             
     362! weta_in=2.0E-07 (default)
     363! wetb_in=0.36 (default)
     364! wetc_in=0.9 (default)
     365
     366
     367              cl=2E-7*prec**0.36
     368              if (dquer(ks).gt.0) then ! is particle
     369                S_i=0.9/cl
     370              else ! is gas
     371                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
     372                S_i=1/cle
     373              endif
     374              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
     375            else ! PS: no cloud diagnosed, old scheme,
     376!CPS          using with fixed a,b for simplicity, one may wish to change!!
     377              wetscav = 1.e-4*prec**0.62
     378            endif
     379
     380
     381            wetdeposit(ks)=xmass1(jpart,ks)*    &
     382            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
     383             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS)
     384            restmass = xmass1(jpart,ks)-wetdeposit(ks)
     385            if (ioutputforeachrelease.eq.1) then
     386              kp=npoint(jpart)
     387            else
     388              kp=1
     389            endif
     390            if (restmass .gt. smallnum) then
     391              xmass1(jpart,ks)=restmass
     392!cccccccccccccccc depostatistic
     393!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
     394!cccccccccccccccc depostatistic
     395            else
     396              xmass1(jpart,ks)=0.
     397            endif
     398!C Correct deposited mass to the last time step when radioactive decay of
     399!C gridded deposited mass was calculated
     400            if (decay(ks).gt.0.) then
     401              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
     402            endif
     403          else  ! weta(k)<0
     404             wetdeposit(ks)=0.
     405          endif
     406
     407     endif !on scheme
     408
     409
     410
    291411    end do
    292412
    … …  
    295415  !*****************************************************************************
    296416
    297     if (ldirect.eq.1) then
    298     call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
    299          real(ytra1(jpart)),nage,kp)
    300     if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
    301          wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
    302          nage,kp)
    303     endif
     417!    if (ldirect.eq.1) then
     418!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
     419!         real(ytra1(jpart)),nage,kp)
     420!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
     421!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
     422!         nage,kp)
     423!    endif
     424
     425   !PS
     426        if (ldirect.eq.1) then
     427          call wetdepokernel(nclass(jpart),wetdeposit,     &
     428             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
     429          if (nested_output.eq.1)                          &
     430           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
     431             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp) 
     432        endif
     433
    304434
    30543520  continue
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG