source: flexpart.git/src/readreleases.f90 @ b7ae015

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025bugfixes+enhancementsdevfp9.3.1-20161214-nc4grib2nc4_repairrelease-10release-10.4.1scaling-bugunivie
Last change on this file since b7ae015 was b7ae015, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 9 years ago

clean up runtime messages and adjust verbosity level
  • Property mode set to 100755
File size: 25.4 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 below-cloud scavenging     *
60  ! weta_in, wetb_in    parameters to determine the in-cloud scavenging        *
61  ! wetc_in, wetd_in    parameters to determine the in-cloud scavenging        *
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  !                                                                            *
66  !*****************************************************************************
67
68  use point_mod
69  use xmass_mod
70  use par_mod
71  use com_mod
72
73  implicit none
74
75  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
76  integer,parameter :: num_min_discrete=100
77  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
78  real(kind=dp) :: jul1,jul2,julm,juldate
79  character(len=50) :: line
80  logical :: old
81
82  ! help variables for namelist reading
83  integer :: numpoints, parts, readerror
84  integer*2 :: zkind
85  integer :: idate1, itime1, idate2, itime2
86  real :: lon1,lon2,lat1,lat2,z1,z2
87  character*40 :: comment
88  integer,parameter :: unitreleasesout=2
89  real,allocatable, dimension (:) :: mass
90  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
91
92  ! declare namelists
93  namelist /releases_ctrl/ &
94    nspec, &
95    specnum_rel
96
97  namelist /release/ &
98    idate1, itime1, &
99    idate2, itime2, &
100    lon1, lon2, &
101    lat1, lat2, &
102    z1, z2, &
103    zkind, &
104    mass, &
105    parts, &
106    comment
107
108  numpoint=0
109
110  ! allocate with maxspec for first input loop
111  allocate(mass(maxspec),stat=stat)
112  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
113  allocate(specnum_rel(maxspec),stat=stat)
114  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
115
116  ! presetting namelist releases_ctrl
117  nspec = -1  ! use negative value to determine failed namelist input
118  specnum_rel = 0
119
120  !sec, read release to find how many releasepoints should be allocated
121  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
122
123  ! check if namelist input provided
124  read(unitreleases,releases_ctrl,iostat=readerror)
125
126  ! prepare namelist output if requested
127  if (nmlout.eqv..true.) then
128    !open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
129    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',err=1000)
130  endif
131
132  if ((readerror.ne.0).or.(nspec.lt.0)) then
133
134    ! no namelist format, close file and allow reopening in old format
135    close(unitreleases)
136    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
137
138    readerror=1 ! indicates old format
139
140    ! Check the format of the RELEASES file (either in free format,
141    ! or using a formatted mask)
142    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
143    !**************************************************************************
144
145    call skplin(12,unitreleases)
146    read (unitreleases,901) line
147901 format (a)
148    if (index(line,'Total') .eq. 0) then
149      old = .false.
150    else
151      old = .true.
152    endif
153    rewind(unitreleases)
154
155
156    ! Skip first 11 lines (file header)
157    !**********************************
158
159    call skplin(11,unitreleases)
160
161    read(unitreleases,*,err=998) nspec
162    if (old) call skplin(2,unitreleases)
163    do i=1,nspec
164      read(unitreleases,*,err=998) specnum_rel(i)
165      if (old) call skplin(2,unitreleases)
166    end do
167
168    numpoint=0
169100 numpoint=numpoint+1
170    read(unitreleases,*,end=25)
171    read(unitreleases,*,err=998,end=25) idum,idum
172    if (old) call skplin(2,unitreleases)
173    read(unitreleases,*,err=998) idum,idum
174    if (old) call skplin(2,unitreleases)
175    read(unitreleases,*,err=998) xdum
176    if (old) call skplin(2,unitreleases)
177    read(unitreleases,*,err=998) xdum
178    if (old) call skplin(2,unitreleases)
179    read(unitreleases,*,err=998) xdum
180    if (old) call skplin(2,unitreleases)
181    read(unitreleases,*,err=998) xdum
182    if (old) call skplin(2,unitreleases)
183    read(unitreleases,*,err=998) idum
184    if (old) call skplin(2,unitreleases)
185    read(unitreleases,*,err=998) xdum
186    if (old) call skplin(2,unitreleases)
187    read(unitreleases,*,err=998) xdum
188    if (old) call skplin(2,unitreleases)
189    read(unitreleases,*,err=998) idum
190    if (old) call skplin(2,unitreleases)
191    do i=1,nspec
192      read(unitreleases,*,err=998) xdum
193      if (old) call skplin(2,unitreleases)
194    end do
195    !save compoint only for the first 1000 release points
196    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
197    if (old) call skplin(1,unitreleases)
198
199    goto 100
200
20125  numpoint=numpoint-1
202
203  else
204
205    readerror=0
206    do while (readerror.eq.0)
207      idate1=-1
208      read(unitreleases,release,iostat=readerror)
209      if ((idate1.lt.0).or.(readerror.ne.0)) then
210        readerror=1
211      else
212        numpoint=numpoint+1
213      endif
214    end do
215    readerror=0
216  endif ! if namelist input
217
218  rewind(unitreleases)
219
220  ! allocate arrays of matching size for number of species (namelist output)
221  deallocate(mass)
222  allocate(mass(nspec),stat=stat)
223  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
224  allocate(specnum_rel2(nspec),stat=stat)
225  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
226  specnum_rel2=specnum_rel(1:nspec)
227  deallocate(specnum_rel)
228  allocate(specnum_rel(nspec),stat=stat)
229  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
230  specnum_rel=specnum_rel2
231  deallocate(specnum_rel2)
232
233  !allocate memory for numpoint releaspoints
234  allocate(ireleasestart(numpoint),stat=stat)
235  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
236  allocate(ireleaseend(numpoint),stat=stat)
237  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
238  allocate(xpoint1(numpoint),stat=stat)
239  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
240  allocate(xpoint2(numpoint),stat=stat)
241  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
242  allocate(ypoint1(numpoint),stat=stat)
243  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
244  allocate(ypoint2(numpoint),stat=stat)
245  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
246  allocate(zpoint1(numpoint),stat=stat)
247  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
248  allocate(zpoint2(numpoint),stat=stat)
249  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
250  allocate(kindz(numpoint),stat=stat)
251  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
252  allocate(xmass(numpoint,maxspec),stat=stat)
253  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
254  allocate(rho_rel(numpoint),stat=stat)
255  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
256  allocate(npart(numpoint),stat=stat)
257  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
258  allocate(xmasssave(numpoint),stat=stat)
259  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
260
261  if (verbosity.gt.0) then
262    write (*,*) 'readreleases> Releasepoints : ', numpoint
263  endif
264
265  do i=1,numpoint
266    xmasssave(i)=0.
267  end do
268
269  !now save the information
270  DEP=.false.
271  DRYDEP=.false.
272  WETDEP=.false.
273  OHREA=.false.
274  do i=1,maxspec
275    DRYDEPSPEC(i)=.false.
276  end do
277
278  if (readerror.ne.0) then
279    ! Skip first 11 lines (file header)
280    !**********************************
281
282    call skplin(11,unitreleases)
283
284    ! Assign species-specific parameters needed for physical processes
285    !*************************************************************************
286
287    read(unitreleases,*,err=998) nspec
288    if (nspec.gt.maxspec) goto 994
289    if (old) call skplin(2,unitreleases)
290  endif
291
292  ! namelist output
293  if (nmlout.eqv..true.) then
294    write(unitreleasesout,nml=releases_ctrl)
295  endif
296
297  do i=1,nspec
298    if (verbosity.gt.0) then
299      print*, 'readreleases> call readspecies', i
300    endif
301 
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      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    if (weta(i).gt.0.)  then
369      WETDEP=.true.
370      write (*,*) 'Below-cloud scavenging: ON'
371      if (verbosity.gt.0) then
372      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
373      endif
374    else
375      if (verbosity.gt.0) then
376      write (*,*) 'Below-cloud scavenging: OFF'
377      endif
378    endif
379   
380    ! NIK 31.01.2013 + 10.12.2013
381    if (weta_in(i).gt.0.)  then
382      WETDEP=.true.
383      write (*,*) 'In-cloud scavenging: ON'
384      if (verbosity.gt.0) then
385      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
386      endif 
387  else
388      if (verbosity.gt.0) then
389      write (*,*) 'In-cloud scavenging: OFF'
390      endif
391    endif
392
393    if (ohreact(i).gt.0) then
394      OHREA=.true.
395      write (*,*) 'OHreaction: ON (',ohreact(i),i,')'
396    endif
397
398    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
399      DRYDEP=.true.
400      DRYDEPSPEC(i)=.true.
401    endif
402
403  end do
404
405  if (WETDEP.or.DRYDEP) DEP=.true.
406
407  ! Read specifications for each release point
408  !*******************************************
409  numpoints=numpoint
410  numpoint=0
411  numpartmax=0
412  releaserate=0.
413101   numpoint=numpoint+1
414
415  if (readerror.lt.1) then ! reading namelist format
416
417    if (numpoint.gt.numpoints) goto 250
418    zkind = 1
419    mass = 0
420    parts = 0
421    comment = ' '
422    read(unitreleases,release,iostat=readerror)
423    id1=idate1
424    it1=itime1
425    id2=idate2
426    it2=itime2
427    xpoint1(numpoint)=lon1
428    xpoint2(numpoint)=lon2
429    ypoint1(numpoint)=lat1
430    ypoint2(numpoint)=lat2
431    zpoint1(numpoint)=z1
432    zpoint2(numpoint)=z2
433    kindz(numpoint)=zkind
434    do i=1,nspec
435      xmass(numpoint,i)=mass(i)
436    end do
437    npart(numpoint)=parts
438    compoint(min(1001,numpoint))=comment
439
440    ! namelist output
441    if (nmlout.eqv..true.) then
442      write(unitreleasesout,nml=release)
443    endif
444
445  else
446
447    read(unitreleases,*,end=250)
448    read(unitreleases,*,err=998,end=250) id1,it1
449    if (old) call skplin(2,unitreleases)
450    read(unitreleases,*,err=998) id2,it2
451    if (old) call skplin(2,unitreleases)
452    read(unitreleases,*,err=998) xpoint1(numpoint)
453    if (old) call skplin(2,unitreleases)
454    read(unitreleases,*,err=998) ypoint1(numpoint)
455    if (old) call skplin(2,unitreleases)
456    read(unitreleases,*,err=998) xpoint2(numpoint)
457    if (old) call skplin(2,unitreleases)
458    read(unitreleases,*,err=998) ypoint2(numpoint)
459    if (old) call skplin(2,unitreleases)
460    read(unitreleases,*,err=998) kindz(numpoint)
461    if (old) call skplin(2,unitreleases)
462    read(unitreleases,*,err=998) zpoint1(numpoint)
463    if (old) call skplin(2,unitreleases)
464    read(unitreleases,*,err=998) zpoint2(numpoint)
465    if (old) call skplin(2,unitreleases)
466    read(unitreleases,*,err=998) npart(numpoint)
467    if (old) call skplin(2,unitreleases)
468    do i=1,nspec
469      read(unitreleases,*,err=998) xmass(numpoint,i)
470      if (old) call skplin(2,unitreleases)
471      mass(i)=xmass(numpoint,i)
472    end do
473    !save compoint only for the first 1000 release points
474    if (numpoint.le.1000) then
475      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
476      comment=compoint(numpoint)(1:40)
477    else
478      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
479      comment=compoint(1001)(1:40)
480    endif
481    if (old) call skplin(1,unitreleases)
482
483    ! namelist output
484    if (nmlout.eqv..true.) then
485      idate1=id1
486      itime1=it1
487      idate2=id2
488      itime2=it2
489      lon1=xpoint1(numpoint)
490      lon2=xpoint2(numpoint)
491      lat1=ypoint1(numpoint)
492      lat2=ypoint2(numpoint)
493      z1=zpoint1(numpoint)
494      z2=zpoint2(numpoint)
495      zkind=kindz(numpoint)
496      parts=npart(numpoint)
497      write(unitreleasesout,nml=release)
498    endif
499
500    if (numpoint.le.1000) then
501      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
502           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
503           (compoint(numpoint)(1:8).eq.'        ')) goto 250
504    else
505      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
506           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
507    endif
508
509  endif ! if namelist format
510
511
512  if (verbosity.gt.1 .and. numpoint.eq.1) then ! verbosity 2 or larger
513    write(*,*) 'numpoint=', numpoint
514    print*,  id1,it1
515    print*,  id2,it2
516    print*,  xpoint1(numpoint)
517    print*,  ypoint1(numpoint)
518    print*,  xpoint2(numpoint)
519    print*,  ypoint2(numpoint)
520    print*,  'kindz=' , kindz(numpoint)
521    print*,  zpoint1(numpoint)
522    print*,  zpoint2(numpoint)
523    print*,  npart(numpoint)
524    do i=1,nspec
525      !mass(i)=
526      print*, 'xmass=', xmass(numpoint,i)
527    end do
528    print*, compoint(numpoint)
529  endif
530
531
532  ! If a release point contains no particles, stop and issue error message
533  !***********************************************************************
534
535  if (npart(numpoint).eq.0) then
536    write(*,*) 'FLEXPART MODEL ERROR'
537    write(*,*) 'RELEASES file is corrupt.'
538    write(*,*) 'At least for one release point, there are zero'
539    write(*,*) 'particles released. Make changes to RELEASES.'
540    stop
541  endif
542
543  ! Check whether x coordinates of release point are within model domain
544  !*********************************************************************
545
546   if (xpoint1(numpoint).lt.xlon0) &
547        xpoint1(numpoint)=xpoint1(numpoint)+360.
548   if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
549        xpoint1(numpoint)=xpoint1(numpoint)-360.
550   if (xpoint2(numpoint).lt.xlon0) &
551        xpoint2(numpoint)=xpoint2(numpoint)+360.
552   if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
553        xpoint2(numpoint)=xpoint2(numpoint)-360.
554
555  ! Determine relative beginning and ending times of particle release
556  !******************************************************************
557
558  jul1=juldate(id1,it1)
559  jul2=juldate(id2,it2)
560  julm=(jul1+jul2)/2.
561  if (jul1.gt.jul2) then
562    write(*,*) 'FLEXPART MODEL ERROR'
563    write(*,*) 'Release stops before it begins.'
564    write(*,*) 'Make changes to file RELEASES.'
565    stop
566  endif
567  if (mdomainfill.eq.0) then   ! no domain filling
568    if (ldirect.eq.1) then
569      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
570        write(*,*) 'FLEXPART MODEL ERROR'
571        write(*,*) 'Release starts before simulation begins or ends'
572        write(*,*) 'after simulation stops.'
573        write(*,*) 'Make files COMMAND and RELEASES consistent.'
574        write(*,*) jul1, ' < ' , bdate
575        write(*,*) ' .or. '
576        write(*,*) jul2 , ' > ', edate
577       
578        stop
579      endif
580      if (npart(numpoint).gt.num_min_discrete) then
581        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
582        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
583      else
584        ireleasestart(numpoint)=int((julm-bdate)*86400.)
585        ireleaseend(numpoint)=int((julm-bdate)*86400.)
586      endif
587    else if (ldirect.eq.-1) then
588      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
589        write(*,*) 'FLEXPART MODEL ERROR'
590        write(*,*) 'Release starts before simulation begins or ends'
591        write(*,*) 'after simulation stops.'
592        write(*,*) 'Make files COMMAND and RELEASES consistent.'
593        stop
594      endif
595      if (npart(numpoint).gt.num_min_discrete) then
596        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
597        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
598      else
599        ireleasestart(numpoint)=int((julm-bdate)*86400.)
600        ireleaseend(numpoint)=int((julm-bdate)*86400.)
601      endif
602    endif
603  endif
604
605  if (verbosity.gt.1 .and. numpoint.eq.1) then ! verbosity 2 or larger
606    print*, 'ireleasestart(',numpoint,')', ireleasestart(numpoint)
607    print*, 'ireleaseend(',numpoint,')', ireleaseend(numpoint)
608  endif
609
610  ! Determine the release rate (particles per second) and total number
611  ! of particles released during the simulation
612  !*******************************************************************
613
614  if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
615    releaserate=releaserate+real(npart(numpoint))/ &
616         real(ireleaseend(numpoint)-ireleasestart(numpoint))
617  else
618    releaserate=99999999
619  endif
620  numpartmax=numpartmax+npart(numpoint)
621  goto 101
622
623250   close(unitreleases)
624
625  if (nmlout.eqv..true.) then
626    close(unitreleasesout)
627  endif
628
629  write (*,*) 'Particles allocated (maxpart)  : ',maxpart
630  write (*,*) 'Particles released (numpartmax): ',numpartmax
631  numpoint=numpoint-1
632
633  if (ioutputforeachrelease.eq.1) then
634    maxpointspec_act=numpoint
635  else
636    maxpointspec_act=1
637  endif
638
639  ! Check, whether the total number of particles may exceed totally allowed
640  ! number of particles at some time during the simulation
641  !************************************************************************
642
643  if (releaserate.gt. &
644       0.99*real(maxpart)/real(lage(nageclass))) then
645    if (numpartmax.gt.maxpart) then
646  write(*,*) '#####################################################'
647  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
648  write(*,*) '####                                             ####'
649  write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
650  write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
651  write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
652  write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
653  write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
654  write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
655  write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
656  write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
657  write(*,*) '#####################################################'
658      write(*,*) 'Maximum release rate may be: ',releaserate, &
659           ' particles per second'
660      write(*,*) 'Maximum allowed release rate is: ', &
661           real(maxpart)/real(lage(nageclass)),' particles per second'
662      write(*,*) &
663           'Total number of particles released during the simulation is: ', &
664           numpartmax
665      write(*,*) 'Maximum allowed number of particles is: ',maxpart
666    endif
667  endif
668
669  return
670
671994   write(*,*) '#####################################################'
672  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
673  write(*,*) '####                                             ####'
674  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
675  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
676  write(*,*) '#####################################################'
677  stop
678
679998   write(*,*) '#####################################################'
680  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
681  write(*,*) '####                                             ####'
682  write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS            ####'
683  write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR       ####'
684  write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"-           ####'
685  write(*,*) '#### FILE ...                                    ####'
686  write(*,*) '#####################################################'
687  stop
688
689
690999   write(*,*) '#####################################################'
691  write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
692  write(*,*)
693  write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
694  write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
695  write(*,*) 'PERMITTED FOR ANY ACCESS'
696  write(*,*) '#####################################################'
697  stop
698
6991000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
700  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
701  write(*,'(a)') path(2)(1:length(2))
702  stop
703
704end subroutine readreleases
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG