move version 9.1.8 form branches to trunk. Contributions from HSO, saeck,...

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

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/trunk@20 ef8cc7e1-21b7-489e-abab-c1baa636049d

src/FLEXPART.f90

@@ -48,6 +48,8 @@ program flexpart
 
   integer :: i,j,ix,jy,inest
   integer :: idummy = -320
+  character(len=256) :: inline_options  !pathfile, flexversion, arg2
+
 
   ! Generate a large number of random numbers
   !******************************************
@@ -57,56 +59,164 @@ program flexpart
   end do
   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.'
+            
+  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 pathnames where input/output files are stored
-  !*******************************************************
-
-  call readpaths
 
   ! 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     
+ 
+
 
   ! Read, which wind fields are available within the modelling period
   !******************************************************************
 
+  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
 
 
   ! Read the output grid specifications
   !************************************
 
+  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
   !*************************************************
   !SEC: now only needed SPECIES are read in readreleases.f
@@ -116,12 +226,18 @@ program flexpart
   ! Read the landuse inventory
   !***************************
 
+  if (verbosity.gt.0) then
+    print*,'call readlanduse'
+  endif
   call readlanduse
 
 
   ! Assign fractional cover of landuse classes to each ECMWF grid point
   !********************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'call assignland'
+  endif
   call assignland
 
 
@@ -129,23 +245,35 @@ program flexpart
   ! Read the coordinates of the release locations
   !**********************************************
 
+  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
 
 
   ! Initialize all particles to non-existent
   !*****************************************
 
+  if (verbosity.gt.0) then
+    print*,'Initialize all particles to non-existent'
+  endif
   do j=1,maxpart
     itra1(j)=-999999999
   end do
@@ -154,8 +282,14 @@ program flexpart
   !*************************************************************
 
   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
   endif
@@ -165,6 +299,10 @@ program flexpart
   ! Allocate fluxes and OHfield if necessary
   !***************************************************************
 
+
+  if (verbosity.gt.0) then
+    print*,'call outgrid_init'
+  endif
   call outgrid_init
   if (nested_output.eq.1) call outgrid_init_nest
 
@@ -172,16 +310,38 @@ program flexpart
   ! Read the OH field
   !******************
 
-  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"
   ! and open files that are to be kept open throughout the simulation
   !******************************************************************
 
+
+  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
 
@@ -189,6 +349,9 @@ program flexpart
   ! Releases can only start and end at discrete times (multiples of lsynctime)
   !***************************************************************************
 
+  if (verbosity.gt.0) then
+    print*,'discretize release times'
+  endif
   do i=1,numpoint
     ireleasestart(i)=nint(real(ireleasestart(i))/ &
          real(lsynctime))*lsynctime
@@ -200,6 +363,10 @@ program flexpart
   ! 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
     do ix=0,nxmin1
       cbaseflux(ix,jy)=0.
@@ -217,6 +384,18 @@ program flexpart
   ! Calculate particle trajectories
   !********************************
 
+  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
 
 

src/com_mod.f90

@@ -6,7 +6,7 @@
 !                                                                              *
 !        June 1996                                                             *
 !                                                                              *
-!        Last update: 9 August 2000                                            *
+!        Last update:15 August 2013 IP                                         *
 !                                                                              *
 !*******************************************************************************
 
@@ -24,10 +24,13 @@ module com_mod
 
   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
 
   !********************************************************
   ! Variables defining the general model run specifications
@@ -64,7 +67,7 @@ module com_mod
   real :: ctl,fine
   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
 
   ! ctl      factor, by which time step must be smaller than Lagrangian time scale
@@ -90,6 +93,8 @@ module com_mod
   !     0=no, 1=mass unit, 2=mass mixing ratio unit
   ! 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
 
@@ -138,6 +143,8 @@ module com_mod
   !real xmass(maxpoint,maxspec)
   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)
   real :: vsetaver(maxspec),cunningham(maxspec),weightmolar(maxspec)
@@ -171,7 +178,9 @@ module com_mod
   ! decay                   decay constant of radionuclide
 
   ! 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
   ! reldiff                 diffusivitiy of species relative to diff. of H2O
@@ -330,8 +339,12 @@ module com_mod
   real :: tth(0:nxmax-1,0:nymax-1,nuvzmax,2)
   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
   ! uupol,vvpol [m/s]    wind components in polar stereographic projection
@@ -345,6 +358,10 @@ module com_mod
   !      cloud, no precipitation      1
   !      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
 
   ! 2d fields
@@ -661,7 +678,6 @@ module com_mod
   ! ohreact             OH reaction rate
 
 
-
   !********************
   ! Random number field
   !********************
@@ -670,5 +686,12 @@ module com_mod
 
   ! 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

src/concoutput.f90

@@ -107,8 +107,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
   call caldate(jul,jjjjmmdd,ihmmss)
   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,
   ! for backward simulations, output fields have dimension MAXPOINT.
@@ -157,7 +158,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
         xl=outlon0+real(ix)*dxout
         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)
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &

src/concoutput_surf.f90

+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
+       drygridtotalunc)
+  !                        i     i          o             o
+  !       o
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the concentration grid and the receptor concentrations.      *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     24 May 1995                                                            *
+  !                                                                            *
+  !     13 April 1999, Major update: if output size is smaller, dump output    *
+  !                    in sparse matrix format; additional output of           *
+  !                    uncertainty                                             *
+  !                                                                            *
+  !     05 April 2000, Major update: output of age classes; output for backward*
+  !                    runs is time spent in grid cell times total mass of     *
+  !                    species.                                                *
+  !                                                                            *
+  !     17 February 2002, Appropriate dimensions for backward and forward runs *
+  !                       are now specified in file par_mod                    *
+  !                                                                            *
+  !     June 2006, write grid in sparse matrix with a single write command     *
+  !                in order to save disk space                                 *
+  !                                                                            *
+  !     2008 new sparse matrix format                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! outnum          number of samples                                          *
+  ! ncells          number of cells with non-zero concentrations               *
+  ! sparse          .true. if in sparse matrix format, else .false.            *
+  ! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use unc_mod
+  use point_mod
+  use outg_mod
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
+  integer :: sp_count_i,sp_count_r
+  real :: sp_fact
+  real :: outnum,densityoutrecept(maxreceptor),xl,yl
+
+  !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
+  !    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
+  !    +    maxageclass)
+  !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
+  !    +       maxageclass)
+  !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
+  !    +       maxpointspec_act,maxageclass)
+  !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
+  !    +       maxpointspec_act,maxageclass),
+  !    +     drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+  !    +     maxpointspec_act,maxageclass),
+  !    +     wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+  !    +     maxpointspec_act,maxageclass)
+  !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
+  !real sparse_dump_r(numxgrid*numygrid*numzgrid)
+  !integer sparse_dump_i(numxgrid*numygrid*numzgrid)
+
+  !real sparse_dump_u(numxgrid*numygrid*numzgrid)
+  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
+  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+  real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
+  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+  real,parameter :: weightair=28.97
+  logical :: sp_zer
+  character :: adate*8,atime*6
+  character(len=3) :: anspec
+
+
+  if (verbosity.eq.1) then
+     print*,'inside concoutput_surf '
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+  !write(unitdates,'(a)') adate//atime
+
+    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,
+  ! for backward simulations, output fields have dimension MAXPOINT.
+  ! Thus, make loops either about nspec, or about numpoint
+  !*****************************************************************
+
+
+  if (ldirect.eq.1) then
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=1
+      end do
+    end do
+  else
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=xmass(kp,ks)
+      end do
+    end do
+  endif
+
+
+  if (verbosity.eq.1) then
+     print*,'concoutput_surf 2'
+     CALL SYSTEM_CLOCK(count_clock)
+     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
+  endif
+
+  !*******************************************************************
+  ! Compute air density: sufficiently accurate to take it
+  ! from coarse grid at some time