source: branches/FP_AI/src/readreleases.f90 @ 23

Last change on this file since 23 was 23, checked in by igpis, 10 years ago

start tracking test environment directory FP_AI

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