Context Navigation

conccalc.f90 @ e52967c

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

Last change on this file since e52967c was e52967c, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago
Quick fix for segmentation fault when particle is at north pole
Property mode set to `100644`
File size: 18.6 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
	22	subroutine conccalc(itime,weight)
[462f74b]	23	! i i
[e200b7a]	24	!*****************************************************************************
	25	! *
	26	! Calculation of the concentrations on a regular grid using volume *
	27	! sampling *
	28	! *
	29	! Author: A. Stohl *
	30	! *
	31	! 24 May 1996 *
	32	! *
	33	! April 2000: Update to calculate age spectra *
	34	! Bug fix to avoid negative conc. at the domain boundaries, *
	35	! as suggested by Petra Seibert *
	36	! *
	37	! 2 July 2002: re-order if-statements in order to optimize CPU time *
	38	! *
	39	! *
	40	!*****************************************************************************
	41	! *
	42	! Variables: *
	43	! nspeciesdim = nspec for forward runs, 1 for backward runs *
	44	! *
	45	!*****************************************************************************
	46
	47	use unc_mod
	48	use outg_mod
	49	use par_mod
	50	use com_mod
	51
	52	implicit none
	53
	54	integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
	55	integer :: il,ind,indz,indzp,nrelpointer
	56	real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
	57	real :: weight,hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy
	58	real :: rhoprof(2),rhoi
	59	real :: xl,yl,wx,wy,w
	60	real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
[9287c01]	61	! integer xscav_count
[e200b7a]	62
	63	! For forward simulations, make a loop over the number of species;
	64	! for backward simulations, make an additional loop over the
	65	! releasepoints
	66	!***************************************************************************
[9287c01]	67	! xscav_count=0
[e200b7a]	68	do i=1,numpart
	69	if (itra1(i).ne.itime) goto 20
	70
	71	! Determine age class of the particle
	72	itage=abs(itra1(i)-itramem(i))
	73	do nage=1,nageclass
	74	if (itage.lt.lage(nage)) goto 33
	75	end do
	76	33 continue
	77
[9287c01]	78	! if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
[33279f7]	79
[e200b7a]	80	! For special runs, interpolate the air density to the particle position
	81	!************************************************************************
	82	!***********************************************************************
	83	!AF IND_SOURCE switches between different units for concentrations at the source
	84	!Af NOTE that in backward simulations the release of particles takes place
	85	!Af at the receptor and the sampling at the source.
	86	!Af 1="mass"
	87	!Af 2="mass mixing ratio"
	88	!Af IND_RECEPTOR switches between different units for concentrations at the receptor
	89	!Af 1="mass"
	90	!Af 2="mass mixing ratio"
	91
	92	!Af switches for the conccalcfile:
	93	!AF IND_SAMP = 0 : xmass * 1
	94	!Af IND_SAMP = -1 : xmass / rho
	95
	96	!Af ind_samp is defined in readcommand.f
	97
	98	if ( ind_samp .eq. -1 ) then
	99
	100	ix=int(xtra1(i))
	101	jy=int(ytra1(i))
	102	ixp=ix+1
	103	jyp=jy+1
	104	ddx=xtra1(i)-real(ix)
	105	ddy=ytra1(i)-real(jy)
	106	rddx=1.-ddx
	107	rddy=1.-ddy
	108	p1=rddx*rddy
	109	p2=ddx*rddy
	110	p3=rddx*ddy
	111	p4=ddx*ddy
	112
[e52967c]	113	! eso: Temporary fix for particle exactly at north pole
	114	if (jyp >= nymax) then
	115	! write(,) 'WARNING: conccalc.f90 jyp >= nymax'
	116	jyp=jyp-1
	117	end if
	118
[e200b7a]	119	do il=2,nz
	120	if (height(il).gt.ztra1(i)) then
	121	indz=il-1
	122	indzp=il
	123	goto 6
	124	endif
	125	end do
	126	6 continue
	127
	128	dz1=ztra1(i)-height(indz)
	129	dz2=height(indzp)-ztra1(i)
	130	dz=1./(dz1+dz2)
	131
	132	! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
	133	!*****************************************************************************
	134	do ind=indz,indzp
[16b61a5]	135	rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,memind(2)) &
[e200b7a]	136	+p2*rho(ixp,jy ,ind,2) &
	137	+p3*rho(ix ,jyp,ind,2) &
	138	+p4*rho(ixp,jyp,ind,2)
	139	end do
	140	rhoi=(dz1rhoprof(2)+dz2rhoprof(1))*dz
	141	elseif (ind_samp.eq.0) then
	142	rhoi = 1.
	143	endif
	144
	145
	146	!****************************************************************************
	147	! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
	148	!****************************************************************************
	149
	150
	151	! For backward simulations, look from which release point the particle comes from
	152	! For domain-filling trajectory option, npoint contains a consecutive particle
	153	! number, not the release point information. Therefore, nrelpointer is set to 1
	154	! for the domain-filling option.
	155	!*****************************************************************************
	156
	157	if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
	158	nrelpointer=1
	159	else
	160	nrelpointer=npoint(i)
	161	endif
	162
	163	do kz=1,numzgrid ! determine height of cell
	164	if (outheight(kz).gt.ztra1(i)) goto 21
	165	end do
	166	21 continue
	167	if (kz.le.numzgrid) then ! inside output domain
	168
	169
	170	!********************************
	171	! Do everything for mother domain
	172	!********************************
	173
	174	xl=(xtra1(i)*dx+xoutshift)/dxout
	175	yl=(ytra1(i)*dy+youtshift)/dyout
	176	ix=int(xl)
	177	if (xl.lt.0.) ix=ix-1
	178	jy=int(yl)
	179	if (yl.lt.0.) jy=jy-1
	180
	181
	182
	183	! For particles aged less than 3 hours, attribute particle mass to grid cell
	184	! it resides in rather than use the kernel, in order to avoid its smoothing effect.
	185	! For older particles, use the uniform kernel.
	186	! If a particle is close to the domain boundary, do not use the kernel either.
	187	!*****************************************************************************
	188
[4c64400]	189	if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
[e200b7a]	190	(xl.gt.real(numxgrid-1)-0.5).or. &
[2bec33e]	191	(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
[e200b7a]	192	if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
	193	(jy.le.numygrid-1)) then
[54cbd6c]	194	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	195	do ks=1,nspec
	196	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	197	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	198	xmass1(i,ks)/rhoiweightmax(xscav_frac1(i,ks),0.0)
[462f74b]	199	end do
	200	else
	201	do ks=1,nspec
	202	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	203	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	204	xmass1(i,ks)/rhoi*weight
[462f74b]	205	end do
	206	endif
[e200b7a]	207	endif
	208
[92a74b2]	209	else ! attribution via uniform kernel
[e200b7a]	210
	211	ddx=xl-real(ix) ! distance to left cell border
	212	ddy=yl-real(jy) ! distance to lower cell border
	213	if (ddx.gt.0.5) then
	214	ixp=ix+1
	215	wx=1.5-ddx
	216	else
	217	ixp=ix-1
	218	wx=0.5+ddx
	219	endif
	220
	221	if (ddy.gt.0.5) then
	222	jyp=jy+1
	223	wy=1.5-ddy
	224	else
	225	jyp=jy-1
	226	wy=0.5+ddy
	227	endif
	228
	229
	230	! Determine mass fractions for four grid points
	231	!**********************************************
	232
	233	if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
	234	if ((jy.ge.0).and.(jy.le.numygrid-1)) then
	235	w=wx*wy
[54cbd6c]	236	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	237	do ks=1,nspec
	238	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	239	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	240	xmass1(i,ks)/rhoiwweight*max(xscav_frac1(i,ks),0.0)
[462f74b]	241	end do
	242	else
	243	do ks=1,nspec
	244	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	245	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	246	xmass1(i,ks)/rhoiweightw
[462f74b]	247	end do
	248	endif
[e200b7a]	249	endif
	250
	251	if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
	252	w=wx*(1.-wy)
[54cbd6c]	253	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	254	do ks=1,nspec
	255	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
	256	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	257	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	258	end do
	259	else
	260	do ks=1,nspec
	261	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	262	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
	263	xmass1(i,ks)/rhoiweightw
[462f74b]	264	end do
	265	endif
[e200b7a]	266	endif
[462f74b]	267	endif !ix ge 0
[e200b7a]	268
	269
	270	if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
	271	if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
	272	w=(1.-wx)*(1.-wy)
[54cbd6c]	273	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	274	do ks=1,nspec
	275	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
	276	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	277	xmass1(i,ks)/rhoiwweight*max(xscav_frac1(i,ks),0.0)
[462f74b]	278	end do
	279	else
	280	do ks=1,nspec
	281	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	282	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
	283	xmass1(i,ks)/rhoiweightw
[462f74b]	284	end do
	285	endif
[e200b7a]	286	endif
	287
	288	if ((jy.ge.0).and.(jy.le.numygrid-1)) then
	289	w=(1.-wx)*wy
[54cbd6c]	290	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	291	do ks=1,nspec
	292	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	293	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	294	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	295	end do
	296	else
	297	do ks=1,nspec
	298	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	299	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	300	xmass1(i,ks)/rhoiweightw
[462f74b]	301	end do
	302	endif
[e200b7a]	303	endif
[462f74b]	304	endif !ixp ge 0
	305	endif
[e200b7a]	306
	307	!************************************
	308	! Do everything for the nested domain
	309	!************************************
	310
	311	if (nested_output.eq.1) then
	312	xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
	313	yl=(ytra1(i)*dy+youtshiftn)/dyoutn
	314	ix=int(xl)
	315	if (xl.lt.0.) ix=ix-1
	316	jy=int(yl)
	317	if (yl.lt.0.) jy=jy-1
	318
	319
	320	! For particles aged less than 3 hours, attribute particle mass to grid cell
	321	! it resides in rather than use the kernel, in order to avoid its smoothing effect.
	322	! For older particles, use the uniform kernel.
	323	! If a particle is close to the domain boundary, do not use the kernel either.
	324	!*****************************************************************************
	325
	326	if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
	327	(xl.gt.real(numxgridn-1)-0.5).or. &
[2bec33e]	328	(yl.gt.real(numygridn-1)-0.5).or.(lnokernel)) then ! no kernel, direct attribution to grid cell
[e200b7a]	329	if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
	330	(jy.le.numygridn-1)) then
[54cbd6c]	331	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	332	do ks=1,nspec
	333	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	334	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	335	xmass1(i,ks)/rhoiweightmax(xscav_frac1(i,ks),0.0)
[462f74b]	336	end do
	337	else
	338	do ks=1,nspec
	339	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	340	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	341	xmass1(i,ks)/rhoi*weight
[462f74b]	342	end do
	343	endif
[e200b7a]	344	endif
	345
	346	else ! attribution via uniform kernel
	347
	348	ddx=xl-real(ix) ! distance to left cell border
	349	ddy=yl-real(jy) ! distance to lower cell border
	350	if (ddx.gt.0.5) then
	351	ixp=ix+1
	352	wx=1.5-ddx
	353	else
	354	ixp=ix-1
	355	wx=0.5+ddx
	356	endif
	357
	358	if (ddy.gt.0.5) then
	359	jyp=jy+1
	360	wy=1.5-ddy
	361	else
	362	jyp=jy-1
	363	wy=0.5+ddy
	364	endif
	365
	366
	367	! Determine mass fractions for four grid points
	368	!**********************************************
	369
	370	if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
	371	if ((jy.ge.0).and.(jy.le.numygridn-1)) then
	372	w=wx*wy
[54cbd6c]	373	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	374	do ks=1,nspec
	375	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	376	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[88929bf]	377	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	378	end do
	379	else
	380	do ks=1,nspec
	381	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	382	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	383	xmass1(i,ks)/rhoiweightw
[462f74b]	384	end do
	385	endif
[e200b7a]	386	endif
	387
	388	if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
	389	w=wx*(1.-wy)
[54cbd6c]	390	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	391	do ks=1,nspec
	392	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
	393	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
[6473ad3]	394	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	395	end do
	396	else
	397	do ks=1,nspec
	398	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	399	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
	400	xmass1(i,ks)/rhoiweightw
[462f74b]	401	end do
	402	endif
[e200b7a]	403	endif
	404	endif
	405
	406
	407	if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
	408	if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
	409	w=(1.-wx)*(1.-wy)
[54cbd6c]	410	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	411	do ks=1,nspec
	412	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
	413	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
[6473ad3]	414	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	415	end do
	416	else
	417	do ks=1,nspec
	418	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	419	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
	420	xmass1(i,ks)/rhoiweightw
[462f74b]	421	end do
	422	endif
[e200b7a]	423	endif
	424
	425	if ((jy.ge.0).and.(jy.le.numygridn-1)) then
	426	w=(1.-wx)*wy
[54cbd6c]	427	if (DRYBKDEP.or.WETBKDEP) then
[462f74b]	428	do ks=1,nspec
	429	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
	430	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
[6473ad3]	431	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
[462f74b]	432	end do
	433	else
	434	do ks=1,nspec
	435	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
[e200b7a]	436	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
	437	xmass1(i,ks)/rhoiweightw
[462f74b]	438	end do
	439	endif
[e200b7a]	440	endif
	441	endif
	442	endif
	443	endif
	444	endif
	445	20 continue
	446	end do
[9287c01]	447	! write(,) 'xscav count:',xscav_count
[e200b7a]	448
	449	!***********************************************************************
	450	! 2. Evaluate concentrations at receptor points, using the kernel method
	451	!***********************************************************************
	452
	453	do n=1,numreceptor
	454
	455
	456	! Reset concentrations
	457	!*********************
	458
	459	do ks=1,nspec
	460	c(ks)=0.
	461	end do
	462
	463
	464	! Estimate concentration at receptor
	465	!***********************************
	466
	467	do i=1,numpart
	468
	469	if (itra1(i).ne.itime) goto 40
	470	itage=abs(itra1(i)-itramem(i))
	471
	472	hz=min(50.+0.3*sqrt(real(itage)),hzmax)
	473	zd=ztra1(i)/hz
	474	if (zd.gt.1.) goto 40 ! save computing time, leave loop
	475
	476	hx=min((0.29+2.222e-3sqrt(real(itage)))dx+ &
	477	real(itage)*1.2e-5,hxmax) ! 80 km/day
	478	xd=(xtra1(i)-xreceptor(n))/hx
	479	if (xd*xd.gt.1.) goto 40 ! save computing time, leave loop
	480
	481	hy=min((0.18+1.389e-3sqrt(real(itage)))dy+ &
	482	real(itage)*7.5e-6,hymax) ! 80 km/day
	483	yd=(ytra1(i)-yreceptor(n))/hy
	484	if (yd*yd.gt.1.) goto 40 ! save computing time, leave loop
	485	h=hxhyhz
	486
	487	r2=xdxd+ydyd+zd*zd
	488	if (r2.lt.1.) then
	489	xkern=factor*(1.-r2)
	490	do ks=1,nspec
	491	c(ks)=c(ks)+xmass1(i,ks)*xkern/h
	492	end do
	493	endif
	494	40 continue
	495	end do
	496
	497	do ks=1,nspec
	498	creceptor(n,ks)=creceptor(n,ks)+2.weightc(ks)/receptorarea(n)
	499	end do
	500	end do
	501
	502	end subroutine conccalc

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

Download in other formats:

Original Format