readspecies.f90 16.2 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
subroutine readspecies(id_spec,pos_spec)

  !*****************************************************************************
  !                                                                            *
  !     This routine reads names and physical constants of chemical species/   *
  !     radionuclides given in the parameter pos_spec                          *
  !                                                                            *
  !   Author: A. Stohl                                                         *
  !                                                                            *
  !   11 July 1996                                                             *
14
  !                                                                            *
15
16
17
  !   Changes:                                                                 *
  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
  !                                                                            *
18
19
20
  !   HSO, 13 August 2013
  !   added optional namelist input
  !                                                                            *
Matthias Langer's avatar
 
Matthias Langer committed
21
22
23
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
24
25
26
27
28
29
30
31
32
33
  ! decaytime(maxtable)   half time for radiological decay                     *
  ! specname(maxtable)    names of chemical species, radionuclides             *
  ! weta_gas, wetb_gas    Parameters for below-cloud scavenging of gasses      *
  ! crain_aero,csnow_aero Parameters for below-cloud scavenging of aerosols    *
  ! ccn_aero,in_aero      Parameters for in-cloud scavenging of aerosols       *
  ! ohcconst              OH reaction rate constant C                          *
  ! ohdconst              OH reaction rate constant D                          *
  ! ohnconst              OH reaction rate constant n                          *
  ! id_spec               SPECIES number as referenced in RELEASE file         *
  ! id_pos                position where SPECIES data shall be stored          *
Matthias Langer's avatar
 
Matthias Langer committed
34
35
36
37
38
39
40
41
42
43
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod

  implicit none

44
  integer :: i, pos_spec,j,icheck_dow_hour
Matthias Langer's avatar
 
Matthias Langer committed
45
46
47
48
  integer :: idow,ihour,id_spec
  character(len=3) :: aspecnumb
  logical :: spec_found

49
  character(len=16) :: pspecies
50
  real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
51
  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
52
  real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
Sabine's avatar
Sabine committed
53
  real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
54
  integer :: readerror
55

56
! declare namelist
57
  namelist /species_params/ &
58
59
       pspecies, pdecay, pweta_gas, pwetb_gas, &
       pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
60
       preldiff, phenry, pf0, pdensity, pdquer, &
Sabine's avatar
Sabine committed
61
62
       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
       parea_dow, parea_hour, ppoint_dow, ppoint_hour
63

64
  pspecies="" ! read failure indicator value
65
  pdecay=-999.9
66
67
68
69
70
71
  pweta_gas=-9.9E-09
  pwetb_gas=0.0
  pcrain_aero=-9.9E-09
  pcsnow_aero=-9.9E-09
  pccn_aero=-9.9E-09
  pin_aero=-9.9E-09
72
73
74
75
76
77
78
  preldiff=-9.9
  phenry=0.0
  pf0=0.0
  pdensity=-9.9E09
  pdquer=0.0
  pdsigma=0.0
  pdryvel=-9.99
79
80
  pohcconst=-9.99
  pohdconst=-9.9E-09
81
  pohnconst=2.0
82
  pweightmolar=-999.9
83
84
85
86
87
  parea_dow=-999.9
  parea_hour=-999.9
  ppoint_dow=-999.9
  ppoint_hour=-999.9

88

Sabine's avatar
Sabine committed
89
  do j=1,24           ! initialize everything to no variation
90
91
    parea_hour(j)=1.
    ppoint_hour(j)=1.
Sabine's avatar
Sabine committed
92
93
94
95
    area_hour(pos_spec,j)=1.
    point_hour(pos_spec,j)=1.
  end do
  do j=1,7
96
97
    parea_dow(j)=1.
    ppoint_dow(j)=1.
Sabine's avatar
Sabine committed
98
99
100
101
    area_dow(pos_spec,j)=1.
    point_dow(pos_spec,j)=1.
  end do

102
103
! Open the SPECIES file and read species names and properties
!************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
104
105
  specnum(pos_spec)=id_spec
  write(aspecnumb,'(i3.3)') specnum(pos_spec)
106
  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
107
!write(*,*) 'reading SPECIES',specnum(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
108
109
110

  ASSSPEC=.FALSE.

111
! try namelist input
112
113
  read(unitspecies,species_params,iostat=readerror)
  close(unitspecies)
114

115
116
117
  if ((len(trim(pspecies)).eq.0).or.(readerror.ne.0)) then ! no namelist found
    if (lroot) write(*,*) "SPECIES file not in NAMELIST format, attempting to &
         &read as fixed format"
118
119
120
121
122
123
124
125

    readerror=1

    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)

    do i=1,6
      read(unitspecies,*)
    end do
Matthias Langer's avatar
 
Matthias Langer committed
126
127

    read(unitspecies,'(a10)',end=22) species(pos_spec)
128
!  write(*,*) species(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
129
    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
130
!  write(*,*) decay(pos_spec)
131
132
133
134
135
136
137
138
    read(unitspecies,'(e18.1)',end=22) weta_gas(pos_spec)
!  write(*,*) weta_gas(pos_spec)
    read(unitspecies,'(f18.2)',end=22) wetb_gas(pos_spec)
!  write(*,*) wetb_gas(pos_spec)
    read(unitspecies,'(e18.1)',end=22) crain_aero(pos_spec)
!  write(*,*) crain_aero(pos_spec)
    read(unitspecies,'(f18.2)',end=22) csnow_aero(pos_spec)
!  write(*,*) csnow_aero(pos_spec)
139
!*** NIK 31.01.2013: including in-cloud scavening parameters
140
141
142
143
    read(unitspecies,'(e18.1)',end=22) ccn_aero(pos_spec)
!  write(*,*) ccn_aero(pos_spec)
    read(unitspecies,'(f18.2)',end=22) in_aero(pos_spec)
!  write(*,*) in_aero(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
144
    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
145
!  write(*,*) reldiff(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
146
    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
147
!  write(*,*) henry(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
148
    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
149
!  write(*,*) f0(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
150
    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
151
!  write(*,*) density(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
152
    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
153
!  write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
154
    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
155
!  write(*,*) dsigma(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
156
    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
157
!  write(*,*) dryvel(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
158
    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
159
160
161
162
163
!  write(*,*) weightmolar(pos_spec)
    read(unitspecies,'(e18.2)',end=22) ohcconst(pos_spec)
!  write(*,*) ohcconst(pos_spec)
    read(unitspecies,'(f8.2)',end=22) ohdconst(pos_spec)
!  write(*,*) ohdconst(pos_spec)
164
165
    read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
!  write(*,*) ohnconst(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
166

Sabine's avatar
Sabine committed
167
168
169
170
171
172
173
174
175
176
177
178
! Read in daily and day-of-week variation of emissions, if available
!*******************************************************************

    read(unitspecies,*,end=22)
    do j=1,24     ! 24 hours, starting with 0-1 local time
      read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
    end do
    read(unitspecies,*)
    do j=1,7      ! 7 days of the week, starting with Monday
      read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
    end do

179
180
    pspecies=species(pos_spec)
    pdecay=decay(pos_spec)
181
182
183
184
185
186
    pweta_gas=weta_gas(pos_spec)
    pwetb_gas=wetb_gas(pos_spec)
    pcrain_aero=crain_aero(pos_spec)
    pcsnow_aero=csnow_aero(pos_spec)
    pccn_aero=ccn_aero(pos_spec)
    pin_aero=in_aero(pos_spec)
187
188
189
190
191
192
193
194
    preldiff=reldiff(pos_spec)
    phenry=henry(pos_spec)
    pf0=f0(pos_spec)
    pdensity=density(pos_spec)
    pdquer=dquer(pos_spec)
    pdsigma=dsigma(pos_spec)
    pdryvel=dryvel(pos_spec)
    pweightmolar=weightmolar(pos_spec)
195
196
    pohcconst=ohcconst(pos_spec)
    pohdconst=ohdconst(pos_spec)
197
    pohnconst=ohnconst(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
198

Sabine's avatar
Sabine committed
199
200
201
202
203
204
205
206
207
208
209

    do j=1,24     ! 24 hours, starting with 0-1 local time
      parea_hour(j)=area_hour(pos_spec,j)
      ppoint_hour(j)=point_hour(pos_spec,j)
    end do
    do j=1,7      ! 7 days of the week, starting with Monday
      parea_dow(j)=area_dow(pos_spec,j)
      ppoint_dow(j)=point_dow(pos_spec,j)
    end do

  else ! namelist available
210
211
212

    species(pos_spec)=pspecies
    decay(pos_spec)=pdecay
213
214
    weta_gas(pos_spec)=pweta_gas
    wetb_gas(pos_spec)=pwetb_gas
215
216
    crain_aero(pos_spec)=pcrain_aero
    csnow_aero(pos_spec)=pcsnow_aero
217
218
    ccn_aero(pos_spec)=pccn_aero
    in_aero(pos_spec)=pin_aero
219
220
221
222
223
224
225
226
    reldiff(pos_spec)=preldiff
    henry(pos_spec)=phenry
    f0(pos_spec)=pf0
    density(pos_spec)=pdensity
    dquer(pos_spec)=pdquer
    dsigma(pos_spec)=pdsigma
    dryvel(pos_spec)=pdryvel
    weightmolar(pos_spec)=pweightmolar
227
228
    ohcconst(pos_spec)=pohcconst
    ohdconst(pos_spec)=pohdconst
229
    ohnconst(pos_spec)=pohnconst
Sabine's avatar
Sabine committed
230

231
232

    icheck_dow_hour=0
Sabine's avatar
Sabine committed
233
234
235
    do j=1,24     ! 24 hours, starting with 0-1 local time
      area_hour(pos_spec,j)=parea_hour(j)
      point_hour(pos_spec,j)=ppoint_hour(j)
236
      if (parea_hour(j).ne.1 .or. ppoint_hour(j).ne.1) icheck_dow_hour=1
Sabine's avatar
Sabine committed
237
238
239
240
    end do
    do j=1,7      ! 7 days of the week, starting with Monday
      area_dow(pos_spec,j)=parea_dow(j)
      point_dow(pos_spec,j)=ppoint_dow(j)
241
      if (parea_dow(j).ne.1 .or. ppoint_dow(j).ne.1) icheck_dow_hour=1
Sabine's avatar
Sabine committed
242
243
    end do

244
  endif
Matthias Langer's avatar
 
Matthias Langer committed
245

246
  i=pos_spec
Matthias Langer's avatar
 
Matthias Langer committed
247

248
249
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
250

251
  if (lroot) then
252
253
254
255
256
257
258
! ZHG 2016.04.07 Start of changes
    write(*,*) ' '
    if (dquer(pos_spec) .gt.0)  write(*,'(a,i3,a,a,a)')       ' SPECIES: ', &
         &id_spec,'  ', species(pos_spec),'  (AEROSOL) '
    if (dquer(pos_spec) .le.0)  write(*,'(a,i3,a,a,a)')       ' SPECIES: ', &
         &id_spec,'  ', species(pos_spec),'  (GAS) '

259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
! Particles
!**********
    if (dquer(pos_spec).gt.0) then
      if (ccn_aero(pos_spec) .gt. 0) then
        write(*,'(a,f5.2)') '  Particle CCN  efficiency (CCNeff):', ccn_aero(pos_spec)
      else 
        write(*,'(a)')      '  Particle CCN  efficiency (CCNeff):   OFF'
      endif
      if (in_aero(pos_spec) .gt. 0) then
        write(*,'(a,f5.2)') '  Particle  IN  efficiency (INeff) :', in_aero(pos_spec)
      else
        write(*,'(a)')      '  Particle  IN  efficiency (INeff) :   OFF'
      endif
      if (crain_aero(pos_spec) .gt. 0) then
        write(*,'(a,f5.2)') '  Particle Rain efficiency (Crain) :', crain_aero(pos_spec)
      else
        write(*,'(a)')      '  Particle Rain efficiency (Crain) :   OFF'
      endif
      if (csnow_aero(pos_spec) .gt. 0) then
        write(*,'(a,f5.2)') '  Particle Snow efficiency (Csnow) :', csnow_aero(pos_spec)
      else
        write(*,'(a)')      '  Particle Snow efficiency (Csnow) :   OFF'
Espen Sollum's avatar
Espen Sollum committed
281
      end if
282
283
      if (density(pos_spec) .gt. 0) then
        write(*,'(a)') '  Dry deposition is turned         :   ON'
284
285
286
        if (reldiff(pos_spec).gt.0) then
           stop 'density>0 (SPECIES is a particle) implies reldiff <=0  '
        endif
287
      else
288
289
290
        if (reldiff(pos_spec).le.0) then
           stop 'density<=0 (SPECIES is a gas) implies reldiff >0  '
        endif      
291
        write(*,'(a)') '  Dry deposition is (density<0)    :   OFF'
Espen Sollum's avatar
Espen Sollum committed
292
      end if
293
294
      if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. &
           &ccn_aero(pos_spec).gt.1.0 .or. in_aero(pos_spec).gt.1.0) then
Espen Sollum's avatar
Espen Sollum committed
295
        write(*,*) '*******************************************'
296
297
298
299
300
        write(*,*) ' WARNING: Particle Scavenging parameter likely out of range '
        write(*,*) '       Likely   range for Crain    0.0-10'
        write(*,*) '       Likely   range for Csnow    0.0-10'
        write(*,*) '       Physical range for CCNeff   0.0-1'
        write(*,*) '       Physical range for INeff    0.0-1'
Espen Sollum's avatar
Espen Sollum committed
301
302
        write(*,*) '*******************************************'
      end if
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
    else
! Gas
!****
      if (weta_gas(pos_spec) .gt. 0 .and. wetb_gas(pos_spec).gt.0) then
        write(*,*)          '  Wet removal for gases      is turned: ON'
        write(*,*)          '  Gas below-cloud scavenging parameter A  ', &
             &weta_gas(pos_spec)
        write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', &
             &wetb_gas(pos_spec)
      else
        write(*,*)          '  Wet removal for gases      is turned: OFF '
      end if
      if (reldiff(i).gt.0.) then
        write(*,*)          '  Dry deposition for gases   is turned: ON '
      else
        write(*,*)          '  Dry deposition for gases   is turned: OFF '
      end if
      if (weta_gas(pos_spec).gt.0.) then !if wet deposition is turned on
        if (weta_gas(pos_spec).gt.1E-04 .or. weta_gas(pos_spec).lt.1E-09 .or. &
             &wetb_gas(pos_spec).gt.0.8 .or. wetb_gas(pos_spec).lt.0.4) then
Espen Sollum's avatar
Espen Sollum committed
323
          write(*,*) '*******************************************'
324
325
326
          write(*,*) ' WARNING: Gas below-cloud scavengig is out of likely range'
          write(*,*) '          Likely range for A is 1E-04 to 1E-08'
          write(*,*) '          Likely range for B is 0.60  to 0.80 ' 
Espen Sollum's avatar
Espen Sollum committed
327
328
329
          write(*,*) '*******************************************'
        end if
      endif
330
331
332
333

      if (((weta_gas(pos_spec).gt.0).or.(wetb_gas(pos_spec).gt.0)).and.&
           &(henry(pos_spec).le.0)) then
        if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
334
335
      endif
    end if
336
  end if
337

338
  !  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
339
  if (dquer(i).gt.0 .and. dsigma(i).le.1.) then !dsigma(i)=1.0001   ! avoid floating exception
340
341
342
343
344
345
346
347
348
349
350
351
352
    !write(*,*) '#### FLEXPART MODEL ERROR!                      ####'
    write(*,*) '#### FLEXPART MODEL WARNING                     ####'
    write(*,*) '#### in SPECIES_',aspecnumb, '                             ####'
    write(*,*) '#### from v10.4 dsigma has to be larger than 1  ####'  
    write(*,*) '#### to adapt older SPECIES files,              ####' 
    write(*,*) '#### if dsigma was < 1                          ####' 
    write(*,*) '#### use the reciprocal of the old dsigma       ####'  
    if (.not.debug_mode) then 
       stop
    else
       write(*,*) 'debug mode: continue'
    endif
  endif
353
354
355
356
357
358
359
360

  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
    write(*,*) '#### PARTICLE AND GAS.                       ####'
    write(*,*) '#### SPECIES NUMBER',aspecnumb
    stop
  endif
361
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
362

363
364
365
366
367
  if ((icheck_dow_hour.eq.1).and.(ldirect.lt.0)) then
    write(*,*) '#### FLEXPART MODEL WARNING                     ####'
    write(*,*) '#### The variation for an emission release      ####'
    write(*,*) '#### will have no effect in backward mode       ####'  
  endif
Matthias Langer's avatar
 
Matthias Langer committed
368

369
22 close(unitspecies)
Matthias Langer's avatar
 
Matthias Langer committed
370

371
372
373
! namelist output if requested
  if (nmlout.and.lroot) then
    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000)
374
375
376
377
378
    write(unitspecies,nml=species_params)
    close(unitspecies)
  endif

  return
Matthias Langer's avatar
 
Matthias Langer committed
379

380
996 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
381
382
383
384
385
386
387
388
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
  write(*,*) '#### CONSTANT IS SET                            ####'
  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
  write(*,*) '#####################################################'
  stop


389
997 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
390
391
392
393
394
395
396
397
398
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### THE ASSSOCIATED SPECIES HAS TO BE DEFINED  #### '
  write(*,*) '#### BEFORE THE ONE WHICH POINTS AT IT          #### '
  write(*,*) '#### PLEASE CHANGE ORDER IN RELEASES OR ADD     #### '
  write(*,*) '#### THE ASSOCIATED SPECIES IN RELEASES         #### '
  write(*,*) '#####################################################'
  stop


399
998 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
400
401
402
403
404
405
406
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### THE SPECIES FILE FOR SPECIES ', id_spec
  write(*,*) '#### CANNOT BE FOUND: CREATE FILE'
  write(*,*) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb
  write(*,*) '#####################################################'
  stop

407
408
409
410
1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
  write(*,'(a)') path(2)(1:length(2))
  stop
Matthias Langer's avatar
 
Matthias Langer committed
411
412

end subroutine readspecies