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
RevLine 
[e200b7a]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                                                             *
[8a65cb0]11  !                                                                            *
[4fbe7a5]12  !   Changes:                                                                 *
13  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
14  !                                                                            *
[b4d29ce]15  !   HSO, 13 August 2013
16  !   added optional namelist input
17  !                                                                            *
[e200b7a]18  !*****************************************************************************
19  !                                                                            *
20  ! Variables:                                                                 *
[05cf28d]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          *
[e200b7a]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
[b4d29ce]46  character(len=16) :: pspecies
[341f4b7]47  real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
[e9e0f06]48  real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst
[341f4b7]49  real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
[d6a245b]50  real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
[e9e0f06]51  integer :: readerror
[b4d29ce]52
[8a65cb0]53! declare namelist
[b4d29ce]54  namelist /species_params/ &
[341f4b7]55       pspecies, pdecay, pweta_gas, pwetb_gas, &
56       pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
[8a65cb0]57       preldiff, phenry, pf0, pdensity, pdquer, &
[d6a245b]58       pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pohnconst, &
59       parea_dow, parea_hour, ppoint_dow, ppoint_hour
[b4d29ce]60
[341f4b7]61  pspecies="" ! read failure indicator value
[b4d29ce]62  pdecay=-999.9
[341f4b7]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
[b4d29ce]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
[8a65cb0]76  pohcconst=-9.99
77  pohdconst=-9.9E-09
[78e62dc]78  pohnconst=2.0
[341f4b7]79  pweightmolar=-999.9
[46f93d5]80  parea_dow=-999.9
81  parea_hour=-999.9
82  ppoint_dow=-999.9
83  ppoint_hour=-999.9
84
[b4d29ce]85
[d6a245b]86  do j=1,24           ! initialize everything to no variation
[6a713e8]87    parea_hour(j)=1.
88    ppoint_hour(j)=1.
[d6a245b]89    area_hour(pos_spec,j)=1.
90    point_hour(pos_spec,j)=1.
91  end do
92  do j=1,7
[6a713e8]93    parea_dow(j)=1.
94    ppoint_dow(j)=1.
[d6a245b]95    area_dow(pos_spec,j)=1.
96    point_dow(pos_spec,j)=1.
97  end do
98
[8a65cb0]99! Open the SPECIES file and read species names and properties
100!************************************************************
[e200b7a]101  specnum(pos_spec)=id_spec
102  write(aspecnumb,'(i3.3)') specnum(pos_spec)
[b4d29ce]103  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
[8a65cb0]104!write(*,*) 'reading SPECIES',specnum(pos_spec)
[e200b7a]105
106  ASSSPEC=.FALSE.
107
[8a65cb0]108! try namelist input
[b4d29ce]109  read(unitspecies,species_params,iostat=readerror)
110  close(unitspecies)
[8a65cb0]111
[05cf28d]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"
[b4d29ce]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
[e200b7a]123
124    read(unitspecies,'(a10)',end=22) species(pos_spec)
[8a65cb0]125!  write(*,*) species(pos_spec)
[e200b7a]126    read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
[8a65cb0]127!  write(*,*) decay(pos_spec)
[341f4b7]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)
[8a65cb0]136!*** NIK 31.01.2013: including in-cloud scavening parameters
[341f4b7]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)
[e200b7a]141    read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
[8a65cb0]142!  write(*,*) reldiff(pos_spec)
[e200b7a]143    read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
[8a65cb0]144!  write(*,*) henry(pos_spec)
[e200b7a]145    read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
[8a65cb0]146!  write(*,*) f0(pos_spec)
[e200b7a]147    read(unitspecies,'(e18.1)',end=22) density(pos_spec)
[8a65cb0]148!  write(*,*) density(pos_spec)
[e200b7a]149    read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
[05cf28d]150!  write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
[e200b7a]151    read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
[8a65cb0]152!  write(*,*) dsigma(pos_spec)
[e200b7a]153    read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
[8a65cb0]154!  write(*,*) dryvel(pos_spec)
[e200b7a]155    read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
[8a65cb0]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)
[78e62dc]161    read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
162!  write(*,*) ohnconst(pos_spec)
[e200b7a]163
[d6a245b]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
[b4d29ce]176    pspecies=species(pos_spec)
177    pdecay=decay(pos_spec)
[341f4b7]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)
[b4d29ce]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)
[8a65cb0]192    pohcconst=ohcconst(pos_spec)
193    pohdconst=ohdconst(pos_spec)
[78e62dc]194    pohnconst=ohnconst(pos_spec)
[e200b7a]195
[d6a245b]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
[b4d29ce]207
208    species(pos_spec)=pspecies
209    decay(pos_spec)=pdecay
[341f4b7]210    weta_gas(pos_spec)=pweta_gas
211    wetb_gas(pos_spec)=pwetb_gas
[c8fc724]212    crain_aero(pos_spec)=pcrain_aero
213    csnow_aero(pos_spec)=pcsnow_aero
[341f4b7]214    ccn_aero(pos_spec)=pccn_aero
215    in_aero(pos_spec)=pin_aero
[b4d29ce]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
[8a65cb0]224    ohcconst(pos_spec)=pohcconst
225    ohdconst(pos_spec)=pohdconst
[78e62dc]226    ohnconst(pos_spec)=pohnconst
[d6a245b]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
[b4d29ce]237  endif
[e200b7a]238
[b4d29ce]239  i=pos_spec
[e200b7a]240
[8a65cb0]241!NIK 16.02.2015
242! Check scavenging parameters given in SPECIES file
[9484483]243
[05cf28d]244  if (lroot) then
[62e65c7]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
[05cf28d]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'
[5f9d14a]274      end if
[05cf28d]275      if (density(pos_spec) .gt. 0) then
276        write(*,'(a)') '  Dry deposition is turned         :   ON'
[07c3e71]277        if (reldiff(pos_spec).gt.0) then
278           stop 'density>0 (SPECIES is a particle) implies reldiff <=0  '
279        endif
[05cf28d]280      else
[07c3e71]281        if (reldiff(pos_spec).le.0) then
282           stop 'density<=0 (SPECIES is a gas) implies reldiff >0  '
283        endif     
[05cf28d]284        write(*,'(a)') '  Dry deposition is (density<0)    :   OFF'
[5f9d14a]285      end if
[05cf28d]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
[5f9d14a]288        write(*,*) '*******************************************'
[05cf28d]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'
[5f9d14a]294        write(*,*) '*******************************************'
295      end if
[05cf28d]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
[5f9d14a]316          write(*,*) '*******************************************'
[05cf28d]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 '
[5f9d14a]320          write(*,*) '*******************************************'
321        end if
322      endif
[05cf28d]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
[8a65cb0]327      endif
328    end if
[05cf28d]329  end if
[b4d29ce]330
[07c3e71]331  !  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
[f04992e]332  if (dquer(i).gt.0 .and. dsigma(i).le.1.) then !dsigma(i)=1.0001   ! avoid floating exception
[07c3e71]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
[b4d29ce]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
[8a65cb0]35420 continue
[e200b7a]355
356
[b4d29ce]35722 close(unitspecies)
[e200b7a]358
[8a65cb0]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)
[b4d29ce]362    write(unitspecies,nml=species_params)
363    close(unitspecies)
364  endif
365
366  return
[e200b7a]367
[8a65cb0]368996 write(*,*) '#####################################################'
[e200b7a]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
[8a65cb0]377997 write(*,*) '#####################################################'
[e200b7a]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
[8a65cb0]387998 write(*,*) '#####################################################'
[e200b7a]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
[b4d29ce]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
[e200b7a]399
400end subroutine readspecies
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG