3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
6 use cam_history, only: outfld, addfld, phys_decomp
8 use module_cam_support, only: outfld, addfld, phys_decomp
10 use module_state_description, only: CAMUWPBLSCHEME, MYJPBLSCHEME
12 use module_state_description, only: CAMUWPBLSCHEME, MYJPBLSCHEME, &
16 use error_function, only: erfc
18 use cam_logfile, only: iulog
19 use ppgrid, only: pcols, pver, pverp
20 use abortutils, only: endrun
21 use spmd_utils, only: masterproc
23 use module_cam_support, only: iulog, pcols, pver, pverp, endrun, masterproc
37 integer , parameter :: r8 = selected_real_kind(12) ! 8 byte real
38 real(r8), parameter :: unset_r8 = huge(1.0_r8)
39 real(r8) :: xlv ! Latent heat of vaporization
40 real(r8) :: xlf ! Latent heat of fusion
41 real(r8) :: xls ! Latent heat of sublimation = xlv + xlf
42 real(r8) :: cp ! Specific heat of dry air
43 real(r8) :: zvir ! rh2o/rair - 1
44 real(r8) :: r ! Gas constant for dry air
45 real(r8) :: g ! Gravitational constant
46 real(r8) :: ep2 ! mol wgt water vapor / mol wgt dry air
47 real(r8) :: p00 ! Reference pressure for exner function
48 real(r8) :: rovcp ! R/cp
50 ! Tuning parameters set via namelist
51 real(r8) :: rpen ! For penetrative entrainment efficiency
53 !===============================================================================
55 !===============================================================================
57 real(r8) function exnf(pressure)
58 real(r8), intent(in) :: pressure
59 exnf = (pressure/p00)**rovcp
63 !===============================================================================
65 subroutine uwshcu_readnl(nlfile)
67 use namelist_utils, only: find_group_name
68 use units, only: getunit, freeunit
72 character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
76 integer :: unitn, ierr
77 character(len=*), parameter :: subname = 'uwshcu_readnl'
80 real(r8) :: uwshcu_rpen = unset_r8 ! For penetrative entrainment efficiency
82 namelist /uwshcu_nl/ uwshcu_rpen
83 !-----------------------------------------------------------------------------
87 open( unitn, file=trim(nlfile), status='old' )
88 call find_group_name(unitn, 'uwshcu_nl', status=ierr)
90 read(unitn, uwshcu_nl, iostat=ierr)
92 call endrun(subname // ':: ERROR reading namelist')
100 ! Broadcast namelist variables
101 call mpibcast(uwshcu_rpen, 1, mpir8, 0, mpicom)
106 !Balwinder.Singh@pnnl.gov: rpen value has been obtained from previous incarnation
107 !of this module (see CESM101)
111 end subroutine uwshcu_readnl
113 !===============================================================================
115 subroutine init_uwshcu( kind, xlv_in, cp_in, xlf_in, zvir_in, r_in, g_in, ep2_in &
117 , rushten, rvshten, rthshten, rqvshten, &
118 rqcshten, rqrshten, rqishten, rqsshten, rqgshten, &
119 p_qc, p_qr, p_qi, p_qs, p_qg, &
120 bl_pbl_physics, param_first_scalar, restart, &
121 ids, ide, jds, jde, kds, kde, &
122 ims, ime, jms, jme, kms, kme, &
123 its, ite, jts, jte, kts, kte )
127 !~add in zeroing of things like rliq that may come from ZM or may not
129 !------------------------------------------------------------- !
131 ! Initialize key constants for the shallow convection package. !
132 !------------------------------------------------------------- !
134 use cam_history, only: outfld, addfld, phys_decomp
135 use ppgrid, only: pcols, pver, pverp
137 use module_cam_support, only: outfld, addfld, phys_decomp, pcols, pver, pverp
140 integer , intent(in) :: kind ! kind of reals being passed in
141 real(r8), intent(in) :: xlv_in ! Latent heat of vaporization
142 real(r8), intent(in) :: xlf_in ! Latent heat of fusion
143 real(r8), intent(in) :: cp_in ! Specific heat of dry air
144 real(r8), intent(in) :: zvir_in ! rh2o/rair - 1
145 real(r8), intent(in) :: r_in ! Gas constant for dry air
146 real(r8), intent(in) :: g_in ! Gravitational constant
147 real(r8), intent(in) :: ep2_in ! mol wgt water vapor / mol wgt dry air
149 character(len=*), parameter :: subname = 'init_uwshcu'
151 LOGICAL , INTENT(IN) :: restart
152 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
153 ims, ime, jms, jme, kms, kme, &
154 its, ite, jts, jte, kts, kte, &
155 p_qc, p_qr, p_qi, p_qs, p_qg, &
156 param_first_scalar, &
159 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: &
170 integer :: i, itf, j, jtf, k, ktf
177 !Balwinder.Singh@pnnl.gov - Call to assign rpen
178 call uwshcu_readnl('dummyString')
181 ! The CAM UW shallow Cu scheme only works with PBL schemes that generate
182 ! TKE. So, for example, YSU would be an invalid choice to combine with
183 ! the UW scheme. For now, we are coupling into the tke_myj version of
184 ! TKE used by MYJ and CAMUWPBL. At some point, it would be great to unify
185 ! all the different TKEs output so that things work together properly.
186 ! For now, we are not trying to deal with the cacophony of variables,
187 ! i.e. tke & tke_pbl.
189 select case(bl_pbl_physics)
190 #if ( NMM_CORE == 1 )
191 case (CAMUWPBLSCHEME, MYJPBLSCHEME)
193 case (CAMUWPBLSCHEME, MYJPBLSCHEME, MYNNPBLSCHEME)
195 !These are acceptable PBL schemes.
197 call wrf_error_fatal("The CAMUWSHCU scheme requires CAMUWPBLSCHEME, MYJPBLSCHEME or MYNN.")
200 ! Initialize module_cam_support variables...
201 ! This could be moved to a master "cam-init" subroutine once multiple CAM
202 ! parameterizations are implemented in WRF. For now, it doesn't hurt to
203 ! have these possibly initialized more than once since they are
204 ! essentially constant.
209 ! ------------------------- !
210 ! Internal Output Variables !
211 ! ------------------------- !
213 call addfld( 'qtflx_Cu' , 'kg/m2/s' , pverp , 'A' , 'Convective qt flux' , phys_decomp )
214 call addfld( 'slflx_Cu' , 'J/m2/s' , pverp , 'A' , 'Convective sl flux' , phys_decomp )
215 call addfld( 'uflx_Cu' , 'kg/m/s2' , pverp , 'A' , 'Convective u flux' , phys_decomp )
216 call addfld( 'vflx_Cu' , 'kg/m/s2' , pverp , 'A' , 'Convective v flux' , phys_decomp )
218 call addfld( 'qtten_Cu' , 'kg/kg/s' , pver , 'A' , 'qt tendency by convection' , phys_decomp )
219 call addfld( 'slten_Cu' , 'J/kg/s' , pver , 'A' , 'sl tendency by convection' , phys_decomp )
220 call addfld( 'uten_Cu' , 'm/s2' , pver , 'A' , ' u tendency by convection' , phys_decomp )
221 call addfld( 'vten_Cu' , 'm/s2' , pver , 'A' , ' v tendency by convection' , phys_decomp )
222 call addfld( 'qvten_Cu' , 'kg/kg/s' , pver , 'A' , 'qv tendency by convection' , phys_decomp )
223 call addfld( 'qlten_Cu' , 'kg/kg/s' , pver , 'A' , 'ql tendency by convection' , phys_decomp )
224 call addfld( 'qiten_Cu' , 'kg/kg/s' , pver , 'A' , 'qi tendency by convection' , phys_decomp )
226 call addfld( 'cbmf_Cu' , 'kg/m2/s' , 1 , 'A' , 'Cumulus base mass flux' , phys_decomp )
227 call addfld( 'ufrcinvbase_Cu' , 'fraction', 1 , 'A' , 'Cumulus fraction at PBL top' , phys_decomp )
228 call addfld( 'ufrclcl_Cu' , 'fraction', 1 , 'A' , 'Cumulus fraction at LCL' , phys_decomp )
229 call addfld( 'winvbase_Cu' , 'm/s' , 1 , 'A' , 'Cumulus vertical velocity at PBL top' , phys_decomp )
230 call addfld( 'wlcl_Cu' , 'm/s' , 1 , 'A' , 'Cumulus vertical velocity at LCL' , phys_decomp )
231 call addfld( 'plcl_Cu' , 'Pa' , 1 , 'A' , 'LCL of source air' , phys_decomp )
232 call addfld( 'pinv_Cu' , 'Pa' , 1 , 'A' , 'PBL top pressure' , phys_decomp )
233 call addfld( 'plfc_Cu' , 'Pa' , 1 , 'A' , 'LFC of source air' , phys_decomp )
234 call addfld( 'pbup_Cu' , 'Pa' , 1 , 'A' , 'Highest interface level of positive cumulus buoyancy', phys_decomp )
235 call addfld( 'ppen_Cu' , 'Pa' , 1 , 'A' , 'Highest level where cumulus w is 0' , phys_decomp )
236 call addfld( 'qtsrc_Cu' , 'kg/kg' , 1 , 'A' , 'Cumulus source air qt' , phys_decomp )
237 call addfld( 'thlsrc_Cu' , 'K' , 1 , 'A' , 'Cumulus source air thl' , phys_decomp )
238 call addfld( 'thvlsrc_Cu' , 'K' , 1 , 'A' , 'Cumulus source air thvl' , phys_decomp )
239 call addfld( 'emfkbup_Cu' , 'kg/m2/s' , 1 , 'A' , 'Penetrative mass flux at kbup' , phys_decomp )
240 call addfld( 'cin_Cu' , 'J/kg' , 1 , 'A' , 'CIN upto LFC' , phys_decomp )
241 call addfld( 'cinlcl_Cu' , 'J/kg' , 1 , 'A' , 'CIN upto LCL' , phys_decomp )
242 call addfld( 'cbmflimit_Cu' , 'kg/m2/s' , 1 , 'A' , 'cbmflimiter' , phys_decomp )
243 call addfld( 'tkeavg_Cu' , 'm2/s2' , 1 , 'A' , 'Average tke within PBL for convection scheme' , phys_decomp )
244 call addfld( 'zinv_Cu' , 'm' , 1 , 'A' , 'PBL top height' , phys_decomp )
245 call addfld( 'rcwp_Cu' , 'kg/m2' , 1 , 'A' , 'Cumulus LWP+IWP' , phys_decomp )
246 call addfld( 'rlwp_Cu' , 'kg/m2' , 1 , 'A' , 'Cumulus LWP' , phys_decomp )
247 call addfld( 'riwp_Cu' , 'kg/m2' , 1 , 'A' , 'Cumulus IWP' , phys_decomp )
248 call addfld( 'tophgt_Cu' , 'm' , 1 , 'A' , 'Cumulus top height' , phys_decomp )
250 call addfld( 'wu_Cu' , 'm/s' , pverp , 'A' , 'Convective updraft vertical velocity' , phys_decomp )
251 call addfld( 'ufrc_Cu' , 'fraction', pverp , 'A' , 'Convective updraft fractional area' , phys_decomp )
252 call addfld( 'qtu_Cu' , 'kg/kg' , pverp , 'A' , 'Cumulus updraft qt' , phys_decomp )
253 call addfld( 'thlu_Cu' , 'K' , pverp , 'A' , 'Cumulus updraft thl' , phys_decomp )
254 call addfld( 'thvu_Cu' , 'K' , pverp , 'A' , 'Cumulus updraft thv' , phys_decomp )
255 call addfld( 'uu_Cu' , 'm/s' , pverp , 'A' , 'Cumulus updraft uwnd' , phys_decomp )
256 call addfld( 'vu_Cu' , 'm/s' , pverp , 'A' , 'Cumulus updraft vwnd' , phys_decomp )
257 call addfld( 'qtu_emf_Cu' , 'kg/kg' , pverp , 'A' , 'qt of penatratively entrained air' , phys_decomp )
258 call addfld( 'thlu_emf_Cu' , 'K' , pverp , 'A' , 'thl of penatratively entrained air' , phys_decomp )
259 call addfld( 'uu_emf_Cu' , 'm/s' , pverp , 'A' , 'uwnd of penatratively entrained air' , phys_decomp )
260 call addfld( 'vu_emf_Cu' , 'm/s' , pverp , 'A' , 'vwnd of penatratively entrained air' , phys_decomp )
261 call addfld( 'umf_Cu' , 'kg/m2/s' , pverp , 'A' , 'Cumulus updraft mass flux' , phys_decomp )
262 call addfld( 'uemf_Cu' , 'kg/m2/s' , pverp , 'A' , 'Cumulus net ( updraft + entrainment ) mass flux' , phys_decomp )
263 call addfld( 'qcu_Cu' , 'kg/kg' , pver , 'A' , 'Cumulus updraft LWC+IWC' , phys_decomp )
264 call addfld( 'qlu_Cu' , 'kg/kg' , pver , 'A' , 'Cumulus updraft LWC' , phys_decomp )
265 call addfld( 'qiu_Cu' , 'kg/kg' , pver , 'A' , 'Cumulus updraft IWC' , phys_decomp )
266 call addfld( 'cufrc_Cu' , 'fraction', pver , 'A' , 'Cumulus cloud fraction' , phys_decomp )
267 call addfld( 'fer_Cu' , '1/m' , pver , 'A' , 'Cumulus lateral fractional entrainment rate' , phys_decomp )
268 call addfld( 'fdr_Cu' , '1/m' , pver , 'A' , 'Cumulus lateral fractional detrainment Rate' , phys_decomp )
270 call addfld( 'dwten_Cu' , 'kg/kg/s' , pver , 'A' , 'Expellsion rate of cumulus cloud water to env.' , phys_decomp )
271 call addfld( 'diten_Cu' , 'kg/kg/s' , pver , 'A' , 'Expellsion rate of cumulus ice water to env.' , phys_decomp )
272 call addfld( 'qrten_Cu' , 'kg/kg/s' , pver , 'A' , 'Production rate of rain by cumulus' , phys_decomp )
273 call addfld( 'qsten_Cu' , 'kg/kg/s' , pver , 'A' , 'Production rate of snow by cumulus' , phys_decomp )
274 call addfld( 'flxrain_Cu' , 'kg/m2/s' , pverp , 'A' , 'Rain flux induced by Cumulus' , phys_decomp )
275 call addfld( 'flxsnow_Cu' , 'kg/m2/s' , pverp , 'A' , 'Snow flux induced by Cumulus' , phys_decomp )
276 call addfld( 'ntraprd_Cu' , 'kg/kg/s' , pver , 'A' , 'Net production rate of rain by Cumulus' , phys_decomp )
277 call addfld( 'ntsnprd_Cu' , 'kg/kg/s' , pver , 'A' , 'Net production rate of snow by Cumulus' , phys_decomp )
279 call addfld( 'excessu_Cu' , 'no' , pver , 'A' , 'Updraft saturation excess' , phys_decomp )
280 call addfld( 'excess0_Cu' , 'no' , pver , 'A' , 'Environmental saturation excess' , phys_decomp )
281 call addfld( 'xc_Cu' , 'no' , pver , 'A' , 'Critical mixing ratio' , phys_decomp )
282 call addfld( 'aquad_Cu' , 'no' , pver , 'A' , 'aquad' , phys_decomp )
283 call addfld( 'bquad_Cu' , 'no' , pver , 'A' , 'bquad' , phys_decomp )
284 call addfld( 'cquad_Cu' , 'no' , pver , 'A' , 'cquad' , phys_decomp )
285 call addfld( 'bogbot_Cu' , 'no' , pver , 'A' , 'Cloud buoyancy at the bottom interface' , phys_decomp )
286 call addfld( 'bogtop_Cu' , 'no' , pver , 'A' , 'Cloud buoyancy at the top interface' , phys_decomp )
288 call addfld('exit_UWCu_Cu' , 'no' , 1 , 'A' , 'exit_UWCu' , phys_decomp )
289 call addfld('exit_conden_Cu' , 'no' , 1 , 'A' , 'exit_conden' , phys_decomp )
290 call addfld('exit_klclmkx_Cu' , 'no' , 1 , 'A' , 'exit_klclmkx' , phys_decomp )
291 call addfld('exit_klfcmkx_Cu' , 'no' , 1 , 'A' , 'exit_klfcmkx' , phys_decomp )
292 call addfld('exit_ufrc_Cu' , 'no' , 1 , 'A' , 'exit_ufrc' , phys_decomp )
293 call addfld('exit_wtw_Cu' , 'no' , 1 , 'A' , 'exit_wtw' , phys_decomp )
294 call addfld('exit_drycore_Cu' , 'no' , 1 , 'A' , 'exit_drycore' , phys_decomp )
295 call addfld('exit_wu_Cu' , 'no' , 1 , 'A' , 'exit_wu' , phys_decomp )
296 call addfld('exit_cufilter_Cu', 'no' , 1 , 'A' , 'exit_cufilter' , phys_decomp )
297 call addfld('exit_kinv1_Cu' , 'no' , 1 , 'A' , 'exit_kinv1' , phys_decomp )
298 call addfld('exit_rei_Cu' , 'no' , 1 , 'A' , 'exit_rei' , phys_decomp )
300 call addfld('limit_shcu_Cu' , 'no' , 1 , 'A' , 'limit_shcu' , phys_decomp )
301 call addfld('limit_negcon_Cu' , 'no' , 1 , 'A' , 'limit_negcon' , phys_decomp )
302 call addfld('limit_ufrc_Cu' , 'no' , 1 , 'A' , 'limit_ufrc' , phys_decomp )
303 call addfld('limit_ppen_Cu' , 'no' , 1 , 'A' , 'limit_ppen' , phys_decomp )
304 call addfld('limit_emf_Cu' , 'no' , 1 , 'A' , 'limit_emf' , phys_decomp )
305 call addfld('limit_cinlcl_Cu' , 'no' , 1 , 'A' , 'limit_cinlcl' , phys_decomp )
306 call addfld('limit_cin_Cu' , 'no' , 1 , 'A' , 'limit_cin' , phys_decomp )
307 call addfld('limit_cbmf_Cu' , 'no' , 1 , 'A' , 'limit_cbmf' , phys_decomp )
308 call addfld('limit_rei_Cu' , 'no' , 1 , 'A' , 'limit_rei' , phys_decomp )
309 call addfld('ind_delcin_Cu' , 'no' , 1 , 'A' , 'ind_delcin' , phys_decomp )
311 if( kind .ne. r8 ) then
312 write(iulog,*) subname//': ERROR -- real KIND does not match internal specification.'
314 call wrf_message(iulog)
316 call endrun(subname//': ERROR -- real KIND does not match internal specification.')
330 if (rpen == unset_r8) then
331 call endrun(subname//': uwshcu_rpen must be set in the namelist')
334 if ( masterproc ) then
335 write(iulog,*) subname//': tuning parameters: rpen=',rpen
337 call wrf_debug(1,iulog)
343 ! Set initial values for WRF arrays dependent on restart conditions
353 if( p_qc > param_first_scalar ) rqcshten(i,k,j) = 0.
354 if( p_qr > param_first_scalar ) rqrshten(i,k,j) = 0.
355 if( p_qi > param_first_scalar ) rqishten(i,k,j) = 0.
356 if( p_qs > param_first_scalar ) rqsshten(i,k,j) = 0.
357 if( p_qg > param_first_scalar ) rqgshten(i,k,j) = 0.
364 end subroutine init_uwshcu
366 subroutine compute_uwshcu_inv( mix , mkx , iend , ncnst , dt , &
367 ps0_inv , zs0_inv , p0_inv , z0_inv , dp0_inv , &
368 u0_inv , v0_inv , qv0_inv , ql0_inv , qi0_inv , &
369 t0_inv , s0_inv , tr0_inv , &
370 tke_inv , cldfrct_inv, concldfrct_inv, pblh , cush , &
371 umf_inv , slflx_inv , qtflx_inv , &
372 flxprc1_inv, flxsnow1_inv, &
373 qvten_inv, qlten_inv , qiten_inv , &
374 sten_inv , uten_inv , vten_inv , trten_inv , &
375 qrten_inv, qsten_inv , precip , snow , evapc_inv, &
376 cufrc_inv, qcu_inv , qlu_inv , qiu_inv , &
377 cbmf , qc_inv , rliq , &
378 cnt_inv , cnb_inv , qsat , lchnk , dpdry0_inv )
381 integer , intent(in) :: lchnk
382 integer , intent(in) :: mix
383 integer , intent(in) :: mkx
384 integer , intent(in) :: iend
385 integer , intent(in) :: ncnst
386 real(r8), intent(in) :: dt ! Time step : 2*delta_t [ s ]
387 real(r8), intent(in) :: ps0_inv(mix,mkx+1) ! Environmental pressure at the interfaces [ Pa ]
388 real(r8), intent(in) :: zs0_inv(mix,mkx+1) ! Environmental height at the interfaces [ m ]
389 real(r8), intent(in) :: p0_inv(mix,mkx) ! Environmental pressure at the layer mid-point [ Pa ]
390 real(r8), intent(in) :: z0_inv(mix,mkx) ! Environmental height at the layer mid-point [ m ]
391 real(r8), intent(in) :: dp0_inv(mix,mkx) ! Environmental layer pressure thickness [ Pa ] > 0.
392 real(r8), intent(in) :: dpdry0_inv(mix,mkx) ! Environmental dry layer pressure thickness [ Pa ]
393 real(r8), intent(in) :: u0_inv(mix,mkx) ! Environmental zonal wind [ m/s ]
394 real(r8), intent(in) :: v0_inv(mix,mkx) ! Environmental meridional wind [ m/s ]
395 real(r8), intent(in) :: qv0_inv(mix,mkx) ! Environmental water vapor specific humidity [ kg/kg ]
396 real(r8), intent(in) :: ql0_inv(mix,mkx) ! Environmental liquid water specific humidity [ kg/kg ]
397 real(r8), intent(in) :: qi0_inv(mix,mkx) ! Environmental ice specific humidity [ kg/kg ]
398 real(r8), intent(in) :: t0_inv(mix,mkx) ! Environmental temperature [ K ]
399 real(r8), intent(in) :: s0_inv(mix,mkx) ! Environmental dry static energy [ J/kg ]
400 real(r8), intent(in) :: tr0_inv(mix,mkx,ncnst) ! Environmental tracers [ #, kg/kg ]
401 real(r8), intent(in) :: tke_inv(mix,mkx+1) ! Turbulent kinetic energy at the interfaces [ m2/s2 ]
402 real(r8), intent(in) :: cldfrct_inv(mix,mkx) ! Total cloud fraction at the previous time step [ fraction ]
403 real(r8), intent(in) :: concldfrct_inv(mix,mkx) ! Total convective ( shallow + deep ) cloud fraction at the previous time step [ fraction ]
404 real(r8), intent(in) :: pblh(mix) ! Height of PBL [ m ]
405 real(r8), intent(inout) :: cush(mix) ! Convective scale height [ m ]
406 real(r8), intent(out) :: umf_inv(mix,mkx+1) ! Updraft mass flux at the interfaces [ kg/m2/s ]
407 real(r8), intent(out) :: qvten_inv(mix,mkx) ! Tendency of water vapor specific humidity [ kg/kg/s ]
408 real(r8), intent(out) :: qlten_inv(mix,mkx) ! Tendency of liquid water specific humidity [ kg/kg/s ]
409 real(r8), intent(out) :: qiten_inv(mix,mkx) ! Tendency of ice specific humidity [ kg/kg/s ]
410 real(r8), intent(out) :: sten_inv(mix,mkx) ! Tendency of dry static energy [ J/kg/s ]
411 real(r8), intent(out) :: uten_inv(mix,mkx) ! Tendency of zonal wind [ m/s2 ]
412 real(r8), intent(out) :: vten_inv(mix,mkx) ! Tendency of meridional wind [ m/s2 ]
413 real(r8), intent(out) :: trten_inv(mix,mkx,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ]
414 real(r8), intent(out) :: qrten_inv(mix,mkx) ! Tendency of rain water specific humidity [ kg/kg/s ]
415 real(r8), intent(out) :: qsten_inv(mix,mkx) ! Tendency of snow specific humidity [ kg/kg/s ]
416 real(r8), intent(out) :: precip(mix) ! Precipitation ( rain + snow ) flux at the surface [ m/s ]
417 real(r8), intent(out) :: snow(mix) ! Snow flux at the surface [ m/s ]
418 real(r8), intent(out) :: evapc_inv(mix,mkx) ! Evaporation of precipitation [ kg/kg/s ]
419 real(r8), intent(out) :: rliq(mix) ! Vertical integral of tendency of detrained cloud condensate qc [ m/s ]
420 real(r8), intent(out) :: slflx_inv(mix,mkx+1) ! Updraft liquid static energy flux [ J/kg * kg/m2/s ]
421 real(r8), intent(out) :: qtflx_inv(mix,mkx+1) ! Updraft total water flux [ kg/kg * kg/m2/s ]
422 real(r8), intent(out) :: flxprc1_inv(mix,mkx+1) ! uw grid-box mean rain+snow flux (kg m^-2 s^-1) for physics buffer calls in convect_shallow.F90
423 real(r8), intent(out) :: flxsnow1_inv(mix,mkx+1) ! uw grid-box mean snow flux (kg m^-2 s^-1) for physics buffer calls in convect_shallow.F90
425 real(r8), intent(out) :: cufrc_inv(mix,mkx) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ]
426 real(r8), intent(out) :: qcu_inv(mix,mkx) ! Liquid+ice specific humidity within cumulus updraft [ kg/kg ]
427 real(r8), intent(out) :: qlu_inv(mix,mkx) ! Liquid water specific humidity within cumulus updraft [ kg/kg ]
428 real(r8), intent(out) :: qiu_inv(mix,mkx) ! Ice specific humidity within cumulus updraft [ kg/kg ]
429 real(r8), intent(out) :: qc_inv(mix,mkx) ! Tendency of cumulus condensate detrained into the environment [ kg/kg/s ]
430 real(r8), intent(out) :: cbmf(mix) ! Cumulus base mass flux [ kg/m2/s ]
431 real(r8), intent(out) :: cnt_inv(mix) ! Cumulus top interface index, cnt = kpen [ no ]
432 real(r8), intent(out) :: cnb_inv(mix) ! Cumulus base interface index, cnb = krel - 1 [ no ]
433 integer , external :: qsat ! Function pointer to sat vap pressure function
435 real(r8) :: ps0(mix,0:mkx) ! Environmental pressure at the interfaces [ Pa ]
436 real(r8) :: zs0(mix,0:mkx) ! Environmental height at the interfaces [ m ]
437 real(r8) :: p0(mix,mkx) ! Environmental pressure at the layer mid-point [ Pa ]
438 real(r8) :: z0(mix,mkx) ! Environmental height at the layer mid-point [ m ]
439 real(r8) :: dp0(mix,mkx) ! Environmental layer pressure thickness [ Pa ] > 0.
440 real(r8) :: dpdry0(mix,mkx) ! Environmental dry layer pressure thickness [ Pa ]
441 real(r8) :: u0(mix,mkx) ! Environmental zonal wind [ m/s ]
442 real(r8) :: v0(mix,mkx) ! Environmental meridional wind [ m/s ]
443 real(r8) :: tke(mix,0:mkx) ! Turbulent kinetic energy at the interfaces [ m2/s2 ]
444 real(r8) :: cldfrct(mix,mkx) ! Total cloud fraction at the previous time step [ fraction ]
445 real(r8) :: concldfrct(mix,mkx) ! Total convective ( shallow + deep ) cloud fraction at the previous time step [ fraction ]
446 real(r8) :: qv0(mix,mkx) ! Environmental water vapor specific humidity [ kg/kg ]
447 real(r8) :: ql0(mix,mkx) ! Environmental liquid water specific humidity [ kg/kg ]
448 real(r8) :: qi0(mix,mkx) ! Environmental ice specific humidity [ kg/kg ]
449 real(r8) :: t0(mix,mkx) ! Environmental temperature [ K ]
450 real(r8) :: s0(mix,mkx) ! Environmental dry static energy [ J/kg ]
451 real(r8) :: tr0(mix,mkx,ncnst) ! Environmental tracers [ #, kg/kg ]
452 real(r8) :: umf(mix,0:mkx) ! Updraft mass flux at the interfaces [ kg/m2/s ]
453 real(r8) :: qvten(mix,mkx) ! Tendency of water vapor specific humidity [ kg/kg/s ]
454 real(r8) :: qlten(mix,mkx) ! Tendency of liquid water specific humidity [ kg/kg/s ]
455 real(r8) :: qiten(mix,mkx) ! tendency of ice specific humidity [ kg/kg/s ]
456 real(r8) :: sten(mix,mkx) ! Tendency of static energy [ J/kg/s ]
457 real(r8) :: uten(mix,mkx) ! Tendency of zonal wind [ m/s2 ]
458 real(r8) :: vten(mix,mkx) ! Tendency of meridional wind [ m/s2 ]
459 real(r8) :: trten(mix,mkx,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ]
460 real(r8) :: qrten(mix,mkx) ! Tendency of rain water specific humidity [ kg/kg/s ]
461 real(r8) :: qsten(mix,mkx) ! Tendency of snow speficif humidity [ kg/kg/s ]
462 real(r8) :: evapc(mix,mkx) ! Tendency of evaporation of precipitation [ kg/kg/s ]
463 real(r8) :: slflx(mix,0:mkx) ! Updraft liquid static energy flux [ J/kg * kg/m2/s ]
464 real(r8) :: qtflx(mix,0:mkx) ! Updraft total water flux [ kg/kg * kg/m2/s ]
465 real(r8) :: flxprc1(mix,0:mkx) ! uw grid-box mean rain+snow flux (kg m^-2 s^-1) for physics buffer calls in convect_shallow.F90
466 real(r8) :: flxsnow1(mix,0:mkx) ! uw grid-box mean snow flux (kg m^-2 s^-1) for physics buffer calls in convect_shallow.F90
467 real(r8) :: cufrc(mix,mkx) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ]
468 real(r8) :: qcu(mix,mkx) ! Condensate water specific humidity within cumulus updraft at the layer mid-point [ kg/kg ]
469 real(r8) :: qlu(mix,mkx) ! Liquid water specific humidity within cumulus updraft at the layer mid-point [ kg/kg ]
470 real(r8) :: qiu(mix,mkx) ! Ice specific humidity within cumulus updraft at the layer mid-point [ kg/kg ]
471 real(r8) :: qc(mix,mkx) ! Tendency of cumulus condensate detrained into the environment [ kg/kg/s ]
472 real(r8) :: cnt(mix) ! Cumulus top interface index, cnt = kpen [ no ]
473 real(r8) :: cnb(mix) ! Cumulus base interface index, cnb = krel - 1 [ no ]
474 integer :: k ! Vertical index for local fields [ no ]
475 integer :: k_inv ! Vertical index for incoming fields [ no ]
476 integer :: m ! Tracer index [ no ]
480 p0(:iend,k) = p0_inv(:iend,k_inv)
481 u0(:iend,k) = u0_inv(:iend,k_inv)
482 v0(:iend,k) = v0_inv(:iend,k_inv)
483 z0(:iend,k) = z0_inv(:iend,k_inv)
484 dp0(:iend,k) = dp0_inv(:iend,k_inv)
485 dpdry0(:iend,k) = dpdry0_inv(:iend,k_inv)
486 qv0(:iend,k) = qv0_inv(:iend,k_inv)
487 ql0(:iend,k) = ql0_inv(:iend,k_inv)
488 qi0(:iend,k) = qi0_inv(:iend,k_inv)
489 t0(:iend,k) = t0_inv(:iend,k_inv)
490 s0(:iend,k) = s0_inv(:iend,k_inv)
491 cldfrct(:iend,k) = cldfrct_inv(:iend,k_inv)
492 concldfrct(:iend,k) = concldfrct_inv(:iend,k_inv)
494 tr0(:iend,k,m) = tr0_inv(:iend,k_inv,m)
500 ps0(:iend,k) = ps0_inv(:iend,k_inv)
501 zs0(:iend,k) = zs0_inv(:iend,k_inv)
502 tke(:iend,k) = tke_inv(:iend,k_inv)
505 call compute_uwshcu( mix , mkx , iend , ncnst , dt , &
506 ps0 , zs0 , p0 , z0 , dp0 , &
507 u0 , v0 , qv0 , ql0 , qi0 , &
509 tke , cldfrct, concldfrct, pblh , cush , &
510 umf , slflx , qtflx , &
511 flxprc1 , flxsnow1 , &
512 qvten, qlten , qiten , &
513 sten , uten , vten , trten , &
514 qrten, qsten , precip , snow , evapc, &
515 cufrc, qcu , qlu , qiu , &
517 cnt , cnb , qsat , lchnk , dpdry0 )
519 ! Reverse cloud top/base interface indices
521 cnt_inv(:iend) = mkx + 1 - cnt(:iend)
522 cnb_inv(:iend) = mkx + 1 - cnb(:iend)
526 umf_inv(:iend,k_inv) = umf(:iend,k)
527 slflx_inv(:iend,k_inv) = slflx(:iend,k)
528 qtflx_inv(:iend,k_inv) = qtflx(:iend,k)
529 flxprc1_inv(:iend,k_inv) = flxprc1(:iend,k) ! reversed for output to cam
530 flxsnow1_inv(:iend,k_inv) = flxsnow1(:iend,k) ! ""
535 qvten_inv(:iend,k_inv) = qvten(:iend,k)
536 qlten_inv(:iend,k_inv) = qlten(:iend,k)
537 qiten_inv(:iend,k_inv) = qiten(:iend,k)
538 sten_inv(:iend,k_inv) = sten(:iend,k)
539 uten_inv(:iend,k_inv) = uten(:iend,k)
540 vten_inv(:iend,k_inv) = vten(:iend,k)
541 qrten_inv(:iend,k_inv) = qrten(:iend,k)
542 qsten_inv(:iend,k_inv) = qsten(:iend,k)
543 evapc_inv(:iend,k_inv) = evapc(:iend,k)
544 cufrc_inv(:iend,k_inv) = cufrc(:iend,k)
545 qcu_inv(:iend,k_inv) = qcu(:iend,k)
546 qlu_inv(:iend,k_inv) = qlu(:iend,k)
547 qiu_inv(:iend,k_inv) = qiu(:iend,k)
548 qc_inv(:iend,k_inv) = qc(:iend,k)
550 trten_inv(:iend,k_inv,m) = trten(:iend,k,m)
555 end subroutine compute_uwshcu_inv
557 subroutine compute_uwshcu( mix , mkx , iend , ncnst , dt , &
558 ps0_in , zs0_in , p0_in , z0_in , dp0_in , &
559 u0_in , v0_in , qv0_in , ql0_in , qi0_in , &
560 t0_in , s0_in , tr0_in , &
561 tke_in , cldfrct_in, concldfrct_in, pblh_in , cush_inout, &
562 umf_out , slflx_out , qtflx_out , &
563 flxprc1_out , flxsnow1_out , &
564 qvten_out, qlten_out , qiten_out , &
565 sten_out , uten_out , vten_out , trten_out, &
566 qrten_out, qsten_out , precip_out , snow_out , evapc_out , &
567 cufrc_out, qcu_out , qlu_out , qiu_out , &
568 cbmf_out , qc_out , rliq_out , &
569 cnt_out , cnb_out , qsat , lchnk , dpdry0_in )
571 ! ------------------------------------------------------------ !
573 ! University of Washington Shallow Convection Scheme !
575 ! Described in Park and Bretherton. 2008. J. Climate : !
577 ! 'The University of Washington shallow convection and !
578 ! moist turbulent schemes and their impact on climate !
579 ! simulations with the Community Atmosphere Model' !
581 ! Coded by Sungsu Park. Oct.2005. !
583 ! For questions, send an email to sungsup@ucar.edu or !
584 ! sungsu@atmos.washington.edu !
586 ! ------------------------------------------------------------ !
588 use cam_history, only : outfld, addfld, phys_decomp
590 use module_cam_support, only : outfld, addfld, phys_decomp
592 use constituents, only : qmin, cnst_get_type_byind, cnst_get_ind
594 use modal_aero_data, only : ntot_amode, numptr_amode
599 ! ---------------------- !
600 ! Input-Output Variables !
601 ! ---------------------- !
603 integer , intent(in) :: lchnk
604 integer , intent(in) :: mix
605 integer , intent(in) :: mkx
606 integer , intent(in) :: iend
607 integer , intent(in) :: ncnst
608 real(r8), intent(in) :: dt ! Time step : 2*delta_t [ s ]
609 real(r8), intent(in) :: ps0_in(mix,0:mkx) ! Environmental pressure at the interfaces [ Pa ]
610 real(r8), intent(in) :: zs0_in(mix,0:mkx) ! Environmental height at the interfaces [ m ]
611 real(r8), intent(in) :: p0_in(mix,mkx) ! Environmental pressure at the layer mid-point [ Pa ]
612 real(r8), intent(in) :: z0_in(mix,mkx) ! Environmental height at the layer mid-point [ m ]
613 real(r8), intent(in) :: dp0_in(mix,mkx) ! Environmental layer pressure thickness [ Pa ] > 0.
614 real(r8), intent(in) :: dpdry0_in(mix,mkx) ! Environmental dry layer pressure thickness [ Pa ]
615 real(r8), intent(in) :: u0_in(mix,mkx) ! Environmental zonal wind [ m/s ]
616 real(r8), intent(in) :: v0_in(mix,mkx) ! Environmental meridional wind [ m/s ]
617 real(r8), intent(in) :: qv0_in(mix,mkx) ! Environmental water vapor specific humidity [ kg/kg ]
618 real(r8), intent(in) :: ql0_in(mix,mkx) ! Environmental liquid water specific humidity [ kg/kg ]
619 real(r8), intent(in) :: qi0_in(mix,mkx) ! Environmental ice specific humidity [ kg/kg ]
620 real(r8), intent(in) :: t0_in(mix,mkx) ! Environmental temperature [ K ]
621 real(r8), intent(in) :: s0_in(mix,mkx) ! Environmental dry static energy [ J/kg ]
622 real(r8), intent(in) :: tr0_in(mix,mkx,ncnst) ! Environmental tracers [ #, kg/kg ]
623 real(r8), intent(in) :: tke_in(mix,0:mkx) ! Turbulent kinetic energy at the interfaces [ m2/s2 ]
624 real(r8), intent(in) :: cldfrct_in(mix,mkx) ! Total cloud fraction at the previous time step [ fraction ]
625 real(r8), intent(in) :: concldfrct_in(mix,mkx) ! Total convective cloud fraction at the previous time step [ fraction ]
626 real(r8), intent(in) :: pblh_in(mix) ! Height of PBL [ m ]
627 real(r8), intent(inout) :: cush_inout(mix) ! Convective scale height [ m ]
629 real(r8) tw0_in(mix,mkx) ! Wet bulb temperature [ K ]
630 real(r8) qw0_in(mix,mkx) ! Wet-bulb specific humidity [ kg/kg ]
632 real(r8), intent(out) :: umf_out(mix,0:mkx) ! Updraft mass flux at the interfaces [ kg/m2/s ]
633 real(r8), intent(out) :: qvten_out(mix,mkx) ! Tendency of water vapor specific humidity [ kg/kg/s ]
634 real(r8), intent(out) :: qlten_out(mix,mkx) ! Tendency of liquid water specific humidity [ kg/kg/s ]
635 real(r8), intent(out) :: qiten_out(mix,mkx) ! Tendency of ice specific humidity [ kg/kg/s ]
636 real(r8), intent(out) :: sten_out(mix,mkx) ! Tendency of dry static energy [ J/kg/s ]
637 real(r8), intent(out) :: uten_out(mix,mkx) ! Tendency of zonal wind [ m/s2 ]
638 real(r8), intent(out) :: vten_out(mix,mkx) ! Tendency of meridional wind [ m/s2 ]
639 real(r8), intent(out) :: trten_out(mix,mkx,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ]
640 real(r8), intent(out) :: qrten_out(mix,mkx) ! Tendency of rain water specific humidity [ kg/kg/s ]
641 real(r8), intent(out) :: qsten_out(mix,mkx) ! Tendency of snow specific humidity [ kg/kg/s ]
642 real(r8), intent(out) :: precip_out(mix) ! Precipitation ( rain + snow ) rate at surface [ m/s ]
643 real(r8), intent(out) :: snow_out(mix) ! Snow rate at surface [ m/s ]
644 real(r8), intent(out) :: evapc_out(mix,mkx) ! Tendency of evaporation of precipitation [ kg/kg/s ]
645 real(r8), intent(out) :: slflx_out(mix,0:mkx) ! Updraft/pen.entrainment liquid static energy flux [ J/kg * kg/m2/s ]
646 real(r8), intent(out) :: qtflx_out(mix,0:mkx) ! updraft/pen.entrainment total water flux [ kg/kg * kg/m2/s ]
647 real(r8), intent(out) :: flxprc1_out(mix,0:mkx) ! precip (rain+snow) flux
648 real(r8), intent(out) :: flxsnow1_out(mix,0:mkx) ! snow flux
649 real(r8), intent(out) :: cufrc_out(mix,mkx) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ]
650 real(r8), intent(out) :: qcu_out(mix,mkx) ! Condensate water specific humidity within cumulus updraft [ kg/kg ]
651 real(r8), intent(out) :: qlu_out(mix,mkx) ! Liquid water specific humidity within cumulus updraft [ kg/kg ]
652 real(r8), intent(out) :: qiu_out(mix,mkx) ! Ice specific humidity within cumulus updraft [ kg/kg ]
653 real(r8), intent(out) :: cbmf_out(mix) ! Cloud base mass flux [ kg/m2/s ]
654 real(r8), intent(out) :: qc_out(mix,mkx) ! Tendency of detrained cumulus condensate into the environment [ kg/kg/s ]
655 real(r8), intent(out) :: rliq_out(mix) ! Vertical integral of qc_out [ m/s ]
656 real(r8), intent(out) :: cnt_out(mix) ! Cumulus top interface index, cnt = kpen [ no ]
657 real(r8), intent(out) :: cnb_out(mix) ! Cumulus base interface index, cnb = krel - 1 [ no ]
660 ! Internal Output Variables
663 integer , external :: qsat
664 real(r8) qtten_out(mix,mkx) ! Tendency of qt [ kg/kg/s ]
665 real(r8) slten_out(mix,mkx) ! Tendency of sl [ J/kg/s ]
666 real(r8) ufrc_out(mix,0:mkx) ! Updraft fractional area at the interfaces [ fraction ]
667 real(r8) uflx_out(mix,0:mkx) ! Updraft/pen.entrainment zonal momentum flux [ m/s/m2/s ]
668 real(r8) vflx_out(mix,0:mkx) ! Updraft/pen.entrainment meridional momentum flux [ m/s/m2/s ]
669 real(r8) fer_out(mix,mkx) ! Fractional lateral entrainment rate [ 1/Pa ]
670 real(r8) fdr_out(mix,mkx) ! Fractional lateral detrainment rate [ 1/Pa ]
671 real(r8) cinh_out(mix) ! Convective INhibition upto LFC (CIN) [ J/kg ]
672 real(r8) trflx_out(mix,0:mkx,ncnst) ! Updraft/pen.entrainment tracer flux [ #/m2/s, kg/kg/m2/s ]
674 ! -------------------------------------------- !
675 ! One-dimensional variables at each grid point !
676 ! -------------------------------------------- !
680 real(r8) ps0(0:mkx) ! Environmental pressure at the interfaces [ Pa ]
681 real(r8) zs0(0:mkx) ! Environmental height at the interfaces [ m ]
682 real(r8) p0(mkx) ! Environmental pressure at the layer mid-point [ Pa ]
683 real(r8) z0(mkx) ! Environmental height at the layer mid-point [ m ]
684 real(r8) dp0(mkx) ! Environmental layer pressure thickness [ Pa ] > 0.
685 real(r8) dpdry0(mkx) ! Environmental dry layer pressure thickness [ Pa ]
686 real(r8) u0(mkx) ! Environmental zonal wind [ m/s ]
687 real(r8) v0(mkx) ! Environmental meridional wind [ m/s ]
688 real(r8) tke(0:mkx) ! Turbulent kinetic energy at the interfaces [ m2/s2 ]
689 real(r8) cldfrct(mkx) ! Total cloud fraction at the previous time step [ fraction ]
690 real(r8) concldfrct(mkx) ! Total convective cloud fraction at the previous time step [ fraction ]
691 real(r8) qv0(mkx) ! Environmental water vapor specific humidity [ kg/kg ]
692 real(r8) ql0(mkx) ! Environmental liquid water specific humidity [ kg/kg ]
693 real(r8) qi0(mkx) ! Environmental ice specific humidity [ kg/kg ]
694 real(r8) t0(mkx) ! Environmental temperature [ K ]
695 real(r8) s0(mkx) ! Environmental dry static energy [ J/kg ]
696 real(r8) pblh ! Height of PBL [ m ]
697 real(r8) cush ! Convective scale height [ m ]
698 real(r8) tr0(mkx,ncnst) ! Environmental tracers [ #, kg/kg ]
700 ! 2. Environmental variables directly derived from the input variables
702 real(r8) qt0(mkx) ! Environmental total specific humidity [ kg/kg ]
703 real(r8) thl0(mkx) ! Environmental liquid potential temperature [ K ]
704 real(r8) thvl0(mkx) ! Environmental liquid virtual potential temperature [ K ]
705 real(r8) ssqt0(mkx) ! Linear internal slope of environmental total specific humidity [ kg/kg/Pa ]
706 real(r8) ssthl0(mkx) ! Linear internal slope of environmental liquid potential temperature [ K/Pa ]
707 real(r8) ssu0(mkx) ! Linear internal slope of environmental zonal wind [ m/s/Pa ]
708 real(r8) ssv0(mkx) ! Linear internal slope of environmental meridional wind [ m/s/Pa ]
709 real(r8) thv0bot(mkx) ! Environmental virtual potential temperature at the bottom of each layer [ K ]
710 real(r8) thv0top(mkx) ! Environmental virtual potential temperature at the top of each layer [ K ]
711 real(r8) thvl0bot(mkx) ! Environmental liquid virtual potential temperature at the bottom of each layer [ K ]
712 real(r8) thvl0top(mkx) ! Environmental liquid virtual potential temperature at the top of each layer [ K ]
713 real(r8) exn0(mkx) ! Exner function at the layer mid points [ no ]
714 real(r8) exns0(0:mkx) ! Exner function at the interfaces [ no ]
715 real(r8) sstr0(mkx,ncnst) ! Linear slope of environmental tracers [ #/Pa, kg/kg/Pa ]
717 ! 2-1. For preventing negative condensate at the provisional time step
719 real(r8) qv0_star(mkx) ! Environmental water vapor specific humidity [ kg/kg ]
720 real(r8) ql0_star(mkx) ! Environmental liquid water specific humidity [ kg/kg ]
721 real(r8) qi0_star(mkx) ! Environmental ice specific humidity [ kg/kg ]
722 real(r8) t0_star(mkx) ! Environmental temperature [ K ]
723 real(r8) s0_star(mkx) ! Environmental dry static energy [ J/kg ]
725 ! 3. Variables associated with cumulus convection
727 real(r8) umf(0:mkx) ! Updraft mass flux at the interfaces [ kg/m2/s ]
728 real(r8) emf(0:mkx) ! Penetrative entrainment mass flux at the interfaces [ kg/m2/s ]
729 real(r8) qvten(mkx) ! Tendency of water vapor specific humidity [ kg/kg/s ]
730 real(r8) qlten(mkx) ! Tendency of liquid water specific humidity [ kg/kg/s ]
731 real(r8) qiten(mkx) ! Tendency of ice specific humidity [ kg/kg/s ]
732 real(r8) sten(mkx) ! Tendency of dry static energy [ J/kg ]
733 real(r8) uten(mkx) ! Tendency of zonal wind [ m/s2 ]
734 real(r8) vten(mkx) ! Tendency of meridional wind [ m/s2 ]
735 real(r8) qrten(mkx) ! Tendency of rain water specific humidity [ kg/kg/s ]
736 real(r8) qsten(mkx) ! Tendency of snow specific humidity [ kg/kg/s ]
737 real(r8) precip ! Precipitation rate ( rain + snow) at the surface [ m/s ]
738 real(r8) snow ! Snow rate at the surface [ m/s ]
739 real(r8) evapc(mkx) ! Tendency of evaporation of precipitation [ kg/kg/s ]
740 real(r8) slflx(0:mkx) ! Updraft/pen.entrainment liquid static energy flux [ J/kg * kg/m2/s ]
741 real(r8) qtflx(0:mkx) ! Updraft/pen.entrainment total water flux [ kg/kg * kg/m2/s ]
742 real(r8) uflx(0:mkx) ! Updraft/pen.entrainment flux of zonal momentum [ m/s/m2/s ]
743 real(r8) vflx(0:mkx) ! Updraft/pen.entrainment flux of meridional momentum [ m/s/m2/s ]
744 real(r8) cufrc(mkx) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ]
745 real(r8) qcu(mkx) ! Condensate water specific humidity within convective updraft [ kg/kg ]
746 real(r8) qlu(mkx) ! Liquid water specific humidity within convective updraft [ kg/kg ]
747 real(r8) qiu(mkx) ! Ice specific humidity within convective updraft [ kg/kg ]
748 real(r8) dwten(mkx) ! Detrained water tendency from cumulus updraft [ kg/kg/s ]
749 real(r8) diten(mkx) ! Detrained ice tendency from cumulus updraft [ kg/kg/s ]
750 real(r8) fer(mkx) ! Fractional lateral entrainment rate [ 1/Pa ]
751 real(r8) fdr(mkx) ! Fractional lateral detrainment rate [ 1/Pa ]
752 real(r8) uf(mkx) ! Zonal wind at the provisional time step [ m/s ]
753 real(r8) vf(mkx) ! Meridional wind at the provisional time step [ m/s ]
754 real(r8) qc(mkx) ! Tendency due to detrained 'cloud water + cloud ice' (without rain-snow contribution) [ kg/kg/s ]
755 real(r8) qc_l(mkx) ! Tendency due to detrained 'cloud water' (without rain-snow contribution) [ kg/kg/s ]
756 real(r8) qc_i(mkx) ! Tendency due to detrained 'cloud ice' (without rain-snow contribution) [ kg/kg/s ]
767 real(r8) rliq ! Vertical integral of qc [ m/s ]
768 real(r8) cnt ! Cumulus top interface index, cnt = kpen [ no ]
769 real(r8) cnb ! Cumulus base interface index, cnb = krel - 1 [ no ]
770 real(r8) qtten(mkx) ! Tendency of qt [ kg/kg/s ]
771 real(r8) slten(mkx) ! Tendency of sl [ J/kg/s ]
772 real(r8) ufrc(0:mkx) ! Updraft fractional area [ fraction ]
773 real(r8) trten(mkx,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ]
774 real(r8) trflx(0:mkx,ncnst) ! Flux of tracers due to convection [ # * kg/m2/s, kg/kg * kg/m2/s ]
775 real(r8) trflx_d(0:mkx) ! Adjustive downward flux of tracers to prevent negative tracers
776 real(r8) trflx_u(0:mkx) ! Adjustive upward flux of tracers to prevent negative tracers
777 real(r8) trmin ! Minimum concentration of tracers allowed
780 !----- Variables used for the calculation of condensation sink associated with compensating subsidence
781 ! In the current code, this 'sink' tendency is simply set to be zero.
783 real(r8) uemf(0:mkx) ! Net updraft mass flux at the interface ( emf + umf ) [ kg/m2/s ]
784 real(r8) comsub(mkx) ! Compensating subsidence at the layer mid-point ( unit of mass flux, umf ) [ kg/m2/s ]
785 real(r8) qlten_sink(mkx) ! Liquid condensate tendency by compensating subsidence/upwelling [ kg/kg/s ]
786 real(r8) qiten_sink(mkx) ! Ice condensate tendency by compensating subsidence/upwelling [ kg/kg/s ]
787 real(r8) nlten_sink(mkx) ! Liquid droplets # tendency by compensating subsidence/upwelling [ kg/kg/s ]
788 real(r8) niten_sink(mkx) ! Ice droplets # tendency by compensating subsidence/upwelling [ kg/kg/s ]
789 real(r8) thlten_sub, qtten_sub ! Tendency of conservative scalars by compensating subsidence/upwelling
790 real(r8) qlten_sub, qiten_sub ! Tendency of ql0, qi0 by compensating subsidence/upwelling
791 real(r8) nlten_sub, niten_sub ! Tendency of nl0, ni0 by compensating subsidence/upwelling
792 real(r8) thl_prog, qt_prog ! Prognosed 'thl, qt' by compensating subsidence/upwelling
794 !----- Variables describing cumulus updraft
796 real(r8) wu(0:mkx) ! Updraft vertical velocity at the interface [ m/s ]
797 real(r8) thlu(0:mkx) ! Updraft liquid potential temperature at the interface [ K ]
798 real(r8) qtu(0:mkx) ! Updraft total specific humidity at the interface [ kg/kg ]
799 real(r8) uu(0:mkx) ! Updraft zonal wind at the interface [ m/s ]
800 real(r8) vu(0:mkx) ! Updraft meridional wind at the interface [ m/s ]
801 real(r8) thvu(0:mkx) ! Updraft virtual potential temperature at the interface [ m/s ]
802 real(r8) rei(mkx) ! Updraft fractional mixing rate with the environment [ 1/Pa ]
803 real(r8) tru(0:mkx,ncnst) ! Updraft tracers [ #, kg/kg ]
805 !----- Variables describing conservative scalars of entraining downdrafts at the
806 ! entraining interfaces, i.e., 'kbup <= k < kpen-1'. At the other interfaces,
807 ! belows are simply set to equal to those of updraft for simplicity - but it
808 ! does not influence numerical calculation.
810 real(r8) thlu_emf(0:mkx) ! Penetrative downdraft liquid potential temperature at entraining interfaces [ K ]
811 real(r8) qtu_emf(0:mkx) ! Penetrative downdraft total water at entraining interfaces [ kg/kg ]
812 real(r8) uu_emf(0:mkx) ! Penetrative downdraft zonal wind at entraining interfaces [ m/s ]
813 real(r8) vu_emf(0:mkx) ! Penetrative downdraft meridional wind at entraining interfaces [ m/s ]
814 real(r8) tru_emf(0:mkx,ncnst) ! Penetrative Downdraft tracers at entraining interfaces [ #, kg/kg ]
816 !----- Variables associated with evaporations of convective 'rain' and 'snow'
818 real(r8) flxrain(0:mkx) ! Downward rain flux at each interface [ kg/m2/s ]
819 real(r8) flxsnow(0:mkx) ! Downward snow flux at each interface [ kg/m2/s ]
820 real(r8) ntraprd(mkx) ! Net production ( production - evaporation + melting ) rate of rain in each layer [ kg/kg/s ]
821 real(r8) ntsnprd(mkx) ! Net production ( production - evaporation + freezing ) rate of snow in each layer [ kg/kg/s ]
822 real(r8) flxsntm ! Downward snow flux at the top of each layer after melting [ kg/m2/s ]
823 real(r8) snowmlt ! Snow melting tendency [ kg/kg/s ]
824 real(r8) subsat ! Sub-saturation ratio (1-qv/qs) [ no unit ]
825 real(r8) evprain ! Evaporation rate of rain [ kg/kg/s ]
826 real(r8) evpsnow ! Evaporation rate of snow [ kg/kg/s ]
827 real(r8) evplimit ! Limiter of 'evprain + evpsnow' [ kg/kg/s ]
828 real(r8) evplimit_rain ! Limiter of 'evprain' [ kg/kg/s ]
829 real(r8) evplimit_snow ! Limiter of 'evpsnow' [ kg/kg/s ]
830 real(r8) evpint_rain ! Vertically-integrated evaporative flux of rain [ kg/m2/s ]
831 real(r8) evpint_snow ! Vertically-integrated evaporative flux of snow [ kg/m2/s ]
832 real(r8) kevp ! Evaporative efficiency [ complex unit ]
834 !----- Other internal variables
836 integer kk, mm, k, i, m, kp1, km1
837 integer iter_scaleh, iter_xc
838 integer id_check, status
839 integer klcl ! Layer containing LCL of source air
840 integer kinv ! Inversion layer with PBL top interface as a lower interface
841 integer krel ! Release layer where buoyancy sorting mixing occurs for the first time
842 integer klfc ! LFC layer of cumulus source air
843 integer kbup ! Top layer in which cloud buoyancy is positive at the top interface
844 integer kpen ! Highest layer with positive updraft vertical velocity - top layer cumulus can reach
846 logical forcedCu ! If 'true', cumulus updraft cannot overcome the buoyancy barrier just above the PBL top.
847 real(r8) thlsrc, qtsrc, usrc, vsrc, thvlsrc ! Updraft source air properties
848 real(r8) PGFc, uplus, vplus
849 real(r8) trsrc(ncnst), tre(ncnst)
850 real(r8) plcl, plfc, prel, wrel
852 real(r8) ee2, ud2, wtw, wtwb, wtwh
854 real(r8) cldhgt, scaleh, tscaleh, cridis, rle, rkm
855 real(r8) rkfre, sigmaw, epsvarw, tkeavg, dpsum, dpi, thvlmin
856 real(r8) thlxsat, qtxsat, thvxsat, x_cu, x_en, thv_x0, thv_x1
857 real(r8) thj, qvj, qlj, qij, thvj, tj, thv0j, rho0j, rhos0j, qse
859 real(r8) pe, dpe, exne, thvebot, thle, qte, ue, ve, thlue, qtue, wue
860 real(r8) mu, mumin0, mumin1, mumin2, mulcl, mulclstar
861 real(r8) cbmf, wcrit, winv, wlcl, ufrcinv, ufrclcl, rmaxfrac
862 real(r8) criqc, exql, exqi, ppen
863 real(r8) thl0top, thl0bot, qt0bot, qt0top, thvubot, thvutop
864 real(r8) thlu_top, qtu_top, qlu_top, qiu_top, qlu_mid, qiu_mid, exntop
865 real(r8) thl0lcl, qt0lcl, thv0lcl, thv0rel, rho0inv, autodet
866 real(r8) aquad, bquad, cquad, xc1, xc2, excessu, excess0, xsat, xs1, xs2
867 real(r8) bogbot, bogtop, delbog, drage, expfac, rbuoy, rdrag
868 real(r8) rcwp, rlwp, riwp, qcubelow, qlubelow, qiubelow
869 real(r8) rainflx, snowflx
872 real(r8) gam(1) ! (L/cp)*dqs/dT
874 real(r8) xsrc, xmean, xtop, xbot, xflx(0:mkx)
877 !----- Some diagnostic internal output variables
879 real(r8) ufrcinvbase_out(mix) ! Cumulus updraft fraction at the PBL top [ fraction ]
880 real(r8) ufrclcl_out(mix) ! Cumulus updraft fraction at the LCL ( or PBL top when LCL is below PBL top ) [ fraction ]
881 real(r8) winvbase_out(mix) ! Cumulus updraft velocity at the PBL top [ m/s ]
882 real(r8) wlcl_out(mix) ! Cumulus updraft velocity at the LCL ( or PBL top when LCL is below PBL top ) [ m/s ]
883 real(r8) plcl_out(mix) ! LCL of source air [ Pa ]
884 real(r8) pinv_out(mix) ! PBL top pressure [ Pa ]
885 real(r8) plfc_out(mix) ! LFC of source air [ Pa ]
886 real(r8) pbup_out(mix) ! Highest interface level of positive buoyancy [ Pa ]
887 real(r8) ppen_out(mix) ! Highest interface evel where Cu w = 0 [ Pa ]
888 real(r8) qtsrc_out(mix) ! Sourse air qt [ kg/kg ]
889 real(r8) thlsrc_out(mix) ! Sourse air thl [ K ]
890 real(r8) thvlsrc_out(mix) ! Sourse air thvl [ K ]
891 real(r8) emfkbup_out(mix) ! Penetrative downward mass flux at 'kbup' interface [ kg/m2/s ]
892 real(r8) cinlclh_out(mix) ! Convective INhibition upto LCL (CIN) [ J/kg = m2/s2 ]
893 real(r8) tkeavg_out(mix) ! Average tke over the PBL [ m2/s2 ]
894 real(r8) cbmflimit_out(mix) ! Cloud base mass flux limiter [ kg/m2/s ]
895 real(r8) zinv_out(mix) ! PBL top height [ m ]
896 real(r8) rcwp_out(mix) ! Layer mean Cumulus LWP+IWP [ kg/m2 ]
897 real(r8) rlwp_out(mix) ! Layer mean Cumulus LWP [ kg/m2 ]
898 real(r8) riwp_out(mix) ! Layer mean Cumulus IWP [ kg/m2 ]
899 real(r8) wu_out(mix,0:mkx) ! Updraft vertical velocity ( defined from the release level to 'kpen-1' interface )
900 real(r8) qtu_out(mix,0:mkx) ! Updraft qt [ kg/kg ]
901 real(r8) thlu_out(mix,0:mkx) ! Updraft thl [ K ]
902 real(r8) thvu_out(mix,0:mkx) ! Updraft thv [ K ]
903 real(r8) uu_out(mix,0:mkx) ! Updraft zonal wind [ m/s ]
904 real(r8) vu_out(mix,0:mkx) ! Updraft meridional wind [ m/s ]
905 real(r8) qtu_emf_out(mix,0:mkx) ! Penetratively entrained qt [ kg/kg ]
906 real(r8) thlu_emf_out(mix,0:mkx) ! Penetratively entrained thl [ K ]
907 real(r8) uu_emf_out(mix,0:mkx) ! Penetratively entrained u [ m/s ]
908 real(r8) vu_emf_out(mix,0:mkx) ! Penetratively entrained v [ m/s ]
909 real(r8) uemf_out(mix,0:mkx) ! Net upward mass flux including penetrative entrainment (umf+emf) [ kg/m2/s ]
910 real(r8) tru_out(mix,0:mkx,ncnst) ! Updraft tracers [ #, kg/kg ]
911 real(r8) tru_emf_out(mix,0:mkx,ncnst) ! Penetratively entrained tracers [ #, kg/kg ]
913 real(r8) wu_s(0:mkx) ! Same as above but for implicit CIN
914 real(r8) qtu_s(0:mkx)
915 real(r8) thlu_s(0:mkx)
916 real(r8) thvu_s(0:mkx)
919 real(r8) qtu_emf_s(0:mkx)
920 real(r8) thlu_emf_s(0:mkx)
921 real(r8) uu_emf_s(0:mkx)
922 real(r8) vu_emf_s(0:mkx)
923 real(r8) uemf_s(0:mkx)
924 real(r8) tru_s(0:mkx,ncnst)
925 real(r8) tru_emf_s(0:mkx,ncnst)
927 real(r8) dwten_out(mix,mkx)
928 real(r8) diten_out(mix,mkx)
929 real(r8) flxrain_out(mix,0:mkx)
930 real(r8) flxsnow_out(mix,0:mkx)
931 real(r8) ntraprd_out(mix,mkx)
932 real(r8) ntsnprd_out(mix,mkx)
934 real(r8) dwten_s(mkx)
935 real(r8) diten_s(mkx)
936 real(r8) flxrain_s(0:mkx)
937 real(r8) flxsnow_s(0:mkx)
938 real(r8) ntraprd_s(mkx)
939 real(r8) ntsnprd_s(mkx)
941 real(r8) excessu_arr_out(mix,mkx)
942 real(r8) excessu_arr(mkx)
943 real(r8) excessu_arr_s(mkx)
944 real(r8) excess0_arr_out(mix,mkx)
945 real(r8) excess0_arr(mkx)
946 real(r8) excess0_arr_s(mkx)
947 real(r8) xc_arr_out(mix,mkx)
949 real(r8) xc_arr_s(mkx)
950 real(r8) aquad_arr_out(mix,mkx)
951 real(r8) aquad_arr(mkx)
952 real(r8) aquad_arr_s(mkx)
953 real(r8) bquad_arr_out(mix,mkx)
954 real(r8) bquad_arr(mkx)
955 real(r8) bquad_arr_s(mkx)
956 real(r8) cquad_arr_out(mix,mkx)
957 real(r8) cquad_arr(mkx)
958 real(r8) cquad_arr_s(mkx)
959 real(r8) bogbot_arr_out(mix,mkx)
960 real(r8) bogbot_arr(mkx)
961 real(r8) bogbot_arr_s(mkx)
962 real(r8) bogtop_arr_out(mix,mkx)
963 real(r8) bogtop_arr(mkx)
964 real(r8) bogtop_arr_s(mkx)
966 real(r8) exit_UWCu(mix)
967 real(r8) exit_conden(mix)
968 real(r8) exit_klclmkx(mix)
969 real(r8) exit_klfcmkx(mix)
970 real(r8) exit_ufrc(mix)
971 real(r8) exit_wtw(mix)
972 real(r8) exit_drycore(mix)
973 real(r8) exit_wu(mix)
974 real(r8) exit_cufilter(mix)
975 real(r8) exit_kinv1(mix)
976 real(r8) exit_rei(mix)
978 real(r8) limit_shcu(mix)
979 real(r8) limit_negcon(mix)
980 real(r8) limit_ufrc(mix)
981 real(r8) limit_ppen(mix)
982 real(r8) limit_emf(mix)
983 real(r8) limit_cinlcl(mix)
984 real(r8) limit_cin(mix)
985 real(r8) limit_cbmf(mix)
986 real(r8) limit_rei(mix)
987 real(r8) ind_delcin(mix)
989 real(r8) :: ufrcinvbase_s, ufrclcl_s, winvbase_s, wlcl_s, plcl_s, pinv_s, plfc_s, &
990 qtsrc_s, thlsrc_s, thvlsrc_s, emfkbup_s, cinlcl_s, pbup_s, ppen_s, cbmflimit_s, &
991 tkeavg_s, zinv_s, rcwp_s, rlwp_s, riwp_s
992 real(r8) :: ufrcinvbase, winvbase, pinv, zinv, emfkbup, cbmflimit, rho0rel
994 !----- Variables for implicit CIN computation
996 real(r8), dimension(mkx) :: qv0_s , ql0_s , qi0_s , s0_s , u0_s , &
997 v0_s , t0_s , qt0_s , thl0_s , thvl0_s , qvten_s , &
998 qlten_s, qiten_s , qrten_s , qsten_s , sten_s , evapc_s , &
999 uten_s , vten_s , cufrc_s , qcu_s , qlu_s , qiu_s , &
1000 fer_s , fdr_s , qc_s , qtten_s , slten_s
1001 real(r8), dimension(0:mkx) :: umf_s , slflx_s , qtflx_s , ufrc_s , uflx_s , vflx_s
1002 real(r8) :: cush_s , precip_s, snow_s , cin_s , rliq_s, cbmf_s, cnt_s, cnb_s
1003 real(r8) :: cin_i,cin_f,del_CIN,ke,alpha,thlj
1004 real(r8) :: cinlcl_i,cinlcl_f,del_cinlcl
1007 real(r8), dimension(mkx,ncnst) :: tr0_s, trten_s
1008 real(r8), dimension(0:mkx,ncnst) :: trflx_s
1010 !----- Variables for temporary storages
1012 real(r8), dimension(mkx) :: qv0_o, ql0_o, qi0_o, t0_o, s0_o, u0_o, v0_o
1013 real(r8), dimension(mkx) :: qt0_o , thl0_o , thvl0_o , &
1014 qvten_o , qlten_o , qiten_o , qrten_o , qsten_o , &
1015 sten_o , uten_o , vten_o , qcu_o , qlu_o , &
1016 qiu_o , cufrc_o , evapc_o , &
1017 thv0bot_o, thv0top_o, thvl0bot_o, thvl0top_o, &
1018 ssthl0_o , ssqt0_o , ssu0_o , ssv0_o , qc_o , &
1020 real(r8), dimension(0:mkx) :: umf_o , slflx_o , qtflx_o , ufrc_o
1021 real(r8), dimension(mix) :: cush_o , precip_o , snow_o , rliq_o, cbmf_o, cnt_o, cnb_o
1022 real(r8), dimension(0:mkx) :: uflx_o , vflx_o
1023 real(r8) :: tkeavg_o , thvlmin_o, qtsrc_o , thvlsrc_o, thlsrc_o , &
1024 usrc_o , vsrc_o , plcl_o , plfc_o , &
1026 integer :: kinv_o , klcl_o , klfc_o
1028 real(r8), dimension(mkx,ncnst) :: tr0_o
1029 real(r8), dimension(mkx,ncnst) :: trten_o, sstr0_o
1030 real(r8), dimension(0:mkx,ncnst) :: trflx_o
1031 real(r8), dimension(ncnst) :: trsrc_o
1032 integer :: ixnumliq, ixnumice
1034 ! ------------------ !
1036 ! Define Parameters !
1038 ! ------------------ !
1040 ! ------------------------ !
1041 ! Iterative xc calculation !
1042 ! ------------------------ !
1044 integer , parameter :: niter_xc = 2
1046 ! ----------------------------------------------------------- !
1047 ! Choice of 'CIN = cin' (.true.) or 'CIN = cinlcl' (.false.). !
1048 ! ----------------------------------------------------------- !
1050 logical , parameter :: use_CINcin = .true.
1052 ! --------------------------------------------------------------- !
1053 ! Choice of 'explicit' ( 1 ) or 'implicit' ( 2 ) CIN. !
1055 ! When choose 'CIN = cinlcl' above, it is recommended not to use !
1056 ! implicit CIN, i.e., do 'NOT' choose simultaneously : !
1057 ! [ 'use_CINcin=.false. & 'iter_cin=2' ] !
1058 ! since 'cinlcl' will be always set to zero whenever LCL is below !
1059 ! the PBL top interface in the current code. So, averaging cinlcl !
1060 ! of two iter_cin steps is likely not so good. Except that, all !
1061 ! the other combinations of 'use_CINcin' & 'iter_cin' are OK. !
1063 ! Feb 2007, Bundy: Note that use_CINcin = .false. will try to use !
1064 ! a variable (del_cinlcl) that is not currently set !
1066 ! --------------------------------------------------------------- !
1068 integer , parameter :: iter_cin = 2
1070 ! ---------------------------------------------------------------- !
1071 ! Choice of 'self-detrainment' by negative buoyancy in calculating !
1072 ! cumulus updraft mass flux at the top interface in each layer. !
1073 ! ---------------------------------------------------------------- !
1075 logical , parameter :: use_self_detrain = .false.
1077 ! --------------------------------------------------------- !
1078 ! Cumulus momentum flux : turn-on (.true.) or off (.false.) !
1079 ! --------------------------------------------------------- !
1081 logical , parameter :: use_momenflx = .true.
1083 ! ----------------------------------------------------------------------------------------- !
1084 ! Penetrative Entrainment : Cumulative ( .true. , original ) or Non-Cumulative ( .false. ) !
1085 ! This option ( .false. ) is designed to reduce the sensitivity to the vertical resolution. !
1086 ! ----------------------------------------------------------------------------------------- !
1088 logical , parameter :: use_cumpenent = .true.
1090 ! --------------------------------------------------------------------------------------------------------------- !
1091 ! Computation of the grid-mean condensate tendency. !
1092 ! use_expconten = .true. : explcitly compute tendency by condensate detrainment and compensating subsidence !
1093 ! use_expconten = .false. : use the original proportional condensate tendency equation. ( original ) !
1094 ! --------------------------------------------------------------------------------------------------------------- !
1096 logical , parameter :: use_expconten = .true.
1098 ! --------------------------------------------------------------------------------------------------------------- !
1099 ! Treatment of reserved condensate !
1100 ! use_unicondet = .true. : detrain condensate uniformly over the environment ( original ) !
1101 ! use_unicondet = .false. : detrain condensate into the pre-existing stratus !
1102 ! --------------------------------------------------------------------------------------------------------------- !
1104 logical , parameter :: use_unicondet = .false.
1106 ! ----------------------- !
1107 ! For lateral entrainment !
1108 ! ----------------------- !
1110 parameter (rle = 0.1_r8) ! For critical stopping distance for lateral entrainment [no unit]
1111 ! parameter (rkm = 16.0_r8) ! Determine the amount of air that is involved in buoyancy-sorting [no unit]
1112 parameter (rkm = 14.0_r8) ! Determine the amount of air that is involved in buoyancy-sorting [no unit]
1114 parameter (rkfre = 1.0_r8) ! Vertical velocity variance as fraction of tke.
1115 parameter (rmaxfrac = 0.10_r8) ! Maximum allowable 'core' updraft fraction
1116 parameter (mumin1 = 0.906_r8) ! Normalized CIN ('mu') corresponding to 'rmaxfrac' at the PBL top
1117 ! obtaind by inverting 'rmaxfrac = 0.5*erfc(mumin1)'.
1118 ! [ rmaxfrac:mumin1 ] = [ 0.05:1.163, 0.075:1.018, 0.1:0.906, 0.15:0.733, 0.2:0.595, 0.25:0.477 ]
1119 parameter (rbuoy = 1.0_r8) ! For nonhydrostatic pressure effects on updraft [no unit]
1120 parameter (rdrag = 1.0_r8) ! Drag coefficient [no unit]
1122 parameter (epsvarw = 5.e-4_r8) ! Variance of w at PBL top by meso-scale component [m2/s2]
1123 parameter (PGFc = 0.7_r8) ! This is used for calculating vertical variations cumulus
1124 ! 'u' & 'v' by horizontal PGF during upward motion [no unit]
1126 ! ---------------------------------------- !
1127 ! Bulk microphysics controlling parameters !
1128 ! --------------------------------------------------------------------------- !
1129 ! criqc : Maximum condensate that can be hold by cumulus updraft [kg/kg] !
1130 ! frc_rasn : Fraction of precipitable condensate in the expelled cloud water !
1131 ! from cumulus updraft. The remaining fraction ('1-frc_rasn') is !
1132 ! 'suspended condensate'. !
1133 ! 0 : all expelled condensate is 'suspended condensate' !
1134 ! 1 : all expelled condensate is 'precipitable condensate' !
1135 ! kevp : Evaporative efficiency !
1136 ! noevap_krelkpen : No evaporation from 'krel' to 'kpen' layers !
1137 ! --------------------------------------------------------------------------- !
1139 parameter ( criqc = 0.7e-3_r8 )
1140 parameter ( frc_rasn = 1.0_r8 )
1141 parameter ( kevp = 2.e-6_r8 )
1142 logical, parameter :: noevap_krelkpen = .false.
1144 !------------------------!
1146 ! Start Main Calculation !
1148 !------------------------!
1150 call cnst_get_ind( 'NUMLIQ', ixnumliq )
1151 call cnst_get_ind( 'NUMICE', ixnumice )
1153 ! ------------------------------------------------------- !
1154 ! Initialize output variables defined for all grid points !
1155 ! ------------------------------------------------------- !
1157 umf_out(:iend,0:mkx) = 0.0_r8
1158 slflx_out(:iend,0:mkx) = 0.0_r8
1159 qtflx_out(:iend,0:mkx) = 0.0_r8
1160 flxprc1_out(:iend,0:mkx) = 0.0_r8
1161 flxsnow1_out(:iend,0:mkx) = 0.0_r8
1162 qvten_out(:iend,:mkx) = 0.0_r8
1163 qlten_out(:iend,:mkx) = 0.0_r8
1164 qiten_out(:iend,:mkx) = 0.0_r8
1165 sten_out(:iend,:mkx) = 0.0_r8
1166 uten_out(:iend,:mkx) = 0.0_r8
1167 vten_out(:iend,:mkx) = 0.0_r8
1168 qrten_out(:iend,:mkx) = 0.0_r8
1169 qsten_out(:iend,:mkx) = 0.0_r8
1170 precip_out(:iend) = 0.0_r8
1171 snow_out(:iend) = 0.0_r8
1172 evapc_out(:iend,:mkx) = 0.0_r8
1173 cufrc_out(:iend,:mkx) = 0.0_r8
1174 qcu_out(:iend,:mkx) = 0.0_r8
1175 qlu_out(:iend,:mkx) = 0.0_r8
1176 qiu_out(:iend,:mkx) = 0.0_r8
1177 fer_out(:iend,:mkx) = 0.0_r8
1178 fdr_out(:iend,:mkx) = 0.0_r8
1179 cinh_out(:iend) = -1.0_r8
1180 cinlclh_out(:iend) = -1.0_r8
1181 cbmf_out(:iend) = 0.0_r8
1182 qc_out(:iend,:mkx) = 0.0_r8
1183 rliq_out(:iend) = 0.0_r8
1184 cnt_out(:iend) = real(mkx, r8)
1185 cnb_out(:iend) = 0.0_r8
1186 qtten_out(:iend,:mkx) = 0.0_r8
1187 slten_out(:iend,:mkx) = 0.0_r8
1188 ufrc_out(:iend,0:mkx) = 0.0_r8
1190 uflx_out(:iend,0:mkx) = 0.0_r8
1191 vflx_out(:iend,0:mkx) = 0.0_r8
1193 trten_out(:iend,:mkx,:ncnst) = 0.0_r8
1194 trflx_out(:iend,0:mkx,:ncnst)= 0.0_r8
1196 ufrcinvbase_out(:iend) = 0.0_r8
1197 ufrclcl_out(:iend) = 0.0_r8
1198 winvbase_out(:iend) = 0.0_r8
1199 wlcl_out(:iend) = 0.0_r8
1200 plcl_out(:iend) = 0.0_r8
1201 pinv_out(:iend) = 0.0_r8
1202 plfc_out(:iend) = 0.0_r8
1203 pbup_out(:iend) = 0.0_r8
1204 ppen_out(:iend) = 0.0_r8
1205 qtsrc_out(:iend) = 0.0_r8
1206 thlsrc_out(:iend) = 0.0_r8
1207 thvlsrc_out(:iend) = 0.0_r8
1208 emfkbup_out(:iend) = 0.0_r8
1209 cbmflimit_out(:iend) = 0.0_r8
1210 tkeavg_out(:iend) = 0.0_r8
1211 zinv_out(:iend) = 0.0_r8
1212 rcwp_out(:iend) = 0.0_r8
1213 rlwp_out(:iend) = 0.0_r8
1214 riwp_out(:iend) = 0.0_r8
1216 wu_out(:iend,0:mkx) = 0.0_r8
1217 qtu_out(:iend,0:mkx) = 0.0_r8
1218 thlu_out(:iend,0:mkx) = 0.0_r8
1219 thvu_out(:iend,0:mkx) = 0.0_r8
1220 uu_out(:iend,0:mkx) = 0.0_r8
1221 vu_out(:iend,0:mkx) = 0.0_r8
1222 qtu_emf_out(:iend,0:mkx) = 0.0_r8
1223 thlu_emf_out(:iend,0:mkx) = 0.0_r8
1224 uu_emf_out(:iend,0:mkx) = 0.0_r8
1225 vu_emf_out(:iend,0:mkx) = 0.0_r8
1226 uemf_out(:iend,0:mkx) = 0.0_r8
1228 tru_out(:iend,0:mkx,:ncnst) = 0.0_r8
1229 tru_emf_out(:iend,0:mkx,:ncnst) = 0.0_r8
1231 dwten_out(:iend,:mkx) = 0.0_r8
1232 diten_out(:iend,:mkx) = 0.0_r8
1233 flxrain_out(:iend,0:mkx) = 0.0_r8
1234 flxsnow_out(:iend,0:mkx) = 0.0_r8
1235 ntraprd_out(:iend,mkx) = 0.0_r8
1236 ntsnprd_out(:iend,mkx) = 0.0_r8
1238 excessu_arr_out(:iend,:mkx) = 0.0_r8
1239 excess0_arr_out(:iend,:mkx) = 0.0_r8
1240 xc_arr_out(:iend,:mkx) = 0.0_r8
1241 aquad_arr_out(:iend,:mkx) = 0.0_r8
1242 bquad_arr_out(:iend,:mkx) = 0.0_r8
1243 cquad_arr_out(:iend,:mkx) = 0.0_r8
1244 bogbot_arr_out(:iend,:mkx) = 0.0_r8
1245 bogtop_arr_out(:iend,:mkx) = 0.0_r8
1247 exit_UWCu(:iend) = 0.0_r8
1248 exit_conden(:iend) = 0.0_r8
1249 exit_klclmkx(:iend) = 0.0_r8
1250 exit_klfcmkx(:iend) = 0.0_r8
1251 exit_ufrc(:iend) = 0.0_r8
1252 exit_wtw(:iend) = 0.0_r8
1253 exit_drycore(:iend) = 0.0_r8
1254 exit_wu(:iend) = 0.0_r8
1255 exit_cufilter(:iend) = 0.0_r8
1256 exit_kinv1(:iend) = 0.0_r8
1257 exit_rei(:iend) = 0.0_r8
1259 limit_shcu(:iend) = 0.0_r8
1260 limit_negcon(:iend) = 0.0_r8
1261 limit_ufrc(:iend) = 0.0_r8
1262 limit_ppen(:iend) = 0.0_r8
1263 limit_emf(:iend) = 0.0_r8
1264 limit_cinlcl(:iend) = 0.0_r8
1265 limit_cin(:iend) = 0.0_r8
1266 limit_cbmf(:iend) = 0.0_r8
1267 limit_rei(:iend) = 0.0_r8
1269 ind_delcin(:iend) = 0.0_r8
1271 !--------------------------------------------------------------!
1273 ! Start the column i loop where i is a horozontal column index !
1275 !--------------------------------------------------------------!
1277 ! Compute wet-bulb temperature and specific humidity
1278 ! for treating evaporation of precipitation.
1280 call findsp( lchnk, iend, qv0_in, t0_in, p0_in, tw0_in, qw0_in )
1286 ! -------------------------------------------- !
1287 ! Define 1D input variables at each grid point !
1288 ! -------------------------------------------- !
1290 ps0(0:mkx) = ps0_in(i,0:mkx)
1291 zs0(0:mkx) = zs0_in(i,0:mkx)
1292 p0(:mkx) = p0_in(i,:mkx)
1293 z0(:mkx) = z0_in(i,:mkx)
1294 dp0(:mkx) = dp0_in(i,:mkx)
1295 dpdry0(:mkx) = dpdry0_in(i,:mkx)
1296 u0(:mkx) = u0_in(i,:mkx)
1297 v0(:mkx) = v0_in(i,:mkx)
1298 qv0(:mkx) = qv0_in(i,:mkx)
1299 ql0(:mkx) = ql0_in(i,:mkx)
1300 qi0(:mkx) = qi0_in(i,:mkx)
1301 t0(:mkx) = t0_in(i,:mkx)
1302 s0(:mkx) = s0_in(i,:mkx)
1303 tke(0:mkx) = tke_in(i,0:mkx)
1304 cldfrct(:mkx) = cldfrct_in(i,:mkx)
1305 concldfrct(:mkx) = concldfrct_in(i,:mkx)
1307 cush = cush_inout(i)
1309 tr0(:mkx,m) = tr0_in(i,:mkx,m)
1312 ! --------------------------------------------------------- !
1313 ! Compute other basic thermodynamic variables directly from !
1314 ! the input variables at each grid point !
1315 ! --------------------------------------------------------- !
1317 !----- 1. Compute internal environmental variables
1319 exn0(:mkx) = (p0(:mkx)/p00)**rovcp
1320 exns0(0:mkx) = (ps0(0:mkx)/p00)**rovcp
1321 qt0(:mkx) = (qv0(:mkx) + ql0(:mkx) + qi0(:mkx))
1322 thl0(:mkx) = (t0(:mkx) - xlv*ql0(:mkx)/cp - xls*qi0(:mkx)/cp)/exn0(:mkx)
1323 thvl0(:mkx) = (1._r8 + zvir*qt0(:mkx))*thl0(:mkx)
1325 !----- 2. Compute slopes of environmental variables in each layer
1326 ! Dimension of ssthl0(:mkx) is implicit.
1328 ssthl0 = slope(mkx,thl0,p0)
1329 ssqt0 = slope(mkx,qt0 ,p0)
1330 ssu0 = slope(mkx,u0 ,p0)
1331 ssv0 = slope(mkx,v0 ,p0)
1333 sstr0(:mkx,m) = slope(mkx,tr0(:mkx,m),p0)
1336 !----- 3. Compute "thv0" and "thvl0" at the top/bottom interfaces in each layer
1337 ! There are computed from the reconstructed thl, qt at the top/bottom.
1341 thl0bot = thl0(k) + ssthl0(k)*(ps0(k-1) - p0(k))
1342 qt0bot = qt0(k) + ssqt0(k) *(ps0(k-1) - p0(k))
1343 call conden(ps0(k-1),thl0bot,qt0bot,thj,qvj,qlj,qij,qse,id_check,qsat)
1344 if( id_check .eq. 1 ) then
1345 exit_conden(i) = 1._r8
1349 thv0bot(k) = thj*(1._r8 + zvir*qvj - qlj - qij)
1350 thvl0bot(k) = thl0bot*(1._r8 + zvir*qt0bot)
1352 thl0top = thl0(k) + ssthl0(k)*(ps0(k) - p0(k))
1353 qt0top = qt0(k) + ssqt0(k) *(ps0(k) - p0(k))
1354 call conden(ps0(k),thl0top,qt0top,thj,qvj,qlj,qij,qse,id_check,qsat)
1355 if( id_check .eq. 1 ) then
1356 exit_conden(i) = 1._r8
1360 thv0top(k) = thj*(1._r8 + zvir*qvj - qlj - qij)
1361 thvl0top(k) = thl0top*(1._r8 + zvir*qt0top)
1365 ! ------------------------------------------------------------ !
1366 ! Save input and related environmental thermodynamic variables !
1367 ! for use at "iter_cin=2" when "del_CIN >= 0" !
1368 ! ------------------------------------------------------------ !
1370 qv0_o(:mkx) = qv0(:mkx)
1371 ql0_o(:mkx) = ql0(:mkx)
1372 qi0_o(:mkx) = qi0(:mkx)
1373 t0_o(:mkx) = t0(:mkx)
1374 s0_o(:mkx) = s0(:mkx)
1375 u0_o(:mkx) = u0(:mkx)
1376 v0_o(:mkx) = v0(:mkx)
1377 qt0_o(:mkx) = qt0(:mkx)
1378 thl0_o(:mkx) = thl0(:mkx)
1379 thvl0_o(:mkx) = thvl0(:mkx)
1380 ssthl0_o(:mkx) = ssthl0(:mkx)
1381 ssqt0_o(:mkx) = ssqt0(:mkx)
1382 thv0bot_o(:mkx) = thv0bot(:mkx)
1383 thv0top_o(:mkx) = thv0top(:mkx)
1384 thvl0bot_o(:mkx) = thvl0bot(:mkx)
1385 thvl0top_o(:mkx) = thvl0top(:mkx)
1386 ssu0_o(:mkx) = ssu0(:mkx)
1387 ssv0_o(:mkx) = ssv0(:mkx)
1389 tr0_o(:mkx,m) = tr0(:mkx,m)
1390 sstr0_o(:mkx,m) = sstr0(:mkx,m)
1393 ! ---------------------------------------------- !
1394 ! Initialize output variables at each grid point !
1395 ! ---------------------------------------------- !
1399 slflx(0:mkx) = 0.0_r8
1400 qtflx(0:mkx) = 0.0_r8
1401 uflx(0:mkx) = 0.0_r8
1402 vflx(0:mkx) = 0.0_r8
1403 qvten(:mkx) = 0.0_r8
1404 qlten(:mkx) = 0.0_r8
1405 qiten(:mkx) = 0.0_r8
1409 qrten(:mkx) = 0.0_r8
1410 qsten(:mkx) = 0.0_r8
1411 dwten(:mkx) = 0.0_r8
1412 diten(:mkx) = 0.0_r8
1415 evapc(:mkx) = 0.0_r8
1416 cufrc(:mkx) = 0.0_r8
1430 qtten(:mkx) = 0.0_r8
1431 slten(:mkx) = 0.0_r8
1432 ufrc(0:mkx) = 0.0_r8
1434 thlu(0:mkx) = 0.0_r8
1439 thvu(0:mkx) = 0.0_r8
1440 thlu_emf(0:mkx) = 0.0_r8
1441 qtu_emf(0:mkx) = 0.0_r8
1442 uu_emf(0:mkx) = 0.0_r8
1443 vu_emf(0:mkx) = 0.0_r8
1445 ufrcinvbase = 0.0_r8
1451 excessu_arr(:mkx) = 0.0_r8
1452 excess0_arr(:mkx) = 0.0_r8
1453 xc_arr(:mkx) = 0.0_r8
1454 aquad_arr(:mkx) = 0.0_r8
1455 bquad_arr(:mkx) = 0.0_r8
1456 cquad_arr(:mkx) = 0.0_r8
1457 bogbot_arr(:mkx) = 0.0_r8
1458 bogtop_arr(:mkx) = 0.0_r8
1460 uemf(0:mkx) = 0.0_r8
1461 comsub(:mkx) = 0.0_r8
1462 qlten_sink(:mkx) = 0.0_r8
1463 qiten_sink(:mkx) = 0.0_r8
1464 nlten_sink(:mkx) = 0.0_r8
1465 niten_sink(:mkx) = 0.0_r8
1468 trflx(0:mkx,m) = 0.0_r8
1469 trten(:mkx,m) = 0.0_r8
1470 tru(0:mkx,m) = 0.0_r8
1471 tru_emf(0:mkx,m) = 0.0_r8
1474 !-----------------------------------------------!
1475 ! Below 'iter' loop is for implicit CIN closure !
1476 !-----------------------------------------------!
1478 ! ----------------------------------------------------------------------------- !
1479 ! It is important to note that this iterative cin loop is located at the outest !
1480 ! shell of the code. Thus, source air properties can also be changed during the !
1481 ! iterative cin calculation, because cumulus convection induces non-zero fluxes !
1482 ! even at interfaces below PBL top height through 'fluxbelowinv' subroutine. !
1483 ! ----------------------------------------------------------------------------- !
1485 do iter = 1, iter_cin
1487 ! ---------------------------------------------------------------------- !
1488 ! Cumulus scale height !
1489 ! In contrast to the premitive code, cumulus scale height is iteratively !
1490 ! calculated at each time step, and at each iterative cin step. !
1491 ! It is not clear whether I should locate below two lines within or out !
1492 ! of the iterative cin loop. !
1493 ! ---------------------------------------------------------------------- !
1498 ! ----------------------------------------------------------------------- !
1499 ! Find PBL top height interface index, 'kinv-1' where 'kinv' is the layer !
1500 ! index with PBLH in it. When PBLH is exactly at interface, 'kinv' is the !
1501 ! layer index having PBLH as a lower interface. !
1502 ! In the previous code, I set the lower limit of 'kinv' by 2 in order to !
1503 ! be consistent with the other parts of the code. However in the modified !
1504 ! code, I allowed 'kinv' to be 1 & if 'kinv = 1', I just exit the program !
1505 ! without performing cumulus convection. This new approach seems to be !
1506 ! more reasonable: if PBL height is within 'kinv=1' layer, surface is STL !
1507 ! interface (bflxs <= 0) and interface just above the surface should be !
1508 ! either non-turbulent (Ri>0.19) or stably turbulent (0<=Ri<0.19 but this !
1509 ! interface is identified as a base external interface of upperlying CL. !
1510 ! Thus, when 'kinv=1', PBL scheme guarantees 'bflxs <= 0'. For this case !
1511 ! it is reasonable to assume that cumulus convection does not happen. !
1512 ! When these is SBCL, PBL height from the PBL scheme is likely to be very !
1513 ! close at 'kinv-1' interface, but not exactly, since 'zi' information is !
1514 ! changed between two model time steps. In order to ensure correct identi !
1515 ! fication of 'kinv' for general case including SBCL, I imposed an offset !
1516 ! of 5 [m] in the below 'kinv' finding block. !
1517 ! ----------------------------------------------------------------------- !
1519 do k = mkx - 1, 1, -1
1520 if( (pblh + 5._r8 - zs0(k))*(pblh + 5._r8 - zs0(k+1)) .lt. 0._r8 ) then
1528 if( kinv .le. 1 ) then
1529 exit_kinv1(i) = 1._r8
1533 ! From here, it must be 'kinv >= 2'.
1535 ! -------------------------------------------------------------------------- !
1536 ! Find PBL averaged tke ('tkeavg') and minimum 'thvl' ('thvlmin') in the PBL !
1537 ! In the current code, 'tkeavg' is obtained by averaging all interfacial TKE !
1538 ! within the PBL. However, in order to be conceptually consistent with PBL !
1539 ! scheme, 'tkeavg' should be calculated by considering surface buoyancy flux.!
1540 ! If surface buoyancy flux is positive ( bflxs >0 ), surface interfacial TKE !
1541 ! should be included in calculating 'tkeavg', while if bflxs <= 0, surface !
1542 ! interfacial TKE should not be included in calculating 'tkeavg'. I should !
1543 ! modify the code when 'bflxs' is available as an input of cumulus scheme. !
1544 ! 'thvlmin' is a minimum 'thvl' within PBL obtained by comparing top & base !
1545 ! interface values of 'thvl' in each layers within the PBL. !
1546 ! -------------------------------------------------------------------------- !
1551 do k = 0, kinv - 1 ! Here, 'k' is an interfacial layer index.
1553 dpi = ps0(0) - p0(1)
1554 elseif( k .eq. (kinv-1) ) then
1555 dpi = p0(kinv-1) - ps0(kinv-1)
1557 dpi = p0(k) - p0(k+1)
1560 tkeavg = tkeavg + dpi*tke(k)
1561 if( k .ne. 0 ) thvlmin = min(thvlmin,min(thvl0bot(k),thvl0top(k)))
1563 tkeavg = tkeavg/dpsum
1565 ! ------------------------------------------------------------------ !
1566 ! Find characteristics of cumulus source air: qtsrc,thlsrc,usrc,vsrc !
1567 ! Note that 'thlsrc' was con-cocked using 'thvlsrc' and 'qtsrc'. !
1568 ! 'qtsrc' is defined as the lowest layer mid-point value; 'thlsrc' !
1569 ! is from 'qtsrc' and 'thvlmin=thvlsrc'; 'usrc' & 'vsrc' are defined !
1570 ! as the values just below the PBL top interface. !
1571 ! ------------------------------------------------------------------ !
1575 thlsrc = thvlsrc / ( 1._r8 + zvir * qtsrc )
1576 usrc = u0(kinv-1) + ssu0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
1577 vsrc = v0(kinv-1) + ssv0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
1582 ! ------------------------------------------------------------------ !
1583 ! Find LCL of the source air and a layer index containing LCL (klcl) !
1584 ! When the LCL is exactly at the interface, 'klcl' is a layer index !
1585 ! having 'plcl' as the lower interface similar to the 'kinv' case. !
1586 ! In the previous code, I assumed that if LCL is located within the !
1587 ! lowest model layer ( 1 ) or the top model layer ( mkx ), then no !
1588 ! convective adjustment is performed and just exited. However, in !
1589 ! the revised code, I relaxed the first constraint and even though !
1590 ! LCL is at the lowest model layer, I allowed cumulus convection to !
1591 ! be initiated. For this case, cumulus convection should be started !
1592 ! from the PBL top height, as shown in the following code. !
1593 ! When source air is already saturated even at the surface, klcl is !
1595 ! ------------------------------------------------------------------ !
1597 plcl = qsinvert(qtsrc,thlsrc,ps0(0),qsat)
1599 if( ps0(k) .lt. plcl ) then
1608 if( plcl .lt. 30000._r8 ) then
1609 ! if( klcl .eq. mkx ) then
1610 exit_klclmkx(i) = 1._r8
1615 ! ------------------------------------------------------------- !
1616 ! Calculate environmental virtual potential temperature at LCL, !
1617 !'thv0lcl' which is solely used in the 'cin' calculation. Note !
1618 ! that 'thv0lcl' is calculated first by calculating 'thl0lcl' !
1619 ! and 'qt0lcl' at the LCL, and performing 'conden' afterward, !
1620 ! in fully consistent with the other parts of the code. !
1621 ! ------------------------------------------------------------- !
1623 thl0lcl = thl0(klcl) + ssthl0(klcl) * ( plcl - p0(klcl) )
1624 qt0lcl = qt0(klcl) + ssqt0(klcl) * ( plcl - p0(klcl) )
1625 call conden(plcl,thl0lcl,qt0lcl,thj,qvj,qlj,qij,qse,id_check,qsat)
1626 if( id_check .eq. 1 ) then
1627 exit_conden(i) = 1._r8
1631 thv0lcl = thj * ( 1._r8 + zvir * qvj - qlj - qij )
1633 ! ------------------------------------------------------------------------ !
1634 ! Compute Convective Inhibition, 'cin' & 'cinlcl' [J/kg]=[m2/s2] TKE unit. !
1636 ! 'cin' (cinlcl) is computed from the PBL top interface to LFC (LCL) using !
1637 ! piecewisely reconstructed environmental profiles, assuming environmental !
1638 ! buoyancy profile within each layer ( or from LCL to upper interface in !
1639 ! each layer ) is simply a linear profile. For the purpose of cin (cinlcl) !
1640 ! calculation, we simply assume that lateral entrainment does not occur in !
1641 ! updrafting cumulus plume, i.e., cumulus source air property is conserved.!
1642 ! Below explains some rules used in the calculations of cin (cinlcl). In !
1643 ! general, both 'cin' and 'cinlcl' are calculated from a PBL top interface !
1644 ! to LCL and LFC, respectively : !
1645 ! 1. If LCL is lower than the PBL height, cinlcl = 0 and cin is calculated !
1646 ! from PBL height to LFC. !
1647 ! 2. If LCL is higher than PBL height, 'cinlcl' is calculated by summing !
1648 ! both positive and negative cloud buoyancy up to LCL using 'single_cin'!
1649 ! From the LCL to LFC, however, only negative cloud buoyancy is counted !
1650 ! to calculate final 'cin' upto LFC. !
1651 ! 3. If either 'cin' or 'cinlcl' is negative, they are set to be zero. !
1652 ! In the below code, 'klfc' is the layer index containing 'LFC' similar to !
1653 ! 'kinv' and 'klcl'. !
1654 ! ------------------------------------------------------------------------ !
1661 ! ------------------------------------------------------------------------- !
1662 ! Case 1. LCL height is higher than PBL interface ( 'pLCL <= ps0(kinv-1)' ) !
1663 ! ------------------------------------------------------------------------- !
1665 if( klcl .ge. kinv ) then
1667 do k = kinv, mkx - 1
1668 if( k .lt. klcl ) then
1671 cin = cin + single_cin(ps0(k-1),thv0bot(k),ps0(k),thv0top(k),thvubot,thvutop)
1672 elseif( k .eq. klcl ) then
1673 !----- Bottom to LCL
1676 cin = cin + single_cin(ps0(k-1),thv0bot(k),plcl,thv0lcl,thvubot,thvutop)
1677 if( cin .lt. 0._r8 ) limit_cinlcl(i) = 1._r8
1678 cinlcl = max(cin,0._r8)
1682 call conden(ps0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check,qsat)
1683 if( id_check .eq. 1 ) then
1684 exit_conden(i) = 1._r8
1688 thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij )
1689 call getbuoy(plcl,thv0lcl,ps0(k),thv0top(k),thvubot,thvutop,plfc,cin)
1690 if( plfc .gt. 0._r8 ) then
1696 call conden(ps0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check,qsat)
1697 if( id_check .eq. 1 ) then
1698 exit_conden(i) = 1._r8
1702 thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij )
1703 call getbuoy(ps0(k-1),thv0bot(k),ps0(k),thv0top(k),thvubot,thvutop,plfc,cin)
1704 if( plfc .gt. 0._r8 ) then
1711 ! ----------------------------------------------------------------------- !
1712 ! Case 2. LCL height is lower than PBL interface ( 'pLCL > ps0(kinv-1)' ) !
1713 ! ----------------------------------------------------------------------- !
1717 do k = kinv, mkx - 1
1718 call conden(ps0(k-1),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check,qsat)
1719 if( id_check .eq. 1 ) then
1720 exit_conden(i) = 1._r8
1724 thvubot = thj * ( 1._r8 + zvir*qvj - qlj - qij )
1725 call conden(ps0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check,qsat)
1726 if( id_check .eq. 1 ) then
1727 exit_conden(i) = 1._r8
1731 thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij )
1732 call getbuoy(ps0(k-1),thv0bot(k),ps0(k),thv0top(k),thvubot,thvutop,plfc,cin)
1733 if( plfc .gt. 0._r8 ) then
1738 endif ! End of CIN case selection
1741 if( cin .lt. 0._r8 ) limit_cin(i) = 1._r8
1742 cin = max(0._r8,cin)
1743 if( klfc .ge. mkx ) then
1745 ! write(iulog,*) 'klfc >= mkx'
1746 exit_klfcmkx(i) = 1._r8
1751 ! ---------------------------------------------------------------------- !
1752 ! In order to calculate implicit 'cin' (or 'cinlcl'), save the initially !
1753 ! calculated 'cin' and 'cinlcl', and other related variables. These will !
1754 ! be restored after calculating implicit CIN. !
1755 ! ---------------------------------------------------------------------- !
1757 if( iter .eq. 1 ) then
1760 ke = rbuoy / ( rkfre * tkeavg + epsvarw )
1775 trsrc_o(m) = trsrc(m)
1779 ! Modification : If I impose w = max(0.1_r8, w) up to the top interface of
1780 ! klfc, I should only use cinlfc. That is, if I want to
1781 ! use cinlcl, I should not impose w = max(0.1_r8, w).
1782 ! Using cinlcl is equivalent to treating only 'saturated'
1783 ! moist convection. Note that in this sense, I should keep
1784 ! the functionality of both cinlfc and cinlcl.
1785 ! However, the treatment of penetrative entrainment level becomes
1786 ! ambiguous if I choose 'cinlcl'. Thus, the best option is to use
1789 ! -------------------------------------------------------------------------- !
1790 ! Calculate implicit 'cin' by averaging initial and final cins. Note that !
1791 ! implicit CIN is adopted only when cumulus convection stabilized the system,!
1792 ! i.e., only when 'del_CIN >0'. If 'del_CIN<=0', just use explicit CIN. Note !
1793 ! also that since 'cinlcl' is set to zero whenever LCL is below the PBL top, !
1794 ! (see above CIN calculation part), the use of 'implicit CIN=cinlcl' is not !
1795 ! good. Thus, when using implicit CIN, always try to only use 'implicit CIN= !
1796 ! cin', not 'implicit CIN=cinlcl'. However, both 'CIN=cin' and 'CIN=cinlcl' !
1797 ! are good when using explicit CIN. !
1798 ! -------------------------------------------------------------------------- !
1800 if( iter .ne. 1 ) then
1804 if( use_CINcin ) then
1805 del_CIN = cin_f - cin_i
1807 del_CIN = cinlcl_f - cinlcl_i
1810 if( del_CIN .gt. 0._r8 ) then
1812 ! -------------------------------------------------------------- !
1813 ! Calculate implicit 'cin' and 'cinlcl'. Note that when we chose !
1814 ! to use 'implicit CIN = cin', choose 'cinlcl = cinlcl_i' below: !
1815 ! because iterative CIN only aims to obtain implicit CIN, once !
1816 ! we obtained 'implicit CIN=cin', it is good to use the original !
1817 ! profiles information for all the other variables after that. !
1818 ! Note 'cinlcl' will be explicitly used in calculating 'wlcl' & !
1819 ! 'ufrclcl' after calculating 'winv' & 'ufrcinv' at the PBL top !
1820 ! interface later, after calculating 'cbmf'. !
1821 ! -------------------------------------------------------------- !
1823 alpha = compute_alpha( del_CIN, ke )
1824 cin = cin_i + alpha * del_CIN
1825 if( use_CINcin ) then
1828 cinlcl = cinlcl_i + alpha * del_cinlcl
1831 ! ----------------------------------------------------------------- !
1832 ! Restore the original values from the previous 'iter_cin' step (1) !
1833 ! to compute correct tendencies for (n+1) time step by implicit CIN !
1834 ! ----------------------------------------------------------------- !
1850 trsrc(m) = trsrc_o(m)
1853 qv0(:mkx) = qv0_o(:mkx)
1854 ql0(:mkx) = ql0_o(:mkx)
1855 qi0(:mkx) = qi0_o(:mkx)
1856 t0(:mkx) = t0_o(:mkx)
1857 s0(:mkx) = s0_o(:mkx)
1858 u0(:mkx) = u0_o(:mkx)
1859 v0(:mkx) = v0_o(:mkx)
1860 qt0(:mkx) = qt0_o(:mkx)
1861 thl0(:mkx) = thl0_o(:mkx)
1862 thvl0(:mkx) = thvl0_o(:mkx)
1863 ssthl0(:mkx) = ssthl0_o(:mkx)
1864 ssqt0(:mkx) = ssqt0_o(:mkx)
1865 thv0bot(:mkx) = thv0bot_o(:mkx)
1866 thv0top(:mkx) = thv0top_o(:mkx)
1867 thvl0bot(:mkx) = thvl0bot_o(:mkx)
1868 thvl0top(:mkx) = thvl0top_o(:mkx)
1869 ssu0(:mkx) = ssu0_o(:mkx)
1870 ssv0(:mkx) = ssv0_o(:mkx)
1872 tr0(:mkx,m) = tr0_o(:mkx,m)
1873 sstr0(:mkx,m) = sstr0_o(:mkx,m)
1876 ! ------------------------------------------------------ !
1877 ! Initialize all fluxes, tendencies, and other variables !
1878 ! in association with cumulus convection. !
1879 ! ------------------------------------------------------ !
1883 slflx(0:mkx) = 0.0_r8
1884 qtflx(0:mkx) = 0.0_r8
1885 uflx(0:mkx) = 0.0_r8
1886 vflx(0:mkx) = 0.0_r8
1887 qvten(:mkx) = 0.0_r8
1888 qlten(:mkx) = 0.0_r8
1889 qiten(:mkx) = 0.0_r8
1893 qrten(:mkx) = 0.0_r8
1894 qsten(:mkx) = 0.0_r8
1895 dwten(:mkx) = 0.0_r8
1896 diten(:mkx) = 0.0_r8
1899 evapc(:mkx) = 0.0_r8
1900 cufrc(:mkx) = 0.0_r8
1913 qtten(:mkx) = 0.0_r8
1914 slten(:mkx) = 0.0_r8
1915 ufrc(0:mkx) = 0.0_r8
1917 thlu(0:mkx) = 0.0_r8
1922 thvu(0:mkx) = 0.0_r8
1923 thlu_emf(0:mkx) = 0.0_r8
1924 qtu_emf(0:mkx) = 0.0_r8
1925 uu_emf(0:mkx) = 0.0_r8
1926 vu_emf(0:mkx) = 0.0_r8
1929 trflx(0:mkx,m) = 0.0_r8
1930 trten(:mkx,m) = 0.0_r8
1931 tru(0:mkx,m) = 0.0_r8
1932 tru_emf(0:mkx,m) = 0.0_r8
1935 ! -------------------------------------------------- !
1936 ! Below are diagnostic output variables for detailed !
1937 ! analysis of cumulus scheme. !
1938 ! -------------------------------------------------- !
1940 ufrcinvbase = 0.0_r8
1946 excessu_arr(:mkx) = 0.0_r8
1947 excess0_arr(:mkx) = 0.0_r8
1948 xc_arr(:mkx) = 0.0_r8
1949 aquad_arr(:mkx) = 0.0_r8
1950 bquad_arr(:mkx) = 0.0_r8
1951 cquad_arr(:mkx) = 0.0_r8
1952 bogbot_arr(:mkx) = 0.0_r8
1953 bogtop_arr(:mkx) = 0.0_r8
1955 else ! When 'del_CIN < 0', use explicit CIN instead of implicit CIN.
1957 ! ----------------------------------------------------------- !
1958 ! Identifier showing whether explicit or implicit CIN is used !
1959 ! ----------------------------------------------------------- !
1961 ind_delcin(i) = 1._r8
1963 ! --------------------------------------------------------- !
1964 ! Restore original output values of "iter_cin = 1" and exit !
1965 ! --------------------------------------------------------- !
1967 umf_out(i,0:mkx) = umf_s(0:mkx)
1968 qvten_out(i,:mkx) = qvten_s(:mkx)
1969 qlten_out(i,:mkx) = qlten_s(:mkx)
1970 qiten_out(i,:mkx) = qiten_s(:mkx)
1971 sten_out(i,:mkx) = sten_s(:mkx)
1972 uten_out(i,:mkx) = uten_s(:mkx)
1973 vten_out(i,:mkx) = vten_s(:mkx)
1974 qrten_out(i,:mkx) = qrten_s(:mkx)
1975 qsten_out(i,:mkx) = qsten_s(:mkx)
1976 precip_out(i) = precip_s
1977 snow_out(i) = snow_s
1978 evapc_out(i,:mkx) = evapc_s(:mkx)
1979 cush_inout(i) = cush_s
1980 cufrc_out(i,:mkx) = cufrc_s(:mkx)
1981 slflx_out(i,0:mkx) = slflx_s(0:mkx)
1982 qtflx_out(i,0:mkx) = qtflx_s(0:mkx)
1983 qcu_out(i,:mkx) = qcu_s(:mkx)
1984 qlu_out(i,:mkx) = qlu_s(:mkx)
1985 qiu_out(i,:mkx) = qiu_s(:mkx)
1986 cbmf_out(i) = cbmf_s
1987 qc_out(i,:mkx) = qc_s(:mkx)
1988 rliq_out(i) = rliq_s
1992 trten_out(i,:mkx,m) = trten_s(:mkx,m)
1995 ! ------------------------------------------------------------------------------ !
1996 ! Below are diagnostic output variables for detailed analysis of cumulus scheme. !
1997 ! The order of vertical index is reversed for this internal diagnostic output. !
1998 ! ------------------------------------------------------------------------------ !
2000 fer_out(i,mkx:1:-1) = fer_s(:mkx)
2001 fdr_out(i,mkx:1:-1) = fdr_s(:mkx)
2003 cinlclh_out(i) = cinlcl_s
2004 qtten_out(i,mkx:1:-1) = qtten_s(:mkx)
2005 slten_out(i,mkx:1:-1) = slten_s(:mkx)
2006 ufrc_out(i,mkx:0:-1) = ufrc_s(0:mkx)
2007 uflx_out(i,mkx:0:-1) = uflx_s(0:mkx)
2008 vflx_out(i,mkx:0:-1) = vflx_s(0:mkx)
2010 ufrcinvbase_out(i) = ufrcinvbase_s
2011 ufrclcl_out(i) = ufrclcl_s
2012 winvbase_out(i) = winvbase_s
2013 wlcl_out(i) = wlcl_s
2014 plcl_out(i) = plcl_s
2015 pinv_out(i) = pinv_s
2016 plfc_out(i) = plfc_s
2017 pbup_out(i) = pbup_s
2018 ppen_out(i) = ppen_s
2019 qtsrc_out(i) = qtsrc_s
2020 thlsrc_out(i) = thlsrc_s
2021 thvlsrc_out(i) = thvlsrc_s
2022 emfkbup_out(i) = emfkbup_s
2023 cbmflimit_out(i) = cbmflimit_s
2024 tkeavg_out(i) = tkeavg_s
2025 zinv_out(i) = zinv_s
2026 rcwp_out(i) = rcwp_s
2027 rlwp_out(i) = rlwp_s
2028 riwp_out(i) = riwp_s
2030 wu_out(i,mkx:0:-1) = wu_s(0:mkx)
2031 qtu_out(i,mkx:0:-1) = qtu_s(0:mkx)
2032 thlu_out(i,mkx:0:-1) = thlu_s(0:mkx)
2033 thvu_out(i,mkx:0:-1) = thvu_s(0:mkx)
2034 uu_out(i,mkx:0:-1) = uu_s(0:mkx)
2035 vu_out(i,mkx:0:-1) = vu_s(0:mkx)
2036 qtu_emf_out(i,mkx:0:-1) = qtu_emf_s(0:mkx)
2037 thlu_emf_out(i,mkx:0:-1) = thlu_emf_s(0:mkx)
2038 uu_emf_out(i,mkx:0:-1) = uu_emf_s(0:mkx)
2039 vu_emf_out(i,mkx:0:-1) = vu_emf_s(0:mkx)
2040 uemf_out(i,mkx:0:-1) = uemf_s(0:mkx)
2042 dwten_out(i,mkx:1:-1) = dwten_s(:mkx)
2043 diten_out(i,mkx:1:-1) = diten_s(:mkx)
2044 flxrain_out(i,mkx:0:-1) = flxrain_s(0:mkx)
2045 flxsnow_out(i,mkx:0:-1) = flxsnow_s(0:mkx)
2046 ntraprd_out(i,mkx:1:-1) = ntraprd_s(:mkx)
2047 ntsnprd_out(i,mkx:1:-1) = ntsnprd_s(:mkx)
2049 excessu_arr_out(i,mkx:1:-1) = excessu_arr_s(:mkx)
2050 excess0_arr_out(i,mkx:1:-1) = excess0_arr_s(:mkx)
2051 xc_arr_out(i,mkx:1:-1) = xc_arr_s(:mkx)
2052 aquad_arr_out(i,mkx:1:-1) = aquad_arr_s(:mkx)
2053 bquad_arr_out(i,mkx:1:-1) = bquad_arr_s(:mkx)
2054 cquad_arr_out(i,mkx:1:-1) = cquad_arr_s(:mkx)
2055 bogbot_arr_out(i,mkx:1:-1) = bogbot_arr_s(:mkx)
2056 bogtop_arr_out(i,mkx:1:-1) = bogtop_arr_s(:mkx)
2059 trflx_out(i,mkx:0:-1,m) = trflx_s(0:mkx,m)
2060 tru_out(i,mkx:0:-1,m) = tru_s(0:mkx,m)
2061 tru_emf_out(i,mkx:0:-1,m) = tru_emf_s(0:mkx,m)
2071 ! ------------------------------------------------------------------ !
2072 ! Define a release level, 'prel' and release layer, 'krel'. !
2073 ! 'prel' is the lowest level from which buoyancy sorting occurs, and !
2074 ! 'krel' is the layer index containing 'prel' in it, similar to the !
2075 ! previous definitions of 'kinv', 'klcl', and 'klfc'. In order to !
2076 ! ensure that only PBL scheme works within the PBL, if LCL is below !
2077 ! PBL top height, then 'krel = kinv', while if LCL is above PBL top !
2078 ! height, then 'krel = klcl'. Note however that regardless of the !
2079 ! definition of 'krel', cumulus convection induces fluxes within PBL !
2080 ! through 'fluxbelowinv'. We can make cumulus convection start from !
2081 ! any level, even within the PBL by appropriately defining 'krel' & !
2082 ! 'prel' here. Then it must be accompanied by appropriate definition !
2083 ! of source air properties, CIN, and re-setting of 'fluxbelowinv', & !
2084 ! many other stuffs. !
2085 ! Note that even when 'prel' is located above the PBL top height, we !
2086 ! still have cumulus convection between PBL top height and 'prel': !
2087 ! we simply assume that no lateral mixing occurs in this range. !
2088 ! ------------------------------------------------------------------ !
2090 if( klcl .lt. kinv ) then
2093 thv0rel = thv0bot(krel)
2100 ! --------------------------------------------------------------------------- !
2101 ! Calculate cumulus base mass flux ('cbmf'), fractional area ('ufrcinv'), and !
2102 ! and mean vertical velocity (winv) of cumulus updraft at PBL top interface. !
2103 ! Also, calculate updraft fractional area (ufrclcl) and vertical velocity at !
2104 ! the LCL (wlcl). When LCL is below PBLH, cinlcl = 0 and 'ufrclcl = ufrcinv', !
2105 ! and 'wlcl = winv. !
2106 ! Only updrafts strong enough to overcome CIN can rise over PBL top interface.!
2107 ! Thus, in order to calculate cumulus mass flux at PBL top interface, 'cbmf',!
2108 ! we need to know 'CIN' ( the strength of potential energy barrier ) and !
2109 ! 'sigmaw' ( a standard deviation of updraft vertical velocity at the PBL top !
2110 ! interface, a measure of turbulentce strength in the PBL ). Naturally, the !
2111 ! ratio of these two variables, 'mu' - normalized CIN by TKE- is key variable !
2112 ! controlling 'cbmf'. If 'mu' becomes large, only small fraction of updrafts !
2113 ! with very strong TKE can rise over the PBL - both 'cbmf' and 'ufrc' becomes !
2114 ! small, but 'winv' becomes large ( this can be easily understood by PDF of w !
2115 ! at PBL top ). If 'mu' becomes small, lots of updraft can rise over the PBL !
2116 ! top - both 'cbmf' and 'ufrc' becomes large, but 'winv' becomes small. Thus, !
2117 ! all of the key variables associated with cumulus convection at the PBL top !
2118 ! - 'cbmf', 'ufrc', 'winv' where 'cbmf = rho*ufrc*winv' - are a unique functi !
2119 ! ons of 'mu', normalized CIN. Although these are uniquely determined by 'mu',!
2120 ! we usually impose two comstraints on 'cbmf' and 'ufrc': (1) because we will !
2121 ! simply assume that subsidence warming and drying of 'kinv-1' layer in assoc !
2122 ! iation with 'cbmf' at PBL top interface is confined only in 'kinv-1' layer, !
2123 ! cbmf must not be larger than the mass within the 'kinv-1' layer. Otherwise, !
2124 ! instability will occur due to the breaking of stability con. If we consider !
2125 ! semi-Lagrangian vertical advection scheme and explicitly consider the exten !
2126 ! t of vertical movement of each layer in association with cumulus mass flux, !
2127 ! we don't need to impose this constraint. However, using a semi-Lagrangian !
2128 ! scheme is a future research subject. Note that this constraint should be ap !
2129 ! plied for all interfaces above PBL top as well as PBL top interface. As a !
2130 ! result, this 'cbmf' constraint impose a 'lower' limit on mu - 'mumin0'. (2) !
2131 ! in order for mass flux parameterization - rho*(w'a')= M*(a_c-a_e) - to be !
2132 ! valid, cumulus updraft fractional area should be much smaller than 1. In !
2133 ! current code, we impose 'rmaxfrac = 0.1 ~ 0.2' through the whole vertical !
2134 ! layers where cumulus convection occurs. At the PBL top interface, the same !
2135 ! constraint is made by imposing another lower 'lower' limit on mu, 'mumin1'. !
2136 ! After that, also limit 'ufrclcl' to be smaller than 'rmaxfrac' by 'mumin2'. !
2137 ! --------------------------------------------------------------------------- !
2139 ! --------------------------------------------------------------------------- !
2140 ! Calculate normalized CIN, 'mu' satisfying all the three constraints imposed !
2141 ! on 'cbmf'('mumin0'), 'ufrc' at the PBL top - 'ufrcinv' - ( by 'mumin1' from !
2142 ! a parameter sentence), and 'ufrc' at the LCL - 'ufrclcl' ( by 'mumin2'). !
2143 ! Note that 'cbmf' does not change between PBL top and LCL because we assume !
2144 ! that buoyancy sorting does not occur when cumulus updraft is unsaturated. !
2145 ! --------------------------------------------------------------------------- !
2147 if( use_CINcin ) then
2148 wcrit = sqrt( 2._r8 * cin * rbuoy )
2150 wcrit = sqrt( 2._r8 * cinlcl * rbuoy )
2152 sigmaw = sqrt( rkfre * tkeavg + epsvarw )
2153 mu = wcrit/sigmaw/1.4142_r8
2154 if( mu .ge. 3._r8 ) then
2155 ! write(iulog,*) 'mu >= 3'
2159 rho0inv = ps0(kinv-1)/(r*thv0top(kinv-1)*exns0(kinv-1))
2160 cbmf = (rho0inv*sigmaw/2.5066_r8)*exp(-mu**2)
2161 ! 1. 'cbmf' constraint
2162 cbmflimit = 0.9_r8*dp0(kinv-1)/g/dt
2164 if( cbmf .gt. cbmflimit ) mumin0 = sqrt(-log(2.5066_r8*cbmflimit/rho0inv/sigmaw))
2165 ! 2. 'ufrcinv' constraint
2166 mu = max(max(mu,mumin0),mumin1)
2167 ! 3. 'ufrclcl' constraint
2168 mulcl = sqrt(2._r8*cinlcl*rbuoy)/1.4142_r8/sigmaw
2169 mulclstar = sqrt(max(0._r8,2._r8*(exp(-mu**2)/2.5066_r8)**2*(1._r8/erfc(mu)**2-0.25_r8/rmaxfrac**2)))
2170 if( mulcl .gt. 1.e-8_r8 .and. mulcl .gt. mulclstar ) then
2171 mumin2 = compute_mumin2(mulcl,rmaxfrac,mu)
2172 if( mu .gt. mumin2 ) then
2173 write(iulog,*) 'Critical error in mu calculation in UW_ShCu'
2175 call wrf_message(iulog)
2180 if( mu .eq. mumin2 ) limit_ufrc(i) = 1._r8
2182 if( mu .eq. mumin0 ) limit_cbmf(i) = 1._r8
2183 if( mu .eq. mumin1 ) limit_ufrc(i) = 1._r8
2185 ! ------------------------------------------------------------------- !
2186 ! Calculate final ['cbmf','ufrcinv','winv'] at the PBL top interface. !
2187 ! Note that final 'cbmf' here is obtained in such that 'ufrcinv' and !
2188 ! 'ufrclcl' are smaller than ufrcmax with no instability. !
2189 ! ------------------------------------------------------------------- !
2191 cbmf = (rho0inv*sigmaw/2.5066_r8)*exp(-mu**2)
2192 winv = sigmaw*(2._r8/2.5066_r8)*exp(-mu**2)/erfc(mu)
2193 ufrcinv = cbmf/winv/rho0inv
2195 ! ------------------------------------------------------------------- !
2196 ! Calculate ['ufrclcl','wlcl'] at the LCL. When LCL is below PBL top, !
2197 ! it automatically becomes 'ufrclcl = ufrcinv' & 'wlcl = winv', since !
2198 ! it was already set to 'cinlcl=0' if LCL is below PBL top interface. !
2199 ! Note 'cbmf' at the PBL top is the same as 'cbmf' at the LCL. Note !
2200 ! also that final 'cbmf' here is obtained in such that 'ufrcinv' and !
2201 ! 'ufrclcl' are smaller than ufrcmax and there is no instability. !
2202 ! By construction, it must be 'wlcl > 0' but for assurance, I checked !
2203 ! this again in the below block. If 'ufrclcl < 0.1%', just exit. !
2204 ! ------------------------------------------------------------------- !
2206 wtw = winv * winv - 2._r8 * cinlcl * rbuoy
2207 if( wtw .le. 0._r8 ) then
2208 ! write(iulog,*) 'wlcl < 0 at the LCL'
2214 ufrclcl = cbmf/wlcl/rho0inv
2216 if( ufrclcl .le. 0.0001_r8 ) then
2217 ! write(iulog,*) 'ufrclcl <= 0.0001'
2218 exit_ufrc(i) = 1._r8
2222 ufrc(krel-1) = ufrclcl
2224 ! ----------------------------------------------------------------------- !
2225 ! Below is just diagnostic output for detailed analysis of cumulus scheme !
2226 ! ----------------------------------------------------------------------- !
2228 ufrcinvbase = ufrcinv
2230 umf(kinv-1:krel-1) = cbmf
2231 wu(kinv-1:krel-1) = winv
2233 ! -------------------------------------------------------------------------- !
2234 ! Define updraft properties at the level where buoyancy sorting starts to be !
2235 ! happening, i.e., by definition, at 'prel' level within the release layer. !
2236 ! Because no lateral entrainment occurs upto 'prel', conservative scalars of !
2237 ! cumulus updraft at release level is same as those of source air. However, !
2238 ! horizontal momentums of source air are modified by horizontal PGF forcings !
2239 ! from PBL top interface to 'prel'. For this case, we should add additional !
2240 ! horizontal momentum from PBL top interface to 'prel' as will be done below !
2241 ! to 'usrc' and 'vsrc'. Note that below cumulus updraft properties - umf, wu,!
2242 ! thlu, qtu, thvu, uu, vu - are defined all interfaces not at the layer mid- !
2243 ! point. From the index notation of cumulus scheme, wu(k) is the cumulus up- !
2244 ! draft vertical velocity at the top interface of k layer. !
2245 ! Diabatic horizontal momentum forcing should be treated as a kind of 'body' !
2246 ! forcing without actual mass exchange between convective updraft and !
2247 ! environment, but still taking horizontal momentum from the environment to !
2248 ! the convective updrafts. Thus, diabatic convective momentum transport !
2249 ! vertically redistributes environmental horizontal momentum. !
2250 ! -------------------------------------------------------------------------- !
2255 thlu(krel-1) = thlsrc
2257 call conden(prel,thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check,qsat)
2258 if( id_check .eq. 1 ) then
2259 exit_conden(i) = 1._r8
2263 thvu(krel-1) = thj * ( 1._r8 + zvir*qvj - qlj - qij )
2267 if( krel .eq. kinv ) then
2268 uplus = PGFc * ssu0(kinv) * ( prel - ps0(kinv-1) )
2269 vplus = PGFc * ssv0(kinv) * ( prel - ps0(kinv-1) )
2271 do k = kinv, max(krel-1,kinv)
2272 uplus = uplus + PGFc * ssu0(k) * ( ps0(k) - ps0(k-1) )
2273 vplus = vplus + PGFc * ssv0(k) * ( ps0(k) - ps0(k-1) )
2275 uplus = uplus + PGFc * ssu0(krel) * ( prel - ps0(krel-1) )
2276 vplus = vplus + PGFc * ssv0(krel) * ( prel - ps0(krel-1) )
2278 uu(krel-1) = usrc + uplus
2279 vu(krel-1) = vsrc + vplus
2282 tru(krel-1,m) = trsrc(m)
2285 ! -------------------------------------------------------------------------- !
2286 ! Define environmental properties at the level where buoyancy sorting occurs !
2287 ! ('pe', normally, layer midpoint except in the 'krel' layer). In the 'krel' !
2288 ! layer where buoyancy sorting starts to occur, however, 'pe' is defined !
2289 ! differently because LCL is regarded as lower interface for mixing purpose. !
2290 ! -------------------------------------------------------------------------- !
2292 pe = 0.5_r8 * ( prel + ps0(krel) )
2293 dpe = prel - ps0(krel)
2296 thle = thl0(krel) + ssthl0(krel) * ( pe - p0(krel) )
2297 qte = qt0(krel) + ssqt0(krel) * ( pe - p0(krel) )
2298 ue = u0(krel) + ssu0(krel) * ( pe - p0(krel) )
2299 ve = v0(krel) + ssv0(krel) * ( pe - p0(krel) )
2301 tre(m) = tr0(krel,m) + sstr0(krel,m) * ( pe - p0(krel) )
2304 !-------------------------!
2305 ! Buoyancy-Sorting Mixing !
2306 !-------------------------!------------------------------------------------ !
2308 ! In order to complete buoyancy-sorting mixing at layer mid-point, and so !
2309 ! calculate 'updraft mass flux, updraft w velocity, conservative scalars' !
2310 ! at the upper interface of each layer, we need following 3 information. !
2312 ! 1. Pressure where mixing occurs ('pe'), and temperature at 'pe' which is !
2313 ! necessary to calculate various thermodynamic coefficients at pe. This !
2314 ! temperature is obtained by undiluted cumulus properties lifted to pe. !
2315 ! 2. Undiluted updraft properties at pe - conservative scalar and vertical !
2316 ! velocity -which are assumed to be the same as the properties at lower !
2317 ! interface only for calculation of fractional lateral entrainment and !
2318 ! detrainment rate ( fer(k) and fdr(k) [Pa-1] ), respectively. Final !
2319 ! values of cumulus conservative scalars and w at the top interface are !
2320 ! calculated afterward after obtaining fer(k) & fdr(k). !
2321 ! 3. Environmental properties at pe. !
2322 ! ------------------------------------------------------------------------- !
2324 ! ------------------------------------------------------------------------ !
2325 ! Define cumulus scale height. !
2326 ! Cumulus scale height is defined as the maximum height cumulus can reach. !
2327 ! In case of premitive code, cumulus scale height ('cush') at the current !
2328 ! time step was assumed to be the same as 'cush' of previous time step. !
2329 ! However, I directly calculated cush at each time step using an iterative !
2330 ! method. Note that within the cumulus scheme, 'cush' information is used !
2331 ! only at two places during buoyancy-sorting process: !
2332 ! (1) Even negatively buoyancy mixtures with strong vertical velocity !
2333 ! enough to rise up to 'rle*scaleh' (rle = 0.1) from pe are entrained !
2334 ! into cumulus updraft, !
2335 ! (2) The amount of mass that is involved in buoyancy-sorting mixing !
2336 ! process at pe is rei(k) = rkm/scaleh/rho*g [Pa-1] !
2337 ! In terms of (1), I think critical stopping distance might be replaced by !
2338 ! layer thickness. In future, we will use rei(k) = (0.5*rkm/z0(k)/rho/g). !
2339 ! In the premitive code, 'scaleh' was largely responsible for the jumping !
2340 ! variation of precipitation amount. !
2341 ! ------------------------------------------------------------------------ !
2344 if( tscaleh .lt. 0.0_r8 ) scaleh = 1000._r8
2346 ! Save time : Set iter_scaleh = 1. This will automatically use 'cush' from the previous time step
2347 ! at the first implicit iteration. At the second implicit iteration, it will use
2348 ! the updated 'cush' by the first implicit cin. So, this updating has an effect of
2349 ! doing one iteration for cush calculation, which is good.
2350 ! So, only this setting of 'iter_scaleh = 1' is sufficient-enough to save computation time.
2353 do iter_scaleh = 1, 3
2355 ! ---------------------------------------------------------------- !
2356 ! Initialization of 'kbup' and 'kpen' !
2357 ! ---------------------------------------------------------------- !
2358 ! 'kbup' is the top-most layer in which cloud buoyancy is positive !
2359 ! both at the top and bottom interface of the layer. 'kpen' is the !
2360 ! layer upto which cumulus panetrates ,i.e., cumulus w at the base !
2361 ! interface is positive, but becomes negative at the top interface.!
2362 ! Here, we initialize 'kbup' and 'kpen'. These initializations are !
2363 ! not trivial but important, expecially in calculating turbulent !
2364 ! fluxes without confliction among several physics as explained in !
2365 ! detail in the part of turbulent fluxes calculation later. Note !
2366 ! that regardless of whether 'kbup' and 'kpen' are updated or not !
2367 ! during updraft motion, penetrative entrainments are dumped down !
2368 ! across the top interface of 'kbup' later. More specifically,!
2369 ! penetrative entrainment heat and moisture fluxes are calculated !
2370 ! from the top interface of 'kbup' layer to the base interface of !
2371 ! 'kpen' layer. Because of this, initialization of 'kbup' & 'kpen' !
2372 ! influence the convection system when there are not updated. The !
2373 ! below initialization of 'kbup = krel' assures that penetrative !
2374 ! entrainment fluxes always occur at interfaces above the PBL top !
2375 ! interfaces (i.e., only at interfaces k >=kinv ), which seems to !
2376 ! be attractable considering that the most correct fluxes at the !
2377 ! PBL top interface can be ontained from the 'fluxbelowinv' using !
2378 ! reconstructed PBL height. !
2379 ! The 'kbup = krel'(after going through the whole buoyancy sorting !
2380 ! proces during updraft motion) implies that cumulus updraft from !
2381 ! the PBL top interface can not reach to the LFC,so that 'kbup' is !
2382 ! not updated during upward. This means that cumulus updraft did !
2383 ! not fully overcome the buoyancy barrier above just the PBL top. !
2384 ! If 'kpen' is not updated either ( i.e., cumulus cannot rise over !
2385 ! the top interface of release layer),penetrative entrainment will !
2386 ! not happen at any interfaces. If cumulus updraft can rise above !
2387 ! the release layer but cannot fully overcome the buoyancy barrier !
2388 ! just above PBL top interface, penetratve entrainment occurs at !
2389 ! several above interfaces, including the top interface of release !
2390 ! layer. In the latter case, warming and drying tendencies will be !
2391 ! be initiated in 'krel' layer. Note current choice of 'kbup=krel' !
2392 ! is completely compatible with other flux physics without double !
2393 ! or miss counting turbulent fluxes at any interface. However, the !
2394 ! alternative choice of 'kbup=krel-1' also has itw own advantage - !
2395 ! when cumulus updraft cannot overcome buoyancy barrier just above !
2396 ! PBL top, entrainment warming and drying are concentrated in the !
2397 ! 'kinv-1' layer instead of 'kinv' layer for this case. This might !
2398 ! seems to be more dynamically reasonable, but I will choose the !
2399 ! 'kbup = krel' choice since it is more compatible with the other !
2400 ! parts of the code, expecially, when we chose ' use_emf=.false. ' !
2401 ! as explained in detail in turbulent flux calculation part. !
2402 ! ---------------------------------------------------------------- !
2407 ! ------------------------------------------------------------ !
2408 ! Since 'wtw' is continuously updated during vertical motion, !
2409 ! I need below initialization command within this 'iter_scaleh'!
2410 ! do loop. Similarily, I need initializations of environmental !
2411 ! properties at 'krel' layer as below. !
2412 ! ------------------------------------------------------------ !
2415 pe = 0.5_r8 * ( prel + ps0(krel) )
2416 dpe = prel - ps0(krel)
2419 thle = thl0(krel) + ssthl0(krel) * ( pe - p0(krel) )
2420 qte = qt0(krel) + ssqt0(krel) * ( pe - p0(krel) )
2421 ue = u0(krel) + ssu0(krel) * ( pe - p0(krel) )
2422 ve = v0(krel) + ssv0(krel) * ( pe - p0(krel) )
2424 tre(m) = tr0(krel,m) + sstr0(krel,m) * ( pe - p0(krel) )
2427 ! ----------------------------------------------------------------------- !
2428 ! Cumulus rises upward from 'prel' ( or base interface of 'krel' layer ) !
2429 ! until updraft vertical velocity becomes zero. !
2430 ! Buoyancy sorting is performed via two stages. (1) Using cumulus updraft !
2431 ! properties at the base interface of each layer,perform buoyancy sorting !
2432 ! at the layer mid-point, 'pe', and update cumulus properties at the top !
2433 ! interface, and then (2) by averaging updated cumulus properties at the !
2434 ! top interface and cumulus properties at the base interface, calculate !
2435 ! cumulus updraft properties at pe that will be used in buoyancy sorting !
2436 ! mixing - thlue, qtue and, wue. Using this averaged properties, perform !
2437 ! buoyancy sorting again at pe, and re-calculate fer(k) and fdr(k). Using !
2438 ! this recalculated fer(k) and fdr(k), finally calculate cumulus updraft !
2439 ! properties at the top interface - thlu, qtu, thvu, uu, vu. In the below,!
2440 ! 'iter_xc = 1' performs the first stage, while 'iter_xc= 2' performs the !
2441 ! second stage. We can increase the number of iterations, 'nter_xc'.as we !
2442 ! want, but a sample test indicated that about 3 - 5 iterations produced !
2443 ! satisfactory converent solution. Finally, identify 'kbup' and 'kpen'. !
2444 ! ----------------------------------------------------------------------- !
2446 do k = krel, mkx - 1 ! Here, 'k' is a layer index.
2455 do iter_xc = 1, niter_xc
2457 wtw = wu(km1) * wu(km1)
2459 ! ---------------------------------------------------------------- !
2460 ! Calculate environmental and cumulus saturation 'excess' at 'pe'. !
2461 ! Note that in order to calculate saturation excess, we should use !
2462 ! liquid water temperature instead of temperature as the argument !
2463 ! of "qsat". But note normal argument of "qsat" is temperature. !
2464 ! ---------------------------------------------------------------- !
2466 call conden(pe,thle,qte,thj,qvj,qlj,qij,qse,id_check,qsat)
2467 if( id_check .eq. 1 ) then
2468 exit_conden(i) = 1._r8
2472 thv0j = thj * ( 1._r8 + zvir*qvj - qlj - qij )
2473 rho0j = pe / ( r * thv0j * exne )
2474 qsat_arg = thle*exne
2475 status = qsat(qsat_arg,pe,es(1),qs(1),gam(1),1)
2476 excess0 = qte - qs(1)
2478 call conden(pe,thlue,qtue,thj,qvj,qlj,qij,qse,id_check,qsat)
2479 if( id_check .eq. 1 ) then
2480 exit_conden(i) = 1._r8
2484 ! ----------------------------------------------------------------- !
2485 ! Detrain excessive condensate larger than 'criqc' from the cumulus !
2486 ! updraft before performing buoyancy sorting. All I should to do is !
2487 ! to update 'thlue' & 'que' here. Below modification is completely !
2488 ! compatible with the other part of the code since 'thule' & 'qtue' !
2489 ! are used only for buoyancy sorting. I found that as long as I use !
2490 ! 'niter_xc >= 2', detraining excessive condensate before buoyancy !
2491 ! sorting has negligible influence on the buoyancy sorting results. !
2492 ! ----------------------------------------------------------------- !
2493 if( (qlj + qij) .gt. criqc ) then
2494 exql = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij )
2495 exqi = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij )
2496 qtue = qtue - exql - exqi
2497 thlue = thlue + (xlv/cp/exne)*exql + (xls/cp/exne)*exqi
2499 call conden(pe,thlue,qtue,thj,qvj,qlj,qij,qse,id_check,qsat)
2500 if( id_check .eq. 1 ) then
2501 exit_conden(i) = 1._r8
2505 thvj = thj * ( 1._r8 + zvir * qvj - qlj - qij )
2506 tj = thj * exne ! This 'tj' is used for computing thermo. coeffs. below
2507 qsat_arg = thlue*exne
2508 status = qsat(qsat_arg,pe,es(1),qs(1),gam(1),1)
2509 excessu = qtue - qs(1)
2511 ! ------------------------------------------------------------------- !
2512 ! Calculate critical mixing fraction, 'xc'. Mixture with mixing ratio !
2513 ! smaller than 'xc' will be entrained into cumulus updraft. Both the !
2514 ! saturated updrafts with 'positive buoyancy' or 'negative buoyancy + !
2515 ! strong vertical velocity enough to rise certain threshold distance' !
2516 ! are kept into the updraft in the below program. If the core updraft !
2517 ! is unsaturated, we can set 'xc = 0' and let the cumulus convection !
2518 ! still works or we may exit. !
2519 ! Current below code does not entrain unsaturated mixture. However it !
2520 ! should be modified such that it also entrain unsaturated mixture. !
2521 ! ------------------------------------------------------------------- !
2523 ! ----------------------------------------------------------------- !
2524 ! cridis : Critical stopping distance for buoyancy sorting purpose. !
2525 ! scaleh is only used here. !
2526 ! ----------------------------------------------------------------- !
2528 cridis = rle*scaleh ! Original code
2529 ! cridis = 1._r8*(zs0(k) - zs0(k-1)) ! New code
2531 ! ---------------- !
2532 ! Buoyancy Sorting !
2533 ! ---------------- !
2535 ! ----------------------------------------------------------------- !
2536 ! Case 1 : When both cumulus and env. are unsaturated or saturated. !
2537 ! ----------------------------------------------------------------- !
2539 if( ( excessu .le. 0._r8 .and. excess0 .le. 0._r8 ) .or. ( excessu .ge. 0._r8 .and. excess0 .ge. 0._r8 ) ) then
2540 xc = min(1._r8,max(0._r8,1._r8-2._r8*rbuoy*g*cridis/wue**2._r8*(1._r8-thvj/thv0j)))
2541 ! Below 3 lines are diagnostic output not influencing
2542 ! numerical calculations.
2547 ! -------------------------------------------------- !
2548 ! Case 2 : When either cumulus or env. is saturated. !
2549 ! -------------------------------------------------- !
2550 xsat = excessu / ( excessu - excess0 );
2551 thlxsat = thlue + xsat * ( thle - thlue );
2552 qtxsat = qtue + xsat * ( qte - qtue );
2553 call conden(pe,thlxsat,qtxsat,thj,qvj,qlj,qij,qse,id_check,qsat)
2554 if( id_check .eq. 1 ) then
2555 exit_conden(i) = 1._r8
2559 thvxsat = thj * ( 1._r8 + zvir * qvj - qlj - qij )
2560 ! -------------------------------------------------- !
2561 ! kk=1 : Cumulus Segment, kk=2 : Environment Segment !
2562 ! -------------------------------------------------- !
2564 if( kk .eq. 1 ) then
2566 thv_x1 = ( 1._r8 - 1._r8/xsat ) * thvj + ( 1._r8/xsat ) * thvxsat;
2569 thv_x0 = ( xsat / ( xsat - 1._r8 ) ) * thv0j + ( 1._r8/( 1._r8 - xsat ) ) * thvxsat;
2572 bquad = 2._r8*rbuoy*g*cridis*(thv_x1 - thv_x0)/thv0j - 2._r8*wue**2;
2573 cquad = 2._r8*rbuoy*g*cridis*(thv_x0 - thv0j)/thv0j + wue**2;
2574 if( kk .eq. 1 ) then
2575 if( ( bquad**2-4._r8*aquad*cquad ) .ge. 0._r8 ) then
2576 call roots(aquad,bquad,cquad,xs1,xs2,status)
2577 x_cu = min(1._r8,max(0._r8,min(xsat,min(xs1,xs2))))
2582 if( ( bquad**2-4._r8*aquad*cquad) .ge. 0._r8 ) then
2583 call roots(aquad,bquad,cquad,xs1,xs2,status)
2584 x_en = min(1._r8,max(0._r8,max(xsat,min(xs1,xs2))))
2590 if( x_cu .eq. xsat ) then
2591 xc = max(x_cu, x_en);
2597 ! ------------------------------------------------------------------------ !
2598 ! Compute fractional lateral entrainment & detrainment rate in each layers.!
2599 ! The unit of rei(k), fer(k), and fdr(k) is [Pa-1]. Alternative choice of !
2600 ! 'rei(k)' is also shown below, where coefficient 0.5 was from approximate !
2601 ! tuning against the BOMEX case. !
2602 ! In order to prevent the onset of instability in association with cumulus !
2603 ! induced subsidence advection, cumulus mass flux at the top interface in !
2604 ! any layer should be smaller than ( 90% of ) total mass within that layer.!
2605 ! I imposed limits on 'rei(k)' as below, in such that stability condition !
2606 ! is always satisfied. !
2607 ! Below limiter of 'rei(k)' becomes negative for some cases, causing error.!
2608 ! So, for the time being, I came back to the original limiter. !
2609 ! ------------------------------------------------------------------------ !
2611 ud2 = 1._r8 - 2._r8*xc + xc**2
2612 ! rei(k) = ( rkm / scaleh / g / rho0j ) ! Default.
2613 rei(k) = ( 0.5_r8 * rkm / z0(k) / g /rho0j ) ! Alternative.
2614 if( xc .gt. 0.5_r8 ) rei(k) = min(rei(k),0.9_r8*log(dp0(k)/g/dt/umf(km1) + 1._r8)/dpe/(2._r8*xc-1._r8))
2615 fer(k) = rei(k) * ee2
2616 fdr(k) = rei(k) * ud2
2618 ! ------------------------------------------------------------------------------ !
2619 ! Iteration Start due to 'maxufrc' constraint [ ****************************** ] !
2620 ! ------------------------------------------------------------------------------ !
2622 ! -------------------------------------------------------------------------- !
2623 ! Calculate cumulus updraft mass flux and penetrative entrainment mass flux. !
2624 ! Note that non-zero penetrative entrainment mass flux will be asigned only !
2625 ! to interfaces from the top interface of 'kbup' layer to the base interface !
2626 ! of 'kpen' layer as will be shown later. !
2627 ! -------------------------------------------------------------------------- !
2629 umf(k) = umf(km1) * exp( dpe * ( fer(k) - fdr(k) ) )
2632 ! --------------------------------------------------------- !
2633 ! Compute cumulus updraft properties at the top interface. !
2634 ! Also use Tayler expansion in order to treat limiting case !
2635 ! --------------------------------------------------------- !
2637 if( fer(k)*dpe .lt. 1.e-4_r8 ) then
2638 thlu(k) = thlu(km1) + ( thle + ssthl0(k) * dpe / 2._r8 - thlu(km1) ) * fer(k) * dpe
2639 qtu(k) = qtu(km1) + ( qte + ssqt0(k) * dpe / 2._r8 - qtu(km1) ) * fer(k) * dpe
2640 uu(k) = uu(km1) + ( ue + ssu0(k) * dpe / 2._r8 - uu(km1) ) * fer(k) * dpe - PGFc * ssu0(k) * dpe
2641 vu(k) = vu(km1) + ( ve + ssv0(k) * dpe / 2._r8 - vu(km1) ) * fer(k) * dpe - PGFc * ssv0(k) * dpe
2643 tru(k,m) = tru(km1,m) + ( tre(m) + sstr0(k,m) * dpe / 2._r8 - tru(km1,m) ) * fer(k) * dpe
2646 thlu(k) = ( thle + ssthl0(k) / fer(k) - ssthl0(k) * dpe / 2._r8 ) - &
2647 ( thle + ssthl0(k) * dpe / 2._r8 - thlu(km1) + ssthl0(k) / fer(k) ) * exp(-fer(k) * dpe)
2648 qtu(k) = ( qte + ssqt0(k) / fer(k) - ssqt0(k) * dpe / 2._r8 ) - &
2649 ( qte + ssqt0(k) * dpe / 2._r8 - qtu(km1) + ssqt0(k) / fer(k) ) * exp(-fer(k) * dpe)
2650 uu(k) = ( ue + ( 1._r8 - PGFc ) * ssu0(k) / fer(k) - ssu0(k) * dpe / 2._r8 ) - &
2651 ( ue + ssu0(k) * dpe / 2._r8 - uu(km1) + ( 1._r8 - PGFc ) * ssu0(k) / fer(k) ) * exp(-fer(k) * dpe)
2652 vu(k) = ( ve + ( 1._r8 - PGFc ) * ssv0(k) / fer(k) - ssv0(k) * dpe / 2._r8 ) - &
2653 ( ve + ssv0(k) * dpe / 2._r8 - vu(km1) + ( 1._r8 - PGFc ) * ssv0(k) / fer(k) ) * exp(-fer(k) * dpe)
2655 tru(k,m) = ( tre(m) + sstr0(k,m) / fer(k) - sstr0(k,m) * dpe / 2._r8 ) - &
2656 ( tre(m) + sstr0(k,m) * dpe / 2._r8 - tru(km1,m) + sstr0(k,m) / fer(k) ) * exp(-fer(k) * dpe)
2660 !------------------------------------------------------------------- !
2661 ! Expel some of cloud water and ice from cumulus updraft at the top !
2662 ! interface. Note that this is not 'detrainment' term but a 'sink' !
2663 ! term of cumulus updraft qt ( or one part of 'source' term of mean !
2664 ! environmental qt ). At this stage, as the most simplest choice, if !
2665 ! condensate amount within cumulus updraft is larger than a critical !
2666 ! value, 'criqc', expels the surplus condensate from cumulus updraft !
2667 ! to the environment. A certain fraction ( e.g., 'frc_sus' ) of this !
2668 ! expelled condesnate will be in a form that can be suspended in the !
2669 ! layer k where it was formed, while the other fraction, '1-frc_sus' !
2670 ! will be in a form of precipitatble (e.g.,can potentially fall down !
2671 ! across the base interface of layer k ). In turn we should describe !
2672 ! subsequent falling of precipitable condensate ('1-frc_sus') across !
2673 ! the base interface of the layer k, & evaporation of precipitating !
2674 ! water in the below layer k-1 and associated evaporative cooling of !
2675 ! the later, k-1, and falling of 'non-evaporated precipitating water !
2676 ! ( which was initially formed in layer k ) and a newly-formed preci !
2677 ! pitable water in the layer, k-1', across the base interface of the !
2678 ! lower layer k-1. Cloud microphysics should correctly describe all !
2679 ! of these process. In a near future, I should significantly modify !
2680 ! this cloud microphysics, including precipitation-induced downdraft !
2682 ! ------------------------------------------------------------------ !
2684 call conden(ps0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check,qsat)
2685 if( id_check .eq. 1 ) then
2686 exit_conden(i) = 1._r8
2690 if( (qlj + qij) .gt. criqc ) then
2691 exql = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij )
2692 exqi = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij )
2693 ! ---------------------------------------------------------------- !
2694 ! It is very important to re-update 'qtu' and 'thlu' at the upper !
2695 ! interface after expelling condensate from cumulus updraft at the !
2696 ! top interface of the layer. As mentioned above, this is a 'sink' !
2697 ! of cumulus qt (or equivalently, a 'source' of environmentasl qt),!
2698 ! not a regular convective'detrainment'. !
2699 ! ---------------------------------------------------------------- !
2700 qtu(k) = qtu(k) - exql - exqi
2701 thlu(k) = thlu(k) + (xlv/cp/exns0(k))*exql + (xls/cp/exns0(k))*exqi
2702 ! ---------------------------------------------------------------- !
2703 ! Expelled cloud condensate into the environment from the updraft. !
2704 ! After all the calculation later, 'dwten' and 'diten' will have a !
2705 ! unit of [ kg/kg/s ], because it is a tendency of qt. Restoration !
2706 ! of 'dwten' and 'diten' to this correct unit through multiplying !
2707 ! 'umf(k)*g/dp0(k)' will be performed later after finally updating !
2708 ! 'umf' using a 'rmaxfrac' constraint near the end of this updraft !
2709 ! buoyancy sorting loop. !
2710 ! ---------------------------------------------------------------- !
2717 ! ----------------------------------------------------------------- !
2718 ! Update 'thvu(k)' after detraining condensate from cumulus updraft.!
2719 ! ----------------------------------------------------------------- !
2720 call conden(ps0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check,qsat)
2721 if( id_check .eq. 1 ) then
2722 exit_conden(i) = 1._r8
2726 thvu(k) = thj * ( 1._r8 + zvir * qvj - qlj - qij )
2728 ! ----------------------------------------------------------- !
2729 ! Calculate updraft vertical velocity at the upper interface. !
2730 ! In order to calculate 'wtw' at the upper interface, we use !
2731 ! 'wtw' at the lower interface. Note 'wtw' is continuously !
2732 ! updated as cumulus updraft rises. !
2733 ! ----------------------------------------------------------- !
2735 bogbot = rbuoy * ( thvu(km1) / thvebot - 1._r8 ) ! Cloud buoyancy at base interface
2736 bogtop = rbuoy * ( thvu(k) / thv0top(k) - 1._r8 ) ! Cloud buoyancy at top interface
2738 delbog = bogtop - bogbot
2739 drage = fer(k) * ( 1._r8 + rdrag )
2740 expfac = exp(-2._r8*drage*dpe)
2743 if( drage*dpe .gt. 1.e-3_r8 ) then
2744 wtw = wtw*expfac + (delbog + (1._r8-expfac)*(bogbot + delbog/(-2._r8*drage*dpe)))/(rho0j*drage)
2746 wtw = wtw + dpe * ( bogbot + bogtop ) / rho0j
2749 ! Force the plume rise at least to klfc of the undiluted plume.
2750 ! Because even the below is not complete, I decided not to include this.
2752 ! if( k .le. klfc ) then
2753 ! wtw = max( 1.e-2_r8, wtw )
2756 ! -------------------------------------------------------------- !
2757 ! Repeat 'iter_xc' iteration loop until 'iter_xc = niter_xc'. !
2758 ! Also treat the case even when wtw < 0 at the 'kpen' interface. !
2759 ! -------------------------------------------------------------- !
2761 if( wtw .gt. 0._r8 ) then
2762 thlue = 0.5_r8 * ( thlu(km1) + thlu(k) )
2763 qtue = 0.5_r8 * ( qtu(km1) + qtu(k) )
2764 wue = 0.5_r8 * sqrt( max( wtwb + wtw, 0._r8 ) )
2769 enddo ! End of 'iter_xc' loop
2773 ! --------------------------------------------------------------------------- !
2774 ! Add the contribution of self-detrainment to vertical variations of cumulus !
2775 ! updraft mass flux. The reason why we are trying to include self-detrainment !
2776 ! is as follows. In current scheme, vertical variation of updraft mass flux !
2777 ! is not fully consistent with the vertical variation of updraft vertical w. !
2778 ! For example, within a given layer, let's assume that cumulus w is positive !
2779 ! at the base interface, while negative at the top interface. This means that !
2780 ! cumulus updraft cannot reach to the top interface of the layer. However, !
2781 ! cumulus updraft mass flux at the top interface is not zero according to the !
2782 ! vertical tendency equation of cumulus mass flux. Ideally, cumulus updraft !
2783 ! mass flux at the top interface should be zero for this case. In order to !
2784 ! assures that cumulus updraft mass flux goes to zero when cumulus updraft !
2785 ! vertical velocity goes to zero, we are imposing self-detrainment term as !
2786 ! below by considering layer-mean cloud buoyancy and cumulus updraft vertical !
2787 ! velocity square at the top interface. Use of auto-detrainment term will be !
2788 ! determined by setting 'use_self_detrain=.true.' in the parameter sentence. !
2789 ! --------------------------------------------------------------------------- !
2791 if( use_self_detrain ) then
2792 autodet = min( 0.5_r8*g*(bogbot+bogtop)/(max(wtw,0._r8)+1.e-4_r8), 0._r8 )
2793 umf(k) = umf(k) * exp( 0.637_r8*(dpe/rho0j/g) * autodet )
2795 if( umf(k) .eq. 0._r8 ) wtw = -1._r8
2797 ! -------------------------------------- !
2798 ! Below block is just a dignostic output !
2799 ! -------------------------------------- !
2801 excessu_arr(k) = excessu
2802 excess0_arr(k) = excess0
2804 aquad_arr(k) = aquad
2805 bquad_arr(k) = bquad
2806 cquad_arr(K) = cquad
2807 bogbot_arr(k) = bogbot
2808 bogtop_arr(k) = bogtop
2810 ! ------------------------------------------------------------------- !
2811 ! 'kbup' is the upper most layer in which cloud buoyancy is positive !
2812 ! both at the base and top interface. 'kpen' is the upper most layer !
2813 ! up to cumulus can reach. Usually, 'kpen' is located higher than the !
2814 ! 'kbup'. Note we initialized these by 'kbup = krel' & 'kpen = krel'. !
2815 ! As explained before, it is possible that only 'kpen' is updated, !
2816 ! while 'kbup' keeps its initialization value. For this case, current !
2817 ! scheme will simply turns-off penetrative entrainment fluxes and use !
2818 ! normal buoyancy-sorting fluxes for 'kbup <= k <= kpen-1' interfaces,!
2819 ! in order to describe shallow continental cumulus convection. !
2820 ! ------------------------------------------------------------------- !
2822 ! if( bogbot .gt. 0._r8 .and. bogtop .gt. 0._r8 ) then
2823 ! if( bogtop .gt. 0._r8 ) then
2824 if( bogtop .gt. 0._r8 .and. wtw .gt. 0._r8 ) then
2828 if( wtw .le. 0._r8 ) then
2834 if( wu(k) .gt. 100._r8 ) then
2840 ! ---------------------------------------------------------------------------- !
2841 ! Iteration end due to 'rmaxfrac' constraint [ ***************************** ] !
2842 ! ---------------------------------------------------------------------------- !
2844 ! ---------------------------------------------------------------------- !
2845 ! Calculate updraft fractional area at the upper interface and set upper !
2846 ! limit to 'ufrc' by 'rmaxfrac'. In order to keep the consistency among !
2847 ! ['ufrc','umf','wu (or wtw)'], if ufrc is limited by 'rmaxfrac', either !
2848 ! 'umf' or 'wu' should be changed. Although both 'umf' and 'wu (wtw)' at !
2849 ! the current upper interface are used for updating 'umf' & 'wu' at the !
2850 ! next upper interface, 'umf' is a passive variable not influencing the !
2851 ! buoyancy sorting process in contrast to 'wtw'. This is a reason why we !
2852 ! adjusted 'umf' instead of 'wtw'. In turn we updated 'fdr' here instead !
2853 ! of 'fer', which guarantees that all previously updated thermodynamic !
2854 ! variables at the upper interface before applying 'rmaxfrac' constraint !
2855 ! are already internally consistent, even though 'ufrc' is limited by !
2856 ! 'rmaxfrac'. Thus, we don't need to go through interation loop again.If !
2857 ! If we update 'fer' however, we should go through above iteration loop. !
2858 ! ---------------------------------------------------------------------- !
2860 rhos0j = ps0(k) / ( r * 0.5_r8 * ( thv0bot(k+1) + thv0top(k) ) * exns0(k) )
2861 ufrc(k) = umf(k) / ( rhos0j * wu(k) )
2862 if( ufrc(k) .gt. rmaxfrac ) then
2863 limit_ufrc(i) = 1._r8
2865 umf(k) = rmaxfrac * rhos0j * wu(k)
2866 fdr(k) = fer(k) - log( umf(k) / umf(km1) ) / dpe
2869 ! ------------------------------------------------------------ !
2870 ! Update environmental properties for at the mid-point of next !
2871 ! upper layer for use in buoyancy sorting. !
2872 ! ------------------------------------------------------------ !
2877 thvebot = thv0bot(k+1)
2886 end do ! End of cumulus updraft loop from the 'krel' layer to 'kpen' layer.
2888 ! ------------------------------------------------------------------------------- !
2889 ! Up to this point, we finished all of buoyancy sorting processes from the 'krel' !
2890 ! layer to 'kpen' layer: at the top interface of individual layers, we calculated !
2891 ! updraft and penetrative mass fluxes [ umf(k) & emf(k) = 0 ], updraft fractional !
2892 ! area [ ufrc(k) ], updraft vertical velocity [ wu(k) ], updraft thermodynamic !
2893 ! variables [thlu(k),qtu(k),uu(k),vu(k),thvu(k)]. In the layer,we also calculated !
2894 ! fractional entrainment-detrainment rate [ fer(k), fdr(k) ], and detrainment ten !
2895 ! dency of water and ice from cumulus updraft [ dwten(k), diten(k) ]. In addition,!
2896 ! we updated and identified 'krel' and 'kpen' layer index, if any. In the 'kpen' !
2897 ! layer, we calculated everything mentioned above except the 'wu(k)' and 'ufrc(k)'!
2898 ! since a real value of updraft vertical velocity is not defined at the kpen top !
2899 ! interface (note 'ufrc' at the top interface of layer is calculated from 'umf(k)'!
2900 ! and 'wu(k)'). As mentioned before, special treatment is required when 'kbup' is !
2901 ! not updated and so 'kbup = krel'. !
2902 ! ------------------------------------------------------------------------------- !
2904 ! ------------------------------------------------------------------------------ !
2905 ! During the 'iter_scaleh' iteration loop, non-physical ( with non-zero values ) !
2906 ! values can remain in the variable arrays above (also 'including' in case of wu !
2907 ! and ufrc at the top interface) the 'kpen' layer. This can happen when the kpen !
2908 ! layer index identified from the 'iter_scaleh = 1' iteration loop is located at !
2909 ! above the kpen layer index identified from 'iter_scaleh = 3' iteration loop. !
2910 ! Thus, in the following calculations, we should only use the values in each !
2911 ! variables only up to finally identified 'kpen' layer & 'kpen' interface except !
2912 ! 'wu' and 'ufrc' at the top interface of 'kpen' layer. Note that in order to !
2913 ! prevent any problems due to these non-physical values, I re-initialized the !
2914 ! values of [ umf(kpen:mkx), emf(kpen:mkx), dwten(kpen+1:mkx), diten(kpen+1:mkx),!
2915 ! fer(kpen:mkx), fdr(kpen+1:mkx), ufrc(kpen:mkx) ] to be zero after 'iter_scaleh'!
2917 ! ------------------------------------------------------------------------------ !
2921 ! ------------------------------------------------------------------------------ !
2922 ! Calculate 'ppen( < 0 )', updarft penetrative distance from the lower interface !
2923 ! of 'kpen' layer. Note that bogbot & bogtop at the 'kpen' layer either when fer !
2924 ! is zero or non-zero was already calculated above. !
2925 ! It seems that below qudarature solving formula is valid only when bogbot < 0. !
2926 ! Below solving equation is clearly wrong ! I should revise this ! !
2927 ! ------------------------------------------------------------------------------ !
2929 if( drage .eq. 0._r8 ) then
2930 aquad = ( bogtop - bogbot ) / ( ps0(kpen) - ps0(kpen-1) )
2931 bquad = 2._r8 * bogbot
2932 cquad = -wu(kpen-1)**2 * rho0j
2933 call roots(aquad,bquad,cquad,xc1,xc2,status)
2934 if( status .eq. 0 ) then
2935 if( xc1 .le. 0._r8 .and. xc2 .le. 0._r8 ) then
2936 ppen = max( xc1, xc2 )
2937 ppen = min( 0._r8,max( -dp0(kpen), ppen ) )
2938 elseif( xc1 .gt. 0._r8 .and. xc2 .gt. 0._r8 ) then
2940 write(iulog,*) 'Warning : UW-Cumulus penetrates upto kpen interface'
2942 call wrf_message(iulog)
2945 ppen = min( xc1, xc2 )
2946 ppen = min( 0._r8,max( -dp0(kpen), ppen ) )
2950 write(iulog,*) 'Warning : UW-Cumulus penetrates upto kpen interface'
2952 call wrf_message(iulog)
2956 ppen = compute_ppen(wtwb,drage,bogbot,bogtop,rho0j,dp0(kpen))
2958 if( ppen .eq. -dp0(kpen) .or. ppen .eq. 0._r8 ) limit_ppen(i) = 1._r8
2960 ! -------------------------------------------------------------------- !
2961 ! Re-calculate the amount of expelled condensate from cloud updraft !
2962 ! at the cumulus top. This is necessary for refined calculations of !
2963 ! bulk cloud microphysics at the cumulus top. Note that ppen < 0._r8 !
2964 ! In the below, I explicitly calculate 'thlu_top' & 'qtu_top' by !
2965 ! using non-zero 'fer(kpen)'. !
2966 ! -------------------------------------------------------------------- !
2968 if( fer(kpen)*(-ppen) .lt. 1.e-4_r8 ) then
2969 thlu_top = thlu(kpen-1) + ( thl0(kpen) + ssthl0(kpen) * (-ppen) / 2._r8 - thlu(kpen-1) ) * fer(kpen) * (-ppen)
2970 qtu_top = qtu(kpen-1) + ( qt0(kpen) + ssqt0(kpen) * (-ppen) / 2._r8 - qtu(kpen-1) ) * fer(kpen) * (-ppen)
2972 thlu_top = ( thl0(kpen) + ssthl0(kpen) / fer(kpen) - ssthl0(kpen) * (-ppen) / 2._r8 ) - &
2973 ( thl0(kpen) + ssthl0(kpen) * (-ppen) / 2._r8 - thlu(kpen-1) + ssthl0(kpen) / fer(kpen) ) * &
2974 exp(-fer(kpen) * (-ppen))
2975 qtu_top = ( qt0(kpen) + ssqt0(kpen) / fer(kpen) - ssqt0(kpen) * (-ppen) / 2._r8 ) - &
2976 ( qt0(kpen) + ssqt0(kpen) * (-ppen) / 2._r8 - qtu(kpen-1) + ssqt0(kpen) / fer(kpen) ) * &
2977 exp(-fer(kpen) * (-ppen))
2980 call conden(ps0(kpen-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check,qsat)
2981 if( id_check .eq. 1 ) then
2982 exit_conden(i) = 1._r8
2986 exntop = ((ps0(kpen-1)+ppen)/p00)**rovcp
2987 if( (qlj + qij) .gt. criqc ) then
2988 dwten(kpen) = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij )
2989 diten(kpen) = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij )
2990 qtu_top = qtu_top - dwten(kpen) - diten(kpen)
2991 thlu_top = thlu_top + (xlv/cp/exntop)*dwten(kpen) + (xls/cp/exntop)*diten(kpen)
2997 ! ----------------------------------------------------------------------- !
2998 ! Calculate cumulus scale height as the top height that cumulus can reach.!
2999 ! ----------------------------------------------------------------------- !
3001 rhos0j = ps0(kpen-1)/(r*0.5_r8*(thv0bot(kpen)+thv0top(kpen-1))*exns0(kpen-1))
3002 cush = zs0(kpen-1) - ppen/rhos0j/g
3005 end do ! End of 'iter_scaleh' loop.
3007 ! -------------------------------------------------------------------- !
3008 ! The 'forcedCu' is logical identifier saying whether cumulus updraft !
3009 ! overcome the buoyancy barrier just above the PBL top. If it is true, !
3010 ! cumulus did not overcome the barrier - this is a shallow convection !
3011 ! with negative cloud buoyancy, mimicking shallow continental cumulus !
3012 ! convection. Depending on 'forcedCu' parameter, treatment of heat & !
3013 ! moisture fluxes at the entraining interfaces, 'kbup <= k < kpen - 1' !
3014 ! will be set up in a different ways, as will be shown later. !
3015 ! -------------------------------------------------------------------- !
3017 if( kbup .eq. krel ) then
3019 limit_shcu(i) = 1._r8
3022 limit_shcu(i) = 0._r8
3025 ! ------------------------------------------------------------------ !
3026 ! Filtering of unerasonable cumulus adjustment here. This is a very !
3027 ! important process which should be done cautiously. Various ways of !
3028 ! filtering are possible depending on cases mainly using the indices !
3029 ! of key layers - 'klcl','kinv','krel','klfc','kbup','kpen'. At this !
3030 ! stage, the followings are all possible : 'kinv >= 2', 'klcl >= 1', !
3031 ! 'krel >= kinv', 'kbup >= krel', 'kpen >= krel'. I must design this !
3032 ! filtering very cautiously, in such that none of realistic cumulus !
3033 ! convection is arbitrarily turned-off. Potentially, I might turn-off!
3034 ! cumulus convection if layer-mean 'ql > 0' in the 'kinv-1' layer,in !
3035 ! order to suppress cumulus convection growing, based at the Sc top. !
3036 ! This is one of potential future modifications. Note that ppen < 0. !
3037 ! ------------------------------------------------------------------ !
3039 cldhgt = ps0(kpen-1) + ppen
3041 ! write(iulog,*) 'forcedCu - did not overcome initial buoyancy barrier'
3042 exit_cufilter(i) = 1._r8
3046 ! Limit 'additional shallow cumulus' for DYCOMS simulation.
3047 ! if( cldhgt.ge.88000._r8 ) then
3052 ! ------------------------------------------------------------------------------ !
3053 ! Re-initializing some key variables above the 'kpen' layer in order to suppress !
3054 ! the influence of non-physical values above 'kpen', in association with the use !
3055 ! of 'iter_scaleh' loop. Note that umf, emf, ufrc are defined at the interfaces !
3056 ! (0:mkx), while 'dwten','diten', 'fer', 'fdr' are defined at layer mid-points. !
3057 ! Initialization of 'fer' and 'fdr' is for correct writing purpose of diagnostic !
3058 ! output. Note that we set umf(kpen)=emf(kpen)=ufrc(kpen)=0, in consistent with !
3059 ! wtw < 0 at the top interface of 'kpen' layer. However, we still have non-zero !
3060 ! expelled cloud condensate in the 'kpen' layer. !
3061 ! ------------------------------------------------------------------------------ !
3063 umf(kpen:mkx) = 0._r8
3064 emf(kpen:mkx) = 0._r8
3065 ufrc(kpen:mkx) = 0._r8
3066 dwten(kpen+1:mkx) = 0._r8
3067 diten(kpen+1:mkx) = 0._r8
3068 fer(kpen+1:mkx) = 0._r8
3069 fdr(kpen+1:mkx) = 0._r8
3071 ! ------------------------------------------------------------------------ !
3072 ! Calculate downward penetrative entrainment mass flux, 'emf(k) < 0', and !
3073 ! thermodynamic properties of penetratively entrained airs at entraining !
3074 ! interfaces. emf(k) is defined from the top interface of the layer kbup !
3075 ! to the bottom interface of the layer 'kpen'. Note even when kbup = krel,!
3076 ! i.e.,even when 'kbup' was not updated in the above buoyancy sorting do !
3077 ! loop (i.e., 'kbup' remains as the initialization value), below do loop !
3078 ! of penetrative entrainment flux can be performed without any conceptual !
3079 ! or logical problems, because we have already computed all the variables !
3080 ! necessary for performing below penetrative entrainment block. !
3081 ! In the below 'do' loop, 'k' is an interface index at which non-zero 'emf'!
3082 ! (penetrative entrainment mass flux) is calculated. Since cumulus updraft !
3083 ! is negatively buoyant in the layers between the top interface of 'kbup' !
3084 ! layer (interface index, kbup) and the top interface of 'kpen' layer, the !
3085 ! fractional lateral entrainment, fer(k) within these layers will be close !
3086 ! to zero - so it is likely that only strong lateral detrainment occurs in !
3087 ! thses layers. Under this situation,we can easily calculate the amount of !
3088 ! detrainment cumulus air into these negatively buoyanct layers by simply !
3089 ! comparing cumulus updraft mass fluxes between the base and top interface !
3090 ! of each layer: emf(k) = emf(k-1)*exp(-fdr(k)*dp0(k)) !
3091 ! ~ emf(k-1)*(1-rei(k)*dp0(k)) !
3092 ! emf(k-1)-emf(k) ~ emf(k-1)*rei(k)*dp0(k) !
3093 ! Current code assumes that about 'rpen~10' times of these detrained mass !
3094 ! are penetratively re-entrained down into the 'k-1' interface. And all of !
3095 ! these detrained masses are finally dumped down into the top interface of !
3096 ! 'kbup' layer. Thus, the amount of penetratively entrained air across the !
3097 ! top interface of 'kbup' layer with 'rpen~10' becomes too large. !
3098 ! Note that this penetrative entrainment part can be completely turned-off !
3099 ! and we can simply use normal buoyancy-sorting involved turbulent fluxes !
3100 ! by modifying 'penetrative entrainment fluxes' part below. !
3101 ! ------------------------------------------------------------------------ !
3103 ! -----------------------------------------------------------------------!
3104 ! Calculate entrainment mass flux and conservative scalars of entraining !
3105 ! free air at interfaces of 'kbup <= k < kpen - 1' !
3106 ! ---------------------------------------------------------------------- !
3109 thlu_emf(k) = thlu(k)
3114 tru_emf(k,m) = tru(k,m)
3118 do k = kpen - 1, kbup, -1 ! Here, 'k' is an interface index at which
3119 ! penetrative entrainment fluxes are calculated.
3121 rhos0j = ps0(k) / ( r * 0.5_r8 * ( thv0bot(k+1) + thv0top(k) ) * exns0(k) )
3123 if( k .eq. kpen - 1 ) then
3125 ! ------------------------------------------------------------------------ !
3126 ! Note that 'ppen' has already been calculated in the above 'iter_scaleh' !
3127 ! loop assuming zero lateral entrainmentin the layer 'kpen'. !
3128 ! ------------------------------------------------------------------------ !
3130 ! -------------------------------------------------------------------- !
3131 ! Calculate returning mass flux, emf ( < 0 ) !
3132 ! Current penetrative entrainment rate with 'rpen~10' is too large and !
3133 ! future refinement is necessary including the definition of 'thl','qt'!
3134 ! of penetratively entrained air. Penetratively entrained airs across !
3135 ! the 'kpen-1' interface is assumed to have the properties of the base !
3136 ! interface of 'kpen' layer. Note that 'emf ~ - umf/ufrc = - w * rho'. !
3137 ! Thus, below limit sets an upper limit of |emf| to be ~ 10cm/s, which !
3138 ! is very loose constraint. Here, I used more restricted constraint on !
3139 ! the limit of emf, assuming 'emf' cannot exceed a net mass within the !
3140 ! layer above the interface. Similar to the case of warming and drying !
3141 ! due to cumulus updraft induced compensating subsidence, penetrative !
3142 ! entrainment induces compensating upwelling - in order to prevent !
3143 ! numerical instability in association with compensating upwelling, we !
3144 ! should similarily limit the amount of penetrative entrainment at the !
3145 ! interface by the amount of masses within the layer just above the !
3146 ! penetratively entraining interface. !
3147 ! -------------------------------------------------------------------- !
3149 if( ( umf(k)*ppen*rei(kpen)*rpen ) .lt. -0.1_r8*rhos0j ) limit_emf(i) = 1._r8
3150 if( ( umf(k)*ppen*rei(kpen)*rpen ) .lt. -0.9_r8*dp0(kpen)/g/dt ) limit_emf(i) = 1._r8
3152 emf(k) = max( max( umf(k)*ppen*rei(kpen)*rpen, -0.1_r8*rhos0j), -0.9_r8*dp0(kpen)/g/dt)
3153 thlu_emf(k) = thl0(kpen) + ssthl0(kpen) * ( ps0(k) - p0(kpen) )
3154 qtu_emf(k) = qt0(kpen) + ssqt0(kpen) * ( ps0(k) - p0(kpen) )
3155 uu_emf(k) = u0(kpen) + ssu0(kpen) * ( ps0(k) - p0(kpen) )
3156 vu_emf(k) = v0(kpen) + ssv0(kpen) * ( ps0(k) - p0(kpen) )
3158 tru_emf(k,m) = tr0(kpen,m) + sstr0(kpen,m) * ( ps0(k) - p0(kpen) )
3161 else ! if(k.lt.kpen-1).
3163 ! --------------------------------------------------------------------------- !
3164 ! Note we are coming down from the higher interfaces to the lower interfaces. !
3165 ! Also note that 'emf < 0'. So, below operation is a summing not subtracting. !
3166 ! In order to ensure numerical stability, I imposed a modified correct limit !
3167 ! of '-0.9*dp0(k+1)/g/dt' on emf(k). !
3168 ! --------------------------------------------------------------------------- !
3170 if( use_cumpenent ) then ! Original Cumulative Penetrative Entrainment
3172 if( ( emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.1_r8*rhos0j ) limit_emf(i) = 1
3173 if( ( emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.9_r8*dp0(k+1)/g/dt ) limit_emf(i) = 1
3174 emf(k) = max(max(emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen, -0.1_r8*rhos0j), -0.9_r8*dp0(k+1)/g/dt )
3175 if( abs(emf(k)) .gt. abs(emf(k+1)) ) then
3176 thlu_emf(k) = ( thlu_emf(k+1) * emf(k+1) + thl0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k)
3177 qtu_emf(k) = ( qtu_emf(k+1) * emf(k+1) + qt0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k)
3178 uu_emf(k) = ( uu_emf(k+1) * emf(k+1) + u0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k)
3179 vu_emf(k) = ( vu_emf(k+1) * emf(k+1) + v0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k)
3181 tru_emf(k,m) = ( tru_emf(k+1,m) * emf(k+1) + tr0(k+1,m) * ( emf(k) - emf(k+1) ) ) / emf(k)
3184 thlu_emf(k) = thl0(k+1)
3185 qtu_emf(k) = qt0(k+1)
3189 tru_emf(k,m) = tr0(k+1,m)
3193 else ! Alternative Non-Cumulative Penetrative Entrainment
3195 if( ( -umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.1_r8*rhos0j ) limit_emf(i) = 1
3196 if( ( -umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.9_r8*dp0(k+1)/g/dt ) limit_emf(i) = 1
3197 emf(k) = max(max(-umf(k)*dp0(k+1)*rei(k+1)*rpen, -0.1_r8*rhos0j), -0.9_r8*dp0(k+1)/g/dt )
3198 thlu_emf(k) = thl0(k+1)
3199 qtu_emf(k) = qt0(k+1)
3203 tru_emf(k,m) = tr0(k+1,m)
3210 ! ---------------------------------------------------------------------------- !
3211 ! In this GCM modeling framework, all what we should do is to calculate heat !
3212 ! and moisture fluxes at the given geometrically-fixed height interfaces - we !
3213 ! don't need to worry about movement of material height surface in association !
3214 ! with compensating subsidence or unwelling, in contrast to the bulk modeling. !
3215 ! In this geometrically fixed height coordinate system, heat and moisture flux !
3216 ! at the geometrically fixed height handle everything - a movement of material !
3217 ! surface is implicitly treated automatically. Note that in terms of turbulent !
3218 ! heat and moisture fluxes at model interfaces, both the cumulus updraft mass !
3219 ! flux and penetratively entraining mass flux play the same role -both of them !
3220 ! warms and dries the 'kbup' layer, cools and moistens the 'kpen' layer, and !
3221 ! cools and moistens any intervening layers between 'kbup' and 'kpen' layers. !
3222 ! It is important to note these identical roles on turbulent heat and moisture !
3223 ! fluxes of 'umf' and 'emf'. !
3224 ! When 'kbup' is a stratocumulus-topped PBL top interface, increase of 'rpen' !
3225 ! is likely to strongly diffuse stratocumulus top interface, resulting in the !
3226 ! reduction of cloud fraction. In this sense, the 'kbup' interface has a very !
3227 ! important meaning and role : across the 'kbup' interface, strong penetrative !
3228 ! entrainment occurs, thus any sharp gradient properties across that interface !
3229 ! are easily diffused through strong mass exchange. Thus, an initialization of !
3230 ! 'kbup' (and also 'kpen') should be done very cautiously as mentioned before. !
3231 ! In order to prevent this stron diffusion for the shallow cumulus convection !
3232 ! based at the Sc top, it seems to be good to initialize 'kbup = krel', rather !
3233 ! that 'kbup = krel-1'. !
3234 ! ---------------------------------------------------------------------------- !
3238 !------------------------------------------------------------------ !
3240 ! Compute turbulent heat, moisture, momentum flux at all interfaces !
3242 !------------------------------------------------------------------ !
3243 ! It is very important to note that in calculating turbulent fluxes !
3244 ! below, we must not double count turbulent flux at any interefaces.!
3245 ! In the below, turbulent fluxes at the interfaces (interface index !
3246 ! k) are calculated by the following 4 blocks in consecutive order: !
3248 ! (1) " 0 <= k <= kinv - 1 " : PBL fluxes. !
3249 ! From 'fluxbelowinv' using reconstructed PBL height. Currently,!
3250 ! the reconstructed PBLs are independently calculated for each !
3251 ! individual conservative scalar variables ( qt, thl, u, v ) in !
3252 ! each 'fluxbelowinv', instead of being uniquely calculated by !
3253 ! using thvl. Turbulent flux at the surface is assumed to be 0. !
3254 ! (2) " kinv <= k <= krel - 1 " : Non-buoyancy sorting fluxes !
3255 ! Assuming cumulus mass flux and cumulus updraft thermodynamic !
3256 ! properties (except u, v which are modified by the PGFc during !
3257 ! upward motion) are conserved during a updraft motion from the !
3258 ! PBL top interface to the release level. If these layers don't !
3259 ! exist (e,g, when 'krel = kinv'), then current routine do not !
3260 ! perform this routine automatically. So I don't need to modify !
3262 ! (3) " krel <= k <= kbup - 1 " : Buoyancy sorting fluxes !
3263 ! From laterally entraining-detraining buoyancy sorting plumes. !
3264 ! (4) " kbup <= k < kpen-1 " : Penetrative entrainment fluxes !
3265 ! From penetratively entraining plumes, !
3267 ! In case of normal situation, turbulent interfaces in each groups !
3268 ! are mutually independent of each other. Thus double flux counting !
3269 ! or ambiguous flux counting requiring the choice among the above 4 !
3270 ! groups do not occur normally. However, in case that cumulus plume !
3271 ! could not completely overcome the buoyancy barrier just above the !
3272 ! PBL top interface and so 'kbup = krel' (.forcedCu=.true.) ( here, !
3273 ! it can be either 'kpen = krel' as the initialization, or ' kpen > !
3274 ! krel' if cumulus updraft just penetrated over the top of release !
3275 ! layer ). If this happens, we should be very careful in organizing !
3276 ! the sequence of the 4 calculation routines above - note that the !
3277 ! routine located at the later has the higher priority. Additional !
3278 ! feature I must consider is that when 'kbup = kinv - 1' (this is a !
3279 ! combined situation of 'kbup=krel-1' & 'krel = kinv' when I chose !
3280 ! 'kbup=krel-1' instead of current choice of 'kbup=krel'), a strong !
3281 ! penetrative entrainment fluxes exists at the PBL top interface, & !
3282 ! all of these fluxes are concentrated (deposited) within the layer !
3283 ! just below PBL top interface (i.e., 'kinv-1' layer). On the other !
3284 ! hand, in case of 'fluxbelowinv', only the compensating subsidence !
3285 ! effect is concentrated in the 'kinv-1' layer and 'pure' turbulent !
3286 ! heat and moisture fluxes ( 'pure' means the fluxes not associated !
3287 ! with compensating subsidence) are linearly distributed throughout !
3288 ! the whole PBL. Thus different choice of the above flux groups can !
3289 ! produce very different results. Output variable should be written !
3290 ! consistently to the choice of computation sequences. !
3291 ! When the case of 'kbup = krel(-1)' happens,another way to dealing !
3292 ! with this case is to simply ' exit ' the whole cumulus convection !
3293 ! calculation without performing any cumulus convection. We can !
3294 ! choose this approach by specifying a condition in the 'Filtering !
3295 ! of unreasonable cumulus adjustment' just after 'iter_scaleh'. But !
3296 ! this seems not to be a good choice (although this choice was used !
3297 ! previous code ), since it might arbitrary damped-out the shallow !
3298 ! cumulus convection over the continent land, where shallow cumulus !
3299 ! convection tends to be negatively buoyant. !
3300 ! ----------------------------------------------------------------- !
3302 ! --------------------------------------------------- !
3303 ! 1. PBL fluxes : 0 <= k <= kinv - 1 !
3304 ! All the information necessary to reconstruct PBL !
3305 ! height are passed to 'fluxbelowinv'. !
3306 ! --------------------------------------------------- !
3310 xtop = qt0(kinv+1) + ssqt0(kinv+1) * ( ps0(kinv) - p0(kinv+1) )
3311 xbot = qt0(kinv-1) + ssqt0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
3312 call fluxbelowinv( cbmf, ps0(0:mkx), mkx, kinv, dt, xsrc, xmean, xtop, xbot, xflx )
3313 qtflx(0:kinv-1) = xflx(0:kinv-1)
3317 xtop = thl0(kinv+1) + ssthl0(kinv+1) * ( ps0(kinv) - p0(kinv+1) )
3318 xbot = thl0(kinv-1) + ssthl0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
3319 call fluxbelowinv( cbmf, ps0(0:mkx), mkx, kinv, dt, xsrc, xmean, xtop, xbot, xflx )
3320 slflx(0:kinv-1) = cp * exns0(0:kinv-1) * xflx(0:kinv-1)
3324 xtop = u0(kinv+1) + ssu0(kinv+1) * ( ps0(kinv) - p0(kinv+1) )
3325 xbot = u0(kinv-1) + ssu0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
3326 call fluxbelowinv( cbmf, ps0(0:mkx), mkx, kinv, dt, xsrc, xmean, xtop, xbot, xflx )
3327 uflx(0:kinv-1) = xflx(0:kinv-1)
3331 xtop = v0(kinv+1) + ssv0(kinv+1) * ( ps0(kinv) - p0(kinv+1) )
3332 xbot = v0(kinv-1) + ssv0(kinv-1) * ( ps0(kinv-1) - p0(kinv-1) )
3333 call fluxbelowinv( cbmf, ps0(0:mkx), mkx, kinv, dt, xsrc, xmean, xtop, xbot, xflx )
3334 vflx(0:kinv-1) = xflx(0:kinv-1)
3339 xtop = tr0(kinv+1,m) + sstr0(kinv+1,m) * ( ps0(kinv) - p0(kinv+1) )
3340 xbot = tr0(kinv-1,m) + sstr0(kinv-1,m) * ( ps0(kinv-1) - p0(kinv-1) )
3341 call fluxbelowinv( cbmf, ps0(0:mkx), mkx, kinv, dt, xsrc, xmean, xtop, xbot, xflx )
3342 trflx(0:kinv-1,m) = xflx(0:kinv-1)
3345 ! -------------------------------------------------------------- !
3346 ! 2. Non-buoyancy sorting fluxes : kinv <= k <= krel - 1 !
3347 ! Note that when 'krel = kinv', below block is never executed !
3348 ! as in a desirable, expected way ( but I must check if this !
3349 ! is the case ). The non-buoyancy sorting fluxes are computed !
3350 ! only when 'krel > kinv'. !
3351 ! -------------------------------------------------------------- !
3355 do k = kinv, krel - 1
3357 qtflx(k) = cbmf * ( qtsrc - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3358 slflx(k) = cbmf * ( thlsrc - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) ) * cp * exns0(k)
3359 uplus = uplus + PGFc * ssu0(k) * ( ps0(k) - ps0(k-1) )
3360 vplus = vplus + PGFc * ssv0(k) * ( ps0(k) - ps0(k-1) )
3361 uflx(k) = cbmf * ( usrc + uplus - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3362 vflx(k) = cbmf * ( vsrc + vplus - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3364 trflx(k,m) = cbmf * ( trsrc(m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) )
3368 ! ------------------------------------------------------------------------ !
3369 ! 3. Buoyancy sorting fluxes : krel <= k <= kbup - 1 !
3370 ! In case that 'kbup = krel - 1 ' ( or even in case 'kbup = krel' ), !
3371 ! buoyancy sorting fluxes are not calculated, which is consistent, !
3372 ! desirable feature. !
3373 ! ------------------------------------------------------------------------ !
3375 do k = krel, kbup - 1
3377 slflx(k) = cp * exns0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3378 qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3379 uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3380 vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) )
3382 trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) )
3386 ! ------------------------------------------------------------------------- !
3387 ! 4. Penetrative entrainment fluxes : kbup <= k <= kpen - 1 !
3388 ! The only confliction that can happen is when 'kbup = kinv-1'. For this !
3389 ! case, turbulent flux at kinv-1 is calculated both from 'fluxbelowinv' !
3390 ! and here as penetrative entrainment fluxes. Since penetrative flux is !
3391 ! calculated later, flux at 'kinv - 1 ' will be that of penetrative flux.!
3392 ! However, turbulent flux calculated at 'kinv - 1' from penetrative entr.!
3393 ! is less attractable, since more reasonable turbulent flux at 'kinv-1' !
3394 ! should be obtained from 'fluxbelowinv', by considering re-constructed !
3395 ! inversion base height. This conflicting problem can be solved if we can!
3396 ! initialize 'kbup = krel', instead of kbup = krel - 1. This choice seems!
3397 ! to be more reasonable since it is not conflicted with 'fluxbelowinv' in!
3398 ! calculating fluxes at 'kinv - 1' ( for this case, flux at 'kinv-1' is !
3399 ! always from 'fluxbelowinv' ), and flux at 'krel-1' is calculated from !
3400 ! the non-buoyancy sorting flux without being competed with penetrative !
3401 ! entrainment fluxes. Even when we use normal cumulus flux instead of !
3402 ! penetrative entrainment fluxes at 'kbup <= k <= kpen-1' interfaces, !
3403 ! the initialization of kbup=krel perfectly works without any conceptual !
3404 ! confliction. Thus it seems to be much better to choose 'kbup = krel' !
3405 ! initialization of 'kbup', which is current choice. !
3406 ! Note that below formula uses conventional updraft cumulus fluxes for !
3407 ! shallow cumulus which did not overcome the first buoyancy barrier above!
3408 ! PBL top while uses penetrative entrainment fluxes for the other cases !
3409 ! 'kbup <= k <= kpen-1' interfaces. Depending on cases, however, I can !
3410 ! selelct different choice. !
3411 ! ------------------------------------------------------------------------------------------------------------------ !
3412 ! if( forcedCu ) then !
3413 ! slflx(k) = cp * exns0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) ) !
3414 ! qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) ) !
3415 ! uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) ) !
3416 ! vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) ) !
3418 ! trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) ) !
3421 ! slflx(k) = cp * exns0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) ) ) ) !
3422 ! qtflx(k) = emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) ) ) ) !
3423 ! uflx(k) = emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( ps0(k) - p0(k) ) ) ) !
3424 ! vflx(k) = emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( ps0(k) - p0(k) ) ) ) !
3426 ! trflx(k,m) = emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( ps0(k) - p0(k) ) ) ) !
3430 ! if( use_uppenent ) then ! Combined Updraft + Penetrative Entrainment Flux !
3431 ! slflx(k) = cp * exns0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & !
3432 ! cp * exns0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) ) ) ) !
3433 ! qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & !
3434 ! emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) ) ) ) !
3435 ! uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & !
3436 ! emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( ps0(k) - p0(k) ) ) ) !
3437 ! vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & !
3438 ! emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( ps0(k) - p0(k) ) ) ) !
3440 ! trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) ) + & !
3441 ! emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( ps0(k) - p0(k) ) ) ) !
3443 ! ------------------------------------------------------------------------------------------------------------------ !
3445 do k = kbup, kpen - 1
3447 slflx(k) = cp * exns0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) ) ) )
3448 qtflx(k) = emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) ) ) )
3449 uflx(k) = emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( ps0(k) - p0(k) ) ) )
3450 vflx(k) = emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( ps0(k) - p0(k) ) ) )
3452 trflx(k,m) = emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( ps0(k) - p0(k) ) ) )
3456 ! ------------------------------------------- !
3457 ! Turn-off cumulus momentum flux as an option !
3458 ! ------------------------------------------- !
3460 if( .not. use_momenflx ) then
3465 ! -------------------------------------------------------- !
3466 ! Condensate tendency by compensating subsidence/upwelling !
3467 ! -------------------------------------------------------- !
3470 do k = 0, kinv - 2 ! Assume linear updraft mass flux within the PBL.
3471 uemf(k) = cbmf * ( ps0(0) - ps0(k) ) / ( ps0(0) - ps0(kinv-1) )
3473 uemf(kinv-1:krel-1) = cbmf
3474 uemf(krel:kbup-1) = umf(krel:kbup-1)
3475 uemf(kbup:kpen-1) = emf(kbup:kpen-1) ! Only use penetrative entrainment flux consistently.
3477 comsub(1:mkx) = 0._r8
3479 comsub(k) = 0.5_r8 * ( uemf(k) + uemf(k-1) )
3483 if( comsub(k) .ge. 0._r8 ) then
3484 if( k .eq. mkx ) then
3492 thlten_sub = g * comsub(k) * ( thl0(k+1) - thl0(k) ) / ( p0(k) - p0(k+1) )
3493 qtten_sub = g * comsub(k) * ( qt0(k+1) - qt0(k) ) / ( p0(k) - p0(k+1) )
3494 qlten_sub = g * comsub(k) * ( ql0(k+1) - ql0(k) ) / ( p0(k) - p0(k+1) )
3495 qiten_sub = g * comsub(k) * ( qi0(k+1) - qi0(k) ) / ( p0(k) - p0(k+1) )
3496 nlten_sub = g * comsub(k) * ( tr0(k+1,ixnumliq) - tr0(k,ixnumliq) ) / ( p0(k) - p0(k+1) )
3497 niten_sub = g * comsub(k) * ( tr0(k+1,ixnumice) - tr0(k,ixnumice) ) / ( p0(k) - p0(k+1) )
3508 thlten_sub = g * comsub(k) * ( thl0(k) - thl0(k-1) ) / ( p0(k-1) - p0(k) )
3509 qtten_sub = g * comsub(k) * ( qt0(k) - qt0(k-1) ) / ( p0(k-1) - p0(k) )
3510 qlten_sub = g * comsub(k) * ( ql0(k) - ql0(k-1) ) / ( p0(k-1) - p0(k) )
3511 qiten_sub = g * comsub(k) * ( qi0(k) - qi0(k-1) ) / ( p0(k-1) - p0(k) )
3512 nlten_sub = g * comsub(k) * ( tr0(k,ixnumliq) - tr0(k-1,ixnumliq) ) / ( p0(k-1) - p0(k) )
3513 niten_sub = g * comsub(k) * ( tr0(k,ixnumice) - tr0(k-1,ixnumice) ) / ( p0(k-1) - p0(k) )
3516 thl_prog = thl0(k) + thlten_sub * dt
3517 qt_prog = max( qt0(k) + qtten_sub * dt, 1.e-12_r8 )
3518 call conden(p0(k),thl_prog,qt_prog,thj,qvj,qlj,qij,qse,id_check,qsat)
3519 if( id_check .eq. 1 ) then
3523 ! qlten_sink(k) = ( qlj - ql0(k) ) / dt
3524 ! qiten_sink(k) = ( qij - qi0(k) ) / dt
3525 qlten_sink(k) = max( qlten_sub, - ql0(k) / dt ) ! For consistency with prognostic macrophysics scheme
3526 qiten_sink(k) = max( qiten_sub, - qi0(k) / dt ) ! For consistency with prognostic macrophysics scheme
3527 nlten_sink(k) = max( nlten_sub, - tr0(k,ixnumliq) / dt )
3528 niten_sink(k) = max( niten_sub, - tr0(k,ixnumice) / dt )
3531 ! --------------------------------------------- !
3533 ! Calculate convective tendencies at each layer !
3535 ! --------------------------------------------- !
3537 ! ----------------- !
3538 ! Momentum tendency !
3539 ! ----------------- !
3543 uten(k) = ( uflx(km1) - uflx(k) ) * g / dp0(k)
3544 vten(k) = ( vflx(km1) - vflx(k) ) * g / dp0(k)
3545 uf(k) = u0(k) + uten(k) * dt
3546 vf(k) = v0(k) + vten(k) * dt
3548 ! trten(k,m) = ( trflx(km1,m) - trflx(k,m) ) * g / dp0(k)
3549 ! ! Limit trten(k,m) such that negative value is not developed.
3550 ! ! This limitation does not conserve grid-mean tracers and future
3551 ! ! refinement is required for tracer-conserving treatment.
3552 ! trten(k,m) = max(trten(k,m),-tr0(k,m)/dt)
3556 ! ----------------------------------------------------------------- !
3557 ! Tendencies of thermodynamic variables. !
3558 ! This part requires a careful treatment of bulk cloud microphysics.!
3559 ! Relocations of 'precipitable condensates' either into the surface !
3560 ! or into the tendency of 'krel' layer will be performed just after !
3561 ! finishing the below 'do-loop'. !
3562 ! ----------------------------------------------------------------- !
3572 ! ------------------------------------------------------------------------------ !
3573 ! Compute 'slten', 'qtten', 'qvten', 'qlten', 'qiten', and 'sten' !
3575 ! Key assumptions made in this 'cumulus scheme' are : !
3576 ! 1. Cumulus updraft expels condensate into the environment at the top interface !
3577 ! of each layer. Note that in addition to this expel process ('source' term), !
3578 ! cumulus updraft can modify layer mean condensate through normal detrainment !
3579 ! forcing or compensating subsidence. !
3580 ! 2. Expelled water can be either 'sustaining' or 'precipitating' condensate. By !
3581 ! definition, 'suataining condensate' will remain in the layer where it was !
3582 ! formed, while 'precipitating condensate' will fall across the base of the !
3583 ! layer where it was formed. !
3584 ! 3. All precipitating condensates are assumed to fall into the release layer or !
3585 ! ground as soon as it was formed without being evaporated during the falling !
3586 ! process down to the desinated layer ( either release layer of surface ). !
3587 ! ------------------------------------------------------------------------------ !
3589 ! ------------------------------------------------------------------------- !
3590 ! 'dwten(k)','diten(k)' : Production rate of condensate within the layer k !
3591 ! [ kg/kg/s ] by the expels of condensate from cumulus updraft. !
3592 ! It is important to note that in terms of moisture tendency equation, this !
3593 ! is a 'source' term of enviromental 'qt'. More importantly, these source !
3594 ! are already counted in the turbulent heat and moisture fluxes we computed !
3595 ! until now, assuming all the expelled condensate remain in the layer where !
3596 ! it was formed. Thus, in calculation of 'qtten' and 'slten' below, we MUST !
3597 ! NOT add or subtract these terms explicitly in order not to double or miss !
3598 ! count, unless some expelled condensates fall down out of the layer. Note !
3599 ! this falling-down process ( i.e., precipitation process ) and associated !
3600 ! 'qtten' and 'slten' and production of surface precipitation flux will be !
3601 ! treated later in 'zm_conv_evap' in 'convect_shallow_tend' subroutine. !
3602 ! In below, we are converting expelled cloud condensate into correct unit. !
3603 ! I found that below use of '0.5 * (umf(k-1) + umf(k))' causes conservation !
3604 ! errors at some columns in global simulation. So, I returned to originals. !
3605 ! This will cause no precipitation flux at 'kpen' layer since umf(kpen)=0. !
3606 ! ------------------------------------------------------------------------- !
3608 dwten(k) = dwten(k) * 0.5_r8 * ( umf(k-1) + umf(k) ) * g / dp0(k) ! [ kg/kg/s ]
3609 diten(k) = diten(k) * 0.5_r8 * ( umf(k-1) + umf(k) ) * g / dp0(k) ! [ kg/kg/s ]
3611 ! dwten(k) = dwten(k) * umf(k) * g / dp0(k) ! [ kg/kg/s ]
3612 ! diten(k) = diten(k) * umf(k) * g / dp0(k) ! [ kg/kg/s ]
3614 ! --------------------------------------------------------------------------- !
3615 ! 'qrten(k)','qsten(k)' : Production rate of rain and snow within the layer k !
3616 ! [ kg/kg/s ] by cumulus expels of condensates to the environment.!
3617 ! This will be falled-out of the layer where it was formed and will be dumped !
3618 ! dumped into the release layer assuming that there is no evaporative cooling !
3619 ! while precipitable condensate moves to the relaes level. This is reasonable !
3620 ! assumtion if cumulus is purely vertical and so the path along which precita !
3621 ! ble condensate falls is fully saturared. This 're-allocation' process of !
3622 ! precipitable condensate into the release layer is fully described in this !
3623 ! convection scheme. After that, the dumped water into the release layer will !
3624 ! falling down across the base of release layer ( or LCL, if exact treatment !
3625 ! is required ) and will be allowed to be evaporated in layers below release !
3626 ! layer, and finally non-zero surface precipitation flux will be calculated. !
3627 ! This latter process will be separately treated 'zm_conv_evap' routine. !
3628 ! --------------------------------------------------------------------------- !
3630 qrten(k) = frc_rasn * dwten(k)
3631 qsten(k) = frc_rasn * diten(k)
3633 ! ----------------------------------------------------------------------- !
3634 ! 'rainflx','snowflx' : Cumulative rain and snow flux integrated from the !
3635 ! [ kg/m2/s ] release leyer to the 'kpen' layer. Note that even !
3636 ! though wtw(kpen) < 0 (and umf(kpen) = 0) at the top interface of 'kpen' !
3637 ! layer, 'dwten(kpen)' and diten(kpen) were calculated after calculating !
3638 ! explicit cloud top height. Thus below calculation of precipitation flux !
3639 ! is correct. Note that precipitating condensates are formed only in the !
3640 ! layers from 'krel' to 'kpen', including the two layers. !
3641 ! ----------------------------------------------------------------------- !
3643 rainflx = rainflx + qrten(k) * dp0(k) / g
3644 snowflx = snowflx + qsten(k) * dp0(k) / g
3646 ! ------------------------------------------------------------------------ !
3647 ! 'slten(k)','qtten(k)' !
3648 ! Note that 'slflx(k)' and 'qtflx(k)' we have calculated already included !
3649 ! all the contributions of (1) expels of condensate (dwten(k), diten(k)), !
3650 ! (2) mass detrainment ( delta * umf * ( qtu - qt ) ), & (3) compensating !
3651 ! subsidence ( M * dqt / dz ). Thus 'slflx(k)' and 'qtflx(k)' we computed !
3652 ! is a hybrid turbulent flux containing one part of 'source' term - expel !
3653 ! of condensate. In order to calculate 'slten' and 'qtten', we should add !
3654 ! additional 'source' term, if any. If the expelled condensate falls down !
3655 ! across the base of the layer, it will be another sink (negative source) !
3656 ! term. Note also that we included frictional heating terms in the below !
3657 ! calculation of 'slten'. !
3658 ! ------------------------------------------------------------------------ !
3660 slten(k) = ( slflx(km1) - slflx(k) ) * g / dp0(k)
3662 slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( &
3663 uflx(k)*(uf(k+1) - uf(k) + u0(k+1) - u0(k)) + &
3664 vflx(k)*(vf(k+1) - vf(k) + v0(k+1) - v0(k)))
3665 elseif( k .ge. 2 .and. k .le. kpen-1 ) then
3666 slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( &
3667 uflx(k)*(uf(k+1) - uf(k) + u0(k+1) - u0(k)) + &
3668 uflx(k-1)*(uf(k) - uf(k-1) + u0(k) - u0(k-1)) + &
3669 vflx(k)*(vf(k+1) - vf(k) + v0(k+1) - v0(k)) + &
3670 vflx(k-1)*(vf(k) - vf(k-1) + v0(k) - v0(k-1)))
3671 elseif( k .eq. kpen ) then
3672 slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( &
3673 uflx(k-1)*(uf(k) - uf(k-1) + u0(k) - u0(k-1)) + &
3674 vflx(k-1)*(vf(k) - vf(k-1) + v0(k) - v0(k-1)))
3676 qtten(k) = ( qtflx(km1) - qtflx(k) ) * g / dp0(k)
3678 ! ---------------------------------------------------------------------------- !
3679 ! Compute condensate tendency, including reserved condensate !
3680 ! We assume that eventual detachment and detrainment occurs in kbup layer due !
3681 ! to downdraft buoyancy sorting. In the layer above the kbup, only penetrative !
3682 ! entrainment exists. Penetrative entrained air is assumed not to contain any !
3684 ! ---------------------------------------------------------------------------- !
3686 ! Compute in-cumulus condensate at the layer mid-point.
3688 if( k .lt. krel .or. k .gt. kpen ) then
3693 elseif( k .eq. krel ) then
3694 call conden(prel,thlu(krel-1),qtu(krel-1),thj,qvj,qlj,qij,qse,id_check,qsat)
3695 if( id_check .eq. 1 ) then
3696 exit_conden(i) = 1._r8
3702 call conden(ps0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check,qsat)
3703 if( id_check .eq. 1 ) then
3704 exit_conden(i) = 1._r8
3708 qlu_mid = 0.5_r8 * ( qlubelow + qlj ) * ( prel - ps0(k) )/( ps0(k-1) - ps0(k) )
3709 qiu_mid = 0.5_r8 * ( qiubelow + qij ) * ( prel - ps0(k) )/( ps0(k-1) - ps0(k) )
3710 elseif( k .eq. kpen ) then
3711 call conden(ps0(k-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check,qsat)
3712 if( id_check .eq. 1 ) then
3713 exit_conden(i) = 1._r8
3717 qlu_mid = 0.5_r8 * ( qlubelow + qlj ) * ( -ppen ) /( ps0(k-1) - ps0(k) )
3718 qiu_mid = 0.5_r8 * ( qiubelow + qij ) * ( -ppen ) /( ps0(k-1) - ps0(k) )
3722 call conden(ps0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check,qsat)
3723 if( id_check .eq. 1 ) then
3724 exit_conden(i) = 1._r8
3728 qlu_mid = 0.5_r8 * ( qlubelow + qlj )
3729 qiu_mid = 0.5_r8 * ( qiubelow + qij )
3734 ! 1. Sustained Precipitation
3736 qc_l(k) = ( 1._r8 - frc_rasn ) * dwten(k) ! [ kg/kg/s ]
3737 qc_i(k) = ( 1._r8 - frc_rasn ) * diten(k) ! [ kg/kg/s ]
3739 ! 2. Detrained Condensate
3741 if( k .le. kbup ) then
3742 qc_l(k) = qc_l(k) + g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qlu_mid ! [ kg/kg/s ]
3743 qc_i(k) = qc_i(k) + g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qiu_mid ! [ kg/kg/s ]
3744 qc_lm = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * ql0(k)
3745 qc_im = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qi0(k)
3746 ! Below 'nc_lm', 'nc_im' should be used only when frc_rasn = 1.
3747 nc_lm = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * tr0(k,ixnumliq)
3748 nc_im = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * tr0(k,ixnumice)
3756 ! 3. Detached Updraft
3758 if( k .eq. kbup ) then
3759 qc_l(k) = qc_l(k) + g * umf(k) * qlj / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3760 qc_i(k) = qc_i(k) + g * umf(k) * qij / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3761 qc_lm = qc_lm - g * umf(k) * ql0(k) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3762 qc_im = qc_im - g * umf(k) * qi0(k) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3763 nc_lm = nc_lm - g * umf(k) * tr0(k,ixnumliq) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3764 nc_im = nc_im - g * umf(k) * tr0(k,ixnumice) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3767 ! 4. Cumulative Penetrative entrainment detrained in the 'kbup' layer
3768 ! Explicitly compute the properties detrained penetrative entrained airs in k = kbup layer.
3770 if( k .eq. kbup ) then
3771 call conden(p0(k),thlu_emf(k),qtu_emf(k),thj,qvj,ql_emf_kbup,qi_emf_kbup,qse,id_check,qsat)
3772 if( id_check .eq. 1 ) then
3776 if( ql_emf_kbup .gt. 0._r8 ) then
3777 nl_emf_kbup = tru_emf(k,ixnumliq)
3781 if( qi_emf_kbup .gt. 0._r8 ) then
3782 ni_emf_kbup = tru_emf(k,ixnumice)
3786 qc_lm = qc_lm - g * emf(k) * ( ql_emf_kbup - ql0(k) ) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3787 qc_im = qc_im - g * emf(k) * ( qi_emf_kbup - qi0(k) ) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3788 nc_lm = nc_lm - g * emf(k) * ( nl_emf_kbup - tr0(k,ixnumliq) ) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3789 nc_im = nc_im - g * emf(k) * ( ni_emf_kbup - tr0(k,ixnumice) ) / ( ps0(k-1) - ps0(k) ) ! [ kg/kg/s ]
3792 qlten_det = qc_l(k) + qc_lm
3793 qiten_det = qc_i(k) + qc_im
3795 ! --------------------------------------------------------------------------------- !
3796 ! 'qlten(k)','qiten(k)','qvten(k)','sten(k)' !
3797 ! Note that falling of precipitation will be treated later. !
3798 ! The prevension of negative 'qv,ql,qi' will be treated later in positive_moisture. !
3799 ! --------------------------------------------------------------------------------- !
3801 if( use_expconten ) then
3802 if( use_unicondet ) then
3805 qlten(k) = frc_rasn * dwten(k) + qlten_sink(k) + qlten_det
3806 qiten(k) = frc_rasn * diten(k) + qiten_sink(k) + qiten_det
3808 qlten(k) = qc_l(k) + frc_rasn * dwten(k) + ( max( 0._r8, ql0(k) + ( qc_lm + qlten_sink(k) ) * dt ) - ql0(k) ) / dt
3809 qiten(k) = qc_i(k) + frc_rasn * diten(k) + ( max( 0._r8, qi0(k) + ( qc_im + qiten_sink(k) ) * dt ) - qi0(k) ) / dt
3810 trten(k,ixnumliq) = max( nc_lm + nlten_sink(k), - tr0(k,ixnumliq) / dt )
3811 trten(k,ixnumice) = max( nc_im + niten_sink(k), - tr0(k,ixnumice) / dt )
3814 if( use_unicondet ) then
3818 qlten(k) = dwten(k) + ( qtten(k) - dwten(k) - diten(k) ) * ( ql0(k) / qt0(k) )
3819 qiten(k) = diten(k) + ( qtten(k) - dwten(k) - diten(k) ) * ( qi0(k) / qt0(k) )
3822 qvten(k) = qtten(k) - qlten(k) - qiten(k)
3823 sten(k) = slten(k) + xlv * qlten(k) + xls * qiten(k)
3825 ! -------------------------------------------------------------------------- !
3826 ! 'rliq' : Verticall-integrated 'suspended cloud condensate' !
3827 ! [m/s] This is so called 'reserved liquid water' in other subroutines !
3828 ! of CAM3, since the contribution of this term should not be included into !
3829 ! the tendency of each layer or surface flux (precip) within this cumulus !
3830 ! scheme. The adding of this term to the layer tendency will be done inthe !
3831 ! 'stratiform_tend', just after performing sediment process there. !
3832 ! The main problem of these rather going-back-and-forth and stupid-seeming !
3833 ! approach is that the sediment process of suspendened condensate will not !
3834 ! be treated at all in the 'stratiform_tend'. !
3835 ! Note that 'precip' [m/s] is vertically-integrated total 'rain+snow' formed !
3836 ! from the cumulus updraft. Important : in the below, 1000 is rhoh2o ( water !
3837 ! density ) [ kg/m^3 ] used for unit conversion from [ kg/m^2/s ] to [ m/s ] !
3838 ! for use in stratiform.F90. !
3839 ! -------------------------------------------------------------------------- !
3841 qc(k) = qc_l(k) + qc_i(k)
3842 rliq = rliq + qc(k) * dp0(k) / g / 1000._r8 ! [ m/s ]
3846 precip = rainflx + snowflx ! [ kg/m2/s ]
3847 snow = snowflx ! [ kg/m2/s ]
3849 ! ---------------------------------------------------------------- !
3850 ! Now treats the 'evaporation' and 'melting' of rain ( qrten ) and !
3851 ! snow ( qsten ) during falling process. Below algorithms are from !
3852 ! 'zm_conv_evap' but with some modification, which allows separate !
3853 ! treatment of 'rain' and 'snow' condensates. Note that I included !
3854 ! the evaporation dynamics into the convection scheme for complete !
3855 ! development of cumulus scheme especially in association with the !
3856 ! implicit CIN closure. In compatible with this internal treatment !
3857 ! of evaporation, I should modify 'convect_shallow', in such that !
3858 ! 'zm_conv_evap' is not performed when I choose UW PBL-Cu schemes. !
3859 ! ---------------------------------------------------------------- !
3863 flxrain(0:mkx) = 0._r8
3864 flxsnow(0:mkx) = 0._r8
3865 ntraprd(:mkx) = 0._r8
3866 ntsnprd(:mkx) = 0._r8
3868 do k = mkx, 1, -1 ! 'k' is a layer index : 'mkx'('1') is the top ('bottom') layer
3870 ! ----------------------------------------------------------------------------- !
3871 ! flxsntm [kg/m2/s] : Downward snow flux at the top of each layer after melting.!
3872 ! snowmlt [kg/kg/s] : Snow melting tendency. !
3873 ! Below allows melting of snow when it goes down into the warm layer below. !
3874 ! ----------------------------------------------------------------------------- !
3876 if( t0(k) .gt. 273.16_r8 ) then
3877 snowmlt = max( 0._r8, flxsnow(k) * g / dp0(k) )
3882 ! ----------------------------------------------------------------- !
3883 ! Evaporation rate of 'rain' and 'snow' in the layer k, [ kg/kg/s ] !
3884 ! where 'rain' and 'snow' are coming down from the upper layers. !
3885 ! I used the same evaporative efficiency both for 'rain' and 'snow'.!
3886 ! Note that evaporation is not allowed in the layers 'k >= krel' by !
3887 ! assuming that inside of cumulus cloud, across which precipitation !
3888 ! is falling down, is fully saturated. !
3889 ! The asumptions in association with the 'evplimit_rain(snow)' are !
3890 ! 1. Do not allow evaporation to supersate the layer !
3891 ! 2. Do not evaporate more than the flux falling into the layer !
3892 ! 3. Total evaporation cannot exceed the input total surface flux !
3893 ! ----------------------------------------------------------------- !
3895 status = qsat(t0(k),p0(k),es(1),qs(1),gam(1), 1)
3896 subsat = max( ( 1._r8 - qv0(k)/qs(1) ), 0._r8 )
3897 if( noevap_krelkpen ) then
3898 if( k .ge. krel ) subsat = 0._r8
3901 evprain = kevp * subsat * sqrt(flxrain(k)+snowmlt*dp0(k)/g)
3902 evpsnow = kevp * subsat * sqrt(max(flxsnow(k)-snowmlt*dp0(k)/g,0._r8))
3904 evplimit = max( 0._r8, ( qw0_in(i,k) - qv0(k) ) / dt )
3906 evplimit_rain = min( evplimit, ( flxrain(k) + snowmlt * dp0(k) / g ) * g / dp0(k) )
3907 evplimit_rain = min( evplimit_rain, ( rainflx - evpint_rain ) * g / dp0(k) )
3908 evprain = max(0._r8,min( evplimit_rain, evprain ))
3910 evplimit_snow = min( evplimit, max( flxsnow(k) - snowmlt * dp0(k) / g , 0._r8 ) * g / dp0(k) )
3911 evplimit_snow = min( evplimit_snow, ( snowflx - evpint_snow ) * g / dp0(k) )
3912 evpsnow = max(0._r8,min( evplimit_snow, evpsnow ))
3914 if( ( evprain + evpsnow ) .gt. evplimit ) then
3915 tmp1 = evprain * evplimit / ( evprain + evpsnow )
3916 tmp2 = evpsnow * evplimit / ( evprain + evpsnow )
3921 evapc(k) = evprain + evpsnow
3923 ! ------------------------------------------------------------- !
3924 ! Vertically-integrated evaporative fluxes of 'rain' and 'snow' !
3925 ! ------------------------------------------------------------- !
3927 evpint_rain = evpint_rain + evprain * dp0(k) / g
3928 evpint_snow = evpint_snow + evpsnow * dp0(k) / g
3930 ! -------------------------------------------------------------- !
3931 ! Net 'rain' and 'snow' production rate in the layer [ kg/kg/s ] !
3932 ! -------------------------------------------------------------- !
3934 ntraprd(k) = qrten(k) - evprain + snowmlt
3935 ntsnprd(k) = qsten(k) - evpsnow - snowmlt
3937 ! -------------------------------------------------------------------------------- !
3938 ! Downward fluxes of 'rain' and 'snow' fluxes at the base of the layer [ kg/m2/s ] !
3939 ! Note that layer index increases with height. !
3940 ! -------------------------------------------------------------------------------- !
3942 flxrain(k-1) = flxrain(k) + ntraprd(k) * dp0(k) / g
3943 flxsnow(k-1) = flxsnow(k) + ntsnprd(k) * dp0(k) / g
3944 flxrain(k-1) = max( flxrain(k-1), 0._r8 )
3945 if( flxrain(k-1) .eq. 0._r8 ) ntraprd(k) = -flxrain(k) * g / dp0(k)
3946 flxsnow(k-1) = max( flxsnow(k-1), 0._r8 )
3947 if( flxsnow(k-1) .eq. 0._r8 ) ntsnprd(k) = -flxsnow(k) * g / dp0(k)
3949 ! ---------------------------------- !
3950 ! Calculate thermodynamic tendencies !
3951 ! --------------------------------------------------------------------------- !
3952 ! Note that equivalently, we can write tendency formula of 'sten' and 'slten' !
3953 ! by 'sten(k) = sten(k) - xlv*evprain - xls*evpsnow - (xls-xlv)*snowmlt' & !
3954 ! 'slten(k) = sten(k) - xlv*qlten(k) - xls*qiten(k)'. !
3955 ! The above formula is equivalent to the below formula. However below formula !
3956 ! is preferred since we have already imposed explicit constraint on 'ntraprd' !
3957 ! and 'ntsnprd' in case that flxrain(k-1) < 0 & flxsnow(k-1) < 0._r8 !
3958 ! Note : In future, I can elborate the limiting of 'qlten','qvten','qiten' !
3959 ! such that that energy and moisture conservation error is completely !
3961 ! Re-storation to the positive condensate will be performed later below !
3962 ! --------------------------------------------------------------------------- !
3964 qlten(k) = qlten(k) - qrten(k)
3965 qiten(k) = qiten(k) - qsten(k)
3966 qvten(k) = qvten(k) + evprain + evpsnow
3967 qtten(k) = qlten(k) + qiten(k) + qvten(k)
3968 if( ( qv0(k) + qvten(k)*dt ) .lt. qmin(1) .or. &
3969 ( ql0(k) + qlten(k)*dt ) .lt. qmin(2) .or. &
3970 ( qi0(k) + qiten(k)*dt ) .lt. qmin(3) ) then
3971 limit_negcon(i) = 1._r8
3973 sten(k) = sten(k) - xlv*evprain - xls*evpsnow - (xls-xlv)*snowmlt
3974 slten(k) = sten(k) - xlv*qlten(k) - xls*qiten(k)
3976 ! slten(k) = slten(k) + xlv * ntraprd(k) + xls * ntsnprd(k)
3977 ! sten(k) = slten(k) + xlv * qlten(k) + xls * qiten(k)
3981 ! ------------------------------------------------------------- !
3982 ! Calculate final surface flux of precipitation, rain, and snow !
3983 ! Convert unit to [m/s] for use in 'check_energy_chng'. !
3984 ! ------------------------------------------------------------- !
3986 precip = ( flxrain(0) + flxsnow(0) ) / 1000._r8
3987 snow = flxsnow(0) / 1000._r8
3989 ! --------------------------------------------------------------------------- !
3990 ! Until now, all the calculations are done completely in this shallow cumulus !
3991 ! scheme. If you want to use this cumulus scheme other than CAM3, then do not !
3992 ! perform below block. However, for compatible use with the other subroutines !
3993 ! in CAM3, I should subtract the effect of 'qc(k)' ('rliq') from the tendency !
3994 ! equation in each layer, since this effect will be separately added later in !
3995 ! in 'stratiform_tend' just after performing sediment process there. In order !
3996 ! to be consistent with 'stratiform_tend', just subtract qc(k) from tendency !
3997 ! equation of each layer, but do not add it to the 'precip'. Apprently, this !
3998 ! will violate energy and moisture conservations. However, when performing !
3999 ! conservation check in 'tphysbc.F90' just after 'convect_shallow_tend', we !
4000 ! will add 'qc(k)' ( rliq ) to the surface flux term just for the purpose of !
4001 ! passing the energy-moisture conservation check. Explicit adding-back of 'qc'!
4002 ! to the individual layer tendency equation will be done in 'stratiform_tend' !
4003 ! after performing sediment process there. Simply speaking, in 'tphysbc' just !
4004 ! after 'convect_shallow_tend', we will dump 'rliq' into surface as a 'rain' !
4005 ! in order to satisfy energy and moisture conservation, and in the following !
4006 ! 'stratiform_tend', we will restore it back to 'qlten(k)' ( 'ice' will go to !
4007 ! 'water' there) from surface precipitation. This is a funny but conceptually !
4008 ! entertaining procedure. One concern I have for this complex process is that !
4009 ! output-writed stratiform precipitation amount will be underestimated due to !
4010 ! arbitrary subtracting of 'rliq' in stratiform_tend, where !
4011 ! ' prec_str = prec_sed + prec_pcw - rliq' and 'rliq' is not real but fake. !
4012 ! However, as shown in 'srfxfer.F90', large scale precipitation amount (PRECL)!
4013 ! that is writed-output is corrected written since in 'srfxfer.F90', PRECL = !
4014 ! 'prec_sed + prec_pcw', without including 'rliq'. So current code is correct.!
4015 ! Note also in 'srfxfer.F90', convective precipitation amount is 'PRECC = !
4016 ! prec_zmc(i) + prec_cmf(i)' which is also correct. !
4017 ! --------------------------------------------------------------------------- !
4020 qtten(k) = qtten(k) - qc(k)
4021 qlten(k) = qlten(k) - qc_l(k)
4022 qiten(k) = qiten(k) - qc_i(k)
4023 slten(k) = slten(k) + ( xlv * qc_l(k) + xls * qc_i(k) )
4024 ! ---------------------------------------------------------------------- !
4025 ! Since all reserved condensates will be treated as liquid water in the !
4026 ! 'check_energy_chng' & 'stratiform_tend' without an explicit conversion !
4027 ! algorithm, I should consider explicitly the energy conversions between !
4028 ! 'ice' and 'liquid' - i.e., I should convert 'ice' to 'liquid' and the !
4029 ! necessary energy for this conversion should be subtracted from 'sten'. !
4030 ! Without this conversion here, energy conservation error come out. Note !
4031 ! that there should be no change of 'qvten(k)'. !
4032 ! ---------------------------------------------------------------------- !
4033 sten(k) = sten(k) - ( xls - xlv ) * qc_i(k)
4036 ! --------------------------------------------------------------- !
4037 ! Prevent the onset-of negative condensate at the next time step !
4038 ! Potentially, this block can be moved just in front of the above !
4040 ! --------------------------------------------------------------- !
4042 ! Modification : I should check whether this 'positive_moisture_single' routine is
4043 ! consistent with the one used in UW PBL and cloud macrophysics schemes.
4044 ! Modification : Below may overestimate resulting 'ql, qi' if we use the new 'qc_l', 'qc_i'
4045 ! in combination with the original computation of qlten, qiten. However,
4046 ! if we use new 'qlten,qiten', there is no problem.
4048 qv0_star(:mkx) = qv0(:mkx) + qvten(:mkx) * dt
4049 ql0_star(:mkx) = ql0(:mkx) + qlten(:mkx) * dt
4050 qi0_star(:mkx) = qi0(:mkx) + qiten(:mkx) * dt
4051 s0_star(:mkx) = s0(:mkx) + sten(:mkx) * dt
4052 call positive_moisture_single( xlv, xls, mkx, dt, qmin(1), qmin(2), qmin(3), dp0, &
4053 qv0_star, ql0_star, qi0_star, s0_star, qvten, qlten, qiten, sten )
4054 qtten(:mkx) = qvten(:mkx) + qlten(:mkx) + qiten(:mkx)
4055 slten(:mkx) = sten(:mkx) - xlv * qlten(:mkx) - xls * qiten(:mkx)
4057 ! --------------------- !
4058 ! Tendencies of tracers !
4059 ! --------------------- !
4063 if( m .ne. ixnumliq .and. m .ne. ixnumice ) then
4067 do mm = 1, ntot_amode
4068 if( m .eq. numptr_amode(mm) ) then
4075 trflx_d(0:mkx) = 0._r8
4076 trflx_u(0:mkx) = 0._r8
4078 if( cnst_get_type_byind(m) .eq. 'wet' ) then
4084 dum = ( tr0(k,m) - trmin ) * pdelx / g / dt + trflx(km1,m) - trflx(k,m) + trflx_d(km1)
4085 trflx_d(k) = min( 0._r8, dum )
4088 if( cnst_get_type_byind(m) .eq. 'wet' ) then
4094 dum = ( tr0(k,m) - trmin ) * pdelx / g / dt + trflx(km1,m) - trflx(k,m) + &
4095 trflx_d(km1) - trflx_d(k) - trflx_u(k)
4096 trflx_u(km1) = max( 0._r8, -dum )
4099 if( cnst_get_type_byind(m) .eq. 'wet' ) then
4105 ! Check : I should re-check whether '_u', '_d' are correctly ordered in
4106 ! the below tendency computation.
4107 trten(k,m) = ( trflx(km1,m) - trflx(k,m) + &
4108 trflx_d(km1) - trflx_d(k) + &
4109 trflx_u(km1) - trflx_u(k) ) * g / pdelx
4116 ! ---------------------------------------------------------------- !
4117 ! Cumpute default diagnostic outputs !
4118 ! Note that since 'qtu(krel-1:kpen-1)' & 'thlu(krel-1:kpen-1)' has !
4119 ! been adjusted after detraining cloud condensate into environment !
4120 ! during cumulus updraft motion, below calculations will exactly !
4121 ! reproduce in-cloud properties as shown in the output analysis. !
4122 ! ---------------------------------------------------------------- !
4124 call conden(prel,thlu(krel-1),qtu(krel-1),thj,qvj,qlj,qij,qse,id_check,qsat)
4125 if( id_check .eq. 1 ) then
4126 exit_conden(i) = 1._r8
4130 qcubelow = qlj + qij
4137 ! --------------------------------------------------------------------- !
4138 ! In the below calculations, I explicitly considered cloud base ( LCL ) !
4139 ! and cloud top height ( ps0(kpen-1) + ppen ) !
4140 ! ----------------------------------------------------------------------!
4141 do k = krel, kpen ! This is a layer index
4142 ! ------------------------------------------------------------------ !
4143 ! Calculate cumulus condensate at the upper interface of each layer. !
4144 ! Note 'ppen < 0' and at 'k=kpen' layer, I used 'thlu_top'&'qtu_top' !
4145 ! which explicitly considered zero or non-zero 'fer(kpen)'. !
4146 ! ------------------------------------------------------------------ !
4147 if( k .eq. kpen ) then
4148 call conden(ps0(k-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check,qsat)
4150 call conden(ps0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check,qsat)
4152 if( id_check .eq. 1 ) then
4153 exit_conden(i) = 1._r8
4157 ! ---------------------------------------------------------------- !
4158 ! Calculate in-cloud mean LWC ( qlu(k) ), IWC ( qiu(k) ), & layer !
4159 ! mean cumulus fraction ( cufrc(k) ), vertically-integrated layer !
4160 ! mean LWP and IWP. Expel some of in-cloud condensate at the upper !
4161 ! interface if it is largr than criqc. Note cumulus cloud fraction !
4162 ! is assumed to be twice of core updraft fractional area. Thus LWP !
4163 ! and IWP will be twice of actual value coming from our scheme. !
4164 ! ---------------------------------------------------------------- !
4165 qcu(k) = 0.5_r8 * ( qcubelow + qlj + qij )
4166 qlu(k) = 0.5_r8 * ( qlubelow + qlj )
4167 qiu(k) = 0.5_r8 * ( qiubelow + qij )
4168 cufrc(k) = ( ufrc(k-1) + ufrc(k) )
4169 if( k .eq. krel ) then
4170 cufrc(k) = ( ufrclcl + ufrc(k) )*( prel - ps0(k) )/( ps0(k-1) - ps0(k) )
4171 else if( k .eq. kpen ) then
4172 cufrc(k) = ( ufrc(k-1) + 0._r8 )*( -ppen ) /( ps0(k-1) - ps0(k) )
4173 if( (qlj + qij) .gt. criqc ) then
4174 qcu(k) = 0.5_r8 * ( qcubelow + criqc )
4175 qlu(k) = 0.5_r8 * ( qlubelow + criqc * qlj / ( qlj + qij ) )
4176 qiu(k) = 0.5_r8 * ( qiubelow + criqc * qij / ( qlj + qij ) )
4179 rcwp = rcwp + ( qlu(k) + qiu(k) ) * ( ps0(k-1) - ps0(k) ) / g * cufrc(k)
4180 rlwp = rlwp + qlu(k) * ( ps0(k-1) - ps0(k) ) / g * cufrc(k)
4181 riwp = riwp + qiu(k) * ( ps0(k-1) - ps0(k) ) / g * cufrc(k)
4182 qcubelow = qlj + qij
4186 ! ------------------------------------ !
4187 ! Cloud top and base interface indices !
4188 ! ------------------------------------ !
4189 cnt = real( kpen, r8 )
4190 cnb = real( krel - 1, r8 )
4192 ! ------------------------------------------------------------------------- !
4193 ! End of formal calculation. Below blocks are for implicit CIN calculations !
4194 ! with re-initialization and save variables at iter_cin = 1._r8 !
4195 ! ------------------------------------------------------------------------- !
4197 ! --------------------------------------------------------------- !
4198 ! Adjust the original input profiles for implicit CIN calculation !
4199 ! --------------------------------------------------------------- !
4201 if( iter .ne. iter_cin ) then
4203 ! ------------------------------------------------------------------- !
4204 ! Save the output from "iter_cin = 1" !
4205 ! These output will be writed-out if "iter_cin = 1" was not performed !
4206 ! for some reasons. !
4207 ! ------------------------------------------------------------------- !
4209 qv0_s(:mkx) = qv0(:mkx) + qvten(:mkx) * dt
4210 ql0_s(:mkx) = ql0(:mkx) + qlten(:mkx) * dt
4211 qi0_s(:mkx) = qi0(:mkx) + qiten(:mkx) * dt
4212 s0_s(:mkx) = s0(:mkx) + sten(:mkx) * dt
4213 u0_s(:mkx) = u0(:mkx) + uten(:mkx) * dt
4214 v0_s(:mkx) = v0(:mkx) + vten(:mkx) * dt
4215 qt0_s(:mkx) = qv0_s(:mkx) + ql0_s(:mkx) + qi0_s(:mkx)
4216 t0_s(:mkx) = t0(:mkx) + sten(:mkx) * dt / cp
4218 tr0_s(:mkx,m) = tr0(:mkx,m) + trten(:mkx,m) * dt
4221 umf_s(0:mkx) = umf(0:mkx)
4222 qvten_s(:mkx) = qvten(:mkx)
4223 qlten_s(:mkx) = qlten(:mkx)
4224 qiten_s(:mkx) = qiten(:mkx)
4225 sten_s(:mkx) = sten(:mkx)
4226 uten_s(:mkx) = uten(:mkx)
4227 vten_s(:mkx) = vten(:mkx)
4228 qrten_s(:mkx) = qrten(:mkx)
4229 qsten_s(:mkx) = qsten(:mkx)
4232 evapc_s(:mkx) = evapc(:mkx)
4234 cufrc_s(:mkx) = cufrc(:mkx)
4235 slflx_s(0:mkx) = slflx(0:mkx)
4236 qtflx_s(0:mkx) = qtflx(0:mkx)
4237 qcu_s(:mkx) = qcu(:mkx)
4238 qlu_s(:mkx) = qlu(:mkx)
4239 qiu_s(:mkx) = qiu(:mkx)
4240 fer_s(:mkx) = fer(:mkx)
4241 fdr_s(:mkx) = fdr(:mkx)
4246 qc_s(:mkx) = qc(:mkx)
4249 qtten_s(:mkx) = qtten(:mkx)
4250 slten_s(:mkx) = slten(:mkx)
4251 ufrc_s(0:mkx) = ufrc(0:mkx)
4253 uflx_s(0:mkx) = uflx(0:mkx)
4254 vflx_s(0:mkx) = vflx(0:mkx)
4256 ufrcinvbase_s = ufrcinvbase
4258 winvbase_s = winvbase
4261 pinv_s = ps0(kinv-1)
4264 ppen_s = ps0(kpen-1) + ppen
4268 emfkbup_s = emf(kbup)
4269 cbmflimit_s = cbmflimit
4271 zinv_s = zs0(kinv-1)
4276 wu_s(0:mkx) = wu(0:mkx)
4277 qtu_s(0:mkx) = qtu(0:mkx)
4278 thlu_s(0:mkx) = thlu(0:mkx)
4279 thvu_s(0:mkx) = thvu(0:mkx)
4280 uu_s(0:mkx) = uu(0:mkx)
4281 vu_s(0:mkx) = vu(0:mkx)
4282 qtu_emf_s(0:mkx) = qtu_emf(0:mkx)
4283 thlu_emf_s(0:mkx) = thlu_emf(0:mkx)
4284 uu_emf_s(0:mkx) = uu_emf(0:mkx)
4285 vu_emf_s(0:mkx) = vu_emf(0:mkx)
4286 uemf_s(0:mkx) = uemf(0:mkx)
4288 dwten_s(:mkx) = dwten(:mkx)
4289 diten_s(:mkx) = diten(:mkx)
4290 flxrain_s(0:mkx) = flxrain(0:mkx)
4291 flxsnow_s(0:mkx) = flxsnow(0:mkx)
4292 ntraprd_s(:mkx) = ntraprd(:mkx)
4293 ntsnprd_s(:mkx) = ntsnprd(:mkx)
4295 excessu_arr_s(:mkx) = excessu_arr(:mkx)
4296 excess0_arr_s(:mkx) = excess0_arr(:mkx)
4297 xc_arr_s(:mkx) = xc_arr(:mkx)
4298 aquad_arr_s(:mkx) = aquad_arr(:mkx)
4299 bquad_arr_s(:mkx) = bquad_arr(:mkx)
4300 cquad_arr_s(:mkx) = cquad_arr(:mkx)
4301 bogbot_arr_s(:mkx) = bogbot_arr(:mkx)
4302 bogtop_arr_s(:mkx) = bogtop_arr(:mkx)
4305 trten_s(:mkx,m) = trten(:mkx,m)
4306 trflx_s(0:mkx,m) = trflx(0:mkx,m)
4307 tru_s(0:mkx,m) = tru(0:mkx,m)
4308 tru_emf_s(0:mkx,m) = tru_emf(0:mkx,m)
4311 ! ----------------------------------------------------------------------------- !
4312 ! Recalculate environmental variables for new cin calculation at "iter_cin = 2" !
4313 ! using the updated state variables. Perform only for variables necessary for !
4314 ! the new cin calculation. !
4315 ! ----------------------------------------------------------------------------- !
4317 qv0(:mkx) = qv0_s(:mkx)
4318 ql0(:mkx) = ql0_s(:mkx)
4319 qi0(:mkx) = qi0_s(:mkx)
4320 s0(:mkx) = s0_s(:mkx)
4321 t0(:mkx) = t0_s(:mkx)
4323 qt0(:mkx) = (qv0(:mkx) + ql0(:mkx) + qi0(:mkx))
4324 thl0(:mkx) = (t0(:mkx) - xlv*ql0(:mkx)/cp - xls*qi0(:mkx)/cp)/exn0(:mkx)
4325 thvl0(:mkx) = (1._r8 + zvir*qt0(:mkx))*thl0(:mkx)
4327 ssthl0 = slope(mkx,thl0,p0) ! Dimension of ssthl0(:mkx) is implicit
4328 ssqt0 = slope(mkx,qt0 ,p0)
4329 ssu0 = slope(mkx,u0 ,p0)
4330 ssv0 = slope(mkx,v0 ,p0)
4332 sstr0(:mkx,m) = slope(mkx,tr0(:mkx,m),p0)
4337 thl0bot = thl0(k) + ssthl0(k) * ( ps0(k-1) - p0(k) )
4338 qt0bot = qt0(k) + ssqt0(k) * ( ps0(k-1) - p0(k) )
4339 call conden(ps0(k-1),thl0bot,qt0bot,thj,qvj,qlj,qij,qse,id_check,qsat)
4340 if( id_check .eq. 1 ) then
4341 exit_conden(i) = 1._r8
4345 thv0bot(k) = thj * ( 1._r8 + zvir*qvj - qlj - qij )
4346 thvl0bot(k) = thl0bot * ( 1._r8 + zvir*qt0bot )
4348 thl0top = thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) )
4349 qt0top = qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) )
4350 call conden(ps0(k),thl0top,qt0top,thj,qvj,qlj,qij,qse,id_check,qsat)
4351 if( id_check .eq. 1 ) then
4352 exit_conden(i) = 1._r8
4356 thv0top(k) = thj * ( 1._r8 + zvir*qvj - qlj - qij )
4357 thvl0top(k) = thl0top * ( 1._r8 + zvir*qt0top )
4361 endif ! End of 'if(iter .ne. iter_cin)' if sentence.
4363 end do ! End of implicit CIN loop (cin_iter)
4365 ! ----------------------- !
4366 ! Update Output Variables !
4367 ! ----------------------- !
4369 umf_out(i,0:mkx) = umf(0:mkx)
4370 slflx_out(i,0:mkx) = slflx(0:mkx)
4371 qtflx_out(i,0:mkx) = qtflx(0:mkx)
4372 !the indices are not reversed, these variables go into compute_mcshallow_inv, this is why they are called "flxprc1" and "flxsnow1".
4373 flxprc1_out(i,0:mkx) = flxrain(0:mkx) + flxsnow(0:mkx)
4374 flxsnow1_out(i,0:mkx) = flxsnow(0:mkx)
4375 qvten_out(i,:mkx) = qvten(:mkx)
4376 qlten_out(i,:mkx) = qlten(:mkx)
4377 qiten_out(i,:mkx) = qiten(:mkx)
4378 sten_out(i,:mkx) = sten(:mkx)
4379 uten_out(i,:mkx) = uten(:mkx)
4380 vten_out(i,:mkx) = vten(:mkx)
4381 qrten_out(i,:mkx) = qrten(:mkx)
4382 qsten_out(i,:mkx) = qsten(:mkx)
4383 precip_out(i) = precip
4385 evapc_out(i,:mkx) = evapc(:mkx)
4386 cufrc_out(i,:mkx) = cufrc(:mkx)
4387 qcu_out(i,:mkx) = qcu(:mkx)
4388 qlu_out(i,:mkx) = qlu(:mkx)
4389 qiu_out(i,:mkx) = qiu(:mkx)
4390 cush_inout(i) = cush
4393 qc_out(i,:mkx) = qc(:mkx)
4398 trten_out(i,:mkx,m) = trten(:mkx,m)
4401 ! ------------------------------------------------- !
4402 ! Below are specific diagnostic output for detailed !
4403 ! analysis of cumulus scheme !
4404 ! ------------------------------------------------- !
4406 fer_out(i,mkx:1:-1) = fer(:mkx)
4407 fdr_out(i,mkx:1:-1) = fdr(:mkx)
4409 cinlclh_out(i) = cinlcl
4410 qtten_out(i,mkx:1:-1) = qtten(:mkx)
4411 slten_out(i,mkx:1:-1) = slten(:mkx)
4412 ufrc_out(i,mkx:0:-1) = ufrc(0:mkx)
4413 uflx_out(i,mkx:0:-1) = uflx(0:mkx)
4414 vflx_out(i,mkx:0:-1) = vflx(0:mkx)
4416 ufrcinvbase_out(i) = ufrcinvbase
4417 ufrclcl_out(i) = ufrclcl
4418 winvbase_out(i) = winvbase
4421 pinv_out(i) = ps0(kinv-1)
4423 pbup_out(i) = ps0(kbup)
4424 ppen_out(i) = ps0(kpen-1) + ppen
4425 qtsrc_out(i) = qtsrc
4426 thlsrc_out(i) = thlsrc
4427 thvlsrc_out(i) = thvlsrc
4428 emfkbup_out(i) = emf(kbup)
4429 cbmflimit_out(i) = cbmflimit
4430 tkeavg_out(i) = tkeavg
4431 zinv_out(i) = zs0(kinv-1)
4436 wu_out(i,mkx:0:-1) = wu(0:mkx)
4437 qtu_out(i,mkx:0:-1) = qtu(0:mkx)
4438 thlu_out(i,mkx:0:-1) = thlu(0:mkx)
4439 thvu_out(i,mkx:0:-1) = thvu(0:mkx)
4440 uu_out(i,mkx:0:-1) = uu(0:mkx)
4441 vu_out(i,mkx:0:-1) = vu(0:mkx)
4442 qtu_emf_out(i,mkx:0:-1) = qtu_emf(0:mkx)
4443 thlu_emf_out(i,mkx:0:-1) = thlu_emf(0:mkx)
4444 uu_emf_out(i,mkx:0:-1) = uu_emf(0:mkx)
4445 vu_emf_out(i,mkx:0:-1) = vu_emf(0:mkx)
4446 uemf_out(i,mkx:0:-1) = uemf(0:mkx)
4448 dwten_out(i,mkx:1:-1) = dwten(:mkx)
4449 diten_out(i,mkx:1:-1) = diten(:mkx)
4450 flxrain_out(i,mkx:0:-1) = flxrain(0:mkx)
4451 flxsnow_out(i,mkx:0:-1) = flxsnow(0:mkx)
4452 ntraprd_out(i,mkx:1:-1) = ntraprd(:mkx)
4453 ntsnprd_out(i,mkx:1:-1) = ntsnprd(:mkx)
4455 excessu_arr_out(i,mkx:1:-1) = excessu_arr(:mkx)
4456 excess0_arr_out(i,mkx:1:-1) = excess0_arr(:mkx)
4457 xc_arr_out(i,mkx:1:-1) = xc_arr(:mkx)
4458 aquad_arr_out(i,mkx:1:-1) = aquad_arr(:mkx)
4459 bquad_arr_out(i,mkx:1:-1) = bquad_arr(:mkx)
4460 cquad_arr_out(i,mkx:1:-1) = cquad_arr(:mkx)
4461 bogbot_arr_out(i,mkx:1:-1) = bogbot_arr(:mkx)
4462 bogtop_arr_out(i,mkx:1:-1) = bogtop_arr(:mkx)
4465 trflx_out(i,mkx:0:-1,m) = trflx(0:mkx,m)
4466 tru_out(i,mkx:0:-1,m) = tru(0:mkx,m)
4467 tru_emf_out(i,mkx:0:-1,m) = tru_emf(0:mkx,m)
4470 333 if(id_exit) then ! Exit without cumulus convection
4472 exit_UWCu(i) = 1._r8
4474 ! --------------------------------------------------------------------- !
4475 ! Initialize output variables when cumulus convection was not performed.!
4476 ! --------------------------------------------------------------------- !
4478 umf_out(i,0:mkx) = 0._r8
4479 slflx_out(i,0:mkx) = 0._r8
4480 qtflx_out(i,0:mkx) = 0._r8
4481 qvten_out(i,:mkx) = 0._r8
4482 qlten_out(i,:mkx) = 0._r8
4483 qiten_out(i,:mkx) = 0._r8
4484 sten_out(i,:mkx) = 0._r8
4485 uten_out(i,:mkx) = 0._r8
4486 vten_out(i,:mkx) = 0._r8
4487 qrten_out(i,:mkx) = 0._r8
4488 qsten_out(i,:mkx) = 0._r8
4489 precip_out(i) = 0._r8
4491 evapc_out(i,:mkx) = 0._r8
4492 cufrc_out(i,:mkx) = 0._r8
4493 qcu_out(i,:mkx) = 0._r8
4494 qlu_out(i,:mkx) = 0._r8
4495 qiu_out(i,:mkx) = 0._r8
4496 cush_inout(i) = -1._r8
4499 qc_out(i,:mkx) = 0._r8
4501 cnb_out(i) = real(mkx, r8)
4503 fer_out(i,mkx:1:-1) = 0._r8
4504 fdr_out(i,mkx:1:-1) = 0._r8
4505 cinh_out(i) = -1._r8
4506 cinlclh_out(i) = -1._r8
4507 qtten_out(i,mkx:1:-1) = 0._r8
4508 slten_out(i,mkx:1:-1) = 0._r8
4509 ufrc_out(i,mkx:0:-1) = 0._r8
4510 uflx_out(i,mkx:0:-1) = 0._r8
4511 vflx_out(i,mkx:0:-1) = 0._r8
4513 ufrcinvbase_out(i) = 0._r8
4514 ufrclcl_out(i) = 0._r8
4515 winvbase_out(i) = 0._r8
4522 qtsrc_out(i) = 0._r8
4523 thlsrc_out(i) = 0._r8
4524 thvlsrc_out(i) = 0._r8
4525 emfkbup_out(i) = 0._r8
4526 cbmflimit_out(i) = 0._r8
4527 tkeavg_out(i) = 0._r8
4533 wu_out(i,mkx:0:-1) = 0._r8
4534 qtu_out(i,mkx:0:-1) = 0._r8
4535 thlu_out(i,mkx:0:-1) = 0._r8
4536 thvu_out(i,mkx:0:-1) = 0._r8
4537 uu_out(i,mkx:0:-1) = 0._r8
4538 vu_out(i,mkx:0:-1) = 0._r8
4539 qtu_emf_out(i,mkx:0:-1) = 0._r8
4540 thlu_emf_out(i,mkx:0:-1) = 0._r8
4541 uu_emf_out(i,mkx:0:-1) = 0._r8
4542 vu_emf_out(i,mkx:0:-1) = 0._r8
4543 uemf_out(i,mkx:0:-1) = 0._r8
4545 dwten_out(i,mkx:1:-1) = 0._r8
4546 diten_out(i,mkx:1:-1) = 0._r8
4547 flxrain_out(i,mkx:0:-1) = 0._r8
4548 flxsnow_out(i,mkx:0:-1) = 0._r8
4549 ntraprd_out(i,mkx:1:-1) = 0._r8
4550 ntsnprd_out(i,mkx:1:-1) = 0._r8
4552 excessu_arr_out(i,mkx:1:-1) = 0._r8
4553 excess0_arr_out(i,mkx:1:-1) = 0._r8
4554 xc_arr_out(i,mkx:1:-1) = 0._r8
4555 aquad_arr_out(i,mkx:1:-1) = 0._r8
4556 bquad_arr_out(i,mkx:1:-1) = 0._r8
4557 cquad_arr_out(i,mkx:1:-1) = 0._r8
4558 bogbot_arr_out(i,mkx:1:-1) = 0._r8
4559 bogtop_arr_out(i,mkx:1:-1) = 0._r8
4562 trten_out(i,:mkx,m) = 0._r8
4563 trflx_out(i,mkx:0:-1,m) = 0._r8
4564 tru_out(i,mkx:0:-1,m) = 0._r8
4565 tru_emf_out(i,mkx:0:-1,m) = 0._r8
4570 end do ! end of big i loop for each column.
4572 ! ---------------------------------------- !
4573 ! Writing main diagnostic output variables !
4574 ! ---------------------------------------- !
4576 call outfld( 'qtflx_Cu' , qtflx_out(:,mkx:0:-1), mix, lchnk )
4577 call outfld( 'slflx_Cu' , slflx_out(:,mkx:0:-1), mix, lchnk )
4578 call outfld( 'uflx_Cu' , uflx_out, mix, lchnk )
4579 call outfld( 'vflx_Cu' , vflx_out, mix, lchnk )
4581 call outfld( 'qtten_Cu' , qtten_out, mix, lchnk )
4582 call outfld( 'slten_Cu' , slten_out, mix, lchnk )
4583 call outfld( 'uten_Cu' , uten_out(:,mkx:1:-1), mix, lchnk )
4584 call outfld( 'vten_Cu' , vten_out(:,mkx:1:-1), mix, lchnk )
4585 call outfld( 'qvten_Cu' , qvten_out(:,mkx:1:-1), mix, lchnk )
4586 call outfld( 'qlten_Cu' , qlten_out(:,mkx:1:-1), mix, lchnk )
4587 call outfld( 'qiten_Cu' , qiten_out(:,mkx:1:-1), mix, lchnk )
4589 call outfld( 'cbmf_Cu' , cbmf_out, mix, lchnk )
4590 call outfld( 'ufrcinvbase_Cu' , ufrcinvbase_out, mix, lchnk )
4591 call outfld( 'ufrclcl_Cu' , ufrclcl_out, mix, lchnk )
4592 call outfld( 'winvbase_Cu' , winvbase_out, mix, lchnk )
4593 call outfld( 'wlcl_Cu' , wlcl_out, mix, lchnk )
4594 call outfld( 'plcl_Cu' , plcl_out, mix, lchnk )
4595 call outfld( 'pinv_Cu' , pinv_out, mix, lchnk )
4596 call outfld( 'plfc_Cu' , plfc_out, mix, lchnk )
4597 call outfld( 'pbup_Cu' , pbup_out, mix, lchnk )
4598 call outfld( 'ppen_Cu' , ppen_out, mix, lchnk )
4599 call outfld( 'qtsrc_Cu' , qtsrc_out, mix, lchnk )
4600 call outfld( 'thlsrc_Cu' , thlsrc_out, mix, lchnk )
4601 call outfld( 'thvlsrc_Cu' , thvlsrc_out, mix, lchnk )
4602 call outfld( 'emfkbup_Cu' , emfkbup_out, mix, lchnk )
4603 call outfld( 'cin_Cu' , cinh_out, mix, lchnk )
4604 call outfld( 'cinlcl_Cu' , cinlclh_out, mix, lchnk )
4605 call outfld( 'cbmflimit_Cu' , cbmflimit_out, mix, lchnk )
4606 call outfld( 'tkeavg_Cu' , tkeavg_out, mix, lchnk )
4607 call outfld( 'zinv_Cu' , zinv_out, mix, lchnk )
4608 call outfld( 'rcwp_Cu' , rcwp_out, mix, lchnk )
4609 call outfld( 'rlwp_Cu' , rlwp_out, mix, lchnk )
4610 call outfld( 'riwp_Cu' , riwp_out, mix, lchnk )
4611 call outfld( 'tophgt_Cu' , cush_inout, mix, lchnk )
4613 call outfld( 'wu_Cu' , wu_out, mix, lchnk )
4614 call outfld( 'ufrc_Cu' , ufrc_out, mix, lchnk )
4615 call outfld( 'qtu_Cu' , qtu_out, mix, lchnk )
4616 call outfld( 'thlu_Cu' , thlu_out, mix, lchnk )
4617 call outfld( 'thvu_Cu' , thvu_out, mix, lchnk )
4618 call outfld( 'uu_Cu' , uu_out, mix, lchnk )
4619 call outfld( 'vu_Cu' , vu_out, mix, lchnk )
4620 call outfld( 'qtu_emf_Cu' , qtu_emf_out, mix, lchnk )
4621 call outfld( 'thlu_emf_Cu' , thlu_emf_out, mix, lchnk )
4622 call outfld( 'uu_emf_Cu' , uu_emf_out, mix, lchnk )
4623 call outfld( 'vu_emf_Cu' , vu_emf_out, mix, lchnk )
4624 call outfld( 'umf_Cu' , umf_out(:,mkx:0:-1), mix, lchnk )
4625 call outfld( 'uemf_Cu' , uemf_out, mix, lchnk )
4626 call outfld( 'qcu_Cu' , qcu_out(:,mkx:1:-1), mix, lchnk )
4627 call outfld( 'qlu_Cu' , qlu_out(:,mkx:1:-1), mix, lchnk )
4628 call outfld( 'qiu_Cu' , qiu_out(:,mkx:1:-1), mix, lchnk )
4629 call outfld( 'cufrc_Cu' , cufrc_out(:,mkx:1:-1), mix, lchnk )
4630 call outfld( 'fer_Cu' , fer_out, mix, lchnk )
4631 call outfld( 'fdr_Cu' , fdr_out, mix, lchnk )
4633 call outfld( 'dwten_Cu' , dwten_out, mix, lchnk )
4634 call outfld( 'diten_Cu' , diten_out, mix, lchnk )
4635 call outfld( 'qrten_Cu' , qrten_out(:,mkx:1:-1), mix, lchnk )
4636 call outfld( 'qsten_Cu' , qsten_out(:,mkx:1:-1), mix, lchnk )
4637 call outfld( 'flxrain_Cu' , flxrain_out, mix, lchnk )
4638 call outfld( 'flxsnow_Cu' , flxsnow_out, mix, lchnk )
4639 call outfld( 'ntraprd_Cu' , ntraprd_out, mix, lchnk )
4640 call outfld( 'ntsnprd_Cu' , ntsnprd_out, mix, lchnk )
4642 call outfld( 'excessu_Cu' , excessu_arr_out, mix, lchnk )
4643 call outfld( 'excess0_Cu' , excess0_arr_out, mix, lchnk )
4644 call outfld( 'xc_Cu' , xc_arr_out, mix, lchnk )
4645 call outfld( 'aquad_Cu' , aquad_arr_out, mix, lchnk )
4646 call outfld( 'bquad_Cu' , bquad_arr_out, mix, lchnk )
4647 call outfld( 'cquad_Cu' , cquad_arr_out, mix, lchnk )
4648 call outfld( 'bogbot_Cu' , bogbot_arr_out, mix, lchnk )
4649 call outfld( 'bogtop_Cu' , bogtop_arr_out, mix, lchnk )
4651 call outfld( 'exit_UWCu_Cu' , exit_UWCu, mix, lchnk )
4652 call outfld( 'exit_conden_Cu' , exit_conden, mix, lchnk )
4653 call outfld( 'exit_klclmkx_Cu' , exit_klclmkx, mix, lchnk )
4654 call outfld( 'exit_klfcmkx_Cu' , exit_klfcmkx, mix, lchnk )
4655 call outfld( 'exit_ufrc_Cu' , exit_ufrc, mix, lchnk )
4656 call outfld( 'exit_wtw_Cu' , exit_wtw, mix, lchnk )
4657 call outfld( 'exit_drycore_Cu' , exit_drycore, mix, lchnk )
4658 call outfld( 'exit_wu_Cu' , exit_wu, mix, lchnk )
4659 call outfld( 'exit_cufilter_Cu', exit_cufilter, mix, lchnk )
4660 call outfld( 'exit_kinv1_Cu' , exit_kinv1, mix, lchnk )
4661 call outfld( 'exit_rei_Cu' , exit_rei, mix, lchnk )
4663 call outfld( 'limit_shcu_Cu' , limit_shcu, mix, lchnk )
4664 call outfld( 'limit_negcon_Cu' , limit_negcon, mix, lchnk )
4665 call outfld( 'limit_ufrc_Cu' , limit_ufrc, mix, lchnk )
4666 call outfld( 'limit_ppen_Cu' , limit_ppen, mix, lchnk )
4667 call outfld( 'limit_emf_Cu' , limit_emf, mix, lchnk )
4668 call outfld( 'limit_cinlcl_Cu' , limit_cinlcl, mix, lchnk )
4669 call outfld( 'limit_cin_Cu' , limit_cin, mix, lchnk )
4670 call outfld( 'limit_cbmf_Cu' , limit_cbmf, mix, lchnk )
4671 call outfld( 'limit_rei_Cu' , limit_rei, mix, lchnk )
4672 call outfld( 'ind_delcin_Cu' , ind_delcin, mix, lchnk )
4676 end subroutine compute_uwshcu
4678 ! ------------------------------ !
4680 ! Beginning of subroutine blocks !
4682 ! ------------------------------ !
4684 subroutine getbuoy(pbot,thv0bot,ptop,thv0top,thvubot,thvutop,plfc,cin)
4685 ! ----------------------------------------------------------- !
4686 ! Subroutine to calculate integrated CIN [ J/kg = m2/s2 ] and !
4687 ! 'cinlcl, plfc' if any. Assume 'thv' is linear in each layer !
4688 ! both for cumulus and environment. Note that this subroutine !
4689 ! only include positive CIN in calculation - if there are any !
4690 ! negative CIN, it is assumed to be zero. This is slightly !
4691 ! different from 'single_cin' below, where both positive and !
4692 ! negative CIN are included. !
4693 ! ----------------------------------------------------------- !
4694 real(r8) pbot,thv0bot,ptop,thv0top,thvubot,thvutop,plfc,cin,frc
4696 if( thvubot .gt. thv0bot .and. thvutop .gt. thv0top ) then
4699 elseif( thvubot .le. thv0bot .and. thvutop .le. thv0top ) then
4700 cin = cin - ( (thvubot/thv0bot - 1._r8) + (thvutop/thv0top - 1._r8)) * (pbot - ptop) / &
4701 ( pbot/(r*thv0bot*exnf(pbot)) + ptop/(r*thv0top*exnf(ptop)) )
4702 elseif( thvubot .gt. thv0bot .and. thvutop .le. thv0top ) then
4703 frc = ( thvutop/thv0top - 1._r8 ) / ( (thvutop/thv0top - 1._r8) - (thvubot/thv0bot - 1._r8) )
4704 cin = cin - ( thvutop/thv0top - 1._r8 ) * ( (ptop + frc*(pbot - ptop)) - ptop ) / &
4705 ( pbot/(r*thv0bot*exnf(pbot)) + ptop/(r*thv0top*exnf(ptop)) )
4707 frc = ( thvubot/thv0bot - 1._r8 ) / ( (thvubot/thv0bot - 1._r8) - (thvutop/thv0top - 1._r8) )
4708 plfc = pbot - frc * ( pbot - ptop )
4709 cin = cin - ( thvubot/thv0bot - 1._r8)*(pbot - plfc)/ &
4710 ( pbot/(r*thv0bot*exnf(pbot)) + ptop/(r*thv0top * exnf(ptop)))
4714 end subroutine getbuoy
4716 function single_cin(pbot,thv0bot,ptop,thv0top,thvubot,thvutop)
4717 ! ------------------------------------------------------- !
4718 ! Function to calculate a single layer CIN by summing all !
4719 ! positive and negative CIN. !
4720 ! ------------------------------------------------------- !
4721 real(r8) :: single_cin
4722 real(r8) pbot,thv0bot,ptop,thv0top,thvubot,thvutop
4724 single_cin = ( (1._r8 - thvubot/thv0bot) + (1._r8 - thvutop/thv0top)) * ( pbot - ptop ) / &
4725 ( pbot/(r*thv0bot*exnf(pbot)) + ptop/(r*thv0top*exnf(ptop)) )
4727 end function single_cin
4730 subroutine conden(p,thl,qt,th,qv,ql,qi,rvls,id_check,qsat)
4731 ! --------------------------------------------------------------------- !
4732 ! Calculate thermodynamic properties from a given set of ( p, thl, qt ) !
4733 ! --------------------------------------------------------------------- !
4735 real(r8), intent(in) :: p
4736 real(r8), intent(in) :: thl
4737 real(r8), intent(in) :: qt
4738 real(r8), intent(out) :: th
4739 real(r8), intent(out) :: qv
4740 real(r8), intent(out) :: ql
4741 real(r8), intent(out) :: qi
4742 real(r8), intent(out) :: rvls
4743 integer , intent(out) :: id_check
4744 integer , external :: qsat
4745 real(r8) :: tc,temps,t
4746 real(r8) :: leff, nu, qc
4747 integer :: iteration
4748 real(r8) :: es(1) ! Saturation vapor pressure
4749 real(r8) :: qs(1) ! Saturation spec. humidity
4750 real(r8) :: gam(1) ! (L/cp)*dqs/dT
4751 integer :: status ! Return status of qsat call
4754 ! Modification : In order to be compatible with the dlf treatment in stratiform.F90,
4755 ! we may use ( 268.15, 238.15 ) with 30K ramping instead of 20 K,
4756 ! in computing ice fraction below.
4757 ! Note that 'cldwat_fice' uses ( 243.15, 263.15 ) with 20K ramping for stratus.
4758 nu = max(min((268._r8 - tc)/20._r8,1.0_r8),0.0_r8) ! Fraction of ice in the condensate.
4759 leff = (1._r8 - nu)*xlv + nu*xls ! This is an estimate that hopefully speeds convergence
4761 ! --------------------------------------------------------------------------- !
4762 ! Below "temps" and "rvls" are just initial guesses for iteration loop below. !
4763 ! Note that the output "temps" from the below iteration loop is "temperature" !
4764 ! NOT "liquid temperature". !
4765 ! --------------------------------------------------------------------------- !
4768 status = qsat(temps,p,es(1),qs(1),gam(1), 1)
4771 if( qs(1) .ge. qt ) then
4779 do iteration = 1, 10
4780 temps = temps + ( (tc-temps)*cp/leff + qt - rvls )/( cp/leff + ep2*leff*rvls/r/temps/temps )
4781 status = qsat(temps,p,es(1),qs(1),gam(1),1)
4784 qc = max(qt - qs(1),0._r8)
4786 ql = qc*(1._r8 - nu)
4789 if( abs((temps-(leff/cp)*qc)-tc) .ge. 1._r8 ) then
4797 end subroutine conden
4799 subroutine roots(a,b,c,r1,r2,status)
4800 ! --------------------------------------------------------- !
4801 ! Subroutine to solve the second order polynomial equation. !
4802 ! I should check this subroutine later. !
4803 ! --------------------------------------------------------- !
4804 real(r8), intent(in) :: a
4805 real(r8), intent(in) :: b
4806 real(r8), intent(in) :: c
4807 real(r8), intent(out) :: r1
4808 real(r8), intent(out) :: r2
4809 integer , intent(out) :: status
4814 if( a .eq. 0._r8 ) then ! Form b*x + c = 0
4815 if( b .eq. 0._r8 ) then ! Failure: c = 0
4822 if( b .eq. 0._r8 ) then ! Form a*x**2 + c = 0
4823 if( a*c .gt. 0._r8 ) then ! Failure: x**2 = -c/a < 0
4829 else ! Form a*x**2 + b*x + c = 0
4830 if( (b**2 - 4._r8*a*c) .lt. 0._r8 ) then ! Failure, no real roots
4833 q = -0.5_r8*(b + sign(1.0_r8,b)*sqrt(b**2 - 4._r8*a*c))
4841 end subroutine roots
4843 function slope(mkx,field,p0)
4844 ! ------------------------------------------------------------------ !
4845 ! Function performing profile reconstruction of conservative scalars !
4846 ! in each layer. This is identical to profile reconstruction used in !
4847 ! UW-PBL scheme but from bottom to top layer here. At the lowest !
4848 ! layer near to surface, slope is defined using the two lowest layer !
4849 ! mid-point values. I checked this subroutine and it is correct. !
4850 ! ------------------------------------------------------------------ !
4851 real(r8) :: slope(mkx)
4852 integer, intent(in) :: mkx
4853 real(r8), intent(in) :: field(mkx)
4854 real(r8), intent(in) :: p0(mkx)
4860 below = ( field(2) - field(1) ) / ( p0(2) - p0(1) )
4862 above = ( field(k) - field(k-1) ) / ( p0(k) - p0(k-1) )
4863 if( above .gt. 0._r8 ) then
4864 slope(k-1) = max(0._r8,min(above,below))
4866 slope(k-1) = min(0._r8,max(above,below))
4870 slope(mkx) = slope(mkx-1)
4875 function qsinvert(qt,thl,psfc,qsat)
4876 ! ----------------------------------------------------------------- !
4877 ! Function calculating saturation pressure ps (or pLCL) from qt and !
4878 ! thl ( liquid potential temperature, NOT liquid virtual potential !
4879 ! temperature) by inverting Bolton formula. I should check later if !
4880 ! current use of 'leff' instead of 'xlv' here is reasonable or not. !
4881 ! ----------------------------------------------------------------- !
4882 real(r8) :: qsinvert
4883 real(r8) qt, thl, psfc
4884 real(r8) ps, Pis, Ts, err, dlnqsdT, dTdPis
4885 real(r8) dPisdps, dlnqsdps, derrdps, dps
4886 real(r8) Ti, rhi, TLCL, PiLCL, psmin, dpsmax
4888 integer, external :: qsat
4889 real(r8) :: es(1) ! saturation vapor pressure
4890 real(r8) :: qs(1) ! saturation spec. humidity
4891 real(r8) :: gam(1) ! (L/cp)*dqs/dT
4892 integer :: status ! return status of qsat call
4893 real(r8) :: leff, nu
4895 psmin = 100._r8*100._r8 ! Default saturation pressure [Pa] if iteration does not converge
4896 dpsmax = 1._r8 ! Tolerance [Pa] for convergence of iteration
4898 ! ------------------------------------ !
4899 ! Calculate best initial guess of pLCL !
4900 ! ------------------------------------ !
4902 Ti = thl*(psfc/p00)**rovcp
4903 status = qsat(Ti,psfc,es(1),qs(1),gam(1),1)
4905 if( rhi .le. 0.01_r8 ) then
4906 write(iulog,*) 'Source air is too dry and pLCL is set to psmin in uwshcu.F90'
4908 call wrf_message(iulog)
4913 TLCL = 55._r8 + 1._r8/(1._r8/(Ti-55._r8)-log(rhi)/2840._r8); ! Bolton's formula. MWR.1980.Eq.(22)
4915 ps = p00*(PiLCL)**(1._r8/rovcp)
4918 Pis = (ps/p00)**rovcp
4920 status = qsat(Ts,ps,es(1),qs(1),gam(1),1)
4922 nu = max(min((268._r8 - Ts)/20._r8,1.0_r8),0.0_r8)
4923 leff = (1._r8 - nu)*xlv + nu*xls
4924 dlnqsdT = gam(1)*(cp/leff)/qs(1)
4926 dPisdps = rovcp*Pis/ps
4927 dlnqsdps = -1._r8/(ps - (1._r8 - ep2)*es(1))
4928 derrdps = -qs(1)*(dlnqsdT * dTdPis * dPisdps + dlnqsdps)
4931 if( ps .lt. 0._r8 ) then
4932 write(iulog,*) 'pLCL iteration is negative and set to psmin in uwshcu.F90', qt, thl, psfc
4934 call wrf_message(iulog)
4939 if( abs(dps) .le. dpsmax ) then
4944 write(iulog,*) 'pLCL does not converge and is set to psmin in uwshcu.F90', qt, thl, psfc
4946 call wrf_message(iulog)
4950 end function qsinvert
4952 real(r8) function compute_alpha(del_CIN,ke)
4953 ! ------------------------------------------------ !
4954 ! Subroutine to compute proportionality factor for !
4955 ! implicit CIN calculation. !
4956 ! ------------------------------------------------ !
4957 real(r8) :: del_CIN, ke
4960 integer :: iteration
4963 do iteration = 1, 10
4964 x1 = x0 - (exp(-x0*ke*del_CIN) - x0)/(-ke*del_CIN*exp(-x0*ke*del_CIN) - 1._r8)
4971 end function compute_alpha
4973 real(r8) function compute_mumin2(mulcl,rmaxfrac,mulow)
4974 ! --------------------------------------------------------- !
4975 ! Subroutine to compute critical 'mu' (normalized CIN) such !
4976 ! that updraft fraction at the LCL is equal to 'rmaxfrac'. !
4977 ! --------------------------------------------------------- !
4978 real(r8) :: mulcl, rmaxfrac, mulow
4979 real(r8) :: x0, x1, ex, ef, exf, f, fs
4980 integer :: iteration
4983 do iteration = 1, 10
4986 ! if(x0.ge.3._r8) then
4987 ! compute_mumin2 = 3._r8
4991 f = 0.5_r8*exf**2 - 0.5_r8*(ex/2._r8/rmaxfrac)**2 - (mulcl*2.5066_r8/2._r8)**2
4992 fs = (2._r8*exf**2)*(exf/sqrt(3.141592_r8)-x0) + (0.5_r8*x0*ex**2)/(rmaxfrac**2)
5000 end function compute_mumin2
5002 real(r8) function compute_ppen(wtwb,D,bogbot,bogtop,rho0j,dpen)
5003 ! ----------------------------------------------------------- !
5004 ! Subroutine to compute critical 'ppen[Pa]<0' ( pressure dis. !
5005 ! from 'ps0(kpen-1)' to the cumulus top where cumulus updraft !
5006 ! vertical velocity is exactly zero ) by considering exact !
5007 ! non-zero fer(kpen). !
5008 ! ----------------------------------------------------------- !
5009 real(r8) :: wtwb, D, bogbot, bogtop, rho0j, dpen
5010 real(r8) :: x0, x1, f, fs, SB, s00
5011 integer :: iteration
5014 SB = ( bogtop - bogbot ) / dpen
5015 ! Sign of slope, 'f' at x = 0
5016 ! If 's00>0', 'w' increases with height.
5017 s00 = bogbot / rho0j - D * wtwb
5019 if( D*dpen .lt. 1.e-8 ) then
5020 if( s00 .ge. 0._r8 ) then
5023 x0 = max(0._r8,min(dpen,-0.5_r8*wtwb/s00))
5026 if( s00 .ge. 0._r8 ) then
5032 f = exp(-2._r8*D*x0)*(wtwb-(bogbot-SB/(2._r8*D))/(D*rho0j)) + &
5033 (SB*x0+bogbot-SB/(2._r8*D))/(D*rho0j)
5034 fs = -2._r8*D*exp(-2._r8*D*x0)*(wtwb-(bogbot-SB/(2._r8*D))/(D*rho0j)) + &
5036 if( fs .ge. 0._r8 ) then
5037 fs = max(fs, 1.e-10_r8)
5039 fs = min(fs,-1.e-10_r8)
5047 compute_ppen = -max(0._r8,min(dpen,x0))
5049 end function compute_ppen
5051 subroutine fluxbelowinv(cbmf,ps0,mkx,kinv,dt,xsrc,xmean,xtopin,xbotin,xflx)
5052 ! ------------------------------------------------------------------------- !
5053 ! Subroutine to calculate turbulent fluxes at and below 'kinv-1' interfaces.!
5054 ! Check in the main program such that input 'cbmf' should not be zero. !
5055 ! If the reconstructed inversion height does not go down below the 'kinv-1' !
5056 ! interface, then turbulent flux at 'kinv-1' interface is simply a product !
5057 ! of 'cmbf' and 'qtsrc-xbot' where 'xbot' is the value at the top interface !
5058 ! of 'kinv-1' layer. This flux is linearly interpolated down to the surface !
5059 ! assuming turbulent fluxes at surface are zero. If reconstructed inversion !
5060 ! height goes down below the 'kinv-1' interface, subsidence warming &drying !
5061 ! measured by 'xtop-xbot', where 'xtop' is the value at the base interface !
5062 ! of 'kinv+1' layer, is added ONLY to the 'kinv-1' layer, using appropriate !
5063 ! mass weighting ( rpinv and rcbmf, or rr = rpinv / rcbmf ) between current !
5064 ! and next provisional time step. Also impose a limiter to enforce outliers !
5065 ! of thermodynamic variables in 'kinv' layer to come back to normal values !
5066 ! at the next step. !
5067 ! ------------------------------------------------------------------------- !
5068 integer, intent(in) :: mkx, kinv
5069 real(r8), intent(in) :: cbmf, dt, xsrc, xmean, xtopin, xbotin
5070 real(r8), intent(in), dimension(0:mkx) :: ps0
5071 real(r8), intent(out), dimension(0:mkx) :: xflx
5073 real(r8) rcbmf, rpeff, dp, rr, pinv_eff, xtop, xbot, pinv, xtop_ori, xbot_ori
5076 dp = ps0(kinv-1) - ps0(kinv)
5078 if( abs(xbotin-xtopin) .le. 1.e-13_r8 ) then
5079 xbot = xbotin - 1.e-13_r8
5080 xtop = xtopin + 1.e-13_r8
5089 ! -------------------------------------- !
5090 ! Compute reconstructed inversion height !
5091 ! -------------------------------------- !
5094 rcbmf = ( cbmf * g * dt ) / dp ! Can be larger than 1 : 'OK'
5096 rpeff = ( xmean - xtop ) / ( xbot - xtop )
5098 rpeff = ( xmean - xtop ) / max(1.e-13_r8, (xbot - xtop) )
5101 rpeff = min( max(0._r8,rpeff), 1._r8 ) ! As of this, 0<= rpeff <= 1
5102 if( rpeff .eq. 0._r8 .or. rpeff .eq. 1._r8 ) then
5106 ! Below two commented-out lines are the old code replacing the above 'if' block.
5107 ! if(rpeff.eq.1) xbot = xmean
5108 ! if(rpeff.eq.0) xtop = xmean
5110 pinv = ps0(kinv-1) - rpeff * dp ! "pinv" before detraining mass
5111 pinv_eff = ps0(kinv-1) + ( rcbmf - rpeff ) * dp ! Effective "pinv" after detraining mass
5112 ! ----------------------------------------------------------------------- !
5113 ! Compute turbulent fluxes. !
5114 ! Below two cases exactly converges at 'kinv-1' interface when rr = 1._r8 !
5115 ! ----------------------------------------------------------------------- !
5117 xflx(k) = cbmf * ( xsrc - xbot ) * ( ps0(0) - ps0(k) ) / ( ps0(0) - pinv )
5119 if( rr .le. 1._r8 ) then
5120 xflx(kinv-1) = xflx(kinv-1) - ( 1._r8 - rr ) * cbmf * ( xtop_ori - xbot_ori )
5124 end subroutine fluxbelowinv
5126 subroutine positive_moisture_single( xlv, xls, mkx, dt, qvmin, qlmin, qimin, dp, qv, ql, qi, s, qvten, qlten, qiten, sten )
5127 ! ------------------------------------------------------------------------------- !
5128 ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
5129 ! force them to be larger than minimum value by (1) condensating water vapor !
5130 ! into liquid or ice, and (2) by transporting water vapor from the very lower !
5131 ! layer. '2._r8' is multiplied to the minimum values for safety. !
5132 ! Update final state variables and tendencies associated with this correction. !
5133 ! If any condensation happens, update (s,t) too. !
5134 ! Note that (qv,ql,qi,s) are final state variables after applying corresponding !
5135 ! input tendencies and corrective tendencies !
5136 ! ------------------------------------------------------------------------------- !
5138 integer, intent(in) :: mkx
5139 real(r8), intent(in) :: xlv, xls
5140 real(r8), intent(in) :: dt, qvmin, qlmin, qimin
5141 real(r8), intent(in) :: dp(mkx)
5142 real(r8), intent(inout) :: qv(mkx), ql(mkx), qi(mkx), s(mkx)
5143 real(r8), intent(inout) :: qvten(mkx), qlten(mkx), qiten(mkx), sten(mkx)
5145 real(r8) dql, dqi, dqv, sum, aa, dum
5147 do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface
5148 dql = max(0._r8,1._r8*qlmin-ql(k))
5149 dqi = max(0._r8,1._r8*qimin-qi(k))
5150 qlten(k) = qlten(k) + dql/dt
5151 qiten(k) = qiten(k) + dqi/dt
5152 qvten(k) = qvten(k) - (dql+dqi)/dt
5153 sten(k) = sten(k) + xlv * (dql/dt) + xls * (dqi/dt)
5156 qv(k) = qv(k) - dql - dqi
5157 s(k) = s(k) + xlv * dql + xls * dqi
5158 dqv = max(0._r8,1._r8*qvmin-qv(k))
5159 qvten(k) = qvten(k) + dqv/dt
5162 qv(k-1) = qv(k-1) - dqv*dp(k)/dp(k-1)
5163 qvten(k-1) = qvten(k-1) - dqv*dp(k)/dp(k-1)/dt
5165 qv(k) = max(qv(k),qvmin)
5166 ql(k) = max(ql(k),qlmin)
5167 qi(k) = max(qi(k),qimin)
5169 ! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
5170 ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
5171 ! preserves column moisture.
5172 if( dqv .gt. 1.e-20_r8 ) then
5175 if( qv(k) .gt. 2._r8*qvmin ) sum = sum + qv(k)*dp(k)
5177 aa = dqv*dp(1)/max(1.e-20_r8,sum)
5178 if( aa .lt. 0.5_r8 ) then
5180 if( qv(k) .gt. 2._r8*qvmin ) then
5183 qvten(k) = qvten(k) - dum/dt
5187 write(iulog,*) 'Full positive_moisture is impossible in uwshcu'
5189 call wrf_message(iulog)
5195 end subroutine positive_moisture_single
5197 subroutine findsp (lchnk, ncol, q, t, p, tsp, qsp)
5199 !-----------------------------------------------------------------------
5202 ! find the wet bulb temperature for a given t and q
5203 ! in a longitude height section
5204 ! wet bulb temp is the temperature and spec humidity that is
5205 ! just saturated and has the same enthalpy
5206 ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
5207 ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
5210 ! a Newton method is used
5211 ! first guess uses an algorithm provided by John Petch from the UKMO
5212 ! we exclude points where the physical situation is unrealistic
5213 ! e.g. where the temperature is outside the range of validity for the
5214 ! saturation vapor pressure, or where the water vapor pressure
5215 ! exceeds the ambient pressure, or the saturation specific humidity is
5220 !-----------------------------------------------------------------------
5222 use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, &
5228 integer, intent(in) :: lchnk ! chunk identifier
5229 integer, intent(in) :: ncol ! number of atmospheric columns
5231 real(r8), intent(in) :: q(pcols,pver) ! water vapor (kg/kg)
5232 real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
5233 real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa)
5237 real(r8), intent(out) :: tsp(pcols,pver) ! saturation temp (K)
5238 real(r8), intent(out) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg)
5242 integer i ! work variable
5243 integer k ! work variable
5244 logical lflg ! work variable
5245 integer iter ! work variable
5246 integer l ! work variable
5247 logical :: error_found
5249 real(r8) omeps ! 1 minus epsilon
5250 real(r8) trinv ! work variable
5251 real(r8) es ! sat. vapor pressure
5252 real(r8) desdt ! change in sat vap pressure wrt temperature
5253 ! real(r8) desdp ! change in sat vap pressure wrt pressure
5254 real(r8) dqsdt ! change in sat spec. hum. wrt temperature
5255 real(r8) dgdt ! work variable
5256 real(r8) g ! work variable
5257 real(r8) weight(pcols) ! work variable
5258 real(r8) hlatsb ! (sublimation)
5259 real(r8) hlatvp ! (vaporization)
5260 real(r8) hltalt(pcols,pver) ! lat. heat. of vap.
5261 real(r8) tterm ! work var.
5262 real(r8) qs ! spec. hum. of water vapor
5263 real(r8) tc ! crit temp of transition to ice
5267 real(r8) t1, q1, dt, dq
5269 real(r8) qvd, a1, tmp
5271 real(r8) r1b, c1, c2, c3
5276 real(r8) enin(pcols), enout(pcols)
5277 real(r8) tlim(pcols)
5279 omeps = 1.0_r8 - epsqs
5280 trinv = 1.0_r8/ttrice
5281 a1 = 7.5_r8*log(10._r8)
5284 dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei
5285 dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei
5286 dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
5287 dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
5288 tt0 = 273.15_r8 ! Freezing temperature
5289 ! tmin = 173.16 ! the coldest temperature we can deal with
5291 ! max number of times to iterate the calculation
5297 ! first guess on the wet bulb temperature
5302 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
5305 call wrf_message(iulog)
5307 write(iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
5309 call wrf_message(iulog)
5313 ! limit the temperature range to that relevant to the sat vap pres tables
5314 #if ( ! defined WACCM_MOZART )
5315 tlim(i) = min(max(t(i,k),173._r8),373._r8)
5317 tlim(i) = min(max(t(i,k),128._r8),373._r8)
5319 es = estblf(tlim(i))
5320 denom = p(i,k) - omeps*es
5324 ! make sure a meaningful calculation is possible
5325 if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
5327 ! Saturation specific humidity
5329 qs = min(epsqs*es/denom,1._r8)
5331 ! "generalized" analytic expression for t derivative of es
5332 ! accurate to within 1 percent for 173.16 < t < 373.16
5334 ! Weighting of hlat accounts for transition from water to ice
5335 ! polynomial expression approximates difference between es over
5336 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
5337 ! -40): required for accurate estimate of es derivative in transition
5338 ! range from ice to water also accounting for change of hlatv with t
5339 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
5342 lflg = (tc >= -ttrice .and. tc < 0.0_r8)
5343 weight(i) = min(-tc*trinv,1.0_r8)
5344 hlatsb = hlatv + weight(i)*hlatf
5345 hlatvp = hlatv - 2369.0_r8*tc
5346 if (tlim(i) < tt0) then
5347 hltalt(i,k) = hlatsb
5349 hltalt(i,k) = hlatvp
5351 enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)
5353 ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
5356 c2 = (tlim(i) + 36._r8)**2
5357 r1b = c2/(c2 + c1*qs)
5359 tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
5361 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
5362 write(iulog,*) ' relative humidity ', q(i,k)/qs
5364 call wrf_message(iulog)
5366 write(iulog,*) ' first guess ', tsp(i,k)
5368 call wrf_message(iulog)
5372 es = estblf(tsp(i,k))
5373 qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
5382 ! now iterate on first guess
5388 if (doit(i) == 0) then
5389 es = estblf(tsp(i,k))
5391 ! Saturation specific humidity
5393 qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
5395 ! "generalized" analytic expression for t derivative of es
5396 ! accurate to within 1 percent for 173.16 < t < 373.16
5398 ! Weighting of hlat accounts for transition from water to ice
5399 ! polynomial expression approximates difference between es over
5400 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
5401 ! -40): required for accurate estimate of es derivative in transition
5402 ! range from ice to water also accounting for change of hlatv with t
5403 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
5406 lflg = (tc >= -ttrice .and. tc < 0.0_r8)
5407 weight(i) = min(-tc*trinv,1.0_r8)
5408 hlatsb = hlatv + weight(i)*hlatf
5409 hlatvp = hlatv - 2369.0_r8*tc
5410 if (tsp(i,k) < tt0) then
5411 hltalt(i,k) = hlatsb
5413 hltalt(i,k) = hlatvp
5416 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
5420 desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
5421 dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
5422 ! g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
5423 g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
5424 dgdt = -(cp + hltalt(i,k)*dqsdt)
5425 t1 = tsp(i,k) - g/dgdt
5426 dt = abs(t1 - tsp(i,k))/t1
5427 tsp(i,k) = max(t1,tmin)
5428 es = estblf(tsp(i,k))
5429 q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
5430 dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
5433 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
5434 write(iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
5436 call wrf_message(iulog)
5442 ! if converged at this point, exclude it from more iterations
5443 if (dt < dttol .and. dq < dqtol) then
5446 enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
5447 ! bail out if we are too near the end of temp range
5448 #if ( ! defined WACCM_MOZART )
5449 if (tsp(i,k) < 174.16_r8) then
5451 if (tsp(i,k) < 130.16_r8) then
5457 end do ! do i = 1,ncol
5459 if (dtm < dttol .and. dqm < dqtol) then
5463 end do ! do l = 1,iter
5466 error_found = .false.
5467 if (dtm > dttol .or. dqm > dqtol) then
5469 if (doit(i) == 0) error_found = .true.
5471 if (error_found) then
5473 if (doit(i) == 0) then
5474 write(iulog,*) ' findsp not converging at point i, k ', i, k
5476 call wrf_message(iulog)
5478 write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
5480 call wrf_message(iulog)
5482 write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
5484 call wrf_message(iulog)
5486 call endrun ('FINDSP')
5492 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
5493 error_found = .true.
5496 if (error_found) then
5498 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
5499 write(iulog,*) ' the enthalpy is not conserved for point ', &
5500 i, k, enin(i), enout(i)
5502 call wrf_message(iulog)
5504 write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
5506 call wrf_message(iulog)
5508 write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
5510 call wrf_message(iulog)
5512 call endrun ('FINDSP')
5517 end do ! level loop (k=1,pver)
5520 end subroutine findsp
5522 ! ------------------------ !
5524 ! End of subroutine blocks !
5526 ! ------------------------ !