source: flexpart.git/src/readspecies.f90 @ 3481cc1

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

move license from headers to a different file

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