| blob | 561b461089ec6a949d34d479e882a234ae59341a |
1 #define WRF_PORT
2 #define MODAL_AERO
3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
5 module eddy_diff
7 !--------------------------------------------------------------------------------- !
8 ! !
9 ! The University of Washington Moist Turbulence Scheme to compute eddy diffusion !
10 ! coefficients associated with dry and moist turbulences in the whole !
11 ! atmospheric layers. !
12 ! !
13 ! For detailed description of the code and its performances, see !
14 ! !
15 ! 1.'A new moist turbulence parametrization in the Community Atmosphere Model' !
16 ! by Christopher S. Bretherton and Sungsu Park. J. Climate. 2009. 22. 3422-3448 !
17 ! 2.'The University of Washington shallow convection and moist turbulence schemes !
18 ! and their impact on climate simulations with the Community Atmosphere Model' !
19 ! by Sungsu Park and Christopher S. Bretherton. J. Climate. 2009. 22. 3449-3469 !
20 ! !
21 ! For questions on the scheme and code, send an email to !
22 ! Sungsu Park at sungsup@ucar.edu (tel: 303-497-1375) !
23 ! Chris Bretherton at breth@washington.edu !
24 ! !
25 ! Developed by Chris Bretherton at the University of Washington, Seattle, WA. !
26 ! Sungsu Park at the CGD/NCAR, Boulder, CO. !
27 ! Last coded on May.2006, Dec.2009 by Sungsu Park. !
28 ! !
29 !--------------------------------------------------------------------------------- !
31 use diffusion_solver, only : vdiff_selector
32 #ifndef WRF_PORT
33 use cam_history, only : outfld, addfld, phys_decomp
34 use cam_logfile, only : iulog
35 use ppgrid, only : pver
36 #else
37 use module_cam_support, only: iulog, pver,outfld, addfld, phys_decomp
38 #endif
40 implicit none
41 private
42 save
44 public init_eddy_diff
45 public compute_eddy_diff
47 type(vdiff_selector) :: fieldlist_wet ! Logical switches for moist mixing ratio diffusion
48 type(vdiff_selector) :: fieldlist_dry ! Logical switches for dry mixing ratio diffusion
49 integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
51 ! --------------------------------- !
52 ! PBL Parameters used in the UW PBL !
53 ! --------------------------------- !
55 character, parameter :: sftype = 'l' ! Method for calculating saturation fraction
57 character(len=4), parameter :: choice_evhc = 'maxi' ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_radf
58 character(len=6), parameter :: choice_radf = 'maxi' ! 'orig', 'ramp', 'maxi' : recommended to be used with choice_evhc
59 character(len=6), parameter :: choice_SRCL = 'nonamb' ! 'origin', 'remove', 'nonamb'
61 character(len=6), parameter :: choice_tunl = 'rampcl' ! 'origin', 'rampsl'(Sungsu), 'rampcl'(Chris)
62 real(r8), parameter :: ctunl = 2._r8 ! Maximum asympt leng = ctunl*tunl when choice_tunl = 'rampsl(cl)' [ no unit ]
63 character(len=6), parameter :: choice_leng = 'origin' ! 'origin', 'takemn'
64 real(r8), parameter :: cleng = 3._r8 ! Order of 'leng' when choice_leng = 'origin' [ no unit ]
65 character(len=6), parameter :: choice_tkes = 'ibprod' ! 'ibprod' (include tkes in computing bprod), 'ebprod'(exclude)
67 ! Parameters for 'sedimenttaion-entrainment feedback' for liquid stratus
68 ! If .false., no sedimentation entrainment feedback ( i.e., use default evhc )
70 logical, parameter :: id_sedfact = .false.
71 real(r8), parameter :: ased = 9._r8 ! Valid only when id_sedfact = .true.
73 ! --------------------------------------------------------------------------------------------------- !
74 ! Parameters governing entrainment efficiency A = a1l(i)*evhc, evhc = 1 + a2l * a3l * L * ql / jt2slv !
75 ! Here, 'ql' is cloud-top LWC and 'jt2slv' is the jump in 'slv' across !
76 ! the cloud-top entrainment zone ( across two grid layers to consider full mixture ) !
77 ! --------------------------------------------------------------------------------------------------- !
79 real(r8), parameter :: a1l = 0.10_r8 ! Dry entrainment efficiency for TKE closure
80 ! a1l = 0.2*tunl*erat^-1.5, where erat = <e>/wstar^2 for dry CBL = 0.3.
82 real(r8), parameter :: a1i = 0.2_r8 ! Dry entrainment efficiency for wstar closure
83 real(r8), parameter :: ccrit = 0.5_r8 ! Minimum allowable sqrt(tke)/wstar. Used in solving cubic equation for 'ebrk'
84 real(r8), parameter :: wstar3factcrit = 0.5_r8 ! 1/wstar3factcrit is the maximally allowed enhancement of 'wstar3' due to entrainment.
86 real(r8), parameter :: a2l = 30._r8 ! Moist entrainment enhancement param (recommended range : 10~30 )
87 real(r8), parameter :: a3l = 0.8_r8 ! Approximation to a complicated thermodynamic parameters
89 real(r8), parameter :: jbumin = .001_r8 ! Minimum buoyancy jump at an entrainment jump, [m/s2]
90 real(r8), parameter :: evhcmax = 10._r8 ! Upper limit of evaporative enhancement factor
92 real(r8), parameter :: ustar_min = 0.01_r8 ! Minimum permitted value of ustar [ m/s ]
93 real(r8), parameter :: onet = 1._r8/3._r8 ! 1/3 power in wind gradient expression [ no unit ]
94 #ifndef WRF_PORT
95 !pver is not a parameter in cam_support module in WRF, therefore ncvmax cannot be equated to pver as parameter
96 integer, parameter :: ncvmax = pver ! Max numbers of CLs (good to set to 'pver')
97 #else
98 integer :: ncvmax
99 #endif
100 real(r8), parameter :: qmin = 1.e-5_r8 ! Minimum grid-mean LWC counted as clouds [kg/kg]
101 real(r8), parameter :: ntzero = 1.e-12_r8 ! Not zero (small positive number used in 's2')
102 real(r8), parameter :: b1 = 5.8_r8 ! TKE dissipation D = e^3/(b1*leng), e = b1*W.
103 real(r8) :: b123 ! b1**(2/3)
104 real(r8), parameter :: tunl = 0.085_r8 ! Asympt leng = tunl*(turb lay depth)
105 real(r8), parameter :: alph1 = 0.5562_r8 ! alph1~alph5 : Galperin instability function parameters
106 real(r8), parameter :: alph2 = -4.3640_r8 ! These coefficients are used to calculate
107 real(r8), parameter :: alph3 = -34.6764_r8 ! 'sh' and 'sm' from 'gh'.
108 real(r8), parameter :: alph4 = -6.1272_r8 !
109 real(r8), parameter :: alph5 = 0.6986_r8 !
110 real(r8), parameter :: ricrit = 0.19_r8 ! Critical Richardson number for turbulence. Can be any value >= 0.19.
111 real(r8), parameter :: ae = 1._r8 ! TKE transport efficiency [no unit]
112 real(r8), parameter :: rinc = -0.04_r8 ! Minimum W/<W> used for CL merging test
113 real(r8), parameter :: wpertmin = 1.e-6_r8 ! Minimum PBL eddy vertical velocity perturbation
114 real(r8), parameter :: wfac = 1._r8 ! Ratio of 'wpert' to sqrt(tke) for CL.
115 real(r8), parameter :: tfac = 1._r8 ! Ratio of 'tpert' to (w't')/wpert for CL. Same ratio also used for q
116 real(r8), parameter :: fak = 8.5_r8 ! Constant in surface temperature excess for stable STL. [ no unit ]
117 real(r8), parameter :: rcapmin = 0.1_r8 ! Minimum allowable e/<e> in a CL
118 real(r8), parameter :: rcapmax = 2.0_r8 ! Maximum allowable e/<e> in a CL
119 real(r8), parameter :: tkemax = 20._r8 ! TKE is capped at tkemax [m2/s2]
120 real(r8), parameter :: lambda = 0.5_r8 ! Under-relaxation factor ( 0 < lambda =< 1 )
122 logical, parameter :: use_kvf = .false. ! .true. (.false.) : initialize kvh/kvm = kvf ( 0. )
123 logical, parameter :: use_dw_surf = .true. ! Used in 'zisocl'. Default is 'true'
124 ! If 'true', surface interfacial energy does not contribute to the CL mean
125 ! stbility functions after finishing merging. For this case,
126 ! 'dl2n2_surf' is only used for a merging test based on 'l2n2'
127 ! If 'false',surface interfacial enery explicitly contribute to CL mean
128 ! stability functions after finishing merging. For this case,
129 ! 'dl2n2_surf' and 'dl2s2_surf' are directly used for calculating
130 ! surface interfacial layer energetics
132 logical, parameter :: set_qrlzero = .false. ! .true. ( .false.) : turning-off ( on) radiative-turbulence interaction by setting qrl = 0.
134 ! ------------------------------------- !
135 ! PBL Parameters not used in the UW PBL !
136 ! ------------------------------------- !
138 real(r8), parameter :: pblmaxp = 4.e4_r8 ! PBL max depth in pressure units.
139 real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri).
140 real(r8), parameter :: betam = 15.0_r8 ! Constant in wind gradient expression.
141 real(r8), parameter :: betas = 5.0_r8 ! Constant in surface layer gradient expression.
142 real(r8), parameter :: betah = 15.0_r8 ! Constant in temperature gradient expression.
143 real(r8), parameter :: fakn = 7.2_r8 ! Constant in turbulent prandtl number.
144 real(r8), parameter :: ricr = 0.3_r8 ! Critical richardson number.
145 real(r8), parameter :: sffrac = 0.1_r8 ! Surface layer fraction of boundary layer
146 real(r8), parameter :: binm = betam*sffrac ! betam * sffrac
147 real(r8), parameter :: binh = betah*sffrac ! betah * sffrac
149 ! ------------------------------------------------------- !
150 ! PBL constants set using values from other parts of code !
151 ! ------------------------------------------------------- !
153 real(r8) :: cpair ! Specific heat of dry air
154 real(r8) :: rair ! Gas const for dry air
155 real(r8) :: zvir ! rh2o/rair - 1
156 real(r8) :: latvap ! Latent heat of vaporization
157 real(r8) :: latice ! Latent heat of fusion
158 real(r8) :: latsub ! Latent heat of sublimation
159 real(r8) :: g ! Gravitational acceleration
160 real(r8) :: vk ! Von Karman's constant
161 real(r8) :: ccon ! fak * sffrac * vk
163 integer :: ntop_turb ! Top interface level to which turbulent vertical diffusion is applied ( = 1 )
164 integer :: nbot_turb ! Bottom interface level to which turbulent vertical diff is applied ( = pver )
166 real(r8), allocatable :: ml2(:) ! Mixing lengths squared. Not used in the UW PBL. Used for computing free air diffusivity.
168 CONTAINS
170 !============================================================================ !
171 ! !
172 !============================================================================ !
174 subroutine init_eddy_diff( kind, pver, gravx, cpairx, rairx, zvirx, &
175 latvapx, laticex, ntop_eddy, nbot_eddy, hypm, vkx )
176 !---------------------------------------------------------------- !
177 ! Purpose: !
178 ! Initialize time independent constants/variables of PBL package. !
179 !---------------------------------------------------------------- !
180 use diffusion_solver, only: init_vdiff, vdiff_select
181 #ifndef WRF_PORT
182 use cam_history, only: outfld, addfld, phys_decomp
183 #else
184 use module_cam_support, only: outfld, addfld, phys_decomp
185 #endif
186 implicit none
187 ! --------- !
188 ! Arguments !
189 ! --------- !
190 integer, intent(in) :: kind ! Kind of reals being passed in
191 integer, intent(in) :: pver ! Number of vertical layers
192 integer, intent(in) :: ntop_eddy ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
193 integer, intent(in) :: nbot_eddy ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )
194 real(r8), intent(in) :: gravx ! Acceleration of gravity
195 real(r8), intent(in) :: cpairx ! Specific heat of dry air
196 real(r8), intent(in) :: rairx ! Gas constant for dry air
197 real(r8), intent(in) :: zvirx ! rh2o/rair - 1
198 real(r8), intent(in) :: latvapx ! Latent heat of vaporization
199 real(r8), intent(in) :: laticex ! Latent heat of fusion
200 real(r8), intent(in) :: hypm(pver) ! Reference pressures at midpoints
201 real(r8), intent(in) :: vkx ! Von Karman's constant
203 character(128) :: errstring ! Error status for init_vdiff
204 integer :: k ! Vertical loop index
205 #ifdef WRF_PORT
206 !As pver is not parameter in the cam_support module in WRF, a value to ncvmax is given here
207 ncvmax = pver ! Max numbers of CLs (good to set to 'pver')
208 #endif
209 if( kind .ne. r8 ) then
210 write(iulog,*) 'wrong KIND of reals passed to init_diffusvity -- exiting.'
211 #ifdef WRF_PORT
212 call wrf_message(iulog)
213 #endif
214 stop 'init_eddy_diff'
215 endif
217 ! --------------- !
218 ! Basic constants !
219 ! --------------- !
221 cpair = cpairx
222 rair = rairx
223 g = gravx
224 zvir = zvirx
225 latvap = latvapx
226 latice = laticex
227 latsub = latvap + latice
228 vk = vkx
229 ccon = fak*sffrac*vk
230 ntop_turb = ntop_eddy
231 nbot_turb = nbot_eddy
232 b123 = b1**(2._r8/3._r8)
234 ! Set the square of the mixing lengths. Only for CAM3 HB PBL scheme.
235 ! Not used for UW moist PBL. Used for free air eddy diffusivity.
236 #ifndef WRF_PORT
237 ! ml2 is never used currently as use_kvf is set to false
238 allocate(ml2(pver+1))
239 ml2(1:ntop_turb) = 0._r8
240 do k = ntop_turb + 1, nbot_turb
241 ml2(k) = 30.0_r8**2
242 end do
243 ml2(nbot_turb+1:pver+1) = 0._r8
244 #endif
246 ! Initialize diffusion solver module
248 call init_vdiff(r8, 1, rair, g, fieldlist_wet, fieldlist_dry, errstring)
250 ! Select the fields which will be diffused
252 if(vdiff_select(fieldlist_wet,'s').ne.'') write(iulog,*) 'error: ', vdiff_select(fieldlist_wet,'s')
253 #ifdef WRF_PORT
254 call wrf_message(iulog)
255 #endif
256 if(vdiff_select(fieldlist_wet,'q',1).ne.'') write(iulog,*) 'error: ', vdiff_select(fieldlist_wet,'q',1)
257 #ifdef WRF_PORT
258 call wrf_message(iulog)
259 #endif
260 if(vdiff_select(fieldlist_wet,'u').ne.'') write(iulog,*) 'error: ', vdiff_select(fieldlist_wet,'u')
261 #ifdef WRF_PORT
262 call wrf_message(iulog)
263 #endif
264 if(vdiff_select(fieldlist_wet,'v').ne.'') write(iulog,*) 'error: ', vdiff_select(fieldlist_wet,'v')
265 #ifdef WRF_PORT
266 call wrf_message(iulog)
267 #endif
269 ! ------------------------------------------------------------------- !
270 ! Writing outputs for detailed analysis of UW moist turbulence scheme !
271 ! ------------------------------------------------------------------- !
273 call addfld('UW_errorPBL', 'm2/s', 1, 'A', 'Error function of UW PBL', phys_decomp )
274 call addfld('UW_n2', 's-2', pver, 'A', 'Buoyancy Frequency, LI', phys_decomp )
275 call addfld('UW_s2', 's-2', pver, 'A', 'Shear Frequency, LI', phys_decomp )
276 call addfld('UW_ri', 'no', pver, 'A', 'Interface Richardson Number, I', phys_decomp )
277 call addfld('UW_sfuh', 'no', pver, 'A', 'Upper-Half Saturation Fraction, L', phys_decomp )
278 call addfld('UW_sflh', 'no', pver, 'A', 'Lower-Half Saturation Fraction, L', phys_decomp )
279 call addfld('UW_sfi', 'no', pver+1, 'A', 'Interface Saturation Fraction, I', phys_decomp )
280 call addfld('UW_cldn', 'no', pver, 'A', 'Cloud Fraction, L', phys_decomp )
281 call addfld('UW_qrl', 'g*W/m2', pver, 'A', 'LW cooling rate, L', phys_decomp )
282 call addfld('UW_ql', 'kg/kg', pver, 'A', 'ql(LWC), L', phys_decomp )
283 call addfld('UW_chu', 'g*kg/J', pver+1, 'A', 'Buoyancy Coefficient, chu, I', phys_decomp )
284 call addfld('UW_chs', 'g*kg/J', pver+1, 'A', 'Buoyancy Coefficient, chs, I', phys_decomp )
285 call addfld('UW_cmu', 'g/kg/kg', pver+1, 'A', 'Buoyancy Coefficient, cmu, I', phys_decomp )
286 call addfld('UW_cms', 'g/kg/kg', pver+1, 'A', 'Buoyancy Coefficient, cms, I', phys_decomp )
287 call addfld('UW_tke', 'm2/s2', pver+1, 'A', 'TKE, I', phys_decomp )
288 call addfld('UW_wcap', 'm2/s2', pver+1, 'A', 'Wcap, I', phys_decomp )
289 call addfld('UW_bprod', 'm2/s3', pver+1, 'A', 'Buoyancy production, I', phys_decomp )
290 call addfld('UW_sprod', 'm2/s3', pver+1, 'A', 'Shear production, I', phys_decomp )
291 call addfld('UW_kvh', 'm2/s', pver+1, 'A', 'Eddy diffusivity of heat, I', phys_decomp )
292 call addfld('UW_kvm', 'm2/s', pver+1, 'A', 'Eddy diffusivity of uv, I', phys_decomp )
293 call addfld('UW_pblh', 'm', 1, 'A', 'PBLH, 1', phys_decomp )
294 call addfld('UW_pblhp', 'Pa', 1, 'A', 'PBLH pressure, 1', phys_decomp )
295 call addfld('UW_tpert', 'K', 1, 'A', 'Convective T excess, 1', phys_decomp )
296 call addfld('UW_qpert', 'kg/kg', 1, 'A', 'Convective qt excess, I', phys_decomp )
297 call addfld('UW_wpert', 'm/s', 1, 'A', 'Convective W excess, I', phys_decomp )
298 call addfld('UW_ustar', 'm/s', 1, 'A', 'Surface Frictional Velocity, 1', phys_decomp )
299 call addfld('UW_tkes', 'm2/s2', 1, 'A', 'Surface TKE, 1', phys_decomp )
300 call addfld('UW_minpblh', 'm', 1, 'A', 'Minimum PBLH, 1', phys_decomp )
301 call addfld('UW_turbtype', 'no', pver+1, 'A', 'Interface Turbulence Type, I', phys_decomp )
302 call addfld('UW_kbase_o', 'no', ncvmax, 'A', 'Initial CL Base Exterbal Interface Index, CL', phys_decomp )
303 call addfld('UW_ktop_o', 'no', ncvmax, 'A', 'Initial Top Exterbal Interface Index, CL', phys_decomp )
304 call addfld('UW_ncvfin_o', '#', 1, 'A', 'Initial Total Number of CL regimes, CL', phys_decomp )
305 call addfld('UW_kbase_mg', 'no', ncvmax, 'A', 'kbase after merging, CL', phys_decomp )
306 call addfld('UW_ktop_mg', 'no', ncvmax, 'A', 'ktop after merging, CL', phys_decomp )
307 call addfld('UW_ncvfin_mg', '#', 1, 'A', 'ncvfin after merging, CL', phys_decomp )
308 call addfld('UW_kbase_f', 'no', ncvmax, 'A', 'Final kbase with SRCL, CL', phys_decomp )
309 call addfld('UW_ktop_f', 'no', ncvmax, 'A', 'Final ktop with SRCL, CL', phys_decomp )
310 call addfld('UW_ncvfin_f', '#', 1, 'A', 'Final ncvfin with SRCL, CL', phys_decomp )
311 call addfld('UW_wet', 'm/s', ncvmax, 'A', 'Entrainment rate at CL top, CL', phys_decomp )
312 call addfld('UW_web', 'm/s', ncvmax, 'A', 'Entrainment rate at CL base, CL', phys_decomp )
313 call addfld('UW_jtbu', 'm/s2', ncvmax, 'A', 'Buoyancy jump across CL top, CL', phys_decomp )
314 call addfld('UW_jbbu', 'm/s2', ncvmax, 'A', 'Buoyancy jump across CL base, CL', phys_decomp )
315 call addfld('UW_evhc', 'no', ncvmax, 'A', 'Evaporative enhancement factor, CL', phys_decomp )
316 call addfld('UW_jt2slv', 'J/kg', ncvmax, 'A', 'slv jump for evhc, CL', phys_decomp )
317 call addfld('UW_n2ht', 's-2', ncvmax, 'A', 'n2 at just below CL top interface, CL', phys_decomp )
318 call addfld('UW_n2hb', 's-2', ncvmax, 'A', 'n2 at just above CL base interface', phys_decomp )
319 call addfld('UW_lwp', 'kg/m2', ncvmax, 'A', 'LWP in the CL top layer, CL', phys_decomp )
320 call addfld('UW_optdepth', 'no', ncvmax, 'A', 'Optical depth of the CL top layer, CL', phys_decomp )
321 call addfld('UW_radfrac', 'no', ncvmax, 'A', 'Fraction of radiative cooling confined in the CL top', phys_decomp )
322 call addfld('UW_radf', 'm2/s3', ncvmax, 'A', 'Buoyancy production at the CL top by radf, I', phys_decomp )
323 call addfld('UW_wstar', 'm/s', ncvmax, 'A', 'Convective velocity, Wstar, CL', phys_decomp )
324 call addfld('UW_wstar3fact', 'no', ncvmax, 'A', 'Enhancement of wstar3 due to entrainment, CL', phys_decomp )
325 call addfld('UW_ebrk', 'm2/s2', ncvmax, 'A', 'CL-averaged TKE, CL', phys_decomp )
326 call addfld('UW_wbrk', 'm2/s2', ncvmax, 'A', 'CL-averaged W, CL', phys_decomp )
327 call addfld('UW_lbrk', 'm', ncvmax, 'A', 'CL internal thickness, CL', phys_decomp )
328 call addfld('UW_ricl', 'no', ncvmax, 'A', 'CL-averaged Ri, CL', phys_decomp )
329 call addfld('UW_ghcl', 'no', ncvmax, 'A', 'CL-averaged gh, CL', phys_decomp )
330 call addfld('UW_shcl', 'no', ncvmax, 'A', 'CL-averaged sh, CL', phys_decomp )
331 call addfld('UW_smcl', 'no', ncvmax, 'A', 'CL-averaged sm, CL', phys_decomp )
332 call addfld('UW_gh', 'no', pver+1, 'A', 'gh at all interfaces, I', phys_decomp )
333 call addfld('UW_sh', 'no', pver+1, 'A', 'sh at all interfaces, I', phys_decomp )
334 call addfld('UW_sm', 'no', pver+1, 'A', 'sm at all interfaces, I', phys_decomp )
335 call addfld('UW_ria', 'no', pver+1, 'A', 'ri at all interfaces, I', phys_decomp )
336 call addfld('UW_leng', 'm/s', pver+1, 'A', 'Turbulence length scale, I', phys_decomp )
338 return
340 end subroutine init_eddy_diff
342 !=============================================================================== !
343 ! !
344 !=============================================================================== !
346 subroutine compute_eddy_diff( lchnk , &
347 pcols , pver , ncol , t , qv , ztodt , &
348 ql , qi , s , rpdel , cldn , qrl , wsedl , &
349 z , zi , pmid , pi , u , v , &
350 taux , tauy , shflx , qflx , wstarent , nturb , &
351 ustar , pblh , kvm_in , kvh_in , kvm_out , kvh_out , kvq , &
352 cgh , cgs , tpert , qpert , wpert , tke , bprod , &
353 sprod , sfi , qsat , kvinit , &
354 tauresx, tauresy, ksrftms , &
355 ipbl , kpblh , wstarPBL , turbtype, sm_aw )
357 !-------------------------------------------------------------------- !
358 ! Purpose: Interface to compute eddy diffusivities. !
359 ! Eddy diffusivities are calculated in a fully implicit way !
360 ! through iteration process. !
361 ! Author: Sungsu Park. August. 2006. !
362 ! May. 2008. !
363 !-------------------------------------------------------------------- !
364 use diffusion_solver, only: compute_vdiff
365 #ifndef WRF_PORT
366 use cam_history, only: outfld, addfld, phys_decomp
367 ! use physics_types, only: physics_state
368 use phys_debug_util, only: phys_debug_col
369 use time_manager, only: is_first_step, get_nstep
370 #else
371 use module_cam_support, only: outfld, addfld, phys_decomp
372 #endif
374 implicit none
376 ! type(physics_state) :: state ! Physics state variables
378 ! --------------- !
379 ! Input Variables !
380 ! --------------- !
382 integer, intent(in) :: lchnk
383 integer, intent(in) :: pcols ! Number of atmospheric columns [ # ]
384 integer, intent(in) :: pver ! Number of atmospheric layers [ # ]
385 integer, intent(in) :: ncol ! Number of atmospheric columns [ # ]
386 integer, intent(in) :: nturb ! Number of iteration steps for calculating eddy diffusivity [ # ]
387 logical, intent(in) :: wstarent ! .true. means use the 'wstar' entrainment closure.
388 logical, intent(in) :: kvinit ! 'true' means time step = 1 : used for initializing kvh, kvm (uses kvf or zero)
389 real(r8), intent(in) :: ztodt ! Physics integration time step 2 delta-t [ s ]
390 real(r8), intent(in) :: t(pcols,pver) ! Temperature [K]
391 real(r8), intent(in) :: qv(pcols,pver) ! Water vapor specific humidity [ kg/kg ]
392 real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ]
393 real(r8), intent(in) :: qi(pcols,pver) ! Ice specific humidity [ kg/kg ]
394 real(r8), intent(in) :: s(pcols,pver) ! Dry static energy [ J/kg ]
395 real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel where 'pdel' is thickness of the layer [ Pa ]
396 real(r8), intent(in) :: cldn(pcols,pver) ! Stratiform cloud fraction [ fraction ]
397 real(r8), intent(in) :: qrl(pcols,pver) ! LW cooling rate
398 real(r8), intent(in) :: wsedl(pcols,pver) ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
399 real(r8), intent(in) :: z(pcols,pver) ! Layer mid-point height above surface [ m ]
400 real(r8), intent(in) :: zi(pcols,pver+1) ! Interface height above surface [ m ]
401 real(r8), intent(in) :: pmid(pcols,pver) ! Layer mid-point pressure [ Pa ]
402 real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressure [ Pa ]
403 real(r8), intent(in) :: u(pcols,pver) ! Zonal velocity [ m/s ]
404 real(r8), intent(in) :: v(pcols,pver) ! Meridional velocity [ m/s ]
405 real(r8), intent(in) :: taux(pcols) ! Zonal wind stress at surface [ N/m2 ]
406 real(r8), intent(in) :: tauy(pcols) ! Meridional wind stress at surface [ N/m2 ]
407 real(r8), intent(in) :: shflx(pcols) ! Sensible heat flux at surface [ unit ? ]
408 real(r8), intent(in) :: qflx(pcols) ! Water vapor flux at surface [ unit ? ]
409 real(r8), intent(in) :: kvm_in(pcols,pver+1) ! kvm saved from last timestep [ m2/s ]
410 real(r8), intent(in) :: kvh_in(pcols,pver+1) ! kvh saved from last timestep [ m2/s ]
411 real(r8), intent(in) :: ksrftms(pcols) ! Surface drag coefficient of turbulent mountain stress [ unit ? ]
413 ! ---------------- !
414 ! Output Variables !
415 ! ---------------- !
417 real(r8), intent(out) :: kvm_out(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
418 real(r8), intent(out) :: kvh_out(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
419 real(r8), intent(out) :: kvq(pcols,pver+1) ! Eddy diffusivity for constituents, moisture and tracers [ m2/s ] (note not having '_out')
420 real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ]
421 real(r8), intent(out) :: pblh(pcols) ! PBL top height [ m ]
422 real(r8), intent(out) :: cgh(pcols,pver+1) ! Counter-gradient term for heat [ J/kg/m ]
423 real(r8), intent(out) :: cgs(pcols,pver+1) ! Counter-gradient star [ cg/flux ]
424 real(r8), intent(out) :: tpert(pcols) ! Convective temperature excess [ K ]
425 real(r8), intent(out) :: qpert(pcols) ! Convective humidity excess [ kg/kg ]
426 real(r8), intent(out) :: wpert(pcols) ! Turbulent velocity excess [ m/s ]
427 real(r8), intent(out) :: tke(pcols,pver+1) ! Turbulent kinetic energy [ m2/s2 ]
428 real(r8), intent(out) :: bprod(pcols,pver+1) ! Buoyancy production [ m2/s3 ]
429 real(r8), intent(out) :: sprod(pcols,pver+1) ! Shear production [ m2/s3 ]
430 real(r8), intent(out) :: sfi(pcols,pver+1) ! Interfacial layer saturation fraction [ fraction ]
431 real(r8), intent(out) :: turbtype(pcols,pver+1) ! Turbulence type identifier at all interfaces [ no unit ]
432 real(r8), intent(out) :: sm_aw(pcols,pver+1) ! Normalized Galperin instability function for momentum [ no unit ]
433 ! This is 1 when neutral condition (Ri=0), 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19.
434 real(r8), intent(out) :: ipbl(pcols) ! If 1, PBL is CL, while if 0, PBL is STL.
435 real(r8), intent(out) :: kpblh(pcols) ! Layer index containing PBL top within or at the base interface
436 real(r8), intent(out) :: wstarPBL(pcols) ! Convective velocity within PBL [ m/s ]
438 ! ---------------------- !
439 ! Input-Output Variables !
440 ! ---------------------- !
442 real(r8), intent(inout) :: tauresx(pcols) ! Residual stress to be added in vdiff to correct for turb
443 real(r8), intent(inout) :: tauresy(pcols) ! Stress mismatch between sfc and atm accumulated in prior timesteps
445 ! --------------- !
446 ! Local Variables !
447 ! --------------- !
449 integer icol
450 integer i, k, iturb, status
451 integer, external :: qsat
452 character(128) :: errstring ! Error status for compute_vdiff
454 real(r8) :: tautotx(pcols) ! Total stress including tms
455 real(r8) :: tautoty(pcols) ! Total stress including tms
456 real(r8) :: kvf(pcols,pver+1) ! Free atmospheric eddy diffusivity [ m2/s ]
457 real(r8) :: kvm(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
458 real(r8) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
459 real(r8) :: kvm_preo(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
460 real(r8) :: kvh_preo(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
461 real(r8) :: kvm_pre(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
462 real(r8) :: kvh_pre(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
463 real(r8) :: errorPBL(pcols) ! Error function showing whether PBL produced convergent solution or not. [ unit ? ]
464 real(r8) :: s2(pcols,pver) ! Shear squared, defined at interfaces except surface [ s-2 ]
465 real(r8) :: n2(pcols,pver) ! Buoyancy frequency, defined at interfaces except surface [ s-2 ]
466 real(r8) :: ri(pcols,pver) ! Richardson number, 'n2/s2', defined at interfaces except surface [ s-2 ]
467 real(r8) :: pblhp(pcols) ! PBL top pressure [ Pa ]
468 real(r8) :: minpblh(pcols) ! Minimum PBL height based on surface stress
470 real(r8) :: qt(pcols,pver) ! Total specific humidity [ kg/kg ]
471 real(r8) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ]
472 real(r8) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ]
473 real(r8) :: sl(pcols,pver) ! Liquid water static energy [ J/kg ]
474 real(r8) :: slv(pcols,pver) ! Liquid water virtual static energy [ J/kg ]
475 real(r8) :: slslope(pcols,pver) ! Slope of 'sl' in each layer
476 real(r8) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer
477 real(r8) :: rrho(pcols) ! Density at the lowest layer
478 real(r8) :: qvfd(pcols,pver) ! Specific humidity for diffusion [ kg/kg ]
479 real(r8) :: tfd(pcols,pver) ! Temperature for diffusion [ K ]
480 real(r8) :: slfd(pcols,pver) ! Liquid static energy [ J/kg ]
481 real(r8) :: qtfd(pcols,pver) ! Total specific humidity [ kg/kg ]
482 real(r8) :: qlfd(pcols,pver) ! Liquid water specific humidity for diffusion [ kg/kg ]
483 real(r8) :: ufd(pcols,pver) ! U-wind for diffusion [ m/s ]
484 real(r8) :: vfd(pcols,pver) ! V-wind for diffusion [ m/s ]
486 ! Buoyancy coefficients : w'b' = ch * w'sl' + cm * w'qt'
488 real(r8) :: chu(pcols,pver+1) ! Heat buoyancy coef for dry states, defined at each interface, finally.
489 real(r8) :: chs(pcols,pver+1) ! Heat buoyancy coef for sat states, defined at each interface, finally.
490 real(r8) :: cmu(pcols,pver+1) ! Moisture buoyancy coef for dry states, defined at each interface, finally.
491 real(r8) :: cms(pcols,pver+1) ! Moisture buoyancy coef for sat states, defined at each interface, finally.
493 real(r8) :: jnk1d(pcols)
494 real(r8) :: jnk2d(pcols,pver+1)
495 real(r8) :: zero(pcols)
496 real(r8) :: zero2d(pcols,pver+1)
497 real(r8) :: es(1) ! Saturation vapor pressure
498 real(r8) :: qs(1) ! Saturation specific humidity
499 real(r8) :: gam(1) ! (L/cp)*dqs/dT
500 real(r8) :: ep2, templ, temps
502 ! ------------------------------- !
503 ! Variables for diagnostic output !
504 ! ------------------------------- !
506 real(r8) :: tkes(pcols) ! TKE at surface interface [ m2/s2 ]
507 real(r8) :: kbase_o(pcols,ncvmax) ! Original external base interface index of CL from 'exacol'
508 real(r8) :: ktop_o(pcols,ncvmax) ! Original external top interface index of CL from 'exacol'
509 real(r8) :: ncvfin_o(pcols) ! Original number of CLs from 'exacol'
510 real(r8) :: kbase_mg(pcols,ncvmax) ! 'kbase' after extending-merging from 'zisocl'
511 real(r8) :: ktop_mg(pcols,ncvmax) ! 'ktop' after extending-merging from 'zisocl'
512 real(r8) :: ncvfin_mg(pcols) ! 'ncvfin' after extending-merging from 'zisocl'
513 real(r8) :: kbase_f(pcols,ncvmax) ! Final 'kbase' after extending-merging & including SRCL
514 real(r8) :: ktop_f(pcols,ncvmax) ! Final 'ktop' after extending-merging & including SRCL
515 real(r8) :: ncvfin_f(pcols) ! Final 'ncvfin' after extending-merging & including SRCL
516 real(r8) :: wet(pcols,ncvmax) ! Entrainment rate at the CL top [ m/s ]
517 real(r8) :: web(pcols,ncvmax) ! Entrainment rate at the CL base [ m/s ]. Set to zero if CL is based at surface.
518 real(r8) :: jtbu(pcols,ncvmax) ! Buoyancy jump across the CL top [ m/s2 ]
519 real(r8) :: jbbu(pcols,ncvmax) ! Buoyancy jump across the CL base [ m/s2 ]
520 real(r8) :: evhc(pcols,ncvmax) ! Evaporative enhancement factor at the CL top
521 real(r8) :: jt2slv(pcols,ncvmax) ! Jump of slv ( across two layers ) at CL top used only for evhc [ J/kg ]
522 real(r8) :: n2ht(pcols,ncvmax) ! n2 defined at the CL top interface but using sfuh(kt) instead of sfi(kt) [ s-2 ]
523 real(r8) :: n2hb(pcols,ncvmax) ! n2 defined at the CL base interface but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
524 real(r8) :: lwp(pcols,ncvmax) ! LWP in the CL top layer [ kg/m2 ]
525 real(r8) :: opt_depth(pcols,ncvmax) ! Optical depth of the CL top layer
526 real(r8) :: radinvfrac(pcols,ncvmax) ! Fraction of radiative cooling confined in the top portion of CL top layer
527 real(r8) :: radf(pcols,ncvmax) ! Buoyancy production at the CL top due to LW radiative cooling [ m2/s3 ]
528 real(r8) :: wstar(pcols,ncvmax) ! Convective velocity in each CL [ m/s ]
529 real(r8) :: wstar3fact(pcols,ncvmax) ! Enhancement of 'wstar3' due to entrainment (inverse) [ no unit ]
530 real(r8) :: ebrk(pcols,ncvmax) ! Net mean TKE of CL including entrainment effect [ m2/s2 ]
531 real(r8) :: wbrk(pcols,ncvmax) ! Net mean normalized TKE (W) of CL, 'ebrk/b1' including entrainment effect [ m2/s2 ]
532 real(r8) :: lbrk(pcols,ncvmax) ! Energetic internal thickness of CL [m]
533 real(r8) :: ricl(pcols,ncvmax) ! CL internal mean Richardson number
534 real(r8) :: ghcl(pcols,ncvmax) ! Half of normalized buoyancy production of CL
535 real(r8) :: shcl(pcols,ncvmax) ! Galperin instability function of heat-moisture of CL
536 real(r8) :: smcl(pcols,ncvmax) ! Galperin instability function of mementum of CL
537 real(r8) :: ghi(pcols,pver+1) ! Half of normalized buoyancy production at all interfaces
538 real(r8) :: shi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces
539 real(r8) :: smi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces
540 real(r8) :: rii(pcols,pver+1) ! Interfacial Richardson number defined at all interfaces
541 real(r8) :: lengi(pcols,pver+1) ! Turbulence length scale at all interfaces [ m ]
542 real(r8) :: wcap(pcols,pver+1) ! Normalized TKE at all interfaces [ m2/s2 ]
544 ! ---------- !
545 ! Initialize !
546 ! ---------- !
548 zero(:) = 0._r8
549 zero2d(:,:) = 0._r8
551 ! ----------------------- !
552 ! Main Computation Begins !
553 ! ----------------------- !
555 ufd(:ncol,:) = u(:ncol,:)
556 vfd(:ncol,:) = v(:ncol,:)
557 tfd(:ncol,:) = t(:ncol,:)
558 qvfd(:ncol,:) = qv(:ncol,:)
559 qlfd(:ncol,:) = ql(:ncol,:)
561 do iturb = 1, nturb
563 ! Compute total stress by including 'tms'.
564 ! Here, in computing 'tms', we can use either iteratively changed 'ufd,vfd' or the
565 ! initially given 'u,v' to the PBL scheme. Note that normal stress, 'taux, tauy'
566 ! are not changed by iteration. In order to treat 'tms' in a fully implicit way,
567 ! I am using updated wind, here.
569 tautotx(:ncol) = taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver)
570 tautoty(:ncol) = tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver)
572 ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)
574 call trbintd( &
575 pcols , pver , ncol , z , ufd , vfd , tfd , pmid , &
576 tautotx , tautoty , ustar , rrho , s2 , n2 , ri , zi , &
577 pi , cldn , qtfd , qvfd , qlfd , qi , sfi , sfuh , &
578 sflh , slfd , slv , slslope , qtslope , chs , chu , cms , &
579 cmu , minpblh , qsat )
581 ! Save initial (i.e., before iterative diffusion) profile of (qt,sl) at each iteration.
582 ! Only necessary for (qt,sl) not (u,v) because (qt,sl) are newly calculated variables.
584 if( iturb .eq. 1 ) then
585 qt(:ncol,:) = qtfd(:ncol,:)
586 sl(:ncol,:) = slfd(:ncol,:)
587 endif
589 ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
590 #ifndef WRF_PORT
591 call austausch_atm( pcols, pver, ncol, ri, s2, kvf )
592 #else
593 if(use_kvf)call austausch_atm( pcols, pver, ncol, ri, s2, kvf )
594 #endif
596 ! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number
597 ! This is necessary for 'wstar-based' entrainment closure.
599 if( iturb .eq. 1 ) then
600 if( kvinit ) then
601 ! First iteration of first model timestep : Use free tropospheric value or zero.
602 if( use_kvf ) then
603 kvh(:ncol,:) = kvf(:ncol,:)
604 kvm(:ncol,:) = kvf(:ncol,:)
605 else
606 kvh(:ncol,:) = 0._r8
607 kvm(:ncol,:) = 0._r8
608 endif
609 else
610 ! First iteration on any model timestep except the first : Use value from previous timestep
611 kvh(:ncol,:) = kvh_in(:ncol,:)
612 kvm(:ncol,:) = kvm_in(:ncol,:)
613 endif
614 else
615 ! Not the first iteration : Use from previous iteration
616 kvh(:ncol,:) = kvh_out(:ncol,:)
617 kvm(:ncol,:) = kvm_out(:ncol,:)
618 endif
620 ! Calculate eddy diffusivity (kvh_out,kvm_out) and (tke,bprod,sprod) using
621 ! a given (kvh,kvm) which are used only for initializing (bprod,sprod) at
622 ! the first part of caleddy. (bprod,sprod) are fully updated at the end of
623 ! caleddy after calculating (kvh_out,kvm_out)
625 call caleddy( pcols , pver , ncol , &
626 slfd , qtfd , qlfd , slv ,ufd , &
627 vfd , pi , z , zi , &
628 qflx , shflx , slslope , qtslope , &
629 chu , chs , cmu , cms ,sfuh , &
630 sflh , n2 , s2 , ri ,rrho , &
631 pblh , ustar , &
632 kvh , kvm , kvh_out , kvm_out , &
633 tpert , qpert , qrl , kvf , tke , &
634 wstarent , bprod , sprod , minpblh , wpert , &
635 tkes , turbtype , sm_aw , &
636 kbase_o , ktop_o , ncvfin_o , &
637 kbase_mg , ktop_mg , ncvfin_mg , &
638 kbase_f , ktop_f , ncvfin_f , &
639 wet , web , jtbu , jbbu , &
640 evhc , jt2slv , n2ht , n2hb , &
641 lwp , opt_depth , radinvfrac, radf , &
642 wstar , wstar3fact, &
643 ebrk , wbrk , lbrk , ricl , ghcl , &
644 shcl , smcl , ghi , shi , smi , &
645 rii , lengi , wcap , pblhp , cldn , &
646 ipbl , kpblh , wsedl )
648 ! Calculate errorPBL to check whether PBL produced convergent solutions or not.
650 if( iturb .eq. nturb ) then
651 do i = 1, ncol
652 errorPBL(i) = 0._r8
653 do k = 1, pver
654 errorPBL(i) = errorPBL(i) + ( kvh(i,k) - kvh_out(i,k) )**2
655 end do
656 errorPBL(i) = sqrt(errorPBL(i)/pver)
657 end do
658 end if
660 ! Eddy diffusivities which will be used for the initialization of (bprod,
661 ! sprod) in 'caleddy' at the next iteration step.
663 if( iturb .gt. 1 .and. iturb .lt. nturb ) then
664 kvm_out(:ncol,:) = lambda * kvm_out(:ncol,:) + ( 1._r8 - lambda ) * kvm(:ncol,:)
665 kvh_out(:ncol,:) = lambda * kvh_out(:ncol,:) + ( 1._r8 - lambda ) * kvh(:ncol,:)
666 endif
668 ! Set nonlocal terms to zero for flux diagnostics, since not used by caleddy.
670 cgh(:ncol,:) = 0._r8
671 cgs(:ncol,:) = 0._r8
673 if( iturb .lt. nturb ) then
675 ! Each time we diffuse the original state
677 slfd(:ncol,:) = sl(:ncol,:)
678 qtfd(:ncol,:) = qt(:ncol,:)
679 ufd(:ncol,:) = u(:ncol,:)
680 vfd(:ncol,:) = v(:ncol,:)
682 ! Diffuse initial profile of each time step using a given (kvh_out,kvm_out)
683 ! In the below 'compute_vdiff', (slfd,qtfd,ufd,vfd) are 'inout' variables.
685 call compute_vdiff( lchnk , &
686 pcols , pver , 1 , ncol , pmid , &
687 pi , rpdel , t , ztodt , taux , &
688 tauy , shflx , qflx , ntop_turb , nbot_turb , &
689 kvh_out , kvm_out , kvh_out , cgs , cgh , &
690 zi , ksrftms , zero , fieldlist_wet, &
691 ufd , vfd , qtfd , slfd , &
692 jnk1d , jnk1d , jnk2d , jnk1d , errstring , &
693 tauresx , tauresy , 0 , .false. )
695 ! Retrieve (tfd,qvfd,qlfd) from (slfd,qtfd) in order to
696 ! use 'trbintd' at the next iteration.
698 do k = 1, pver
699 do i = 1, ncol
700 ! ----------------------------------------------------- !
701 ! Compute the condensate 'qlfd' in the updated profiles !
702 ! ----------------------------------------------------- !
703 ! Option.1 : Assume grid-mean condensate is homogeneously diffused by the moist turbulence scheme.
704 ! This should bs used if 'pseudodiff = .false.' in vertical_diffusion.F90.
705 ! Modification : Need to be check whether below is correct in the presence of ice, qi.
706 ! I should understand why the variation of ice, qi is neglected during diffusion.
707 templ = ( slfd(i,k) - g*z(i,k) ) / cpair
708 status = qsat( templ, pmid(i,k), es(1), qs(1), gam(1), 1 )
709 ep2 = .622_r8
710 temps = templ + ( qtfd(i,k) - qs(1) ) / ( cpair / latvap + latvap * qs(1) / ( rair * templ**2 ) )
711 status = qsat( temps, pmid(i,k), es(1), qs(1), gam(1), 1 )
712 qlfd(i,k) = max( qtfd(i,k) - qi(i,k) - qs(1) ,0._r8 )
713 ! Option.2 : Assume condensate is not diffused by the moist turbulence scheme.
714 ! This should bs used if 'pseudodiff = .true.' in vertical_diffusion.F90.
715 ! qlfd(i,k) = ql(i,k)
716 ! ----------------------------- !
717 ! Compute the other 'qvfd, tfd' !
718 ! ----------------------------- !
719 qvfd(i,k) = max( 0._r8, qtfd(i,k) - qi(i,k) - qlfd(i,k) )
720 tfd(i,k) = ( slfd(i,k) + latvap * qlfd(i,k) + latsub * qi(i,k) - g*z(i,k)) / cpair
721 end do
722 end do
723 endif
725 ! Debug
726 ! icol = phys_debug_col(lchnk)
727 ! if( icol > 0 .and. get_nstep() .ge. 1 ) then
728 ! write(iulog,*) ' '
729 ! write(iulog,*) 'eddy_diff debug at the end of iteration'
730 ! write(iulog,*) 't, qv, ql, cld, u, v'
731 ! do k = pver-3, pver
732 ! write (iulog,*) k, tfd(icol,k), qvfd(icol,k), qlfd(icol,k), cldn(icol,k), ufd(icol,k), vfd(icol,k)
733 ! end do
734 ! endif
735 ! Debug
737 end do ! End of 'iturb' iteration
739 kvq(:ncol,:) = kvh_out(:ncol,:)
741 ! Compute 'wstar' within the PBL for use in the future convection scheme.
743 do i = 1, ncol
744 if( ipbl(i) .eq. 1._r8 ) then
745 wstarPBL(i) = max( 0._r8, wstar(i,1) )
746 else
747 wstarPBL(i) = 0._r8
748 endif
749 end do
751 ! --------------------------------------------------------------- !
752 ! Writing for detailed diagnostic analysis of UW moist PBL scheme !
753 ! --------------------------------------------------------------- !
755 call outfld( 'UW_errorPBL', errorPBL, pcols, lchnk )
757 call outfld( 'UW_n2', n2, pcols, lchnk )
758 call outfld( 'UW_s2', s2, pcols, lchnk )
759 call outfld( 'UW_ri', ri, pcols, lchnk )
761 call outfld( 'UW_sfuh', sfuh, pcols, lchnk )
762 call outfld( 'UW_sflh', sflh, pcols, lchnk )
763 call outfld( 'UW_sfi', sfi, pcols, lchnk )
765 call outfld( 'UW_cldn', cldn, pcols, lchnk )
766 call outfld( 'UW_qrl', qrl, pcols, lchnk )
767 call outfld( 'UW_ql', qlfd, pcols, lchnk )
769 call outfld( 'UW_chu', chu, pcols, lchnk )
770 call outfld( 'UW_chs', chs, pcols, lchnk )
771 call outfld( 'UW_cmu', cmu, pcols, lchnk )
772 call outfld( 'UW_cms', cms, pcols, lchnk )
774 call outfld( 'UW_tke', tke, pcols, lchnk )
775 call outfld( 'UW_wcap', wcap, pcols, lchnk )
776 call outfld( 'UW_bprod', bprod, pcols, lchnk )
777 call outfld( 'UW_sprod', sprod, pcols, lchnk )
779 call outfld( 'UW_kvh', kvh_out, pcols, lchnk )
780 call outfld( 'UW_kvm', kvm_out, pcols, lchnk )
782 call outfld( 'UW_pblh', pblh, pcols, lchnk )
783 call outfld( 'UW_pblhp', pblhp, pcols, lchnk )
784 call outfld( 'UW_tpert', tpert, pcols, lchnk )
785 call outfld( 'UW_qpert', qpert, pcols, lchnk )
786 call outfld( 'UW_wpert', wpert, pcols, lchnk )
788 call outfld( 'UW_ustar', ustar, pcols, lchnk )
789 call outfld( 'UW_tkes', tkes, pcols, lchnk )
790 call outfld( 'UW_minpblh', minpblh, pcols, lchnk )
792 call outfld( 'UW_turbtype', turbtype, pcols, lchnk )
794 call outfld( 'UW_kbase_o', kbase_o, pcols, lchnk )
795 call outfld( 'UW_ktop_o', ktop_o, pcols, lchnk )
796 call outfld( 'UW_ncvfin_o', ncvfin_o, pcols, lchnk )
798 call outfld( 'UW_kbase_mg', kbase_mg, pcols, lchnk )
799 call outfld( 'UW_ktop_mg', ktop_mg, pcols, lchnk )
800 call outfld( 'UW_ncvfin_mg', ncvfin_mg, pcols, lchnk )
802 call outfld( 'UW_kbase_f', kbase_f, pcols, lchnk )
803 call outfld( 'UW_ktop_f', ktop_f, pcols, lchnk )
804 call outfld( 'UW_ncvfin_f', ncvfin_f, pcols, lchnk )
806 call outfld( 'UW_wet', wet, pcols, lchnk )
807 call outfld( 'UW_web', web, pcols, lchnk )
808 call outfld( 'UW_jtbu', jtbu, pcols, lchnk )
809 call outfld( 'UW_jbbu', jbbu, pcols, lchnk )
810 call outfld( 'UW_evhc', evhc, pcols, lchnk )
811 call outfld( 'UW_jt2slv', jt2slv, pcols, lchnk )
812 call outfld( 'UW_n2ht', n2ht, pcols, lchnk )
813 call outfld( 'UW_n2hb', n2hb, pcols, lchnk )
814 call outfld( 'UW_lwp', lwp, pcols, lchnk )
815 call outfld( 'UW_optdepth', opt_depth, pcols, lchnk )
816 call outfld( 'UW_radfrac', radinvfrac, pcols, lchnk )
817 call outfld( 'UW_radf', radf, pcols, lchnk )
818 call outfld( 'UW_wstar', wstar, pcols, lchnk )
819 call outfld( 'UW_wstar3fact', wstar3fact, pcols, lchnk )
820 call outfld( 'UW_ebrk', ebrk, pcols, lchnk )
821 call outfld( 'UW_wbrk', wbrk, pcols, lchnk )
822 call outfld( 'UW_lbrk', lbrk, pcols, lchnk )
823 call outfld( 'UW_ricl', ricl, pcols, lchnk )
824 call outfld( 'UW_ghcl', ghcl, pcols, lchnk )
825 call outfld( 'UW_shcl', shcl, pcols, lchnk )
826 call outfld( 'UW_smcl', smcl, pcols, lchnk )
828 call outfld( 'UW_gh', ghi, pcols, lchnk )
829 call outfld( 'UW_sh', shi, pcols, lchnk )
830 call outfld( 'UW_sm', smi, pcols, lchnk )
831 call outfld( 'UW_ria', rii, pcols, lchnk )
832 call outfld( 'UW_leng', lengi, pcols, lchnk )
834 return
836 end subroutine compute_eddy_diff
838 !=============================================================================== !
839 ! !
840 !=============================================================================== !
842 subroutine sfdiag( pcols , pver , ncol , qt , ql , sl , &
843 pi , pm , zi , cld , sfi , sfuh , &
844 sflh , slslope , qtslope , qsat )
845 !----------------------------------------------------------------------- !
846 ! !
847 ! Purpose: Interface for calculating saturation fractions at upper and !
848 ! lower-half layers, & interfaces for use by turbulence scheme !
849 ! !
850 ! Method : Various but 'l' should be chosen for consistency. !
851 ! !
852 ! Author : B. Stevens and C. Bretherton (August 2000) !
853 ! Sungsu Park. August 2006. !
854 ! May. 2008. !
855 ! !
856 ! S.Park : The computed saturation fractions are repeatedly !
857 ! used to compute buoyancy coefficients in'trbintd' & 'caleddy'.!
858 !----------------------------------------------------------------------- !
860 implicit none
862 ! --------------- !
863 ! Input arguments !
864 ! --------------- !
866 integer, external :: qsat
868 integer, intent(in) :: pcols ! Number of atmospheric columns
869 integer, intent(in) :: pver ! Number of atmospheric layers
870 integer, intent(in) :: ncol ! Number of atmospheric columns
872 real(r8), intent(in) :: sl(pcols,pver) ! Liquid water static energy [ J/kg ]
873 real(r8), intent(in) :: qt(pcols,pver) ! Total water specific humidity [ kg/kg ]
874 real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ]
875 real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressures [ Pa ]
876 real(r8), intent(in) :: pm(pcols,pver) ! Layer mid-point pressures [ Pa ]
877 real(r8), intent(in) :: zi(pcols,pver+1) ! Interface heights [ m ]
878 real(r8), intent(in) :: cld(pcols,pver) ! Stratiform cloud fraction [ fraction ]
879 real(r8), intent(in) :: slslope(pcols,pver) ! Slope of 'sl' in each layer
880 real(r8), intent(in) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer
882 ! ---------------- !
883 ! Output arguments !
884 ! ---------------- !
886 real(r8), intent(out) :: sfi(pcols,pver+1) ! Interfacial layer saturation fraction [ fraction ]
887 real(r8), intent(out) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ]
888 real(r8), intent(out) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ]
890 ! --------------- !
891 ! Local Variables !
892 ! --------------- !
894 integer :: i ! Longitude index
895 integer :: k ! Vertical index
896 integer :: km1 ! k-1
897 integer :: status ! Status returned by function calls
898 real(r8) :: sltop, slbot ! sl at top/bot of grid layer
899 real(r8) :: qttop, qtbot ! qt at top/bot of grid layer
900 real(r8) :: tltop(1), tlbot(1) ! Liquid water temperature at top/bot of grid layer
901 real(r8) :: qxtop, qxbot ! Sat excess at top/bot of grid layer
902 real(r8) :: qxm ! Sat excess at midpoint
903 real(r8) :: es(1) ! Saturation vapor pressure
904 real(r8) :: qs(1) ! Saturation spec. humidity
905 real(r8) :: gam(1) ! (L/cp)*dqs/dT
906 real(r8) :: cldeff(pcols,pver) ! Effective Cloud Fraction [ fraction ]
908 ! ----------------------- !
909 ! Main Computation Begins !
910 ! ----------------------- !
912 sfi(1:ncol,:) = 0._r8
913 sfuh(1:ncol,:) = 0._r8
914 sflh(1:ncol,:) = 0._r8
915 cldeff(1:ncol,:) = 0._r8
917 select case (sftype)
918 case ('d')
919 ! ----------------------------------------------------------------------- !
920 ! Simply use the given stratus fraction ('horizontal' cloud partitioning) !
921 ! ----------------------------------------------------------------------- !
922 do k = ntop_turb + 1, nbot_turb
923 km1 = k - 1
924 do i = 1, ncol
925 sfuh(i,k) = cld(i,k)
926 sflh(i,k) = cld(i,k)
927 sfi(i,k) = 0.5_r8 * ( sflh(i,km1) + min( sflh(i,km1), sfuh(i,k) ) )
928 end do
929 end do
930 do i = 1, ncol
931 sfi(i,pver+1) = sflh(i,pver)
932 end do
933 case ('l')
934 ! ------------------------------------------ !
935 ! Use modified stratus fraction partitioning !
936 ! ------------------------------------------ !
937 do k = ntop_turb + 1, nbot_turb
938 km1 = k - 1
939 do i = 1, ncol
940 cldeff(i,k) = cld(i,k)
941 sfuh(i,k) = cld(i,k)
942 sflh(i,k) = cld(i,k)
943 if( ql(i,k) .lt. qmin ) then
944 sfuh(i,k) = 0._r8
945 sflh(i,k) = 0._r8
946 end if
947 ! Modification : The contribution of ice should be carefully considered.
948 if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then
949 cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
950 sfuh(i,k) = cldeff(i,k)
951 sflh(i,k) = cldeff(i,k)
952 elseif( choice_evhc .eq. 'maxi' .or. choice_radf .eq. 'maxi' ) then
953 cldeff(i,k) = cld(i,k)
954 sfuh(i,k) = cldeff(i,k)
955 sflh(i,k) = cldeff(i,k)
956 endif
957 ! At the stratus top, take the minimum interfacial saturation fraction
958 sfi(i,k) = 0.5_r8 * ( sflh(i,km1) + min( sfuh(i,k), sflh(i,km1) ) )
959 ! Modification : Currently sfi at the top and surface interfaces are set to be zero.
960 ! Also, sfuh and sflh in the top model layer is set to be zero.
961 ! However, I may need to set
962 ! do i = 1, ncol
963 ! sfi(i,pver+1) = sflh(i,pver)
964 ! end do
965 ! for treating surface-based fog.
966 ! OK. I added below block similar to the other cases.
967 end do
968 end do
969 do i = 1, ncol
970 sfi(i,pver+1) = sflh(i,pver)
971 end do
972 case ('u')
973 ! ------------------------------------------------------------------------- !
974 ! Use unsaturated buoyancy - since sfi, sfuh, sflh have already been zeroed !
975 ! nothing more need be done for this case. !
976 ! ------------------------------------------------------------------------- !
977 case ('z')
978 ! ------------------------------------------------------------------------- !
979 ! Calculate saturation fraction based on whether the air just above or just !
980 ! below the interface is saturated, i.e. with vertical cloud partitioning. !
981 ! The saturation fraction of the interfacial layer between mid-points k and !
982 ! k+1 is computed by averaging the saturation fraction of the half-layers !
983 ! above and below the interface, with a special provision for cloud tops !
984 ! (more cloud in the half-layer below than in the half-layer above).In each !
985 ! half-layer, vertical partitioning of cloud based on the slopes diagnosed !
986 ! above is used. Loop down through the layers, computing the saturation !
987 ! fraction in each half-layer (sfuh for upper half, sflh for lower half). !
988 ! Once sfuh(i,k) is computed, use with sflh(i,k-1) to determine saturation !
989 ! fraction sfi(i,k) for interfacial layer k-0.5. !
990 ! This is 'not' chosen for full consistent treatment of stratus fraction in !
991 ! all physics schemes. !
992 ! ------------------------------------------------------------------------- !
993 do k = ntop_turb + 1, nbot_turb
994 km1 = k - 1
995 do i = 1, ncol
996 ! Compute saturation excess at the mid-point of layer k
997 sltop = sl(i,k) + slslope(i,k) * ( pi(i,k) - pm(i,k) )
998 qttop = qt(i,k) + qtslope(i,k) * ( pi(i,k) - pm(i,k) )
999 tltop(1) = ( sltop - g * zi(i,k) ) / cpair
1000 status = qsat( tltop(1), pi(i,k), es(1), qs(1), gam(1), 1 )
1001 qxtop = qttop - qs(1)
1002 slbot = sl(i,k) + slslope(i,k) * ( pi(i,k+1) - pm(i,k) )
1003 qtbot = qt(i,k) + qtslope(i,k) * ( pi(i,k+1) - pm(i,k) )
1004 tlbot(1) = ( slbot - g * zi(i,k+1) ) / cpair
1005 status = qsat( tlbot(1), pi(i,k+1), es(1), qs(1), gam(1), 1 )
1006 qxbot = qtbot - qs(1)
1007 qxm = qxtop + ( qxbot - qxtop ) * ( pm(i,k) - pi(i,k) ) / ( pi(i,k+1) - pi(i,k) )
1008 ! Find the saturation fraction sfuh(i,k) of the upper half of layer k.
1009 if( ( qxtop .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
1010 sfuh(i,k) = 0._r8
1011 else if( ( qxtop .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
1012 sfuh(i,k) = 1._r8
1013 else ! Either qxm < 0 and qxtop > 0 or vice versa
1014 sfuh(i,k) = max( qxtop, qxm ) / abs( qxtop - qxm )
1015 end if
1016 ! Combine with sflh(i) (still for layer k-1) to get interfac layer saturation fraction
1017 sfi(i,k) = 0.5_r8 * ( sflh(i,k-1) + min( sflh(i,k-1), sfuh(i,k) ) )
1018 ! Update sflh to be for the lower half of layer k.
1019 if( ( qxbot .lt. 0._r8 ) .and. ( qxm .lt. 0._r8 ) ) then
1020 sflh(i,k) = 0._r8
1021 else if( ( qxbot .gt. 0._r8 ) .and. ( qxm .gt. 0._r8 ) ) then
1022 sflh(i,k) = 1._r8
1023 else ! Either qxm < 0 and qxbot > 0 or vice versa
1024 sflh(i,k) = max( qxbot, qxm ) / abs( qxbot - qxm )
1025 end if
1026 end do ! i
1027 end do ! k
1028 do i = 1, ncol
1029 sfi(i,pver+1) = sflh(i,pver) ! Saturation fraction in the lowest half-layer.
1030 end do
1031 end select
1033 return
1034 end subroutine sfdiag
1036 !=============================================================================== !
1037 ! !
1038 !=============================================================================== !
1040 subroutine trbintd( pcols , pver , ncol , &
1041 z , u , v , &
1042 t , pmid , taux , &
1043 tauy , ustar , rrho , &
1044 s2 , n2 , ri , &
1045 zi , pi , cld , &
1046 qt , qv , ql , qi , sfi , sfuh , &
1047 sflh , sl , slv , slslope , qtslope , &
1048 chs , chu , cms , cmu , minpblh , qsat )
1049 !----------------------------------------------------------------------- !
1050 ! Purpose: Calculate buoyancy coefficients at all interfaces including !
1051 ! surface. Also, computes the profiles of ( sl,qt,n2,s2,ri ). !
1052 ! Note that (n2,s2,ri) are defined at each interfaces except !
1053 ! surface. !
1054 ! !
1055 ! Author: B. Stevens ( Extracted from pbldiff, August, 2000 ) !
1056 ! Sungsu Park ( August 2006, May. 2008 ) !
1057 !----------------------------------------------------------------------- !
1059 implicit none
1061 ! --------------- !
1062 ! Input arguments !
1063 ! --------------- !
1065 integer, intent(in) :: pcols ! Number of atmospheric columns
1066 integer, intent(in) :: pver ! Number of atmospheric layers
1067 integer, intent(in) :: ncol ! Number of atmospheric columns
1068 real(r8), intent(in) :: z(pcols,pver) ! Layer mid-point height above surface [ m ]
1069 real(r8), intent(in) :: u(pcols,pver) ! Layer mid-point u [ m/s ]
1070 real(r8), intent(in) :: v(pcols,pver) ! Layer mid-point v [ m/s ]
1071 real(r8), intent(in) :: t(pcols,pver) ! Layer mid-point temperature [ K ]
1072 real(r8), intent(in) :: pmid(pcols,pver) ! Layer mid-point pressure [ Pa ]
1073 real(r8), intent(in) :: taux(pcols) ! Surface u stress [ N/m2 ]
1074 real(r8), intent(in) :: tauy(pcols) ! Surface v stress [ N/m2 ]
1075 real(r8), intent(in) :: zi(pcols,pver+1) ! Interface height [ m ]
1076 real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressure [ Pa ]
1077 real(r8), intent(in) :: cld(pcols,pver) ! Stratus fraction
1078 real(r8), intent(in) :: qv(pcols,pver) ! Water vapor specific humidity [ kg/kg ]
1079 real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ]
1080 real(r8), intent(in) :: qi(pcols,pver) ! Ice water specific humidity [ kg/kg ]
1081 integer, external :: qsat
1083 ! ---------------- !
1084 ! Output arguments !
1085 ! ---------------- !
1087 real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ]
1088 real(r8), intent(out) :: s2(pcols,pver) ! Interfacial ( except surface ) shear squared [ s-2 ]
1089 real(r8), intent(out) :: n2(pcols,pver) ! Interfacial ( except surface ) buoyancy frequency [ s-2 ]
1090 real(r8), intent(out) :: ri(pcols,pver) ! Interfacial ( except surface ) Richardson number, 'n2/s2'
1092 real(r8), intent(out) :: qt(pcols,pver) ! Total specific humidity [ kg/kg ]
1093 real(r8), intent(out) :: sfi(pcols,pver+1) ! Interfacial layer saturation fraction [ fraction ]
1094 real(r8), intent(out) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ]
1095 real(r8), intent(out) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ]
1096 real(r8), intent(out) :: sl(pcols,pver) ! Liquid water static energy [ J/kg ]
1097 real(r8), intent(out) :: slv(pcols,pver) ! Liquid water virtual static energy [ J/kg ]
1099 real(r8), intent(out) :: chu(pcols,pver+1) ! Heat buoyancy coef for dry states at all interfaces, finally. [ unit ? ]
1100 real(r8), intent(out) :: chs(pcols,pver+1) ! heat buoyancy coef for sat states at all interfaces, finally. [ unit ? ]
1101 real(r8), intent(out) :: cmu(pcols,pver+1) ! Moisture buoyancy coef for dry states at all interfaces, finally. [ unit ? ]
1102 real(r8), intent(out) :: cms(pcols,pver+1) ! Moisture buoyancy coef for sat states at all interfaces, finally. [ unit ? ]
1103 real(r8), intent(out) :: slslope(pcols,pver) ! Slope of 'sl' in each layer
1104 real(r8), intent(out) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer
1105 real(r8), intent(out) :: rrho(pcols) ! 1./bottom level density [ m3/kg ]
1106 real(r8), intent(out) :: minpblh(pcols) ! Minimum PBL height based on surface stress [ m ]
1108 ! --------------- !
1109 ! Local Variables !
1110 ! --------------- !
1112 integer :: i ! Longitude index
1113 integer :: k, km1 ! Level index
1114 integer :: status ! Status returned by function calls
1116 real(r8) :: qs(pcols,pver) ! Saturation specific humidity
1117 real(r8) :: es(pcols,pver) ! Saturation vapor pressure
1118 real(r8) :: gam(pcols,pver) ! (l/cp)*(d(qs)/dT)
1119 real(r8) :: rdz ! 1 / (delta z) between midpoints
1120 real(r8) :: dsldz ! 'delta sl / delta z' at interface
1121 real(r8) :: dqtdz ! 'delta qt / delta z' at interface
1122 real(r8) :: ch ! 'sfi' weighted ch at the interface
1123 real(r8) :: cm ! 'sfi' weighted cm at the interface
1124 real(r8) :: bfact ! Buoyancy factor in n2 calculations
1125 real(r8) :: product ! Intermediate vars used to find slopes
1126 real(r8) :: dsldp_a, dqtdp_a ! Slopes across interface above
1127 real(r8) :: dsldp_b(pcols), dqtdp_b(pcols) ! Slopes across interface below
1129 ! ----------------------- !
1130 ! Main Computation Begins !
1131 ! ----------------------- !
1133 ! Compute ustar, and kinematic surface fluxes from surface energy fluxes
1135 do i = 1, ncol
1136 rrho(i) = rair * t(i,pver) / pmid(i,pver)
1137 ustar(i) = max( sqrt( sqrt( taux(i)**2 + tauy(i)**2 ) * rrho(i) ), ustar_min )
1138 minpblh(i) = 100.0_r8 * ustar(i) ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.
1139 end do
1141 ! Calculate conservative scalars (qt,sl,slv) and buoyancy coefficients at the layer mid-points.
1142 ! Note that 'ntop_turb = 1', 'nbot_turb = pver'
1144 do k = ntop_turb, nbot_turb
1145 status = qsat( t(1,k), pmid(1,k), es(1,k), qs(1,k), gam(1,k), ncol )
1146 do i = 1, ncol
1147 qt(i,k) = qv(i,k) + ql(i,k) + qi(i,k)
1148 sl(i,k) = cpair * t(i,k) + g * z(i,k) - latvap * ql(i,k) - latsub * qi(i,k)
1149 slv(i,k) = sl(i,k) * ( 1._r8 + zvir * qt(i,k) )
1150 ! Thermodynamic coefficients for buoyancy flux - in this loop these are
1151 ! calculated at mid-points; later, they will be averaged to interfaces,
1152 ! where they will ultimately be used. At the surface, the coefficients
1153 ! are taken from the lowest mid point.
1154 bfact = g / ( t(i,k) * ( 1._r8 + zvir * qv(i,k) - ql(i,k) - qi(i,k) ) )
1155 chu(i,k) = ( 1._r8 + zvir * qt(i,k) ) * bfact / cpair
1156 chs(i,k) = ( ( 1._r8 + ( 1._r8 + zvir ) * gam(i,k) * cpair * t(i,k) / latvap ) / ( 1._r8 + gam(i,k) ) ) * bfact / cpair
1157 cmu(i,k) = zvir * bfact * t(i,k)
1158 cms(i,k) = latvap * chs(i,k) - bfact * t(i,k)
1159 end do
1160 end do
1162 do i = 1, ncol
1163 chu(i,pver+1) = chu(i,pver)
1164 chs(i,pver+1) = chs(i,pver)
1165 cmu(i,pver+1) = cmu(i,pver)
1166 cms(i,pver+1) = cms(i,pver)
1167 end do
1169 ! Compute slopes of conserved variables sl, qt within each layer k.
1170 ! 'a' indicates the 'above' gradient from layer k-1 to layer k and
1171 ! 'b' indicates the 'below' gradient from layer k to layer k+1.
1172 ! We take a smaller (in absolute value) of these gradients as the
1173 ! slope within layer k. If they have opposite signs, gradient in
1174 ! layer k is taken to be zero. I should re-consider whether this
1175 ! profile reconstruction is the best or not.
1176 ! This is similar to the profile reconstruction used in the UWShCu.
1178 do i = 1, ncol
1179 ! Slopes at endpoints determined by extrapolation
1180 slslope(i,pver) = ( sl(i,pver) - sl(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
1181 qtslope(i,pver) = ( qt(i,pver) - qt(i,pver-1) ) / ( pmid(i,pver) - pmid(i,pver-1) )
1182 slslope(i,1) = ( sl(i,2) - sl(i,1) ) / ( pmid(i,2) - pmid(i,1) )
1183 qtslope(i,1) = ( qt(i,2) - qt(i,1) ) / ( pmid(i,2) - pmid(i,1) )
1184 dsldp_b(i) = slslope(i,1)
1185 dqtdp_b(i) = qtslope(i,1)
1186 end do
1188 do k = 2, pver - 1
1189 do i = 1, ncol
1190 dsldp_a = dsldp_b(i)
1191 dqtdp_a = dqtdp_b(i)
1192 dsldp_b(i) = ( sl(i,k+1) - sl(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
1193 dqtdp_b(i) = ( qt(i,k+1) - qt(i,k) ) / ( pmid(i,k+1) - pmid(i,k) )
1194 product = dsldp_a * dsldp_b(i)
1195 if( product .le. 0._r8 ) then
1196 slslope(i,k) = 0._r8
1197 else if( product .gt. 0._r8 .and. dsldp_a .lt. 0._r8 ) then
1198 slslope(i,k) = max( dsldp_a, dsldp_b(i) )
1199 else if( product .gt. 0._r8 .and. dsldp_a .gt. 0._r8 ) then
1200 slslope(i,k) = min( dsldp_a, dsldp_b(i) )
1201 end if
1202 product = dqtdp_a*dqtdp_b(i)
1203 if( product .le. 0._r8 ) then
1204 qtslope(i,k) = 0._r8
1205 else if( product .gt. 0._r8 .and. dqtdp_a .lt. 0._r8 ) then
1206 qtslope(i,k) = max( dqtdp_a, dqtdp_b(i) )
1207 else if( product .gt. 0._r8 .and. dqtdp_a .gt. 0._r8 ) then
1208 qtslope(i,k) = min( dqtdp_a, dqtdp_b(i) )
1209 end if
1210 end do ! i
1211 end do ! k
1213 ! Compute saturation fraction at the interfacial layers for use in buoyancy
1214 ! flux computation.
1216 call sfdiag( pcols , pver , ncol , qt , ql , sl , &
1217 pi , pmid , zi , cld , sfi , sfuh , &
1218 sflh , slslope , qtslope , qsat )
1220 ! Calculate buoyancy coefficients at all interfaces (1:pver+1) and (n2,s2,ri)
1221 ! at all interfaces except surface. Note 'nbot_turb = pver', 'ntop_turb = 1'.
1222 ! With the previous definition of buoyancy coefficients at the surface, the
1223 ! resulting buoyancy coefficients at the top and surface interfaces becomes
1224 ! identical to the buoyancy coefficients at the top and bottom layers. Note
1225 ! that even though the dimension of (s2,n2,ri) is 'pver', they are defined
1226 ! at interfaces ( not at the layer mid-points ) except the surface.
1228 do k = nbot_turb, ntop_turb + 1, -1
1229 km1 = k - 1
1230 do i = 1, ncol
1231 rdz = 1._r8 / ( z(i,km1) - z(i,k) )
1232 dsldz = ( sl(i,km1) - sl(i,k) ) * rdz
1233 dqtdz = ( qt(i,km1) - qt(i,k) ) * rdz
1234 chu(i,k) = ( chu(i,km1) + chu(i,k) ) * 0.5_r8
1235 chs(i,k) = ( chs(i,km1) + chs(i,k) ) * 0.5_r8
1236 cmu(i,k) = ( cmu(i,km1) + cmu(i,k) ) * 0.5_r8
1237 cms(i,k) = ( cms(i,km1) + cms(i,k) ) * 0.5_r8
1238 ch = chu(i,k) * ( 1._r8 - sfi(i,k) ) + chs(i,k) * sfi(i,k)
1239 cm = cmu(i,k) * ( 1._r8 - sfi(i,k) ) + cms(i,k) * sfi(i,k)
1240 n2(i,k) = ch * dsldz + cm * dqtdz
1241 s2(i,k) = ( ( u(i,km1) - u(i,k) )**2 + ( v(i,km1) - v(i,k) )**2) * rdz**2
1242 s2(i,k) = max( ntzero, s2(i,k) )
1243 ri(i,k) = n2(i,k) / s2(i,k)
1244 end do
1245 end do
1246 do i = 1, ncol
1247 n2(i,1) = n2(i,2)
1248 s2(i,1) = s2(i,2)
1249 ri(i,1) = ri(i,2)
1250 end do
1252 return
1254 end subroutine trbintd
1256 !=============================================================================== !
1257 ! !
1258 !=============================================================================== !
1260 subroutine austausch_atm( pcols, pver, ncol, ri, s2, kvf )
1262 !---------------------------------------------------------------------- !
1263 ! !
1264 ! Purpose: Computes exchange coefficients for free turbulent flows. !
1265 ! This is not used in the UW moist turbulence scheme. !
1266 ! !
1267 ! Method: !
1268 ! !
1269 ! The free atmosphere diffusivities are based on standard mixing length !
1270 ! forms for the neutral diffusivity multiplied by functns of Richardson !
1271 ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for !
1272 ! momentum, potential temperature, and constitutents. !
1273 ! !
1274 ! The stable Richardson num function (Ri>0) is taken from Holtslag and !
1275 ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri)) !
1276 ! The unstable Richardson number function (Ri<0) is taken from CCM1. !
1277 ! f = sqrt(1 - 18*Ri) !
1278 ! !
1279 ! Author: B. Stevens (rewrite, August 2000) !
1280 ! !
1281 !---------------------------------------------------------------------- !
1282 implicit none
1284 ! --------------- !
1285 ! Input arguments !
1286 ! --------------- !
1288 integer, intent(in) :: pcols ! Number of atmospheric columns
1289 integer, intent(in) :: pver ! Number of atmospheric layers
1290 integer, intent(in) :: ncol ! Number of atmospheric columns
1292 real(r8), intent(in) :: s2(pcols,pver) ! Shear squared
1293 real(r8), intent(in) :: ri(pcols,pver) ! Richardson no
1295 ! ---------------- !
1296 ! Output arguments !
1297 ! ---------------- !
1299 real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers
1301 ! --------------- !
1302 ! Local Variables !
1303 ! --------------- !
1305 real(r8) :: fofri ! f(ri)
1306 real(r8) :: kvn ! Neutral Kv
1308 integer :: i ! Longitude index
1309 integer :: k ! Vertical index
1311 ! ----------------------- !
1312 ! Main Computation Begins !
1313 ! ----------------------- !
1315 kvf(:ncol,:) = 0.0_r8
1316 kvf(:ncol,pver+1) = 0.0_r8
1317 kvf(:ncol,1:ntop_turb) = 0.0_r8
1319 ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
1321 do k = ntop_turb + 1, nbot_turb
1322 do i = 1, ncol
1323 if( ri(i,k) < 0.0_r8 ) then
1324 fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
1325 else
1326 fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) )
1327 end if
1328 kvn = ml2(k) * sqrt(s2(i,k))
1329 kvf(i,k) = max( zkmin, kvn * fofri )
1330 end do
1331 end do
1333 return
1335 end subroutine austausch_atm
1337 ! ---------------------------------------------------------------------------- !
1338 ! !
1339 ! The University of Washington Moist Turbulence Scheme !
1340 ! !
1341 ! Authors : Chris Bretherton at the University of Washington, Seattle, WA !
1342 ! Sungsu Park at the CGD/NCAR, Boulder, CO !
1343 ! !
1344 ! ---------------------------------------------------------------------------- !
1346 subroutine caleddy( pcols , pver , ncol , &
1347 sl , qt , ql , slv , u , &
1348 v , pi , z , zi , &
1349 qflx , shflx , slslope , qtslope , &
1350 chu , chs , cmu , cms , sfuh , &
1351 sflh , n2 , s2 , ri , rrho , &
1352 pblh , ustar , &
1353 kvh_in , kvm_in , kvh , kvm , &
1354 tpert , qpert , qrlin , kvf , tke , &
1355 wstarent , bprod , sprod , minpblh , wpert , &
1356 tkes , turbtype_f , sm_aw , &
1357 kbase_o , ktop_o , ncvfin_o , &
1358 kbase_mg , ktop_mg , ncvfin_mg , &
1359 kbase_f , ktop_f , ncvfin_f , &
1360 wet_CL , web_CL , jtbu_CL , jbbu_CL , &
1361 evhc_CL , jt2slv_CL , n2ht_CL , n2hb_CL , lwp_CL , &
1362 opt_depth_CL , radinvfrac_CL, radf_CL , wstar_CL , wstar3fact_CL, &
1363 ebrk , wbrk , lbrk , ricl , ghcl , &
1364 shcl , smcl , &
1365 gh_a , sh_a , sm_a , ri_a , leng , &
1366 wcap , pblhp , cld , ipbl , kpblh , &
1367 wsedl )
1369 !--------------------------------------------------------------------------------- !
1370 ! !
1371 ! Purpose : This is a driver routine to compute eddy diffusion coefficients !
1372 ! for heat (sl), momentum (u, v), moisture (qt), and other trace !
1373 ! constituents. This scheme uses first order closure for stable !
1374 ! turbulent layers (STL). For convective layers (CL), entrainment !
1375 ! closure is used at the CL external interfaces, which is coupled !
1376 ! to the diagnosis of a CL regime mean TKE from the instantaneous !
1377 ! thermodynamic and velocity profiles. The CLs are diagnosed by !
1378 ! extending original CL layers of moist static instability into !
1379 ! adjacent weakly stably stratified interfaces, stopping if the !
1380 ! stability is too strong. This allows a realistic depiction of !
1381 ! dry convective boundary layers with a downgradient approach. !
1382 ! !
1383 ! NOTE: This routine currently assumes ntop_turb = 1, nbot_turb = pver !
1384 ! ( turbulent diffusivities computed at all interior interfaces ) !
1385 ! and will require modification to handle a different ntop_turb. !
1386 ! !
1387 ! Authors: Sungsu Park and Chris Bretherton. 08/2006, 05/2008. !
1388 ! !
1389 ! For details, see !
1390 ! !
1391 ! 1. 'A new moist turbulence parametrization in the Community Atmosphere Model' !
1392 ! by Christopher S. Bretherton & Sungsu Park. J. Climate. 22. 3422-3448. 2009. !
1393 ! !
1394 ! 2. 'The University of Washington shallow convection and moist turbulence schemes !
1395 ! and their impact on climate simulations with the Community Atmosphere Model' !
1396 ! by Sungsu Park & Christopher S. Bretherton. J. Climate. 22. 3449-3469. 2009. !
1397 ! !
1398 ! For questions on the scheme and code, send an email to !
1399 ! sungsup@ucar.edu or breth@washington.edu !
1400 ! !
1401 !--------------------------------------------------------------------------------- !
1403 ! ---------------- !
1404 ! Inputs variables !
1405 ! ---------------- !
1407 implicit none
1409 integer, intent(in) :: pcols ! Number of atmospheric columns
1410 integer, intent(in) :: pver ! Number of atmospheric layers
1411 integer, intent(in) :: ncol ! Number of atmospheric columns
1412 real(r8), intent(in) :: u(pcols,pver) ! U wind [ m/s ]
1413 real(r8), intent(in) :: v(pcols,pver) ! V wind [ m/s ]
1414 real(r8), intent(in) :: sl(pcols,pver) ! Liquid water static energy, cp * T + g * z - Lv * ql - Ls * qi [ J/kg ]
1415 real(r8), intent(in) :: slv(pcols,pver) ! Liquid water virtual static energy, sl * ( 1 + 0.608 * qt ) [ J/kg ]
1416 real(r8), intent(in) :: qt(pcols,pver) ! Total speccific humidity qv + ql + qi [ kg/kg ]
1417 real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ]
1418 real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressures [ Pa ]
1419 real(r8), intent(in) :: z(pcols,pver) ! Layer midpoint height above surface [ m ]
1420 real(r8), intent(in) :: zi(pcols,pver+1) ! Interface height above surface, i.e., zi(pver+1) = 0 all over the globe [ m ]
1421 real(r8), intent(in) :: chu(pcols,pver+1) ! Buoyancy coeffi. unsaturated sl (heat) coef. at all interfaces. [ unit ? ]
1422 real(r8), intent(in) :: chs(pcols,pver+1) ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces. [ unit ? ]
1423 real(r8), intent(in) :: cmu(pcols,pver+1) ! Buoyancy coeffi. unsaturated qt (moisture) coef. at all interfaces [ unit ? ]
1424 real(r8), intent(in) :: cms(pcols,pver+1) ! Buoyancy coeffi. saturated qt (moisture) coef. at all interfaces [ unit ? ]
1425 real(r8), intent(in) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ]
1426 real(r8), intent(in) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ]
1427 real(r8), intent(in) :: n2(pcols,pver) ! Interfacial (except surface) moist buoyancy frequency [ s-2 ]
1428 real(r8), intent(in) :: s2(pcols,pver) ! Interfacial (except surface) shear frequency [ s-2 ]
1429 real(r8), intent(in) :: ri(pcols,pver) ! Interfacial (except surface) Richardson number
1430 real(r8), intent(in) :: qflx(pcols) ! Kinematic surface constituent ( water vapor ) flux [ kg/m2/s ]
1431 real(r8), intent(in) :: shflx(pcols) ! Kinematic surface heat flux [ unit ? ]
1432 real(r8), intent(in) :: slslope(pcols,pver) ! Slope of 'sl' in each layer [ J/kg/Pa ]
1433 real(r8), intent(in) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer [ kg/kg/Pa ]
1434 real(r8), intent(in) :: qrlin(pcols,pver) ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
1435 real(r8), intent(in) :: wsedl(pcols,pver) ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
1436 real(r8), intent(in) :: ustar(pcols) ! Surface friction velocity [ m/s ]
1437 real(r8), intent(in) :: rrho(pcols) ! 1./bottom mid-point density. Specific volume [ m3/kg ]
1438 real(r8), intent(in) :: kvf(pcols,pver+1) ! Free atmosphere eddy diffusivity [ m2/s ]
1439 logical, intent(in) :: wstarent ! Switch for choosing wstar3 entrainment parameterization
1440 real(r8), intent(in) :: minpblh(pcols) ! Minimum PBL height based on surface stress [ m ]
1441 real(r8), intent(in) :: kvh_in(pcols,pver+1) ! kvh saved from last timestep or last iterative step [ m2/s ]
1442 real(r8), intent(in) :: kvm_in(pcols,pver+1) ! kvm saved from last timestep or last iterative step [ m2/s ]
1443 real(r8), intent(in) :: cld(pcols,pver) ! Stratus Cloud Fraction [ fraction ]
1445 ! ---------------- !
1446 ! Output variables !
1447 ! ---------------- !
1449 real(r8), intent(out) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat, moisture, and tracers [ m2/s ]
1450 real(r8), intent(out) :: kvm(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
1451 real(r8), intent(out) :: pblh(pcols) ! PBL top height [ m ]
1452 real(r8), intent(out) :: pblhp(pcols) ! PBL top height pressure [ Pa ]
1453 real(r8), intent(out) :: tpert(pcols) ! Convective temperature excess [ K ]
1454 real(r8), intent(out) :: qpert(pcols) ! Convective humidity excess [ kg/kg ]
1455 real(r8), intent(out) :: wpert(pcols) ! Turbulent velocity excess [ m/s ]
1456 real(r8), intent(out) :: tke(pcols,pver+1) ! Turbulent kinetic energy [ m2/s2 ], 'tkes' at surface, pver+1.
1457 real(r8), intent(out) :: bprod(pcols,pver+1) ! Buoyancy production [ m2/s3 ], 'bflxs' at surface, pver+1.
1458 real(r8), intent(out) :: sprod(pcols,pver+1) ! Shear production [ m2/s3 ], (ustar(i)**3)/(vk*z(i,pver)) at surface, pver+1.
1459 real(r8), intent(out) :: turbtype_f(pcols,pver+1) ! Turbulence type at each interface:
1460 ! 0. = Non turbulence interface
1461 ! 1. = Stable turbulence interface
1462 ! 2. = CL interior interface ( if bflxs > 0, surface is this )
1463 ! 3. = Bottom external interface of CL
1464 ! 4. = Top external interface of CL.
1465 ! 5. = Double entraining CL external interface
1466 real(r8), intent(out) :: sm_aw(pcols,pver+1) ! Galperin instability function of momentum for use in the microphysics [ no unit ]
1467 real(r8), intent(out) :: ipbl(pcols) ! If 1, PBL is CL, while if 0, PBL is STL.
1468 real(r8), intent(out) :: kpblh(pcols) ! Layer index containing PBL within or at the base interface
1470 ! --------------------------- !
1471 ! Diagnostic output variables !
1472 ! --------------------------- !
1474 real(r8) :: tkes(pcols) ! TKE at surface [ m2/s2 ]
1475 real(r8) :: kbase_o(pcols,ncvmax) ! Original external base interface index of CL just after 'exacol'
1476 real(r8) :: ktop_o(pcols,ncvmax) ! Original external top interface index of CL just after 'exacol'
1477 real(r8) :: ncvfin_o(pcols) ! Original number of CLs just after 'exacol'
1478 real(r8) :: kbase_mg(pcols,ncvmax) ! kbase just after extending-merging (after 'zisocl') but without SRCL
1479 real(r8) :: ktop_mg(pcols,ncvmax) ! ktop just after extending-merging (after 'zisocl') but without SRCL
1480 real(r8) :: ncvfin_mg(pcols) ! ncvfin just after extending-merging (after 'zisocl') but without SRCL
1481 real(r8) :: kbase_f(pcols,ncvmax) ! Final kbase after adding SRCL
1482 real(r8) :: ktop_f(pcols,ncvmax) ! Final ktop after adding SRCL
1483 real(r8) :: ncvfin_f(pcols) ! Final ncvfin after adding SRCL
1484 real(r8) :: wet_CL(pcols,ncvmax) ! Entrainment rate at the CL top [ m/s ]
1485 real(r8) :: web_CL(pcols,ncvmax) ! Entrainment rate at the CL base [ m/s ]
1486 real(r8) :: jtbu_CL(pcols,ncvmax) ! Buoyancy jump across the CL top [ m/s2 ]
1487 real(r8) :: jbbu_CL(pcols,ncvmax) ! Buoyancy jump across the CL base [ m/s2 ]
1488 real(r8) :: evhc_CL(pcols,ncvmax) ! Evaporative enhancement factor at the CL top
1489 real(r8) :: jt2slv_CL(pcols,ncvmax) ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
1490 real(r8) :: n2ht_CL(pcols,ncvmax) ! n2 defined at the CL top interface but using sfuh(kt) instead of sfi(kt) [ s-2 ]
1491 real(r8) :: n2hb_CL(pcols,ncvmax) ! n2 defined at the CL base interface but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
1492 real(r8) :: lwp_CL(pcols,ncvmax) ! LWP in the CL top layer [ kg/m2 ]
1493 real(r8) :: opt_depth_CL(pcols,ncvmax) ! Optical depth of the CL top layer
1494 real(r8) :: radinvfrac_CL(pcols,ncvmax) ! Fraction of LW radiative cooling confined in the top portion of CL
1495 real(r8) :: radf_CL(pcols,ncvmax) ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
1496 real(r8) :: wstar_CL(pcols,ncvmax) ! Convective velocity of CL including entrainment contribution finally [ m/s ]
1497 real(r8) :: wstar3fact_CL(pcols,ncvmax) ! "wstar3fact" of CL. Entrainment enhancement of wstar3 (inverse)
1499 real(r8) :: gh_a(pcols,pver+1) ! Half of normalized buoyancy production, -l2n2/2e. [ no unit ]
1500 real(r8) :: sh_a(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces [ no unit ]
1501 real(r8) :: sm_a(pcols,pver+1) ! Galperin instability function of momentum at all interfaces [ no unit ]
1502 real(r8) :: ri_a(pcols,pver+1) ! Interfacial Richardson number at all interfaces [ no unit ]
1504 real(r8) :: ebrk(pcols,ncvmax) ! Net CL mean TKE [ m2/s2 ]
1505 real(r8) :: wbrk(pcols,ncvmax) ! Net CL mean normalized TKE [ m2/s2 ]
1506 real(r8) :: lbrk(pcols,ncvmax) ! Net energetic integral thickness of CL [ m ]
1507 real(r8) :: ricl(pcols,ncvmax) ! Mean Richardson number of CL ( l2n2/l2s2 )
1508 real(r8) :: ghcl(pcols,ncvmax) ! Half of normalized buoyancy production of CL
1509 real(r8) :: shcl(pcols,ncvmax) ! Instability function of heat and moisture of CL
1510 real(r8) :: smcl(pcols,ncvmax) ! Instability function of momentum of CL
1512 real(r8) :: leng(pcols,pver+1) ! Turbulent length scale [ m ], 0 at the surface.
1513 real(r8) :: wcap(pcols,pver+1) ! Normalized TKE [m2/s2], 'tkes/b1' at the surface and 'tke/b1' at
1514 ! the top/bottom entrainment interfaces of CL assuming no transport.
1515 ! ------------------------ !
1516 ! Local Internal Variables !
1517 ! ------------------------ !
1519 logical :: belongcv(pcols,pver+1) ! True for interfaces in a CL (both interior and exterior are included)
1520 logical :: belongst(pcols,pver+1) ! True for stable turbulent layer interfaces (STL)
1521 logical :: in_CL ! True if interfaces k,k+1 both in same CL.
1522 logical :: extend ! True when CL is extended in zisocl
1523 logical :: extend_up ! True when CL is extended upward in zisocl
1524 logical :: extend_dn ! True when CL is extended downward in zisocl
1526 integer :: i ! Longitude index
1527 integer :: k ! Vertical index
1528 integer :: ks ! Vertical index
1529 integer :: ncvfin(pcols) ! Total number of CL in column
1530 integer :: ncvf ! Total number of CL in column prior to adding SRCL
1531 integer :: ncv ! Index of current CL
1532 integer :: ncvnew ! Index of added SRCL appended after regular CLs from 'zisocl'
1533 integer :: ncvsurf ! If nonzero, CL index based on surface (usually 1, but can be > 1 when SRCL is based at sfc)
1534 integer :: kbase(pcols,ncvmax) ! Vertical index of CL base interface
1535 integer :: ktop(pcols,ncvmax) ! Vertical index of CL top interface
1536 integer :: kb, kt ! kbase and ktop for current CL
1537 integer :: ktblw ! ktop of the CL located at just below the current CL
1538 integer :: turbtype(pcols,pver+1) ! Interface turbulence type :
1539 ! 0 = Non turbulence interface
1540 ! 1 = Stable turbulence interface
1541 ! 2 = CL interior interface ( if bflxs > 0, sfc is this )
1542 ! 3 = Bottom external interface of CL
1543 ! 4 = Top external interface of CL
1544 ! 5 = Double entraining CL external interface
1545 integer :: ktopbl(pcols) ! PBL top height or interface index
1546 real(r8) :: bflxs(pcols) ! Surface buoyancy flux [ m2/s3 ]
1547 real(r8) :: rcap ! 'tke/ebrk' at all interfaces of CL. Set to 1 at the CL entrainment interfaces
1548 real(r8) :: jtzm ! Interface layer thickness of CL top interface [ m ]
1549 real(r8) :: jtsl ! Jump of s_l across CL top interface [ J/kg ]
1550 real(r8) :: jtqt ! Jump of q_t across CL top interface [ kg/kg ]
1551 real(r8) :: jtbu ! Jump of buoyancy across CL top interface [ m/s2 ]
1552 real(r8) :: jtu ! Jump of u across CL top interface [ m/s ]
1553 real(r8) :: jtv ! Jump of v across CL top interface [ m/s ]
1554 real(r8) :: jt2slv ! Jump of slv ( across two layers ) at CL top for use only in evhc [ J/kg ]
1555 real(r8) :: radf ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
1556 real(r8) :: jbzm ! Interface layer thickness of CL base interface [ m ]
1557 real(r8) :: jbsl ! Jump of s_l across CL base interface [ J/kg ]
1558 real(r8) :: jbqt ! Jump of q_t across CL top interface [ kg/kg ]
1559 real(r8) :: jbbu ! Jump of buoyancy across CL base interface [ m/s2 ]
1560 real(r8) :: jbu ! Jump of u across CL base interface [ m/s ]
1561 real(r8) :: jbv ! Jump of v across CL base interface [ m/s ]
1562 real(r8) :: ch ! Buoyancy coefficients defined at the CL top and base interfaces using CL internal
1563 real(r8) :: cm ! sfuh(kt) and sflh(kb-1) instead of sfi(kt) and sfi(kb), respectively. These are
1564 ! used for entrainment calculation at CL external interfaces and SRCL identification.
1565 real(r8) :: n2ht ! n2 defined at the CL top interface but using sfuh(kt) instead of sfi(kt) [ s-2 ]
1566 real(r8) :: n2hb ! n2 defined at the CL base interface but using sflh(kb-1) instead of sfi(kb) [ s-2 ]
1567 real(r8) :: n2htSRCL ! n2 defined at the upper-half layer of SRCL. This is used only for identifying SRCL.
1568 ! n2htSRCL use SRCL internal slope sl and qt as well as sfuh(kt) instead of sfi(kt) [ s-2 ]
1569 real(r8) :: gh ! Half of normalized buoyancy production ( -l2n2/2e ) [ no unit ]
1570 real(r8) :: sh ! Galperin instability function for heat and moisture
1571 real(r8) :: sm ! Galperin instability function for momentum
1572 real(r8) :: lbulk ! Depth of turbulent layer, Master length scale (not energetic length)
1573 real(r8) :: dzht ! Thickness of top half-layer [ m ]
1574 real(r8) :: dzhb ! Thickness of bottom half-layer [ m ]
1575 real(r8) :: rootp ! Sqrt(net CL-mean TKE including entrainment contribution) [ m/s ]
1576 real(r8) :: evhc ! Evaporative enhancement factor: (1+E) with E = evap. cool. efficiency [ no unit ]
1577 real(r8) :: kentr ! Effective entrainment diffusivity 'wet*dz', 'web*dz' [ m2/s ]
1578 real(r8) :: lwp ! Liquid water path in the layer kt [ kg/m2 ]
1579 real(r8) :: opt_depth ! Optical depth of the layer kt [ no unit ]
1580 real(r8) :: radinvfrac ! Fraction of LW cooling in the layer kt concentrated at the CL top [ no unit ]
1581 real(r8) :: wet ! CL top entrainment rate [ m/s ]
1582 real(r8) :: web ! CL bot entrainment rate [ m/s ]. Set to zero if CL is based at surface.
1583 real(r8) :: vyt ! n2ht/n2 at the CL top interface
1584 real(r8) :: vyb ! n2hb/n2 at the CL base interface
1585 real(r8) :: vut ! Inverse Ri (=s2/n2) at the CL top interface
1586 real(r8) :: vub ! Inverse Ri (=s2/n2) at the CL base interface
1587 real(r8) :: fact ! Factor relating TKE generation to entrainment [ no unit ]
1588 real(r8) :: trma ! Intermediate variables used for solving quadratic ( for gh from ri )
1589 real(r8) :: trmb ! and cubic equations ( for ebrk: the net CL mean TKE )
1590 real(r8) :: trmc !
1591 real(r8) :: trmp !
1592 real(r8) :: trmq !
1593 real(r8) :: qq !
1594 real(r8) :: det !
1595 real(r8) :: gg ! Intermediate variable used for calculating stability functions of
1596 ! SRCL or SBCL based at the surface with bflxs > 0.
1597 real(r8) :: dzhb5 ! Half thickness of the bottom-most layer of current CL regime
1598 real(r8) :: dzht5 ! Half thickness of the top-most layer of adjacent CL regime just below current CL
1599 real(r8) :: qrlw(pcols,pver) ! Local grid-mean LW heating rate : [K/s] * cpair * dp = [ W/kg*Pa ]
1601 real(r8) :: cldeff(pcols,pver) ! Effective stratus fraction
1602 real(r8) :: qleff ! Used for computing evhc
1603 real(r8) :: tunlramp ! Ramping tunl
1604 real(r8) :: leng_imsi ! For Kv = max(Kv_STL, Kv_entrain)
1605 real(r8) :: tke_imsi !
1606 real(r8) :: kvh_imsi !
1607 real(r8) :: kvm_imsi !
1608 real(r8) :: alph4exs ! For extended stability function in the stable regime
1609 real(r8) :: ghmin !
1611 real(r8) :: sedfact ! For 'sedimentation-entrainment feedback'
1613 ! Local variables specific for 'wstar' entrainment closure
1615 real(r8) :: cet ! Proportionality coefficient between wet and wstar3
1616 real(r8) :: ceb ! Proportionality coefficient between web and wstar3
1617 real(r8) :: wstar ! Convective velocity for CL [ m/s ]
1618 real(r8) :: wstar3 ! Cubed convective velocity for CL [ m3/s3 ]
1619 real(r8) :: wstar3fact ! 1/(relative change of wstar^3 by entrainment)
1620 real(r8) :: rmin ! sqrt(p)
1621 real(r8) :: fmin ! f(rmin), where f(r) = r^3 - 3*p*r - 2q
1622 real(r8) :: rcrit ! ccrit*wstar
1623 real(r8) :: fcrit ! f(rcrit)
1624 logical noroot ! True if f(r) has no root r > rcrit
1626 !-----------------------!
1627 ! Start of Main Program !
1628 !-----------------------!
1630 ! Option: Turn-off LW radiative-turbulence interaction in PBL scheme
1631 ! by setting qrlw = 0. Logical parameter 'set_qrlzero' was
1632 ! defined in the first part of 'eddy_diff.F90' module.
1634 if( set_qrlzero ) then
1635 qrlw(:,:) = 0._r8
1636 else
1637 qrlw(:ncol,:pver) = qrlin(:ncol,:pver)
1638 endif
1640 ! Define effective stratus fraction using the grid-mean ql.
1641 ! Modification : The contribution of ice should be carefully considered.
1642 ! This should be done in combination with the 'qrlw' and
1643 ! overlapping assumption of liquid and ice stratus.
1645 do k = 1, pver
1646 do i = 1, ncol
1647 if( choice_evhc .eq. 'ramp' .or. choice_radf .eq. 'ramp' ) then
1648 cldeff(i,k) = cld(i,k) * min( ql(i,k) / qmin, 1._r8 )
1649 else
1650 cldeff(i,k) = cld(i,k)
1651 endif
1652 end do
1653 end do
1655 ! For an extended stability function in the stable regime, re-define
1656 ! alph4exe and ghmin. This is for future work.
1658 if( ricrit .eq. 0.19_r8 ) then
1659 alph4exs = alph4
1660 ghmin = -3.5334_r8
1661 elseif( ricrit .gt. 0.19_r8 ) then
1662 alph4exs = -2._r8 * b1 * alph2 / ( alph3 - 2._r8 * b1 * alph5 ) / ricrit
1663 ghmin = -1.e10_r8
1664 else
1665 write(iulog,*) 'Error : ricrit should be larger than 0.19 in UW PBL'
1666 #ifdef WRF_PORT
1667 call wrf_message(iulog)
1668 #endif
1669 stop
1670 endif
1672 !
1673 ! Initialization of Diagnostic Output
1674 !
1676 do i = 1, ncol
1677 wet_CL(i,:ncvmax) = 0._r8
1678 web_CL(i,:ncvmax) = 0._r8
1679 jtbu_CL(i,:ncvmax) = 0._r8
1680 jbbu_CL(i,:ncvmax) = 0._r8
1681 evhc_CL(i,:ncvmax) = 0._r8
1682 jt2slv_CL(i,:ncvmax) = 0._r8
1683 n2ht_CL(i,:ncvmax) = 0._r8
1684 n2hb_CL(i,:ncvmax) = 0._r8
1685 lwp_CL(i,:ncvmax) = 0._r8
1686 opt_depth_CL(i,:ncvmax) = 0._r8
1687 radinvfrac_CL(i,:ncvmax) = 0._r8
1688 radf_CL(i,:ncvmax) = 0._r8
1689 wstar_CL(i,:ncvmax) = 0._r8
1690 wstar3fact_CL(i,:ncvmax) = 0._r8
1691 ricl(i,:ncvmax) = 0._r8
1692 ghcl(i,:ncvmax) = 0._r8
1693 shcl(i,:ncvmax) = 0._r8
1694 smcl(i,:ncvmax) = 0._r8
1695 ebrk(i,:ncvmax) = 0._r8
1696 wbrk(i,:ncvmax) = 0._r8
1697 lbrk(i,:ncvmax) = 0._r8
1698 gh_a(i,:pver+1) = 0._r8
1699 sh_a(i,:pver+1) = 0._r8
1700 sm_a(i,:pver+1) = 0._r8
1701 ri_a(i,:pver+1) = 0._r8
1702 sm_aw(i,:pver+1) = 0._r8
1703 ipbl(i) = 0._r8
1704 kpblh(i) = real(pver,r8)
1705 end do
1707 ! kvh and kvm are stored over timesteps in 'vertical_diffusion.F90' and
1708 ! passed in as kvh_in and kvm_in. However, at the first timestep they
1709 ! need to be computed and these are done just before calling 'caleddy'.
1710 ! kvm and kvh are also stored over iterative time step in the first part
1711 ! of 'eddy_diff.F90'
1713 do k = 1, pver + 1
1714 do i = 1, ncol
1715 ! Initialize kvh and kvm to zero or kvf
1716 if( use_kvf ) then
1717 kvh(i,k) = kvf(i,k)
1718 kvm(i,k) = kvf(i,k)
1719 else
1720 kvh(i,k) = 0._r8
1721 kvm(i,k) = 0._r8
1722 end if
1723 ! Zero diagnostic quantities for the new diffusion step.
1724 wcap(i,k) = 0._r8
1725 leng(i,k) = 0._r8
1726 tke(i,k) = 0._r8
1727 turbtype(i,k) = 0
1728 end do
1729 end do
1731 ! Initialize 'bprod' [ m2/s3 ] and 'sprod' [ m2/s3 ] at all interfaces.
1732 ! Note this initialization is a hybrid initialization since 'n2' [s-2] and 's2' [s-2]
1733 ! are calculated from the given current initial profile, while 'kvh_in' [m2/s] and
1734 ! 'kvm_in' [m2/s] are from the previous iteration or previous time step.
1735 ! This initially guessed 'bprod' and 'sprod' will be updated at the end of this
1736 ! 'caleddy' subroutine for diagnostic output.
1737 ! This computation of 'brpod,sprod' below is necessary for wstar-based entrainment closure.
1739 do k = 2, pver
1740 do i = 1, ncol
1741 bprod(i,k) = -kvh_in(i,k) * n2(i,k)
1742 sprod(i,k) = kvm_in(i,k) * s2(i,k)
1743 end do
1744 end do
1746 ! Set 'bprod' and 'sprod' at top and bottom interface.
1747 ! In calculating 'surface' (actually lowest half-layer) buoyancy flux,
1748 ! 'chu' at surface is defined to be the same as 'chu' at the mid-point
1749 ! of lowest model layer (pver) at the end of 'trbind'. The same is for
1750 ! the other buoyancy coefficients. 'sprod(i,pver+1)' is defined in a
1751 ! consistent way as the definition of 'tkes' in the original code.
1752 ! ( Important Option ) If I want to isolate surface buoyancy flux from
1753 ! the other parts of CL regimes energetically even though bflxs > 0,
1754 ! all I should do is to re-define 'bprod(i,pver+1)=0' in the below 'do'
1755 ! block. Additionally for merging test of extending SBCL based on 'l2n2'
1756 ! in 'zisocl', I should use 'l2n2 = - wint / sh' for similar treatment
1757 ! as previous code. All other parts of the code are fully consistently
1758 ! treated by these change only.
1759 ! My future general convection scheme will use bflxs(i).
1761 do i = 1, ncol
1762 bprod(i,1) = 0._r8 ! Top interface
1763 sprod(i,1) = 0._r8 ! Top interface
1764 ch = chu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + chs(i,pver+1) * sflh(i,pver)
1765 cm = cmu(i,pver+1) * ( 1._r8 - sflh(i,pver) ) + cms(i,pver+1) * sflh(i,pver)
1766 bflxs(i) = ch * shflx(i) * rrho(i) + cm * qflx(i) * rrho(i)
1767 if( choice_tkes .eq. 'ibprod' ) then
1768 bprod(i,pver+1) = bflxs(i)
1769 else
1770 bprod(i,pver+1) = 0._r8
1771 endif
1772 sprod(i,pver+1) = (ustar(i)**3)/(vk*z(i,pver))
1773 end do
1775 ! Initially identify CL regimes in 'exacol'
1776 ! ktop : Interface index of the CL top external interface
1777 ! kbase : Interface index of the CL base external interface
1778 ! ncvfin: Number of total CLs
1779 ! Note that if surface buoyancy flux is positive ( bflxs = bprod(i,pver+1) > 0 ),
1780 ! surface interface is identified as an internal interface of CL. However, even
1781 ! though bflxs <= 0, if 'pver' interface is a CL internal interface (ri(pver)<0),
1782 ! surface interface is identified as an external interface of CL. If bflxs =< 0
1783 ! and ri(pver) >= 0, then surface interface is identified as a stable turbulent
1784 ! intereface (STL) as shown at the end of 'caleddy'. Even though a 'minpblh' is
1785 ! passed into 'exacol', it is not used in the 'exacol'.
1787 call exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin )
1789 ! Diagnostic output of CL interface indices before performing 'extending-merging'
1790 ! of CL regimes in 'zisocl'
1791 do i = 1, ncol
1792 do k = 1, ncvmax
1793 kbase_o(i,k) = real(kbase(i,k),r8)
1794 ktop_o(i,k) = real(ktop(i,k),r8)
1795 ncvfin_o(i) = real(ncvfin(i),r8)
1796 end do
1797 end do
1799 ! ----------------------------------- !
1800 ! Perform calculation for each column !
1801 ! ----------------------------------- !
1803 do i = 1, ncol
1805 ! Define Surface Interfacial Layer TKE, 'tkes'.
1806 ! In the current code, 'tkes' is used as representing TKE of surface interfacial
1807 ! layer (low half-layer of surface-based grid layer). In the code, when bflxs>0,
1808 ! surface interfacial layer is assumed to be energetically coupled to the other
1809 ! parts of the CL regime based at the surface. In this sense, it is conceptually
1810 ! more reasonable to include both 'bprod' and 'sprod' in the definition of 'tkes'.
1811 ! Since 'tkes' cannot be negative, it is lower bounded by small positive number.
1812 ! Note that inclusion of 'bprod' in the definition of 'tkes' may increase 'ebrk'
1813 ! and 'wstar3', and eventually, 'wet' at the CL top, especially when 'bflxs>0'.
1814 ! This might help to solve the problem of too shallow PBLH over the overcast Sc
1815 ! regime. If I want to exclude 'bprod(i,pver+1)' in calculating 'tkes' even when
1816 ! bflxs > 0, all I should to do is to set 'bprod(i,pver+1) = 0' in the above
1817 ! initialization 'do' loop (explained above), NOT changing the formulation of
1818 ! tkes(i) in the below block. This is because for consistent treatment in the
1819 ! other parts of the code also.
1821 ! tkes(i) = (b1*vk*z(i,pver)*sprod(i,pver+1))**(2._r8/3._r8)
1822 tkes(i) = max(b1*vk*z(i,pver)*(bprod(i,pver+1)+sprod(i,pver+1)), 1.e-7_r8)**(2._r8/3._r8)
1823 tkes(i) = min(tkes(i), tkemax)
1824 tke(i,pver+1) = tkes(i)
1825 wcap(i,pver+1) = tkes(i)/b1
1827 ! Extend and merge the initially identified CLs, relabel the CLs, and calculate
1828 ! CL internal mean energetics and stability functions in 'zisocl'.
1829 ! The CL nearest to the surface is CL(1) and the CL index, ncv, increases
1830 ! with height. The following outputs are from 'zisocl'. Here, the dimension
1831 ! of below outputs are (pcols,ncvmax) (except the 'ncvfin(pcols)' and
1832 ! 'belongcv(pcols,pver+1)) and 'ncv' goes from 1 to 'ncvfin'.
1833 ! For 'ncv = ncvfin+1, ncvmax', below output are already initialized to be zero.
1834 ! ncvfin : Total number of CLs
1835 ! kbase(ncv) : Base external interface index of CL
1836 ! ktop : Top external interface index of CL
1837 ! belongcv : True if the interface (either internal or external) is CL
1838 ! ricl : Mean Richardson number of internal CL
1839 ! ghcl : Normalized buoyancy production '-l2n2/2e' [no unit] of internal CL
1840 ! shcl : Galperin instability function of heat-moisture of internal CL
1841 ! smcl : Galperin instability function of momentum of internal CL
1842 ! lbrk, <l>int : Thickness of (energetically) internal CL (lint, [m])
1843 ! wbrk, <W>int : Mean normalized TKE of internal CL ([m2/s2])
1844 ! ebrk, <e>int : Mean TKE of internal CL (b1*wbrk,[m2/s2])
1845 ! The ncvsurf is an identifier saying which CL regime is based at the surface.
1846 ! If 'ncvsurf=1', then the first CL regime is based at the surface. If surface
1847 ! interface is not a part of CL (neither internal nor external), 'ncvsurf = 0'.
1848 ! After identifying and including SRCLs into the normal CL regimes (where newly
1849 ! identified SRCLs are simply appended to the normal CL regimes using regime
1850 ! indices of 'ncvfin+1','ncvfin+2' (as will be shown in the below SRCL part),..
1851 ! where 'ncvfin' is the final CL regime index produced after extending-merging
1852 ! in 'zisocl' but before adding SRCLs), if any newly identified SRCL (e.g.,
1853 ! 'ncvfin+1') is based at surface, then 'ncvsurf = ncvfin+1'. Thus 'ncvsurf' can
1854 ! be 0, 1, or >1. 'ncvsurf' can be a useful diagnostic output.
1856 ncvsurf = 0
1857 if( ncvfin(i) .gt. 0 ) then
1858 call zisocl( pcols , pver , i , &
1859 z , zi , n2 , s2 , &
1860 bprod , sprod , bflxs , tkes , &
1861 ncvfin , kbase , ktop , belongcv, &
1862 ricl , ghcl , shcl , smcl , &
1863 lbrk , wbrk , ebrk , &
1864 extend , extend_up, extend_dn )
1865 if( kbase(i,1) .eq. pver + 1 ) ncvsurf = 1
1866 else
1867 belongcv(i,:) = .false.
1868 endif
1870 ! Diagnostic output after finishing extending-merging process in 'zisocl'
1871 ! Since we are adding SRCL additionally, we need to print out these here.
1873 do k = 1, ncvmax
1874 kbase_mg(i,k) = real(kbase(i,k))
1875 ktop_mg(i,k) = real(ktop(i,k))
1876 ncvfin_mg(i) = real(ncvfin(i))
1877 end do
1879 ! ----------------------- !
1880 ! Identification of SRCLs !
1881 ! ----------------------- !
1883 ! Modification : This cannot identify the 'cirrus' layer due to the condition of
1884 ! ql(i,k) .gt. qmin. This should be modified in future to identify
1885 ! a single thin cirrus layer.
1886 ! Instead of ql, we may use cldn in future, including ice
1887 ! contribution.
1889 ! ------------------------------------------------------------------------------ !
1890 ! Find single-layer radiatively-driven cloud-topped convective layers (SRCLs). !
1891 ! SRCLs extend through a single model layer k, with entrainment at the top and !
1892 ! bottom interfaces, unless bottom interface is the surface. !
1893 ! The conditions for an SRCL is identified are: !
1894 ! !
1895 ! 1. Cloud in the layer, k : ql(i,k) .gt. qmin = 1.e-5 [ kg/kg ] !
1896 ! 2. No cloud in the above layer (else assuming that some fraction of the LW !
1897 ! flux divergence in layer k is concentrated at just below top interface !
1898 ! of layer k is invalid). Then, this condition might be sensitive to the !
1899 ! vertical resolution of grid. !
1900 ! 3. LW radiative cooling (SW heating is assumed uniformly distributed through !
1901 ! layer k, so not relevant to buoyancy production) in the layer k. However, !
1902 ! SW production might also contribute, which may be considered in a future. !
1903 ! 4. Internal stratification 'n2ht' of upper-half layer should be unstable. !
1904 ! The 'n2ht' is pure internal stratification of upper half layer, obtained !
1905 ! using internal slopes of sl, qt in layer k (in contrast to conventional !
1906 ! interfacial slope) and saturation fraction in the upper-half layer, !
1907 ! sfuh(k) (in contrast to sfi(k)). !
1908 ! 5. Top and bottom interfaces not both in the same existing convective layer. !
1909 ! If SRCL is within the previouisly identified CL regimes, we don't define !
1910 ! a new SRCL. !
1911 ! 6. k >= ntop_turb + 1 = 2 !
1912 ! 7. Ri at the top interface > ricrit = 0.19 (otherwise turbulent mixing will !
1913 ! broadly distribute the cloud top in the vertical, preventing localized !
1914 ! radiative destabilization at the top interface). !
1915 ! !
1916 ! Note if 'k = pver', it identifies a surface-based single fog layer, possibly, !
1917 ! warm advection fog. Note also the CL regime index of SRCLs itself increases !
1918 ! with height similar to the regular CLs indices identified from 'zisocl'. !
1919 ! ------------------------------------------------------------------------------ !
1921 ncv = 1
1922 ncvf = ncvfin(i)
1924 if( choice_SRCL .eq. 'remove' ) goto 222
1926 do k = nbot_turb, ntop_turb + 1, -1 ! 'k = pver, 2, -1' is a layer index.
1928 if( ql(i,k) .gt. qmin .and. ql(i,k-1) .lt. qmin .and. qrlw(i,k) .lt. 0._r8 &
1929 .and. ri(i,k) .ge. ricrit ) then
1931 ! In order to avoid any confliction with the treatment of ambiguous layer,
1932 ! I need to impose an additional constraint that ambiguous layer cannot be
1933 ! SRCL. So, I added constraint that 'k+1' interface (base interface of k
1934 ! layer) should not be a part of previously identified CL. Since 'belongcv'
1935 ! is even true for external entrainment interfaces, below constraint is
1936 ! fully sufficient.
1938 if( choice_SRCL .eq. 'nonamb' .and. belongcv(i,k+1) ) then
1939 go to 220
1940 endif
1942 ch = ( 1._r8 - sfuh(i,k) ) * chu(i,k) + sfuh(i,k) * chs(i,k)
1943 cm = ( 1._r8 - sfuh(i,k) ) * cmu(i,k) + sfuh(i,k) * cms(i,k)
1945 n2htSRCL = ch * slslope(i,k) + cm * qtslope(i,k)
1947 if( n2htSRCL .le. 0._r8 ) then
1949 ! Test if bottom and top interfaces are part of the pre-existing CL.
1950 ! If not, find appropriate index for the new SRCL. Note that this
1951 ! calculation makes use of 'ncv set' obtained from 'zisocl'. The
1952 ! 'in_CL' is a parameter testing whether the new SRCL is already
1953 ! within the pre-existing CLs (.true.) or not (.false.).
1955 in_CL = .false.
1957 do while ( ncv .le. ncvf )
1958 if( ktop(i,ncv) .le. k ) then
1959 if( kbase(i,ncv) .gt. k ) then
1960 in_CL = .true.
1961 endif
1962 exit ! Exit from 'do while' loop if SRCL is within the CLs.
1963 else
1964 ncv = ncv + 1 ! Go up one CL
1965 end if
1966 end do ! ncv
1968 if( .not. in_CL ) then ! SRCL is not within the pre-existing CLs.
1970 ! Identify a new SRCL and add it to the pre-existing CL regime group.
1972 ncvfin(i) = ncvfin(i) + 1
1973 ncvnew = ncvfin(i)
1974 ktop(i,ncvnew) = k
1975 kbase(i,ncvnew) = k+1
1976 belongcv(i,k) = .true.
1977 belongcv(i,k+1) = .true.
1979 ! Calculate internal energy of SRCL. There is no internal energy if
1980 ! SRCL is elevated from the surface. Also, we simply assume neutral
1981 ! stability function. Note that this assumption of neutral stability
1982 ! does not influence numerical calculation- stability functions here
1983 ! are just for diagnostic output. In general SRCLs other than a SRCL
1984 ! based at surface with bflxs <= 0, there is no other way but to use
1985 ! neutral stability function. However, in case of SRCL based at the
1986 ! surface, we can explicitly calculate non-zero stability functions
1987 ! in a consistent way. Even though stability functions of SRCL are
1988 ! just diagnostic outputs not influencing numerical calculations, it
1989 ! would be informative to write out correct reasonable values rather
1990 ! than simply assuming neutral stability. I am doing this right now.
1991 ! Similar calculations were done for the SBCL and when surface inter
1992 ! facial layer was merged by overlying CL in 'ziscol'.
1994 if( k .lt. pver ) then
1996 wbrk(i,ncvnew) = 0._r8
1997 ebrk(i,ncvnew) = 0._r8
1998 lbrk(i,ncvnew) = 0._r8
1999 ghcl(i,ncvnew) = 0._r8
2000 shcl(i,ncvnew) = 0._r8
2001 smcl(i,ncvnew) = 0._r8
2002 ricl(i,ncvnew) = 0._r8
2004 else ! Surface-based fog
2006 if( bflxs(i) .gt. 0._r8 ) then ! Incorporate surface TKE into CL interior energy
2007 ! It is likely that this case cannot exist since
2008 ! if surface buoyancy flux is positive, it would
2009 ! have been identified as SBCL in 'zisocl' ahead.
2010 ebrk(i,ncvnew) = tkes(i)
2011 lbrk(i,ncvnew) = z(i,pver)
2012 wbrk(i,ncvnew) = tkes(i) / b1
2014 write(iulog,*) 'Major mistake in SRCL: bflxs > 0 for surface-based SRCL'
2015 #ifdef WRF_PORT
2016 call wrf_message(iulog)
2017 #endif
2018 write(iulog,*) 'bflxs = ', bflxs(i)
2019 #ifdef WRF_PORT
2020 call wrf_message(iulog)
2021 #endif
2022 write(iulog,*) 'ncvfin_o = ', ncvfin_o(i)
2023 #ifdef WRF_PORT
2024 call wrf_message(iulog)
2025 #endif
2026 write(iulog,*) 'ncvfin_mg = ', ncvfin_mg(i)
2027 #ifdef WRF_PORT
2028 call wrf_message(iulog)
2029 #endif
2030 do ks = 1, ncvmax
2031 write(iulog,*) 'ncv =', ks, ' ', kbase_o(i,ks), ktop_o(i,ks), kbase_mg(i,ks), ktop_mg(i,ks)
2032 #ifdef WRF_PORT
2033 call wrf_message(iulog)
2034 #endif
2035 end do
2036 stop
2038 else ! Don't incorporate surface interfacial TKE into CL interior energy
2040 ebrk(i,ncvnew) = 0._r8
2041 lbrk(i,ncvnew) = 0._r8
2042 wbrk(i,ncvnew) = 0._r8
2044 endif
2046 ! Calculate stability functions (ghcl, shcl, smcl, ricl) explicitly
2047 ! using an reverse procedure starting from tkes(i). Note that it is
2048 ! possible to calculate stability functions even when bflxs < 0.
2049 ! Previous code just assumed neutral stability functions. Note that
2050 ! since alph5 = 0.7 > 0, alph3 = -35 < 0, the denominator of gh is
2051 ! always positive if bflxs > 0. However, if bflxs < 0, denominator
2052 ! can be zero. For this case, we provide a possible maximum negative
2053 ! value (the most stable state) to gh. Note also tkes(i) is always a
2054 ! positive value by a limiter. Also, sprod(i,pver+1) > 0 by limiter.
2056 gg = 0.5_r8 * vk * z(i,pver) * bprod(i,pver+1) / ( tkes(i)**(3._r8/2._r8) )
2057 if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
2058 ! gh = -0.28_r8
2059 ! gh = -3.5334_r8
2060 gh = ghmin
2061 else
2062 gh = gg / ( alph5 - gg * alph3 )
2063 end if
2064 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
2065 ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
2066 gh = min(max(gh,ghmin),0.0233_r8)
2067 ghcl(i,ncvnew) = gh
2068 shcl(i,ncvnew) = max(0._r8,alph5/(1._r8+alph3*gh))
2069 smcl(i,ncvnew) = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
2070 ricl(i,ncvnew) = -(smcl(i,ncvnew)/shcl(i,ncvnew))*(bprod(i,pver+1)/sprod(i,pver+1))
2072 ! 'ncvsurf' is CL regime index based at the surface. If there is no
2073 ! such regime, then 'ncvsurf = 0'.
2075 ncvsurf = ncvnew
2077 end if
2079 end if
2081 end if
2083 end if
2085 220 continue
2087 end do ! End of 'k' loop where 'k' is a grid layer index running from 'pver' to 2
2089 222 continue
2091 ! -------------------------------------------------------------------------- !
2092 ! Up to this point, we identified all kinds of CL regimes : !
2093 ! 1. A SBCL. By construction, 'bflxs > 0' for SBCL. !
2094 ! 2. Surface-based CL with multiple layers and 'bflxs =< 0' !
2095 ! 3. Surface-based CL with multiple layers and 'bflxs > 0' !
2096 ! 4. Regular elevated CL with two entraining interfaces !
2097 ! 5. SRCLs. If SRCL is based at surface, it will be bflxs < 0. !
2098 ! '1-4' were identified from 'zisocl' while '5' were identified separately !
2099 ! after performing 'zisocl'. CL regime index of '1-4' increases with height !
2100 ! ( e.g., CL = 1 is the CL regime nearest to the surface ) while CL regime !
2101 ! index of SRCL is simply appended after the final index of CL regimes from !
2102 ! 'zisocl'. However, CL regime indices of SRCLs itself increases with height !
2103 ! when there are multiple SRCLs, similar to the regular CLs from 'zisocl'. !
2104 ! -------------------------------------------------------------------------- !
2106 ! Diagnostic output of final CL regimes indices
2108 do k = 1, ncvmax
2109 kbase_f(i,k) = real(kbase(i,k))
2110 ktop_f(i,k) = real(ktop(i,k))
2111 ncvfin_f(i) = real(ncvfin(i))
2112 end do
2114 ! ---------------------------------------- !
2115 ! Perform do loop for individual CL regime !
2116 ! ---------------------------------------- ! -------------------------------- !
2117 ! For individual CLs, compute !
2118 ! 1. Entrainment rates at the CL top and (if any) base interfaces using !
2119 ! appropriate entrainment closure (current code use 'wstar' closure). !
2120 ! 2. Net CL mean (i.e., including entrainment contribution) TKE (ebrk) !
2121 ! and normalized TKE (wbrk). !
2122 ! 3. TKE (tke) and normalized TKE (wcap) profiles at all CL interfaces. !
2123 ! 4. ( kvm, kvh ) profiles at all CL interfaces. !
2124 ! 5. ( bprod, sprod ) profiles at all CL interfaces. !
2125 ! Also calculate !
2126 ! 1. PBL height as the top external interface of surface-based CL, if any. !
2127 ! 2. Characteristic excesses of convective 'updraft velocity (wpert)', !
2128 ! 'temperature (tpert)', and 'moisture (qpert)' in the surface-based CL, !
2129 ! if any, for use in the separate convection scheme. !
2130 ! If there is no surface-based CL, 'PBL height' and 'convective excesses' are !
2131 ! calculated later from surface-based STL (Stable Turbulent Layer) properties.!
2132 ! --------------------------------------------------------------------------- !
2134 ktblw = 0
2135 do ncv = 1, ncvfin(i)
2137 kt = ktop(i,ncv)
2138 kb = kbase(i,ncv)
2139 ! Check whether surface interface is energetically interior or not.
2140 if( kb .eq. (pver+1) .and. bflxs(i) .le. 0._r8 ) then
2141 lbulk = zi(i,kt) - z(i,pver)
2142 else
2143 lbulk = zi(i,kt) - zi(i,kb)
2144 end if
2146 ! Calculate 'turbulent length scale (leng)' and 'normalized TKE (wcap)'
2147 ! at all CL interfaces except the surface. Note that below 'wcap' at
2148 ! external interfaces are not correct. However, it does not influence
2149 ! numerical calculation and correct normalized TKE at the entraining
2150 ! interfaces will be re-calculated at the end of this 'do ncv' loop.
2152 do k = min(kb,pver), kt, -1
2153 if( choice_tunl .eq. 'rampcl' ) then
2154 ! In order to treat the case of 'ricl(i,ncv) >> 0' of surface-based SRCL
2155 ! with 'bflxs(i) < 0._r8', I changed ricl(i,ncv) -> min(0._r8,ricl(i,ncv))
2156 ! in the below exponential. This is necessary to prevent the model crash
2157 ! by too large values (e.g., 700) of ricl(i,ncv)
2158 tunlramp = ctunl*tunl*(1._r8-(1._r8-1._r8/ctunl)*exp(min(0._r8,ricl(i,ncv))))
2159 tunlramp = min(max(tunlramp,tunl),ctunl*tunl)
2160 elseif( choice_tunl .eq. 'rampsl' ) then
2161 tunlramp = ctunl*tunl
2162 ! tunlramp = 0.765_r8
2163 else
2164 tunlramp = tunl
2165 endif
2166 if( choice_leng .eq. 'origin' ) then
2167 leng(i,k) = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
2168 ! leng(i,k) = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
2169 else
2170 leng(i,k) = min( vk*zi(i,k), tunlramp*lbulk )
2171 endif
2172 wcap(i,k) = (leng(i,k)**2) * (-shcl(i,ncv)*n2(i,k)+smcl(i,ncv)*s2(i,k))
2173 end do ! k
2175 ! Calculate basic cross-interface variables ( jump condition ) across the
2176 ! base external interface of CL.
2178 if( kb .lt. pver+1 ) then
2180 jbzm = z(i,kb-1) - z(i,kb) ! Interfacial layer thickness [m]
2181 jbsl = sl(i,kb-1) - sl(i,kb) ! Interfacial jump of 'sl' [J/kg]
2182 jbqt = qt(i,kb-1) - qt(i,kb) ! Interfacial jump of 'qt' [kg/kg]
2183 jbbu = n2(i,kb) * jbzm ! Interfacial buoyancy jump [m/s2] considering saturation ( > 0 )
2184 jbbu = max(jbbu,jbumin) ! Set minimum buoyancy jump, jbumin = 1.e-3
2185 jbu = u(i,kb-1) - u(i,kb) ! Interfacial jump of 'u' [m/s]
2186 jbv = v(i,kb-1) - v(i,kb) ! Interfacial jump of 'v' [m/s]
2187 ch = (1._r8 -sflh(i,kb-1))*chu(i,kb) + sflh(i,kb-1)*chs(i,kb) ! Buoyancy coefficient just above the base interface
2188 cm = (1._r8 -sflh(i,kb-1))*cmu(i,kb) + sflh(i,kb-1)*cms(i,kb) ! Buoyancy coefficient just above the base interface
2189 n2hb = (ch*jbsl + cm*jbqt)/jbzm ! Buoyancy frequency [s-2] just above the base interface
2190 vyb = n2hb*jbzm/jbbu ! Ratio of 'n2hb/n2' at 'kb' interface
2191 vub = min(1._r8,(jbu**2+jbv**2)/(jbbu*jbzm) ) ! Ratio of 's2/n2 = 1/Ri' at 'kb' interface
2193 else
2195 ! Below setting is necessary for consistent treatment when 'kb' is at the surface.
2196 jbbu = 0._r8
2197 n2hb = 0._r8
2198 vyb = 0._r8
2199 vub = 0._r8
2200 web = 0._r8
2202 end if
2204 ! Calculate basic cross-interface variables ( jump condition ) across the
2205 ! top external interface of CL. The meanings of variables are similar to
2206 ! the ones at the base interface.
2208 jtzm = z(i,kt-1) - z(i,kt)
2209 jtsl = sl(i,kt-1) - sl(i,kt)
2210 jtqt = qt(i,kt-1) - qt(i,kt)
2211 jtbu = n2(i,kt)*jtzm ! Note : 'jtbu' is guaranteed positive by definition of CL top.
2212 jtbu = max(jtbu,jbumin) ! But threshold it anyway to be sure.
2213 jtu = u(i,kt-1) - u(i,kt)
2214 jtv = v(i,kt-1) - v(i,kt)
2215 ch = (1._r8 -sfuh(i,kt))*chu(i,kt) + sfuh(i,kt)*chs(i,kt)
2216 cm = (1._r8 -sfuh(i,kt))*cmu(i,kt) + sfuh(i,kt)*cms(i,kt)
2217 n2ht = (ch*jtsl + cm*jtqt)/jtzm
2218 vyt = n2ht*jtzm/jtbu
2219 vut = min(1._r8,(jtu**2+jtv**2)/(jtbu*jtzm))
2221 ! Evaporative enhancement factor of entrainment rate at the CL top interface, evhc.
2222 ! We take the full inversion strength to be 'jt2slv = slv(i,kt-2)-slv(i,kt)'
2223 ! where 'kt-1' is in the ambiguous layer. However, for a cloud-topped CL overlain
2224 ! by another CL, it is possible that 'slv(i,kt-2) < slv(i,kt)'. To avoid negative
2225 ! or excessive evhc, we lower-bound jt2slv and upper-bound evhc. Note 'jtslv' is
2226 ! used only for calculating 'evhc' : when calculating entrainment rate, we will
2227 ! use normal interfacial buoyancy jump across CL top interface.
2229 evhc = 1._r8
2230 jt2slv = 0._r8
2232 ! Modification : I should check whether below 'jbumin' produces reasonable limiting value.
2233 ! In addition, our current formulation does not consider ice contribution.
2235 if( choice_evhc .eq. 'orig' ) then
2237 if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
2238 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2239 jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
2240 evhc = 1._r8 + a2l * a3l * latvap * ql(i,kt) / jt2slv
2241 evhc = min( evhc, evhcmax )
2242 end if
2244 elseif( choice_evhc .eq. 'ramp' ) then
2246 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2247 jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
2248 evhc = 1._r8 + max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * a2l * a3l * latvap * ql(i,kt) / jt2slv
2249 evhc = min( evhc, evhcmax )
2251 elseif( choice_evhc .eq. 'maxi' ) then
2253 qleff = max( ql(i,kt-1), ql(i,kt) )
2254 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2255 jt2slv = max( jt2slv, jbumin*slv(i,kt-1)/g )
2256 evhc = 1._r8 + a2l * a3l * latvap * qleff / jt2slv
2257 evhc = min( evhc, evhcmax )
2259 endif
2261 ! Calculate cloud-top radiative cooling contribution to buoyancy production.
2262 ! Here, 'radf' [m2/s3] is additional buoyancy flux at the CL top interface
2263 ! associated with cloud-top LW cooling being mainly concentrated near the CL
2264 ! top interface ( just below CL top interface ). Contribution of SW heating
2265 ! within the cloud is not included in this radiative buoyancy production
2266 ! since SW heating is more broadly distributed throughout the CL top layer.
2268 lwp = 0._r8
2269 opt_depth = 0._r8
2270 radinvfrac = 0._r8
2271 radf = 0._r8
2273 if( choice_radf .eq. 'orig' ) then
2275 if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
2277 lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
2278 opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
2280 ! Approximate LW cooling fraction concentrated at the inversion by using
2281 ! polynomial approx to exact formula 1-2/opt_depth+2/(exp(opt_depth)-1))
2283 radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
2284 radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling = [ W/kg ]
2285 radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
2286 ! We can disable cloud LW cooling contribution to turbulence by uncommenting:
2287 ! radf = 0._r8
2289 end if
2291 elseif( choice_radf .eq. 'ramp' ) then
2293 lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
2294 opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
2295 radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
2296 radinvfrac = max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * radinvfrac
2297 radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling [W/kg]
2298 radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
2300 elseif( choice_radf .eq. 'maxi' ) then
2302 ! Radiative flux divergence both in 'kt' and 'kt-1' layers are included
2303 ! 1. From 'kt' layer
2304 lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
2305 opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
2306 radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
2307 radf = max( radinvfrac * qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 )
2308 ! 2. From 'kt-1' layer and add the contribution from 'kt' layer
2309 lwp = ql(i,kt-1) * ( pi(i,kt) - pi(i,kt-1) ) / g
2310 opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
2311 radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth) + opt_depth**2 )
2312 radf = radf + max( radinvfrac * qrlw(i,kt-1) / ( pi(i,kt-1) - pi(i,kt) ) * ( zi(i,kt-1) - zi(i,kt) ), &
2313 0._r8 )
2314 radf = max( radf, 0._r8 ) * chs(i,kt)
2316 endif
2318 ! ------------------------------------------------------------------- !
2319 ! Calculate 'wstar3' by summing buoyancy productions within CL from !
2320 ! 1. Interior buoyancy production ( bprod: fcn of TKE ) !
2321 ! 2. Cloud-top radiative cooling !
2322 ! 3. Surface buoyancy flux contribution only when bflxs > 0. !
2323 ! Note that master length scale, lbulk, has already been !
2324 ! corrctly defined at the first part of this 'do ncv' loop !
2325 ! considering the sign of bflxs. !
2326 ! This 'wstar3' is used for calculation of entrainment rate. !
2327 ! Note that this 'wstar3' formula does not include shear production !
2328 ! and the effect of drizzle, which should be included later. !
2329 ! Q : Strictly speaking, in calculating interior buoyancy production, !
2330 ! the use of 'bprod' is not correct, since 'bprod' is not correct !
2331 ! value but initially guessed value. More reasonably, we should !
2332 ! use '-leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)' instead !
2333 ! of 'bprod(i,k)', although this is still an approximation since !
2334 ! tke(i,k) is not exactly 'b1*wcap(i,k)' due to a transport term.!
2335 ! However since iterative calculation will be performed after all,!
2336 ! below might also be OK. But I should test this alternative. !
2337 ! ------------------------------------------------------------------- !
2339 dzht = zi(i,kt) - z(i,kt) ! Thickness of CL top half-layer
2340 dzhb = z(i,kb-1) - zi(i,kb) ! Thickness of CL bot half-layer
2341 wstar3 = radf * dzht
2342 do k = kt + 1, kb - 1 ! If 'kt = kb - 1', this loop will not be performed.
2343 wstar3 = wstar3 + bprod(i,k) * ( z(i,k-1) - z(i,k) )
2344 ! Below is an alternative which may speed up convergence.
2345 ! However, for interfaces merged into original CL, it can
2346 ! be 'wcap(i,k)<0' since 'n2(i,k)>0'. Thus, I should use
2347 ! the above original one.
2348 ! wstar3 = wstar3 - leng(i,k)*sqrt(b1*wcap(i,k))*shcl(i,ncv)*n2(i,k)* &
2349 ! (z(i,k-1) - z(i,k))
2350 end do
2351 if( kb .eq. (pver+1) .and. bflxs(i) .gt. 0._r8 ) then
2352 wstar3 = wstar3 + bflxs(i) * dzhb
2353 ! wstar3 = wstar3 + bprod(i,pver+1) * dzhb
2354 end if
2355 wstar3 = max( 2.5_r8 * wstar3, 0._r8 )
2357 ! -------------------------------------------------------------- !
2358 ! Below single block is for 'sedimentation-entrainment feedback' !
2359 ! -------------------------------------------------------------- !
2361 if( id_sedfact ) then
2362 ! wsed = 7.8e5_r8*(ql(i,kt)/ncliq(i,kt))**(2._r8/3._r8)
2363 sedfact = exp(-ased*wsedl(i,kt)/(wstar3**(1._r8/3._r8)+1.e-6))
2364 if( choice_evhc .eq. 'orig' ) then
2365 if (ql(i,kt).gt.qmin .and. ql(i,kt-1).lt.qmin) then
2366 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2367 jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
2368 evhc = 1._r8+sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
2369 evhc = min(evhc,evhcmax)
2370 end if
2371 elseif( choice_evhc .eq. 'ramp' ) then
2372 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2373 jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
2374 evhc = 1._r8+max(cldeff(i,kt)-cldeff(i,kt-1),0._r8)*sedfact*a2l*a3l*latvap*ql(i,kt) / jt2slv
2375 evhc = min(evhc,evhcmax)
2376 elseif( choice_evhc .eq. 'maxi' ) then
2377 qleff = max(ql(i,kt-1),ql(i,kt))
2378 jt2slv = slv(i,max(kt-2,1)) - slv(i,kt)
2379 jt2slv = max(jt2slv, jbumin*slv(i,kt-1)/g)
2380 evhc = 1._r8+sedfact*a2l*a3l*latvap*qleff / jt2slv
2381 evhc = min(evhc,evhcmax)
2382 endif
2383 endif
2385 ! -------------------------------------------------------------------------- !
2386 ! Now diagnose CL top and bottom entrainment rates (and the contribution of !
2387 ! top/bottom entrainments to wstar3) using entrainment closures of the form !
2388 ! !
2389 ! wet = cet*wstar3, web = ceb*wstar3 !
2390 ! !
2391 ! where cet and ceb depend on the entrainment interface jumps, ql, etc. !
2392 ! No entrainment is diagnosed unless the wstar3 > 0. Note '1/wstar3fact' is !
2393 ! a factor indicating the enhancement of wstar3 due to entrainment process. !
2394 ! Q : Below setting of 'wstar3fact = max(..,0.5)'might prevent the possible !
2395 ! case when buoyancy consumption by entrainment is stronger than cloud !
2396 ! top radiative cooling production. Is that OK ? No. According to bulk !
2397 ! modeling study, entrainment buoyancy consumption was always a certain !
2398 ! fraction of other net productions, rather than a separate sum. Thus, !
2399 ! below max limit of wstar3fact is correct. 'wstar3fact = max(.,0.5)' !
2400 ! prevents unreasonable enhancement of CL entrainment rate by cloud-top !
2401 ! entrainment instability, CTEI. !
2402 ! Q : Use of the same dry entrainment coefficient, 'a1i' both at the CL top !
2403 ! and base interfaces may result in too small 'wstar3' and 'ebrk' below, !
2404 ! as was seen in my generalized bulk modeling study. This should be re- !
2405 ! considered later !
2406 ! -------------------------------------------------------------------------- !
2408 if( wstar3 .gt. 0._r8 ) then
2409 cet = a1i * evhc / ( jtbu * lbulk )
2410 if( kb .eq. pver + 1 ) then
2411 wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht, wstar3factcrit )
2412 else
2413 ceb = a1i / ( jbbu * lbulk )
2414 wstar3fact = max( 1._r8 + 2.5_r8 * cet * n2ht * jtzm * dzht &
2415 + 2.5_r8 * ceb * n2hb * jbzm * dzhb, wstar3factcrit )
2416 end if
2417 wstar3 = wstar3 / wstar3fact
2418 else ! wstar3 == 0
2419 wstar3fact = 0._r8 ! This is just for dianostic output
2420 cet = 0._r8
2421 ceb = 0._r8
2422 end if
2424 ! ---------------------------------------------------------------------------- !
2425 ! Calculate net CL mean TKE including entrainment contribution by solving a !
2426 ! canonical cubic equation. The solution of cubic equ. is 'rootp**2 = ebrk' !
2427 ! where 'ebrk' originally (before solving cubic eq.) was interior CL mean TKE, !
2428 ! but after solving cubic equation, it is replaced by net CL mean TKE in the !
2429 ! same variable 'ebrk'. !
2430 ! ---------------------------------------------------------------------------- !
2431 ! Solve cubic equation (canonical form for analytic solution) !
2432 ! r^3 - 3*trmp*r - 2*trmq = 0, r = sqrt<e> !
2433 ! to estimate <e> for CL, derived from layer-mean TKE balance: !
2434 ! !
2435 ! <e>^(3/2)/(b_1*<l>) \approx <B + S> (*) !
2436 ! <B+S> = (<B+S>_int * l_int + <B+S>_et * dzt + <B+S>_eb * dzb)/lbulk !
2437 ! <B+S>_int = <e>^(1/2)/(b_1*<l>)*<e>_int !
2438 ! <B+S>_et = (-vyt+vut)*wet*jtbu + radf !
2439 ! <B+S>_eb = (-vyb+vub)*web*jbbu !
2440 ! !
2441 ! where: !
2442 ! <> denotes a vertical avg (over the whole CL unless indicated) !
2443 ! l_int (called lbrk below) is aggregate thickness of interior CL layers !
2444 ! dzt = zi(i,kt)-z(i,kt) is thickness of top entrainment layer !
2445 ! dzb = z(i,kb-1)-zi(i,kb) is thickness of bot entrainment layer !
2446 ! <e>_int (called ebrk below) is the CL-mean TKE if only interior !
2447 ! interfaces contributed. !
2448 ! wet, web are top. bottom entrainment rates !
2449 ! !
2450 ! For a single-level radiatively-driven convective layer, there are no !
2451 ! interior interfaces so 'ebrk' = 'lbrk' = 0. If the CL goes to the !
2452 ! surface, 'vyb' and 'vub' are set to zero before and 'ebrk' and 'lbrk' !
2453 ! have already incorporated the surface interfacial layer contribution, !
2454 ! so the same formulas still apply. !
2455 ! !
2456 ! In the original formulation based on TKE, !
2457 ! wet*jtbu = a1l*evhc*<e>^3/2/leng(i,kt) !
2458 ! web*jbbu = a1l*<e>^3/2/leng(i,kt) !
2459 ! !
2460 ! In the wstar formulation !
2461 ! wet*jtbu = a1i*evhc*wstar3/lbulk !
2462 ! web*jbbu = a1i*wstar3/lbulk, !
2463 ! ---------------------------------------------------------------------------- !
2465 fact = ( evhc * ( -vyt + vut ) * dzht + ( -vyb + vub ) * dzhb * leng(i,kb) / leng(i,kt) ) / lbulk
2467 if( wstarent ) then
2469 ! (Option 1) 'wstar' entrainment formulation
2470 ! Here trmq can have either sign, and will usually be nonzero even for non-
2471 ! cloud topped CLs. If trmq > 0, there will be two positive roots r; we take
2472 ! the larger one. Why ? If necessary, we limit entrainment and wstar to prevent
2473 ! a solution with r < ccrit*wstar ( Why ? ) where we take ccrit = 0.5.
2475 trma = 1._r8
2476 trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / 3._r8 + ntzero
2477 trmq = 0.5_r8 * b1 * ( leng(i,kt) / lbulk ) * ( radf * dzht + a1i * fact * wstar3 )
2479 ! Check if there is an acceptable root with r > rcrit = ccrit*wstar.
2480 ! To do this, first find local minimum fmin of the cubic f(r) at sqrt(p),
2481 ! and value fcrit = f(rcrit).
2483 rmin = sqrt(trmp)
2484 fmin = rmin * ( rmin * rmin - 3._r8 * trmp ) - 2._r8 * trmq
2485 wstar = wstar3**onet
2486 rcrit = ccrit * wstar
2487 fcrit = rcrit * ( rcrit * rcrit - 3._r8 * trmp ) - 2._r8 * trmq
2489 ! No acceptable root exists (noroot = .true.) if either:
2490 ! 1) rmin < rcrit (in which case cubic is monotone increasing for r > rcrit)
2491 ! and f(rcrit) > 0.
2492 ! or 2) rmin > rcrit (in which case min of f(r) in r > rcrit is at rmin)
2493 ! and f(rmin) > 0.
2494 ! In this case, we reduce entrainment and wstar3 such that r/wstar = ccrit;
2495 ! this changes the coefficients of the cubic. It might be informative to
2496 ! check when and how many 'noroot' cases occur, since when 'noroot', we
2497 ! will impose arbitrary limit on 'wstar3, wet, web, and ebrk' using ccrit.
2499 noroot = ( ( rmin .lt. rcrit ) .and. ( fcrit .gt. 0._r8 ) ) &
2500 .or. ( ( rmin .ge. rcrit ) .and. ( fmin .gt. 0._r8 ) )
2501 if( noroot ) then ! Solve cubic for r
2502 trma = 1._r8 - b1 * ( leng(i,kt) / lbulk ) * a1i * fact / ccrit**3
2503 trma = max( trma, 0.5_r8 ) ! Limit entrainment enhancement of ebrk
2504 trmp = trmp / trma
2505 trmq = 0.5_r8 * b1 * ( leng(i,kt) / lbulk ) * radf * dzht / trma
2506 end if ! noroot
2508 ! Solve the cubic equation
2510 qq = trmq**2 - trmp**3
2511 if( qq .ge. 0._r8 ) then
2512 rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
2513 else
2514 rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
2515 end if
2517 ! Adjust 'wstar3' only if there is 'noroot'.
2518 ! And calculate entrainment rates at the top and base interfaces.
2520 if( noroot ) wstar3 = ( rootp / ccrit )**3 ! Adjust wstar3
2521 wet = cet * wstar3 ! Find entrainment rates
2522 if( kb .lt. pver + 1 ) web = ceb * wstar3 ! When 'kb.eq.pver+1', it was set to web=0.
2524 else !
2526 ! (Option.2) wstarentr = .false. Use original entrainment formulation.
2527 ! trmp > 0 if there are interior interfaces in CL, trmp = 0 otherwise.
2528 ! trmq > 0 if there is cloudtop radiative cooling, trmq = 0 otherwise.
2530 trma = 1._r8 - b1 * a1l * fact
2531 trma = max( trma, 0.5_r8 ) ! Prevents runaway entrainment instability
2532 trmp = ebrk(i,ncv) * ( lbrk(i,ncv) / lbulk ) / ( 3._r8 * trma )
2533 trmq = 0.5_r8 * b1 * ( leng(i,kt) / lbulk ) * radf * dzht / trma
2535 qq = trmq**2 - trmp**3
2536 if( qq .ge. 0._r8 ) then
2537 rootp = ( trmq + sqrt(qq) )**(1._r8/3._r8) + ( max( trmq - sqrt(qq), 0._r8 ) )**(1._r8/3._r8)
2538 else ! Also part of case 3
2539 rootp = 2._r8 * sqrt(trmp) * cos( acos( trmq / sqrt(trmp**3) ) / 3._r8 )
2540 end if ! qq
2542 ! Find entrainment rates and limit them by free-entrainment values a1l*sqrt(e)
2544 wet = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kt) * jtbu ), 1._r8 )
2545 if( kb .lt. pver + 1 ) web = a1l * rootp * min( evhc * rootp**2 / ( leng(i,kb) * jbbu ), 1._r8 )
2547 end if ! wstarentr
2549 ! ---------------------------------------------------- !
2550 ! Finally, get the net CL mean TKE and normalized TKE !
2551 ! ---------------------------------------------------- !
2553 ebrk(i,ncv) = rootp**2
2554 ebrk(i,ncv) = min(ebrk(i,ncv),tkemax) ! Limit CL-avg TKE used for entrainment
2555 wbrk(i,ncv) = ebrk(i,ncv)/b1
2557 ! The only way ebrk = 0 is for SRCL which are actually radiatively cooled
2558 ! at top interface. In this case, we remove 'convective' label from the
2559 ! interfaces around this layer. This case should now be impossible, so
2560 ! we flag it. Q: I can't understand why this case is impossible now. Maybe,
2561 ! due to various limiting procedures used in solving cubic equation ?
2562 ! In case of SRCL, 'ebrk' should be positive due to cloud top LW radiative
2563 ! cooling contribution, although 'ebrk(internal)' of SRCL before including
2564 ! entrainment contribution (which include LW cooling contribution also) is
2565 ! zero.
2567 if( ebrk(i,ncv) .le. 0._r8 ) then
2568 write(iulog,*) 'CALEDDY: Warning, CL with zero TKE, i, kt, kb ', i, kt, kb
2569 #ifdef WRF_PORT
2570 call wrf_message(iulog)
2571 #endif
2572 belongcv(i,kt) = .false.
2573 belongcv(i,kb) = .false.
2574 end if
2576 ! ----------------------------------------------------------------------- !
2577 ! Calculate complete TKE profiles at all CL interfaces, capped by tkemax. !
2578 ! We approximate TKE = <e> at entrainment interfaces. However when CL is !
2579 ! based at surface, correct 'tkes' will be inserted to tke(i,pver+1). !
2580 ! Note that this approximation at CL external interfaces do not influence !
2581 ! numerical calculation since 'e' at external interfaces are not used in !
2582 ! actual numerical calculation afterward. In addition in order to extract !
2583 ! correct TKE averaged over the PBL in the cumulus scheme,it is necessary !
2584 ! to set e = <e> at the top entrainment interface. Since net CL mean TKE !
2585 ! 'ebrk' obtained by solving cubic equation already includes tkes ( tkes !
2586 ! is included when bflxs > 0 but not when bflxs <= 0 into internal ebrk ),!
2587 ! 'tkes' should be written to tke(i,pver+1) !
2588 ! ----------------------------------------------------------------------- !
2590 ! 1. At internal interfaces
2591 do k = kb - 1, kt + 1, -1
2592 rcap = ( b1 * ae + wcap(i,k) / wbrk(i,ncv) ) / ( b1 * ae + 1._r8 )
2593 rcap = min( max(rcap,rcapmin), rcapmax )
2594 tke(i,k) = ebrk(i,ncv) * rcap
2595 tke(i,k) = min( tke(i,k), tkemax )
2596 kvh(i,k) = leng(i,k) * sqrt(tke(i,k)) * shcl(i,ncv)
2597 kvm(i,k) = leng(i,k) * sqrt(tke(i,k)) * smcl(i,ncv)
2598 bprod(i,k) = -kvh(i,k) * n2(i,k)
2599 sprod(i,k) = kvm(i,k) * s2(i,k)
2600 turbtype(i,k) = 2 ! CL interior interfaces.
2601 sm_aw(i,k) = smcl(i,ncv)/alph1 ! Diagnostic output for microphysics
2602 end do
2604 ! 2. At CL top entrainment interface
2605 kentr = wet * jtzm
2606 kvh(i,kt) = kentr
2607 kvm(i,kt) = kentr
2608 bprod(i,kt) = -kentr * n2ht + radf ! I must use 'n2ht' not 'n2'
2609 sprod(i,kt) = kentr * s2(i,kt)
2610 turbtype(i,kt) = 4 ! CL top entrainment interface
2611 trmp = -b1 * ae / ( 1._r8 + b1 * ae )
2612 trmq = -(bprod(i,kt)+sprod(i,kt))*b1*leng(i,kt)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
2613 rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
2614 rcap = min( max(rcap,rcapmin), rcapmax )
2615 tke(i,kt) = ebrk(i,ncv) * rcap
2616 tke(i,kt) = min( tke(i,kt), tkemax )
2617 sm_aw(i,kt) = smcl(i,ncv) / alph1 ! Diagnostic output for microphysics
2619 ! 3. At CL base entrainment interface and double entraining interfaces
2620 ! When current CL base is also the top interface of CL regime below,
2621 ! simply add the two contributions for calculating eddy diffusivity
2622 ! and buoyancy/shear production. Below code correctly works because
2623 ! we (CL regime index) always go from surface upward.
2625 if( kb .lt. pver + 1 ) then
2627 kentr = web * jbzm
2629 if( kb .ne. ktblw ) then
2631 kvh(i,kb) = kentr
2632 kvm(i,kb) = kentr
2633 bprod(i,kb) = -kvh(i,kb)*n2hb ! I must use 'n2hb' not 'n2'
2634 sprod(i,kb) = kvm(i,kb)*s2(i,kb)
2635 turbtype(i,kb) = 3 ! CL base entrainment interface
2636 trmp = -b1*ae/(1._r8+b1*ae)
2637 trmq = -(bprod(i,kb)+sprod(i,kb))*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
2638 rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
2639 rcap = min( max(rcap,rcapmin), rcapmax )
2640 tke(i,kb) = ebrk(i,ncv) * rcap
2641 tke(i,kb) = min( tke(i,kb),tkemax )
2643 else
2645 kvh(i,kb) = kvh(i,kb) + kentr
2646 kvm(i,kb) = kvm(i,kb) + kentr
2647 ! dzhb5 : Half thickness of the lowest layer of current CL regime
2648 ! dzht5 : Half thickness of the highest layer of adjacent CL regime just below current CL.
2649 dzhb5 = z(i,kb-1) - zi(i,kb)
2650 dzht5 = zi(i,kb) - z(i,kb)
2651 bprod(i,kb) = ( dzht5*bprod(i,kb) - dzhb5*kentr*n2hb ) / ( dzhb5 + dzht5 )
2652 sprod(i,kb) = ( dzht5*sprod(i,kb) + dzhb5*kentr*s2(i,kb) ) / ( dzhb5 + dzht5 )
2653 trmp = -b1*ae/(1._r8+b1*ae)
2654 trmq = -kentr*(s2(i,kb)-n2hb)*b1*leng(i,kb)/(1._r8+b1*ae)/(ebrk(i,ncv)**(3._r8/2._r8))
2655 rcap = compute_cubic(0._r8,trmp,trmq)**2._r8
2656 rcap = min( max(rcap,rcapmin), rcapmax )
2657 tke_imsi = ebrk(i,ncv) * rcap
2658 tke_imsi = min( tke_imsi, tkemax )
2659 tke(i,kb) = ( dzht5*tke(i,kb) + dzhb5*tke_imsi ) / ( dzhb5 + dzht5 )
2660 tke(i,kb) = min(tke(i,kb),tkemax)
2661 turbtype(i,kb) = 5 ! CL double entraining interface
2663 end if
2665 else
2667 ! If CL base interface is surface, compute similarly using wcap(i,kb)=tkes/b1
2668 ! Even when bflx < 0, use the same formula in order to impose consistency of
2669 ! tke(i,kb) at bflx = 0._r8
2671 rcap = (b1*ae + wcap(i,kb)/wbrk(i,ncv))/(b1*ae + 1._r8)
2672 rcap = min( max(rcap,rcapmin), rcapmax )
2673 tke(i,kb) = ebrk(i,ncv) * rcap
2674 tke(i,kb) = min( tke(i,kb),tkemax )
2676 end if
2678 ! For double entraining interface, simply use smcl(i,ncv) of the overlying CL.
2679 ! Below 'sm_aw' is a diagnostic output for use in the microphysics.
2680 ! When 'kb' is surface, 'sm' will be over-written later below.
2682 sm_aw(i,kb) = smcl(i,ncv)/alph1
2684 ! Calculate wcap at all interfaces of CL. Put a minimum threshold on TKE
2685 ! to prevent possible division by zero. 'wcap' at CL internal interfaces
2686 ! are already calculated in the first part of 'do ncv' loop correctly.
2687 ! When 'kb.eq.pver+1', below formula produces the identical result to the
2688 ! 'tkes(i)/b1' if leng(i,kb) is set to vk*z(i,pver). Note wcap(i,pver+1)
2689 ! is already defined as 'tkes(i)/b1' at the first part of caleddy.
2691 wcap(i,kt) = (bprod(i,kt)+sprod(i,kt))*leng(i,kt)/sqrt(max(tke(i,kt),1.e-6_r8))
2692 if( kb .lt. pver + 1 ) then
2693 wcap(i,kb) = (bprod(i,kb)+sprod(i,kb))*leng(i,kb)/sqrt(max(tke(i,kb),1.e-6_r8))
2694 end if
2696 ! Save the index of upper external interface of current CL-regime in order to
2697 ! handle the case when this interface is also the lower external interface of
2698 ! CL-regime located just above.
2700 ktblw = kt
2702 ! Diagnostic Output
2704 wet_CL(i,ncv) = wet
2705 web_CL(i,ncv) = web
2706 jtbu_CL(i,ncv) = jtbu
2707 jbbu_CL(i,ncv) = jbbu
2708 evhc_CL(i,ncv) = evhc
2709 jt2slv_CL(i,ncv) = jt2slv
2710 n2ht_CL(i,ncv) = n2ht
2711 n2hb_CL(i,ncv) = n2hb
2712 lwp_CL(i,ncv) = lwp
2713 opt_depth_CL(i,ncv) = opt_depth
2714 radinvfrac_CL(i,ncv) = radinvfrac
2715 radf_CL(i,ncv) = radf
2716 wstar_CL(i,ncv) = wstar
2717 wstar3fact_CL(i,ncv) = wstar3fact
2719 end do ! ncv
2721 ! Calculate PBL height and characteristic cumulus excess for use in the
2722 ! cumulus convection shceme. Also define turbulence type at the surface
2723 ! when the lowest CL is based at the surface. These are just diagnostic
2724 ! outputs, not influencing numerical calculation of current PBL scheme.
2725 ! If the lowest CL is based at the surface, define the PBL depth as the
2726 ! CL top interface. The same rule is applied for all CLs including SRCL.
2728 if( ncvsurf .gt. 0 ) then
2730 ktopbl(i) = ktop(i,ncvsurf)
2731 pblh(i) = zi(i, ktopbl(i))
2732 pblhp(i) = pi(i, ktopbl(i))
2733 wpert(i) = max(wfac*sqrt(ebrk(i,ncvsurf)),wpertmin)
2734 tpert(i) = max(abs(shflx(i)*rrho(i)/cpair)*tfac/wpert(i),0._r8)
2735 qpert(i) = max(abs(qflx(i)*rrho(i))*tfac/wpert(i),0._r8)
2737 if( bflxs(i) .gt. 0._r8 ) then
2738 turbtype(i,pver+1) = 2 ! CL interior interface
2739 else
2740 turbtype(i,pver+1) = 3 ! CL external base interface
2741 endif
2743 ipbl(i) = 1._r8
2744 kpblh(i) = ktopbl(i) - 1._r8
2746 end if ! End of the calculationf of te properties of surface-based CL.
2748 ! -------------------------------------------- !
2749 ! Treatment of Stable Turbulent Regime ( STL ) !
2750 ! -------------------------------------------- !
2752 ! Identify top and bottom most (internal) interfaces of STL except surface.
2753 ! Also, calculate 'turbulent length scale (leng)' at each STL interfaces.
2755 belongst(i,1) = .false. ! k = 1 (top interface) is assumed non-turbulent
2756 do k = 2, pver ! k is an interface index
2757 belongst(i,k) = ( ri(i,k) .lt. ricrit ) .and. ( .not. belongcv(i,k) )
2758 if( belongst(i,k) .and. ( .not. belongst(i,k-1) ) ) then
2759 kt = k ! Top interface index of STL
2760 elseif( .not. belongst(i,k) .and. belongst(i,k-1) ) then
2761 kb = k - 1 ! Base interface index of STL
2762 lbulk = z(i,kt-1) - z(i,kb)
2763 do ks = kt, kb
2764 if( choice_tunl .eq. 'rampcl' ) then
2765 tunlramp = tunl
2766 elseif( choice_tunl .eq. 'rampsl' ) then
2767 tunlramp = max( 1.e-3_r8, ctunl * tunl * exp(-log(ctunl)*ri(i,ks)/ricrit) )
2768 ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
2769 else
2770 tunlramp = tunl
2771 endif
2772 if( choice_leng .eq. 'origin' ) then
2773 leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
2774 ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
2775 else
2776 leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )
2777 endif
2778 end do
2779 end if
2780 end do ! k
2782 ! Now look whether STL extends to ground. If STL extends to surface,
2783 ! re-define master length scale,'lbulk' including surface interfacial
2784 ! layer thickness, and re-calculate turbulent length scale, 'leng' at
2785 ! all STL interfaces again. Note that surface interface is assumed to
2786 ! always be STL if it is not CL.
2788 belongst(i,pver+1) = .not. belongcv(i,pver+1)
2790 if( belongst(i,pver+1) ) then ! kb = pver+1 (surface STL)
2792 turbtype(i,pver+1) = 1 ! Surface is STL interface
2794 if( belongst(i,pver) ) then ! STL includes interior
2795 ! 'kt' already defined above as the top interface of STL
2796 lbulk = z(i,kt-1)
2797 else ! STL with no interior turbulence
2798 kt = pver+1
2799 lbulk = z(i,kt-1)
2800 end if
2802 ! PBL height : Layer mid-point just above the highest STL interface
2803 ! Note in contrast to the surface based CL regime where PBL height
2804 ! was defined at the top external interface, PBL height of surface
2805 ! based STL is defined as the layer mid-point.
2807 ktopbl(i) = kt - 1
2808 pblh(i) = z(i,ktopbl(i))
2809 pblhp(i) = 0.5_r8 * ( pi(i,ktopbl(i)) + pi(i,ktopbl(i)+1) )
2811 ! Re-calculate turbulent length scale including surface interfacial
2812 ! layer contribution to lbulk.
2814 do ks = kt, pver
2815 if( choice_tunl .eq. 'rampcl' ) then
2816 tunlramp = tunl
2817 elseif( choice_tunl .eq. 'rampsl' ) then
2818 tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,ks)/ricrit))
2819 ! tunlramp = 0.065_r8 + 0.7_r8 * exp(-20._r8*ri(i,ks))
2820 else
2821 tunlramp = tunl
2822 endif
2823 if( choice_leng .eq. 'origin' ) then
2824 leng(i,ks) = ( (vk*zi(i,ks))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
2825 ! leng(i,ks) = vk*zi(i,ks) / (1._r8+vk*zi(i,ks)/(tunlramp*lbulk))
2826 else
2827 leng(i,ks) = min( vk*zi(i,ks), tunlramp*lbulk )
2828 endif
2829 end do ! ks
2831 ! Characteristic cumulus excess of surface-based STL.
2832 ! We may be able to use ustar for wpert.
2834 wpert(i) = 0._r8
2835 tpert(i) = max(shflx(i)*rrho(i)/cpair*fak/ustar(i),0._r8) ! CCM stable-layer forms
2836 qpert(i) = max(qflx(i)*rrho(i)*fak/ustar(i),0._r8)
2838 ipbl(i) = 0._r8
2839 kpblh(i) = ktopbl(i)
2841 end if
2843 ! Calculate stability functions and energetics at the STL interfaces
2844 ! except the surface. Note that tke(i,pver+1) and wcap(i,pver+1) are
2845 ! already calculated in the first part of 'caleddy', kvm(i,pver+1) &
2846 ! kvh(i,pver+1) were already initialized to be zero, bprod(i,pver+1)
2847 ! & sprod(i,pver+1) were direcly calculated from the bflxs and ustar.
2848 ! Note transport term is assumed to be negligible at STL interfaces.
2850 do k = 2, pver
2852 if( belongst(i,k) ) then
2854 turbtype(i,k) = 1 ! STL interfaces
2855 trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
2856 trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
2857 trmc = ri(i,k)
2858 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
2859 ! Sanity Check
2860 if( det .lt. 0._r8 ) then
2861 write(iulog,*) 'The det < 0. for the STL in UW eddy_diff'
2862 #ifdef WRF_PORT
2863 call wrf_message(iulog)
2864 #endif
2865 stop
2866 end if
2867 gh = (-trmb + sqrt(det))/(2._r8*trma)
2868 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
2869 ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
2870 gh = min(max(gh,ghmin),0.0233_r8)
2871 sh = max(0._r8,alph5/(1._r8+alph3*gh))
2872 sm = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
2874 tke(i,k) = b1*(leng(i,k)**2)*(-sh*n2(i,k)+sm*s2(i,k))
2875 tke(i,k) = min(tke(i,k),tkemax)
2876 wcap(i,k) = tke(i,k)/b1
2877 kvh(i,k) = leng(i,k) * sqrt(tke(i,k)) * sh
2878 kvm(i,k) = leng(i,k) * sqrt(tke(i,k)) * sm
2879 bprod(i,k) = -kvh(i,k) * n2(i,k)
2880 sprod(i,k) = kvm(i,k) * s2(i,k)
2882 sm_aw(i,k) = sm/alph1 ! This is diagnostic output for use in the microphysics
2884 end if
2886 end do ! k
2888 ! --------------------------------------------------- !
2889 ! End of treatment of Stable Turbulent Regime ( STL ) !
2890 ! --------------------------------------------------- !
2892 ! --------------------------------------------------------------- !
2893 ! Re-computation of eddy diffusivity at the entrainment interface !
2894 ! assuming that it is purely STL (0<Ri<0.19). Note even Ri>0.19, !
2895 ! turbulent can exist at the entrainment interface since 'Sh,Sm' !
2896 ! do not necessarily go to zero even when Ri>0.19. Since Ri can !
2897 ! be fairly larger than 0.19 at the entrainment interface, I !
2898 ! should set minimum value of 'tke' to be 0. in order to prevent !
2899 ! sqrt(tke) from being imaginary. !
2900 ! --------------------------------------------------------------- !
2902 ! goto 888
2904 do k = 2, pver
2906 if( ( turbtype(i,k) .eq. 3 ) .or. ( turbtype(i,k) .eq. 4 ) .or. &
2907 ( turbtype(i,k) .eq. 5 ) ) then
2909 trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
2910 trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
2911 trmc = ri(i,k)
2912 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
2913 gh = (-trmb + sqrt(det))/(2._r8*trma)
2914 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
2915 ! gh = min(max(gh,-3.5334_r8),0.0233_r8)
2916 gh = min(max(gh,ghmin),0.0233_r8)
2917 sh = max(0._r8,alph5/(1._r8+alph3*gh))
2918 sm = max(0._r8,(alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
2920 lbulk = z(i,k-1) - z(i,k)
2922 if( choice_tunl .eq. 'rampcl' ) then
2923 tunlramp = tunl
2924 elseif( choice_tunl .eq. 'rampsl' ) then
2925 tunlramp = max(1.e-3_r8,ctunl*tunl*exp(-log(ctunl)*ri(i,k)/ricrit))
2926 ! tunlramp = 0.065_r8 + 0.7_r8*exp(-20._r8*ri(i,k))
2927 else
2928 tunlramp = tunl
2929 endif
2930 if( choice_leng .eq. 'origin' ) then
2931 leng_imsi = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
2932 ! leng_imsi = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
2933 else
2934 leng_imsi = min( vk*zi(i,k), tunlramp*lbulk )
2935 endif
2937 tke_imsi = b1*(leng_imsi**2)*(-sh*n2(i,k)+sm*s2(i,k))
2938 tke_imsi = min(max(tke_imsi,0._r8),tkemax)
2939 kvh_imsi = leng_imsi * sqrt(tke_imsi) * sh
2940 kvm_imsi = leng_imsi * sqrt(tke_imsi) * sm
2942 if( kvh(i,k) .lt. kvh_imsi ) then
2943 kvh(i,k) = kvh_imsi
2944 kvm(i,k) = kvm_imsi
2945 leng(i,k) = leng_imsi
2946 tke(i,k) = tke_imsi
2947 wcap(i,k) = tke_imsi / b1
2948 bprod(i,k) = -kvh_imsi * n2(i,k)
2949 sprod(i,k) = kvm_imsi * s2(i,k)
2950 sm_aw(i,k) = sm/alph1 ! This is diagnostic output for use in the microphysics
2951 turbtype(i,k) = 1 ! This was added on Dec.10.2009 for use in microphysics.
2952 endif
2954 end if
2956 end do
2958 ! 888 continue
2960 ! ------------------------------------------------------------------ !
2961 ! End of recomputation of eddy diffusivity at entrainment interfaces !
2962 ! ------------------------------------------------------------------ !
2964 ! As an option, we can impose a certain minimum back-ground diffusivity.
2966 ! do k = 1, pver+1
2967 ! kvh(i,k) = max(0.01_r8,kvh(i,k))
2968 ! kvm(i,k) = max(0.01_r8,kvm(i,k))
2969 ! enddo
2971 ! --------------------------------------------------------------------- !
2972 ! Diagnostic Output !
2973 ! Just for diagnostic purpose, calculate stability functions at each !
2974 ! interface including surface. Instead of assuming neutral stability, !
2975 ! explicitly calculate stability functions using an reverse procedure !
2976 ! starting from tkes(i) similar to the case of SRCL and SBCL in zisocl. !
2977 ! Note that it is possible to calculate stability functions even when !
2978 ! bflxs < 0. Note that this inverse method allows us to define Ri even !
2979 ! at the surface. Note also tkes(i) and sprod(i,pver+1) are always !
2980 ! positive values by limiters (e.g., ustar_min = 0.01). !
2981 ! Dec.12.2006 : Also just for diagnostic output, re-set !
2982 ! 'bprod(i,pver+1)= bflxs(i)' here. Note that this setting does not !
2983 ! influence numerical calculation at all - it is just for diagnostic !
2984 ! output. !
2985 ! --------------------------------------------------------------------- !
2987 bprod(i,pver+1) = bflxs(i)
2989 gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
2990 if( abs(alph5-gg*alph3) .le. 1.e-7_r8 ) then
2991 ! gh = -0.28_r8
2992 if( bprod(i,pver+1) .gt. 0._r8 ) then
2993 gh = -3.5334_r8
2994 else
2995 gh = ghmin
2996 endif
2997 else
2998 gh = gg/(alph5-gg*alph3)
2999 end if
3001 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
3002 if( bprod(i,pver+1) .gt. 0._r8 ) then
3003 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3004 else
3005 gh = min(max(gh,ghmin),0.0233_r8)
3006 endif
3008 gh_a(i,pver+1) = gh
3009 sh_a(i,pver+1) = max(0._r8,alph5/(1._r8+alph3*gh))
3010 if( bprod(i,pver+1) .gt. 0._r8 ) then
3011 sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
3012 else
3013 sm_a(i,pver+1) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
3014 endif
3015 sm_aw(i,pver+1) = sm_a(i,pver+1)/alph1
3016 ri_a(i,pver+1) = -(sm_a(i,pver+1)/sh_a(i,pver+1))*(bprod(i,pver+1)/sprod(i,pver+1))
3018 do k = 1, pver
3019 if( ri(i,k) .lt. 0._r8 ) then
3020 trma = alph3*alph4*ri(i,k) + 2._r8*b1*(alph2-alph4*alph5*ri(i,k))
3021 trmb = (alph3+alph4)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
3022 trmc = ri(i,k)
3023 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
3024 gh = (-trmb + sqrt(det))/(2._r8*trma)
3025 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3026 gh_a(i,k) = gh
3027 sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
3028 sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh))
3029 ri_a(i,k) = ri(i,k)
3030 else
3031 if( ri(i,k) .gt. ricrit ) then
3032 gh_a(i,k) = ghmin
3033 sh_a(i,k) = 0._r8
3034 sm_a(i,k) = 0._r8
3035 ri_a(i,k) = ri(i,k)
3036 else
3037 trma = alph3*alph4exs*ri(i,k) + 2._r8*b1*(alph2-alph4exs*alph5*ri(i,k))
3038 trmb = (alph3+alph4exs)*ri(i,k) + 2._r8*b1*(-alph5*ri(i,k)+alph1)
3039 trmc = ri(i,k)
3040 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
3041 gh = (-trmb + sqrt(det))/(2._r8*trma)
3042 gh = min(max(gh,ghmin),0.0233_r8)
3043 gh_a(i,k) = gh
3044 sh_a(i,k) = max(0._r8,alph5/(1._r8+alph3*gh))
3045 sm_a(i,k) = max(0._r8,(alph1+alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4exs*gh))
3046 ri_a(i,k) = ri(i,k)
3047 endif
3048 endif
3050 end do
3052 do k = 1, pver + 1
3053 turbtype_f(i,k) = real(turbtype(i,k))
3054 end do
3056 end do ! End of column index loop, i
3058 return
3060 end subroutine caleddy
3062 !============================================================================== !
3063 ! !
3064 !============================================================================== !
3066 subroutine exacol( pcols, pver, ncol, ri, bflxs, minpblh, zi, ktop, kbase, ncvfin )
3068 ! ---------------------------------------------------------------------------- !
3069 ! Object : Find unstable CL regimes and determine the indices !
3070 ! kbase, ktop which delimit these unstable layers : !
3071 ! ri(kbase) > 0 and ri(ktop) > 0, but ri(k) < 0 for ktop < k < kbase. !
3072 ! Author : Chris Bretherton 08/2000, !
3073 ! Sungsu Park 08/2006, 11/2008 !
3074 !----------------------------------------------------------------------------- !
3076 implicit none
3078 ! --------------- !
3079 ! Input variables !
3080 ! --------------- !
3082 integer, intent(in) :: pcols ! Number of atmospheric columns
3083 integer, intent(in) :: pver ! Number of atmospheric vertical layers
3084 integer, intent(in) :: ncol ! Number of atmospheric columns
3086 real(r8), intent(in) :: ri(pcols,pver) ! Moist gradient Richardson no.
3087 real(r8), intent(in) :: bflxs(pcols) ! Buoyancy flux at surface
3088 real(r8), intent(in) :: minpblh(pcols) ! Minimum PBL height based on surface stress
3089 real(r8), intent(in) :: zi(pcols,pver+1) ! Interface heights
3091 ! ---------------- !
3092 ! Output variables !
3093 ! ---------------- !
3095 integer, intent(out) :: kbase(pcols,ncvmax) ! External interface index of CL base
3096 integer, intent(out) :: ktop(pcols,ncvmax) ! External interface index of CL top
3097 integer, intent(out) :: ncvfin(pcols) ! Total number of CLs
3099 ! --------------- !
3100 ! Local variables !
3101 ! --------------- !
3103 integer :: i
3104 integer :: k
3105 integer :: ncv
3106 real(r8) :: rimaxentr
3107 real(r8) :: riex(pver+1) ! Column Ri profile extended to surface
3109 ! ----------------------- !
3110 ! Main Computation Begins !
3111 ! ----------------------- !
3113 do i = 1, ncol
3114 ncvfin(i) = 0
3115 do ncv = 1, ncvmax
3116 ktop(i,ncv) = 0
3117 kbase(i,ncv) = 0
3118 end do
3119 end do
3121 ! ------------------------------------------------------ !
3122 ! Find CL regimes starting from the surface going upward !
3123 ! ------------------------------------------------------ !
3125 rimaxentr = 0._r8
3127 do i = 1, ncol
3129 riex(2:pver) = ri(i,2:pver)
3131 ! Below allows consistent treatment of surface and other interfaces.
3132 ! Simply, if surface buoyancy flux is positive, Ri of surface is set to be negative.
3134 riex(pver+1) = rimaxentr - bflxs(i)
3136 ncv = 0
3137 k = pver + 1 ! Work upward from surface interface
3139 do while ( k .gt. ntop_turb + 1 )
3141 ! Below means that if 'bflxs > 0' (do not contain '=' sign), surface
3142 ! interface is energetically interior surface.
3144 if( riex(k) .lt. rimaxentr ) then
3146 ! Identify a new CL
3148 ncv = ncv + 1
3150 ! First define 'kbase' as the first interface below the lower-most unstable interface
3151 ! Thus, Richardson number at 'kbase' is positive.
3153 kbase(i,ncv) = min(k+1,pver+1)
3155 ! Decrement k until top unstable level
3157 do while( riex(k) .lt. rimaxentr .and. k .gt. ntop_turb + 1 )
3158 k = k - 1
3159 end do
3161 ! ktop is the first interface above upper-most unstable interface
3162 ! Thus, Richardson number at 'ktop' is positive.
3164 ktop(i,ncv) = k
3166 else
3168 ! Search upward for a CL.
3170 k = k - 1
3172 end if
3174 end do ! End of CL regime finding for each atmospheric column
3176 ncvfin(i) = ncv
3178 end do ! End of atmospheric column do loop
3180 return
3182 end subroutine exacol
3184 !============================================================================== !
3185 ! !
3186 !============================================================================== !
3188 subroutine zisocl( pcols , pver , long , &
3189 z , zi , n2 , s2 , &
3190 bprod , sprod , bflxs, tkes , &
3191 ncvfin , kbase , ktop , belongcv, &
3192 ricl , ghcl , shcl , smcl , &
3193 lbrk , wbrk , ebrk , extend , extend_up, extend_dn )
3195 !------------------------------------------------------------------------ !
3196 ! Object : This 'zisocl' vertically extends original CLs identified from !
3197 ! 'exacol' using a merging test based on either 'wint' or 'l2n2' !
3198 ! and identify new CL regimes. Similar to the case of 'exacol', !
3199 ! CL regime index increases with height. After identifying new !
3200 ! CL regimes ( kbase, ktop, ncvfin ),calculate CL internal mean !
3201 ! energetics (lbrk : energetic thickness integral, wbrk, ebrk ) !
3202 ! and stability functions (ricl, ghcl, shcl, smcl) by including !
3203 ! surface interfacial layer contribution when bflxs > 0. Note !
3204 ! that there are two options in the treatment of the energetics !
3205 ! of surface interfacial layer (use_dw_surf= 'true' or 'false') !
3206 ! Author : Sungsu Park 08/2006, 11/2008 !
3207 !------------------------------------------------------------------------ !
3209 implicit none
3211 ! --------------- !
3212 ! Input variables !
3213 ! --------------- !
3215 integer, intent(in) :: long ! Longitude of the column
3216 integer, intent(in) :: pcols ! Number of atmospheric columns
3217 integer, intent(in) :: pver ! Number of atmospheric vertical layers
3218 real(r8), intent(in) :: z(pcols, pver) ! Layer mid-point height [ m ]
3219 real(r8), intent(in) :: zi(pcols, pver+1) ! Interface height [ m ]
3220 real(r8), intent(in) :: n2(pcols, pver) ! Buoyancy frequency at interfaces except surface [ s-2 ]
3221 real(r8), intent(in) :: s2(pcols, pver) ! Shear frequency at interfaces except surface [ s-2 ]
3222 real(r8), intent(in) :: bprod(pcols,pver+1) ! Buoyancy production [ m2/s3 ]. bprod(i,pver+1) = bflxs
3223 real(r8), intent(in) :: sprod(pcols,pver+1) ! Shear production [ m2/s3 ]. sprod(i,pver+1) = usta**3/(vk*z(i,pver))
3224 real(r8), intent(in) :: bflxs(pcols) ! Surface buoyancy flux [ m2/s3 ]. bprod(i,pver+1) = bflxs
3225 real(r8), intent(in) :: tkes(pcols) ! TKE at the surface [ s2/s2 ]
3227 ! ---------------------- !
3228 ! Input/output variables !
3229 ! ---------------------- !
3231 integer, intent(inout) :: kbase(pcols,ncvmax) ! Base external interface index of CL
3232 integer, intent(inout) :: ktop(pcols,ncvmax) ! Top external interface index of CL
3233 integer, intent(inout) :: ncvfin(pcols) ! Total number of CLs
3235 ! ---------------- !
3236 ! Output variables !
3237 ! ---------------- !
3239 logical, intent(out) :: belongcv(pcols,pver+1) ! True if interface is in a CL ( either internal or external )
3240 real(r8), intent(out) :: ricl(pcols,ncvmax) ! Mean Richardson number of internal CL
3241 real(r8), intent(out) :: ghcl(pcols,ncvmax) ! Half of normalized buoyancy production of internal CL
3242 real(r8), intent(out) :: shcl(pcols,ncvmax) ! Galperin instability function of heat-moisture of internal CL
3243 real(r8), intent(out) :: smcl(pcols,ncvmax) ! Galperin instability function of momentum of internal CL
3244 real(r8), intent(out) :: lbrk(pcols,ncvmax) ! Thickness of (energetically) internal CL ( lint, [m] )
3245 real(r8), intent(out) :: wbrk(pcols,ncvmax) ! Mean normalized TKE of internal CL [ m2/s2 ]
3246 real(r8), intent(out) :: ebrk(pcols,ncvmax) ! Mean TKE of internal CL ( b1*wbrk, [m2/s2] )
3248 ! ------------------ !
3249 ! Internal variables !
3250 ! ------------------ !
3252 logical :: extend ! True when CL is extended in zisocl
3253 logical :: extend_up ! True when CL is extended upward in zisocl
3254 logical :: extend_dn ! True when CL is extended downward in zisocl
3255 logical :: bottom ! True when CL base is at surface ( kb = pver + 1 )
3257 integer :: i ! Local index for the longitude
3258 integer :: ncv ! CL Index increasing with height
3259 integer :: incv
3260 integer :: k
3261 integer :: kb ! Local index for kbase
3262 integer :: kt ! Local index for ktop
3263 integer :: ncvinit ! Value of ncv at routine entrance
3264 integer :: cntu ! Number of merged CLs during upward extension of individual CL
3265 integer :: cntd ! Number of merged CLs during downward extension of individual CL
3266 integer :: kbinc ! Index for incorporating underlying CL
3267 integer :: ktinc ! Index for incorporating overlying CL
3269 real(r8) :: wint ! Normalized TKE of internal CL
3270 real(r8) :: dwinc ! Normalized TKE of CL external interfaces
3271 real(r8) :: dw_surf ! Normalized TKE of surface interfacial layer
3272 real(r8) :: dzinc
3273 real(r8) :: gh
3274 real(r8) :: sh
3275 real(r8) :: sm
3276 real(r8) :: gh_surf ! Half of normalized buoyancy production in surface interfacial layer
3277 real(r8) :: sh_surf ! Galperin instability function in surface interfacial layer
3278 real(r8) :: sm_surf ! Galperin instability function in surface interfacial layer
3279 real(r8) :: l2n2 ! Vertical integral of 'l^2N^2' over CL. Include thickness product
3280 real(r8) :: l2s2 ! Vertical integral of 'l^2S^2' over CL. Include thickness product
3281 real(r8) :: dl2n2 ! Vertical integration of 'l^2*N^2' of CL external interfaces
3282 real(r8) :: dl2s2 ! Vertical integration of 'l^2*S^2' of CL external interfaces
3283 real(r8) :: dl2n2_surf ! 'dl2n2' defined in the surface interfacial layer
3284 real(r8) :: dl2s2_surf ! 'dl2s2' defined in the surface interfacial layer
3285 real(r8) :: lint ! Thickness of (energetically) internal CL
3286 real(r8) :: dlint ! Interfacial layer thickness of CL external interfaces
3287 real(r8) :: dlint_surf ! Surface interfacial layer thickness
3288 real(r8) :: lbulk ! Master Length Scale : Whole CL thickness from top to base external interface
3289 real(r8) :: lz ! Turbulent length scale
3290 real(r8) :: ricll ! Mean Richardson number of internal CL
3291 real(r8) :: trma
3292 real(r8) :: trmb
3293 real(r8) :: trmc
3294 real(r8) :: det
3295 real(r8) :: zbot ! Height of CL base
3296 real(r8) :: l2rat ! Square of ratio of actual to initial CL (not used)
3297 real(r8) :: gg ! Intermediate variable used for calculating stability functions of SBCL
3298 real(r8) :: tunlramp ! Ramping tunl
3300 ! ----------------------- !
3301 ! Main Computation Begins !
3302 ! ----------------------- !
3304 i = long
3306 ! Initialize main output variables
3308 do k = 1, ncvmax
3309 ricl(i,k) = 0._r8
3310 ghcl(i,k) = 0._r8
3311 shcl(i,k) = 0._r8
3312 smcl(i,k) = 0._r8
3313 lbrk(i,k) = 0._r8
3314 wbrk(i,k) = 0._r8
3315 ebrk(i,k) = 0._r8
3316 end do
3317 extend = .false.
3318 extend_up = .false.
3319 extend_dn = .false.
3321 ! ----------------------------------------------------------- !
3322 ! Loop over each CL to see if any of them need to be extended !
3323 ! ----------------------------------------------------------- !
3325 ncv = 1
3327 do while( ncv .le. ncvfin(i) )
3329 ncvinit = ncv
3330 cntu = 0
3331 cntd = 0
3332 kb = kbase(i,ncv)
3333 kt = ktop(i,ncv)
3335 ! ---------------------------------------------------------------------------- !
3336 ! Calculation of CL interior energetics including surface before extension !
3337 ! ---------------------------------------------------------------------------- !
3338 ! Note that the contribution of interior interfaces (not surface) to 'wint' is !
3339 ! accounted by using '-sh*l2n2 + sm*l2s2' while the contribution of surface is !
3340 ! accounted by using 'dwsurf = tkes/b1' when bflxs > 0. This approach is fully !
3341 ! reasonable. Another possible alternative, which seems to be also consistent !
3342 ! is to calculate 'dl2n2_surf' and 'dl2s2_surf' of surface interfacial layer !
3343 ! separately, and this contribution is explicitly added by initializing 'l2n2' !
3344 ! 'l2s2' not by zero, but by 'dl2n2_surf' and 'ds2n2_surf' below. At the same !
3345 ! time, 'dwsurf' should be excluded in 'wint' calculation below. The only diff.!
3346 ! between two approaches is that in case of the latter approach, contributions !
3347 ! of surface interfacial layer to the CL mean stability function (ri,gh,sh,sm) !
3348 ! are explicitly included while the first approach is not. In this sense, the !
3349 ! second approach seems to be more conceptually consistent, but currently, I !
3350 ! (Sungsu) will keep the first default approach. There is a switch !
3351 ! 'use_dw_surf' at the first part of eddy_diff.F90 chosing one of !
3352 ! these two options. !
3353 ! ---------------------------------------------------------------------------- !
3355 ! ------------------------------------------------------ !
3356 ! Step 0: Calculate surface interfacial layer energetics !
3357 ! ------------------------------------------------------ !
3359 lbulk = zi(i,kt) - zi(i,kb)
3360 dlint_surf = 0._r8
3361 dl2n2_surf = 0._r8
3362 dl2s2_surf = 0._r8
3363 dw_surf = 0._r8
3364 if( kb .eq. pver+1 ) then
3366 if( bflxs(i) .gt. 0._r8 ) then
3368 ! Calculate stability functions of surface interfacial layer
3369 ! from the given 'bprod(i,pver+1)' and 'sprod(i,pver+1)' using
3370 ! inverse approach. Since alph5>0 and alph3<0, denominator of
3371 ! gg is always positive if bprod(i,pver+1)>0.
3373 gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
3374 gh = gg/(alph5-gg*alph3)
3375 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
3376 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3377 sh = alph5/(1._r8+alph3*gh)
3378 sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
3379 ricll = min(-(sm/sh)*(bprod(i,pver+1)/sprod(i,pver+1)),ricrit)
3381 ! Calculate surface interfacial layer contribution to CL internal
3382 ! energetics. By construction, 'dw_surf = -dl2n2_surf + ds2n2_surf'
3383 ! is exactly satisfied, which corresponds to assuming turbulent
3384 ! length scale of surface interfacial layer = vk * z(i,pver). Note
3385 ! 'dl2n2_surf','dl2s2_surf','dw_surf' include thickness product.
3387 dlint_surf = z(i,pver)
3388 dl2n2_surf = -vk*(z(i,pver)**2)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
3389 dl2s2_surf = vk*(z(i,pver)**2)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
3390 dw_surf = (tkes(i)/b1)*z(i,pver)
3392 else
3394 ! Note that this case can happen when surface is an external
3395 ! interface of CL.
3396 lbulk = zi(i,kt) - z(i,pver)
3398 end if
3400 end if
3402 ! ------------------------------------------------------ !
3403 ! Step 1: Include surface interfacial layer contribution !
3404 ! ------------------------------------------------------ !
3406 lint = dlint_surf
3407 l2n2 = dl2n2_surf
3408 l2s2 = dl2s2_surf
3409 wint = dw_surf
3410 if( use_dw_surf ) then
3411 l2n2 = 0._r8
3412 l2s2 = 0._r8
3413 else
3414 wint = 0._r8
3415 end if
3417 ! --------------------------------------------------------------------------------- !
3418 ! Step 2. Include the contribution of 'pure internal interfaces' other than surface !
3419 ! --------------------------------------------------------------------------------- !
3421 if( kt .lt. kb - 1 ) then ! The case of non-SBCL.
3423 do k = kb - 1, kt + 1, -1
3424 if( choice_tunl .eq. 'rampcl' ) then
3425 ! Modification : I simply used the average tunlramp between the two limits.
3426 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3427 elseif( choice_tunl .eq. 'rampsl' ) then
3428 tunlramp = ctunl*tunl
3429 ! tunlramp = 0.765_r8
3430 else
3431 tunlramp = tunl
3432 endif
3433 if( choice_leng .eq. 'origin' ) then
3434 lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3435 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
3436 else
3437 lz = min( vk*zi(i,k), tunlramp*lbulk )
3438 endif
3439 dzinc = z(i,k-1) - z(i,k)
3440 l2n2 = l2n2 + lz*lz*n2(i,k)*dzinc
3441 l2s2 = l2s2 + lz*lz*s2(i,k)*dzinc
3442 lint = lint + dzinc
3443 end do
3445 ! Calculate initial CL stability functions (gh,sh,sm) and net
3446 ! internal energy of CL including surface contribution if any.
3448 ! Modification : It seems that below cannot be applied when ricrit > 0.19.
3449 ! May need future generalization.
3451 ricll = min(l2n2/max(l2s2,ntzero),ricrit) ! Mean Ri of internal CL
3452 trma = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
3453 trmb = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
3454 trmc = ricll
3455 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
3456 gh = (-trmb + sqrt(det))/2._r8/trma
3457 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
3458 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3459 sh = alph5/(1._r8+alph3*gh)
3460 sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
3461 wint = wint - sh*l2n2 + sm*l2s2
3463 else ! The case of SBCL
3465 ! If there is no pure internal interface, use only surface interfacial
3466 ! values. However, re-set surface interfacial values such that it can
3467 ! be used in the merging tests (either based on 'wint' or 'l2n2') and
3468 ! in such that surface interfacial energy is not double-counted.
3469 ! Note that regardless of the choise of 'use_dw_surf', below should be
3470 ! kept as it is below, for consistent merging test of extending SBCL.
3472 lint = dlint_surf
3473 l2n2 = dl2n2_surf
3474 l2s2 = dl2s2_surf
3475 wint = dw_surf
3477 ! Aug.29.2006 : Only for the purpose of merging test of extending SRCL
3478 ! based on 'l2n2', re-define 'l2n2' of surface interfacial layer using
3479 ! 'wint'. This part is designed for similar treatment of merging as in
3480 ! the original 'eddy_diff.F90' code, where 'l2n2' of SBCL was defined
3481 ! as 'l2n2 = - wint / sh'. Note that below block is used only when (1)
3482 ! surface buoyancy production 'bprod(i,pver+1)' is NOT included in the
3483 ! calculation of surface TKE in the initialization of 'bprod(i,pver+1)'
3484 ! in the main subroutine ( even though bflxs > 0 ), and (2) to force
3485 ! current scheme be similar to the previous scheme in the treatment of
3486 ! extending-merging test of SBCL based on 'l2n2'. Otherwise below line
3487 ! must be commented out. Note at this stage, correct non-zero value of
3488 ! 'sh' has been already computed.
3490 if( choice_tkes .eq. 'ebprod' ) then
3491 l2n2 = - wint / sh
3492 endif
3494 endif
3496 ! Set consistent upper limits on 'l2n2' and 'l2s2'. Below limits are
3497 ! reasonable since l2n2 of CL interior interface is always negative.
3499 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3500 l2s2 = min( l2s2, tkemax*lint/(b1*sm))
3502 ! Note that at this stage, ( gh, sh, sm ) are the values of surface
3503 ! interfacial layer if there is no pure internal interface, while if
3504 ! there is pure internal interface, ( gh, sh, sm ) are the values of
3505 ! pure CL interfaces or the values that include both the CL internal
3506 ! interfaces and surface interfaces, depending on the 'use_dw_surf'.
3508 ! ----------------------------------------------------------------------- !
3509 ! Perform vertical extension-merging process !
3510 ! ----------------------------------------------------------------------- !
3511 ! During the merging process, we assumed ( lbulk, sh, sm ) of CL external !
3512 ! interfaces are the same as the ones of the original merging CL. This is !
3513 ! an inevitable approximation since we don't know ( sh, sm ) of external !
3514 ! interfaces at this stage. Note that current default merging test is !
3515 ! purely based on buoyancy production without including shear production, !
3516 ! since we used 'l2n2' instead of 'wint' as a merging parameter. However, !
3517 ! merging test based on 'wint' maybe conceptually more attractable. !
3518 ! Downward CL merging process is identical to the upward merging process, !
3519 ! but when the base of extended CL reaches to the surface, surface inter !
3520 ! facial layer contribution to the energetic of extended CL must be done !
3521 ! carefully depending on the sign of surface buoyancy flux. The contribu !
3522 ! tion of surface interfacial layer energetic is included to the internal !
3523 ! energetics of merging CL only when bflxs > 0. !
3524 ! ----------------------------------------------------------------------- !
3526 ! ---------------------------- !
3527 ! Step 1. Extend the CL upward !
3528 ! ---------------------------- !
3530 extend = .false. ! This will become .true. if CL top or base is extended
3532 ! Calculate contribution of potentially incorporable CL top interface
3534 if( choice_tunl .eq. 'rampcl' ) then
3535 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3536 elseif( choice_tunl .eq. 'rampsl' ) then
3537 tunlramp = ctunl*tunl
3538 ! tunlramp = 0.765_r8
3539 else
3540 tunlramp = tunl
3541 endif
3542 if( choice_leng .eq. 'origin' ) then
3543 lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3544 ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
3545 else
3546 lz = min( vk*zi(i,kt), tunlramp*lbulk )
3547 endif
3549 dzinc = z(i,kt-1)-z(i,kt)
3550 dl2n2 = lz*lz*n2(i,kt)*dzinc
3551 dl2s2 = lz*lz*s2(i,kt)*dzinc
3552 dwinc = -sh*dl2n2 + sm*dl2s2
3554 ! ------------ !
3555 ! Merging Test !
3556 ! ------------ !
3558 ! do while ( dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) ) ! Merging test based on wint
3559 ! do while ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) ) ! Merging test based on l2n2
3560 #ifndef WRF_PORT
3561 do while ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) ) ! Integral merging test
3562 #else
3563 do while ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) .and. kt-1 .gt. ntop_turb ) ! Integral merging test
3564 #endif
3566 ! Add contribution of top external interface to interior energy.
3567 ! Note even when we chose 'use_dw_surf='true.', the contribution
3568 ! of surface interfacial layer to 'l2n2' and 'l2s2' are included
3569 ! here. However it is not double counting of surface interfacial
3570 ! energy : surface interfacial layer energy is counted in 'wint'
3571 ! formula and 'l2n2' is just used for performing merging test in
3572 ! this 'do while' loop.
3574 lint = lint + dzinc
3575 l2n2 = l2n2 + dl2n2
3576 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3577 l2s2 = l2s2 + dl2s2
3578 wint = wint + dwinc
3580 ! Extend top external interface of CL upward after merging
3582 kt = kt - 1
3583 extend = .true.
3584 extend_up = .true.
3585 if( kt .eq. ntop_turb ) then
3586 write(iulog,*) 'zisocl: Error: Tried to extend CL to the model top'
3587 #ifdef WRF_PORT
3588 call wrf_message(iulog)
3589 #endif
3590 stop
3591 end if
3593 ! If the top external interface of extending CL is the same as the
3594 ! top interior interface of the overlying CL, overlying CL will be
3595 ! automatically merged. Then,reduce total number of CL regime by 1.
3596 ! and increase 'cntu'(number of merged CLs during upward extension)
3597 ! by 1.
3599 ktinc = kbase(i,ncv+cntu+1) - 1 ! Lowest interior interface of overlying CL
3601 if( kt .eq. ktinc ) then
3603 do k = kbase(i,ncv+cntu+1) - 1, ktop(i,ncv+cntu+1) + 1, -1
3605 if( choice_tunl .eq. 'rampcl' ) then
3606 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3607 elseif( choice_tunl .eq. 'rampsl' ) then
3608 tunlramp = ctunl*tunl
3609 ! tunlramp = 0.765_r8
3610 else
3611 tunlramp = tunl
3612 endif
3613 if( choice_leng .eq. 'origin' ) then
3614 lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3615 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
3616 else
3617 lz = min( vk*zi(i,k), tunlramp*lbulk )
3618 endif
3620 dzinc = z(i,k-1)-z(i,k)
3621 dl2n2 = lz*lz*n2(i,k)*dzinc
3622 dl2s2 = lz*lz*s2(i,k)*dzinc
3623 dwinc = -sh*dl2n2 + sm*dl2s2
3625 lint = lint + dzinc
3626 l2n2 = l2n2 + dl2n2
3627 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3628 l2s2 = l2s2 + dl2s2
3629 wint = wint + dwinc
3631 end do
3633 kt = ktop(i,ncv+cntu+1)
3634 ncvfin(i) = ncvfin(i) - 1
3635 cntu = cntu + 1
3637 end if
3639 ! Again, calculate the contribution of potentially incorporatable CL
3640 ! top external interface of CL regime.
3642 if( choice_tunl .eq. 'rampcl' ) then
3643 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3644 elseif( choice_tunl .eq. 'rampsl' ) then
3645 tunlramp = ctunl*tunl
3646 ! tunlramp = 0.765_r8
3647 else
3648 tunlramp = tunl
3649 endif
3650 if( choice_leng .eq. 'origin' ) then
3651 lz = ( (vk*zi(i,kt))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3652 ! lz = vk*zi(i,kt) / (1._r8+vk*zi(i,kt)/(tunlramp*lbulk))
3653 else
3654 lz = min( vk*zi(i,kt), tunlramp*lbulk )
3655 endif
3657 dzinc = z(i,kt-1)-z(i,kt)
3658 dl2n2 = lz*lz*n2(i,kt)*dzinc
3659 dl2s2 = lz*lz*s2(i,kt)*dzinc
3660 dwinc = -sh*dl2n2 + sm*dl2s2
3662 end do ! End of upward merging test 'do while' loop
3664 ! Update CL interface indices appropriately if any CL was merged.
3665 ! Note that below only updated the interface index of merged CL,
3666 ! not the original merging CL. Updates of 'kbase' and 'ktop' of
3667 ! the original merging CL will be done after finishing downward
3668 ! extension also later.
3670 if( cntu .gt. 0 ) then
3671 do incv = 1, ncvfin(i) - ncv
3672 kbase(i,ncv+incv) = kbase(i,ncv+cntu+incv)
3673 ktop(i,ncv+incv) = ktop(i,ncv+cntu+incv)
3674 end do
3675 end if
3677 ! ------------------------------ !
3678 ! Step 2. Extend the CL downward !
3679 ! ------------------------------ !
3681 if( kb .ne. pver + 1 ) then
3683 ! Calculate contribution of potentially incorporable CL base interface
3685 if( choice_tunl .eq. 'rampcl' ) then
3686 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3687 elseif( choice_tunl .eq. 'rampsl' ) then
3688 tunlramp = ctunl*tunl
3689 ! tunlramp = 0.765_r8
3690 else
3691 tunlramp = tunl
3692 endif
3693 if( choice_leng .eq. 'origin' ) then
3694 lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3695 ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
3696 else
3697 lz = min( vk*zi(i,kb), tunlramp*lbulk )
3698 endif
3700 dzinc = z(i,kb-1)-z(i,kb)
3701 dl2n2 = lz*lz*n2(i,kb)*dzinc
3702 dl2s2 = lz*lz*s2(i,kb)*dzinc
3703 dwinc = -sh*dl2n2 + sm*dl2s2
3705 ! ------------ !
3706 ! Merging test !
3707 ! ------------ !
3709 ! In the below merging tests, I must keep '.and.(kb.ne.pver+1)',
3710 ! since 'kb' is continuously updated within the 'do while' loop
3711 ! whenever CL base is merged.
3713 ! do while( ( dwinc .gt. ( rinc*dzinc*wint/(lint+(1._r8-rinc)*dzinc)) ) & ! Merging test based on wint
3714 ! do while( ( -dl2n2 .gt. (-rinc*dzinc*l2n2/(lint+(1._r8-rinc)*dzinc)) ) & ! Merging test based on l2n2
3715 ! .and.(kb.ne.pver+1))
3716 do while( ( -dl2n2 .gt. (-rinc*l2n2/(1._r8-rinc)) ) & ! Integral merging test
3717 .and.(kb.ne.pver+1))
3719 ! Add contributions from interfacial layer kb to CL interior
3721 lint = lint + dzinc
3722 l2n2 = l2n2 + dl2n2
3723 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3724 l2s2 = l2s2 + dl2s2
3725 wint = wint + dwinc
3727 ! Extend the base external interface of CL downward after merging
3729 kb = kb + 1
3730 extend = .true.
3731 extend_dn = .true.
3733 ! If the base external interface of extending CL is the same as the
3734 ! base interior interface of the underlying CL, underlying CL will
3735 ! be automatically merged. Then, reduce total number of CL by 1.
3736 ! For a consistent treatment with 'upward' extension, I should use
3737 ! 'kbinc = kbase(i,ncv-1) - 1' instead of 'ktop(i,ncv-1) + 1' below.
3738 ! However, it seems that these two methods produce the same results.
3739 ! Note also that in contrast to upward merging, the decrease of ncv
3740 ! should be performed here.
3741 ! Note that below formula correctly works even when upperlying CL
3742 ! regime incorporates below SBCL.
3744 kbinc = 0
3745 if( ncv .gt. 1 ) kbinc = ktop(i,ncv-1) + 1
3746 if( kb .eq. kbinc ) then
3748 do k = ktop(i,ncv-1) + 1, kbase(i,ncv-1) - 1
3750 if( choice_tunl .eq. 'rampcl' ) then
3751 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3752 elseif( choice_tunl .eq. 'rampsl' ) then
3753 tunlramp = ctunl*tunl
3754 ! tunlramp = 0.765_r8
3755 else
3756 tunlramp = tunl
3757 endif
3758 if( choice_leng .eq. 'origin' ) then
3759 lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3760 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
3761 else
3762 lz = min( vk*zi(i,k), tunlramp*lbulk )
3763 endif
3765 dzinc = z(i,k-1)-z(i,k)
3766 dl2n2 = lz*lz*n2(i,k)*dzinc
3767 dl2s2 = lz*lz*s2(i,k)*dzinc
3768 dwinc = -sh*dl2n2 + sm*dl2s2
3770 lint = lint + dzinc
3771 l2n2 = l2n2 + dl2n2
3772 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3773 l2s2 = l2s2 + dl2s2
3774 wint = wint + dwinc
3776 end do
3778 ! We are incorporating interior of CL ncv-1, so merge
3779 ! this CL into the current CL.
3781 kb = kbase(i,ncv-1)
3782 ncv = ncv - 1
3783 ncvfin(i) = ncvfin(i) -1
3784 cntd = cntd + 1
3786 end if
3788 ! Calculate the contribution of potentially incorporatable CL
3789 ! base external interface. Calculate separately when the base
3790 ! of extended CL is surface and non-surface.
3792 if( kb .eq. pver + 1 ) then
3794 if( bflxs(i) .gt. 0._r8 ) then
3795 ! Calculate stability functions of surface interfacial layer
3796 gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
3797 gh_surf = gg/(alph5-gg*alph3)
3798 ! gh_surf = min(max(gh_surf,-0.28_r8),0.0233_r8)
3799 gh_surf = min(max(gh_surf,-3.5334_r8),0.0233_r8)
3800 sh_surf = alph5/(1._r8+alph3*gh_surf)
3801 sm_surf = (alph1 + alph2*gh_surf)/(1._r8+alph3*gh_surf)/(1._r8+alph4*gh_surf)
3802 ! Calculate surface interfacial layer contribution. By construction,
3803 ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'
3804 dlint_surf = z(i,pver)
3805 dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh_surf*sqrt(tkes(i)))
3806 dl2s2_surf = vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm_surf*sqrt(tkes(i)))
3807 dw_surf = (tkes(i)/b1)*z(i,pver)
3808 else
3809 dlint_surf = 0._r8
3810 dl2n2_surf = 0._r8
3811 dl2s2_surf = 0._r8
3812 dw_surf = 0._r8
3813 end if
3814 ! If (kb.eq.pver+1), updating of CL internal energetics should be
3815 ! performed here inside of 'do while' loop, since 'do while' loop
3816 ! contains the constraint of '.and.(kb.ne.pver+1)',so updating of
3817 ! CL internal energetics cannot be performed within this do while
3818 ! loop when kb.eq.pver+1. Even though I updated all 'l2n2','l2s2',
3819 ! 'wint' below, only the updated 'wint' is used in the following
3820 ! numerical calculation.
3821 lint = lint + dlint_surf
3822 l2n2 = l2n2 + dl2n2_surf
3823 l2n2 = -min(-l2n2, tkemax*lint/(b1*sh))
3824 l2s2 = l2s2 + dl2s2_surf
3825 wint = wint + dw_surf
3827 else
3829 if( choice_tunl .eq. 'rampcl' ) then
3830 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3831 elseif( choice_tunl .eq. 'rampsl' ) then
3832 tunlramp = ctunl*tunl
3833 ! tunlramp = 0.765_r8
3834 else
3835 tunlramp = tunl
3836 endif
3837 if( choice_leng .eq. 'origin' ) then
3838 lz = ( (vk*zi(i,kb))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3839 ! lz = vk*zi(i,kb) / (1._r8+vk*zi(i,kb)/(tunlramp*lbulk))
3840 else
3841 lz = min( vk*zi(i,kb), tunlramp*lbulk )
3842 endif
3844 dzinc = z(i,kb-1)-z(i,kb)
3845 dl2n2 = lz*lz*n2(i,kb)*dzinc
3846 dl2s2 = lz*lz*s2(i,kb)*dzinc
3847 dwinc = -sh*dl2n2 + sm*dl2s2
3849 end if
3851 end do ! End of merging test 'do while' loop
3853 if( (kb.eq.pver+1) .and. (ncv.ne.1) ) then
3854 write(iulog,*) 'Major mistake zisocl: the CL based at surface is not indexed 1'
3855 #ifdef WRF_PORT
3856 call wrf_message(iulog)
3857 #endif
3858 stop
3859 end if
3861 end if ! Done with bottom extension of CL
3863 ! Update CL interface indices appropriately if any CL was merged.
3864 ! Note that below only updated the interface index of merged CL,
3865 ! not the original merging CL. Updates of 'kbase' and 'ktop' of
3866 ! the original merging CL will be done later below. I should
3867 ! check in detail if below index updating is correct or not.
3869 if( cntd .gt. 0 ) then
3870 do incv = 1, ncvfin(i) - ncv
3871 kbase(i,ncv+incv) = kbase(i,ncvinit+incv)
3872 ktop(i,ncv+incv) = ktop(i,ncvinit+incv)
3873 end do
3874 end if
3876 ! Sanity check for positive wint.
3878 if( wint .lt. 0.01_r8 ) then
3879 wint = 0.01_r8
3880 end if
3882 ! -------------------------------------------------------------------------- !
3883 ! Finally update CL mean internal energetics including surface contribution !
3884 ! after finishing all the CL extension-merging process. As mentioned above, !
3885 ! there are two possible ways in the treatment of surface interfacial layer, !
3886 ! either through 'dw_surf' or 'dl2n2_surf and dl2s2_surf' by setting logical !
3887 ! variable 'use_dw_surf' =.true. or .false. In any cases, we should avoid !
3888 ! double counting of surface interfacial layer and one single consistent way !
3889 ! should be used throughout the program. !
3890 ! -------------------------------------------------------------------------- !
3892 if( extend ) then
3894 ktop(i,ncv) = kt
3895 kbase(i,ncv) = kb
3897 ! ------------------------------------------------------ !
3898 ! Step 1: Include surface interfacial layer contribution !
3899 ! ------------------------------------------------------ !
3901 lbulk = zi(i,kt) - zi(i,kb)
3902 dlint_surf = 0._r8
3903 dl2n2_surf = 0._r8
3904 dl2s2_surf = 0._r8
3905 dw_surf = 0._r8
3906 if( kb .eq. pver + 1 ) then
3907 if( bflxs(i) .gt. 0._r8 ) then
3908 ! Calculate stability functions of surface interfacial layer
3909 gg = 0.5_r8*vk*z(i,pver)*bprod(i,pver+1)/(tkes(i)**(3._r8/2._r8))
3910 gh = gg/(alph5-gg*alph3)
3911 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
3912 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3913 sh = alph5/(1._r8+alph3*gh)
3914 sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
3915 ! Calculate surface interfacial layer contribution. By construction,
3916 ! it exactly becomes 'dw_surf = -dl2n2_surf + ds2n2_surf'
3917 dlint_surf = z(i,pver)
3918 dl2n2_surf = -vk*(z(i,pver)**2._r8)*bprod(i,pver+1)/(sh*sqrt(tkes(i)))
3919 dl2s2_surf = vk*(z(i,pver)**2._r8)*sprod(i,pver+1)/(sm*sqrt(tkes(i)))
3920 dw_surf = (tkes(i)/b1)*z(i,pver)
3921 else
3922 lbulk = zi(i,kt) - z(i,pver)
3923 end if
3924 end if
3925 lint = dlint_surf
3926 l2n2 = dl2n2_surf
3927 l2s2 = dl2s2_surf
3928 wint = dw_surf
3929 if( use_dw_surf ) then
3930 l2n2 = 0._r8
3931 l2s2 = 0._r8
3932 else
3933 wint = 0._r8
3934 end if
3936 ! -------------------------------------------------------------- !
3937 ! Step 2. Include the contribution of 'pure internal interfaces' !
3938 ! -------------------------------------------------------------- !
3940 do k = kt + 1, kb - 1
3941 if( choice_tunl .eq. 'rampcl' ) then
3942 tunlramp = 0.5_r8*(1._r8+ctunl)*tunl
3943 elseif( choice_tunl .eq. 'rampsl' ) then
3944 tunlramp = ctunl*tunl
3945 ! tunlramp = 0.765_r8
3946 else
3947 tunlramp = tunl
3948 endif
3949 if( choice_leng .eq. 'origin' ) then
3950 lz = ( (vk*zi(i,k))**(-cleng) + (tunlramp*lbulk)**(-cleng) )**(-1._r8/cleng)
3951 ! lz = vk*zi(i,k) / (1._r8+vk*zi(i,k)/(tunlramp*lbulk))
3952 else
3953 lz = min( vk*zi(i,k), tunlramp*lbulk )
3954 endif
3955 dzinc = z(i,k-1) - z(i,k)
3956 lint = lint + dzinc
3957 l2n2 = l2n2 + lz*lz*n2(i,k)*dzinc
3958 l2s2 = l2s2 + lz*lz*s2(i,k)*dzinc
3959 end do
3961 ricll = min(l2n2/max(l2s2,ntzero),ricrit)
3962 trma = alph3*alph4*ricll+2._r8*b1*(alph2-alph4*alph5*ricll)
3963 trmb = ricll*(alph3+alph4)+2._r8*b1*(-alph5*ricll+alph1)
3964 trmc = ricll
3965 det = max(trmb*trmb-4._r8*trma*trmc,0._r8)
3966 gh = (-trmb + sqrt(det))/2._r8/trma
3967 ! gh = min(max(gh,-0.28_r8),0.0233_r8)
3968 gh = min(max(gh,-3.5334_r8),0.0233_r8)
3969 sh = alph5 / (1._r8+alph3*gh)
3970 sm = (alph1 + alph2*gh)/(1._r8+alph3*gh)/(1._r8+alph4*gh)
3971 ! Even though the 'wint' after finishing merging was positive, it is
3972 ! possible that re-calculated 'wint' here is negative. In this case,
3973 ! correct 'wint' to be a small positive number
3974 wint = max( wint - sh*l2n2 + sm*l2s2, 0.01_r8 )
3976 end if
3978 ! ---------------------------------------------------------------------- !
3979 ! Calculate final output variables of each CL (either has merged or not) !
3980 ! ---------------------------------------------------------------------- !
3982 lbrk(i,ncv) = lint
3983 wbrk(i,ncv) = wint/lint
3984 ebrk(i,ncv) = b1*wbrk(i,ncv)
3985 ebrk(i,ncv) = min(ebrk(i,ncv),tkemax)
3986 ricl(i,ncv) = ricll
3987 ghcl(i,ncv) = gh
3988 shcl(i,ncv) = sh
3989 smcl(i,ncv) = sm
3991 ! Increment counter for next CL. I should check if the increament of 'ncv'
3992 ! below is reasonable or not, since whenever CL is merged during downward
3993 ! extension process, 'ncv' is lowered down continuously within 'do' loop.
3994 ! But it seems that below 'ncv = ncv + 1' is perfectly correct.
3996 ncv = ncv + 1
3998 end do ! End of loop over each CL regime, ncv.
4000 ! ---------------------------------------------------------- !
4001 ! Re-initialize external interface indices which are not CLs !
4002 ! ---------------------------------------------------------- !
4004 do ncv = ncvfin(i) + 1, ncvmax
4005 ktop(i,ncv) = 0
4006 kbase(i,ncv) = 0
4007 end do
4009 ! ------------------------------------------------ !
4010 ! Update CL interface identifiers, 'belongcv' !
4011 ! CL external interfaces are also identified as CL !
4012 ! ------------------------------------------------ !
4014 do k = 1, pver + 1
4015 belongcv(i,k) = .false.
4016 end do
4018 do ncv = 1, ncvfin(i)
4019 do k = ktop(i,ncv), kbase(i,ncv)
4020 belongcv(i,k) = .true.
4021 end do
4022 end do
4024 return
4026 end subroutine zisocl
4028 real(r8) function compute_cubic(a,b,c)
4029 ! ------------------------------------------------------------------------- !
4030 ! Solve canonical cubic : x^3 + a*x^2 + b*x + c = 0, x = sqrt(e)/sqrt(<e>) !
4031 ! Set x = max(xmin,x) at the end !
4032 ! ------------------------------------------------------------------------- !
4033 implicit none
4034 real(r8), intent(in) :: a, b, c
4035 real(r8) qq, rr, dd, theta, aa, bb, x1, x2, x3
4036 real(r8), parameter :: xmin = 1.e-2_r8
4038 qq = (a**2-3._r8*b)/9._r8
4039 rr = (2._r8*a**3 - 9._r8*a*b + 27._r8*c)/54._r8
4041 dd = rr**2 - qq**3
4042 if( dd .le. 0._r8 ) then
4043 theta = acos(rr/qq**(3._r8/2._r8))
4044 x1 = -2._r8*sqrt(qq)*cos(theta/3._r8) - a/3._r8
4045 x2 = -2._r8*sqrt(qq)*cos((theta+2._r8*3.141592)/3._r8) - a/3._r8
4046 x3 = -2._r8*sqrt(qq)*cos((theta-2._r8*3.141592)/3._r8) - a/3._r8
4047 compute_cubic = max(max(max(x1,x2),x3),xmin)
4048 return
4049 else
4050 if( rr .ge. 0._r8 ) then
4051 aa = -(sqrt(rr**2-qq**3)+rr)**(1._r8/3._r8)
4052 else
4053 aa = (sqrt(rr**2-qq**3)-rr)**(1._r8/3._r8)
4054 endif
4055 if( aa .eq. 0._r8 ) then
4056 bb = 0._r8
4057 else
4058 bb = qq/aa
4059 endif
4060 compute_cubic = max((aa+bb)-a/3._r8,xmin)
4061 return
4062 endif
4064 return
4065 end function compute_cubic
4067 END MODULE eddy_diff