Ignore:
Files:
12 added
11 deleted
24 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/FLEXPART.f90

    r14 r14  
    4949  integer :: i,j,ix,jy,inest
    5050  integer :: idummy = -320
    51   character(len=256) :: pathfile
    5251
    5352  ! Generate a large number of random numbers
    … …  
    6160  ! Print the GPL License statement
    6261  !*******************************************************
    63   ! print*,'Welcome to FLEXPART Version 9.1 (Build 20121029)'
    64   print*,'Welcome to FLEXPART Version 9.1.3 (Build 2013158)'
     62  print*,'Welcome to FLEXPART Version 9.0'
    6563  print*,'FLEXPART is free software released under the GNU Genera'// &
    6664       'l Public License.'
    … …  
    6967  !*******************************************************
    7068
    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)
     69  call readpaths
    8470
    8571  ! Read the user specifications for the current model run

  • trunk/src/com_mod.f90

    • Property svn:executable set to *
    r14 r14  
    77!        June 1996                                                             *
    88!                                                                              *
    9 !        Last update:15 August 2013 IP                                         *
     9!        Last update: 9 August 2000                                            *
    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)
    143141  real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
    144142  real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
    … …  
    174172
    175173  ! WET DEPOSITION
    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
     174  ! weta, wetb              parameters for determining wet scavenging coefficients
    179175
    180176  ! GAS DEPOSITION
    … …  
    335331  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    336332  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
    337   !scavenging NIK, PS
    338333  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    339334  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 
    343335
    344336  ! uu,vv,ww [m/2]       wind components in x,y and z direction
    … …  
    354346  !      rainout  conv/lsp dominated  2/3
    355347  !      washout  conv/lsp dominated  4/5
    356 ! PS 2013
    357 !c icloudbot (m)        cloud bottom height
    358 !c icloudthck (m)       cloud thickness     
    359 
    360348  ! pplev for the GFS version
    361349

  • trunk/src/concoutput.f90

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

  • trunk/src/erf.f90

    • Property svn:executable set to *
    r14 r14  
    4242
    4343  integer :: j
    44   real(kind=dp) :: x,tmp,ser,xx,gammln
    45   real(KIND=dp) :: cof(6) = (/ &
     44  real :: x,tmp,ser,xx,gammln
     45  real :: cof(6) = (/ &
    4646       76.18009173_dp, -86.50532033_dp, 24.01409822_dp, &
    4747       -1.231739516_dp, .120858003e-2_dp, -.536382e-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
     48  real :: stp = 2.50662827465_dp
     49  real :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
    5050
    5151  x=xx-one
    … …  
    6262function gammp(a,x)
    6363
    64   use par_mod, only: dp
    65 
    66   implicit none
    67 
    68   real(KIND=dp) :: a, x, gln, gamser, gammp, gammcf
     64  implicit none
     65
     66  real :: a, x, gln, gamser, gammp, gammcf
    6967
    7068  if(x .lt. 0. .or. a .le. 0.) then
    … …  
    8381function gammq(a,x)
    8482
    85   use par_mod, only: dp
    86 
    87   implicit none
    88 
    89   real(KIND=dp) :: a, x, gln, gamser, gammq, gammcf
     83  implicit none
     84
     85  real :: a, x, gln, gamser, gammq, gammcf
    9086
    9187  if(x.lt.0..or.a.le.0.) then
    … …  
    104100subroutine gser(gamser,a,x,gln)
    105101
    106   use par_mod, only: dp
    107 
    108102  implicit none
    109103
    110104  integer :: n
    111   real(KIND=dp) :: gamser, a, x, gln, ap, summ, del
    112   real(KIND=dp), external :: gammln
     105  real :: gamser, a, x, gln, ap, summ, del
     106  real, external :: gammln
    113107
    114108  integer,parameter :: itmax=100
    … …  
    140134subroutine gcf(gammcf,a,x,gln)
    141135
    142   use par_mod, only: dp
    143 
    144136  implicit none
    145137
    146138  integer :: n
    147   real(KIND=4) :: gammcf, x, a, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
    148   real(KIND=dp), external :: gammln
     139  real :: gammcf, a, x, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
     140  real, external :: gammln
    149141
    150142  integer,parameter :: itmax=100
    … …  
    180172function erf(x)
    181173
    182   use par_mod, only: dp
    183 
    184   implicit none
    185 
    186   real(KIND=dp) :: x, erf
    187   real(KIND=dp), external :: gammp
     174  implicit none
     175
     176  real :: x, erf
     177  real, external :: gammp
    188178
    189179  if(x.lt.0.)then
    … …  
    196186function erfc(x)
    197187
    198   use par_mod, only: dp
    199 
    200   implicit none
    201 
    202   real(KIND=dp) :: x, erfc
    203   real(KIND=dp), external :: gammp, gammq
     188  implicit none
     189
     190  real :: x, erfc
     191  real, external :: gammp, gammq
    204192
    205193  if(x.lt.0.)then
    … …  
    212200function erfcc(x)
    213201
    214   use par_mod, only: dp
    215 
    216   implicit none
    217 
    218   real(KIND=dp) :: x, z, t, erfcc
     202  implicit none
     203
     204  real :: x, z, t, erfcc
    219205
    220206  z=abs(x)

  • trunk/src/getvdep.f90

    • Property svn:executable set to *

  • trunk/src/gridcheck.f90

    • Property svn:executable set to *
    r14 r14  
    9898
    9999  integer :: isec1(56),isec2(22+nxmax+nymax)
    100   !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
    101   real(kind=4) :: zsec2(184),zsec4(jpunp)
     100  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
    102101  character(len=1) :: opt
    103102
    … …  
    467466    k=nlev_ec+1+numskip+i
    468467    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

  • trunk/src/gridcheck_gfs.f90

    • Property svn:executable set to *

  • trunk/src/gridcheck_nests.f90

    • Property svn:executable set to *

  • trunk/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
    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 
     22subroutine 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
    3326  !****************************************************************************
    3427  !                                                                           *
    … …  
    4740  !                                                                           *
    4841  !     30 August 1996                                                        *
    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  *
     42  !                                                                           *
    5343  !****************************************************************************
    5444  !                                                                           *
    … …  
    8070
    8171  integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
    82   !integer :: itime,itime1,itime2,level,indexh
    83   integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
    84   integer :: intiy1,intiy2,ipsum,icmv 
     72  integer :: itime,itime1,itime2,level,indexh
    8573  real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
    8674  real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
    8775  real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
    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
     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
    9378
    9479
    … …  
    142127         + p3*yy3(ix ,jyp,level,indexh) &
    143128         + 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 
    183129  end do
    184130
    … …  
    187133  ! 2.) Temporal interpolation (linear)
    188134  !************************************
    189 
    190       if (itime .lt. itime1) then
    191         print*,'interpol_rain.f90'
    192         print*,itime,itime1,itime2
    193         stop 'ITIME PROBLEM'
    194       endif
    195 
    196135
    197136  dt1=real(itime-itime1)
    … …  
    204143
    205144
    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 
    223145end subroutine interpol_rain

  • trunk/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)
    1112FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
    1213#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)

  • trunk/src/par_mod.f90

    • Property svn:executable set to *
    r14 r14  
    2828!        1997                                                                  *
    2929!                                                                              *
    30 !        Last update 15 August 2013 IP                                         *
     30!        Last update 10 August 2000                                            *
    3131!                                                                              *
    3232!*******************************************************************************
    … …  
    122122
    123123  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
    124   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
     124  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
    125125  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    126   !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
    127   integer,parameter :: nxshift=359 ! for ECMWF
     126  integer,parameter :: nxshift=359  ! for ECMWF
    128127  !integer,parameter :: nxshift=0     ! for GFS
    129128
    … …  
    151150  !*********************************************
    152151
    153   !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    154   integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
     152  integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
     153  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
    155154
    156155  ! maxnests                maximum number of nested grids
    … …  
    195194  !**************************************************
    196195
    197   integer,parameter :: maxpart=150000
    198   integer,parameter :: maxspec=1
     196  integer,parameter :: maxpart=2000000
     197  integer,parameter :: maxspec=4
    199198
    200199
    … …  
    259258  integer,parameter :: unitboundcond=89
    260259
    261 !******************************************************
    262 ! integer code for missing values, used in wet scavenging (PS, 2012)
    263 !******************************************************
    264 
    265       integer icmv
    266       parameter(icmv=-9999)
    267 
    268260end module par_mod

  • trunk/src/partdep.f90

    • Property svn:executable set to *

  • trunk/src/readageclasses.f90

    • Property svn:executable set to *

  • trunk/src/readavailable.f90

    • Property svn:executable set to *
    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))
    127125    open(unitavailab,file=path(numpath+2*(k-1)+2) &
    128126         (1:length(numpath+2*(k-1)+2)),status='old',err=998)
    … …  
    277275  return
    278276
    279 998   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE   #### '
     277998   write(*,*) ' #### FLEXPART MODEL ERROR! FILE   #### '
    280278  write(*,'(a)') '     '//path(numpath+2*(k-1)+2) &
    281279       (1:length(numpath+2*(k-1)+2))
    … …  
    283281  stop
    284282
    285 999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
     283999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
    286284  write(*,'(a)') '     '//path(4)(1:length(4))
    287285  write(*,*) ' #### CANNOT BE OPENED           #### '

  • trunk/src/readcommand.f90

    r14 r14  
    8080  character(len=50) :: line
    8181  logical :: old
    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
     82
    13883
    13984  ! Open the command file and read user options
    140   ! Namelist input first: try to read as namelist file
     85  !********************************************
     86
     87
     88  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'
    14194  !**************************************************************************
    142   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    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)
     95
     96  call skplin(9,unitcommand)
     97  read (unitcommand,901) line
     98901   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
    153159  close(unitcommand)
    154160
    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 
    240161  ifine=max(ifine,1)
     162
    241163
    242164  ! Determine how Markov chain is formulated (for w or for w/sigw)

  • trunk/src/readdepo.f90

    • Property svn:executable set to *

  • trunk/src/readpartpositions.f90

    r14 r14  
    7777
    7878  read(unitpartin) numpointin
    79   if (numpointin.ne.numpoint) goto 995
     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
    8086  do i=1,numpointin
    8187    read(unitpartin)
    … …  
    146152  stop
    147153
    148 995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
    149   write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
    150   write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
    151   stop
    152154
    153155996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '

  • trunk/src/readpaths.f90

    • Property svn:executable set to *
    r14 r14  
    2020!**********************************************************************
    2121
    22 subroutine readpaths(pathfile)
     22subroutine readpaths
    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                   *
    3532  !                                                                            *
    3633  !*****************************************************************************
    … …  
    5047  implicit none
    5148
    52   integer   :: i
    53   character(256) :: pathfile
     49  integer :: i
    5450
    5551  ! Read the pathname information stored in unitpath
    5652  !*************************************************
    5753
    58   open(unitpath,file=trim(pathfile),status='old',err=999)
     54
     55  open(unitpath,file='pathnames',status='old',err=999)
    5956
    6057  do i=1,numpath
    … …  
    7370    length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
    7471  end do
    75   print*,length(5),length(6)
    7672
    7773

  • trunk/src/readreleases.f90

    • Property svn:executable set to *
    r14 r14  
    3333  !     Update: 29 January 2001                                                *
    3434  !     Release altitude can be either in magl or masl                         *
    35   !     Update: Aug 2013 IP                                                                        *
     35  !                                                                            *
    3636  !*****************************************************************************
    3737  !                                                                            *
    … …  
    5555  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
    5656  !                     area                                                   *
    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)    * 
     57  ! weta, wetb          parameters to determine the wet scavenging coefficient *
    6158  ! 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                     *
    6459  !                                                                            *
    6560  !*****************************************************************************
    … …  
    7368
    7469  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
    75   integer,parameter :: num_min_discrete=100
    7670  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    77   real(kind=dp) :: jul1,jul2,julm,juldate
     71  real(kind=dp) :: jul1,jul2,juldate
    7872  character(len=50) :: line
    7973  logical :: old
    … …  
    283277    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
    284278
    285 
    286 
    287279  ! Check if wet deposition or OH reaction shall be calculated
    288280  !***********************************************************
    289281    if (weta(i).gt.0.)  then
    290282      WETDEP=.true.
    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 
     283      write (*,*) 'Wetdeposition switched on: ',weta(i),i
     284    endif
    301285    if (ohreact(i).gt.0) then
    302286      OHREA=.true.
    … …  
    392376  jul1=juldate(id1,it1)
    393377  jul2=juldate(id2,it2)
    394   julm=(jul1+jul2)/2.
    395378  if (jul1.gt.jul2) then
    396379    write(*,*) 'FLEXPART MODEL ERROR'
    … …  
    408391        stop
    409392      endif
    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
     393      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     394      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
    417395    else if (ldirect.eq.-1) then
    418396      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
    … …  
    423401        stop
    424402      endif
    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
     403      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     404      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
    432405    endif
    433406  endif

  • trunk/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   !                                                                            *
    3633  !*****************************************************************************
    3734  !                                                                            *
    … …  
    3936  ! decaytime(maxtable)  half time for radiological decay                      *
    4037  ! specname(maxtable)   names of chemical species, radionuclides              *
    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         *
     38  ! wetscava, wetscavb   Parameters for determining scavenging coefficient     *
    4639  ! ohreact              OH reaction rate                                      *
    4740  ! id_spec              SPECIES number as referenced in RELEASE file          *
    … …  
    8578    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
    8679  !  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 
    9880    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
    9981  !  write(*,*) reldiff(pos_spec)

  • trunk/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
    108103
    109104  hflswitch=.false.
    … …  
    370365  endif
    371366
    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 
    388367  ! For global fields, assign the leftmost data column also to the rightmost
    389368  ! data column; if required, shift whole grid by nxshift grid points

  • trunk/src/readwind_gfs.f90

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

  • trunk/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
    5351  !*****************************************************************************
    5452  !                                                                            *
    … …  
    7270
    7371  integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
    74  integer :: rain_cloud_above,kz_inv !SE
    75   integer icloudtop !PS
     72  integer :: rain_cloud_above,kz_inv
    7673  real :: f_qvsat,pressure
    77  !real :: rh,lsp,convp
    78   real :: rh,lsp,convp,prec,rhmin 
     74  real :: rh,lsp,convp
    7975  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
    8076  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
    8177  real :: xlon,ylat,xlonr,dzdx,dzdy
    82   real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
     78  real :: dzdx1,dzdx2,dzdy1,dzdy2
    8379  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
    8480  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    … …  
    8783  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8884  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
    89   logical lconvectprec
    9085  real,parameter :: const=r_air/ga
    91   parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
    9286
    9387  logical :: init = .true.
    … …  
    551545  !   create a cloud and rainout/washout field, clouds occur where rh>80%
    552546  !   total cloudheight is stored at level 0
    553 
    554 
    555 
    556547  do jy=0,nymin1
    557548    do ix=0,nxmin1
    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
    611 98        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
     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
    620570            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
     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
    643578            endif
    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 
     579         endif
     580      end do
    657581    end do
    658582  end do
    … …  
    681605  !104   continue
    682606  ! close(4)
    683 
    684 
    685607end subroutine verttransform

  • trunk/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
    4441  !*****************************************************************************
    4542  !                                                                            *
    … …  
    7471
    7572  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
    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
     73  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
     74  integer :: ks, kp
    8075  real :: S_i, act_temp, cl, cle ! in cloud scavenging
    8176  real :: clouds_h ! cloud height for the specific grid point
    82   real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
     77  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
    8378  real :: wetdeposit(maxspec),restmass
    8479  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    … …  
    151146    interp_time=nint(itime-0.5*ltsample)
    152147
    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 
     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
    181159
    182160  ! get the level were the actual particle is in
    183161      do il=2,nz
    184162        if (height(il).gt.ztra1(jpart)) then
    185           !hz=il-1
    186           kz=il-1
     163          hz=il-1
    187164          goto 26
    188165        endif
    … …  
    197174  ! scavenging is done
    198175
    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 
     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'
    226187
    227188  ! 1) Parameterization of the the area fraction of the grid cell where the
    … …  
    267228  !**********************************************************
    268229
    269     do ks=1,nspec            ! loop over species
     230    do ks=1,nspec                                  ! loop over species
    270231      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 
    277232      if (weta(ks).gt.0.) then
    278233        if (clouds_v.ge.4) then
    … …  
    281236          wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
    282237  !        write(*,*) 'bel. wetscav: ',wetscav
    283 
    284238        else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
    285239  !        IN CLOUD SCAVENGING
    … …  
    291245             act_temp=tt(ix,jy,hz,n)
    292246          endif
    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)
     247          cl=2E-7*prec**0.36
    300248          if (dquer(ks).gt.0) then ! is particle
    301             S_i=wetc_in(ks)/cl
     249            S_i=0.9/cl
    302250           else ! is gas
    303251            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    304252            S_i=1/cle
    305253           endif
    306            wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
     254           wetscav=S_i*prec/3.6E6/clouds_h
    307255  !         write(*,*) 'in. wetscav:'
    308256  !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
    … …  
    341289         wetdeposit(ks)=0.
    342290      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 
    411291    end do
    412292
    … …  
    415295  !*****************************************************************************
    416296
    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 
     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
    434304
    43530520  continue
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG