Changeset 20

Timestamp:

Dec 23, 2013, 6:23:38 PM (11 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:

• FLEXPART.f90 (modified) (15 diffs)
• com_mod.f90 (modified) (10 diffs)
• concoutput.f90 (modified) (2 diffs)
• concoutput_surf.f90 (added)
• concoutput_surf_nest.f90 (added)
• gridcheck.f90 (modified) (2 diffs)
• gridcheck_fnl.f90 (added)
• interpol_rain.f90 (modified) (6 diffs)
• makefile (added)
• makefile.ecmwf_absoft (deleted)
• makefile.ecmwf_gfortran (deleted)
• makefile.ecmwf_gfortran_jfb (deleted)
• makefile.gfs_absoft (deleted)
• makefile.gfs_gfortran (deleted)
• makefile.henrik (deleted)
• makefile.test (deleted)
• openouttraj.f90 (modified) (1 diff)
• outg_mod.f90 (modified) (1 diff)
• outgrid_init.f90 (modified) (1 diff)
• par_mod.f90 (modified) (5 diffs)
• readavailable.f90 (modified) (3 diffs)
• readcommand.f90 (modified) (3 diffs)
• readoutgrid.f90 (modified) (1 diff)
• readpartpositions.f90 (modified) (2 diffs)
• readpaths.f90 (modified) (4 diffs)
• readreleases.f90 (modified) (5 diffs)
• readwind.f90 (modified) (2 diffs)
• readwind_gfs.f90 (modified) (1 diff)
• timemanager.f90 (modified) (5 diffs)
• verttransform.f90 (modified) (5 diffs)
• wetdepo.f90 (modified) (9 diffs)
• writeheader.f90 (modified) (1 diff)
• writeheader_nest_surf.f90 (added)
• writeheader_surf.f90 (added)
• writeheader_txt.f90 (added)

trunk/src/FLEXPART.f90

@@ -49 +49 @@
   integer :: i,j,ix,jy,inest
   integer :: idummy = -320
+  character(len=256) :: inline_options  !pathfile, flexversion, arg2
+
 
   ! Generate a large number of random numbers
@@ -58 +60 @@
   call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
 
+  ! 
+  flexversion='Version 9.1.8  (2013-12-08)'
+  !verbosity=0
+  ! Read the pathnames where input/output files are stored
+  !*******************************************************
+
+  inline_options='none'
+  select case (iargc())
+  case (2)
+    call getarg(1,arg1)
+    pathfile=arg1
+    call getarg(2,arg2)
+    inline_options=arg2
+  case (1)
+    call getarg(1,arg1)
+    pathfile=arg1
+    verbosity=0
+    if (arg1(1:1).eq.'-') then
+        write(pathfile,'(a11)') './pathnames'
+        inline_options=arg1 
+    endif
+  case (0)
+    write(pathfile,'(a11)') './pathnames'
+    verbosity=0
+  end select
+  
+    if (inline_options(1:1).eq.'-') then
+      print*, 'inline options=', inline_options       
+      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
+         print*, 'verbose mode 1: additional information will be displayed'
+         verbosity=1
+      endif
+      if (trim(inline_options).eq.'-v2') then
+         print*, 'verbose mode 2: additional information will be displayed'
+         verbosity=2
+      endif
+      if (trim(inline_options).eq.'-i') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=1
+         info_flag=1
+      endif
+      if (trim(inline_options).eq.'-i2') then
+         print*, 'info mode: will provide run specific information and stop'
+         verbosity=2
+         info_flag=1
+      endif
+    endif
+
+
   ! Print the GPL License statement
   !*******************************************************
-  print*,'Welcome to FLEXPART Version 9.0'
+  print*,'Welcome to FLEXPART', trim(flexversion)
   print*,'FLEXPART is free software released under the GNU Genera'// &
        'l Public License.'
-
-  ! Read the pathnames where input/output files are stored
-  !*******************************************************
-
-  call readpaths
+            
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readpaths'
+  endif 
+  call readpaths(pathfile)
+   
+ 
+  if (verbosity.gt.1) then !show clock info 
+     !print*,'length(4)',length(4)
+     !count=0,count_rate=1000
+     CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
+     !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
+     !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
+  endif
+
 
   ! Read the user specifications for the current model run
   !*******************************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readcommand'
+  endif
   call readcommand
+  if (verbosity.gt.0) then
+        WRITE(*,*) '    ldirect=', ldirect
+        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
+        WRITE(*,*) '    iedate,ietime=', iedate,ietime
+        if (verbosity.gt.1) then   
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+        endif     
+  endif
 
 
   ! Read the age classes to be used
   !********************************
-
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readageclasses'
+  endif
   call readageclasses
+
+  if (verbosity.gt.1) then   
+                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+  endif     
+ 
 
 
@@ -84 +167 @@
   !******************************************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readavailable'
+  endif  
   call readavailable
-
 
   ! Read the model grid specifications,
   ! both for the mother domain and eventual nests
   !**********************************************
+ 
+  if (verbosity.gt.0) then
+     WRITE(*,*) 'call gridcheck'
+  endif
 
   call gridcheck
+
+  if (verbosity.gt.1) then   
+     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+  endif      
+  
+
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call gridcheck_nests'
+  endif  
   call gridcheck_nests
 
@@ -98 +197 @@
   !************************************
 
+  if (verbosity.gt.0) then
+        WRITE(*,*) 'call readoutgrid'
+  endif
+
   call readoutgrid
-  if (nested_output.eq.1) call readoutgrid_nest
-
+
+  if (nested_output.eq.1) then
+          call readoutgrid_nest
+    if (verbosity.gt.0) then
+        WRITE(*,*) '# readoutgrid_nest'
+    endif
+  endif
 
   ! Read the receptor points for which extra concentrations are to be calculated
   !*****************************************************************************
 
+  if (verbosity.eq.1) then
+     print*,'call readreceptors'
+  endif
   call readreceptors
-
 
   ! Read the physico-chemical species property table
@@ -117 +227 @@
   !***************************
 
+  if (verbosity.gt.0) then
+    print*,'call readlanduse'
+  endif
   call readlanduse
 
@@ -123 +236 @@
   !********************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'call assignland'
+  endif
   call assignland
 
@@ -130 +246 @@
   !**********************************************
 
+  if (verbosity.gt.0) then
+    print*,'call readreleases'
+  endif
   call readreleases
+
 
   ! Read and compute surface resistances to dry deposition of gases
   !****************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'call readdepo'
+  endif
   call readdepo
-
 
   ! Convert the release point coordinates from geografical to grid coordinates
   !***************************************************************************
 
-  call coordtrafo
+  call coordtrafo  
+  if (verbosity.gt.0) then
+    print*,'call coordtrafo'
+  endif
 
 
@@ -147 +272 @@
   !*****************************************
 
+  if (verbosity.gt.0) then
+    print*,'Initialize all particles to non-existent'
+  endif
   do j=1,maxpart
     itra1(j)=-999999999
@@ -155 +283 @@
 
   if (ipin.eq.1) then
+    if (verbosity.gt.0) then
+      print*,'call readpartpositions'
+    endif
     call readpartpositions
   else
+    if (verbosity.gt.0) then
+      print*,'numpart=0, numparticlecount=0'
+    endif    
     numpart=0
     numparticlecount=0
@@ -166 +300 @@
   !***************************************************************
 
+
+  if (verbosity.gt.0) then
+    print*,'call outgrid_init'
+  endif
   call outgrid_init
   if (nested_output.eq.1) call outgrid_init_nest
@@ -173 +311 @@
   !******************
 
-  if (OHREA.eqv..TRUE.) &
+  if (OHREA.eqv..TRUE.) then
+    if (verbosity.gt.0) then
+       print*,'call readOHfield'
+    endif
        call readOHfield
+  endif
 
   ! Write basic information on the simulation to a file "header"
@@ -180 +322 @@
   !******************************************************************
 
+
+  if (verbosity.gt.0) then
+    print*,'call writeheader'
+  endif
+
   call writeheader
-  if (nested_output.eq.1) call writeheader_nest
-  open(unitdates,file=path(2)(1:length(2))//'dates')
+  ! FLEXPART 9.2 ticket ?? write header in ASCII format 
+  call writeheader_txt
+  !if (nested_output.eq.1) call writeheader_nest
+  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
+
+  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
+  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
+
+
+
+  !open(unitdates,file=path(2)(1:length(2))//'dates')
+
+  if (verbosity.gt.0) then
+    print*,'call openreceptors'
+  endif
   call openreceptors
   if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
@@ -190 +350 @@
   !***************************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'discretize release times'
+  endif
   do i=1,numpoint
     ireleasestart(i)=nint(real(ireleasestart(i))/ &
@@ -200 +363 @@
   ! Initialize cloud-base mass fluxes for the convection scheme
   !************************************************************
+
+  if (verbosity.gt.0) then
+    print*,'Initialize cloud-base mass fluxes for the convection scheme'
+  endif
 
   do jy=0,nymin1
@@ -218 +385 @@
   !********************************
 
+  if (verbosity.gt.0) then
+     if (verbosity.gt.1) then   
+       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+     endif
+     if (info_flag.eq.1) then
+         print*, 'info only mode (stop)'    
+         stop
+     endif
+     print*,'call timemanager'
+  endif
+
   call timemanager
 

trunk/src/com_mod.f90

@@ -7 +7 @@
 !        June 1996                                                             *
 !                                                                              *
-!        Last update: 9 August 2000                                            *
+!        Last update:15 August 2013 IP                                         *
 !                                                                              *
 !*******************************************************************************
@@ -25 +25 @@
   character :: path(numpath+2*maxnests)*120
   integer :: length(numpath+2*maxnests)
-
+  character(len=256) :: pathfile, flexversion, arg1, arg2
+  
   ! path                    path names needed for trajectory model
   ! length                  length of path names needed for trajectory model
-
+  ! pathfile                file where pathnames are stored
+  ! flexversion             version of flexpart
+  ! arg                     input arguments from launch at command line
 
   !********************************************************
@@ -65 +68 @@
   integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
   integer :: mquasilag,nested_output,ind_source,ind_receptor
-  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond
+  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
   logical :: turbswitch
 
@@ -91 +94 @@
   ! mquasilag 0: normal run
   !      1: Particle position output is produced in a condensed format and particles are numbered
+  ! surf_only   switch output in grid_time files for surface only or full vertical resolution
+  !      0=no (full vertical resolution), 1=yes (surface only)
   ! nested_output: 0 no, 1 yes
   ! turbswitch              determines how the Markov chain is formulated
@@ -139 +144 @@
   real :: decay(maxspec)
   real :: weta(maxspec),wetb(maxspec)
+! NIK: 31.01.2013- parameters for in-cloud scavening
+  real :: weta_in(maxspec), wetb_in(maxspec), wetc_in(maxspec), wetd_in(maxspec)
   real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
   real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
@@ -172 +179 @@
 
   ! WET DEPOSITION
-  ! weta, wetb              parameters for determining wet scavenging coefficients
+  ! weta, wetb              parameters for determining below-cloud wet scavenging coefficients
+  ! weta_in, wetb_in       parameters for determining in-cloud wet scavenging coefficients
+  ! wetc_in, wetd_in       parameters for determining in-cloud wet scavenging coefficients
 
   ! GAS DEPOSITION
@@ -331 +340 @@
   real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
   real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
+  !scavenging NIK, PS
   integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
   integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
+      integer icloudbot(0:nxmax-1,0:nymax-1,2)
+      integer icloudthck(0:nxmax-1,0:nymax-1,2)
+
 
   ! uu,vv,ww [m/2]       wind components in x,y and z direction
@@ -346 +359 @@
   !      rainout  conv/lsp dominated  2/3
   !      washout  conv/lsp dominated  4/5
+! PS 2013
+!c icloudbot (m)        cloud bottom height
+!c icloudthck (m)       cloud thickness     
+
   ! pplev for the GFS version
 
@@ -662 +679 @@
 
 
-
   !********************
   ! Random number field
@@ -671 +687 @@
   ! rannumb                 field of normally distributed random numbers
 
+  !********************
+  ! Verbosity, testing flags
+  !********************   
+  integer :: verbosity=0
+  integer :: info_flag=0
+  INTEGER :: count_clock, count_clock0,  count_rate, count_max
+   
 
 end module com_mod

trunk/src/concoutput.f90

@@ -108 +108 @@
   write(adate,'(i8.8)') jjjjmmdd
   write(atime,'(i6.6)') ihmmss
+  open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
   write(unitdates,'(a)') adate//atime
-
+  close(unitdates)  
 
   ! For forward simulations, output fields have dimension MAXSPEC,
@@ -158 +159 @@
         yl=outlat0+real(jy)*dyout
         xl=(xl-xlon0)/dx
-        yl=(yl-ylat0)/dx
+        yl=(yl-ylat0)/dy !v9.1.1 
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)

trunk/src/gridcheck.f90

@@ -98 +98 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  real(kind=4) :: zsec2(184),zsec4(jpunp)
   character(len=1) :: opt
 
@@ -466 +467 @@
     k=nlev_ec+1+numskip+i
     akm(nwz-i+1)=zsec2(j)
-  !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
     bkm(nwz-i+1)=zsec2(k)
   end do

trunk/src/interpol_rain.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
-       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
-  !                          i   i   i    i    i     i   i
-  !i    i    i  i    i     i      i      i     o     o     o
+!subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
+!       ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+!  !                          i   i   i    i    i     i   i
+!  !i    i    i  i    i     i      i      i     o     o     o
+
+      subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
+     ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+     intiy1,intiy2,icmv)
+!                               i   i   i    i    i     i   i 
+!     i    i    i  i    i     i      i      i     o     o     o
+
   !****************************************************************************
   !                                                                           *
@@ -40 +47 @@
   !                                                                           *
   !     30 August 1996                                                        *
-  !                                                                           *
+  !
+  !*  Petra Seibert, 2011/2012:
+  !*  Add interpolation of cloud bottom and cloud thickness
+  !*  for fix to SE's new wet scavenging scheme  *
   !****************************************************************************
   !                                                                           *
@@ -70 +80 @@
 
   integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
-  integer :: itime,itime1,itime2,level,indexh
+  !integer :: itime,itime1,itime2,level,indexh
+  integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4 
+  integer :: intiy1,intiy2,ipsum,icmv  
   real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
   real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
-  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
-  real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
+  integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
+  real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
+  !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+  real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4  
+  !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
 
 
@@ -127 +142 @@
          + p3*yy3(ix ,jyp,level,indexh) &
          + p4*yy3(ixp,jyp,level,indexh)
+
+!CPS clouds:
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi1(m)=icmv
+        else
+          yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh)    &
+             + ip2*p2*iy1(ixp,jy ,indexh)       &
+             + ip3*p3*iy1(ix ,jyp,indexh)       &
+             + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
+        endif
+        
+        ip1=1
+        ip2=1
+        ip3=1
+        ip4=1
+        if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
+        if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
+        if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
+        if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
+        ipsum= ip1+ip2+ip3+ip4
+        if (ipsum .eq. 0) then
+          yi2(m)=icmv
+        else
+          yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh)  &
+             + ip2*p2*iy2(ixp,jy ,indexh)     &
+             + ip3*p3*iy2(ix ,jyp,indexh)     &
+             + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
+        endif
+!CPS end clouds
+
   end do
 
@@ -134 +188 @@
   !************************************
 
+      if (abs(itime) .lt. abs(itime1)) then
+        print*,'interpol_rain.f90'
+        print*,itime,itime1,itime2
+        stop 'ITIME PROBLEM'
+      endif
+
+
   dt1=real(itime-itime1)
   dt2=real(itime2-itime)
@@ -143 +204 @@
 
 
+!PS clouds:
+      intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
+      if (yi1(1) .eq. float(icmv)) intiy1=yi1(2) 
+      if (yi1(2) .eq. float(icmv)) intiy1=yi1(1) 
+
+      intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
+      if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)  
+      if (yi2(2) .eq. float(icmv)) intiy2=yi2(1) 
+      
+      if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
+        intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
+      else
+        intiy1=icmv
+        intiy2=icmv
+      endif
+!PS end clouds
+
 end subroutine interpol_rain

trunk/src/openouttraj.f90

@@ -54 +54 @@
 
   if (ldirect.eq.1) then
-  write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime,'FLEXPART V8.2'
+  write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime, trim(flexversion)
   else
-  write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime,'FLEXPART V8.2'
+  write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime, trim(flexversion)
   endif
   write(unitouttraj,*) method,lsubgrid,lconvection

trunk/src/outg_mod.f90

@@ -43 +43 @@
   real,allocatable, dimension (:,:) :: wetgridsigma
   real,allocatable, dimension (:) :: sparse_dump_r
+  real,allocatable, dimension (:) :: sparse_dump_u
   integer,allocatable, dimension (:) :: sparse_dump_i
 

trunk/src/outgrid_init.f90

@@ -247 +247 @@
        max(numygrid,numygridn)*numzgrid),stat=stat)
     if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
+   allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
+       max(numygrid,numygridn)*numzgrid),stat=stat)
+        if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
   allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
        max(numygrid,numygridn)*numzgrid),stat=stat)

trunk/src/par_mod.f90

@@ -28 +28 @@
 !        1997                                                                  *
 !                                                                              *
-!        Last update 10 August 2000                                            *
+!        Last update 15 August 2013 IP                                         *
 !                                                                              *
 !*******************************************************************************
@@ -120 +120 @@
   ! Maximum dimensions of the input mother grids
   !*********************************************
-
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
+  
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
   !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
-  integer,parameter :: nxshift=359  ! for ECMWF
-  !integer,parameter :: nxshift=0     ! for GFS
+  !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
+
+  integer,parameter :: nxshift=359 ! for ECMWF
+  !integer,parameter :: nxshift=0     ! for GFS or FNL (XF)
 
   integer,parameter :: nconvlevmax = nuvzmax-1
@@ -150 +153 @@
   !*********************************************
 
-  integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
-  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
-
+  !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
+  !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
+  integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
   ! maxnests                maximum number of nested grids
   ! nxmaxn,nymaxn           maximum dimension of nested wind fields in
@@ -194 +197 @@
   !**************************************************
 
-  integer,parameter :: maxpart=2000000
-  integer,parameter :: maxspec=4
+  !integer,parameter :: maxpart=4000000
+  integer,parameter :: maxpart=2000
+  integer,parameter :: maxspec=6
 
 
@@ -255 +259 @@
   integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
   integer,parameter :: unitOH=1
-  integer,parameter :: unitdates=94, unitheader=90, unitshortpart=95
+  integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
   integer,parameter :: unitboundcond=89
 
+!******************************************************
+! integer code for missing values, used in wet scavenging (PS, 2012)
+!******************************************************
+
+    !  integer icmv
+    !  parameter(icmv=-9999)
+      integer,parameter ::  icmv=-9999
+
+! Parameters for testing
+!*******************************************
+!  integer :: verbosity=0
+
 end module par_mod

trunk/src/readavailable.f90

@@ -123 +123 @@
 
   do k=1,numbnests
+  print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
+  print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     open(unitavailab,file=path(numpath+2*(k-1)+2) &
          (1:length(numpath+2*(k-1)+2)),status='old',err=998)
@@ -275 +277 @@
   return
 
-998   write(*,*) ' #### FLEXPART MODEL ERROR! FILE   #### '
+998   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE   #### '
   write(*,'(a)') '     '//path(numpath+2*(k-1)+2) &
        (1:length(numpath+2*(k-1)+2))
@@ -281 +283 @@
   stop
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
+999   write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
   write(*,'(a)') '     '//path(4)(1:length(4))
   write(*,*) ' #### CANNOT BE OPENED           #### '

trunk/src/readcommand.f90

@@ -80 +80 @@
   character(len=50) :: line
   logical :: old
-
+  logical :: nmlout=.true. !.false.
+  integer :: readerror
+
+  namelist /command/ &
+    ldirect, &
+    ibdate,ibtime, &
+    iedate,ietime, &
+    loutstep, &
+    loutaver, &
+    loutsample, &
+    itsplit, &
+    lsynctime, &
+    ctl, &
+    ifine, &
+    iout, &
+    ipout, &
+    lsubgrid, &
+    lconvection, &
+    lagespectra, &
+    ipin, &
+    ioutputforeachrelease, &
+    iflux, &
+    mdomainfill, &
+    ind_source, &
+    ind_receptor, &
+    mquasilag, &
+    nested_output, &
+    linit_cond, &
+    surf_only    
+
+  ! Presetting namelist command
+  ldirect=1
+  ibdate=20000101
+  ibtime=0
+  iedate=20000102
+  ietime=0
+  loutstep=10800
+  loutaver=10800
+  loutsample=900
+  itsplit=999999999
+  lsynctime=900
+  ctl=-5.0
+  ifine=4
+  iout=3
+  ipout=0
+  lsubgrid=1
+  lconvection=1
+  lagespectra=0
+  ipin=1
+  ioutputforeachrelease=1
+  iflux=1
+  mdomainfill=0
+  ind_source=1
+  ind_receptor=1
+  mquasilag=0
+  nested_output=0
+  linit_cond=0
+  surf_only=0 
 
   ! Open the command file and read user options
-  !********************************************
-
-
+  ! Namelist input first: try to read as namelist file
+  !**************************************************************************
   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-       err=999)
-
-  ! Check the format of the COMMAND file (either in free format,
-  ! or using formatted mask)
-  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
-  !**************************************************************************
-
-  call skplin(9,unitcommand)
-  read (unitcommand,901) line
-901   format (a)
-  if (index(line,'LDIRECT') .eq. 0) then
-    old = .false.
-  else
-    old = .true.
-  endif
-  rewind(unitcommand)
-
-  ! Read parameters
-  !****************
-
-  call skplin(7,unitcommand)
-  if (old) call skplin(1,unitcommand)
-
-  read(unitcommand,*) ldirect
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ibdate,ibtime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iedate,ietime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutstep
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutaver
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) loutsample
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) itsplit
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lsynctime
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ctl
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ifine
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iout
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ipout
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lsubgrid
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lconvection
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) lagespectra
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ipin
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ioutputforeachrelease
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) iflux
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) mdomainfill
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ind_source
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) ind_receptor
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) mquasilag
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) nested_output
-  if (old) call skplin(3,unitcommand)
-  read(unitcommand,*) linit_cond
+         form='formatted',iostat=readerror)
+  ! If fail, check if file does not exist
+  if (readerror.ne.0) then
+    
+    print*,'***ERROR: file COMMAND not found in ' 
+    print*, path(1)(1:length(1))//'COMMAND'
+    print*, 'Check your pathnames file.'
+    stop
+
+  endif
+
+  read(unitcommand,command,iostat=readerror)
   close(unitcommand)
 
+  ! If error in namelist format, try to open with old input code
+  if (readerror.ne.0) then
+
+    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
+         err=999)
+
+    ! Check the format of the COMMAND file (either in free format,
+    ! or using formatted mask)
+    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+    !**************************************************************************
+
+    call skplin(9,unitcommand)
+    read (unitcommand,901) line
+  901   format (a)
+    if (index(line,'LDIRECT') .eq. 0) then
+      old = .false.
+    else
+      old = .true.
+    endif
+    rewind(unitcommand)
+
+    ! Read parameters
+    !****************
+
+    call skplin(7,unitcommand)
+    if (old) call skplin(1,unitcommand)
+
+    read(unitcommand,*) ldirect
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ibdate,ibtime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iedate,ietime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutstep
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutaver
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) loutsample
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) itsplit
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lsynctime
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ctl
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ifine
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iout
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ipout
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lsubgrid
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lconvection
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) lagespectra
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ipin
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ioutputforeachrelease
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) iflux
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) mdomainfill
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ind_source
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) ind_receptor
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) mquasilag
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) nested_output
+    if (old) call skplin(3,unitcommand)
+    read(unitcommand,*) linit_cond
+    close(unitcommand)
+
+  endif ! input format
+
+  ! write command file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
+    open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
+    write(unitcommand,nml=command)
+    close(unitcommand)
+     ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
+     ! write(unitheader,NML=COMMAND)
+     !close(unitheader) 
+  endif
+
   ifine=max(ifine,1)
-
 
   ! Determine how Markov chain is formulated (for w or for w/sigw)
@@ -371 +457 @@
   endif
 
-  if(lsubgrid.ne.1) then
+  if(lsubgrid.ne.1.and.verbosity.eq.0) then
     write(*,*) '             ----------------               '
     write(*,*) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
@@ -505 +591 @@
   stop
 
+1000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
+       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+        write(*,'(a)') path(2)(1:length(1))
+        stop
 end subroutine readcommand

trunk/src/readoutgrid.f90

@@ -91 +91 @@
   if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
        .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
+    write(*,*) 'outlon0,outlat0:'
     write(*,*) outlon0,outlat0
+    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
     write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'

trunk/src/readpartpositions.f90

@@ -77 +77 @@
 
   read(unitpartin) numpointin
-  if (numpointin.ne.numpoint) then
-     write(*,*) ' ####              WARNING                     #### '
-     write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
-     write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
-     write(*,*) ' #### THIS IS JUST A CHECK TO SEE IF YOU READ  #### '
-     write(*,*) ' #### THE CORRECT PARTICLE DUMP FILE.          #### '
-  endif
+  if (numpointin.ne.numpoint) goto 995
+999 continue 
   do i=1,numpointin
     read(unitpartin)
@@ -152 +147 @@
   stop
 
+!995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
+995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
+  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
+  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
+!  stop
+  goto 999 
 
 996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '

trunk/src/readpaths.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine readpaths
+subroutine readpaths !(pathfile)
 
   !*****************************************************************************
@@ -30 +30 @@
   !                                                                            *
   !     1 February 1994                                                        *
+  !     last modified                                                          *
+  !     HS, 7.9.2012                                                           *
+  !     option to give pathnames file as command line option                   *
   !                                                                            *
   !*****************************************************************************
@@ -47 +50 @@
   implicit none
 
-  integer :: i
+  integer   :: i
+  character(256) :: string_test 
+  character(1) :: character_test 
 
   ! Read the pathname information stored in unitpath
   !*************************************************
 
-
-  open(unitpath,file='pathnames',status='old',err=999)
+  open(unitpath,file=trim(pathfile),status='old',err=999)
 
   do i=1,numpath
     read(unitpath,'(a)',err=998) path(i)
     length(i)=index(path(i),' ')-1
+
+    
+    string_test = path(i)
+    character_test = string_test(length(i):length(i))
+    !print*, 'character_test,  string_test ', character_test,  string_test 
+      if ((character_test .NE. '/') .AND. (i .LT. 4))  then
+         print*, 'WARNING: path not ending in /' 
+         print*, path(i)
+         path(i) = string_test(1:length(i)) // '/'
+         length(i)=length(i)+1
+         print*, 'fix: padded with /' 
+         print*, path(i)
+         print*, 'length(i) increased 1' 
+      endif
   end do
 
@@ -70 +88 @@
     length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
   end do
+  print*,length(5),length(6)
 
 

trunk/src/readreleases.f90

@@ -57 +57 @@
   ! weta, wetb          parameters to determine the wet scavenging coefficient *
   ! zpoint1,zpoint2     height range, over which release takes place           *
+  ! num_min_discrete    if less, release cannot be randomized and happens at   *
+  !                     time mid-point of release interval                     *
   !                                                                            *
   !*****************************************************************************
@@ -68 +70 @@
 
   integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
-  real(kind=dp) :: jul1,jul2,juldate
+  real(kind=dp) :: jul1,jul2,julm,juldate
   character(len=50) :: line
   logical :: old
@@ -376 +379 @@
   jul1=juldate(id1,it1)
   jul2=juldate(id2,it2)
+  julm=(jul1+jul2)/2.
   if (jul1.gt.jul2) then
     write(*,*) 'FLEXPART MODEL ERROR'
@@ -391 +395 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     else if (ldirect.eq.-1) then
       if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
@@ -401 +410 @@
         stop
       endif
-      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      if (npart(numpoint).gt.num_min_discrete) then
+        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+      else
+        ireleasestart(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+      endif
     endif
   endif

trunk/src/readwind.f90

@@ -101 +101 @@
   character(len=24) :: gribErrorMsg = 'Error reading grib file'
   character(len=20) :: gribFunction = 'readwind'
+
+  !HSO conversion of ECMWF etadot to etadot*dp/deta
+  logical :: etacon=.false.
+  real,parameter :: p00=101325.
+  real :: dak,dbk
 
   hflswitch=.false.
@@ -365 +370 @@
   endif
 
+  ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
+  if(etacon.eqv..true.) then
+    do k=1,nwzmax
+      dak=akm(k+1)-akm(k)
+      dbk=bkm(k+1)-bkm(k)
+      do i=0,nxmin1
+        do j=0,nymin1 
+          wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
+          if (k.gt.1) then 
+            wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
+          endif
+        end do
+      end do
+    end do
+  endif
+
   ! For global fields, assign the leftmost data column also to the rightmost
   ! data column; if required, shift whole grid by nxshift grid points

trunk/src/readwind_gfs.f90

@@ -703 +703 @@
   endif
 
-
   if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
   if(iumax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'

trunk/src/timemanager.f90

@@ -128 +128 @@
   !**********************************************************************
 
-  !write (*,*) 'starting simulation'
+
+  if (verbosity.gt.0) then
+    write (*,*) 'timemanager> starting simulation'
+    if (verbosity.gt.1) then
+      CALL SYSTEM_CLOCK(count_clock)
+      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+    endif     
+  endif
+
   do itime=0,ideltas,lsynctime
 
@@ -142 +150 @@
   !********************************************************************
 
-    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) &
+    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+        if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call wetdepo'
+        endif     
          call wetdepo(itime,lsynctime,loutnext)
+    endif
 
     if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
@@ -156 +168 @@
   !*************************************
 
-      if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) &
-           call convmix(itime)
+   if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
+        if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call convmix -- backward'
+        endif         
+      call convmix(itime)
+        if (verbosity.gt.1) then
+          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+          CALL SYSTEM_CLOCK(count_clock)
+          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+        endif 
+   endif
 
   ! Get necessary wind fields if not available
   !*******************************************
-
+    if (verbosity.gt.0) then
+           write (*,*) 'timemanager> call getfields'
+    endif 
     call getfields(itime,nstop1)
+        if (verbosity.gt.1) then
+          CALL SYSTEM_CLOCK(count_clock)
+          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+        endif 
     if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
+
   ! Release particles
   !******************
 
+    if (verbosity.gt.0) then
+           write (*,*) 'timemanager>  Release particles'
+    endif 
+
     if (mdomainfill.ge.1) then
       if (itime.eq.0) then
+        if (verbosity.gt.0) then
+          write (*,*) 'timemanager>  call init_domainfill'
+        endif       
         call init_domainfill
       else
+        if (verbosity.gt.0) then
+          write (*,*) 'timemanager>  call boundcond_domainfill'
+        endif   
         call boundcond_domainfill(itime,loutend)
       endif
     else
+      if (verbosity.gt.0) then
+        print*,'call releaseparticles'  
+      endif
       call releaseparticles(itime)
+      if (verbosity.gt.1) then
+        CALL SYSTEM_CLOCK(count_clock)
+        WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+      endif 
     endif
 
@@ -182 +227 @@
   !**************************************************************
 
-      if ((ldirect.eq.1).and.(lconvection.eq.1)) &
-           call convmix(itime)
-
+   if ((ldirect.eq.1).and.(lconvection.eq.1)) then
+     if (verbosity.gt.0) then
+       write (*,*) 'timemanager> call convmix -- forward'
+     endif    
+     call convmix(itime)
+   endif
 
   ! If middle of averaging period of output fields is reached, accumulated
@@ -299 +347 @@
       if ((itime.eq.loutend).and.(outnum.gt.0.)) then
         if ((iout.le.3.).or.(iout.eq.5)) then
+          if (surf_only.ne.1) then 
           call concoutput(itime,outnum,gridtotalunc, &
                wetgridtotalunc,drygridtotalunc)
-          if (nested_output.eq.1) call concoutput_nest(itime,outnum)
+          else  
+  if (verbosity.eq.1) then
+     print*,'call concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+          call concoutput_surf(itime,outnum,gridtotalunc, &
+               wetgridtotalunc,drygridtotalunc)
+  if (verbosity.eq.1) then
+     print*,'called concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+          endif
+
+          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
+          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
           outnum=0.
         endif

trunk/src/verttransform.f90

@@ -49 +49 @@
   ! Sabine Eckhardt, March 2007
   ! added the variable cloud for use with scavenging - descr. in com_mod
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
+  ! note that also other subroutines are affected by the fix
   !*****************************************************************************
   !                                                                            *
@@ -70 +72 @@
 
   integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
-  integer :: rain_cloud_above,kz_inv
+ integer :: rain_cloud_above,kz_inv !SE
+  integer icloudtop !PS
   real :: f_qvsat,pressure
-  real :: rh,lsp,convp
+ !real :: rh,lsp,convp
+  real :: rh,lsp,convp,prec,rhmin  
   real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
   real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
   real :: xlon,ylat,xlonr,dzdx,dzdy
-  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
   real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
   real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -83 +87 @@
   real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
   real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+  logical lconvectprec
   real,parameter :: const=r_air/ga
+  parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
 
   logical :: init = .true.
@@ -545 +551 @@
   !   create a cloud and rainout/washout field, clouds occur where rh>80%
   !   total cloudheight is stored at level 0
+
+
+
   do jy=0,nymin1
     do ix=0,nxmin1
-      rain_cloud_above=0
-      lsp=lsprec(ix,jy,1,n)
-      convp=convprec(ix,jy,1,n)
-      cloudsh(ix,jy,n)=0
-      do kz_inv=1,nz-1
-         kz=nz-kz_inv+1
-         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
-         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
-         clouds(ix,jy,kz,n)=0
-         if (rh.gt.0.8) then ! in cloud
-            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
-               rain_cloud_above=1
-               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
-                    height(kz)-height(kz-1)
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
-               else
-                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
-               endif
-            else ! no precipitation
-                  clouds(ix,jy,kz,n)=1 ! cloud
+
+
+
+  !    rain_cloud_above=0
+  !    lsp=lsprec(ix,jy,1,n)
+  !    convp=convprec(ix,jy,1,n)
+  !    cloudsh(ix,jy,n)=0
+  !    do kz_inv=1,nz-1
+  !       kz=nz-kz_inv+1
+  !       pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+  !       rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+  !       clouds(ix,jy,kz,n)=0
+  !       if (rh.gt.0.8) then ! in cloud
+  !          if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+  !             rain_cloud_above=1
+  !             cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+  !                  height(kz)-height(kz-1)
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+  !             else
+  !                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+  !             endif
+  !          else ! no precipitation
+  !                clouds(ix,jy,kz,n)=1 ! cloud
+  !          endif
+  !       else ! no cloud
+  !          if (rain_cloud_above.eq.1) then ! scavenging
+  !             if (lsp.ge.convp) then
+  !                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+  !             else
+  !                clouds(ix,jy,kz,n)=4 ! convp dominated washout
+  !             endif
+  !          endif
+  !       endif
+  !    end do
+
+
+   ! PS 3012
+
+             lsp=lsprec(ix,jy,1,n)
+          convp=convprec(ix,jy,1,n)
+          prec=lsp+convp
+          if (lsp.gt.convp) then !  prectype='lsp'
+            lconvectprec = .false.
+          else ! prectype='cp '
+            lconvectprec = .true.
+          endif
+          rhmin = 0.90 ! standard condition for presence of clouds
+!PS       note that original by Sabine Eckhart was 80%
+!PS       however, for T<-20 C we consider saturation over ice
+!PS       so I think 90% should be enough          
+          icloudbot(ix,jy,n)=icmv
+          icloudtop=icmv ! this is just a local variable
+98        do kz=1,nz
+            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
+            if (rh .gt. rhmin) then
+              if (icloudbot(ix,jy,n) .eq. icmv) then
+                icloudbot(ix,jy,n)=nint(height(kz))
+              endif
+              icloudtop=nint(height(kz)) ! use int to save memory
             endif
-         else ! no cloud
-            if (rain_cloud_above.eq.1) then ! scavenging
-               if (lsp.ge.convp) then
-                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
-               else
-                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
-               endif
+          enddo
+
+!PS try to get a cloud thicker than 50 m 
+!PS if there is at least .01 mm/h  - changed to 0.002 and put into
+!PS parameter precpmin        
+          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
+              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
+              prec .gt. precmin) then
+            rhmin = rhmin - 0.05
+            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
+          endif
+!PS implement a rough fix for badly represented convection
+!PS is based on looking at a limited set of comparison data
+          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
+             prec .gt. precmin) then 
+
+            if (convp .lt. 0.1) then
+              icloudbot(ix,jy,n) = 500
+              icloudtop =         8000
+            else
+              icloudbot(ix,jy,n) = 0
+              icloudtop =      10000
             endif
-         endif
-      end do
+          endif
+          if (icloudtop .ne. icmv) then
+            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+          else
+            icloudthck(ix,jy,n) = icmv
+          endif
+!PS  get rid of too thin clouds      
+          if (icloudthck(ix,jy,n) .lt. 50) then
+            icloudbot(ix,jy,n)=icmv
+            icloudthck(ix,jy,n)=icmv
+          endif
+
+
     end do
   end do
@@ -605 +681 @@
   !104   continue
   ! close(4)
+
+
 end subroutine verttransform

trunk/src/wetdepo.f90

@@ -39 +39 @@
   ! Code may not be correct for decay of deposition!                           *
   !                                                                            *
+  ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
+  ! scheme, not dated
+  ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
   !*****************************************************************************
   !                                                                            *
@@ -71 +74 @@
 
   integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
-  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
-  integer :: ks, kp
+  integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
+  integer :: kz !PS scheme
+  integer :: ks, kp, n1,n2, icbot,ictop, indcloud
+  integer :: scheme_number ! NIK==1, PS ==2
   real :: S_i, act_temp, cl, cle ! in cloud scavenging
   real :: clouds_h ! cloud height for the specific grid point
-  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+  real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -146 +151 @@
     interp_time=nint(itime-0.5*ltsample)
 
-    if (ngrid.eq.0) then
-      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
-           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
-           memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    else
-      call interpol_rain_nests(lsprecn,convprecn,tccn, &
-           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
-           memtime(1),memtime(2),interp_time,lsp,convp,cc)
-    endif
-
-    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+!    PS nest case still needs to be implemented!!    
+!    if (ngrid.eq.0) then
+!      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+!           1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+           call interpol_rain(lsprec,convprec,tcc,     &
+            icloudbot,icloudthck,nxmax,nymax,1,nx,ny,         &
+            memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
+            memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv)  
+!    else
+!      call interpol_rain_nests(lsprecn,convprecn,tccn, &
+!           nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+!           memtime(1),memtime(2),interp_time,lsp,convp,cc)
+!    endif
+
+
+ !   if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+ !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
+        prec = lsp+convp
+        precsub = 0.01
+        if (prec .lt. precsub) then
+          goto 20
+        else
+          f = (prec-precsub)/prec
+          lsp = f*lsp
+          convp = f*convp
+        endif
+
 
   ! get the level were the actual particle is in
       do il=2,nz
         if (height(il).gt.ztra1(jpart)) then
-          hz=il-1
+          !hz=il-1
+          kz=il-1
           goto 26
         endif
@@ -174 +197 @@
   ! scavenging is done
 
-  if (ngrid.eq.0) then
-     clouds_v=clouds(ix,jy,hz,n)
-     clouds_h=cloudsh(ix,jy,n)
-  else
-     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
-     clouds_h=cloudsnh(ix,jy,n,ngrid)
-  endif
-  !write(*,*) 'there is
-  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
-  if (clouds_v.le.1) goto 20
-  !write (*,*) 'there is scavenging'
+!old scheme
+!  if (ngrid.eq.0) then
+!     clouds_v=clouds(ix,jy,hz,n)
+!     clouds_h=cloudsh(ix,jy,n)
+!  else
+!     clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+!     clouds_h=cloudsnh(ix,jy,n,ngrid)
+!  endif
+!  !write(*,*) 'there is
+!  !    + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
+!  if (clouds_v.le.1) goto 20
+!  !write (*,*) 'there is scavenging'
+
+  ! PS: part of 2011/2012 fix 
+
+        if (ztra1(jpart) .le. float(ictop)) then
+          if (ztra1(jpart) .gt. float(icbot)) then
+            indcloud = 2 ! in-cloud
+          else
+            indcloud = 1 ! below-cloud
+          endif
+        elseif (ictop .eq. icmv) then
+          indcloud = 0 ! no cloud found, use old scheme
+        else
+          goto 20 ! above cloud
+        endif
+
 
   ! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228 +267 @@
   !**********************************************************
 
-    do ks=1,nspec                                  ! loop over species
+    do ks=1,nspec            ! loop over species
       wetdeposit(ks)=0.
+
+    
+    !conflicting changes to the same routine: 1=NIK 2 =PS
+    scheme_number=2   
+    if (scheme_number.eq.1) then !NIK 
+
       if (weta(ks).gt.0.) then
         if (clouds_v.ge.4) then
@@ -236 +281 @@
           wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
   !        write(*,*) 'bel. wetscav: ',wetscav
+
         else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
   !        IN CLOUD SCAVENGING
@@ -245 +291 @@
              act_temp=tt(ix,jy,hz,n)
           endif
-          cl=2E-7*prec**0.36
+
+! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening 
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
+          cl=weta_in(ks)*prec**wetb_in(ks)
           if (dquer(ks).gt.0) then ! is particle
-            S_i=0.9/cl
+            S_i=wetc_in(ks)/cl
            else ! is gas
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
             S_i=1/cle
            endif
-           wetscav=S_i*prec/3.6E6/clouds_h
+           wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
   !         write(*,*) 'in. wetscav:'
   !    +          ,wetscav,cle,cl,act_temp,prec,clouds_h
@@ -289 +341 @@
          wetdeposit(ks)=0.
       endif ! weta(k)
+
+     elseif (scheme_number.eq.2) then ! PS
+
+!PS          indcloud=0 ! Use this for FOR TESTING, 
+!PS          will skip the new in/below cloud method !!!             
+
+          if (weta(ks).gt.0.) then
+            if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
+!C               for aerosols and not highliy soluble substances weta=5E-6
+              wetscav=weta(ks)*prec**wetb(ks)                ! scavenging coeff.
+!c             write(*,*) 'bel. wetscav: ',wetscav
+            elseif (indcloud .eq. 2) then !  IN CLOUD SCAVENGING
+              if (ngrid.gt.0) then
+                 act_temp=ttn(ix,jy,kz,n,ngrid)
+              else
+                 act_temp=tt(ix,jy,kz,n)
+              endif
+
+! from NIK              
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+
+
+              cl=2E-7*prec**0.36
+              if (dquer(ks).gt.0) then ! is particle
+                S_i=0.9/cl
+              else ! is gas
+                cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+                S_i=1/cle
+              endif
+              wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+            else ! PS: no cloud diagnosed, old scheme, 
+!CPS          using with fixed a,b for simplicity, one may wish to change!!
+              wetscav = 1.e-4*prec**0.62
+            endif
+
+
+            wetdeposit(ks)=xmass1(jpart,ks)*    &
+            ! (1.-exp(-wetscav*abs(ltsample)))*fraction  ! wet deposition
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction  ! fraction = grfraction (PS) 
+            restmass = xmass1(jpart,ks)-wetdeposit(ks)
+            if (ioutputforeachrelease.eq.1) then
+              kp=npoint(jpart)
+            else
+              kp=1
+            endif
+            if (restmass .gt. smallnum) then
+              xmass1(jpart,ks)=restmass
+!cccccccccccccccc depostatistic
+!c            wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!cccccccccccccccc depostatistic
+            else
+              xmass1(jpart,ks)=0.
+            endif
+!C Correct deposited mass to the last time step when radioactive decay of 
+!C gridded deposited mass was calculated
+            if (decay(ks).gt.0.) then
+              wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
+            endif
+          else  ! weta(k)<0
+             wetdeposit(ks)=0.
+          endif 
+
+     endif !on scheme
+
+
+
     end do
 
@@ -295 +415 @@
   !*****************************************************************************
 
-    if (ldirect.eq.1) then
-    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-         real(ytra1(jpart)),nage,kp)
-    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
-         nage,kp)
-    endif
+!    if (ldirect.eq.1) then
+!    call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
+!         real(ytra1(jpart)),nage,kp)
+!    if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
+!         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
+!         nage,kp)
+!    endif
+
+   !PS
+        if (ldirect.eq.1) then
+          call wetdepokernel(nclass(jpart),wetdeposit,     &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+          if (nested_output.eq.1)                          &
+           call wetdepokernel_nest(nclass(jpart),wetdeposit, &
+             sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)  
+        endif
+
 
 20  continue

trunk/src/writeheader.f90

@@ -67 +67 @@
 
   if (ldirect.eq.1) then
-    write(unitheader) ibdate,ibtime,'FLEXPART V9.0'
+    write(unitheader) ibdate,ibtime, trim(flexversion)
   else
-    write(unitheader) iedate,ietime,'FLEXPART V9.0'
+    write(unitheader) iedate,ietime, trim(flexversion)
   endif