source: flexpart.git/src/readspecies.f90 @ 46f93d5

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 46f93d5 was 46f93d5, checked in by Ignacio Pisso <ip@…>, 5 years ago

default values for point, area for dow, hour

  • Property mode set to 100644
File size: 16.1 KB
Line 
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                                                             *
32  !                                                                            *
33  !   Changes:                                                                 *
34  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
35  !                                                                            *
36  !   HSO, 13 August 2013
37  !   added optional namelist input
38  !                                                                            *
39  !*****************************************************************************
40  !                                                                            *
41  ! Variables:                                                                 *
42  ! decaytime(maxtable)   half time for radiological decay                     *
43  ! specname(maxtable)    names of chemical species, radionuclides             *
44  ! weta_gas, wetb_gas    Parameters for below-cloud scavenging of gasses      *
45  ! crain_aero,csnow_aero Parameters for below-cloud scavenging of aerosols    *
46  ! ccn_aero,in_aero      Parameters for in-cloud scavenging of aerosols       *
47  ! ohcconst              OH reaction rate constant C                          *
48  ! ohdconst              OH reaction rate constant D                          *
49  ! ohnconst              OH reaction rate constant n                          *
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
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
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       pspecies, pdecay, pweta_gas, pwetb_gas, &
77       pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
78       preldiff, phenry, pf0, pdensity, pdquer, &
79       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
80       parea_dow, parea_hour, ppoint_dow, ppoint_hour
81
82  pspecies="" ! read failure indicator value
83  pdecay=-999.9
84  pweta_gas=-9.9E-09
85  pwetb_gas=0.0
86  pcrain_aero=-9.9E-09
87  pcsnow_aero=-9.9E-09
88  pccn_aero=-9.9E-09
89  pin_aero=-9.9E-09
90  preldiff=-9.9
91  phenry=0.0
92  pf0=0.0
93  pdensity=-9.9E09
94  pdquer=0.0
95  pdsigma=0.0
96  pdryvel=-9.99
97  pohcconst=-9.99
98  pohdconst=-9.9E-09
99  pohnconst=2.0
100  pweightmolar=-999.9
101  parea_dow=-999.9
102  parea_hour=-999.9
103  ppoint_dow=-999.9
104  ppoint_hour=-999.9
105
106
107  do j=1,24           ! initialize everything to no variation
108    parea_hour(j)=1.
109    ppoint_hour(j)=1.
110    area_hour(pos_spec,j)=1.
111    point_hour(pos_spec,j)=1.
112  end do
113  do j=1,7
114    parea_dow(j)=1.
115    ppoint_dow(j)=1.
116    area_dow(pos_spec,j)=1.
117    point_dow(pos_spec,j)=1.
118  end do
119
120! Open the SPECIES file and read species names and properties
121!************************************************************
122  specnum(pos_spec)=id_spec
123  write(aspecnumb,'(i3.3)') specnum(pos_spec)
124  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
125!write(*,*) 'reading SPECIES',specnum(pos_spec)
126
127  ASSSPEC=.FALSE.
128
129! try namelist input
130  read(unitspecies,species_params,iostat=readerror)
131  close(unitspecies)
132
133  if ((len(trim(pspecies)).eq.0).or.(readerror.ne.0)) then ! no namelist found
134    if (lroot) write(*,*) "SPECIES file not in NAMELIST format, attempting to &
135         &read as fixed format"
136
137    readerror=1
138
139    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
140
141    do i=1,6
142      read(unitspecies,*)
143    end do
144
145    read(unitspecies,'(a10)',end=22) species(pos_spec)
146!  write(*,*) species(pos_spec)
147    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
148!  write(*,*) decay(pos_spec)
149    read(unitspecies,'(e18.1)',end=22) weta_gas(pos_spec)
150!  write(*,*) weta_gas(pos_spec)
151    read(unitspecies,'(f18.2)',end=22) wetb_gas(pos_spec)
152!  write(*,*) wetb_gas(pos_spec)
153    read(unitspecies,'(e18.1)',end=22) crain_aero(pos_spec)
154!  write(*,*) crain_aero(pos_spec)
155    read(unitspecies,'(f18.2)',end=22) csnow_aero(pos_spec)
156!  write(*,*) csnow_aero(pos_spec)
157!*** NIK 31.01.2013: including in-cloud scavening parameters
158    read(unitspecies,'(e18.1)',end=22) ccn_aero(pos_spec)
159!  write(*,*) ccn_aero(pos_spec)
160    read(unitspecies,'(f18.2)',end=22) in_aero(pos_spec)
161!  write(*,*) in_aero(pos_spec)
162    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
163!  write(*,*) reldiff(pos_spec)
164    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
165!  write(*,*) henry(pos_spec)
166    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
167!  write(*,*) f0(pos_spec)
168    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
169!  write(*,*) density(pos_spec)
170    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
171!  write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
172    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
173!  write(*,*) dsigma(pos_spec)
174    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
175!  write(*,*) dryvel(pos_spec)
176    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
177!  write(*,*) weightmolar(pos_spec)
178    read(unitspecies,'(e18.2)',end=22) ohcconst(pos_spec)
179!  write(*,*) ohcconst(pos_spec)
180    read(unitspecies,'(f8.2)',end=22) ohdconst(pos_spec)
181!  write(*,*) ohdconst(pos_spec)
182    read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
183!  write(*,*) ohnconst(pos_spec)
184
185! Read in daily and day-of-week variation of emissions, if available
186!*******************************************************************
187
188    read(unitspecies,*,end=22)
189    do j=1,24     ! 24 hours, starting with 0-1 local time
190      read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
191    end do
192    read(unitspecies,*)
193    do j=1,7      ! 7 days of the week, starting with Monday
194      read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
195    end do
196
197    pspecies=species(pos_spec)
198    pdecay=decay(pos_spec)
199    pweta_gas=weta_gas(pos_spec)
200    pwetb_gas=wetb_gas(pos_spec)
201    pcrain_aero=crain_aero(pos_spec)
202    pcsnow_aero=csnow_aero(pos_spec)
203    pccn_aero=ccn_aero(pos_spec)
204    pin_aero=in_aero(pos_spec)
205    preldiff=reldiff(pos_spec)
206    phenry=henry(pos_spec)
207    pf0=f0(pos_spec)
208    pdensity=density(pos_spec)
209    pdquer=dquer(pos_spec)
210    pdsigma=dsigma(pos_spec)
211    pdryvel=dryvel(pos_spec)
212    pweightmolar=weightmolar(pos_spec)
213    pohcconst=ohcconst(pos_spec)
214    pohdconst=ohdconst(pos_spec)
215    pohnconst=ohnconst(pos_spec)
216
217
218    do j=1,24     ! 24 hours, starting with 0-1 local time
219      parea_hour(j)=area_hour(pos_spec,j)
220      ppoint_hour(j)=point_hour(pos_spec,j)
221    end do
222    do j=1,7      ! 7 days of the week, starting with Monday
223      parea_dow(j)=area_dow(pos_spec,j)
224      ppoint_dow(j)=point_dow(pos_spec,j)
225    end do
226
227  else ! namelist available
228
229    species(pos_spec)=pspecies
230    decay(pos_spec)=pdecay
231    weta_gas(pos_spec)=pweta_gas
232    wetb_gas(pos_spec)=pwetb_gas
233    crain_aero(pos_spec)=pcrain_aero
234    csnow_aero(pos_spec)=pcsnow_aero
235    ccn_aero(pos_spec)=pccn_aero
236    in_aero(pos_spec)=pin_aero
237    reldiff(pos_spec)=preldiff
238    henry(pos_spec)=phenry
239    f0(pos_spec)=pf0
240    density(pos_spec)=pdensity
241    dquer(pos_spec)=pdquer
242    dsigma(pos_spec)=pdsigma
243    dryvel(pos_spec)=pdryvel
244    weightmolar(pos_spec)=pweightmolar
245    ohcconst(pos_spec)=pohcconst
246    ohdconst(pos_spec)=pohdconst
247    ohnconst(pos_spec)=pohnconst
248
249    do j=1,24     ! 24 hours, starting with 0-1 local time
250      area_hour(pos_spec,j)=parea_hour(j)
251      point_hour(pos_spec,j)=ppoint_hour(j)
252    end do
253    do j=1,7      ! 7 days of the week, starting with Monday
254      area_dow(pos_spec,j)=parea_dow(j)
255      point_dow(pos_spec,j)=ppoint_dow(j)
256    end do
257
258  endif
259
260  i=pos_spec
261
262!NIK 16.02.2015
263! Check scavenging parameters given in SPECIES file
264
265  if (lroot) then
266! ZHG 2016.04.07 Start of changes
267    write(*,*) ' '
268    if (dquer(pos_spec) .gt.0)  write(*,'(a,i3,a,a,a)')       ' SPECIES: ', &
269         &id_spec,'  ', species(pos_spec),'  (AEROSOL) '
270    if (dquer(pos_spec) .le.0)  write(*,'(a,i3,a,a,a)')       ' SPECIES: ', &
271         &id_spec,'  ', species(pos_spec),'  (GAS) '
272
273! Particles
274!**********
275    if (dquer(pos_spec).gt.0) then
276      if (ccn_aero(pos_spec) .gt. 0) then
277        write(*,'(a,f5.2)') '  Particle CCN  efficiency (CCNeff):', ccn_aero(pos_spec)
278      else
279        write(*,'(a)')      '  Particle CCN  efficiency (CCNeff):   OFF'
280      endif
281      if (in_aero(pos_spec) .gt. 0) then
282        write(*,'(a,f5.2)') '  Particle  IN  efficiency (INeff) :', in_aero(pos_spec)
283      else
284        write(*,'(a)')      '  Particle  IN  efficiency (INeff) :   OFF'
285      endif
286      if (crain_aero(pos_spec) .gt. 0) then
287        write(*,'(a,f5.2)') '  Particle Rain efficiency (Crain) :', crain_aero(pos_spec)
288      else
289        write(*,'(a)')      '  Particle Rain efficiency (Crain) :   OFF'
290      endif
291      if (csnow_aero(pos_spec) .gt. 0) then
292        write(*,'(a,f5.2)') '  Particle Snow efficiency (Csnow) :', csnow_aero(pos_spec)
293      else
294        write(*,'(a)')      '  Particle Snow efficiency (Csnow) :   OFF'
295      end if
296      if (density(pos_spec) .gt. 0) then
297        write(*,'(a)') '  Dry deposition is turned         :   ON'
298      else
299        write(*,'(a)') '  Dry deposition is (density<0)    :   OFF'
300      end if
301      if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. &
302           &ccn_aero(pos_spec).gt.1.0 .or. in_aero(pos_spec).gt.1.0) then
303        write(*,*) '*******************************************'
304        write(*,*) ' WARNING: Particle Scavenging parameter likely out of range '
305        write(*,*) '       Likely   range for Crain    0.0-10'
306        write(*,*) '       Likely   range for Csnow    0.0-10'
307        write(*,*) '       Physical range for CCNeff   0.0-1'
308        write(*,*) '       Physical range for INeff    0.0-1'
309        write(*,*) '*******************************************'
310      end if
311    else
312! Gas
313!****
314      if (weta_gas(pos_spec) .gt. 0 .and. wetb_gas(pos_spec).gt.0) then
315        write(*,*)          '  Wet removal for gases      is turned: ON'
316        write(*,*)          '  Gas below-cloud scavenging parameter A  ', &
317             &weta_gas(pos_spec)
318        write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', &
319             &wetb_gas(pos_spec)
320      else
321        write(*,*)          '  Wet removal for gases      is turned: OFF '
322      end if
323      if (reldiff(i).gt.0.) then
324        write(*,*)          '  Dry deposition for gases   is turned: ON '
325      else
326        write(*,*)          '  Dry deposition for gases   is turned: OFF '
327      end if
328      if (weta_gas(pos_spec).gt.0.) then !if wet deposition is turned on
329        if (weta_gas(pos_spec).gt.1E-04 .or. weta_gas(pos_spec).lt.1E-09 .or. &
330             &wetb_gas(pos_spec).gt.0.8 .or. wetb_gas(pos_spec).lt.0.4) then
331          write(*,*) '*******************************************'
332          write(*,*) ' WARNING: Gas below-cloud scavengig is out of likely range'
333          write(*,*) '          Likely range for A is 1E-04 to 1E-08'
334          write(*,*) '          Likely range for B is 0.60  to 0.80 '
335          write(*,*) '*******************************************'
336        end if
337      endif
338
339      if (((weta_gas(pos_spec).gt.0).or.(wetb_gas(pos_spec).gt.0)).and.&
340           &(henry(pos_spec).le.0)) then
341        if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
342      endif
343    end if
344  end if
345
346  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
347
348  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
349    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
350    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
351    write(*,*) '#### PARTICLE AND GAS.                       ####'
352    write(*,*) '#### SPECIES NUMBER',aspecnumb
353    stop
354  endif
35520 continue
356
357
35822 close(unitspecies)
359
360! namelist output if requested
361  if (nmlout.and.lroot) then
362    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000)
363    write(unitspecies,nml=species_params)
364    close(unitspecies)
365  endif
366
367  return
368
369996 write(*,*) '#####################################################'
370  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
371  write(*,*) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS  #### '
372  write(*,*) '#### CONSTANT IS SET                            ####'
373  write(*,*) '#### PLEASE MODIFY SPECIES DESCR. FILE!        #### '
374  write(*,*) '#####################################################'
375  stop
376
377
378997 write(*,*) '#####################################################'
379  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
380  write(*,*) '#### THE ASSSOCIATED SPECIES HAS TO BE DEFINED  #### '
381  write(*,*) '#### BEFORE THE ONE WHICH POINTS AT IT          #### '
382  write(*,*) '#### PLEASE CHANGE ORDER IN RELEASES OR ADD     #### '
383  write(*,*) '#### THE ASSOCIATED SPECIES IN RELEASES         #### '
384  write(*,*) '#####################################################'
385  stop
386
387
388998 write(*,*) '#####################################################'
389  write(*,*) '#### FLEXPART MODEL ERROR!                      #### '
390  write(*,*) '#### THE SPECIES FILE FOR SPECIES ', id_spec
391  write(*,*) '#### CANNOT BE FOUND: CREATE FILE'
392  write(*,*) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb
393  write(*,*) '#####################################################'
394  stop
395
3961000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
397  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
398  write(*,'(a)') path(2)(1:length(2))
399  stop
400
401end subroutine readspecies
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG