Context Navigation

timemanager.f90 @ d6a0977

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

Last change on this file since d6a0977 was d6a0977, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago
Updates to Henrik's wet depo scheme
Property mode set to `100644`
File size: 28.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 timemanager
23
24	!*****************************************************************************
25	! *
26	! Handles the computation of trajectories, i.e. determines which *
27	! trajectories have to be computed at what time. *
28	! Manages dry+wet deposition routines, radioactive decay and the computation *
29	! of concentrations. *
30	! *
31	! Author: A. Stohl *
32	! *
33	! 20 May 1996 *
34	! *
35	!*****************************************************************************
36	! Changes, Bernd C. Krueger, Feb. 2001: *
37	! Call of convmix when new windfield is read *
38	!------------------------------------ *
39	! Changes Petra Seibert, Sept 2002 *
40	! fix wet scavenging problem *
41	! Code may not be correct for decay of deposition! *
42	! Changes Petra Seibert, Nov 2002 *
43	! call convection BEFORE new fields are read in BWD mode *
44	! Changes Caroline Forster, Feb 2005 *
45	! new interface between flexpart and convection scheme *
46	! Emanuel's latest subroutine convect43c.f is used *
47	! Changes Stefan Henne, Harald Sodemann, 2013-2014 *
48	! added netcdf output code *
49	! Changes Espen Sollum 2014 *
50	! For compatibility with MPI version, *
51	! variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod *
52	!*****************************************************************************
53	! *
54	! Variables: *
55	! DEP .true. if either wet or dry deposition is switched on *
56	! decay(maxspec) [1/s] decay constant for radioactive decay *
57	! DRYDEP .true. if dry deposition is switched on *
58	! ideltas [s] modelling period *
59	! itime [s] actual temporal position of calculation *
60	! ldeltat [s] time since computation of radioact. decay of depositions*
61	! loutaver [s] averaging period for concentration calculations *
62	! loutend [s] end of averaging for concentration calculations *
63	! loutnext [s] next time at which output fields shall be centered *
64	! loutsample [s] sampling interval for averaging of concentrations *
65	! loutstart [s] start of averaging for concentration calculations *
66	! loutstep [s] time interval for which concentrations shall be *
67	! calculated *
68	! npoint(maxpart) index, which starting point the trajectory has *
69	! starting positions of trajectories *
70	! nstop serves as indicator for fate of particles *
71	! in the particle loop *
72	! nstop1 serves as indicator for wind fields (see getfields) *
73	! outnum number of samples for each concentration calculation *
74	! outnum number of samples for each concentration calculation *
75	! prob probability of absorption at ground due to dry *
76	! deposition *
77	! WETDEP .true. if wet deposition is switched on *
78	! weight weight for each concentration sample (1/2 or 1) *
79	! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to *
80	! turbulence *
81	! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to inter- *
82	! polation *
83	! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = *
84	! spatial positions of trajectories *
85	! *
86	! Constants: *
87	! maxpart maximum number of trajectories *
88	! *
89	!*****************************************************************************
90
91	use unc_mod
92	use point_mod
93	use xmass_mod
94	use flux_mod
95	use outg_mod
96	use oh_mod
97	use par_mod
98	use com_mod
99	use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
100	&concoutput_surf_netcdf,concoutput_surf_nest_netcdf
101
102	implicit none
103
104	integer :: j,ks,kp,l,n,itime,nstop,nstop1
105	! integer :: ksp
106	integer :: loutnext,loutstart,loutend
107	integer :: ix,jy,ldeltat,itage,nage
108	integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
109	real :: outnum,weight,prob(maxspec),decfact
110	! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
111	! real :: us(maxpart),vs(maxpart),ws(maxpart)
112	! integer(kind=2) :: cbt(maxpart)
113	real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
114	real :: drygridtotalunc,xold,yold,zold,xmassfract
115	!double precision xm(maxspec,maxpointspec_act),
116	! + xm_depw(maxspec,maxpointspec_act),
117	! + xm_depd(maxspec,maxpointspec_act)
118
119
120	!open(88,file='TEST.dat')
121
122	! First output for time 0
123	!************************
124
125	loutnext=loutstep/2
126	outnum=0.
127	loutstart=loutnext-loutaver/2
128	loutend=loutnext+loutaver/2
129
130	! open(127,file=path(2)(1:length(2))//'depostat.dat'
131	! + ,form='unformatted')
132	!write (,) 'writing deposition statistics depostat.dat!'
133
134	!**********************************************************************
135	! Loop over the whole modelling period in time steps of mintime seconds
136	!**********************************************************************
137
138	!ZHG 2015
139	!CGZ-lifetime: set lifetime to 0
140	checklifetime(:,:)=0
141	species_lifetime(:,:)=0
142	print*, 'Initialized lifetime'
143	!CGZ-lifetime: set lifetime to 0
144
145
146
147	if (verbosity.gt.0) then
148	write (,) 'timemanager> starting simulation'
149	if (verbosity.gt.1) then
150	CALL SYSTEM_CLOCK(count_clock)
151	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
152	endif
153	endif
154
155	do itime=0,ideltas,lsynctime
156
157	! Computation of wet deposition, OH reaction and mass transfer
158	! between two species every lsynctime seconds
159	! maybe wet depo frequency can be relaxed later but better be on safe side
160	! wetdepo must be called BEFORE new fields are read in but should not
161	! be called in the very beginning before any fields are loaded, or
162	! before particles are in the system
163	! Code may not be correct for decay of deposition
164	! changed by Petra Seibert 9/02
165	!********************************************************************
166
167	if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
168	if (verbosity.gt.0) then
169	write (,) 'timemanager> call wetdepo'
170	endif
171	call wetdepo(itime,lsynctime,loutnext)
172	endif
173
174	if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
175	call ohreaction(itime,lsynctime,loutnext)
176
177	if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
178	stop 'associated species not yet implemented!'
179	! call transferspec(itime,lsynctime,loutnext)
180	endif
181
182	! compute convection for backward runs
183	!*************************************
184
185	if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
186	if (verbosity.gt.0) then
187	write (,) 'timemanager> call convmix -- backward'
188	endif
189	call convmix(itime)
190	if (verbosity.gt.1) then
191	!CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
192	CALL SYSTEM_CLOCK(count_clock)
193	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
194	endif
195	endif
196
197	! Get necessary wind fields if not available
198	!*******************************************
199	if (verbosity.gt.0) then
200	write (,) 'timemanager> call getfields'
201	endif
202	call getfields(itime,nstop1)
203	if (verbosity.gt.1) then
204	CALL SYSTEM_CLOCK(count_clock)
205	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
206	endif
207	if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
208
209	! Get hourly OH fields if not available
210	!****************************************************
211	if (OHREA) then
212	if (verbosity.gt.0) then
213	write (,) 'timemanager> call gethourlyOH'
214	endif
215	call gethourlyOH(itime)
216	if (verbosity.gt.1) then
217	CALL SYSTEM_CLOCK(count_clock)
218	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
219	endif
220	endif
221
222	! Release particles
223	!******************
224
225	if (verbosity.gt.0) then
226	write (,) 'timemanager> Release particles'
227	endif
228
229	if (mdomainfill.ge.1) then
230	if (itime.eq.0) then
231	if (verbosity.gt.0) then
232	write (,) 'timemanager> call init_domainfill'
233	endif
234	call init_domainfill
235	else
236	if (verbosity.gt.0) then
237	write (,) 'timemanager> call boundcond_domainfill'
238	endif
239	call boundcond_domainfill(itime,loutend)
240	endif
241	else
242	if (verbosity.gt.0) then
243	print*,'call releaseparticles'
244	endif
245	call releaseparticles(itime)
246	if (verbosity.gt.1) then
247	CALL SYSTEM_CLOCK(count_clock)
248	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
249	endif
250	endif
251
252
253	! Compute convective mixing for forward runs
254	! for backward runs it is done before next windfield is read in
255	!**************************************************************
256
257	if ((ldirect.eq.1).and.(lconvection.eq.1)) then
258	if (verbosity.gt.0) then
259	write (,) 'timemanager> call convmix -- forward'
260	endif
261	call convmix(itime)
262	endif
263
264	! If middle of averaging period of output fields is reached, accumulated
265	! deposited mass radioactively decays
266	!***********************************************************************
267
268	if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
269	do ks=1,nspec
270	do kp=1,maxpointspec_act
271	if (decay(ks).gt.0.) then
272	do nage=1,nageclass
273	do l=1,nclassunc
274	! Mother output grid
275	do jy=0,numygrid-1
276	do ix=0,numxgrid-1
277	wetgridunc(ix,jy,ks,kp,l,nage)= &
278	wetgridunc(ix,jy,ks,kp,l,nage)* &
279	exp(-1.outstepdecay(ks))
280	drygridunc(ix,jy,ks,kp,l,nage)= &
281	drygridunc(ix,jy,ks,kp,l,nage)* &
282	exp(-1.outstepdecay(ks))
283	end do
284	end do
285	! Nested output grid
286	if (nested_output.eq.1) then
287	do jy=0,numygridn-1
288	do ix=0,numxgridn-1
289	wetgriduncn(ix,jy,ks,kp,l,nage)= &
290	wetgriduncn(ix,jy,ks,kp,l,nage)* &
291	exp(-1.outstepdecay(ks))
292	drygriduncn(ix,jy,ks,kp,l,nage)= &
293	drygriduncn(ix,jy,ks,kp,l,nage)* &
294	exp(-1.outstepdecay(ks))
295	end do
296	end do
297	endif
298	end do
299	end do
300	endif
301	end do
302	end do
303	endif
304
305	!!! CHANGE: These lines may be switched on to check the conservation
306	!!! of mass within FLEXPART
307	! if (itime.eq.loutnext) then
308	! do 247 ksp=1, nspec
309	! do 247 kp=1, maxpointspec_act
310	!47 xm(ksp,kp)=0.
311
312	! do 249 ksp=1, nspec
313	! do 249 j=1,numpart
314	! if (ioutputforeachrelease.eq.1) then
315	! kp=npoint(j)
316	! else
317	! kp=1
318	! endif
319	! if (itra1(j).eq.itime) then
320	! xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
321	! write(,) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
322	! endif
323	!49 continue
324	! do 248 ksp=1,nspec
325	! do 248 kp=1,maxpointspec_act
326	! xm_depw(ksp,kp)=0.
327	! xm_depd(ksp,kp)=0.
328	! do 248 nage=1,nageclass
329	! do 248 ix=0,numxgrid-1
330	! do 248 jy=0,numygrid-1
331	! do 248 l=1,nclassunc
332	! xm_depw(ksp,kp)=xm_depw(ksp,kp)
333	! + +wetgridunc(ix,jy,ksp,kp,l,nage)
334	!48 xm_depd(ksp,kp)=xm_depd(ksp,kp)
335	! + +drygridunc(ix,jy,ksp,kp,l,nage)
336	! do 246 ksp=1,nspec
337	!46 write(88,'(2i10,3e12.3)')
338	! + itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
339	! + (xm_depw(ksp,kp),kp=1,maxpointspec_act),
340	! + (xm_depd(ksp,kp),kp=1,maxpointspec_act)
341	! endif
342	!!! CHANGE
343
344
345
346	! Check whether concentrations are to be calculated
347	!**************************************************
348
349	if ((ldirectitime.ge.ldirectloutstart).and. &
350	(ldirectitime.le.ldirectloutend)) then ! add to grid
351	if (mod(itime-loutstart,loutsample).eq.0) then
352
353	! If we are exactly at the start or end of the concentration averaging interval,
354	! give only half the weight to this sample
355	!*****************************************************************************
356
357	if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
358	weight=0.5
359	else
360	weight=1.0
361	endif
362	outnum=outnum+weight
363	call conccalc(itime,weight)
364	endif
365
366
367	if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
368	call partoutput_short(itime) ! dump particle positions in extremely compressed format
369
370
371	! Output and reinitialization of grid
372	! If necessary, first sample of new grid is also taken
373	!*****************************************************
374
375	if ((itime.eq.loutend).and.(outnum.gt.0.)) then
376	if ((iout.le.3.).or.(iout.eq.5)) then
377	if (surf_only.ne.1) then
378	if (lnetcdfout.eq.1) then
379	call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
380	else
381	call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
382	endif
383	else
384	if (verbosity.eq.1) then
385	print*,'call concoutput_surf '
386	call system_clock(count_clock)
387	write(,) 'system clock',count_clock - count_clock0
388	endif
389	if (lnetcdfout.eq.1) then
390	call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
391	else
392	call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
393	if (verbosity.eq.1) then
394	print*,'called concoutput_surf '
395	call system_clock(count_clock)
396	write(,) 'system clock',count_clock - count_clock0
397	endif
398	endif
399	endif
400
401	if (nested_output .eq. 1) then
402	if (lnetcdfout.eq.0) then
403	if (surf_only.ne.1) then
404	call concoutput_nest(itime,outnum)
405	else
406	call concoutput_surf_nest(itime,outnum)
407	endif
408	else
409	if (surf_only.ne.1) then
410	call concoutput_nest_netcdf(itime,outnum)
411	else
412	call concoutput_surf_nest_netcdf(itime,outnum)
413	endif
414	endif
415	endif
416	outnum=0.
417	endif
418	if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
419	if (iflux.eq.1) call fluxoutput(itime)
420	write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
421
422	!CGZ-lifetime: output species lifetime
423	!ZHG
424	write(,) 'Overview species lifetime in days', &
425	real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
426	write(,) 'all info:',species_lifetime
427	!ZHG
428	!CGZ-lifetime: output species lifetime
429
430	!write(*,46) float(itime)/3600,itime,numpart
431	45 format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES: Uncertainty: ',3f7.3)
432	46 format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
433	if (ipout.ge.1) call partoutput(itime) ! dump particle positions
434	loutnext=loutnext+loutstep
435	loutstart=loutnext-loutaver/2
436	loutend=loutnext+loutaver/2
437	if (itime.eq.loutstart) then
438	weight=0.5
439	outnum=outnum+weight
440	call conccalc(itime,weight)
441	endif
442
443
444	! Check, whether particles are to be split:
445	! If so, create new particles and attribute all information from the old
446	! particles also to the new ones; old and new particles both get half the
447	! mass of the old ones
448	!************************************************************************
449
450	if (ldirectitime.ge.ldirectitsplit) then
451	n=numpart
452	do j=1,numpart
453	if (ldirectitime.ge.ldirectitrasplit(j)) then
454	if (n.lt.maxpart) then
455	n=n+1
456	itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
457	itrasplit(n)=itrasplit(j)
458	itramem(n)=itramem(j)
459	itra1(n)=itra1(j)
460	idt(n)=idt(j)
461	npoint(n)=npoint(j)
462	nclass(n)=nclass(j)
463	xtra1(n)=xtra1(j)
464	ytra1(n)=ytra1(j)
465	ztra1(n)=ztra1(j)
466	uap(n)=uap(j)
467	ucp(n)=ucp(j)
468	uzp(n)=uzp(j)
469	us(n)=us(j)
470	vs(n)=vs(j)
471	ws(n)=ws(j)
472	cbt(n)=cbt(j)
473	do ks=1,nspec
474	xmass1(j,ks)=xmass1(j,ks)/2.
475	xmass1(n,ks)=xmass1(j,ks)
476	end do
477	endif
478	endif
479	end do
480	numpart=n
481	endif
482	endif
483	endif
484
485
486	if (itime.eq.ideltas) exit ! almost finished
487
488	! Compute interval since radioactive decay of deposited mass was computed
489	!************************************************************************
490
491	if (itime.lt.loutnext) then
492	ldeltat=itime-(loutnext-loutstep)
493	else ! first half of next interval
494	ldeltat=itime-loutnext
495	endif
496
497
498	! Loop over all particles
499	!************************
500	! Various variables for testing reason of CBL scheme, by mc
501	well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
502	well_mixed_norm=0. !erase normalization to test well mixed condition: modified by mc
503	avg_ol=0.
504	avg_wst=0.
505	avg_h=0.
506	avg_air_dens=0. !erase vector to obtain air density at particle positions: modified by mc
507	!-----------------------------------------------------------------------------
508	do j=1,numpart
509
510
511	! If integration step is due, do it
512	!**********************************
513
514	if (itra1(j).eq.itime) then
515
516	if (ioutputforeachrelease.eq.1) then
517	kp=npoint(j)
518	else
519	kp=1
520	endif
521	! Determine age class of the particle
522	itage=abs(itra1(j)-itramem(j))
523	do nage=1,nageclass
524	if (itage.lt.lage(nage)) exit
525	end do
526
527	! Initialize newly released particle
528	!***********************************
529
530	if ((itramem(j).eq.itime).or.(itime.eq.0)) &
531	call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
532	us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
533
534	! Memorize particle positions
535	!****************************
536
537	xold=xtra1(j)
538	yold=ytra1(j)
539	zold=ztra1(j)
540
541	! Integrate Lagevin equation for lsynctime seconds
542	!*************************************************
543
544	call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
545	us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
546	cbt(j))
547
548	! Calculate the gross fluxes across layer interfaces
549	!***************************************************
550
551	if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
552
553
554	! Determine, when next time step is due
555	! If trajectory is terminated, mark it
556	!**************************************
557
558	if (nstop.gt.1) then
559	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
560	itra1(j)=-999999999
561	else
562	itra1(j)=itime+lsynctime
563
564
565	! Dry deposition and radioactive decay for each species
566	! Also check maximum (of all species) of initial mass remaining on the particle;
567	! if it is below a threshold value, terminate particle
568	!*****************************************************************************
569
570	xmassfract=0.
571	do ks=1,nspec
572	if (decay(ks).gt.0.) then ! radioactive decay
573	decfact=exp(-real(abs(lsynctime))*decay(ks))
574	else
575	decfact=1.
576	endif
577
578	if (DRYDEPSPEC(ks)) then ! dry deposition
579	drydeposit(ks)=xmass1(j,ks)prob(ks)decfact
580	xmass1(j,ks)=xmass1(j,ks)(1.-prob(ks))decfact
581	if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
582	drydeposit(ks)=drydeposit(ks)* &
583	exp(real(abs(ldeltat))*decay(ks))
584	endif
585	else ! no dry deposition
586	xmass1(j,ks)=xmass1(j,ks)*decfact
587	endif
588
589
590	if (mdomainfill.eq.0) then
591	if (xmass(npoint(j),ks).gt.0.) &
592	xmassfract=max(xmassfract,real(npart(npoint(j)))* &
593	xmass1(j,ks)/xmass(npoint(j),ks))
594	!ZHG 2015
595	!CGZ-lifetime: Check mass fraction left/save lifetime
596	if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
597	!Mass below 1% of initial >register lifetime
598	checklifetime(j,ks)=abs(itra1(j)-itramem(j))
599	species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
600	species_lifetime(ks,2)= species_lifetime(ks,2)+1
601	endif
602	!CGZ-lifetime: Check mass fraction left/save lifetime
603	!ZHG 2015
604	else
605	xmassfract=1.
606	endif
607	end do
608
609	if (xmassfract.lt.0.0001) then ! terminate all particles carrying less mass
610	itra1(j)=-999999999
611	if (verbosity.gt.0) then
612	print*,'terminated particle ',j,' for small mass'
613	endif
614	endif
615
616	! Sabine Eckhardt, June 2008
617	! don't create depofield for backward runs
618	if (DRYDEP.AND.(ldirect.eq.1)) then
619	call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
620	real(ytra1(j)),nage,kp)
621	if (nested_output.eq.1) call drydepokernel_nest( &
622	nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
623	nage,kp)
624	endif
625
626	! Terminate trajectories that are older than maximum allowed age
627	!***************************************************************
628
629	if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
630	if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
631	itra1(j)=-999999999
632	if (verbosity.gt.0) then
633	print*,'terminated particle ',j,' for age'
634	endif
635	endif
636	endif
637
638	endif
639
640	end do !loop over particles
641
642	! Counter of "unstable" particle velocity during a time scale of
643	! maximumtl=20 minutes (defined in com_mod)
644	!***************************************************************
645
646	total_nan_intl=0
647	i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
648	sum_nan_count(i_nan)=nan_count
649	if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
650	do ii_nan=1, (maxtl/lsynctime)
651	total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
652	end do
653	! Output to keep track of the numerical instabilities in CBL simulation and if
654	! they are compromising the final result (or not)
655	if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl
656
657	!!------------------------------------------------------------------------------
658	! these lines below to test the well-mixed condition, modified by mc, not to be
659	! included in final release:
660	! if (itime.eq.0) then
661	! open(551,file='/home/mc/test_cbl/out/WELLMIXEDTEST_CBL_lonlat_9_33_100_3hours_3htp_cd.DAT')
662	! open(552,file='/home/mc/test_cbl/out/avg_ol_h_wst_lonlat_9_33_100_3hours_3htp_cd.DAT')
663	! end if
664	! write(552,'(5F16.7)')itime*1./3600.,avg_wst/well_mixed_norm,avg_ol/well_mixed_norm,avg_h/well_mixed_norm
665	! do j=1,25
666	! !write(551,))itime1.,h_well/50.j,well_mixed_vector(j)/well_mixed_norm50.
667	! avg_air_dens(j)=avg_air_dens(j)/well_mixed_vector(j)
668	! write(551,'(5F16.7)')itime1.,h_well/25.j,well_mixed_vector(j)/well_mixed_norm25.,avg_air_dens(j),0.04j
669	! end do
670	!!------------------------------------------------------------------------------
671
672	end do
673
674
675	! Complete the calculation of initial conditions for particles not yet terminated
676	!*****************************************************************************
677
678	do j=1,numpart
679	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
680	end do
681
682	if (ipout.eq.2) call partoutput(itime) ! dump particle positions
683
684	if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
685
686	!close(104)
687
688	! De-allocate memory and end
689	!***************************
690
691	if (iflux.eq.1) then
692	deallocate(flux)
693	endif
694	if (OHREA) then
695	deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
696	endif
697	if (ldirect.gt.0) then
698	deallocate(drygridunc,wetgridunc)
699	endif
700	deallocate(gridunc)
701	deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
702	deallocate(ireleasestart,ireleaseend,npart,kindz)
703	deallocate(xmasssave)
704	if (nested_output.eq.1) then
705	deallocate(orooutn, arean, volumen)
706	if (ldirect.gt.0) then
707	deallocate(griduncn,drygriduncn,wetgriduncn)
708	endif
709	endif
710	deallocate(outheight,outheighthalf)
711	deallocate(oroout, area, volume)
712
713	end subroutine timemanager
714

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

Download in other formats:

Original Format