Context Navigation

readreleases.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 11 years ago
move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others
Property svn:executable set to ``*
File size: 19.6 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 readreleases
23
24	!*****************************************************************************
25	! *
26	! This routine reads the release point specifications for the current *
27	! model run. Several release points can be used at the same time. *
28	! *
29	! Author: A. Stohl *
30	! *
31	! 18 May 1996 *
32	! *
33	! Update: 29 January 2001 *
34	! Release altitude can be either in magl or masl *
35	! *
36	!*****************************************************************************
37	! *
38	! Variables: *
39	! decay decay constant of species *
40	! dquer [um] mean particle diameters *
41	! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
42	! are between 0.1dquer and 10dquer *
43	! ireleasestart, ireleaseend [s] starting time and ending time of each *
44	! release *
45	! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
46	! is in hPa *
47	! npart number of particles to be released *
48	! nspec number of species to be released *
49	! density [kg/m3] density of the particles *
50	! rm [s/m] Mesophyll resistance *
51	! species name of species *
52	! xmass total mass of each species *
53	! xpoint1,ypoint1 geograf. coordinates of lower left corner of release *
54	! area *
55	! xpoint2,ypoint2 geograf. coordinates of upper right corner of release *
56	! area *
57	! weta, wetb parameters to determine the wet scavenging coefficient *
58	! zpoint1,zpoint2 height range, over which release takes place *
59	! num_min_discrete if less, release cannot be randomized and happens at *
60	! time mid-point of release interval *
61	! *
62	!*****************************************************************************
63
64	use point_mod
65	use xmass_mod
66	use par_mod
67	use com_mod
68
69	implicit none
70
71	integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
72	integer,parameter :: num_min_discrete=100
73	real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
74	real(kind=dp) :: jul1,jul2,julm,juldate
75	character(len=50) :: line
76	logical :: old
77
78	!sec, read release to find how many releasepoints should be allocated
79
80	open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
81	err=999)
82
83	! Check the format of the RELEASES file (either in free format,
84	! or using a formatted mask)
85	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
86	!**************************************************************************
87
88	call skplin(12,unitreleases)
89	read (unitreleases,901) line
90	901 format (a)
91	if (index(line,'Total') .eq. 0) then
92	old = .false.
93	else
94	old = .true.
95	endif
96	rewind(unitreleases)
97
98
99	! Skip first 11 lines (file header)
100	!**********************************
101
102	call skplin(11,unitreleases)
103
104
105	read(unitreleases,*,err=998) nspec
106	if (old) call skplin(2,unitreleases)
107	do i=1,nspec
108	read(unitreleases,*,err=998) specnum_rel
109	if (old) call skplin(2,unitreleases)
110	end do
111
112	numpoint=0
113	100 numpoint=numpoint+1
114	read(unitreleases,*,end=25)
115	read(unitreleases,*,err=998,end=25) idum,idum
116	if (old) call skplin(2,unitreleases)
117	read(unitreleases,*,err=998) idum,idum
118	if (old) call skplin(2,unitreleases)
119	read(unitreleases,*,err=998) xdum
120	if (old) call skplin(2,unitreleases)
121	read(unitreleases,*,err=998) xdum
122	if (old) call skplin(2,unitreleases)
123	read(unitreleases,*,err=998) xdum
124	if (old) call skplin(2,unitreleases)
125	read(unitreleases,*,err=998) xdum
126	if (old) call skplin(2,unitreleases)
127	read(unitreleases,*,err=998) idum
128	if (old) call skplin(2,unitreleases)
129	read(unitreleases,*,err=998) xdum
130	if (old) call skplin(2,unitreleases)
131	read(unitreleases,*,err=998) xdum
132	if (old) call skplin(2,unitreleases)
133	read(unitreleases,*,err=998) idum
134	if (old) call skplin(2,unitreleases)
135	do i=1,nspec
136	read(unitreleases,*,err=998) xdum
137	if (old) call skplin(2,unitreleases)
138	end do
139	!save compoint only for the first 1000 release points
140	read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
141	if (old) call skplin(1,unitreleases)
142
143	goto 100
144
145	25 numpoint=numpoint-1
146
147	!allocate memory for numpoint releaspoint
148	allocate(ireleasestart(numpoint) &
149	,stat=stat)
150	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
151	allocate(ireleaseend(numpoint) &
152	,stat=stat)
153	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
154	allocate(xpoint1(numpoint) &
155	,stat=stat)
156	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
157	allocate(xpoint2(numpoint) &
158	,stat=stat)
159	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
160	allocate(ypoint1(numpoint) &
161	,stat=stat)
162	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
163	allocate(ypoint2(numpoint) &
164	,stat=stat)
165	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
166	allocate(zpoint1(numpoint) &
167	,stat=stat)
168	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
169	allocate(zpoint2(numpoint) &
170	,stat=stat)
171	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
172	allocate(kindz(numpoint) &
173	,stat=stat)
174	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
175	allocate(xmass(numpoint,maxspec) &
176	,stat=stat)
177	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
178	allocate(rho_rel(numpoint) &
179	,stat=stat)
180	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
181	allocate(npart(numpoint) &
182	,stat=stat)
183	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
184	allocate(xmasssave(numpoint) &
185	,stat=stat)
186	if (stat.ne.0) write(,)'ERROR: could not allocate RELEASPOINT'
187
188	write (,) ' Releasepoints allocated: ', numpoint
189
190	do i=1,numpoint
191	xmasssave(i)=0.
192	end do
193
194	!now save the information
195	DEP=.false.
196	DRYDEP=.false.
197	WETDEP=.false.
198	OHREA=.false.
199	do i=1,maxspec
200	DRYDEPSPEC(i)=.false.
201	end do
202
203	rewind(unitreleases)
204
205
206	! Skip first 11 lines (file header)
207	!**********************************
208
209	call skplin(11,unitreleases)
210
211
212	! Assign species-specific parameters needed for physical processes
213	!*************************************************************************
214
215	read(unitreleases,*,err=998) nspec
216	if (nspec.gt.maxspec) goto 994
217	if (old) call skplin(2,unitreleases)
218	do i=1,nspec
219	read(unitreleases,*,err=998) specnum_rel
220	if (old) call skplin(2,unitreleases)
221	call readspecies(specnum_rel,i)
222
223	! For backward runs, only 1 species is allowed
224	!*********************************************
225
226	!if ((ldirect.lt.0).and.(nspec.gt.1)) then
227	!write(,) '#####################################################'
228	!write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
229	!write(,) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
230	!write(,) '#####################################################'
231	! stop
232	!endif
233
234	! Molecular weight
235	!*****************
236
237	if (((iout.eq.2).or.(iout.eq.3)).and. &
238	(weightmolar(i).lt.0.)) then
239	write(,) 'For mixing ratio output, valid molar weight'
240	write(,) 'must be specified for all simulated species.'
241	write(,) 'Check table SPECIES or choose concentration'
242	write(,) 'output instead if molar weight is not known.'
243	stop
244	endif
245
246
247	! Radioactive decay
248	!******************
249
250	decay(i)=0.693147/decay(i) !conversion half life to decay constant
251
252
253	! Dry deposition of gases
254	!************************
255
256	if (reldiff(i).gt.0.) &
257	rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance
258
259	! Dry deposition of particles
260	!****************************
261
262	vsetaver(i)=0.
263	cunningham(i)=0.
264	dquer(i)=dquer(i)*1000000. ! Conversion m to um
265	if (density(i).gt.0.) then ! Additional parameters
266	call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
267	do j=1,ni
268	fract(i,j)=fracth(j)
269	schmi(i,j)=schmih(j)
270	vset(i,j)=vsh(j)
271	cunningham(i)=cunningham(i)+cun*fract(i,j)
272	vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
273	end do
274	write(,) 'Average setting velocity: ',i,vsetaver(i)
275	endif
276
277	! Dry deposition for constant deposition velocity
278	!************************************************
279
280	dryvel(i)=dryvel(i)*0.01 ! conversion to m/s
281
282	! Check if wet deposition or OH reaction shall be calculated
283	!***********************************************************
284	if (weta(i).gt.0.) then
285	WETDEP=.true.
286	write (,) 'Wetdeposition switched on: ',weta(i),i
287	endif
288	if (ohreact(i).gt.0) then
289	OHREA=.true.
290	write (,) 'OHreaction switched on: ',ohreact(i),i
291	endif
292
293
294	if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
295	(dryvel(i).gt.0.)) then
296	DRYDEP=.true.
297	DRYDEPSPEC(i)=.true.
298	endif
299
300	end do
301
302	if (WETDEP.or.DRYDEP) DEP=.true.
303
304	! Read specifications for each release point
305	!*******************************************
306
307	numpoint=0
308	numpartmax=0
309	releaserate=0.
310	1000 numpoint=numpoint+1
311	read(unitreleases,*,end=250)
312	read(unitreleases,*,err=998,end=250) id1,it1
313	if (old) call skplin(2,unitreleases)
314	read(unitreleases,*,err=998) id2,it2
315	if (old) call skplin(2,unitreleases)
316	read(unitreleases,*,err=998) xpoint1(numpoint)
317	if (old) call skplin(2,unitreleases)
318	read(unitreleases,*,err=998) ypoint1(numpoint)
319	if (old) call skplin(2,unitreleases)
320	read(unitreleases,*,err=998) xpoint2(numpoint)
321	if (old) call skplin(2,unitreleases)
322	read(unitreleases,*,err=998) ypoint2(numpoint)
323	if (old) call skplin(2,unitreleases)
324	read(unitreleases,*,err=998) kindz(numpoint)
325	if (old) call skplin(2,unitreleases)
326	read(unitreleases,*,err=998) zpoint1(numpoint)
327	if (old) call skplin(2,unitreleases)
328	read(unitreleases,*,err=998) zpoint2(numpoint)
329	if (old) call skplin(2,unitreleases)
330	read(unitreleases,*,err=998) npart(numpoint)
331	if (old) call skplin(2,unitreleases)
332	do i=1,nspec
333	read(unitreleases,*,err=998) xmass(numpoint,i)
334	if (old) call skplin(2,unitreleases)
335	end do
336	!save compoint only for the first 1000 release points
337	if (numpoint.le.1000) then
338	read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
339	else
340	read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
341	endif
342	if (old) call skplin(1,unitreleases)
343	if (numpoint.le.1000) then
344	if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
345	(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
346	(compoint(numpoint)(1:8).eq.' ')) goto 250
347	else
348	if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
349	(xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
350	endif
351
352
353	! If a release point contains no particles, stop and issue error message
354	!***********************************************************************
355
356	if (npart(numpoint).eq.0) then
357	write(,) 'FLEXPART MODEL ERROR'
358	write(,) 'RELEASES file is corrupt.'
359	write(,) 'At least for one release point, there are zero'
360	write(,) 'particles released. Make changes to RELEASES.'
361	stop
362	endif
363
364	! Check whether x coordinates of release point are within model domain
365	!*********************************************************************
366
367	if (xpoint1(numpoint).lt.xlon0) &
368	xpoint1(numpoint)=xpoint1(numpoint)+360.
369	if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
370	xpoint1(numpoint)=xpoint1(numpoint)-360.
371	if (xpoint2(numpoint).lt.xlon0) &
372	xpoint2(numpoint)=xpoint2(numpoint)+360.
373	if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
374	xpoint2(numpoint)=xpoint2(numpoint)-360.
375
376	! Determine relative beginning and ending times of particle release
377	!******************************************************************
378
379	jul1=juldate(id1,it1)
380	jul2=juldate(id2,it2)
381	julm=(jul1+jul2)/2.
382	if (jul1.gt.jul2) then
383	write(,) 'FLEXPART MODEL ERROR'
384	write(,) 'Release stops before it begins.'
385	write(,) 'Make changes to file RELEASES.'
386	stop
387	endif
388	if (mdomainfill.eq.0) then ! no domain filling
389	if (ldirect.eq.1) then
390	if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
391	write(,) 'FLEXPART MODEL ERROR'
392	write(,) 'Release starts before simulation begins or ends'
393	write(,) 'after simulation stops.'
394	write(,) 'Make files COMMAND and RELEASES consistent.'
395	stop
396	endif
397	if (npart(numpoint).gt.num_min_discrete) then
398	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
399	ireleaseend(numpoint)=int((jul2-bdate)*86400.)
400	else
401	ireleasestart(numpoint)=int((julm-bdate)*86400.)
402	ireleaseend(numpoint)=int((julm-bdate)*86400.)
403	endif
404	else if (ldirect.eq.-1) then
405	if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
406	write(,) 'FLEXPART MODEL ERROR'
407	write(,) 'Release starts before simulation begins or ends'
408	write(,) 'after simulation stops.'
409	write(,) 'Make files COMMAND and RELEASES consistent.'
410	stop
411	endif
412	if (npart(numpoint).gt.num_min_discrete) then
413	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
414	ireleaseend(numpoint)=int((jul2-bdate)*86400.)
415	else
416	ireleasestart(numpoint)=int((julm-bdate)*86400.)
417	ireleaseend(numpoint)=int((julm-bdate)*86400.)
418	endif
419	endif
420	endif
421
422
423	! Check, whether the total number of particles may exceed totally allowed
424	! number of particles at some time during the simulation
425	!************************************************************************
426
427	! Determine the release rate (particles per second) and total number
428	! of particles released during the simulation
429	!*******************************************************************
430
431	if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
432	releaserate=releaserate+real(npart(numpoint))/ &
433	real(ireleaseend(numpoint)-ireleasestart(numpoint))
434	else
435	releaserate=99999999
436	endif
437	numpartmax=numpartmax+npart(numpoint)
438	goto 1000
439
440
441	250 close(unitreleases)
442
443	write (,) ' Particles allocated for this run: ',maxpart, ', released in simulation: ', numpartmax
444	numpoint=numpoint-1
445
446	if (ioutputforeachrelease.eq.1) then
447	maxpointspec_act=numpoint
448	else
449	maxpointspec_act=1
450	endif
451
452	if (releaserate.gt. &
453	0.99*real(maxpart)/real(lage(nageclass))) then
454	if (numpartmax.gt.maxpart) then
455	write(,) '#####################################################'
456	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
457	write(,) '#### ####'
458	write(,) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
459	write(,) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
460	write(,) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
461	write(,) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
462	write(,) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
463	write(,) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
464	write(,) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
465	write(,) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
466	write(,) '#####################################################'
467	write(,) 'Maximum release rate may be: ',releaserate, &
468	' particles per second'
469	write(,) 'Maximum allowed release rate is: ', &
470	real(maxpart)/real(lage(nageclass)),' particles per second'
471	write(,) &
472	'Total number of particles released during the simulation is: ', &
473	numpartmax
474	write(,) 'Maximum allowed number of particles is: ',maxpart
475	endif
476	endif
477
478	return
479
480	994 write(,) '#####################################################'
481	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
482	write(,) '#### ####'
483	write(,) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
484	write(,) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
485	write(,) '#####################################################'
486	stop
487
488	998 write(,) '#####################################################'
489	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
490	write(,) '#### ####'
491	write(,) '#### FATAL ERROR - FILE "RELEASES" IS ####'
492	write(,) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
493	write(,) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
494	write(,) '#### FILE ... ####'
495	write(,) '#####################################################'
496	stop
497
498
499	999 write(,) '#####################################################'
500	write(,) ' FLEXPART MODEL SUBROUTINE READRELEASES: '
501	write(,)
502	write(,) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
503	write(,) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
504	write(,) 'PERMITTED FOR ANY ACCESS'
505	write(,) '#####################################################'
506	stop
507
508	end subroutine readreleases

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

Download in other formats:

Original Format