readspecies.f90 14.7 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
42
43
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! decaytime(maxtable)  half time for radiological decay                      *
  ! specname(maxtable)   names of chemical species, radionuclides              *
44
45
46
  ! weta, wetb           Parameters for determining below-cloud scavenging     *
  ! weta_in              Parameter for determining in-cloud scavenging         *
  ! wetb_in              Parameter for determining in-cloud scavenging         *
47
48
  ! ohcconst             OH reaction rate constant C                           *
  ! ohdconst             OH reaction rate constant D                           *
Matthias Langer's avatar
 
Matthias Langer committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
  ! id_spec              SPECIES number as referenced in RELEASE file          *
  ! id_pos               position where SPECIES data shall be stored           *
  !                                                                            *
  ! 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

66
67
  character(len=16) :: pspecies
  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
68
  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
69
70
71
  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
  integer :: readerror

72
! declare namelist
73
  namelist /species_params/ &
74
75
76
77
       pspecies, pdecay, pweta, pwetb, &
       pweta_in, pwetb_in, pwetc_in, pwetd_in, &
       preldiff, phenry, pf0, pdensity, pdquer, &
       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
78
79
80
81
82
83
84

  pspecies=" "
  pdecay=-999.9
  pweta=-9.9E-09
  pwetb=0.0
  pweta_in=-9.9E-09
  pwetb_in=-9.9E-09
85
86
!  pwetc_in=-9.9E-09
!  pwetd_in=-9.9E-09
87
88
89
90
91
92
93
  preldiff=-9.9
  phenry=0.0
  pf0=0.0
  pdensity=-9.9E09
  pdquer=0.0
  pdsigma=0.0
  pdryvel=-9.99
94
95
  pohcconst=-9.99
  pohdconst=-9.9E-09
96
97
98
99
  pspec_ass=-9
  pkao=-99.99
  pweightmolar=-789.0 ! read failure indicator value

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

  ASSSPEC=.FALSE.

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

113
114
115
116
117
118
119
120
121
  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found

    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
122
123

    read(unitspecies,'(a10)',end=22) species(pos_spec)
124
!  write(*,*) species(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
125
    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
126
!  write(*,*) decay(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
127
    read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
128
!  write(*,*) weta(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
129
    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
130
131
132
133
134
135
136
137
138
139
140
!  write(*,*) wetb(pos_spec)

!*** NIK 31.01.2013: including in-cloud scavening parameters
    read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
!  write(*,*) weta_in(pos_spec)
    read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
!  write(*,*) wetb_in(pos_spec)
! read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
!  write(*,*) wetc_in(pos_spec)
! read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
!  write(*,*) wetd_in(pos_spec)
141

Matthias Langer's avatar
 
Matthias Langer committed
142
    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
143
!  write(*,*) reldiff(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
144
    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
145
!  write(*,*) henry(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
146
    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
147
!  write(*,*) f0(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
148
    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
149
!  write(*,*) density(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
150
    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
151
!  write(*,*) dquer(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
152
    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
153
!  write(*,*) dsigma(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
154
    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
155
!  write(*,*) dryvel(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
156
    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
157
158
159
160
161
!  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)
Matthias Langer's avatar
 
Matthias Langer committed
162
    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
163
!  write(*,*) spec_ass(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
164
    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
165
!       write(*,*) kao(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
166

167
168
169
170
171
172
    pspecies=species(pos_spec)
    pdecay=decay(pos_spec)
    pweta=weta(pos_spec)
    pwetb=wetb(pos_spec)
    pweta_in=weta_in(pos_spec)
    pwetb_in=wetb_in(pos_spec)
173
174
!    pwetc_in=wetc_in(pos_spec)
!    pwetd_in=wetd_in(pos_spec)
175
176
177
178
179
180
181
182
    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)
183
184
    pohcconst=ohcconst(pos_spec)
    pohdconst=ohdconst(pos_spec)
185
186
    pspec_ass=spec_ass(pos_spec)
    pkao=kao(pos_spec)
Matthias Langer's avatar
 
Matthias Langer committed
187

188
189
190
191
192
193
194
195
  else

    species(pos_spec)=pspecies
    decay(pos_spec)=pdecay
    weta(pos_spec)=pweta
    wetb(pos_spec)=pwetb
    weta_in(pos_spec)=pweta_in
    wetb_in(pos_spec)=pwetb_in
196
197
!    wetc_in(pos_spec)=pwetc_in
!    wetd_in(pos_spec)=pwetd_in
198
199
200
201
202
203
204
205
    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
206
207
    ohcconst(pos_spec)=pohcconst
    ohdconst(pos_spec)=pohdconst
208
209
210
211
    spec_ass(pos_spec)=pspec_ass
    kao(pos_spec)=pkao

  endif
Matthias Langer's avatar
 
Matthias Langer committed
212

213
  i=pos_spec
Matthias Langer's avatar
 
Matthias Langer committed
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
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
  if (dquer(pos_spec).gt.0) then !is particle
    write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
         &(Rain collection efficiency)  ', weta(pos_spec)
    write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
         &(Snow collection efficiency)  ', wetb(pos_spec)
    if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
      write(*,*) '*******************************************'
      write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
           &is out of likely range'
      write(*,*) '          Likely range is 0.0-1.0'
      write(*,*) '*******************************************'
    endif
    write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
         &(CCN efficiency)  ', weta_in(pos_spec)
    write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
         &(IN efficiency)  ', wetb_in(pos_spec)
    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7 )then
      write(*,*) '*******************************************'
      write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
      write(*,*) '          The default value is 0.9 for CCN '
      write(*,*) '*******************************************'
    endif
    if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
      write(*,*) '*******************************************'
      write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
      write(*,*) '          The default value is 0.1 for IN '
      write(*,*) '*******************************************'
    endif

  else !is gas
    write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
    write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
    write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
    if (wetb(pos_spec).gt.0.) then !if wet deposition is turned on
      if (wetb(pos_spec).gt.0.8 .or. wetb(pos_spec).lt.0.6) then
        write(*,*) '*******************************************'
        write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
        write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
        write(*,*) '*******************************************'
      endif
    end if
    if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
      if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
        write(*,*) '*******************************************'
        write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
        write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
        write(*,*) '*******************************************'
      endif
    end if
  endif


269
  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
270
    if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
271
272
273
274
275
276
277
278
279
280
281
282
283
  endif

  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
    if (spec_found.eqv..FALSE.) then
      goto 997
Matthias Langer's avatar
 
Matthias Langer committed
284
    endif
285
286
287
288
289
290
291
292
293
294
295
296
  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
297
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
298
299


300
301
302
! 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
303

304
305
306
307
308
309
310
311
312
313
  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
314
315
316
317
318
319
320
321
322
323

    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

324
325
326
  endif

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

328
329
330
! 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)
331
332
333
334
335
    write(unitspecies,nml=species_params)
    close(unitspecies)
  endif

  return
Matthias Langer's avatar
 
Matthias Langer committed
336

337
996 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
338
339
340
341
342
343
344
345
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
  write(*,*) '#### CONSTANT IS SET                            ####'
  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
  write(*,*) '#####################################################'
  stop


346
997 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
347
348
349
350
351
352
353
354
355
  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


356
998 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
357
358
359
360
361
362
363
  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

364
365
366
367
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
368
369

end subroutine readspecies