Changeset 20 for trunk


Ignore:
Timestamp:
Dec 23, 2013, 6:23:38 PM (5 years ago)
Author:
igpis
Message:

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

Location:
trunk/src
Files:
7 added
7 deleted
21 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/FLEXPART.f90

    r4 r20  
    4949  integer :: i,j,ix,jy,inest
    5050  integer :: idummy = -320
     51  character(len=256) :: inline_options  !pathfile, flexversion, arg2
     52
    5153
    5254  ! Generate a large number of random numbers
     
    5860  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
    5961
     62  !
     63  flexversion='Version 9.1.8  (2013-12-08)'
     64  !verbosity=0
     65  ! Read the pathnames where input/output files are stored
     66  !*******************************************************
     67
     68  inline_options='none'
     69  select case (iargc())
     70  case (2)
     71    call getarg(1,arg1)
     72    pathfile=arg1
     73    call getarg(2,arg2)
     74    inline_options=arg2
     75  case (1)
     76    call getarg(1,arg1)
     77    pathfile=arg1
     78    verbosity=0
     79    if (arg1(1:1).eq.'-') then
     80        write(pathfile,'(a11)') './pathnames'
     81        inline_options=arg1
     82    endif
     83  case (0)
     84    write(pathfile,'(a11)') './pathnames'
     85    verbosity=0
     86  end select
     87 
     88    if (inline_options(1:1).eq.'-') then
     89      print*, 'inline options=', inline_options       
     90      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
     91         print*, 'verbose mode 1: additional information will be displayed'
     92         verbosity=1
     93      endif
     94      if (trim(inline_options).eq.'-v2') then
     95         print*, 'verbose mode 2: additional information will be displayed'
     96         verbosity=2
     97      endif
     98      if (trim(inline_options).eq.'-i') then
     99         print*, 'info mode: will provide run specific information and stop'
     100         verbosity=1
     101         info_flag=1
     102      endif
     103      if (trim(inline_options).eq.'-i2') then
     104         print*, 'info mode: will provide run specific information and stop'
     105         verbosity=2
     106         info_flag=1
     107      endif
     108    endif
     109
     110
    60111  ! Print the GPL License statement
    61112  !*******************************************************
    62   print*,'Welcome to FLEXPART Version 9.0'
     113  print*,'Welcome to FLEXPART', trim(flexversion)
    63114  print*,'FLEXPART is free software released under the GNU Genera'// &
    64115       'l Public License.'
    65 
    66   ! Read the pathnames where input/output files are stored
    67   !*******************************************************
    68 
    69   call readpaths
     116           
     117  if (verbosity.gt.0) then
     118        WRITE(*,*) 'call readpaths'
     119  endif
     120  call readpaths(pathfile)
     121   
     122 
     123  if (verbosity.gt.1) then !show clock info
     124     !print*,'length(4)',length(4)
     125     !count=0,count_rate=1000
     126     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
     127     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
     128     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
     129     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
     130     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
     131  endif
     132
    70133
    71134  ! Read the user specifications for the current model run
    72135  !*******************************************************
    73136
     137  if (verbosity.gt.0) then
     138        WRITE(*,*) 'call readcommand'
     139  endif
    74140  call readcommand
     141  if (verbosity.gt.0) then
     142        WRITE(*,*) '    ldirect=', ldirect
     143        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
     144        WRITE(*,*) '    iedate,ietime=', iedate,ietime
     145        if (verbosity.gt.1) then   
     146                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     147                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     148        endif     
     149  endif
    75150
    76151
    77152  ! Read the age classes to be used
    78153  !********************************
    79 
     154  if (verbosity.gt.0) then
     155        WRITE(*,*) 'call readageclasses'
     156  endif
    80157  call readageclasses
     158
     159  if (verbosity.gt.1) then   
     160                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     161                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     162  endif     
     163 
    81164
    82165
     
    84167  !******************************************************************
    85168
     169  if (verbosity.gt.0) then
     170        WRITE(*,*) 'call readavailable'
     171  endif 
    86172  call readavailable
    87 
    88173
    89174  ! Read the model grid specifications,
    90175  ! both for the mother domain and eventual nests
    91176  !**********************************************
     177 
     178  if (verbosity.gt.0) then
     179     WRITE(*,*) 'call gridcheck'
     180  endif
    92181
    93182  call gridcheck
     183
     184  if (verbosity.gt.1) then   
     185     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     186     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     187  endif     
     188 
     189
     190  if (verbosity.gt.0) then
     191        WRITE(*,*) 'call gridcheck_nests'
     192  endif 
    94193  call gridcheck_nests
    95194
     
    98197  !************************************
    99198
     199  if (verbosity.gt.0) then
     200        WRITE(*,*) 'call readoutgrid'
     201  endif
     202
    100203  call readoutgrid
    101   if (nested_output.eq.1) call readoutgrid_nest
    102 
     204
     205  if (nested_output.eq.1) then
     206          call readoutgrid_nest
     207    if (verbosity.gt.0) then
     208        WRITE(*,*) '# readoutgrid_nest'
     209    endif
     210  endif
    103211
    104212  ! Read the receptor points for which extra concentrations are to be calculated
    105213  !*****************************************************************************
    106214
     215  if (verbosity.eq.1) then
     216     print*,'call readreceptors'
     217  endif
    107218  call readreceptors
    108 
    109219
    110220  ! Read the physico-chemical species property table
     
    117227  !***************************
    118228
     229  if (verbosity.gt.0) then
     230    print*,'call readlanduse'
     231  endif
    119232  call readlanduse
    120233
     
    123236  !********************************************************************
    124237
     238  if (verbosity.gt.0) then
     239    print*,'call assignland'
     240  endif
    125241  call assignland
    126242
     
    130246  !**********************************************
    131247
     248  if (verbosity.gt.0) then
     249    print*,'call readreleases'
     250  endif
    132251  call readreleases
     252
    133253
    134254  ! Read and compute surface resistances to dry deposition of gases
    135255  !****************************************************************
    136256
     257  if (verbosity.gt.0) then
     258    print*,'call readdepo'
     259  endif
    137260  call readdepo
    138 
    139261
    140262  ! Convert the release point coordinates from geografical to grid coordinates
    141263  !***************************************************************************
    142264
    143   call coordtrafo
     265  call coordtrafo 
     266  if (verbosity.gt.0) then
     267    print*,'call coordtrafo'
     268  endif
    144269
    145270
     
    147272  !*****************************************
    148273
     274  if (verbosity.gt.0) then
     275    print*,'Initialize all particles to non-existent'
     276  endif
    149277  do j=1,maxpart
    150278    itra1(j)=-999999999
     
    155283
    156284  if (ipin.eq.1) then
     285    if (verbosity.gt.0) then
     286      print*,'call readpartpositions'
     287    endif
    157288    call readpartpositions
    158289  else
     290    if (verbosity.gt.0) then
     291      print*,'numpart=0, numparticlecount=0'
     292    endif   
    159293    numpart=0
    160294    numparticlecount=0
     
    166300  !***************************************************************
    167301
     302
     303  if (verbosity.gt.0) then
     304    print*,'call outgrid_init'
     305  endif
    168306  call outgrid_init
    169307  if (nested_output.eq.1) call outgrid_init_nest
     
    173311  !******************
    174312
    175   if (OHREA.eqv..TRUE.) &
     313  if (OHREA.eqv..TRUE.) then
     314    if (verbosity.gt.0) then
     315       print*,'call readOHfield'
     316    endif
    176317       call readOHfield
     318  endif
    177319
    178320  ! Write basic information on the simulation to a file "header"
     
    180322  !******************************************************************
    181323
     324
     325  if (verbosity.gt.0) then
     326    print*,'call writeheader'
     327  endif
     328
    182329  call writeheader
    183   if (nested_output.eq.1) call writeheader_nest
    184   open(unitdates,file=path(2)(1:length(2))//'dates')
     330  ! FLEXPART 9.2 ticket ?? write header in ASCII format
     331  call writeheader_txt
     332  !if (nested_output.eq.1) call writeheader_nest
     333  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
     334
     335  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
     336  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
     337
     338
     339
     340  !open(unitdates,file=path(2)(1:length(2))//'dates')
     341
     342  if (verbosity.gt.0) then
     343    print*,'call openreceptors'
     344  endif
    185345  call openreceptors
    186346  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
     
    190350  !***************************************************************************
    191351
     352  if (verbosity.gt.0) then
     353    print*,'discretize release times'
     354  endif
    192355  do i=1,numpoint
    193356    ireleasestart(i)=nint(real(ireleasestart(i))/ &
     
    200363  ! Initialize cloud-base mass fluxes for the convection scheme
    201364  !************************************************************
     365
     366  if (verbosity.gt.0) then
     367    print*,'Initialize cloud-base mass fluxes for the convection scheme'
     368  endif
    202369
    203370  do jy=0,nymin1
     
    218385  !********************************
    219386
     387  if (verbosity.gt.0) then
     388     if (verbosity.gt.1) then   
     389       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     390       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     391     endif
     392     if (info_flag.eq.1) then
     393         print*, 'info only mode (stop)'   
     394         stop
     395     endif
     396     print*,'call timemanager'
     397  endif
     398
    220399  call timemanager
    221400
  • trunk/src/com_mod.f90

    r4 r20  
    77!        June 1996                                                             *
    88!                                                                              *
    9 !        Last update: 9 August 2000                                            *
     9!        Last update:15 August 2013 IP                                         *
    1010!                                                                              *
    1111!*******************************************************************************
     
    2525  character :: path(numpath+2*maxnests)*120
    2626  integer :: length(numpath+2*maxnests)
    27 
     27  character(len=256) :: pathfile, flexversion, arg1, arg2
     28 
    2829  ! path                    path names needed for trajectory model
    2930  ! length                  length of path names needed for trajectory model
    30 
     31  ! pathfile                file where pathnames are stored
     32  ! flexversion             version of flexpart
     33  ! arg                     input arguments from launch at command line
    3134
    3235  !********************************************************
     
    6568  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
    6669  integer :: mquasilag,nested_output,ind_source,ind_receptor
    67   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond
     70  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
    6871  logical :: turbswitch
    6972
     
    9194  ! mquasilag 0: normal run
    9295  !      1: Particle position output is produced in a condensed format and particles are numbered
     96  ! surf_only   switch output in grid_time files for surface only or full vertical resolution
     97  !      0=no (full vertical resolution), 1=yes (surface only)
    9398  ! nested_output: 0 no, 1 yes
    9499  ! turbswitch              determines how the Markov chain is formulated
     
    139144  real :: decay(maxspec)
    140145  real :: weta(maxspec),wetb(maxspec)
     146! NIK: 31.01.2013- parameters for in-cloud scavening
     147  real :: weta_in(maxspec), wetb_in(maxspec), wetc_in(maxspec), wetd_in(maxspec)
    141148  real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
    142149  real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
     
    172179
    173180  ! WET DEPOSITION
    174   ! weta, wetb              parameters for determining wet scavenging coefficients
     181  ! weta, wetb              parameters for determining below-cloud wet scavenging coefficients
     182  ! weta_in, wetb_in       parameters for determining in-cloud wet scavenging coefficients
     183  ! wetc_in, wetd_in       parameters for determining in-cloud wet scavenging coefficients
    175184
    176185  ! GAS DEPOSITION
     
    331340  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
    332341  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
     342  !scavenging NIK, PS
    333343  integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
    334344  integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
     345      integer icloudbot(0:nxmax-1,0:nymax-1,2)
     346      integer icloudthck(0:nxmax-1,0:nymax-1,2)
     347
    335348
    336349  ! uu,vv,ww [m/2]       wind components in x,y and z direction
     
    346359  !      rainout  conv/lsp dominated  2/3
    347360  !      washout  conv/lsp dominated  4/5
     361! PS 2013
     362!c icloudbot (m)        cloud bottom height
     363!c icloudthck (m)       cloud thickness     
     364
    348365  ! pplev for the GFS version
    349366
     
    662679
    663680
    664 
    665681  !********************
    666682  ! Random number field
     
    671687  ! rannumb                 field of normally distributed random numbers
    672688
     689  !********************
     690  ! Verbosity, testing flags
     691  !********************   
     692  integer :: verbosity=0
     693  integer :: info_flag=0
     694  INTEGER :: count_clock, count_clock0,  count_rate, count_max
     695   
    673696
    674697end module com_mod
  • trunk/src/concoutput.f90

    r4 r20  
    108108  write(adate,'(i8.8)') jjjjmmdd
    109109  write(atime,'(i6.6)') ihmmss
     110  open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
    110111  write(unitdates,'(a)') adate//atime
    111 
     112  close(unitdates) 
    112113
    113114  ! For forward simulations, output fields have dimension MAXSPEC,
     
    158159        yl=outlat0+real(jy)*dyout
    159160        xl=(xl-xlon0)/dx
    160         yl=(yl-ylat0)/dx
     161        yl=(yl-ylat0)/dy !v9.1.1
    161162        iix=max(min(nint(xl),nxmin1),0)
    162163        jjy=max(min(nint(yl),nymin1),0)
  • trunk/src/gridcheck.f90

    r4 r20  
    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
  • trunk/src/interpol_rain.f90

    r4 r20  
    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 (abs(itime) .lt. abs(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
  • trunk/src/openouttraj.f90

    r4 r20  
    5454
    5555  if (ldirect.eq.1) then
    56   write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime,'FLEXPART V8.2'
     56  write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime, trim(flexversion)
    5757  else
    58   write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime,'FLEXPART V8.2'
     58  write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime, trim(flexversion)
    5959  endif
    6060  write(unitouttraj,*) method,lsubgrid,lconvection
  • trunk/src/outg_mod.f90

    r4 r20  
    4343  real,allocatable, dimension (:,:) :: wetgridsigma
    4444  real,allocatable, dimension (:) :: sparse_dump_r
     45  real,allocatable, dimension (:) :: sparse_dump_u
    4546  integer,allocatable, dimension (:) :: sparse_dump_i
    4647
  • trunk/src/outgrid_init.f90

    r4 r20  
    247247       max(numygrid,numygridn)*numzgrid),stat=stat)
    248248    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
     249
     250   allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
     251       max(numygrid,numygridn)*numzgrid),stat=stat)
     252        if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
     253
    249254  allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
    250255       max(numygrid,numygridn)*numzgrid),stat=stat)
  • trunk/src/par_mod.f90

    r4 r20  
    2828!        1997                                                                  *
    2929!                                                                              *
    30 !        Last update 10 August 2000                                            *
     30!        Last update 15 August 2013 IP                                         *
    3131!                                                                              *
    3232!*******************************************************************************
     
    120120  ! Maximum dimensions of the input mother grids
    121121  !*********************************************
    122 
    123   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
     122 
     123  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
     124  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
     125  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    125126  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
    126   integer,parameter :: nxshift=359  ! for ECMWF
    127   !integer,parameter :: nxshift=0     ! for GFS
     127  !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
     128
     129  integer,parameter :: nxshift=359 ! for ECMWF
     130  !integer,parameter :: nxshift=0     ! for GFS or FNL (XF)
    128131
    129132  integer,parameter :: nconvlevmax = nuvzmax-1
     
    150153  !*********************************************
    151154
    152   integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
    153   !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
    154 
     155  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
     156  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
     157  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
    155158  ! maxnests                maximum number of nested grids
    156159  ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
     
    194197  !**************************************************
    195198
    196   integer,parameter :: maxpart=2000000
    197   integer,parameter :: maxspec=4
     199  !integer,parameter :: maxpart=4000000
     200  integer,parameter :: maxpart=2000
     201  integer,parameter :: maxspec=6
    198202
    199203
     
    255259  integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
    256260  integer,parameter :: unitOH=1
    257   integer,parameter :: unitdates=94, unitheader=90, unitshortpart=95
     261  integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
    258262  integer,parameter :: unitboundcond=89
    259263
     264!******************************************************
     265! integer code for missing values, used in wet scavenging (PS, 2012)
     266!******************************************************
     267
     268    !  integer icmv
     269    !  parameter(icmv=-9999)
     270      integer,parameter ::  icmv=-9999
     271
     272! Parameters for testing
     273!*******************************************
     274!  integer :: verbosity=0
     275
    260276end module par_mod
  • trunk/src/readavailable.f90

    r4 r20  
    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           #### '
  • trunk/src/readcommand.f90

    r4 r20  
    8080  character(len=50) :: line
    8181  logical :: old
    82 
     82  logical :: nmlout=.true. !.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    surf_only   
     111
     112  ! Presetting namelist command
     113  ldirect=1
     114  ibdate=20000101
     115  ibtime=0
     116  iedate=20000102
     117  ietime=0
     118  loutstep=10800
     119  loutaver=10800
     120  loutsample=900
     121  itsplit=999999999
     122  lsynctime=900
     123  ctl=-5.0
     124  ifine=4
     125  iout=3
     126  ipout=0
     127  lsubgrid=1
     128  lconvection=1
     129  lagespectra=0
     130  ipin=1
     131  ioutputforeachrelease=1
     132  iflux=1
     133  mdomainfill=0
     134  ind_source=1
     135  ind_receptor=1
     136  mquasilag=0
     137  nested_output=0
     138  linit_cond=0
     139  surf_only=0
    83140
    84141  ! Open the command file and read user options
    85   !********************************************
    86 
    87 
     142  ! Namelist input first: try to read as namelist file
     143  !**************************************************************************
    88144  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
     145         form='formatted',iostat=readerror)
     146  ! If fail, check if file does not exist
     147  if (readerror.ne.0) then
     148   
     149    print*,'***ERROR: file COMMAND not found in '
     150    print*, path(1)(1:length(1))//'COMMAND'
     151    print*, 'Check your pathnames file.'
     152    stop
     153
     154  endif
     155
     156  read(unitcommand,command,iostat=readerror)
    159157  close(unitcommand)
    160158
     159  ! If error in namelist format, try to open with old input code
     160  if (readerror.ne.0) then
     161
     162    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
     163         err=999)
     164
     165    ! Check the format of the COMMAND file (either in free format,
     166    ! or using formatted mask)
     167    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
     168    !**************************************************************************
     169
     170    call skplin(9,unitcommand)
     171    read (unitcommand,901) line
     172  901   format (a)
     173    if (index(line,'LDIRECT') .eq. 0) then
     174      old = .false.
     175    else
     176      old = .true.
     177    endif
     178    rewind(unitcommand)
     179
     180    ! Read parameters
     181    !****************
     182
     183    call skplin(7,unitcommand)
     184    if (old) call skplin(1,unitcommand)
     185
     186    read(unitcommand,*) ldirect
     187    if (old) call skplin(3,unitcommand)
     188    read(unitcommand,*) ibdate,ibtime
     189    if (old) call skplin(3,unitcommand)
     190    read(unitcommand,*) iedate,ietime
     191    if (old) call skplin(3,unitcommand)
     192    read(unitcommand,*) loutstep
     193    if (old) call skplin(3,unitcommand)
     194    read(unitcommand,*) loutaver
     195    if (old) call skplin(3,unitcommand)
     196    read(unitcommand,*) loutsample
     197    if (old) call skplin(3,unitcommand)
     198    read(unitcommand,*) itsplit
     199    if (old) call skplin(3,unitcommand)
     200    read(unitcommand,*) lsynctime
     201    if (old) call skplin(3,unitcommand)
     202    read(unitcommand,*) ctl
     203    if (old) call skplin(3,unitcommand)
     204    read(unitcommand,*) ifine
     205    if (old) call skplin(3,unitcommand)
     206    read(unitcommand,*) iout
     207    if (old) call skplin(3,unitcommand)
     208    read(unitcommand,*) ipout
     209    if (old) call skplin(3,unitcommand)
     210    read(unitcommand,*) lsubgrid
     211    if (old) call skplin(3,unitcommand)
     212    read(unitcommand,*) lconvection
     213    if (old) call skplin(3,unitcommand)
     214    read(unitcommand,*) lagespectra
     215    if (old) call skplin(3,unitcommand)
     216    read(unitcommand,*) ipin
     217    if (old) call skplin(3,unitcommand)
     218    read(unitcommand,*) ioutputforeachrelease
     219    if (old) call skplin(3,unitcommand)
     220    read(unitcommand,*) iflux
     221    if (old) call skplin(3,unitcommand)
     222    read(unitcommand,*) mdomainfill
     223    if (old) call skplin(3,unitcommand)
     224    read(unitcommand,*) ind_source
     225    if (old) call skplin(3,unitcommand)
     226    read(unitcommand,*) ind_receptor
     227    if (old) call skplin(3,unitcommand)
     228    read(unitcommand,*) mquasilag
     229    if (old) call skplin(3,unitcommand)
     230    read(unitcommand,*) nested_output
     231    if (old) call skplin(3,unitcommand)
     232    read(unitcommand,*) linit_cond
     233    close(unitcommand)
     234
     235  endif ! input format
     236
     237  ! write command file in namelist format to output directory if requested
     238  if (nmlout.eqv..true.) then
     239    !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
     240    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
     241    write(unitcommand,nml=command)
     242    close(unitcommand)
     243     ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
     244     ! write(unitheader,NML=COMMAND)
     245     !close(unitheader)
     246  endif
     247
    161248  ifine=max(ifine,1)
    162 
    163249
    164250  ! Determine how Markov chain is formulated (for w or for w/sigw)
     
    371457  endif
    372458
    373   if(lsubgrid.ne.1) then
     459  if(lsubgrid.ne.1.and.verbosity.eq.0) then
    374460    write(*,*) '             ----------------               '
    375461    write(*,*) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
     
    505591  stop
    506592
     5931000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
     594       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     595        write(*,'(a)') path(2)(1:length(1))
     596        stop
    507597end subroutine readcommand
  • trunk/src/readoutgrid.f90

    r4 r20  
    9191  if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
    9292       .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
     93    write(*,*) 'outlon0,outlat0:'
    9394    write(*,*) outlon0,outlat0
     95    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
    9496    write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
    9597    write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
  • trunk/src/readpartpositions.f90

    r4 r20  
    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
     80999 continue
    8681  do i=1,numpointin
    8782    read(unitpartin)
     
    152147  stop
    153148
     149!995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
     150995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
     151  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
     152  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
     153!  stop
     154  goto 999
    154155
    155156996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
  • trunk/src/readpaths.f90

    r4 r20  
    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) :: string_test
     54  character(1) :: character_test
    5055
    5156  ! Read the pathname information stored in unitpath
    5257  !*************************************************
    5358
    54 
    55   open(unitpath,file='pathnames',status='old',err=999)
     59  open(unitpath,file=trim(pathfile),status='old',err=999)
    5660
    5761  do i=1,numpath
    5862    read(unitpath,'(a)',err=998) path(i)
    5963    length(i)=index(path(i),' ')-1
     64
     65   
     66    string_test = path(i)
     67    character_test = string_test(length(i):length(i))
     68    !print*, 'character_test,  string_test ', character_test,  string_test
     69      if ((character_test .NE. '/') .AND. (i .LT. 4))  then
     70         print*, 'WARNING: path not ending in /'
     71         print*, path(i)
     72         path(i) = string_test(1:length(i)) // '/'
     73         length(i)=length(i)+1
     74         print*, 'fix: padded with /'
     75         print*, path(i)
     76         print*, 'length(i) increased 1'
     77      endif
    6078  end do
    6179
     
    7088    length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
    7189  end do
     90  print*,length(5),length(6)
    7291
    7392
  • trunk/src/readreleases.f90

    r4 r20  
    5757  ! weta, wetb          parameters to determine the wet scavenging coefficient *
    5858  ! zpoint1,zpoint2     height range, over which release takes place           *
     59  ! num_min_discrete    if less, release cannot be randomized and happens at   *
     60  !                     time mid-point of release interval                     *
    5961  !                                                                            *
    6062  !*****************************************************************************
     
    6870
    6971  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
     72  integer,parameter :: num_min_discrete=100
    7073  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    71   real(kind=dp) :: jul1,jul2,juldate
     74  real(kind=dp) :: jul1,jul2,julm,juldate
    7275  character(len=50) :: line
    7376  logical :: old
     
    376379  jul1=juldate(id1,it1)
    377380  jul2=juldate(id2,it2)
     381  julm=(jul1+jul2)/2.
    378382  if (jul1.gt.jul2) then
    379383    write(*,*) 'FLEXPART MODEL ERROR'
     
    391395        stop
    392396      endif
    393       ireleasestart(numpoint)=int((jul1-bdate)*86400.)
    394       ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     397      if (npart(numpoint).gt.num_min_discrete) then
     398        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     399        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     400      else
     401        ireleasestart(numpoint)=int((julm-bdate)*86400.)
     402        ireleaseend(numpoint)=int((julm-bdate)*86400.)
     403      endif
    395404    else if (ldirect.eq.-1) then
    396405      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
     
    401410        stop
    402411      endif
    403       ireleasestart(numpoint)=int((jul1-bdate)*86400.)
    404       ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     412      if (npart(numpoint).gt.num_min_discrete) then
     413        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
     414        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
     415      else
     416        ireleasestart(numpoint)=int((julm-bdate)*86400.)
     417        ireleaseend(numpoint)=int((julm-bdate)*86400.)
     418      endif
    405419    endif
    406420  endif
  • trunk/src/readwind.f90

    r4 r20  
    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
  • trunk/src/readwind_gfs.f90

    r4 r20  
    703703  endif
    704704
    705 
    706705  if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
    707706  if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
  • trunk/src/timemanager.f90

    r4 r20  
    128128  !**********************************************************************
    129129
    130   !write (*,*) 'starting simulation'
     130
     131  if (verbosity.gt.0) then
     132    write (*,*) 'timemanager> starting simulation'
     133    if (verbosity.gt.1) then
     134      CALL SYSTEM_CLOCK(count_clock)
     135      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
     136    endif     
     137  endif
     138
    131139  do itime=0,ideltas,lsynctime
    132140
     
    142150  !********************************************************************
    143151
    144     if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) &
     152    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
     153        if (verbosity.gt.0) then
     154           write (*,*) 'timemanager> call wetdepo'
     155        endif     
    145156         call wetdepo(itime,lsynctime,loutnext)
     157    endif
    146158
    147159    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
     
    156168  !*************************************
    157169
    158       if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) &
    159            call convmix(itime)
     170   if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
     171        if (verbosity.gt.0) then
     172           write (*,*) 'timemanager> call convmix -- backward'
     173        endif         
     174      call convmix(itime)
     175        if (verbosity.gt.1) then
     176          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     177          CALL SYSTEM_CLOCK(count_clock)
     178          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
     179        endif
     180   endif
    160181
    161182  ! Get necessary wind fields if not available
    162183  !*******************************************
    163 
     184    if (verbosity.gt.0) then
     185           write (*,*) 'timemanager> call getfields'
     186    endif
    164187    call getfields(itime,nstop1)
     188        if (verbosity.gt.1) then
     189          CALL SYSTEM_CLOCK(count_clock)
     190          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
     191        endif
    165192    if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
     193
    166194  ! Release particles
    167195  !******************
    168196
     197    if (verbosity.gt.0) then
     198           write (*,*) 'timemanager>  Release particles'
     199    endif
     200
    169201    if (mdomainfill.ge.1) then
    170202      if (itime.eq.0) then
     203        if (verbosity.gt.0) then
     204          write (*,*) 'timemanager>  call init_domainfill'
     205        endif       
    171206        call init_domainfill
    172207      else
     208        if (verbosity.gt.0) then
     209          write (*,*) 'timemanager>  call boundcond_domainfill'
     210        endif   
    173211        call boundcond_domainfill(itime,loutend)
    174212      endif
    175213    else
     214      if (verbosity.gt.0) then
     215        print*,'call releaseparticles' 
     216      endif
    176217      call releaseparticles(itime)
     218      if (verbosity.gt.1) then
     219        CALL SYSTEM_CLOCK(count_clock)
     220        WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
     221      endif
    177222    endif
    178223
     
    182227  !**************************************************************
    183228
    184       if ((ldirect.eq.1).and.(lconvection.eq.1)) &
    185            call convmix(itime)
    186 
     229   if ((ldirect.eq.1).and.(lconvection.eq.1)) then
     230     if (verbosity.gt.0) then
     231       write (*,*) 'timemanager> call convmix -- forward'
     232     endif   
     233     call convmix(itime)
     234   endif
    187235
    188236  ! If middle of averaging period of output fields is reached, accumulated
     
    299347      if ((itime.eq.loutend).and.(outnum.gt.0.)) then
    300348        if ((iout.le.3.).or.(iout.eq.5)) then
     349          if (surf_only.ne.1) then
    301350          call concoutput(itime,outnum,gridtotalunc, &
    302351               wetgridtotalunc,drygridtotalunc)
    303           if (nested_output.eq.1) call concoutput_nest(itime,outnum)
     352          else 
     353  if (verbosity.eq.1) then
     354     print*,'call concoutput_surf '
     355     CALL SYSTEM_CLOCK(count_clock)
     356     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
     357  endif
     358          call concoutput_surf(itime,outnum,gridtotalunc, &
     359               wetgridtotalunc,drygridtotalunc)
     360  if (verbosity.eq.1) then
     361     print*,'called concoutput_surf '
     362     CALL SYSTEM_CLOCK(count_clock)
     363     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
     364  endif
     365          endif
     366
     367          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
     368          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
    304369          outnum=0.
    305370        endif
  • trunk/src/verttransform.f90

    r4 r20  
    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
  • trunk/src/wetdepo.f90

    r4 r20  
    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
  • trunk/src/writeheader.f90

    r4 r20  
    6767
    6868  if (ldirect.eq.1) then
    69     write(unitheader) ibdate,ibtime,'FLEXPART V9.0'
     69    write(unitheader) ibdate,ibtime, trim(flexversion)
    7070  else
    71     write(unitheader) iedate,ietime,'FLEXPART V9.0'
     71    write(unitheader) iedate,ietime, trim(flexversion)
    7272  endif
    7373
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG