Context Navigation

timemanager.f90 @ f3054ea

GFS_025dev

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

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

Download in other formats:

Original Format