Context Navigation

calcpar_nests.f90

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

Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago
add SPDX-License-Identifier to all .f90 files
Property mode set to `100644`
File size: 8.6 KB

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

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

Download in other formats:

Original Format