Context Navigation

timemanager_mpi.f90 @ b1f28c3

dev

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

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

Download in other formats:

Original Format