readspecies.f90 15.6 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, pspec_ass, pkao
70
  real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
71
72
  integer :: readerror

73
! declare namelist
74
  namelist /species_params/ &
75
76
       pspecies, pdecay, pweta_gas, pwetb_gas, &
       pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
77
       preldiff, phenry, pf0, pdensity, pdquer, &
78
       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
79

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

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
    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
167
!  write(*,*) spec_ass(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
168
    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
169
!       write(*,*) kao(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
170

171
172
    pspecies=species(pos_spec)
    pdecay=decay(pos_spec)
173
174
175
176
177
178
    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)
179
180
181
182
183
184
185
186
    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)
187
188
    pohcconst=ohcconst(pos_spec)
    pohdconst=ohdconst(pos_spec)
189
    pohnconst=ohnconst(pos_spec)
190
191
    pspec_ass=spec_ass(pos_spec)
    pkao=kao(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
192

193
194
195
196
  else

    species(pos_spec)=pspecies
    decay(pos_spec)=pdecay
197
198
199
200
201
202
    weta_gas(pos_spec)=pweta_gas
    wetb_gas(pos_spec)=pwetb_gas
    crain_aero=pcrain_aero
    csnow_aero=pcsnow_aero
    ccn_aero(pos_spec)=pccn_aero
    in_aero(pos_spec)=pin_aero
203
204
205
206
207
208
209
210
    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
211
212
    ohcconst(pos_spec)=pohcconst
    ohdconst(pos_spec)=pohdconst
213
    ohnconst(pos_spec)=pohnconst
214
215
216
217
    spec_ass(pos_spec)=pspec_ass
    kao(pos_spec)=pkao

  endif
Matthias Langer's avatar
 
Matthias Langer committed
218

219
  i=pos_spec
Matthias Langer's avatar
 
Matthias Langer committed
220

221
222
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
223

224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
  if (lroot) then
! 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
247
      end if
248
249
250
251
      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
252
      end if
253
254
      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
255
        write(*,*) '*******************************************'
256
257
258
259
260
        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
261
262
        write(*,*) '*******************************************'
      end if
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
    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
283
          write(*,*) '*******************************************'
284
285
286
          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
287
288
289
          write(*,*) '*******************************************'
        end if
      endif
290
291
292
293

      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
294
295
      endif
    end if
296
  end if
297
298
299
300
301
302
303
304
305
306

  if (spec_ass(pos_spec).gt.0) then
    spec_found=.FALSE.
    do j=1,pos_spec-1
      if (spec_ass(pos_spec).eq.specnum(j)) then
        spec_ass(pos_spec)=j
        spec_found=.TRUE.
        ASSSPEC=.TRUE.
      endif
    end do
307
    if (spec_found.eqv..false.) then
308
      goto 997
Matthias Langer's avatar
 
Matthias Langer committed
309
    endif
310
311
312
313
314
315
316
317
318
319
320
321
  endif

  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
  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
322
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
323
324


325
326
327
! Read in daily and day-of-week variation of emissions, if available
!*******************************************************************
! HSO: This is not yet implemented as namelist parameters
Matthias Langer's avatar
 
Matthias Langer committed
328

329
330
331
332
333
334
335
336
337
338
  do j=1,24           ! initialize everything to no variation
    area_hour(i,j)=1.
    point_hour(i,j)=1.
  end do
  do j=1,7
    area_dow(i,j)=1.
    point_dow(i,j)=1.
  end do

  if (readerror.ne.0) then ! text format input
Matthias Langer's avatar
 
Matthias Langer committed
339
340
341
342
343
344
345
346
347
348

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

349
350
351
  endif

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

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

  return
Matthias Langer's avatar
 
Matthias Langer committed
361

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


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


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

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

end subroutine readspecies