Context Navigation

timemanager.f90 @ 37

Last change on this file since 37 was 37, checked in by pesei, 9 years ago
Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
File size: 24.4 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	!*****************************************************************************
50	! *
51	! Variables: *
52	! dep .true. if either wet or dry deposition is switched on *
53	! decay(maxspec) [1/s] decay constant for radioactive decay *
54	! drydep .true. if dry deposition is switched on *
55	! ideltas [s] modelling period *
56	! itime [s] actual temporal position of calculation *
57	! ldeltat [s] time since computation of radioact. decay of depositions*
58	! loutaver [s] averaging period for concentration calculations *
59	! loutend [s] end of averaging for concentration calculations *
60	! loutnext [s] next time at which output fields shall be centered *
61	! loutsample [s] sampling interval for averaging of concentrations *
62	! loutstart [s] start of averaging for concentration calculations *
63	! loutstep [s] time interval for which concentrations shall be *
64	! calculated *
65	! npoint(maxpart) index, which starting point the trajectory has *
66	! starting positions of trajectories *
67	! nstop serves as indicator for fate of particles *
68	! in the particle loop *
69	! nstop1 serves as indicator for wind fields (see getfields) *
70	! outnum number of samples for each concentration calculation *
71	! outnum number of samples for each concentration calculation *
72	! prob probability of absorption at ground due to dry *
73	! deposition *
74	! wetdep .true. if wet deposition is switched on *
75	! weight weight for each concentration sample (1/2 or 1) *
76	! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to *
77	! turbulence *
78	! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to inter- *
79	! polation *
80	! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = *
81	! spatial positions of trajectories *
82	! *
83	! Constants: *
84	! maxpart maximum number of trajectories *
85	! *
86	!*****************************************************************************
87
88	use unc_mod
89	use point_mod
90	use xmass_mod
91	use flux_mod
92	use outg_mod
93	use oh_mod
94	use par_mod
95	use com_mod
96	#ifdef NETCDF_OUTPUT
97	use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
98	#endif
99
100	implicit none
101
102	integer :: j,ks,kp,l,n,itime,nstop,nstop1
103	! integer :: ksp
104	integer :: loutnext,loutstart,loutend
105	integer :: ix,jy,ldeltat,itage,nage
106	real :: outnum,weight,prob(maxspec)
107	real :: uap(maxpart),ucp(maxpart),uzp(maxpart),decfact
108	real :: us(maxpart),vs(maxpart),ws(maxpart)
109	integer(kind=2) :: cbt(maxpart)
110	real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
111	real :: drygridtotalunc,xold,yold,zold,xmassfract
112	!double precision xm(maxspec,maxpointspec_act),
113	! + xm_depw(maxspec,maxpointspec_act),
114	! + xm_depd(maxspec,maxpointspec_act)
115
116
117	!open(88,file='TEST.dat')
118
119	! First output for time 0
120	!************************
121
122	loutnext=loutstep/2
123	outnum=0.
124	loutstart=loutnext-loutaver/2
125	loutend=loutnext+loutaver/2
126
127	! open(127,file=path(2)(1:length(2))//'depostat.dat'
128	! + ,form='unformatted')
129	!write (,) 'writing deposition statistics depostat.dat!'
130
131	!**********************************************************************
132	! Loop over the whole modelling period in time steps of mintime seconds
133	!**********************************************************************
134
135
136	itime=0
137	!write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
138	write(*,46) float(itime)/3600,itime,numpart
139	call verboutput(verbosity,1,'timemanager> starting simulation')
140	if (verbosity.gt.1) then
141	CALL SYSTEM_CLOCK(count_clock)
142	WRITE(,) 'timemanager> SYSTEM CLOCK', (count_clock-count_clock0) / real(count_rate)
143	endif
144
145	do itime=0,ideltas,lsynctime
146
147	! Computation of wet deposition, OH reaction and mass transfer
148	! between two species every lsynctime seconds
149	! maybe wet depo frequency can be relaxed later but better be on safe side
150	! wetdepo must be called BEFORE new fields are read in but should not
151	! be called in the very beginning before any fields are loaded, or
152	! before particles are in the system
153	! Code may not be correct for decay of deposition
154	! changed by Petra Seibert 9/02
155	!********************************************************************
156
157	if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
158	if (verbosity.gt.0) then
159	write (,) 'timemanager> call wetdepo'
160	endif
161	call wetdepo(itime,lsynctime,loutnext)
162	endif
163
164	if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
165	call ohreaction(itime,lsynctime,loutnext)
166
167	if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
168	stop 'associated species not yet implemented!'
169	! call transferspec(itime,lsynctime,loutnext)
170	endif
171
172	! compute convection for backward runs
173	!*************************************
174
175	if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
176	if (verbosity.gt.0) then
177	write (,) 'timemanager> call convmix -- backward'
178	endif
179	call convmix(itime)
180	if (verbosity.gt.1) then
181	!CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
182	CALL SYSTEM_CLOCK(count_clock)
183	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
184	endif
185	endif
186
187	! Get necessary wind fields if not available
188	!*******************************************
189	if (verbosity.gt.0) then
190	write (,) 'timemanager> call getfields'
191	endif
192	call getfields(itime,nstop1)
193	if (verbosity.gt.1) then
194	CALL SYSTEM_CLOCK(count_clock)
195	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
196	endif
197	if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
198
199	! Release particles
200	!******************
201
202	if (verbosity.gt.0) then
203	write (,) 'timemanager> Release particles'
204	endif
205
206	if (mdomainfill.ge.1) then
207	if (itime.eq.0) then
208	if (verbosity.gt.0) then
209	write (,) 'timemanager> call init_domainfill'
210	endif
211	call init_domainfill
212	else
213	if (verbosity.gt.0) then
214	write (,) 'timemanager> call boundcond_domainfill'
215	endif
216	call boundcond_domainfill(itime,loutend)
217	endif
218	else
219	if (verbosity.gt.0) then
220	print*,'call releaseparticles'
221	endif
222	call releaseparticles(itime)
223	if (verbosity.gt.1) then
224	CALL SYSTEM_CLOCK(count_clock)
225	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
226	endif
227	endif
228
229
230	! Compute convective mixing for forward runs
231	! for backward runs it is done before next windfield is read in
232	!**************************************************************
233
234	if ((ldirect.eq.1).and.(lconvection.eq.1)) then
235	if (verbosity.gt.0) then
236	write (,) 'timemanager> call convmix -- forward'
237	endif
238	call convmix(itime)
239	endif
240
241	! If middle of averaging period of output fields is reached, accumulated
242	! deposited mass radioactively decays
243	!***********************************************************************
244
245	if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
246	do ks=1,nspec
247	do kp=1,maxpointspec_act
248	if (decay(ks).gt.0.) then
249	do nage=1,nageclass
250	do l=1,nclassunc
251	! Mother output grid
252	do jy=0,numygrid-1
253	do ix=0,numxgrid-1
254	wetgridunc(ix,jy,ks,kp,l,nage)= &
255	wetgridunc(ix,jy,ks,kp,l,nage)* &
256	exp(-1.outstepdecay(ks))
257	drygridunc(ix,jy,ks,kp,l,nage)= &
258	drygridunc(ix,jy,ks,kp,l,nage)* &
259	exp(-1.outstepdecay(ks))
260	end do
261	end do
262	! Nested output grid
263	if (nested_output.eq.1) then
264	do jy=0,numygridn-1
265	do ix=0,numxgridn-1
266	wetgriduncn(ix,jy,ks,kp,l,nage)= &
267	wetgriduncn(ix,jy,ks,kp,l,nage)* &
268	exp(-1.outstepdecay(ks))
269	drygriduncn(ix,jy,ks,kp,l,nage)= &
270	drygriduncn(ix,jy,ks,kp,l,nage)* &
271	exp(-1.outstepdecay(ks))
272	end do
273	end do
274	endif
275	end do
276	end do
277	endif
278	end do
279	end do
280	endif
281
282	!!! CHANGE: These lines may be switched on to check the conservation
283	!!! of mass within FLEXPART
284	! if (itime.eq.loutnext) then
285	! do 247 ksp=1, nspec
286	! do 247 kp=1, maxpointspec_act
287	!47 xm(ksp,kp)=0.
288
289	! do 249 ksp=1, nspec
290	! do 249 j=1,numpart
291	! if (ioutputforeachrelease.eq.1) then
292	! kp=npoint(j)
293	! else
294	! kp=1
295	! endif
296	! if (itra1(j).eq.itime) then
297	! xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
298	! write(,) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
299	! endif
300	!49 continue
301	! do 248 ksp=1,nspec
302	! do 248 kp=1,maxpointspec_act
303	! xm_depw(ksp,kp)=0.
304	! xm_depd(ksp,kp)=0.
305	! do 248 nage=1,nageclass
306	! do 248 ix=0,numxgrid-1
307	! do 248 jy=0,numygrid-1
308	! do 248 l=1,nclassunc
309	! xm_depw(ksp,kp)=xm_depw(ksp,kp)
310	! + +wetgridunc(ix,jy,ksp,kp,l,nage)
311	!48 xm_depd(ksp,kp)=xm_depd(ksp,kp)
312	! + +drygridunc(ix,jy,ksp,kp,l,nage)
313	! do 246 ksp=1,nspec
314	!46 write(88,'(2i10,3e12.3)')
315	! + itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
316	! + (xm_depw(ksp,kp),kp=1,maxpointspec_act),
317	! + (xm_depd(ksp,kp),kp=1,maxpointspec_act)
318	! endif
319	!!! CHANGE
320
321
322
323	! Check whether concentrations are to be calculated
324	!**************************************************
325
326	if ((ldirectitime.ge.ldirectloutstart).and. &
327	(ldirectitime.le.ldirectloutend)) then ! add to grid
328	if (mod(itime-loutstart,loutsample).eq.0) then
329
330	! If we are exactly at the start or end of the concentration averaging interval,
331	! give only half the weight to this sample
332	!*****************************************************************************
333
334	if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
335	weight=0.5
336	else
337	weight=1.0
338	endif
339	outnum=outnum+weight
340	call conccalc(itime,weight)
341	endif
342
343
344	if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
345	call partoutput_short(itime) ! dump particle positions in extremely compressed format
346
347
348	! Output and reinitialization of grid
349	! If necessary, first sample of new grid is also taken
350	!*****************************************************
351
352	if ((itime.eq.loutend).and.(outnum.gt.0.)) then
353	if ((iout.le.3.).or.(iout.eq.5)) then
354	if (surf_only.ne.1) then
355	if (lnetcdfout.eq.1) then
356	#ifdef NETCDF_OUTPUT
357	call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
358	#endif
359	else
360	call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
361	endif
362	else
363	if (verbosity.eq.1) then
364	print*,'call concoutput_surf '
365	call system_clock(count_clock)
366	write(,) 'system clock',count_clock - count_clock0
367	endif
368	if (lnetcdfout.eq.1) then
369	#ifdef NETCDF_OUTPUT
370	call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
371	#endif
372	else
373	call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
374	if (verbosity.eq.1) then
375	print*,'called concoutput_surf '
376	call system_clock(count_clock)
377	write(,) 'system clock',count_clock - count_clock0
378	endif
379	endif
380	endif
381
382	if (nested_output .eq. 1) then
383	if (lnetcdfout.eq.0) then
384	if (surf_only.ne.1) then
385	call concoutput_nest(itime,outnum)
386	else
387	call concoutput_surf_nest(itime,outnum)
388	endif
389	else
390	#ifdef NETCDF_OUTPUT
391	if (surf_only.ne.1) then
392	call concoutput_nest_netcdf(itime,outnum)
393	else
394	call concoutput_surf_nest_netcdf(itime,outnum)
395	endif
396	#endif
397	endif
398	endif
399	outnum=0.
400	endif
401	if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
402	if (iflux.eq.1) call fluxoutput(itime)
403	!write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
404	write(*,46) float(itime)/3600,itime,numpart
405	45 format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES: Uncertainty: ',3f7.3)
406	46 format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
407	if (ipout.ge.1) call partoutput(itime) ! dump particle positions
408	loutnext=loutnext+loutstep
409	loutstart=loutnext-loutaver/2
410	loutend=loutnext+loutaver/2
411	if (itime.eq.loutstart) then
412	weight=0.5
413	outnum=outnum+weight
414	call conccalc(itime,weight)
415	endif
416
417
418	! Check, whether particles are to be split:
419	! If so, create new particles and attribute all information from the old
420	! particles also to the new ones; old and new particles both get half the
421	! mass of the old ones
422	!************************************************************************
423
424	if (ldirectitime.ge.ldirectitsplit) then
425	n=numpart
426	do j=1,numpart
427	if (ldirectitime.ge.ldirectitrasplit(j)) then
428	if (n.lt.maxpart) then
429	n=n+1
430	itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
431	itrasplit(n)=itrasplit(j)
432	itramem(n)=itramem(j)
433	itra1(n)=itra1(j)
434	idt(n)=idt(j)
435	npoint(n)=npoint(j)
436	nclass(n)=nclass(j)
437	xtra1(n)=xtra1(j)
438	ytra1(n)=ytra1(j)
439	ztra1(n)=ztra1(j)
440	uap(n)=uap(j)
441	ucp(n)=ucp(j)
442	uzp(n)=uzp(j)
443	us(n)=us(j)
444	vs(n)=vs(j)
445	ws(n)=ws(j)
446	cbt(n)=cbt(j)
447	do ks=1,nspec
448	xmass1(j,ks)=xmass1(j,ks)/2.
449	xmass1(n,ks)=xmass1(j,ks)
450	end do
451	endif
452	endif
453	end do
454	numpart=n
455	endif
456	endif
457	endif
458
459
460	if (itime.eq.ideltas) exit ! almost finished
461
462	! Compute interval since radioactive decay of deposited mass was computed
463	!************************************************************************
464
465	if (itime.lt.loutnext) then
466	ldeltat=itime-(loutnext-loutstep)
467	else ! first half of next interval
468	ldeltat=itime-loutnext
469	endif
470
471
472	! Loop over all particles
473	!************************
474
475	do j=1,numpart
476
477
478	! If integration step is due, do it
479	!**********************************
480
481	if (itra1(j).eq.itime) then
482
483	if (ioutputforeachrelease.eq.1) then
484	kp=npoint(j)
485	else
486	kp=1
487	endif
488	! Determine age class of the particle
489	itage=abs(itra1(j)-itramem(j))
490	do nage=1,nageclass
491	if (itage.lt.lage(nage)) exit
492	end do
493
494	! Initialize newly released particle
495	!***********************************
496
497	if ((itramem(j).eq.itime).or.(itime.eq.0)) &
498	call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
499	us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
500
501	! Memorize particle positions
502	!****************************
503
504	xold=real(xtra1(j))
505	yold=real(ytra1(j))
506	zold=ztra1(j)
507
508	! Integrate Lagevin equation for lsynctime seconds
509	!*************************************************
510
511	call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
512	us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
513	cbt(j))
514
515	! Calculate the gross fluxes across layer interfaces
516	!***************************************************
517
518	if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
519
520
521	! Determine, when next time step is due
522	! If trajectory is terminated, mark it
523	!**************************************
524
525	if (nstop.gt.1) then
526	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
527	itra1(j)=-999999999
528	else
529	itra1(j)=itime+lsynctime
530
531
532	! Dry deposition and radioactive decay for each species
533	! Also check maximum (of all species) of initial mass remaining on the particle;
534	! if it is below a threshold value, terminate particle
535	!*****************************************************************************
536
537	xmassfract=0.
538	do ks=1,nspec
539	if (decay(ks).gt.0.) then ! radioactive decay
540	decfact=exp(-real(abs(lsynctime))*decay(ks))
541	else
542	decfact=1.
543	endif
544
545	if (drydepspec(ks)) then ! dry deposition
546	drydeposit(ks)=xmass1(j,ks)prob(ks)decfact
547	xmass1(j,ks)=xmass1(j,ks)(1.-prob(ks))decfact
548	if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
549	drydeposit(ks)=drydeposit(ks)* &
550	exp(real(abs(ldeltat))*decay(ks))
551	endif
552	else ! no dry deposition
553	xmass1(j,ks)=xmass1(j,ks)*decfact
554	endif
555
556	if (mdomainfill.eq.0) then
557	if (xmass(npoint(j),ks).gt.0.) &
558	xmassfract=max(xmassfract,real(npart(npoint(j)))* &
559	xmass1(j,ks)/xmass(npoint(j),ks))
560	else
561	xmassfract=1.
562	endif
563	end do
564
565	if (xmassfract.lt.0.0001) then ! terminate all particles carrying less mass
566	itra1(j)=-999999999
567	if (verbosity.gt.0) then
568	print*,'terminated particle ',j,' for small mass'
569	endif
570	endif
571
572	! Sabine Eckhardt, June 2008
573	! don't create depofield for backward runs
574	if (drydep.AND.(ldirect.eq.1)) then
575	call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
576	if (nested_output.eq.1) then
577	call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
578	endif
579	endif
580
581	! Terminate trajectories that are older than maximum allowed age
582	!***************************************************************
583
584	if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
585	if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
586	itra1(j)=-999999999
587	if (verbosity.gt.0) then
588	print*,'terminated particle ',j,' for age'
589	endif
590	endif
591	endif
592
593	endif
594
595	end do
596
597	end do
598
599
600	! Complete the calculation of initial conditions for particles not yet terminated
601	!*****************************************************************************
602
603	do j=1,numpart
604	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
605	end do
606
607	if (ipout.eq.2) call partoutput(itime) ! dump particle positions
608
609	if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
610
611	close(104)
612
613	! De-allocate memory and end
614	!***************************
615
616	if (iflux.eq.1) then
617	deallocate(flux)
618	endif
619	if (ohrea.eqv..TRUE.) then
620	deallocate(OH_field,OH_field_height)
621	endif
622	if (ldirect.gt.0) then
623	deallocate(drygridunc,wetgridunc)
624	endif
625	deallocate(gridunc)
626	deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
627	deallocate(ireleasestart,ireleaseend,npart,kindz)
628	deallocate(xmasssave)
629	if (nested_output.eq.1) then
630	deallocate(orooutn, arean, volumen)
631	if (ldirect.gt.0) then
632	deallocate(griduncn,drygriduncn,wetgriduncn)
633	endif
634	endif
635	deallocate(outheight,outheighthalf)
636	deallocate(oroout,area,volume)
637
638	end subroutine timemanager
639

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

Download in other formats:

Original Format