readreleases.f90 25.1 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine readreleases

  !*****************************************************************************
  !                                                                            *
  !     This routine reads the release point specifications for the current    *
  !     model run. Several release points can be used at the same time.        *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     18 May 1996                                                            *
  !                                                                            *
  !     Update: 29 January 2001                                                *
  !     Release altitude can be either in magl or masl                         *
35
36
  !     HSO, 12 August 2013
  !     Added optional namelist input
Matthias Langer's avatar
 
Matthias Langer committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! decay               decay constant of species                              *
  ! dquer [um]          mean particle diameters                                *
  ! dsigma              e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
  !                     are between 0.1*dquer and 10*dquer                     *
  ! ireleasestart, ireleaseend [s] starting time and ending time of each       *
  !                     release                                                *
  ! kindz               1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
  !                     is in hPa                                              *
  ! npart               number of particles to be released                     *
  ! nspec               number of species to be released                       *
  ! density [kg/m3]     density of the particles                               *
  ! rm [s/m]            Mesophyll resistance                                   *
  ! species             name of species                                        *
  ! xmass               total mass of each species                             *
  ! xpoint1,ypoint1     geograf. coordinates of lower left corner of release   *
  !                     area                                                   *
  ! xpoint2,ypoint2     geograf. coordinates of upper right corner of release  *
  !                     area                                                   *
59
60
61
  ! weta_gas, wetb_gas  parameters for below-cloud scavenging (gas)            *
  ! crain_aero, csnow_aero  parameters for below-cloud scavenging (aerosol)    *
  ! ccn_aero, in_aero   parameters for in-cloud scavenging (aerosol)           *
Matthias Langer's avatar
 
Matthias Langer committed
62
  ! zpoint1,zpoint2     height range, over which release takes place           *
63
64
  ! num_min_discrete    if less, release cannot be randomized and happens at   *
  !                     time mid-point of release interval                     *
65
  ! lroot               true if serial version, or if MPI and root process     *
Matthias Langer's avatar
 
Matthias Langer committed
66
67
68
69
70
71
72
73
74
75
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use xmass_mod
  use par_mod
  use com_mod

  implicit none

76
  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
77
  integer,parameter :: num_min_discrete=100
Matthias Langer's avatar
 
Matthias Langer committed
78
  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
79
  real(kind=dp) :: jul1,jul2,julm,juldate
Matthias Langer's avatar
 
Matthias Langer committed
80
81
82
  character(len=50) :: line
  logical :: old

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  ! help variables for namelist reading
  integer :: numpoints, parts, readerror
  integer*2 :: zkind
  integer :: idate1, itime1, idate2, itime2
  real :: lon1,lon2,lat1,lat2,z1,z2
  character*40 :: comment
  integer,parameter :: unitreleasesout=2
  real,allocatable, dimension (:) :: mass
  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2

  ! declare namelists
  namelist /releases_ctrl/ &
    nspec, &
    specnum_rel

  namelist /release/ &
    idate1, itime1, &
    idate2, itime2, &
    lon1, lon2, &
    lat1, lat2, &
    z1, z2, &
    zkind, &
    mass, &
    parts, &
    comment
Matthias Langer's avatar
 
Matthias Langer committed
108

109
  numpoint=0
Matthias Langer's avatar
 
Matthias Langer committed
110

111
112
113
114
115
  ! allocate with maxspec for first input loop
  allocate(mass(maxspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
  allocate(specnum_rel(maxspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
Matthias Langer's avatar
 
Matthias Langer committed
116

117
118
119
120
121
122
123
124
125
126
127
  ! presetting namelist releases_ctrl
  nspec = -1  ! use negative value to determine failed namelist input
  specnum_rel = 0

  !sec, read release to find how many releasepoints should be allocated
  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)

  ! check if namelist input provided
  read(unitreleases,releases_ctrl,iostat=readerror)

  ! prepare namelist output if requested
128
129
  if (nmlout.and.lroot) then
    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='replace',err=1000)
Matthias Langer's avatar
 
Matthias Langer committed
130
131
  endif

132
  if ((readerror.ne.0).or.(nspec.lt.0)) then
Matthias Langer's avatar
 
Matthias Langer committed
133

134
135
136
    ! no namelist format, close file and allow reopening in old format
    close(unitreleases)
    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
Matthias Langer's avatar
 
Matthias Langer committed
137

138
    readerror=1 ! indicates old format
Matthias Langer's avatar
 
Matthias Langer committed
139

140
141
142
143
    ! Check the format of the RELEASES file (either in free format,
    ! or using a formatted mask)
    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
    !**************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
144

145
146
147
148
149
150
151
152
153
    call skplin(12,unitreleases)
    read (unitreleases,901) line
901 format (a)
    if (index(line,'Total') .eq. 0) then
      old = .false.
    else
      old = .true.
    endif
    rewind(unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
154

155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178

    ! Skip first 11 lines (file header)
    !**********************************

    call skplin(11,unitreleases)

    read(unitreleases,*,err=998) nspec
    if (old) call skplin(2,unitreleases)
    do i=1,nspec
      read(unitreleases,*,err=998) specnum_rel(i)
      if (old) call skplin(2,unitreleases)
    end do

    numpoint=0
100 numpoint=numpoint+1
    read(unitreleases,*,end=25)
    read(unitreleases,*,err=998,end=25) idum,idum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) idum,idum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
179
180
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) idum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) idum
    if (old) call skplin(2,unitreleases)
    do i=1,nspec
      read(unitreleases,*,err=998) xdum
      if (old) call skplin(2,unitreleases)
    end do
    !save compoint only for the first 1000 release points
    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
    if (old) call skplin(1,unitreleases)

    goto 100

25  numpoint=numpoint-1

  else 

    readerror=0
    do while (readerror.eq.0) 
      idate1=-1
      read(unitreleases,release,iostat=readerror)
      if ((idate1.lt.0).or.(readerror.ne.0)) then
        readerror=1
      else
        numpoint=numpoint+1
      endif
    end do
    readerror=0
  endif ! if namelist input

  rewind(unitreleases)

  ! allocate arrays of matching size for number of species (namelist output)
  deallocate(mass)
  allocate(mass(nspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
  allocate(specnum_rel2(nspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
  specnum_rel2=specnum_rel(1:nspec)
227
228
229
  deallocate(specnum_rel) 
! eso: BUG, crashes here for nspec=12 and maxspec=6,
! TODO: catch error and exit
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
  allocate(specnum_rel(nspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
  specnum_rel=specnum_rel2
  deallocate(specnum_rel2)

  !allocate memory for numpoint releaspoints
  allocate(ireleasestart(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
  allocate(ireleaseend(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
  allocate(xpoint1(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
  allocate(xpoint2(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
  allocate(ypoint1(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
  allocate(ypoint2(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
  allocate(zpoint1(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
  allocate(zpoint2(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
  allocate(kindz(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
  allocate(xmass(numpoint,maxspec),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
  allocate(rho_rel(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
  allocate(npart(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
  allocate(xmasssave(numpoint),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'

263
  if (lroot) write (*,*) 'Releasepoints : ', numpoint
264
265
266

  do i=1,numpoint
    xmasssave(i)=0.
Matthias Langer's avatar
 
Matthias Langer committed
267
268
269
270
271
272
273
274
275
276
277
  end do

  !now save the information
  DEP=.false.
  DRYDEP=.false.
  WETDEP=.false.
  OHREA=.false.
  do i=1,maxspec
    DRYDEPSPEC(i)=.false.
  end do

278
279
280
  if (readerror.ne.0) then
    ! Skip first 11 lines (file header)
    !**********************************
Matthias Langer's avatar
 
Matthias Langer committed
281

282
    call skplin(11,unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
283

284
285
    ! Assign species-specific parameters needed for physical processes
    !*************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
286

287
288
289
290
    read(unitreleases,*,err=998) nspec
    if (nspec.gt.maxspec) goto 994
    if (old) call skplin(2,unitreleases)
  endif
Matthias Langer's avatar
 
Matthias Langer committed
291

292
  ! namelist output
293
  if (nmlout.and.lroot) then
294
295
    write(unitreleasesout,nml=releases_ctrl)
  endif
Matthias Langer's avatar
 
Matthias Langer committed
296
297

  do i=1,nspec
298
299
300
301
302
303
304
    if (readerror.ne.0) then
      read(unitreleases,*,err=998) specnum_rel(i)
      if (old) call skplin(2,unitreleases)
      call readspecies(specnum_rel(i),i)
    else
      call readspecies(specnum_rel(i),i)
    endif
Matthias Langer's avatar
 
Matthias Langer committed
305

306
307
    ! For backward runs, only 1 species is allowed
    !*********************************************
Matthias Langer's avatar
 
Matthias Langer committed
308

309
310
311
312
313
314
315
    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
    !write(*,*) '#####################################################'
    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
    !write(*,*) '#####################################################'
    !  stop
    !endif
Matthias Langer's avatar
 
Matthias Langer committed
316

317
318
    ! Molecular weight
    !*****************
Matthias Langer's avatar
 
Matthias Langer committed
319

320
    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
Matthias Langer's avatar
 
Matthias Langer committed
321
322
323
324
325
326
327
      write(*,*) 'For mixing ratio output, valid molar weight'
      write(*,*) 'must be specified for all simulated species.'
      write(*,*) 'Check table SPECIES or choose concentration'
      write(*,*) 'output instead if molar weight is not known.'
      stop
    endif

328
329
    ! Radioactive decay
    !******************
Matthias Langer's avatar
 
Matthias Langer committed
330
331
332
333

    decay(i)=0.693147/decay(i) !conversion half life to decay constant


334
335
    ! Dry deposition of gases
    !************************
Matthias Langer's avatar
 
Matthias Langer committed
336

337
    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
Matthias Langer's avatar
 
Matthias Langer committed
338

339
340
    ! Dry deposition of particles
    !****************************
Matthias Langer's avatar
 
Matthias Langer committed
341
342
343
344

    vsetaver(i)=0.
    cunningham(i)=0.
    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
345
346
    if (density(i).gt.0.) then         ! Additional parameters
      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
Matthias Langer's avatar
 
Matthias Langer committed
347
348
349
350
351
352
353
      do j=1,ni
        fract(i,j)=fracth(j)
        schmi(i,j)=schmih(j)
        vset(i,j)=vsh(j)
        cunningham(i)=cunningham(i)+cun*fract(i,j)
        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
      end do
354
      if (lroot) write(*,*) 'Average settling velocity: ',i,vsetaver(i)
Matthias Langer's avatar
 
Matthias Langer committed
355
356
    endif

357
358
    ! Dry deposition for constant deposition velocity
    !************************************************
Matthias Langer's avatar
 
Matthias Langer committed
359
360
361

    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s

362
363
  ! Check if wet deposition or OH reaction shall be calculated
  !***********************************************************
364
365
366
367

  ! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
    if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
         &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.)))  then
Matthias Langer's avatar
 
Matthias Langer committed
368
      WETDEP=.true.
369
      if (lroot) then
370
371
        write (*,*) '  Below-cloud scavenging: ON'
!  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
372
      end if
373
    else
374
      if (lroot) write (*,*) '  Below-cloud scavenging: OFF'
Matthias Langer's avatar
 
Matthias Langer committed
375
    endif
376
    
377
! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
378
    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
379
      WETDEP=.true.
380
      if (lroot) then
381
        write (*,*) '  In-cloud scavenging: ON'
flexpart's avatar
flexpart committed
382
!        write (*,*) 'In-cloud scavenging coefficients: ',&
383
!           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
384
385
      end if
    else
386
      if (lroot) write (*,*) '  In-cloud scavenging: OFF' 
387
388
    endif

389
    if (ohcconst(i).gt.0.) then
Matthias Langer's avatar
 
Matthias Langer committed
390
      OHREA=.true.
391
      if (lroot) write (*,*) '  OHreaction switched on: ',ohcconst(i),i
Matthias Langer's avatar
 
Matthias Langer committed
392
393
    endif

394
    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
Matthias Langer's avatar
 
Matthias Langer committed
395
396
397
398
399
400
      DRYDEP=.true.
      DRYDEPSPEC(i)=.true.
    endif

  end do

401
  if (WETDEP.or.DRYDEP) DEP=.true.
Matthias Langer's avatar
 
Matthias Langer committed
402
403
404

  ! Read specifications for each release point
  !*******************************************
405
  numpoints=numpoint
Matthias Langer's avatar
 
Matthias Langer committed
406
407
408
  numpoint=0
  numpartmax=0
  releaserate=0.
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
101   numpoint=numpoint+1

  if (readerror.lt.1) then ! reading namelist format

    if (numpoint.gt.numpoints) goto 250
    zkind = 1
    mass = 0
    parts = 0
    comment = ' '
    read(unitreleases,release,iostat=readerror)
    id1=idate1
    it1=itime1
    id2=idate2
    it2=itime2
    xpoint1(numpoint)=lon1
    xpoint2(numpoint)=lon2
    ypoint1(numpoint)=lat1
    ypoint2(numpoint)=lat2
    zpoint1(numpoint)=z1
    zpoint2(numpoint)=z2
    kindz(numpoint)=zkind
    do i=1,nspec
      xmass(numpoint,i)=mass(i)
    end do
    npart(numpoint)=parts
    compoint(min(1001,numpoint))=comment

    ! namelist output
437
    if (nmlout.and.lroot) then
438
439
440
      write(unitreleasesout,nml=release)
    endif

Matthias Langer's avatar
 
Matthias Langer committed
441
442
  else

443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    read(unitreleases,*,end=250)
    read(unitreleases,*,err=998,end=250) id1,it1
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) id2,it2
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xpoint1(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) ypoint1(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) xpoint2(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) ypoint2(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) kindz(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) zpoint1(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) zpoint2(numpoint)
    if (old) call skplin(2,unitreleases)
    read(unitreleases,*,err=998) npart(numpoint)
    if (old) call skplin(2,unitreleases)
    do i=1,nspec
      read(unitreleases,*,err=998) xmass(numpoint,i)
      if (old) call skplin(2,unitreleases)
      mass(i)=xmass(numpoint,i)
    end do
    !save compoint only for the first 1000 release points
    if (numpoint.le.1000) then
      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
      comment=compoint(numpoint)(1:40)
    else
      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
      comment=compoint(1001)(1:40)
    endif
    if (old) call skplin(1,unitreleases)

    ! namelist output
480
    if (nmlout.and.lroot) then
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
      idate1=id1
      itime1=it1
      idate2=id2
      itime2=it2
      lon1=xpoint1(numpoint)
      lon2=xpoint2(numpoint)
      lat1=ypoint1(numpoint)
      lat2=ypoint2(numpoint)
      z1=zpoint1(numpoint)
      z2=zpoint2(numpoint)
      zkind=kindz(numpoint)
      parts=npart(numpoint)
      write(unitreleasesout,nml=release)
    endif

    if (numpoint.le.1000) then
      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
           (compoint(numpoint)(1:8).eq.'        ')) goto 250
    else
      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
    endif

  endif ! if namelist format
Matthias Langer's avatar
 
Matthias Langer committed
506
507
508
509
510
511
512
513
514
515
516
517

  ! If a release point contains no particles, stop and issue error message
  !***********************************************************************

  if (npart(numpoint).eq.0) then
    write(*,*) 'FLEXPART MODEL ERROR'
    write(*,*) 'RELEASES file is corrupt.'
    write(*,*) 'At least for one release point, there are zero'
    write(*,*) 'particles released. Make changes to RELEASES.'
    stop
  endif

518
519
520
521
522
523
524
525
526
527
528
529
530
531
  ! If FLEXPART is run for backward deposition force zpoint
  !*********************************************************************
  if (WETBKDEP) then
      zpoint1(numpoint)=0.
      zpoint2(numpoint)=20000.
      kindz(numpoint)=1
  endif
  if (DRYBKDEP) then
      zpoint1(numpoint)=0.
      zpoint2(numpoint)=2.*href
      kindz(numpoint)=1
  endif


Matthias Langer's avatar
 
Matthias Langer committed
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
  ! Check whether x coordinates of release point are within model domain
  !*********************************************************************

   if (xpoint1(numpoint).lt.xlon0) &
        xpoint1(numpoint)=xpoint1(numpoint)+360.
   if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
        xpoint1(numpoint)=xpoint1(numpoint)-360.
   if (xpoint2(numpoint).lt.xlon0) &
        xpoint2(numpoint)=xpoint2(numpoint)+360.
   if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
        xpoint2(numpoint)=xpoint2(numpoint)-360.

  ! Determine relative beginning and ending times of particle release
  !******************************************************************

  jul1=juldate(id1,it1)
  jul2=juldate(id2,it2)
549
  julm=(jul1+jul2)/2.
Matthias Langer's avatar
 
Matthias Langer committed
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
  if (jul1.gt.jul2) then
    write(*,*) 'FLEXPART MODEL ERROR'
    write(*,*) 'Release stops before it begins.'
    write(*,*) 'Make changes to file RELEASES.'
    stop
  endif
  if (mdomainfill.eq.0) then   ! no domain filling
    if (ldirect.eq.1) then
      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
        write(*,*) 'FLEXPART MODEL ERROR'
        write(*,*) 'Release starts before simulation begins or ends'
        write(*,*) 'after simulation stops.'
        write(*,*) 'Make files COMMAND and RELEASES consistent.'
        stop
      endif
565
566
567
568
569
570
571
      if (npart(numpoint).gt.num_min_discrete) then
        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
      else
        ireleasestart(numpoint)=int((julm-bdate)*86400.)
        ireleaseend(numpoint)=int((julm-bdate)*86400.)
      endif
Matthias Langer's avatar
 
Matthias Langer committed
572
573
574
575
576
577
578
579
    else if (ldirect.eq.-1) then
      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
        write(*,*) 'FLEXPART MODEL ERROR'
        write(*,*) 'Release starts before simulation begins or ends'
        write(*,*) 'after simulation stops.'
        write(*,*) 'Make files COMMAND and RELEASES consistent.'
        stop
      endif
580
581
582
583
584
585
586
      if (npart(numpoint).gt.num_min_discrete) then
        ireleasestart(numpoint)=int((jul1-bdate)*86400.)
        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
      else
        ireleasestart(numpoint)=int((julm-bdate)*86400.)
        ireleaseend(numpoint)=int((julm-bdate)*86400.)
      endif
Matthias Langer's avatar
 
Matthias Langer committed
587
588
589
590
591
592
593
594
595
596
597
598
599
600
    endif
  endif

  ! Determine the release rate (particles per second) and total number
  ! of particles released during the simulation
  !*******************************************************************

  if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
    releaserate=releaserate+real(npart(numpoint))/ &
         real(ireleaseend(numpoint)-ireleasestart(numpoint))
  else
    releaserate=99999999
  endif
  numpartmax=numpartmax+npart(numpoint)
601
  goto 101
Matthias Langer's avatar
 
Matthias Langer committed
602
603
604

250   close(unitreleases)

605
  if (nmlout.and.lroot) then
606
607
608
    close(unitreleasesout)
  endif

609
610
  if (lroot) write (*,*) 'Particles allocated (maxpart)  : ',maxpart
  if (lroot) write (*,*) 'Particles released (numpartmax): ',numpartmax
Matthias Langer's avatar
 
Matthias Langer committed
611
612
613
614
615
616
617
618
  numpoint=numpoint-1

  if (ioutputforeachrelease.eq.1) then
    maxpointspec_act=numpoint
  else
    maxpointspec_act=1
  endif

619
620
621
622
  ! Check, whether the total number of particles may exceed totally allowed
  ! number of particles at some time during the simulation
  !************************************************************************

Matthias Langer's avatar
 
Matthias Langer committed
623
624
  if (releaserate.gt. &
       0.99*real(maxpart)/real(lage(nageclass))) then
625
    if (numpartmax.gt.maxpart.and.lroot) then
Matthias Langer's avatar
 
Matthias Langer committed
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
  write(*,*) '#####################################################'
  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
  write(*,*) '####                                             ####'
  write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
  write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
  write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED   ####'
  write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
  write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
  write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP          ####'
  write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE  ####'
  write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT.   ####'
  write(*,*) '#####################################################'
      write(*,*) 'Maximum release rate may be: ',releaserate, &
           ' particles per second'
      write(*,*) 'Maximum allowed release rate is: ', &
           real(maxpart)/real(lage(nageclass)),' particles per second'
      write(*,*) &
           'Total number of particles released during the simulation is: ', &
           numpartmax
      write(*,*) 'Maximum allowed number of particles is: ',maxpart
    endif
  endif

  return

994   write(*,*) '#####################################################'
  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
  write(*,*) '####                                             ####'
  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
  write(*,*) '#####################################################'
  stop

998   write(*,*) '#####################################################'
  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
  write(*,*) '####                                             ####'
  write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS            ####'
  write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR       ####'
  write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"-           ####'
  write(*,*) '#### FILE ...                                    ####'
  write(*,*) '#####################################################'
  stop


999   write(*,*) '#####################################################'
  write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
  write(*,*)
  write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
  write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
  write(*,*) 'PERMITTED FOR ANY ACCESS'
  write(*,*) '#####################################################'
  stop

679
680
681
682
683
1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
  write(*,'(a)') path(2)(1:length(2))
  stop

Matthias Langer's avatar
 
Matthias Langer committed
684
end subroutine readreleases