Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

timemanager.f90 @ 14

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

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

Download in other formats:

Original Format