Context Navigation

readspecies.f90 @ 8a65cb0

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

Last change on this file since 8a65cb0 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago
Added code, makefile for dev branch
Property mode set to `100644`
File size: 14.7 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
22	subroutine 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, wetb Parameters for determining below-cloud scavenging *
45	! weta_in Parameter for determining in-cloud scavenging *
46	! wetb_in Parameter for determining in-cloud scavenging *
47	! ohcconst OH reaction rate constant C *
48	! ohdconst OH reaction rate constant D *
49	! id_spec SPECIES number as referenced in RELEASE file *
50	! id_pos position where SPECIES data shall be stored *
51	! *
52	! Constants: *
53	! *
54	!*****************************************************************************
55
56	use par_mod
57	use com_mod
58
59	implicit none
60
61	integer :: i, pos_spec,j
62	integer :: idow,ihour,id_spec
63	character(len=3) :: aspecnumb
64	logical :: spec_found
65
66	character(len=16) :: pspecies
67	real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
68	real :: pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
69	real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
70	integer :: readerror
71
72	! declare namelist
73	namelist /species_params/ &
74	pspecies, pdecay, pweta, pwetb, &
75	pweta_in, pwetb_in, pwetc_in, pwetd_in, &
76	preldiff, phenry, pf0, pdensity, pdquer, &
77	pdsigma, pdryvel, pweightmolar, pohcconst, pohdconst, pspec_ass, pkao
78
79	pspecies=" "
80	pdecay=-999.9
81	pweta=-9.9E-09
82	pwetb=0.0
83	pweta_in=-9.9E-09
84	pwetb_in=-9.9E-09
85	! pwetc_in=-9.9E-09
86	! pwetd_in=-9.9E-09
87	preldiff=-9.9
88	phenry=0.0
89	pf0=0.0
90	pdensity=-9.9E09
91	pdquer=0.0
92	pdsigma=0.0
93	pdryvel=-9.99
94	pohcconst=-9.99
95	pohdconst=-9.9E-09
96	pspec_ass=-9
97	pkao=-99.99
98	pweightmolar=-789.0 ! read failure indicator value
99
100	! Open the SPECIES file and read species names and properties
101	!************************************************************
102	specnum(pos_spec)=id_spec
103	write(aspecnumb,'(i3.3)') specnum(pos_spec)
104	open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
105	!write(,) 'reading SPECIES',specnum(pos_spec)
106
107	ASSSPEC=.FALSE.
108
109	! try namelist input
110	read(unitspecies,species_params,iostat=readerror)
111	close(unitspecies)
112
113	if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
114
115	readerror=1
116
117	open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
118
119	do i=1,6
120	read(unitspecies,*)
121	end do
122
123	read(unitspecies,'(a10)',end=22) species(pos_spec)
124	! write(,) species(pos_spec)
125	read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
126	! write(,) decay(pos_spec)
127	read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
128	! write(,) weta(pos_spec)
129	read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
130	! write(,) wetb(pos_spec)
131
132	!*** NIK 31.01.2013: including in-cloud scavening parameters
133	read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
134	! write(,) weta_in(pos_spec)
135	read(unitspecies,'(f18.2)',end=22) wetb_in(pos_spec)
136	! write(,) wetb_in(pos_spec)
137	! read(unitspecies,'(f18.2)',end=22) wetc_in(pos_spec)
138	! write(,) wetc_in(pos_spec)
139	! read(unitspecies,'(f18.2)',end=22) wetd_in(pos_spec)
140	! write(,) wetd_in(pos_spec)
141
142	read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
143	! write(,) reldiff(pos_spec)
144	read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
145	! write(,) henry(pos_spec)
146	read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
147	! write(,) f0(pos_spec)
148	read(unitspecies,'(e18.1)',end=22) density(pos_spec)
149	! write(,) density(pos_spec)
150	read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
151	! write(,) dquer(pos_spec)
152	read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
153	! write(,) dsigma(pos_spec)
154	read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
155	! write(,) dryvel(pos_spec)
156	read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
157	! write(,) weightmolar(pos_spec)
158	read(unitspecies,'(e18.2)',end=22) ohcconst(pos_spec)
159	! write(,) ohcconst(pos_spec)
160	read(unitspecies,'(f8.2)',end=22) ohdconst(pos_spec)
161	! write(,) ohdconst(pos_spec)
162	read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
163	! write(,) spec_ass(pos_spec)
164	read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
165	! write(,) kao(pos_spec)
166
167	pspecies=species(pos_spec)
168	pdecay=decay(pos_spec)
169	pweta=weta(pos_spec)
170	pwetb=wetb(pos_spec)
171	pweta_in=weta_in(pos_spec)
172	pwetb_in=wetb_in(pos_spec)
173	! pwetc_in=wetc_in(pos_spec)
174	! pwetd_in=wetd_in(pos_spec)
175	preldiff=reldiff(pos_spec)
176	phenry=henry(pos_spec)
177	pf0=f0(pos_spec)
178	pdensity=density(pos_spec)
179	pdquer=dquer(pos_spec)
180	pdsigma=dsigma(pos_spec)
181	pdryvel=dryvel(pos_spec)
182	pweightmolar=weightmolar(pos_spec)
183	pohcconst=ohcconst(pos_spec)
184	pohdconst=ohdconst(pos_spec)
185	pspec_ass=spec_ass(pos_spec)
186	pkao=kao(pos_spec)
187
188	else
189
190	species(pos_spec)=pspecies
191	decay(pos_spec)=pdecay
192	weta(pos_spec)=pweta
193	wetb(pos_spec)=pwetb
194	weta_in(pos_spec)=pweta_in
195	wetb_in(pos_spec)=pwetb_in
196	! wetc_in(pos_spec)=pwetc_in
197	! wetd_in(pos_spec)=pwetd_in
198	reldiff(pos_spec)=preldiff
199	henry(pos_spec)=phenry
200	f0(pos_spec)=pf0
201	density(pos_spec)=pdensity
202	dquer(pos_spec)=pdquer
203	dsigma(pos_spec)=pdsigma
204	dryvel(pos_spec)=pdryvel
205	weightmolar(pos_spec)=pweightmolar
206	ohcconst(pos_spec)=pohcconst
207	ohdconst(pos_spec)=pohdconst
208	spec_ass(pos_spec)=pspec_ass
209	kao(pos_spec)=pkao
210
211	endif
212
213	i=pos_spec
214
215	!NIK 16.02.2015
216	! Check scavenging parameters given in SPECIES file
217	if (dquer(pos_spec).gt.0) then !is particle
218	write(*,'(a,f5.2)') ' Particle below-cloud scavenging parameter A &
219	&(Rain collection efficiency) ', weta(pos_spec)
220	write(*,'(a,f5.2)') ' Particle below-cloud scavenging parameter B &
221	&(Snow collection efficiency) ', wetb(pos_spec)
222	if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
223	write(,) '*******************************************'
224	write(,) ' WARNING: Particle below-cloud scavenging parameter A or B &
225	&is out of likely range'
226	write(,) ' Likely range is 0.0-1.0'
227	write(,) '*******************************************'
228	endif
229	write(*,'(a,f5.2)') ' Particle in-cloud scavenging parameter Ai &
230	&(CCN efficiency) ', weta_in(pos_spec)
231	write(*,'(a,f5.2)') ' Particle in-cloud scavenging parameter Bi &
232	&(IN efficiency) ', wetb_in(pos_spec)
233	if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7 )then
234	write(,) '*******************************************'
235	write(,) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
236	write(,) ' The default value is 0.9 for CCN '
237	write(,) '*******************************************'
238	endif
239	if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
240	write(,) '*******************************************'
241	write(,) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
242	write(,) ' The default value is 0.1 for IN '
243	write(,) '*******************************************'
244	endif
245
246	else !is gas
247	write(*,'(a,f5.2)') ' Gas below-cloud scavenging parameter A ', weta(pos_spec)
248	write(*,'(a,f5.2)') ' Gas below-cloud scavenging parameter B ', wetb(pos_spec)
249	write(,) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
250	if (wetb(pos_spec).gt.0.) then !if wet deposition is turned on
251	if (wetb(pos_spec).gt.0.8 .or. wetb(pos_spec).lt.0.6) then
252	write(,) '*******************************************'
253	write(,) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
254	write(,) ' Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
255	write(,) '*******************************************'
256	endif
257	end if
258	if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
259	if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
260	write(,) '*******************************************'
261	write(,) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
262	write(,) ' Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
263	write(,) '*******************************************'
264	endif
265	end if
266	endif
267
268
269	if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
270	if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
271	endif
272
273	if (spec_ass(pos_spec).gt.0) then
274	spec_found=.FALSE.
275	do j=1,pos_spec-1
276	if (spec_ass(pos_spec).eq.specnum(j)) then
277	spec_ass(pos_spec)=j
278	spec_found=.TRUE.
279	ASSSPEC=.TRUE.
280	endif
281	end do
282	if (spec_found.eqv..FALSE.) then
283	goto 997
284	endif
285	endif
286
287	if (dsigma(i).eq.1.) dsigma(i)=1.0001 ! avoid floating exception
288	if (dsigma(i).eq.0.) dsigma(i)=1.0001 ! avoid floating exception
289
290	if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
291	write(,) '#### FLEXPART MODEL ERROR! FILE "SPECIES" ####'
292	write(,) '#### IS CORRUPT. SPECIES CANNOT BE BOTH ####'
293	write(,) '#### PARTICLE AND GAS. ####'
294	write(,) '#### SPECIES NUMBER',aspecnumb
295	stop
296	endif
297	20 continue
298
299
300	! Read in daily and day-of-week variation of emissions, if available
301	!*******************************************************************
302	! HSO: This is not yet implemented as namelist parameters
303
304	do j=1,24 ! initialize everything to no variation
305	area_hour(i,j)=1.
306	point_hour(i,j)=1.
307	end do
308	do j=1,7
309	area_dow(i,j)=1.
310	point_dow(i,j)=1.
311	end do
312
313	if (readerror.ne.0) then ! text format input
314
315	read(unitspecies,*,end=22)
316	do j=1,24 ! 24 hours, starting with 0-1 local time
317	read(unitspecies,*) ihour,area_hour(i,j),point_hour(i,j)
318	end do
319	read(unitspecies,*)
320	do j=1,7 ! 7 days of the week, starting with Monday
321	read(unitspecies,*) idow,area_dow(i,j),point_dow(i,j)
322	end do
323
324	endif
325
326	22 close(unitspecies)
327
328	! namelist output if requested
329	if (nmlout.and.lroot) then
330	open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='replace',err=1000)
331	write(unitspecies,nml=species_params)
332	close(unitspecies)
333	endif
334
335	return
336
337	996 write(,) '#####################################################'
338	write(,) '#### FLEXPART MODEL ERROR! #### '
339	write(,) '#### WET DEPOSITION SWITCHED ON, BUT NO HENRYS #### '
340	write(,) '#### CONSTANT IS SET ####'
341	write(,) '#### PLEASE MODIFY SPECIES DESCR. FILE! #### '
342	write(,) '#####################################################'
343	stop
344
345
346	997 write(,) '#####################################################'
347	write(,) '#### FLEXPART MODEL ERROR! #### '
348	write(,) '#### THE ASSSOCIATED SPECIES HAS TO BE DEFINED #### '
349	write(,) '#### BEFORE THE ONE WHICH POINTS AT IT #### '
350	write(,) '#### PLEASE CHANGE ORDER IN RELEASES OR ADD #### '
351	write(,) '#### THE ASSOCIATED SPECIES IN RELEASES #### '
352	write(,) '#####################################################'
353	stop
354
355
356	998 write(,) '#####################################################'
357	write(,) '#### FLEXPART MODEL ERROR! #### '
358	write(,) '#### THE SPECIES FILE FOR SPECIES ', id_spec
359	write(,) '#### CANNOT BE FOUND: CREATE FILE'
360	write(,) '#### ',path(1)(1:length(1)),'SPECIES/SPECIES_',aspecnumb
361	write(,) '#####################################################'
362	stop
363
364	1000 write(,) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
365	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
366	write(*,'(a)') path(2)(1:length(2))
367	stop
368
369	end subroutine readspecies

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

Download in other formats:

Original Format