source: flexpart.git/src/readspecies.f90 @ 26f6039

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 26f6039 was 0f20c31, checked in by flexpart <>, 9 years ago

Updated SPECIES files
  • Property mode set to 100644
File size: 15.3 KB
RevLine 
[e200b7a]1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine readspecies(id_spec,pos_spec)
23
24  !*****************************************************************************
25  !                                                                            *
26  !     This routine reads names and physical constants of chemical species/   *
27  !     radionuclides given in the parameter pos_spec                          *
28  !                                                                            *
29  !   Author: A. Stohl                                                         *
30  !                                                                            *
31  !   11 July 1996                                                             *
[8a65cb0]32  !                                                                            *
[4fbe7a5]33  !   Changes:                                                                 *
34  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
35  !                                                                            *
[b4d29ce]36  !   HSO, 13 August 2013
37  !   added optional namelist input
38  !                                                                            *
[e200b7a]39  !*****************************************************************************
40  !                                                                            *
41  ! Variables:                                                                 *
42  ! decaytime(maxtable)  half time for radiological decay                      *
43  ! specname(maxtable)   names of chemical species, radionuclides              *
[4fbe7a5]44  ! weta, wetb           Parameters for determining below-cloud scavenging     *
45  ! weta_in              Parameter for determining in-cloud scavenging         *
46  ! wetb_in              Parameter for determining in-cloud scavenging         *
[8a65cb0]47  ! ohcconst             OH reaction rate constant C                           *
48  ! ohdconst             OH reaction rate constant D                           *
[78e62dc]49  ! ohnconst             OH reaction rate constant n                           *
[e200b7a]50  ! id_spec              SPECIES number as referenced in RELEASE file          *
51  ! id_pos               position where SPECIES data shall be stored           *
52  !                                                                            *
53  ! Constants:                                                                 *
54  !                                                                            *
55  !*****************************************************************************
56
57  use par_mod
58  use com_mod
59
60  implicit none
61
62  integer :: i, pos_spec,j
63  integer :: idow,ihour,id_spec
64  character(len=3) :: aspecnumb
65  logical :: spec_found
66
[b4d29ce]67  character(len=16) :: pspecies
68  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
[78e62dc]69  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
[b4d29ce]70  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
71  integer :: readerror
72
[8a65cb0]73! declare namelist
[b4d29ce]74  namelist /species_params/ &
[8a65cb0]75       pspecies, pdecay, pweta, pwetb, &
76       pweta_in, pwetb_in, pwetc_in, pwetd_in, &
77       preldiff, phenry, pf0, pdensity, pdquer, &
[78e62dc]78       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, pspec_ass, pkao
[b4d29ce]79
80  pspecies=" "
81  pdecay=-999.9
82  pweta=-9.9E-09
83  pwetb=0.0
84  pweta_in=-9.9E-09
85  pwetb_in=-9.9E-09
[8a65cb0]86!  pwetc_in=-9.9E-09
87!  pwetd_in=-9.9E-09
[b4d29ce]88  preldiff=-9.9
89  phenry=0.0
90  pf0=0.0
91  pdensity=-9.9E09
92  pdquer=0.0
93  pdsigma=0.0
94  pdryvel=-9.99
[8a65cb0]95  pohcconst=-9.99
96  pohdconst=-9.9E-09
[78e62dc]97  pohnconst=2.0
[b4d29ce]98  pspec_ass=-9
99  pkao=-99.99
100  pweightmolar=-789.0 ! read failure indicator value
101
[8a65cb0]102! Open the SPECIES file and read species names and properties
103!************************************************************
[e200b7a]104  specnum(pos_spec)=id_spec
105  write(aspecnumb,'(i3.3)') specnum(pos_spec)
[b4d29ce]106  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
[8a65cb0]107!write(*,*) 'reading SPECIES',specnum(pos_spec)
[e200b7a]108
109  ASSSPEC=.FALSE.
110
[8a65cb0]111! try namelist input
[b4d29ce]112  read(unitspecies,species_params,iostat=readerror)
113  close(unitspecies)
[8a65cb0]114
[b4d29ce]115  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
116
117    readerror=1
118
119    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
120
121    do i=1,6
122      read(unitspecies,*)
123    end do
[e200b7a]124
125    read(unitspecies,'(a10)',end=22) species(pos_spec)
[8a65cb0]126!  write(*,*) species(pos_spec)
[e200b7a]127    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
[8a65cb0]128!  write(*,*) decay(pos_spec)
[e200b7a]129    read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
[8a65cb0]130!  write(*,*) weta(pos_spec)
[e200b7a]131    read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
[8a65cb0]132!  write(*,*) wetb(pos_spec)
133
134!*** NIK 31.01.2013: including in-cloud scavening parameters
135    read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
136!  write(*,*) weta_in(pos_spec)
137    read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
138!  write(*,*) wetb_in(pos_spec)
139! read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
140!  write(*,*) wetc_in(pos_spec)
141! read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
142!  write(*,*) wetd_in(pos_spec)
[4fbe7a5]143
[e200b7a]144    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
[8a65cb0]145!  write(*,*) reldiff(pos_spec)
[e200b7a]146    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
[8a65cb0]147!  write(*,*) henry(pos_spec)
[e200b7a]148    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
[8a65cb0]149!  write(*,*) f0(pos_spec)
[e200b7a]150    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
[8a65cb0]151!  write(*,*) density(pos_spec)
[e200b7a]152    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
[8a65cb0]153!  write(*,*) dquer(pos_spec)
[e200b7a]154    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
[8a65cb0]155!  write(*,*) dsigma(pos_spec)
[e200b7a]156    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
[8a65cb0]157!  write(*,*) dryvel(pos_spec)
[e200b7a]158    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
[8a65cb0]159!  write(*,*) weightmolar(pos_spec)
160    read(unitspecies,'(e18.2)',end=22) ohcconst(pos_spec)
161!  write(*,*) ohcconst(pos_spec)
162    read(unitspecies,'(f8.2)',end=22) ohdconst(pos_spec)
163!  write(*,*) ohdconst(pos_spec)
[78e62dc]164    read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
165!  write(*,*) ohnconst(pos_spec)
[e200b7a]166    read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
[8a65cb0]167!  write(*,*) spec_ass(pos_spec)
[e200b7a]168    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
[8a65cb0]169!       write(*,*) kao(pos_spec)
[e200b7a]170
[b4d29ce]171    pspecies=species(pos_spec)
172    pdecay=decay(pos_spec)
173    pweta=weta(pos_spec)
174    pwetb=wetb(pos_spec)
175    pweta_in=weta_in(pos_spec)
176    pwetb_in=wetb_in(pos_spec)
[8a65cb0]177!    pwetc_in=wetc_in(pos_spec)
178!    pwetd_in=wetd_in(pos_spec)
[b4d29ce]179    preldiff=reldiff(pos_spec)
180    phenry=henry(pos_spec)
181    pf0=f0(pos_spec)
182    pdensity=density(pos_spec)
183    pdquer=dquer(pos_spec)
184    pdsigma=dsigma(pos_spec)
185    pdryvel=dryvel(pos_spec)
186    pweightmolar=weightmolar(pos_spec)
[8a65cb0]187    pohcconst=ohcconst(pos_spec)
188    pohdconst=ohdconst(pos_spec)
[78e62dc]189    pohnconst=ohnconst(pos_spec)
[b4d29ce]190    pspec_ass=spec_ass(pos_spec)
191    pkao=kao(pos_spec)
[e200b7a]192
[b4d29ce]193  else
194
195    species(pos_spec)=pspecies
196    decay(pos_spec)=pdecay
197    weta(pos_spec)=pweta
198    wetb(pos_spec)=pwetb
199    weta_in(pos_spec)=pweta_in
200    wetb_in(pos_spec)=pwetb_in
[8a65cb0]201!    wetc_in(pos_spec)=pwetc_in
202!    wetd_in(pos_spec)=pwetd_in
[b4d29ce]203    reldiff(pos_spec)=preldiff
204    henry(pos_spec)=phenry
205    f0(pos_spec)=pf0
206    density(pos_spec)=pdensity
207    dquer(pos_spec)=pdquer
208    dsigma(pos_spec)=pdsigma
209    dryvel(pos_spec)=pdryvel
210    weightmolar(pos_spec)=pweightmolar
[8a65cb0]211    ohcconst(pos_spec)=pohcconst
212    ohdconst(pos_spec)=pohdconst
[78e62dc]213    ohnconst(pos_spec)=pohnconst
[b4d29ce]214    spec_ass(pos_spec)=pspec_ass
215    kao(pos_spec)=pkao
216
217  endif
[e200b7a]218
[b4d29ce]219  i=pos_spec
[e200b7a]220
[8a65cb0]221!NIK 16.02.2015
222! Check scavenging parameters given in SPECIES file
223  if (dquer(pos_spec).gt.0) then !is particle
[5f9d14a]224    if (lroot) then
225      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
226           &(Rain collection efficiency)  ', weta(pos_spec)
227      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
228           &(Snow collection efficiency)  ', wetb(pos_spec)
229    end if
[8a65cb0]230    if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
[5f9d14a]231      if (lroot) then
232        write(*,*) '*******************************************'
233        write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
234             &is out of likely range'
235        write(*,*) '          Likely range is 0.0-1.0'
236        write(*,*) '*******************************************'
237      end if
[8a65cb0]238    endif
[5f9d14a]239    if (lroot) then
240      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
241           &(CCN efficiency)  ', weta_in(pos_spec)
242      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
243           &(IN efficiency)  ', wetb_in(pos_spec)
244    end if
245    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7) then
246      if (lroot) then
247        write(*,*) '*******************************************'
248        write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
249        write(*,*) '          The default value is 0.9 for CCN '
250        write(*,*) '*******************************************'
251      end if
[8a65cb0]252    endif
253    if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
[5f9d14a]254      if (lroot) then
255        write(*,*) '*******************************************'
256        write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
257        write(*,*) '          The default value is 0.1 for IN '
258        write(*,*) '*******************************************'
259      end if
[8a65cb0]260    endif
261
262  else !is gas
[5f9d14a]263    if (lroot) then
[0f20c31]264      write(*,*) '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
[5f9d14a]265      write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
266      write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
[8a65cb0]267    end if
268    if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
269      if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
[5f9d14a]270        if (lroot) then
271          write(*,*) '*******************************************'
272          write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
273          write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
274          write(*,*) '*******************************************'
275        end if
276      endif
277    end if
278    if (wetb(pos_spec).gt.0.) then !if wet deposition is turned on
279      if (wetb(pos_spec).gt.0.8 .or. wetb(pos_spec).lt.0.6) then
280        if (lroot) then
281          write(*,*) '*******************************************'
282          write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
283          write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
284          write(*,*) '*******************************************'
285        end if
[8a65cb0]286      endif
287    end if
288  endif
289
290
[b4d29ce]291  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
[8a65cb0]292    if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
[b4d29ce]293  endif
294
295  if (spec_ass(pos_spec).gt.0) then
296    spec_found=.FALSE.
297    do j=1,pos_spec-1
298      if (spec_ass(pos_spec).eq.specnum(j)) then
299        spec_ass(pos_spec)=j
300        spec_found=.TRUE.
301        ASSSPEC=.TRUE.
302      endif
303    end do
304    if (spec_found.eqv..FALSE.) then
305      goto 997
[e200b7a]306    endif
[b4d29ce]307  endif
308
309  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
310  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
311
312  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
313    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
314    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
315    write(*,*) '#### PARTICLE AND GAS.                       ####'
316    write(*,*) '#### SPECIES NUMBER',aspecnumb
317    stop
318  endif
[8a65cb0]31920 continue
[e200b7a]320
321
[8a65cb0]322! Read in daily and day-of-week variation of emissions, if available
323!*******************************************************************
324! HSO: This is not yet implemented as namelist parameters
[e200b7a]325
[8a65cb0]326  do j=1,24           ! initialize everything to no variation
327    area_hour(i,j)=1.
328    point_hour(i,j)=1.
329  end do
330  do j=1,7
331    area_dow(i,j)=1.
332    point_dow(i,j)=1.
333  end do
334
335  if (readerror.ne.0) then ! text format input
[e200b7a]336
337    read(unitspecies,*,end=22)
338    do j=1,24     ! 24 hours, starting with 0-1 local time
339      read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
340    end do
341    read(unitspecies,*)
342    do j=1,7      ! 7 days of the week, starting with Monday
343      read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
344    end do
345
[b4d29ce]346  endif
347
34822 close(unitspecies)
[e200b7a]349
[8a65cb0]350! namelist output if requested
351  if (nmlout.and.lroot) then
352    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000)
[b4d29ce]353    write(unitspecies,nml=species_params)
354    close(unitspecies)
355  endif
356
357  return
[e200b7a]358
[8a65cb0]359996 write(*,*) '#####################################################'
[e200b7a]360  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
361  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
362  write(*,*) '#### CONSTANT IS SET                            ####'
363  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
364  write(*,*) '#####################################################'
365  stop
366
367
[8a65cb0]368997 write(*,*) '#####################################################'
[e200b7a]369  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
370  write(*,*) '#### THE ASSSOCIATED SPECIES HAS TO BE DEFINED  #### '
371  write(*,*) '#### BEFORE THE ONE WHICH POINTS AT IT          #### '
372  write(*,*) '#### PLEASE CHANGE ORDER IN RELEASES OR ADD     #### '
373  write(*,*) '#### THE ASSOCIATED SPECIES IN RELEASES         #### '
374  write(*,*) '#####################################################'
375  stop
376
377
[8a65cb0]378998 write(*,*) '#####################################################'
[e200b7a]379  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
380  write(*,*) '#### THE SPECIES FILE FOR SPECIES ', id_spec
381  write(*,*) '#### CANNOT BE FOUND: CREATE FILE'
382  write(*,*) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb
383  write(*,*) '#####################################################'
384  stop
385
[b4d29ce]3861000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
387  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
388  write(*,'(a)') path(2)(1:length(2))
389  stop
[e200b7a]390
391end subroutine readspecies
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG