Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cam_bl_eddy_diff.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 224.4 KB

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

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

Download in other formats:

Original Format