Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

calcpar.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: 13.6 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 calcpar(n,uuh,vvh,pvh)
24	! i i i o
25	!*******************************************************************************
26	! *
27	! Computation of several boundary layer parameters needed for the *
28	! dispersion calculation and calculation of dry deposition velocities. *
29	! All parameters are calculated over the entire grid. *
30	! *
31	! Note: This is the FLEXPART_WRF version of subroutine calcpar. *
32	! *
33	! Author: A. Stohl *
34	! *
35	! 21 May 1995 *
36	! *
37	! ------------------------------------------------------------------ *
38	! Petra Seibert, Feb 2000: *
39	! convection scheme: *
40	! new variables in call to richardson *
41	! *
42	! Changes, Bernd C. Krueger, Feb. 2001:
43	! Variables tth and qvh (on eta coordinates) in common block
44	! *
45	! 17 Oct 2005 - R. Easter - added ierr in call to richardson *
46	! 18 Oct 2005 - J. Fast - limit ustar to < 5.0 m/s *
47	! -- Oct 2005 - R. Easter - use xy_to_ll_wrf to get latitude *
48	! use pph for calculating zlev *
49	! pass level-2 pph directly to obukhov *
50	! 11 Nov 2005 - R. Easter - changed name of "xy to latlon" routine *
51	! 15 Nov 2005 - R. Easter - pass pplev to richardson instead of akz,bkz *
52	! July 2012: J. Brioude: coded in fortran 90 and parallelized *
53	!*******************************************************************************
54	! *
55	! Variables: *
56	! n temporal index for meteorological fields (1 to 3) *
57	! *
58	! Constants: *
59	! *
60	! *
61	! Functions: *
62	! scalev computation of ustar *
63	! obukhov computatio of Obukhov length *
64	! *
65	!*******************************************************************************
66
67	use par_mod
68	use com_mod
69
70	implicit none
71
72	integer :: n,ix,jy,i,kz,lz,kzmin,ierr
73	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
74	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
75	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
76	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
77	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
78
79	! real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
80	! real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
81	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
82	real,parameter :: const=r_air/ga
83
84	real :: xlon,dumx,dumy,dumxb,dumyb,pplev(nuvzmax),hmixdummy
85
86	! Loop over entire grid
87	!**********************
88	! ientry = ientry + 1
89
90	!$OMP PARALLEL DEFAULT(SHARED) &
91	!$OMP PRIVATE(i,ix,jy,kz,lz,kzmin,tvold,pold,zold,zlev,tv,pint, &
92	!$OMP rh,ierr,subsceff,ulev,vlev,pplev,ttlev,qvlev,ol,altmin,xlon,ylat )
93	!$OMP DO
94	do jy=0,nymin1
95
96	do ix=0,nxmin1
97
98	! Set minimum height for tropopause
99	!**********************************
100
101	! FLEXPART_WRF - use this routine to get lat,lon
102	! ylat=ylat0+real(jy)*dy
103	call xyindex_to_ll_wrf( 0, real(ix), real(jy), xlon, ylat )
104
105	! if ( ((ix.eq.0) .or. (ix.eq.nxmin1) .or.
106	! & (ix.eq.nxmin1/2)) .and.
107	! & ((jy.eq.0) .or. (jy.eq.nymin1) .or.
108	! & (jy.eq.nymin1/2)) ) then
109	! if (ientry .eq. 1) then
110	! write(*,'(a,2i4,2f12.5)')
111	! & 'calcpar i,j, xlon,ylat', ix, jy, xlon, ylat
112	! write(*,'(a, 8x,2f12.5)')
113	! & ' dlon,dlat',
114	! & (xlon-xlon2d(ix,jy)), (ylat-ylat2d(ix,jy))
115	! call ll_to_xyindex_wrf(
116	! & xlon2d(ix,jy), ylat2d(ix,jy), dumx, dumy )
117	! write(*,'(a, 8x,2f12.5)')
118	! & ' dxkm,dykm',
119	! & ((dumx-ix)dx1.0e-3), ((dumy-jy)dy1.0e-3)
120	!
121	! if ((ix .eq. 0) .and. (jy .eq. 0)) then
122	! dumxb = 2.33
123	! dumyb = 3.44
124	! call xyindex_to_ll_wrf( 0, dumxb, dumyb, dumx, dumy )
125	! call ll_to_xyindex_wrf( dumx, dumy, dumx, dumy )
126	! write(*,'(a,2f5.2,2f12.5)')
127	! & 'xi,yj, dxkm,dykm', dumxb, dumyb,
128	! & ((dumx-dumxb)dx1.0e-3), ((dumy-dumyb)dy1.0e-3)
129	! dumxb = 4.55
130	! dumyb = 6.77
131	! call xyindex_to_ll_wrf( 0, dumxb, dumyb, dumx, dumy )
132	! call ll_to_xyindex_wrf( dumx, dumy, dumx, dumy )
133	! write(*,'(a,2f5.2,2f12.5)')
134	! & 'xi,yj, dxkm,dykm', dumxb, dumyb,
135	! & ((dumx-dumxb)dx1.0e-3), ((dumy-dumyb)dy1.0e-3)
136	! end if
137	!
138	! end if
139	! end if
140
141	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
142	altmin = 5000.
143	else
144	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
145	altmin=2500.+(40.-ylat)*125.
146	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
147	altmin=2500.+(40.+ylat)*125.
148	else
149	altmin=2500.
150	endif
151	endif
152
153	! 1) Calculation of friction velocity
154	!************************************
155	if ( (.not.strswitch)) then
156	ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
157	td2(ix,jy,1,n),surfstr(ix,jy,1,n))
158	endif
159	if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
160	! FLEXPART_WRF - limit ustar
161	if (ustar(ix,jy,1,n).ge.5.0) ustar(ix,jy,1,n)=5.0
162
163	! 2) Calculation of inverse Obukhov length scale
164	!***********************************************
165
166	! FLEXPART_WRF - pass k=2 pressure directly
167	! ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n),
168	! + tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
169	ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
170	tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n), &
171	pph(ix,jy,2,n) )
172	if (ol.ne.0.) then
173	oli(ix,jy,1,n)=1./ol
174	else
175	oli(ix,jy,1,n)=99999.
176	endif
177
178
179	! 3) Calculation of convective velocity scale and mixing height
180	!**************************************************************
181
182	do i=1,nuvz
183	ulev(i) =uuh(ix,jy,i)
184	vlev(i) =vvh(ix,jy,i)
185	pplev(i)=pph(ix,jy,i,n)
186	ttlev(i)=tth(ix,jy,i,n)
187	qvlev(i)=qvh(ix,jy,i,n)
188	zlev(i)=0.5*(zzh(ix,jy,i+1,n)+zzh(ix,jy,i,n))-zzh(ix,jy,1,n)
189	enddo
190	! FLEXPART_WRF - use & check ierr argument
191	! FLEXPART_WRF - pass pplev instead of akz,bkz
192	! call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev,
193	! + ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n),
194	! + td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
195	call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
196	ulev,vlev,nuvz, pplev,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
197	td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus, &
198	! td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus, &
199	ierr,sfc_option )
200	!JB
201	! no reflec
202	! hmix(ix,jy,1,n)=5000.
203
204	if (ierr .gt. 0) then
205	write(*,9500) 'warning', ix, jy
206	else if (ierr .lt. 0) then
207	write(*,9500) 'failure', ix, jy
208	stop
209	end if
210	9500 format( 'calcpar - richardson ', a, ' - ix,jy=', 2i5 )
211
212
213	if(lsubgrid.eq.1) then
214	subsceff=min(excessoro(ix,jy),hmixplus)
215	! subsceff=hmixplus
216	else
217	subsceff=0.
218	endif
219	!
220	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
221	!
222	hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
223	! print*,'hmix',hmix(ix,jy,1,n),subsceff
224	hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
225	hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
226
227	! 4) Calculation of dry deposition velocities
228	!********************************************
229
230	if (DRYDEP) then
231	z0(4)=0.016ustar(ix,jy,1,n)ustar(ix,jy,1,n)/ga
232	z0(9)=0.016ustar(ix,jy,1,n)ustar(ix,jy,1,n)/ga
233
234	! Calculate relative humidity at surface
235	!***************************************
236	rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
237
238	call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
239	tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
240	ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n),vd)
241
242	do i=1,nspec
243	vdep(ix,jy,i,n)=vd(i)
244	enddo
245	endif
246
247	!******************************************************
248	! Calculate height of thermal tropopause (Hoinka, 1997)
249	!******************************************************
250
251	! 1) Calculate altitudes of ECMWF model levels
252	!*********************************************
253
254	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
255	ps(ix,jy,1,n))
256	pold=ps(ix,jy,1,n)
257	zold=0.
258	! FLEXPART_WRF - set zlev(1)
259	zlev(1)=zold
260	do kz=2,nuvz
261	! FLEXPART_WRF - use pph for pressure
262	! pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
263	pint=pph(ix,jy,kz,n) ! pressure on model layers
264	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
265
266	if (abs(tv-tvold).gt.0.2) then
267	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
268	else
269	zlev(kz)=zold+constlog(pold/pint)tv
270	endif
271	tvold=tv
272	pold=pint
273	zold=zlev(kz)
274	enddo
275
276	! 2) Define a minimum level kzmin, from which upward the tropopause is
277	! searched for. This is to avoid inversions in the lower troposphere
278	! to be identified as the tropopause
279	!************************************************************************
280
281	do kz=1,nuvz
282	if (zlev(kz).ge.altmin) then
283	kzmin=kz
284	goto 45
285	endif
286	end do
287	45 continue
288
289	! 3) Search for first stable layer above minimum height that fulfills the
290	! thermal tropopause criterion
291	!************************************************************************
292
293	do kz=kzmin,nuvz
294	do lz=kz+1,nuvz
295	if ((zlev(lz)-zlev(kz)).gt.2000.) then
296	if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
297	(zlev(lz)-zlev(kz))).lt.0.002) then
298	tropopause(ix,jy,1,n)=zlev(kz)
299	goto 51
300	endif
301	goto 50
302	endif
303	end do
304	50 continue
305	end do
306	51 continue
307
308
309	end do
310	end do
311	!$OMP END DO
312	!$OMP END PARALLEL
313
314	! Calculation of potential vorticity on 3-d grid, if plume trajectory mode is used
315	!*********************************************************************************
316
317	if ((iout.eq.4).or.(iout.eq.5).or.(mdomainfill.eq.2)) then
318	call calcpv(n,uuh,vvh,pvh)
319	endif
320
321
322	end subroutine calcpar
323

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

Download in other formats:

Original Format