Select Git revision
readreleases.f90 20.76 KiB
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine 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, wetb parameters to determine the wet scavenging coefficient *
! zpoint1,zpoint2 height range, over which release takes place *
! *
!*****************************************************************************
use point_mod
use xmass_mod
use par_mod
use com_mod
implicit none
integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
integer :: specnum_rel(maxspec)
real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
real(kind=dp) :: jul1,jul2,juldate
character(len=50) :: line
logical :: old
! help variables for namelist reading
integer :: numpoints, parts, readerror
integer*2 :: zkind
integer :: idate1, itime1, idate2, itime2
real :: lon1,lon2,lat1,lat2,z1,z2 !,mass(nspec)
real :: mass(45)
character*40 :: comment
! 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
! presetting namelist releases_ctrl
nspec = -1 ! use negative value to determine failed namelist input
specnum_rel(1) = 0
!sec, read release to find how many releasepoints should be allocated
open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
err=999)
! check if namelist input provided
read(unitreleases,releases_ctrl,iostat=readerror)
if ((readerror.ne.0).or.(nspec.lt.0)) then
! no namelist format, close file and allow reopening in old format
rewind(unitreleases)
readerror=1 ! indicates old format
! Check the format of the RELEASES file (either in free format,
! or using a formatted mask)
! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
!**************************************************************************
call skplin(12,unitreleases)
read (unitreleases,901) line
901 format (a)
if (index(line,'Total') .eq. 0) then
old = .false.
else
old = .true.
endif
rewind(unitreleases)
! Skip first 11 lines (file header)
!**********************************
call skplin(11,unitreleases)
read(unitreleases,*,err=998) nspec
if (old) call skplin(2,unitreleases)
do i=1,nspec
read(unitreleases,*,err=998) specnum_rel(1)
if (old) call skplin(2,unitreleases)
end do
100 numpoint=numpoint+1
read(unitreleases,*,end=25)
read(unitreleases,*,err=998,end=25) idum,idum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) idum,idum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) idum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) idum
if (old) call skplin(2,unitreleases)
do i=1,nspec
read(unitreleases,*,err=998) xdum
if (old) call skplin(2,unitreleases)
end do
!save compoint only for the first 1000 release points
read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
if (old) call skplin(1,unitreleases)
goto 100
25 numpoint=numpoint-1
else
readerror=0
do while (readerror.eq.0)
idate1=-1
read(unitreleases,release,iostat=readerror)
if ((idate1.lt.0).or.(readerror.ne.0)) then
readerror=1
else
numpoint=numpoint+1
endif
end do
readerror=0
endif ! if namelist input
rewind(unitreleases)
!allocate memory for numpoint releaspoint
allocate(ireleasestart(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(ireleaseend(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(xpoint1(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(xpoint2(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(ypoint1(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(ypoint2(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(zpoint1(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(zpoint2(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(kindz(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(xmass(numpoint,maxspec),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(rho_rel(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(npart(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
allocate(xmasssave(numpoint),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
write (*,*) numpoint,' releasepoints allocated.'
do i=1,numpoint
xmasssave(i)=0.
end do
!now save the information
DEP=.false.
DRYDEP=.false.
WETDEP=.false.
OHREA=.false.
do i=1,maxspec
DRYDEPSPEC(i)=.false.
end do
if (readerror.ne.0) then
! Skip first 11 lines (file header)
!**********************************
call skplin(11,unitreleases)
! Assign species-specific parameters needed for physical processes
!*************************************************************************
read(unitreleases,*,err=998) nspec
if (nspec.gt.maxspec) goto 994
if (old) call skplin(2,unitreleases)
endif
do i=1,nspec
if (readerror.ne.0) then
read(unitreleases,*,err=998) specnum_rel(1)
if (old) call skplin(2,unitreleases)
call readspecies(specnum_rel(1),i)
else
call readspecies(specnum_rel(i),i)
endif
! For backward runs, only 1 species is allowed
!*********************************************
!if ((ldirect.lt.0).and.(nspec.gt.1)) then
!write(*,*) '#####################################################'
!write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
!write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
!write(*,*) '#####################################################'
! stop
!endif
! 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.'
stop
endif
! Radioactive decay
!******************
decay(i)=0.693147/decay(i) !conversion half life to decay constant
! 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),fracth,schmih,cun,vsh)
do j=1,ni
fract(i,j)=fracth(j)
schmi(i,j)=schmih(j)
vset(i,j)=vsh(j)
cunningham(i)=cunningham(i)+cun*fract(i,j)
vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
end do
write(*,*) 'Average setting 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 or OH reaction shall be calculated
!***********************************************************
if (weta(i).gt.0.) then
WETDEP=.true.
write (*,*) 'Wetdeposition switched on: ',weta(i),i
endif
if (ohreact(i).gt.0) then
OHREA=.true.
write (*,*) 'OHreaction switched on: ',ohreact(i),i
endif
if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
(dryvel(i).gt.0.)) then
DRYDEP=.true.
DRYDEPSPEC(i)=.true.
endif
end do
if (WETDEP.or.DRYDEP) DEP=.true.
! Read specifications for each release point
!*******************************************
numpoints=numpoint
numpoint=0
numpartmax=0
releaserate=0.
1000 numpoint=numpoint+1
if (readerror.eq.0) then ! reading namelist format
if (numpoint.gt.numpoints) goto 250
zkind = 1
mass = 0
parts = 0
comment = ' '
read(unitreleases,release,iostat=readerror)
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
else
read(unitreleases,*,end=250)
read(unitreleases,*,err=998,end=250) id1,it1
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) id2,it2
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xpoint1(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) ypoint1(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) xpoint2(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) ypoint2(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) kindz(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) zpoint1(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) zpoint2(numpoint)
if (old) call skplin(2,unitreleases)
read(unitreleases,*,err=998) npart(numpoint)
if (old) call skplin(2,unitreleases)
do i=1,nspec
read(unitreleases,*,err=998) xmass(numpoint,i)
if (old) call skplin(2,unitreleases)
end do
!save compoint only for the first 1000 release points
if (numpoint.le.1000) then
read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
else
read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
endif
if (old) call skplin(1,unitreleases)
if (numpoint.le.1000) then
if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
(compoint(numpoint)(1:8).eq.' ')) goto 250
else
if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
endif
endif ! if namelist format
! 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.'
stop
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)
if (jul1.gt.jul2) then
write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release stops before it begins.'
write(*,*) 'Make changes to file RELEASES.'
stop
endif
if (mdomainfill.eq.0) then ! no domain filling
if (ldirect.eq.1) then
if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release starts before simulation begins or ends'
write(*,*) 'after simulation stops.'
write(*,*) 'Make files COMMAND and RELEASES consistent.'
stop
endif
ireleasestart(numpoint)=int((jul1-bdate)*86400.)
ireleaseend(numpoint)=int((jul2-bdate)*86400.)
else if (ldirect.eq.-1) then
if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release starts before simulation begins or ends'
write(*,*) 'after simulation stops.'
write(*,*) 'Make files COMMAND and RELEASES consistent.'
stop
endif
ireleasestart(numpoint)=int((jul1-bdate)*86400.)
ireleaseend(numpoint)=int((jul2-bdate)*86400.)
endif
endif
! Check, whether the total number of particles may exceed totally allowed
! number of particles at some time during the simulation
!************************************************************************
! 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 1000
250 close(unitreleases)
write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ', numpartmax
numpoint=numpoint-1
if (ioutputforeachrelease.eq.1) then
maxpointspec_act=numpoint
else
maxpointspec_act=1
endif
if (releaserate.gt. &
0.99*real(maxpart)/real(lage(nageclass))) then
if (numpartmax.gt.maxpart) 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
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(*,*) '#####################################################'
stop
998 write(*,*) '#####################################################'
write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
write(*,*) '#### ####'
write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS ####'
write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
write(*,*) '#### FILE ... ####'
write(*,*) '#####################################################'
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(*,*) '#####################################################'
stop
end subroutine readreleases