Context Navigation

conccalc.f90 @ 420423c

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

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

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

Download in other formats:

Original Format