readspecies.f90 15.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
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
!NIK 16.02.2015
! Check scavenging parameters given in SPECIES file
  if (dquer(pos_spec).gt.0) then !is particle
Espen Sollum's avatar
Espen Sollum committed
218
219
220
221
222
223
    if (lroot) then
      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)
    end if
224
    if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
Espen Sollum's avatar
Espen Sollum committed
225
226
227
228
229
230
231
      if (lroot) 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(*,*) '*******************************************'
      end if
232
    endif
Espen Sollum's avatar
Espen Sollum committed
233
234
235
236
237
238
239
240
241
242
243
244
245
    if (lroot) then
      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)
    end if
    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7) then
      if (lroot) 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(*,*) '*******************************************'
      end if
246
247
    endif
    if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
Espen Sollum's avatar
Espen Sollum committed
248
249
250
251
252
253
      if (lroot) 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(*,*) '*******************************************'
      end if
254
255
256
    endif

  else !is gas
Espen Sollum's avatar
Espen Sollum committed
257
258
259
260
    if (lroot) then
      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'
261
262
263
    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
Espen Sollum's avatar
Espen Sollum committed
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
        if (lroot) 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(*,*) '*******************************************'
        end if
      endif
    end if
    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
        if (lroot) 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(*,*) '*******************************************'
        end if
280
281
282
283
284
      endif
    end if
  endif


285
  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
286
    if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
287
288
289
290
291
292
293
294
295
296
297
298
299
  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
300
    endif
301
302
303
304
305
306
307
308
309
310
311
312
  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
313
20 continue
Matthias Langer's avatar
 
Matthias Langer committed
314
315


316
317
318
! 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
319

320
321
322
323
324
325
326
327
328
329
  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
330
331
332
333
334
335
336
337
338
339

    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

340
341
342
  endif

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

344
345
346
! 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)
347
348
349
350
351
    write(unitspecies,nml=species_params)
    close(unitspecies)
  endif

  return
Matthias Langer's avatar
 
Matthias Langer committed
352

353
996 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
354
355
356
357
358
359
360
361
  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
  write(*,*) '#### CONSTANT IS SET                            ####'
  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
  write(*,*) '#####################################################'
  stop


362
997 write(*,*) '#####################################################'
Matthias Langer's avatar
 
Matthias Langer committed
363
364
365
366
367
368
369
370
371
  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


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

380
381
382
383
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
384
385

end subroutine readspecies