Context Navigation

readreleases.f90 @ d7935de

univie

Last change on this file since d7935de was d7935de, checked in by pesei <petra seibert at univie ac at>, 6 years ago

modify most input read subroutines

changed some variable names (mostly for I-N reasons)
includes two names appearing also in timemanager, com_mod
corrected a few mistakes
simplified some parts of code
changed options/RELEASES which is in nml fmt correspondingly

Property mode set to 100644

File size: 26.1 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	! HSO, 12 August 2013: Added optional namelist input *
36	! PS, 2/2015: access= -> position= (ported from 9.2 review) *
37	! PS, 8/2018: check for uninitialised mass (ticket:204), *
38	! rename int/real variables, a bit of code layout improvement *
39	! *
40	!*****************************************************************************
41	! *
42	! Variables: *
43	! decay decay constant of species *
44	! dquer [um] mean particle diameters *
45	! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass*
46	! are between 0.1dquer and 10dquer *
47	! ireleasestart, ireleaseend [s] starting time and ending time of each *
48	! release *
49	! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint*
50	! is in hPa *
51	! npart number of particles to be released *
52	! nspec number of species to be released *
53	! density [kg/m3] density of the particles *
54	! rm [s/m] Mesophyll resistance *
55	! species name of species *
56	! xmass total mass of each species *
57	! xpoint1,ypoint1 geograf. coordinates of lower left corner of release *
58	! area *
59	! xpoint2,ypoint2 geograf. coordinates of upper right corner of release *
60	! area *
61	! weta_gas, wetb_gas parameters for below-cloud scavenging (gas) *
62	! crain_aero, csnow_aero parameters for below-cloud scavenging (aerosol) *
63	! ccn_aero, in_aero parameters for in-cloud scavenging (aerosol) *
64	! zpoint1,zpoint2 height range, over which release takes place *
65	! num_min_discrete if less, release cannot be randomized and happens at *
66	! time mid-point of release interval *
67	! lroot true if serial version, or if MPI and root process *
68	! *
69	!*****************************************************************************
70
71	use point_mod
72	use xmass_mod
73	use par_mod
74	use com_mod
75
76	implicit none
77
78	integer :: numpartmax,i,j,id1,it1,id2,it2,idum,iallocstat,irel,ispc,nsettle
79	integer,parameter :: num_min_discrete=100
80	real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
81	real(kind=dp) :: jul1,jul2,julm,juldate
82	character(len=50) :: line
83	logical :: old
84
85	! help variables for namelist reading
86	integer :: numpoints, nparts, ios
87	integer*2 :: ikindz
88	integer :: idate1, itime1, idate2, itime2
89	real :: rlon1,rlon2,rlat1,rlat2,z1,z2
90	character*40 :: comment
91	integer,parameter :: unitreleasesout=2
92	real,allocatable, dimension (:) :: rmass
93	integer,allocatable, dimension (:) :: numspec_rel,numspec_rel2
94
95	! declare namelists
96	namelist /releases_ctrl/ &
97	nspec, &
98	numspec_rel
99
100	namelist /release/ &
101	idate1, itime1, &
102	idate2, itime2, &
103	rlon1, rlon2, &
104	rlat1, rlat2, &
105	z1, z2, &
106	ikindz, &
107	rmass, &
108	nparts, &
109	comment
110	! initialiase mass
111
112	numpoint=0
113
114	! allocate with maxspec for first input loop
115	allocate(rmass(maxspec),stat=iallocstat)
116	if (iallocstat.ne.0) write(,)'ERROR: could not allocate rmass'
117	allocate(numspec_rel(maxspec),stat=iallocstat)
118	if (iallocstat.ne.0) write(,)'ERROR: could not allocate numspec_rel'
119
120	! presetting namelist releases_ctrl
121	nspec = -1 ! use negative value to determine failed namelist input
122	numspec_rel = 0
123
124	!SEC: read release to find how many releasepoints should be allocated
125	open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
126	form='formatted',err=999)
127
128	! prepare namelist output if requested
129	if (nmlout.and.lroot) then
130	open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist', &
131	position='append',status='replace',err=1000)
132	endif
133
134	! check whether namelist input is available
135	read(unitreleases,releases_ctrl,iostat=ios)
136
137	if (ios .ne. 0 .or. nspec .lt. 0 ) then
138
139	! not in namelist format, close file and reopen as simple text file
140
141	close(unitreleases)
142	open(unitreleases,file=path(1)(1:length(1))//'RELEASES', status='old', &
143	err=999)
144
145	ios=1 ! indicates old format
146
147	! Check the format of the RELEASES file
148	! (either in free format, or using a formatted mask)
149	! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
150	!**************************************************************************
151	! Note PS: something has changed, now Total is checked. Is this compatible??
152
153	call skplin(12,unitreleases)
154	read (unitreleases,901) line
155	901 format (a)
156	if (index(line,'Total') .eq. 0) then
157	old = .false.
158	else
159	old = .true.
160	endif
161	rewind(unitreleases)
162
163
164	! Skip first 11 lines (file header)
165	!**********************************
166
167	call skplin(11,unitreleases)
168
169	read(unitreleases,*,err=998) nspec
170	if (old) call skplin(2,unitreleases)
171	do i=1,nspec
172	read(unitreleases,*,err=998) numspec_rel(i)
173	if (old) call skplin(2,unitreleases)
174	end do
175
176	numpoint=0
177	100 numpoint=numpoint+1
178	read(unitreleases,*,end=25)
179	read(unitreleases,*,err=998,end=25) idum,idum
180	if (old) call skplin(2,unitreleases)
181	read(unitreleases,*,err=998) idum,idum
182	if (old) call skplin(2,unitreleases)
183	read(unitreleases,*,err=998) xdum
184	if (old) call skplin(2,unitreleases)
185	read(unitreleases,*,err=998) xdum
186	if (old) call skplin(2,unitreleases)
187	read(unitreleases,*,err=998) xdum
188	if (old) call skplin(2,unitreleases)
189	read(unitreleases,*,err=998) xdum
190	if (old) call skplin(2,unitreleases)
191	read(unitreleases,*,err=998) idum
192	if (old) call skplin(2,unitreleases)
193	read(unitreleases,*,err=998) xdum
194	if (old) call skplin(2,unitreleases)
195	read(unitreleases,*,err=998) xdum
196	if (old) call skplin(2,unitreleases)
197	read(unitreleases,*,err=998) idum
198	if (old) call skplin(2,unitreleases)
199	do i=1,nspec
200	read(unitreleases,*,err=998) xdum
201	if (old) call skplin(2,unitreleases)
202	end do
203	!save compoint only for the first 1000 release points
204	read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
205	if (old) call skplin(1,unitreleases)
206
207	goto 100
208
209	25 numpoint=numpoint-1
210
211	else ! read namelist input
212
213	ios=0
214	do while (ios.eq.0) ! loop over release points
215	idate1=-1
216	read(unitreleases,release,iostat=ios)
217	if (idate1.lt.0 .or. ios.ne.0) then
218	ios=1
219	else
220	numpoint=numpoint+1
221	endif
222	end do
223	ios=0
224
225	endif
226
227	rewind(unitreleases)
228
229	if (nspec.gt.maxspec) goto 994
230
231	! allocate arrays of matching size for number of species (namelist output)
232	deallocate(rmass)
233	allocate(rmass(nspec),stat=iallocstat)
234	if (iallocstat.ne.0) write(,)'ERROR: could not allocate rmass'
235	allocate(numspec_rel2(nspec),stat=iallocstat)
236	if (iallocstat.ne.0) write(,)'ERROR: could not allocate numspec_rel2'
237
238	numspec_rel2 = numspec_rel(1:nspec)
239
240	deallocate(numspec_rel)
241	! eso: BUG, crashes here for nspec=12 and maxspec=6,
242	! TODO: catch error and exit
243	allocate(numspec_rel(nspec),stat=iallocstat)
244	if (iallocstat.ne.0) write(,)'ERROR: could not allocate numspec_rel'
245
246	numspec_rel = numspec_rel2
247
248	deallocate(numspec_rel2)
249
250	! allocate memory for numpoint releaspoints
251	allocate ( &
252	ireleasestart(numpoint), ireleaseend(numpoint),xpoint1(numpoint), &
253	xpoint2(numpoint), ypoint1(numpoint), ypoint2(numpoint), &
254	zpoint1(numpoint), zpoint2(numpoint), kindz(numpoint), &
255	xmass(numpoint,maxspec), rho_rel(numpoint), npart(numpoint), &
256	xmasssave(numpoint), &
257	stat=iallocstat )
258	if (iallocstat.ne.0) write(,)'ERROR: could not allocate variables' // &
259	' dimensioned numpoint =', numpoint
260
261	if (lroot) write (,) 'Release points : ', numpoint
262
263	xmasssave(:)=0.
264
265	!now save the information
266	dep =.false.
267	drydep=.false.
268	wetdep=.false.
269	ohrea =.false.
270	drydepspec(:)=.false.
271	wetdepspec(:)=.false.
272
273	if (ios.ne.0) then
274
275	! Skip first 11 lines (file header)
276	!**********************************
277
278	call skplin(11,unitreleases)
279
280	! Assign species-specific parameters needed for physical processes
281	!*************************************************************************
282
283	read(unitreleases,*,err=998) nspec
284	if (nspec.gt.maxspec) goto 994
285	if (old) call skplin(2,unitreleases)
286
287	endif
288
289	! namelist output
290	if (nmlout.and.lroot) write(unitreleasesout,nml=releases_ctrl)
291
292	! species loop
293	do i=1,nspec
294
295	if (ios.ne.0) then
296	read(unitreleases,*,err=998) numspec_rel(i)
297	if (old) call skplin(2,unitreleases)
298	call readspecies(numspec_rel(i),i)
299	else
300	call readspecies(numspec_rel(i),i)
301	endif
302
303	! For backward runs, only 1 species is allowed
304	!*********************************************
305
306	!if ((ldirect.lt.0).and.(nspec.gt.1)) then
307	!write(,) '#####################################################'
308	!write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
309	!write(,) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
310	!write(,) '#####################################################'
311	! stop
312	!endif
313
314	! Molecular weight
315	!*****************
316
317	if ( (iout.eq.2 .or. iout.eq.3) .and. weightmolar(i).lt.0. ) then
318	write(,) 'For mixing ratio output, valid molar weight'
319	write(,) 'must be specified for all simulated species.'
320	write(,) 'Check table SPECIES or choose concentration'
321	write(,) 'output instead if molar weight is not known.'
322	stop
323	endif
324
325	! Radioactive decay
326	!******************
327
328	decay(i)=0.693147/decay(i) !conversion half life to decay constant
329	! Note PS: this should be replace by = alog(2.)/decay(i)
330	! keep it for the moment so that results are not affect in test phase of v10
331
332	! Dry deposition of gases
333	!************************
334
335	if (reldiff(i).gt.0.) rm(i)=1./( henry(i)/3000. + 100.*f0(i) )! mesophyll resistance
336
337	! Dry deposition of particles
338	!****************************
339
340	vsetaver(i) = 0.
341	cunningham(i) = 0.
342	dquer(i) = dquer(i)*1.e6 ! Conversion m to um
343	if (density(i).gt.0.) then ! Additional parameters
344	call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
345	do j=1,ni
346	fract(i,j) = fracth(j)
347	schmi(i,j) = schmih(j)
348	vset(i,j) = vsh(j)
349	cunningham(i) = cunningham(i) + cun*fract(i,j)
350	vsetaver(i) = vsetaver(i) - vset(i,j)*fract(i,j)
351	end do
352	if (lroot) write(,) 'Average settling velocity: ',i,vsetaver(i)
353	endif
354
355	! Dry deposition for constant deposition velocity
356	!************************************************
357
358	dryvel(i)=dryvel(i)*0.01 ! conversion to m/s
359
360	! Check if wet deposition or OH reaction shall be calculated
361	!***********************************************************
362
363	! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
364	if ( &
365	dquer(i).le.0. .and. (weta_gas(i) .gt.0. .or. wetb_gas(i) .gt.0.) .or. &
366	dquer(i).gt.0. .and. (crain_aero(i).gt.0. .or. csnow_aero(i).gt.0.) &
367	) then
368	wetdep = .true.
369	wetdepspec(i) = .true.
370	if (lroot) then
371	write (,) ' Below-cloud scavenging: ON'
372	! write (,) 'Below-cloud scavenging coefficients: ',weta(i),i
373	end if
374	else
375	if (lroot) write (,) ' Below-cloud scavenging: OFF'
376	endif
377
378	! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
379	if (dquer(i).gt.0. .and. (ccn_aero(i).gt.0. .or. in_aero(i).gt.0.)) then
380	wetdep = .true.
381	wetdepspec(i) = .true.
382	if (lroot) then
383	write (,) ' In-cloud scavenging: ON'
384	! write (,) 'In-cloud scavenging coefficients: ',&
385	! &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
386	end if
387	else
388	if (lroot) write (,) ' In-cloud scavenging: OFF'
389	endif
390
391	if (ohcconst(i).gt.0.) then
392	ohrea = .true.
393	if (lroot) write (,) ' OHreaction switched on: ',ohcconst(i),i
394	endif
395
396	if (reldiff(i).gt.0. .or. density(i).gt.0. .or. dryvel(i).gt.0.) then
397	drydep = .true.
398	drydepspec(i) = .true.
399	endif
400
401	end do ! end loop over species
402
403	if (WETDEP.or.DRYDEP) DEP=.true.
404
405	! Read specifications for each release point
406	!*******************************************
407	numpoints=numpoint
408	numpoint=0
409	numpartmax=0
410	releaserate=0.
411	101 numpoint=numpoint+1
412
413	if (ios.lt.1) then ! reading namelist format
414
415	if (numpoint.gt.numpoints) goto 250
416	ikindz = 1
417	rmass = -1. ! initialise as negative to catch if not given in input
418	nparts = 0
419	comment = ' '
420	read(unitreleases,release,iostat=ios)
421	id1=idate1
422	it1=itime1
423	id2=idate2
424	it2=itime2
425	xpoint1(numpoint)=rlon1
426	xpoint2(numpoint)=rlon2
427	ypoint1(numpoint)=rlat1
428	ypoint2(numpoint)=rlat2
429	zpoint1(numpoint)=z1
430	zpoint2(numpoint)=z2
431	kindz(numpoint)=ikindz
432	do i=1,nspec
433	if (rmass(i) .lt. 0.) goto 995 ! write error and stop
434	xmass(numpoint,i)=rmass(i)
435	end do
436	npart(numpoint)=nparts
437	compoint(min(1001,numpoint))=comment
438
439	! namelist output
440	if (nmlout.and.lroot) then
441	write(unitreleasesout,nml=release)
442	endif
443
444	else ! simple text format
445
446	read(unitreleases,*,end=250)
447	read(unitreleases,*,err=998,end=250) id1,it1
448	if (old) call skplin(2,unitreleases)
449	read(unitreleases,*,err=998) id2,it2
450	if (old) call skplin(2,unitreleases)
451	read(unitreleases,*,err=998) xpoint1(numpoint)
452	if (old) call skplin(2,unitreleases)
453	read(unitreleases,*,err=998) ypoint1(numpoint)
454	if (old) call skplin(2,unitreleases)
455	read(unitreleases,*,err=998) xpoint2(numpoint)
456	if (old) call skplin(2,unitreleases)
457	read(unitreleases,*,err=998) ypoint2(numpoint)
458	if (old) call skplin(2,unitreleases)
459	read(unitreleases,*,err=998) kindz(numpoint)
460	if (old) call skplin(2,unitreleases)
461	read(unitreleases,*,err=998) zpoint1(numpoint)
462	if (old) call skplin(2,unitreleases)
463	read(unitreleases,*,err=998) zpoint2(numpoint)
464	if (old) call skplin(2,unitreleases)
465	read(unitreleases,*,err=998) npart(numpoint)
466	if (old) call skplin(2,unitreleases)
467	do i=1,nspec
468	read(unitreleases,*,err=998) xmass(numpoint,i)
469	if (old) call skplin(2,unitreleases)
470	rmass(i)=xmass(numpoint,i)
471	end do
472	!save compoint only for the first 1000 release points
473	if (numpoint.le.1000) then
474	read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
475	comment=compoint(numpoint)(1:40)
476	else
477	read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
478	comment=compoint(1001)(1:40)
479	endif
480	if (old) call skplin(1,unitreleases)
481
482	! namelist output
483	if (nmlout.and.lroot) then
484	idate1=id1
485	itime1=it1
486	idate2=id2
487	itime2=it2
488	rlon1=xpoint1(numpoint)
489	rlon2=xpoint2(numpoint)
490	rlat1=ypoint1(numpoint)
491	rlat2=ypoint2(numpoint)
492	z1=zpoint1(numpoint)
493	z2=zpoint2(numpoint)
494	ikindz=kindz(numpoint)
495	nparts=npart(numpoint)
496	write(unitreleasesout,nml=release)
497	endif
498
499	if (numpoint.le.1000) then
500	if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
501	xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0. .and. &
502	compoint(numpoint)(1:8).eq.' ') goto 250
503	else
504	if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
505	xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0.) goto 250
506	endif
507
508	endif ! if namelist format
509
510	! If a release point contains no particles, stop and issue error message
511	!***********************************************************************
512
513	if (npart(numpoint).eq.0) then
514	write(,) 'FLEXPART MODEL ERROR'
515	write(,) 'RELEASES file is corrupt.'
516	write(,) 'At least for one release point, there are zero'
517	write(,) 'particles released. Make changes to RELEASES.'
518	stop
519	endif
520
521	! If FLEXPART is run for backward deposition force zpoint
522	!*********************************************************************
523	if (WETBKDEP) then
524	zpoint1(numpoint)=0.
525	zpoint2(numpoint)=20000.
526	kindz(numpoint)=1
527	endif
528	if (DRYBKDEP) then
529	zpoint1(numpoint)=0.
530	zpoint2(numpoint)=2.*href
531	kindz(numpoint)=1
532	endif
533
534
535	! Check whether x coordinates of release point are within model domain
536	!*********************************************************************
537
538	if (xpoint1(numpoint) .lt. xlon0) &
539	xpoint1(numpoint) = xpoint1(numpoint) + 360.
540	if (xpoint1(numpoint) .gt. xlon0+nxmin1*dx) &
541	xpoint1(numpoint) = xpoint1(numpoint) - 360.
542	if (xpoint2(numpoint) .lt. xlon0) &
543	xpoint2(numpoint) = xpoint2(numpoint) + 360.
544	if (xpoint2(numpoint) .gt. xlon0+nxmin1*dx) &
545	xpoint2(numpoint) = xpoint2(numpoint) - 360.
546
547	! Determine relative beginning and ending times of particle release
548	!******************************************************************
549
550	jul1=juldate(id1,it1)
551	jul2=juldate(id2,it2)
552	julm=(jul1+jul2)/2.
553	if (jul1.gt.jul2) then
554	write(,) 'FLEXPART MODEL ERROR'
555	write(,) 'Release stops before it begins.'
556	write(,) 'Make changes to file RELEASES.'
557	stop
558	endif
559	if (mdomainfill.eq.0) then ! no domain filling
560	if (ldirect.eq.1) then
561	if (jul1.lt.bdate .or. jul2.gt.edate) then
562	write(,) 'FLEXPART MODEL ERROR'
563	write(,) 'Release starts before simulation begins or ends'
564	write(,) 'after simulation stops.'
565	write(,) 'Make files COMMAND and RELEASES consistent.'
566	stop
567	endif
568	if (npart(numpoint).gt.num_min_discrete) then
569	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
570	ireleaseend(numpoint)=int((jul2-bdate)*86400.)
571	else
572	ireleasestart(numpoint)=int((julm-bdate)*86400.)
573	ireleaseend(numpoint)=int((julm-bdate)*86400.)
574	endif
575	else if (ldirect.eq.-1) then
576	if (jul1.lt.edate .or. jul2.gt.bdate) then
577	write(,) 'FLEXPART MODEL ERROR'
578	write(,) 'Release starts before simulation begins or ends'
579	write(,) 'after simulation stops.'
580	write(,) 'Make files COMMAND and RELEASES consistent.'
581	stop
582	endif
583	if (npart(numpoint).gt.num_min_discrete) then
584	ireleasestart(numpoint)=int((jul1-bdate)*86400.)
585	ireleaseend(numpoint) = int((jul2-bdate)*86400.)
586	else
587	ireleasestart(numpoint)=int((julm-bdate)*86400.)
588	ireleaseend(numpoint) = int((julm-bdate)*86400.)
589	endif
590	endif
591	endif
592
593
594	! Determine the release rate (particles per second) and total number
595	! of particles released during the simulation
596	!*******************************************************************
597
598	if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
599	releaserate = releaserate + real( npart(numpoint) ) / &
600	real( ireleaseend(numpoint) - ireleasestart(numpoint) )
601	else
602	releaserate=99999999
603	endif
604	numpartmax=numpartmax+npart(numpoint)
605	goto 101
606
607	250 close(unitreleases)
608
609	if (nmlout.and.lroot) then
610	close(unitreleasesout)
611	endif
612
613	if (lroot) write (,) 'Particles allocated (maxpart) : ',maxpart
614	if (lroot) write (,) 'Particles released (numpartmax): ',numpartmax
615	numpoint=numpoint-1
616
617	if (ioutputforeachrelease.eq.1) then
618	maxpointspec_act=numpoint
619	else
620	maxpointspec_act=1
621	endif
622
623	! Disable settling if more than 1 species at any release point
624	! or if MQUASILAG and more than one species
625	!*************************************************************
626
627	if (mquasilag.ne.0) then
628	if (nspec.gt.1) lsettling=.false.
629	else
630	do irel=1,numpoint
631	nsettle=0
632	do ispc=1,nspec
633	if (xmass(irel,ispc).gt.tiny(1.)) nsettle=nsettle+1
634	end do
635	if (nsettle.gt.1) lsettling=.false.
636	end do
637	end if
638
639	if (lroot) then
640	if (.not.lsettling) &
641	write(,) 'WARNING: more than 1 species per release point,' &
642	// ' settling disabled'
643	end if
644
645	! Check, whether the total number of particles may exceed totally allowed
646	! number of particles at some time during the simulation
647	!************************************************************************
648
649	if (releaserate .gt. 0.99*real(maxpart)/real(lage(nageclass))) then
650	if (numpartmax.gt.maxpart.and.lroot) then
651	write(,) '#####################################################'
652	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
653	write(,) '#### ####'
654	write(,) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
655	write(,) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
656	write(,) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
657	write(,) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
658	write(,) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
659	write(,) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
660	write(,) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
661	write(,) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
662	write(,) '#####################################################'
663	write(,) 'Maximum release rate may be: ',releaserate, &
664	' particles per second'
665	write(,) 'Maximum allowed release rate is: ', &
666	real(maxpart)/real(lage(nageclass)),' particles per second'
667	write(,) 'Total number of particles released during the simulation is:'&
668	,numpartmax
669	write(,) 'Maximum allowed number of particles is: ',maxpart
670	endif
671	endif
672
673
674	if (lroot) then
675	write(*,FMT='(A,ES14.7)') ' Total mass released:', &
676	sum(xmass(1:numpoint,1:nspec))
677	end if
678
679	return
680
681	995 write(,) '#####################################################'
682	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
683	write(,) '#### ####'
684	write(,) '#### ERROR - mass for species',i,' not given ####'
685	write(,) '#####################################################'
686	stop 'readreleases - undefined mass'
687
688	994 write(,) '######################################################'
689	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
690	write(,) '#### ####'
691	write(,) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS ####'
692	write(,) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
693	write(,) '######################################################'
694	stop
695
696	998 write(,) '#####################################################'
697	write(,) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
698	write(,) '#### ####'
699	write(,) '#### FATAL ERROR - FILE "RELEASES" IS ####'
700	write(,) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
701	write(,) '#### MISTAKES OR GET A NEW "RELEASES"- ####'
702	write(,) '#### FILE ... ####'
703	write(,) '#####################################################'
704	stop
705
706	999 write(,) '#####################################################'
707	write(,) ' FLEXPART MODEL SUBROUTINE READRELEASES: '
708	write(,)
709	write(,) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
710	write(,) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
711	write(,) 'PERMITTED FOR ANY ACCESS'
712	write(,) '#####################################################'
713	stop
714
715	1000 write(,) ' #### FLEXPART MODEL ERROR! FILE "RELEASES" #### '
716	write(,) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
717	write(*,'(a)') path(2)(1:length(2))
718	stop
719
720	end subroutine readreleases

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

Download in other formats:

Original Format