Context Navigation

timemanager.f90 @ 4138764

dev

Last change on this file since 4138764 was 4138764, checked in by Sabine <sabine.eckhardt@…>, 3 years ago
Merge remote-tracking branch 'refs/remotes/origin/dev' into dev
Property mode set to `100644`
File size: 29.4 KB

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

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

Download in other formats:

Original Format