Context Navigation

timemanager.f90 @ c0884a8

univie

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

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

Download in other formats:

Original Format