Changeset d7935de in flexpart.git

Timestamp:

Sep 11, 2018, 6:06:31 PM (6 years ago)

Author:

pesei <petra seibert at univie ac at>

Branches:

univie

Children:

93786a1

Parents:

34f1452

Message:

modify most input read subroutines

changed some variable names (mostly for I-N reasons)
includes two names appearing also in timemanager, com_mod
corrected a few mistakes
simplified some parts of code
changed options/RELEASES which is in nml fmt correspondingly

Files:

• options/RELEASES (modified) (2 diffs)
• src/advance.f90 (modified) (2 diffs)
• src/com_mod.f90 (modified) (4 diffs)
• src/initialize.f90 (modified) (1 diff)
• src/readOHfield.f90 (modified) (1 diff)
• src/readageclasses.f90 (modified) (5 diffs)
• src/readcommand.f90 (modified) (15 diffs)
• src/readoutgrid.f90 (modified) (6 diffs)
• src/readoutgrid_nest.f90 (modified) (6 diffs)
• src/readreceptors.f90 (modified) (9 diffs)
• src/readreleases.f90 (modified) (27 diffs)
• src/timemanager.f90 (modified) (1 diff)
• src/timemanager_mpi.f90 (modified) (1 diff)

options/RELEASES

@@ -1 +1 @@
 &RELEASES_CTRL
  NSPEC=          1,
- SPECNUM_REL=        24,
+ NUMSPEC_REL=        24,
  /
 &RELEASE
@@ -8 +8 @@
  IDATE2=   20170102,
  ITIME2=     090000,
- LON1=    0.000    ,
- LON2=    0.000    ,
- LAT1=   0.000    ,
- LAT2=   0.000    ,
+ RLON1=    0.000    ,
+ RLON2=    0.000    ,
+ RLAT1=   0.000    ,
+ RLAT2=   0.000    ,
  Z1=    50.000       ,
  Z2=    50.000       ,
- ZKIND=     1,
- MASS=  1.0000E8    ,
- PARTS=   10000 
+ IKINDZ=     1,
+ RMASS=  1.0000E8    ,
+ NPARTS=   10000 
  COMMENT="TEST1",
  /

src/advance.f90

@@ -421 +421 @@
   !*************************************************************
   !************** CBL options added by mc see routine cblf90 ***
-            if (cblflag.eq.1) then  !modified by mc
+            if (iflagcbl.eq.1) then  !modified by mc
               if (-h/ol.gt.5) then  !modified by mc
               !if (ol.lt.0.) then   !modified by mc  
@@ -515 +515 @@
 
       end do
-      if (cblflag.ne.1) nrand=nrand+i
+      if (iflagcbl.ne.1) nrand=nrand+i
 
   ! Determine time step for next integration

src/com_mod.f90

@@ -26 +26 @@
   integer :: length(numpath+2*maxnests)
   character(len=256) :: pathfile, flexversion, flexversion_major, arg1, arg2
-  character(len=256) :: ohfields_path
+  character(len=256) :: path_ohfields
   
   ! path                    path names needed for trajectory model
@@ -34 +34 @@
   ! flexversion_major       version of flexpart (major version number)
   ! arg                     input arguments from launch at command line
-  ! ohfields_path           path to binary files for OH fields
+  ! path_ohfields           path to binary files for OH fields
 
   !********************************************************
@@ -73 +73 @@
   integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
   logical :: turbswitch
-  integer :: cblflag !added by mc for cbl
+  integer :: iflagcbl !added by mc for cbl
 
   ! ctl      factor, by which time step must be smaller than Lagrangian time scale
@@ -104 +104 @@
 
   ! ind_rel and ind_samp  are used within the code to change between mass and mass-mix (see readcommand.f)
-  ! cblflag !: 1 activate cbl skewed pdf routines with bi-gaussina pdf whan OL<0 added by mc
+  ! iflagcbl !: 1 activate cbl skewed pdf routines with bi-gaussina pdf whan OL<0 added by mc
 
 

src/initialize.f90

@@ -158 +158 @@
     if (.not.turbswitch) then     ! modified by mc
       wp=wp*sigw
-    else if (cblflag.eq.1) then   ! modified by mc
+    else if (iflagcbl.eq.1) then   ! modified by mc
       if(-h/ol.gt.5) then
 !if (ol.lt.0.) then

src/readOHfield.f90

@@ -60 +60 @@
 
 
-  open(unitOH,file=trim(ohfields_path) &
+  open(unitOH,file=trim(path_ohfields) &
        //'OH_FIELDS/OH_variables.bin',status='old', &
        form='UNFORMATTED', iostat=ierr, convert='little_endian')
 
   if(ierr.ne.0) then
-    write(*,*) 'Cannot read binary OH fields in ',trim(ohfields_path)//'OH_FIELDS/OH_variables.bin'
+    write(*,*) 'Cannot read binary OH fields in ',trim(path_ohfields)//'OH_FIELDS/OH_variables.bin'
     stop
   endif

src/readageclasses.f90

@@ -29 +29 @@
   !     Author: A. Stohl                                                       *
   !     20 March 2000                                                          *
+  !
   !     HSO, 1 July 2014                                                       *
   !     Added optional namelist input                                          *
+  !                                                                            *
+  !     PS, 6/2015-9/2018 some variable names changed as in readreleases.f90   *
+  !       catch nageclass>maxage properly                                      *
   !                                                                            *
   !*****************************************************************************
@@ -47 +51 @@
   integer :: i
 
-  ! namelist help variables
-  integer :: readerror
+  ! namelist aux variables
+  integer :: ios
 
   ! namelist declaration
-  namelist /ageclass/ &
+  namelist /nml_ageclass/ &
     nageclass, &
     lage
@@ -71 +75 @@
   !************************************************
 
-  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
+  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
+    form='formatted',status='old',err=999)
 
   ! try to read in as a namelist
-  read(unitageclasses,ageclass,iostat=readerror)
+  read(unitageclasses, nml_ageclass, iostat=ios)
   close(unitageclasses)
 
-  if ((nageclass.lt.0).or.(readerror.ne.0)) then
-    open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',status='old',err=999)
-    do i=1,13
-      read(unitageclasses,*)
-    end do
-    read(unitageclasses,*) nageclass
-    read(unitageclasses,*) lage(1)
-    do i=2,nageclass
-      read(unitageclasses,*) lage(i)
+  if (ios .ne. 0) then ! failed to read nml, assume simple text
+    open(unitageclasses,file=path(1)(1:length(1))// &
+      'AGECLASSES',status='old',err=999)
+    call skplin(13,unitageclasses)
+    read(unitageclasses,*) nageclass ! number of classes
+    if (nageclass.gt.maxageclass) goto 1001
+    do i=1,nageclass
+      read(unitageclasses,*) lage(i) ! max age per classes
     end do
     close(unitageclasses)
   endif
 
+  if (nageclass.lt.1) goto 1002
+
   ! write ageclasses file in namelist format to output directory if requested
   if (nmlout.and.lroot) then
-    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist',err=1000)
-    write(unitageclasses,nml=ageclass)
+    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist', &
+      err=1000)
+    write(unitageclasses,nml=nml_ageclass)
     close(unitageclasses)
-  endif
-
-  if (nageclass.gt.maxageclass) then
-    write(*,*) ' #### FLEXPART MODEL ERROR! NUMBER OF AGE     #### '
-    write(*,*) ' #### CLASSES GREATER THAN MAXIMUM ALLOWED.   #### '
-    write(*,*) ' #### CHANGE SETTINGS IN FILE AGECLASSES OR   #### '
-    write(*,*) ' #### RECOMPILE WITH LARGER MAXAGECLASS IN    #### '
-    write(*,*) ' #### FILE PAR_MOD.                        #### '
-    stop
   endif
 
@@ -126 +124 @@
 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
+  write(*,'(a)') trim(path(1))
   stop
 
@@ -134 +132 @@
   stop
 
+1001 continue
+    write(*,*) ' #### FLEXPART MODEL ERROR! NUMBER OF AGE     #### '
+    write(*,*) ' #### CLASSES GREATER THAN MAXIMUM ALLOWED.   #### '
+    write(*,*) ' #### CHANGE SETTINGS IN FILE AGECLASSES OR   #### '
+    write(*,*) ' #### RECOMPILE WITH LARGER MAXAGECLASS IN    #### '
+    write(*,*) ' #### FILE PAR_MOD.                           #### '
+    stop
+
+1002 continue
+    write(*,*) ' #### FLEXPART MODEL ERROR! NUMBER OF AGE     #### '
+    write(*,*) ' #### CLASSES < 1                             #### '
+    write(*,*) ' #### CHANGE SETTINGS IN FILE AGECLASSES      #### '
+    stop
 
 end subroutine readageclasses

src/readcommand.f90

@@ -33 +33 @@
   !     Unknown, unknown: various                                              *
   !     Petra Seibert, 2018-06-08: improve error msgs                          *
+  !     PS 6/2015: Minor changes in variable names and layout                  *
   !                                                                            *
   !*****************************************************************************
@@ -84 +85 @@
   character(len=50) :: line
   logical :: old
-  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, &
-  lnetcdfout, &
-  surf_only, &
-  cblflag, &
-  ohfields_path
-
-  ! Presetting namelist command
-  ldirect=0
-  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
-  lnetcdfout=0
-  surf_only=0 
-  cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
-  ohfields_path="../../flexin/"
+  integer :: ios, icmdstat
+
+  namelist /nml_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, &
+    lnetcdfout, &
+    surf_only, &
+    iflagcbl, &
+    path_ohfields
+
+! Set default values for namelist
+    ldirect=0
+    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
+    lnetcdfout=0
+    surf_only=0 
+    iflagcbl=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
+    path_ohfields="../../flexin/"
 
   !Af set release-switch
-  WETBKDEP=.false.
-  DRYBKDEP=.false.
-
+    wetbkdep=.false.
+    drybkdep=.false.
+ 
   ! Open the command file and read user options
   ! Namelist input first: try to read as namelist file
@@ -158 +159 @@
     form='formatted',err=999)
 
-  ! try namelist input (default)
-  read(unitcommand,command,iostat=readerror)
+! try namelist input
+  read(unitcommand, nml_command, iostat=ios)
   close(unitcommand)
 
   ! distinguish namelist from fixed text input
-  if ((readerror.ne.0).or.(ldirect.eq.0)) then ! parse as text file format
+
+  if (ios .ne. 0) then ! simple text file format
  
     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)
+    ! 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)
+    read (unitcommand,900) line
     if (index(line,'LDIRECT') .eq. 0) then
       old = .false.
@@ -243 +244 @@
     ! Removed for backwards compatibility.
     ! if (old) call skplin(3,unitcommand)  !added by mc
-    ! read(unitcommand,*) cblflag          !added by mc
+    ! read(unitcommand,*) iflagcbl          !added by mc
 
     close(unitcommand)
@@ -252 +253 @@
   if (nmlout.and.lroot) then
     open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=998)
-    write(unitcommand,nml=command)
+    write(unitcommand,nml=nml_command)
     close(unitcommand)
   endif
@@ -260 +261 @@
   ! Determine how Markov chain is formulated (for w or for w/sigw)
   !***************************************************************
-  if (cblflag.eq.1) then !---- added by mc to properly set parameters for CBL simulations 
+  if (iflagcbl.eq.1) then !added by mc to set parameters for CBL simulations 
     turbswitch=.true.
-    if (lsynctime>maxtl) lsynctime=maxtl  !maxtl defined in com_mod.f90
+    if (lsynctime.gt.maxtl) lsynctime=maxtl  !maxtl defined in com_mod.f90
     if (ctl.lt.5) then
-      print *,'WARNING: CBL flag active the ratio of TLu/dt has been set to 5'
+      print*,'WARNING: CBL flag active; ctl (TLu/dt) has been set to 5'
       ctl=5.
-    end if
-    if (ifine*ctl.lt.50) then
+    endif
+    if (ifine*ctl.lt.50.) then
       ifine=int(50./ctl)+1
-
-      print *,'WARNING: CBL flag active the ratio of TLW/dt was < 50, ifine has been re-set to',ifine
-!pause
+      print *,'WARNING: CBL flag active; ctl (TLW/dt) was < 50,'// &
+        ' ifine has been re-set to', ifine
     endif
-    print *,'WARNING: CBL flag active the ratio of TLW/dt is ',ctl*ifine
-    print *,'WARNING: CBL flag active lsynctime is ',lsynctime
+    print*,'WARNING: CBL flag active; reduced ctl is ',ctl*ifine
+    print*,'WARNING: CBL flag active; lsynctime is ',lsynctime
   else                    !added by mc
+  ! note PS: shouldn't we print some msg as above also in the ntext case?
     if (ctl.ge.0.1) then
       turbswitch=.true.
@@ -286 +287 @@
   ctl=1./ctl
 
-  ! Set the switches required for the various options for input/output units
-  !*************************************************************************
-  !AF Set the switches IND_REL and IND_SAMP for the release and sampling
-  !Af switches for the releasefile:
-  !Af IND_REL =  1 : xmass * rho
-  !Af IND_REL =  0 : xmass * 1
-
-  !Af switches for the conccalcfile:
-  !AF IND_SAMP =  0 : xmass * 1
-  !Af IND_SAMP = -1 : xmass / rho
-
-  !AF IND_SOURCE switches between different units for concentrations at the source
-  !Af   NOTE that in backward simulations the release of computational particles
-  !Af   takes place at the "receptor" and the sampling of p[articles at the "source".
-  !Af          1 = mass units
-  !Af          2 = mass mixing ratio units
-  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
-  !Af          1 = mass units
-  !Af          2 = mass mixing ratio units
-  !            3 = wet deposition in outputfield
-  !            4 = dry deposition in outputfield
+! Set the switches required for the various options for input/output units
+!*************************************************************************
+!AF Set the switches IND_REL and IND_SAMP for the release and sampling
+!Af switches for the releasefile:
+!Af IND_REL =  1 : xmass * rho
+!Af IND_REL =  0 : xmass * 1
+
+!Af switches for the conccalcfile:
+!AF IND_SAMP =  0 : xmass * 1
+!Af IND_SAMP = -1 : xmass / rho
+
+!AF IND_SOURCE switches between different units for concentrations at the source
+!Af   NOTE that in backward simulations the release of computational particles
+!Af   takes place at the "receptor" and the sampling of p[articles at the "source".
+!Af          1 = mass units
+!Af          2 = mass mixing ratio units
+!Af IND_RECEPTOR switches between different units for concentrations at the receptor
+!Af          1 = mass units
+!Af          2 = mass mixing ratio units
+!            3 = wet deposition in outputfield
+!            4 = dry deposition in outputfield
 
   if ( ldirect .eq. 1 ) then  ! FWD-Run
@@ -328 +329 @@
         ind_samp = 0
      endif
+! note PS: why do we suddenly switch to CASE syntax??
+! not helpful.
      select case (ind_receptor)
      case (1)  !  1 .. concentration at receptor
@@ -337 +340 @@
         if (lroot) then
           write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
-          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
+          write(*,*) ' #### Release height is forced to 0 - 20 km    #### '
           write(*,*) ' #### Release is performed above ground lev    #### '
         end if
-         WETBKDEP=.true.
+         wetbkdep=.true.
          allocate(xscav_frac1(maxpart,maxspec))
      case (4)  ! 4 .. dry deposition in outputfield
@@ -349 +352 @@
            write(*,*) ' #### Release is performed above ground lev    #### '
          end if
-         DRYBKDEP=.true.
+         drybkdep=.true.
          allocate(xscav_frac1(maxpart,maxspec))
      end select
@@ -399 +402 @@
   endif
 
-! Check for netcdf output switch
-!*******************************
+!  check for netcdf output switch (use for non-namelist input only!)
   if (iout.ge.8) then
      lnetcdfout = 1
      iout = iout - 8
 #ifndef USE_NCF
-     write(*,*) 'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
-     write(*,*) 'Please recompile with netcdf library (`make [...] ncf=yes`) or use standard output format.'
+     write(*,*) 'ERROR: netcdf output not activated during compile '// &
+       'time but used in COMMAND file!'
+     write(*,*) 'Please recompile with netcdf library (`make [...] ncf=yes`)'//&
+       ' or use standard output format.'
      stop
 #endif
@@ -414 +418 @@
   !**********************************************************************
 
-  if ((iout.lt.1).or.(iout.gt.6)) then
+  if (iout.lt.1 .or. iout.gt.6) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
     write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
@@ -423 +427 @@
 
   !AF check consistency between units and volume mixing ratio
-  if ( ((iout.eq.2).or.(iout.eq.3)).and. &
-       (ind_source.gt.1 .or.ind_receptor.gt.1) ) then
+  if ( (iout.eq.2       .or. iout.eq.3) .and. &
+       (ind_source.gt.1 .or. ind_receptor.gt.1) ) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
     write(*,*) ' #### VOLUME MIXING RATIO ONLY SUPPORTED      #### '
@@ -479 +483 @@
 
   ! For domain-filling trajectories, a plume centroid trajectory makes no sense,
-  ! For backward runs, only residence time output (iout=1) or plume trajectories (iout=4),
-  ! or both (iout=5) makes sense; other output options are "forbidden"
+  ! For backward runs, only residence time output (iout=1) or plume trajectories
+  ! (iout=4), or both (iout=5) makes sense; other output options are "forbidden"
   !*****************************************************************************
 
@@ -649 +653 @@
     path(2)(1:length(2))//'COMMAND.namelist'
   stop 'stopped in readcommand'
+
 999 write(*,900) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"'
   write(*,900)   ' #### CANNOT OPEN '//path(1)(1:length(1))//'COMMAND'
   stop 'stopped in readcommand'
-900 format (a)
+
+900 format(a)
+
 end subroutine readcommand

src/readoutgrid.f90

@@ -29 +29 @@
   !                                                                            *
   !     4 June 1996                                                            *
-  !     HSO, 1 July 2014
-  !     Added optional namelist input
+  !     HSO, 1 July 2014: Add optional namelist input                          *
+  !     PS, 6/2015-9/2018: read regular input with free format                 *
+  !       simplify code and rename some variables                              *
   !                                                                            *
   !*****************************************************************************
@@ -51 +52 @@
   implicit none
 
-  integer :: i,j,stat
-  real :: outhelp,xr,xr1,yr,yr1
+  integer :: i,kz,istat
+  real :: xr,xr1,yr,yr1
   real,parameter :: eps=1.e-4
 
   ! namelist variables
   integer, parameter :: maxoutlev=500
-  integer :: readerror
-  real,allocatable, dimension (:) :: outheights
+  integer :: ios
+  real,allocatable, dimension (:) :: outheights,outaux
+  logical :: lnml
 
   ! declare namelist
-  namelist /outgrid/ &
+  namelist /nml_outgrid/ &
     outlon0,outlat0, &
     numxgrid,numygrid, &
     dxout,dyout, &
-    outheights
-
-  ! allocate large array for reading input
-  allocate(outheights(maxoutlev),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
-
-  ! helps identifying failed namelist input
-  dxout=-1.0
-  outheights=-1.0
+    outheight
+
+! allocate outheights for nml read with max dimension
+  allocate(outheights(maxoutlev),outaux(maxoutlev),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate outheights'
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
+  outheight(:) = -999. ! initialise for later finding #valid levels
+  open(unitoutgrid,file=trim(path(1))//'OUTGRID',status='old',&
+    form='formatted',err=999)
 
   ! try namelist input
-  read(unitoutgrid,outgrid,iostat=readerror)
+  read(unitoutgrid,nml_outgrid,iostat=ios)
   close(unitoutgrid)
 
-  if ((dxout.le.0).or.(readerror.ne.0)) then
-
-    readerror=1
-
-    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
-
+  if (ios .eq. 0) then ! namelist works
+
+    lnml = .true.
+
+  else ! read as regular text
+
+    lnml = .false.
+
+    open(unitoutgrid,file=trim(path(1))//'OUTGRID',status='old',err=999)
     call skplin(5,unitoutgrid)
 
-    ! 1.  Read horizontal grid specifications
+   ! Read horizontal grid specifications
     !****************************************
 
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlon0
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlat0
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numxgrid
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numygrid
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dxout
-    call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dyout
-
-  endif
+    read(unitoutgrid,*) outlon0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,*) outlat0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,*) numxgrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,*) numygrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,*) dxout
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,*) dyout
+
+  endif ! read OUTGRID file
 
   ! Check validity of output grid (shall be within model domain)
@@ -124 +127 @@
     write(*,*) ' #### GRID IS OUTSIDE MODEL DOMAIN. CHANGE    ####'
     write(*,*) ' #### FILE OUTGRID IN DIRECTORY               ####'
-    write(*,'(a)') path(1)(1:length(1))
+    write(*,'(a)') trim(path(1))
     stop
   endif
 
-  ! 2. Count Vertical levels of output grid
-  !****************************************
-
-  if (readerror.ne.0) then
-    j=0
-100 j=j+1
-    do i=1,3
-      read(unitoutgrid,*,end=99)
-    end do
-    read(unitoutgrid,'(4x,f7.1)',end=99) outhelp
-    if (outhelp.eq.0.) goto 99
-    goto 100
-99  numzgrid=j-1
-  else
-    do i=1,maxoutlev
-      if (outheights(i).lt.0) exit
-    end do
-    numzgrid=i-1
-  end if
-
-  allocate(outheight(numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
-  allocate(outheighthalf(numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
-
-  ! 2. Vertical levels of output grid
-  !**********************************
-
-  if (readerror.ne.0) then
-
-    rewind(unitoutgrid)
-    call skplin(29,unitoutgrid)
-
-    do j=1,numzgrid
-      do i=1,3
-        read(unitoutgrid,*)
-      end do
-      read(unitoutgrid,'(4x,f7.1)') outhelp
-      outheight(j)=outhelp
-      outheights(j)=outhelp
-    end do
-    close(unitoutgrid)
-
-  else
-
-    do j=1,numzgrid
-      outheight(j)=outheights(j)
-    end do
-
+! Read (if .not. lmnl) and count vertical levels of output grid 
+!**************************************************************
+
+  do kz = 1,maxoutlev
+    if (lnml) then ! we have read them already
+      if (outheight(kz) .lt. 0.) exit ! 1st nondefined level
+    else
+      call skplin(3,unitoutgrid)
+      read(unitoutgrid,*,end=10) outheight(kz)
+    endif
+  end do
+10 continue  
+
+  numzgrid = kz - 1 ! number of outgrid levels
+
+! allocate the required length only, shuffle data  
+  outaux = outheight ! shuffle
+
+  deallocate(outheights)
+  allocate(outheight(numzgrid),outheighthalf(numzgrid),stat=istat)
+  if (istat .ne. 0) then
+    write(*,*) 'ERROR: could not allocate outheight and outheighthalf'
+    stop 'readoutgrid error'
   endif
 
-  ! write outgrid file in namelist format to output directory if requested
-  if (nmlout.and.lroot) then
-    ! reallocate outheights with actually required dimension for namelist writing
-    deallocate(outheights)
-    allocate(outheights(numzgrid),stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
-
-    do j=1,numzgrid
-      outheights(j)=outheight(j)
-    end do
-
-    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
-    write(unitoutgrid,nml=outgrid)
+  outheight=outaux(1:numzgrid) ! shuffle back
+  deallocate (outaux) 
+
+! write outgrid file in namelist format to output directory if requested
+  if (nmlout) then
+    open(unitoutgrid,file=trim(path(2))//'OUTGRID.namelist',err=1000)
+    write(unitoutgrid,nml=nml_outgrid)
     close(unitoutgrid)
   endif
@@ -198 +169 @@
   !***************************************************************
 
-  do j=2,numzgrid
-    if (outheight(j).le.outheight(j-1)) then
-    write(*,*) ' #### FLEXPART MODEL ERROR! YOUR SPECIFICATION#### '
-    write(*,*) ' #### OF OUTPUT LEVELS IS CORRUPT AT LEVEL    #### '
-    write(*,*) ' #### ',j,'                              #### '
-    write(*,*) ' #### PLEASE MAKE CHANGES IN FILE OUTGRID.    #### '
-    endif
+  do kz=2,numzgrid
+    if (outheight(kz) .le. outheight(kz-1)) goto 998
   end do
 
@@ -210 +176 @@
   !*****************************************************************
 
-  outheighthalf(1)=outheight(1)/2.
-  do j=2,numzgrid
-    outheighthalf(j)=(outheight(j-1)+outheight(j))/2.
+  outheighthalf(1) = 0.5*outheight(1)
+  do kz = 2,numzgrid
+    outheighthalf(kz) = 0.5*(outheight(kz-1)+outheight(kz))
   end do
 
@@ -218 +184 @@
   youtshift=ylat0-outlat0
 
-  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
-  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
-  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
-  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
-  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
+  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate oroout'
+  allocate(area(0:numxgrid-1,0:numygrid-1),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate area'
+  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate volume'
+  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate areaeast'
+  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=istat)
+  if (istat .ne. 0) write(*,*)'ERROR: could not allocate areanorth'
+  
   return
 
+998 continue
+  write(*,*) ' #### FLEXPART MODEL ERROR! YOUR SPECIFICATION#### '
+  write(*,*) ' #### OF OUTPUT LEVELS NOT INCREASING AT LEVEL#### '
+  write(*,*) ' #### ',kz,'                                  #### '
+  write(*,*) ' #### PLEASE MAKE CHANGES IN FILE OUTGRID.    #### '
+  STOP 'readoutgrid error'
 
 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
+  write(*,'(a)') trim(path(1))
   stop
 
 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
+  write(*,'(a)') trim(path(2))
   stop
 

src/readoutgrid_nest.f90

@@ -29 +29 @@
   !                                                                            *
   !     4 June 1996                                                            *
+  !     HSO, 1 July 2014: Add optional namelist input                          *
+  !     PS, 6/2015-9/2018: read regular input with free format                 *
+  !       and rename some variables                                            *
   !                                                                            *
   !*****************************************************************************
@@ -53 +56 @@
   real,parameter :: eps=1.e-4
 
-  integer :: readerror
+  integer :: ios
 
-  ! declare namelist
-  namelist /outgridn/ &
-    outlon0n,outlat0n, &
-    numxgridn,numygridn, &
-    dxoutn,dyoutn
+! declare namelist
+  namelist /nml_outgridn/ &
+    outlon0n,outlat0n,numxgridn,numygridn,dxoutn,dyoutn
 
   ! helps identifying failed namelist input
@@ -67 +68 @@
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',&
+    status='old',err=999)
 
   ! try namelist input
-  read(unitoutgrid,outgridn,iostat=readerror)
+  read(unitoutgrid,nml_outgridn,iostat=ios)
   close(unitoutgrid)
 
-  if ((dxoutn.le.0).or.(readerror.ne.0)) then
+  if (dxoutn.le.0 .or.ios.ne.0) then
 
-    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
+    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',&
+      err=999)
     call skplin(5,unitoutgrid)
 
-    ! 1.  Read horizontal grid specifications
-    !****************************************
+  ! Read horizontal grid specifications
+  ! ***********************************
 
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlon0n
+    read(unitoutgrid,*) outlon0n
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f11.4)') outlat0n
+    read(unitoutgrid,*) outlat0n
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numxgridn
+    read(unitoutgrid,*) numxgridn
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,i5)') numygridn
+    read(unitoutgrid,*) numygridn
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dxoutn
+    read(unitoutgrid,*) dxoutn
     call skplin(3,unitoutgrid)
-    read(unitoutgrid,'(4x,f12.5)') dyoutn
+    read(unitoutgrid,*) dyoutn
 
     close(unitoutgrid)
@@ -99 +102 @@
   ! write outgrid_nest file in namelist format to output directory if requested
   if (nmlout.and.lroot) then
-    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',err=1000)
-    write(unitoutgrid,nml=outgridn)
+    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',&
+      err=1000)
+    write(unitoutgrid,nml=nml_outgridn)
     close(unitoutgrid)
   endif
 
-  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
-  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
-  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
+  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=ios)
+  if (ios.ne.0) write(*,*)'ERROR: could not allocate orooutn'
+  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=ios)
+  if (ios.ne.0) write(*,*)'ERROR: could not allocate arean'
+  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=ios)
+  if (ios.ne.0) write(*,*)'ERROR: could not allocate volumen'
 
   ! Check validity of output grid (shall be within model domain)
@@ -118 +122 @@
   xr1=xlon0+real(nxmin1)*dx
   yr1=ylat0+real(nymin1)*dy
-  if ((outlon0n+eps.lt.xlon0).or.(outlat0n+eps.lt.ylat0) &
-       .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
+  if (outlon0n+eps.lt.xlon0 .or. outlat0n+eps.lt.ylat0 &
+    .or. xr.gt.xr1+eps .or. yr.gt.yr1+eps) then
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
     write(*,*) ' #### NEST IS OUTSIDE MODEL DOMAIN. CHANGE    ####'
     write(*,*) ' #### FILE OUTGRID IN DIRECTORY               ####'
-    write(*,'(a)') path(1)(1:length(1))
+    write(*,'(a)') trim(path(1))
     stop
   endif
@@ -133 +137 @@
 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
+  write(*,'(a)') trim(path(1))
   stop
 
 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
+  write(*,'(a)') trim(path(2))
   stop
 

src/readreceptors.f90

@@ -27 +27 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !     1 August 1996                                                          *
-  !     HSO, 14 August 2013
-  !     Added optional namelist input
+  !     1 August 1996        
+  !                                                  *
+  !     HSO, 14 August 2013: Added optional namelist input
+  !     PS, 2/2015: access= -> position=
+  !     PS, 6/2015: variable names, simplify code
   !                                                                            *
   !*****************************************************************************
@@ -52 +54 @@
   character(len=16) :: receptor
 
-  integer :: readerror
-  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
+  integer :: ios
+  real :: xlon,ylat   ! for namelist input, lon/lat are used instead of x,y
   integer,parameter :: unitreceptorout=2
 
   ! declare namelist
-  namelist /receptors/ &
-    receptor, lon, lat
+  namelist /nml_receptors/ receptor, xlon, ylat
 
-  lon=-999.9
-
-  ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
+!CPS I comment this out - why should we not have receptor output in bwd runs? 
+  ! For backward runs, do not allow receptor output. Thus, set number of
+  ! receptors to zero
   !*****************************************************************************
 
-  if (ldirect.lt.0) then
-    numreceptor=0
-    return
-  endif
+!  if (ldirect.lt.0) then
+!    numreceptor=0
+!    return
+!  endif
 
 
@@ -74 +75 @@
   !************************************************************
 
-  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
+  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',&
+    status='old',err=999)
 
   ! try namelist input
-  read(unitreceptor,receptors,iostat=readerror)
+  read(unitreceptor,nml_receptors,iostat=ios)
 
   ! prepare namelist output if requested
-  if (nmlout.and.lroot) then
-    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',&
-         &access='append',status='replace',err=1000)
-  endif
+  if (nmlout) open(unitreceptorout,file=path(2)(1:length(2))// &
+    'RECEPTORS.namelist', position='append',status='new',err=1000)
 
-  if ((lon.lt.-900).or.(readerror.ne.0)) then
+  if (ios .ne. 0) then ! read as regular text file
 
     close(unitreceptor)
-    open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
+    open(unitreceptor,file=path(1)(1:length(1))// &
+     'RECEPTORS',status='old',err=999)
+ 
     call skplin(5,unitreceptor)
 
@@ -101 +103 @@
     read(unitreceptor,'(4x,a16)',end=99) receptor
     call skplin(3,unitreceptor)
-    read(unitreceptor,'(4x,f11.4)',end=99) x
+      read(unitreceptor,'(4x,f11.4)',end=99) xlon
     call skplin(3,unitreceptor)
-    read(unitreceptor,'(4x,f11.4)',end=99) y
-    if ((x.eq.0.).and.(y.eq.0.).and. &
+      read(unitreceptor,'(4x,f11.4)',end=99) ylat
+      if ((xlon.eq.0.).and.(ylat.eq.0.).and. &
          (receptor.eq.'                ')) then
       j=j-1
@@ -116 +118 @@
     endif
     receptorname(j)=receptor
-    xreceptor(j)=(x-xlon0)/dx       ! transform to grid coordinates
-    yreceptor(j)=(y-ylat0)/dy
-    xm=r_earth*cos(y*pi/180.)*dx/180.*pi
+      xreceptor(j)=(xlon-xlon0)/dx       ! transform to grid coordinates
+      yreceptor(j)=(ylat-ylat0)/dy
+      xm=r_earth*cos(ylat*pi/180.)*dx/180.*pi
     ym=r_earth*dy/180.*pi
     receptorarea(j)=xm*ym
-
-    ! write receptors file in namelist format to output directory if requested
-    if (nmlout.and.lroot) then
-      lon=x
-      lat=y
-      write(unitreceptorout,nml=receptors)
-    endif
+  ! write receptors file in namelist format to output directory if requested
+    if (nmlout.and.lroot) write(unitreceptorout,nml=nml_receptors)
 
     goto 100
@@ -136 +133 @@
 
     j=0
-    do while (readerror.eq.0)
+    do while (ios .eq. 0)
       j=j+1
-      lon=-999.9
-      read(unitreceptor,receptors,iostat=readerror)
-      if ((lon.lt.-900).or.(readerror.ne.0)) then
-        readerror=1
-      else
+      read(unitreceptor,nml_receptors,iostat=ios)
+      if (ios .eq. 0) then
         if (j.gt.maxreceptor) then
           write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
@@ -150 +144 @@
         endif
         receptorname(j)=receptor
-        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
-        yreceptor(j)=(lat-ylat0)/dy
-        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
+        xreceptor(j)=(xlon-xlon0)/dx       ! transform to grid coordinates
+        yreceptor(j)=(ylat-ylat0)/dy
+        xm=r_earth*cos(ylat*pi/180.)*dx/180.*pi
         ym=r_earth*dy/180.*pi
         receptorarea(j)=xm*ym
       endif
 
-      ! write receptors file in namelist format to output directory if requested
-      if (nmlout.and.lroot) then
-        write(unitreceptorout,nml=receptors)
-      endif
+    ! write receptors file in namelist format to output directory if requested
+      if (nmlout.and.lroot) write(unitreceptorout,nml=nml_receptors)
 
     end do
@@ -166 +158 @@
     close(unitreceptor)
 
-  endif
+  endif ! end no-nml / nml bloc
 
-  if (nmlout.and.lroot) then
-    close(unitreceptorout)
-  endif
+  if (nmlout .and. lroot) close(unitreceptorout)
 
   return
@@ -176 +166 @@
 
 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(1)(1:length(1))
+  write(*,*)   ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') trim(path(1))
   stop
 
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,'(a)') path(2)(1:length(2))
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
+  write(*,*)    ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') trim(path(2))
   stop
 

src/readreleases.f90

@@ -33 +33 @@
   !     Update: 29 January 2001                                                *
   !     Release altitude can be either in magl or masl                         *
-  !     HSO, 12 August 2013
-  !     Added optional namelist input
+  !     HSO, 12 August 2013: Added optional namelist input                     *
+  !     PS, 2/2015: access= -> position= (ported from 9.2 review)              *
+  !     PS, 8/2018: check for uninitialised mass (ticket:204),                 *
+  !          rename int/real variables, a bit of code layout improvement       *
   !                                                                            *
   !*****************************************************************************
@@ -74 +76 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,iallocstat,irel,ispc,nsettle
   integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
   real(kind=dp) :: jul1,jul2,julm,juldate
-  real,parameter :: eps2=1.e-9
   character(len=50) :: line
   logical :: old
 
   ! help variables for namelist reading
-  integer :: numpoints, parts, readerror
-  integer*2 :: zkind
+  integer :: numpoints, nparts, ios
+  integer*2 :: ikindz
   integer :: idate1, itime1, idate2, itime2
-  real :: lon1,lon2,lat1,lat2,z1,z2
+  real :: rlon1,rlon2,rlat1,rlat2,z1,z2
   character*40 :: comment
   integer,parameter :: unitreleasesout=2
-  real,allocatable, dimension (:) :: mass
-  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
+  real,allocatable, dimension (:) :: rmass
+  integer,allocatable, dimension (:) :: numspec_rel,numspec_rel2
 
   ! declare namelists
   namelist /releases_ctrl/ &
        nspec, &
-       specnum_rel
+       numspec_rel
 
   namelist /release/ &
        idate1, itime1, &
        idate2, itime2, &
-       lon1, lon2, &
-       lat1, lat2, &
+       rlon1, rlon2, &
+       rlat1, rlat2, &
        z1, z2, &
-       zkind, &
-       mass, &
-       parts, &
+       ikindz, &
+       rmass, &
+       nparts, &
        comment
+! initialiase mass
 
   numpoint=0
 
-  ! allocate with maxspec for first input loop
-  allocate(mass(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-
-  ! presetting namelist releases_ctrl
+! allocate with maxspec for first input loop
+  allocate(rmass(maxspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate rmass'
+  allocate(numspec_rel(maxspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel'
+
+! presetting namelist releases_ctrl
   nspec = -1  ! use negative value to determine failed namelist input
-  specnum_rel = 0
-
-  !sec, read release to find how many releasepoints should be allocated
-  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
-
-  ! check if namelist input provided
-  read(unitreleases,releases_ctrl,iostat=readerror)
-
-  ! prepare namelist output if requested
+  numspec_rel = 0
+
+!SEC: read release to find how many releasepoints should be allocated
+  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
+    form='formatted',err=999)
+
+! prepare namelist output if requested
   if (nmlout.and.lroot) then
-    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='replace',err=1000)
-  endif
-
-  if ((readerror.ne.0).or.(nspec.lt.0)) then
-
-  ! no namelist format, close file and allow reopening in old format
+    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist', &
+      position='append',status='replace',err=1000)
+  endif
+
+! check whether namelist input is available
+  read(unitreleases,releases_ctrl,iostat=ios)
+
+  if (ios .ne. 0 .or. nspec .lt. 0 ) then
+
+  ! not in namelist format, close file and reopen as simple text file
+
     close(unitreleases)
-    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
-
-    readerror=1 ! indicates old format
-
-  ! Check the format of the RELEASES file (either in free format,
-  ! or using a formatted mask)
+    open(unitreleases,file=path(1)(1:length(1))//'RELEASES', status='old', &
+      err=999)
+
+    ios=1 ! indicates old format
+
+  ! Check the format of the RELEASES file
+  ! (either in free format, or using a formatted mask)
   ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
   !**************************************************************************
+! Note PS: something has changed, now Total is checked. Is this compatible??
 
     call skplin(12,unitreleases)
@@ -163 +170 @@
     if (old) call skplin(2,unitreleases)
     do i=1,nspec
-      read(unitreleases,*,err=998) specnum_rel(i)
+      read(unitreleases,*,err=998) numspec_rel(i)
       if (old) call skplin(2,unitreleases)
     end do
@@ -202 +209 @@
 25  numpoint=numpoint-1
 
-  else 
-
-    readerror=0
-    do while (readerror.eq.0) 
+  else ! read namelist input
+
+    ios=0
+    do while (ios.eq.0)  ! loop over release points
       idate1=-1
-      read(unitreleases,release,iostat=readerror)
-      if ((idate1.lt.0).or.(readerror.ne.0)) then
-        readerror=1
+      read(unitreleases,release,iostat=ios)
+      if (idate1.lt.0 .or. ios.ne.0) then
+        ios=1
       else
         numpoint=numpoint+1
       endif
     end do
-    readerror=0
-  endif ! if namelist input
+    ios=0
+
+  endif 
 
   rewind(unitreleases)
@@ -221 +229 @@
   if (nspec.gt.maxspec) goto 994
 
-  ! allocate arrays of matching size for number of species (namelist output)
-  deallocate(mass)
-  allocate(mass(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel2(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
-  specnum_rel2=specnum_rel(1:nspec)
-  deallocate(specnum_rel) 
+! allocate arrays of matching size for number of species (namelist output)
+  deallocate(rmass)
+  allocate(rmass(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate rmass'
+  allocate(numspec_rel2(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel2'
+
+  numspec_rel2 = numspec_rel(1:nspec)
+
+  deallocate(numspec_rel)
   ! eso: BUG, crashes here for nspec=12 and maxspec=6,
   ! TODO: catch error and exit
-  allocate(specnum_rel(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-  specnum_rel=specnum_rel2
-  deallocate(specnum_rel2)
-
-  !allocate memory for numpoint releaspoints
-  allocate(ireleasestart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
-  allocate(ireleaseend(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
-  allocate(xpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
-  allocate(xpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
-  allocate(ypoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
-  allocate(ypoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
-  allocate(zpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
-  allocate(zpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
-  allocate(kindz(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
-  allocate(xmass(numpoint,maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
-  allocate(rho_rel(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
-  allocate(npart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
-  allocate(xmasssave(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
-
-  if (lroot) write (*,*) 'Releasepoints : ', numpoint
-
-  do i=1,numpoint
-    xmasssave(i)=0.
-  end do
+  allocate(numspec_rel(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel'
+
+  numspec_rel = numspec_rel2
+
+  deallocate(numspec_rel2)
+
+! allocate memory for numpoint releaspoints
+  allocate ( &
+    ireleasestart(numpoint), ireleaseend(numpoint),xpoint1(numpoint), &
+    xpoint2(numpoint), ypoint1(numpoint), ypoint2(numpoint), &
+    zpoint1(numpoint), zpoint2(numpoint), kindz(numpoint), &
+    xmass(numpoint,maxspec), rho_rel(numpoint), npart(numpoint), &
+    xmasssave(numpoint), &
+      stat=iallocstat )
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate variables' // &
+    ' dimensioned numpoint =', numpoint
+
+  if (lroot) write (*,*) 'Release points : ', numpoint
+
+  xmasssave(:)=0.
 
   !now save the information
-  DEP=.false.
-  DRYDEP=.false.
-  WETDEP=.false.
-  OHREA=.false.
-  do i=1,maxspec
-    DRYDEPSPEC(i)=.false.
-    WETDEPSPEC(i)=.false.
-  end do
-
-  if (readerror.ne.0) then
+  dep   =.false.
+  drydep=.false.
+  wetdep=.false.
+  ohrea =.false.
+  drydepspec(:)=.false.
+  wetdepspec(:)=.false.
+
+  if (ios.ne.0) then
+
   ! Skip first 11 lines (file header)
   !**********************************
@@ -292 +284 @@
     if (nspec.gt.maxspec) goto 994
     if (old) call skplin(2,unitreleases)
-  endif
-
-  ! namelist output
-  if (nmlout.and.lroot) then
-    write(unitreleasesout,nml=releases_ctrl)
-  endif
-
+
+  endif
+
+! namelist output
+  if (nmlout.and.lroot) write(unitreleasesout,nml=releases_ctrl)
+
+! species loop
   do i=1,nspec
-    if (readerror.ne.0) then
-      read(unitreleases,*,err=998) specnum_rel(i)
+  
+    if (ios.ne.0) then
+      read(unitreleases,*,err=998) numspec_rel(i)
       if (old) call skplin(2,unitreleases)
-      call readspecies(specnum_rel(i),i)
+      call readspecies(numspec_rel(i),i)
     else
-      call readspecies(specnum_rel(i),i)
+      call readspecies(numspec_rel(i),i)
     endif
 
@@ -322 +315 @@
   !*****************
 
-    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
+    if ( (iout.eq.2 .or. iout.eq.3) .and. weightmolar(i).lt.0. ) then
       write(*,*) 'For mixing ratio output, valid molar weight'
       write(*,*) 'must be specified for all simulated species.'
@@ -334 +327 @@
 
     decay(i)=0.693147/decay(i) !conversion half life to decay constant
-
+! Note PS: this should be replace by = alog(2.)/decay(i)
+! keep it for the moment so that results are not affect in test phase of v10
 
   ! Dry deposition of gases
   !************************
 
-    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
+    if (reldiff(i).gt.0.) rm(i)=1./( henry(i)/3000. + 100.*f0(i) )! mesophyll resistance
 
   ! Dry deposition of particles
   !****************************
 
-    vsetaver(i)=0.
-    cunningham(i)=0.
-    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
+    vsetaver(i) = 0.
+    cunningham(i) = 0.
+    dquer(i) = dquer(i)*1.e6 ! Conversion m to um
     if (density(i).gt.0.) then         ! Additional parameters
       call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
       do j=1,ni
-        fract(i,j)=fracth(j)
-        schmi(i,j)=schmih(j)
-        vset(i,j)=vsh(j)
-        cunningham(i)=cunningham(i)+cun*fract(i,j)
-        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
+        fract(i,j) = fracth(j)
+        schmi(i,j) = schmih(j)
+        vset(i,j) = vsh(j)
+        cunningham(i) = cunningham(i) + cun*fract(i,j)
+        vsetaver(i) = vsetaver(i) - vset(i,j)*fract(i,j)
       end do
       if (lroot) write(*,*) 'Average settling velocity: ',i,vsetaver(i)
@@ -368 +362 @@
 
   ! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
-    if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
-         &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.)))  then
-      WETDEP=.true.
-      WETDEPSPEC(i)=.true.
+    if ( &
+      dquer(i).le.0. .and. (weta_gas(i)  .gt.0. .or. wetb_gas(i)  .gt.0.) .or. &
+      dquer(i).gt.0. .and. (crain_aero(i).gt.0. .or. csnow_aero(i).gt.0.)  &
+       )  then
+      wetdep        = .true.
+      wetdepspec(i) = .true.
       if (lroot) then
         write (*,*) '  Below-cloud scavenging: ON'
@@ -381 +377 @@
 
   ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
-    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
-      WETDEP=.true.
-      WETDEPSPEC(i)=.true.
+    if (dquer(i).gt.0. .and. (ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
+      wetdep        = .true.
+      wetdepspec(i) = .true.
       if (lroot) then
         write (*,*) '  In-cloud scavenging: ON'
@@ -394 +390 @@
 
     if (ohcconst(i).gt.0.) then
-      OHREA=.true.
+      ohrea = .true.
       if (lroot) write (*,*) '  OHreaction switched on: ',ohcconst(i),i
     endif
 
-    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
-      DRYDEP=.true.
-      DRYDEPSPEC(i)=.true.
+    if (reldiff(i).gt.0. .or. density(i).gt.0. .or. dryvel(i).gt.0.) then
+      drydep        = .true.
+      drydepspec(i) = .true.
     endif
 
@@ -415 +411 @@
 101 numpoint=numpoint+1
 
-  if (readerror.lt.1) then ! reading namelist format
+  if (ios.lt.1) then ! reading namelist format
 
     if (numpoint.gt.numpoints) goto 250
-    zkind = 1
-    mass = 0
-    parts = 0
+    ikindz = 1
+    rmass = -1. ! initialise as negative to catch if not given in input
+    nparts = 0
     comment = ' '
-    read(unitreleases,release,iostat=readerror)
+    read(unitreleases,release,iostat=ios)
     id1=idate1
     it1=itime1
     id2=idate2
     it2=itime2
-    xpoint1(numpoint)=lon1
-    xpoint2(numpoint)=lon2
-    ypoint1(numpoint)=lat1
-    ypoint2(numpoint)=lat2
+    xpoint1(numpoint)=rlon1
+    xpoint2(numpoint)=rlon2
+    ypoint1(numpoint)=rlat1
+    ypoint2(numpoint)=rlat2
     zpoint1(numpoint)=z1
     zpoint2(numpoint)=z2
-    kindz(numpoint)=zkind
+    kindz(numpoint)=ikindz
     do i=1,nspec
-      xmass(numpoint,i)=mass(i)
+      if (rmass(i) .lt. 0.) goto 995 ! write error and stop
+      xmass(numpoint,i)=rmass(i)
     end do
-    npart(numpoint)=parts
+    npart(numpoint)=nparts
     compoint(min(1001,numpoint))=comment
 
@@ -445 +442 @@
     endif
 
-  else
+  else ! simple text format
 
     read(unitreleases,*,end=250)
@@ -471 +468 @@
       read(unitreleases,*,err=998) xmass(numpoint,i)
       if (old) call skplin(2,unitreleases)
-      mass(i)=xmass(numpoint,i)
+      rmass(i)=xmass(numpoint,i)
     end do
   !save compoint only for the first 1000 release points
@@ -489 +486 @@
       idate2=id2
       itime2=it2
-      lon1=xpoint1(numpoint)
-      lon2=xpoint2(numpoint)
-      lat1=ypoint1(numpoint)
-      lat2=ypoint2(numpoint)
+      rlon1=xpoint1(numpoint)
+      rlon2=xpoint2(numpoint)
+      rlat1=ypoint1(numpoint)
+      rlat2=ypoint2(numpoint)
       z1=zpoint1(numpoint)
       z2=zpoint2(numpoint)
-      zkind=kindz(numpoint)
-      parts=npart(numpoint)
+      ikindz=kindz(numpoint)
+      nparts=npart(numpoint)
       write(unitreleasesout,nml=release)
     endif
 
     if (numpoint.le.1000) then
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
-           (compoint(numpoint)(1:8).eq.'        ')) goto 250
+      if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
+          xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0. .and. &
+          compoint(numpoint)(1:8).eq.'        ') goto 250
     else
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
+      if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
+          xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0.) goto 250
     endif
 
@@ -539 +536 @@
   !*********************************************************************
 
-  if (xpoint1(numpoint).lt.xlon0) &
-       xpoint1(numpoint)=xpoint1(numpoint)+360.
-  if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
-       xpoint1(numpoint)=xpoint1(numpoint)-360.
-  if (xpoint2(numpoint).lt.xlon0) &
-       xpoint2(numpoint)=xpoint2(numpoint)+360.
-  if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
-       xpoint2(numpoint)=xpoint2(numpoint)-360.
+  if (xpoint1(numpoint) .lt. xlon0) &
+      xpoint1(numpoint) = xpoint1(numpoint) + 360.
+  if (xpoint1(numpoint) .gt. xlon0+nxmin1*dx) &
+      xpoint1(numpoint) = xpoint1(numpoint) - 360.
+  if (xpoint2(numpoint) .lt. xlon0) &
+      xpoint2(numpoint) = xpoint2(numpoint) + 360.
+  if (xpoint2(numpoint) .gt. xlon0+nxmin1*dx) &
+      xpoint2(numpoint) = xpoint2(numpoint) - 360.
 
   ! Determine relative beginning and ending times of particle release
@@ -562 +559 @@
   if (mdomainfill.eq.0) then   ! no domain filling
     if (ldirect.eq.1) then
-      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
+      if (jul1.lt.bdate .or. jul2.gt.edate) then
         write(*,*) 'FLEXPART MODEL ERROR'
         write(*,*) 'Release starts before simulation begins or ends'
@@ -577 +574 @@
       endif
     else if (ldirect.eq.-1) then
-      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
+      if (jul1.lt.edate .or. jul2.gt.bdate) then
         write(*,*) 'FLEXPART MODEL ERROR'
         write(*,*) 'Release starts before simulation begins or ends'
@@ -586 +583 @@
       if (npart(numpoint).gt.num_min_discrete) then
         ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+        ireleaseend(numpoint) = int((jul2-bdate)*86400.)
       else
         ireleasestart(numpoint)=int((julm-bdate)*86400.)
-        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint) = int((julm-bdate)*86400.)
       endif
     endif
@@ -600 +597 @@
 
   if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
-    releaserate=releaserate+real(npart(numpoint))/ &
-         real(ireleaseend(numpoint)-ireleasestart(numpoint))
+    releaserate = releaserate + real( npart(numpoint) ) / &
+      real( ireleaseend(numpoint) - ireleasestart(numpoint) )
   else
     releaserate=99999999
@@ -634 +631 @@
       nsettle=0
       do ispc=1,nspec
-        if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
+        if (xmass(irel,ispc).gt.tiny(1.)) nsettle=nsettle+1
       end do
       if (nsettle.gt.1) lsettling=.false.
@@ -641 +638 @@
 
   if (lroot) then
-    if (.not.lsettling) then 
-      write(*,*) 'WARNING: more than 1 species per release point, settling &
-           &disabled'
-    end if
+    if (.not.lsettling) &
+     write(*,*) 'WARNING: more than 1 species per release point,' &
+           // ' settling disabled'
   end if
 
@@ -651 +647 @@
   !************************************************************************
 
-  if (releaserate.gt. &
-       0.99*real(maxpart)/real(lage(nageclass))) then
+  if (releaserate .gt. 0.99*real(maxpart)/real(lage(nageclass))) then
     if (numpartmax.gt.maxpart.and.lroot) then
       write(*,*) '#####################################################'
@@ -670 +665 @@
       write(*,*) 'Maximum allowed release rate is: ', &
            real(maxpart)/real(lage(nageclass)),' particles per second'
-      write(*,*) &
-           'Total number of particles released during the simulation is: ', &
-           numpartmax
+      write(*,*) 'Total number of particles released during the simulation is:'&
+        ,numpartmax
       write(*,*) 'Maximum allowed number of particles is: ',maxpart
     endif
@@ -679 +673 @@
 
   if (lroot) then
-    write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
+    write(*,FMT='(A,ES14.7)') ' Total mass released:', &
+      sum(xmass(1:numpoint,1:nspec))
   end if
 
   return
 
-994 write(*,*) '#####################################################'
+995 write(*,*) '#####################################################'
   write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
   write(*,*) '####                                             ####'
-  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
-  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
+  write(*,*) '#### ERROR - mass for species',i,' not given     ####'
   write(*,*) '#####################################################'
+  stop 'readreleases - undefined mass'
+
+994 write(*,*) '######################################################'
+  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:      ####'
+  write(*,*) '####                                              ####'
+  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS ####'
+  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES.  ####'
+  write(*,*) '######################################################'
   stop
 
@@ -702 +704 @@
   stop
 
-
 999 write(*,*) '#####################################################'
   write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
@@ -713 +714 @@
 
 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,*)    ' #### CANNOT BE OPENED IN THE DIRECTORY        #### '
   write(*,'(a)') path(2)(1:length(2))
   stop

src/timemanager.f90

@@ -719 +719 @@
   ! Output to keep track of the numerical instabilities in CBL simulation and if
   ! they are compromising the final result (or not)
-    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl  
+    if (iflagcbl.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl  
           
   end do

src/timemanager_mpi.f90

@@ -876 +876 @@
 ! Output to keep track of the numerical instabilities in CBL simulation
 ! and if they are compromising the final result (or not):
-    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl  
+    if (iflagcbl.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl