Maintenance is scheduled between 16:00 and 23:59 CEST (14:00 and 21:59 UTC) on Thursday 2021-10-28. The system may be unavailable at any time during this timeframe. Please plan accordingly.

Commit d7935deb authored by pesei's avatar pesei
Browse files

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
parent 34f14525
&RELEASES_CTRL
NSPEC= 1,
SPECNUM_REL= 24,
NUMSPEC_REL= 24,
/
&RELEASE
IDATE1= 20170102,
ITIME1= 090000,
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",
/
......@@ -420,7 +420,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
if (dtftlw.lt..5) then
!*************************************************************
!************** 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
!if (ol.gt.0.) then !modified by mc : for test
......@@ -514,7 +514,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
endif
end do
if (cblflag.ne.1) nrand=nrand+i
if (iflagcbl.ne.1) nrand=nrand+i
! Determine time step for next integration
!*****************************************
......
......@@ -25,7 +25,7 @@ module com_mod
character :: path(numpath+2*maxnests)*120
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
! length length of path names needed for trajectory model
......@@ -33,7 +33,7 @@ module com_mod
! flexversion version of flexpart (descriptive long string)
! 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
!********************************************************
! Variables defining the general model run specifications
......@@ -72,7 +72,7 @@ module com_mod
integer :: mquasilag,nested_output,ind_source,ind_receptor
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
! ifine reduction factor for time step used for vertical wind
......@@ -103,7 +103,7 @@ module com_mod
! turbswitch determines how the Markov chain is formulated
! 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
integer :: mintime,itsplit
......
......@@ -157,7 +157,7 @@ subroutine initialize(itime,ldt,up,vp,wp, &
wp=rannumb(nrand+2)
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
!if (ol.gt.0.) then !by mc : only for test correct is lt.0
......
......@@ -59,12 +59,12 @@ subroutine readOHfield
real, parameter :: scalehgt=7000. ! scale height in metres
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
......
......@@ -28,9 +28,13 @@ subroutine readageclasses
! *
! 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 *
! *
!*****************************************************************************
! *
! Variables: *
......@@ -46,11 +50,11 @@ subroutine readageclasses
integer :: i
! namelist help variables
integer :: readerror
! namelist aux variables
integer :: ios
! namelist declaration
namelist /ageclass/ &
namelist /nml_ageclass/ &
nageclass, &
lage
......@@ -70,41 +74,35 @@ subroutine readageclasses
! open the AGECLASSSES file and read user options
!************************************************
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
if (lage(1).le.0) then
write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST #### '
write(*,*) ' #### CLASS MUST BE GREATER THAN ZERO. CHANGE #### '
......@@ -125,7 +123,7 @@ subroutine readageclasses
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
1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
......@@ -133,5 +131,18 @@ subroutine readageclasses
write(*,'(a)') path(2)(1:length(2))
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
......@@ -32,6 +32,7 @@ subroutine readcommand
! HSO, 1 July 2014: Added optional namelist input *
! Unknown, unknown: various *
! Petra Seibert, 2018-06-08: improve error msgs *
! PS 6/2015: Minor changes in variable names and layout *
! *
!*****************************************************************************
! *
......@@ -83,97 +84,97 @@ subroutine readcommand
real(kind=dp) :: juldate
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
!**************************************************************************
open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
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.
if (lroot) write(*,*) 'COMMAND in old short format, &
......@@ -242,7 +243,7 @@ subroutine readcommand
read(unitcommand,*) surf_only
! 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)
......@@ -251,7 +252,7 @@ subroutine readcommand
! write command file in namelist format to output directory if requested
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
......@@ -259,22 +260,22 @@ subroutine readcommand
! 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.
else
......@@ -285,27 +286,27 @@ subroutine readcommand
fine=1./real(ifine)
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
!Af set release-switch
......@@ -327,6 +328,8 @@ subroutine readcommand
else ! mass mix
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
ind_rel = 1
......@@ -336,10 +339,10 @@ subroutine readcommand
ind_rel = 3
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
ind_rel = 4
......@@ -348,7 +351,7 @@ subroutine readcommand
write(*,*) ' #### Releaseheight is forced to 0 - 2*href #### '
write(*,*) ' #### Release is performed above ground lev #### '
end if
DRYBKDEP=.true.
drybkdep=.true.
allocate(xscav_frac1(maxpart,maxspec))
end select
endif
......@@ -398,14 +401,15 @@ subroutine readcommand
mintime=lsynctime
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
endif
......@@ -413,7 +417,7 @@ subroutine readcommand
! Check whether a valid option for gridded model output has been chosen
!**********************************************************************
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 #### '
write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### '
......@@ -422,8 +426,8 @@ subroutine readcommand
endif
!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 #### '
write(*,*) ' #### FOR MASS UNITS (at the moment) #### '
......@@ -478,8 +482,8 @@ subroutine readcommand
! 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"
!*****************************************************************************
if (ldirect.lt.0) then
......@@ -648,8 +652,11 @@ subroutine readcommand
write(*,900) ' #### CANNOT WRITE TO '// &
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
......@@ -28,8 +28,9 @@ subroutine readoutgrid
! Author: A. Stohl *
! *
! 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 *
! *
!*****************************************************************************
! *
......@@ -50,64 +51,66 @@ subroutine readoutgrid
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=