CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
shoc.F90
1
3
5module shoc
6 use machine, only: kind_phys
7
8 implicit none
9
10 private
11
12 public shoc_run, shoc_init
13 integer, parameter :: kp = kind_phys
14
15contains
16
17subroutine shoc_init (do_shoc, errmsg, errflg)
18 implicit none
19 logical, intent(in) :: do_shoc
20 character(len=*), intent(out) :: errmsg
21 integer, intent(out) :: errflg
22
23! Initialize CCPP error handling variables
24 errmsg = ''
25 errflg = 0
26
27! Consistency checks
28 if (.not. do_shoc) then
29 errflg = 1
30 write(errmsg,'(*(a))') 'Logic error: do_shoc == .false.'
31 return
32 end if
33end subroutine shoc_init
34
38subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, &
39 con_pi, con_fvirt, con_eps, dtp, prsl, delp, phii, phil, u, v, omega, rhc, &
40 supice, pcrit, cefac, cesfac, tkef1, dis_opt, hflx, evap, prnum, &
41 gt0, gq0, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc, &
42 cld_sgs, tke, tkh, wthv_sec, errmsg, errflg)
43
44 implicit none
45
46 integer, intent(in) :: nx, nzm, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc
47 real(kind=kind_phys), intent(in) :: tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, &
48 con_rd, con_pi, con_fvirt, con_eps, &
49 dtp, supice, pcrit, cefac, cesfac, tkef1, dis_opt
50 !
51 real(kind=kind_phys), intent(in), dimension(:) :: hflx, evap
52 real(kind=kind_phys), intent(in), dimension(:,:) :: prsl, delp, phil, u, v, omega, rhc, prnum
53 real(kind=kind_phys), intent(in), dimension(:,:) :: phii
54 !
55 real(kind=kind_phys), intent(inout), dimension(:,:) :: gt0, tke, tkh, wthv_sec
56 real(kind=kind_phys), intent(inout), dimension(:,:), optional :: cld_sgs
57 real(kind=kind_phys), intent(inout), dimension(:,:,:) :: gq0
58
59 character(len=*), intent(out) :: errmsg
60 integer, intent(out) :: errflg
61
62 real(kind=kind_phys), parameter :: epsq = 1.0e-20_kp, zero=0.0_kp, one=1.0_kp
63
64 integer :: i, k
65
66 real(kind=kind_phys) :: tem
67 real(kind=kind_phys), dimension(nx,nzm) :: qi ! local array of suspended cloud ice
68 real(kind=kind_phys), dimension(nx,nzm) :: qc ! local array of suspended cloud water
69 real(kind=kind_phys), dimension(nx,nzm) :: qsnw ! local array of suspended snowq
70 real(kind=kind_phys), dimension(nx,nzm) :: qrn ! local array of suepended rain
71 real(kind=kind_phys), dimension(nx,nzm) :: qgl ! local array of suspended graupel
72 real(kind=kind_phys), dimension(nx,nzm) :: ncpl ! local array of cloud water number concentration
73 real(kind=kind_phys), dimension(nx,nzm) :: ncpi ! local array of cloud ice number concentration
74
75! Initialize CCPP error handling variables
76
77 errmsg = ''
78 errflg = 0
79
80 if (ntiw < 0) then ! this is valid only for Zhao-Carr scheme
81 do k=1,nzm
82 do i=1,nx
83 qc(i,k) = gq0(i,k,ntcw)
84 if (abs(qc(i,k)) < epsq) then
85 qc(i,k) = zero
86 endif
87 tem = qc(i,k) * max(zero, min(one, (tcr-gt0(i,k))*tcrf))
88 qi(i,k) = tem ! ice
89 qc(i,k) = qc(i,k) - tem ! water
90 qrn(i,k) = zero
91 qsnw(i,k) = zero
92 ncpl(i,k) = zero
93 ncpi(i,k) = zero
94 enddo
95 enddo
96 else
97 if (ntgl > 0) then ! graupel exists - combine graupel with snow
98 do k=1,nzm
99 do i=1,nx
100 qc(i,k) = gq0(i,k,ntcw)
101 qi(i,k) = gq0(i,k,ntiw)
102 qrn(i,k) = gq0(i,k,ntrw)
103 qsnw(i,k) = gq0(i,k,ntsw) + gq0(i,k,ntgl)
104 enddo
105 enddo
106 else ! no graupel
107 do k=1,nzm
108 do i=1,nx
109 qc(i,k) = gq0(i,k,ntcw)
110 qi(i,k) = gq0(i,k,ntiw)
111 qrn(i,k) = gq0(i,k,ntrw)
112 qsnw(i,k) = gq0(i,k,ntsw)
113 enddo
114 enddo
115 endif
116 if (ntlnc > 0) then
117 do k=1,nzm
118 do i=1,nx
119 ncpl(i,k) = gq0(i,k,ntlnc)
120 ncpi(i,k) = gq0(i,k,ntinc)
121 enddo
122 enddo
123 endif
124 endif
125
126 ! phy_f3d(1,1,ntot3d-2) - shoc determined sgs clouds
127 ! phy_f3d(1,1,ntot3d-1) - shoc determined diffusion coefficients
128 ! phy_f3d(1,1,ntot3d ) - shoc determined w'theta'
129
130 call shoc_work (nx, nx, nzm, nzm+1, dtp, prsl, delp, &
131 phii, phil, u, v, omega, gt0, gq0(:,:,1), qi, qc, qsnw, qrn, &
132 rhc, supice, pcrit, cefac, cesfac, tkef1, dis_opt, &
133 cld_sgs, tke, hflx, evap, prnum, tkh, wthv_sec, &
134 ntlnc, ncpl, ncpi, &
135 con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, &
136 con_fvirt, con_eps)
137
138 if (ntiw < 0) then ! this is valid only for Zhao-Carr scheme
139 do k=1,nzm
140 do i=1,nx
141 gq0(i,k,ntcw) = qc(i,k) + qi(i,k)
142 enddo
143 enddo
144 else
145 do k=1,nzm
146 do i=1,nx
147 gq0(i,k,ntcw) = qc(i,k)
148 gq0(i,k,ntiw) = qi(i,k)
149 enddo
150 enddo
151 if (ntlnc > 0) then
152 do k=1,nzm
153 do i=1,nx
154 gq0(i,k,ntlnc) = ncpl(i,k)
155 gq0(i,k,ntinc) = ncpi(i,k)
156 enddo
157 enddo
158 endif
159 endif
160
161end subroutine shoc_run
162
163 ! Implementation of the Simplified High Order Closure (SHOC) scheme
164 ! of Bogenschutz and Krueger (2013), J. Adv. Model. Earth Syst, 5, 195-211,
165 ! doi: 10.1002/jame.200118. (further referred to as BK13)
166 ! in a single column form suitable for use in a GCM physics package.
167 ! Alex Belochitski, heavily based on the code of Peter Bogenschutz.
168 ! S Moorthi - optimization, cleanup, improve and customize for gsm
169 ! - improved solution for sgs-tke equation
170 ! S Moorthi - 05-11-17 - modified shear production term to eliminate
171 ! spurious tke ove Antarctica.
172 ! S Moorthi - 01-12-17 - added extra pressure dependent tke dissipation at
173 ! pressures below a critical value pcrit
174 ! S Moorthi - 04-12-17 - fixed a bug in the definition of hl on input
175 ! replacing fac_fus by fac_sub
176 ! S.Moorthi - 00-00-17 - added an alternate option for near boundary cek following
177 ! Scipion et. al., from U. Oklahoma.
178 subroutine shoc_work (ix, nx, nzm, nz, dtn, &
179 prsl, delp, phii, phil, u, v, omega, tabs, &
180 qwv, qi, qc, qpi, qpl, rhc, supice, &
181 pcrit, cefac, cesfac, tkef1, dis_opt, &
182 cld_sgs, tke, hflx, evap, prnum, tkh, &
183 wthv_sec, ntlnc, ncpl, ncpi, &
184 cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps)
185
186 use funcphys , only : fpvsl, fpvsi, fpvs ! saturation vapor pressure for water & ice
187
188 implicit none
189
190 real, intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps
191 integer, intent(in) :: ix ! max number of points in the physics window in the x
192 integer, intent(in) :: nx ! Number of points in the physics window in the x
193
194 integer, intent(in) :: nzm ! Number of vertical layers
195 integer, intent(in) :: nz ! Number of layer interfaces (= nzm + 1)
196 integer, intent(in) :: ntlnc ! index of liquid water number concentration
197 real, intent(in) :: dtn ! Physics time step, s
198
199 real, intent(in) :: pcrit ! pressure in Pa below which additional tke dissipation is applied
200 real, intent(in) :: cefac ! tunable multiplier to dissipation term
201 real, intent(in) :: cesfac ! tunable multiplier to dissipation term for bottom level
202 real, intent(in) :: tkef1 ! uncentering terms in implicit tke integration
203 real, intent(in) :: dis_opt ! when > 0 use different formula for near surface dissipation
204
205 real, intent(in) :: hflx(nx)
206 real, intent(in) :: evap(nx)
207
208! The interface is talored to GFS in a sense that input variables are 2D
209
210 real, intent(in) :: prsl (ix,nzm) ! mean layer presure
211 real, intent(in) :: delp (ix,nzm) ! layer presure depth
212 real, intent(in) :: phii (ix,nz ) ! interface geopotential height
213 real, intent(in) :: phil (ix,nzm) ! layer geopotential height
214 real, intent(in) :: u (ix,nzm) ! u-wind, m/s
215 real, intent(in) :: v (ix,nzm) ! v-wind, m/s
216 real, intent(in) :: omega (ix,nzm) ! omega, Pa/s
217 real, intent(inout) :: tabs (ix,nzm) ! temperature, K
218 real, intent(inout) :: qwv (ix,nzm) ! water vapor mixing ratio, kg/kg
219 real, intent(inout) :: qc (ix,nzm) ! cloud water mixing ratio, kg/kg
220 real, intent(inout) :: qi (ix,nzm) ! cloud ice mixing ratio, kg/kg
221! Anning Cheng 03/11/2016 SHOC feedback to number concentration
222 real, intent(inout) :: ncpl (nx,nzm) ! cloud water number concentration,/m^3
223 real, intent(inout) :: ncpi (nx,nzm) ! cloud ice number concentration,/m^3
224 real, intent(in) :: qpl (nx,nzm) ! rain mixing ratio, kg/kg
225 real, intent(in) :: qpi (nx,nzm) ! snow mixing ratio, kg/kg
226
227 real, intent(in) :: rhc (nx,nzm) ! critical relative humidity
228 real, intent(in) :: supice ! ice supersaturation parameter
229 real, intent(out) :: cld_sgs(ix,nzm) ! sgs cloud fraction
230! real, intent(inout) :: cld_sgs(nx,nzm) ! sgs cloud fraction
231 real, intent(inout) :: tke (ix,nzm) ! turbulent kinetic energy. m**2/s**2
232! real, intent(inout) :: tk (nx,nzm) ! eddy viscosity
233 real, intent(inout) :: tkh (ix,nzm) ! eddy diffusivity
234 real, intent(in) :: prnum (nx,nzm) ! turbulent Prandtl number
235 real, intent(inout) :: wthv_sec (ix,nzm) ! Buoyancy flux, K*m/s
236
237 real, parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, &
238 three=3.0_kp, oneb3=one/three, twoby3=two/three, fourb3=twoby3+twoby3
239 real, parameter :: sqrt2 = sqrt(two), twoby15 = two / 15.0_kp, &
240 nmin = 1.0_kp, ri_cub = 6.4e-14_kp, rl_cub = 1.0e-15_kp, &
241 skew_facw=1.2_kp, skew_fact=0.0_kp, &
242 tkhmax=300.0_kp, qcmin=1.0e-9_kp
243 real :: lsub, fac_cond, fac_fus, cpolv, fac_sub, ggri, kapa, gocp, &
244 rog, sqrtpii, epsterm, onebeps, onebrvcp
245
246! SHOC tunable parameters
247
248 real, parameter :: lambda = 0.04_kp
249! real, parameter :: min_tke = 1.0e-6_kp ! Minumum TKE value, m**2/s**2
250 real, parameter :: min_tke = 1.0e-4_kp ! Minumum TKE value, m**2/s**2
251! real, parameter :: max_tke = 100.0_kp ! Maximum TKE value, m**2/s**2
252 real, parameter :: max_tke = 40.0_kp ! Maximum TKE value, m**2/s**2
253! Maximum turbulent eddy length scale, m
254! real, parameter :: max_eddy_length_scale = 2000.0_kp
255 real, parameter :: max_eddy_length_scale = 1000.0_kp
256! Maximum "return-to-isotropy" time scale, s
257 real, parameter :: max_eddy_dissipation_time_scale = 2000.0_kp
258 real, parameter :: Pr = 1.0_kp ! Prandtl number
259
260! Constants for the TKE dissipation term based on Deardorff (1980)
261 real, parameter :: pt19=0.19_kp, pt51=0.51_kp, pt01=0.01_kp, atmin=0.01_kp, atmax=one-atmin
262 real, parameter :: Cs = 0.15_kp, epsln=1.0e-6_kp
263! real, parameter :: Ck = 0.2_kp ! Coeff in the eddy diffusivity - TKE relationship, see Eq. 7 in BK13
264 real, parameter :: Ck = 0.1_kp ! Coeff in the eddy diffusivity - TKE relationship, see Eq. 7 in BK13
265
266! real, parameter :: Ce = Ck**3/(0.7*Cs**4)
267! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 2.2
268! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 3.0 , Ces = Ce
269! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 2.5 , Ces = Ce * 3.0 / 2.5
270! real, parameter :: Ces = Ce/0.7*3.0
271
272! real, parameter :: Ce = Ck**3/(0.7*Cs**4), Ces = Ce*3.0/0.7 ! Commented Moor
273
274 real, parameter :: Ce = ck**3/cs**4, ces = ce
275! real, parameter :: Ce = Ck**3/Cs**4, Ces = Ce*3.0/0.7
276
277! real, parameter :: vonk=0.35 ! Von Karman constant
278 real, parameter :: vonk=0.4_kp ! Von Karman constant Moorthi - as in GFS
279 real, parameter :: tscale=400.0_kp ! time scale set based off of similarity results of BK13, s
280 real, parameter :: w_tol_sqd = 4.0e-04_kp ! Min vlaue of second moment of w
281! real, parameter :: w_tol_sqd = 1.0e-04_kp ! Min vlaue of second moment of w
282 real, parameter :: w_thresh = 0.0_kp, thresh = 0.0_kp
283 real, parameter :: w3_tol = 1.0e-20_kp ! Min vlaue of third moment of w
284
285
286! These parameters are a tie-in with a microphysical scheme
287! Double check their values for the Zhao-Carr scheme.
288 real, parameter :: tbgmin = 233.16_kp ! Minimum temperature for cloud water., K (ZC)
289! real, parameter :: tbgmin = 258.16_kp ! Minimum temperature for cloud water., K (ZC)
290! real, parameter :: tbgmin = 253.16_kp ! Minimum temperature for cloud water., K
291 real, parameter :: tbgmax = 273.16_kp ! Maximum temperature for cloud ice, K
292 real, parameter :: a_bg = one/(tbgmax-tbgmin)
293!
294! Parameters to tune the second order moments- No tuning is performed currently
295
296! real, parameter :: thl2tune = 2.0_kp, qw2tune = 2.0_kp, qwthl2tune = 2.0_kp, &
297 real, parameter :: thl2tune = 1.0_kp, qw2tune = 1.0_kp, qwthl2tune = 1.0_kp, &
298! thl_tol = 1.0e-4_kp, rt_tol = 1.0e-8_kp, basetemp = 300.0_kp
299 thl_tol = 1.0e-2_kp, rt_tol = 1.0e-4_kp
300
301 integer, parameter :: nitr=6
302
303! Local variables. Note that pressure is in millibars in the SHOC code.
304
305 real zl (nx,nzm) ! height of the pressure levels above surface, m
306 real zi (nx,nz) ! height of the interface levels, m
307 real adzl (nx,nzm) ! layer thickness i.e. zi(k+1)-zi(k) - defined at levels
308 real adzi (nx,nz) ! level thickness i.e. zl(k)-zl(k-1) - defined at interface
309
310 real hl (nx,nzm) ! liquid/ice water static energy , K
311 real qv (nx,nzm) ! water vapor, kg/kg
312 real qcl (nx,nzm) ! liquid water (condensate), kg/kg
313 real qci (nx,nzm) ! ice water (condensate), kg/kg
314 real w (nx,nzm) ! z-wind, m/s
315 real bet (nx,nzm) ! ggr/tv0
316 real gamaz (nx,nzm) ! ggr/cp*z
317
318! Moments of the trivariate double Gaussian PDF for the SGS total water mixing ratio
319! SGS liquid/ice static energy, and vertical velocity
320
321 real qw_sec (nx,nzm) ! Second moment total water mixing ratio, kg^2/kg^2
322 real thl_sec (nx,nzm) ! Second moment liquid/ice static energy, K^2
323 real qwthl_sec(nx,nzm) ! Covariance tot. wat. mix. ratio and static energy, K*kg/kg
324 real wqw_sec (nx,nzm) ! Turbulent flux of tot. wat. mix., kg/kg*m/s
325 real wthl_sec (nx,nzm) ! Turbulent flux of liquid/ice static energy, K*m/s
326 real w_sec (nx,nzm) ! Second moment of vertical velocity, m**2/s**2
327 real w3 (nx,nzm) ! Third moment of vertical velocity, m**3/s**3
328 real wqp_sec (nx,nzm) ! Turbulent flux of precipitation, kg/kg*m/s
329
330! Eddy length formulation
331 real smixt (nx,nzm) ! Turbulent length scale, m
332 real isotropy (nx,nzm) ! "Return-to-isotropy" eddy dissipation time scale, s
333! real isotropy_debug (nx,nzm) ! Return to isotropy scale, s without artificial limits
334 real brunt (nx,nzm) ! Moist Brunt-Vaisalla frequency, s^-1
335 real conv_vel2(nx,nzm) ! Convective velocity scale cubed, m^3/s^3
336
337 real cek(nx)
338
339! Output of SHOC
340 real diag_frac, diag_qn, diag_qi, diag_ql
341
342! real diag_frac(nx,nzm) ! SGS cloud fraction
343! real diag_qn (nx,nzm) ! SGS cloud+ice condensate, kg/kg
344! real diag_qi (nx,nzm) ! SGS ice condensate, kg/kg
345! real diag_ql (nx,nzm) ! SGS liquid condensate, kg/kg
346
347
348! Horizontally averaged variables
349! real conv_vel(nzm) ! Convective velocity scale cubed, m^3/s^3
350 real wqlsb (nzm) ! liquid water flux, kg/kg/ m/s
351 real wqisb (nzm) ! ice flux, kg/kg m/s
352! real thlv (nzm) ! Grid-scale level-average virtual potential temperature
353! (not used)
354
355
356! Local variables
357
358! real, dimension(nx,nzm) :: tkesbbuoy, tkesbshear, tkesbdiss, tkesbbuoy_debug &
359! tkebuoy_sgs, total_water, tscale1_debug, brunt2
360
361 real, dimension(nx,nzm) :: total_water, brunt2, thv, tkesbdiss
362 real, dimension(nx,nzm) :: def2
363 real, dimension(nx) :: denom, numer, l_inf, cldarr, thedz, thedz2
364
365 real lstarn, depth, omn, betdz, bbb, term, qsatt, dqsat, &
366 conv_var, tkes, skew_w, skew_qw, aterm, w1_1, w1_2, w2_1, &
367 w2_2, w3var, thl1_1, thl1_2, thl2_1, thl2_2, qw1_1, qw1_2, qw2_1, &
368 qw2_2, ql1, ql2, w_ql1, w_ql2, &
369 r_qwthl_1, r_wqw_1, r_wthl_1, testvar, s1, s2, std_s1, std_s2, C1, C2, &
370 thl_first, qw_first, w_first, Tl1_1, Tl1_2, betatest, pval, pkap, &
371 w2thl, w2qw,w2ql, w2ql_1, w2ql_2, &
372 thec, thlsec, qwsec, qwthlsec, wqwsec, wthlsec, thestd,dum, &
373 cqt1, cthl1, cqt2, cthl2, qn1, qn2, qi1, qi2, omn1, omn2, &
374 basetemp2, beta1, beta2, qs1, qs2, &
375 esval, esval2, om1, om2, epss, &
376 lstarn1, lstarn2, sqrtw2, sqrtthl, sqrtqt, &
377 sqrtstd1, sqrtstd2, tsign, tvar, sqrtw2t, wqls, wqis, &
378 sqrtqw2_1, sqrtqw2_2, sqrtthl2_1, sqrtthl2_2, sm, prespot, &
379 corrtest1, corrtest2, wrk, wrk1, wrk2, wrk3, onema, pfac
380
381
382 integer i,k,km1,ku,kd,ka,kb
383
384
385!calculate derived constants
386 lsub = lcond+lfus
387 fac_cond = lcond/cp
388 fac_fus = lfus/cp
389 cpolv = cp/lcond
390 fac_sub = lsub/cp
391 ggri = one/ggr
392 kapa = rgas/cp
393 gocp = ggr/cp
394 rog = rgas*ggri
395 sqrtpii = one/sqrt(pi+pi)
396 epsterm = rgas/rv
397 onebeps = one/epsterm
398 onebrvcp = one/(rv*cp)
399 epss = eps * supice
400
401! Map GFS variables to those of SHOC - SHOC operates on 3D fields
402! Here a Y-dimension is added to the input variables, along with some unit conversions
403
404 do k=1,nz
405 do i=1,nx
406 zi(i,k) = phii(i,k) * ggri
407 enddo
408 enddo
409
410!
411! move water from vapor to condensate if the condensate is negative
412!
413 do k=1,nzm
414 do i=1,nx
415 if (qc(i,k) < zero) then
416 qwv(i,k) = qwv(i,k) + qc(i,k)
417 tabs(i,k) = tabs(i,k) - fac_cond * qc(i,k)
418 qc(i,k) = zero
419 endif
420 if (qi(i,k) < zero) then
421 qwv(i,k) = qwv(i,k) + qi(i,k)
422 tabs(i,k) = tabs(i,k) - fac_sub * qi(i,k)
423 qi(i,k) = zero
424 endif
425!
426! testing removal of ice when too warm to sustain ice
427!
428! if (qi(i,k) > zero .and. tabs(i,k) > 273.16) then
429! wrk = (tabs(i,k) - 273.16) / fac_sub
430! if (wrk < qi(i,k)) then
431! wrk = qi(i,k) - wrk
432! qi(i,k) = wrk
433! qwv(i,k) = qwv(i,k) + wrk
434! tabs(i,k) = 273.16
435! else
436! tabs(i,k) = tabs(i,k) - qi(i,k) / fac_sub
437! qwv(i,k) = qwv(i,k) + qi(i,k)
438! qi(i,k) = 0.0
439! endif
440! endif
441
442 enddo
443 enddo
444! fill negative water vapor from below
445 do k=nzm,2,-1
446 km1 = k - 1
447 do i=1,nx
448 if (qwv(i,k) < zero) then
449 qwv(i,k) = qwv(i,km1) + qwv(i,k) * delp(i,k) / delp(i,km1)
450 endif
451 enddo
452 enddo
453
454 do k=1,nzm
455 do i=1,nx
456 zl(i,k) = phil(i,k) * ggri
457 wrk = one / prsl(i,k)
458 qv(i,k) = max(qwv(i,k), zero)
459 thv(i,k) = tabs(i,k) * (one+epsv*qv(i,k))
460 w(i,k) = - rog * omega(i,k) * thv(i,k) * wrk
461 qcl(i,k) = max(qc(i,k), zero)
462 qci(i,k) = max(qi(i,k), zero)
463!
464! qpl(i,k) = zero ! comment or remove when using with prognostic rain/snow
465! qpi(i,k) = zero ! comment or remove when using with prognostic rain/snow
466
467 wqp_sec(i,k) = zero ! Turbulent flux of precipiation
468!
469 total_water(i,k) = qcl(i,k) + qci(i,k) + qv(i,k)
470
471 prespot = (100000.0_kp*wrk) ** kapa ! Exner function
472 bet(i,k) = ggr/(tabs(i,k)*prespot) ! Moorthi
473 thv(i,k) = thv(i,k)*prespot ! Moorthi
474!
475! Lapse rate * height = reference temperature
476 gamaz(i,k) = gocp * zl(i,k)
477
478! Liquid/ice water static energy - ! Note the the units are degrees K
479 hl(i,k) = tabs(i,k) + gamaz(i,k) - fac_cond*(qcl(i,k)+qpl(i,k)) &
480 - fac_sub *(qci(i,k)+qpi(i,k))
481 w3(i,k) = zero
482 enddo
483 enddo
484
485! Define vertical grid increments for later use in the vertical differentiation
486
487 do k=2,nzm
488 km1 = k - 1
489 do i=1,nx
490 adzi(i,k) = zl(i,k) - zl(i,km1)
491 adzl(i,km1) = zi(i,k) - zi(i,km1)
492 enddo
493 enddo
494 do i=1,nx
495 adzi(i,1) = (zl(i,1)-zi(i,1)) ! unused in the code
496 adzi(i,nz) = adzi(i,nzm) ! at the top - probably unused
497 adzl(i,nzm) = zi(i,nz) - zi(i,nzm)
498!
499 wthl_sec(i,1) = hflx(i)
500 wqw_sec(i,1) = evap(i)
501 enddo
502
503
504 call tke_shoc() ! Integrate prognostic TKE equation forward in time
505
506
507! diagnose second order moments of the subgrid PDF following
508! Redelsperger J.L., and G. Sommeria, 1986, JAS, 43, 2619-2635 sans the use of stabilty
509! weighting functions - Result is in global variables w_sec, thl_sec, qw_sec, and qwthl_sec
510
511! call diag_moments(total_water,tke,tkh)
512
513! Second moment of vertical velocity.
514! Note that Eq 6 in BK13 gives a different expression that is dependent on
515! vertical gradient of grid scale vertical velocity
516
517 do k=1,nzm
518 ku = k+1
519 kd = k-1
520 ka = ku
521 kb = k
522 if (k == 1) then
523 kd = k
524 kb = ka
525 elseif (k == nzm) then
526 ku = k
527 ka = kb
528 endif
529 do i=1,nx
530 if (tke(i,k) > zero) then
531! wrk = half*(tkh(i,ka)+tkh(i,kb))*(w(i,ku) - w(i,kd)) &
532 wrk = half*(tkh(i,ka)*prnum(i,ka)+tkh(i,kb)*prnum(i,kb))*(w(i,ku) - w(i,kd)) &
533 * sqrt(tke(i,k)) / (zl(i,ku) - zl(i,kd))
534 w_sec(i,k) = max(twoby3 * tke(i,k) - twoby15 * wrk, zero)
535! w_sec(i,k) = max(twoby3 * tke(i,k), zero)
536 else
537 w_sec(i,k) = zero
538 endif
539 enddo
540 enddo
541
542 do k=2,nzm
543
544 km1 = k - 1
545 do i=1,nx
546
547! Use backward difference in the vertical, use averaged values of "return-to-isotropy"
548! time scale and diffusion coefficient
549
550 wrk1 = one / adzi(i,k) ! adzi(k) = (zl(k)-zl(km1))
551! wrk3 = max(tkh(i,k),pt01) * wrk1
552 wrk3 = max(tkh(i,k),epsln) * wrk1
553
554 sm = half*(isotropy(i,k)+isotropy(i,km1))*wrk1*wrk3 ! Tau*Kh/dz^2
555
556! SGS vertical flux liquid/ice water static energy. Eq 1 in BK13
557! No rain, snow or graupel in pdf (Annig, 08/29/2018)
558
559 wrk1 = hl(i,k) - hl(i,km1) &
560 + (qpl(i,k) - qpl(i,km1)) * fac_cond &
561 + (qpi(i,k) - qpi(i,km1)) * fac_sub
562 wthl_sec(i,k) = - wrk3 * wrk1
563
564! SGS vertical flux of total water. Eq 2 in BK13
565
566 wrk2 = total_water(i,k) - total_water(i,km1)
567 wqw_sec(i,k) = - wrk3 * wrk2
568
569! Second moment of liquid/ice water static energy. Eq 4 in BK13
570
571 thl_sec(i,k) = thl2tune * sm * wrk1 * wrk1
572
573! Second moment of total water mixing ratio. Eq 3 in BK13
574
575 qw_sec(i,k) = qw2tune * sm * wrk2 * wrk2
576
577! Covariance of total water mixing ratio and liquid/ice water static energy.
578! Eq 5 in BK13
579
580 qwthl_sec(i,k) = qwthl2tune * sm * wrk1 * wrk2
581
582 enddo ! i loop
583 enddo ! k loop
584
585! These would be at the surface - do we need them?
586 do i=1,nx
587! wthl_sec(i,1) = wthl_sec(i,2)
588! wqw_sec(i,1) = wqw_sec(i,2)
589 thl_sec(i,1) = thl_sec(i,2)
590 qw_sec(i,1) = qw_sec(i,2)
591 qwthl_sec(i,1) = qwthl_sec(i,2)
592 enddo
593
594! Diagnose the third moment of SGS vertical velocity
595
596 call canuto()
597
598! Recover parameters of the subgrid PDF using diagnosed moments
599! and calculate SGS cloudiness, condensation and it's effects on temeperature
600! and moisture variables
601
602 call assumed_pdf()
603
604contains
605
606 subroutine tke_shoc()
607
608! This subroutine solves the TKE equation,
609! Heavily based on SAM's tke_full.f90 by Marat Khairoutdinov
610
611 real grd,betdz,Cee,lstarn, lstarp, bbb, omn, omp,qsatt,dqsat, smix, &
612 buoy_sgs,ratio,a_prod_sh,a_prod_bu,a_diss,a_prod_bu_debug, buoy_sgs_debug, &
613 tscale1, wrk, wrk1, wtke, wtk2, rdtn, tkef2
614 integer i,k,ku,kd,itr,k1
615
616 rdtn = one / dtn
617
618 call tke_shear_prod(def2) ! Calculate shear production of TKE
619
620! Ensure values of TKE are reasonable
621
622 do k=1,nzm
623 do i=1,nx
624 tke(i,k) = max(min_tke,tke(i,k))
625 tkesbdiss(i,k) = zero
626! tkesbshear(i,k) = zero
627! tkesbbuoy(i,k) = zero
628 enddo
629 enddo
630
631 call eddy_length() ! Find turbulent mixing length
632 call check_eddy() ! Make sure it's reasonable
633
634 tkef2 = one - tkef1
635 do k=1,nzm
636 ku = k+1
637 kd = k
638
639! Cek = Ce * cefac
640
641 if(k == 1) then
642 ku = 2
643 kd = 2
644! Cek = Ces
645 elseif(k == nzm) then
646 ku = k
647 kd = k
648! Cek = Ces
649 endif
650
651 if (dis_opt > 0) then
652 do i=1,nx
653 wrk = (zl(i,k)-zi(i,1)) / adzl(i,1) + 1.5_kp
654 cek(i) = (one + two / max((wrk*wrk - 3.3_kp), 0.5_kp)) * cefac
655 enddo
656 else
657 if (k == 1) then
658 cek = ces * cesfac
659 else
660 cek = ce * cefac
661 endif
662 endif
663
664 do i=1,nx
665 grd = adzl(i,k) ! adzl(k) = zi(k+1)-zi(k)
666
667
668! TKE boyancy production term. wthv_sec (buoyancy flux) is calculated in
669! assumed_pdf(). The value used here is from the previous time step
670
671 a_prod_bu = ggr / thv(i,k) * wthv_sec(i,k)
672
673! If wthv_sec from subgrid PDF is not available use Brunt-Vaisalla frequency from eddy_length()
674
675!Obtain Brunt-Vaisalla frequency from diagnosed SGS buoyancy flux
676!Presumably it is more precise than BV freq. calculated in eddy_length()?
677
678 buoy_sgs = - (a_prod_bu+a_prod_bu) / (tkh(i,ku)+tkh(i,kd) + 0.0001_kp) ! tkh is eddy thermal diffussivity
679
680
681!Compute $c_k$ (variable Cee) for the TKE dissipation term following Deardorff (1980)
682
683 if (buoy_sgs <= zero) then
684 smix = grd
685 else
686 smix = min(grd,max(0.1_kp*grd, 0.76_kp*sqrt(tke(i,k)/(buoy_sgs+1.0e-10_kp))))
687 endif
688
689 ratio = smix/grd
690 cee = cek(i) * (pt19 + pt51*ratio) * max(one, sqrt(pcrit/prsl(i,k)))
691
692! TKE shear production term
693 a_prod_sh = half*(def2(i,ku)*tkh(i,ku)*prnum(i,ku) &
694 + def2(i,kd)*tkh(i,kd)*prnum(i,kd))
695
696
697! smixt (turb. mixing lenght) is calculated in eddy_length()
698! Explicitly integrate TKE equation forward in time
699! a_diss = Cee/smixt(i,k)*tke(i,k)**1.5 ! TKE dissipation term
700! tke(i,k) = max(zero,tke(i,k)+dtn*(max(zero,a_prod_sh+a_prod_bu)-a_diss))
701
702! Semi-implicitly integrate TKE equation forward in time
703
704 wtke = tke(i,k)
705 wtk2 = wtke
706! wrk = (dtn*Cee)/smixt(i,k)
707 wrk = (dtn*cee) / smixt(i,k)
708 wrk1 = wtke + dtn*(a_prod_sh+a_prod_bu)
709
710 do itr=1,nitr ! iterate for implicit solution
711 wtke = min(max(min_tke, wtke), max_tke)
712 a_diss = wrk*sqrt(wtke) ! Coefficient in the TKE dissipation term
713 wtke = wrk1 / (one+a_diss)
714 wtke = tkef1*wtke + tkef2*wtk2 ! tkef1+tkef2 = 1.0
715 wtk2 = wtke
716 enddo
717
718 tke(i,k) = min(max(min_tke, wtke), max_tke)
719 a_diss = wrk*sqrt(tke(i,k))
720
721 tscale1 = (dtn+dtn) / a_diss ! corrected Eq 8 in BK13 -- tau = 2*tke/eps
722
723 tkesbdiss(i,k) = rdtn*a_diss*tke(i,k) ! TKE dissipation term, epsilon
724
725
726! Calculate "return-to-isotropy" eddy dissipation time scale, see Eq. 8 in BK13
727
728 if (buoy_sgs <= zero) then
729 isotropy(i,k) = min(max_eddy_dissipation_time_scale, tscale1)
730 else
731 isotropy(i,k) = min(max_eddy_dissipation_time_scale, &
732 tscale1/(one+lambda*buoy_sgs*tscale1*tscale1))
733 endif
734
735! TKE budget terms
736
737! tkesbdiss(i,k) = a_diss
738! tkesbshear(i,k) = a_prod_sh
739! tkesbbuoy(i,k) = a_prod_bu
740! tkesbbuoy_debug(i,k) = a_prod_bu_debug
741! tkebuoy_sgs(i,k) = buoy_sgs
742
743 enddo ! i loop
744 enddo ! k loop
745 wrk = half * ck
746 do k=2,nzm
747 k1 = k - 1
748 do i=1,nx
749 tkh(i,k) = min(tkhmax, wrk * (isotropy(i,k) * tke(i,k) &
750 + isotropy(i,k1) * tke(i,k1))) ! Eddy thermal diffusivity
751 enddo ! i
752 enddo ! k
753
754
755 end subroutine tke_shoc
756
757
758 subroutine tke_shear_prod(def2)
759
760! Calculate TKE shear production term
761
762 real, intent(out) :: def2(nx,nzm)
763
764 real rdzw, wrku, wrkv, wrkw
765 integer i,k,k1
766
767! Calculate TKE shear production term at layer interface
768
769 do k=2,nzm
770 k1 = k - 1
771 do i=1,nx
772 rdzw = one / adzi(i,k)
773 wrku = (u(i,k)-u(i,k1)) * rdzw
774 wrkv = (v(i,k)-v(i,k1)) * rdzw
775! wrkw = (w(i,k)-w(i,k1)) * rdzw
776 def2(i,k) = wrku*wrku + wrkv*wrkv !+ 2*wrkw(1) * wrkw(1)
777 enddo
778 enddo ! k loop
779 do i=1,nx
780! def2(i,1) = def2(i,2)
781 def2(i,1) = (u(i,1)*u(i,1) + v(i,1)*v(i,1)) / (zl(i,1)*zl(i,1))
782 enddo
783
784 end subroutine tke_shear_prod
785
786 subroutine eddy_length()
787
788! This subroutine computes the turbulent length scale based on a new
789! formulation described in BK13
790
791! Local variables
792 real wrk, wrk1, wrk2, wrk3
793 integer i, k, kk, kl, ku, kb, kc, kli, kui
794
795 do i=1,nx
796 cldarr(i) = zero
797 numer(i) = zero
798 denom(i) = zero
799 enddo
800
801! Find the length scale outside of clouds, that includes boundary layers.
802
803 do k=1,nzm
804 do i=1,nx
805
806! Reinitialize the mixing length related arrays to zero
807! smixt(i,k) = one ! shoc_mod module variable smixt
808 smixt(i,k) = epsln ! shoc_mod module variable smixt
809 brunt(i,k) = zero
810
811!Eq. 11 in BK13 (Eq. 4.13 in Pete's dissertation)
812!Outside of cloud, integrate from the surface to the cloud base
813!Should the 'if' below check if the cloud liquid < a small constant instead?
814
815 if (qcl(i,k)+qci(i,k) <= qcmin) then
816 tkes = sqrt(tke(i,k)) * adzl(i,k)
817 numer(i) = numer(i) + tkes*zl(i,k) ! Numerator in Eq. 11 in BK13
818 denom(i) = denom(i) + tkes ! Denominator in Eq. 11 in BK13
819 else
820 cldarr(i) = one ! Take note of columns containing cloud.
821 endif
822 enddo
823 enddo
824
825! Calculate the measure of PBL depth, Eq. 11 in BK13 (Is this really PBL depth?)
826 do i=1,nx
827 if (denom(i) > zero .and. numer(i) > zero) then
828 l_inf(i) = min(0.1_kp * (numer(i)/denom(i)), 100.0_kp)
829 else
830 l_inf(i) = 100.0_kp
831 endif
832 enddo
833
834!Calculate length scale outside of cloud, Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation)
835 do k=1,nzm
836
837 kb = k-1
838 kc = k+1
839 if (k == 1) then
840 kb = 1
841 kc = 2
842 thedz(:) = adzi(:,kc)
843 elseif (k == nzm) then
844 kb = nzm-1
845 kc = nzm
846 thedz(:) = adzi(:,k)
847 else
848 thedz(:) = adzi(:,kc) + adzi(:,k) ! = (z(k+1)-z(k-1))
849 endif
850
851 do i=1,nx
852
853! vars module variable bet (=ggr/tv0) ; grid module variable adzi
854
855 betdz = bet(i,k) / thedz(i)
856
857 tkes = sqrt(tke(i,k))
858
859! Compute local Brunt-Vaisalla frequency
860
861 wrk = qcl(i,k) + qci(i,k)
862 if (wrk > zero) then ! If in the cloud
863
864! Find the in-cloud Brunt-Vaisalla frequency
865
866 omn = qcl(i,k) / (wrk+1.0e-20_kp) ! Ratio of liquid water to total water
867
868! Latent heat of phase transformation based on relative water phase content
869! fac_cond = lcond/cp, fac_fus = lfus/cp
870
871 lstarn = fac_cond + (one-omn)*fac_fus
872
873! Derivative of saturation mixing ratio over water/ice wrt temp. based on relative water phase content
874 dqsat = omn * dtqsatw(tabs(i,k),prsl(i,k)) &
875 + (one-omn) * dtqsati(tabs(i,k),prsl(i,k))
876
877! Saturation mixing ratio over water/ice wrt temp based on relative water phase content
878
879 qsatt = omn * qsatw(tabs(i,k),prsl(i,k)) &
880 + (one-omn) * qsati(tabs(i,k),prsl(i,k))
881
882! liquid/ice moist static energy static energy divided by cp?
883
884 bbb = (one + epsv*qsatt-wrk-qpl(i,k)-qpi(i,k) &
885 + 1.61_kp*tabs(i,k)*dqsat) / (one+lstarn*dqsat)
886
887! Calculate Brunt-Vaisalla frequency using centered differences in the vertical
888
889 brunt(i,k) = betdz*(bbb*(hl(i,kc)-hl(i,kb)) &
890 + (bbb*lstarn - (one+lstarn*dqsat)*tabs(i,k)) &
891 * (total_water(i,kc)-total_water(i,kb)) &
892 + (bbb*fac_cond - (one+fac_cond*dqsat)*tabs(i,k))*(qpl(i,kc)-qpl(i,kb)) &
893 + (bbb*fac_sub - (one+fac_sub*dqsat)*tabs(i,k))*(qpi(i,kc)-qpi(i,kb)) )
894
895 else ! outside of cloud
896
897! Find outside-of-cloud Brunt-Vaisalla frequency
898! Only unsaturated air, rain and snow contribute to virt. pot. temp.
899! liquid/ice moist static energy divided by cp?
900
901 bbb = one + epsv*qv(i,k) - qpl(i,k) - qpi(i,k)
902 brunt(i,k) = betdz*( bbb*(hl(i,kc)-hl(i,kb)) &
903 + epsv*tabs(i,k)*(total_water(i,kc)-total_water(i,kb)) &
904 + (bbb*fac_cond-tabs(i,k))*(qpl(i,kc)-qpl(i,kb)) &
905 + (bbb*fac_sub -tabs(i,k))*(qpi(i,kc)-qpi(i,kb)) )
906 endif
907
908! Reduction of mixing length in the stable regions (where B.-V. freq. > 0) is required.
909! Here we find regions of Brunt-Vaisalla freq. > 0 for later use.
910
911 if (brunt(i,k) >= zero) then
912 brunt2(i,k) = brunt(i,k)
913 else
914 brunt2(i,k) = zero
915 endif
916
917! Calculate turbulent length scale in the boundary layer.
918! See Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation)
919
920! Keep the length scale adequately small near the surface following Blackadar (1984)
921! Note that this is not documented in BK13 and was added later for SP-CAM runs
922
923! if (k == 1) then
924! term = 600.*tkes
925! smixt(i,k) = term + (0.4*zl(i,k)-term)*exp(-zl(i,k)*0.01)
926! else
927
928! tscale is the eddy turnover time scale in the boundary layer and is
929! an empirically derived constant
930
931 if (tkes > zero .and. l_inf(i) > zero) then
932 wrk1 = one / (tscale*tkes*vonk*zl(i,k))
933 wrk2 = one / (tscale*tkes*l_inf(i))
934 wrk1 = wrk1 + wrk2 + pt01 * brunt2(i,k) / tke(i,k)
935 wrk1 = sqrt(one / max(wrk1,1.0e-8_kp)) * (one/0.3_kp)
936! smixt(i,k) = min(max_eddy_length_scale, 2.8284*sqrt(wrk1)/0.3)
937 smixt(i,k) = min(max_eddy_length_scale, wrk1)
938
939! smixt(i,k) = min(max_eddy_length_scale,(2.8284*sqrt(1./((1./(tscale*tkes*vonk*zl(i,k))) &
940! + (1./(tscale*tkes*l_inf(i)))+0.01*(brunt2(i,k)/tke(i,k)))))/0.3)
941! else
942! smixt(i,k) = zero
943 endif
944
945! endif
946
947
948 enddo
949 enddo
950
951! Now find the in-cloud turbulence length scale
952! See Eq. 13 in BK13 (Eq. 4.18 in Pete's disseration)
953
954!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
955! Remove after coupling to subgrid PDF.
956!wthv_sec = -300/ggr*brunt*tk
957!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
958
959! determine cubed convective velocity scale (conv_vel2) inside the cloud
960
961! call conv_scale() ! inlining the relevant code
962
963! do i=1,nx
964! conv_vel2(i,1) = zero ! Convective velocity scale cubed
965! enddo
966 ! Integrate velocity scale in the vertical
967! do k=2,nzm
968! do i=1,nx
969! conv_vel2(i,k) = conv_vel2(i,k-1) &
970! + 2.5*adzi(i,k)*bet(i,k)*wthv_sec(i,k)
971! enddo
972! enddo
973
974 do i=1,nx
975
976 if (cldarr(i) == 1) then ! If there's a cloud in this column
977
978 kl = 0
979 ku = 0
980 do k=2,nzm-3
981
982! Look for the cloud base in this column
983! thresh (=0) is a variable local to eddy_length(). Should be a module constant.
984 wrk = qcl(i,k) + qci(i,k)
985 if (wrk > qcmin) then
986 if (kl == 0) then
987 kl = k
988 endif
989
990! Look for the cloud top in this column
991 if (qcl(i,k+1)+qci(i,k+1) <= qcmin) then
992 ku = k
993! conv_vel2 (Cubed convective velocity scale) is calculated in conv_scale()
994! Use the value of conv_vel2 at the top of the cloud.
995! conv_var = conv_vel2(i,k)** oneb3
996 endif
997 endif
998
999! Compute the mixing length scale for the cloud layer that we just found
1000! if (kl > 0 .and. ku > 0 .and. ku-kl > 1) then
1001! if (kl > 0 .and. ku > 0 .and. ku-kl > 0) then
1002 if (kl > 0 .and. ku >= kl) then
1003! The calculation below finds the integral in the Eq. 10 in BK13 for the current cloud
1004 conv_var = zero
1005 do kk=kl,ku
1006 conv_var = conv_var+ 2.5_kp*adzi(i,kk)*bet(i,kk)*wthv_sec(i,kk)
1007 enddo
1008 conv_var = conv_var ** oneb3
1009
1010 if (conv_var > zero) then ! If convective vertical velocity scale > 0
1011
1012 depth = (zl(i,ku)-zl(i,kl)) + adzl(i,kl)
1013
1014 do kk=kl,ku
1015! in-cloud turbulence length scale, Eq. 13 in BK13 (Eq. 4.18)
1016
1017! wrk = conv_var/(depth*sqrt(tke(i,kk)))
1018! wrk = wrk * wrk + pt01*brunt2(i,kk)/tke(i,kk)
1019
1020 wrk = conv_var/(depth*depth*sqrt(tke(i,kk))) &
1021 + pt01*brunt2(i,kk)/tke(i,kk)
1022
1023 smixt(i,kk) = min(max_eddy_length_scale, (one/0.3_kp)*sqrt(one/wrk))
1024
1025 enddo
1026
1027 endif ! If convective vertical velocity scale > 0
1028 kl = zero
1029 ku = zero
1030 endif ! if inside the cloud layer
1031
1032 enddo ! k=2,nzm-3
1033 endif ! if in the cloudy column
1034 enddo ! i=1,nx
1035
1036
1037 end subroutine eddy_length
1038
1039
1040 subroutine conv_scale()
1041
1042! This subroutine calculates the cubed convective velocity scale needed
1043! for the definition of the length scale in clouds
1044! See Eq. 16 in BK13 (Eq. 4.21 in Pete's dissertation)
1045
1046 integer i, k
1047
1048!!!!!!!!!
1049!! A bug in formulation of conv_vel
1050! Obtain it by averaging conv_vel2 in the horizontal
1051!!!!!!!!!!
1052
1053! conv_vel(1)=zero ! Horizontally averaged convective velocity scale cubed
1054 do i=1,nx
1055 conv_vel2(i,1) = zero ! Convective velocity scale cubed
1056 enddo
1057! Integrate velocity scale in the vertical
1058 do k=2,nzm
1059! conv_vel(k)=conv_vel(k-1)
1060 do i=1,nx
1061!**********************************************************************
1062!Do not include grid-scale contribution to convective velocity scale in GCM applications
1063! conv_vel(k)=conv_vel(k-1)+2.5*adzi(k)*bet(k)*(tvwle(k)+tvws(k))
1064! conv_vel(k)=conv_vel(k)+2.5*adzi(i,k)*bet(i,k)*(tvws(k))
1065!Do not include grid-scale contribution to convective velocity scale in GCM applications
1066! conv_vel2(i,k)=conv_vel2(i,k-1)+2.5*adzi(k)*bet(k)*(tvwle(k)+wthv_sec(i,k))
1067!**********************************************************************
1068
1069 conv_vel2(i,k) = conv_vel2(i,k-1) &
1070 + 2.5_kp*adzi(i,k)*bet(i,k)*wthv_sec(i,k)
1071 enddo
1072 enddo
1073
1074 end subroutine conv_scale
1075
1076
1077 subroutine check_eddy()
1078
1079! This subroutine checks eddy length values
1080
1081 integer i, k, kb, ks, zend
1082 real wrk
1083! real zstart, zthresh, qthresh
1084
1085! Temporary kludge for marine stratocumulus under very strong inversions at coarse resolution
1086! Placement until some explicity PBL top is put in
1087! Not used.
1088! zthresh = 100.
1089! qthresh = -6.0
1090
1091 do k=1,nzm
1092
1093 if (k == nzm) then
1094 kb = k
1095 else
1096 kb = k+1
1097 endif
1098
1099 do i=1,nx
1100
1101 wrk = 0.1_kp*adzl(i,k)
1102 ! Minimum 0.1 of local dz
1103 smixt(i,k) = max(wrk, min(max_eddy_length_scale,smixt(i,k)))
1104
1105! If chracteristic grid dimension in the horizontal< 1000m, set lengthscale to
1106! be not larger that that.
1107! if (sqrt(dx*dy) .le. 1000.) smixt(i,k)=min(sqrt(dx*dy),smixt(i,k))
1108
1109 if (qcl(i,kb) == zero .and. qcl(i,k) > zero .and. brunt(i,k) > 1.0e-4_kp) then
1110!If just above the cloud top and atmosphere is stable, set to 0.1 of local dz
1111 smixt(i,k) = wrk
1112 endif
1113
1114 enddo ! i
1115 enddo ! k
1116
1117 end subroutine check_eddy
1118
1119 subroutine canuto()
1120
1121! Subroutine impements an analytic expression for the third moment of SGS vertical velocity
1122! based on Canuto et at, 2001, JAS, 58, 1169-1172 (further referred to as C01)
1123! This allows to avoid having a prognostic equation for the third moment.
1124! Result is returned in a global variable w3 defined at the interface levels.
1125
1126! Local variables
1127 integer i, k, kb, kc
1128
1129 real bet2, f0, f1, f2, f3, f4, f5, iso, isosqr, &
1130 omega0, omega1, omega2, X0, Y0, X1, Y1, AA0, AA1, buoy_sgs2, &
1131 wrk, wrk1, wrk2, wrk3, avew
1132! cond, wrk, wrk1, wrk2, wrk3, avew
1133!
1134! See Eq. 7 in C01 (B.7 in Pete's dissertation)
1135 real, parameter :: c=7.0_kp, a0=0.52_kp/(c*c*(c-2.0_kp)), a1=0.87_kp/(c*c), &
1136 a2=0.5_kp/c, a3=0.6_kp/(c*(c-2.0_kp)), a4=2.4_kp/(3.0_kp*c+5.0_kp), &
1137 a5=0.6_kp/(c*(3.0_kp*c+5.0_kp))
1138!Moorthi a5=0.6_kp/(c*(3.0_kp+5.0_kp*c))
1139
1140! do k=1,nzm
1141 do k=2,nzm
1142
1143 kb = k-1
1144 kc = k+1
1145
1146! if(k == 1) then
1147! kb = 1
1148! kc = 2
1149! do i=1,nx
1150! thedz(i) = one / adzl(i,kc)
1151! thedz2(i) = thedz(i)
1152! enddo
1153! elseif(k == nzm) then
1154 if(k == nzm) then
1155 kb = nzm-1
1156 kc = nzm
1157 do i=1,nx
1158 thedz(i) = one / adzi(i,k)
1159 thedz2(i) = one / adzl(i,kb)
1160 enddo
1161 else
1162 do i=1,nx
1163 thedz(i) = one / adzi(i,k)
1164 thedz2(i) = one / (adzl(i,k)+adzl(i,kb))
1165 enddo
1166 endif
1167
1168 do i=1,nx
1169
1170 iso = half*(isotropy(i,k)+isotropy(i,kb))
1171 isosqr = iso*iso ! Two-level average of "return-to-isotropy" time scale squared
1172 buoy_sgs2 = isosqr*half*(brunt(i,k)+brunt(i,kb))
1173 bet2 = half*(bet(i,k)+bet(i,kb)) !Two-level average of BV frequency squared
1174
1175
1176! Compute functions f0-f5, see Eq, 8 in C01 (B.8 in Pete's dissertation)
1177
1178
1179 avew = half*(w_sec(i,k)+w_sec(i,kb))
1180
1181
1182! This is not a bug, but an algorithmical change.
1183! The line below calculates cond_w ,an estimate of the maximum allowed value of the third moment.
1184! It is used at the end of this subroutine to limit the value of w3.
1185! Here the second moment is interpolated from the layer centers to the interface, where w3 is defined.
1186! In the presence of strong vertical gradients of w2, the value interpolated to the interface can
1187! be as much as twice as as large (or as small) as the value on in layer center. When the skewness
1188! of W PDF is calculated in assumed_pdf(), the code there uses w2 on the layer center, and the value
1189! of w3 interpolated from the interfaces to the layer center. The errors introduced due to dual
1190! interpolation are amplified by exponentiation during the calculation of skewness
1191! and result in (ususally negative) values
1192! of skewness of W PDF that are too large ( < -10). The resulting PDF consists of two delta
1193! functions set far apart from each other on the W axis, and produces unrealistically
1194! high values of diagnosed SGS fluxes.
1195
1196! The solution is to move clipping of the third moment from canuto() to right before the
1197! calculation of skewness in assumed_pdf()
1198! cond_w = 1.2*sqrt(max(1.0e-20, 2.*avew*avew*avew))
1200!
1201
1202 wrk1 = bet2*iso
1203 wrk2 = thedz2(i)*wrk1*wrk1*iso
1204 wrk3 = thl_sec(i,kc) - thl_sec(i,kb)
1205
1206 f0 = wrk2 * wrk1 * wthl_sec(i,k) * wrk3
1207
1208 wrk = wthl_sec(i,kc) - wthl_sec(i,kb)
1209
1210 f1 = wrk2 * (wrk*wthl_sec(i,k) + half*avew*wrk3)
1211
1212 wrk1 = bet2*isosqr
1213 f2 = thedz(i)*wrk1*wthl_sec(i,k)*(w_sec(i,k)-w_sec(i,kb)) &
1214 + (thedz2(i)+thedz2(i))*bet(i,k)*isosqr*wrk
1215
1216 f3 = thedz2(i)*wrk1*wrk + thedz(i)*bet2*isosqr*(wthl_sec(i,k)*(tke(i,k)-tke(i,kb)))
1217
1218 wrk1 = thedz(i)*iso*avew
1219 f4 = wrk1*(w_sec(i,k)-w_sec(i,kb) + tke(i,k)-tke(i,kb))
1220
1221 f5 = wrk1*(w_sec(i,k)-w_sec(i,kb))
1222
1223
1224! Compute the "omega" terms, see Eq. 6 in C01 (B.6 in Pete's dissertation)
1225
1226 omega0 = a4 / (one-a5*buoy_sgs2)
1227 omega1 = omega0 / (c+c)
1228 omega2 = omega1*f3+(5.0_kp/4.0_kp)*omega0*f4
1229
1230! Compute the X0, Y0, X1, Y1 terms, see Eq. 5 a-b in C01 (B.5 in Pete's dissertation)
1231
1232 wrk1 = one / (one-(a1+a3)*buoy_sgs2)
1233 wrk2 = one / (one-a3*buoy_sgs2)
1234 x0 = wrk1 * (a2*buoy_sgs2*(one-a3*buoy_sgs2))
1235 y0 = wrk2 * (two*a2*buoy_sgs2*x0)
1236 x1 = wrk1 * (a0*f0+a1*f1+a2*(one-a3*buoy_sgs2)*f2)
1237 y1 = wrk2 * (two*a2*(buoy_sgs2*x1+(a0/a1)*f0+f1))
1238
1239! Compute the A0, A1 terms, see Eq. 5d in C01 (B.5 in Pete's dissertation)
1240
1241 aa0 = omega0*x0 + omega1*y0
1242 aa1 = omega0*x1 + omega1*y1 + omega2
1243
1244! Finally, we have the third moment of w, see Eq. 4c in C01 (B.4 in Pete's dissertation)
1245! cond_w is an estimate of third moment from second oment - If the third moment is larger
1246! than the estimate - limit w3.
1247
1248
1249! Move clipping of w3 to assumed_pdf()
1250! w3(i,k) = max(-cond_w, min(cond_w, (AA1-1.2*X1-1.5*f5)/(c-1.2*X0+AA0)))
1251 w3(i,k) = (aa1-1.2_kp*x1-1.5_kp*f5)/(c-1.2_kp*x0+aa0)
1253
1254! Implemetation of the C01 approach in this subroutine is nearly complete
1255! (the missing part are Eqs. 5c and 5e which are very simple)
1256! therefore it's easy to diagnose other third order moments obtained in C01 using this code.
1257
1258 enddo
1259 enddo
1260 do i=1,nx
1261 w3(i,1) = w3(i,2)
1262 enddo
1263
1264 end subroutine canuto
1265
1266 subroutine assumed_pdf()
1267
1268! Compute SGS buoyancy flux, SGS cloud fraction, and SGS condensation
1269! using assumed analytic double-gaussian PDF for SGS vertical velocity,
1270! moisture, and liquid/ice water static energy, based on the
1271! general approach of Larson et al 2002, JAS, 59, 3519-3539,
1272! and Golaz et al 2002, JAS, 59, 3540-3551
1273! References in the comments in this code are given to
1274! the Appendix A of Pete Bogenschutz's dissertation.
1275
1276! Local variables
1277
1278 integer i,k,ku,kd
1279 real wrk, wrk1, wrk2, wrk3, wrk4, bastoeps, eps_ss1, eps_ss2, cond_w
1280
1281! bastoeps = basetemp / epsterm
1282
1283
1284! Initialize for statistics
1285 do k=1,nzm
1286 wqlsb(k) = zero
1287 wqisb(k) = zero
1288 enddo
1289
1290 DO k=1,nzm
1291
1292 kd = k
1293 ku = k + 1
1294! if (k == nzm) ku = k
1295
1296 DO i=1,nx
1297
1298! Initialize cloud variables to zero
1299 diag_qn = zero
1300 diag_frac = zero
1301 diag_ql = zero
1302 diag_qi = zero
1303
1304 pval = prsl(i,k)
1305 pfac = pval * 1.0e-5_kp
1306 pkap = pfac ** kapa
1307
1308! Read in liquid/ice static energy, total water mixing ratio,
1309! and vertical velocity to variables PDF needs
1310 thl_first = hl(i,k) + fac_cond*qpl(i,k) + fac_sub*qpi(i,k)
1311 qw_first = total_water(i,k)
1312! w_first = half*(w(i,kd)+w(i,ku))
1313 w_first = w(i,k)
1314
1315
1316! GET ALL INPUT VARIABLES ON THE SAME GRID
1317! Points to be computed with relation to thermo point
1318! Read in points that need to be averaged
1319
1320 if (k < nzm) then
1321 w3var = half*(w3(i,kd)+w3(i,ku))
1322 thlsec = max(zero, half*(thl_sec(i,kd)+thl_sec(i,ku)) )
1323 qwsec = max(zero, half*(qw_sec(i,kd)+qw_sec(i,ku)) )
1324 qwthlsec = half * (qwthl_sec(i,kd) + qwthl_sec(i,ku))
1325 wqwsec = half * (wqw_sec(i,kd) + wqw_sec(i,ku))
1326 wthlsec = half * (wthl_sec(i,kd) + wthl_sec(i,ku))
1327 else ! at the model top assuming zeros
1328 w3var = half*w3(i,k)
1329 thlsec = max(zero, half*thl_sec(i,k))
1330 qwsec = max(zero, half*qw_sec(i,k))
1331 qwthlsec = half * qwthl_sec(i,k)
1332 wqwsec = half * wqw_sec(i,k)
1333 wthlsec = half * wthl_sec(i,k)
1334 endif
1335
1336! w3var = w3(i,k)
1337! thlsec = max(zero,thl_sec(i,k))
1338! qwsec = max(zero,qw_sec(i,k))
1339! qwthlsec = qwthl_sec(i,k)
1340! wqwsec = wqw_sec(i,k)
1341! wthlsec = wthl_sec(i,k)
1342
1343! Compute square roots of some variables so we don't have to do it again
1344 if (w_sec(i,k) > zero) then
1345 sqrtw2 = sqrt(w_sec(i,k))
1346 else
1347 sqrtw2 = zero
1348 endif
1349 if (thlsec > zero) then
1350 sqrtthl = sqrt(thlsec)
1351 else
1352 sqrtthl = zero
1353 endif
1354 if (qwsec > zero) then
1355 sqrtqt = sqrt(qwsec)
1356 else
1357 sqrtqt = zero
1358 endif
1359
1360
1361! Find parameters of the double Gaussian PDF of vertical velocity
1362
1363! Skewness of vertical velocity
1364! Skew_w = w3var / w_sec(i,k)**(3./2.)
1365! Skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2) ! Moorthi
1366
1367 IF (w_sec(i,k) <= w_tol_sqd) THEN ! If variance of w is too small then
1368 ! PDF is a sum of two delta functions
1369 skew_w = zero
1370 w1_1 = w_first
1371 w1_2 = w_first
1372 w2_1 = zero
1373 w2_2 = zero
1374 aterm = half
1375 onema = half
1376 ELSE
1377
1378! Clip w3
1379 cond_w = 1.2_kp*sqrt2*max(w3_tol, sqrtw2*sqrtw2*sqrtw2)
1380 w3var = max(-cond_w, min(cond_w, w3var))
1382
1383 skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2) ! Moorthi
1384! Proportionality coefficients between widths of each vertical velocity
1385! gaussian and the sqrt of the second moment of w
1386 w2_1 = 0.4_kp
1387 w2_2 = 0.4_kp
1388
1389! Compute realtive weight of the first PDF "plume"
1390! See Eq A4 in Pete's dissertaion - Ensure 0.01 < a < 0.99
1391
1392 wrk = one - w2_1
1393 aterm = max(atmin,min(half*(one-skew_w*sqrt(one/(4.0_kp*wrk*wrk*wrk+skew_w*skew_w))),atmax))
1394 onema = one - aterm
1395
1396 sqrtw2t = sqrt(wrk)
1397
1398! Eq. A.5-A.6
1399 wrk = sqrt(onema/aterm)
1400 w1_1 = sqrtw2t * wrk
1401 w1_2 = - sqrtw2t / wrk
1402
1403 w2_1 = w2_1 * w_sec(i,k)
1404 w2_2 = w2_2 * w_sec(i,k)
1405
1406 ENDIF
1407
1408! Find parameters of the PDF of liquid/ice static energy
1409
1410 IF (thlsec <= thl_tol*thl_tol .or. abs(w1_2-w1_1) <= w_thresh) THEN
1411 thl1_1 = thl_first
1412 thl1_2 = thl_first
1413 thl2_1 = zero
1414 thl2_2 = zero
1415 sqrtthl2_1 = zero
1416 sqrtthl2_2 = zero
1417 ELSE
1418
1419 corrtest1 = max(-one,min(one,wthlsec/(sqrtw2*sqrtthl)))
1420
1421 thl1_1 = -corrtest1 / w1_2 ! A.7
1422 thl1_2 = -corrtest1 / w1_1 ! A.8
1423
1424 wrk1 = thl1_1 * thl1_1
1425 wrk2 = thl1_2 * thl1_2
1426 wrk3 = three * (one - aterm*wrk1 - onema*wrk2)
1427 wrk4 = -skew_facw*skew_w - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2 ! testing - Moorthi
1428! wrk4 = -skew_fact*Skew_w - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2 ! testing - Moorthi
1429! wrk4 = - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2
1430 wrk = three * (thl1_2-thl1_1)
1431 if (wrk /= zero) then
1432 thl2_1 = thlsec * min(100.0_kp,max(zero,(thl1_2*wrk3-wrk4)/(aterm*wrk))) ! A.10
1433 thl2_2 = thlsec * min(100.0_kp,max(zero,(-thl1_1*wrk3+wrk4)/(onema*wrk))) ! A.11
1434 else
1435 thl2_1 = zero
1436 thl2_2 = zero
1437 endif
1438!
1439 thl1_1 = thl1_1*sqrtthl + thl_first
1440 thl1_2 = thl1_2*sqrtthl + thl_first
1441
1442 sqrtthl2_1 = sqrt(thl2_1)
1443 sqrtthl2_2 = sqrt(thl2_2)
1444
1445 ENDIF
1446
1447! FIND PARAMETERS FOR TOTAL WATER MIXING RATIO
1448
1449 IF (qwsec <= rt_tol*rt_tol .or. abs(w1_2-w1_1) <= w_thresh) THEN
1450 qw1_1 = qw_first
1451 qw1_2 = qw_first
1452 qw2_1 = zero
1453 qw2_2 = zero
1454 sqrtqw2_1 = zero
1455 sqrtqw2_2 = zero
1456 ELSE
1457
1458 corrtest2 = max(-one,min(one,wqwsec/(sqrtw2*sqrtqt)))
1459
1460 qw1_1 = - corrtest2 / w1_2 ! A.7
1461 qw1_2 = - corrtest2 / w1_1 ! A.8
1462
1463 tsign = abs(qw1_2-qw1_1)
1464
1465! Skew_qw = skew_facw*Skew_w
1466
1467 IF (tsign > 0.4_kp) THEN
1468 skew_qw = skew_facw*skew_w
1469 ELSEIF (tsign <= 0.2_kp) THEN
1470 skew_qw = zero
1471 ELSE
1472 skew_qw = (skew_facw/0.2_kp) * skew_w * (tsign-0.2_kp)
1473 ENDIF
1474
1475 wrk1 = qw1_1 * qw1_1
1476 wrk2 = qw1_2 * qw1_2
1477 wrk3 = three * (one - aterm*wrk1 - onema*wrk2)
1478 wrk4 = skew_qw - aterm*wrk1*qw1_1 - onema*wrk2*qw1_2
1479 wrk = three * (qw1_2-qw1_1)
1480
1481 if (wrk /= zero) then
1482 qw2_1 = qwsec * min(100.0_kp,max(zero,( qw1_2*wrk3-wrk4)/(aterm*wrk))) ! A.10
1483 qw2_2 = qwsec * min(100.0_kp,max(zero,(-qw1_1*wrk3+wrk4)/(onema*wrk))) ! A.11
1484 else
1485 qw2_1 = zero
1486 qw2_2 = zero
1487 endif
1488!
1489 qw1_1 = qw1_1*sqrtqt + qw_first
1490 qw1_2 = qw1_2*sqrtqt + qw_first
1491
1492 sqrtqw2_1 = sqrt(qw2_1)
1493 sqrtqw2_2 = sqrt(qw2_2)
1494
1495 ENDIF
1496
1497! CONVERT FROM TILDA VARIABLES TO "REAL" VARIABLES
1498
1499 w1_1 = w1_1*sqrtw2 + w_first
1500 w1_2 = w1_2*sqrtw2 + w_first
1501
1502! FIND WITHIN-PLUME CORRELATIONS
1503
1504 testvar = aterm*sqrtqw2_1*sqrtthl2_1 + onema*sqrtqw2_2*sqrtthl2_2
1505
1506 IF (testvar == zero) THEN
1507 r_qwthl_1 = zero
1508 ELSE
1509 r_qwthl_1 = max(-one,min(one,(qwthlsec-aterm*(qw1_1-qw_first)*(thl1_1-thl_first) &
1510 -onema*(qw1_2-qw_first)*(thl1_2-thl_first))/testvar)) ! A.12
1511 ENDIF
1512
1513! BEGIN TO COMPUTE CLOUD PROPERTY STATISTICS
1514
1515! wrk1 = gamaz(i,k) - fac_cond*qpl(i,k) - fac_sub*qpi(i,k)
1516! Tl1_1 = thl1_1 - wrk1
1517! Tl1_2 = thl1_2 - wrk1
1518
1519 tl1_1 = thl1_1 - gamaz(i,k)
1520 tl1_2 = thl1_2 - gamaz(i,k)
1521
1522! Now compute qs
1523
1524! Partition based on temperature for the first plume
1525
1526 IF (tl1_1 >= tbgmax) THEN
1527 lstarn1 = lcond
1528 esval = min(fpvsl(tl1_1), pval)
1529 qs1 = eps * esval / (pval-0.378_kp*esval)
1530 ELSE IF (tl1_1 <= tbgmin) THEN
1531 lstarn1 = lsub
1532 esval = min(fpvsi(tl1_1), pval)
1533 qs1 = epss * esval / (pval-0.378_kp*esval)
1534 ELSE
1535 om1 = max(zero, min(one, a_bg*(tl1_1-tbgmin)))
1536 lstarn1 = lcond + (one-om1)*lfus
1537 esval = min(fpvsl(tl1_1), pval)
1538 esval2 = min(fpvsi(tl1_1), pval)
1539 qs1 = om1 * eps * esval / (pval-0.378_kp*esval) &
1540 + (one-om1) * epss * esval2 / (pval-0.378_kp*esval2)
1541 ENDIF
1542
1543! beta1 = (rgas/rv)*(lstarn1/(rgas*Tl1_1))*(lstarn1/(cp*Tl1_1))
1544! beta1 = (lstarn1*lstarn1*onebrvcp) / (Tl1_1*Tl1_1) ! A.18
1545
1546 beta1 = lstarn1 / tl1_1
1547 beta1 = beta1 * beta1 * onebrvcp
1548
1549
1550! Are the two plumes equal? If so then set qs and beta
1551! in each column to each other to save computation
1552 IF (tl1_1 == tl1_2) THEN
1553 qs2 = qs1
1554 beta2 = beta1
1555 ELSE
1556 IF (tl1_2 >= tbgmax) THEN
1557 lstarn2 = lcond
1558 esval = min(fpvsl(tl1_2), pval)
1559 qs2 = eps * esval / (pval-0.378_kp*esval)
1560 ELSE IF (tl1_2 <= tbgmin) THEN
1561 lstarn2 = lsub
1562 esval = min(fpvsi(tl1_2), pval)
1563 qs2 = epss * esval / (pval-0.378_kp*esval)
1564 ELSE
1565 om2 = max(zero, min(one, a_bg*(tl1_2-tbgmin)))
1566 lstarn2 = lcond + (one-om2)*lfus
1567 esval = min(fpvsl(tl1_2), pval)
1568 esval2 = min(fpvsi(tl1_2), pval)
1569 qs2 = om2 * eps * esval / (pval-0.378_kp*esval) &
1570 + (one-om2) * epss * esval2 / (pval-0.378_kp*esval2)
1571 ENDIF
1572
1573! beta2 = (rgas/rv)*(lstarn2/(rgas*Tl1_2))*(lstarn2/(cp*Tl1_2)) ! A.18
1574! beta2 = (lstarn2*lstarn2*onebrvcp) / (Tl1_2*Tl1_2) ! A.18
1575
1576 beta2 = lstarn2 / tl1_2
1577 beta2 = beta2 * beta2 * onebrvcp
1578
1579
1580 ENDIF
1581
1582 qs1 = qs1 * rhc(i,k)
1583 qs2 = qs2 * rhc(i,k)
1584
1585! Now compute cloud stuff - compute s term
1586
1587 cqt1 = one / (one+beta1*qs1) ! A.19
1588 wrk = qs1 * (one+beta1*qw1_1) * cqt1
1589 s1 = qw1_1 - wrk ! A.17
1590 cthl1 = cqt1*wrk*cpolv*beta1*pkap ! A.20
1591
1592 wrk1 = cthl1 * cthl1
1593 wrk2 = cqt1 * cqt1
1594! std_s1 = sqrt(max(zero,wrk1*thl2_1+wrk2*qw2_1-2.*cthl1*sqrtthl2_1*cqt1*sqrtqw2_1*r_qwthl_1))
1595 std_s1 = sqrt(max(zero, wrk1*thl2_1+wrk2*qw2_1 &
1596 - two*cthl1*sqrtthl2_1*cqt1*sqrtqw2_1*r_qwthl_1))
1597
1598 qn1 = zero
1599 c1 = zero
1600
1601 IF (std_s1 > zero) THEN
1602 wrk = s1 / (std_s1*sqrt2)
1603 c1 = max(zero, min(one, half*(one+erf(wrk)))) ! A.15
1604
1605 IF (c1 > zero) qn1 = s1*c1 + (std_s1*sqrtpii)*exp(-wrk*wrk) ! A.16
1606 ELSEIF (s1 >= qcmin) THEN
1607 c1 = one
1608 qn1 = s1
1609 ENDIF
1610
1611! now compute non-precipitating cloud condensate
1612
1613! If two plumes exactly equal, then just set many of these
1614! variables to themselves to save on computation.
1615 IF (qw1_1 == qw1_2 .and. thl2_1 == thl2_2 .and. qs1 == qs2) THEN
1616 s2 = s1
1617 cthl2 = cthl1
1618 cqt2 = cqt1
1619 std_s2 = std_s1
1620 c2 = c1
1621 qn2 = qn1
1622 ELSE
1623
1624 cqt2 = one / (one+beta2*qs2)
1625 wrk = qs2 * (one+beta2*qw1_2) * cqt2
1626 s2 = qw1_2 - wrk
1627 cthl2 = wrk*cqt2*cpolv*beta2*pkap
1628 wrk1 = cthl2 * cthl2
1629 wrk2 = cqt2 * cqt2
1630! std_s2 = sqrt(max(zero,wrk1*thl2_2+wrk2*qw2_2-2.*cthl2*sqrtthl2_2*cqt2*sqrtqw2_2*r_qwthl_1))
1631 std_s2 = sqrt(max(zero, wrk1*thl2_2+wrk2*qw2_2 &
1632 - two*cthl2*sqrtthl2_2*cqt2*sqrtqw2_2*r_qwthl_1))
1633
1634 qn2 = zero
1635 c2 = zero
1636
1637 IF (std_s2 > zero) THEN
1638 wrk = s2 / (std_s2*sqrt2)
1639 c2 = max(zero, min(one, half*(one+erf(wrk))))
1640 IF (c2 > zero) qn2 = s2*c2 + (std_s2*sqrtpii)*exp(-wrk*wrk)
1641 ELSEIF (s2 >= qcmin) THEN
1642 c2 = one
1643 qn2 = s2
1644 ENDIF
1645
1646 ENDIF
1647
1648! finally, compute the SGS cloud fraction
1649 diag_frac = aterm*c1 + onema*c2
1650
1651 om1 = max(zero, min(one, (tl1_1-tbgmin)*a_bg))
1652 om2 = max(zero, min(one, (tl1_2-tbgmin)*a_bg))
1653
1654 qn1 = min(qn1,qw1_1)
1655 qn2 = min(qn2,qw1_2)
1656
1657 ql1 = qn1*om1
1658 ql2 = qn2*om2
1659
1660 qi1 = qn1 - ql1
1661 qi2 = qn2 - ql2
1662
1663 diag_qn = min(max(zero, aterm*qn1 + onema*qn2), total_water(i,k))
1664 diag_ql = min(max(zero, aterm*ql1 + onema*ql2), diag_qn)
1665 diag_qi = max(zero, diag_qn - diag_ql)
1666
1667! Update temperature variable based on diagnosed cloud properties
1668 om1 = max(zero, min(one, (tabs(i,k)-tbgmin)*a_bg))
1669 lstarn1 = lcond + (one-om1)*lfus
1670 tabs(i,k) = hl(i,k) - gamaz(i,k) + fac_cond*(diag_ql+qpl(i,k)) &
1671 + fac_sub *(diag_qi+qpi(i,k)) &
1672 + tkesbdiss(i,k) * (dtn/cp) ! tke dissipative heating
1673
1674! Update ncpl and ncpi Anning Cheng 03/11/2016
1675! ncpl(i,k) = diag_ql/max(qc(i,k),1.e-10)*ncpl(i,k)
1676
1677! Update ncpl and ncpi Moorthi 12/12/2018
1678 if (ntlnc > 0) then ! liquid and ice number concentrations predicted
1679 if (ncpl(i,k) > nmin) then
1680 ncpl(i,k) = diag_ql/max(qc(i,k),1.0e-10_kp)*ncpl(i,k)
1681 else
1682 ncpl(i,k) = max(diag_ql/(fourb3*pi*rl_cub*997.0_kp), nmin)
1683 endif
1684 if (ncpi(i,k) > nmin) then
1685 ncpi(i,k) = diag_qi/max(qi(i,k),1.0e-10_kp)*ncpi(i,k)
1686 else
1687 ncpi(i,k) = max(diag_qi/(fourb3*pi*ri_cub*500.0_kp), nmin)
1688 endif
1689 endif
1690
1691! Update moisture fields
1692 qc(i,k) = diag_ql
1693 qi(i,k) = diag_qi
1694 qwv(i,k) = max(zero, total_water(i,k) - diag_qn)
1695 cld_sgs(i,k) = diag_frac
1696
1697! Compute the liquid water flux
1698 wqls = aterm * ((w1_1-w_first)*ql1) + onema * ((w1_2-w_first)*ql2)
1699 wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2)
1700
1701! Compute statistics for the fluxes so we don't have to save these variables
1702 wqlsb(k) = wqlsb(k) + wqls
1703 wqisb(k) = wqisb(k) + wqis
1704
1705! diagnostic buoyancy flux. Includes effects from liquid water, ice
1706! condensate, liquid & ice precipitation
1707! wrk = epsv * basetemp
1708 wrk = epsv * thv(i,k)
1709
1710 bastoeps = onebeps * thv(i,k)
1711
1712 if (k < nzm) then
1713 wthv_sec(i,k) = wthlsec + wrk*wqwsec &
1714 + (fac_cond-bastoeps)*wqls &
1715 + (fac_sub-bastoeps) *wqis &
1716 + ((lstarn1/cp)-thv(i,k))*half*(wqp_sec(i,kd)+wqp_sec(i,ku))
1717 else
1718 wthv_sec(i,k) = wthlsec + wrk*wqwsec &
1719 + (fac_cond-bastoeps)*wqls &
1720 + (fac_sub-bastoeps) *wqis &
1721 + ((lstarn1/cp)-thv(i,k))*half*wqp_sec(i,k)
1722 endif
1723
1724! wthv_sec(i,k) = wthlsec + wrk*wqwsec &
1725! + (fac_cond-bastoeps)*wqls &
1726! + (fac_sub-bastoeps)*wqis &
1727! + ((lstarn1/cp)-basetemp)*half*(wqp_sec(i,kd)+wqp_sec(i,ku))
1728
1729 ENDDO
1730 ENDDO
1731
1732
1733 end subroutine assumed_pdf
1734
1735
1736! Saturation vapor pressure and mixing ratio subroutines
1737! Based on Flatau et al (1992), J. App. Met., 31, 1507-1513
1738! Code by Marat Khairoutdinov
1739
1740
1741 real function esatw(t)
1742 real t ! temperature (K)
1743 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1744 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1745 6.11239921, 0.443987641, 0.142986287e-1, &
1746 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
1747 0.640689451e-10, -0.952447341e-13,-0.976195544e-15/
1748 real dt
1749 dt = max(-80.,t-273.16)
1750 esatw = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1751 end function esatw
1752
1753 real function qsatw(t,p)
1754! implicit none
1755 real t ! temperature (K)
1756 real p ! pressure (Pa)
1757 real esat
1758! esat = fpvs(t)
1759 esat = fpvsl(t)
1760 qsatw = 0.622 * esat/max(esat,p-0.378*esat)
1761! esat = esatw(t)
1762! qsatw = 0.622 * esat/max(esat,p-esat)
1763 end function qsatw
1764
1765
1766 real function esati(t)
1767 real t ! temperature (K)
1768 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1769 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1770 6.11147274, 0.503160820, 0.188439774e-1, &
1771 0.420895665e-3, 0.615021634e-5, 0.602588177e-7, &
1772 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
1773 real dt
1774! real esatw
1775 if(t > 273.15) then
1776 esati = esatw(t)
1777 else if(t.gt.185.) then
1778 dt = t-273.16
1779 esati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1780 else ! use some additional interpolation below 184K
1781 dt = max(-100.,t-273.16)
1782 esati = 0.00763685 + dt*(0.000151069+dt*7.48215e-07)
1783 endif
1784 end function esati
1785
1786 real function qsati(t,p)
1787 real t ! temperature (K)
1788 real p ! pressure (Pa)
1789 real esat !,esati
1790! esat = fpvs(t)
1791 esat = fpvsi(t)
1792 qsati = 0.622 * esat/max(esat,p-0.378*esat)
1793! esat = esati(t)
1794! qsati = 0.622 * esat/max(esat,p-esat)
1795 end function qsati
1796
1797 real function dtesatw(t)
1798 real t ! temperature (K)
1799 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1800 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1801 0.443956472, 0.285976452e-1, 0.794747212e-3, &
1802 0.121167162e-4, 0.103167413e-6, 0.385208005e-9, &
1803 -0.604119582e-12, -0.792933209e-14, -0.599634321e-17/
1804 real dt
1805 dt = max(-80.,t-273.16)
1806 dtesatw = a0 + dt* (a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1807 end function dtesatw
1808
1809 real function dtqsatw(t,p)
1810 real t ! temperature (K)
1811 real p ! pressure (Pa)
1812! real dtesatw
1813 dtqsatw = 100.0*0.622*dtesatw(t)/p
1814 end function dtqsatw
1815
1816 real function dtesati(t)
1817 real t ! temperature (K)
1818 real a0,a1,a2,a3,a4,a5,a6,a7,a8
1819 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / &
1820 0.503223089, 0.377174432e-1, 0.126710138e-2, &
1821 0.249065913e-4, 0.312668753e-6, 0.255653718e-8, &
1822 0.132073448e-10, 0.390204672e-13, 0.497275778e-16/
1823 real dt
1824! real dtesatw
1825 if(t > 273.15) then
1826 dtesati = dtesatw(t)
1827 else if(t > 185.) then
1828 dt = t-273.16
1829 dtesati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
1830 else ! use additional interpolation below 185K
1831 dt = max(-100.,t-273.16)
1832 dtesati = 0.0013186 + dt*(2.60269e-05+dt*1.28676e-07)
1833 endif
1834 end function dtesati
1835
1836
1837 real function dtqsati(t,p)
1838 real t ! temperature (K)
1839 real p ! pressure (Pa)
1840! real dtesati
1841 dtqsati = 100.0*0.622*dtesati(t)/p
1842 end function dtqsati
1843
1844end subroutine shoc_work
1845
1846end module shoc
subroutine esat(t, esw, esi, desw, desi)
use polynomials to calculate saturation vapor pressure and derivative with respect to temperature: ov...
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module contains the CCPP-compliant SHOC scheme.
Definition shoc.F90:5