Context Navigation

calcpar_nests.f90 @ c0884a8

univie

Last change on this file since c0884a8 was c0884a8, checked in by pesei <petra seibert at univie ac at>, 6 years ago
replace CTBTO code for checking type of GRIB
Property mode set to `100644`
File size: 10.1 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,id_centre)
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	! Unified ECMWF and GFS builds *
45	! Marian Harustak, 12.5.2017 *
46	! - Added passing of id_centre as it was needed by called routines *
47	! *
48	! Petra Seibert, 2018-06-26: simplified version met data format detection *
49	! *
50	!*****************************************************************************
51	! *
52	! Variables: *
53	! n temporal index for meteorological fields (1 to 3) *
54	! id_centre format of metdata (ecmwf/gfs) *
55	! *
56	! Constants: *
57	! *
58	! *
59	! Functions: *
60	! scalev computation of ustar *
61	! obukhov computatio of Obukhov length *
62	! *
63	!*****************************************************************************
64
65	use par_mod
66	use com_mod
67
68	implicit none
69
70	integer :: id_centre
71	integer :: n,ix,jy,i,l,kz,lz,kzmin
72	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
73	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
74	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
75	real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
76	real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
77	real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
78	real,parameter :: const=r_air/ga
79
80
81	! Loop over all nests
82	!********************
83
84	do l=1,numbnests
85
86	! Loop over entire grid
87	!**********************
88
89	do jy=0,nyn(l)-1
90
91	! Set minimum height for tropopause
92	!**********************************
93
94	ylat=ylat0n(l)+real(jy)*dyn(l)
95	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
96	altmin = 5000.
97	else
98	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
99	altmin=2500.+(40.-ylat)*125.
100	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
101	altmin=2500.+(40.+ylat)*125.
102	else
103	altmin=2500.
104	endif
105	endif
106
107	do ix=0,nxn(l)-1
108
109	! 1) Calculation of friction velocity
110	!************************************
111
112	ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
113	td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
114	if (ustarn(ix,jy,1,n,l).le.1.e-8) ustarn(ix,jy,1,n,l)=1.e-8
115
116	! 2) Calculation of inverse Obukhov length scale
117	!***********************************************
118
119	ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
120	td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
121	sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,id_centre)
122	if (ol.ne.0.) then
123	olin(ix,jy,1,n,l)=1./ol
124	else
125	olin(ix,jy,1,n,l)=99999.
126	endif
127
128
129	! 3) Calculation of convective velocity scale and mixing height
130	!**************************************************************
131
132	do i=1,nuvz
133	ulev(i)=uuhn(ix,jy,i,l)
134	vlev(i)=vvhn(ix,jy,i,l)
135	ttlev(i)=tthn(ix,jy,i,n,l)
136	qvlev(i)=qvhn(ix,jy,i,n,l)
137	end do
138
139	call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
140	qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
141	tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
142	wstarn(ix,jy,1,n,l),hmixplus,id_centre)
143
144	if(lsubgrid.eq.1) then
145	subsceff=min(excessoron(ix,jy,l),hmixplus)
146	else
147	subsceff=0.0
148	endif
149	!
150	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
151	!
152	hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
153	hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
154	hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
155
156
157	! 4) Calculation of dry deposition velocities
158	!********************************************
159
160	if (DRYDEP) then
161	! z0(4)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
162	! z0(9)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
163	z0(7)=0.016ustarn(ix,jy,1,n,l)ustarn(ix,jy,1,n,l)/ga
164
165	! Calculate relative humidity at surface
166	!***************************************
167	rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
168
169	call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
170	tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
171	ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
172	convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
173
174	do i=1,nspec
175	vdepn(ix,jy,i,n,l)=vd(i)
176	end do
177
178	endif
179
180	!******************************************************
181	! Calculate height of thermal tropopause (Hoinka, 1997)
182	!******************************************************
183
184	! 1) Calculate altitudes of ECMWF model levels
185	!*********************************************
186
187	tvold=tt2n(ix,jy,1,n,l)(1.+0.378ew(td2n(ix,jy,1,n,l))/ &
188	psn(ix,jy,1,n,l))
189	pold=psn(ix,jy,1,n,l)
190	zold=0.
191	do kz=2,nuvz
192	pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) ! pressure on model layers
193	tv=tthn(ix,jy,kz,n,l)(1.+0.608qvhn(ix,jy,kz,n,l))
194
195	if (abs(tv-tvold).gt.0.2) then
196	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
197	else
198	zlev(kz)=zold+constlog(pold/pint)tv
199	endif
200	tvold=tv
201	pold=pint
202	zold=zlev(kz)
203	end do
204
205	! 2) Define a minimum level kzmin, from which upward the tropopause is
206	! searched for. This is to avoid inversions in the lower troposphere
207	! to be identified as the tropopause
208	!************************************************************************
209
210	do kz=1,nuvz
211	if (zlev(kz).ge.altmin) then
212	kzmin=kz
213	goto 45
214	endif
215	end do
216	45 continue
217
218	! 3) Search for first stable layer above minimum height that fulfills the
219	! thermal tropopause criterion
220	!************************************************************************
221
222	do kz=kzmin,nuvz
223	do lz=kz+1,nuvz
224	if ((zlev(lz)-zlev(kz)).gt.2000.) then
225	if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
226	(zlev(lz)-zlev(kz))).lt.0.002) then
227	tropopausen(ix,jy,1,n,l)=zlev(kz)
228	goto 51
229	endif
230	goto 50
231	endif
232	end do
233	50 continue
234	end do
235	51 continue
236
237
238	end do
239	end do
240
241	! Calculation of potential vorticity on 3-d grid
242	!***********************************************
243
244	call calcpv_nests(l,n,uuhn,vvhn,pvhn)
245
246	end do
247
248
249	end subroutine calcpar_nests

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

Download in other formats:

Original Format