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

  use point_mod
  use xmass_mod
  use par_mod
  use com_mod

  implicit none

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

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

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

89
  numpoint=0
Matthias Langer's avatar
 
Matthias Langer committed
90

91
92
93
94
95
  ! 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
96

97
98
99
100
101
102
103
104
105
106
107
  ! 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
108
109
  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
110
111
  endif

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

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

118
    readerror=1 ! indicates old format
Matthias Langer's avatar
 
Matthias Langer committed
119

120
121
122
123
  ! 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
124

125
126
127
128
129
130
131
132
133
    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
134

135

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

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

200
201
  if (nspec.gt.maxspec) goto 994

202
203
204
205
206
207
208
  ! 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)
209
  deallocate(specnum_rel) 
210
211
  ! eso: BUG, crashes here for nspec=12 and maxspec=6,
  ! TODO: catch error and exit
212
213
214
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
  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'

245
  if (lroot) write (*,*) 'Releasepoints : ', numpoint
246
247
248

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

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

261
  if (readerror.ne.0) then
262
263
  ! Skip first 11 lines (file header)
  !**********************************
Matthias Langer's avatar
 
Matthias Langer committed
264

265
    call skplin(11,unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
266

267
268
  ! Assign species-specific parameters needed for physical processes
  !*************************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
269

270
271
272
273
    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
274

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

  do i=1,nspec
281
282
283
284
285
286
287
    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
288

289
290
  ! For backward runs, only 1 species is allowed
  !*********************************************
Matthias Langer's avatar
 
Matthias Langer committed
291

292
293
294
295
296
297
298
  !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
299

300
301
  ! Molecular weight
  !*****************
Matthias Langer's avatar
 
Matthias Langer committed
302

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

311
312
  ! Radioactive decay
  !******************
Matthias Langer's avatar
 
Matthias Langer committed
313
314
315
316

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


317
318
  ! Dry deposition of gases
  !************************
Matthias Langer's avatar
 
Matthias Langer committed
319

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

322
323
  ! Dry deposition of particles
  !****************************
Matthias Langer's avatar
 
Matthias Langer committed
324
325
326
327

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

340
341
  ! Dry deposition for constant deposition velocity
  !************************************************
Matthias Langer's avatar
 
Matthias Langer committed
342
343
344

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

345
346
  ! Check if wet deposition or OH reaction shall be calculated
  !***********************************************************
347
348
349
350

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

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

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

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

384
  end do ! end loop over species
Matthias Langer's avatar
 
Matthias Langer committed
385

386
  if (WETDEP.or.DRYDEP) DEP=.true.
Matthias Langer's avatar
 
Matthias Langer committed
387
388
389

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

  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

421
  ! namelist output
422
    if (nmlout.and.lroot) then
423
424
425
      write(unitreleasesout,nml=release)
    endif

Matthias Langer's avatar
 
Matthias Langer committed
426
427
  else

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

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

  ! 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

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


Matthias Langer's avatar
 
Matthias Langer committed
517
518
519
  ! Check whether x coordinates of release point are within model domain
  !*********************************************************************

520
521
522
523
524
525
526
527
  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
528
529
530
531
532
533

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

  jul1=juldate(id1,it1)
  jul2=juldate(id2,it2)
534
  julm=(jul1+jul2)/2.
Matthias Langer's avatar
 
Matthias Langer committed
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
  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
550
551
552
553
554
555
556
      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
557
558
559
560
561
562
563
564
    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
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
    endif
  endif

575

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

589
250 close(unitreleases)
Matthias Langer's avatar
 
Matthias Langer committed
590

591
  if (nmlout.and.lroot) then
592
593
594
    close(unitreleasesout)
  endif

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

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

605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
  ! 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

628
629
630
631
  ! 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
632
633
  if (releaserate.gt. &
       0.99*real(maxpart)/real(lage(nageclass))) then
634
    if (numpartmax.gt.maxpart.and.lroot) then
635
636
637
638
639
640
641
642
643
644
645
646
      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
647
648
649
650
651
652
653
654
655
656
657
      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

658

659
660
661
662
  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
663
664
  return

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

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


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

693
694
695
696
697
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
698
end subroutine readreleases