Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_temf.F @ 1

Last change on this file since 1 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: 64.8 KB

Line
1	!wrf:model_layer:physics
2	!
3	!
4	!
5	!
6	module module_bl_temf
7	contains
8	!
9	!-------------------------------------------------------------------
10	!
11	subroutine temfpbl(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d,rho, &
12	rublten,rvblten,rthblten, &
13	rqvblten,rqcblten,rqiblten,flag_qi, &
14	g,cp,rcp,r_d,r_v,cpv, &
15	z,xlv,psfc, &
16	mut,p_top, &
17	znt,ht,ust,zol,hol,hpbl,psim,psih, &
18	xland,hfx,qfx,tsk,qsfc,gz1oz0,wspd,br, &
19	dt,dtmin,kpbl2d, &
20	svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt, &
21	kh_temf,km_temf, &
22	u10,v10,t2, &
23	te_temf,shf_temf,qf_temf,uw_temf,vw_temf, &
24	wupd_temf,mf_temf,thup_temf,qtup_temf,qlup_temf, &
25	cf3d_temf,cfm_temf, &
26	hd_temf,lcl_temf,hct_temf, &
27	flhc,flqc,exch_temf, &
28	fCor, &
29	ids,ide, jds,jde, kds,kde, &
30	ims,ime, jms,jme, kms,kme, &
31	its,ite, jts,jte, kts,kte &
32	)
33	!-------------------------------------------------------------------
34	implicit none
35	!-------------------------------------------------------------------
36	! New variables for TEMF
37	!-- te_temf Total energy from this scheme
38	!-- shf_temf Sensible heat flux profile from this scheme (kinematic)
39	!-- qf_temf Moisture flux profile from this scheme (kinematic)
40	!-- uw_temf U momentum flux component from this scheme
41	!-- vw_temf V momentum flux component from this scheme
42	!-- kh_temf Exchange coefficient for heat (3D)
43	!-- km_temf Exchange coefficient for momentum (3D)
44	!-- wupd_temf Updraft velocity from TEMF BL scheme
45	!-- mf_temf Mass flux from TEMF BL scheme
46	!-- thup_temf Updraft thetal from TEMF BL scheme
47	!-- qtup_temf Updraft qt from TEMF BL scheme
48	!-- qlup_temf Updraft ql from TEMF BL scheme
49	!-- cf3d_temf 3D cloud fraction from TEMF BL scheme
50	!-- cfm_temf Column cloud fraction from TEMF BL scheme
51	!-- exch_temf Surface exchange coefficient (as for moisture) from TEMF surface layer scheme
52	!-- flhc Surface exchange coefficient for heat (needed by surface scheme)
53	!-- flqc Surface exchange coefficient for moisture (including moisture availablity)
54	!-- fCor Coriolis parameter (from grid%f)
55	!
56	!-- u3d 3d u-velocity interpolated to theta points (m/s)
57	!-- v3d 3d v-velocity interpolated to theta points (m/s)
58	!-- th3d 3d potential temperature (k)
59	!-- t3d temperature (k)
60	!-- qv3d 3d water vapor mixing ratio (kg/kg)
61	!-- qc3d 3d cloud mixing ratio (kg/kg)
62	!-- qi3d 3d ice mixing ratio (kg/kg)
63	! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
64	!-- p3d 3d pressure (pa)
65	!-- p3di 3d pressure (pa) at interface level
66	!-- pi3d 3d exner function (dimensionless)
67	!-- rho 3d dry air density (kg/m^3)
68	!-- rublten u tendency due to
69	! pbl parameterization (m/s/s)
70	!-- rvblten v tendency due to
71	! pbl parameterization (m/s/s)
72	!-- rthblten theta tendency due to
73	! pbl parameterization (K/s)
74	!-- rqvblten qv tendency due to
75	! pbl parameterization (kg/kg/s)
76	!-- rqcblten qc tendency due to
77	! pbl parameterization (kg/kg/s)
78	!-- rqiblten qi tendency due to
79	! pbl parameterization (kg/kg/s)
80	!-- cp heat capacity at constant pressure for dry air (j/kg/k)
81	!-- g acceleration due to gravity (m/s^2)
82	!-- rovcp r/cp
83	!-- r_d gas constant for dry air (j/kg/k)
84	!-- rovg r/g
85	!-- z height above sea level (m)
86	!-- xlv latent heat of vaporization (j/kg)
87	!-- r_v gas constant for water vapor (j/kg/k)
88	!-- psfc pressure at the surface (pa)
89	!-- znt roughness length (m)
90	!-- ht terrain height ASL (m)
91	!-- ust u* in similarity theory (m/s)
92	!-- zol z/l height over monin-obukhov length
93	!-- hol pbl height over monin-obukhov length
94	!-- hpbl pbl height (m)
95	!-- psim similarity stability function for momentum
96	!-- psih similarity stability function for heat
97	!-- xland land mask (1 for land, 2 for water)
98	!-- hfx upward heat flux at the surface (w/m^2)
99	!-- qfx upward moisture flux at the surface (kg/m^2/s)
100	!-- tsk surface temperature (k)
101	!-- qsfc surface specific humidity (kg/kg)
102	!-- gz1oz0 log(z/z0) where z0 is roughness length
103	!-- wspd wind speed at lowest model level (m/s)
104	!-- u10 u-wind speed at 10 m (m/s)
105	!-- v10 v-wind speed at 10 m (m/s)
106	!-- br bulk richardson number in surface layer
107	!-- dt time step (s)
108	!-- dtmin time step (minute)
109	!-- rvovrd r_v divided by r_d (dimensionless)
110	!-- svp1 constant for saturation vapor pressure (kpa)
111	!-- svp2 constant for saturation vapor pressure (dimensionless)
112	!-- svp3 constant for saturation vapor pressure (k)
113	!-- svpt0 constant for saturation vapor pressure (k)
114	!-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
115	!-- ep2 constant for specific humidity calculation
116	!-- karman von karman constant
117	!-- eomeg angular velocity of earths rotation (rad/s)
118	!-- stbolt stefan-boltzmann constant (w/m^2/k^4)
119	!-- ids start index for i in domain
120	!-- ide end index for i in domain
121	!-- jds start index for j in domain
122	!-- jde end index for j in domain
123	!-- kds start index for k in domain
124	!-- kde end index for k in domain
125	!-- ims start index for i in memory
126	!-- ime end index for i in memory
127	!-- jms start index for j in memory
128	!-- jme end index for j in memory
129	!-- kms start index for k in memory
130	!-- kme end index for k in memory
131	!-- its start index for i in tile
132	!-- ite end index for i in tile
133	!-- jts start index for j in tile
134	!-- jte end index for j in tile
135	!-- kts start index for k in tile
136	!-- kte end index for k in tile
137	!-------------------------------------------------------------------
138	! Arguments
139	!
140	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
141	ims,ime, jms,jme, kms,kme, &
142	its,ite, jts,jte, kts,kte
143	!
144	real, intent(in ) :: dt,dtmin,g,cp,rcp,r_d,r_v,xlv,cpv
145	!
146	real, intent(in ) :: svp1,svp2,svp3,svpt0
147	real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
148	!
149	real, dimension( ims:ime, kms:kme, jms:jme ) , &
150	intent(in ) :: qv3d, qc3d, qi3d, &
151	p3d, pi3d, th3d, t3d, &
152	z, rho
153	!
154	real, dimension( ims:ime, kms:kme, jms:jme ) , &
155	intent(inout) :: te_temf
156	real, dimension( ims:ime, kms:kme, jms:jme ) , &
157	intent( out) :: shf_temf, qf_temf, uw_temf, vw_temf , &
158	wupd_temf, mf_temf, thup_temf, qtup_temf, &
159	qlup_temf,cf3d_temf
160	real, dimension( ims:ime, jms:jme ) , &
161	intent(inout) :: flhc, flqc, exch_temf
162	real, dimension( ims:ime, jms:jme ) , &
163	intent(in ) :: fCor
164	real, dimension( ims:ime, jms:jme ) , &
165	intent( out) :: hd_temf, lcl_temf, hct_temf, cfm_temf
166	!
167	real, dimension( ims:ime, kms:kme, jms:jme ) , &
168	intent(in ) :: p3di
169	!
170	real, dimension( ims:ime, kms:kme, jms:jme ) , &
171	intent(inout) :: rublten, rvblten, &
172	rthblten, &
173	rqvblten, rqcblten, rqiblten
174	!
175	real, dimension( ims:ime, kms:kme, jms:jme ) , &
176	intent(inout) :: kh_temf, km_temf
177	real, dimension( ims:ime, jms:jme ) , &
178	intent(inout) :: u10, v10, t2
179	!
180	real, dimension( ims:ime, jms:jme ) , &
181	intent(in ) :: xland, &
182	psim, psih, gz1oz0, br, &
183	psfc, tsk, qsfc
184	!
185	real, dimension( ims:ime, jms:jme ) , &
186	intent(inout) :: hfx, qfx
187	real, dimension( ims:ime, jms:jme ) , &
188	intent(inout) :: hol, ust, hpbl, znt, wspd, zol
189	real, dimension( ims:ime, jms:jme ) , &
190	intent(in ) :: ht
191	!
192	real, dimension( ims:ime, kms:kme, jms:jme ) , &
193	intent(in ) :: u3d, v3d
194	!
195	integer, dimension( ims:ime, jms:jme ) , &
196	intent(out ) :: kpbl2d
197	!
198	logical, intent(in) :: flag_qi
199	!
200	! real, dimension( ims:ime, kms:kme, jms:jme ), &
201	! optional , &
202	! intent(inout) :: rqiblten
203	!
204	real, dimension( ims:ime, jms:jme ) , &
205	optional , &
206	intent(in ) :: mut
207	!
208	real, optional, intent(in ) :: p_top
209	!
210	!-------------------------------------------------------
211	! Local variables
212	integer :: j
213
214	do j = jts,jte
215	call temf2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
216	,tx=t3d(ims,kms,j),thx=th3d(ims,kms,j) &
217	,qvx=qv3d(ims,kms,j),qcx=qc3d(ims,kms,j) &
218	,qix=qi3d(ims,kms,j) &
219	,p2d=p3d(ims,kms,j),p2di=p3di(ims,kms,j) &
220	,pi2d=pi3d(ims,kms,j),rho=rho(ims,kms,j) &
221	,rubltenx=rublten(ims,kms,j),rvbltenx=rvblten(ims,kms,j) &
222	,rthbltenx=rthblten(ims,kms,j),rqvbltenx=rqvblten(ims,kms,j) &
223	,rqcbltenx=rqcblten(ims,kms,j),rqibltenx=rqiblten(ims,kms,j) &
224	,g=g,cp=cp,rcp=rcp,r_d=r_d,r_v=r_v,cpv=cpv &
225	,z2d=z(ims,kms,j) &
226	,xlv=xlv &
227	,psfcpa=psfc(ims,j),znt=znt(ims,j),zsrf=ht(ims,j),ust=ust(ims,j) &
228	,zol=zol(ims,j),hol=hol(ims,j),hpbl=hpbl(ims,j) &
229	,psim=psim(ims,j) &
230	,psih=psih(ims,j),xland=xland(ims,j) &
231	,hfx=hfx(ims,j),qfx=qfx(ims,j) &
232	,tsk=tsk(ims,j),qsfc=qsfc(ims,j),gz1oz0=gz1oz0(ims,j) &
233	,wspd=wspd(ims,j),br=br(ims,j) &
234	,dt=dt,dtmin=dtmin,kpbl1d=kpbl2d(ims,j) &
235	,svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0 &
236	,ep1=ep1,ep2=ep2,karman=karman,eomeg=eomeg &
237	,stbolt=stbolt &
238	,kh_temfx=kh_temf(ims,kms,j),km_temfx=km_temf(ims,kms,j) &
239	,u10=u10(ims,j),v10=v10(ims,j),t2=t2(ims,j) &
240	,te_temfx=te_temf(ims,kms,j) &
241	,shf_temfx=shf_temf(ims,kms,j),qf_temfx=qf_temf(ims,kms,j) &
242	,uw_temfx=uw_temf(ims,kms,j),vw_temfx=vw_temf(ims,kms,j) &
243	,wupd_temfx=wupd_temf(ims,kms,j),mf_temfx=mf_temf(ims,kms,j) &
244	,thup_temfx=thup_temf(ims,kms,j),qtup_temfx=qtup_temf(ims,kms,j) &
245	,qlup_temfx=qlup_temf(ims,kms,j) &
246	,cf3d_temfx=cf3d_temf(ims,kms,j),cfm_temfx=cfm_temf(ims,j) &
247	,hd_temfx=hd_temf(ims,j),lcl_temfx=lcl_temf(ims,j) &
248	,hct_temfx=hct_temf(ims,j),exch_temfx=exch_temf(ims,j) &
249	,flhc=flhc(ims,j),flqc=flqc(ims,j) &
250	,fCor=fCor(ims,j) &
251	,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
252	,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
253	,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
254	enddo
255	!
256	end subroutine temfpbl
257	!
258	!-------------------------------------------------------------------
259	!
260	subroutine temf2d(j,ux,vx,tx,thx,qvx,qcx,qix,p2d,p2di,pi2d,rho, &
261	rubltenx,rvbltenx,rthbltenx, &
262	rqvbltenx,rqcbltenx,rqibltenx, &
263	g,cp,rcp,r_d,r_v,cpv, &
264	z2d, &
265	xlv,psfcpa, &
266	znt,zsrf,ust,zol,hol,hpbl,psim,psih, &
267	xland,hfx,qfx,tsk,qsfc,gz1oz0,wspd,br, &
268	dt,dtmin,kpbl1d, &
269	svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt, &
270	kh_temfx,km_temfx, &
271	u10,v10,t2, &
272	te_temfx,shf_temfx,qf_temfx,uw_temfx,vw_temfx, &
273	wupd_temfx,mf_temfx,thup_temfx,qtup_temfx,qlup_temfx, &
274	cf3d_temfx,cfm_temfx, &
275	hd_temfx,lcl_temfx,hct_temfx,exch_temfx, &
276	flhc,flqc, &
277	fCor, &
278	ids,ide, jds,jde, kds,kde, &
279	ims,ime, jms,jme, kms,kme, &
280	its,ite, jts,jte, kts,kte &
281	)
282	!-------------------------------------------------------------------
283	implicit none
284	!-------------------------------------------------------------------
285	!
286	! This is the Total Energy - Mass Flux (TEMF) PBL scheme.
287	! Initial implementation 2010 by Wayne Angevine, CIRES/NOAA ESRL.
288	! References:
289	! Angevine et al., 2010, MWR
290	! Angevine, 2005, JAM
291	! Mauritsen et al., 2007, JAS
292	!
293	!-------------------------------------------------------------------
294	!
295	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
296	ims,ime, jms,jme, kms,kme, &
297	its,ite, jts,jte, kts,kte, j
298	!
299	real, intent(in ) :: dt,dtmin,g,cp,rcp,r_d,r_v,cpv,xlv
300	!
301	real, intent(in ) :: svp1,svp2,svp3,svpt0
302	real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
303	!
304	real, dimension( ims:ime, kms:kme ), &
305	intent(in) :: z2d
306	!
307	real, dimension( ims:ime, kms:kme ) , &
308	intent(in ) :: ux, vx
309	real, dimension( ims:ime, kms:kme ) , &
310	intent(inout) :: te_temfx
311	real, dimension( ims:ime, kms:kme ) , &
312	intent( out) :: shf_temfx, qf_temfx, uw_temfx, vw_temfx , &
313	wupd_temfx, mf_temfx,thup_temfx, &
314	qtup_temfx, qlup_temfx, cf3d_temfx
315	real, dimension( ims:ime ) , &
316	intent( out) :: hd_temfx, lcl_temfx, hct_temfx, cfm_temfx
317	real, dimension( ims:ime ) , &
318	intent(in ) :: fCor
319	real, dimension( ims:ime ) , &
320	intent(inout) :: flhc, flqc, exch_temfx
321	real, dimension( ims:ime, kms:kme ) , &
322	intent(in ) :: tx, thx, qvx, qcx, qix, pi2d, rho
323	real, dimension( ims:ime, kms:kme ) , &
324	intent(in ) :: p2di, p2d
325	!
326	real, dimension( ims:ime, kms:kme ) , &
327	intent(inout) :: rubltenx, rvbltenx, rthbltenx, &
328	rqvbltenx, rqcbltenx, rqibltenx
329	!
330	real, dimension( ims:ime ) , &
331	intent(inout) :: hol, ust, hpbl, znt
332	real, dimension( ims:ime ) , &
333	intent(in ) :: xland, zsrf
334	real, dimension( ims:ime ) , &
335	intent(inout) :: hfx, qfx
336	!
337	real, dimension( ims:ime ), intent(inout) :: wspd
338	real, dimension( ims:ime ), intent(in ) :: br
339	!
340	real, dimension( ims:ime ), intent(in ) :: psim, psih
341	real, dimension( ims:ime ), intent(in ) :: gz1oz0
342	!
343	real, dimension( ims:ime ), intent(in ) :: psfcpa
344	real, dimension( ims:ime ), intent(in ) :: tsk, qsfc
345	real, dimension( ims:ime ), intent(inout) :: zol
346	integer, dimension( ims:ime ), intent(out ) :: kpbl1d
347	real, dimension( ims:ime, kms:kme ) , &
348	intent(inout) :: kh_temfx, km_temfx
349	!
350	real, dimension( ims:ime ) , &
351	intent(inout) :: u10, v10, t2
352	!
353	!
354	!-----------------------------------------------------------
355	! Local variables
356	!
357	! TE model constants
358	logical, parameter :: MFopt = .true. ! Use mass flux or not
359	! real, parameter :: visc_temf = 1.57e-5
360	! real, parameter :: conduc_temf = 1.57e-5 / 0.733
361	real, parameter :: visc_temf = 1.57e-4 ! WA TEST bigger minimum K
362	real, parameter :: conduc_temf = 1.57e-4 / 0.733
363	real, parameter :: Pr_temf = 0.733
364	real, parameter :: TEmin = 1e-3
365	real, parameter :: ftau0 = 0.17
366	real, parameter :: fth0 = 0.145
367	! real, parameter :: fth0 = 0.12 ! WA 10/13/10 to make PrT0 ~= 1
368	real, parameter :: critRi = 0.25
369	real, parameter :: Cf = 0.185
370	real, parameter :: CN = 2.0
371	! real, parameter :: Ceps = ftau0**1.5
372	real, parameter :: Ceps = 0.070
373	real, parameter :: Cgamma = Ceps
374	real, parameter :: Cphi = Ceps
375	! real, parameter :: PrT0 = Cphi/Ceps * ftau02. / 2 / fth02.
376	real, parameter :: PrT0 = Cphi/Ceps * ftau02 / 2. / fth02
377	! EDMF constants
378	real, parameter :: CM = 0.03 ! Proportionality constant for subcloud MF
379	real, parameter :: Cdelt = 0.006 ! Prefactor for detrainment rate
380	real, parameter :: Cw = 0.5 ! Prefactor for surface wUPD
381	real, parameter :: Cc = 3.0 ! Prefactor for convective length scale
382	real, parameter :: lasymp = 200.0 ! Asymptotic length scale WA 11/20/09
383	real, parameter :: hmax = 4000.0 ! Max hd,hct WA 11/20/09
384	!
385	integer :: i, k, kt ! Loop variable
386	integer, dimension( its:ite) :: h0idx
387	real, dimension( its:ite) :: h0
388	real, dimension( its:ite) :: wstr, ang, wm
389	real, dimension( its:ite) :: hd,lcl,hct,ht
390	real, dimension( its:ite) :: convection_TKE_surface_src, sfcFTE
391	real, dimension( its:ite) :: sfcTHVF
392	real, dimension( its:ite) :: z0t
393	integer, dimension( its:ite) :: hdidx,lclidx,hctidx,htidx
394	integer, dimension( its:ite) :: tval
395	! real, dimension( its:ite ) :: sfcHF, sfcQF
396	real, dimension( its:ite, kts:kte) :: thetal, qt
397	real, dimension( its:ite, kts:kte) :: u_temf, v_temf
398	real, dimension( its:ite, kts:kte) :: rv, rl, rt
399	real, dimension( its:ite, kts:kte) :: chi_poisson, gam
400	real, dimension( its:ite, kts:kte) :: dthdz, dqtdz, dudz, dvdz
401	real, dimension( its:ite, kts:kte) :: lepsmin
402	real, dimension( its:ite, kts:kte) :: thetav
403	real, dimension( its:ite, kts:kte) :: MFCth, MFCq, MFCu, MFCv
404	real, dimension( its:ite, kts:kte) :: MFCql, MFCthv, MFCTE
405	real, dimension( its:ite, kts:kte) :: epsmf, deltmf, dMdz
406	real, dimension( its:ite, kts:kte) :: UUPD, VUPD
407	real, dimension( its:ite, kts:kte) :: thetavUPD, qlUPD, TEUPD
408	real, dimension( its:ite, kts:kte) :: thetavUPDmoist, wupd_dry
409	real, dimension( its:ite, kts:kte) :: B, Bmoist
410	real, dimension( its:ite, kts:kte) :: zm, zt, dzm, dzt
411	real, dimension( its:ite, kts:kte) :: dthUPDdz, dqtup_temfxdz, dwUPDdz
412	real, dimension( its:ite, kts:kte) :: dwUPDmoistdz
413	real, dimension( its:ite, kts:kte) :: dUUPDdz, dVUPDdz, dTEUPDdz
414	real, dimension( its:ite, kts:kte) :: TUPD, rstUPD, rUPD, rlUPD, qstUPD
415	real, dimension( its:ite, kts:kte) :: N2, S, Ri, beta, ftau, fth, ratio
416	real, dimension( its:ite, kts:kte) :: TKE, TE2
417	real, dimension( its:ite, kts:kte) :: ustrtilde, linv, leps
418	real, dimension( its:ite, kts:kte) :: km, kh
419	real, dimension( its:ite, kts:kte) :: Fz, QFK, uwk, vwk
420	real, dimension( its:ite, kts:kte) :: km_conv, kh_conv, lconv
421	real, dimension( its:ite, kts:kte) :: alpha2, beta2 ! For thetav flux calculation
422	real, dimension( its:ite, kts:kte) :: THVF, buoy_src, srcs
423	real, dimension( its:ite, kts:kte) :: u_new, v_new
424	real, dimension( its:ite, kts:kte) :: thx_new, qvx_new, qcx_new
425	real, dimension( its:ite, kts:kte) :: thup_new, qvup_new
426	real, dimension( its:ite, kts:kte) :: beta1 ! For saturation humidity calculations
427	real Cepsmf ! Prefactor for entrainment rate
428	real red_fact ! WA TEST for reducing MF components
429	logical is_convective
430	! Vars for cloud fraction calculation
431	real, dimension( its:ite, kts:kte) :: au, sigq, qst, satdef
432	real sigq2, rst
433
434	!----------------------------------------------------------------------
435	! Grid staggering: Matlab version has mass and turbulence levels.
436	! WRF has full levels (with w) and half levels (u,v,theta,q*). Both
437	! sets of levels use the same indices (kts:kte). See pbl_driver or
438	! WRF Physics doc for (a few) details.
439	! So mass levels correspond to half levels.
440	! WRF full levels are ignored, we define our own turbulence levels
441	! in order to put the first one below the first half level.
442	! Another difference is that
443	! the Matlab version (and the Mauritsen et al. paper) consider the
444	! first mass level to be at z0 (effectively the surface). WRF considers
445	! the first half level to be above the effective surface. The first half
446	! level, at k=1, has nonzero values of u,v for example. Here we convert
447	! all incoming variables to internal ones with the correct indexing
448	! in order to make the code consistent with the Matlab version. We
449	! already had to do this for thetal and qt anyway, so the only additional
450	! overhead is for u and v.
451	! I use suffixes m for mass and t for turbulence as in Matlab for things
452	! like indices.
453	! Note that zsrf is the terrain height ASL, from Registry variable ht.
454	! Translations (Matlab to WRF):
455	! dzt -> calculated below
456	! dzm -> not supplied, calculated below
457	! k -> karman
458	! z0 -> znt
459	! z0t -> not in WRF, calculated below
460	! zt -> calculated below
461	! zm -> (z2d - zsrf) but NOTE zm(1) is now z0 (znt) and zm(2) is
462	! z2d(1) - zsrf
463	!
464	! WA I take the temperature at z0 to be
465	! TSK. This isn't exactly robust. Also I pass out the surface
466	! exchange coefficients flhc, flqc for the surface scheme to use in the
467	! next timestep.
468	! WA 2/16/11 removed calculation of flhc, flqc which are not needed here.
469	! These should be removed from the calling sequence someday.
470	!
471	! Other notes:
472	! - I have often used 1 instead of kts below, because the scheme demands
473	! to know where the surface is. It won't work if kts .NE. 1.
474
475
476	do i = its,ite ! Main loop
477
478	! Get incoming surface theta from TSK (WA for now)
479	thetal(i,1) = tsk(i) / pi2d(i,1) ! WA really should use Exner func. at z0
480	if (exch_temfx(i) > 1.0e-12) then
481	qt(i,1) = qfx(i) / exch_temfx(i) + qvx(i,1) ! WA assumes no liquid at z0
482	else
483	qt(i,1) = qvx(i,1)
484	end if
485	rv(i,1) = qt(i,1) / (1.-qt(i,1)) ! Water vapor
486	rl(i,1) = 0.
487	rt(i,1) = rv(i,1) + rl(i,1) ! Total water (without ice)
488	chi_poisson(i,1) = rcp * (1.+rv(i,1)/ep2) / (1.+rv(i,1)*cpv/cp)
489	gam(i,1) = rv(i,1) * r_v / (cp + rv(i,1)*cpv)
490	! thetav(i,1) = thetal(i,1) * (1. + 0.608*qt(i,1)) ! WA Assumes ql(env)=0, what if it isn't?
491	thetav(i,1) = thetal(i,1) * (1. + 0.608*qt(i,1) - qcx(i,1)) ! WA 4/6/10 allow environment liquid
492	! WA TEST (R5) set z0t = z0
493	! z0t(i) = znt(i) / 10.0 ! WA this is hard coded in Matlab version
494	z0t(i) = znt(i)
495
496	! Convert incoming theta to thetal and qv,qc to qt
497	! NOTE this is where the indexing gets changed from WRF to TEMF basis
498	do k = kts+1,kte
499	! Convert specific humidities to mixing ratios
500	rv(i,k) = qvx(i,k-1) / (1.-qvx(i,k-1)) ! Water vapor
501	rl(i,k) = qcx(i,k-1) / (1.-qcx(i,k-1)) ! Liquid water
502	rt(i,k) = rv(i,k) + rl(i,k) ! Total water (without ice)
503	chi_poisson(i,k) = rcp * (1.+rv(i,k)/ep2) / (1.+rv(i,k)*cpv/cp)
504	gam(i,k) = rt(i,k) * r_v / (cp + rt(i,k)*cpv)
505	thetal(i,k) = thx(i,k-1) * ((ep2+rv(i,k))/(ep2+rt(i,k)))*chi_poisson(i,k) (rv(i,k)/rt(i,k))*(-gam(i,k)) exp( -xlvrl(i,k) / ((cp + rt(i,k)cpv) * tx(i,k)))
506	qt(i,k) = qvx(i,k-1) + qcx(i,k-1)
507	! thetav(i,k) = thetal(i,k) * (1. + 0.608*qt(i,k)) ! WA Assumes ql(env)=0, what if it isn't?
508	thetav(i,k) = thetal(i,k) * (1. + 0.608*qt(i,k) - qcx(i,k-1)) ! WA 4/6/10 allow environment liquid
509	end do
510
511	! Convert incoming u,v to internal u_temf, v_temf
512	! NOTE this is where the indexing gets changed from WRF to TEMF basis
513	u_temf(i,1) = 0. ! zero winds at z0
514	v_temf(i,1) = 0.
515	do k = kts+1,kte
516	u_temf(i,k) = ux(i,k-1)
517	v_temf(i,k) = vx(i,k-1)
518	end do
519
520	! Get delta height at half (mass) levels
521	zm(i,1) = znt(i)
522	dzt(i,1) = z2d(i,1) - zsrf(i) - zm(i,1)
523	! Get height and delta at turbulence levels
524	zt(i,1) = (z2d(i,1) - zsrf(i) - znt(i)) / 2.
525	do kt = kts+1,kte
526	zm(i,kt) = z2d(i,kt-1) - zsrf(i) ! Convert indexing from WRF to TEMF
527	zt(i,kt) = (zm(i,kt) + z2d(i,kt) - zsrf(i)) / 2.
528	dzm(i,kt) = zt(i,kt) - zt(i,kt-1)
529	dzt(i,kt) = z2d(i,kt+1) - z2d(i,kt)
530	end do
531	dzm(i,1) = dzm(i,2) ! WA why?
532	dzt(i,kte) = dzt(i,kte-1) ! WA 12/23/09
533
534	! Gradients at first level
535	dthdz(i,1) = (thetal(i,2)-thetal(i,1)) / (zt(i,1) * log10(zm(i,2)/z0t(i)))
536	dqtdz(i,1) = (qt(i,2)-qt(i,1)) / (zt(i,1) * log10(zm(i,2)/z0t(i)))
537	dudz(i,1) = (u_temf(i,2)-u_temf(i,1)) / (zt(i,1) * log10(zm(i,2)/znt(i)))
538	dvdz(i,1) = (v_temf(i,2)-v_temf(i,1)) / (zt(i,1) * log10(zm(i,2)/znt(i)))
539
540	! Surface thetaV flux from Stull p.147
541	sfcTHVF(i) = hfx(i)/(rho(i,1)cp) (1.+0.608(qvx(i,1)+qcx(i,1))) + 0.608thetav(i,1)*qf_temfx(i,1)
542
543	! WA use hd_temf to calculate w* instead of finding h0 here????
544	! Watch initialization!
545	h0idx(i) = 1
546	h0(i) = zm(i,1)
547
548	! WA TEST (R4) remove lower limit on leps
549	! lepsmin(i,kts) = min(0.4*zt(i,kts), 5.)
550	lepsmin(i,kts) = 0.
551
552	do k = kts+1,kte-1
553	! WA TEST (R4) remove lower limit on leps
554	! lepsmin(i,k) = min(0.4*zt(i,k), 5.)
555	lepsmin(i,k) = 0.
556	! lepsmin(i,k) = min(zt(i,k), 20.) ! WA to deal with runaway
557
558	! Mean gradients
559	! dthdz(i,k) = (thx(i,k) - thx(i,k-1)) / dzt(i,k) ! WA 1/12/10
560	dthdz(i,k) = (thetal(i,k+1) - thetal(i,k)) / dzt(i,k)
561	dqtdz(i,k) = (qt(i,k+1) - qt(i,k)) / dzt(i,k)
562	dudz(i,k) = (u_temf(i,k+1) - u_temf(i,k)) / dzt(i,k)
563	dvdz(i,k) = (v_temf(i,k+1) - v_temf(i,k)) / dzt(i,k)
564
565	! Find h0 (should eventually be interpolated for smoothness)
566	if (thetav(i,k) > thetav(i,1) .AND. h0idx(i) .EQ. 1) then
567	h0idx(i) = k
568	h0(i) = zm(i,k)
569	end if
570	end do
571
572	! Gradients at top level
573
574	dthdz(i,kte) = dthdz(i,kte-1)
575	dqtdz(i,kte) = dqtdz(i,kte-1)
576	dudz(i,kte) = dudz(i,kte-1)
577	dvdz(i,kte) = dvdz(i,kte-1)
578
579	if ( hfx(i) > 0.) then
580	! wstr(i) = (g * h0(i) / thetav(i,2) * shf_temfx(i,1) ) ** (1./3.)
581	wstr(i) = (g * h0(i) / thetav(i,2) * hfx(i)/(rho(i,1)cp) ) * (1./3.)
582	else
583	wstr(i) = 0.
584	end if
585
586
587	! Set flag convective or not for use below
588	is_convective = wstr(i) > 0. .AND. MFopt .AND. dthdz(i,1)<0. .AND. dthdz(i,2)<0. ! WA 12/16/09 require two levels of negative (unstable) gradient
589
590	! Find stability parameters and length scale (on turbulence levels)
591	do kt = 1,kte-1
592	N2(i,kt) = 2. * g / (thetav(i,kt) + thetav(i,kt+1))*dthdz(i,kt)
593	S(i,kt) = sqrt(dudz(i,kt)2. + dvdz(i,kt)2.)
594	Ri(i,kt) = N2(i,kt) / S(i,kt)**2.
595	if (S(i,kt) < 1e-15) then
596	if (N2(i,kt) >= 0) then
597	Ri(i,kt) = 10.
598	else
599	Ri(i,kt) = -1.
600	end if
601	end if
602	beta(i,kt) = 2. * g / (thetav(i,kt)+thetav(i,kt+1))
603	if (Ri(i,kt) > 0) then
604	ratio(i,kt) = Ri(i,kt)/(Cphi*2.ftau0*2./(2.Ceps*2.fth0*2.)+3.Ri(i,kt))
605	ftau(i,kt) = ftau0 * ((3./4.) / (1.+4.*Ri(i,kt)) + 1./4.)
606	fth(i,kt) = fth0 / (1.+4.*Ri(i,kt))
607	TE2(i,kt) = 2. * te_temfx(i,kt) * ratio(i,kt) * N2(i,kt) / beta(i,kt)**2.
608	else
609	ratio(i,kt) = Ri(i,kt)/(Cphi*2.ftau0*2./(-2.Ceps*2.fth0*2.)+2.Ri(i,kt))
610	ftau(i,kt) = ftau0
611	fth(i,kt) = fth0
612	TE2(i,kt) = 0.
613	end if
614	TKE(i,kt) = te_temfx(i,kt) * (1. - ratio(i,kt))
615	ustrtilde(i,kt) = sqrt(ftau(i,kt) * TKE(i,kt))
616	if (N2(i,kt) > 0.) then
617	linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / (Cfustrtilde(i,kt)) + sqrt(N2(i,kt))/(CNustrtilde(i,kt)) + 1./lasymp ! WA Test 11/20/09
618	else
619	linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / (Cf*ustrtilde(i,kt)) + 1./lasymp ! WA Test 11/20/09
620	end if
621	leps(i,kt) = 1./linv(i,kt)
622	leps(i,kt) = max(leps(i,kt),lepsmin(i,kt))
623	end do
624	S(i,kte) = 0.0
625	N2(i,kte) = 0.0
626	TKE(i,kte) = 0.0
627	linv(i,kte) = linv(i,kte-1)
628	leps(i,kte) = leps(i,kte-1)
629
630
631	! Find diffusion coefficients
632	! First use basic formulae for stable and neutral cases,
633	! then for convective conditions, and finally choose the larger
634	! WA 12/23/09 use convective form up to hd/2 always
635	! WA 12/28/09 after restructuring, this block is above MF block,
636	! so hd is not yet available for this timestep, must use h0,
637	! or use hd from previous timestep but be careful about initialization.
638	do kt = 1,kte-1 ! WA 12/22/09
639	! WA 4/8/10 remove beta term to avoid negative and huge values
640	! of km due to very small denominator. This is an interim fix
641	! until we find something better (more theoretically sound).
642	! km(i,kt) = TKE(i,kt)*1.5 ftau(i,kt)*2. / (-beta(i,kt) fth(i,kt) * sqrt(TE2(i,kt)) + Ceps * sqrt(TKE(i,kt)*te_temfx(i,kt)) / leps(i,kt))
643	km(i,kt) = TKE(i,kt)*1.5 ftau(i,kt)*2. / (Ceps sqrt(TKE(i,kt)*te_temfx(i,kt)) / leps(i,kt))
644	kh(i,kt) = 2. * leps(i,kt) * fth(i,kt)*2. TKE(i,kt) / sqrt(te_temfx(i,kt)) / Cphi
645	if ( is_convective) then
646	if (kt <= h0idx(i)) then
647	lconv(i,kt) = 1. / (1. / (karmanzt(i,kt)) + Cc / (karman (h0(i) - zt(i,kt))))
648	else
649	lconv(i,kt) = 0.
650	end if
651	! WA 12/15/09 use appropriate coeffs to match kh_conv and kh at neutral
652	kh_conv(i,kt) = ftau0*2. / Ceps / PrT0 sqrt(TKE(i,kt)) * lconv(i,kt)
653	if (kh_conv(i,kt) < 0.) then
654	kh_conv(i,kt) = 0.
655	end if
656	km_conv(i,kt) = PrT0 * kh_conv(i,kt)
657	if (zt(i,kt) <= h0(i)/2.) then
658	km(i,kt) = km_conv(i,kt)
659	kh(i,kt) = kh_conv(i,kt)
660	end if
661	! WA TEST 1/11/10 go back to max in upper BL
662	if (zt(i,kt) > h0(i)/2. .AND. kt <= h0idx(i)) then
663	km(i,kt) = max(km(i,kt),km_conv(i,kt),visc_temf)
664	kh(i,kt) = max(kh(i,kt),kh_conv(i,kt),conduc_temf)
665	end if
666	end if ! is_convective
667	km(i,kt) = max(km(i,kt),visc_temf)
668	kh(i,kt) = max(kh(i,kt),conduc_temf)
669	Fz(i,kt) = -kh(i,kt) * dthdz(i,kt) ! Diffusive heat flux
670	end do
671	km(i,kte) = km(i,kte-1) ! WA 12/22/09
672	kh(i,kte) = kh(i,kte-1)
673	Fz(i,kte) = 0.0 ! WA 4/2/10
674
675
676	!* Mass flux block starts here *
677
678	if ( is_convective) then
679
680	Cepsmf = 2. / max(200.,h0(i))
681	Cepsmf = max(Cepsmf,0.002) ! WA 7/20/10
682	do k = kts,kte
683	! Calculate lateral entrainment fraction for subcloud layer
684	! epsilon and delta are defined on mass grid (half levels)
685	epsmf(i,k) = Cepsmf
686	end do
687
688	! Initialize updraft
689	thup_temfx(i,1) = thetal(i,1) ! No excess
690	qtup_temfx(i,1) = qt(i,1) ! No excess
691	rUPD(i,1) = qtup_temfx(i,1) / (1. - qtup_temfx(i,1))
692	wupd_temfx(i,1) = Cw * wstr(i)
693	wupd_dry(i,1) = Cw * wstr(i)
694	UUPD(i,1) = u_temf(i,1)
695	VUPD(i,1) = v_temf(i,1)
696	thetavUPD(i,1) = thup_temfx(i,1) * (1. + 0.608*qtup_temfx(i,1)) ! WA Assumes no liquid
697	thetavUPDmoist(i,1) = thup_temfx(i,1) * (1. + 0.608*qtup_temfx(i,1)) ! WA Assumes no liquid
698	TEUPD(i,1) = te_temfx(i,1) + g / thetav(i,1) * sfcTHVF(i)
699	! qlUPD(i,1) = 0.
700	qlUPD(i,1) = qcx(i,1) ! WA allow environment liquid
701	TUPD(i,1) = thup_temfx(i,1) * pi2d(i,1)
702	rstUPD(i,1) = rsat(p2d(i,1),TUPD(i,1),ep2)
703	rlUPD(i,1) = 0.
704
705	! Calculate updraft parameters counting up
706	do k = 2,kte
707	dthUPDdz(i,k-1) = -epsmf(i,k) * (thup_temfx(i,k-1) - thetal(i,k-1))
708	thup_temfx(i,k) = thup_temfx(i,k-1) + dthUPDdz(i,k-1) * dzm(i,k-1)
709	dqtup_temfxdz(i,k-1) = -epsmf(i,k) * (qtup_temfx(i,k-1) - qt(i,k-1))
710	qtup_temfx(i,k) = qtup_temfx(i,k-1) + dqtup_temfxdz(i,k-1) * dzm(i,k-1)
711	thetavUPD(i,k) = thup_temfx(i,k) * (1. + 0.608*qtup_temfx(i,k)) ! WA Assumes no liquid
712	B(i,k-1) = g * (thetavUPD(i,k) - thetav(i,k)) / thetav(i,k)
713	if ( wupd_dry(i,k-1) < 1e-15 ) then
714	wupd_dry(i,k) = 0.
715	else
716	dwUPDdz(i,k-1) = -2. epsmf(i,k)wupd_dry(i,k-1) + 0.33*B(i,k-1)/wupd_dry(i,k-1)
717	wupd_dry(i,k) = wupd_dry(i,k-1) + dwUPDdz(i,k-1) * dzm(i,k-1)
718	end if
719	dUUPDdz(i,k-1) = -epsmf(i,k) * (UUPD(i,k-1) - u_temf(i,k-1))
720	UUPD(i,k) = UUPD(i,k-1) + dUUPDdz(i,k-1) * dzm(i,k-1)
721	dVUPDdz(i,k-1) = -epsmf(i,k) * (VUPD(i,k-1) - v_temf(i,k-1))
722	VUPD(i,k) = VUPD(i,k-1) + dVUPDdz(i,k-1) * dzm(i,k-1)
723	dTEUPDdz(i,k-1) = -epsmf(i,k) * (TEUPD(i,k-1) - te_temfx(i,k-1))
724	TEUPD(i,k) = TEUPD(i,k-1) + dTEUPDdz(i,k-1) * dzm(i,k-1)
725	! Alternative updraft velocity based on moist thetav
726	! Need thetavUPDmoist, qlUPD
727	rUPD(i,k) = qtup_temfx(i,k) / (1. - qtup_temfx(i,k))
728	! WA Updraft temperature assuming no liquid
729	TUPD(i,k) = thup_temfx(i,k) * pi2d(i,k)
730	! Updraft saturation mixing ratio
731	! rstUPD(i,k) = rsat(p2d(i,k),TUPD(i,k),ep2) ! WA 4/19/10
732	rstUPD(i,k) = rsat(p2d(i,k-1),TUPD(i,k),ep2)
733	! Correct to actual temperature (Sommeria & Deardorff 1977)
734	beta1(i,k) = 0.622 * (xlv/(r_dTUPD(i,k))) (xlv/(cp*TUPD(i,k)))
735	rstUPD(i,k) = rstUPD(i,k) * (1.0+beta1(i,k)rUPD(i,k)) / (1.0+beta1(i,k)rstUPD(i,k))
736	qstUPD(i,k) = rstUPD(i,k) / (1. + rstUPD(i,k))
737	if (rUPD(i,k) > rstUPD(i,k)) then
738	rlUPD(i,k) = rUPD(i,k) - rstUPD(i,k)
739	qlUPD(i,k) = rlUPD(i,k) / (1. + rlUPD(i,k))
740	thetavUPDmoist(i,k) = (thup_temfx(i,k) + ((xlv/cp)qlUPD(i,k)/pi2d(i,k))) (1. + 0.608*qstUPD(i,k) - qlUPD(i,k))
741	else
742	rlUPD(i,k) = 0.
743	! qlUPD(i,k) = 0.
744	qlUPD(i,k) = qcx(i,k-1) ! WA 4/6/10 allow environment liquid
745	! WA does this make sense? Should be covered above?
746	thetavUPDmoist(i,k) = thup_temfx(i,k) * (1. + 0.608*qtup_temfx(i,k))
747	end if
748	Bmoist(i,k-1) = g * (thetavUPDmoist(i,k) - thetav(i,k)) / thetav(i,k)
749	if ( wupd_temfx(i,k-1) < 1e-15 ) then
750	wupd_temfx(i,k) = 0.
751	else
752	dwUPDmoistdz(i,k-1) = -2. epsmf(i,k)wupd_temfx(i,k-1) + 0.33*Bmoist(i,k-1)/wupd_temfx(i,k-1)
753	wupd_temfx(i,k) = wupd_temfx(i,k-1) + dwUPDmoistdz(i,k-1) * dzm(i,k-1)
754	end if
755	end do
756
757	! Find hd based on wUPD
758	if (wupd_dry(i,1) == 0.) then
759	hdidx(i) = 1
760	else
761	hdidx(i) = kte ! In case wUPD <= 0 not found
762	do k = 2,kte
763	! if (wupd_dry(i,k) <= 0.) then
764	if (wupd_dry(i,k) <= 0. .OR. zm(i,k) > hmax) then ! WA Test
765	hdidx(i) = k
766	goto 100 ! FORTRAN made me do it!
767	end if
768	end do
769	end if
770	100 hd(i) = zm(i,hdidx(i))
771	! kpbl1d(i) = hd(i) ! WA not sure if this is what I want for diagnostic out to larger WRF universe....and it's not right if not convective
772	kpbl1d(i) = hdidx(i) ! WA 5/11/10 kpbl should be index
773	hpbl(i) = hd(i) ! WA 5/11/10 hpbl is height. Should still be replaced by something that works whether convective or not.
774
775	! Find LCL, hct, and ht
776	lclidx(i) = kte ! In case LCL not found
777	do k = kts,kte
778	if (rUPD(i,k) > rstUPD(i,k)) then
779	lclidx(i) = k
780	goto 200
781	end if
782	end do
783	200 lcl(i) = zm(i,lclidx(i))
784
785	if (hd(i) > lcl(i)) then ! Forced cloud (at least) occurs
786	! Find hct based on wUPDmoist
787	if (wupd_temfx(i,1) == 0.) then
788	hctidx(i) = 1
789	else
790	hctidx(i) = kte ! In case wUPD <= 0 not found
791	do k = 2,kte
792	if (wupd_temfx(i,k) <= 0. .OR. zm(i,k) > hmax) then ! WA Test
793	hctidx(i) = k
794	goto 300 ! FORTRAN made me do it!
795	end if
796	end do
797	end if
798	300 hct(i) = zm(i,hctidx(i))
799	if (hctidx(i) <= hdidx(i)+1) then ! No active cloud
800	hct(i) = hd(i)
801	hctidx(i) = hdidx(i)
802	else
803	end if
804	else ! No cloud
805	hct(i) = hd(i)
806	hctidx(i) = hdidx(i)
807	end if
808	ht(i) = max(hd(i),hct(i))
809	htidx(i) = max(hdidx(i),hctidx(i))
810
811	! Now truncate updraft at ht with taper
812	do k = 1,kte
813	if (zm(i,k) < 0.9*ht(i)) then ! Below taper region
814	tval(i) = 1
815	else if (zm(i,k) >= 0.9ht(i) .AND. zm(i,k) <= 1.0ht(i)) then
816	! Within taper region
817	tval(i) = 1. - ((zm(i,k) - 0.9ht(i)) / (1.0ht(i) - 0.9*ht(i)))
818	else ! Above taper region
819	tval(i) = 0.
820	end if
821	thup_temfx(i,k) = tval(i) * thup_temfx(i,k) + (1-tval(i))*thetal(i,k)
822	thetavUPD(i,k) = tval(i) * thetavUPD(i,k) + (1-tval(i))*thetav(i,k)
823	qtup_temfx(i,k) = tval(i) * qtup_temfx(i,k) + (1-tval(i)) * qt(i,k)
824	qlUPD(i,k) = tval(i) * qlUPD(i,k) + (1-tval(i)) * qcx(i,k-1)
825	UUPD(i,k) = tval(i) * UUPD(i,k) + (1-tval(i)) * u_temf(i,k)
826	VUPD(i,k) = tval(i) * VUPD(i,k) + (1-tval(i)) * v_temf(i,k)
827	TEUPD(i,k) = tval(i) * TEUPD(i,k) + (1-tval(i)) * te_temfx(i,k)
828	if (zm(i,k) > ht(i)) then ! WA this is just for cleanliness
829	wupd_temfx(i,k) = 0.
830	dwUPDmoistdz(i,k) = 0.
831	wupd_dry(i,k) = 0.
832	dwUPDdz(i,k) = 0.
833	end if
834	end do
835
836	! Calculate lateral detrainment rate for cloud layer
837	deltmf(i,1) = Cepsmf
838	do k = 2,kte-1
839	if (hctidx(i) > hdidx(i)+1) then ! Some cloud
840	deltmf(i,k) = 0.9 * Cepsmf + Cdelt * (atan((zm(i,k)-(lcl(i)+(hct(i)-lcl(i))/1.5))/((hct(i)-lcl(i))/8))+(3.1415926/2))/3.1415926
841	else if (k < hdidx(i)) then ! No cloud, below hd
842	deltmf(i,k) = Cepsmf + 0.05 * 1. / (hd(i) - zm(i,k))
843	else if (k >= hdidx(i)) then ! No cloud, above hd
844	deltmf(i,k) = deltmf(i,k-1)
845	end if
846	end do
847
848	! Calculate mass flux (defined on turbulence levels)
849	mf_temfx(i,1) = CM * wstr(i)
850	do kt = 2,kte-1
851	dMdz(i,kt) = (epsmf(i,kt) - deltmf(i,kt)) * mf_temfx(i,kt-1) * dzt(i,kt)
852	mf_temfx(i,kt) = mf_temfx(i,kt-1) + dMdz(i,kt)
853	end do
854
855	! WA 12/28/09 If mass flux component > diffusive
856	! component at second level,
857	! reduce M to prevent a stable layer
858	MFCth(i,2) = mf_temfx(i,2) * (thup_temfx(i,2)-thetal(i,2) + thup_temfx(i,3)-thetal(i,3)) / 2.
859	if (MFCth(i,2) > Fz(i,2)) then
860	red_fact = Fz(i,2) / MFCth(i,2)
861	do kt = 1,kte
862	mf_temfx(i,kt) = mf_temfx(i,kt) * red_fact
863	end do
864	end if ! Reduce M to prevent stable layer at second level
865
866	! Calculate mass flux contributions to fluxes (defined on turb levels)
867	! Use log interpolation at first level
868	MFCth(i,1) = mf_temfx(i,1) * (thup_temfx(i,1)-thetal(i,1) + (thup_temfx(i,2)-thetal(i,2) - (thup_temfx(i,1)-thetal(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i)))
869	MFCq(i,1) = mf_temfx(i,1) * (qtup_temfx(i,1)-qt(i,1) + (qtup_temfx(i,2)-qt(i,2) - (qtup_temfx(i,1)-qt(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i)))
870	MFCu(i,1) = mf_temfx(i,1) * (UUPD(i,1)-u_temf(i,1) + (UUPD(i,2)-u_temf(i,2) - (UUPD(i,1)-u_temf(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i)))
871	MFCv(i,1) = mf_temfx(i,1) * (VUPD(i,1)-v_temf(i,1) + (VUPD(i,2)-v_temf(i,2) - (VUPD(i,1)-v_temf(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i)))
872	MFCql(i,1) = mf_temfx(i,1) * (qlUPD(i,1)-qcx(i,1) + (qlUPD(i,2)-qcx(i,2) - (qlUPD(i,1)-qcx(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i)))
873	MFCTE(i,1) = mf_temfx(i,1) * (TEUPD(i,1)-te_temfx(i,1) + (TEUPD(i,2)-te_temfx(i,2) - (TEUPD(i,1)-te_temfx(i,1))) * log(zt(i,1)/znt(i))/log(zm(i,2)/znt(i))) ! WA Check this
874	do kt = 2,kte-1
875	MFCth(i,kt) = mf_temfx(i,kt) * (thup_temfx(i,kt)-thetal(i,kt) + thup_temfx(i,kt+1)-thetal(i,kt+1)) / 2.
876	MFCq(i,kt) = mf_temfx(i,kt) * (qtup_temfx(i,kt)-qt(i,kt) + qtup_temfx(i,kt+1)-qt(i,kt+1)) / 2.
877	MFCu(i,kt) = mf_temfx(i,kt) * (UUPD(i,kt)-u_temf(i,kt) + UUPD(i,kt+1)-u_temf(i,kt+1)) / 2.
878	MFCv(i,kt) = mf_temfx(i,kt) * (VUPD(i,kt)-v_temf(i,kt) + VUPD(i,kt+1)-v_temf(i,kt+1)) / 2.
879	MFCql(i,kt) = mf_temfx(i,kt) * (qlUPD(i,kt)-qcx(i,kt-1) + qlUPD(i,kt+1)-qcx(i,kt)) / 2.
880	MFCTE(i,kt) = mf_temfx(i,kt) * (TEUPD(i,kt)-te_temfx(i,kt)) ! TE is on turb levels
881	end do
882	MFCth(i,kte) = 0
883	MFCq(i,kte) = 0
884	MFCu(i,kte) = 0
885	MFCv(i,kte) = 0
886	MFCql(i,kte) = 0
887	MFCTE(i,kte) = 0
888
889	! Calculate cloud fraction (on mass levels)
890	cf3d_temfx(i,1) = 0.0
891	cfm_temfx(i) = 0.0
892	do k = 2,kte
893	! if (wupd_temfx(i,k-1) >= 1.0e-15 .AND. wupd_temfx(i,k) >= 1.0e-15 .AND. .NOT. isnan(wupd_temfx(i,k-1)) .AND. .NOT. isnan(wupd_temfx(i,k))) then
894	if (wupd_temfx(i,k-1) >= 1.0e-15 .AND. wupd_temfx(i,k) >= 1.0e-15) then
895	au(i,k) = ((mf_temfx(i,k-1)+mf_temfx(i,k))/2.0) / ((wupd_temfx(i,k-1)+wupd_temfx(i,k))/2.0) ! WA average before divide, is that best?
896	else
897	au(i,k) = 0.0
898	end if
899	sigq2 = au(i,k) * (qtup_temfx(i,k)-qt(i,k))
900	if (sigq2 > 0.0) then
901	sigq(i,k) = sqrt(sigq2)
902	else
903	sigq(i,k) = 0.0
904	end if
905	! rst = rsat(p2d(i,k),thx(i,k)*pi2d(i,k),ep2)
906	rst = rsat(p2d(i,k-1),thx(i,k-1)*pi2d(i,k-1),ep2)
907	qst(i,k) = rst / (1. + rst)
908	satdef(i,k) = qt(i,k) - qst(i,k)
909	if (satdef(i,k) <= 0.0) then
910	if (sigq(i,k) > 1.0e-15) then
911	cf3d_temfx(i,k) = max(0.5 + 0.36 * atan(1.55*(satdef(i,k)/sigq(i,k))),0.0)
912	else
913	cf3d_temfx(i,k) = 0.0
914	end if
915	else
916	cf3d_temfx(i,k) = 1.0
917	end if
918	if (zm(i,k) < lcl(i)) then
919	cf3d_temfx(i,k) = 0.0
920	end if
921	! Put max value so far into cfm
922	if (zt(i,k) <= hmax) then
923	cfm_temfx(i) = max(cf3d_temfx(i,k),cfm_temfx(i))
924	end if
925	end do
926
927	else ! not is_convective, no MF components
928	do kt = 1,kte
929	MFCth(i,kt) = 0
930	MFCq(i,kt) = 0
931	MFCu(i,kt) = 0
932	MFCv(i,kt) = 0
933	MFCql(i,kt) = 0
934	MFCTE(i,kt) = 0
935	end do
936	lcl(i) = zm(i,kte-1)
937	hct(i) = zm(i,1)
938	hctidx(i) = 1
939	hd(i) = zm(i,1)
940	hdidx(i) = 1
941	ht(i) = hd(i)
942	! Cloud fraction calculations
943	cf3d_temfx(i,1) = 0.0
944	cfm_temfx(i) = 0.0
945	do k = 2,kte
946	if (qcx(i,k-1) > 1.0e-15) then
947	cf3d_temfx(i,k) = 1.0
948	else
949	cf3d_temfx(i,k) = 0.0
950	end if
951	! Put max value so far into cfm
952	if (zt(i,k) <= hmax) then
953	cfm_temfx(i) = max(cf3d_temfx(i,k),cfm_temfx(i))
954	end if
955	end do
956
957	end if ! MF components or not
958	cf3d_temfx(i,kte) = 0.0
959	! Mass flux block ends here
960
961	! Flux profiles
962	do kt = 2,kte
963	! Fz(i,kt) = -kh(i,kt) * dthdz(i,kt)
964	shf_temfx(i,kt) = Fz(i,kt) + MFCth(i,kt)
965	QFK(i,kt) = -kh(i,kt) * dqtdz(i,kt)
966	qf_temfx(i,kt) = QFK(i,kt) + MFCq(i,kt)
967	uwk(i,kt) = -km(i,kt) * dudz(i,kt)
968	uw_temfx(i,kt) = uwk(i,kt) + MFCu(i,kt)
969	vwk(i,kt) = -km(i,kt) * dvdz(i,kt)
970	vw_temfx(i,kt) = vwk(i,kt) + MFCv(i,kt)
971	end do
972
973	! Surface momentum fluxes
974	ust(i) = sqrt(ftau(i,1)/ftau0) * sqrt(u_temf(i,2)2. + v_temf(i,2)2.) * leps(i,1) / log(zm(i,2)/znt(i)) / zt(i,1)
975	ang(i) = atan2(v_temf(i,2),u_temf(i,2))
976	uw_temfx(i,1) = -cos(ang(i)) * ust(i)**2.
977	vw_temfx(i,1) = -sin(ang(i)) * ust(i)**2.
978
979	! Calculate mixed scaling velocity (Moeng & Sullivan 1994 JAS p.1021)
980	! Replaces ust everywhere (WA need to reconsider?)
981	! wm(i) = (1./5. * (wstr(i)*3. + 5. ust(i)3.)) (1./3.)
982	! WA TEST (R2,R11) 7/23/10 reduce velocity scale to fix excessive fluxes
983	wm(i) = 0.5 * (1./5. * (wstr(i)*3. + 5. ust(i)3.)) (1./3.)
984	! WA TEST 2/14/11 limit contribution of w*
985	! wm(i) = 0.5 * (1./5. * (min(0.8,wstr(i))*3. + 5. ust(i)3.)) (1./3.)
986	! WA TEST (R3-R11) 7/23/10 wm = u*
987	! wm(i) = ust(i)
988
989	! Specified flux versions (flux is modified by land surface)
990	shf_temfx(i,1) = hfx(i)/(rho(i,1)cp) + (shf_temfx(i,2) - hfx(i)/(rho(i,1)cp)) * (zt(i,2)-zt(i,1)) / (zt(i,2)-znt(i))
991	qf_temfx(i,1) = qfx(i)/rho(i,1) + (qf_temfx(i,2)-qfx(i)/rho(i,1)) * (zt(i,2)-zt(i,1)) / (zt(i,2)-znt(i))
992	Fz(i,1) = shf_temfx(i,1) - MFCth(i,1)
993	QFK(i,1) = qf_temfx(i,1) - MFCq(i,1)
994
995
996	! Calculate thetav and its flux
997	! From Lewellen 2004 eq. 3
998	! WA The constants and combinations below should probably be replaced
999	! by something more consistent with the WRF system, but for now
1000	! I don't want to take the chance of breaking something.
1001	do kt = 2,kte-1
1002	alpha2(i,kt) = 0.61 * (thetal(i,kt) + thetal(i,kt+1)) / 2.
1003	beta2(i,kt) = (100000. / p2di(i,kt))*0.286 2.45e-6 / 1004.67 - 1.61 * (thetal(i,kt) + thetal(i,kt+1)) / 2.
1004	end do
1005	alpha2(i,1) = 0.61 * (thetal(i,1) + (thetal(i,2)-thetal(i,1)) * (zt(i,2) - znt(i)) / (zm(i,2) - znt(i)))
1006	alpha2(i,kte) = 0.61 * thetal(i,kte)
1007	beta2(i,1) = (100000. / p2di(i,1))*0.286 2.45e-6 / 1004.67 - 1.61 * (thetal(i,1) + (thetal(i,2) - thetal(i,1)) * (zt(i,2) - znt(i)) / (zm(i,2) - znt(i)))
1008	beta2(i,kte) = (100000. / p2di(i,kte))*0.286 2.45e-6 / 1004.67 - 1.61 * thetal(i,kte)
1009	if ( is_convective ) then ! Activate EDMF components
1010	do kt = 1,kte-1
1011	MFCthv(i,kt) = (1. + 0.61 * (qtup_temfx(i,kt)+qtup_temfx(i,kt+1))) / 2. * MFCth(i,kt) + alpha2(i,kt) * MFCq(i,kt) + beta2(i,kt) * MFCql(i,kt)
1012	end do
1013	MFCthv(i,kte) = 0.
1014	else ! No MF components
1015	do kt = 1,kte
1016	MFCthv(i,kt) = 0.
1017	end do
1018	end if
1019
1020	do kt = 1,kte
1021	THVF(i,kt) = (1. + 0.61 * qt(i,kt)) * Fz(i,kt) + alpha2(i,kt) * QFK(i,kt) + MFCthv(i,kt)
1022	end do
1023
1024	! Update mean variables:
1025	! This is done with implicit solver for diffusion part followed
1026	! by explicit solution for MF terms.
1027	! Note that Coriolis terms that were source terms for u and v
1028	! in Matlab code are now handled by WRF outside this PBL context.
1029
1030	u_new(i,:) = u_temf(i,:)
1031	call solve_implicit_temf(km(i,kts:kte-1),u_new(i,kts+1:kte),uw_temfx(i,1),dzm(i,kts:kte-1),dzt(i,kts:kte-1),kts,kte-1,dt,.FALSE.)
1032	do k = 2,kte-1
1033	u_new(i,k) = u_new(i,k) + dt * (-(MFCu(i,k)-MFCu(i,k-1))) / dzm(i,k)
1034	end do
1035
1036	v_new(i,:) = v_temf(i,:)
1037	call solve_implicit_temf(km(i,kts:kte-1),v_new(i,kts+1:kte),vw_temfx(i,1),dzm(i,kts:kte-1),dzt(i,kts:kte-1),kts,kte-1,dt,.FALSE.)
1038	do k = 2,kte-1
1039	v_new(i,k) = v_new(i,k) + dt * (-(MFCv(i,k)-MFCv(i,k-1))) / dzm(i,k)
1040	end do
1041
1042	call solve_implicit_temf(kh(i,kts:kte-1),thetal(i,kts+1:kte),Fz(i,1),dzm(i,kts:kte-1),dzt(i,kts:kte-1),kts,kte-1,dt,.FALSE.)
1043	do k = 2,kte-1
1044	thetal(i,k) = thetal(i,k) + dt * (-(MFCth(i,k)-MFCth(i,k-1))) / dzm(i,k)
1045	end do
1046
1047	call solve_implicit_temf(kh(i,kts:kte-1),qt(i,kts+1:kte),QFK(i,1),dzm(i,kts:kte-1),dzt(i,kts:kte-1),kts,kte-1,dt,.FALSE.)
1048	do k = 2,kte-1
1049	qt(i,k) = qt(i,k) + dt * (-(MFCq(i,k)-MFCq(i,k-1))) / dzm(i,k)
1050	end do
1051
1052	! Stepping TE forward is a bit more complicated
1053	te_temfx(i,1) = ust(i)*2. / ftau(i,1) (1. + ratio(i,1))
1054	if ( is_convective ) then
1055	! WA currently disabled if MFopt=false, is that right?
1056	convection_TKE_surface_src(i) = 2. * beta(i,1) * shf_temfx(i,1)
1057	else
1058	convection_TKE_surface_src(i) = 0.
1059	end if
1060	te_temfx(i,1) = max(te_temfx(i,1), (leps(i,1) / Cgamma * (ust(i)*2. S(i,1) + convection_TKE_surface_src(i)))**(2./3.))
1061	sfcFTE(i) = -(km(i,1)+km(i,2)) / 2. * (te_temfx(i,2)-te_temfx(i,1)) / dzm(i,2)
1062
1063	do kt = 1,kte
1064	if (THVF(i,kt) >= 0) then
1065	buoy_src(i,kt) = 2. * g / thetav(i,kt) * THVF(i,kt)
1066	else
1067	buoy_src(i,kt) = 0. ! Cancel buoyancy term when locally stable
1068	end if
1069	srcs(i,kt) = -uw_temfx(i,kt) * dudz(i,kt) - vw_temfx(i,kt) * dvdz(i,kt) - Cgamma * te_temfx(i,kt)*1.5 linv(i,kt) + buoy_src(i,kt)
1070	end do
1071	call solve_implicit_temf((km(i,kts:kte-1)+km(i,kts+1:kte))/2.0,te_temfx(i,kts+1:kte),sfcFTE(i),dzt(i,kts+1:kte),dzt(i,kts:kte-1),kts,kte-1,dt,.false.)
1072	do kt = 2,kte-1
1073	te_temfx(i,kt) = te_temfx(i,kt) + dt * srcs(i,kt)
1074	te_temfx(i,kt) = te_temfx(i,kt) + dt * (-(MFCTE(i,kt)-MFCTE(i,kt-1))) / dzt(i,kt)
1075	if (te_temfx(i,kt) < TEmin) te_temfx(i,kt) = TEmin
1076	end do
1077	te_temfx(i,kte) = 0.0 ! WA 4/2/10
1078	do kt = 2,kte-1
1079	if (te_temfx(i,kt) > 30.0) then ! WA DEBUG
1080	te_temfx(i,kt) = 30.0 ! WA TEST limit max TE
1081	end if
1082	end do
1083
1084	! Done with updates, now convert internal variables back to WRF vars
1085	do k = kts,kte
1086	! Populate kh_temfx, km_temfx from kh, km
1087	kh_temfx(i,k) = kh(i,k)
1088	km_temfx(i,k) = km(i,k)
1089	end do
1090
1091	! Convert thetal, qt back to theta, qv, qc
1092	! See opposite conversion at top of subroutine
1093	! WA this accounts for offset of indexing between
1094	! WRF and TEMF, see notes at top of this routine.
1095	call thlqt2thqvqc(thetal(i,kts+1:kte),qt(i,kts+1:kte),thx_new(i,kts:kte-1),qvx_new(i,kts:kte-1),qcx_new(i,kts:kte-1),p2d(i,kts:kte-1),pi2d(i,kts:kte-1),kts,kte-1,ep2,xlv,cp)
1096
1097	do k = kts,kte-1
1098	! Calculate tendency terms
1099	! WA this accounts for offset of indexing between
1100	! WRF and TEMF, see notes at top of this routine.
1101	rubltenx(i,k) = (u_new(i,k+1) - u_temf(i,k+1)) / dt
1102	rvbltenx(i,k) = (v_new(i,k+1) - v_temf(i,k+1)) / dt
1103	rthbltenx(i,k) = (thx_new(i,k) - thx(i,k)) / dt
1104	rqvbltenx(i,k) = (qvx_new(i,k) - qvx(i,k)) / dt
1105	rqcbltenx(i,k) = (qcx_new(i,k) - qcx(i,k)) / dt
1106	end do
1107	rubltenx(i,kte) = 0.
1108	rvbltenx(i,kte) = 0.
1109	rthbltenx(i,kte) = 0.
1110	rqvbltenx(i,kte) = 0.
1111	rqcbltenx(i,kte) = 0.
1112
1113	! Populate surface exchange coefficient variables to go back out
1114	! for next time step of surface scheme
1115	! Unit specifications in SLAB and sfclay are conflicting and probably
1116	! incorrect. This will give a dynamic heat flux (W/m^2) or moisture
1117	! flux (kg(water)/(m^2*s)) when multiplied by a difference.
1118	! These formulae are the same as what's used above to get surface
1119	! flux from surface temperature and specific humidity.
1120	! WA 2/16/11 removed, not needed in BL
1121	! flhc(i) = rho(i,1) * cp * fth(i,1)/fth0 * wm(i) * leps(i,1) / PrT0 / log(zm(i,2)/z0t(i)) / zt(i,1)
1122	! flqc(i) = rho(i,1) * fth(i,1)/fth0 * wm(i) * leps(i,1) / PrT0 / log(zm(i,2)/z0t(i)) / zt(i,1)
1123	! WA Must exchange coeffs be limited to avoid runaway in some
1124	! (convective?) conditions? Something like this is done in sfclay.
1125	! Doing nothing for now.
1126
1127	! Populate 10 m winds and 2 m temp
1128	! WA Note this only works if first mass level is above 10 m
1129	u10(i) = u_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i))
1130	v10(i) = v_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i))
1131	t2(i) = (tsk(i)/pi2d(i,1) + (thx_new(i,1) - tsk(i)/pi2d(i,1)) * log(2.0/z0t(i)) / log(zm(i,2)/z0t(i))) * pi2d(i,1) ! WA this should also use pi at z0
1132
1133	! Populate diagnostic variables
1134	hd_temfx(i) = hd(i)
1135	lcl_temfx(i) = lcl(i)
1136	hct_temfx(i) = hct(i)
1137
1138	! Send updraft liquid water back
1139	if ( is_convective) then
1140	do k = kts,kte-1
1141	qlup_temfx(i,k) = qlUPD(i,k)
1142	end do
1143	else
1144	qlup_temfx(i,1) = qcx(i,1)
1145	do k = kts+1,kte-1
1146	qlup_temfx(i,k) = qcx(i,k-1)
1147	end do
1148	end if
1149	qlup_temfx(i,kte) = qcx(i,kte)
1150
1151	end do ! Main (i) loop
1152
1153	end subroutine temf2d
1154	!
1155	!--------------------------------------------------------------------
1156	!
1157	subroutine thlqt2thqvqc(thetal,qt,theta,qv,qc,p,piex,kbot,ktop,ep2,Lv,Cp)
1158
1159	! Calculates theta, qv, qc from thetal, qt.
1160	! Originally from RAMS (subroutine satadjst) by way of Hongli Jiang.
1161
1162	implicit none
1163	integer, intent(in ) :: kbot, ktop
1164	real, dimension( kbot:ktop ), intent(in ) :: thetal, qt
1165	real, dimension( kbot:ktop ), intent( out) :: theta, qv, qc
1166	real, dimension( kbot:ktop ), intent(in ) :: p, piex
1167	real, intent(in ) :: ep2, Lv, Cp
1168	!
1169	! Local variables
1170	integer :: k, iterate
1171	real :: T1, Tt
1172	real, dimension( kbot:ktop) :: rst
1173	real, dimension( kbot:ktop) :: Tair, rc, rt, rv
1174	!
1175	do k = kbot,ktop
1176	T1 = thetal(k) * piex(k) ! First guess T is just thetal converted to T
1177	Tair(k) = T1
1178	rt(k) = qt(k) / (1. - qt(k))
1179
1180	do iterate = 1,20
1181	rst(k) = rsat(p(k),Tair(k),ep2)
1182	rc(k) = max(rt(k) - rst(k), 0.)
1183	Tt = 0.7Tair(k) + 0.3T1 * (1.+Lvrc(k) / (Cpmax(Tair(k),253.)))
1184	if ( abs(Tt - Tair(k)) < 0.001) GOTO 100
1185	Tair(k) = Tt
1186	end do
1187	100 continue
1188	rv(k) = rt(k) - rc(k)
1189	qv(k) = rv(k) / (1. + rv(k))
1190	qc(k) = rc(k) / (1. + rc(k))
1191	theta(k) = Tair(k) / piex(k)
1192	end do ! k loop
1193	return
1194	end subroutine thlqt2thqvqc
1195	!
1196	!--------------------------------------------------------------------
1197	!
1198	subroutine findhct_te( thetavenv,thetaparin,qpar, &
1199	rpar,hdidx,paridx,zm,hct,hctidx,p,piex,ep2,kbot,ktop)
1200
1201	! Calculates the cloud top height (limit of convection) for the
1202	! updraft properties. hct is the level at which a parcel lifted
1203	! at the moist adiabatic rate where it is saturated and at the dry
1204	! adiabatic rate otherwise, first has thetav cooler than the environment.
1205	! Loops start at LCL (paridx).
1206
1207	implicit none
1208	integer, intent(in) :: kbot, ktop
1209	integer, intent(in) :: paridx, hdidx
1210	real, intent(in) :: ep2
1211	real, dimension( kbot:ktop), intent(in) :: thetavenv
1212	real, dimension( kbot:ktop), intent(in) :: thetaparin
1213	real, dimension( kbot:ktop), intent(in) :: qpar, rpar, zm, p, piex
1214	real, intent(out) :: hct
1215	integer, intent(out) :: hctidx
1216	! Local variables
1217	integer k
1218	real, dimension( kbot:ktop) :: thetapar, thetavpar, qlpar, Tpar, rsatpar
1219	real, dimension( kbot:ktop) :: qsatpar
1220	real :: gammas, TparC
1221
1222	thetapar(paridx) = thetaparin(paridx)
1223	Tpar(paridx) = thetapar(paridx) * piex(paridx)
1224	hctidx = ktop ! In case hct not found
1225	do k = paridx+1,ktop
1226	! Find saturation mixing ratio at parcel temperature
1227	rsatpar(k) = rsat(p(k-1),Tpar(k-1),ep2)
1228	qsatpar(k) = rsatpar(k) / (1. + rsatpar(k))
1229
1230	! When parcel is unsaturated, its temperature changes
1231	! at dry adiabatic rate (in other words, theta is constant).
1232	if (rpar(k) <= rsatpar(k)) then
1233	thetapar(k) = thetapar(k-1)
1234	Tpar(k) = thetapar(k) * piex(k)
1235	thetavpar(k) = thetapar(k) * (1.+0.608*qpar(k))
1236	else
1237	! When parcel is saturated, its temperature changes at
1238	! moist adiabatic rate
1239	! Calculate moist adiabatic lapse rate according to Gill A4.12
1240	! Requires T in deg.C
1241	TparC = Tpar(k-1) - 273.15
1242	gammas = 6.4 - 0.12 * TparC + 2.5e-5 * TparC*3. + (-2.4 + 1.e-3 (TparC-5.)*2.) (1. - p(k-1)/100000.)
1243	Tpar(k) = Tpar(k-1) - gammas/1000. * (zm(k)-zm(k-1))
1244	thetapar(k) = Tpar(k) / piex(k)
1245	qlpar(k) = qpar(k) - qsatpar(k) ! Liquid water in parcel
1246	thetavpar(k) = thetapar(k) * (1. + 0.608 * qsatpar(k) - qlpar(k))
1247	end if
1248	if (thetavenv(k) >= thetavpar(k)) then
1249	hctidx = k
1250	goto 1000
1251	end if
1252	end do
1253	1000 hct = zm(hctidx)
1254	return
1255	end subroutine findhct_te
1256	!
1257	!--------------------------------------------------------------------
1258	!
1259	real function rsat(p,T,ep2)
1260
1261	! Calculates the saturation mixing ratio with respect to liquid water
1262	! Arguments are pressure (Pa) and absolute temperature (K)
1263	! Uses the formula from the ARM intercomparison setup.
1264	! Converted from Matlab by WA 7/28/08
1265
1266	implicit none
1267	real p, T, ep2
1268	real temp, x
1269	real, parameter :: c0 = 0.6105851e+3
1270	real, parameter :: c1 = 0.4440316e+2
1271	real, parameter :: c2 = 0.1430341e+1
1272	real, parameter :: c3 = 0.2641412e-1
1273	real, parameter :: c4 = 0.2995057e-3
1274	real, parameter :: c5 = 0.2031998e-5
1275	real, parameter :: c6 = 0.6936113e-8
1276	real, parameter :: c7 = 0.2564861e-11
1277	real, parameter :: c8 = -0.3704404e-13
1278
1279	temp = T - 273.15
1280
1281	x =c0+temp(c1+temp(c2+temp(c3+temp(c4+temp(c5+temp(c6+temp(c7+tempc8)))))))
1282	rsat = ep2*x/(p-x)
1283
1284	return
1285	end function rsat
1286	!
1287	!--------------------------------------------------------------------
1288	!
1289	subroutine solve_implicit_temf(Khlf,psi_n,srf_flux,dzm,dzt,kbot,ktop,dt,print_flag)
1290
1291	! Implicit solution of vertical diffusion for conserved variable
1292	! psi given diffusivity Khlf on turbulence levels,
1293	! and surface flux srf_flux.
1294	! dzm is delta_z of mass levels, dzt is delta_z of turbulence levels.
1295	! dt is timestep (s).
1296
1297	implicit none
1298	integer :: kbot, ktop
1299	logical :: print_flag
1300	real :: srf_flux, dt
1301	real, dimension( kbot:ktop ), intent(in ) :: Khlf
1302	real, dimension( kbot:ktop ), intent(in ) :: dzm, dzt
1303	real, dimension( kbot:ktop ), intent(inout) :: psi_n
1304	!
1305	! Local variables
1306	integer :: k
1307	real, dimension( kbot:ktop ) :: AU, BU, CU, YU
1308	!
1309	AU(kbot) = Khlf(kbot) / (dzm(kbot)*dzt(kbot))
1310	BU(kbot) = -1.0/dt - Khlf(kbot+1)/(dzm(kbot+1)*dzt(kbot+1))
1311	CU(kbot) = Khlf(kbot+1)/(dzm(kbot)*dzt(kbot+1))
1312	YU(kbot) = -psi_n(kbot)/dt - srf_flux/dzm(kbot)
1313
1314	do k = kbot+1,ktop-1
1315	! Subdiagonal (A) vector
1316	AU(k) = Khlf(k) / (dzm(k) * dzt(k))
1317	! Main diagonal (B) vector
1318	BU(k) = -1.0/dt - (Khlf(k)/dzt(k) + Khlf(k+1)/dzt(k+1)) / dzm(k)
1319	! Superdiagonal (C) vector
1320	CU(k) = Khlf(k+1) / (dzm(k)*dzt(k+1))
1321	! Result vector
1322	YU(k) = -psi_n(k)/dt
1323	end do ! k loop
1324
1325	AU(ktop) = 0.
1326	BU(ktop) = -1.0 / dt
1327	YU(ktop) = -psi_n(ktop) / dt
1328
1329	! Compute result with tridiagonal routine
1330	psi_n = trid(AU,BU,CU,YU,kbot,ktop)
1331
1332	return
1333	end subroutine solve_implicit_temf
1334	!
1335	!--------------------------------------------------------------------
1336	!
1337	function trid(a,b,c,r,kbot,ktop)
1338
1339	! Solves tridiagonal linear system.
1340	! Inputs are subdiagonal vector a, main diagonal b, superdiagonal c,
1341	! result r, column top and bottom indices kbot and ktop.
1342	! Originally from Numerical Recipes section 2.4.
1343
1344	implicit none
1345	real, dimension( kbot:ktop ) :: trid
1346	integer :: kbot, ktop
1347	real, dimension( kbot:ktop ), intent(in ) :: a, b, c, r
1348	!
1349	! Local variables
1350	integer :: k
1351	real :: bet
1352	real, dimension( kbot:ktop ) :: gam, u
1353	!
1354	bet = b(kbot)
1355	u(kbot) = r(kbot) / bet
1356
1357	do k = kbot+1,ktop
1358	gam(k) = c(k-1) / bet
1359	bet = b(k) - a(k)*gam(k)
1360	u(k) = (r(k) - a(k)*u(k-1)) / bet
1361	end do
1362
1363	do k = ktop-1,kbot,-1
1364	u(k) = u(k) - gam(k+1)*u(k+1)
1365	end do
1366
1367	trid = u
1368
1369	return
1370	end function trid
1371	!
1372	!--------------------------------------------------------------------
1373	!
1374	subroutine temfinit(rublten,rvblten,rthblten,rqvblten, &
1375	rqcblten,rqiblten,p_qi,p_first_scalar, &
1376	restart, allowed_to_read, &
1377	te_temf, cf3d_temf, &
1378	ids, ide, jds, jde, kds, kde, &
1379	ims, ime, jms, jme, kms, kme, &
1380	its, ite, jts, jte, kts, kte )
1381	!-------------------------------------------------------------------
1382	implicit none
1383	!-------------------------------------------------------------------
1384	!
1385	logical , intent(in) :: restart, allowed_to_read
1386	integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1387	ims, ime, jms, jme, kms, kme, &
1388	its, ite, jts, jte, kts, kte
1389	integer , intent(in) :: p_qi,p_first_scalar
1390	real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
1391	rublten, &
1392	rvblten, &
1393	rthblten, &
1394	rqvblten, &
1395	rqcblten, &
1396	rqiblten, &
1397	te_temf, &
1398	cf3d_temf
1399	! Local variables
1400	integer :: i, j, k, itf, jtf, ktf
1401	real, parameter :: TEmin = 1e-3
1402	!
1403	jtf = min0(jte,jde-1)
1404	ktf = min0(kte,kde-1)
1405	itf = min0(ite,ide-1)
1406	!
1407	if(.not.restart)then
1408	do j = jts,jtf
1409	do k = kts,ktf
1410	do i = its,itf
1411	rublten(i,k,j) = 0.
1412	rvblten(i,k,j) = 0.
1413	rthblten(i,k,j) = 0.
1414	rqvblten(i,k,j) = 0.
1415	rqcblten(i,k,j) = 0.
1416	te_temf(i,k,j) = TEmin
1417	cf3d_temf(i,k,j) = 0.
1418	enddo
1419	enddo
1420	enddo
1421	endif
1422	!
1423	if (p_qi .ge. p_first_scalar .and. .not.restart) then
1424	do j = jts,jtf
1425	do k = kts,ktf
1426	do i = its,itf
1427	rqiblten(i,k,j) = 0.
1428	enddo
1429	enddo
1430	enddo
1431	endif
1432	!
1433	end subroutine temfinit
1434	!-------------------------------------------------------------------
1435	end module module_bl_temf

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

Download in other formats:

Original Format