source: flexpart.git/src/readreleases.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 5184a7c, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Enable settling with multiple species if from separate releases

  • Property mode set to 100644
File size: 26.0 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_gas, wetb_gas  parameters for below-cloud scavenging (gas)            *
60  ! crain_aero, csnow_aero  parameters for below-cloud scavenging (aerosol)    *
61  ! ccn_aero, in_aero   parameters for in-cloud scavenging (aerosol)           *
62  ! zpoint1,zpoint2     height range, over which release takes place           *
63  ! num_min_discrete    if less, release cannot be randomized and happens at   *
64  !                     time mid-point of release interval                     *
65  ! lroot               true if serial version, or if MPI and root process     *
66  !                                                                            *
67  !*****************************************************************************
68
69  use point_mod
70  use xmass_mod
71  use par_mod
72  use com_mod
73
74  implicit none
75
76  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
77  integer,parameter :: num_min_discrete=100
78  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
79  real(kind=dp) :: jul1,jul2,julm,juldate
80  real,parameter :: eps2=1.e-9
81  character(len=50) :: line
82  logical :: old
83
84  ! help variables for namelist reading
85  integer :: numpoints, parts, readerror
86  integer*2 :: zkind
87  integer :: idate1, itime1, idate2, itime2
88  real :: lon1,lon2,lat1,lat2,z1,z2
89  character*40 :: comment
90  integer,parameter :: unitreleasesout=2
91  real,allocatable, dimension (:) :: mass
92  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
93
94  ! declare namelists
95  namelist /releases_ctrl/ &
96       nspec, &
97       specnum_rel
98
99  namelist /release/ &
100       idate1, itime1, &
101       idate2, itime2, &
102       lon1, lon2, &
103       lat1, lat2, &
104       z1, z2, &
105       zkind, &
106       mass, &
107       parts, &
108       comment
109
110  numpoint=0
111
112  ! allocate with maxspec for first input loop
113  allocate(mass(maxspec),stat=stat)
114  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
115  allocate(specnum_rel(maxspec),stat=stat)
116  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
117
118  ! presetting namelist releases_ctrl
119  nspec = -1  ! use negative value to determine failed namelist input
120  specnum_rel = 0
121
122  !sec, read release to find how many releasepoints should be allocated
123  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
124
125  ! check if namelist input provided
126  read(unitreleases,releases_ctrl,iostat=readerror)
127
128  ! prepare namelist output if requested
129  if (nmlout.and.lroot) then
130    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='replace',err=1000)
131  endif
132
133  if ((readerror.ne.0).or.(nspec.lt.0)) then
134
135  ! no namelist format, close file and allow reopening in old format
136    close(unitreleases)
137    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
138
139    readerror=1 ! indicates old format
140
141  ! Check the format of the RELEASES file (either in free format,
142  ! or using a formatted mask)
143  ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
144  !**************************************************************************
145
146    call skplin(12,unitreleases)
147    read (unitreleases,901) line
148901 format (a)
149    if (index(line,'Total') .eq. 0) then
150      old = .false.
151    else
152      old = .true.
153    endif
154    rewind(unitreleases)
155
156
157  ! Skip first 11 lines (file header)
158  !**********************************
159
160    call skplin(11,unitreleases)
161
162    read(unitreleases,*,err=998) nspec
163    if (old) call skplin(2,unitreleases)
164    do i=1,nspec
165      read(unitreleases,*,err=998) specnum_rel(i)
166      if (old) call skplin(2,unitreleases)
167    end do
168
169    numpoint=0
170100 numpoint=numpoint+1
171    read(unitreleases,*,end=25)
172    read(unitreleases,*,err=998,end=25) idum,idum
173    if (old) call skplin(2,unitreleases)
174    read(unitreleases,*,err=998) idum,idum
175    if (old) call skplin(2,unitreleases)
176    read(unitreleases,*,err=998) xdum
177    if (old) call skplin(2,unitreleases)
178    read(unitreleases,*,err=998) xdum
179    if (old) call skplin(2,unitreleases)
180    read(unitreleases,*,err=998) xdum
181    if (old) call skplin(2,unitreleases)
182    read(unitreleases,*,err=998) xdum
183    if (old) call skplin(2,unitreleases)
184    read(unitreleases,*,err=998) idum
185    if (old) call skplin(2,unitreleases)
186    read(unitreleases,*,err=998) xdum
187    if (old) call skplin(2,unitreleases)
188    read(unitreleases,*,err=998) xdum
189    if (old) call skplin(2,unitreleases)
190    read(unitreleases,*,err=998) idum
191    if (old) call skplin(2,unitreleases)
192    do i=1,nspec
193      read(unitreleases,*,err=998) xdum
194      if (old) call skplin(2,unitreleases)
195    end do
196  !save compoint only for the first 1000 release points
197    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
198    if (old) call skplin(1,unitreleases)
199
200    goto 100
201
20225  numpoint=numpoint-1
203
204  else
205
206    readerror=0
207    do while (readerror.eq.0)
208      idate1=-1
209      read(unitreleases,release,iostat=readerror)
210      if ((idate1.lt.0).or.(readerror.ne.0)) then
211        readerror=1
212      else
213        numpoint=numpoint+1
214      endif
215    end do
216    readerror=0
217  endif ! if namelist input
218
219  rewind(unitreleases)
220
221  if (nspec.gt.maxspec) goto 994
222
223  ! allocate arrays of matching size for number of species (namelist output)
224  deallocate(mass)
225  allocate(mass(nspec),stat=stat)
226  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
227  allocate(specnum_rel2(nspec),stat=stat)
228  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
229  specnum_rel2=specnum_rel(1:nspec)
230  deallocate(specnum_rel)
231  ! eso: BUG, crashes here for nspec=12 and maxspec=6,
232  ! TODO: catch error and exit
233  allocate(specnum_rel(nspec),stat=stat)
234  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
235  specnum_rel=specnum_rel2
236  deallocate(specnum_rel2)
237
238  !allocate memory for numpoint releaspoints
239  allocate(ireleasestart(numpoint),stat=stat)
240  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
241  allocate(ireleaseend(numpoint),stat=stat)
242  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
243  allocate(xpoint1(numpoint),stat=stat)
244  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
245  allocate(xpoint2(numpoint),stat=stat)
246  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
247  allocate(ypoint1(numpoint),stat=stat)
248  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
249  allocate(ypoint2(numpoint),stat=stat)
250  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
251  allocate(zpoint1(numpoint),stat=stat)
252  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
253  allocate(zpoint2(numpoint),stat=stat)
254  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
255  allocate(kindz(numpoint),stat=stat)
256  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
257  allocate(xmass(numpoint,maxspec),stat=stat)
258  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
259  allocate(rho_rel(numpoint),stat=stat)
260  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
261  allocate(npart(numpoint),stat=stat)
262  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
263  allocate(xmasssave(numpoint),stat=stat)
264  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
265
266  if (lroot) write (*,*) 'Releasepoints : ', numpoint
267
268  do i=1,numpoint
269    xmasssave(i)=0.
270  end do
271
272  !now save the information
273  DEP=.false.
274  DRYDEP=.false.
275  WETDEP=.false.
276  OHREA=.false.
277  do i=1,maxspec
278    DRYDEPSPEC(i)=.false.
279    WETDEPSPEC(i)=.false.
280  end do
281
282  if (readerror.ne.0) then
283  ! Skip first 11 lines (file header)
284  !**********************************
285
286    call skplin(11,unitreleases)
287
288  ! Assign species-specific parameters needed for physical processes
289  !*************************************************************************
290
291    read(unitreleases,*,err=998) nspec
292    if (nspec.gt.maxspec) goto 994
293    if (old) call skplin(2,unitreleases)
294  endif
295
296  ! namelist output
297  if (nmlout.and.lroot) then
298    write(unitreleasesout,nml=releases_ctrl)
299  endif
300
301  do i=1,nspec
302    if (readerror.ne.0) then
303      read(unitreleases,*,err=998) specnum_rel(i)
304      if (old) call skplin(2,unitreleases)
305      call readspecies(specnum_rel(i),i)
306    else
307      call readspecies(specnum_rel(i),i)
308    endif
309
310  ! For backward runs, only 1 species is allowed
311  !*********************************************
312
313  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
314  !write(*,*) '#####################################################'
315  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
316  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
317  !write(*,*) '#####################################################'
318  !  stop
319  !endif
320
321  ! Molecular weight
322  !*****************
323
324    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
325      write(*,*) 'For mixing ratio output, valid molar weight'
326      write(*,*) 'must be specified for all simulated species.'
327      write(*,*) 'Check table SPECIES or choose concentration'
328      write(*,*) 'output instead if molar weight is not known.'
329      stop
330    endif
331
332  ! Radioactive decay
333  !******************
334
335    decay(i)=0.693147/decay(i) !conversion half life to decay constant
336
337
338  ! Dry deposition of gases
339  !************************
340
341    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
342
343  ! Dry deposition of particles
344  !****************************
345
346    vsetaver(i)=0.
347    cunningham(i)=0.
348    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
349    if (density(i).gt.0.) then         ! Additional parameters
350      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
351      do j=1,ni
352        fract(i,j)=fracth(j)
353        schmi(i,j)=schmih(j)
354        vset(i,j)=vsh(j)
355        cunningham(i)=cunningham(i)+cun*fract(i,j)
356        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
357      end do
358      if (lroot) write(*,*) 'Average settling velocity: ',i,vsetaver(i)
359    endif
360
361  ! Dry deposition for constant deposition velocity
362  !************************************************
363
364    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
365
366  ! Check if wet deposition or OH reaction shall be calculated
367  !***********************************************************
368
369  ! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
370    if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
371         &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.)))  then
372      WETDEP=.true.
373      WETDEPSPEC(i)=.true.
374      if (lroot) then
375        write (*,*) '  Below-cloud scavenging: ON'
376  !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
377      end if
378    else
379      if (lroot) write (*,*) '  Below-cloud scavenging: OFF'
380    endif
381
382  ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
383    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
384      WETDEP=.true.
385      WETDEPSPEC(i)=.true.
386      if (lroot) then
387        write (*,*) '  In-cloud scavenging: ON'
388  !        write (*,*) 'In-cloud scavenging coefficients: ',&
389  !           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
390      end if
391    else
392      if (lroot) write (*,*) '  In-cloud scavenging: OFF'
393    endif
394
395    if (ohcconst(i).gt.0.) then
396      OHREA=.true.
397      if (lroot) write (*,*) '  OHreaction switched on: ',ohcconst(i),i
398    endif
399
400    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
401      DRYDEP=.true.
402      DRYDEPSPEC(i)=.true.
403    endif
404
405  end do ! end loop over species
406
407  if (WETDEP.or.DRYDEP) DEP=.true.
408
409  ! Read specifications for each release point
410  !*******************************************
411  numpoints=numpoint
412  numpoint=0
413  numpartmax=0
414  releaserate=0.
415101 numpoint=numpoint+1
416
417  if (readerror.lt.1) then ! reading namelist format
418
419    if (numpoint.gt.numpoints) goto 250
420    zkind = 1
421    mass = 0
422    parts = 0
423    comment = ' '
424    read(unitreleases,release,iostat=readerror)
425    id1=idate1
426    it1=itime1
427    id2=idate2
428    it2=itime2
429    xpoint1(numpoint)=lon1
430    xpoint2(numpoint)=lon2
431    ypoint1(numpoint)=lat1
432    ypoint2(numpoint)=lat2
433    zpoint1(numpoint)=z1
434    zpoint2(numpoint)=z2
435    kindz(numpoint)=zkind
436    do i=1,nspec
437      xmass(numpoint,i)=mass(i)
438    end do
439    npart(numpoint)=parts
440    compoint(min(1001,numpoint))=comment
441
442  ! namelist output
443    if (nmlout.and.lroot) then
444      write(unitreleasesout,nml=release)
445    endif
446
447  else
448
449    read(unitreleases,*,end=250)
450    read(unitreleases,*,err=998,end=250) id1,it1
451    if (old) call skplin(2,unitreleases)
452    read(unitreleases,*,err=998) id2,it2
453    if (old) call skplin(2,unitreleases)
454    read(unitreleases,*,err=998) xpoint1(numpoint)
455    if (old) call skplin(2,unitreleases)
456    read(unitreleases,*,err=998) ypoint1(numpoint)
457    if (old) call skplin(2,unitreleases)
458    read(unitreleases,*,err=998) xpoint2(numpoint)
459    if (old) call skplin(2,unitreleases)
460    read(unitreleases,*,err=998) ypoint2(numpoint)
461    if (old) call skplin(2,unitreleases)
462    read(unitreleases,*,err=998) kindz(numpoint)
463    if (old) call skplin(2,unitreleases)
464    read(unitreleases,*,err=998) zpoint1(numpoint)
465    if (old) call skplin(2,unitreleases)
466    read(unitreleases,*,err=998) zpoint2(numpoint)
467    if (old) call skplin(2,unitreleases)
468    read(unitreleases,*,err=998) npart(numpoint)
469    if (old) call skplin(2,unitreleases)
470    do i=1,nspec
471      read(unitreleases,*,err=998) xmass(numpoint,i)
472      if (old) call skplin(2,unitreleases)
473      mass(i)=xmass(numpoint,i)
474    end do
475  !save compoint only for the first 1000 release points
476    if (numpoint.le.1000) then
477      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
478      comment=compoint(numpoint)(1:40)
479    else
480      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
481      comment=compoint(1001)(1:40)
482    endif
483    if (old) call skplin(1,unitreleases)
484
485  ! namelist output
486    if (nmlout.and.lroot) then
487      idate1=id1
488      itime1=it1
489      idate2=id2
490      itime2=it2
491      lon1=xpoint1(numpoint)
492      lon2=xpoint2(numpoint)
493      lat1=ypoint1(numpoint)
494      lat2=ypoint2(numpoint)
495      z1=zpoint1(numpoint)
496      z2=zpoint2(numpoint)
497      zkind=kindz(numpoint)
498      parts=npart(numpoint)
499      write(unitreleasesout,nml=release)
500    endif
501
502    if (numpoint.le.1000) then
503      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
504           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
505           (compoint(numpoint)(1:8).eq.'        ')) goto 250
506    else
507      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
508           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
509    endif
510
511  endif ! if namelist format
512
513  ! If a release point contains no particles, stop and issue error message
514  !***********************************************************************
515
516  if (npart(numpoint).eq.0) then
517    write(*,*) 'FLEXPART MODEL ERROR'
518    write(*,*) 'RELEASES file is corrupt.'
519    write(*,*) 'At least for one release point, there are zero'
520    write(*,*) 'particles released. Make changes to RELEASES.'
521    stop
522  endif
523
524  ! If FLEXPART is run for backward deposition force zpoint
525  !*********************************************************************
526  if (WETBKDEP) then
527    zpoint1(numpoint)=0.
528    zpoint2(numpoint)=20000.
529    kindz(numpoint)=1
530  endif
531  if (DRYBKDEP) then
532    zpoint1(numpoint)=0.
533    zpoint2(numpoint)=2.*href
534    kindz(numpoint)=1
535  endif
536
537
538  ! Check whether x coordinates of release point are within model domain
539  !*********************************************************************
540
541  if (xpoint1(numpoint).lt.xlon0) &
542       xpoint1(numpoint)=xpoint1(numpoint)+360.
543  if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
544       xpoint1(numpoint)=xpoint1(numpoint)-360.
545  if (xpoint2(numpoint).lt.xlon0) &
546       xpoint2(numpoint)=xpoint2(numpoint)+360.
547  if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
548       xpoint2(numpoint)=xpoint2(numpoint)-360.
549
550  ! Determine relative beginning and ending times of particle release
551  !******************************************************************
552
553  jul1=juldate(id1,it1)
554  jul2=juldate(id2,it2)
555  julm=(jul1+jul2)/2.
556  if (jul1.gt.jul2) then
557    write(*,*) 'FLEXPART MODEL ERROR'
558    write(*,*) 'Release stops before it begins.'
559    write(*,*) 'Make changes to file RELEASES.'
560    stop
561  endif
562  if (mdomainfill.eq.0) then   ! no domain filling
563    if (ldirect.eq.1) then
564      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
565        write(*,*) 'FLEXPART MODEL ERROR'
566        write(*,*) 'Release starts before simulation begins or ends'
567        write(*,*) 'after simulation stops.'
568        write(*,*) 'Make files COMMAND and RELEASES consistent.'
569        stop
570      endif
571      if (npart(numpoint).gt.num_min_discrete) then
572        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
573        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
574      else
575        ireleasestart(numpoint)=int((julm-bdate)*86400.)
576        ireleaseend(numpoint)=int((julm-bdate)*86400.)
577      endif
578    else if (ldirect.eq.-1) then
579      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
580        write(*,*) 'FLEXPART MODEL ERROR'
581        write(*,*) 'Release starts before simulation begins or ends'
582        write(*,*) 'after simulation stops.'
583        write(*,*) 'Make files COMMAND and RELEASES consistent.'
584        stop
585      endif
586      if (npart(numpoint).gt.num_min_discrete) then
587        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
588        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
589      else
590        ireleasestart(numpoint)=int((julm-bdate)*86400.)
591        ireleaseend(numpoint)=int((julm-bdate)*86400.)
592      endif
593    endif
594  endif
595
596
597  ! Determine the release rate (particles per second) and total number
598  ! of particles released during the simulation
599  !*******************************************************************
600
601  if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
602    releaserate=releaserate+real(npart(numpoint))/ &
603         real(ireleaseend(numpoint)-ireleasestart(numpoint))
604  else
605    releaserate=99999999
606  endif
607  numpartmax=numpartmax+npart(numpoint)
608  goto 101
609
610250 close(unitreleases)
611
612  if (nmlout.and.lroot) then
613    close(unitreleasesout)
614  endif
615
616  if (lroot) write (*,*) 'Particles allocated (maxpart)  : ',maxpart
617  if (lroot) write (*,*) 'Particles released (numpartmax): ',numpartmax
618  numpoint=numpoint-1
619
620  if (ioutputforeachrelease.eq.1) then
621    maxpointspec_act=numpoint
622  else
623    maxpointspec_act=1
624  endif
625
626  ! Disable settling if more than 1 species at any release point
627  ! or if MQUASILAG and more than one species
628  !*************************************************************
629
630  if (mquasilag.ne.0) then
631    if (nspec.gt.1) lsettling=.false.
632  else
633    do irel=1,numpoint
634      nsettle=0
635      do ispc=1,nspec
636        if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
637      end do
638      if (nsettle.gt.1) lsettling=.false.
639    end do
640  end if
641
642  if (lroot) then
643    if (.not.lsettling) then
644      write(*,*) 'WARNING: more than 1 species per release point, settling &
645           &disabled'
646    end if
647  end if
648
649  ! Check, whether the total number of particles may exceed totally allowed
650  ! number of particles at some time during the simulation
651  !************************************************************************
652
653  if (releaserate.gt. &
654       0.99*real(maxpart)/real(lage(nageclass))) then
655    if (numpartmax.gt.maxpart.and.lroot) then
656      write(*,*) '#####################################################'
657      write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
658      write(*,*) '####                                             ####'
659      write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
660      write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
661      write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
662      write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
663      write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
664      write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
665      write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
666      write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
667      write(*,*) '#####################################################'
668      write(*,*) 'Maximum release rate may be: ',releaserate, &
669           ' particles per second'
670      write(*,*) 'Maximum allowed release rate is: ', &
671           real(maxpart)/real(lage(nageclass)),' particles per second'
672      write(*,*) &
673           'Total number of particles released during the simulation is: ', &
674           numpartmax
675      write(*,*) 'Maximum allowed number of particles is: ',maxpart
676    endif
677  endif
678
679
680  if (lroot) then
681    write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
682  end if
683
684  return
685
686994 write(*,*) '#####################################################'
687  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
688  write(*,*) '####                                             ####'
689  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
690  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
691  write(*,*) '#####################################################'
692  stop
693
694998 write(*,*) '#####################################################'
695  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
696  write(*,*) '####                                             ####'
697  write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS            ####'
698  write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR       ####'
699  write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"-           ####'
700  write(*,*) '#### FILE ...                                    ####'
701  write(*,*) '#####################################################'
702  stop
703
704
705999 write(*,*) '#####################################################'
706  write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
707  write(*,*)
708  write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
709  write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
710  write(*,*) 'PERMITTED FOR ANY ACCESS'
711  write(*,*) '#####################################################'
712  stop
713
7141000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
715  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
716  write(*,'(a)') path(2)(1:length(2))
717  stop
718
719end subroutine readreleases
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG