Context Navigation

timemanager_mpi.f90

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

Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago
add SPDX-License-Identifier to all .f90 files
Property mode set to `100644`
File size: 34.0 KB

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

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

Download in other formats:

Original Format