| blob | f869115147751be878deca6ff27fef8555855735 |
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
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 * ftau0**2. / 2 / fth0**2.
376 real, parameter :: PrT0 = Cphi/Ceps * ftau0**2 / 2. / fth0**2
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
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.
476 do i = its,ite ! Main loop
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)
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( -xlv*rl(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
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
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
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)))
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.608*thetav(i,1)*qf_temfx(i,1)
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)
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.
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
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)
565 ! Find h0 (should eventually be interpolated for smoothness)
566 if (thetav(i,k) > thetav(i,1) .AND. h0idx(i) .EQ. 1) then
567 ! WA 9/28/11 limit h0 as for hd and hct
568 if (zm(i,k) < hmax) then
569 h0idx(i) = k
570 h0(i) = zm(i,k)
571 else
572 h0idx(i) = k
573 h0(i) = hmax
574 end if
575 end if
576 end do
578 ! Gradients at top level
580 dthdz(i,kte) = dthdz(i,kte-1)
581 dqtdz(i,kte) = dqtdz(i,kte-1)
582 dudz(i,kte) = dudz(i,kte-1)
583 dvdz(i,kte) = dvdz(i,kte-1)
585 if ( hfx(i) > 0.) then
586 ! wstr(i) = (g * h0(i) / thetav(i,2) * shf_temfx(i,1) ) ** (1./3.)
587 wstr(i) = (g * h0(i) / thetav(i,2) * hfx(i)/(rho(i,1)*cp) ) ** (1./3.)
588 else
589 wstr(i) = 0.
590 end if
593 ! Set flag convective or not for use below
594 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
596 ! Find stability parameters and length scale (on turbulence levels)
597 do kt = 1,kte-1
598 N2(i,kt) = 2. * g / (thetav(i,kt) + thetav(i,kt+1))*dthdz(i,kt)
599 S(i,kt) = sqrt(dudz(i,kt)**2. + dvdz(i,kt)**2.)
600 Ri(i,kt) = N2(i,kt) / S(i,kt)**2.
601 if (S(i,kt) < 1e-15) then
602 if (N2(i,kt) >= 0) then
603 Ri(i,kt) = 10.
604 else
605 Ri(i,kt) = -1.
606 end if
607 end if
608 beta(i,kt) = 2. * g / (thetav(i,kt)+thetav(i,kt+1))
609 if (Ri(i,kt) > 0) then
610 ratio(i,kt) = Ri(i,kt)/(Cphi**2.*ftau0**2./(2.*Ceps**2.*fth0**2.)+3.*Ri(i,kt))
611 ftau(i,kt) = ftau0 * ((3./4.) / (1.+4.*Ri(i,kt)) + 1./4.)
612 fth(i,kt) = fth0 / (1.+4.*Ri(i,kt))
613 TE2(i,kt) = 2. * te_temfx(i,kt) * ratio(i,kt) * N2(i,kt) / beta(i,kt)**2.
614 else
615 ratio(i,kt) = Ri(i,kt)/(Cphi**2.*ftau0**2./(-2.*Ceps**2.*fth0**2.)+2.*Ri(i,kt))
616 ftau(i,kt) = ftau0
617 fth(i,kt) = fth0
618 TE2(i,kt) = 0.
619 end if
620 TKE(i,kt) = te_temfx(i,kt) * (1. - ratio(i,kt))
621 ustrtilde(i,kt) = sqrt(ftau(i,kt) * TKE(i,kt))
622 if (N2(i,kt) > 0.) then
623 linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / (Cf*ustrtilde(i,kt)) + sqrt(N2(i,kt))/(CN*ustrtilde(i,kt)) + 1./lasymp ! WA Test 11/20/09
624 else
625 linv(i,kt) = 1./karman / zt(i,kt) + abs(fCor(i)) / (Cf*ustrtilde(i,kt)) + 1./lasymp ! WA Test 11/20/09
626 end if
627 leps(i,kt) = 1./linv(i,kt)
628 leps(i,kt) = max(leps(i,kt),lepsmin(i,kt))
629 end do
630 S(i,kte) = 0.0
631 N2(i,kte) = 0.0
632 TKE(i,kte) = 0.0
633 linv(i,kte) = linv(i,kte-1)
634 leps(i,kte) = leps(i,kte-1)
637 ! Find diffusion coefficients
638 ! First use basic formulae for stable and neutral cases,
639 ! then for convective conditions, and finally choose the larger
640 ! WA 12/23/09 use convective form up to hd/2 always
641 ! WA 12/28/09 after restructuring, this block is above MF block,
642 ! so hd is not yet available for this timestep, must use h0,
643 ! or use hd from previous timestep but be careful about initialization.
644 do kt = 1,kte-1 ! WA 12/22/09
645 ! WA 4/8/10 remove beta term to avoid negative and huge values
646 ! of km due to very small denominator. This is an interim fix
647 ! until we find something better (more theoretically sound).
648 ! 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))
649 km(i,kt) = TKE(i,kt)**1.5 * ftau(i,kt)**2. / (Ceps * sqrt(TKE(i,kt)*te_temfx(i,kt)) / leps(i,kt))
650 kh(i,kt) = 2. * leps(i,kt) * fth(i,kt)**2. * TKE(i,kt) / sqrt(te_temfx(i,kt)) / Cphi
651 if ( is_convective) then
652 if (kt <= h0idx(i)) then
653 lconv(i,kt) = 1. / (1. / (karman*zt(i,kt)) + Cc / (karman * (h0(i) - zt(i,kt))))
654 else
655 lconv(i,kt) = 0.
656 end if
657 ! WA 12/15/09 use appropriate coeffs to match kh_conv and kh at neutral
658 kh_conv(i,kt) = ftau0**2. / Ceps / PrT0 * sqrt(TKE(i,kt)) * lconv(i,kt)
659 if (kh_conv(i,kt) < 0.) then
660 kh_conv(i,kt) = 0.
661 end if
662 km_conv(i,kt) = PrT0 * kh_conv(i,kt)
663 if (zt(i,kt) <= h0(i)/2.) then
664 km(i,kt) = km_conv(i,kt)
665 kh(i,kt) = kh_conv(i,kt)
666 end if
667 ! WA TEST 1/11/10 go back to max in upper BL
668 if (zt(i,kt) > h0(i)/2. .AND. kt <= h0idx(i)) then
669 km(i,kt) = max(km(i,kt),km_conv(i,kt),visc_temf)
670 kh(i,kt) = max(kh(i,kt),kh_conv(i,kt),conduc_temf)
671 end if
672 end if ! is_convective
673 km(i,kt) = max(km(i,kt),visc_temf)
674 kh(i,kt) = max(kh(i,kt),conduc_temf)
675 Fz(i,kt) = -kh(i,kt) * dthdz(i,kt) ! Diffusive heat flux
676 end do
677 km(i,kte) = km(i,kte-1) ! WA 12/22/09
678 kh(i,kte) = kh(i,kte-1)
679 Fz(i,kte) = 0.0 ! WA 4/2/10
682 !*** Mass flux block starts here ***
684 if ( is_convective) then
686 Cepsmf = 2. / max(200.,h0(i))
687 Cepsmf = max(Cepsmf,0.002) ! WA 7/20/10
688 do k = kts,kte
689 ! Calculate lateral entrainment fraction for subcloud layer
690 ! epsilon and delta are defined on mass grid (half levels)
691 epsmf(i,k) = Cepsmf
692 end do
694 ! Initialize updraft
695 thup_temfx(i,1) = thetal(i,1) ! No excess
696 qtup_temfx(i,1) = qt(i,1) ! No excess
697 rUPD(i,1) = qtup_temfx(i,1) / (1. - qtup_temfx(i,1))
698 wupd_temfx(i,1) = Cw * wstr(i)
699 wupd_dry(i,1) = Cw * wstr(i)
700 UUPD(i,1) = u_temf(i,1)
701 VUPD(i,1) = v_temf(i,1)
702 thetavUPD(i,1) = thup_temfx(i,1) * (1. + 0.608*qtup_temfx(i,1)) ! WA Assumes no liquid
703 thetavUPDmoist(i,1) = thup_temfx(i,1) * (1. + 0.608*qtup_temfx(i,1)) ! WA Assumes no liquid
704 TEUPD(i,1) = te_temfx(i,1) + g / thetav(i,1) * sfcTHVF(i)
705 ! qlUPD(i,1) = 0.
706 qlUPD(i,1) = qcx(i,1) ! WA allow environment liquid
707 TUPD(i,1) = thup_temfx(i,1) * pi2d(i,1)
708 rstUPD(i,1) = rsat(p2d(i,1),TUPD(i,1),ep2)
709 rlUPD(i,1) = 0.
711 ! Calculate updraft parameters counting up
712 do k = 2,kte
713 dthUPDdz(i,k-1) = -epsmf(i,k) * (thup_temfx(i,k-1) - thetal(i,k-1))
714 thup_temfx(i,k) = thup_temfx(i,k-1) + dthUPDdz(i,k-1) * dzm(i,k-1)
715 dqtup_temfxdz(i,k-1) = -epsmf(i,k) * (qtup_temfx(i,k-1) - qt(i,k-1))
716 qtup_temfx(i,k) = qtup_temfx(i,k-1) + dqtup_temfxdz(i,k-1) * dzm(i,k-1)
717 thetavUPD(i,k) = thup_temfx(i,k) * (1. + 0.608*qtup_temfx(i,k)) ! WA Assumes no liquid
718 B(i,k-1) = g * (thetavUPD(i,k) - thetav(i,k)) / thetav(i,k)
719 if ( wupd_dry(i,k-1) < 1e-15 ) then
720 wupd_dry(i,k) = 0.
721 else
722 dwUPDdz(i,k-1) = -2. *epsmf(i,k)*wupd_dry(i,k-1) + 0.33*B(i,k-1)/wupd_dry(i,k-1)
723 wupd_dry(i,k) = wupd_dry(i,k-1) + dwUPDdz(i,k-1) * dzm(i,k-1)
724 end if
725 dUUPDdz(i,k-1) = -epsmf(i,k) * (UUPD(i,k-1) - u_temf(i,k-1))
726 UUPD(i,k) = UUPD(i,k-1) + dUUPDdz(i,k-1) * dzm(i,k-1)
727 dVUPDdz(i,k-1) = -epsmf(i,k) * (VUPD(i,k-1) - v_temf(i,k-1))
728 VUPD(i,k) = VUPD(i,k-1) + dVUPDdz(i,k-1) * dzm(i,k-1)
729 dTEUPDdz(i,k-1) = -epsmf(i,k) * (TEUPD(i,k-1) - te_temfx(i,k-1))
730 TEUPD(i,k) = TEUPD(i,k-1) + dTEUPDdz(i,k-1) * dzm(i,k-1)
731 ! Alternative updraft velocity based on moist thetav
732 ! Need thetavUPDmoist, qlUPD
733 rUPD(i,k) = qtup_temfx(i,k) / (1. - qtup_temfx(i,k))
734 ! WA Updraft temperature assuming no liquid
735 TUPD(i,k) = thup_temfx(i,k) * pi2d(i,k)
736 ! Updraft saturation mixing ratio
737 ! rstUPD(i,k) = rsat(p2d(i,k),TUPD(i,k),ep2) ! WA 4/19/10
738 rstUPD(i,k) = rsat(p2d(i,k-1),TUPD(i,k),ep2)
739 ! Correct to actual temperature (Sommeria & Deardorff 1977)
740 beta1(i,k) = 0.622 * (xlv/(r_d*TUPD(i,k))) * (xlv/(cp*TUPD(i,k)))
741 rstUPD(i,k) = rstUPD(i,k) * (1.0+beta1(i,k)*rUPD(i,k)) / (1.0+beta1(i,k)*rstUPD(i,k))
742 qstUPD(i,k) = rstUPD(i,k) / (1. + rstUPD(i,k))
743 if (rUPD(i,k) > rstUPD(i,k)) then
744 rlUPD(i,k) = rUPD(i,k) - rstUPD(i,k)
745 qlUPD(i,k) = rlUPD(i,k) / (1. + rlUPD(i,k))
746 thetavUPDmoist(i,k) = (thup_temfx(i,k) + ((xlv/cp)*qlUPD(i,k)/pi2d(i,k))) * (1. + 0.608*qstUPD(i,k) - qlUPD(i,k))
747 else
748 rlUPD(i,k) = 0.
749 ! qlUPD(i,k) = 0.
750 qlUPD(i,k) = qcx(i,k-1) ! WA 4/6/10 allow environment liquid
751 ! WA does this make sense? Should be covered above?
752 thetavUPDmoist(i,k) = thup_temfx(i,k) * (1. + 0.608*qtup_temfx(i,k))
753 end if
754 Bmoist(i,k-1) = g * (thetavUPDmoist(i,k) - thetav(i,k)) / thetav(i,k)
755 if ( wupd_temfx(i,k-1) < 1e-15 ) then
756 wupd_temfx(i,k) = 0.
757 else
758 dwUPDmoistdz(i,k-1) = -2. *epsmf(i,k)*wupd_temfx(i,k-1) + 0.33*Bmoist(i,k-1)/wupd_temfx(i,k-1)
759 wupd_temfx(i,k) = wupd_temfx(i,k-1) + dwUPDmoistdz(i,k-1) * dzm(i,k-1)
760 end if
761 end do
763 ! Find hd based on wUPD
764 if (wupd_dry(i,1) == 0.) then
765 hdidx(i) = 1
766 else
767 hdidx(i) = kte ! In case wUPD <= 0 not found
768 do k = 2,kte
769 ! if (wupd_dry(i,k) <= 0.) then
770 if (wupd_dry(i,k) <= 0. .OR. zm(i,k) > hmax) then ! WA Test
771 hdidx(i) = k
772 goto 100 ! FORTRAN made me do it!
773 end if
774 end do
775 end if
776 100 hd(i) = zm(i,hdidx(i))
777 ! 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
778 kpbl1d(i) = hdidx(i) ! WA 5/11/10 kpbl should be index
779 hpbl(i) = hd(i) ! WA 5/11/10 hpbl is height. Should still be replaced by something that works whether convective or not.
781 ! Find LCL, hct, and ht
782 lclidx(i) = kte ! In case LCL not found
783 do k = kts,kte
784 if (rUPD(i,k) > rstUPD(i,k)) then
785 lclidx(i) = k
786 goto 200
787 end if
788 end do
789 200 lcl(i) = zm(i,lclidx(i))
791 if (hd(i) > lcl(i)) then ! Forced cloud (at least) occurs
792 ! Find hct based on wUPDmoist
793 if (wupd_temfx(i,1) == 0.) then
794 hctidx(i) = 1
795 else
796 hctidx(i) = kte ! In case wUPD <= 0 not found
797 do k = 2,kte
798 if (wupd_temfx(i,k) <= 0. .OR. zm(i,k) > hmax) then ! WA Test
799 hctidx(i) = k
800 goto 300 ! FORTRAN made me do it!
801 end if
802 end do
803 end if
804 300 hct(i) = zm(i,hctidx(i))
805 if (hctidx(i) <= hdidx(i)+1) then ! No active cloud
806 hct(i) = hd(i)
807 hctidx(i) = hdidx(i)
808 else
809 end if
810 else ! No cloud
811 hct(i) = hd(i)
812 hctidx(i) = hdidx(i)
813 end if
814 ht(i) = max(hd(i),hct(i))
815 htidx(i) = max(hdidx(i),hctidx(i))
817 ! Now truncate updraft at ht with taper
818 do k = 1,kte
819 if (zm(i,k) < 0.9*ht(i)) then ! Below taper region
820 tval(i) = 1
821 else if (zm(i,k) >= 0.9*ht(i) .AND. zm(i,k) <= 1.0*ht(i)) then
822 ! Within taper region
823 tval(i) = 1. - ((zm(i,k) - 0.9*ht(i)) / (1.0*ht(i) - 0.9*ht(i)))
824 else ! Above taper region
825 tval(i) = 0.
826 end if
827 thup_temfx(i,k) = tval(i) * thup_temfx(i,k) + (1-tval(i))*thetal(i,k)
828 thetavUPD(i,k) = tval(i) * thetavUPD(i,k) + (1-tval(i))*thetav(i,k)
829 qtup_temfx(i,k) = tval(i) * qtup_temfx(i,k) + (1-tval(i)) * qt(i,k)
830 qlUPD(i,k) = tval(i) * qlUPD(i,k) + (1-tval(i)) * qcx(i,k-1)
831 UUPD(i,k) = tval(i) * UUPD(i,k) + (1-tval(i)) * u_temf(i,k)
832 VUPD(i,k) = tval(i) * VUPD(i,k) + (1-tval(i)) * v_temf(i,k)
833 TEUPD(i,k) = tval(i) * TEUPD(i,k) + (1-tval(i)) * te_temfx(i,k)
834 if (zm(i,k) > ht(i)) then ! WA this is just for cleanliness
835 wupd_temfx(i,k) = 0.
836 dwUPDmoistdz(i,k) = 0.
837 wupd_dry(i,k) = 0.
838 dwUPDdz(i,k) = 0.
839 end if
840 end do
842 ! Calculate lateral detrainment rate for cloud layer
843 deltmf(i,1) = Cepsmf
844 do k = 2,kte-1
845 if (hctidx(i) > hdidx(i)+1) then ! Some cloud
846 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
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) + (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)))
875 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)))
876 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)))
877 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)))
878 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)))
879 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
880 do kt = 2,kte-1
881 MFCth(i,kt) = mf_temfx(i,kt) * (thup_temfx(i,kt)-thetal(i,kt) + thup_temfx(i,kt+1)-thetal(i,kt+1)) / 2.
882 MFCq(i,kt) = mf_temfx(i,kt) * (qtup_temfx(i,kt)-qt(i,kt) + qtup_temfx(i,kt+1)-qt(i,kt+1)) / 2.
883 MFCu(i,kt) = mf_temfx(i,kt) * (UUPD(i,kt)-u_temf(i,kt) + UUPD(i,kt+1)-u_temf(i,kt+1)) / 2.
884 MFCv(i,kt) = mf_temfx(i,kt) * (VUPD(i,kt)-v_temf(i,kt) + VUPD(i,kt+1)-v_temf(i,kt+1)) / 2.
885 MFCql(i,kt) = mf_temfx(i,kt) * (qlUPD(i,kt)-qcx(i,kt-1) + qlUPD(i,kt+1)-qcx(i,kt)) / 2.
886 MFCTE(i,kt) = mf_temfx(i,kt) * (TEUPD(i,kt)-te_temfx(i,kt)) ! TE is on turb levels
887 end do
888 MFCth(i,kte) = 0
889 MFCq(i,kte) = 0
890 MFCu(i,kte) = 0
891 MFCv(i,kte) = 0
892 MFCql(i,kte) = 0
893 MFCTE(i,kte) = 0
895 ! Calculate cloud fraction (on mass levels)
896 cf3d_temfx(i,1) = 0.0
897 cfm_temfx(i) = 0.0
898 do k = 2,kte
899 ! 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
900 if (wupd_temfx(i,k-1) >= 1.0e-15 .AND. wupd_temfx(i,k) >= 1.0e-15) then
901 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?
902 else
903 au(i,k) = 0.0
904 end if
905 sigq2 = au(i,k) * (qtup_temfx(i,k)-qt(i,k))
906 if (sigq2 > 0.0) then
907 sigq(i,k) = sqrt(sigq2)
908 else
909 sigq(i,k) = 0.0
910 end if
911 ! rst = rsat(p2d(i,k),thx(i,k)*pi2d(i,k),ep2)
912 rst = rsat(p2d(i,k-1),thx(i,k-1)*pi2d(i,k-1),ep2)
913 qst(i,k) = rst / (1. + rst)
914 satdef(i,k) = qt(i,k) - qst(i,k)
915 if (satdef(i,k) <= 0.0) then
916 if (sigq(i,k) > 1.0e-15) then
917 cf3d_temfx(i,k) = max(0.5 + 0.36 * atan(1.55*(satdef(i,k)/sigq(i,k))),0.0)
918 else
919 cf3d_temfx(i,k) = 0.0
920 end if
921 else
922 cf3d_temfx(i,k) = 1.0
923 end if
924 if (zm(i,k) < lcl(i)) then
925 cf3d_temfx(i,k) = 0.0
926 end if
927 ! Put max value so far into cfm
928 if (zt(i,k) <= hmax) then
929 cfm_temfx(i) = max(cf3d_temfx(i,k),cfm_temfx(i))
930 end if
931 end do
933 else ! not is_convective, no MF components
934 do kt = 1,kte
935 MFCth(i,kt) = 0
936 MFCq(i,kt) = 0
937 MFCu(i,kt) = 0
938 MFCv(i,kt) = 0
939 MFCql(i,kt) = 0
940 MFCTE(i,kt) = 0
941 end do
942 lcl(i) = zm(i,kte-1)
943 hct(i) = zm(i,1)
944 hctidx(i) = 1
945 hd(i) = zm(i,1)
946 hdidx(i) = 1
947 ht(i) = hd(i)
948 ! Cloud fraction calculations
949 cf3d_temfx(i,1) = 0.0
950 cfm_temfx(i) = 0.0
951 do k = 2,kte
952 if (qcx(i,k-1) > 1.0e-15) then
953 cf3d_temfx(i,k) = 1.0
954 else
955 cf3d_temfx(i,k) = 0.0
956 end if
957 ! Put max value so far into cfm
958 if (zt(i,k) <= hmax) then
959 cfm_temfx(i) = max(cf3d_temfx(i,k),cfm_temfx(i))
960 end if
961 end do
963 end if ! MF components or not
964 cf3d_temfx(i,kte) = 0.0
965 ! Mass flux block ends here
967 ! Flux profiles
968 do kt = 2,kte
969 ! Fz(i,kt) = -kh(i,kt) * dthdz(i,kt)
970 shf_temfx(i,kt) = Fz(i,kt) + MFCth(i,kt)
971 QFK(i,kt) = -kh(i,kt) * dqtdz(i,kt)
972 qf_temfx(i,kt) = QFK(i,kt) + MFCq(i,kt)
973 uwk(i,kt) = -km(i,kt) * dudz(i,kt)
974 uw_temfx(i,kt) = uwk(i,kt) + MFCu(i,kt)
975 vwk(i,kt) = -km(i,kt) * dvdz(i,kt)
976 vw_temfx(i,kt) = vwk(i,kt) + MFCv(i,kt)
977 end do
979 ! Surface momentum fluxes
980 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)
981 ang(i) = atan2(v_temf(i,2),u_temf(i,2))
982 uw_temfx(i,1) = -cos(ang(i)) * ust(i)**2.
983 vw_temfx(i,1) = -sin(ang(i)) * ust(i)**2.
985 ! Calculate mixed scaling velocity (Moeng & Sullivan 1994 JAS p.1021)
986 ! Replaces ust everywhere (WA need to reconsider?)
987 ! wm(i) = (1./5. * (wstr(i)**3. + 5. * ust(i)**3.)) ** (1./3.)
988 ! WA TEST (R2,R11) 7/23/10 reduce velocity scale to fix excessive fluxes
989 wm(i) = 0.5 * (1./5. * (wstr(i)**3. + 5. * ust(i)**3.)) ** (1./3.)
990 ! WA TEST 2/14/11 limit contribution of w*
991 ! wm(i) = 0.5 * (1./5. * (min(0.8,wstr(i))**3. + 5. * ust(i)**3.)) ** (1./3.)
992 ! WA TEST (R3-R11) 7/23/10 wm = u*
993 ! wm(i) = ust(i)
995 ! Specified flux versions (flux is modified by land surface)
996 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))
997 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))
998 Fz(i,1) = shf_temfx(i,1) - MFCth(i,1)
999 QFK(i,1) = qf_temfx(i,1) - MFCq(i,1)
1002 ! Calculate thetav and its flux
1003 ! From Lewellen 2004 eq. 3
1004 ! WA The constants and combinations below should probably be replaced
1005 ! by something more consistent with the WRF system, but for now
1006 ! I don't want to take the chance of breaking something.
1007 do kt = 2,kte-1
1008 alpha2(i,kt) = 0.61 * (thetal(i,kt) + thetal(i,kt+1)) / 2.
1009 beta2(i,kt) = (100000. / p2di(i,kt))**0.286 * 2.45e-6 / 1004.67 - 1.61 * (thetal(i,kt) + thetal(i,kt+1)) / 2.
1010 end do
1011 alpha2(i,1) = 0.61 * (thetal(i,1) + (thetal(i,2)-thetal(i,1)) * (zt(i,2) - znt(i)) / (zm(i,2) - znt(i)))
1012 alpha2(i,kte) = 0.61 * thetal(i,kte)
1013 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)))
1014 beta2(i,kte) = (100000. / p2di(i,kte))**0.286 * 2.45e-6 / 1004.67 - 1.61 * thetal(i,kte)
1015 if ( is_convective ) then ! Activate EDMF components
1016 do kt = 1,kte-1
1017 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)
1018 end do
1019 MFCthv(i,kte) = 0.
1020 else ! No MF components
1021 do kt = 1,kte
1022 MFCthv(i,kt) = 0.
1023 end do
1024 end if
1026 do kt = 1,kte
1027 THVF(i,kt) = (1. + 0.61 * qt(i,kt)) * Fz(i,kt) + alpha2(i,kt) * QFK(i,kt) + MFCthv(i,kt)
1028 end do
1030 ! Update mean variables:
1031 ! This is done with implicit solver for diffusion part followed
1032 ! by explicit solution for MF terms.
1033 ! Note that Coriolis terms that were source terms for u and v
1034 ! in Matlab code are now handled by WRF outside this PBL context.
1036 u_new(i,:) = u_temf(i,:)
1037 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.)
1038 do k = 2,kte-1
1039 u_new(i,k) = u_new(i,k) + dt * (-(MFCu(i,k)-MFCu(i,k-1))) / dzm(i,k)
1040 end do
1042 v_new(i,:) = v_temf(i,:)
1043 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.)
1044 do k = 2,kte-1
1045 v_new(i,k) = v_new(i,k) + dt * (-(MFCv(i,k)-MFCv(i,k-1))) / dzm(i,k)
1046 end do
1048 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.)
1049 do k = 2,kte-1
1050 thetal(i,k) = thetal(i,k) + dt * (-(MFCth(i,k)-MFCth(i,k-1))) / dzm(i,k)
1051 end do
1053 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.)
1054 do k = 2,kte-1
1055 qt(i,k) = qt(i,k) + dt * (-(MFCq(i,k)-MFCq(i,k-1))) / dzm(i,k)
1056 end do
1058 ! Stepping TE forward is a bit more complicated
1059 te_temfx(i,1) = ust(i)**2. / ftau(i,1) * (1. + ratio(i,1))
1060 if ( is_convective ) then
1061 ! WA currently disabled if MFopt=false, is that right?
1062 convection_TKE_surface_src(i) = 2. * beta(i,1) * shf_temfx(i,1)
1063 else
1064 convection_TKE_surface_src(i) = 0.
1065 end if
1066 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.))
1067 if (te_temfx(i,1) > 20.0) then
1068 te_temfx(i,1) = 20.0 ! WA 9/28/11 limit max TE
1069 end if
1070 sfcFTE(i) = -(km(i,1)+km(i,2)) / 2. * (te_temfx(i,2)-te_temfx(i,1)) / dzm(i,2)
1072 do kt = 1,kte
1073 if (THVF(i,kt) >= 0) then
1074 buoy_src(i,kt) = 2. * g / thetav(i,kt) * THVF(i,kt)
1075 else
1076 buoy_src(i,kt) = 0. ! Cancel buoyancy term when locally stable
1077 end if
1078 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)
1079 end do
1080 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.)
1081 do kt = 2,kte-1
1082 te_temfx(i,kt) = te_temfx(i,kt) + dt * srcs(i,kt)
1083 te_temfx(i,kt) = te_temfx(i,kt) + dt * (-(MFCTE(i,kt)-MFCTE(i,kt-1))) / dzt(i,kt)
1084 if (te_temfx(i,kt) < TEmin) te_temfx(i,kt) = TEmin
1085 end do
1086 te_temfx(i,kte) = 0.0 ! WA 4/2/10
1087 do kt = 2,kte-1
1088 if (te_temfx(i,kt) > 20.0) then
1089 te_temfx(i,kt) = 20.0 ! WA 9/29/11 reduce limit max TE from 30
1090 end if
1091 end do
1093 ! Done with updates, now convert internal variables back to WRF vars
1094 do k = kts,kte
1095 ! Populate kh_temfx, km_temfx from kh, km
1096 kh_temfx(i,k) = kh(i,k)
1097 km_temfx(i,k) = km(i,k)
1098 end do
1100 ! Convert thetal, qt back to theta, qv, qc
1101 ! See opposite conversion at top of subroutine
1102 ! WA this accounts for offset of indexing between
1103 ! WRF and TEMF, see notes at top of this routine.
1104 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)
1106 do k = kts,kte-1
1107 ! Calculate tendency terms
1108 ! WA this accounts for offset of indexing between
1109 ! WRF and TEMF, see notes at top of this routine.
1110 rubltenx(i,k) = (u_new(i,k+1) - u_temf(i,k+1)) / dt
1111 rvbltenx(i,k) = (v_new(i,k+1) - v_temf(i,k+1)) / dt
1112 rthbltenx(i,k) = (thx_new(i,k) - thx(i,k)) / dt
1113 rqvbltenx(i,k) = (qvx_new(i,k) - qvx(i,k)) / dt
1114 rqcbltenx(i,k) = (qcx_new(i,k) - qcx(i,k)) / dt
1115 end do
1116 rubltenx(i,kte) = 0.
1117 rvbltenx(i,kte) = 0.
1118 rthbltenx(i,kte) = 0.
1119 rqvbltenx(i,kte) = 0.
1120 rqcbltenx(i,kte) = 0.
1122 ! Populate surface exchange coefficient variables to go back out
1123 ! for next time step of surface scheme
1124 ! Unit specifications in SLAB and sfclay are conflicting and probably
1125 ! incorrect. This will give a dynamic heat flux (W/m^2) or moisture
1126 ! flux (kg(water)/(m^2*s)) when multiplied by a difference.
1127 ! These formulae are the same as what's used above to get surface
1128 ! flux from surface temperature and specific humidity.
1129 ! WA 2/16/11 removed, not needed in BL
1130 ! 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)
1131 ! flqc(i) = rho(i,1) * fth(i,1)/fth0 * wm(i) * leps(i,1) / PrT0 / log(zm(i,2)/z0t(i)) / zt(i,1)
1132 ! WA Must exchange coeffs be limited to avoid runaway in some
1133 ! (convective?) conditions? Something like this is done in sfclay.
1134 ! Doing nothing for now.
1136 ! Populate 10 m winds and 2 m temp
1137 ! WA Note this only works if first mass level is above 10 m
1138 u10(i) = u_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i))
1139 v10(i) = v_new(i,2) * log(10.0/znt(i)) / log(zm(i,2)/znt(i))
1140 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
1142 ! Populate diagnostic variables
1143 hd_temfx(i) = hd(i)
1144 lcl_temfx(i) = lcl(i)
1145 hct_temfx(i) = hct(i)
1147 ! Send updraft liquid water back
1148 if ( is_convective) then
1149 do k = kts,kte-1
1150 qlup_temfx(i,k) = qlUPD(i,k)
1151 end do
1152 else
1153 qlup_temfx(i,1) = qcx(i,1)
1154 do k = kts+1,kte-1
1155 qlup_temfx(i,k) = qcx(i,k-1)
1156 end do
1157 end if
1158 qlup_temfx(i,kte) = qcx(i,kte)
1160 end do ! Main (i) loop
1162 end subroutine temf2d
1163 !
1164 !--------------------------------------------------------------------
1165 !
1166 subroutine thlqt2thqvqc(thetal,qt,theta,qv,qc,p,piex,kbot,ktop,ep2,Lv,Cp)
1168 ! Calculates theta, qv, qc from thetal, qt.
1169 ! Originally from RAMS (subroutine satadjst) by way of Hongli Jiang.
1171 implicit none
1172 integer, intent(in ) :: kbot, ktop
1173 real, dimension( kbot:ktop ), intent(in ) :: thetal, qt
1174 real, dimension( kbot:ktop ), intent( out) :: theta, qv, qc
1175 real, dimension( kbot:ktop ), intent(in ) :: p, piex
1176 real, intent(in ) :: ep2, Lv, Cp
1177 !
1178 ! Local variables
1179 integer :: k, iterate
1180 real :: T1, Tt
1181 real, dimension( kbot:ktop) :: rst
1182 real, dimension( kbot:ktop) :: Tair, rc, rt, rv
1183 !
1184 do k = kbot,ktop
1185 T1 = thetal(k) * piex(k) ! First guess T is just thetal converted to T
1186 Tair(k) = T1
1187 rt(k) = qt(k) / (1. - qt(k))
1189 do iterate = 1,20
1190 rst(k) = rsat(p(k),Tair(k),ep2)
1191 rc(k) = max(rt(k) - rst(k), 0.)
1192 Tt = 0.7*Tair(k) + 0.3*T1 * (1.+Lv*rc(k) / (Cp*max(Tair(k),253.)))
1193 if ( abs(Tt - Tair(k)) < 0.001) GOTO 100
1194 Tair(k) = Tt
1195 end do
1196 100 continue
1197 rv(k) = rt(k) - rc(k)
1198 qv(k) = rv(k) / (1. + rv(k))
1199 qc(k) = rc(k) / (1. + rc(k))
1200 theta(k) = Tair(k) / piex(k)
1201 end do ! k loop
1202 return
1203 end subroutine thlqt2thqvqc
1204 !
1205 !--------------------------------------------------------------------
1206 !
1207 subroutine findhct_te( thetavenv,thetaparin,qpar, &
1208 rpar,hdidx,paridx,zm,hct,hctidx,p,piex,ep2,kbot,ktop)
1210 ! Calculates the cloud top height (limit of convection) for the
1211 ! updraft properties. hct is the level at which a parcel lifted
1212 ! at the moist adiabatic rate where it is saturated and at the dry
1213 ! adiabatic rate otherwise, first has thetav cooler than the environment.
1214 ! Loops start at LCL (paridx).
1216 implicit none
1217 integer, intent(in) :: kbot, ktop
1218 integer, intent(in) :: paridx, hdidx
1219 real, intent(in) :: ep2
1220 real, dimension( kbot:ktop), intent(in) :: thetavenv
1221 real, dimension( kbot:ktop), intent(in) :: thetaparin
1222 real, dimension( kbot:ktop), intent(in) :: qpar, rpar, zm, p, piex
1223 real, intent(out) :: hct
1224 integer, intent(out) :: hctidx
1225 ! Local variables
1226 integer k
1227 real, dimension( kbot:ktop) :: thetapar, thetavpar, qlpar, Tpar, rsatpar
1228 real, dimension( kbot:ktop) :: qsatpar
1229 real :: gammas, TparC
1231 thetapar(paridx) = thetaparin(paridx)
1232 Tpar(paridx) = thetapar(paridx) * piex(paridx)
1233 hctidx = ktop ! In case hct not found
1234 do k = paridx+1,ktop
1235 ! Find saturation mixing ratio at parcel temperature
1236 rsatpar(k) = rsat(p(k-1),Tpar(k-1),ep2)
1237 qsatpar(k) = rsatpar(k) / (1. + rsatpar(k))
1239 ! When parcel is unsaturated, its temperature changes
1240 ! at dry adiabatic rate (in other words, theta is constant).
1241 if (rpar(k) <= rsatpar(k)) then
1242 thetapar(k) = thetapar(k-1)
1243 Tpar(k) = thetapar(k) * piex(k)
1244 thetavpar(k) = thetapar(k) * (1.+0.608*qpar(k))
1245 else
1246 ! When parcel is saturated, its temperature changes at
1247 ! moist adiabatic rate
1248 ! Calculate moist adiabatic lapse rate according to Gill A4.12
1249 ! Requires T in deg.C
1250 TparC = Tpar(k-1) - 273.15
1251 gammas = 6.4 - 0.12 * TparC + 2.5e-5 * TparC**3. + (-2.4 + 1.e-3 * (TparC-5.)**2.) * (1. - p(k-1)/100000.)
1252 Tpar(k) = Tpar(k-1) - gammas/1000. * (zm(k)-zm(k-1))
1253 thetapar(k) = Tpar(k) / piex(k)
1254 qlpar(k) = qpar(k) - qsatpar(k) ! Liquid water in parcel
1255 thetavpar(k) = thetapar(k) * (1. + 0.608 * qsatpar(k) - qlpar(k))
1256 end if
1257 if (thetavenv(k) >= thetavpar(k)) then
1258 hctidx = k
1259 goto 1000
1260 end if
1261 end do
1262 1000 hct = zm(hctidx)
1263 return
1264 end subroutine findhct_te
1265 !
1266 !--------------------------------------------------------------------
1267 !
1268 real function rsat(p,T,ep2)
1270 ! Calculates the saturation mixing ratio with respect to liquid water
1271 ! Arguments are pressure (Pa) and absolute temperature (K)
1272 ! Uses the formula from the ARM intercomparison setup.
1273 ! Converted from Matlab by WA 7/28/08
1275 implicit none
1276 real p, T, ep2
1277 real temp, x
1278 real, parameter :: c0 = 0.6105851e+3
1279 real, parameter :: c1 = 0.4440316e+2
1280 real, parameter :: c2 = 0.1430341e+1
1281 real, parameter :: c3 = 0.2641412e-1
1282 real, parameter :: c4 = 0.2995057e-3
1283 real, parameter :: c5 = 0.2031998e-5
1284 real, parameter :: c6 = 0.6936113e-8
1285 real, parameter :: c7 = 0.2564861e-11
1286 real, parameter :: c8 = -0.3704404e-13
1288 temp = T - 273.15
1290 x =c0+temp*(c1+temp*(c2+temp*(c3+temp*(c4+temp*(c5+temp*(c6+temp*(c7+temp*c8)))))))
1291 rsat = ep2*x/(p-x)
1293 return
1294 end function rsat
1295 !
1296 !--------------------------------------------------------------------
1297 !
1298 subroutine solve_implicit_temf(Khlf,psi_n,srf_flux,dzm,dzt,kbot,ktop,dt,print_flag)
1300 ! Implicit solution of vertical diffusion for conserved variable
1301 ! psi given diffusivity Khlf on turbulence levels,
1302 ! and surface flux srf_flux.
1303 ! dzm is delta_z of mass levels, dzt is delta_z of turbulence levels.
1304 ! dt is timestep (s).
1306 implicit none
1307 integer :: kbot, ktop
1308 logical :: print_flag
1309 real :: srf_flux, dt
1310 real, dimension( kbot:ktop ), intent(in ) :: Khlf
1311 real, dimension( kbot:ktop ), intent(in ) :: dzm, dzt
1312 real, dimension( kbot:ktop ), intent(inout) :: psi_n
1313 !
1314 ! Local variables
1315 integer :: k
1316 real, dimension( kbot:ktop ) :: AU, BU, CU, YU
1317 !
1318 AU(kbot) = Khlf(kbot) / (dzm(kbot)*dzt(kbot))
1319 BU(kbot) = -1.0/dt - Khlf(kbot+1)/(dzm(kbot+1)*dzt(kbot+1))
1320 CU(kbot) = Khlf(kbot+1)/(dzm(kbot)*dzt(kbot+1))
1321 YU(kbot) = -psi_n(kbot)/dt - srf_flux/dzm(kbot)
1323 do k = kbot+1,ktop-1
1324 ! Subdiagonal (A) vector
1325 AU(k) = Khlf(k) / (dzm(k) * dzt(k))
1326 ! Main diagonal (B) vector
1327 BU(k) = -1.0/dt - (Khlf(k)/dzt(k) + Khlf(k+1)/dzt(k+1)) / dzm(k)
1328 ! Superdiagonal (C) vector
1329 CU(k) = Khlf(k+1) / (dzm(k)*dzt(k+1))
1330 ! Result vector
1331 YU(k) = -psi_n(k)/dt
1332 end do ! k loop
1334 AU(ktop) = 0.
1335 BU(ktop) = -1.0 / dt
1336 YU(ktop) = -psi_n(ktop) / dt
1338 ! Compute result with tridiagonal routine
1339 psi_n = trid(AU,BU,CU,YU,kbot,ktop)
1341 return
1342 end subroutine solve_implicit_temf
1343 !
1344 !--------------------------------------------------------------------
1345 !
1346 function trid(a,b,c,r,kbot,ktop)
1348 ! Solves tridiagonal linear system.
1349 ! Inputs are subdiagonal vector a, main diagonal b, superdiagonal c,
1350 ! result r, column top and bottom indices kbot and ktop.
1351 ! Originally from Numerical Recipes section 2.4.
1353 implicit none
1354 real, dimension( kbot:ktop ) :: trid
1355 integer :: kbot, ktop
1356 real, dimension( kbot:ktop ), intent(in ) :: a, b, c, r
1357 !
1358 ! Local variables
1359 integer :: k
1360 real :: bet
1361 real, dimension( kbot:ktop ) :: gam, u
1362 !
1363 bet = b(kbot)
1364 u(kbot) = r(kbot) / bet
1366 do k = kbot+1,ktop
1367 gam(k) = c(k-1) / bet
1368 bet = b(k) - a(k)*gam(k)
1369 u(k) = (r(k) - a(k)*u(k-1)) / bet
1370 end do
1372 do k = ktop-1,kbot,-1
1373 u(k) = u(k) - gam(k+1)*u(k+1)
1374 end do
1376 trid = u
1378 return
1379 end function trid
1380 !
1381 !--------------------------------------------------------------------
1382 !
1383 subroutine temfinit(rublten,rvblten,rthblten,rqvblten, &
1384 rqcblten,rqiblten,p_qi,p_first_scalar, &
1385 restart, allowed_to_read, &
1386 te_temf, cf3d_temf, &
1387 ids, ide, jds, jde, kds, kde, &
1388 ims, ime, jms, jme, kms, kme, &
1389 its, ite, jts, jte, kts, kte )
1390 !-------------------------------------------------------------------
1391 implicit none
1392 !-------------------------------------------------------------------
1393 !
1394 logical , intent(in) :: restart, allowed_to_read
1395 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1396 ims, ime, jms, jme, kms, kme, &
1397 its, ite, jts, jte, kts, kte
1398 integer , intent(in) :: p_qi,p_first_scalar
1399 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
1400 rublten, &
1401 rvblten, &
1402 rthblten, &
1403 rqvblten, &
1404 rqcblten, &
1405 rqiblten, &
1406 te_temf, &
1407 cf3d_temf
1408 ! Local variables
1409 integer :: i, j, k, itf, jtf, ktf
1410 real, parameter :: TEmin = 1e-3
1411 !
1412 jtf = min0(jte,jde-1)
1413 ktf = min0(kte,kde-1)
1414 itf = min0(ite,ide-1)
1415 !
1416 if(.not.restart)then
1417 do j = jts,jtf
1418 do k = kts,ktf
1419 do i = its,itf
1420 rublten(i,k,j) = 0.
1421 rvblten(i,k,j) = 0.
1422 rthblten(i,k,j) = 0.
1423 rqvblten(i,k,j) = 0.
1424 rqcblten(i,k,j) = 0.
1425 te_temf(i,k,j) = TEmin
1426 cf3d_temf(i,k,j) = 0.
1427 enddo
1428 enddo
1429 enddo
1430 endif
1431 !
1432 if (p_qi .ge. p_first_scalar .and. .not.restart) then
1433 do j = jts,jtf
1434 do k = kts,ktf
1435 do i = its,itf
1436 rqiblten(i,k,j) = 0.
1437 enddo
1438 enddo
1439 enddo
1440 endif
1441 !
1442 end subroutine temfinit
1443 !-------------------------------------------------------------------
1444 end module module_bl_temf