readspecies.f90 16 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
103
104
105
106
107
108
109
110
111
  do j=1,24           ! initialize everything to no variation
    area_hour(pos_spec,j)=1.
    point_hour(pos_spec,j)=1.
  end do
  do j=1,7
    area_dow(pos_spec,j)=1.
    point_dow(pos_spec,j)=1.
  end do

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

  ASSSPEC=.FALSE.

121
! try namelist input
122
123
  read(unitspecies,species_params,iostat=readerror)
  close(unitspecies)
124

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

    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
136
137

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

Sabine's avatar
Sabine committed
177
178
179
180
181
182
183
184
185
186
187
188
! 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

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

Sabine's avatar
Sabine committed
209
210
211
212
213
214
215
216
217
218
219

    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
220
221
222

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

    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

250
  endif
Matthias Langer's avatar
 
Matthias Langer committed
251

252
  i=pos_spec
Matthias Langer's avatar
 
Matthias Langer committed
253

254
255
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
256

257
  if (lroot) then
258
259
260
261
262
263
264
! 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) '

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
! 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
287
      end if
288
289
290
291
      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
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
339
340
341
342
343
344
345
346

  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
347
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
348
349


350
351
352
  endif

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

354
355
356
! 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)
357
358
359
360
361
    write(unitspecies,nml=species_params)
    close(unitspecies)
  endif

  return
Matthias Langer's avatar
 
Matthias Langer committed
362

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


372
997 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
373
374
375
376
377
378
379
380
381
  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


382
998 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
383
384
385
386
387
388
389
  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

390
391
392
393
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
394
395

end subroutine readspecies