Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

calcpv.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: 14.4 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 calcpv(n,uuh,vvh,pvh)
24	! i i i o
25	!*******************************************************************************
26	! *
27	! Note: This is the FLEXPART_WRF version of subroutine calcpv. *
28	! The computational grid is the WRF x-y grid rather than lat-lon. *
29	! *
30	! Calculation of potential vorticity on 3-d grid. *
31	! *
32	! Author: P. James *
33	! 3 February 2000 *
34	! *
35	! Adaptation to FLEXPART, A. Stohl, 1 May 2000 *
36	! *
37	! 26 Oct 2005, R. Easter - changes associated with WRF horizontal grid. *
38	! For pressure use pph instead of (akz + bkzps)
39	! *** Note -- see ??? comments below regarding the pvh calculation. *
40	! 11 Nov 2005, R. Easter - fixed error involving xy to latlon *
41	! *
42	!*******************************************************************************
43	! *
44	! Variables: *
45	! n temporal index for meteorological fields (1 to 2) *
46	! *
47	! Constants: *
48	! *
49	!*******************************************************************************
50
51	! include 'includepar'
52	! include 'includecom'
53	!
54	! integer n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
55	! integer jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
56	! integer nlck
57	! real vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
58	! real theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
59	! real thup,thdn,eps,p0
60	! parameter(eps=1.e-5,p0=101325)
61	! real uuh(0:nxmax-1,0:nymax-1,nuvzmax)
62	! real vvh(0:nxmax-1,0:nymax-1,nuvzmax)
63	! real pvh(0:nxmax-1,0:nymax-1,nuvzmax)
64	use par_mod
65	use com_mod
66
67	implicit none
68
69	integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
70	integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
71	integer :: nlck
72	real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
73	real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
74	real :: pvavr,ppml(nuvzmax)
75	real :: thup,thdn
76	real,parameter :: eps=1.e-5, p0=101325
77	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
78	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
79	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
80
81
82
83	real :: dumlon,dumlat
84	real :: thh(0:nxmax-1,0:nymax-1,nuvzmax)
85	real :: altit(nuvzmax)
86
87	! Set number of levels to check for adjacent theta
88	nlck=nuvz/3
89	! FLEXPART_WRF -- altit is never used, so don't calculate it
90	! do 5 k=1,nuvz
91	! altit(k)=akz(k)/p0+bkz(k)
92	!5 continue
93	! *** Precalculate all theta levels for efficiency
94	do jy=0,nymin1
95	do kl=1,nuvz
96	do ix=0,nxmin1
97	! FLEXPART_WRF -- use pph here
98	! ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
99	ppmk=pph(ix,jy,kl,n)
100	thh(ix,jy,kl)=tth(ix,jy,kl,n)(100000./ppmk)*kappa
101	enddo
102	enddo
103	enddo
104	!
105	! Loop over entire grid
106	!**********************
107	do jy=0,nymin1
108	if (sglobal.and.jy.eq.0) goto 10
109	if (nglobal.and.jy.eq.nymin1) goto 10
110
111	! for FLEXPART_WRF, x & y coords are in meters
112	! and true latitude varies with both i and j
113	! phi = (ylat0 + jy * dy) * pi / 180.
114	! f = 0.00014585 * sin(phi)
115	! tanphi = tan(phi)
116	! cosphi = cos(phi)
117
118	! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
119	jyvp=jy+1
120	jyvm=jy-1
121	if (jy.eq.0) jyvm=0
122	if (jy.eq.nymin1) jyvp=nymin1
123	! Define absolute gap length
124	jumpy=2
125	if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
126	if (sglobal.and.jy.eq.1) then
127	jyvm=1
128	jumpy=1
129	end if
130	if (nglobal.and.jy.eq.ny-2) then
131	jyvp=ny-2
132	jumpy=1
133	end if
134	juy=jumpy
135	!
136	do ix=0,nxmin1
137
138	! for FLEXPART_WRF, x & y coords are in meters,
139	! and true latitude varies with both i and j
140	call xyindex_to_ll_wrf( &
141	0, real(ix), real(jy), dumlon, dumlat )
142	phi = dumlat * pi / 180.
143	f = 0.00014585 * sin(phi)
144	tanphi = tan(phi)
145	cosphi = cos(phi)
146
147	! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
148	ixvp=ix+1
149	ixvm=ix-1
150	jumpx=2
151	if (xglobal) then
152	ivrp=ixvp
153	ivrm=ixvm
154	if (ixvm.lt.0) ivrm=ixvm+nxmin1
155	if (ixvp.ge.nx) ivrp=ixvp-nx+1
156	else
157	if (ix.eq.0) ixvm=0
158	if (ix.eq.nxmin1) ixvp=nxmin1
159	ivrp=ixvp
160	ivrm=ixvm
161	! Define absolute gap length
162	if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
163	end if
164	jux=jumpx
165	! Precalculate pressure values for efficiency
166	do kl=1,nuvz
167	! FLEXPART_WRF -- use pph here
168	! ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
169	ppml(kl)=pph(ix,jy,kl,n)
170	enddo
171	!
172	! Loop over the vertical
173	!***********************
174
175	do kl=1,nuvz
176	theta=thh(ix,jy,kl)
177	klvrp=kl+1
178	klvrm=kl-1
179	klpt=kl
180	! If top or bottom level, dthetadp is evaluated between the current
181	! level and the level inside, otherwise between level+1 and level-1
182	!
183	if (klvrp.gt.nuvz) klvrp=nuvz
184	if (klvrm.lt.1) klvrm=1
185	thetap=thh(ix,jy,klvrp)
186	thetam=thh(ix,jy,klvrm)
187	dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
188
189	! Compute vertical position at pot. temperature surface on subgrid
190	! and the wind at that position
191	!*****************************************************************
192	! a) in x direction
193	ii=0
194	do i=ixvm,ixvp,jumpx
195	ivr=i
196	if (xglobal) then
197	if (i.lt.0) ivr=ivr+nxmin1
198	if (i.ge.nx) ivr=ivr-nx+1
199	end if
200	ii=ii+1
201	! Search adjacent levels for current theta value
202	! Spiral out from current level for efficiency
203	kup=klpt-1
204	kdn=klpt
205	kch=0
206	40 continue
207	! Upward branch
208	kup=kup+1
209	if (kch.ge.nlck) goto 21 ! No more levels to check,
210	! ! and no values found
211	if (kup.ge.nuvz) goto 41
212	kch=kch+1
213	k=kup
214	thdn=thh(ivr,jy,k)
215	thup=thh(ivr,jy,k+1)
216	if (((thdn.ge.theta).and.(thup.le.theta)).or. &
217	((thdn.le.theta).and.(thup.ge.theta))) then
218	dt1=abs(theta-thdn)
219	dt2=abs(theta-thup)
220	dt=dt1+dt2
221	if (dt.lt.eps) then ! Avoid division by zero error
222	dt1=0.5 ! G.W., 10.4.1996
223	dt2=0.5
224	dt=1.0
225	endif
226	vx(ii)=(vvh(ivr,jy,k)dt2+vvh(ivr,jy,k+1)dt1)/dt
227	goto 20
228	endif
229	41 continue
230	! Downward branch
231	kdn=kdn-1
232	if (kdn.lt.1) goto 40
233	kch=kch+1
234	k=kdn
235	thdn=thh(ivr,jy,k)
236	thup=thh(ivr,jy,k+1)
237	if (((thdn.ge.theta).and.(thup.le.theta)).or. &
238	((thdn.le.theta).and.(thup.ge.theta))) then
239	dt1=abs(theta-thdn)
240	dt2=abs(theta-thup)
241	dt=dt1+dt2
242	if (dt.lt.eps) then ! Avoid division by zero error
243	dt1=0.5 ! G.W., 10.4.1996
244	dt2=0.5
245	dt=1.0
246	endif
247	vx(ii)=(vvh(ivr,jy,k)dt2+vvh(ivr,jy,k+1)dt1)/dt
248	goto 20
249	endif
250	goto 40
251	! This section used when no values were found
252	21 continue
253	! Must use vv at current level and long. jux becomes smaller by 1
254	vx(ii)=vvh(ix,jy,kl)
255	jux=jux-1
256	! Otherwise OK
257	20 continue
258	end do
259
260	if (jux.gt.0) then
261	! for FLEXPART_WRF, dx & dy are in meters.
262	! dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
263	dvdx=(vx(2)-vx(1))/real(jux)/dx
264	else
265	dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
266	! dvdx=dvdx/real(jumpx)/(dx*pi/180.)
267	dvdx=dvdx/real(jumpx)/dx
268	! Only happens if no equivalent theta value
269	! can be found on either side, hence must use values
270	! from either side, same pressure level.
271	end if
272
273	! b) in y direction
274
275	jj=0
276	do j=jyvm,jyvp,jumpy
277	jj=jj+1
278	! Search adjacent levels for current theta value
279	! Spiral out from current level for efficiency
280	kup=klpt-1
281	kdn=klpt
282	kch=0
283	70 continue
284	! Upward branch
285	kup=kup+1
286	if (kch.ge.nlck) goto 51 ! No more levels to check,
287	! ! and no values found
288	if (kup.ge.nuvz) goto 71
289	kch=kch+1
290	k=kup
291	thdn=thh(ix,j,k)
292	thup=thh(ix,j,k+1)
293	if (((thdn.ge.theta).and.(thup.le.theta)).or. &
294	((thdn.le.theta).and.(thup.ge.theta))) then
295	dt1=abs(theta-thdn)
296	dt2=abs(theta-thup)
297	dt=dt1+dt2
298	if (dt.lt.eps) then ! Avoid division by zero error
299	dt1=0.5 ! G.W., 10.4.1996
300	dt2=0.5
301	dt=1.0
302	endif
303	uy(jj)=(uuh(ix,j,k)dt2+uuh(ix,j,k+1)dt1)/dt
304	goto 50
305	endif
306	71 continue
307	! Downward branch
308	kdn=kdn-1
309	if (kdn.lt.1) goto 70
310	kch=kch+1
311	k=kdn
312	thdn=thh(ix,j,k)
313	thup=thh(ix,j,k+1)
314	if (((thdn.ge.theta).and.(thup.le.theta)).or. &
315	((thdn.le.theta).and.(thup.ge.theta))) then
316	dt1=abs(theta-thdn)
317	dt2=abs(theta-thup)
318	dt=dt1+dt2
319	if (dt.lt.eps) then ! Avoid division by zero error
320	dt1=0.5 ! G.W., 10.4.1996
321	dt2=0.5
322	dt=1.0
323	endif
324	uy(jj)=(uuh(ix,j,k)dt2+uuh(ix,j,k+1)dt1)/dt
325	goto 50
326	endif
327	goto 70
328	! This section used when no values were found
329	51 continue
330	! Must use uu at current level and lat. juy becomes smaller by 1
331	uy(jj)=uuh(ix,jy,kl)
332	juy=juy-1
333	! Otherwise OK
334	50 continue
335	end do
336
337	if (juy.gt.0) then
338	! for FLEXPART_WRF, dx & dy are in meters.
339	! dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
340	dudy=(uy(2)-uy(1))/real(juy)/dy
341	else
342	dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
343	! dudy=dudy/real(jumpy)/(dy*pi/180.)
344	dudy=dudy/real(jumpy)/dy
345	end if
346	!
347
348	! for FLEXPART_WRF, dx & dy are in meters.
349	! don't need to divide by r_earth when doing d/dy
350	! don't need to divide by r_earth*cosphi when doing d/dx
351	! ??? I don't understand the uuh*tanphi term, but leave it in for now ???
352	! ??? What is the "-1.e6" factor ???
353	!
354	! pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy
355	! + +uuh(ix,jy,kl)tanphi)/r_earth)(-1.e6)*9.81
356	pvh(ix,jy,kl)=dthetadp*( f + dvdx - dudy &
357	+ (uuh(ix,jy,kl)tanphi/r_earth) )(-1.e6)*9.81
358
359	!
360	! Resest jux and juy
361	jux=jumpx
362	juy=jumpy
363	end do
364	end do
365	10 continue
366	end do
367
368	!
369	! Fill in missing PV values on poles, if present
370	! Use mean PV of surrounding latitude ring
371	!
372	if (sglobal) then
373	do kl=1,nuvz
374	pvavr=0.
375	do ix=0,nxmin1
376	pvavr=pvavr+pvh(ix,1,kl)
377	end do
378	pvavr=pvavr/real(nx)
379	jy=0
380	do ix=0,nxmin1
381	pvh(ix,jy,kl)=pvavr
382	end do
383	end do
384	end if
385	if (nglobal) then
386	do kl=1,nuvz
387	pvavr=0.
388	do ix=0,nxmin1
389	pvavr=pvavr+pvh(ix,ny-2,kl)
390	end do
391	pvavr=pvavr/real(nx)
392	jy=nymin1
393	do ix=0,nxmin1
394	pvh(ix,jy,kl)=pvavr
395	end do
396	end do
397	end if
398
399	end subroutine calcpv
400

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

Download in other formats:

Original Format