Context Navigation

readreleases.f90 @ 3481cc1

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

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 24.6 KB

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

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

Download in other formats:

Original Format