readspecies.f90 16.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
!**********************************************************************
! 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 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                                                             *
32
  !                                                                            *
33
34
35
  !   Changes:                                                                 *
  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
  !                                                                            *
36
37
38
  !   HSO, 13 August 2013
  !   added optional namelist input
  !                                                                            *
Matthias Langer's avatar
 
Matthias Langer committed
39
40
41
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
42
43
44
45
46
47
48
49
50
51
  ! 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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod

  implicit none

  integer :: i, pos_spec,j
  integer :: idow,ihour,id_spec
  character(len=3) :: aspecnumb
  logical :: spec_found

67
  character(len=16) :: pspecies
68
  real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
69
  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
70
  real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
Sabine's avatar
Sabine committed
71
  real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
72
  integer :: readerror
73

74
! declare namelist
75
  namelist /species_params/ &
76
77
       pspecies, pdecay, pweta_gas, pwetb_gas, &
       pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
78
       preldiff, phenry, pf0, pdensity, pdquer, &
Sabine's avatar
Sabine committed
79
80
       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
       parea_dow, parea_hour, ppoint_dow, ppoint_hour
81

82
  pspecies="" ! read failure indicator value
83
  pdecay=-999.9
84
85
86
87
88
89
  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
90
91
92
93
94
95
96
  preldiff=-9.9
  phenry=0.0
  pf0=0.0
  pdensity=-9.9E09
  pdquer=0.0
  pdsigma=0.0
  pdryvel=-9.99
97
98
  pohcconst=-9.99
  pohdconst=-9.9E-09
99
  pohnconst=2.0
100
  pweightmolar=-999.9
101

Sabine's avatar
Sabine committed
102
  do j=1,24           ! initialize everything to no variation
103
104
    parea_hour(j)=1.
    ppoint_hour(j)=1.
Sabine's avatar
Sabine committed
105
106
107
108
    area_hour(pos_spec,j)=1.
    point_hour(pos_spec,j)=1.
  end do
  do j=1,7
109
110
    parea_dow(j)=1.
    ppoint_dow(j)=1.
Sabine's avatar
Sabine committed
111
112
113
114
115
    area_dow(pos_spec,j)=1.
    point_dow(pos_spec,j)=1.
  end do

  if (readerror.ne.0) then ! text format input
116
117
! Open the SPECIES file and read species names and properties
!************************************************************
Matthias Langer's avatar
 
Matthias Langer committed
118
119
  specnum(pos_spec)=id_spec
  write(aspecnumb,'(i3.3)') specnum(pos_spec)
120
  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
121
!write(*,*) 'reading SPECIES',specnum(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
122
123
124

  ASSSPEC=.FALSE.

125
! try namelist input
126
127
  read(unitspecies,species_params,iostat=readerror)
  close(unitspecies)
128

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

    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
140
141

    read(unitspecies,'(a10)',end=22) species(pos_spec)
142
!  write(*,*) species(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
143
    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
144
!  write(*,*) decay(pos_spec)
145
146
147
148
149
150
151
152
    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)
153
!*** NIK 31.01.2013: including in-cloud scavening parameters
154
155
156
157
    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
158
    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
159
!  write(*,*) reldiff(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
160
    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
161
!  write(*,*) henry(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
162
    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
163
!  write(*,*) f0(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
164
    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
165
!  write(*,*) density(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
166
    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
167
!  write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
168
    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
169
!  write(*,*) dsigma(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
170
    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
171
!  write(*,*) dryvel(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
172
    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
173
174
175
176
177
!  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)
178
179
    read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
!  write(*,*) ohnconst(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
180

Sabine's avatar
Sabine committed
181
182
183
184
185
186
187
188
189
190
191
192
! 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

193
194
    pspecies=species(pos_spec)
    pdecay=decay(pos_spec)
195
196
197
198
199
200
    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)
201
202
203
204
205
206
207
208
    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)
209
210
    pohcconst=ohcconst(pos_spec)
    pohdconst=ohdconst(pos_spec)
211
    pohnconst=ohnconst(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
212

Sabine's avatar
Sabine committed
213
214
215
216
217
218
219
220
221
222
223

    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
224
225
226

    species(pos_spec)=pspecies
    decay(pos_spec)=pdecay
227
228
    weta_gas(pos_spec)=pweta_gas
    wetb_gas(pos_spec)=pwetb_gas
229
230
    crain_aero(pos_spec)=pcrain_aero
    csnow_aero(pos_spec)=pcsnow_aero
231
232
    ccn_aero(pos_spec)=pccn_aero
    in_aero(pos_spec)=pin_aero
233
234
235
236
237
238
239
240
    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
241
242
    ohcconst(pos_spec)=pohcconst
    ohdconst(pos_spec)=pohdconst
243
    ohnconst(pos_spec)=pohnconst
Sabine's avatar
Sabine committed
244
245
246
247
248
249
250
251
252
253

    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)
    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)
    end do

254
  endif
Matthias Langer's avatar
 
Matthias Langer committed
255

256
  i=pos_spec
Matthias Langer's avatar
 
Matthias Langer committed
257

258
259
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
260

261
  if (lroot) then
262
263
264
265
266
267
268
! 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) '

269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
! 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
291
      end if
292
293
294
295
      if (density(pos_spec) .gt. 0) then
        write(*,'(a)') '  Dry deposition is turned         :   ON'
      else
        write(*,'(a)') '  Dry deposition is (density<0)    :   OFF'
Espen Sollum's avatar
Espen Sollum committed
296
      end if
297
298
      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
299
        write(*,*) '*******************************************'
300
301
302
303
304
        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
305
306
        write(*,*) '*******************************************'
      end if
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
    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
327
          write(*,*) '*******************************************'
328
329
330
          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
331
332
333
          write(*,*) '*******************************************'
        end if
      endif
334
335
336
337

      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
338
339
      endif
    end if
340
  end if
341
342
343
344
345
346
347
348
349
350

  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception

  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
351
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
352
353


354
355
356
  endif

22 close(unitspecies)
Matthias Langer's avatar
 
Matthias Langer committed
357

358
359
360
! 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)
361
362
363
364
365
    write(unitspecies,nml=species_params)
    close(unitspecies)
  endif

  return
Matthias Langer's avatar
 
Matthias Langer committed
366

367
996 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
368
369
370
371
372
373
374
375
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
  write(*,*) '#### CONSTANT IS SET                            ####'
  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
  write(*,*) '#####################################################'
  stop


376
997 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
377
378
379
380
381
382
383
384
385
  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


386
998 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
387
388
389
390
391
392
393
  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

394
395
396
397
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
398
399

end subroutine readspecies