*
L. Bakels 2022: This module contains all subroutines for * reading option files * *
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later !***************************************************************************** ! * ! L. Bakels 2022: This module contains all subroutines for * ! reading option files * ! * !***************************************************************************** module readoptions_mod use par_mod use com_mod use date_mod use point_mod use windfields_mod implicit none contains subroutine readageclasses !***************************************************************************** ! * ! This routine reads the age classes to be used for the current model * ! run. * ! * ! Author: A. Stohl * ! 20 March 2000 * ! HSO, 1 July 2014 * ! Added optional namelist input * ! * !***************************************************************************** ! * ! Variables: * ! * ! Constants: * ! * !***************************************************************************** implicit none integer :: i ! namelist help variables integer :: ios,stat ! namelist declaration namelist /nage/ & nageclass namelist /ageclass/ & lage nageclass=-1 ! preset to negative value to identify failed namelist input ! If age spectra calculation is switched off, set number of age classes ! to 1 and maximum age to a large number !********************************************************************** if (lagespectra.ne.1) then nageclass=1 allocate( lage(nageclass),stat=stat) if (stat.ne.0) error stop "Could not allocate lage" lage(nageclass)=999999999 return endif ! If age spectra claculation is switched on, ! open the AGECLASSSES file and read user options !************************************************ open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', & form='formatted',status='old',err=999) ! try to read in as a namelist read(unitageclasses,nage,iostat=ios) allocate( lage(nageclass),stat=stat) if (stat.ne.0) error stop "Could not allocate lage" read(unitageclasses,ageclass,iostat=ios) close(unitageclasses) if ((nageclass.lt.0).or.(ios.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 allocate( lage(nageclass),stat=stat) if (stat.ne.0) error stop "Could not allocate lage" read(unitageclasses,*) lage(1) do i=2,nageclass read(unitageclasses,*) lage(i) end do close(unitageclasses) endif ! 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=nage) write(unitageclasses,nml=ageclass) close(unitageclasses) endif if (lage(1).le.0) then write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST #### ' write(*,*) ' #### CLASS MUST BE GREATER THAN ZERO. CHANGE #### ' write(*,*) ' #### SETTINGS IN FILE AGECLASSES. #### ' error stop 'First age class must be larger than zero' endif do i=2,nageclass if (lage(i).le.lage(i-1)) then write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES #### ' write(*,*) ' #### MUST BE GIVEN IN TEMPORAL ORDER. #### ' write(*,*) ' #### CHANGE SETTINGS IN FILE AGECLASSES. #### ' error stop 'Age classes must be in temporal order' endif end do return 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(1)(1:length(1)) error stop 'AGECLASSES cannot be opened' 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop 'AGECLASSES cannot be opened' end subroutine readageclasses subroutine readavailable !***************************************************************************** ! * ! This routine reads the dates and times for which windfields are * ! available. * ! * ! Authors: A. Stohl * ! * ! 6 February 1994 * ! 8 February 1999, Use of nested fields, A. Stohl * ! * !***************************************************************************** ! * ! Variables: * ! bdate beginning date as Julian date * ! beg beginning date for windfields * ! endl ending date for windfields * ! fname filename of wind field, help variable * ! ideltas [s] duration of modelling period * ! idiff time difference between 2 wind fields * ! idiffnorm normal time difference between 2 wind fields * ! idiffmax [s] maximum allowable time between 2 wind fields * ! jul julian date, help variable * ! numbwf actual number of wind fields * ! wfname(numbwf) file names of needed wind fields * ! wftime(numbwf) [s]times of wind fields relative to beginning time * ! wfname1,wftime1 = same as above, but only local (help variables) * ! * ! Constants: * ! unitavailab unit connected to file AVAILABLE * ! * !***************************************************************************** implicit none integer :: i,idiff,ldat,ltim,k,stat integer,allocatable,dimension(:) :: wftime1,tmpwftime,numbwfn integer,allocatable,dimension(:,:) :: wftimen,wftime1n,tmpwftimen logical :: lwarntd=.true. real(kind=dp) :: jul,beg,endl character(len=255) :: fname character(len=255),allocatable,dimension(:) :: wfname1, & tmpwfname character(len=255),allocatable,dimension(:,:) :: wfname1n,tmpwfnamen ! Windfields are only used, if they are within the modelling period. ! However, 1 additional day at the beginning and at the end is used for ! interpolation. -> Compute beginning and ending date for the windfields. !************************************************************************ if (ideltas.gt.0) then ! forward trajectories beg=bdate-1._dp endl=bdate+real(ideltas,kind=dp)/86400._dp+real(idiffmax,kind=dp)/ & 86400._dp else ! backward trajectories beg=bdate+real(ideltas,kind=dp)/86400._dp-real(idiffmax,kind=dp)/ & 86400._dp endl=bdate+1._dp endif ! Open the wind field availability file and read available wind fields ! within the modelling period. !********************************************************************* open(unitavailab,file=path(4)(1:length(4)),status='old', & err=999) do i=1,3 read(unitavailab,*) end do numbwf=0 100 read(unitavailab,'(i8,1x,i6,2(6x,a255))',end=99) & ldat,ltim,fname jul=juldate(ldat,ltim) if ((jul.ge.beg).and.(jul.le.endl)) then numbwf=numbwf+1 allocate( tmpwfname(numbwf),tmpwftime(numbwf),stat=stat) if (stat.ne.0) error stop 'ERROR: could not allocate tmpwfname' if (numbwf.gt.1) then tmpwfname(1:numbwf-1)=wfname1 tmpwftime(1:numbwf-1)=wftime1 endif tmpwfname(numbwf)=fname(1:index(fname,' ')) tmpwftime(numbwf)=nint((jul-bdate)*86400._dp) call move_alloc(tmpwfname,wfname1) call move_alloc(tmpwftime,wftime1) endif goto 100 ! next wind field 99 continue close(unitavailab) ! Open the wind field availability file and read available wind fields ! within the modelling period (nested grids) !********************************************************************* if (numbnests.gt.0) then allocate( numbwfn(numbnests),stat=stat) if (stat.ne.0) error stop 'ERROR: could not allocate numwfn' do k=1,numbnests !print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3) !print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2)) open(unitavailab,file=path(numpath+2*(k-1)+2) & (1:length(numpath+2*(k-1)+2)),status='old',err=998) do i=1,3 read(unitavailab,*) end do numbwfn(k)=0 700 read(unitavailab,'(i8,1x,i6,2(6x,a255))',end=699) ldat, & ltim,fname jul=juldate(ldat,ltim) if ((jul.ge.beg).and.(jul.le.endl)) then numbwfn(k)=numbwfn(k)+1 allocate( tmpwfnamen(numbnests,numbwfn(k)),tmpwftimen(numbnests,numbwfn(k)), & stat=stat) if (stat.ne.0) error stop 'ERROR: could not allocate tmpwfnamen' if (numbwfn(k).gt.1) then tmpwfnamen(:,1:numbwfn(k)-1)=wfname1n tmpwftimen(:,1:numbwfn(k)-1)=wftime1n endif tmpwfnamen(k,numbwfn(k))=fname(1:index(fname,' ')) tmpwftimen(k,numbwfn(k))=nint((jul-bdate)*86400._dp) call move_alloc(tmpwfnamen,wfname1n) call move_alloc(tmpwftimen,wftime1n) endif goto 700 ! next wind field 699 continue close(unitavailab) end do endif ! Check wind field times of file AVAILABLE (expected to be in temporal order) !**************************************************************************** if (numbwf.eq.0) then write(*,*) ' #### FLEXPART MODEL ERROR! NO WIND FIELDS #### ' write(*,*) ' #### AVAILABLE FOR SELECTED TIME PERIOD. #### ' error stop 'No wind fields available for selected time period' endif do i=2,numbwf if (wftime1(i).le.wftime1(i-1)) then write(*,*) 'FLEXPART ERROR: FILE AVAILABLE IS CORRUPT.' write(*,*) 'THE WIND FIELDS ARE NOT IN TEMPORAL ORDER.' write(*,*) 'PLEASE CHECK FIELD ',wfname1(i) error stop 'AVAILABLE file error' endif end do ! Check wind field times of file AVAILABLE for the nested fields ! (expected to be in temporal order) !*************************************************************** do k=1,numbnests if (numbwfn(k).eq.0) then write(*,*) '#### FLEXPART MODEL ERROR! NO WIND FIELDS ####' write(*,*) '#### AVAILABLE FOR SELECTED TIME PERIOD. ####' error stop 'No nested wind fields available for selected time period' endif do i=2,numbwfn(k) if (wftime1n(k,i).le.wftime1n(k,i-1)) then write(*,*) 'FLEXPART ERROR: FILE AVAILABLE IS CORRUPT. ' write(*,*) 'THE NESTED WIND FIELDS ARE NOT IN TEMPORAL ORDER.' write(*,*) 'PLEASE CHECK FIELD ',wfname1n(k,i) write(*,*) 'AT NESTING LEVEL ',k error stop 'nested AVAILABLE file error' endif end do end do ! Allocating global fields storing the windfield names and times !*************************************************************** allocate(wfname(numbwf),stat=stat) if (stat.ne.0) error stop "Could not allocate wfname" allocate(wftime(numbwf),stat=stat) if (stat.ne.0) error stop "Could not allocate wftime" if (numbnests.gt.0) then allocate(wfnamen(numbnests,numbwf),stat=stat) if (stat.ne.0) error stop "Could not allocate wfnamen" allocate(wftimen(numbnests,numbwf),stat=stat) if (stat.ne.0) error stop "Could not allocate wftimen" endif ! For backward trajectories, reverse the order of the windfields !*************************************************************** if (ideltas.ge.0) then do i=1,numbwf wfname(i)=wfname1(i) wftime(i)=wftime1(i) end do do k=1,numbnests do i=1,numbwfn(k) wfnamen(k,i)=wfname1n(k,i) wftimen(k,i)=wftime1n(k,i) end do end do else do i=1,numbwf wfname(numbwf-i+1)=wfname1(i) wftime(numbwf-i+1)=wftime1(i) end do do k=1,numbnests do i=1,numbwfn(k) wfnamen(k,numbwfn(k)-i+1)=wfname1n(k,i) wftimen(k,numbwfn(k)-i+1)=wftime1n(k,i) end do end do endif ! Check the time difference between the wind fields. If it is big, ! write a warning message. If it is too big, terminate the trajectory. !********************************************************************* do i=2,numbwf idiff=abs(wftime(i)-wftime(i-1)) if (idiff.gt.idiffmax.and.lroot) then write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO' write(*,*) 'WIND FIELDS IS TOO BIG FOR TRANSPORT CALCULATION.& &' write(*,*) 'THEREFORE, TRAJECTORIES HAVE TO BE SKIPPED.' else if (idiff.gt.idiffnorm.and.lroot.and.lwarntd) then write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO' write(*,*) 'WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION' write(*,*) 'OF SIMULATION QUALITY.' lwarntd=.false. ! only issue this warning once endif end do do k=1,numbnests if (numbwfn(k).ne.numbwf) then write(*,*) 'FLEXPART ERROR: THE AVAILABLE FILES FOR THE' write(*,*) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH' write(*,*) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. ' write(*,*) 'ERROR AT NEST LEVEL: ',k error stop 'Not the same number of nested wind field files as mother domain' endif do i=1,numbwf if (wftimen(k,i).ne.wftime(i)) then write(*,*) 'FLEXPART ERROR: THE AVAILABLE FILES FOR THE' write(*,*) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH' write(*,*) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. ' write(*,*) 'ERROR AT NEST LEVEL: ',k error stop 'Not the same time period of nested wind field files as mother domain' endif end do end do ! Reset the times of the wind fields that are kept in memory to no time !********************************************************************** do i=1,2 memind(i)=i memtime(i)=999999999 end do deallocate( wftime1,wfname1) if (numbnests.gt.0) deallocate(wftime1n,wftimen,wfname1n,numbwfn) return 998 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### ' write(*,'(a)') ' '//path(numpath+2*(k-1)+2) & (1:length(numpath+2*(k-1)+2)) write(*,*) ' #### CANNOT BE OPENED #### ' error stop 'Nested AVAILABLE file cannot be opened' 999 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### ' write(*,'(a)') ' '//path(4)(1:length(4)) write(*,*) ' #### CANNOT BE OPENED #### ' error stop 'AVAILABLE file cannot be opened' end subroutine readavailable subroutine readcommand !***************************************************************************** ! * ! This routine reads the user specifications for the current model run. * ! * ! Author: A. Stohl * ! * ! 18 May 1996 * ! HSO, 1 July 2014 * ! Added optional namelist input * ! * ! * ! January 2024 Rona Thompson * ! Added new variables for LCM * ! * !***************************************************************************** ! * ! Variables: * ! bdate beginning date as Julian date * ! ctl factor by which time step must be smaller than * ! Lagrangian time scale * ! ibdate,ibtime beginnning date and time (YYYYMMDD, HHMISS) * ! ideltas [s] modelling period * ! iedate,ietime ending date and time (YYYYMMDD, HHMISS) * ! ifine reduction factor for vertical wind time step * ! outputforeachrel for forward runs it is possible either to create * ! one outputfield or several for each releasepoint * ! iflux switch to turn on (1)/off (0) flux calculations * ! iout 1 for conc. (residence time for backward runs) output,* ! 2 for mixing ratio output, 3 both, 4 for plume * ! trajectory output, 5 = options 1 and 4 * ! ipin 1 continue simulation with restart.bin file, * ! 2 continue simulaion with dumped particle data, 0 no * ! 3 use self-defined initial conditions in netcdf * ! 4 initial run using option 3, restart from restart.bin* ! ipout 0 no particle dump, 1 every output time, 3 only at end* ! ipoutfac increase particle dump interval by factor (default 1) * ! loutaver [s] concentration output is an average over loutaver * ! seconds * ! loutsample [s] average is computed from samples taken every [s] * ! seconds * ! loutstep [s] time interval of concentration output * ! lrecoutstep [s] time interval of receptor output * ! lrecoutaver [s] receptor output is an average of lrecoutaver seconds * ! lrecoutsample [s] average is computed from samples taken every [s] * ! lsynctime [s] synchronisation time interval for all particles * ! lagespectra switch to turn on (1)/off (0) calculation of age * ! spectra * ! lconvection value of either 0 and 1 indicating mixing by * ! convection * ! = 0 .. no convection * ! + 1 .. parameterisation of mixing by subgrid-scale * ! convection = on * ! lsubgrid switch to turn on (1)/off (0) subgrid topography * ! parameterization * ! method method used to compute the particle pseudovelocities * ! mdomainfill 1 use domain-filling option, 0 not, 2 use strat. O3 * ! * ! Constants: * ! unitcommand unit connected to file COMMAND * ! * !***************************************************************************** implicit none character(len=50) :: line integer :: ios integer :: lturbulence_meso,lcmoutput character(len=50) :: ohfields_path ! deprecated namelist /command/ & ldirect, & ibdate,ibtime, & iedate,ietime, & loutstep, & loutaver, & loutsample, & loutrestart, & lrecoutstep, & lrecoutaver, & lrecoutsample, & lsynctime, & ctl, & ifine, & iout, & ipout, & ipoutfac, & lsubgrid, & lconvection, & lturbulence, & lturbulence_meso, & lagespectra, & ipin, & ioutputforeachrelease, & iflux, & mdomainfill, & ind_source, & ind_receptor, & mquasilag, & nested_output, & linit_cond, & lnetcdfout, & sfc_only, & surf_only, & cblflag, & linversionout, & d_trop, & d_strat, & nxshift, & maxthreadgrid, & maxfilesize, & logvertinterp, & ohfields_path, & lcmoutput, & itsplit ! deprecated: only for IO back compatibility ! Presetting namelist command ldirect=0 ibdate=20000101 ibtime=0 iedate=20000102 ietime=0 loutstep=10800 loutaver=10800 loutsample=900 lrecoutstep=-1 lrecoutaver=-1 lrecoutsample=-1 loutrestart=-1 lsynctime=900 ctl=-5.0 ifine=4 iout=3 ipout=0 ipoutfac=1 lsubgrid=1 lconvection=1 lturbulence=1 lturbulence_meso=0 lagespectra=0 ipin=0 ioutputforeachrelease=1 iflux=1 mdomainfill=0 ind_source=1 ind_receptor=1 mquasilag=0 nested_output=0 linit_cond=0 lnetcdfout=1 sfc_only=0 surf_only=-1 cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine linversionout=0 nxshift=-9999 maxthreadgrid=1 maxfilesize=10000 logvertinterp=0 ohfields_path='' lcmoutput=0 itsplit=999999999 ! deprecated: only for IO back compatibility !Af set release-switch 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=ios) if (ios.ne.0) then backspace(unitcommand) read(unitcommand,fmt='(A)') line if (lroot) write(*,*) & 'Invalid line in COMMAND reads: '//trim(line) end if close(unitcommand) if (lturbulence_meso.ne.0) lmesoscale_turb=.true. if (logvertinterp.ne.0) log_interpol=.true. if (surf_only.ne.-1) then write(*,*) 'WARNING: SURF_ONLY in COMMAND will be deprecated and renamed SFC_ONLY' write(*,*) 'Continuing with SURF_ONLY...' sfc_only=surf_only endif ! distinguish namelist from fixed text input if ((ios.ne.0).or.(ldirect.eq.0)) then ! parse as text file format if (lroot) write(*,*) 'COMMAND either having unrecognised entries, & &or in old format, please update to namelist format.' error stop 'COMMAND has unrecognised entries or is not a namelist' endif ! input format ! 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=1000) write(unitcommand,nml=command) close(unitcommand) endif ifine=max(ifine,1) ! 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 turbswitch=.true. 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' ctl=5. end if 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 endif print *,'WARNING: CBL flag active the ratio of TLW/dt is ',ctl*ifine print *,'WARNING: CBL flag active lsynctime is ',lsynctime else !added by mc if (ctl.ge.0.1) then turbswitch=.true. else turbswitch=.false. ifine=1 endif endif !added by mc 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 particles 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 ! 0 = no receptors !Af 1 = mass units !Af 2 = mass mixing ratio units ! 3 = wet deposition in outputfield ! 4 = dry deposition in outputfield ! Settings for LCM output !************************************************************ ! MDOMAINFILL = 1 | LLCMOUTPUT = true ! IND_SOURCE = 1 | IND_SAMP = 0 ! IND_RECEPTOR = 1 | calculates mass ratio mixing ratio ! IOUT = 2 | as ratio species_mass to airtracer_mass !------------------------------------------------------------ if ( ldirect .eq. 1 ) then ! FWD-Run !Af set release-switch if (ind_source .eq. 1 ) then !mass ind_rel = 0 else ! mass mix ind_rel = 1 endif !Af set sampling switch if (ind_receptor .le. 1) then !mass ind_samp = 0 else ! mass mix ind_samp = -1 endif elseif (ldirect .eq. -1 ) then !BWD-Run !Af set sampling switch if (ind_source .eq. 1 ) then !mass ind_samp = -1 else ! mass mix ind_samp = 0 endif select case (ind_receptor) case (1) ! 1 .. concentration at receptor ind_rel = 1 case (2) ! 2 .. mixing ratio at receptor ind_rel = 0 case (3) ! 3 .. wet deposition in outputfield ind_rel = 3 if (lroot) then write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE #### ' write(*,*) ' #### Releaseheight is forced to 0 - 20km #### ' write(*,*) ' #### Release is performed above ground lev #### ' end if WETBKDEP=.true. !allocate(xscav_frac1(maxpart,maxspec)) case (4) ! 4 .. dry deposition in outputfield ind_rel = 4 if (lroot) then write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### ' write(*,*) ' #### Releaseheight is forced to 0 - 2*href #### ' write(*,*) ' #### Release is performed above ground lev #### ' end if DRYBKDEP=.true. !allocate(xscav_frac1(maxpart,maxspec)) end select endif !************************************************************* ! Check whether valid options have been chosen in file COMMAND !************************************************************* ! Check options for initial condition output: Switch off for forward runs !************************************************************************ if (ldirect.eq.1) linit_cond=0 if ((linit_cond.lt.0).or.(linit_cond.gt.2)) then write(*,*) ' #### FLEXPART MODEL ERROR! INVALID OPTION #### ' write(*,*) ' #### FOR LINIT_COND IN FILE "COMMAND". #### ' error stop 'LINIT_COND in COMMAND is invalid' endif ! Check input dates !****************** if (iedate.lt.ibdate) then write(*,*) ' #### FLEXPART MODEL ERROR! BEGINNING DATE #### ' write(*,*) ' #### IS LARGER THAN ENDING DATE. CHANGE #### ' write(*,*) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### ' write(*,*) ' #### "COMMAND". #### ' error stop 'COMMAND: beginning date is larger than ending date' else if (iedate.eq.ibdate) then if (ietime.lt.ibtime) then write(*,*) ' #### FLEXPART MODEL ERROR! BEGINNING TIME #### ' write(*,*) ' #### IS LARGER THAN ENDING TIME. CHANGE #### ' write(*,*) ' #### EITHER POINT 2 OR POINT 3 IN FILE #### ' write(*,*) ' #### "COMMAND". #### ' error stop 'COMMAND: beginning time is larger than ending time' endif endif ! #ifndef USE_NCF ! if ((loutrestart.ne.-1).or.(ipin.ne.0)) then ! write(*,*) ' WARNING: restart option set with intervals' ! write(*,*) ' LOUTRESTART', loutrestart ! write(*,*) ' not possible when using binary gridded output' ! write(*,*) ' ==> RESTART FUNCTION SWITCHED OFF!' ! endif ! if (ipin.ne.0) then ! write(*,*) ' ERROR: restart option not possible using binary' ! write(*,*) ' output.' ! write(*,*) ' Please only use IPIN>0 when compiling and running using' ! write(*,*) ' netcdf output. ' ! endif ! #else ! if ((surf_only.eq.1).or.(linversionout.eq.1)) then ! write(*,*) ' ERROR: NetCDF output for surface only or for inversions' ! write(*,*) ' is not yet implemented. Please compile without NetCDF.' ! error stop 'Surface only option is not supported for NetCDF' ! endif ! #endif ! Determine kind of dispersion method !************************************ if (ctl.gt.0.) then method=1 mintime=minstep else method=0 mintime=lsynctime endif ! Check for netcdf output switch !******************************* if (iout.gt.8) then lnetcdfout=1 iout = iout -8 endif if (lnetcdfout.eq.1) then #ifndef USE_NCF write(*,*) 'WARNING: netcdf output not activated during compile time & &but switched on in COMMAND or set to default value 1.' write(*,*) 'Please recompile with netcdf library (`make [...] ncf=yes`) & &when requiring NetCDF output.' write(*,*) 'LNETCDFOUT set to 0.' lnetcdfout = 0 #endif else #ifdef USE_NCF write(*,*) 'WARNING: Executable compiled using NetCDF libraries, but & &BINARY output is requested. If this was unintended, please add 8 & &to IOUT or set LOUTNETCDF=1 in the COMMAND file.' #endif endif if ((lnetcdfout.eq.1).and.((sfc_only.eq.1).or.(linversionout.eq.1))) then write(*,*) ' WARNING: NetCDF output for surface only or for inversions' write(*,*) ' is not yet implemented. LNETCDFOUT set to 0.' lnetcdfout=0 endif #ifndef USE_NCF if (ipout.ne.0) then write(*,*) 'ERROR: NETCDF missing! Please recompile with the netcdf' write(*,*) 'library if you want the particle dump or set IPOUT=0.' error stop 'FLEXPART not compiled with NetCDF' endif #endif ! Check whether RECEPTOR commands are given, otherwise give them default values !****************************************************************************** if (lrecoutstep.eq.-1) then write(*,*) 'WARNING: FILE COMMAND LRECOUTSTEP not provided,' write(*,*) 'value of LOUTSTEP will be used if RECEPTORS are' write(*,*) 'required.' lrecoutstep=loutstep endif if (lrecoutaver.eq.-1) then write(*,*) 'WARNING: FILE COMMAND LRECOUTAVER not provided,' write(*,*) 'value of LOUTAVER will be used if RECEPTORS are' write(*,*) 'required.' lrecoutaver=loutaver endif if (lrecoutsample.eq.-1) then write(*,*) 'WARNING: FILE COMMAND LRECOUTSTEP not provided,' write(*,*) 'value of LOUTSAMPLE will be used if RECEPTORS are' write(*,*) 'required.' lrecoutsample=loutsample endif ! Check whether a valid option for gridded model output has been chosen !********************************************************************** if (iout.eq.0) then write(*,*) 'WARNING: IOUT set to zero, no gridded information will be & &written to file' else if ((iout.lt.0).or.(iout.gt.5)) 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 #### ' write(*,*) ' #### FOR NETCDF OUTPUT #### ' error stop 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 write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' write(*,*) ' #### VOLUME MIXING RATIO ONLY SUPPORTED #### ' write(*,*) ' #### FOR MASS UNITS (at the moment) #### ' error stop endif ! For quasilag output for each release is forbidden !***************************************************************************** if ((ioutputforeachrelease.eq.1).and.(mquasilag.eq.1)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### OUTPUTFOREACHRELEASE AND QUASILAGRANGIAN####' write(*,*) '#### MODE IS NOT POSSIBLE ! ####' error stop endif ! For quasilag backward is forbidden !***************************************************************************** if ((ldirect.lt.0).and.(mquasilag.eq.1)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### FOR BACKWARD RUNS, QUASILAGRANGIAN MODE ####' write(*,*) '#### IS NOT POSSIBLE ! ####' error stop endif ! For backward runs one releasefield for all releases makes no sense, ! For quasilag and domainfill ioutputforechrelease is forbidden !***************************************************************************** if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####' write(*,*) '#### MUST BE SET TO ONE! ####' error stop endif ! For backward runs one releasefield for all releases makes no sense, ! and is "forbidden" !***************************************************************************** if ((mdomainfill.eq.1).and.(ioutputforeachrelease.eq.1)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR ####' write(*,*) '#### EACH RELEASE IS FORBIDDEN ! ####' error stop endif ! Inversion output format only for backward runs !***************************************************************************** if ((linversionout.eq.1).and.(ldirect.eq.1)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR ####' write(*,*) '#### BACKWARD RUNS ####' error stop endif ! 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" !***************************************************************************** if (ldirect.lt.0) then if ((iout.eq.2).or.(iout.eq.3)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### FOR BACKWARD RUNS, IOUT MUST BE 1,4,OR 5####' error stop endif endif ! For domain-filling trajectories, a plume centroid trajectory makes no sense, ! and is "forbidden" !***************************************************************************** if (mdomainfill.ge.1) then if ((iout.eq.4).or.(iout.eq.5)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####' write(*,*) '#### FOR DOMAIN-FILLING TRAJECTORY OPTION, ####' write(*,*) '#### IOUT MUST NOT BE SET TO 4 OR 5. ####' error stop endif endif ! Check whether a valid options for particle dump has been chosen !**************************************************************** if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3! #### ' error stop endif ! Check whether input and output settings don't contradict !********************************************************* ! if (((iout.eq.4).or.(iout.eq.5)).and.((ipin.eq.3).or.(ipin.eq.4))) then ! write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' ! write(*,*) ' #### IOUT CANNOT BE 4 or 5 (plume) WHEN #### ' ! write(*,*) ' #### READING FROM part_ic.nc (ipin=4/5) #### ' ! error stop ! endif if(lsubgrid.ne.1.and.verbosity.eq.0) then write(*,*) ' ---------------- ' write(*,*) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS' write(*,*) ' NOT PARAMETERIZED DURING THIS SIMULATION. ' write(*,*) ' ---------------- ' endif ! Check whether convection scheme is either turned on or off !*********************************************************** if ((lconvection.ne.0).and.(lconvection.ne.1)) then write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### ' write(*,*) ' #### LCONVECTION MUST BE SET TO EITHER 1 OR 0#### ' error stop endif ! Check whether synchronisation interval is sufficiently short !************************************************************* if (lsynctime.gt.(idiffnorm/2)) then write(*,*) ' #### FLEXPART MODEL ERROR! SYNCHRONISATION #### ' write(*,*) ' #### TIME IS TOO LONG. MAKE IT SHORTER. #### ' write(*,*) ' #### MINIMUM HAS TO BE: ', idiffnorm/2 error stop endif ! Check consistency of the intervals, sampling periods, etc., for model output !***************************************************************************** if (loutaver.eq.0) then write(*,*) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### ' write(*,*) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### ' write(*,*) ' #### ZERO. #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' error stop endif if (loutaver.gt.loutstep) then write(*,*) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### ' write(*,*) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### ' write(*,*) ' #### GREATER THAN INTERVAL OF OUTPUT. #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' error stop endif if (loutsample.gt.loutaver) then write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### ' write(*,*) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### ' write(*,*) ' #### GREATER THAN TIME AVERAGE OF OUTPUT. #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' error stop endif if (mod(loutaver,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### ' write(*,*) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' error stop endif if ((loutaver/lsynctime).lt.2) then write(*,*) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### ' write(*,*) ' #### CONCENTRATION FIELD MUST BE AT LEAST #### ' write(*,*) ' #### TWICE THE SYNCHRONISATION INTERVAL #### ' error stop endif if (mod(loutstep,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### ' write(*,*) ' #### CONCENTRATION FIELDS MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' error stop endif if ((loutstep/lsynctime).lt.2) then write(*,*) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### ' write(*,*) ' #### CONCENTRATION FIELDS MUST BE AT LEAST #### ' write(*,*) ' #### TWICE THE SYNCHRONISATION INTERVAL #### ' error stop endif if (mod(loutsample,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### ' write(*,*) ' #### CONCENTRATION FIELD MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' error stop endif if ((mquasilag.eq.1).and.(iout.ge.4)) then write(*,*) ' #### FLEXPART MODEL ERROR! CONFLICTING #### ' write(*,*) ' #### OPTIONS: IF MQUASILAG=1, PLUME #### ' write(*,*) ' #### TRAJECTORY OUTPUT IS IMPOSSIBLE. #### ' error stop endif ! Check consistency of the intervals for receptors ! ************************************************ ! only if ldirect=1 if (ldirect.eq.1) then if (lrecoutaver.eq.0) then write(*,*) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST NOT BE ZERO #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' stop endif if (lrecoutaver.gt.lrecoutstep) then write(*,*) ' #### FLEXPART MODEL ERROR! TIME AVERAGE OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST NOT BE #### ' write(*,*) ' #### GREATER THAN INTERVAL OF OUTPUT. #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' stop endif if (lrecoutsample.gt.lrecoutaver) then write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST NOT BE #### ' write(*,*) ' #### GREATER THAN TIME AVERAGE OF OUTPUT. #### ' write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### ' stop endif if (mod(lrecoutaver,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' stop endif if ((lrecoutaver/lsynctime).lt.2) then write(*,*) ' #### FLEXPART MODEL ERROR! AVERAGING TIME OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST BE AT LEAST #### ' write(*,*) ' #### TWICE THE SYNCHRONISATION INTERVAL #### ' stop endif if (mod(lrecoutstep,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' stop endif if ((lrecoutstep/lsynctime).lt.2) then write(*,*) ' #### FLEXPART MODEL ERROR! INTERVAL BETWEEN #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST BE AT LEAST #### ' write(*,*) ' #### TWICE THE SYNCHRONISATION INTERVAL #### ' stop endif if (mod(lrecoutsample,lsynctime).ne.0) then write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### ' write(*,*) ' #### RECEPTOR OUTPUT MUST BE A MULTIPLE #### ' write(*,*) ' #### OF THE SYNCHRONISATION INTERVAL #### ' stop endif endif ! ldirect ! Switch for LCM mode !******************************************************************* if ( (lcmoutput.ne.0) .and. ((.not. ind_source.eq.1).or. & (.not. ind_receptor.eq.1).or.(.not. iout.eq.2).or. & (.not. ldirect.eq.1).or.(.not. mdomainfill.eq.1)) ) then write(*,*) 'LCM output requested, but one of the following options' write(*,*) 'is not correctly set in COMMAND:' write(*,*) 'ind_source =', ind_source, 'should be set to 1' write(*,*) 'ind_receptor =', ind_receptor, 'should be set to 1' write(*,*) 'iout =', iout, 'should be set to 2' write(*,*) 'ldirect =', ldirect, 'should be set to 1' write(*,*) 'mdomainfill =', mdomainfill, 'should be set to 1' error stop endif if (lcmoutput.eq.0) then llcmoutput=.false. else llcmoutput=.true. endif write(*,*) 'Switch for LCM output LCMOUTPUT = ',llcmoutput ! Compute modeling time in seconds and beginning date in Julian date !******************************************************************* outstep=real(abs(loutstep)) if (ldirect.eq.1) then bdate=juldate(ibdate,ibtime) edate=juldate(iedate,ietime) ideltas=nint((edate-bdate)*86400.) else if (ldirect.eq.-1) then loutaver=-1*loutaver loutstep=-1*loutstep loutsample=-1*loutsample lsynctime=-1*lsynctime bdate=juldate(iedate,ietime) edate=juldate(ibdate,ibtime) ideltas=nint((edate-bdate)*86400.) else write(*,*) ' #### FLEXPART MODEL ERROR! DIRECTION IN #### ' write(*,*) ' #### FILE "COMMAND" MUST BE EITHER -1 OR 1. #### ' error stop endif return 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(1)(1:length(1)) error stop 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readcommand subroutine readdepo !***************************************************************************** ! * ! Reads dry deposition parameters needed by the procedure of Wesely (1989). * ! Wesely (1989): Parameterization of surface resistances to gaseous * ! dry deposition in regional-scale numerical models. * ! Atmos. Environ. 23, 1293-1304. * ! * ! * ! AUTHOR: Andreas Stohl, 19 May 1995 * ! * !***************************************************************************** ! * ! Variables: * ! * ! rcl(maxspec,5,9) [s/m] Lower canopy resistance * ! rgs(maxspec,5,9) [s/m] Ground resistance * ! rlu(maxspec,5,9) [s/m] Leaf cuticular resistance * ! rm(maxspec) [s/m] Mesophyll resistance, set in readreleases * ! ri(maxspec) [s/m] Stomatal resistance * ! * ! Constants: * ! * !***************************************************************************** implicit none ! FOR THIS SUBROUTINE, numclass=9 IS ASSUMED !******************************************* real :: rluh(5,numclass),rgssh(5,numclass),rgsoh(5,numclass) real :: rclsh(5,numclass),rcloh(5,numclass) integer :: i,j,ic logical :: ios, ios2 character(12) :: file_sfcdepo ! Read deposition constants related with landuse and seasonal category !********************************************************************* file_sfcdepo='sfcdepo.t' inquire(file=path(1)(1:length(1))//trim(file_sfcdepo),exist=ios) if (.not. ios) then file_sfcdepo='surfdepo.t' inquire(file=path(1)(1:length(1))//trim(file_sfcdepo),exist=ios2) if (ios2) then write(*,*) 'WARNING: surfdepo.t should be renamed to sfcdepo.t' write(*,*) 'Continuing with surfdepo.t' endif endif open(unitwesely,file=path(1)(1:length(1))//trim(file_sfcdepo), & status='old',err=999) do i=1,16 read(unitwesely,*) end do do i=1,5 read(unitwesely,*) read(unitwesely,'(8x,13f8.0)') (ri(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rluh(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rac(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rgssh(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rgsoh(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rclsh(i,j),j=1,numclass) read(unitwesely,'(8x,13f8.0)') (rcloh(i,j),j=1,numclass) end do ! TEST ! do 31 i=1,5 ! ri(i,13)=ri(i,5) ! rluh(i,13)=rluh(i,5) ! rac(i,13)=rac(i,5) ! rgssh(i,13)=rgssh(i,5) ! rgsoh(i,13)=rgsoh(i,5) ! rclsh(i,13)=rclsh(i,5) ! rcloh(i,13)=rcloh(i,5) !31 continue ! TEST ! Sabine Eckhardt, Dec 06, set resistances of 9999 to 'infinite' (1E25) do i=1,5 do j=1,numclass if (ri(i,j).eq.9999.) ri(i,j)=1.E25 if (rluh(i,j).eq.9999.) rluh(i,j)=1.E25 if (rac(i,j).eq.9999.) rac(i,j)=1.E25 if (rgssh(i,j).eq.9999.) rgssh(i,j)=1.E25 if (rgsoh(i,j).eq.9999.) rgsoh(i,j)=1.E25 if (rclsh(i,j).eq.9999.) rclsh(i,j)=1.E25 if (rcloh(i,j).eq.9999.) rcloh(i,j)=1.E25 end do end do do i=1,5 do j=1,numclass ri(i,j)=max(ri(i,j),0.001) rluh(i,j)=max(rluh(i,j),0.001) rac(i,j)=max(rac(i,j),0.001) rgssh(i,j)=max(rgssh(i,j),0.001) rgsoh(i,j)=max(rgsoh(i,j),0.001) rclsh(i,j)=max(rclsh(i,j),0.001) rcloh(i,j)=max(rcloh(i,j),0.001) end do end do close(unitwesely) ! Compute additional parameters !****************************** do ic=1,nspec if (reldiff(ic).gt.0.) then ! gas is dry deposited do i=1,5 do j=1,numclass rlu(ic,i,j)=rluh(i,j)/(1.e-5*henry(ic)+f0(ic)) rgs(ic,i,j)=1./(henry(ic)/(10.e5*rgssh(i,j))+f0(ic)/ & rgsoh(i,j)) rcl(ic,i,j)=1./(henry(ic)/(10.e5*rclsh(i,j))+f0(ic)/ & rcloh(i,j)) end do end do endif end do return 999 write(*,*) '### FLEXPART ERROR! FILE ###' write(*,*) '### sfcdepo.t DOES NOT EXIST. ###' error stop end subroutine readdepo subroutine readlanduse !***************************************************************************** ! * ! Reads the landuse inventory into memory and relates it to Leaf Area * ! Index and roughness length. * ! * ! AUTHOR: Andreas Stohl, 10 January 1994 * ! * !***************************************************************************** ! * ! Variables: * ! i loop indices * ! landinvent(1200,600,13) area fractions of 13 landuse categories * ! LENGTH(numpath) length of the path names * ! PATH(numpath) contains the path names * ! unitland unit connected with landuse inventory * ! * ! ----- * ! Sabine Eckhardt, Dec 06 - new landuse inventary * ! after * ! Belward, A.S., Estes, J.E., and Kline, K.D., 1999, * ! The IGBP-DIS 1-Km Land-Cover Data Set DISCover: * ! A Project Overview: Photogrammetric Engineering and Remote Sensing, * ! v. 65, no. 9, p. 1013-1020 * ! * ! LANDUSE CATEGORIES: * ! * ! 1 Urban land * ! 2 Agricultural land * ! 3 Range land * ! 4 Deciduous forest * ! 5 Coniferous forest * ! 6 Mixed forest including wetland * ! 7 water, both salt and fresh * ! 8 barren land mostly desert * ! 9 nonforested wetland * ! 10 mixed agricultural and range land * ! 11 rocky open areas with low growing shrubs * ! 12 ice * ! 13 rainforest * ! * !***************************************************************************** use drydepo_mod implicit none integer :: ix,jy,i,k,lu_cat,lu_perc integer(kind=1) :: ilr integer(kind=1) :: ilr_buffer(2160000) integer :: il,irecread real :: rlr, r2lr logical :: ios,ios2 character(12) :: file_sfcdata ! Read landuse inventory !*********************** ! The landuse information is saved in a compressed format and written ! out by records of the length of 1 BYTE. Each grid cell consists of 3 ! Bytes, which include 3 landuse categories (val 1-13 and 16 percentage ! categories) So one half byte is used to store the Landusecat the other ! for the percentageclass in 6.25 steps (100/6.25=16) ! e.g. ! 4 3 percentage 4 = 4*6.25 => 25% landuse class 3 ! 2 1 percentage 2 = 2*6.25 => 13% landuse class 1 ! 1 12 percentage 1 = 1*6.26 => 6.25% landuse class 12 open(unitland,file=path(1)(1:length(1))//'IGBP_int1.dat',status='old', & form='UNFORMATTED', err=998, convert='little_endian') read (unitland) (ilr_buffer(i),i=1,2160000) close(unitland) irecread=1 do ix=1,1200 do jy=1,600 ! the 3 most abundant landuse categories in the inventory ! first half byte contains the landuse class ! second half byte contains the respective percentage do k=1,3 ! 1 byte is read ilr=ilr_buffer(irecread) ! ilr=0 irecread=irecread+1 ! as only signed integer values exist an unsigned value is constructed if (ilr.lt.0) then il=ilr+256 else il=ilr endif ! dividing by 16 has the effect to get rid of the right half of the byte ! so just the left half remains, this corresponds to a shift right of 4 ! bits rlr=real(il)/16. lu_cat=int(rlr) ! the left half of the byte is substracted from the whole in order to ! get only the right half of the byte r2lr=rlr-int(rlr) ! shift left by 4 lu_perc=int(r2lr*16.) landinvent(ix,jy,k)=lu_cat landinvent(ix,jy,k+3)=lu_perc ! if ((jy.lt.10).and.(ix.lt.10)) write(*,*) 'reading: ',ix,jy,lu_cat,lu_perc end do end do end do ! Read relation landuse,z0 !***************************** file_sfcdata='sfcdata.t' inquire(file=path(1)(1:length(1))//trim(file_sfcdata),exist=ios) if (.not. ios) then file_sfcdata='surfdata.t' inquire(file=path(1)(1:length(1))//trim(file_sfcdata),exist=ios2) if (ios2) then write(*,*) 'WARNING: surfdata.t should be renamed to sfcdata.t' write(*,*) 'Continuing with surfdata.t' endif endif open(unitsfcdata,file=path(1)(1:length(1))//trim(file_sfcdata), & status='old',err=999) do i=1,4 read(unitsfcdata,*) end do do i=1,numclass read(unitsfcdata,'(45x,f15.3)') z0(i) end do close(unitsfcdata) return ! Issue error messages !********************* 998 write(*,*) ' #### FLEXPART ERROR! FILE ####' write(*,*) ' #### ', path(1)(1:length(1))//'IGBP_int1.dat' write(*,*) " #### (LANDUSE INVENTORY) COULD NOT BE OPENED ####" error stop 999 write(*,*) ' #### FLEXPART ERROR! FILE ####' write(*,*) ' #### ', path(1)(1:length(1))//'sfcdata.t' write(*,*) ' #### DOES NOT EXIST. ####' error stop end subroutine readlanduse subroutine readoutgrid !***************************************************************************** ! * ! This routine reads the user specifications for the output grid. * ! * ! Author: A. Stohl * ! * ! 4 June 1996 * ! HSO, 1 July 2014 ! Added optional namelist input ! * !***************************************************************************** ! * ! Variables: * ! dxout,dyout grid distance * ! numxgrid,numygrid,numzgrid grid dimensions * ! outlon0,outlat0 lower left corner of grid * ! outheight(maxzgrid) height levels of output grid [m] * ! * ! Constants: * ! unitoutgrid unit connected to file OUTGRID * ! * !***************************************************************************** use outgrid_mod implicit none integer :: i,j,stat real :: outhelp,xr,xr1,yr,yr1 real,parameter :: eps=1.e-4 ! namelist variables integer, parameter :: maxoutlev=500 integer :: ios real,allocatable, dimension (:) :: outheights ! declare namelist namelist /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 ! Open the OUTGRID file and read output grid specifications !********************************************************** open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', & form='formatted',err=999) ! try namelist input read(unitoutgrid,outgrid,iostat=ios) close(unitoutgrid) if ((dxout.le.0).or.(ios.ne.0)) then ios=1 open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999) call skplin(5,unitoutgrid) ! 1. 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 ! Check validity of output grid (shall be within model domain) !************************************************************* xr=outlon0+real(numxgrid)*dxout yr=outlat0+real(numygrid)*dyout xr1=xlon0+real(nxmin1)*dx yr1=ylat0+real(nymin1)*dy if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) & .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then write(*,*) outlon0,outlat0 write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT ####' write(*,*) ' #### GRID IS OUTSIDE MODEL DOMAIN. CHANGE ####' write(*,*) ' #### FILE OUTGRID IN DIRECTORY ####' write(*,'(a)') path(1)(1:length(1)) error stop endif ! 2. Count Vertical levels of output grid !**************************************** if (ios.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 (ios.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 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) close(unitoutgrid) endif ! Check whether vertical levels are specified in ascending order !*************************************************************** 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 end do ! Determine the half levels, i.e. middle levels of the output grid !***************************************************************** outheighthalf(1)=outheight(1)*0.5 do j=2,numzgrid outheighthalf(j)=(outheight(j-1)+outheight(j))*0.5 end do xoutshift=xlon0-outlon0 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' return 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(1)(1:length(1)) error stop 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readoutgrid subroutine readoutgrid_nest !***************************************************************************** ! * ! This routine reads the user specifications for the output nest. * ! * ! Author: A. Stohl * ! * ! 4 June 1996 * ! * !***************************************************************************** ! * ! Variables: * ! dxoutn,dyoutn grid distances of output nest * ! numxgridn,numygridn,numzgrid nest dimensions * ! outlon0n,outlat0n lower left corner of nest * ! outheight(maxzgrid) height levels of output grid [m] * ! * ! Constants: * ! unitoutgrid unit connected to file OUTGRID * ! * !***************************************************************************** use outgrid_mod implicit none integer :: stat real :: xr,xr1,yr,yr1 real,parameter :: eps=1.e-4 integer :: ios ! declare namelist namelist /outgridn/ & outlon0n,outlat0n, & numxgridn,numygridn, & dxoutn,dyoutn ! helps identifying failed namelist input dxoutn=-1.0 ! Open the OUTGRID file and read output grid specifications !********************************************************** open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999) ! try namelist input read(unitoutgrid,outgridn,iostat=ios) close(unitoutgrid) if ((dxoutn.le.0).or.(ios.ne.0)) then open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999) call skplin(5,unitoutgrid) ! 1. Read horizontal grid specifications !**************************************** call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,f11.4)') outlon0n call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,f11.4)') outlat0n call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,i5)') numxgridn call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,i5)') numygridn call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,f12.5)') dxoutn call skplin(3,unitoutgrid) read(unitoutgrid,'(4x,f12.5)') dyoutn close(unitoutgrid) endif ! 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) 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' ! Check validity of output grid (shall be within model domain) !************************************************************* xr=outlon0n+real(numxgridn)*dxoutn yr=outlat0n+real(numygridn)*dyoutn 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 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)) error stop endif xoutshiftn=xlon0-outlon0n youtshiftn=ylat0-outlat0n return 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(1)(1:length(1)) error stop 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readoutgrid_nest subroutine readpaths !***************************************************************************** ! * ! Reads the pathnames, where input/output files are expected to be. * ! The file pathnames must be available in the current working directory. * ! * ! Author: A. Stohl * ! * ! 1 February 1994 * ! last modified * ! HS, 7.9.2012 * ! option to give pathnames file as command line option * ! * !***************************************************************************** ! * ! Variables: * ! length(numpath) lengths of the path names * ! path(numpath) pathnames of input/output files * ! * ! Constants: * ! numpath number of pathnames to be read in * ! * !***************************************************************************** implicit none integer :: i character(256) :: string_test character(1) :: character_test ! Read the pathname information stored in unitpath !************************************************* open(unitpath,file=trim(pathfile),status='old',err=999) do i=1,numpath read(unitpath,'(a)',err=998) path(i) length(i)=index(path(i),' ')-1 string_test = path(i) character_test = string_test(length(i):length(i)) !print*, 'character_test, string_test ', character_test, string_test if ((character_test .NE. '/') .AND. (i .LT. 4)) then print*, 'WARNING: path not ending in /' print*, path(i) path(i) = string_test(1:length(i)) // '/' length(i)=length(i)+1 print*, 'fix: padded with /' print*, path(i) print*, 'length(i) increased 1' endif end do ! Check whether any nested subdomains are to be used !*************************************************** do i=1,maxnests ! ESO 2016 Added 'end'/'err' in case user forgot '====' at end of file and ! maxnests > numbnests read(unitpath,'(a)', end=30, err=30) path(numpath+2*(i-1)+1) read(unitpath,'(a)', end=30, err=30) path(numpath+2*(i-1)+2) if (path(numpath+2*(i-1)+1)(1:5).eq.'=====') goto 30 length(numpath+2*(i-1)+1)=index(path(numpath+2*(i-1)+1),' ')-1 length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1 end do ! Determine number of available nested domains !********************************************* 30 numbnests=i-1 close(unitpath) return 998 write(*,*) ' #### TRAJECTORY MODEL ERROR! ERROR WHILE #### ' write(*,*) ' #### READING FILE PATHNAMES. #### ' error stop 999 write(*,*) ' #### TRAJECTORY MODEL ERROR! FILE "pathnames"#### ' write(*,*) ' #### CANNOT BE OPENED IN THE CURRENT WORKING #### ' write(*,*) ' #### DIRECTORY. #### ' error stop end subroutine readpaths subroutine readreceptors !***************************************************************************** ! * ! This routine reads the user specifications for the receptor points. * ! * ! Author: A. Stohl * ! 1 August 1996 * ! * ! HSO, 14 August 2013: Added optional namelist input ! PS, 2/2015: access= -> position= ! PS, 6/2015: variable names, simplify code ! PS, 3/2023: remove position=append, makes no sense for new file * ! * !***************************************************************************** ! * ! Variables: * ! receptorarea area of dx*dy at location of receptor * ! receptorname names of receptors * ! xreceptor,yreceptor coordinates of receptor points * ! * ! Constants: * ! unitreceptor unit connected to file RECEPTORS * ! * !***************************************************************************** implicit none integer :: j real :: xm,ym character(len=16) :: receptor integer :: ios real :: lon,lat,alt ! for namelist input, lon/lat are used instead of x,y real(kind=dp) :: time ! real,allocatable,dimension(:) :: tmpxrec,tmpyrec,tmprecarea ! character(len=16),allocatable,dimension(:) :: tmprecname ! declare namelist namelist /receptors/ & receptor, lon, lat, alt, time numreceptor=0 ! Initialise numreceptor ! For backward runs no receptor output !************************************* if (ldirect.lt.0) then return endif ! Open the RECEPTORS file and read output grid specifications !************************************************************ open (unitreceptor,file=trim(path(1))//'RECEPTORS',form='formatted', & status='old',err=999) lon = -999. lat = -999. time = -999. lrecregular = .false. ! try namelist input read(unitreceptor,receptors,iostat=ios) close (unitreceptor) if ((lon.lt.-900).or.(ios.ne.0)) then go to 999 else ! only namelist input possible ! prepare namelist output if requested if (nmlout) open(unitreceptorout,file=trim(path(2))// & 'RECEPTORS.namelist',status='replace',err=1000) ! Get number of receptors !************************ open (unitreceptor,file=trim(path(1))//'RECEPTORS',status='old',err=999) j=0 do while (ios.eq.0) lon=-999.9 read(unitreceptor,receptors,iostat=ios) if ((lon.lt.-900).or.(ios.ne.0)) exit ! skip receptors for which a timestamp is given but are not in simulation window if ((time.ne.-999.).and.((time.lt.bdate).or.(time.ge.edate))) cycle j=j+1 end do numreceptor=j write(*,*) 'Number of receptors: ',numreceptor close (unitreceptor) ! Allocate arrays !**************** allocate(receptorname(numreceptor),xreceptor(numreceptor),& yreceptor(numreceptor),zreceptor(numreceptor),& treceptor(numreceptor),receptorarea(numreceptor)) ! Read the names and coordinates of the receptors !************************************************ open (unitreceptor,file=trim(path(1))//'RECEPTORS',status='old',iostat=ios) j=0 do while (ios.eq.0) lon=-999.9 read(unitreceptor,receptors,iostat=ios) if ((lon.lt.-900).or.(ios.ne.0)) exit ! read error ! skip receptors for which a timestamp is given but are not in simulation window if ((time.ne.-999.).and.((time.lt.bdate).or.(time.ge.edate))) cycle j=j+1 receptorname(j)=receptor xreceptor(j)=(lon-xlon0)/dx ! transform to grid coordinates yreceptor(j)=(lat-ylat0)/dy zreceptor(j)=alt if (time.ne.-999.) then treceptor(j)=int((time-bdate)*24.*3600.) ! time in sec ! round to nearest 10 seconds treceptor(j)=nint(real(treceptor(j))/10.)*10 else treceptor(j)=-999 endif xm=r_earth*cos(lat*pi/180.)*dx/180.*pi ym=r_earth*dy/180.*pi receptorarea(j)=xm*ym ! write receptors in namelist format to output directory if requested if (nmlout) write(unitreceptorout,nml=receptors) end do close (unitreceptor) if (nmlout) close (unitreceptorout) endif ! if not timestamp given in namelist assume regular output ! according to COMMAND file settings if (.not.any(treceptor.ne.-999)) then lrecregular=.true. endif !! testing ! write(*,*) 'readreceptors: ' ! do j=1,numreceptor ! print*, 'receptorname = ',receptorname(j) ! print*, 'xreceptor, yreceptor, zreceptor = ',xreceptor(j), yreceptor(j), zreceptor(j) ! print*, 'treceptor = ',treceptor(j) ! end do !! return 999 write(*,*) ' #### FLEXPART WARNING: File RECEPTORS cannot be opened #### ' write(*,*) ' #### in directory '//trim(path(1))//' #### ' write(*,*) ' #### continuing without RECEPTOR output #### ' numreceptor=0 return 1000 write(*,*) ' #### FLEXPART MODEL ERROR! File "RECEPTORS" #### ' write(*,*) ' #### cannot be opened in the output directory #### ' write(*,'(a)') ' #### '//trim(path(2)) write(*,*) ' #### either write perm missing or old file exists ####' error stop end subroutine readreceptors subroutine readreleases !***************************************************************************** ! * ! This routine reads the release point specifications for the current * ! model run. Several release points can be used at the same time. * ! * ! Author: A. Stohl * ! * ! 18 May 1996 * ! * ! Update: 29 January 2001 * ! Release altitude can be either in magl or masl * ! HSO, 12 August 2013 ! Added optional namelist input ! * !***************************************************************************** ! * ! Variables: * ! decay decay constant of species * ! dquer [um] mean particle diameters * ! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass* ! are between 0.1*dquer and 10*dquer * ! ireleasestart, ireleaseend [s] starting time and ending time of each * ! release * ! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint* ! is in hPa * ! npart number of particles to be released * ! nspec number of species to be released * ! density [kg/m3] density of the particles * ! rm [s/m] Mesophyll resistance * ! species name of species * ! xmass total mass of each species * ! xpoint1,ypoint1 geograf. coordinates of lower left corner of release * ! area * ! xpoint2,ypoint2 geograf. coordinates of upper right corner of release * ! area * ! weta_gas, wetb_gas parameters for below-cloud scavenging (gas) * ! crain_aero, csnow_aero parameters for below-cloud scavenging (aerosol) * ! ccn_aero, in_aero parameters for in-cloud scavenging (aerosol) * ! zpoint1,zpoint2 height range, over which release takes place * ! num_min_discrete if less, release cannot be randomized and happens at * ! time mid-point of release interval * ! lroot true if serial version, or if MPI and root process * ! * !***************************************************************************** use point_mod use xmass_mod use drydepo_mod implicit none integer :: numpartmax,i,j,id1,it1,id2,it2,stat,irel,ispc,nsettle integer,parameter :: num_min_discrete=100 real :: releaserate,cun real(kind=dp) :: jul1,jul2,julm real,parameter :: eps2=1.e-9 character(len=50) :: line ! help variables for namelist reading integer :: numpoints, parts, ios, nspec_init integer*2 :: zkind integer :: idate1, itime1, idate2, itime2 real :: lon1,lon2,lat1,lat2,z1,z2 character*40 :: comment integer,parameter :: unitreleasesout=2 real,allocatable, dimension (:) :: mass integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2 real,allocatable,dimension(:) :: vsh,fracth,schmih ! declare namelists namelist /releases_ctrl/ & nspec, & specnum_rel namelist /release/ & idate1, itime1, & idate2, itime2, & lon1, lon2, & lat1, lat2, & z1, z2, & zkind, & mass, & parts, & comment numpoint=0 nspec_init=50 ! necessary to allocate specnum_rel, would be cleaner to change RELEASES ! in one extra namelist allocate(specnum_rel(nspec_init),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_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=ios) if (nspec.gt.nspec_init) then write(*,*) 'Requested number of species:',nspec error stop 'More than 50 species at a time is not possible when using & & RELEASES. You can do this with the part_ic.nc option (IPIN=3).' end if ! Allocate fields with maxspec maxspec=nspec call alloc_com() ! allocate with maxspec for first input loop allocate(mass(maxspec),stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not allocate mass' if (ios.ne.0) then backspace(unitreleases) read(unitreleases,fmt='(A)') line if (lroot) write(*,*) & 'Invalid line in RELEASES: '//trim(line) !cgz; Check if the number of species in RELEASES_CTRL is larger than the maximum number of species in par_mod if ((lroot) .and. nspec.gt.maxspec) goto 994 end if ! 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 ((ios.ne.0).or.(nspec.lt.0)) then if (lroot) write(*,*) 'RELEASE either having unrecognised entries, & &or in old format, please update to namelist format.' error stop else if ((ipin.ne.3).and.(ipin.ne.4)) then ! Not necessary to read releases when using part_ic.nc ios=0 do while (ios.eq.0) idate1=-1 read(unitreleases,release,iostat=ios) if ((idate1.lt.0).or.(ios.ne.0)) then ios=1 else numpoint=numpoint+1 endif end do ios=0 else numpoint=1 endif endif rewind(unitreleases) 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) ! 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 !now save the information DEP=.false. DRYDEP=.false. WETDEP=.false. CLREA=.false. LDECAY=.false. do i=1,maxspec DRYDEPSPEC(i)=.false. WETDEPSPEC(i)=.false. end do ! namelist output if (nmlout.and.lroot) then write(unitreleasesout,nml=releases_ctrl) endif do i=1,nspec call readspecies(specnum_rel(i),i) end do ! Allocate fields that depend on ndia call alloc_com_ndia do i=1,nspec ! Allocate temporary memory necessary for the different diameter bins !******************************************************************** allocate( vsh(ndia(i)),fracth(ndia(i)),schmih(ndia(i)),stat=stat) if (stat.ne.0) error stop "Could not allocate vsh,fracth,schmih" ! Molecular weight !***************** 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.' write(*,*) 'Check table SPECIES or choose concentration' write(*,*) 'output instead if molar weight is not known.' error stop endif ! Radioactive decay !****************** if (decay(i).gt.0) then LDECAY=.true. decay(i)=0.693147/decay(i) !conversion half life to decay constant endif ! Dry deposition of gases !************************ 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 if (density(i).gt.0.) then ! Additional parameters call part0(dquer(i),dsigma(i),density(i),ndia(i),fracth,schmih,cun,vsh) do j=1,ndia(i) 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) endif ! Dry deposition for constant deposition velocity !************************************************ dryvel(i)=dryvel(i)*0.01 ! conversion to m/s ! Check if wet deposition shall be calculated !********************************************* ! 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 (lroot) then write (*,*) ' Below-cloud scavenging: ON' ! write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i end if else if (lroot) write (*,*) ' Below-cloud scavenging: OFF' endif ! 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 (lroot) then write (*,*) ' In-cloud scavenging: ON' ! write (*,*) 'In-cloud scavenging coefficients: ',& ! &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i end if else if (lroot) write (*,*) ' In-cloud scavenging: OFF' endif if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then DRYDEP=.true. DRYDEPSPEC(i)=.true. endif deallocate(vsh,fracth,schmih) end do ! end loop over species if (WETDEP.or.DRYDEP) DEP=.true. ! Check if chemical reaction shall be calculated !*********************************************** if (any(reaccconst(:,:).gt.0.)) then CLREA=.true. if (lroot) write (*,*) ' Chemical reactions switched on' endif ! Check if emissions shall be used !********************************* if (any(emis_path(:).ne."")) then LEMIS=.true. if (lroot) write(*,*) ' Emissions switched on' endif ! Not necessary to read releases when using part_ic.nc !***************************************************** if ((ipin.eq.3).or.(ipin.eq.4)) then maxpointspec_act=1 return endif ! Read specifications for each release point !******************************************* numpoints=numpoint numpoint=0 numpartmax=0 releaserate=0. 101 numpoint=numpoint+1 if (numpoint.gt.numpoints) goto 250 zkind = 1 mass = 0 parts = 0 comment = ' ' 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 zpoint1(numpoint)=z1 zpoint2(numpoint)=z2 kindz(numpoint)=zkind do i=1,nspec xmass(numpoint,i)=mass(i) end do npart(numpoint)=parts compoint(min(1001,numpoint))=comment ! namelist output if (nmlout.and.lroot) then write(unitreleasesout,nml=release) endif ! If a release point contains no particles, stop and issue error message !*********************************************************************** if (npart(numpoint).eq.0) then write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'RELEASES file is corrupt.' write(*,*) 'At least for one release point, there are zero' write(*,*) 'particles released. Make changes to RELEASES.' error stop endif ! If FLEXPART is run for backward deposition force zpoint !********************************************************************* if (WETBKDEP) then zpoint1(numpoint)=0. zpoint2(numpoint)=20000. kindz(numpoint)=1 endif if (DRYBKDEP) then zpoint1(numpoint)=0. zpoint2(numpoint)=2.*href kindz(numpoint)=1 endif ! Check whether x coordinates of release point are within model domain !********************************************************************* 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 !****************************************************************** jul1=juldate(id1,it1) jul2=juldate(id2,it2) julm=(jul1+jul2)*0.5 if (jul1.gt.jul2) then write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'Release stops before it begins.' write(*,*) 'Make changes to file RELEASES.' error stop endif if (mdomainfill.eq.0) then ! no domain filling if (ldirect.eq.1) then if (((jul1.lt.bdate).or.(jul2.gt.edate)).and.(ipin.eq.0)) then write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'Release starts before simulation begins or ends' write(*,*) 'after simulation stops.' write(*,*) 'Make files COMMAND and RELEASES consistent.' error stop endif if (npart(numpoint).gt.num_min_discrete) then ireleasestart(numpoint)=int((jul1-bdate)*86400.) ireleaseend(numpoint)=int((jul2-bdate)*86400.) else ireleasestart(numpoint)=int((julm-bdate)*86400.) ireleaseend(numpoint)=int((julm-bdate)*86400.) endif else if (ldirect.eq.-1) then if (((jul1.lt.edate).or.(jul2.gt.bdate)).and.(ipin.eq.0)) then write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'Release starts before simulation begins or ends' write(*,*) 'after simulation stops.' write(*,*) 'Make files COMMAND and RELEASES consistent.' error stop endif if (npart(numpoint).gt.num_min_discrete) then ireleasestart(numpoint)=int((jul1-bdate)*86400.) ireleaseend(numpoint)=int((jul2-bdate)*86400.) else ireleasestart(numpoint)=int((julm-bdate)*86400.) ireleaseend(numpoint)=int((julm-bdate)*86400.) endif endif endif ! Determine the release rate (particles per second) and total number ! of particles released during the simulation !******************************************************************* if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then releaserate=releaserate+real(npart(numpoint))/ & real(ireleaseend(numpoint)-ireleasestart(numpoint)) else releaserate=99999999. endif numpartmax=numpartmax+npart(numpoint) goto 101 250 close(unitreleases) if (nmlout.and.lroot) then close(unitreleasesout) endif !if (lroot) write (*,*) 'Particles allocated (maxpart) : ',maxpart if (lroot) write (*,*) 'Particles released (numpartmax): ',numpartmax numpoint=numpoint-1 if (ioutputforeachrelease.eq.1) then maxpointspec_act=numpoint else maxpointspec_act=1 endif ! Disable settling if more than 1 species at any release point ! or if MQUASILAG and more than one species !************************************************************* if (mquasilag.ne.0) then if (nspec.gt.1) lsettling=.false. else do irel=1,numpoint nsettle=0 do ispc=1,nspec if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1 end do if (nsettle.gt.1) lsettling=.false. end do end if if (lroot) then if (.not.lsettling) then write(*,*) 'WARNING: more than 1 species per release point, settling & &disabled' end if end if ! Check, whether the total number of particles may exceed totally allowed ! number of particles at some time during the simulation !************************************************************************ ! if (releaserate.gt. & ! 0.99*real(maxpart)/real(lage(nageclass))) then ! if (numpartmax.gt.maxpart.and.lroot) then ! write(*,*) '#####################################################' ! write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' ! write(*,*) '#### ####' ! write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####' ! write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####' ! write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####' ! write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####' ! write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####' ! write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####' ! write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####' ! write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####' ! write(*,*) '#####################################################' ! write(*,*) 'Maximum release rate may be: ',releaserate, & ! ' particles per second' ! 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(*,*) 'Maximum allowed number of particles is: ',maxpart ! endif ! endif if (lroot) then write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec)) end if return 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(*,*) '#####################################################' error stop 999 write(*,*) '#####################################################' write(*,*) ' FLEXPART MODEL SUBROUTINE READRELEASES: ' write(*,*) write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS' write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT' write(*,*) 'PERMITTED FOR ANY ACCESS' write(*,*) '#####################################################' error stop 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readreleases subroutine readspecies(id_spec,pos_spec) !***************************************************************************** ! * ! This routine reads names and physical constants of chemical species/ * ! radionuclides given in the parameter pos_spec * ! * ! Author: A. Stohl * ! * ! 11 July 1996 * ! * ! Changes: * ! N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging * ! * ! HSO, 13 August 2013 ! added optional namelist input ! * ! R. Thompson, 18.01.2024 * ! variables for LCM * ! * !***************************************************************************** ! * ! Variables: * ! decaytime(maxtable) half time for radiological decay * ! specname(maxtable) names of chemical species, radionuclides * ! weta_gas, wetb_gas Parameters for below-cloud scavenging of gasses * ! crain_aero,csnow_aero Parameters for below-cloud scavenging of aerosols * ! ccn_aero,in_aero Parameters for in-cloud scavenging of aerosols * ! reaccconst Chemical reaction rate constant C * ! reacdconst Chemical reaction rate constant D * ! reacnconst Chemical reaction rate constant n * ! id_spec SPECIES number as referenced in RELEASE file * ! id_pos position where SPECIES data shall be stored * ! Constants: * ! * !***************************************************************************** implicit none integer :: i, pos_spec,j integer :: idow,ihour,id_spec character(len=3) :: aspecnumb character(len=16) :: pspecies, pemis_name character(len=256) :: pemis_path, pemis_file integer :: pemis_unit real :: pemis_coeff character(len=50) :: line real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer real :: pdsigma, pdryvel, pweightmolar, pdia character(len=10), allocatable, dimension(:) :: preactions real, allocatable, dimension(:) :: pcconst, pdconst, pnconst real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero real :: pohcconst,pohdconst,pohnconst ! deprecated real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24) !integer :: pndia integer :: ios integer :: pshape,porient ! Daria Tatsii: species shape properties real ::pla,pia,psa,f,e,paspectratio ! declare namelist namelist /species_params/ & pspecies, pdecay, pweta_gas, pwetb_gas, & pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, & preldiff, phenry, pf0, pdensity, pdquer, pdia, & pdsigma, pdryvel, pweightmolar, pohnconst, & preactions, pcconst, pdconst, pnconst, pohcconst, pohdconst, & pemis_path, pemis_file, pemis_name, pemis_unit, pemis_coeff, & parea_dow, parea_hour, ppoint_dow, ppoint_hour, & pshape, paspectratio, pla, pia, psa, porient !pndia, ! allocate reaction variables allocate(preactions(maxreagent)) allocate(pcconst(maxreagent)) allocate(pdconst(maxreagent)) allocate(pnconst(maxreagent)) if (.not.allocated(reaccconst)) then allocate(reaccconst(maxreagent,nspec)) allocate(reacdconst(maxreagent,nspec)) allocate(reacnconst(maxreagent,nspec)) reaccconst(:,:)=-9.99e-9 reacdconst(:,:)=-9.99 reacnconst(:,:)=-9.99 allocate(emis_path(nspec)) allocate(emis_file(nspec)) allocate(emis_name(nspec)) allocate(emis_unit(nspec)) allocate(emis_coeff(nspec)) endif pspecies="" ! read failure indicator value pdecay=-999.9 pweta_gas=-9.9E-09 pwetb_gas=0.0 pcrain_aero=-9.9E-09 pcsnow_aero=-9.9E-09 pccn_aero=-9.9E-09 pin_aero=-9.9E-09 preldiff=-9.9 phenry=-9.9 pf0=0.0 pdensity=-9.9E09 pdquer=-9.9 pdia=-9.9 pdsigma=0.0 !pndia=1 pdryvel=-9.99 preactions(:)="" pcconst(:)=-9.99e-9 pdconst(:)=-9.99 pnconst(:)=-9.99 pohcconst=-9.99 pohdconst=-9.99 pohnconst=-9.99 pweightmolar=-999.9 parea_dow=-999.9 parea_hour=-999.9 ppoint_dow=-999.9 ppoint_hour=-999.9 pshape=0 ! 0 for sphere, 1 for other shapes paspectratio=-1. pla=-1. ! longest axis in micrometer pia=-1. ! Intermediate axis psa=-1. ! Smallest axis porient=0 ! 0 for horizontal, 1 for random pemis_path="" ! read failure indicator value pemis_file="" ! read failure indicator value pemis_name="" ! read failure indicator value pemis_unit=0 pemis_coeff=1. do j=1,24 ! initialize everything to no variation parea_hour(j)=1. ppoint_hour(j)=1. area_hour(pos_spec,j)=1. point_hour(pos_spec,j)=1. end do do j=1,7 parea_dow(j)=1. ppoint_dow(j)=1. area_dow(pos_spec,j)=1. point_dow(pos_spec,j)=1. end do ! Open the SPECIES file and read species names and properties !************************************************************ specnum(pos_spec)=id_spec write(aspecnumb,'(i3.3)') specnum(pos_spec) open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb, & status='old',form='formatted',err=998) write(*,*) 'reading SPECIES',specnum(pos_spec) ASSSPEC=.FALSE. read(unitspecies,species_params,iostat=ios) ! check on which line of species file problem occurs if (ios.ne.0) then backspace(unitspecies) read(unitspecies,fmt='(A)') line if (lroot) write(*,*) & 'Invalid line in species: '//trim(line) end if close(unitspecies) if ((len(trim(pspecies)).eq.0).or.(ios.ne.0)) then ! no namelist found if (lroot) then write(*,*) "FLEXPART ERROR: SPECIES file not in NAMELIST format" write(*,*) "fixed format no longer supported" endif error stop endif if ((pohcconst.ne.-9.99).or.(pohdconst.ne.-9.99).or.(pohnconst.ne.-9.99)) then write(*,*) "ERROR: POHCCONST,POHDCONST, and POHNCONST in SPECIES file are deprecated." error stop endif species(pos_spec)=pspecies decay(pos_spec)=pdecay weta_gas(pos_spec)=pweta_gas wetb_gas(pos_spec)=pwetb_gas crain_aero(pos_spec)=pcrain_aero csnow_aero(pos_spec)=pcsnow_aero ccn_aero(pos_spec)=pccn_aero in_aero(pos_spec)=pin_aero reldiff(pos_spec)=preldiff henry(pos_spec)=phenry f0(pos_spec)=pf0 density(pos_spec)=pdensity if (pdia.ne.-9.9) then dquer(pos_spec)=pdia else if (pdquer.ne.-9.9) then write(*,*) 'WARNING: PDQUER will be depricated, please use PDIA instead.' dquer(pos_spec)=pdquer ! For backwards compatibility else dquer(pos_spec)=0.0 endif dsigma(pos_spec)=pdsigma ! ndia(pos_spec)=pndia dryvel(pos_spec)=pdryvel weightmolar(pos_spec)=pweightmolar emis_path(pos_spec)=pemis_path emis_file(pos_spec)=pemis_file emis_name(pos_spec)=pemis_name emis_unit(pos_spec)=pemis_unit emis_coeff(pos_spec)=pemis_coeff ishape(pos_spec)=pshape orient(pos_spec)=porient ! Daria Tatsii 2023: compute particle shape dimensions if (ishape(pos_spec).ge.1) then ! Compute shape according to given axes select case (ishape(pos_spec)) case (1) write(*,*) "Particle shape USER-DEFINED for particle", id_spec if ((psa.le.0.0).or.(pia.le.0.0).or.(pla.le.0.0)) then write(*,*) "#### ERROR: Shape=1 (user-defined) is chosen, & &but no valid axes are provided." write(*,*) "#### SPECIES file requires SA, IA, and LA parameter & &greater than zero." error stop endif write(*,*) "SA,IA,LA:",psa,pia,pla case (2) ! Cylinders (fibers) ! if (paspectratio.le.0.0) then write(*,*) "#### ERROR: Shape=2 cylinder is chosen, but no valid apect ratio is provided." write(*,*) "#### SPECIES file requires ASPECTRATIO parameter greater than zero." error stop endif psa=(((dquer(pos_spec)**3.0)*2.0)/ & (3.0*paspectratio))**(1.0/3.0) pia=psa pla=psa*paspectratio write(*,*) "Particle shape CYLINDER for particle", id_spec write(*,*) "SA,IA,LA:",psa,pia,pla case (3) ! Cubes ! write(*,*) "Particle shape CUBE for particle", id_spec psa=((dquer(pos_spec)**3)*pi/6.0)**(1.0/3.0) pia=(2.0**0.5)*psa pla=(3.0**0.5)*psa if ((psa.le.0.0).or.(pia.le.0.0).or.(pla.le.0.0)) then write(*,*) "#### ERROR: Shape=3 (user-defined) is chosen, but no valid axes are provided." write(*,*) "#### SPECIES file requires SA, IA, and LA parameter greater than zero." error stop endif write(*,*) "SA,IA,LA:",psa,pia,pla case (4) ! Tetrahedrons ! write(*,*) "Particle shape TETRAHEDRON for particle", id_spec pla=((dquer(pos_spec)**3)*pi*2**(0.5))**(1.0/3.0) pia=pla*((3.0/4.0)**(0.5)) psa=pla*((2.0/3.0)**(0.5)) if ((psa.le.0.0).or.(pia.le.0.0).or.(pla.le.0.0)) then write(*,*) "#### ERROR: Shape=4 (user-defined) is chosen, but no valid axes are provided." write(*,*) "#### SPECIES file requires SA, IA, and LA parameter greater than zero." error stop endif write(*,*) "SA,IA,LA:",psa,pia,pla case (5) ! Octahedrons ! write(*,*) "Particle shape OCTAHEDRON for particle", id_spec psa=dquer(pos_spec)*(pi/(2.0*2.0**(0.5)))**3 pia=psa pla=psa*(2.0**(0.5)) if ((psa.le.0.0).or.(pia.le.0.0).or.(pla.le.0.0)) then write(*,*) "#### ERROR: Shape=5 (user-defined) is chosen, but no valid axes are provided." write(*,*) "#### SPECIES file requires SA, IA, and LA parameter greater than zero." error stop endif write(*,*) "SA,IA,LA:",psa,pia,pla case (6) ! Ellipsoids ! write(*,*) "Particle shape ELLIPSOID for particle", id_spec psa=dquer(pos_spec)/(2.0**(1.0/3.0)) pia=psa pla=2*pia if ((psa.le.0.0).or.(pia.le.0.0).or.(pla.le.0.0)) then write(*,*) "#### ERROR: Shape=6 (user-defined) is chosen, but no valid axes are provided." write(*,*) "#### SPECIES file requires SA, IA, and LA parameter greater than zero." error stop endif write(*,*) "SA,IA,LA:",psa,pia,pla end select ! When using the shape option, dquer is the sphere equivalent diameter f=psa/pia e=pia/pla ! Drag coefficient scheme by Bagheri & Bonadonna, 2016 to define settling velocities of other shapes (by D.Tatsii) if ((ishape(pos_spec).eq.2).or.((ishape(pos_spec).eq.1).and. & (pia.eq.psa).and.(pla.ge.20.0*pia))) then Fn(pos_spec)=f*f*e ! simplified equation, validated by experiments with fibers Fs(pos_spec)=f*e**(1.3) ! simplified equation, validated by experiments with fibers else Fn(pos_spec)=f*f*e*((dquer(pos_spec))**3)/(psa*pia*pla) ! Newton's regime Fs(pos_spec)=f*e**(1.3)*(dquer(pos_spec)**3/(psa*pia*pla)) ! Stokes' regime endif ! Pre-compute ks and kn values needed for horizontal and average orientation (B&B Figure 12 k_(s,max)) ks1(pos_spec)=(Fs(pos_spec)**(1./3.) + Fs(pos_spec)**(-1./3.))/2. ks2(pos_spec)=0.5*((Fs(pos_spec)**0.05)+(Fs(pos_spec)**(-0.36))) kn2(pos_spec)=10.**(alpha2*(-log10(Fn(pos_spec)))**beta2) else ! Spheres write(*,*) "Particle shape SPHERE for particle", id_spec endif ! assign chemical reaction rate constants to table !************************************************* print*, 'readspecies: preactions = ',preactions if (any(preactions.ne."")) then do j=1,maxreagent if ( preactions(j).eq."" ) cycle do i=1,maxreagent if (preactions(j).eq.reagents(i)) then reaccconst(i,pos_spec)=pcconst(j) reacdconst(i,pos_spec)=pdconst(j) reacnconst(i,pos_spec)=pnconst(j) exit endif end do if (i.gt.nreagent) then write(*,*) '#### FLEXPART MODEL ERROR ####' write(*,*) '#### REAGENT NOT FOUND FOR ####' write(*,*) '#### REACTION '//trim(preactions(j))//' ####' error stop endif end do endif ! Check reaction rates if (lroot) then write(*,*) 'Reaction rates for ',species(pos_spec),':' do j=1,nreagent if (reaccconst(j,pos_spec).lt.0.) reaccconst(j,pos_spec)=0. if (reacdconst(j,pos_spec).lt.0.) reacdconst(j,pos_spec)=0. if (reacnconst(j,pos_spec).lt.0.) reacnconst(j,pos_spec)=0. write(*,*) reagents(j),': C, D, N = ',reaccconst(j,pos_spec),reacdconst(j,pos_spec),reacnconst(j,pos_spec) end do endif do j=1,24 ! 24 hours, starting with 0-1 local time area_hour(pos_spec,j)=parea_hour(j) point_hour(pos_spec,j)=ppoint_hour(j) end do do j=1,7 ! 7 days of the week, starting with Monday area_dow(pos_spec,j)=parea_dow(j) point_dow(pos_spec,j)=ppoint_dow(j) end do i=pos_spec !NIK 16.02.2015 ! Check scavenging parameters given in SPECIES file if (lroot) then ! ZHG 2016.04.07 Start of changes write(*,*) ' ' if (dquer(pos_spec) .gt.0) write(*,'(a,i3,a,a,a)') ' SPECIES: ', & id_spec,' ', species(pos_spec),' (AEROSOL) ' if (dquer(pos_spec) .le.0) write(*,'(a,i3,a,a,a)') ' SPECIES: ', & id_spec,' ', species(pos_spec),' (GAS) ' ! Particles !********** if (dquer(pos_spec).gt.0) then if (ccn_aero(pos_spec) .gt. 0) then write(*,'(a,f5.2)') ' Particle CCN efficiency (CCNeff):', ccn_aero(pos_spec) else write(*,'(a)') ' Particle CCN efficiency (CCNeff): OFF' endif if (in_aero(pos_spec) .gt. 0) then write(*,'(a,f5.2)') ' Particle IN efficiency (INeff) :', in_aero(pos_spec) else write(*,'(a)') ' Particle IN efficiency (INeff) : OFF' endif if (crain_aero(pos_spec) .gt. 0) then write(*,'(a,f5.2)') ' Particle Rain efficiency (Crain) :', crain_aero(pos_spec) else write(*,'(a)') ' Particle Rain efficiency (Crain) : OFF' endif if (csnow_aero(pos_spec) .gt. 0) then write(*,'(a,f5.2)') ' Particle Snow efficiency (Csnow) :', csnow_aero(pos_spec) else write(*,'(a)') ' Particle Snow efficiency (Csnow) : OFF' end if if (density(pos_spec) .gt. 0) then write(*,'(a)') ' Dry deposition is turned : ON' if (reldiff(pos_spec).gt.0) then error stop 'density>0 (SPECIES is a particle) implies reldiff <=0 ' endif else if (reldiff(pos_spec).le.0) then error stop 'density<=0 (SPECIES is a gas) implies reldiff >0 ' endif write(*,'(a)') ' Dry deposition is (density<0) : OFF' end if if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. & &ccn_aero(pos_spec).gt.1.0 .or. in_aero(pos_spec).gt.1.0) then write(*,*) '*******************************************' write(*,*) ' WARNING: Particle Scavenging parameter likely out of range ' write(*,*) ' Likely range for Crain 0.0-10' write(*,*) ' Likely range for Csnow 0.0-10' write(*,*) ' Physical range for CCNeff 0.0-1' write(*,*) ' Physical range for INeff 0.0-1' write(*,*) '*******************************************' end if else ! Gas !**** if (weta_gas(pos_spec) .gt. 0 .and. wetb_gas(pos_spec).gt.0) then write(*,*) ' Wet removal for gases is turned: ON' write(*,*) ' Gas below-cloud scavenging parameter A ', & &weta_gas(pos_spec) write(*,'(a,f5.2)') ' Gas below-cloud scavenging parameter B ', & &wetb_gas(pos_spec) else write(*,*) ' Wet removal for gases is turned: OFF ' end if if (reldiff(i).gt.0.) then write(*,*) ' Dry deposition for gases is turned: ON ' else write(*,*) ' Dry deposition for gases is turned: OFF ' end if if (weta_gas(pos_spec).gt.0.) then !if wet deposition is turned on if (weta_gas(pos_spec).gt.1E-04 .or. weta_gas(pos_spec).lt.1E-09 .or. & &wetb_gas(pos_spec).gt.0.8 .or. wetb_gas(pos_spec).lt.0.4) then write(*,*) '*******************************************' write(*,*) ' WARNING: Gas below-cloud scavengig is out of likely range' write(*,*) ' Likely range for A is 1E-04 to 1E-08' write(*,*) ' Likely range for B is 0.60 to 0.80 ' write(*,*) '*******************************************' end if endif if (((weta_gas(pos_spec).gt.0).or.(wetb_gas(pos_spec).gt.0)).and.& &(henry(pos_spec).le.0)) then if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set endif end if end if !if (ndia(pos_spec).gt.maxndia) then ! maxndia=ndia(pos_spec) !endif ndia(pos_spec)=maxndia ! Setting all ndia to maxndia (par_mod.f90) ! if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception if (dquer(i).gt.0 .and. dsigma(i).le.1.) then !dsigma(i)=1.0001 ! avoid floating exception !write(*,*) '#### FLEXPART MODEL ERROR! ####' write(*,*) '#### FLEXPART MODEL WARNING ####' write(*,*) '#### in SPECIES_',aspecnumb, ' ####' write(*,*) '#### from v10.4 dsigma has to be larger than 1 ####' write(*,*) '#### to adapt older SPECIES files, ####' write(*,*) '#### if dsigma was < 1 ####' write(*,*) '#### use the reciprocal of the old dsigma ####' if (.not.debug_mode) then error stop else write(*,*) 'debug mode: continue' endif endif if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####' write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####' write(*,*) '#### PARTICLE AND GAS. ####' write(*,*) '#### SPECIES NUMBER',aspecnumb error stop endif 22 close(unitspecies) ! namelist output if requested if (nmlout.and.lroot) then open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000) write(unitspecies,nml=species_params) close(unitspecies) endif return 996 write(*,*) '#####################################################' write(*,*) '#### FLEXPART MODEL ERROR! #### ' write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS #### ' write(*,*) '#### CONSTANT IS SET ####' write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE! #### ' write(*,*) '#####################################################' error stop 998 write(*,*) '#####################################################' write(*,*) '#### FLEXPART MODEL ERROR! #### ' write(*,*) '#### THE SPECIES FILE FOR SPECIES ', id_spec write(*,*) '#### CANNOT BE FOUND: CREATE FILE' write(*,*) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb write(*,*) '#####################################################' error stop 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readspecies subroutine readpartoptions !***************************************************************************** ! * ! This routine reads the age classes to be used for the current model * ! run. * ! * ! Author: A. Stohl * ! 20 March 2000 * ! HSO, 1 July 2014 * ! Added optional namelist input * ! * !***************************************************************************** ! * ! Variables: * ! * ! Constants: * ! * !***************************************************************************** implicit none integer :: np,stat ! namelist help variables integer :: ios logical :: & longitude=.false., & longitude_average=.false., & latitude=.false., & latitude_average=.false., & height=.false., & height_average=.false., & pv=.false., & pv_average=.false., & qv=.false., & qv_average=.false., & density=.false., & density_average=.false., & temperature=.false., & temperature_average=.false., & pressure=.false., & pressure_average=.false., & mixingheight=.false., & mixingheight_average=.false., & tropopause=.false., & tropopause_average=.false., & topography=.false., & topography_average=.false., & mass=.false., & mass_average=.false., & u=.false., & u_average=.false., & v=.false., & v_average=.false., & w=.false., & w_average=.false., & vsettling=.false., & vsettling_average=.false., & wetdeposition=.false., & drydeposition=.false. ! namelist declaration namelist /partoptions/ & longitude, & longitude_average, & latitude, & latitude_average, & height, & height_average, & pv, & pv_average, & qv, & qv_average, & density, & density_average, & temperature, & temperature_average, & pressure, & pressure_average, & mixingheight, & mixingheight_average, & tropopause, & tropopause_average, & topography, & topography_average, & mass, & mass_average, & u, & u_average, & v, & v_average, & w, & w_average, & vsettling, & vsettling_average, & wetdeposition, & drydeposition ! If age spectra claculation is switched on, ! open the AGECLASSSES file and read user options !************************************************ open(unitpartoptions,file=path(1)(1:length(1))//'PARTOPTIONS',form='formatted',status='old',err=9999) ! try to read in as a namelist read(unitpartoptions,partoptions,iostat=ios) close(unitpartoptions) if (ios.ne.0) then write(*,*) 'Namelist error in PARTOPTIONS file', trim(path(1)(1:length(1))//'PARTOPTIONS') error stop endif allocate( partopt(num_partopt),stat=stat) if (stat.ne.0) error stop "Could not allocate partopt" ! Save values in particle options derived type !********************************************* partopt(1)%long_name='longitude' partopt(1)%short_name='lon' partopt(1)%name='LO' partopt(1)%print=longitude partopt(2)%long_name='longitude_average' partopt(2)%short_name='lon_av' partopt(2)%name='lo' partopt(2)%print=longitude_average partopt(2)%average=.true. partopt(3)%long_name='latitude' partopt(3)%short_name='lat' partopt(3)%name='LA' partopt(3)%print=latitude partopt(4)%long_name='latitude_average' partopt(4)%short_name='lat_av' partopt(4)%name='la' partopt(4)%print=latitude_average partopt(4)%average=.true. partopt(5)%long_name='height' partopt(5)%short_name='z' partopt(5)%name='ZZ' partopt(5)%print=height partopt(6)%long_name='height_average' partopt(6)%short_name='z_av' partopt(6)%name='zz' partopt(6)%print=height_average partopt(6)%average=.true. partopt(7)%long_name='potential_vorticity' partopt(7)%short_name='pv' partopt(7)%name='PV' partopt(7)%print=pv partopt(8)%long_name='potential_vorticity_average' partopt(8)%short_name='pv_av' partopt(8)%name='pv' partopt(8)%print=pv_average partopt(8)%average=.true. partopt(9)%long_name='specific_humidity' partopt(9)%short_name='sh' partopt(9)%name='QV' partopt(9)%print=qv partopt(10)%long_name='specific_humidity_average' partopt(10)%short_name='sh_av' partopt(10)%name='qv' partopt(10)%print=qv_average partopt(10)%average=.true. partopt(11)%long_name='density' partopt(11)%short_name='rho' partopt(11)%name='RH' partopt(11)%print=density partopt(12)%long_name='density_average' partopt(12)%short_name='rho_av' partopt(12)%name='rh' partopt(12)%print=density_average partopt(12)%average=.true. partopt(13)%long_name='temperature' partopt(13)%short_name='T' partopt(13)%name='TT' partopt(13)%print=temperature partopt(14)%long_name='temperature_average' partopt(14)%short_name='T_av' partopt(14)%name='tt' partopt(14)%print=temperature_average partopt(14)%average=.true. partopt(15)%long_name='pressure' partopt(15)%short_name='prs' partopt(15)%name='PR' partopt(15)%print=pressure partopt(16)%long_name='pressure_average' partopt(16)%short_name='prs_av' partopt(16)%name='pr' partopt(16)%print=pressure_average partopt(16)%average=.true. partopt(17)%long_name='mixingheight' partopt(17)%short_name='hmix' partopt(17)%name='HM' partopt(17)%print=mixingheight partopt(18)%long_name='mixingheight_average' partopt(18)%short_name='hmix_av' partopt(18)%name='hm' partopt(18)%print=mixingheight_average partopt(18)%average=.true. partopt(19)%long_name='tropopause' partopt(19)%short_name='tro' partopt(19)%name='TR' partopt(19)%print=tropopause partopt(20)%long_name='tropopause_average' partopt(20)%short_name='tro_av' partopt(20)%name='tr' partopt(20)%print=tropopause_average partopt(20)%average=.true. partopt(21)%long_name='topography' partopt(21)%short_name='to' partopt(21)%name='TO' partopt(21)%print=topography partopt(22)%long_name='topography_average' partopt(22)%short_name='to_av' partopt(22)%name='to' partopt(22)%print=topography_average partopt(22)%average=.true. partopt(23)%long_name='mass' partopt(23)%short_name='m' partopt(23)%name='MA' partopt(23)%print=mass partopt(24)%long_name='mass_average' partopt(24)%short_name='m_av' partopt(24)%name='ma' partopt(24)%print=mass_average partopt(24)%average=.true. partopt(25)%long_name='longitudinal_velocity' partopt(25)%short_name='u' partopt(25)%name='UU' partopt(25)%print=u partopt(26)%long_name='longitudinal_velocity_average' partopt(26)%short_name='u_av' partopt(26)%name='uu' partopt(26)%print=u_average partopt(26)%average=.true. partopt(27)%long_name='latitudinal_velocity' partopt(27)%short_name='v' partopt(27)%name='VV' partopt(27)%print=v partopt(28)%long_name='latitudinal_velocity_average' partopt(28)%short_name='v_av' partopt(28)%name='vv' partopt(28)%print=v_average partopt(28)%average=.true. partopt(29)%long_name='vertical_velocity' partopt(29)%short_name='w' partopt(29)%name='WW' partopt(29)%print=w partopt(30)%long_name='vertical_velocity_average' partopt(30)%short_name='w_av' partopt(30)%name='ww' partopt(30)%print=w_average partopt(30)%average=.true. partopt(31)%long_name='settling_velocity' partopt(31)%short_name='vset' partopt(31)%name='VS' partopt(31)%print=vsettling partopt(32)%long_name='settling_velocity_average' partopt(32)%short_name='vset_av' partopt(32)%name='vs' partopt(32)%print=vsettling_average partopt(32)%average=.true. partopt(33)%long_name='wetdeposition' partopt(33)%short_name='wetdepo' partopt(33)%name='WD' partopt(33)%print=wetdeposition partopt(34)%long_name='drydeposition' partopt(34)%short_name='drydepo' partopt(34)%name='DD' partopt(34)%print=drydeposition ! Numbers are assigned to the averaged fields for proper ! allocation and reading in particle_mod and output_mod !****************************************************** n_average=0 do np=1,num_partopt if (partopt(np)%average .and. partopt(np)%print) then n_average=n_average+1 partopt(np)%i_average = n_average if ((partopt(np)%name.eq.'MA') .or. (partopt(np)%name.eq.'ma')) then n_average=n_average + (maxspec-1) endif endif end do ! write partoptions file in namelist format to output directory if requested if (nmlout.and.lroot) then open(unitpartoptions,file=path(2)(1:length(2))//'PARTOPTIONS.namelist',err=10000) write(unitpartoptions,nml=partoptions) close(unitpartoptions) endif ! Restart files, when using in combination with averaged particle output, ! need to be synchronised to prevent false averages in the first step of ! the new run !************************************************************************ if ((ipout.ne.0).and.(n_average.gt.0).and.(loutrestart.ne.-1)) then if (mod(loutrestart,ipoutfac*loutstep).ne.0) then write(*,*) '### FLEXPART MODEL ERROR! FILE COMMAND: ###' write(*,*) '### LOUTRESTART NEEDS TO BE DIVISABLE BY ###' write(*,*) '### LOUTSTEP*IPOUTFAC. ###' error stop endif endif return 9999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "PARTOPTIONS" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(1)(1:length(1)) error stop 10000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "PARTOPTIONS" #### ' write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' write(*,'(a)') path(2)(1:length(2)) error stop end subroutine readpartoptions subroutine skplin(nlines,iunit) ! i i !***************************************************************************** ! * ! This routine reads nlines from unit iunit and discards them * ! * ! Authors: Petra Seibert * ! * ! 31 Dec 1998 * ! * !***************************************************************************** ! * ! Variables: * ! * ! iunit unit number from which lines are to be skipped * ! nlines number of lines to be skipped * ! * !***************************************************************************** implicit none integer :: i,iunit, nlines do i=1,nlines read(iunit,*) end do end subroutine skplin end module readoptions_mod