Context Navigation

timemanager.f90 @ 278f4ed

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

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

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

Download in other formats:

Original Format