Local branch for espen, working with OpenMP version of readwind

options/COMMAND

 &COMMAND
  LDIRECT=          1,
- IBDATE=   20140911,
+ IBDATE=   20070911,
  IBTIME=     050000,
- IEDATE=   20140911,
- IETIME=     230000,
+ IEDATE=   20070912,
+ IETIME=     050000,
  LOUTSTEP=       3600,
  LOUTAVER=       3600,
  LOUTSAMPLE=        900,

options/COMMAND.namelist

+&COMMAND
+ LDIRECT=          1,
+ IBDATE=   20140911,
+ IBTIME=     050000,
+ IEDATE=   20140911,
+ IETIME=     150000,
+ LOUTSTEP=       3600,
+ LOUTAVER=       3600,
+ LOUTSAMPLE=        900,
+ ITSPLIT=  99997200,
+ LSYNCTIME=        900,
+ CTL= -5.0000000    ,
+ IFINE=          4,
+ IOUT=          1,
+ IPOUT=          0,
+ LSUBGRID=          1,
+ LCONVECTION=          1,
+ LAGESPECTRA=          0,
+ IPIN=          0,
+ IOUTPUTFOREACHRELEASE=          0,
+ IFLUX=          0,
+ MDOMAINFILL=          0,
+ IND_SOURCE=          1,
+ IND_RECEPTOR=          1,
+ MQUASILAG=          0,
+ NESTED_OUTPUT=          1,
+ LINIT_COND=          0,
+ SURF_ONLY=          0,
+ CBLFLAG=          0,
+ OHFIELDS_PATH=          "../../flexin/",
+ /

options/COMMAND.oldformat

+********************************************************************************
+*                                                                              *
+*      Input file for the Lagrangian particle dispersion model FLEXPART        *
+*                           Please select your options                         *
+*                                                                              *
+********************************************************************************
+
+1. __                3X, I2
+    1       
+   LDIRECT           1 FOR FORWARD SIMULATION, -1 FOR BACKWARD SIMULATION
+
+2. ________ ______   3X, I8, 1X, I6
+   20140911 050000
+   YYYYMMDD HHMISS   BEGINNING DATE OF SIMULATION
+
+3. ________ ______   3X, I8, 1X, I6
+   20140911 150000
+   YYYYMMDD HHMISS   ENDING DATE OF SIMULATION
+
+4. _____             3X, I5
+    3600 
+   SSSSS             OUTPUT EVERY SSSSS SECONDS
+
+5. _____             3X, I5
+    3600 
+   SSSSS             TIME AVERAGE OF OUTPUT (IN SSSSS SECONDS)
+
+6. _____             3X, I5
+     900 
+   SSSSS             SAMPLING RATE OF OUTPUT (IN SSSSS SECONDS)
+
+7. _________         3X, I9
+   999999999
+   SSSSSSSSS         TIME CONSTANT FOR PARTICLE SPLITTING (IN SECONDS)
+
+8. _____             3X, I5
+     900 
+   SSSSS             SYNCHRONISATION INTERVAL OF FLEXPART (IN SECONDS)
+
+9.  ---.--           4X, F6.4
+     -5.0
+    CTL              FACTOR, BY WHICH TIME STEP MUST BE SMALLER THAN TL
+
+10. ---              4X, I3
+      4
+    IFINE            DECREASE OF TIME STEP FOR VERTICAL MOTION BY FACTOR IFINE
+
+11. -                4X, I1
+    1  
+    IOUT             1 CONCENTRATION (RESIDENCE TIME FOR BACKWARD RUNS) OUTPUT, 2 MIXING RATIO OUTPUT, 3 BOTH,4 PLUME TRAJECT., 5=1+4
+
+12. -                4X, I1
+    0  
+    IPOUT            PARTICLE DUMP: 0 NO, 1 EVERY OUTPUT INTERVAL, 2 ONLY AT END
+
+13. _                4X, I1
+    1
+    LSUBGRID         SUBGRID TERRAIN EFFECT PARAMETERIZATION: 1 YES, 0 NO
+
+14. _                4X, I1
+    1
+    LCONVECTION      CONVECTION: 1 YES, 0 NO
+
+15. _                4X, I1
+    0
+    LAGESPECTRA      AGE SPECTRA: 1 YES, 0 NO
+
+16. _                4X, I1
+    0
+    IPIN             CONTINUE SIMULATION WITH DUMPED PARTICLE DATA: 1 YES, 0 NO
+
+17. _                
+    0                4X,I1 
+    IOFR             IOUTPUTFOREACHREL CREATE AN OUPUT FILE FOR EACH RELEASE LOCATION: 1 YES, 0 NO
+
+18. _                4X, I1
+    0
+    IFLUX            CALCULATE FLUXES: 1 YES, 0 NO
+
+19. _                4X, I1
+    0
+    MDOMAINFILL      DOMAIN-FILLING TRAJECTORY OPTION: 1 YES, 0 NO, 2 STRAT. O3 TRACER
+
+20. _                4X, I1
+    1
+    IND_SOURCE       1=MASS UNIT , 2=MASS MIXING RATIO UNIT 
+
+21. _                4X, I1
+    1
+    IND_RECEPTOR     1=MASS UNIT , 2=MASS MIXING RATIO UNIT 
+
+22. _                4X, I1
+    0
+    MQUASILAG        QUASILAGRANGIAN MODE TO TRACK INDIVIDUAL PARTICLES: 1 YES, 0 NO
+
+23. _                4X, I1
+    1
+    NESTED_OUTPUT    SHALL NESTED OUTPUT BE USED? 1 YES, 0 NO
+
+24. _                4X, I1
+    0
+    LINIT_COND       INITIAL COND. FOR BW RUNS: 0=NO,1=MASS UNIT,2=MASS MIXING RATIO UNIT
+
+25. _                4X, I1
+    0
+    SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
+
+
+1. Simulation direction, 1 for forward, -1 for backward in time 
+	(consult Seibert and Frank, 2004 for backward runs)
+
+2. Beginning date and time of simulation. Must be given in format
+   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,
+   MI is MINUTE and SS is SECOND. Current version utilizes UTC.
+
+3. Ending date and time of simulation. Same format as 3.
+
+4. Average concentrations are calculated every SSSSS seconds.
+
+5. The average concentrations are time averages of SSSSS seconds
+   duration. If SSSSS is 0, instantaneous concentrations are outputted.
+
+6. The concentrations are sampled every SSSSS seconds to calculate the time
+   average concentration. This period must be shorter than the averaging time.
+
+7. Time constant for particle splitting. Particles are split into two
+   after SSSSS seconds, 2xSSSSS seconds, 4xSSSSS seconds, and so on.
+
+8. All processes are synchronized with this time interval (lsynctime).
+   Therefore, all other time constants must be multiples of this value.
+   Output interval and time average of output must be at least twice lsynctime.
+
+9. CTL must be >1 for time steps shorter than the  Lagrangian time scale
+   If CTL<0, a purely random walk simulation is done
+
+10.IFINE=Reduction factor for time step used for vertical wind
+
+11.IOUT determines how the output shall be made: concentration
+   (ng/m3, Bq/m3), mixing ratio (pptv), or both, or plume trajectory mode,
+   or concentration + plume trajectory mode.
+   In plume trajectory mode, output is in the form of average trajectories.
+
+12.IPOUT determines whether particle positions are outputted (in addition
+   to the gridded concentrations or mixing ratios) or not.
+   0=no output, 1 output every output interval, 2 only at end of the
+   simulation
+
+13.Switch on/off subgridscale terrain parameterization (increase of
+   mixing heights due to subgridscale orographic variations)
+
+14.Switch on/off the convection parameterization
+
+15.Switch on/off the calculation of age spectra: if yes, the file AGECLASSES
+   must be available
+
+16. If IPIN=1, a file "partposit_end" from a previous run must be available in
+    the output directory. Particle positions are read in and previous simulation
+    is continued. If IPIN=0, no particles from a previous run are used
+
+17. IF IOUTPUTFOREACHRELEASE is set to 1, one output field for each location
+    in the RLEASE file is created. For backward calculation this should be
+    set to 1. For forward calculation both possibilities are applicable.
+
+18. If IFLUX is set to 1, fluxes of each species through each of the output
+    boxes are calculated. Six fluxes, corresponding to northward, southward,
+    eastward, westward, upward and downward are calculated for each grid cell of
+    the output grid. The control surfaces are placed in the middle of each
+    output grid cell. If IFLUX is set to 0, no fluxes are determined.
+
+19. If MDOMAINFILL is set to 1, the first box specified in file RELEASES is used
+    as the domain where domain-filling trajectory calculations are to be done.
+    Particles are initialized uniformly distributed (according to the air mass
+    distribution) in that domain at the beginning of the simulation, and are
+    created at the boundaries throughout the simulation period.
+
+20. IND_SOURCE switches between different units for concentrations at the source
+    NOTE that in backward simulations the release of computational particles 
+    takes place at the "receptor" and the sampling of particles at the "source".
+          1=mass units (for bwd-runs = concentration)
+          2=mass mixing ratio units 
+21. IND_RECEPTOR switches between different units for concentrations at the receptor
+          1=mass units (concentrations)
+          2=mass mixing ratio units 
+
+22. MQUASILAG indicates whether particles shall be numbered consecutively (1) or
+    with their release location number (0). The first option allows tracking of
+    individual particles using the partposit output files
+
+23. NESTED_OUTPUT decides whether model output shall be made also for a nested
+    output field (normally with higher resolution)
+
+24. LINIT_COND determines whether, for backward runs only, the sensitivity to initial
+    conditions shall be calculated and written to output files
+    0=no output, 1 or 2 determines in which units the initial conditions are provided.
+
+25. SURF_ONLY: When set to 1, concentration/emission sensitivity is written out only
+    for the surface layer; useful for instance when only footprint emission sensitivity is needed
+    but initial conditions are needed on a full 3-D grid
+

options/COMMAND.reference

@@ -105,10 +105,6 @@
     0
     SURF_ONLY        IF THIS IS SET TO 1, OUTPUT IS WRITTEN ONLY OUT FOR LOWEST LAYER
 
-26. _                4X, I1
-    0
-    CBLFLAG          IF THIS IS SET TO 1, USES NEW CONVECTIVE PBL ROUTINES
-
 
 1. Simulation direction, 1 for forward, -1 for backward in time 
 	(consult Seibert and Frank, 2004 for backward runs)
@@ -201,7 +197,3 @@
     for the surface layer; useful for instance when only footprint emission sensitivity is needed
     but initial conditions are needed on a full 3-D grid
 
-26. CBLFLAG: When set to 1, a skewed rather than Gaussian turbulence in the 
-    convective PBL is used. If selected, this option requires shorter time
-    steps (typically values of CTL= 10 and IFINE= 10), and is therefore more 
-    computationally demanding.

options/OUTGRID_NEST

@@ -6,25 +6,25 @@
 ********************************************************************************
 
 1.  ------.----       4X,F11.4
-       -85.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+       -74.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
     OUTLONLEFT        (left boundary of the first grid cell - not its centre)
 
 2.  ------.----       4X,F11.4
-        30.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+        37.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
     OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
 
 3.  -----             4X,I5
-      220             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
+     1000             NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1)
     NUMXGRID
 
 4.  -----             4X,I5
-      120             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
+      200             NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1)
     NUMYGRID
 
 5.  ------.-----      4X,F12.5
-         0.25000      GRID DISTANCE IN X DIRECTION
+         0.01000      GRID DISTANCE IN X DIRECTION
     DXOUTLON
 
 6.  ------.-----      4X,F12.5
-         0.25000      GRID DISTANCE IN Y DIRECTION
+         0.01000      GRID DISTANCE IN Y DIRECTION
     DYOUTLAT

options/RELEASES

@@ -3,34 +3,42 @@
  SPECNUM_REL=         108,
  /
 &RELEASE
- IDATE1=   20140911,
+ IDATE1=   20070911,
  ITIME1=     070000,
- IDATE2=   20140911,
+ IDATE2=   20070911,
  ITIME2=     080000,
- LON1=  -75.250    ,
- LON2=  -75.250    ,
- LAT1=  38.000    ,
- LAT2=  38.000    ,
- Z1=  1500.000000    ,
- Z2=  3000.000000    ,
+ LON1=  -75.255    ,
+ LON2=  -75.255    ,
+ LAT1=  38.005    ,
+ LAT2=  38.005    ,
+ ! Z1=  1500.000000    ,
+ ! Z2=  2500.000000    ,
+ Z1=  2500.000000    ,
+ Z2=  2500.000000    ,
  ZKIND=     1,
  MASS=  2.0000E8    ,
- PARTS=   2000000
+ PARTS=   10000000
  COMMENT="TEST1                          ",
  /
-&RELEASE
- IDATE1=   20140911,
- ITIME1=     070000,
- IDATE2=   20140911,
- ITIME2=     080000,
- LON1=  -79.250   ,
- LON2=  -79.250   ,
- LAT1=  38.100    ,
- LAT2=  38.100    ,
- Z1=  100.000000    ,
- Z2=  1000.000000    ,
- ZKIND=     1,
- MASS=  2.0000E8    ,
- PARTS=   2000000
- COMMENT="TEST2                          ",
- /
+! &RELEASE
+! IDATE1=   20100911,
+! ITIME1=     070000,
+! IDATE2=   20100911,
+! ITIME2=     080000,
+! LON1=  -79.250   ,
+! LON2=  -79.250   ,
+! LAT1=  38.100    ,
+! LAT2=  38.100    ,
+! Z1=  1500.000000    ,
+! Z2=  2500.000000    ,
+! ZKIND=     1,
+! MASS=  2.0000E8    ,
+! PARTS=   2000000
+! COMMENT="TEST2                          ",
+!  /
+
+
+
+
+
+

options/SPECIES/SPECIES_108

@@ -10,15 +10,15 @@ SO2_TEST           Tracer name
 1.0                Below cloud scavenging - B
 0.9                In-cloud scavenging - Ai 
 0.1                In-cloud scavenging - Bi 
--9.9               Dry deposition (gases) - D
-                   Dry deposition (gases) - Henrys const.
-                   Dry deposition (gases) - f0 (reactivity)
-1.9E+03            Dry deposition (particles) - rho
-4.0E-07            Dry deposition (particles) - dquer
-3.0E-01            Dry deposition (particles) - dsig
+2.0               Dry deposition (gases) - D
+1.0E+05            Dry deposition (gases) - Henrys const.
+0.0                Dry deposition (gases) - f0 (reactivity)
+-1.9E+03            Dry deposition (particles) - rho
+-4.0E-07            Dry deposition (particles) - dquer
+-3.0E-01            Dry deposition (particles) - dsig
 -9.99              Alternative: dry deposition velocity
  64.00             molweight
-1.35E-14           OH Reaction rate - C [cm^3/molecule/sec]
+-1.35E-14           OH Reaction rate - C [cm^3/molecule/sec]
 120.0              OH Reaction rate - D [K]
 2.0                OH Reaction rate - N (no unit)
 -9                 number of associated specias (neg. none)

src/com_mod.f90

@@ -731,7 +731,7 @@ module com_mod
   integer :: info_flag=0
   integer :: count_clock, count_clock0,  count_rate, count_max
   real    :: tins
-  logical, parameter :: nmlout=.false.
+  logical, parameter :: nmlout=.true.
 
   ! These variables are used to avoid having separate versions of
   ! files in cases where differences with MPI version is minor (eso)

src/concoutput.f90

@@ -117,7 +117,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 ! Overwrite existing dates file on first call, later append to it
 ! This fixes a bug where the dates file kept growing across multiple runs
 
-! If 'dates' file exists, make a backup
+! If 'dates' file exists in output directory, make a backup
   inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
   if (ldates_file.and.init) then
     open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
@@ -128,7 +128,6 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
       read(unitdates, '(a)', iostat=ierr) dates_char
       if (ierr.ne.0) exit
       write(unit=unittmp, fmt='(a)', iostat=ierr, advance='yes') trim(dates_char)
-!      if (ierr.ne.0) write(*,*) "Write error, ", ierr
     end do
     close(unit=unitdates)
     close(unit=unittmp)

src/makefile

@@ -38,7 +38,7 @@ FLEXPART-GFS 		= FP_gfs_gfortran
 FLEXPART-GFS-MPI 	= FP_gfs_MPI
 
 ifeq ($(gcc), 4.9)
-# Compiled libraries under user user ~flexpart, gfortran v4.9
+# Compiled libraries under users ~flexpart (or ~espen), gfortran v4.9
 	ROOT_DIR = /homevip/flexpart/
 #	ROOT_DIR = /homevip/espen/
 
@@ -69,7 +69,7 @@ O_LEV_DBG = 0 # [0,g]
 ## LIBRARIES
 LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff -llapack  # -fopenmp # -llapack
 
-FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) -mtune=native $(FUSER) # -march=native
+FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -march=native
 
 DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Warray-bounds  -Wall -fcheck=all $(FUSER)  #  -ffpe-trap=invalid,overflow,denormal,underflow,zero -fdump-core
 
@@ -419,7 +419,7 @@ timemanager.o: com_mod.o flux_mod.o netcdf_output_mod.o oh_mod.o outg_mod.o \
 	par_mod.o point_mod.o unc_mod.o xmass_mod.o
 timemanager_mpi.o: com_mod.o flux_mod.o mpi_mod.o oh_mod.o outg_mod.o \
 	par_mod.o point_mod.o unc_mod.o xmass_mod.o
-verttransform.o: cmapf_mod.o com_mod.o par_mod.o
+verttransform.o: cmapf_mod.o com_mod.o par_mod.o 
 verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
 verttransform_nests.o: com_mod.o par_mod.o
 wetdepo.o: com_mod.o par_mod.o point_mod.o

src/mpi_mod.f90

@@ -104,7 +104,7 @@ module mpi_mod
 ! true if only using synchronous MPI send/recv:
 ! If setting this to .false., numwfmem must be set to 3
 !===============================================================================
-  logical :: lmp_sync=.true. 
+  logical :: lmp_sync=.false. 
 !===============================================================================
 
 ! mp_dbg_mode       Used for debugging MPI.
@@ -118,7 +118,7 @@ module mpi_mod
   logical, parameter :: mp_dev_mode = .false.
   logical, parameter :: mp_dbg_out = .false.
   logical, parameter :: mp_time_barrier=.true.
-  logical, parameter :: mp_measure_time=.false.
+  logical, parameter :: mp_measure_time=.true.
   logical, parameter :: mp_exact_numpart=.true.
 
 ! for measuring CPU/Wall time
@@ -213,7 +213,8 @@ contains
       end if
       lmp_sync=.true. ! :DBG: eso fix this...
     end if
-! TODO: Add warnings for unimplemented flexpart features
+
+! TODO: Add more warnings for unimplemented flexpart features
 
 ! Set ID of process that calls getfield/readwind. 
 ! Using the last process in the group ensures statistical identical results
@@ -433,7 +434,7 @@ contains
 !***********************************************************************
 ! Distribute particle variables from pid0 to all processes.
 ! Called from timemanager
-!
+! *NOT IN USE* at the moment, but can be useful for debugging
 !
 !***********************************************************************
     use com_mod
@@ -587,14 +588,14 @@ contains
 ! This can be called from timemanager_mpi, before caclulating                  *
 ! concentrations/writing output grids.                                         *
 !                                                                              *
-! Currently not in use, as each process calculates concentrations              *
-! separately.                                                                  *
+! Currently *not in use*, as each process calculates concentrations            *
+! separately, but can be useful for debugging                                  *
 !                                                                              *
 ! To use this routine (together with mpif_tm_send_vars) to transfer data       *
 ! to the MPI root process (useful for debugging), insert code like this        *
 ! (in timemanager_mpi.f90)                                                     *
 !                                                                              *
-!      if (lroot) tot_numpart=numpart ! root stores total numpart            *
+!      if (lroot) tot_numpart=numpart ! root stores total numpart              *
 !      call mpif_tm_send_dims                                                  *
 !      if (numpart>1) then                                                     *
 !        call mpif_tm_send_vars                                                *
@@ -604,12 +605,12 @@ contains
 !    region:                                                                   *
 !                                                                              *
 !      if (numpart>0) then                                                     *
-!          if (lroot) numpart = tot_numpart                                  *
+!          if (lroot) numpart = tot_numpart                                    *
 !       call mpif_tm_collect_vars                                              *
 !      end if                                                                  *
 !                                                                              *
-!    Then a section that begins with "if (lroot) ..." will behave like       *
-!    serial flexpart                                                           *
+! Then a section that begins with "if (lroot) ..." will behave like            *
+! serial flexpart                                                              *
 !                                                                              *
 !*******************************************************************************
     use com_mod !, only: numpart, nspec, itra1, npoint, nclass
@@ -2065,4 +2066,78 @@ contains
     end subroutine write_data_dbg
 
 
+  subroutine set_fields_synthetic()
+!*******************************************************************************
+! DESCRIPTION
+!   Set all meteorological fields to synthetic (usually constant/homogeneous)
+!   values.
+!   Used for validation and error-checking
+!
+! NOTE
+!   This version uses asynchronious communications.
+!
+! VARIABLES
+!   
+!
+! TODO
+!
+!*******************************************************************************
+    use com_mod
+
+    implicit none
+
+    integer, parameter :: li=1, ui=3 ! wfmem indices (i.e, operate on entire field)
+
+
+! Variables transferred at initialization only
+!*********************************************
+!    readclouds=readclouds_
+    oro(:,:)=0.0
+    excessoro(:,:)=0.0
+    lsm(:,:)=0.0