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