Context Navigation

timemanager.f90 @ 8

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

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

Download in other formats:

Original Format