Context Navigation

timemanager.f90 @ 0a94e13

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

Last change on this file since 0a94e13 was 0a94e13, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago
Added ipout=3 option for time averaged particle output
Property mode set to `100644`
File size: 29.8 KB

Rev	Line
[e200b7a]	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
[6ecb30a]	22	subroutine timemanager(metdata_format)
[e200b7a]	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 *
[8a65cb0]	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	! For compatibility with MPI version, *
	51	! variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod *
[6ecb30a]	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 *
[e200b7a]	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	! outnum number of samples for each concentration calculation *
	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 *
[6ecb30a]	88	! metdata_format format of metdata (ecmwf/gfs) *
[e200b7a]	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
[a9cf4b1]	103	#ifdef USE_NCF
[8a65cb0]	104	use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
	105	&concoutput_surf_netcdf,concoutput_surf_nest_netcdf
[a9cf4b1]	106	#endif
[e200b7a]	107
	108	implicit none
	109
[6ecb30a]	110	integer :: metdata_format
[05cf28d]	111	integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
[e200b7a]	112	! integer :: ksp
	113	integer :: loutnext,loutstart,loutend
[c9cf570]	114	integer :: ix,jy,ldeltat,itage,nage,idummy
[8a65cb0]	115	integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
[a76d954]	116	real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact,wetscav
[8a65cb0]	117	! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
	118	! real :: us(maxpart),vs(maxpart),ws(maxpart)
	119	! integer(kind=2) :: cbt(maxpart)
[6a678e3]	120	real(sp) :: gridtotalunc
	121	real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
	122	real :: xold,yold,zold,xmassfract
[c9cf570]	123	real :: grfraction(3)
[fdc0f03]	124	real, parameter :: e_inv = 1.0/exp(1.0)
[657fa36]	125
[e200b7a]	126	!double precision xm(maxspec,maxpointspec_act),
	127	! + xm_depw(maxspec,maxpointspec_act),
	128	! + xm_depd(maxspec,maxpointspec_act)
	129
	130
	131	!open(88,file='TEST.dat')
	132
	133	! First output for time 0
	134	!************************
	135
	136	loutnext=loutstep/2
	137	outnum=0.
	138	loutstart=loutnext-loutaver/2
	139	loutend=loutnext+loutaver/2
	140
[d56bfae]	141	! open(127,file=path(2)(1:length(2))//'depostat.dat'
	142	! + ,form='unformatted')
	143	!write (,) 'writing deposition statistics depostat.dat!'
[e200b7a]	144
	145	!**********************************************************************
	146	! Loop over the whole modelling period in time steps of mintime seconds
	147	!**********************************************************************
	148
[d6a0977]	149	!ZHG 2015
	150	!CGZ-lifetime: set lifetime to 0
[f75967d]	151	! checklifetime(:,:)=0
	152	! species_lifetime(:,:)=0
	153	! print*, 'Initialized lifetime'
[d6a0977]	154	!CGZ-lifetime: set lifetime to 0
	155
[fe32dca]	156	if (.not.lusekerneloutput) write(,) 'Not using the kernel'
[d1a8707]	157	if (turboff) write(,) 'Turbulence switched off'
	158
[ec7fc72]	159	write(*,46) float(itime)/3600,itime,numpart
[f13406c]	160
	161	if (verbosity.gt.0) then
	162	write (,) 'timemanager> starting simulation'
	163	if (verbosity.gt.1) then
	164	CALL SYSTEM_CLOCK(count_clock)
	165	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
	166	endif
	167	endif
	168
[e200b7a]	169	do itime=0,ideltas,lsynctime
	170
	171	! Computation of wet deposition, OH reaction and mass transfer
	172	! between two species every lsynctime seconds
	173	! maybe wet depo frequency can be relaxed later but better be on safe side
	174	! wetdepo must be called BEFORE new fields are read in but should not
	175	! be called in the very beginning before any fields are loaded, or
	176	! before particles are in the system
	177	! Code may not be correct for decay of deposition
	178	! changed by Petra Seibert 9/02
	179	!********************************************************************
	180
[f13406c]	181	if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
	182	if (verbosity.gt.0) then
	183	write (,) 'timemanager> call wetdepo'
	184	endif
[0581cac]	185	call wetdepo(itime,lsynctime,loutnext)
[f13406c]	186	endif
[e200b7a]	187
	188	if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
	189	call ohreaction(itime,lsynctime,loutnext)
	190
	191	if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
	192	stop 'associated species not yet implemented!'
	193	! call transferspec(itime,lsynctime,loutnext)
	194	endif
	195
	196	! compute convection for backward runs
	197	!*************************************
	198
[f13406c]	199	if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
	200	if (verbosity.gt.0) then
	201	write (,) 'timemanager> call convmix -- backward'
	202	endif
[6ecb30a]	203	call convmix(itime,metdata_format)
[f13406c]	204	if (verbosity.gt.1) then
	205	!CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
	206	CALL SYSTEM_CLOCK(count_clock)
	207	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
	208	endif
	209	endif
[e200b7a]	210
	211	! Get necessary wind fields if not available
	212	!*******************************************
[f13406c]	213	if (verbosity.gt.0) then
	214	write (,) 'timemanager> call getfields'
	215	endif
[6ecb30a]	216	call getfields(itime,nstop1,metdata_format)
[f13406c]	217	if (verbosity.gt.1) then
	218	CALL SYSTEM_CLOCK(count_clock)
	219	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
	220	endif
[e200b7a]	221	if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
[f13406c]	222
[8a65cb0]	223	! Get hourly OH fields if not available
	224	!****************************************************
	225	if (OHREA) then
	226	if (verbosity.gt.0) then
	227	write (,) 'timemanager> call gethourlyOH'
	228	endif
	229	call gethourlyOH(itime)
	230	if (verbosity.gt.1) then
	231	CALL SYSTEM_CLOCK(count_clock)
	232	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
	233	endif
	234	endif
	235
[e200b7a]	236	! Release particles
	237	!******************
	238
[8a65cb0]	239	if (verbosity.gt.0) then
	240	write (,) 'timemanager> Release particles'
	241	endif
[f13406c]	242
[e200b7a]	243	if (mdomainfill.ge.1) then
	244	if (itime.eq.0) then
[f13406c]	245	if (verbosity.gt.0) then
	246	write (,) 'timemanager> call init_domainfill'
	247	endif
[e200b7a]	248	call init_domainfill
	249	else
[f13406c]	250	if (verbosity.gt.0) then
	251	write (,) 'timemanager> call boundcond_domainfill'
	252	endif
[e200b7a]	253	call boundcond_domainfill(itime,loutend)
	254	endif
	255	else
[f13406c]	256	if (verbosity.gt.0) then
[8a65cb0]	257	print*,'call releaseparticles'
[f13406c]	258	endif
[e200b7a]	259	call releaseparticles(itime)
[f13406c]	260	if (verbosity.gt.1) then
	261	CALL SYSTEM_CLOCK(count_clock)
	262	WRITE(,) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
	263	endif
[e200b7a]	264	endif
	265
	266
	267	! Compute convective mixing for forward runs
	268	! for backward runs it is done before next windfield is read in
	269	!**************************************************************
	270
[f13406c]	271	if ((ldirect.eq.1).and.(lconvection.eq.1)) then
	272	if (verbosity.gt.0) then
	273	write (,) 'timemanager> call convmix -- forward'
	274	endif
[6ecb30a]	275	call convmix(itime,metdata_format)
[f13406c]	276	endif
[e200b7a]	277
	278	! If middle of averaging period of output fields is reached, accumulated
	279	! deposited mass radioactively decays
	280	!***********************************************************************
	281
	282	if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
	283	do ks=1,nspec
	284	do kp=1,maxpointspec_act
	285	if (decay(ks).gt.0.) then
	286	do nage=1,nageclass
	287	do l=1,nclassunc
	288	! Mother output grid
	289	do jy=0,numygrid-1
	290	do ix=0,numxgrid-1
	291	wetgridunc(ix,jy,ks,kp,l,nage)= &
	292	wetgridunc(ix,jy,ks,kp,l,nage)* &
	293	exp(-1.outstepdecay(ks))
	294	drygridunc(ix,jy,ks,kp,l,nage)= &
	295	drygridunc(ix,jy,ks,kp,l,nage)* &
	296	exp(-1.outstepdecay(ks))
	297	end do
	298	end do
	299	! Nested output grid
	300	if (nested_output.eq.1) then
	301	do jy=0,numygridn-1
	302	do ix=0,numxgridn-1
	303	wetgriduncn(ix,jy,ks,kp,l,nage)= &
	304	wetgriduncn(ix,jy,ks,kp,l,nage)* &
	305	exp(-1.outstepdecay(ks))
	306	drygriduncn(ix,jy,ks,kp,l,nage)= &
	307	drygriduncn(ix,jy,ks,kp,l,nage)* &
	308	exp(-1.outstepdecay(ks))
	309	end do
	310	end do
	311	endif
	312	end do
	313	end do
	314	endif
	315	end do
	316	end do
	317	endif
	318
	319	!!! CHANGE: These lines may be switched on to check the conservation
	320	!!! of mass within FLEXPART
	321	! if (itime.eq.loutnext) then
	322	! do 247 ksp=1, nspec
[d56bfae]	323	! do 247 kp=1, maxpointspec_act
[e200b7a]	324	!47 xm(ksp,kp)=0.
	325
	326	! do 249 ksp=1, nspec
	327	! do 249 j=1,numpart
	328	! if (ioutputforeachrelease.eq.1) then
	329	! kp=npoint(j)
	330	! else
	331	! kp=1
	332	! endif
	333	! if (itra1(j).eq.itime) then
	334	! xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
	335	! write(,) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
	336	! endif
	337	!49 continue
	338	! do 248 ksp=1,nspec
	339	! do 248 kp=1,maxpointspec_act
	340	! xm_depw(ksp,kp)=0.
	341	! xm_depd(ksp,kp)=0.
	342	! do 248 nage=1,nageclass
	343	! do 248 ix=0,numxgrid-1
	344	! do 248 jy=0,numygrid-1
	345	! do 248 l=1,nclassunc
	346	! xm_depw(ksp,kp)=xm_depw(ksp,kp)
	347	! + +wetgridunc(ix,jy,ksp,kp,l,nage)
	348	!48 xm_depd(ksp,kp)=xm_depd(ksp,kp)
	349	! + +drygridunc(ix,jy,ksp,kp,l,nage)
	350	! do 246 ksp=1,nspec
	351	!46 write(88,'(2i10,3e12.3)')
	352	! + itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
	353	! + (xm_depw(ksp,kp),kp=1,maxpointspec_act),
	354	! + (xm_depd(ksp,kp),kp=1,maxpointspec_act)
	355	! endif
	356	!!! CHANGE
	357
	358
	359
	360	! Check whether concentrations are to be calculated
	361	!**************************************************
	362
	363	if ((ldirectitime.ge.ldirectloutstart).and. &
	364	(ldirectitime.le.ldirectloutend)) then ! add to grid
	365	if (mod(itime-loutstart,loutsample).eq.0) then
	366
	367	! If we are exactly at the start or end of the concentration averaging interval,
	368	! give only half the weight to this sample
	369	!*****************************************************************************
	370
	371	if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
	372	weight=0.5
	373	else
	374	weight=1.0
	375	endif
	376	outnum=outnum+weight
	377	call conccalc(itime,weight)
	378	endif
	379
	380
	381	if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
	382	call partoutput_short(itime) ! dump particle positions in extremely compressed format
	383
	384
	385	! Output and reinitialization of grid
	386	! If necessary, first sample of new grid is also taken
	387	!*****************************************************
	388
	389	if ((itime.eq.loutend).and.(outnum.gt.0.)) then
[657fa36]	390	if ((iout.le.3.).or.(iout.eq.5)) then
[f13406c]	391	if (surf_only.ne.1) then
[8a65cb0]	392	if (lnetcdfout.eq.1) then
[a9cf4b1]	393	#ifdef USE_NCF
[8a65cb0]	394	call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
[a9cf4b1]	395	#endif
[8a65cb0]	396	else
	397	call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
	398	endif
[f13406c]	399	else
[8a65cb0]	400	if (verbosity.eq.1) then
	401	print*,'call concoutput_surf '
	402	call system_clock(count_clock)
	403	write(,) 'system clock',count_clock - count_clock0
	404	endif
[a9cf4b1]	405	if (lnetcdfout.eq.1) then
	406	#ifdef USE_NCF
[8a65cb0]	407	call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
[a9cf4b1]	408	#endif
[8a65cb0]	409	else
	410	call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
	411	if (verbosity.eq.1) then
	412	print*,'called concoutput_surf '
	413	call system_clock(count_clock)
	414	write(,) 'system clock',count_clock - count_clock0
	415	endif
	416	endif
[f13406c]	417	endif
	418
[8a65cb0]	419	if (nested_output .eq. 1) then
	420	if (lnetcdfout.eq.0) then
	421	if (surf_only.ne.1) then
	422	call concoutput_nest(itime,outnum)
	423	else
	424	call concoutput_surf_nest(itime,outnum)
	425	endif
	426	else
[a9cf4b1]	427	#ifdef USE_NCF
[8a65cb0]	428	if (surf_only.ne.1) then
	429	call concoutput_nest_netcdf(itime,outnum)
	430	else
	431	call concoutput_surf_nest_netcdf(itime,outnum)
	432	endif
[a9cf4b1]	433	#endif
[8a65cb0]	434	endif
	435	endif
[e200b7a]	436	outnum=0.
	437	endif
	438	if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
	439	if (iflux.eq.1) call fluxoutput(itime)
[8a65cb0]	440	write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
[d6a0977]	441
	442	!CGZ-lifetime: output species lifetime
	443	!ZHG
[f75967d]	444	! write(,) 'Overview species lifetime in days', &
	445	! real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
	446	! write(,) 'all info:',species_lifetime
[d6a0977]	447	!ZHG
	448	!CGZ-lifetime: output species lifetime
	449
[8a65cb0]	450	!write(*,46) float(itime)/3600,itime,numpart
[c7d1052]	451	45 format(i13,' Seconds simulated: ',i13, ' Particles: Uncertainty: ',3f7.3)
[ec7fc72]	452	46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
[0a94e13]	453	if (ipout.ge.1) then
	454	if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
	455	if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions
	456	endif
[e200b7a]	457	loutnext=loutnext+loutstep
	458	loutstart=loutnext-loutaver/2
	459	loutend=loutnext+loutaver/2
	460	if (itime.eq.loutstart) then
	461	weight=0.5
	462	outnum=outnum+weight
	463	call conccalc(itime,weight)
	464	endif
	465
	466
	467	! Check, whether particles are to be split:
	468	! If so, create new particles and attribute all information from the old
	469	! particles also to the new ones; old and new particles both get half the
	470	! mass of the old ones
	471	!************************************************************************
	472
	473	if (ldirectitime.ge.ldirectitsplit) then
	474	n=numpart
	475	do j=1,numpart
	476	if (ldirectitime.ge.ldirectitrasplit(j)) then
	477	if (n.lt.maxpart) then
	478	n=n+1
	479	itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
	480	itrasplit(n)=itrasplit(j)
	481	itramem(n)=itramem(j)
	482	itra1(n)=itra1(j)
	483	idt(n)=idt(j)
	484	npoint(n)=npoint(j)
	485	nclass(n)=nclass(j)
	486	xtra1(n)=xtra1(j)
	487	ytra1(n)=ytra1(j)
	488	ztra1(n)=ztra1(j)
	489	uap(n)=uap(j)
	490	ucp(n)=ucp(j)
	491	uzp(n)=uzp(j)
	492	us(n)=us(j)
	493	vs(n)=vs(j)
	494	ws(n)=ws(j)
	495	cbt(n)=cbt(j)
	496	do ks=1,nspec
	497	xmass1(j,ks)=xmass1(j,ks)/2.
	498	xmass1(n,ks)=xmass1(j,ks)
	499	end do
	500	endif
	501	endif
	502	end do
	503	numpart=n
	504	endif
	505	endif
	506	endif
	507
	508
	509	if (itime.eq.ideltas) exit ! almost finished
	510
	511	! Compute interval since radioactive decay of deposited mass was computed
	512	!************************************************************************
	513
	514	if (itime.lt.loutnext) then
	515	ldeltat=itime-(loutnext-loutstep)
	516	else ! first half of next interval
	517	ldeltat=itime-loutnext
	518	endif
	519
	520
	521	! Loop over all particles
	522	!************************
[8a65cb0]	523	! Various variables for testing reason of CBL scheme, by mc
	524	well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
	525	well_mixed_norm=0. !erase normalization to test well mixed condition: modified by mc
	526	avg_ol=0.
	527	avg_wst=0.
	528	avg_h=0.
	529	avg_air_dens=0. !erase vector to obtain air density at particle positions: modified by mc
	530	!-----------------------------------------------------------------------------
[e200b7a]	531	do j=1,numpart
	532
	533
	534	! If integration step is due, do it
	535	!**********************************
	536
	537	if (itra1(j).eq.itime) then
	538
	539	if (ioutputforeachrelease.eq.1) then
	540	kp=npoint(j)
	541	else
	542	kp=1
	543	endif
	544	! Determine age class of the particle
	545	itage=abs(itra1(j)-itramem(j))
	546	do nage=1,nageclass
	547	if (itage.lt.lage(nage)) exit
	548	end do
	549
	550	! Initialize newly released particle
	551	!***********************************
	552
	553	if ((itramem(j).eq.itime).or.(itime.eq.0)) &
	554	call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
	555	us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
	556
	557	! Memorize particle positions
	558	!****************************
	559
	560	xold=xtra1(j)
	561	yold=ytra1(j)
	562	zold=ztra1(j)
	563
[d56bfae]	564
	565	! RECEPTOR: dry/wet depovel
	566	!****************************
	567	! Before the particle is moved
	568	! the calculation of the scavenged mass shall only be done once after release
	569	! xscav_frac1 was initialised with a negative value
	570
	571	if (DRYBKDEP) then
	572	do ks=1,nspec
	573	if ((xscav_frac1(j,ks).lt.0)) then
	574	call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
	575	if (DRYDEPSPEC(ks)) then ! dry deposition
	576	xscav_frac1(j,ks)=prob_rec(ks)
	577	else
[b85b020]	578	xmass1(j,ks)=0.
[d56bfae]	579	xscav_frac1(j,ks)=0.
	580	endif
	581	endif
	582	enddo
	583	endif
	584
	585	if (WETBKDEP) then
	586	do ks=1,nspec
	587	if ((xscav_frac1(j,ks).lt.0)) then
	588	call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
[a76d954]	589	if (wetscav.gt.0) then
	590	xscav_frac1(j,ks)=wetscav* &
[d56bfae]	591	(zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
	592	else
	593	xmass1(j,ks)=0.
	594	xscav_frac1(j,ks)=0.
	595	endif
	596	endif
	597	enddo
	598	endif
[54cbd6c]	599
[e200b7a]	600	! Integrate Lagevin equation for lsynctime seconds
	601	!*************************************************
	602
[54cbd6c]	603	if (verbosity.gt.0) then
	604	if (j.eq.1) then
[0581cac]	605	write (,) 'timemanager> call advance'
	606	endif
	607	endif
	608
	609	call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
[5deb48c]	610	us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
	611	cbt(j))
[d56bfae]	612	! write (,) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
[e200b7a]	613
[0a94e13]	614	! Calculate average position for particle dump output
	615	!****************************************************
	616
	617	if (ipout.eq.3) call partpos_average(itime,j)
	618
	619
[e200b7a]	620	! Calculate the gross fluxes across layer interfaces
	621	!***************************************************
	622
	623	if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
	624
	625
	626	! Determine, when next time step is due
	627	! If trajectory is terminated, mark it
	628	!**************************************
	629
	630	if (nstop.gt.1) then
	631	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
	632	itra1(j)=-999999999
	633	else
	634	itra1(j)=itime+lsynctime
	635
	636
	637	! Dry deposition and radioactive decay for each species
	638	! Also check maximum (of all species) of initial mass remaining on the particle;
	639	! if it is below a threshold value, terminate particle
	640	!*****************************************************************************
	641
	642	xmassfract=0.
	643	do ks=1,nspec
	644	if (decay(ks).gt.0.) then ! radioactive decay
	645	decfact=exp(-real(abs(lsynctime))*decay(ks))
	646	else
	647	decfact=1.
	648	endif
	649
	650	if (DRYDEPSPEC(ks)) then ! dry deposition
	651	drydeposit(ks)=xmass1(j,ks)prob(ks)decfact
	652	xmass1(j,ks)=xmass1(j,ks)(1.-prob(ks))decfact
	653	if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
	654	drydeposit(ks)=drydeposit(ks)* &
	655	exp(real(abs(ldeltat))*decay(ks))
	656	endif
	657	else ! no dry deposition
	658	xmass1(j,ks)=xmass1(j,ks)*decfact
	659	endif
	660
[18adf60]	661	! Skip check on mass fraction when npoint represents particle number
	662	if (mdomainfill.eq.0.and.mquasilag.eq.0) then
[e200b7a]	663	if (xmass(npoint(j),ks).gt.0.) &
	664	xmassfract=max(xmassfract,real(npart(npoint(j)))* &
	665	xmass1(j,ks)/xmass(npoint(j),ks))
[d6a0977]	666	!ZHG 2015
	667	!CGZ-lifetime: Check mass fraction left/save lifetime
[6a678e3]	668	! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
[d6a0977]	669	!Mass below 1% of initial >register lifetime
[f75967d]	670	! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
	671	! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
	672	! species_lifetime(ks,2)= species_lifetime(ks,2)+1
	673	! endif
[d6a0977]	674	!CGZ-lifetime: Check mass fraction left/save lifetime
	675	!ZHG 2015
[e200b7a]	676	else
[18adf60]	677	xmassfract=1.0
	678	end if
[e200b7a]	679	end do
	680
[41d8574]	681	if (xmassfract.lt.minmass) then ! terminate all particles carrying less mass
[e200b7a]	682	itra1(j)=-999999999
[8a65cb0]	683	if (verbosity.gt.0) then
	684	print*,'terminated particle ',j,' for small mass'
	685	endif
[e200b7a]	686	endif
	687
	688	! Sabine Eckhardt, June 2008
	689	! don't create depofield for backward runs
	690	if (DRYDEP.AND.(ldirect.eq.1)) then
	691	call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
	692	real(ytra1(j)),nage,kp)
	693	if (nested_output.eq.1) call drydepokernel_nest( &
	694	nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
	695	nage,kp)
	696	endif
	697
	698	! Terminate trajectories that are older than maximum allowed age
	699	!***************************************************************
	700
	701	if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
[8a65cb0]	702	if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
[e200b7a]	703	itra1(j)=-999999999
[8a65cb0]	704	if (verbosity.gt.0) then
	705	print*,'terminated particle ',j,' for age'
	706	endif
[e200b7a]	707	endif
	708	endif
	709
	710	endif
	711
[8a65cb0]	712	end do !loop over particles
	713
	714	! Counter of "unstable" particle velocity during a time scale of
	715	! maximumtl=20 minutes (defined in com_mod)
	716	!***************************************************************
	717
[41d8574]	718	total_nan_intl=0
	719	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)
	720	sum_nan_count(i_nan)=nan_count
	721	if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
	722	do ii_nan=1, (maxtl/lsynctime)
	723	total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
	724	end do
[8a65cb0]	725	! Output to keep track of the numerical instabilities in CBL simulation and if
	726	! they are compromising the final result (or not)
[41d8574]	727	if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl
[8a65cb0]	728
[e200b7a]	729	end do
	730
	731
	732	! Complete the calculation of initial conditions for particles not yet terminated
	733	!*****************************************************************************
	734
	735	do j=1,numpart
	736	if (linit_cond.ge.1) call initial_cond_calc(itime,j)
	737	end do
	738
	739	if (ipout.eq.2) call partoutput(itime) ! dump particle positions
	740
	741	if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
	742
[78e62dc]	743	!close(104)
[e200b7a]	744
	745	! De-allocate memory and end
	746	!***************************
	747
	748	if (iflux.eq.1) then
	749	deallocate(flux)
	750	endif
[8a65cb0]	751	if (OHREA) then
	752	deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
[e200b7a]	753	endif
	754	if (ldirect.gt.0) then
	755	deallocate(drygridunc,wetgridunc)
	756	endif
	757	deallocate(gridunc)
	758	deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
	759	deallocate(ireleasestart,ireleaseend,npart,kindz)
	760	deallocate(xmasssave)
	761	if (nested_output.eq.1) then
	762	deallocate(orooutn, arean, volumen)
	763	if (ldirect.gt.0) then
	764	deallocate(griduncn,drygriduncn,wetgriduncn)
	765	endif
	766	endif
	767	deallocate(outheight,outheighthalf)
	768	deallocate(oroout, area, volume)
	769
	770	end subroutine timemanager
	771

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

Download in other formats:

Original Format