Context Navigation

timemanager_mpi.f90 @ ffbe224

dev

Last change on this file since ffbe224 was ffbe224, checked in by Sabine <sabine.eckhardt@…>, 4 years ago
first version for fluxoutput in netcdf
Property mode set to `100644`
File size: 34.0 KB

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

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

Download in other formats:

Original Format