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