Context Navigation

timemanager.f90 @ c7d1052

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

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

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

Download in other formats:

Original Format