Context Navigation

calcpar_nests.f90 @ 3481cc1

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

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 8.5 KB

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

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

Download in other formats:

Original Format