Context Navigation

timemanager.f90 @ 11d86e7

dev

Last change on this file since 11d86e7 was 11d86e7, checked in by Sabine <sabine.eckhardt@…>, 4 years ago
bugfix in the fluxoutput
Property mode set to `100644`
File size: 29.2 KB

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

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

Download in other formats:

Original Format