source: branches/flexpart91_hasod/src_parallel/readreleases.f90 @ 10

Last change on this file since 10 was 10, checked in by hasod, 7 years ago

ADD: namelist input implemented for all common input files

  • Property svn:executable set to *
File size: 20.7 KB
Line 
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
22subroutine 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  !     HSO, 12 August 2013
36  !     Added optional namelist input
37  !                                                                            *
38  !*****************************************************************************
39  !                                                                            *
40  ! Variables:                                                                 *
41  ! decay               decay constant of species                              *
42  ! dquer [um]          mean particle diameters                                *
43  ! dsigma              e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
44  !                     are between 0.1*dquer and 10*dquer                     *
45  ! ireleasestart, ireleaseend [s] starting time and ending time of each       *
46  !                     release                                                *
47  ! kindz               1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
48  !                     is in hPa                                              *
49  ! npart               number of particles to be released                     *
50  ! nspec               number of species to be released                       *
51  ! density [kg/m3]     density of the particles                               *
52  ! rm [s/m]            Mesophyll resistance                                   *
53  ! species             name of species                                        *
54  ! xmass               total mass of each species                             *
55  ! xpoint1,ypoint1     geograf. coordinates of lower left corner of release   *
56  !                     area                                                   *
57  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
58  !                     area                                                   *
59  ! weta, wetb          parameters to determine the wet scavenging coefficient *
60  ! zpoint1,zpoint2     height range, over which release takes place           *
61  !                                                                            *
62  !*****************************************************************************
63
64  use point_mod
65  use xmass_mod
66  use par_mod
67  use com_mod
68
69  implicit none
70
71  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
72  integer :: specnum_rel(maxspec)
73  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
74  real(kind=dp) :: jul1,jul2,juldate
75  character(len=50) :: line
76  logical :: old
77
78  ! help variables for namelist reading
79  integer :: numpoints, parts, readerror
80  integer*2 :: zkind
81  integer :: idate1, itime1, idate2, itime2
82  real :: lon1,lon2,lat1,lat2,z1,z2,mass(nspec)
83  character*40 :: comment
84
85  ! declare namelists
86  namelist /releases_ctrl/ &
87    nspec, &
88    specnum_rel
89
90  namelist /release/ &
91    idate1, itime1, &
92    idate2, itime2, &
93    lon1, lon2, &
94    lat1, lat2, &
95    z1, z2, &
96    zkind, &
97    mass, &
98    parts, &
99    comment
100
101  numpoint=0
102
103  ! presetting namelist releases_ctrl
104  nspec = -1  ! use negative value to determine failed namelist input
105  specnum_rel(1) = 0
106
107  !sec, read release to find how many releasepoints should be allocated
108  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
109       err=999)
110
111  ! check if namelist input provided
112  read(unitreleases,releases_ctrl,iostat=readerror)
113
114  if ((readerror.ne.0).or.(nspec.lt.0)) then
115
116  ! no namelist format, close file and allow reopening in old format
117  rewind(unitreleases)
118
119  readerror=1 ! indicates old format
120
121  ! Check the format of the RELEASES file (either in free format,
122  ! or using a formatted mask)
123  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
124  !**************************************************************************
125
126  call skplin(12,unitreleases)
127  read (unitreleases,901) line
128901   format (a)
129  if (index(line,'Total') .eq. 0) then
130    old = .false.
131  else
132    old = .true.
133  endif
134  rewind(unitreleases)
135
136  ! Skip first 11 lines (file header)
137  !**********************************
138
139  call skplin(11,unitreleases)
140  read(unitreleases,*,err=998) nspec
141  if (old) call skplin(2,unitreleases)
142  do i=1,nspec
143    read(unitreleases,*,err=998) specnum_rel(1)
144    if (old) call skplin(2,unitreleases)
145  end do
146
147100   numpoint=numpoint+1
148  read(unitreleases,*,end=25)
149  read(unitreleases,*,err=998,end=25) idum,idum
150  if (old) call skplin(2,unitreleases)
151  read(unitreleases,*,err=998) idum,idum
152  if (old) call skplin(2,unitreleases)
153  read(unitreleases,*,err=998) xdum
154  if (old) call skplin(2,unitreleases)
155  read(unitreleases,*,err=998) xdum
156  if (old) call skplin(2,unitreleases)
157  read(unitreleases,*,err=998) xdum
158  if (old) call skplin(2,unitreleases)
159  read(unitreleases,*,err=998) xdum
160  if (old) call skplin(2,unitreleases)
161  read(unitreleases,*,err=998) idum
162  if (old) call skplin(2,unitreleases)
163  read(unitreleases,*,err=998) xdum
164  if (old) call skplin(2,unitreleases)
165  read(unitreleases,*,err=998) xdum
166  if (old) call skplin(2,unitreleases)
167  read(unitreleases,*,err=998) idum
168  if (old) call skplin(2,unitreleases)
169  do i=1,nspec
170    read(unitreleases,*,err=998) xdum
171    if (old) call skplin(2,unitreleases)
172  end do
173  !save compoint only for the first 1000 release points
174  read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
175  if (old) call skplin(1,unitreleases)
176
177  goto 100
178
17925   numpoint=numpoint-1
180
181  else
182    readerror=0
183    do while (readerror.eq.0)
184      idate1=-1
185      read(unitreleases,release,iostat=readerror)
186      if ((idate1.lt.0).or.(readerror.ne.0)) then
187        readerror=1
188      else
189        numpoint=numpoint+1
190      endif
191    end do
192    readerror=0
193  endif ! if namelist input
194
195  rewind(unitreleases)
196
197  !allocate memory for numpoint releaspoint
198  allocate(ireleasestart(numpoint),stat=stat)
199       
200  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
201  allocate(ireleaseend(numpoint),stat=stat)
202  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
203  allocate(xpoint1(numpoint),stat=stat)
204  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
205  allocate(xpoint2(numpoint),stat=stat)
206  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
207  allocate(ypoint1(numpoint),stat=stat)
208  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
209  allocate(ypoint2(numpoint),stat=stat)
210  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
211  allocate(zpoint1(numpoint),stat=stat)
212  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
213  allocate(zpoint2(numpoint),stat=stat)
214  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
215  allocate(kindz(numpoint),stat=stat)
216  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
217  allocate(xmass(numpoint,maxspec),stat=stat)
218  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
219  allocate(rho_rel(numpoint),stat=stat)
220  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
221  allocate(npart(numpoint),stat=stat)
222  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
223  allocate(xmasssave(numpoint),stat=stat)
224  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
225
226  write (*,*) numpoint,' releasepoints allocated.'
227
228  do i=1,numpoint
229    xmasssave(i)=0.
230  end do
231
232  !now save the information
233  DEP=.false.
234  DRYDEP=.false.
235  WETDEP=.false.
236  OHREA=.false.
237  do i=1,maxspec
238    DRYDEPSPEC(i)=.false.
239  end do
240
241  if (readerror.ne.0) then
242  ! Skip first 11 lines (file header)
243  !**********************************
244
245    call skplin(11,unitreleases)
246
247  ! Assign species-specific parameters needed for physical processes
248  !*************************************************************************
249
250    read(unitreleases,*,err=998) nspec
251    if (nspec.gt.maxspec) goto 994
252  if (old) call skplin(2,unitreleases)
253  endif
254  do i=1,nspec
255    if (readerror.ne.0) then
256      read(unitreleases,*,err=998) specnum_rel(1)
257      if (old) call skplin(2,unitreleases)
258      call readspecies(specnum_rel(1),i)
259    else
260      call readspecies(specnum_rel(i),i)
261    endif
262
263  ! For backward runs, only 1 species is allowed
264  !*********************************************
265
266  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
267  !write(*,*) '#####################################################'
268  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
269  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
270  !write(*,*) '#####################################################'
271  !  stop
272  !endif
273
274  ! Molecular weight
275  !*****************
276
277    if (((iout.eq.2).or.(iout.eq.3)).and. &
278         (weightmolar(i).lt.0.)) then
279      write(*,*) 'For mixing ratio output, valid molar weight'
280      write(*,*) 'must be specified for all simulated species.'
281      write(*,*) 'Check table SPECIES or choose concentration'
282      write(*,*) 'output instead if molar weight is not known.'
283      stop
284    endif
285
286
287  ! Radioactive decay
288  !******************
289
290    decay(i)=0.693147/decay(i) !conversion half life to decay constant
291
292
293  ! Dry deposition of gases
294  !************************
295
296    if (reldiff(i).gt.0.) &
297         rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
298
299  ! Dry deposition of particles
300  !****************************
301
302    vsetaver(i)=0.
303    cunningham(i)=0.
304    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
305    if (density(i).gt.0.) then                  ! Additional parameters
306     call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
307      do j=1,ni
308        fract(i,j)=fracth(j)
309        schmi(i,j)=schmih(j)
310        vset(i,j)=vsh(j)
311        cunningham(i)=cunningham(i)+cun*fract(i,j)
312        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
313      end do
314      write(*,*) 'Average setting velocity: ',i,vsetaver(i)
315    endif
316
317  ! Dry deposition for constant deposition velocity
318  !************************************************
319
320    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
321
322  ! Check if wet deposition or OH reaction shall be calculated
323  !***********************************************************
324    if (weta(i).gt.0.)  then
325      WETDEP=.true.
326      write (*,*) 'Wetdeposition switched on: ',weta(i),i
327    endif
328    if (ohreact(i).gt.0) then
329      OHREA=.true.
330      write (*,*) 'OHreaction switched on: ',ohreact(i),i
331    endif
332
333
334    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
335         (dryvel(i).gt.0.)) then
336      DRYDEP=.true.
337      DRYDEPSPEC(i)=.true.
338    endif
339
340  end do
341
342    if (WETDEP.or.DRYDEP) DEP=.true.
343
344  ! Read specifications for each release point
345  !*******************************************
346  numpoints=numpoint
347  numpoint=0
348  numpartmax=0
349  releaserate=0.
3501000   numpoint=numpoint+1
351
352  if (readerror.eq.0) then ! reading namelist format
353
354  if (numpoint.gt.numpoints) goto 250
355  zkind = 1
356  mass = 0
357  parts = 0
358  comment = ' '
359  read(unitreleases,release,iostat=readerror)
360  id1=idate1
361  it1=itime1
362  id2=idate2
363  it2=itime2
364  xpoint1(numpoint)=lon1
365  xpoint2(numpoint)=lon2
366  ypoint1(numpoint)=lat1
367  ypoint2(numpoint)=lat2
368  zpoint1(numpoint)=z1
369  zpoint2(numpoint)=z2
370  kindz(numpoint)=zkind
371  do i=1,nspec
372    xmass(numpoint,i)=mass(i)
373  end do
374  npart(numpoint)=parts
375  compoint(min(1001,numpoint))=comment
376
377  else
378
379  read(unitreleases,*,end=250)
380  read(unitreleases,*,err=998,end=250) id1,it1
381  if (old) call skplin(2,unitreleases)
382  read(unitreleases,*,err=998) id2,it2
383  if (old) call skplin(2,unitreleases)
384  read(unitreleases,*,err=998) xpoint1(numpoint)
385  if (old) call skplin(2,unitreleases)
386  read(unitreleases,*,err=998) ypoint1(numpoint)
387  if (old) call skplin(2,unitreleases)
388  read(unitreleases,*,err=998) xpoint2(numpoint)
389  if (old) call skplin(2,unitreleases)
390  read(unitreleases,*,err=998) ypoint2(numpoint)
391  if (old) call skplin(2,unitreleases)
392  read(unitreleases,*,err=998) kindz(numpoint)
393  if (old) call skplin(2,unitreleases)
394  read(unitreleases,*,err=998) zpoint1(numpoint)
395  if (old) call skplin(2,unitreleases)
396  read(unitreleases,*,err=998) zpoint2(numpoint)
397  if (old) call skplin(2,unitreleases)
398  read(unitreleases,*,err=998) npart(numpoint)
399  if (old) call skplin(2,unitreleases)
400  do i=1,nspec
401    read(unitreleases,*,err=998) xmass(numpoint,i)
402    if (old) call skplin(2,unitreleases)
403  end do
404  !save compoint only for the first 1000 release points
405  if (numpoint.le.1000) then
406    read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
407  else
408    read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
409  endif
410  if (old) call skplin(1,unitreleases)
411  if (numpoint.le.1000) then
412    if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
413         (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
414         (compoint(numpoint)(1:8).eq.'        ')) goto 250
415  else
416    if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
417         (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
418  endif
419
420  endif ! if namelist format
421
422  ! If a release point contains no particles, stop and issue error message
423  !***********************************************************************
424
425  if (npart(numpoint).eq.0) then
426    write(*,*) 'FLEXPART MODEL ERROR'
427    write(*,*) 'RELEASES file is corrupt.'
428    write(*,*) 'At least for one release point, there are zero'
429    write(*,*) 'particles released. Make changes to RELEASES.'
430    stop
431  endif
432
433  ! Check whether x coordinates of release point are within model domain
434  !*********************************************************************
435
436   if (xpoint1(numpoint).lt.xlon0) &
437        xpoint1(numpoint)=xpoint1(numpoint)+360.
438   if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
439        xpoint1(numpoint)=xpoint1(numpoint)-360.
440   if (xpoint2(numpoint).lt.xlon0) &
441        xpoint2(numpoint)=xpoint2(numpoint)+360.
442   if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
443        xpoint2(numpoint)=xpoint2(numpoint)-360.
444
445  ! Determine relative beginning and ending times of particle release
446  !******************************************************************
447
448  jul1=juldate(id1,it1)
449  jul2=juldate(id2,it2)
450  if (jul1.gt.jul2) then
451    write(*,*) 'FLEXPART MODEL ERROR'
452    write(*,*) 'Release stops before it begins.'
453    write(*,*) 'Make changes to file RELEASES.'
454    stop
455  endif
456  if (mdomainfill.eq.0) then   ! no domain filling
457    if (ldirect.eq.1) then
458      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
459        write(*,*) 'FLEXPART MODEL ERROR'
460        write(*,*) 'Release starts before simulation begins or ends'
461        write(*,*) 'after simulation stops.'
462        write(*,*) 'Make files COMMAND and RELEASES consistent.'
463        stop
464      endif
465      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
466      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
467    else if (ldirect.eq.-1) then
468      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
469        write(*,*) 'FLEXPART MODEL ERROR'
470        write(*,*) 'Release starts before simulation begins or ends'
471        write(*,*) 'after simulation stops.'
472        write(*,*) 'Make files COMMAND and RELEASES consistent.'
473        stop
474      endif
475      ireleasestart(numpoint)=int((jul1-bdate)*86400.)
476      ireleaseend(numpoint)=int((jul2-bdate)*86400.)
477    endif
478  endif
479
480
481  ! Check, whether the total number of particles may exceed totally allowed
482  ! number of particles at some time during the simulation
483  !************************************************************************
484
485  ! Determine the release rate (particles per second) and total number
486  ! of particles released during the simulation
487  !*******************************************************************
488
489  if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
490    releaserate=releaserate+real(npart(numpoint))/ &
491         real(ireleaseend(numpoint)-ireleasestart(numpoint))
492  else
493    releaserate=99999999
494  endif
495  numpartmax=numpartmax+npart(numpoint)
496  goto 1000
497
498
499250   close(unitreleases)
500
501  write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
502  numpoint=numpoint-1
503
504  if (ioutputforeachrelease.eq.1) then
505    maxpointspec_act=numpoint
506  else
507    maxpointspec_act=1
508  endif
509
510  if (releaserate.gt. &
511       0.99*real(maxpart)/real(lage(nageclass))) then
512    if (numpartmax.gt.maxpart) then
513  write(*,*) '#####################################################'
514  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
515  write(*,*) '####                                             ####'
516  write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
517  write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
518  write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
519  write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
520  write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
521  write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
522  write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
523  write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
524  write(*,*) '#####################################################'
525      write(*,*) 'Maximum release rate may be: ',releaserate, &
526           ' particles per second'
527      write(*,*) 'Maximum allowed release rate is: ', &
528           real(maxpart)/real(lage(nageclass)),' particles per second'
529      write(*,*) &
530           'Total number of particles released during the simulation is: ', &
531           numpartmax
532      write(*,*) 'Maximum allowed number of particles is: ',maxpart
533    endif
534  endif
535
536  return
537
538994   write(*,*) '#####################################################'
539  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
540  write(*,*) '####                                             ####'
541  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
542  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
543  write(*,*) '#####################################################'
544  stop
545
546998   write(*,*) '#####################################################'
547  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
548  write(*,*) '####                                             ####'
549  write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS            ####'
550  write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR       ####'
551  write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"-           ####'
552  write(*,*) '#### FILE ...                                    ####'
553  write(*,*) '#####################################################'
554  stop
555
556
557999   write(*,*) '#####################################################'
558  write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
559  write(*,*)
560  write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
561  write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
562  write(*,*) 'PERMITTED FOR ANY ACCESS'
563  write(*,*) '#####################################################'
564  stop
565
566end subroutine readreleases
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG