readreleases.f90 24.7 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

Matthias Langer's avatar
 
Matthias Langer committed
4
5
6
7
8
9
10
11
12
13
14
15
16
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                         *
17
18
  !     HSO, 12 August 2013
  !     Added optional namelist input
Matthias Langer's avatar
 
Matthias Langer committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! 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                                                   *
41
42
43
  ! 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
44
  ! zpoint1,zpoint2     height range, over which release takes place           *
45
46
  ! num_min_discrete    if less, release cannot be randomized and happens at   *
  !                     time mid-point of release interval                     *
47
  ! lroot               true if serial version, or if MPI and root process     *
Matthias Langer's avatar
 
Matthias Langer committed
48
49
50
51
52
53
54
55
56
57
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use xmass_mod
  use par_mod
  use com_mod

  implicit none

58
  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
59
  integer,parameter :: num_min_discrete=100
Matthias Langer's avatar
 
Matthias Langer committed
60
  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
61
  real(kind=dp) :: jul1,jul2,julm,juldate
62
  real,parameter :: eps2=1.e-9
Matthias Langer's avatar
 
Matthias Langer committed
63
64
65
  character(len=50) :: line
  logical :: old

66
67
68
69
70
71
72
73
74
75
76
77
  ! 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/ &
78
79
       nspec, &
       specnum_rel
80
81

  namelist /release/ &
82
83
84
85
86
87
88
89
90
       idate1, itime1, &
       idate2, itime2, &
       lon1, lon2, &
       lat1, lat2, &
       z1, z2, &
       zkind, &
       mass, &
       parts, &
       comment
Matthias Langer's avatar
 
Matthias Langer committed
91

92
  numpoint=0
Matthias Langer's avatar
 
Matthias Langer committed
93

94
95
96
97
98
  ! 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
99

100
101
102
103
104
105
106
107
108
109
110
  ! 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
111
112
  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
113
114
  endif

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

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

121
    readerror=1 ! indicates old format
Matthias Langer's avatar
 
Matthias Langer committed
122

123
124
125
126
  ! 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
127

128
129
130
131
132
133
134
135
136
    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
137

138

139
140
  ! Skip first 11 lines (file header)
  !**********************************
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

    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
162
163
    read(unitreleases,*,err=998) xdum
    if (old) call skplin(2,unitreleases)
164
165
166
167
168
169
170
171
172
173
174
175
176
177
    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
178
  !save compoint only for the first 1000 release points
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    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)

203
204
  if (nspec.gt.maxspec) goto 994

205
206
207
208
209
210
211
  ! 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)
212
  deallocate(specnum_rel) 
213
214
  ! eso: BUG, crashes here for nspec=12 and maxspec=6,
  ! TODO: catch error and exit
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
  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'

248
  if (lroot) write (*,*) 'Releasepoints : ', numpoint
249
250
251

  do i=1,numpoint
    xmasssave(i)=0.
Matthias Langer's avatar
 
Matthias Langer committed
252
253
254
255
256
257
258
259
260
  end do

  !now save the information
  DEP=.false.
  DRYDEP=.false.
  WETDEP=.false.
  OHREA=.false.
  do i=1,maxspec
    DRYDEPSPEC(i)=.false.
261
    WETDEPSPEC(i)=.false.
Matthias Langer's avatar
 
Matthias Langer committed
262
263
  end do

264
  if (readerror.ne.0) then
265
266
  ! Skip first 11 lines (file header)
  !**********************************
Matthias Langer's avatar
 
Matthias Langer committed
267

268
    call skplin(11,unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
269

270
271
  ! Assign species-specific parameters needed for physical processes
  !*************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
272

273
274
275
276
    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
277

278
  ! namelist output
279
  if (nmlout.and.lroot) then
280
281
    write(unitreleasesout,nml=releases_ctrl)
  endif
Matthias Langer's avatar
 
Matthias Langer committed
282
283

  do i=1,nspec
284
285
286
287
288
289
290
    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
291

292
293
  ! For backward runs, only 1 species is allowed
  !*********************************************
Matthias Langer's avatar
 
Matthias Langer committed
294

295
296
297
298
299
300
301
  !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
302

303
304
  ! Molecular weight
  !*****************
Matthias Langer's avatar
 
Matthias Langer committed
305

306
    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
Matthias Langer's avatar
 
Matthias Langer committed
307
308
309
310
311
312
313
      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

314
315
  ! Radioactive decay
  !******************
Matthias Langer's avatar
 
Matthias Langer committed
316
317
318
319

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


320
321
  ! Dry deposition of gases
  !************************
Matthias Langer's avatar
 
Matthias Langer committed
322

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

325
326
  ! Dry deposition of particles
  !****************************
Matthias Langer's avatar
 
Matthias Langer committed
327
328
329
330

    vsetaver(i)=0.
    cunningham(i)=0.
    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
331
332
    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
333
334
335
336
337
338
339
      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
340
      if (lroot) write(*,*) 'Average settling velocity: ',i,vsetaver(i)
Matthias Langer's avatar
 
Matthias Langer committed
341
342
    endif

343
344
  ! Dry deposition for constant deposition velocity
  !************************************************
Matthias Langer's avatar
 
Matthias Langer committed
345
346
347

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

348
349
  ! Check if wet deposition or OH reaction shall be calculated
  !***********************************************************
350
351
352
353

  ! 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
354
      WETDEP=.true.
355
      WETDEPSPEC(i)=.true.
356
      if (lroot) then
357
        write (*,*) '  Below-cloud scavenging: ON'
358
  !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
359
      end if
360
    else
361
      if (lroot) write (*,*) '  Below-cloud scavenging: OFF'
Matthias Langer's avatar
 
Matthias Langer committed
362
    endif
363
364

  ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
365
    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
366
      WETDEP=.true.
367
      WETDEPSPEC(i)=.true.
368
      if (lroot) then
369
        write (*,*) '  In-cloud scavenging: ON'
370
371
  !        write (*,*) 'In-cloud scavenging coefficients: ',&
  !           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
372
373
      end if
    else
374
      if (lroot) write (*,*) '  In-cloud scavenging: OFF' 
375
376
    endif

377
    if (ohcconst(i).gt.0.) then
Matthias Langer's avatar
 
Matthias Langer committed
378
      OHREA=.true.
379
      if (lroot) write (*,*) '  OHreaction switched on: ',ohcconst(i),i
Matthias Langer's avatar
 
Matthias Langer committed
380
381
    endif

382
    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
Matthias Langer's avatar
 
Matthias Langer committed
383
384
385
386
      DRYDEP=.true.
      DRYDEPSPEC(i)=.true.
    endif

387
  end do ! end loop over species
Matthias Langer's avatar
 
Matthias Langer committed
388

389
  if (WETDEP.or.DRYDEP) DEP=.true.
Matthias Langer's avatar
 
Matthias Langer committed
390
391
392

  ! Read specifications for each release point
  !*******************************************
393
  numpoints=numpoint
Matthias Langer's avatar
 
Matthias Langer committed
394
395
396
  numpoint=0
  numpartmax=0
  releaserate=0.
397
101 numpoint=numpoint+1
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423

  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

424
  ! namelist output
425
    if (nmlout.and.lroot) then
426
427
428
      write(unitreleasesout,nml=release)
    endif

Matthias Langer's avatar
 
Matthias Langer committed
429
430
  else

431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
    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
457
  !save compoint only for the first 1000 release points
458
459
460
461
462
463
464
465
466
    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)

467
  ! namelist output
468
    if (nmlout.and.lroot) then
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
      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
494
495
496
497
498
499
500
501
502
503
504
505

  ! 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

506
507
508
  ! If FLEXPART is run for backward deposition force zpoint
  !*********************************************************************
  if (WETBKDEP) then
509
510
511
    zpoint1(numpoint)=0.
    zpoint2(numpoint)=20000.
    kindz(numpoint)=1
512
513
  endif
  if (DRYBKDEP) then
514
515
516
    zpoint1(numpoint)=0.
    zpoint2(numpoint)=2.*href
    kindz(numpoint)=1
517
518
519
  endif


Matthias Langer's avatar
 
Matthias Langer committed
520
521
522
  ! Check whether x coordinates of release point are within model domain
  !*********************************************************************

523
524
525
526
527
528
529
530
  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.
Matthias Langer's avatar
 
Matthias Langer committed
531
532
533
534
535
536

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

  jul1=juldate(id1,it1)
  jul2=juldate(id2,it2)
537
  julm=(jul1+jul2)/2.
Matthias Langer's avatar
 
Matthias Langer committed
538
539
540
541
542
543
544
545
546
547
  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'
548
        write(*,*) 'Release starts before simulation begins or ends (1)'
Matthias Langer's avatar
 
Matthias Langer committed
549
        write(*,*) 'after simulation stops.'
550
        write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
Matthias Langer's avatar
 
Matthias Langer committed
551
552
        stop
      endif
553
554
555
556
557
558
559
      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
560
561
562
    else if (ldirect.eq.-1) then
      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
        write(*,*) 'FLEXPART MODEL ERROR'
563
        write(*,*) 'Release starts before simulation begins or ends (2)'
Matthias Langer's avatar
 
Matthias Langer committed
564
        write(*,*) 'after simulation stops.'
565
        write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
Matthias Langer's avatar
 
Matthias Langer committed
566
567
        stop
      endif
568
569
570
571
572
573
574
      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
575
576
577
    endif
  endif

578

Matthias Langer's avatar
 
Matthias Langer committed
579
580
581
582
583
584
585
586
587
588
589
  ! 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)
590
  goto 101
Matthias Langer's avatar
 
Matthias Langer committed
591

592
250 close(unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
593

594
  if (nmlout.and.lroot) then
595
596
597
    close(unitreleasesout)
  endif

598
599
  if (lroot) write (*,*) 'Particles allocated (maxpart)  : ',maxpart
  if (lroot) write (*,*) 'Particles released (numpartmax): ',numpartmax
Matthias Langer's avatar
 
Matthias Langer committed
600
601
602
603
604
605
606
607
  numpoint=numpoint-1

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

608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
  ! Disable settling if more than 1 species at any release point
  ! or if MQUASILAG and more than one species
  !*************************************************************

  if (mquasilag.ne.0) then
    if (nspec.gt.1) lsettling=.false.
  else
    do irel=1,numpoint
      nsettle=0
      do ispc=1,nspec
        if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
      end do
      if (nsettle.gt.1) lsettling=.false.
    end do
  end if

  if (lroot) then
    if (.not.lsettling) then 
      write(*,*) 'WARNING: more than 1 species per release point, settling &
           &disabled'
    end if
  end if

631
632
633
634
  ! 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
635
636
  if (releaserate.gt. &
       0.99*real(maxpart)/real(lage(nageclass))) then
637
    if (numpartmax.gt.maxpart.and.lroot) then
638
639
640
641
642
643
644
645
646
647
648
649
      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(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
650
651
652
653
654
655
656
657
658
659
660
      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

661

662
663
664
665
  if (lroot) then
    write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
  end if

Matthias Langer's avatar
 
Matthias Langer committed
666
667
  return

668
994 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
669
670
671
672
673
674
675
  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
  write(*,*) '####                                             ####'
  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
  write(*,*) '#####################################################'
  stop

676
998 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
677
678
679
680
681
682
683
684
685
686
  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


687
999 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
688
689
690
691
692
693
694
695
  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

696
697
698
699
700
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
701
end subroutine readreleases