Context Navigation

conccalc.f90 @ 47f96e5

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

Last change on this file since 47f96e5 was 08a38b5, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago
Added (hardcoded) option to output number of particles per grid cell.
Property mode set to `100644`
File size: 19.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 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	! integer xscav_count
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	!***************************************************************************
67	! xscav_count=0
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
78	! if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
79
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
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
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
135	rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,memind(2)) &
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
189	if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
190	(xl.gt.real(numxgrid-1)-0.5).or. &
191	(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
192	if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
193	(jy.le.numygrid-1)) then
194	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
198	xmass1(i,ks)/rhoiweightmax(xscav_frac1(i,ks),0.0)
199	end do
200	else
201	if (lparticlecountoutput) then
202	do ks=1,nspec
203	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
204	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
205	end do
206	else
207	do ks=1,nspec
208	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
209	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
210	xmass1(i,ks)/rhoi*weight
211	end do
212	end if
213	endif
214	endif
215
216	else ! attribution via uniform kernel
217
218	ddx=xl-real(ix) ! distance to left cell border
219	ddy=yl-real(jy) ! distance to lower cell border
220	if (ddx.gt.0.5) then
221	ixp=ix+1
222	wx=1.5-ddx
223	else
224	ixp=ix-1
225	wx=0.5+ddx
226	endif
227
228	if (ddy.gt.0.5) then
229	jyp=jy+1
230	wy=1.5-ddy
231	else
232	jyp=jy-1
233	wy=0.5+ddy
234	endif
235
236
237	! Determine mass fractions for four grid points
238	!**********************************************
239
240	if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
241	if ((jy.ge.0).and.(jy.le.numygrid-1)) then
242	w=wx*wy
243	if (DRYBKDEP.or.WETBKDEP) then
244	do ks=1,nspec
245	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
246	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
247	xmass1(i,ks)/rhoiwweight*max(xscav_frac1(i,ks),0.0)
248	end do
249	else
250	do ks=1,nspec
251	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
252	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
253	xmass1(i,ks)/rhoiweightw
254	end do
255	endif
256	endif
257
258	if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
259	w=wx*(1.-wy)
260	if (DRYBKDEP.or.WETBKDEP) then
261	do ks=1,nspec
262	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
263	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
264	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
265	end do
266	else
267	do ks=1,nspec
268	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
269	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
270	xmass1(i,ks)/rhoiweightw
271	end do
272	endif
273	endif
274	endif !ix ge 0
275
276
277	if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
278	if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
279	w=(1.-wx)*(1.-wy)
280	if (DRYBKDEP.or.WETBKDEP) then
281	do ks=1,nspec
282	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
283	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
284	xmass1(i,ks)/rhoiwweight*max(xscav_frac1(i,ks),0.0)
285	end do
286	else
287	do ks=1,nspec
288	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
289	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
290	xmass1(i,ks)/rhoiweightw
291	end do
292	endif
293	endif
294
295	if ((jy.ge.0).and.(jy.le.numygrid-1)) then
296	w=(1.-wx)*wy
297	if (DRYBKDEP.or.WETBKDEP) then
298	do ks=1,nspec
299	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
300	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
301	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
302	end do
303	else
304	do ks=1,nspec
305	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
306	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
307	xmass1(i,ks)/rhoiweightw
308	end do
309	endif
310	endif
311	endif !ixp ge 0
312	endif
313
314	!************************************
315	! Do everything for the nested domain
316	!************************************
317
318	if (nested_output.eq.1) then
319	xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
320	yl=(ytra1(i)*dy+youtshiftn)/dyoutn
321	ix=int(xl)
322	if (xl.lt.0.) ix=ix-1
323	jy=int(yl)
324	if (yl.lt.0.) jy=jy-1
325
326
327	! For particles aged less than 3 hours, attribute particle mass to grid cell
328	! it resides in rather than use the kernel, in order to avoid its smoothing effect.
329	! For older particles, use the uniform kernel.
330	! If a particle is close to the domain boundary, do not use the kernel either.
331	!*****************************************************************************
332
333	if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
334	(xl.gt.real(numxgridn-1)-0.5).or. &
335	(yl.gt.real(numygridn-1)-0.5).or.(lnokernel)) then ! no kernel, direct attribution to grid cell
336	if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
337	(jy.le.numygridn-1)) then
338	if (DRYBKDEP.or.WETBKDEP) then
339	do ks=1,nspec
340	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
341	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
342	xmass1(i,ks)/rhoiweightmax(xscav_frac1(i,ks),0.0)
343	end do
344	else
345	if (lparticlecountoutput) then
346	do ks=1,nspec
347	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
348	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
349	end do
350	else
351	do ks=1,nspec
352	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
353	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
354	xmass1(i,ks)/rhoi*weight
355	end do
356	endif
357	endif
358	endif
359
360	else ! attribution via uniform kernel
361
362	ddx=xl-real(ix) ! distance to left cell border
363	ddy=yl-real(jy) ! distance to lower cell border
364	if (ddx.gt.0.5) then
365	ixp=ix+1
366	wx=1.5-ddx
367	else
368	ixp=ix-1
369	wx=0.5+ddx
370	endif
371
372	if (ddy.gt.0.5) then
373	jyp=jy+1
374	wy=1.5-ddy
375	else
376	jyp=jy-1
377	wy=0.5+ddy
378	endif
379
380
381	! Determine mass fractions for four grid points
382	!**********************************************
383
384	if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
385	if ((jy.ge.0).and.(jy.le.numygridn-1)) then
386	w=wx*wy
387	if (DRYBKDEP.or.WETBKDEP) then
388	do ks=1,nspec
389	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
390	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
391	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
392	end do
393	else
394	do ks=1,nspec
395	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
396	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
397	xmass1(i,ks)/rhoiweightw
398	end do
399	endif
400	endif
401
402	if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
403	w=wx*(1.-wy)
404	if (DRYBKDEP.or.WETBKDEP) then
405	do ks=1,nspec
406	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
407	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
408	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
409	end do
410	else
411	do ks=1,nspec
412	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
413	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
414	xmass1(i,ks)/rhoiweightw
415	end do
416	endif
417	endif
418	endif
419
420
421	if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
422	if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
423	w=(1.-wx)*(1.-wy)
424	if (DRYBKDEP.or.WETBKDEP) then
425	do ks=1,nspec
426	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
427	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
428	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
429	end do
430	else
431	do ks=1,nspec
432	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
433	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
434	xmass1(i,ks)/rhoiweightw
435	end do
436	endif
437	endif
438
439	if ((jy.ge.0).and.(jy.le.numygridn-1)) then
440	w=(1.-wx)*wy
441	if (DRYBKDEP.or.WETBKDEP) then
442	do ks=1,nspec
443	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
444	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
445	xmass1(i,ks)/rhoiweightw*max(xscav_frac1(i,ks),0.0)
446	end do
447	else
448	do ks=1,nspec
449	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
450	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
451	xmass1(i,ks)/rhoiweightw
452	end do
453	endif
454	endif
455	endif
456	endif
457	endif
458	endif
459	20 continue
460	end do
461	! write(,) 'xscav count:',xscav_count
462
463	!***********************************************************************
464	! 2. Evaluate concentrations at receptor points, using the kernel method
465	!***********************************************************************
466
467	do n=1,numreceptor
468
469
470	! Reset concentrations
471	!*********************
472
473	do ks=1,nspec
474	c(ks)=0.
475	end do
476
477
478	! Estimate concentration at receptor
479	!***********************************
480
481	do i=1,numpart
482
483	if (itra1(i).ne.itime) goto 40
484	itage=abs(itra1(i)-itramem(i))
485
486	hz=min(50.+0.3*sqrt(real(itage)),hzmax)
487	zd=ztra1(i)/hz
488	if (zd.gt.1.) goto 40 ! save computing time, leave loop
489
490	hx=min((0.29+2.222e-3sqrt(real(itage)))dx+ &
491	real(itage)*1.2e-5,hxmax) ! 80 km/day
492	xd=(xtra1(i)-xreceptor(n))/hx
493	if (xd*xd.gt.1.) goto 40 ! save computing time, leave loop
494
495	hy=min((0.18+1.389e-3sqrt(real(itage)))dy+ &
496	real(itage)*7.5e-6,hymax) ! 80 km/day
497	yd=(ytra1(i)-yreceptor(n))/hy
498	if (yd*yd.gt.1.) goto 40 ! save computing time, leave loop
499	h=hxhyhz
500
501	r2=xdxd+ydyd+zd*zd
502	if (r2.lt.1.) then
503	xkern=factor*(1.-r2)
504	do ks=1,nspec
505	c(ks)=c(ks)+xmass1(i,ks)*xkern/h
506	end do
507	endif
508	40 continue
509	end do
510
511	do ks=1,nspec
512	creceptor(n,ks)=creceptor(n,ks)+2.weightc(ks)/receptorarea(n)
513	end do
514	end do
515
516	end subroutine conccalc

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

Download in other formats:

Original Format