CTBTO FLEXPART WO4 (2015-10-15)
 All Classes Files Functions Variables
readreleases.f90
Go to the documentation of this file.
1 !**********************************************************************
2 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5 ! *
6 ! This file is part of FLEXPART. *
7 ! *
8 ! FLEXPART is free software: you can redistribute it and/or modify *
9 ! it under the terms of the GNU General Public License as published by*
10 ! the Free Software Foundation, either version 3 of the License, or *
11 ! (at your option) any later version. *
12 ! *
13 ! FLEXPART is distributed in the hope that it will be useful, *
14 ! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16 ! GNU General Public License for more details. *
17 ! *
18 ! You should have received a copy of the GNU General Public License *
19 ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20 !**********************************************************************
21 
22 subroutine readreleases
23 
24  !*****************************************************************************
25  ! *
26  ! This routine reads the release point specifications for the current *
27  ! model run. Several release points can be used at the same time. *
28  ! *
29  ! Author: A. Stohl *
30  ! *
31  ! 18 May 1996 *
32  ! *
33  ! Update: 29 January 2001 *
34  ! Release altitude can be either in magl or masl *
35  ! *
36  !*****************************************************************************
37  ! *
38  ! Variables: *
39  ! decay decay constant of species *
40  ! dquer [um] mean particle diameters *
41  ! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
42  ! are between 0.1*dquer and 10*dquer *
43  ! ireleasestart, ireleaseend [s] starting time and ending time of each *
44  ! release *
45  ! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
46  ! is in hPa *
47  ! npart number of particles to be released *
48  ! nspec number of species to be released *
49  ! density [kg/m3] density of the particles *
50  ! rm [s/m] Mesophyll resistance *
51  ! species name of species *
52  ! xmass total mass of each species *
53  ! xpoint1,ypoint1 geograf. coordinates of lower left corner of release *
54  ! area *
55  ! xpoint2,ypoint2 geograf. coordinates of upper right corner of release *
56  ! area *
57  ! weta, wetb parameters to determine the wet scavenging coefficient *
58  ! zpoint1,zpoint2 height range, over which release takes place *
59  ! *
60  !*****************************************************************************
61 
62  use point_mod
63  use xmass_mod
64  use par_mod
65  use com_mod
66 
67  implicit none
68 
69  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
70  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
71  real(kind=dp) :: jul1,jul2,juldate
72  character(len=50) :: line
73  logical :: old
74 
75  !sec, read release to find how many releasepoints should be allocated
76 
77  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
78  err=999)
79 
80  ! Check the format of the RELEASES file (either in free format,
81  ! or using a formatted mask)
82  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
83  !**************************************************************************
84 
85  call skplin(12,unitreleases)
86  read (unitreleases,901) line
87 901 format (a)
88  if (index(line,'Total') .eq. 0) then
89  old = .false.
90  else
91  old = .true.
92  endif
93  rewind(unitreleases)
94 
95 
96  ! Skip first 11 lines (file header)
97  !**********************************
98 
99  call skplin(11,unitreleases)
100 
101 
102  read(unitreleases,*,err=998) nspec
103  if (old) call skplin(2,unitreleases)
104  do i=1,nspec
105  read(unitreleases,*,err=998) specnum_rel
106  if (old) call skplin(2,unitreleases)
107  end do
108 
109  numpoint=0
110 100 numpoint=numpoint+1
111  read(unitreleases,*,end=25)
112  read(unitreleases,*,err=998,end=25) idum,idum
113  if (old) call skplin(2,unitreleases)
114  read(unitreleases,*,err=998) idum,idum
115  if (old) call skplin(2,unitreleases)
116  read(unitreleases,*,err=998) xdum
117  if (old) call skplin(2,unitreleases)
118  read(unitreleases,*,err=998) xdum
119  if (old) call skplin(2,unitreleases)
120  read(unitreleases,*,err=998) xdum
121  if (old) call skplin(2,unitreleases)
122  read(unitreleases,*,err=998) xdum
123  if (old) call skplin(2,unitreleases)
124  read(unitreleases,*,err=998) idum
125  if (old) call skplin(2,unitreleases)
126  read(unitreleases,*,err=998) xdum
127  if (old) call skplin(2,unitreleases)
128  read(unitreleases,*,err=998) xdum
129  if (old) call skplin(2,unitreleases)
130  read(unitreleases,*,err=998) idum
131  if (old) call skplin(2,unitreleases)
132  do i=1,nspec
133  read(unitreleases,*,err=998) xdum
134  if (old) call skplin(2,unitreleases)
135  end do
136  !save compoint only for the first 1000 release points
137  read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
138  if (old) call skplin(1,unitreleases)
139 
140  goto 100
141 
142 25 numpoint=numpoint-1
143 
144  !allocate memory for numpoint releaspoint
145  allocate(ireleasestart(numpoint) &
146  ,stat=stat)
147  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
148  allocate(ireleaseend(numpoint) &
149  ,stat=stat)
150  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
151  allocate(xpoint1(numpoint) &
152  ,stat=stat)
153  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
154  allocate(xpoint2(numpoint) &
155  ,stat=stat)
156  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
157  allocate(ypoint1(numpoint) &
158  ,stat=stat)
159  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
160  allocate(ypoint2(numpoint) &
161  ,stat=stat)
162  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
163  allocate(zpoint1(numpoint) &
164  ,stat=stat)
165  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
166  allocate(zpoint2(numpoint) &
167  ,stat=stat)
168  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
169  allocate(kindz(numpoint) &
170  ,stat=stat)
171  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
172  allocate(xmass(numpoint,maxspec) &
173  ,stat=stat)
174  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
175  allocate(rho_rel(numpoint) &
176  ,stat=stat)
177  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
178  allocate(npart(numpoint) &
179  ,stat=stat)
180  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
181  allocate(xmasssave(numpoint) &
182  ,stat=stat)
183  if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
184 
185  write (*,*) ' Releasepoints allocated: ', numpoint
186 
187  do i=1,numpoint
188  xmasssave(i)=0.
189  end do
190 
191  !now save the information
192  dep=.false.
193  drydep=.false.
194  wetdep=.false.
195  ohrea=.false.
196  do i=1,maxspec
197  drydepspec(i)=.false.
198  end do
199 
200  rewind(unitreleases)
201 
202 
203  ! Skip first 11 lines (file header)
204  !**********************************
205 
206  call skplin(11,unitreleases)
207 
208 
209  ! Assign species-specific parameters needed for physical processes
210  !*************************************************************************
211 
212  read(unitreleases,*,err=998) nspec
213  if (nspec.gt.maxspec) goto 994
214  if (old) call skplin(2,unitreleases)
215  do i=1,nspec
216  read(unitreleases,*,err=998) specnum_rel
217  if (old) call skplin(2,unitreleases)
218  call readspecies(specnum_rel,i)
219 
220  ! For backward runs, only 1 species is allowed
221  !*********************************************
222 
223  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
224  !write(*,*) '#####################################################'
225  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
226  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
227  !write(*,*) '#####################################################'
228  ! stop
229  !endif
230 
231  ! Molecular weight
232  !*****************
233 
234  if (((iout.eq.2).or.(iout.eq.3)).and. &
235  (weightmolar(i).lt.0.)) then
236  write(*,*) 'For mixing ratio output, valid molar weight'
237  write(*,*) 'must be specified for all simulated species.'
238  write(*,*) 'Check table SPECIES or choose concentration'
239  write(*,*) 'output instead if molar weight is not known.'
240  stop
241  endif
242 
243 
244  ! Radioactive decay
245  !******************
246 
247  decay(i)=0.693147/decay(i) !conversion half life to decay constant
248 
249 
250  ! Dry deposition of gases
251  !************************
252 
253  if (reldiff(i).gt.0.) &
254  rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance
255 
256  ! Dry deposition of particles
257  !****************************
258 
259  vsetaver(i)=0.
260  cunningham(i)=0.
261  dquer(i)=dquer(i)*1000000. ! Conversion m to um
262  if (density(i).gt.0.) then ! Additional parameters
263  call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
264  do j=1,ni
265  fract(i,j)=fracth(j)
266  schmi(i,j)=schmih(j)
267  vset(i,j)=vsh(j)
268  cunningham(i)=cunningham(i)+cun*fract(i,j)
269  vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
270  end do
271  write(*,*) 'Average setting velocity: ',i,vsetaver(i)
272  endif
273 
274  ! Dry deposition for constant deposition velocity
275  !************************************************
276 
277  dryvel(i)=dryvel(i)*0.01 ! conversion to m/s
278 
279  ! Check if wet deposition or OH reaction shall be calculated
280  !***********************************************************
281  if (weta(i).gt.0.) then
282  wetdep=.true.
283  write (*,*) 'Wetdeposition switched on: ',weta(i),i
284  endif
285  if (ohreact(i).gt.0) then
286  ohrea=.true.
287  write (*,*) 'OHreaction switched on: ',ohreact(i),i
288  endif
289 
290 
291  if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
292  (dryvel(i).gt.0.)) then
293  drydep=.true.
294  drydepspec(i)=.true.
295  endif
296 
297  end do
298 
299  if (wetdep.or.drydep) dep=.true.
300 
301  ! Read specifications for each release point
302  !*******************************************
303 
304  numpoint=0
305  numpartmax=0
306  releaserate=0.
307 1000 numpoint=numpoint+1
308  read(unitreleases,*,end=250)
309  read(unitreleases,*,err=998,end=250) id1,it1
310  if (old) call skplin(2,unitreleases)
311  read(unitreleases,*,err=998) id2,it2
312  if (old) call skplin(2,unitreleases)
313  read(unitreleases,*,err=998) xpoint1(numpoint)
314  if (old) call skplin(2,unitreleases)
315  read(unitreleases,*,err=998) ypoint1(numpoint)
316  if (old) call skplin(2,unitreleases)
317  read(unitreleases,*,err=998) xpoint2(numpoint)
318  if (old) call skplin(2,unitreleases)
319  read(unitreleases,*,err=998) ypoint2(numpoint)
320  if (old) call skplin(2,unitreleases)
321  read(unitreleases,*,err=998) kindz(numpoint)
322  if (old) call skplin(2,unitreleases)
323  read(unitreleases,*,err=998) zpoint1(numpoint)
324  if (old) call skplin(2,unitreleases)
325  read(unitreleases,*,err=998) zpoint2(numpoint)
326  if (old) call skplin(2,unitreleases)
327  read(unitreleases,*,err=998) npart(numpoint)
328  if (old) call skplin(2,unitreleases)
329  do i=1,nspec
330  read(unitreleases,*,err=998) xmass(numpoint,i)
331  if (old) call skplin(2,unitreleases)
332  end do
333  !save compoint only for the first 1000 release points
334  if (numpoint.le.1000) then
335  read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
336  else
337  read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
338  endif
339  if (old) call skplin(1,unitreleases)
340  if (numpoint.le.1000) then
341  if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
342  (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
343  (compoint(numpoint)(1:8).eq.' ')) goto 250
344  else
345  if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
346  (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
347  endif
348 
349 
350  ! If a release point contains no particles, stop and issue error message
351  !***********************************************************************
352 
353  if (npart(numpoint).eq.0) then
354  write(*,*) 'FLEXPART MODEL ERROR'
355  write(*,*) 'RELEASES file is corrupt.'
356  write(*,*) 'At least for one release point, there are zero'
357  write(*,*) 'particles released. Make changes to RELEASES.'
358  stop
359  endif
360 
361  ! Check whether x coordinates of release point are within model domain
362  !*********************************************************************
363 
364  if (xpoint1(numpoint).lt.xlon0) &
365  xpoint1(numpoint)=xpoint1(numpoint)+360.
366  if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
367  xpoint1(numpoint)=xpoint1(numpoint)-360.
368  if (xpoint2(numpoint).lt.xlon0) &
369  xpoint2(numpoint)=xpoint2(numpoint)+360.
370  if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
371  xpoint2(numpoint)=xpoint2(numpoint)-360.
372 
373  ! Determine relative beginning and ending times of particle release
374  !******************************************************************
375 
376  jul1=juldate(id1,it1)
377  jul2=juldate(id2,it2)
378  if (jul1.gt.jul2) then
379  write(*,*) 'FLEXPART MODEL ERROR'
380  write(*,*) 'Release stops before it begins.'
381  write(*,*) 'Make changes to file RELEASES.'
382  stop
383  endif
384  if (mdomainfill.eq.0) then ! no domain filling
385  if (ldirect.eq.1) then
386  if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
387  write(*,*) 'FLEXPART MODEL ERROR'
388  write(*,*) 'Release starts before simulation begins or ends'
389  write(*,*) 'after simulation stops.'
390  write(*,*) 'Make files COMMAND and RELEASES consistent.'
391  stop
392  endif
393  ireleasestart(numpoint)=int((jul1-bdate)*86400.)
394  ireleaseend(numpoint)=int((jul2-bdate)*86400.)
395  else if (ldirect.eq.-1) then
396  if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
397  write(*,*) 'FLEXPART MODEL ERROR'
398  write(*,*) 'Release starts before simulation begins or ends'
399  write(*,*) 'after simulation stops.'
400  write(*,*) 'Make files COMMAND and RELEASES consistent.'
401  stop
402  endif
403  ireleasestart(numpoint)=int((jul1-bdate)*86400.)
404  ireleaseend(numpoint)=int((jul2-bdate)*86400.)
405  endif
406  endif
407 
408 
409  ! Check, whether the total number of particles may exceed totally allowed
410  ! number of particles at some time during the simulation
411  !************************************************************************
412 
413  ! Determine the release rate (particles per second) and total number
414  ! of particles released during the simulation
415  !*******************************************************************
416 
417  if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
418  releaserate=releaserate+real(npart(numpoint))/ &
419  real(ireleaseend(numpoint)-ireleasestart(numpoint))
420  else
421  releaserate=99999999
422  endif
423  numpartmax=numpartmax+npart(numpoint)
424  goto 1000
425 
426 
427 250 close(unitreleases)
428 
429  write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ', numpartmax
430  numpoint=numpoint-1
431 
432  if (ioutputforeachrelease.eq.1) then
433  maxpointspec_act=numpoint
434  else
435  maxpointspec_act=1
436  endif
437 
438  if (releaserate.gt. &
439  0.99*real(maxpart)/real(lage(nageclass))) then
440  if (numpartmax.gt.maxpart) then
441  write(*,*) '#####################################################'
442  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
443  write(*,*) '#### ####'
444  write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
445  write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
446  write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
447  write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
448  write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
449  write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
450  write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
451  write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
452  write(*,*) '#####################################################'
453  write(*,*) 'Maximum release rate may be: ',releaserate, &
454  ' particles per second'
455  write(*,*) 'Maximum allowed release rate is: ', &
456  real(maxpart)/real(lage(nageclass)),' particles per second'
457  write(*,*) &
458  'Total number of particles released during the simulation is: ', &
459  numpartmax
460  write(*,*) 'Maximum allowed number of particles is: ',maxpart
461  endif
462  endif
463 
464  return
465 
466 994 write(*,*) '#####################################################'
467  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
468  write(*,*) '#### ####'
469  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
470  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
471  write(*,*) '#####################################################'
472  stop
473 
474 998 write(*,*) '#####################################################'
475  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
476  write(*,*) '#### ####'
477  write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS ####'
478  write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
479  write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
480  write(*,*) '#### FILE ... ####'
481  write(*,*) '#####################################################'
482  stop
483 
484 
485 999 write(*,*) '#####################################################'
486  write(*,*) ' FLEXPART MODEL SUBROUTINE READRELEASES: '
487  write(*,*)
488  write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
489  write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
490  write(*,*) 'PERMITTED FOR ANY ACCESS'
491  write(*,*) '#####################################################'
492  stop
493 
494 end subroutine readreleases
xmass_mod
Definition: xmass_mod.f90:22
skplin
subroutine skplin(nlines, iunit)
Definition: skplin.f90:22
readspecies
subroutine readspecies(id_spec, pos_spec)
Definition: readspecies.f90:22
par_mod
Definition: par_mod.F90:34
readreleases
subroutine readreleases
Definition: readreleases.f90:22
juldate
real(kind=dp) function juldate(yyyymmdd, hhmiss)
Definition: juldate.f90:22
point_mod
Definition: point_mod.f90:22
com_mod
Definition: com_mod.f90:13
part0
subroutine part0(dquer, dsigma, density, fract, schmi, cun, vsh)
Definition: part0.f90:22
Generated on Thu Oct 15 2015 06:20:58 for CTBTO FLEXPART WO4 (2015-10-15) by   doxygen 1.8.6