Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_temfsfclay.F @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 19.1 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3	MODULE module_sf_temfsfclay
4
5	CONTAINS
6
7	!-------------------------------------------------------------------
8	SUBROUTINE temfsfclay(u3d,v3d,th3d,qv3d,p3d,pi3d,rho,z,ht, &
9	cp,g,rovcp,r,xlv,psfc,chs,chs2,cqs2,cpm, &
10	znt,ust,mavail,xland, &
11	hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc, &
12	u10,v10,th2,t2,q2, &
13	svp1,svp2,svp3,svpt0,ep1,ep2, &
14	karman,fCor,te_temf, &
15	hd_temf,exch_temf,wm_temf, &
16	ids,ide, jds,jde, kds,kde, &
17	ims,ime, jms,jme, kms,kme, &
18	its,ite, jts,jte, kts,kte &
19	)
20	!-------------------------------------------------------------------
21	IMPLICIT NONE
22	!-------------------------------------------------------------------
23	!
24	! This is the Total Energy - Mass Flux (TEMF) surface layer scheme.
25	! Initial implementation 2010 by Wayne Angevine, CIRES/NOAA ESRL.
26	! References:
27	! Angevine et al., 2010, MWR
28	! Angevine, 2005, JAM
29	! Mauritsen et al., 2007, JAS
30	!
31	!-------------------------------------------------------------------
32	!-------------------------------------------------------------------
33	!-- u3d 3D u-velocity interpolated to theta points (m/s)
34	!-- v3d 3D v-velocity interpolated to theta points (m/s)
35	!-- th3d potential temperature (K)
36	!-- qv3d 3D water vapor mixing ratio (Kg/Kg)
37	!-- p3d 3D pressure (Pa)
38	!-- cp heat capacity at constant pressure for dry air (J/kg/K)
39	!-- g acceleration due to gravity (m/s^2)
40	!-- rovcp R/CP
41	!-- r gas constant for dry air (J/kg/K)
42	!-- xlv latent heat of vaporization for water (J/kg)
43	!-- psfc surface pressure (Pa)
44	!-- chs heat/moisture exchange coefficient for LSM (m/s)
45	!-- chs2
46	!-- cqs2
47	!-- cpm
48	!-- znt roughness length (m)
49	!-- ust u* in similarity theory (m/s)
50	!-- mavail surface moisture availability (between 0 and 1)
51	!-- xland land mask (1 for land, 2 for water)
52	!-- hfx upward heat flux at the surface (W/m^2)
53	!-- qfx upward moisture flux at the surface (kg/m^2/s)
54	!-- lh net upward latent heat flux at surface (W/m^2)
55	!-- tsk surface temperature (K)
56	!-- flhc exchange coefficient for heat (W/m^2/K)
57	!-- flqc exchange coefficient for moisture (kg/m^2/s)
58	!-- qgh lowest-level saturated mixing ratio
59	!-- qsfc ground saturated mixing ratio
60	!-- u10 diagnostic 10m u wind
61	!-- v10 diagnostic 10m v wind
62	!-- th2 diagnostic 2m theta (K)
63	!-- t2 diagnostic 2m temperature (K)
64	!-- q2 diagnostic 2m mixing ratio (kg/kg)
65	!-- svp1 constant for saturation vapor pressure (kPa)
66	!-- svp2 constant for saturation vapor pressure (dimensionless)
67	!-- svp3 constant for saturation vapor pressure (K)
68	!-- svpt0 constant for saturation vapor pressure (K)
69	!-- ep1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
70	!-- ep2 constant for specific humidity calculation
71	! (R_d/R_v) (dimensionless)
72	!-- karman Von Karman constant
73	!-- fCor Coriolis parameter
74	!-- ids start index for i in domain
75	!-- ide end index for i in domain
76	!-- jds start index for j in domain
77	!-- jde end index for j in domain
78	!-- kds start index for k in domain
79	!-- kde end index for k in domain
80	!-- ims start index for i in memory
81	!-- ime end index for i in memory
82	!-- jms start index for j in memory
83	!-- jme end index for j in memory
84	!-- kms start index for k in memory
85	!-- kme end index for k in memory
86	!-- its start index for i in tile
87	!-- ite end index for i in tile
88	!-- jts start index for j in tile
89	!-- jte end index for j in tile
90	!-- kts start index for k in tile
91	!-- kte end index for k in tile
92	!-------------------------------------------------------------------
93	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
94	ims,ime, jms,jme, kms,kme, &
95	its,ite, jts,jte, kts,kte
96	!
97	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
98	INTENT(IN ) :: u3d, v3d, th3d, qv3d, p3d, pi3d, rho, z
99	REAL, DIMENSION( ims:ime, jms:jme ) , &
100	INTENT(IN ) :: mavail, xland, tsk, fCor, ht, psfc, znt
101	REAL, DIMENSION( ims:ime, jms:jme ) , &
102	INTENT(INOUT) :: hfx, qfx, lh, flhc, flqc
103	REAL, DIMENSION( ims:ime, jms:jme ) , &
104	INTENT(INOUT) :: ust, chs2, cqs2, chs, cpm, qgh, qsfc
105	REAL, DIMENSION( ims:ime, jms:jme ) , &
106	INTENT(OUT ) :: u10, v10, th2, t2, q2
107	REAL, DIMENSION( ims:ime, jms:jme ) , &
108	INTENT(IN ) :: hd_temf
109	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
110	INTENT(INOUT) :: te_temf
111	REAL, DIMENSION( ims:ime, jms:jme ) , &
112	INTENT( OUT) :: exch_temf
113	REAL, DIMENSION( ims:ime, jms:jme ) , &
114	INTENT(INOUT) :: wm_temf
115
116
117	REAL, INTENT(IN ) :: cp,g,rovcp,r,xlv
118	REAL, INTENT(IN ) :: svp1,svp2,svp3,svpt0
119	REAL, INTENT(IN ) :: ep1,ep2,karman
120	!
121	! LOCAL VARS
122
123	INTEGER :: J
124	!
125
126	DO J=jts,jte
127
128	CALL temfsfclay1d(j,u1d=u3d(ims,kms,j),v1d=v3d(ims,kms,j), &
129	th1d=th3d(ims,kms,j),qv1d=qv3d(ims,kms,j),p1d=p3d(ims,kms,j), &
130	pi1d=pi3d(ims,kms,j),rho=rho(ims,kms,j),z=z(ims,kms,j),&
131	zsrf=ht(ims,j), &
132	cp=cp,g=g,rovcp=rovcp,r=r,xlv=xlv,psfc=psfc(ims,j), &
133	chs=chs(ims,j),chs2=chs2(ims,j),cqs2=cqs2(ims,j), &
134	cpm=cpm(ims,j),znt=znt(ims,j),ust=ust(ims,j), &
135	mavail=mavail(ims,j),xland=xland(ims,j), &
136	hfx=hfx(ims,j),qfx=qfx(ims,j),lh=lh(ims,j),tsk=tsk(ims,j), &
137	flhc=flhc(ims,j),flqc=flqc(ims,j),qgh=qgh(ims,j), &
138	qsfc=qsfc(ims,j),u10=u10(ims,j),v10=v10(ims,j), &
139	th2=th2(ims,j),t2=t2(ims,j),q2=q2(ims,j), &
140	svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0, &
141	ep1=ep1,ep2=ep2,karman=karman,fCor=fCor(ims,j), &
142	te_temfx=te_temf(ims,kms,j),hd_temfx=hd_temf(ims,j), &
143	exch_temfx=exch_temf(ims,j),wm_temfx=wm_temf(ims,j), &
144	ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
145	ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
146	its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
147	)
148	ENDDO
149
150	END SUBROUTINE temfsfclay
151
152
153	!-------------------------------------------------------------------
154	SUBROUTINE temfsfclay1d(j,u1d,v1d,th1d,qv1d,p1d, &
155	pi1d,rho,z,zsrf,cp,g,rovcp,r,xlv,psfc, &
156	chs,chs2,cqs2,cpm,znt,ust, &
157	mavail,xland,hfx,qfx,lh,tsk, &
158	flhc,flqc,qgh,qsfc,u10,v10, &
159	th2,t2,q2,svp1,svp2,svp3,svpt0, &
160	ep1,ep2,karman,fCor, &
161	te_temfx,hd_temfx,exch_temfx,wm_temfx, &
162	ids,ide, jds,jde, kds,kde, &
163	ims,ime, jms,jme, kms,kme, &
164	its,ite, jts,jte, kts,kte &
165	)
166	!!-------------------------------------------------------------------
167	IMPLICIT NONE
168	!!-------------------------------------------------------------------
169	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
170	ims,ime, jms,jme, kms,kme, &
171	its,ite, jts,jte, kts,kte, &
172	j
173
174	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: &
175	u1d,v1d,qv1d,p1d,th1d,pi1d,rho,z,zsrf
176	REAL, INTENT(IN ) :: cp,g,rovcp,r,xlv
177	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: psfc,znt
178	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
179	chs,chs2,cqs2,cpm,ust
180	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: mavail,xland
181	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
182	hfx,qfx,lh
183	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: tsk
184	REAL, DIMENSION( ims:ime ), INTENT( OUT) :: &
185	flhc,flqc
186	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: &
187	qgh,qsfc
188	REAL, DIMENSION( ims:ime ), INTENT( OUT) :: &
189	u10,v10,th2,t2,q2
190	REAL, INTENT(IN ) :: svp1,svp2,svp3,svpt0
191	REAL, INTENT(IN ) :: ep1,ep2,karman
192	REAL, DIMENSION( ims:ime ), INTENT(IN ) :: fCor,hd_temfx
193	REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: te_temfx
194	REAL, DIMENSION( ims:ime ), INTENT( OUT) :: exch_temfx, wm_temfx
195	!
196	!! LOCAL VARS
197	! TE model constants
198	real, parameter :: visc_temf = 1.57e-5
199	real, parameter :: conduc_temf = 1.57e-5 / 0.733
200	logical, parameter :: MFopt = .true. ! Use mass flux or not
201	real, parameter :: TEmin = 1e-3
202	real, parameter :: ftau0 = 0.17
203	real, parameter :: fth0 = 0.145
204	! real, parameter :: fth0 = 0.12 ! WA 10/13/10 to make PrT0 ~= 1
205	real, parameter :: Cf = 0.185
206	real, parameter :: CN = 2.0
207	! real, parameter :: Ceps = ftau0**1.5
208	real, parameter :: Ceps = 0.070
209	real, parameter :: Cgamma = Ceps
210	real, parameter :: Cphi = Ceps
211	! real, parameter :: PrT0 = Cphi/Ceps * ftau02. / 2 / fth02.
212	real, parameter :: PrT0 = Cphi/Ceps * ftau02 / 2. / fth02
213	!
214	integer :: i
215	real :: e1
216	real, dimension( its:ite) :: wstr, ang, wm
217	real, dimension( its:ite) :: z0t
218	real, dimension( its:ite) :: dthdz, dqtdz, dudz, dvdz
219	real, dimension( its:ite) :: lepsmin
220	real, dimension( its:ite) :: thetav
221	real, dimension( its:ite) :: zt,zm
222	real, dimension( its:ite) :: N2, S, Ri, beta, ftau, fth, ratio
223	real, dimension( its:ite) :: TKE, TE2
224	real, dimension( its:ite) :: ustrtilde, linv, leps
225	real, dimension( its:ite) :: km, kh
226	real, dimension( its:ite) :: qsfc_air
227	!!-------------------------------------------------------------------
228
229	!!!!!!! ******
230	! WA Known outages: None
231
232	do i = its,ite ! Main loop
233
234	! Calculate surface saturated q and q in air at surface
235	e1=svp1exp(svp2(tsk(i)-svpt0)/(tsk(i)-svp3))
236	qsfc(i)=ep2*e1/((psfc(i)/1000.)-e1)
237	qsfc_air(i) = qsfc(i) * mavail(i)
238	thetav(i) = (tsk(i)/pi1d(i)) * (1. + 0.608*qsfc_air(i)) ! WA Assumes ql(env)=0, what if it isn't?
239	! WA TEST (R5) set z0t = z0
240	! z0t(i) = znt(i) / 10.0 ! WA this is hard coded in Matlab version
241	z0t(i) = znt(i)
242
243	! Get height and delta at turbulence levels and mass levels
244	zt(i) = (z(i) - zsrf(i) - znt(i)) / 2.
245	zm(i) = z(i) - zsrf(i)
246
247	! Gradients at first level
248	dthdz(i) = (th1d(i)-(tsk(i)/pi1d(i))) / (zt(i) * log10(zm(i)/z0t(i)))
249	dqtdz(i) = (qv1d(i)-qsfc_air(i)) / (zt(i) * log10(zm(i)/z0t(i)))
250	dudz(i) = u1d(i) / (zt(i) * log10(zm(i)/znt(i)))
251	dvdz(i) = v1d(i) / (zt(i) * log10(zm(i)/znt(i)))
252
253	! WA doing this because te_temf may not be initialized,
254	! would be better to do it in initialization routine but it's
255	! not available in module_physics_init.
256	if (te_temfx(i) < TEmin) te_temfx(i) = TEmin
257
258	if ( hfx(i) > 0.) then
259	wstr(i) = (g * hd_temfx(i) / thetav(i) * (hfx(i)/(rho(i)cp))) * (1./3.)
260	else
261	wstr(i) = 0.
262	end if
263
264	! Find stability parameters and length scale
265	! WA Calculation of N should really use d(thetaV)/dz not dthdz
266	! WA 7/1/09 allow N to be negative
267	! if ( dthdz(i) >= 0.) then
268	! N(i) = csqrt(g / thetav(i) * dthdz(i))
269	! else
270	! N(i) = 0.
271	! end if
272	N2(i) = g / thetav(i) * dthdz(i)
273	S(i) = sqrt(dudz(i)2. + dvdz(i)2.)
274	! Ri(i) = N(i)2. / S(i)2.
275	Ri(i) = N2(i) / S(i)**2.
276	! if (S(i) < 1e-15) Ri(i) = 1./1e-15
277	if (S(i) < 1e-15) then
278	print *,'In TEMF SFC Limiting Ri,S,N2,Ri,u,v = ',S(i),N2(i),Ri(i),u1d(i),v1d(i)
279	if (N2(i) >= 0) then
280	Ri(i) = 0.2
281	else
282	Ri(i) = -1.
283	end if
284	end if
285	if (Ri(i) > 0.2) then ! WA TEST to prevent runaway
286	Ri(i) = 0.2
287	end if
288	beta(i) = g / thetav(i)
289	! WA 7/1/09 adjust ratio, ftau, fth for Ri>0
290	if (Ri(i) > 0) then
291	ratio(i) = Ri(i)/(Cphi*2.ftau0*2./(2.Ceps*2.fth0*2.)+3.Ri(i))
292	ftau(i) = ftau0 * ((3./4.) / (1.+4.*Ri(i)) + 1./4.)
293	fth(i) = fth0 / (1.+4.*Ri(i))
294	! TE2(i) = 2. * te_temfx(i) * ratio(i) * N(i)2. / beta(i)2.
295	TE2(i) = 2. * te_temfx(i) * ratio(i) * N2(i) / beta(i)**2.
296	else
297	ratio(i) = Ri(i)/(Cphi*2.ftau0*2./(-2.Ceps*2.fth0*2.)+2.Ri(i))
298	ftau(i) = ftau0
299	fth(i) = fth0
300	TE2(i) = 0.
301	end if
302	TKE(i) = te_temfx(i) * (1. - ratio(i))
303	ustrtilde(i) = sqrt(ftau(i) * TKE(i))
304	! linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cfustrtilde(i)) + N(i)/(CNustrtilde(i))
305	if (N2(i) > 0.) then
306	linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cfustrtilde(i)) + sqrt(N2(i))/(CNustrtilde(i))
307	else
308	linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cf*ustrtilde(i))
309	end if
310	leps(i) = 1./linv(i)
311	! WA TEST (R4) remove lower limit on leps
312	! lepsmin(i) = min(0.4*zt(i), 5.)
313	lepsmin(i) = 0.
314	leps(i) = max(leps(i),lepsmin(i))
315
316
317	! Find diffusion coefficients
318	! First use basic formulae for stable and neutral cases,
319	! then for convective conditions, and finally choose the larger
320	km(i) = TKE(i)*1.5 ftau(i)*2. / (-beta(i) fth(i) * sqrt(TE2(i)) + Ceps * sqrt(TKE(i)*te_temfx(i)) / leps(i))
321	kh(i) = 2. * leps(i) * fth(i)*2. TKE(i) / sqrt(te_temfx(i)) / Cphi
322	km(i) = max(km(i),visc_temf)
323	kh(i) = max(kh(i),conduc_temf)
324
325	! Surface fluxes
326	ust(i) = sqrt(ftau(i)/ftau0) * sqrt(u1d(i)2. + v1d(i)2.) * leps(i) / log(zm(i)/znt(i)) / zt(i)
327	ang(i) = atan2(v1d(i),u1d(i))
328
329	! Calculate mixed scaling velocity (Moeng & Sullivan 1994 JAS p.1021)
330	! Replaces ust everywhere (WA need to reconsider?)
331	! WA wm is too large, makes surface flux too big and cools sfc too much
332	! wm(i) = (1./5. * (wstr(i)*3. + 5. ust(i)3.)) (1./3.)
333	! WA TEST (R2,R11) 7/23/10 reduce velocity scale to fix excessive fluxes
334	wm(i) = 0.5 * (1./5. * (wstr(i)*3. + 5. ust(i)3.)) (1./3.)
335	! WA TEST 2/14/11 limit contribution of w*
336	! wm(i) = 0.5 * (1./5. * (min(0.8,wstr(i))*3. + 5. ust(i)3.)) (1./3.)
337	! WA TEST 2/22/11 average with previous value to reduce instability
338	wm(i) = (wm(i) + wm_temfx(i)) / 2.0
339	wm_temfx(i) = wm(i)
340	! WA TEST (R3-R10) 7/23/10 wm = u*
341	! wm(i) = ust(i)
342
343	! Populate surface exchange coefficient variables to go back out
344	! for next time step of surface scheme
345	! Unit specifications in SLAB and sfclay are conflicting and probably
346	! incorrect. This will give a dynamic heat flux (W/m^2) or moisture
347	! flux (kg(water)/(m^2*s)) when multiplied by a difference.
348	! These formulae are the same as what's used above to get surface
349	! flux from surface temperature and specific humidity.
350	flhc(i) = rho(i) * cp * fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(zm(i)/z0t(i)) / zt(i)
351	flqc(i) = rho(i) * fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(zm(i)/z0t(i)) / zt(i) * mavail(i)
352	exch_temfx(i) = flqc(i) / mavail(i)
353	chs(i) = flqc(i) / rho(i) / mavail(i)
354	! WA Must exchange coeffs be limited to avoid runaway in some
355	! (convective?) conditions? Something like this is done in sfclay.
356	! Doing nothing for now.
357
358	! Populate surface heat and moisture fluxes
359	hfx(i) = flhc(i) * (tsk(i) - th1d(i)*pi1d(i))
360	! qfx(i) = flqc(i) * (qsfc_air(i) - qv1d(i)) ! WA 2/16/11
361	qfx(i) = flqc(i) * (qsfc(i) - qv1d(i))
362	qfx(i) = max(qfx(i),0.) ! WA this is done in sfclay, is it right?
363	lh(i)=xlv*qfx(i)
364
365
366	! Populate 10 m winds and 2 m temp and 2 m exchange coeffs
367	! WA Note this only works if first mass level is above 10 m
368	u10(i) = u1d(i) * log(10.0/znt(i)) / log(zm(i)/znt(i))
369	v10(i) = v1d(i) * log(10.0/znt(i)) / log(zm(i)/znt(i))
370	t2(i) = (tsk(i)/pi1d(i) + (th1d(i) - tsk(i)/pi1d(i)) * log(2.0/z0t(i)) / log(zm(i)/z0t(i))) * pi1d(i) ! WA this should also use pi at z0
371	th2(i) = t2(i) / pi1d(i)
372	q2(i) = (qsfc_air(i) + (qv1d(i) - qsfc_air(i)) * log(2.0/znt(i)) / log(zm(i)/znt(i)))
373	! WA are these correct? Difference between chs2 and cqs2 is unclear
374	! At the moment the only difference is z0t vs. znt
375	chs2(i) = fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(2.0/z0t(i)) / zt(i)
376	cqs2(i) = fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(2.0/znt(i)) / zt(i)
377
378	! Calculate qgh (saturated at first-level temp) and cpm
379	e1=svp1exp(svp2((th1d(i)pi1d(i))-svpt0)/((th1d(i)pi1d(i))-svp3))
380	qgh(i)=ep2*e1/((p1d(i)/1000.)-e1)
381	cpm(i)=cp(1.+0.8qv1d(i))
382
383	end do ! Main loop
384
385	END SUBROUTINE temfsfclay1d
386
387	!====================================================================
388	SUBROUTINE temfsfclayinit( restart, allowed_to_read, &
389	wm_temf, &
390	ids, ide, jds, jde, kds, kde, &
391	ims, ime, jms, jme, kms, kme, &
392	its, ite, jts, jte, kts, kte )
393
394	logical , intent(in) :: restart, allowed_to_read
395	REAL, DIMENSION( ims:ime, jms:jme ) , &
396	INTENT( OUT) :: wm_temf
397	integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
398	ims, ime, jms, jme, kms, kme, &
399	its, ite, jts, jte, kts, kte
400
401	! Local variables
402	integer :: i, j, itf, jtf
403	!
404	CALL wrf_debug( 100, 'in temfsfclayinit' )
405	jtf = min0(jte,jde-1)
406	itf = min0(ite,ide-1)
407	!
408	if(.not.restart)then
409	do j = jts,jtf
410	do i = its,itf
411	! do j = jms,jme
412	! do i = ims,ime
413	wm_temf(i,j) = 0.0
414	enddo
415	enddo
416	endif
417
418	END SUBROUTINE temfsfclayinit
419
420	!-------------------------------------------------------------------
421
422	END MODULE module_sf_temfsfclay

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

Download in other formats:

Original Format