Context Navigation

calcpar_nests.f90 @ b398fb6

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair

Last change on this file since b398fb6 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

Property mode set to 100644

File size: 13.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,metdata_format)
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	! Changes Arnold, D. and Morton, D. (2015): *
45	! -- description of local and common variables *
46	!*****************************************************************************
47	! Functions: *
48	! scalev computation of ustar *
49	! obukhov computatio of Obukhov length *
50	! *
51	!*****************************************************************************
52
53	use par_mod
54	use com_mod
55
56	implicit none
57	!***********************************************************************
58	! Subroutine Parameters: *
59	! input *
60	! n temporal index for meteorological fields (1 to 3)*
61	! uuhn,vvhn, wwhn wind components in ECMWF model levels *
62	! metdata_format to identify the met data
63	integer :: n, metdata_format
64	real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
65	real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
66	real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
67	!***********************************************************************
68	! Local variables *
69	! *
70	! ttlev
71	! qvlev
72	! * obukhov_ecmwf subroutine/function to calculate Obukhov length *
73	! * scalev subroutine/function to calculate ustar *
74	! ol obukhov length
75	! hmixplus maximum lifting from availiable kinetic enrgy *
76	! ulev, vlev wind speed at model levels *
77	! ew subroutine/function to calculate saturation *
78	! water vaport for a given air temperature *
79	! rh relative humidity at surface *
80	! vd deposition velocity from all species *
81	! subsceff excess hmin due to subgrid effects *
82	! ylat temporary latitude *
83	! altmin minimum height of the tropopause *
84	! tvoldm pold, zold temporary variables to keep previous values *
85	! pint pressure on model levels *
86	! tv virtual temperature on model levels *
87	! zlev height of model levels *
88	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov_ecmwf,obukhov_gfs,scalev,ol,hmixplus
89	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
90	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
91
92	! Other variables:
93	! ix,jy,kz,i,lz,kzmin loop control indices in each direction *
94	integer :: ix,jy,i,l,kz,lz,kzmin
95	!***********************************************************************
96
97	!***********************************************************************
98	! Local constants *
99	real,parameter :: const=r_air/ga
100	!***********************************************************************
101
102	!***********************************************************************
103	! Global variables *
104	! from par_mod and com_mod *
105	! olin [m] inverse Obukhov length (1/L) *
106	! hmixn [m] mixing height *
107	! wstarn [m/s] convective velocity scale *
108	! ustarn [m/s] friction velocity *
109	! z0n roughness length for the landuse classes *
110	! tropopausen [m] altitude of thermal tropopause *
111	! nxn, nyn actual dimensions of nested wind fields in x and y direction*
112	! dxn, dyn grid distances in x,y direction for the nested grids *
113	! akm, bkm coefficients which regulate vertical discretization of ecmwf*
114	! akz, bkz model discretization coeffizients at the centre of layers *
115	! psn surface pressure for nests *
116	! tt2n 2-m temperature for nests *
117	! tt2dn 2-m dew temperature for nests *
118	! sshfn surface sensible heat flux for nests *
119	! surfstrn surface stress fro nests *
120	! numbnests number of nested grids *
121	! convprecn convective precip on the nested grid *
122	! lsprecn large scale precip on the nested grid *
123	! sdn snow depth on the nested grid *
124	! ssrn surface solar radiation for nests *
125	! xlon0n, ylat0n geographical longitude/latitude of lower left grid *
126	! point of nested wind fields *
127	!
128	!************************************************************************
129
130	!-----------------------------------------------------------------------------
131
132
133	! Loop over all nests
134	!********************
135
136	do l=1,numbnests
137
138	! Loop over entire grid
139	!**********************
140
141	do jy=0,nyn(l)-1
142
143	! Set minimum height for tropopause
144	!**********************************
145
146	ylat=ylat0n(l)+real(jy)*dyn(l)
147	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
148	altmin = 5000.
149	else
150	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
151	altmin=2500.+(40.-ylat)*125.
152	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
153	altmin=2500.+(40.+ylat)*125.
154	else
155	altmin=2500.
156	endif
157	endif
158
159	do ix=0,nxn(l)-1
160
161	! 1) Calculation of friction velocity
162	!************************************
163
164	ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
165	td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
166
167	! 2) Calculation of inverse Obukhov length scale
168	!***********************************************
169
170	if (metdata_format.eq.ecmwf_metdata) then
171	ol=obukhov_ecmwf(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
172	td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
173	sshfn(ix,jy,1,n,l),akm,bkm)
174	endif
175	if (metdata_format.eq.gfs_metdata) then
176	ol=obukhov_gfs(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
177	td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
178	sshfn(ix,jy,1,n,l),akm,bkm)
179	endif
180	if (ol.ne.0.) then
181	olin(ix,jy,1,n,l)=1./ol
182	else
183	olin(ix,jy,1,n,l)=99999.
184	endif
185
186
187	! 3) Calculation of convective velocity scale and mixing height
188	!**************************************************************
189
190	do i=1,nuvz
191	ulev(i)=uuhn(ix,jy,i,l)
192	vlev(i)=vvhn(ix,jy,i,l)
193	ttlev(i)=tthn(ix,jy,i,n,l)
194	qvlev(i)=qvhn(ix,jy,i,n,l)
195	end do
196
197	if ( metdata_format.eq.ecmwf_metdata) then
198	call richardson_ecmwf(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
199	qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
200	tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
201	wstarn(ix,jy,1,n,l),hmixplus)
202	endif
203	if ( metdata_format.eq.gfs_metdata) then
204	call richardson_gfs(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
205	qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
206	tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
207	wstarn(ix,jy,1,n,l),hmixplus)
208	endif
209
210
211	if(lsubgrid.eq.1) then
212	subsceff=min(excessoron(ix,jy,l),hmixplus)
213	else
214	subsceff=0
215	endif
216	!
217	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
218	!
219	hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
220	hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
221	hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
222
223
224	! 4) Calculation of dry deposition velocities
225	!********************************************
226
227	if (DRYDEP) then
228	z0(4)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
229	z0(9)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
230
231	! Calculate relative humidity at surface
232	!***************************************
233	rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
234
235	call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
236	tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
237	ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
238	convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
239
240	do i=1,nspec
241	vdepn(ix,jy,i,n,l)=vd(i)
242	end do
243
244	endif
245
246	!******************************************************
247	! Calculate height of thermal tropopause (Hoinka, 1997)
248	!******************************************************
249
250	! 1) Calculate altitudes of ECMWF model levels
251	!*********************************************
252
253	tvold=tt2n(ix,jy,1,n,l)(1.+0.378ew(td2n(ix,jy,1,n,l))/ &
254	psn(ix,jy,1,n,l))
255	pold=psn(ix,jy,1,n,l)
256	zold=0.
257	do kz=2,nuvz
258	pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) ! pressure on model layers
259	tv=tthn(ix,jy,kz,n,l)(1.+0.608qvhn(ix,jy,kz,n,l))
260
261	if (abs(tv-tvold).gt.0.2) then
262	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
263	else
264	zlev(kz)=zold+constlog(pold/pint)tv
265	endif
266	tvold=tv
267	pold=pint
268	zold=zlev(kz)
269	end do
270
271	! 2) Define a minimum level kzmin, from which upward the tropopause is
272	! searched for. This is to avoid inversions in the lower troposphere
273	! to be identified as the tropopause
274	!************************************************************************
275
276	do kz=1,nuvz
277	if (zlev(kz).ge.altmin) then
278	kzmin=kz
279	goto 45
280	endif
281	end do
282	45 continue
283
284	! 3) Search for first stable layer above minimum height that fulfills the
285	! thermal tropopause criterion
286	!************************************************************************
287
288	do kz=kzmin,nuvz
289	do lz=kz+1,nuvz
290	if ((zlev(lz)-zlev(kz)).gt.2000.) then
291	if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
292	(zlev(lz)-zlev(kz))).lt.0.002) then
293	tropopausen(ix,jy,1,n,l)=zlev(kz)
294	goto 51
295	endif
296	goto 50
297	endif
298	end do
299	50 continue
300	end do
301	51 continue
302
303
304	end do
305	end do
306
307
308	call calcpv_nests(l,n,uuhn,vvhn,pvhn)
309
310	end do
311
312
313	end subroutine calcpar_nests

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

Download in other formats:

Original Format