Context Navigation

calcpar_nests.f90 @ 02095e3

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

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

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

Download in other formats:

Original Format