Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

calcpar_nests.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago
Added parallel version of Flexpart91
File size: 9.7 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 calcpar_nests(n,uuhn,vvhn,pvhn)
23	! i i i o
24	!*****************************************************************************
25	! *
26	! Computation of several boundary layer parameters needed for the *
27	! dispersion calculation and calculation of dry deposition velocities. *
28	! All parameters are calculated over the entire grid. *
29	! This routine is similar to calcpar, but is used for the nested grids. *
30	! *
31	! Author: A. Stohl *
32	! *
33	! 8 February 1999 *
34	! *
35	! ------------------------------------------------------------------ *
36	! Petra Seibert, Feb 2000: *
37	! convection scheme: *
38	! new variables in call to richardson *
39	! *
40	!*****************************************************************************
41	! Changes, Bernd C. Krueger, Feb. 2001:
42	! Variables tth and qvh (on eta coordinates) in common block
43	!*****************************************************************************
44	! *
45	! Variables: *
46	! n temporal index for meteorological fields (1 to 3) *
47	! *
48	! Constants: *
49	! *
50	! *
51	! Functions: *
52	! scalev computation of ustar *
53	! obukhov computatio of Obukhov length *
54	! *
55	!*****************************************************************************
56
57	use par_mod
58	use com_mod
59
60	implicit none
61
62	integer :: n,ix,jy,i,l,kz,lz,kzmin
63	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
64	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
65	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
66	real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
67	real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
68	real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
69	real,parameter :: const=r_air/ga
70
71
72	! Loop over all nests
73	!********************
74
75	nestloop : do l=1,numbnests
76
77	! Loop over entire grid
78	!**********************
79
80	!$OMP PARALLEL DEFAULT(none) &
81	!$OMP PRIVATE(jy, ix, ylat, altmin, ol, ulev, i, vlev, ttlev, qvlev, &
82	!$OMP hmixplus, kz, tvold, pold, zold, tv, pint, zlev, subsceff, rh, vd, &
83	!$OMP kzmin) &
84	!$OMP SHARED(l, nyn, nxn, dyn, ylat0n, n, ustarn, surfstrn, td2n, tt2n, psn, &
85	!$OMP akm, bkm, sshfn, olin, nuvz, tthn, uuhn, vvhn, qvhn, hmixn, wstarn, akz, bkz, &
86	!$OMP lsubgrid, excessoron, DRYDEP, lsprecn, convprecn, vdepn, sdn, ssrn, tropopausen, &
87	!$OMP nspec) &
88	!$OMP COPYIN(z0)
89
90	#if (defined STATIC_SCHED)
91	!$OMP DO SCHEDULE(static)
92	#else
93	!$OMP DO SCHEDULE(dynamic, max(1,(nyn(l)-1)/10))
94	#endif
95
96	latloop : do jy=0,nyn(l)-1
97
98	! Set minimum height for tropopause
99	!**********************************
100
101	ylat=ylat0n(l)+real(jy)*dyn(l)
102	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
103	altmin = 5000.
104	else
105	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
106	altmin=2500.+(40.-ylat)*125.
107	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
108	altmin=2500.+(40.+ylat)*125.
109	else
110	altmin=2500.
111	endif
112	endif
113
114	do ix=0,nxn(l)-1
115
116	! 1) Calculation of friction velocity
117	!************************************
118
119	ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
120	td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
121
122	! 2) Calculation of inverse Obukhov length scale
123	!***********************************************
124
125	ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
126	td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
127	sshfn(ix,jy,1,n,l),akm,bkm)
128	if (ol.ne.0.) then
129	olin(ix,jy,1,n,l)=1./ol
130	else
131	olin(ix,jy,1,n,l)=99999.
132	endif
133
134
135	! 3) Calculation of convective velocity scale and mixing height
136	!**************************************************************
137
138	do i=1,nuvz
139	ulev(i)=uuhn(ix,jy,i,l)
140	vlev(i)=vvhn(ix,jy,i,l)
141	ttlev(i)=tthn(ix,jy,i,n,l)
142	qvlev(i)=qvhn(ix,jy,i,n,l)
143	end do
144
145	call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
146	qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
147	tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
148	wstarn(ix,jy,1,n,l),hmixplus)
149
150	if(lsubgrid.eq.1) then
151	subsceff=min(excessoron(ix,jy,l),hmixplus)
152	else
153	subsceff=0
154	endif
155	!
156	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
157	!
158	hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
159	hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
160	hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
161
162
163	! 4) Calculation of dry deposition velocities
164	!********************************************
165
166	if (DRYDEP) then
167	z0(4)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
168	z0(9)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
169
170	! Calculate relative humidity at surface
171	!***************************************
172	rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
173
174	call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
175	tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
176	ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
177	convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
178
179	do i=1,nspec
180	vdepn(ix,jy,i,n,l)=vd(i)
181	end do
182
183	endif
184
185	!******************************************************
186	! Calculate height of thermal tropopause (Hoinka, 1997)
187	!******************************************************
188
189	! 1) Calculate altitudes of ECMWF model levels
190	!*********************************************
191
192	tvold=tt2n(ix,jy,1,n,l)(1.+0.378ew(td2n(ix,jy,1,n,l))/ &
193	psn(ix,jy,1,n,l))
194	pold=psn(ix,jy,1,n,l)
195	zold=0.
196	do kz=2,nuvz
197	pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) ! pressure on model layers
198	tv=tthn(ix,jy,kz,n,l)(1.+0.608qvhn(ix,jy,kz,n,l))
199
200	if (abs(tv-tvold).gt.0.2) then
201	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
202	else
203	zlev(kz)=zold+constlog(pold/pint)tv
204	endif
205	tvold=tv
206	pold=pint
207	zold=zlev(kz)
208	end do
209
210	! 2) Define a minimum level kzmin, from which upward the tropopause is
211	! searched for. This is to avoid inversions in the lower troposphere
212	! to be identified as the tropopause
213	!************************************************************************
214
215	do kz=1,nuvz
216	if (zlev(kz).ge.altmin) then
217	kzmin=kz
218	goto 45
219	endif
220	end do
221	45 continue
222
223	! 3) Search for first stable layer above minimum height that fulfills the
224	! thermal tropopause criterion
225	!************************************************************************
226
227	do kz=kzmin,nuvz
228	do lz=kz+1,nuvz
229	if ((zlev(lz)-zlev(kz)).gt.2000.) then
230	if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
231	(zlev(lz)-zlev(kz))).lt.0.002) then
232	tropopausen(ix,jy,1,n,l)=zlev(kz)
233	goto 51
234	endif
235	goto 50
236	endif
237	end do
238	50 continue
239	end do
240	51 continue
241
242	end do
243	end do latloop
244	!$OMP END DO
245	!$OMP END PARALLEL
246
247
248	call calcpv_nests(l,n,uuhn,vvhn,pvhn)
249
250	end do nestloop
251
252
253	end subroutine calcpar_nests

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

Download in other formats:

Original Format