Context Navigation

timemanager_mpi.f90 @ f3054ea

GFS_025dev

Last change on this file since f3054ea was f3054ea, checked in by Espen Sollum <eso@…>, 4 years ago
Changed from grib_api to eccodes. MPI: implemented linversionout=1; fixed calculation of grid_initial fields.
Property mode set to `100644`
File size: 36.1 KB

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

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

Download in other formats:

Original Format