Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

timemanager_serial.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago
sources for flexwrf v3.1
File size: 28.7 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23
24	subroutine timemanager(mts)
25
26	!*******************************************************************************
27	! *
28	! Handles the computation of trajectories, i.e. determines which *
29	! trajectories have to be computed at what time. *
30	! Manages dry+wet deposition routines, radioactive decay and the computation *
31	! of concentrations. *
32	! *
33	! Author: A. Stohl *
34	! *
35	! 20 May 1996 *
36	! *
37	! Dec 2005, J. Fast - Only call conccalc & concoutput when (iout.ge.1) *
38	! Aug 2007, W. Wang - call KFeta convection scheme (lconvection=2or3)
39	! Note, backward is unavailabe for lconvection=2
40	! Mar 2012, J. Brioude: modifications to handle openmp and mpi *
41	!*******************************************************************************
42	! Changes, Bernd C. Krueger, Feb. 2001: *
43	! Call of convmix when new windfield is read *
44	!------------------------------------ *
45	! Changes Petra Seibert, Sept 2002 *
46	! fix wet scavenging problem *
47	! Code may not be correct for decay of deposition! *
48	! Changes Petra Seibert, Nov 2002 *
49	! call convection BEFORE new fields are read in BWD mode *
50	! Changes Caroline Forster, Feb 2005
51	! new interface between flexpart and convection scheme
52	! Emanuel's latest subroutine convect43c.f is used
53	!*******************************************************************************
54	! *
55	! Variables: *
56	! DEP .true. if either wet or dry deposition is switched on *
57	! decay(maxspec) [1/s] decay constant for radioactive decay *
58	! DRYDEP .true. if dry deposition is switched on *
59	! ideltas [s] modelling period *
60	! jtime [s] actual temporal position of calculation *
61	! ldeltat [s] time since computation of radioact. decay of depositions *
62	! loutaver [s] averaging period for concentration calculations *
63	! loutend [s] end of averaging for concentration calculations *
64	! loutnext [s] next time at which output fields shall be centered *
65	! loutsample [s] sampling interval for averaging of concentrations *
66	! loutstart [s] start of averaging for concentration calculations *
67	! loutstep [s] time interval for which concentrations shall be calculated*
68	! npoint(maxpart) index, which starting point the trajectory has *
69	! starting positions of trajectories *
70	! nstop serves as indicator for fate of particles *
71	! in the particle loop *
72	! nstop1 serves as indicator for wind fields (see getfields) *
73	! outnum number of samples for each concentration calculation *
74	! outnum number of samples for each concentration calculation *
75	! prob probability of absorption at ground due to dry deposition *
76	! WETDEP .true. if wet deposition is switched on *
77	! weight weight for each concentration sample (1/2 or 1) *
78	! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to turbulence *
79	! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to interpolation *
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	! include 'includepar'
89	! include 'includecom'
90	use unc_mod
91	use point_mod
92	! use xmass_mod
93	use flux_mod
94	use outg_mod
95	use oh_mod
96	use par_mod
97	use com_mod
98	use mt_stream
99
100	! use ran_mod
101	! use interpol_mod
102
103	implicit none
104
105
106	integer :: ix,jy,j,ks,kp,l,n,jtime,nstop,nstop1
107	! integer :: ksp
108	integer :: loutnext,loutstart,loutend,jj,chunksize
109	!,chunksize2
110	integer :: chunksize3,omp_get_num_threads
111	integer :: ldeltat,itage,nage,th_itra1,i
112	real :: outnum,weight,prob(maxspec),nrand,decfact
113	! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
114	! real :: us(maxpart),vs(maxpart),ws(maxpart)
115	! integer(kind=2) :: cbt(maxpart)
116	! real,allocatable, dimension (:) :: uap,ucp,uzp
117	! real,allocatable, dimension (:) :: us,vs,ws
118	! integer(kind=2),allocatable, dimension (:) :: cbt
119	real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
120	real :: drygridtotalunc,xold,yold,zold,xmassfract
121	! integer j,k,l,n,jtime,nstop,nstop1
122	! integer loutnext,loutstart,loutend
123	! integer ix,jy,ldeltat,itage,nage
124	! real outnum,weight,prob(maxspec)
125	! real uap(maxpart),ucp(maxpart),uzp(maxpart),decfact
126	! real us(maxpart),vs(maxpart),ws(maxpart),cbt(maxpart)
127	! real drydeposit(maxspec),gridtotalunc,wetgridtotalunc
128	! real drygridtotalunc,xold,yold,zold
129	! real xm,xm1
130
131
132	integer :: th_npoint,th_idt,th_itramem,jdeb,jfin,stat,th_nclass
133	integer,save :: cpt(maxomp)=0
134	! integer,save :: cpt(24)=0
135	real(kind=dp) :: th_xtra1,th_ytra1
136	real :: th_ztra1,th_uap,th_ucp,th_uzp
137	real :: th_us,th_vs,th_ws,ran3
138	integer(kind=2) :: th_cbt
139	integer :: from
140
141	real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2
142	integer :: ixp,jyp,ngrid,indz,indzp,nbp,jj2,ii,offset
143	logical :: depoindicator(maxspec)
144	logical,save :: indzindicator(nzmax)
145	real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
146	real :: sigw,dsigwdz,dsigw2dz,th_xmass1(maxspec)
147	real :: start, finish
148	real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
149	real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
150	real :: rhoprof(nzmax),rhogradprof(nzmax)
151	real :: tkeprof(nzmax),pttprof(nzmax)
152	real :: u,v,w,usig,vsig,wsig,pvi
153	integer*4 :: now(3)
154	integer :: ttime,cpttra
155	! integer, dimension(MPI_STATUS_SIZE) :: status
156	integer :: myid,ntasks,ierr,islave,tag2,ompid,n_threads,tag3,i_omp
157	type (mt_state) :: mts (0: MAX_STREAM)
158	!************************
159
160	!JB
161	! call MPI_COMM_RANK ( MPI_COMM_WORLD, myid, ierr )
162	! call MPI_COMM_SIZE ( MPI_COMM_WORLD, ntasks, ierr )
163	! myid gives the info on the node id
164	ntasks=1
165	myid=0
166	loutnext=loutstep/2
167	outnum=0.
168	loutstart=loutnext-loutaver/2
169	loutend=loutnext+loutaver/2
170
171	! if (myid.eq.0) then
172	allocate(uap(maxpart) ,stat=stat)
173	allocate(ucp(maxpart) ,stat=stat)
174	allocate(uzp(maxpart) ,stat=stat)
175	allocate(us(maxpart) ,stat=stat)
176	allocate(vs(maxpart) ,stat=stat)
177	allocate(ws(maxpart) ,stat=stat)
178	allocate(cbt(maxpart) ,stat=stat)
179	! endif
180
181	!**********************************************************************
182	! Loop over the whole modelling period in time steps of mintime seconds
183	!**********************************************************************
184
185	! print*,'time',myid,ideltas,lsynctime
186	do jtime=0,ideltas,lsynctime
187
188
189	! print*,'jtime',jtime
190	! Computation of wet deposition every lsynctime seconds
191	! maybe wet depo frequency can be relaxed later but better be on safe side
192	! wetdepo must be called BEFORE new fields are read in but should not
193	! be called in the very beginning before any fields are loaded, or
194	! before particles are in the system
195	! Code may not be correct for decay of deposition
196	! changed by Petra Seibert 9/02
197	!********************************************************************
198
199	if (WETDEP .and. jtime .ne. 0 .and. numpart .gt. 0) &
200	call wetdepo(jtime,lsynctime,loutnext)
201
202	if (OHREA .and. jtime .ne. 0 .and. numpart .gt. 0) &
203	call ohreaction(jtime,lsynctime,loutnext)
204
205	! compute convection for backward runs
206	!*************************************
207
208	! if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(jtime.lt.0))
209	! & call convmix(jtime)
210
211	if ((ldirect.eq.-1).and.(jtime.lt.0)) then
212	if (lconvection .eq. 1) call convmix(jtime)
213	if (lconvection .eq. 2 .or. lconvection .eq. 3) &
214	call convmix_kfeta(jtime)
215	endif
216
217	! Get necessary wind fields if not available
218	!*******************************************
219
220	! call itime(now)
221	! ttime=now(1)3600+now(2)60+now(3)
222	call cpu_time(start)
223	call getfields(jtime,nstop1)
224	if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
225	! call itime(now)
226	! ttime=now(1)3600+now(2)60+now(3)-ttime
227	call cpu_time(finish)
228	! print*,'read wind time',ttime
229
230	! Release particles
231	!******************
232
233	!JB
234	if (myid.eq.0) then ! I let only the master thread releasing the particles and calculate the output
235	! call itime(now)
236	call cpu_time(start)
237	if (mdomainfill.ge.1) then
238	if (jtime.eq.0) then
239	call init_domainfill()
240	else
241	call boundcond_domainfill(jtime,loutend)
242	endif
243	else
244	if (numpoint_option.eq.0) then
245	call releaseparticles_irreg(jtime)
246	elseif (numpoint_option.eq.1) then
247	! print*,'avant release'
248	call releaseparticles_reg(jtime)
249	endif
250	endif
251	! do i=1,numpart
252	! print*,'ipart 2',myid,i,ztra1(i)
253	! enddo
254	! print*,'test rel',npoint(1),npoint(2),npoint(3)
255
256	! print*,'test rel1',npoint(5139),npoint(6002),npoint(100003)
257	! Compute convective mixing for forward runs
258	! for backward runs it is done before next windfield is read in
259	!**************************************************************
260
261	! if ((ldirect.eq.1).and.(lconvection.eq.1)) &
262	! call convmix(jtime)
263
264	if (ldirect.eq.1) then
265	if (lconvection.eq.1)call convmix(jtime)
266	if (lconvection.eq.2 .or. lconvection .eq. 3) &
267	call convmix_kfeta(jtime)
268	endif
269
270	! If middle of averaging period of output fields is reached, accumulated
271	! deposited mass radioactively decays
272	!***********************************************************************
273
274	if (DEP.and.(jtime.eq.loutnext).and.(ldirect.gt.0)) then
275	do ks=1,nspec
276	do kp=1,maxpointspec_act
277	if (decay(ks).gt.0.) then
278	do nage=1,nageclass
279	do l=1,nclassunc
280	! Mother output grid
281	do jy=0,numygrid-1
282	do ix=0,numxgrid-1
283	wetgridunc(ix,jy,ks,kp,l,nage)= &
284	wetgridunc(ix,jy,ks,kp,l,nage)* &
285	exp(-1.outstepdecay(ks))
286	drygridunc(ix,jy,ks,kp,l,nage)= &
287	drygridunc(ix,jy,ks,kp,l,nage)* &
288	exp(-1.outstepdecay(ks))
289	end do
290	end do
291	! Nested output grid
292	if (nested_output.eq.1) then
293	do jy=0,numygridn-1
294	do ix=0,numxgridn-1
295	wetgriduncn(ix,jy,ks,kp,l,nage)= &
296	wetgriduncn(ix,jy,ks,kp,l,nage)* &
297	exp(-1.outstepdecay(ks))
298	drygriduncn(ix,jy,ks,kp,l,nage)= &
299	drygriduncn(ix,jy,ks,kp,l,nage)* &
300	exp(-1.outstepdecay(ks))
301	end do
302	end do
303	endif
304	end do
305	end do
306	endif
307	end do
308	end do
309	endif
310
311	!!! CHANGE: These lines may be switched on to check the conservation
312	!!! of mass within FLEXPART
313
314	! if (mod(jtime,loutsample).eq.0) then
315	! xm=0.
316	! xm1=0.
317	! do 247 j=1,numpart
318	!47 if (itra1(j).eq.jtime) xm1=xm1+xmass1(j,1)
319	! xm=xm1
320	! do 248 nage=1,nageclass
321	! do 248 ix=0,numxgrid-1
322	! do 248 jy=0,numygrid-1
323	! do 248 l=1,nclassunc
324	!48 xm=xm+wetgridunc(ix,jy,1,l,nage)+drygridunc(ix,jy,1,l,nage)
325	! write(*,'(i6,4f10.3)') jtime,xm,xm1
326	! endif
327	!!! CHANGE
328
329
330	! Check whether concentrations are to be calculated
331	!**************************************************
332
333	if ((ldirectjtime.ge.ldirectloutstart).and. &
334	(ldirectjtime.le.ldirectloutend)) then ! add to grid
335	if (mod(jtime-loutstart,loutsample).eq.0) then
336
337	! If we are exactly at the start or end of the concentration averaging interval,
338	! give only half the weight to this sample
339	!*******************************************************************************
340
341	if ((jtime.eq.loutstart).or.(jtime.eq.loutend)) then
342	weight=0.5
343	else
344	weight=1.0
345	endif
346	outnum=outnum+weight
347	if(iout.ge.1) then
348	if (outgrid_option.eq.0) then
349	call conccalc_irreg(jtime,weight)
350	elseif (outgrid_option.eq.1) then
351	call conccalc_reg(jtime,weight)
352	endif
353	endif
354	endif
355
356
357	! if ((mquasilag.eq.1).and.(jtime.eq.(loutstart+loutend)/2)) &
358	! call partoutput_short(jtime) ! dump particle positions in extremely compressed format
359
360
361	! Output and reinitialization of grid
362	! If necessary, first sample of new grid is also taken
363	!*****************************************************
364
365	if ((jtime.eq.loutend).and.(outnum.gt.0.)) then
366	! print*,'iout',iout,ipout,outgrid_option
367	if ((iout.le.3.).or.(iout.eq.5)) then
368	if(iout.ge.1) then
369	if (outgrid_option.eq.0) then
370	call concoutput_irreg(jtime,outnum,gridtotalunc, &
371	wetgridtotalunc,drygridtotalunc)
372	if (nested_output.eq.1) call concoutput_nest_irreg(jtime,outnum)
373	elseif (outgrid_option.eq.1) then
374	call concoutput_reg(jtime,outnum,gridtotalunc, &
375	wetgridtotalunc,drygridtotalunc)
376	if (nested_output.eq.1) call concoutput_nest_reg(jtime,outnum)
377	endif
378	endif
379
380	! if (nested_output.eq.1.and.iout.ge.1)
381	! + call concoutput_nest(jtime,outnum)
382	outnum=0.
383	endif
384	if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(jtime)
385	if (iflux.eq.1) call fluxoutput(jtime)
386	write(*,45) jtime,numpart,gridtotalunc,wetgridtotalunc, &
387	drygridtotalunc
388	45 format(i9,' SECONDS SIMULATED: ',i9, &
389	' PARTICLES: Uncertainty: ',3f7.3)
390	if (ipout.ge.1) call partoutput(jtime) ! dump particle positions
391	loutnext=loutnext+loutstep
392	loutstart=loutnext-loutaver/2
393	loutend=loutnext+loutaver/2
394	if (jtime.eq.loutstart) then
395	weight=0.5
396	outnum=outnum+weight
397	if(iout.ge.1) then
398	if (outgrid_option.eq.0) then
399	call conccalc_irreg(jtime,weight)
400	elseif (outgrid_option.eq.1) then
401	call conccalc_reg(jtime,weight)
402	endif
403	endif
404	endif
405
406
407	! Check, whether particles are to be split:
408	! If so, create new particles and attribute all information from the old
409	! particles also to the new ones; old and new particles both get half the
410	! mass of the old ones
411	!************************************************************************
412
413	if (ldirectjtime.ge.ldirectitsplit) then
414	n=numpart
415	do j=1,numpart
416	if (ldirectjtime.ge.ldirectitrasplit(j)) then
417	if (n.lt.maxpart) then
418	n=n+1
419	itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
420	itrasplit(n)=itrasplit(j)
421	itramem(n)=itramem(j)
422	itra1(n)=itra1(j)
423	idt(n)=idt(j)
424	npoint(n)=npoint(j)
425	nclass(n)=nclass(j)
426	xtra1(n)=xtra1(j)
427	ytra1(n)=ytra1(j)
428	ztra1(n)=ztra1(j)
429	uap(n)=uap(j)
430	ucp(n)=ucp(j)
431	uzp(n)=uzp(j)
432	us(n)=us(j)
433	vs(n)=vs(j)
434	ws(n)=ws(j)
435	cbt(n)=cbt(j)
436	do ks=1,nspec
437	xmass1(j,ks)=xmass1(j,ks)/2.
438	xmass1(n,ks)=xmass1(j,ks)
439	end do
440	endif
441	endif
442	end do
443	numpart=n
444	endif
445	endif
446	endif
447
448
449
450
451	! Loop over all particles
452	!************************
453
454
455	! chunksize=int(numpart/ntasks)+1 !if sent homogeneously
456	! call itime(now)
457	! ttime=now(1)3600+now(2)60+now(3)-ttime
458	call cpu_time(finish)
459
460	! print*,'processing time',ttime
461	endif !over myid
462	!JB
463	! at this stage, I assume that each node has the same shared memory because they run getfields.
464	! now we need to split the trajectories into pieces for each node
465	! if (myid.eq.0) then
466
467	if (jtime.eq.ideltas) exit
468
469	! Compute interval since radioactive decay of deposited mass was computed
470	!************************************************************************
471
472	if (jtime.lt.loutnext) then
473	ldeltat=jtime-(loutnext-loutstep)
474	else ! first half of next interval
475	ldeltat=jtime-loutnext
476	endif
477
478
479	! if (myid.eq.0) then
480	! call itime(now)
481	! ttime=now(1)3600+now(2)60+now(3)
482	! do ii=1,ntasks-1
483	! call MPI_SEND(chunksize,1, MPI_INTEGER, ii,3001, MPI_COMM_WORLD, ierr)
484	! call MPI_SEND(numpart,1, MPI_INTEGER, ii,3002, MPI_COMM_WORLD, ierr)
485	! enddo
486	! else
487	! call MPI_RECV(chunksize,1, MPI_INTEGER, 0,3001, MPI_COMM_WORLD,status, ierr)
488	! call MPI_RECV(numpart,1, MPI_INTEGER, 0,3002, MPI_COMM_WORLD,status, ierr)
489	! endif
490	! print*,'numpart',numpart
491
492	! call itime(now)
493	! ttime=now(1)3600+now(2)60+now(3)
494
495	! initialize the temporary drydeposition grid
496
497	if (DRYDEP.and.ldirect.gt.0) then
498	do ks=1,nspec
499	do kp=1,maxpointspec_act
500	do nage=1,nageclass
501	do jy=0,numygrid-1
502	do ix=0,numxgrid-1
503	do l=1,nclassunc
504	drygridunc2(ix,jy,ks,kp,l,nage)=0.
505	end do
506	end do
507	end do
508	if (nested_output.eq.1) then
509	do jy=0,numygridn-1
510	do ix=0,numxgridn-1
511	do l=1,nclassunc
512	drygriduncn2(ix,jy,ks,kp,l,nage)=0.
513	end do
514	end do
515	end do
516	endif
517	end do
518	end do
519	end do
520	endif
521
522
523	call cpu_time(start)
524	! chunksize3=int(chunksize2/omp_get_num_threads())+1
525	! chunksize3=int(real(chunksize2)/real(omp_get_num_threads())/20.)+1 !more efficient
526
527	! ompid=omp_get_num_threads()
528	ompid=0
529
530	cpttra=0
531	! print*,'chunksi',chunksize2,myid
532	if (numpart.gt.0 ) then
533	! print*,'test rel2',npoint(5139),npoint(6002),npoint(100003)
534	! do jj=1,numpart
535	! do jj=numpart,1,-1
536	! print*,jj
537	do jj=1,numpart
538
539	! If integration step is due, do it
540	!**********************************
541	if (itra1(jj).eq.jtime) then
542	cpttra=cpttra+1
543	if (ioutputforeachrelease.eq.1) then
544	kp=npoint(jj)
545	else
546	kp=1
547	endif
548
549	! Determine age class of the particle
550	itage=abs(itra1(jj)-itramem(jj))
551	do nage=1,nageclass
552	if (itage.lt.lage(nage)) exit
553	enddo
554
555	nbp=1
556	if ((itramem(jj).eq.jtime).or.(jtime.eq.0)) &
557	call initialize(jtime,idt(jj),uap(jj),ucp(jj),uzp(jj), &
558	us(jj),vs(jj),ws(jj),xtra1(jj),ytra1(jj),ztra1(jj),cbt(jj), &
559	! call initialize(jtime,th_idt,th_uap,th_ucp,th_uzp, &
560	! th_us,th_vs,th_ws,th_xtra1,th_ytra1,th_ztra1,th_cbt, &
561	ngrid,depoindicator,indzindicator,cpt(nbp),ompid,myid,1,mts )
562
563	! print*,'after',th_xtra1,th_ytra1,th_ztra1
564	! Memorize particle positions
565	!****************************
566
567	xold=xtra1(jj)
568	yold=ytra1(jj)
569	zold=ztra1(jj)
570	! xold=th_xtra1
571	! yold=th_ytra1
572	! zold=th_ztra1
573	! Integrate Lagevin equation for lsynctime seconds
574	!*************************************************
575	! write(,)'numpart,jtime, particle #=',numpart,jtime,j
576
577	call advance(jtime,npoint(jj),idt(jj),uap(jj),ucp(jj),uzp(jj),us(jj), &
578	vs(jj),ws(jj),nstop,xtra1(jj),ytra1(jj),ztra1(jj),prob,cbt(jj), &
579	! call advance(jtime,th_npoint,th_idt,th_uap,th_ucp,th_uzp, &
580	! th_us,th_vs,th_ws,nstop,th_xtra1,&
581	! th_ytra1,th_ztra1,prob,th_cbt, &
582	ngrid,depoindicator,indzindicator,cpt(nbp),ompid,myid,1,mts )
583	! if (jj.eq.103) print*,'aft',th_xtra1,th_ytra1,th_ztra1
584	! Calculate the gross fluxes across layer interfaces
585	!***************************************************
586
587
588	if (iflux.eq.1) call calcfluxes(nage,jj,xold,yold,zold)
589
590	! if (jj.lt.5) print*,'coord after',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks)
591
592	! Determine, when next time step is due
593	! If trajectory is terminated, mark it
594	!**************************************
595
596	if (nstop.gt.1) then
597	if (linit_cond.ge.1) call initial_cond_calc(jtime,jj)
598	itra1(jj)=-999999999
599	! th_itra1=-999999999
600	else
601	itra1(jj)=jtime+lsynctime
602	! th_itra1=jtime+lsynctime
603
604
605	! if (jj.lt.5) print*,'coord after2',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks)
606	! Dry deposition and radioactive decay for each species
607	!******************************************************
608	xmassfract=0.
609
610	do ks=1,nspec
611	if (decay(ks).gt.0.) then ! radioactive decay
612	decfact=exp(-real(abs(lsynctime))*decay(ks))
613	else
614	decfact=1.
615	endif
616
617	if (DRYDEPSPEC(ks)) then ! dry deposition
618	drydeposit(ks)=xmass1(jj,ks)prob(ks)decfact
619	! drydeposit(ks)=th_xmass1(ks)prob(ks)decfact
620	xmass1(jj,ks)=xmass1(jj,ks)(1.-prob(ks))decfact
621	! th_xmass1(ks)=th_xmass1(ks)(1.-prob(ks))decfact
622	if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
623	drydeposit(ks)=drydeposit(ks)* &
624	exp(real(abs(ldeltat))*decay(ks))
625	endif
626	else ! no dry deposition
627	xmass1(jj,ks)=xmass1(jj,ks)*decfact
628	! th_xmass1(ks)=th_xmass1(ks)*decfact
629	endif
630	! if (jj.lt.5) print*,'coord after3',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks),xmass(th_npoint,1)
631
632	if (mdomainfill.eq.0) then
633	if (xmass(npoint(jj),ks).gt.0.) &
634	xmassfract=max(xmassfract,real(npart(npoint(jj)))* &
635	! xmassfract=max(xmassfract,real(npart(th_npoint))* &
636	xmass1(jj,ks)/xmass(npoint(jj),ks))
637	! th_xmass1(ks)/xmass(th_npoint,ks))
638	else
639	xmassfract=1.
640	endif
641
642	end do
643
644	if (xmassfract.lt.0.000001) then ! terminate all particles carrying less mass
645	itra1(jj)=-999999999
646	! th_itra1=-999999999
647	endif
648
649	! Sabine Eckhardt, June 2008
650	! don't create depofield for backward runs
651	if (DRYDEP.AND.(ldirect.eq.1)) then
652	call drydepokernel(nclass(jj),drydeposit,real(xtra1(jj)), &
653	! call drydepokernel(th_nclass,drydeposit,real(th_xtra1), &
654	real(ytra1(jj)),itage,nage,kp)
655	! real(th_ytra1),itage,nage,kp)
656	if (nested_output.eq.1) call drydepokernel_nest( &
657	nclass(jj),drydeposit,real(xtra1(jj)),real(ytra1(jj)), &
658	! th_nclass,drydeposit,real(th_xtra1),real(th_ytra1), &
659	itage,nage,kp)
660	endif
661
662	! Terminate trajectories that are older than maximum allowed age
663	!***************************************************************
664
665	if (abs(itra1(jj)-itramem(jj)).ge.lage(nageclass)) then
666	! if (abs(th_itra1-th_itramem).ge.lage(nageclass)) then
667	if (linit_cond.ge.1) &
668	call initial_cond_calc(jtime+lsynctime,jj)
669	itra1(jj)=-999999999
670	! th_itra1=-999999999
671	endif
672	endif
673	!! print*,xtra1(j),th_xtra1,OMP_GET_THREAD_NUM()
674
675	endif
676
677	end do !loop over particles
678
679	endif
680
681
682
683	! call itime(now)
684	! ttime=now(1)3600+now(2)60+now(3)-ttime
685	call cpu_time(finish)
686	! print*,'time',ttime,cpttra,myid,OMP_GET_THREAD_NUM()
687	if (option_verbose.eq.1) then
688	print*,'time',finish-start,cpttra,myid,ompid
689	endif
690
691	! update the drydepo
692	if (DRYDEP.and.ldirect.gt.0) then
693	do ks=1,nspec
694	do kp=1,maxpointspec_act
695	do nage=1,nageclass
696
697	do jy=0,numygrid-1
698	do ix=0,numxgrid-1
699	do l=1,nclassunc
700	drygridunc(ix,jy,ks,kp,l,nage)=drygridunc(ix,jy,ks,kp,l,nage) &
701	+drygridunc2(ix,jy,ks,kp,l,nage)
702	end do
703	end do
704	end do
705	if (nested_output.eq.1) then
706	do jy=0,numygridn-1
707	do ix=0,numxgridn-1
708	do l=1,nclassunc
709	drygriduncn(ix,jy,ks,kp,l,nage)=drygriduncn(ix,jy,ks,kp,l,nage) &
710	+drygriduncn2(ix,jy,ks,kp,l,nage)
711	end do
712	end do
713	end do
714	endif
715	end do
716	end do
717	end do
718
719	endif
720
721	end do !loop over time
722
723
724	! Complete the calculation of initial conditions for particles not yet terminated
725	!*****************************************************************************
726
727	do j=1,numpart
728	if (linit_cond.ge.1) call initial_cond_calc(jtime,j)
729	end do
730
731	if (ipout.eq.2) call partoutput(jtime) ! dump particle positions
732
733	if (linit_cond.ge.1) call initial_cond_output(jtime) ! dump initial cond. field
734
735	close(104)
736
737	! De-allocate memory and end
738	!***************************
739
740	if (iflux.eq.1) then
741	deallocate(flux)
742	endif
743	if (OHREA.eqv..TRUE.) then
744	deallocate(OH_field,OH_field_height)
745	endif
746	deallocate(gridunc)
747	deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
748	deallocate(ireleasestart,ireleaseend,npart,kindz)
749	! deallocate(xmasssave)
750	if (myid.eq.0) then
751	if (nested_output.eq.1) then
752	deallocate(orooutn, arean, volumen)
753	if (ldirect.gt.0) then
754	deallocate(griduncn,drygriduncn,wetgriduncn,drygriduncn2)
755	endif
756	endif
757	if (ldirect.gt.0) then
758	if (allocated(drygridunc)) deallocate(drygridunc)
759	if (allocated(wetgridunc)) deallocate(wetgridunc)
760	if (allocated(drygridunc2)) deallocate(drygridunc2)
761	if (allocated(drygriduncn2)) deallocate(drygriduncn2)
762	endif
763	deallocate(outheight,outheighthalf)
764	deallocate(oroout, area, volume)
765	endif
766	end subroutine timemanager
767
768
769

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

Download in other formats:

Original Format