Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

outgrid_init_irreg.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 20.2 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine outgrid_init_irreg
24	!*******************************************************************************
25	! *
26	! Note: This is the FLEXPART_WRF version of subroutine outgrid_init. *
27	! The computational grid is the WRF x-y grid rather than lat-lon. *
28	! *
29	! This routine calculates, for each grid cell of the output grid, the *
30	! volume, the surface area, and the areas of the northward and eastward *
31	! facing surfaces. *
32	! *
33	! Author: A. Stohl *
34	! *
35	! 7 August 2002 *
36	! *
37	! 26 Oct 2005, R. Easter - changes in gridarea, areaeast, areanorth *
38	! associated with WRF horizontal grid. *
39	! Dec 2005, R. Easter - changed names of "lon0" & "lat0" variables *
40	! *
41	!*******************************************************************************
42	! *
43	! Variables: *
44	! *
45	! area surface area of all output grid cells *
46	! areaeast eastward facing wall area of all output grid cells *
47	! areanorth northward facing wall area of all output grid cells *
48	! volume volumes of all output grid cells *
49	! *
50	!*******************************************************************************
51
52	use flux_mod
53	use oh_mod
54	use unc_mod
55	use outg_mod
56	use par_mod
57	use com_mod
58
59	! include 'includepar'
60	! include 'includecom'
61	implicit none
62	integer :: ix,jy,kz,k,i,nage,l,iix,jjy,ixp,jyp,i1,j1,j,ngrid
63	! real ylat,gridarea,ylatp,ylatm,hzone,cosfact,cosfactm,cosfactp
64	real :: ymet,gridarea,xl1,xl2,yl1,yl2,m1,m2,tmpx,tmpy
65	real :: xmet,xl,yl,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn,oroh
66	integer :: ks,kp,stat,ix2,jy2
67	real,parameter :: eps=nxmax/3.e5
68	real :: lon2(4),lat2(4)
69	real ( kind = 8 ) :: sphere01_polygon_area,haversine,area1
70
71
72
73	! Compute surface area and volume of each grid cell: area, volume;
74	! and the areas of the northward and eastward facing walls: areaeast, areanorth
75	!***********************************************************************
76	! do jy=0,10
77	! print*,areamet(1,jy),areamet(2,jy),areamet2(2,jy)
78	! enddo
79	do jy=0,numygrid-1
80
81	! ylat=outlat0+(real(jy)+0.5)*dyout
82	! ylatp=ylat+0.5*dyout
83	! ylatm=ylat-0.5*dyout
84	! if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
85	! hzone=dyoutr_earthpi180
86	! else
87	!
88	!C Calculate area of grid cell with formula M=2piRhdx/360,
89	!C see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
90	!*************************************************************
91	!
92	! cosfact=cos(ylatpi180)r_earth
93	! cosfactp=cos(ylatppi180)r_earth
94	! cosfactm=cos(ylatmpi180)r_earth
95	! if (cosfactp.lt.cosfactm) then
96	! hzone=sqrt(r_earth2-cosfactp2)-
97	! + sqrt(r_earth2-cosfactm2)
98	! else
99	! hzone=sqrt(r_earth2-cosfactm2)-
100	! + sqrt(r_earth2-cosfactp2)
101	! endif
102	! endif
103	!
104	!C Surface are of a grid cell at a latitude ylat
105	!***********************************************
106	!
107	! gridarea=2.pir_earthhzonedxout/360.
108
109	! for FLEXPART_WRF, dx & dy are in meters, and no cos(lat) is needed
110	! ??? maybe should incorporate map factor here,
111	! and also for areaeast & areanorth ???
112
113	do ix=0,numxgrid-1
114	! gridarea=dxout*dyout ! what is needed is the true area based on map factors
115	!JB: find a way to get the area between 2 output grid cell using areamet
116	! xl1=(real(ix)*dxout+out_xm0)/dx
117	! yl1=(real(jy)*dyout+out_ym0)/dy
118	! xl2=(real(ix+1)*dxout+out_xm0)/dx
119	! yl2=(real(jy+1)*dyout+out_ym0)/dy
120	!! xr=out_xm0+real(numxgrid)*dxout
121	! m1=0.5*(m_x(int(xl1),int(yl1),1)+m_x(int(xl2),int(yl1),1))
122	! m2=0.5*(m_y(int(xl1),int(yl1),1)+m_x(int(xl1),int(yl2),1))
123	! area(ix,jy)=dxoutm1dyout*m2
124
125	! A more precise method
126	tmpx=out_xm0+(float(ix))*dxout
127	tmpy=out_ym0+(float(jy))*dyout
128	call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(1),lat2(1))
129	tmpx=out_xm0+(float(ix+1))*dxout
130	tmpy=out_ym0+(float(jy))*dyout
131	call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(2),lat2(2))
132	tmpx=out_xm0+(float(ix+1))*dxout
133	tmpy=out_ym0+(float(jy+1))*dyout
134	call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(3),lat2(3))
135	tmpx=out_xm0+(float(ix))*dxout
136	tmpy=out_ym0+(float(jy+1))*dyout
137	call xymeter_to_ll_wrf_out(tmpx,tmpy,lon2(4),lat2(4))
138	area1=sphere01_polygon_area ( 4, real(lat2,kind=8), real(lon2,kind=8) )
139	area(ix,jy)=real(area1)6370000.6370000./coefdx/coefdx
140	!
141	! Volume = area x box height
142	!***************************
143
144	volume(ix,jy,1)=area(ix,jy)*outheight(1)
145
146	! for FLEXPART_WRF, dx & dy are in meters, and no cos(lat) is needed
147	! areaeast(ix,jy,1)=dyoutr_earthpi180*outheight(1)
148	! areanorth(ix,jy,1)=cos(ylatpi180)dxoutr_earthpi180*
149	! + outheight(1)
150	areaeast(ix,jy,1)=dyout*outheight(1)
151	areanorth(ix,jy,1)=dxout*outheight(1)
152
153	do kz=2,numzgrid
154
155	! areaeast(ix,jy,kz)=dyoutr_earthpi180*
156	! + (outheight(kz)-outheight(kz-1))
157	! areanorth(ix,jy,kz)=cos(ylatpi180)dxoutr_earthpi180*
158	! + (outheight(kz)-outheight(kz-1))
159	areaeast(ix,jy,kz)=dyout*(outheight(kz)-outheight(kz-1))
160	areanorth(ix,jy,kz)=dxout*(outheight(kz)-outheight(kz-1))
161
162	volume(ix,jy,kz)=area(ix,jy)*(outheight(kz)-outheight(kz-1))
163	end do
164	end do
165	end do
166
167
168
169
170	!******************************************************************
171	! Determine average height of model topography in output grid cells
172	!******************************************************************
173
174	! Loop over all output grid cells
175	!********************************
176
177	do jjy=0,numygrid-1
178	do iix=0,numxgrid-1
179	oroh=0.
180
181	! Take 100 samples of the topography in every grid cell
182	!******************************************************
183
184	do j1=1,10
185	! for FLEXPART_WRF, x & y coords are in meters,
186	! and the lon & lat variables below are in meters.
187	ymet=out_ym0+(real(jjy)+real(j1)/10.-0.05)*dyout
188	yl=(ymet-ymet0)/dy
189	do i1=1,10
190	xmet=out_xm0+(real(iix)+real(i1)/10.-0.05)*dxout
191	xl=(xmet-xmet0)/dx
192
193	! Determine the nest we are in
194	!*****************************
195
196	ngrid=0
197	do j=numbnests,1,-1
198	if ((xl.gt.xln(j)).and.(xl.lt.xrn(j)).and. &
199	(yl.gt.yln(j)).and.(yl.lt.yrn(j))) then
200	ngrid=j
201	goto 43
202	endif
203	end do
204	43 continue
205
206	! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
207	!************************************************************************************
208
209	if (ngrid.gt.0) then
210	xtn=(xl-xln(ngrid))*xresoln(ngrid)
211	ytn=(yl-yln(ngrid))*yresoln(ngrid)
212	ix=int(xtn)
213	jy=int(ytn)
214	ddy=ytn-real(jy)
215	ddx=xtn-real(ix)
216	else
217	ix=int(xl)
218	jy=int(yl)
219	ddy=yl-real(jy)
220	ddx=xl-real(ix)
221	endif
222	ixp=ix+1
223	jyp=jy+1
224	rddx=1.-ddx
225	rddy=1.-ddy
226	p1=rddx*rddy
227	p2=ddx*rddy
228	p3=rddx*ddy
229	p4=ddx*ddy
230
231	if (ngrid.gt.0) then
232	oroh=oroh+p1*oron(ix ,jy ,ngrid) &
233	+ p2*oron(ixp,jy ,ngrid) &
234	+ p3*oron(ix ,jyp,ngrid) &
235	+ p4*oron(ixp,jyp,ngrid)
236	else
237	oroh=oroh+p1*oro(ix ,jy) &
238	+ p2*oro(ixp,jy) &
239	+ p3*oro(ix ,jyp) &
240	+ p4*oro(ixp,jyp)
241	endif
242	end do
243	end do
244
245	! Divide by the number of samples taken
246	!**************************************
247
248	oroout(iix,jjy)=oroh/100.
249	end do
250	end do
251
252
253	! if necessary allocate flux fields
254	if (iflux.eq.1) then
255	allocate(flux(6,0:numxgrid-1,0:numygrid-1,numzgrid, &
256	1:nspec,1:maxpointspec_act,1:nageclass),stat=stat)
257	if (stat.ne.0) write(,)'ERROR: could not allocate flux array '
258	endif
259
260	!write (,) 'allocating: in a sec',OHREA
261	if (OHREA.eqv..TRUE.) then
262	! write (,) 'allocating: ',maxxOH,maxyOH,maxzOH
263	allocate(OH_field(12,0:maxxOH-1,0:maxyOH-1,maxzOH) &
264	,stat=stat)
265	if (stat.ne.0) write(,)'ERROR: could not allocate OH array '
266	allocate(OH_field_height(7) &
267	,stat=stat)
268	if (stat.ne.0) write(,)'ERROR: could not allocate OH array '
269	endif
270	! gridunc,griduncn uncertainty of outputted concentrations
271	! print*,'gridunc allocated'
272	allocate(gridunc(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
273	maxpointspec_act,nclassunc,maxageclass),stat=stat)
274	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
275	if (ldirect.gt.0) then
276	allocate(wetgridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
277	maxpointspec_act,nclassunc,maxageclass),stat=stat)
278	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
279	allocate(drygridunc(0:numxgrid-1,0:numygrid-1,maxspec, &
280	maxpointspec_act,nclassunc,maxageclass),stat=stat)
281	allocate(drygridunc2(0:numxgrid-1,0:numygrid-1,maxspec, &
282	maxpointspec_act,nclassunc,maxageclass),stat=stat)
283	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
284	endif
285	!write (,) 'Dimensions for fields', numxgrid,numygrid, &
286	! maxspec,maxpointspec_act,nclassunc,maxageclass
287
288	! print*,'alloc gridunc',numxgrid-1,numygrid-1,numzgrid,maxspec, &
289	! maxpointspec_act,nclassunc,maxageclass
290
291
292	write (,) ' Allocating fields for nested and global output (x,y): ', &
293	max(numxgrid,numxgridn),max(numygrid,numygridn)
294
295	! allocate fields for concoutput with maximum dimension of outgrid
296	! and outgrid_nest
297	allocate(gridsigma(0:max(numxgrid,numxgridn)-1, &
298	0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
299	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
300	allocate(grid(0:max(numxgrid,numxgridn)-1, &
301	0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
302	allocate(grid2(0:max(numxgrid,numxgridn)-1, &
303	0:max(numygrid,numygridn)-1,numzgrid,maxpointspec_act),stat=stat)
304	allocate(grid3(0:max(numxgrid,numxgridn)-1, &
305	0:max(numygrid,numygridn)-1,numzgrid,maxpointspec_act),stat=stat)
306	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
307	allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
308	0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
309	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
310
311	allocate(factor3d(0:max(numxgrid,numxgridn)-1, &
312	0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
313	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
314	allocate(sparse_dump_r(max(numxgrid,numxgridn)* &
315	max(numygrid,numygridn)*numzgrid),stat=stat)
316	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
317	allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
318	max(numygrid,numygridn)*numzgrid),stat=stat)
319	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
320
321	! deposition fields are only allocated for forward runs
322	if (ldirect.gt.0) then
323	allocate(wetgridsigma(0:max(numxgrid,numxgridn)-1, &
324	0:max(numygrid,numygridn)-1),stat=stat)
325	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
326	allocate(drygridsigma(0:max(numxgrid,numxgridn)-1, &
327	0:max(numygrid,numygridn)-1),stat=stat)
328	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
329	allocate(wetgrid(0:max(numxgrid,numxgridn)-1, &
330	0:max(numygrid,numygridn)-1),stat=stat)
331	allocate(wetgrid2(0:max(numxgrid,numxgridn)-1, &
332	0:max(numygrid,numygridn)-1,maxpointspec_act),stat=stat)
333	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
334	allocate(drygrid(0:max(numxgrid,numxgridn)-1, &
335	0:max(numygrid,numygridn)-1),stat=stat)
336	allocate(drygrid2(0:max(numxgrid,numxgridn)-1, &
337	0:max(numygrid,numygridn)-1,maxpointspec_act),stat=stat)
338	if (stat.ne.0) write(,)'ERROR: could not allocate gridunc'
339	endif
340
341	! Initial condition field
342
343	if (linit_cond.gt.0) then
344	allocate(init_cond(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
345	maxpointspec_act),stat=stat)
346	if (stat.ne.0) write(,)'ERROR: could not allocate init_cond'
347	endif
348
349	!************************
350	! Initialize output grids
351	!************************
352
353	do ks=1,nspec
354	do kp=1,maxpointspec_act
355	do i=1,numreceptor
356	! Receptor points
357	creceptor(i,ks)=0.
358	end do
359	do nage=1,nageclass
360	do jy=0,numygrid-1
361	do ix=0,numxgrid-1
362	do l=1,nclassunc
363	! Deposition fields
364	if (ldirect.gt.0) then
365	wetgridunc(ix,jy,ks,kp,l,nage)=0.
366	drygridunc(ix,jy,ks,kp,l,nage)=0.
367	! drygridunc2(ix,jy,ks,kp,l,nage)=0.
368	endif
369	do kz=1,numzgrid
370	if (iflux.eq.1) then
371	! Flux fields
372	do i=1,5
373	flux(i,ix,jy,kz,ks,kp,nage)=0.
374	end do
375	endif
376	! Initial condition field
377	if ((l.eq.1).and.(nage.eq.1).and.(linit_cond.gt.0)) &
378	init_cond(ix,jy,kz,ks,kp)=0.
379	! Concentration fields
380	gridunc(ix,jy,kz,ks,kp,l,nage)=0.
381	end do
382	end do
383	end do
384	end do
385	end do
386	end do
387	end do
388
389
390
391	end subroutine outgrid_init_irreg
392
393	function sphere01_polygon_area ( n, lat, lon )
394
395	!*****************************************************************************80
396	!
397	!! SPHERE01_POLYGON_AREA returns the area of a spherical polygon.
398	!
399	! Discussion:
400	!
401	! On a unit sphere, the area of a spherical polygon with N sides
402	! is equal to the spherical excess:
403	!
404	! E = sum ( interior angles ) - ( N - 2 ) * pi.
405	!
406	! On a sphere with radius R, the area is the spherical excess multiplied
407	! by R * R.
408	!
409	! The code was revised in accordance with suggestions in Carvalho and
410	! Cavalcanti.
411	!
412	! Licensing:
413	!
414	! This code is distributed under the GNU LGPL license.
415	!
416	! Modified:
417	!
418	! 12 August 2005
419	!
420	! Author:
421	!
422	! Original C version by Robert Miller.
423	! FORTRAN90 version by John Burkardt.
424	!
425	! Reference:
426	!
427	! Paulo Cezar Pinto Carvalho, Paulo Roma Cavalcanti,
428	! Point in Polyhedron Testing Using Spherical Polygons,
429	! in Graphics Gems V,
430	! edited by Alan Paeth,
431	! Academic Press, 1995,
432	! ISBN: 0125434553,
433	! LC: T385.G6975.
434	!
435	! Robert Miller,
436	! Computing the Area of a Spherical Polygon,
437	! Graphics Gems, Volume IV, pages 132-138,
438	! Edited by Paul Heckbert,
439	! Academic Press, 1994, T385.G6974.
440	!
441	! Eric Weisstein,
442	! "Spherical Polygon",
443	! CRC Concise Encyclopedia of Mathematics,
444	! CRC Press, 1999.
445	!
446	! Parameters:
447	!
448	! Input, integer ( kind = 4 ) N, the number of vertices.
449	!
450	! Input, real ( kind = 8 ) LAT[N], LON[N], the latitudes and longitudes
451	! of the vertices of the spherical polygon.
452	!
453	! Output, real ( kind = 8 ) SPHERE01_POLYGON_AREA, the area of the
454	! spherical polygon, measured in spherical radians.
455	!
456	implicit none
457
458	integer ( kind = 4 ) n
459
460	real ( kind = 8 ) a
461	real ( kind = 8 ) area
462	real ( kind = 8 ) b
463	real ( kind = 8 ) beta1
464	real ( kind = 8 ) beta2
465	real ( kind = 8 ) c
466	real ( kind = 8 ) cos_b1
467	real ( kind = 8 ) cos_b2
468	real ( kind = 8 ) excess
469	real ( kind = 8 ) hav_a
470	real ( kind = 8 ) haversine
471	integer ( kind = 4 ) j
472	integer ( kind = 4 ) k
473	real ( kind = 8 ) lam
474	real ( kind = 8 ) lam1
475	real ( kind = 8 ) lam2
476	real ( kind = 8 ) lat(n)
477	real ( kind = 8 ) lon(n)
478	real ( kind = 8 ), parameter :: pi_half = 1.5707963267948966192313D+00
479	real ( kind = 8 ) s
480	real ( kind = 8 ) sphere01_polygon_area
481	real ( kind = 8 ) t
482	real ( kind = 8 ),parameter :: degrees_to_radians=3.141592653589793D+00 / 180.0D+00
483
484	area = 0.0D+00
485
486	do j=1,n
487	lon(j)=lon(j)*degrees_to_radians
488	lat(j)=lat(j)*degrees_to_radians
489	enddo
490	do j = 1, n + 1
491	! do j = 1, n
492
493	if ( j == 1 ) then
494	lam1 = lon(j)
495	beta1 = lat(j)
496	lam2 = lon(j+1)
497	beta2 = lat(j+1)
498	cos_b1 = cos ( beta1 )
499	cos_b2 = cos ( beta2 )
500	else
501	! k = mod ( j + 1, n + 1 )
502	! k = mod ( j , n )
503	k=j
504	if (j.gt.n) k=1
505	lam1 = lam2
506	beta1 = beta2
507	lam2 = lon(k)
508	beta2 = lat(k)
509	! print*,'sphere',n,k,lon(k),lat(k)
510	cos_b1 = cos_b2
511	cos_b2 = cos ( beta2 )
512	end if
513
514	if ( lam1 /= lam2 ) then
515
516	hav_a = haversine ( beta2 - beta1 ) &
517	+ cos_b1 * cos_b2 * haversine ( lam2 - lam1 )
518	a = 2.0D+00 * asin ( sqrt ( hav_a ) )
519
520	b = pi_half - beta2
521	c = pi_half - beta1
522	s = 0.5D+00 * ( a + b + c )
523	!
524	! Given the three sides of a spherical triangle, we can use a formula
525	! to find the spherical excess.
526	!
527	t = tan ( s / 2.0D+00 ) * tan ( ( s - a ) / 2.0D+00 ) &
528	* tan ( ( s - b ) / 2.0D+00 ) * tan ( ( s - c ) / 2.0D+00 )
529
530	excess = abs ( 4.0D+00 * atan ( sqrt ( abs ( t ) ) ) )
531
532	if ( lam1 < lam2 ) then
533	lam = lam2 - lam1
534	else
535	lam = lam2 - lam1 + 4.0D+00 * pi_half
536	end if
537
538	if ( 2.0D+00 * pi_half < lam ) then
539	excess = -excess
540	end if
541
542	area = area + excess
543
544	end if
545
546	end do
547
548	sphere01_polygon_area = abs ( area )
549
550	return
551	end
552
553	function haversine ( a )
554
555	!*****************************************************************************80
556	!
557	!! HAVERSINE computes the haversine of an angle.
558	!
559	! Discussion:
560	!
561	! haversine(A) = ( 1 - cos ( A ) ) / 2
562	!
563	! The haversine is useful in spherical trigonometry.
564	!
565	! Licensing:
566	!
567	! This code is distributed under the GNU LGPL license.
568	!
569	! Modified:
570	!
571	! 02 July 2001
572	!
573	! Author:
574	!
575	! John Burkardt
576	!
577	! Parameters:
578	!
579	! Input, real ( kind = 8 ) A, the angle.
580	!
581	! Output, real ( kind = 8 ) HAVERSINE, the haversine of the angle.
582	!
583	implicit none
584
585	real ( kind = 8 ) a
586	real ( kind = 8 ) haversine
587
588	haversine = ( 1.0D+00 - cos ( a ) ) / 2.0D+00
589
590	return
591	end

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

Download in other formats:

Original Format