Context Navigation

conccalc.f90 @ 54cbd6c

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

Last change on this file since 54cbd6c was 54cbd6c, checked in by Sabine <sabine.eckhardt@…>, 8 years ago
all f90 files for dry/wet backward mode
Property mode set to `100644`
File size: 18.6 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	integer :: usekernel
64
65	usekernel=0
66	if (usekernel.ne.1) then
67	write (,) 'NOT USING THE KERNEL!'
68	endif
69	! For forward simulations, make a loop over the number of species;
70	! for backward simulations, make an additional loop over the
71	! releasepoints
72	!***************************************************************************
73
74	do i=1,numpart
75	if (itra1(i).ne.itime) goto 20
76
77	! Determine age class of the particle
78	itage=abs(itra1(i)-itramem(i))
79	do nage=1,nageclass
80	if (itage.lt.lage(nage)) goto 33
81	end do
82	33 continue
83
84
85	! For special runs, interpolate the air density to the particle position
86	!************************************************************************
87	!***********************************************************************
88	!AF IND_SOURCE switches between different units for concentrations at the source
89	!Af NOTE that in backward simulations the release of particles takes place
90	!Af at the receptor and the sampling at the source.
91	!Af 1="mass"
92	!Af 2="mass mixing ratio"
93	!Af IND_RECEPTOR switches between different units for concentrations at the receptor
94	!Af 1="mass"
95	!Af 2="mass mixing ratio"
96
97	!Af switches for the conccalcfile:
98	!AF IND_SAMP = 0 : xmass * 1
99	!Af IND_SAMP = -1 : xmass / rho
100
101	!Af ind_samp is defined in readcommand.f
102
103	if ( ind_samp .eq. -1 ) then
104
105	ix=int(xtra1(i))
106	jy=int(ytra1(i))
107	ixp=ix+1
108	jyp=jy+1
109	ddx=xtra1(i)-real(ix)
110	ddy=ytra1(i)-real(jy)
111	rddx=1.-ddx
112	rddy=1.-ddy
113	p1=rddx*rddy
114	p2=ddx*rddy
115	p3=rddx*ddy
116	p4=ddx*ddy
117
118	do il=2,nz
119	if (height(il).gt.ztra1(i)) then
120	indz=il-1
121	indzp=il
122	goto 6
123	endif
124	end do
125	6 continue
126
127	dz1=ztra1(i)-height(indz)
128	dz2=height(indzp)-ztra1(i)
129	dz=1./(dz1+dz2)
130
131	! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
132	!*****************************************************************************
133	do ind=indz,indzp
134	rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
135	+p2*rho(ixp,jy ,ind,2) &
136	+p3*rho(ix ,jyp,ind,2) &
137	+p4*rho(ixp,jyp,ind,2)
138	end do
139	rhoi=(dz1rhoprof(2)+dz2rhoprof(1))*dz
140	elseif (ind_samp.eq.0) then
141	rhoi = 1.
142	endif
143
144
145	!****************************************************************************
146	! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
147	!****************************************************************************
148
149
150	! For backward simulations, look from which release point the particle comes from
151	! For domain-filling trajectory option, npoint contains a consecutive particle
152	! number, not the release point information. Therefore, nrelpointer is set to 1
153	! for the domain-filling option.
154	!*****************************************************************************
155
156	if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
157	nrelpointer=1
158	else
159	nrelpointer=npoint(i)
160	endif
161
162	do kz=1,numzgrid ! determine height of cell
163	if (outheight(kz).gt.ztra1(i)) goto 21
164	end do
165	21 continue
166	if (kz.le.numzgrid) then ! inside output domain
167
168
169	!********************************
170	! Do everything for mother domain
171	!********************************
172
173	xl=(xtra1(i)*dx+xoutshift)/dxout
174	yl=(ytra1(i)*dy+youtshift)/dyout
175	ix=int(xl)
176	if (xl.lt.0.) ix=ix-1
177	jy=int(yl)
178	if (yl.lt.0.) jy=jy-1
179
180	! if (i.eq.10000) write(,) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
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 (((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)).or.(usekernel.eq.0)) 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)/(1-xscav_frac1(i,ks))/rhoiweightxscav_frac1(i,ks)
199	end do
200	else
201	do ks=1,nspec
202	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
203	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
204	xmass1(i,ks)/rhoi*weight
205	end do
206	endif
207	endif
208
209	else ! attribution via uniform kernel
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
236	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
240	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiwweight*xscav_frac1(i,ks)
241	end do
242	else
243	do ks=1,nspec
244	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
245	gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
246	xmass1(i,ks)/rhoiweightw
247	end do
248	endif
249	endif
250
251	if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
252	w=wx*(1.-wy)
253	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
257	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
258	end do
259	else
260	do ks=1,nspec
261	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
262	gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
263	xmass1(i,ks)/rhoiweightw
264	end do
265	endif
266	endif
267	endif !ix ge 0
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)
273	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
277	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiwweight*xscav_frac1(i,ks)
278	end do
279	else
280	do ks=1,nspec
281	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
282	gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
283	xmass1(i,ks)/rhoiweightw
284	end do
285	endif
286	endif
287
288	if ((jy.ge.0).and.(jy.le.numygrid-1)) then
289	w=(1.-wx)*wy
290	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
294	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
295	end do
296	else
297	do ks=1,nspec
298	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
299	gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
300	xmass1(i,ks)/rhoiweightw
301	end do
302	endif
303	endif
304	endif !ixp ge 0
305	endif
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. &
328	(yl.gt.real(numygridn-1)-0.5)) then ! no kernel, direct attribution to grid cell
329	if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
330	(jy.le.numygridn-1)) then
331	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
335	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightxscav_frac1(i,ks)
336	end do
337	else
338	do ks=1,nspec
339	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
340	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
341	xmass1(i,ks)/rhoi*weight
342	end do
343	endif
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
373	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
377	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
378	end do
379	else
380	do ks=1,nspec
381	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
382	griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
383	xmass1(i,ks)/rhoiweightw
384	end do
385	endif
386	endif
387
388	if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
389	w=wx*(1.-wy)
390	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
394	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
395	end do
396	else
397	do ks=1,nspec
398	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
399	griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
400	xmass1(i,ks)/rhoiweightw
401	end do
402	endif
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)
410	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
414	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
415	end do
416	else
417	do ks=1,nspec
418	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
419	griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
420	xmass1(i,ks)/rhoiweightw
421	end do
422	endif
423	endif
424
425	if ((jy.ge.0).and.(jy.le.numygridn-1)) then
426	w=(1.-wx)*wy
427	if (DRYBKDEP.or.WETBKDEP) then
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)+ &
431	xmass1(i,ks)/(1-xscav_frac1(i,ks))/rhoiweightw*xscav_frac1(i,ks)
432	end do
433	else
434	do ks=1,nspec
435	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
436	griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
437	xmass1(i,ks)/rhoiweightw
438	end do
439	endif
440	endif
441	endif
442	endif
443	endif
444	endif
445	20 continue
446	end do
447
448	!***********************************************************************
449	! 2. Evaluate concentrations at receptor points, using the kernel method
450	!***********************************************************************
451
452	do n=1,numreceptor
453
454
455	! Reset concentrations
456	!*********************
457
458	do ks=1,nspec
459	c(ks)=0.
460	end do
461
462
463	! Estimate concentration at receptor
464	!***********************************
465
466	do i=1,numpart
467
468	if (itra1(i).ne.itime) goto 40
469	itage=abs(itra1(i)-itramem(i))
470
471	hz=min(50.+0.3*sqrt(real(itage)),hzmax)
472	zd=ztra1(i)/hz
473	if (zd.gt.1.) goto 40 ! save computing time, leave loop
474
475	hx=min((0.29+2.222e-3sqrt(real(itage)))dx+ &
476	real(itage)*1.2e-5,hxmax) ! 80 km/day
477	xd=(xtra1(i)-xreceptor(n))/hx
478	if (xd*xd.gt.1.) goto 40 ! save computing time, leave loop
479
480	hy=min((0.18+1.389e-3sqrt(real(itage)))dy+ &
481	real(itage)*7.5e-6,hymax) ! 80 km/day
482	yd=(ytra1(i)-yreceptor(n))/hy
483	if (yd*yd.gt.1.) goto 40 ! save computing time, leave loop
484	h=hxhyhz
485
486	r2=xdxd+ydyd+zd*zd
487	if (r2.lt.1.) then
488	xkern=factor*(1.-r2)
489	do ks=1,nspec
490	c(ks)=c(ks)+xmass1(i,ks)*xkern/h
491	end do
492	endif
493	40 continue
494	end do
495
496	do ks=1,nspec
497	creceptor(n,ks)=creceptor(n,ks)+2.weightc(ks)/receptorarea(n)
498	end do
499	end do
500
501	end subroutine conccalc

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

Download in other formats:

Original Format