Context Navigation

timemanager_mpi.f90 @ 0f7835d

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

Last change on this file since 0f7835d was 0f7835d, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago
Changed a variable in domain-filling option (MPI). Deleted an unused function from mpi_mod.f90
Property mode set to `100644`
File size: 32.3 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	! Changes Espen Sollum 2014 *
50	! MPI version *
51	! Variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod *
52	!*****************************************************************************
53	! *
54	! Variables: *
55	! dep .true. if either wet or dry deposition is switched on *
56	! decay(maxspec) [1/s] decay constant for radioactive decay *
57	! drydep .true. if dry deposition is switched on *
58	! ideltas [s] modelling period *
59	! itime [s] actual temporal position of calculation *
60	! ldeltat [s] time since computation of radioact. decay of depositions*
61	! loutaver [s] averaging period for concentration calculations *
62	! loutend [s] end of averaging for concentration calculations *
63	! loutnext [s] next time at which output fields shall be centered *
64	! loutsample [s] sampling interval for averaging of concentrations *
65	! loutstart [s] start of averaging for concentration calculations *
66	! loutstep [s] time interval for which concentrations shall be *
67	! 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	! memstat additional indicator for wind fields (see getfields) *
74	! outnum number of samples for each concentration calculation *
75	! prob probability of absorption at ground due to dry *
76	! deposition *
77	! wetdep .true. if wet deposition is switched on *
78	! weight weight for each concentration sample (1/2 or 1) *
79	! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to *
80	! turbulence *
81	! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to inter- *
82	! polation *
83	! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = *
84	! spatial positions of trajectories *
85	! *
86	! Constants: *
87	! maxpart maximum number of trajectories *
88	! *
89	!*****************************************************************************
90
91	use unc_mod
92	use point_mod
93	use xmass_mod
94	use flux_mod
95	use outg_mod
96	use oh_mod
97	use par_mod
98	use com_mod
99	use mpi_mod
100	use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
101	&concoutput_surf_netcdf,concoutput_surf_nest_netcdf
102
103	implicit none
104
105	logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
106	integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 !,mind
107	! integer :: ksp
108	integer :: ip
109	integer :: loutnext,loutstart,loutend
110	integer :: ix,jy,ldeltat,itage,nage
111	integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
112	integer :: numpart_tot_mpi ! for summing particles on all processes
113	real :: outnum,weight,prob(maxspec)
114	real :: decfact
115
116	real(sp) :: gridtotalunc
117	real(dep_prec) :: drygridtotalunc=0_dep_prec,wetgridtotalunc=0_dep_prec,&
118	& drydeposit(maxspec)=0_dep_prec
119	real :: xold,yold,zold,xmassfract
120	real, parameter :: e_inv = 1.0/exp(1.0)
121
122	! Measure time spent in timemanager
123	if (mp_measure_time) call mpif_mtime('timemanager',0)
124
125	! First output for time 0
126	!************************
127
128	loutnext=loutstep/2
129	outnum=0.
130	loutstart=loutnext-loutaver/2
131	loutend=loutnext+loutaver/2
132
133
134	!**********************************************************************
135	! Loop over the whole modelling period in time steps of mintime seconds
136	!**********************************************************************
137
138
139	if (lroot.or.mp_dev_mode) then
140	! write(,45) itime,numpartmp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
141	write(,46) float(itime)/3600,itime,numpartmp_partgroup_np
142
143	if (verbosity.gt.0) then
144	write (,) 'timemanager> starting simulation'
145	end if
146	end if ! (lroot)
147
148	!CGZ-lifetime: set lifetime to 0
149	! checklifetime(:,:)=0
150	! species_lifetime(:,:)=0
151	! print*, 'Initialized lifetime'
152	!CGZ-lifetime: set lifetime to 0
153
154
155
156	do itime=0,ideltas,lsynctime
157
158	! Computation of wet deposition, OH reaction and mass transfer
159	! between two species every lsynctime seconds
160	! maybe wet depo frequency can be relaxed later but better be on safe side
161	! wetdepo must be called BEFORE new fields are read in but should not
162	! be called in the very beginning before any fields are loaded, or
163	! before particles are in the system
164	! Code may not be correct for decay of deposition
165	! changed by Petra Seibert 9/02
166	!********************************************************************
167
168	if (mp_dev_mode) write(,) 'myid, itime: ',mp_pid,itime
169
170	if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
171	if (verbosity.gt.0) then
172	write (,) 'timemanager> call wetdepo'
173	endif
174
175	if (mp_measure_time) call mpif_mtime('wetdepo',0)
176
177	! readwind process skips this step
178	if (.not.(lmpreader.and.lmp_use_reader)) then
179	call wetdepo(itime,lsynctime,loutnext)
180	end if
181
182	if (mp_measure_time) call mpif_mtime('wetdepo',1)
183	end if
184
185	if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) then
186	! readwind process skips this step
187	if (.not.(lmpreader.and.lmp_use_reader)) then
188	call ohreaction(itime,lsynctime,loutnext)
189	endif
190	end if
191
192	if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
193	stop 'associated species not yet implemented!'
194	! call transferspec(itime,lsynctime,loutnext)
195	endif
196
197
198	! compute convection for backward runs
199	!*************************************
200
201	if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
202	if (verbosity.gt.0) then
203	write (,) 'timemanager> call convmix -- backward'
204	endif
205	! readwind process skips this step
206	if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime)
207	endif
208
209	! Get necessary wind fields if not available
210	!*******************************************
211	if (verbosity.gt.0 .and. lmpreader) then
212	write (,) 'timemanager> call getfields'
213	endif
214
215	! This time measure includes reading/MPI communication (for the reader process),
216	! or MPI communication time only (for other processes)
217	if (mp_measure_time) call mpif_mtime('getfields',0)
218
219	call getfields(itime,nstop1,memstat)
220
221	if (mp_measure_time) call mpif_mtime('getfields',1)
222
223
224	! Broadcast fields to all MPI processes
225	! Skip if all processes have called getfields or if no new fields
226	!*****************************************************************
227
228	if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',0)
229
230	! Two approaches to MPI getfields is implemented:
231	! Version 1 (lmp_sync=.true.) uses a read-ahead process where send/recv is done
232	! in sync at start of each new field time interval
233	!
234	! Version 2 (lmp_sync=.false.) is for holding three fields in memory. Uses a
235	! read-ahead process where sending/receiving of the 3rd fields is done in
236	! the background in parallel with performing computations with fields 1&2
237	!********************************************************************************
238
239	if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then
240	call mpif_gf_send_vars(memstat)
241	if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
242	! Version 2 (lmp_sync=.false.) is also used whenever 2 new fields are
243	! read (as at first time step), in which case async send/recv is impossible.
244	else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
245	call mpif_gf_send_vars(memstat)
246	if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
247	end if
248
249	if (.not.lmp_sync) then
250
251	! Reader process:
252	if (memstat.gt.0..and.memstat.lt.32.and.lmp_use_reader.and.lmpreader) then
253	if (mp_dev_mode) write(,) 'Reader process: calling mpif_gf_send_vars_async'
254	call mpif_gf_send_vars_async(memstat)
255	if (numbnests>0) call mpif_gf_send_vars_nest_async(memstat)
256	end if
257
258	! Completion check:
259	! Issued at start of each new field period.
260	if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
261	call mpif_gf_request
262	end if
263
264	! Recveiving process(es):
265	! eso TODO: at this point we do not know if clwc/ciwc will be available
266	! at next time step. Issue receive request anyway, cancel at mpif_gf_request
267	if (memstat.gt.0.and.lmp_use_reader.and..not.lmpreader) then
268	if (mp_dev_mode) write(,) 'Receiving process: calling mpif_gf_send_vars_async. PID: ', mp_pid
269	call mpif_gf_recv_vars_async(memstat)
270	if (numbnests>0) call mpif_gf_recv_vars_nest_async(memstat)
271	end if
272
273	end if
274
275	if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',1)
276
277	if (lmpreader.and.nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
278
279	! Reader process goes back to top of time loop (unless simulation end)
280	!*********************************************************************
281
282	if (lmpreader.and.lmp_use_reader) then
283	if (itime.lt.ideltas*ldirect) then
284	cycle
285	else
286	goto 999
287	end if
288	end if
289
290
291	! Get hourly OH fields if not available
292	!****************************************************
293	if (OHREA) then
294	if (verbosity.gt.0) then
295	write (,) 'timemanager> call gethourlyOH'
296	endif
297	call gethourlyOH(itime)
298	endif
299
300
301	! Release particles
302	!******************
303
304	if (verbosity.gt.0.and.lroot) then
305	write (,) 'timemanager> Release particles'
306	endif
307
308	if (mdomainfill.ge.1) then
309	if (itime.eq.0) then
310	if (verbosity.gt.0) then
311	write (,) 'timemanager> call init_domainfill'
312	endif
313	call init_domainfill
314	else
315	if (verbosity.gt.0.and.lroot) then
316	write (,) 'timemanager> call boundcond_domainfill'
317	endif
318	call boundcond_domainfill(itime,loutend)
319	endif
320	else
321	if (verbosity.gt.0.and.lroot) then
322	print*,'call releaseparticles'
323	endif
324	call releaseparticles(itime)
325	endif
326
327
328	! Compute convective mixing for forward runs
329	! for backward runs it is done before next windfield is read in
330	!**************************************************************
331
332	if ((ldirect.eq.1).and.(lconvection.eq.1)) then
333	if (verbosity.gt.0) then
334	write (,) 'timemanager> call convmix -- forward'
335	endif
336	call convmix(itime)
337	endif
338
339	! If middle of averaging period of output fields is reached, accumulated
340	! deposited mass radioactively decays
341	!***********************************************************************
342
343	if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
344	do ks=1,nspec
345	do kp=1,maxpointspec_act
346	if (decay(ks).gt.0.) then ! TODO move this statement up 2 levels
347	do nage=1,nageclass
348	do l=1,nclassunc
349	! Mother output grid
350	do jy=0,numygrid-1
351	do ix=0,numxgrid-1
352	wetgridunc(ix,jy,ks,kp,l,nage)= &
353	wetgridunc(ix,jy,ks,kp,l,nage)* &
354	exp(-1.outstepdecay(ks))
355	drygridunc(ix,jy,ks,kp,l,nage)= &
356	drygridunc(ix,jy,ks,kp,l,nage)* &
357	exp(-1.outstepdecay(ks))
358	end do
359	end do
360	! Nested output grid
361	if (nested_output.eq.1) then
362	do jy=0,numygridn-1
363	do ix=0,numxgridn-1
364	wetgriduncn(ix,jy,ks,kp,l,nage)= &
365	wetgriduncn(ix,jy,ks,kp,l,nage)* &
366	exp(-1.outstepdecay(ks))
367	drygriduncn(ix,jy,ks,kp,l,nage)= &
368	drygriduncn(ix,jy,ks,kp,l,nage)* &
369	exp(-1.outstepdecay(ks))
370	end do
371	end do
372	endif
373	end do
374	end do
375	endif
376	end do
377	end do
378	endif
379
380
381	!!! CHANGE: These lines may be switched on to check the conservation
382	!!! of mass within FLEXPART
383	! if (itime.eq.loutnext) then
384	! do 247 ksp=1, nspec
385	! do 247 kp=1, maxpointspec_act
386	!47 xm(ksp,kp)=0.
387
388	! do 249 ksp=1, nspec
389	! do 249 j=1,numpart
390	! if (ioutputforeachrelease.eq.1) then
391	! kp=npoint(j)
392	! else
393	! kp=1
394	! endif
395	! if (itra1(j).eq.itime) then
396	! xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
397	! write(,) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
398	! endif
399	!49 continue
400	! do 248 ksp=1,nspec
401	! do 248 kp=1,maxpointspec_act
402	! xm_depw(ksp,kp)=0.
403	! xm_depd(ksp,kp)=0.
404	! do 248 nage=1,nageclass
405	! do 248 ix=0,numxgrid-1
406	! do 248 jy=0,numygrid-1
407	! do 248 l=1,nclassunc
408	! xm_depw(ksp,kp)=xm_depw(ksp,kp)
409	! + +wetgridunc(ix,jy,ksp,kp,l,nage)
410	!48 xm_depd(ksp,kp)=xm_depd(ksp,kp)
411	! + +drygridunc(ix,jy,ksp,kp,l,nage)
412	! do 246 ksp=1,nspec
413	!46 write(88,'(2i10,3e12.3)')
414	! + itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
415	! + (xm_depw(ksp,kp),kp=1,maxpointspec_act),
416	! + (xm_depd(ksp,kp),kp=1,maxpointspec_act)
417	! endif
418	!!! CHANGE
419
420
421
422
423	! Check whether concentrations are to be calculated
424	!**************************************************
425
426	if ((ldirectitime.ge.ldirectloutstart).and. &
427	(ldirectitime.le.ldirectloutend)) then ! add to grid
428	if (mod(itime-loutstart,loutsample).eq.0) then
429
430	! If we are exactly at the start or end of the concentration averaging interval,
431	! give only half the weight to this sample
432	!*****************************************************************************
433
434	if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
435	weight=0.5
436	else
437	weight=1.0
438	endif
439	outnum=outnum+weight
440
441	call conccalc(itime,weight)
442
443	endif
444
445	! :TODO: MPI output of particle positions; each process sequentially
446	! access the same file
447	if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
448	call partoutput_short(itime) ! dump particle positions in extremely compressed format
449
450
451	! Output and reinitialization of grid
452	! If necessary, first sample of new grid is also taken
453	!*****************************************************
454
455	if ((itime.eq.loutend).and.(outnum.gt.0.)) then
456	if ((iout.le.3.).or.(iout.eq.5)) then
457
458	! MPI: Root process collects/sums grids
459	!**************************************
460	call mpif_tm_reduce_grid
461
462	if (mp_measure_time) call mpif_mtime('iotime',0)
463	if (surf_only.ne.1) then
464	if (lroot) then
465	if (lnetcdfout.eq.1) then
466	call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
467	&drygridtotalunc)
468	else
469	call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
470	endif
471	else
472	! zero arrays on non-root processes
473	gridunc(:,:,:,:,:,:,:)=0.
474	creceptor(:,:)=0.
475	end if
476	else
477	if (lroot) then
478	if (lnetcdfout.eq.1) then
479	call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
480	&drygridtotalunc)
481	else
482	call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
483	end if
484	else
485	! zero arrays on non-root processes
486	gridunc(:,:,:,:,:,:,:)=0.
487	creceptor(:,:)=0.
488	endif
489	endif
490	if (mp_measure_time) call mpif_mtime('iotime',1)
491
492	if (nested_output.eq.1) then
493
494	! MPI: Root process collects/sums nested grids
495	!*********************************************
496	call mpif_tm_reduce_grid_nest
497
498	if (mp_measure_time) call mpif_mtime('iotime',0)
499
500	if (lnetcdfout.eq.0) then
501	if (surf_only.ne.1) then
502
503	if (lroot) then
504	call concoutput_nest(itime,outnum)
505	else
506	griduncn(:,:,:,:,:,:,:)=0.
507	end if
508
509	else ! :TODO: check for zeroing in the netcdf module
510	call concoutput_surf_nest(itime,outnum)
511
512	end if
513
514	else
515
516	if (surf_only.ne.1) then
517	if (lroot) then
518	call concoutput_nest_netcdf(itime,outnum)
519	else
520	griduncn(:,:,:,:,:,:,:)=0.
521	end if
522	else
523	if (lroot) then
524	call concoutput_surf_nest_netcdf(itime,outnum)
525	else
526	griduncn(:,:,:,:,:,:,:)=0.
527	end if
528	endif
529
530
531	end if
532	end if
533
534
535	outnum=0.
536	endif
537	if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
538	if (iflux.eq.1) call fluxoutput(itime)
539	if (mp_measure_time) call mpif_mtime('iotime',1)
540
541
542	! Decide whether to write an estimate of the number of particles released,
543	! or exact number (require MPI reduce operation)
544	if (mp_dev_mode) then
545	numpart_tot_mpi = numpart
546	else
547	numpart_tot_mpi = numpart*mp_partgroup_np
548	end if
549
550	if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader).and.&
551	&.not.mp_dev_mode) then
552	call MPI_Reduce(numpart, numpart_tot_mpi, 1, MPI_INTEGER, MPI_SUM, id_root, &
553	& mp_comm_used, mp_ierr)
554	endif
555
556	!CGZ-lifetime: output species lifetime
557	if (lroot.or.mp_dev_mode) then
558	! write(,) 'Overview species lifetime in days', &
559	! real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
560	! write(,) 'all info:',species_lifetime
561	write(*,45) itime,numpart_tot_mpi,gridtotalunc,&
562	&wetgridtotalunc,drygridtotalunc
563	! if (verbosity.gt.0) then
564	! write (,) 'timemanager> starting simulation'
565	! end if
566	end if
567
568	45 format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES: Uncertainty: ',3f7.3)
569	46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
570	if (ipout.ge.1) then
571	do ip=0, mp_partgroup_np-1
572	if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
573	call mpif_mpi_barrier
574	end do
575	end if
576
577	loutnext=loutnext+loutstep
578	loutstart=loutnext-loutaver/2
579	loutend=loutnext+loutaver/2
580	if (itime.eq.loutstart) then
581	weight=0.5
582	outnum=outnum+weight
583	call conccalc(itime,weight)
584	endif
585
586
587	! Check, whether particles are to be split:
588	! If so, create new particles and attribute all information from the old
589	! particles also to the new ones; old and new particles both get half the
590	! mass of the old ones
591	!************************************************************************
592
593	if (ldirectitime.ge.ldirectitsplit) then
594	n=numpart
595	do j=1,numpart
596	if (ldirectitime.ge.ldirectitrasplit(j)) then
597	! if (n.lt.maxpart) then
598	if (n.lt.maxpart_mpi) then
599	n=n+1
600	itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
601	itrasplit(n)=itrasplit(j)
602	itramem(n)=itramem(j)
603	itra1(n)=itra1(j)
604	idt(n)=idt(j)
605	npoint(n)=npoint(j)
606	nclass(n)=nclass(j)
607	xtra1(n)=xtra1(j)
608	ytra1(n)=ytra1(j)
609	ztra1(n)=ztra1(j)
610	uap(n)=uap(j)
611	ucp(n)=ucp(j)
612	uzp(n)=uzp(j)
613	us(n)=us(j)
614	vs(n)=vs(j)
615	ws(n)=ws(j)
616	cbt(n)=cbt(j)
617	do ks=1,nspec
618	xmass1(j,ks)=xmass1(j,ks)/2.
619	xmass1(n,ks)=xmass1(j,ks)
620	end do
621	endif
622	endif
623	end do
624	numpart=n
625	endif
626	endif
627	endif
628
629
630	if (itime.eq.ideltas) exit ! almost finished
631
632	! Compute interval since radioactive decay of deposited mass was computed
633	!************************************************************************
634
635	if (itime.lt.loutnext) then
636	ldeltat=itime-(loutnext-loutstep)
637	else ! first half of next interval
638	ldeltat=itime-loutnext
639	endif
640
641
642	! Loop over all particles
643	!************************
644	if (mp_measure_time) call mpif_mtime('partloop1',0)
645
646	!--------------------------------------------------------------------------------
647	! various variables for testing reason of CBL scheme, by mc
648	well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
649	well_mixed_norm=0. !erase normalization to test well mixed condition: modified by mc
650	avg_ol=0.
651	avg_wst=0.
652	avg_h=0.
653	avg_air_dens=0. !erase vector to obtain air density at particle positions: modified by mc
654	!--------------------------------------------------------------------------------
655
656	do j=1,numpart
657
658	! If integration step is due, do it
659	!**********************************
660
661	if (itra1(j).eq.itime) then
662
663	if (ioutputforeachrelease.eq.1) then
664	kp=npoint(j)
665	else
666	kp=1
667	endif
668	! Determine age class of the particle
669	itage=abs(itra1(j)-itramem(j))
670	do nage=1,nageclass
671	if (itage.lt.lage(nage)) exit
672	end do
673
674	! Initialize newly released particle
675	!***********************************
676
677	if ((itramem(j).eq.itime).or.(itime.eq.0)) &
678	call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
679	us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
680
681	! Memorize particle positions
682	!****************************
683
684	xold=xtra1(j)
685	yold=ytra1(j)
686	zold=ztra1(j)
687
688	! Integrate Lagevin equation for lsynctime seconds
689	!*************************************************
690
691	if (mp_measure_time) call mpif_mtime('advance',0)
692
693	call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
694	us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
695	cbt(j))
696
697	if (mp_measure_time) call mpif_mtime('advance',1)
698
699
700	! Calculate the gross fluxes across layer interfaces
701	!***************************************************
702
703	if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
704
705
706	! Determine, when next time step is due
707	! If trajectory is terminated, mark it
708	!**************************************
709
710	if (nstop.gt.1) then
711	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
712	itra1(j)=-999999999
713	else
714	itra1(j)=itime+lsynctime
715
716
717	! Dry deposition and radioactive decay for each species
718	! Also check maximum (of all species) of initial mass remaining on the particle;
719	! if it is below a threshold value, terminate particle
720	!*****************************************************************************
721
722	xmassfract=0.
723	do ks=1,nspec
724	if (decay(ks).gt.0.) then ! radioactive decay
725	decfact=exp(-real(abs(lsynctime))*decay(ks))
726	else
727	decfact=1.
728	endif
729
730	if (DRYDEPSPEC(ks)) then ! dry deposition
731	drydeposit(ks)=xmass1(j,ks)prob(ks)decfact
732	xmass1(j,ks)=xmass1(j,ks)(1.-prob(ks))decfact
733	if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
734	drydeposit(ks)=drydeposit(ks)* &
735	exp(real(abs(ldeltat))*decay(ks))
736	endif
737	else ! no dry deposition
738	xmass1(j,ks)=xmass1(j,ks)*decfact
739	endif
740
741	! Skip check on mass fraction when npoint represents particle number
742	if (mdomainfill.eq.0.and.mquasilag.eq.0) then
743	if (xmass(npoint(j),ks).gt.0.)then
744	xmassfract=max(xmassfract,real(npart(npoint(j)))* &
745	xmass1(j,ks)/xmass(npoint(j),ks))
746
747	!CGZ-lifetime: Check mass fraction left/save lifetime
748	! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
749	!Mass below 1% of initial >register lifetime
750	! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
751
752	! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
753	! species_lifetime(ks,2)= species_lifetime(ks,2)+1
754	! endif
755	!CGZ-lifetime: Check mass fraction left/save lifetime
756
757	endif
758	else
759	xmassfract=1.
760	endif
761	end do
762
763	if (xmassfract.lt.minmass) then ! .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then ! terminate all particles carrying less mass
764	! print,'terminated particle ',j,' for small mass (', sum(real(npart(npoint(j))) &
765	! xmass1(j,:)), ' of ', sum(xmass(npoint(j),:)),')'
766	itra1(j)=-999999999
767	endif
768
769	! Sabine Eckhardt, June 2008
770	! don't create depofield for backward runs
771	if (DRYDEP.AND.(ldirect.eq.1)) then
772	call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
773	real(ytra1(j)),nage,kp)
774	if (nested_output.eq.1) call drydepokernel_nest( &
775	nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
776	nage,kp)
777	endif
778
779	! Terminate trajectories that are older than maximum allowed age
780	!***************************************************************
781
782	if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
783	if (linit_cond.ge.1) &
784	call initial_cond_calc(itime+lsynctime,j)
785	itra1(j)=-999999999
786	!print*, 'terminated particle ',j,'for age'
787	endif
788	endif
789
790	endif
791
792	end do ! j=1, numpart
793
794	if(mp_measure_time) call mpif_mtime('partloop1',1)
795
796
797	! Added by mc: counter of "unstable" particle velocity during a time scale
798	! of maximumtl=20 minutes (defined in com_mod)
799
800	total_nan_intl=0
801	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)
802	sum_nan_count(i_nan)=nan_count
803	if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
804	do ii_nan=1, (maxtl/lsynctime)
805	total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
806	end do
807
808	! Output to keep track of the numerical instabilities in CBL simulation
809	! and if they are compromising the final result (or not):
810	if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl
811
812
813	end do ! itime=0,ideltas,lsynctime
814
815
816	! Complete the calculation of initial conditions for particles not yet terminated
817	!*****************************************************************************
818
819	! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
820	! the tools to do this are already in mpi_mod.f90
821	if (lroot) then
822	do j=1,numpart
823	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
824	end do
825	end if
826
827
828	if (ipout.eq.2) then
829	! MPI process 0 creates the file, the other processes append to it
830	do ip=0, mp_partgroup_np-1
831	if (ip.eq.mp_partid) then
832	!if (mp_dbg_mode) write(,) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
833	call partoutput(itime) ! dump particle positions
834	end if
835	call mpif_mpi_barrier
836	end do
837	end if
838
839	! eso :TODO: MPI
840	if (linit_cond.ge.1.and.lroot) call initial_cond_output(itime) ! dump initial cond. field
841
842
843	! De-allocate memory and end
844	!***************************
845
846	999 if (iflux.eq.1) then
847	deallocate(flux)
848	endif
849	if (OHREA) then
850	deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
851	endif
852	if (ldirect.gt.0) then
853	deallocate(drygridunc,wetgridunc)
854	endif
855	deallocate(gridunc)
856	deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass, checklifetime)
857	deallocate(ireleasestart,ireleaseend,npart,kindz)
858	deallocate(xmasssave)
859	if (nested_output.eq.1) then
860	deallocate(orooutn, arean, volumen)
861	if (ldirect.gt.0) then
862	deallocate(griduncn,drygriduncn,wetgriduncn)
863	endif
864	endif
865	deallocate(outheight,outheighthalf)
866	deallocate(oroout, area, volume)
867
868	if (mp_measure_time) call mpif_mtime('timemanager',1)
869
870	end subroutine timemanager

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

Download in other formats:

Original Format