Context Navigation

readspecies.f90

bugfixes+enhancements

Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 18 months ago

Bugfixes:

options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
src/FLEXPART_MPI.f90: ditto
src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
src/gridcheck_nests.f90: ditto
src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
readwind_ecmwf_mpi.f90: ditto

Code enhancements:

options/OUTGRID: added comments describing contents
options/SPECIES/SPECIES_*: aligned comments
options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
src/FLEXPART.f90: rectified lines indenting, updated date in version string
src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
src/readreleases.f90: write out warning if too few particles are used to randomize release
src/readspecies.f90: write out name of SPECIES file before it is read
src/readwind_*.f90: write out name of input file before opening it
src/writeheader_txt.f90: removed wrong comment

Property mode set to 100644

File size: 15.8 KB

Line
1	! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2	! SPDX-License-Identifier: GPL-3.0-or-later
3
4	subroutine 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	write(,) 'Reading: '// &
107	path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb
108	open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
109	!write(,) 'reading SPECIES',specnum(pos_spec)
110
111	ASSSPEC=.FALSE.
112
113	! try namelist input
114	read(unitspecies,species_params,iostat=readerror)
115	close(unitspecies)
116
117	if ((len(trim(pspecies)).eq.0).or.(readerror.ne.0)) then ! no namelist found
118	if (lroot) write(,) "SPECIES file not in NAMELIST format, attempting to &
119	&read as fixed format"
120
121	readerror=1
122
123	open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
124
125	do i=1,6
126	read(unitspecies,*)
127	end do
128
129	read(unitspecies,'(a10)',end=22) species(pos_spec)
130	! write(,) species(pos_spec)
131	read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
132	! write(,) decay(pos_spec)
133	read(unitspecies,'(e18.1)',end=22) weta_gas(pos_spec)
134	! write(,) weta_gas(pos_spec)
135	read(unitspecies,'(f18.2)',end=22) wetb_gas(pos_spec)
136	! write(,) wetb_gas(pos_spec)
137	read(unitspecies,'(e18.1)',end=22) crain_aero(pos_spec)
138	! write(,) crain_aero(pos_spec)
139	read(unitspecies,'(f18.2)',end=22) csnow_aero(pos_spec)
140	! write(,) csnow_aero(pos_spec)
141	!*** NIK 31.01.2013: including in-cloud scavening parameters
142	read(unitspecies,'(e18.1)',end=22) ccn_aero(pos_spec)
143	! write(,) ccn_aero(pos_spec)
144	read(unitspecies,'(f18.2)',end=22) in_aero(pos_spec)
145	! write(,) in_aero(pos_spec)
146	read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
147	! write(,) reldiff(pos_spec)
148	read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
149	! write(,) henry(pos_spec)
150	read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
151	! write(,) f0(pos_spec)
152	read(unitspecies,'(e18.1)',end=22) density(pos_spec)
153	! write(,) density(pos_spec)
154	read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
155	! write(,) 'dquer(pos_spec):', dquer(pos_spec)
156	read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
157	! write(,) dsigma(pos_spec)
158	read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
159	! write(,) dryvel(pos_spec)
160	read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
161	! write(,) weightmolar(pos_spec)
162	read(unitspecies,'(e18.2)',end=22) ohcconst(pos_spec)
163	! write(,) ohcconst(pos_spec)
164	read(unitspecies,'(f8.2)',end=22) ohdconst(pos_spec)
165	! write(,) ohdconst(pos_spec)
166	read(unitspecies,'(f8.2)',end=22) ohnconst(pos_spec)
167	! write(,) ohnconst(pos_spec)
168
169	! Read in daily and day-of-week variation of emissions, if available
170	!*******************************************************************
171
172	read(unitspecies,*,end=22)
173	do j=1,24 ! 24 hours, starting with 0-1 local time
174	read(unitspecies,*) ihour,area_hour(pos_spec,j),point_hour(pos_spec,j)
175	end do
176	read(unitspecies,*)
177	do j=1,7 ! 7 days of the week, starting with Monday
178	read(unitspecies,*) idow,area_dow(pos_spec,j),point_dow(pos_spec,j)
179	end do
180
181	pspecies=species(pos_spec)
182	pdecay=decay(pos_spec)
183	pweta_gas=weta_gas(pos_spec)
184	pwetb_gas=wetb_gas(pos_spec)
185	pcrain_aero=crain_aero(pos_spec)
186	pcsnow_aero=csnow_aero(pos_spec)
187	pccn_aero=ccn_aero(pos_spec)
188	pin_aero=in_aero(pos_spec)
189	preldiff=reldiff(pos_spec)
190	phenry=henry(pos_spec)
191	pf0=f0(pos_spec)
192	pdensity=density(pos_spec)
193	pdquer=dquer(pos_spec)
194	pdsigma=dsigma(pos_spec)
195	pdryvel=dryvel(pos_spec)
196	pweightmolar=weightmolar(pos_spec)
197	pohcconst=ohcconst(pos_spec)
198	pohdconst=ohdconst(pos_spec)
199	pohnconst=ohnconst(pos_spec)
200
201
202	do j=1,24 ! 24 hours, starting with 0-1 local time
203	parea_hour(j)=area_hour(pos_spec,j)
204	ppoint_hour(j)=point_hour(pos_spec,j)
205	end do
206	do j=1,7 ! 7 days of the week, starting with Monday
207	parea_dow(j)=area_dow(pos_spec,j)
208	ppoint_dow(j)=point_dow(pos_spec,j)
209	end do
210
211	else ! namelist available
212
213	species(pos_spec)=pspecies
214	decay(pos_spec)=pdecay
215	weta_gas(pos_spec)=pweta_gas
216	wetb_gas(pos_spec)=pwetb_gas
217	crain_aero(pos_spec)=pcrain_aero
218	csnow_aero(pos_spec)=pcsnow_aero
219	ccn_aero(pos_spec)=pccn_aero
220	in_aero(pos_spec)=pin_aero
221	reldiff(pos_spec)=preldiff
222	henry(pos_spec)=phenry
223	f0(pos_spec)=pf0
224	density(pos_spec)=pdensity
225	dquer(pos_spec)=pdquer
226	dsigma(pos_spec)=pdsigma
227	dryvel(pos_spec)=pdryvel
228	weightmolar(pos_spec)=pweightmolar
229	ohcconst(pos_spec)=pohcconst
230	ohdconst(pos_spec)=pohdconst
231	ohnconst(pos_spec)=pohnconst
232
233	do j=1,24 ! 24 hours, starting with 0-1 local time
234	area_hour(pos_spec,j)=parea_hour(j)
235	point_hour(pos_spec,j)=ppoint_hour(j)
236	end do
237	do j=1,7 ! 7 days of the week, starting with Monday
238	area_dow(pos_spec,j)=parea_dow(j)
239	point_dow(pos_spec,j)=ppoint_dow(j)
240	end do
241
242	endif
243
244	i=pos_spec
245
246	!NIK 16.02.2015
247	! Check scavenging parameters given in SPECIES file
248
249	if (lroot) then
250	! ZHG 2016.04.07 Start of changes
251	write(,) ' '
252	if (dquer(pos_spec) .gt.0) write(*,'(a,i3,a,a,a)') ' SPECIES: ', &
253	&id_spec,' ', species(pos_spec),' (AEROSOL) '
254	if (dquer(pos_spec) .le.0) write(*,'(a,i3,a,a,a)') ' SPECIES: ', &
255	&id_spec,' ', species(pos_spec),' (GAS) '
256
257	! Particles
258	!**********
259	if (dquer(pos_spec).gt.0) then
260	if (ccn_aero(pos_spec) .gt. 0) then
261	write(*,'(a,f5.2)') ' Particle CCN efficiency (CCNeff):', ccn_aero(pos_spec)
262	else
263	write(*,'(a)') ' Particle CCN efficiency (CCNeff): OFF'
264	endif
265	if (in_aero(pos_spec) .gt. 0) then
266	write(*,'(a,f5.2)') ' Particle IN efficiency (INeff) :', in_aero(pos_spec)
267	else
268	write(*,'(a)') ' Particle IN efficiency (INeff) : OFF'
269	endif
270	if (crain_aero(pos_spec) .gt. 0) then
271	write(*,'(a,f5.2)') ' Particle Rain efficiency (Crain) :', crain_aero(pos_spec)
272	else
273	write(*,'(a)') ' Particle Rain efficiency (Crain) : OFF'
274	endif
275	if (csnow_aero(pos_spec) .gt. 0) then
276	write(*,'(a,f5.2)') ' Particle Snow efficiency (Csnow) :', csnow_aero(pos_spec)
277	else
278	write(*,'(a)') ' Particle Snow efficiency (Csnow) : OFF'
279	end if
280	if (density(pos_spec) .gt. 0) then
281	write(*,'(a)') ' Dry deposition is turned : ON'
282	if (reldiff(pos_spec).gt.0) then
283	stop 'density>0 (SPECIES is a particle) implies reldiff <=0 '
284	endif
285	else
286	if (reldiff(pos_spec).le.0) then
287	stop 'density<=0 (SPECIES is a gas) implies reldiff >0 '
288	endif
289	write(*,'(a)') ' Dry deposition is (density<0) : OFF'
290	end if
291	if (crain_aero(pos_spec).gt.10.0 .or. csnow_aero(pos_spec).gt.10.0 .or. &
292	&ccn_aero(pos_spec).gt.1.0 .or. in_aero(pos_spec).gt.1.0) then
293	write(,) '*******************************************'
294	write(,) ' WARNING: Particle Scavenging parameter likely out of range '
295	write(,) ' Likely range for Crain 0.0-10'
296	write(,) ' Likely range for Csnow 0.0-10'
297	write(,) ' Physical range for CCNeff 0.0-1'
298	write(,) ' Physical range for INeff 0.0-1'
299	write(,) '*******************************************'
300	end if
301	else
302	! Gas
303	!****
304	if (weta_gas(pos_spec) .gt. 0 .and. wetb_gas(pos_spec).gt.0) then
305	write(,) ' Wet removal for gases is turned: ON'
306	write(,) ' Gas below-cloud scavenging parameter A ', &
307	&weta_gas(pos_spec)
308	write(*,'(a,f5.2)') ' Gas below-cloud scavenging parameter B ', &
309	&wetb_gas(pos_spec)
310	else
311	write(,) ' Wet removal for gases is turned: OFF '
312	end if
313	if (reldiff(i).gt.0.) then
314	write(,) ' Dry deposition for gases is turned: ON '
315	else
316	write(,) ' Dry deposition for gases is turned: OFF '
317	end if
318	if (weta_gas(pos_spec).gt.0.) then !if wet deposition is turned on
319	if (weta_gas(pos_spec).gt.1E-04 .or. weta_gas(pos_spec).lt.1E-09 .or. &
320	&wetb_gas(pos_spec).gt.0.8 .or. wetb_gas(pos_spec).lt.0.4) then
321	write(,) '*******************************************'
322	write(,) ' WARNING: Gas below-cloud scavengig is out of likely range'
323	write(,) ' Likely range for A is 1E-04 to 1E-08'
324	write(,) ' Likely range for B is 0.60 to 0.80 '
325	write(,) '*******************************************'
326	end if
327	endif
328
329	if (((weta_gas(pos_spec).gt.0).or.(wetb_gas(pos_spec).gt.0)).and.&
330	&(henry(pos_spec).le.0)) then
331	if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
332	endif
333	end if
334	end if
335
336	! if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception
337	if (dquer(i).gt.0 .and. dsigma(i).le.1.) then !dsigma(i)=1.0001 ! avoid floating exception
338	!write(,) '#### FLEXPART MODEL ERROR! ####'
339	write(,) '#### FLEXPART MODEL WARNING ####'
340	write(,) '#### in SPECIES_',aspecnumb, ' ####'
341	write(,) '#### from v10.4 dsigma has to be larger than 1 ####'
342	write(,) '#### to adapt older SPECIES files, ####'
343	write(,) '#### if dsigma was < 1 ####'
344	write(,) '#### use the reciprocal of the old dsigma ####'
345	if (.not.debug_mode) then
346	stop
347	else
348	write(,) 'debug mode: continue'
349	endif
350	endif
351
352	if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
353	write(,) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####'
354	write(,) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####'
355	write(,) '#### PARTICLE AND GAS. ####'
356	write(,) '#### SPECIES NUMBER',aspecnumb
357	stop
358	endif
359	20 continue
360
361
362	22 close(unitspecies)
363
364	! namelist output if requested
365	if (nmlout.and.lroot) then
366	open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000)
367	write(unitspecies,nml=species_params)
368	close(unitspecies)
369	endif
370
371	return
372
373	996 write(,) '#####################################################'
374	write(,) '#### FLEXPART MODEL ERROR! #### '
375	write(,) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS #### '
376	write(,) '#### CONSTANT IS SET ####'
377	write(,) '#### PLEASE MODIFY SPECIES DESCR. FILE! #### '
378	write(,) '#####################################################'
379	stop
380
381
382	997 write(,) '#####################################################'
383	write(,) '#### FLEXPART MODEL ERROR! #### '
384	write(,) '#### THE ASSSOCIATED SPECIES HAS TO BE DEFINED #### '
385	write(,) '#### BEFORE THE ONE WHICH POINTS AT IT #### '
386	write(,) '#### PLEASE CHANGE ORDER IN RELEASES OR ADD #### '
387	write(,) '#### THE ASSOCIATED SPECIES IN RELEASES #### '
388	write(,) '#####################################################'
389	stop
390
391
392	998 write(,) '#####################################################'
393	write(,) '#### FLEXPART MODEL ERROR! #### '
394	write(,) '#### THE SPECIES FILE FOR SPECIES ', id_spec
395	write(,) '#### CANNOT BE FOUND: CREATE FILE'
396	write(,) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb
397	write(,) '#####################################################'
398	stop
399
400	1000 write(,) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
401	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
402	write(*,'(a)') path(2)(1:length(2))
403	stop
404
405	end subroutine readspecies

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format