source: flexpart.git/src/readspecies.f90 @ 92fab65

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

add SPDX-License-Identifier to all .f90 files

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