Context Navigation

timemanager.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: 29.2 KB

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

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

Download in other formats:

Original Format