Context Navigation

readreleases.f90

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

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

Download in other formats:

Original Format