diff --git a/interfaces/physics_schemes_interface/rose-meta/um-boundary_layer/HEAD/rose-meta.conf b/interfaces/physics_schemes_interface/rose-meta/um-boundary_layer/HEAD/rose-meta.conf index 82f183e99..501081189 100644 --- a/interfaces/physics_schemes_interface/rose-meta/um-boundary_layer/HEAD/rose-meta.conf +++ b/interfaces/physics_schemes_interface/rose-meta/um-boundary_layer/HEAD/rose-meta.conf @@ -257,15 +257,24 @@ help=Relax to SHARPEST: the stable stability function requested relaxes =above the BL top" is selected this option will also smoothly merge =to the asymptotic mixing length by 3km above the surface. = - =Corrected NTML is as the original but corrects the level + =Original with corrected NTML is as the original but corrects the level =(ntml_local) below which the boundary layer value is used. - =The final option extends this version to use interactively diagnosed - =turbulent layer depths throughout the atmospheric column. + = + =Original with diagnosed turbulent layer depths extends this version to + =use interactively diagnosed turbulent layer depths throughout the + =atmospheric column. + = + =Reduce smoothly towards layer boundaries calculates the lengthscale + =within each turbulent layer (including the surface layer) more smoothly, + =by interpolating between sub- and super- critical levels to find the + =base and top heights, and tapering the mixing length down near the base + =and top. ns=namelist/Science/UM Boundary layer sort-key=Panel-A05 value-titles=Relax to SHARPEST,Original with corrected NTML, - =Original with diagnosed turbulent layer depths -values='to_sharp','ntml_corrected','free_trop_layer' + =Original with diagnosed turbulent layer depths, + =Reduce smoothly towards layer boundaries +values='to_sharp','ntml_corrected','free_trop_layer','smooth_to_boundaries' [namelist:blayer=fric_heating] compulsory=true diff --git a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 index c9572fdc4..09273c6c8 100644 --- a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 +++ b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 @@ -47,6 +47,7 @@ module um_physics_init_mod free_atm_mix, free_atm_mix_to_sharp, & free_atm_mix_ntml_corrected, & free_atm_mix_free_trop_layer, & + free_atm_mix_smooth_to_boundaries, & interp_local, interp_local_gradients, & interp_local_cf_dbdz, & new_kcloudtop, p_unstable, & @@ -341,10 +342,11 @@ subroutine um_physics_init() DynDiag_ZL_corrn, blend_allpoints, ng_stress, & BrownGrant97_limited, BrownGrant97_original, lem_std, & lem_adjust, interactive_fluxes, specified_fluxes_only, & - except_disc_inv, ntml_level_corrn, free_trop_layers, sharpest, & - lem_stability, sg_shear_enh_lambda, l_new_kcloudtop, buoy_integ, & - l_reset_dec_thres, DynDiag_ZL_CuOnly, i_interp_local, & - i_interp_local_gradients, l_noice_in_turb, l_use_var_fixes, & + except_disc_inv, ntml_level_corrn, free_trop_layers, & + smooth_to_bdys, sharpest, lem_stability, sg_shear_enh_lambda, & + l_new_kcloudtop, buoy_integ, l_reset_dec_thres, DynDiag_ZL_CuOnly,& + i_interp_local, i_interp_local_gradients, & + l_noice_in_turb, l_use_var_fixes, & i_interp_local_cf_dbdz, tke_diag_fac, a_ent_2, dec_thres_cloud, & dec_thres_cu, near_neut_z_on_l, blend_gridindep_fa, & specified_fluxes_tstar, buoy_integ_low, num_sweeps_bflux, & @@ -737,6 +739,8 @@ subroutine um_physics_init() local_fa = ntml_level_corrn case(free_atm_mix_free_trop_layer) local_fa = free_trop_layers + case(free_atm_mix_smooth_to_boundaries) + local_fa = smooth_to_bdys end select pstb = 2.0_r_um @@ -771,9 +775,9 @@ subroutine um_physics_init() sg_orog_mixing = sg_shear_enh_lambda end select - ! Switch for corrections to variance diagnostics - l_use_var_fixes = .true. + ! TKE scaling parameter and switch for fixes to variance diagnostics tke_diag_fac = 1.0_r_bl + l_use_var_fixes = .true. zhloc_depth_fac = real(zhloc_depth_fac_in, r_bl) if (topography == topography_horizon) then diff --git a/rose-stem/app/lfric_atm/opt/rose-app-comorph_dev.conf b/rose-stem/app/lfric_atm/opt/rose-app-comorph_dev.conf index b0a923929..6bdb06948 100644 --- a/rose-stem/app/lfric_atm/opt/rose-app-comorph_dev.conf +++ b/rose-stem/app/lfric_atm/opt/rose-app-comorph_dev.conf @@ -1,6 +1,7 @@ [namelist:blayer] dzrad_disc_opt='smooth_1p5' entr_smooth_dec='taper_zh' +free_atm_mix='smooth_to_boundaries' kprof_cu='buoy_integ_low' l_converge_ga=.true. l_use_sml_dsc_fixes=.true. diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit.txt index 55f4a4065..51f0bae73 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 46D8336E -Inner product checksum theta = 518BE459 -Inner product checksum u = 6A85DC8A -Inner product checksum mr1 = 3FD15792 -Inner product checksum mr2 = 38051D13 -Inner product checksum mr3 = 37803DF5 -Inner product checksum mr4 = 37C9CDA4 -Inner product checksum mr5 = 36AE26A1 +Inner product checksum rho = 46D83213 +Inner product checksum theta = 518BE374 +Inner product checksum u = 6A8587F4 +Inner product checksum mr1 = 3FD1948D +Inner product checksum mr2 = 37FCE6A0 +Inner product checksum mr3 = 3763F38E +Inner product checksum mr4 = 37BFEBC6 +Inner product checksum mr5 = 368D33E9 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt index 50c43e284..99847fe4f 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ Inner product checksum rho = 42F8B1D8 -Inner product checksum theta = 4D84DAD8 -Inner product checksum u = 601AD43E -Inner product checksum mr1 = 3C82CADA -Inner product checksum mr2 = 3394F0D5 -Inner product checksum mr3 = 31D3641A +Inner product checksum theta = 4D84DAF4 +Inner product checksum u = 601AD45A +Inner product checksum mr1 = 3C82D8FF +Inner product checksum mr2 = 34023FF4 +Inner product checksum mr3 = 31DE067E Inner product checksum mr4 = 0 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt index f13cb866a..6ad9e33a6 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_azspice_gnu_fast-debug-32bit.txt @@ -1,9 +1,9 @@ Inner product checksum rho = 42F50044 -Inner product checksum theta = 4E044DDE -Inner product checksum u = 5FDA7F3A -Inner product checksum mr1 = 3CC8E42B -Inner product checksum mr2 = 32FA1097 -Inner product checksum mr3 = 30126AEC -Inner product checksum mr4 = 2F0F3F10 +Inner product checksum theta = 4E04586C +Inner product checksum u = 5FD80959 +Inner product checksum mr1 = 3CC8E67A +Inner product checksum mr2 = 33312188 +Inner product checksum mr3 = 30913369 +Inner product checksum mr4 = 2FEC44AB Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_comorph_dev-C12_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_comorph_dev-C12_ex1a_cce_fast-debug-32bit.txt index 37f19beb8..4d0f15533 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_comorph_dev-C12_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_comorph_dev-C12_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 46D8350A -Inner product checksum theta = 518BE2E7 -Inner product checksum u = 6A866B87 -Inner product checksum mr1 = 3FD17AA6 -Inner product checksum mr2 = 37FC8F04 -Inner product checksum mr3 = 37768E4E -Inner product checksum mr4 = 37C47143 -Inner product checksum mr5 = 36A12664 +Inner product checksum rho = 46D83118 +Inner product checksum theta = 518BE292 +Inner product checksum u = 6A860958 +Inner product checksum mr1 = 3FD16668 +Inner product checksum mr2 = 37EBD5F2 +Inner product checksum mr3 = 377556E8 +Inner product checksum mr4 = 37C80521 +Inner product checksum mr5 = 368DA046 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt index c79bb3f92..2499ecd08 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_bomex-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ Inner product checksum rho = 42F8B1D6 -Inner product checksum theta = 4D84DAD0 -Inner product checksum u = 601AD442 -Inner product checksum mr1 = 3C82D804 -Inner product checksum mr2 = 33960C14 -Inner product checksum mr3 = 31C57728 +Inner product checksum theta = 4D84DAEF +Inner product checksum u = 601AD464 +Inner product checksum mr1 = 3C82DF3C +Inner product checksum mr2 = 33F9CD32 +Inner product checksum mr3 = 31EBB59C Inner product checksum mr4 = 0 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt index 27b79d6a7..de76ff022 100644 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_scm_comorph_dev_toga-BiP2x2-50000x50000_ex1a_cce_fast-debug-32bit.txt @@ -1,9 +1,9 @@ Inner product checksum rho = 42F50046 -Inner product checksum theta = 4E044F92 -Inner product checksum u = 5FD92014 -Inner product checksum mr1 = 3CC89574 -Inner product checksum mr2 = 32F9FFA6 -Inner product checksum mr3 = 3036128A -Inner product checksum mr4 = 2EFE31A1 +Inner product checksum theta = 4E045583 +Inner product checksum u = 5FD4E2A6 +Inner product checksum mr1 = 3CC87E3A +Inner product checksum mr2 = 32FB9064 +Inner product checksum mr3 = 308473BC +Inner product checksum mr4 = 308C3F3C Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index 3e4176c53..ea4721681 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -82,7 +82,7 @@ subroutine bdy_expl2 ( & nl_bl_levels, local_fa, free_trop_layers, to_sharp_across_1km, & sbl_op, equilibrium_sbl, one_third, two_thirds, blending_option, & blend_except_cu, blend_cth_shcu_only, sg_shear, sg_shear_enh_lambda, & - max_tke, tke_diag_fac, & + max_tke, tke_diag_fac, smooth_to_bdys, & i_interp_local, i_interp_local_gradients, i_interp_local_cf_dbdz, & shallow_cu_maxtop, sc_cftol, near_neut_z_on_l, zero, one, one_half use cloud_inputs_mod, only: i_rhcpt, forced_cu, i_cld_vn, i_pc2_init_method, & @@ -670,12 +670,12 @@ subroutine bdy_expl2 ( & ! boundary layer (metres) as ! determined from the local ! Richardson number profile. -dsldz(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & - bl_levels), & +zdsc_base(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & + ! Height of base of K_top in DSC +dsldz(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & ! TL+gz/cp gradient between ! levels K and K-1 -dsldz_ga(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & - bl_levels), & +dsldz_ga(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & ! TL+gz/cp gradient between ! levels K and K-1, inc gradient adjust dqwdz(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels) @@ -945,7 +945,7 @@ subroutine bdy_expl2 ( & ! heterogeneous land surface this is poorly defined and we can't use Rib ! from the surface scheme as vertically averaging Ri is numerically ! unstable. So, over land, only the average temperature gradient is used -if (.not. l_use_surf_in_ri) then +if ( .not. l_use_surf_in_ri ) then ! if not using surface variables in Ri (l_use_surf_in_ri=false) we ! extrapolate dbdz itself from level 2, with the sl and qw gradients being ! used in the variance calculations and with i_interp_local_cf_dbdz @@ -1859,7 +1859,7 @@ subroutine bdy_expl2 ( & ! INOUT fields ftl,fqw,zhnl,dzh,cumulus,ntml_nl,w,etadot,t1_sd,q1_sd,wtrac_bl, & ! out fields - rhokmz,rhokhz,rhokm_top,rhokh_top,zhsc, & + rhokmz,rhokhz,rhokm_top,rhokh_top,zhsc,zdsc_base, & unstable,dsc,coupled,sml_disc_inv,dsc_disc_inv, & ntdsc,nbdsc,f_ngstress,tke_nl, & grad_t_adj, grad_q_adj, & @@ -1882,7 +1882,7 @@ subroutine bdy_expl2 ( & !$OMP SHARED(pdims,unstable,fb_surf,cumulus,l_shallow,sml_disc_inv,ntpar, & !$OMP ntml_nl,zhnl,grad_t_adj,grad_q_adj,dsc,dsc_disc_inv,ntdsc,nbdsc, & !$OMP zhsc,dzh,coupled,kent,kent_dsc,t_frac,zrzi,we_lim,t_frac_dsc, & -!$OMP zrzi_dsc,we_lim_dsc,kplume) +!$OMP zdsc_base,zrzi_dsc,we_lim_dsc,kplume) !$OMP do SCHEDULE(STATIC) do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end @@ -1903,6 +1903,7 @@ subroutine bdy_expl2 ( & ntdsc(i,j) = 0 nbdsc(i,j) = 0 zhsc(i,j) = zero + zdsc_base(i,j) = zero coupled(i,j) = .false. ! entrainment variables for non-local tracer mixing kent(i,j) = 2 @@ -1988,8 +1989,9 @@ subroutine bdy_expl2 ( & ! in levels/logicals bl_levels,k_log_layr,BL_diag, & ! in fields - sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m_eff_gb,zhpar,ntpar, & - ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq,delta_smag, & + sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m_eff_gb,zhnl,zhpar,zhsc, & + zdsc_base,ntpar,ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq, & + delta_smag, & ! in/out fields cumulus,weight_1dbl, & ! out fields @@ -2025,7 +2027,7 @@ subroutine bdy_expl2 ( & + (weight2/weight1)*weight_1dbl(i,j,k) end if - if (local_fa == free_trop_layers) then + if (local_fa == free_trop_layers .or. local_fa == smooth_to_bdys) then ! elh already included in rhokh_th so no need to calculate ! here, but interpolate elh separately for diagnostic if (BL_diag%l_elh3d) then @@ -2076,7 +2078,7 @@ subroutine bdy_expl2 ( & rhokh(i,j,k) = elh_rho(i,j,k) * rhokh(i,j,k) end if ! test on sbl_op - end if ! test on local_fa = free_trop_layers + end if ! test on local_fa = free_trop_layers or smooth_to_bdys ! Finally multiply RHOKH by dry density if (l_mr_physics) & diff --git a/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 b/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 index faaa5a96a..ed1b099db 100644 --- a/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 +++ b/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 @@ -120,6 +120,10 @@ module bl_option_mod ! free_trop_layers (=3) => as "ntml_level_corrn" but also diagnose ! FA turbulent layer depths integer, parameter :: free_trop_layers = 3 +! smooth_to_bdys (=4) => as "free_trop_layers" but smoothly interpolate +! between sub- and super-critical model-levels to find depth of turbulent +! layers, and taper the mixing-length down near top and bottom of each layer. +integer, parameter :: smooth_to_bdys = 4 ! 10 Switch to keep local mixing in free atmosphere integer :: Keep_Ri_FA = imdi diff --git a/science/physics_schemes/source/boundary_layer/ex_coef.F90 b/science/physics_schemes/source/boundary_layer/ex_coef.F90 index 92de6db4e..feaf14dde 100644 --- a/science/physics_schemes/source/boundary_layer/ex_coef.F90 +++ b/science/physics_schemes/source/boundary_layer/ex_coef.F90 @@ -26,8 +26,8 @@ subroutine ex_coef ( & ! in levels/logicals bl_levels, k_log_layr, BL_diag, & ! in fields - sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m,zhpar,ntpar, & - ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq, delta_smag, & + sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m,zhnl,zhpar,zhsc,zdsc_base, & + ntpar,ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq,delta_smag, & ! in/out fields cumulus,weight_1dbl, & ! out fields @@ -45,7 +45,7 @@ subroutine ex_coef ( & off, on, sharpest, sharp_sea_long_land, sharp_sea_mes_land, & louis_tails, sharp_sea_louis_land, long_tails, mes_tails, ritrans, & neut_cbl, lambda_min_nml, lambda_max_nml, & - lambda_fac, beta_bl, beta_fa, rlinfac, linear0, & + lambda_fac, beta_bl, beta_fa, rlinfac, linear0, smooth_to_bdys, & to_sharp_across_1km, ntml_level_corrn, free_trop_layers, two_thirds, & blending_option, blend_except_cu, blend_gridindep_fa, blend_cth_shcu_only, & extended_tail, zero, one, one_half @@ -70,56 +70,61 @@ subroutine ex_coef ( & integer, intent(in) :: & bl_levels, & - ! in maximum number of boundary layer levels + ! IN maximum number of boundary layer levels k_log_layr - ! in num of levs requiring log-profile correction + ! IN num of levs requiring log-profile correction integer, intent(in) :: & ntml_nl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Number of model layers in the turbulently + ! IN Number of model layers in the turbulently ! mixed layer as determined from the non-local ! scheme. ntdsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Top level of any decoupled Sc + ! IN Top level of any decoupled Sc nbdsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Bottom level of any decoupled Sc layer. + ! IN Bottom level of any decoupled Sc layer. ntpar(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! in Top level of parcel ascent + ! IN Top level of parcel ascent real(kind=r_bl), intent(in) :: & sigma_h(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Standard deviation of subgrid + ! IN Standard deviation of subgrid ! orography (m) rho_wet_tq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels), & - ! in density on theta levels; + ! IN density on theta levels; ! used in RHOKM so wet density rmlmax2(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels), & - ! in Square of asymptotic mixing length for Smagorinsky scheme + ! IN Square of asymptotic mixing length for Smagorinsky scheme z_uv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels+1), & - ! in Z_UV(K) is height of u level k + ! IN Z_UV(K) is height of u level k z_tq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Z_TQ(K) is height of T,Q level k + ! IN Z_TQ(K) is height of T,Q level k ! NOTE: RI(K) is held at Z_TQ(K-1) + zhnl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & + ! IN Height of top of surface-driven non-local mixing zhpar(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Height of top of initial parcel ascent + ! IN Height of top of initial parcel ascent + zhsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & + zdsc_base(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & + ! IN Base and top heights of decoupled Sc-layer z0m(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Roughness length for momentum (m). + ! IN Roughness length for momentum (m). dvdzm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2:bl_levels), & - ! in Modulus of wind shear. + ! IN Modulus of wind shear. ri(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,2:bl_levels), & - ! in Local Richardson number. + ! IN Local Richardson number. flandg(pdims_s%i_start:pdims_s%i_end,pdims_s%j_start:pdims_s%j_end), & - ! in Land fraction on all tiles. + ! IN Land fraction on all tiles. rneutml_sq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Square of the neutral mixing length for Smagorinsky + ! IN Square of the neutral mixing length for Smagorinsky delta_smag(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end) - ! in delta_x used by Smagorinsky + ! IN delta_x used by Smagorinsky logical, intent(in) :: & l_shallow_cth(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! in Flag to indicate shallow convection based on cl-top + ! IN Flag to indicate shallow convection based on cl-top ! Declaration of new BL diagnostics. type (strnewbldiag), intent(in out) :: BL_diag @@ -140,44 +145,44 @@ subroutine ex_coef ( & real(kind=r_bl), intent(out) :: & rhokm(pdims_s%i_start:pdims_s%i_end,pdims_s%j_start:pdims_s%j_end, & bl_levels), & - ! out Layer K-1 - to - layer K exchange coefficient + ! OUT Layer K-1 - to - layer K exchange coefficient ! for momentum, on UV-grid with first and last ! levels set to "missing data" rhokh(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & bl_levels), & - ! out Layer K-1 - to - layer K exchange coefficient + ! OUT Layer K-1 - to - layer K exchange coefficient ! for scalars (but currently on th-levels) ! On out: still to be multiplied by rho(if l_mr_physics) ! and, for Ri-based scheme, interpolated to ! rho levels in BDY_EXPL2 zh_local(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! out Mixing layer height (m). + ! OUT Mixing layer height (m). integer, intent(out) :: & ntml_local(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! out Number of model layers in the turbulently + ! OUT Number of model layers in the turbulently ! mixed layer as determined from the local ! Richardson number profile. real(kind=r_bl), intent(out) :: & lambda_min, & - ! out Min value of length scale LAMBDA. + ! OUT Min value of length scale LAMBDA. fm_3d(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! out stability function for momentum transport. + ! OUT stability function for momentum transport. ! level 1 value is dummy for use in diagnostics fh_3d(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! out stability function for heat and moisture. + ! OUT stability function for heat and moisture. ! level 1 value is dummy for use in diagnostics tke_loc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels), & - ! out Ri-based scheme diagnosed TKE + ! OUT Ri-based scheme diagnosed TKE elm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,2:bl_levels), & - ! out Mixing length for momentum + ! OUT Mixing length for momentum elh(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,2:bl_levels), & - ! out Mixing length for scalars on theta levels + ! OUT Mixing length for scalars on theta levels elh_rho(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2:bl_levels) - ! out Mixing length for scalars on rho levels + ! OUT Mixing length for scalars on rho levels !----------------------------------------------------------------------- ! Local and other symbolic constants :- @@ -227,7 +232,7 @@ subroutine ex_coef ( & ! Fractional weight applied to ! BL function, vs free atmos turb_length(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & - 2:bl_levels), & + 2:bl_levels+1), & ! Turbulent length scale on theta levels, ! indexed as Ri (m) weight_bltop(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) @@ -280,9 +285,6 @@ subroutine ex_coef ( & lambdah, & ! Asymptotic mixing length for turbulent transport ! of heat/moisture. - lambdah_rho, & - ! Asymptotic mixing length for turbulent transport - ! of heat/moisture on rho levels rlambda_fac, & ! reciprocal of lambda_fac turb_length_layer, & @@ -293,7 +295,15 @@ subroutine ex_coef ( & ! top of boundary layer mixing zfa, & ! height to use beta_fa in blendin - zz + zz, & + ! height + zbot, ztop, & + ! temporaries for height of base and top of layer + lambda_min_use, & + ! Value of minimum mixing length used in certain checks + ! (work variable allows them to be optionally overridden) + interp + ! Vertical interpolation weight integer :: & i, & @@ -538,135 +548,302 @@ subroutine ex_coef ( & ! for a cumulus-capped layer, shear driven mixing is allowed to ! dominate (if ISHEAR_BL=1 selected) !----------------------------------------------------------------------- -!$OMP PARALLEL DEFAULT(none) & -!$OMP SHARED( pdims, ishear_bl, ntml_local, ntpar, cumulus, & -!$OMP bl_levels, lambda_min, rlambda_fac, & -!$OMP turb_length, blending_option, rmlmax2) & -!$OMP private( i, k ) -!$OMP do SCHEDULE(STATIC) +!$OMP PARALLEL do DEFAULT(none) SCHEDULE(STATIC) & +!$OMP SHARED( pdims, ishear_bl, ntml_local, ntpar, cumulus ) & +!$OMP private( i ) do i = pdims%i_start, pdims%i_end if ( ishear_bl == 1 .and. ntml_local(i,j) > ntpar(i,j) ) then cumulus(i,j) = .false. end if end do -!$OMP end do +!$OMP end PARALLEL do !----------------------------------------------------------------------- ! 1.3 Search for sub-critical layers above the PBL and set the ! mixing length to scale with these layer depths -!----------------------------------------------------------------------- +!---------------------------------------------------------------------------- +if ( local_fa==smooth_to_bdys ) then + ! Calculation of turbulence length-scale based on sub-critical layer depth, + ! with interpolation to find accurate base and top height of each layer + ! between neighbouring sub- and super- critical model-levels. + ! Also tapers the length-scale down near the top and bottom of each + ! turbulent layer. + +!$OMP PARALLEL DEFAULT(none) & +!$OMP SHARED( pdims, bl_levels, z_uv, z_tq, turb_length, ri, ricrit , & +!$OMP ntpar, ntml_nl, nbdsc, ntdsc, zh_local, zhnl, zhsc, zdsc_base) & +!$OMP private( i, k, kl, kb, kt, zbot, ztop, interp ) + + ! Initialise turb_length to zero; + ! min and max limits will be imposed after calculation + do k = 2, bl_levels+1 !$OMP do SCHEDULE(STATIC) -do k = 2, bl_levels - do i = pdims%i_start, pdims%i_end - turb_length(i,j,k) = lambda_min*rlambda_fac + do i = pdims%i_start, pdims%i_end + turb_length(i,j,k) = zero + end do +!$OMP end do end do -end do -!$OMP end do NOWAIT -if (blending_option == blend_cth_shcu_only) then - ! use Smag mixing length as background length scale if smaller - ! than lambda_min (ie ignore lambda_min for high res simulations) + +!$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + + kb = 0 + kt = 0 + + do k = 2, bl_levels-1 + + ! Find base of sub-critical layer + if ( k==2 .and. ri(i,j,k) <= ricrit(i,j) ) then + zbot = zero + kb = k + end if + if ( ri(i,j,k) > ricrit(i,j) .and. ri(i,j,k+1) <= ricrit(i,j) ) then + interp = ( ri(i,j,k) - ricrit(i,j) ) / ( ri(i,j,k) - ri(i,j,k+1) ) + interp = max( zero, min( one, interp) ) + zbot = (one-interp) * z_tq(i,j,k-1) & + + interp * z_tq(i,j,k) + kb = k+1 + end if + + ! Find top of sub-critical layer + if ( ri(i,j,k) <= ricrit(i,j) .and. ri(i,j,k+1) > ricrit(i,j) ) then + interp = ( ri(i,j,k) - ricrit(i,j) ) / ( ri(i,j,k) - ri(i,j,k+1) ) + interp = max( zero, min( one, interp) ) + ztop = (one-interp) * z_tq(i,j,k-1) & + + interp * z_tq(i,j,k) + kt = k + end if + if ( k==bl_levels-1 .and. ri(i,j,k+1) <= ricrit(i,j) ) then + ztop = z_uv(i,j,k+1) + kt = k+1 + end if + + if ( kb > 0 .and. kt > 0 ) then + ! Once found base and top of sub-critical layer, set + ! lengthscale on all levels within the layer... + + do kl = kb, kt + ! This formula yields turb_length = ztop - zbot + ! at a height half-way between zbot and ztop, + ! but tapers it towards zero near the base and top + ! (such that L varies quadratically with height, + ! to make the diffusivity profile look a bit like + ! the non-local mixing profiles). + turb_length(i,j,kl) & + = ( z_tq(i,j,kl-1) - zbot ) * ( ztop - z_tq(i,j,kl-1) ) & + * 4.0_r_bl / ( ztop - zbot ) + end do + + ! If this sub-critical layer is at the surface + if ( kb == 2 ) then + ! Reset zh_local consistently + zh_local(i,j) = ztop + ! Remove tapering of turb_length at the surface, as elm + ! already gets tapered near the surface after setting from + ! turb_length, so will be double-counting + do kl = kb, kt + if ( z_tq(i,j,kl-1) < one_half * ztop ) then + turb_length(i,j,kl) = ztop + end if + end do + end if + + ! Reset base and top levels ready to search for next + ! sub-critical layer + kb = 0 + kt = 0 + + end if ! ( kb > 0 .AND. kt > 0 ) + end do ! k = 2, bl_levels-1 + + ! Overwrite with length-scales based on non-local mixing depth + ! where appropriate... + + if ( zhnl(i,j) > zh_local(i,j) ) then + ! If surface-driven non-local mixing extends higher than the + ! found surface-sub-critical layer, reset turb_length + ! consistent with + ! zbot = 0, ztop = zhnl + ztop = zhnl(i,j) + do kl = 2, ntml_nl(i,j) + 2 + if ( z_tq(i,j,kl-1) < one_half * ztop ) then + ! No tapering in lower half of surface-layer + turb_length(i,j,kl) = ztop + else if ( z_tq(i,j,kl-1) < ztop ) then + ! Quadratic tapering in upper half + turb_length(i,j,kl) & + = ( z_tq(i,j,kl-1) ) * ( ztop - z_tq(i,j,kl-1) ) & + * 4.0_r_bl / ztop + end if + end do + end if + + if ( ntdsc(i,j) > 0 ) then + ! If there's a decoupled Sc layer, calculate length-scale + ! consistent with + ! zbot = z_bdsc, ztop = z_tdsc + zbot = zdsc_base(i,j) ! z_uv(i,j,nbdsc(i,j)) + ztop = zhsc(i,j) ! z_uv(i,j,ntdsc(i,j)+1) + ! Only use where larger than the existing length-scale + do kl = nbdsc(i,j), ntdsc(i,j) + 2 + turb_length(i,j,kl) = max( turb_length(i,j,kl), & + ( z_tq(i,j,kl-1) - zbot ) * ( ztop - z_tq(i,j,kl-1) ) & + * 4.0_r_bl / ( ztop - zbot ) ) + end do + end if + + end do ! loop over i +!$OMP end do +!$OMP end PARALLEL + +else ! ( .NOT. local_fa==smooth_to_bdys ) + +!$OMP PARALLEL DEFAULT(none) & +!$OMP SHARED( pdims, bl_levels, lambda_min, rlambda_fac, & +!$OMP turb_length, blending_option, rmlmax2) & +!$OMP private( i, k ) !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels do i = pdims%i_start, pdims%i_end - turb_length(i,j,k) = min( turb_length(i,j,k), sqrt(rmlmax2(i,j,k)) ) + turb_length(i,j,k) = lambda_min*rlambda_fac end do end do !$OMP end do NOWAIT -end if + if (blending_option == blend_cth_shcu_only) then + ! use Smag mixing length as background length scale if smaller + ! than lambda_min (ie ignore lambda_min for high res simulations) +!$OMP do SCHEDULE(STATIC) + do k = 2, bl_levels + do i = pdims%i_start, pdims%i_end + turb_length(i,j,k) = min( turb_length(i,j,k), sqrt(rmlmax2(i,j,k)) ) + end do + end do +!$OMP end do NOWAIT + end if !$OMP end PARALLEL -if (local_fa == free_trop_layers) then + if (local_fa == free_trop_layers) then !$OMP PARALLEL do DEFAULT(none) SCHEDULE(STATIC) & !$OMP SHARED( pdims, bl_levels, ntml_local, ri, ricrit, z_uv, & !$OMP turb_length, rlambda_fac, lambda_min ) & !$OMP private( i, k, subcrit, kb, kt, kl, turb_length_layer ) - do i = pdims%i_start, pdims%i_end - subcrit = .false. - do k = 3, bl_levels - - if ( k > ntml_local(i,j)+1 .and. & - ! we know Ri(ntml_local(i,j)+2) > RiCrit - ri(i,j,k) < ricrit(i,j) .and. .not. subcrit ) then - kb = k ! first level of subcritical Ri in layer - subcrit = .true. - end if - if (ri(i,j,k) >= ricrit(i,j) .and. subcrit ) then - kt = k-1 ! last level of subcritical ri - subcrit = .false. - !--------------------------------------------------------- - ! turb_length(k) is held, with Ri(k), on th-level(k-1) - !--------------------------------------------------------- - turb_length_layer = z_uv(i,j,kt) - z_uv(i,j,kb-1) - do kl = kb, kt - turb_length(i,j,kl) = max( turb_length(i,j,kl), & - min(turb_length_layer,lambda_max_nml*rlambda_fac) ) - end do - end if + do i = pdims%i_start, pdims%i_end + subcrit = .false. + do k = 3, bl_levels + + if ( k > ntml_local(i,j)+1 .and. & + ! we know Ri(ntml_local(i,j)+2) > RiCrit + ri(i,j,k) < ricrit(i,j) .and. .not. subcrit ) then + kb = k ! first level of subcritical Ri in layer + subcrit = .true. + end if + if (ri(i,j,k) >= ricrit(i,j) .and. subcrit ) then + kt = k-1 ! last level of subcritical ri + subcrit = .false. + !--------------------------------------------------------- + ! turb_length(k) is held, with Ri(k), on th-level(k-1) + !--------------------------------------------------------- + turb_length_layer = z_uv(i,j,kt) - z_uv(i,j,kb-1) + do kl = kb, kt + turb_length(i,j,kl) = max( turb_length(i,j,kl), & + min(turb_length_layer,lambda_max_nml*rlambda_fac) ) + end do + end if + end do end do - end do !$OMP end PARALLEL do -end if -!----------------------------------------------------------------------- -! When using turb_length, calculate within the BL -! and use the DSC layer depth as the length scale within a DSC layer -! Remember turb_length(k) is held, with Ri(k), on th-level(k-1) -!----------------------------------------------------------------------- -if (blending_option /= off) then + end if + !----------------------------------------------------------------------- + ! When using turb_length, calculate within the BL + ! and use the DSC layer depth as the length scale within a DSC layer + ! Remember turb_length(k) is held, with Ri(k), on th-level(k-1) + !----------------------------------------------------------------------- + if (blending_option /= off) then !$OMP PARALLEL do DEFAULT(none) SCHEDULE(STATIC) & !$OMP SHARED( bl_levels, pdims, ntml_nl, ntml_local, turb_length, z_uv, & !$OMP zh_local, nbdsc, ntdsc ) & !$OMP private( i, k ) - do k = 2, bl_levels - do i = pdims%i_start, pdims%i_end - if ( k-1 <= max(ntml_nl(i,j),ntml_local(i,j)) ) then - turb_length(i,j,k) = max( turb_length(i,j,k), & - max( z_uv(i,j,ntml_nl(i,j)+1), zh_local(i,j) ) ) - end if - if ( k-1 >= nbdsc(i,j) .and. k-1 <= ntdsc(i,j) ) then - turb_length(i,j,k) = max( turb_length(i,j,k), & - ( z_uv(i,j,ntdsc(i,j)+1)-z_uv(i,j,nbdsc(i,j)) ) ) - end if + do k = 2, bl_levels + do i = pdims%i_start, pdims%i_end + if ( k-1 <= max(ntml_nl(i,j),ntml_local(i,j)) ) then + turb_length(i,j,k) = max( turb_length(i,j,k), & + max( z_uv(i,j,ntml_nl(i,j)+1), zh_local(i,j) ) ) + end if + if ( k-1 >= nbdsc(i,j) .and. k-1 <= ntdsc(i,j) ) then + turb_length(i,j,k) = max( turb_length(i,j,k), & + ( z_uv(i,j,ntdsc(i,j)+1)-z_uv(i,j,nbdsc(i,j)) ) ) + end if + end do end do - end do !$OMP end PARALLEL do -end if + end if +end if ! local_fa /= smooth_to_bdys !----------------------------------------------------------------------- ! 2.0 Loop over levels; calculate the mixing lengths !----------------------------------------------------------------------- +! Under the smooth_to_bdys option, we just impose the min limit on lambda +! once at the end; whereas under other options the limiting is done in several +! places, creating potential loop-holes and oddities. +! Setting lambda_min_use to zero disables the earlier checks, +! defering to the limiting imposed at the end. +if ( local_fa==smooth_to_bdys ) then + lambda_min_use = 0.0_r_bl +else + lambda_min_use = lambda_min +end if + do k = 2, bl_levels !$OMP PARALLEL DEFAULT(none) & -!$OMP PRIVATE(z_scale,i,lambdam,lambdah, & -!$OMP lambdah_rho,vkz,f_log,zz,zht,zfa,beta) & +!$OMP PRIVATE(z_scale,i,lambdam,lambdah,vkz,f_log,zz,zht,zfa,beta) & !$OMP SHARED(k,pdims,ri,ricrit,flandg,ntml_local,ntml_nl,z_tq, & !$OMP l_rp2,lambda_min,par_mezcla_rp,zh_local,turb_length,k_log_layr, & !$OMP z_uv,z0m,elm,elh,elh_rho,blending_option,cumulus,l_shallow_cth,zhpar, & -!$OMP ntdsc,weight_1dbl,weight_bltop,delta_smag,rneutml_sq,BL_diag,local_fa) - !----------------------------------------------------------------- - ! 2.1 Calculate asymptotic mixing lengths LAMBDAM and LAMBDAH - !----------------------------------------------------------------- +!$OMP ntdsc,weight_1dbl,weight_bltop,delta_smag,rneutml_sq,BL_diag,local_fa, & +!$OMP lambda_min_use) !$OMP do SCHEDULE(STATIC) do i = pdims%i_start, pdims%i_end + !------------------------------------------------------------------------ + ! 2.1 Calculate asymptotic mixing lengths, LAMBDAM and LAMBDAH + !------------------------------------------------------------------------ + ! First, effectively set lambda for the surface-based turbulent layer, + ! noting that this value will subsequently be overwritten above the BL + !------------------------------------------------------------------------ + if ( local_fa==smooth_to_bdys ) then + ! Include tapering-down of the length-scale near the top of the + ! surface sub-critical layer (tapered version stored in turb_length) + z_scale = turb_length(i,j,k) + else + ! Use the overall depth of the BL + z_scale = zh_local(i,j) + end if if (l_rp2) then - lambdam = max ( lambda_min , par_mezcla_rp(rp_idx)*zh_local(i,j) ) + lambdam = max ( lambda_min_use , par_mezcla_rp(rp_idx)*z_scale ) else - lambdam = max ( lambda_min , lambda_fac*zh_local(i,j) ) + lambdam = max ( lambda_min_use , lambda_fac*z_scale ) end if !----------------------------------------------------------------- - ! Reduce mixing lengths above BL + ! Reduce lambda above the BL !----------------------------------------------------------------- if (k >= ntml_local(i,j)+2) then - lambdam = lambda_min + lambdam = lambda_min_use end if - - lambdah = lambdam - lambdah_rho = lambdah - - if ( local_fa == free_trop_layers ) then + !----------------------------------------------------------------- + ! Potentially (re)inflate lambda in turbulent layers + !----------------------------------------------------------------- + if ( local_fa == free_trop_layers .or. local_fa == smooth_to_bdys ) then lambdam = max( lambdam, lambda_fac*turb_length(i,j,k) ) - lambdah = max( lambdah, lambda_fac*turb_length(i,j,k) ) - ! lambdah_rho does not need to be recalculated under - ! local_fa option "free_trop_layers" as the full KH profile - ! will be interpolated in bdy_expl2 end if + if ( local_fa==smooth_to_bdys ) then + !--------------------------------------------------------------- + ! Finally impose max and min limits on lambda + ! (remembering that lambda_min_use was set to zero for + ! "smooth_to_bdys", so was effectively inactive earlier) + !--------------------------------------------------------------- + lambdam = max( min( lambdam, lambda_max_nml ), lambda_min ) + end if + !----------------------------------------------------------------- + ! Use the same asymptotic mixing length for scalars + !----------------------------------------------------------------- + lambdah = lambdam !----------------------------------------------------------------------- ! 2.2 Calculate mixing lengths ELH, ELM coincident with RI(K) and so ! at Z_TQ(K-1) @@ -676,7 +853,7 @@ subroutine ex_coef ( & ! Note that ELH_RHO is calculated (on rho levels) for direct inclusion ! in RHOKH and also (as elh) on theta levels for the unstable ! stability functions and inclusion in RHOKH before interpolation - ! (under local_fa option "free_trop_layers"). + ! (under local_fa option "free_trop_layers" or "smooth_to_bdys"). ! To save computing logarithms for all K, the values of ELM and ELH ! are unchanged for K > K_LOG_LAYR. @@ -689,13 +866,13 @@ subroutine ex_coef ( & vkz = vkman * ( z_tq(i,j,k) - z_tq(i,j,k-1) ) f_log = log( ( z_tq(i,j,k) + z0m(i,j) ) / & ( z_tq(i,j,k-1) + z0m(i,j) ) ) - elh_rho(i,j,k) = vkz / ( f_log + vkz/lambdah_rho ) + elh_rho(i,j,k) = vkz / ( f_log + vkz/lambdah ) else vkz = vkman * ( z_tq(i,j,k-1) + z0m(i,j) ) elm(i,j,k) = vkz / (one + vkz/lambdam ) elh(i,j,k) = vkz / (one + vkz/lambdah ) vkz = vkman * ( z_uv(i,j,k) + z0m(i,j) ) - elh_rho(i,j,k) = vkz / (one + vkz/lambdah_rho ) + elh_rho(i,j,k) = vkz / (one + vkz/lambdah ) end if end do !$OMP end do @@ -763,8 +940,8 @@ subroutine ex_coef ( & else weight_1dbl(i,j,k) = one end if - if ( local_fa == free_trop_layers .and. & - ri(i,j,k) < ricrit(i,j) ) then + if ( (local_fa==free_trop_layers .or. local_fa==smooth_to_bdys) & + .and. ri(i,j,k) < ricrit(i,j) ) then ! Except in an elevated turbulent layer where we still use ! the standard blending weight z_scale = turb_length(i,j,k) @@ -1201,7 +1378,7 @@ subroutine ex_coef ( & ! 4.0 Calculate exchange coefficients RHO*KM(K), RHO*KH(K) ! both on TH-level K-1 at this stage (RHOKH will be interpolated ! onto uv-levels and then be multiplied by ELH in BDY_EXPL2 if - ! local_fa is not "free_trop_layers") + ! local_fa is not "free_trop_layers" or "smooth_to_bdys") !------------------------------------------------------------------ if (l_subfilter_vert .or. l_subfilter_horiz) then @@ -1218,14 +1395,13 @@ subroutine ex_coef ( & if (l_mr_physics) then ! Note "RHO" here is always wet density (RHO_WET_TQ) so ! save multiplication of RHOKH to after interpolation - rhokh(i,j,k) = elm(i,j,k) * dvdzm(i,j,k) * fh + rhokh(i,j,k) = elm(i,j,k) * dvdzm(i,j,k) * fh else - rhokh(i,j,k) = rho_wet_tq(i,j,k-1) * elm(i,j,k) * dvdzm(i,j,k) & - * fh + rhokh(i,j,k) = rho_wet_tq(i,j,k-1) * elm(i,j,k) * dvdzm(i,j,k) * fh end if ! If using the FA mixing length profile it is simplest to ! interpolate the full KH profile, including elh (in bdy_expl2) - if (local_fa == free_trop_layers) & + if (local_fa == free_trop_layers .or. local_fa == smooth_to_bdys) & rhokh(i,j,k) = rhokh(i,j,k) * elh(i,j,k) if (BL_diag%l_tke) then diff --git a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 index 32fb1559c..21f3980fe 100644 --- a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 @@ -34,7 +34,7 @@ subroutine kmkhz_9c ( & ! INOUT fields ftl,fqw,zh,dzh,cumulus,ntml,w,etadot,t1_sd,q1_sd,wtrac_bl, & ! out fields - rhokm,rhokh,rhokm_top,rhokh_top,zhsc, & + rhokm,rhokh,rhokm_top,rhokh_top,zhsc,zdsc_base, & unstable,dsc,coupled,sml_disc_inv,dsc_disc_inv, & ntdsc,nbdsc,f_ngstress,tke_nl,grad_t_adj, grad_q_adj, & rhof2, rhofsc, ft_nt, fq_nt, ft_nt_dscb, fq_nt_dscb, & @@ -107,18 +107,18 @@ subroutine kmkhz_9c ( & real(kind=r_bl), intent(in) :: & zmaxb_for_dsc, & zmaxt_for_dsc - ! in Max heights to look for DSC cloud + ! IN Max heights to look for DSC cloud ! base and top integer, intent(in) :: & ntpar(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! in Top level of parcel ascent. + ! IN Top level of parcel ascent. ! Used in convection scheme. ! NOTE: CAN BE > BL_LEVELS-1 logical, intent(in) :: & l_shallow(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! in Flag to indicate shallow + ! IN Flag to indicate shallow ! convection ! Water tracer structure containing 'micro_tends' fields @@ -126,76 +126,76 @@ subroutine kmkhz_9c ( & real(kind=r_bl), intent(in) :: & bq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in A buoyancy parameter for clear air + ! IN A buoyancy parameter for clear air ! on p,T,q-levels (full levels). bt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in A buoyancy parameter for clear air + ! IN A buoyancy parameter for clear air ! on p,T,q-levels (full levels). bqm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in A buoyancy parameter for clear air + ! IN A buoyancy parameter for clear air ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. btm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in A buoyancy parameter for clear air + ! IN A buoyancy parameter for clear air ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. bqm_cld(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels), & - ! in A buoyancy parameter for cloudy air + ! IN A buoyancy parameter for cloudy air ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. btm_cld(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels), & - ! in A buoyancy parameter for cloudy air + ! IN A buoyancy parameter for cloudy air ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. a_qs(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Saturated lapse rate factor + ! IN Saturated lapse rate factor ! on p,T,q-levels (full levels). a_qsm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Saturated lapse rate factor + ! IN Saturated lapse rate factor ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. a_dqsdtm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels), & - ! in Saturated lapse rate factor + ! IN Saturated lapse rate factor ! on intermediate levels (half levels): ! (*,K) elements are k+1/2 values. dqsdt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Partial derivative of QSAT w.r.t. + ! IN Partial derivative of QSAT w.r.t. ! temperature. p(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,0:bl_levels), & - ! in P(*,K) is pressure at full level k. + ! IN P(*,K) is pressure at full level k. qw(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Total water content (kg per kg air). + ! IN Total water content (kg per kg air). tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Liquid/frozen water temperature (K). + ! IN Liquid/frozen water temperature (K). t(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Temperature (K). + ! IN Temperature (K). qcf(tdims_l%i_start:tdims_l%i_end,tdims_l%j_start:tdims_l%j_end, & tdims_l%k_start:bl_levels), & - ! in Cloud ice (kg per kg air) + ! IN Cloud ice (kg per kg air) qcl(tdims_l%i_start:tdims_l%i_end,tdims_l%j_start:tdims_l%j_end, & tdims_l%k_start:bl_levels), & - ! in Cloud liquid water + ! IN Cloud liquid water q(tdims_l%i_start:tdims_l%i_end,tdims_l%j_start:tdims_l%j_end, & tdims_l%k_start:bl_levels), & - ! in specific humidity + ! IN specific humidity cf(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Cloud fractions for boundary levs. + ! IN Cloud fractions for boundary levs. z_tq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Z_tq(*,K) is the height of the + ! IN Z_tq(*,K) is the height of the ! k-th full level above the surface. z_uv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels+1), & - ! in Z_uv(*,K) is the height of level + ! IN Z_uv(*,K) is the height of level ! k-1/2 above the surface (m). dzl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Layer depths (m). DZL(,K) is the + ! IN Layer depths (m). DZL(,K) is the ! distance from layer boundary K-1/2 ! to layer boundary K+1/2. For K=1 ! the lower boundary is the surface. rdz(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Reciprocal of distance between + ! IN Reciprocal of distance between ! full levels (m-1). 1/RDZ(,K) is ! the vertical distance from level ! K-1 to level K, except that for @@ -203,44 +203,44 @@ subroutine kmkhz_9c ( & ! lowest atmospheric full level. rho_mix(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & bl_levels), & - ! in density on UV (ie. rho) levels, + ! IN density on UV (ie. rho) levels, ! used in RHOKH so dry density if ! L_mr_physics is true rho_wet_tq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels), & - ! in density on TQ (ie. theta) levels, + ! IN density on TQ (ie. theta) levels, ! used in RHOKM so wet density rho_mix_tq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & bl_levels) - ! in density on TQ (ie. theta) levels, + ! IN density on TQ (ie. theta) levels, ! used in non-turb flux integration ! so dry density if L_mr_physics is true real(kind=r_bl), intent(in) :: & v_s(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Surface friction velocity (m/s) + ! IN Surface friction velocity (m/s) fb_surf(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Surface buoyancy flux over density + ! IN Surface buoyancy flux over density ! (m^2/s^3). rhostar_gb(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Surface air density in kg/m3 + ! IN Surface air density in kg/m3 z_lcl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Height of lifting condensation + ! IN Height of lifting condensation ! level. zhpar(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Height of top of NTPAR + ! IN Height of top of NTPAR ! NOTE: CAN BE ABOVE BL_LEVELS-1 zh_prev(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! in boundary layer height (m) from + ! IN boundary layer height (m) from ! previous timestep real(kind=r_bl), intent(in) :: & rad_hr(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2,bl_levels), & - ! in (LW,SW) radiative heating rates (K/s) + ! IN (LW,SW) radiative heating rates (K/s) micro_tends(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2, bl_levels) - ! in Tendencies from microphysics + ! IN Tendencies from microphysics ! (TL, K/s; QW, kg/kg/s) ! INOUT arrays @@ -281,31 +281,31 @@ subroutine kmkhz_9c ( & ! out arrays integer, intent(out) :: & ntdsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Top level for turb mixing in + ! OUT Top level for turb mixing in ! cloud layer nbdsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Bottom level of any decoupled + ! OUT Bottom level of any decoupled ! turbulently mixed Sc layer sml_disc_inv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Flags for whether discontinuous + ! OUT Flags for whether discontinuous dsc_disc_inv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out inversions are diagnosed + ! OUT inversions are diagnosed kent(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Grid-levels of SML and DSC + ! OUT Grid-levels of SML and DSC kent_dsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out inversions (for tracer mixing) + ! OUT inversions (for tracer mixing) k_plume(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! out Start grid-level for surface-driven plume + ! OUT Start grid-level for surface-driven plume logical, intent(out) :: & unstable(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Flag to indicate an unstable + ! OUT Flag to indicate an unstable ! surface layer. dsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Flag set if decoupled + ! OUT Flag set if decoupled ! stratocumulus layer found coupled(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! out Flag to indicate Sc layer weakly + ! OUT Flag to indicate Sc layer weakly ! decoupled (implies mixing at SML ! top is through K profiles rather ! than entrainment parametrization) @@ -313,80 +313,82 @@ subroutine kmkhz_9c ( & real(kind=r_bl), intent(out) :: & rhokm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2:bl_levels), & - ! out Non-local turbulent mixing + ! OUT Non-local turbulent mixing ! coefficient for momentum. rhokh(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels), & - ! out Non-local turbulent mixing + ! OUT Non-local turbulent mixing ! coefficient for scalars. rhokm_top(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2:bl_levels), & - ! out Top-down turbulent mixing + ! OUT Top-down turbulent mixing ! coefficient for momentum. rhokh_top(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels), & - ! out Top-down turbulent mixing + ! OUT Top-down turbulent mixing ! coefficient for scalars. tke_nl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels), & - ! out Non-local TKE diag (currently times rho) + ! OUT Non-local TKE diag (currently times rho) f_ngstress(pdims_s%i_start:pdims_s%i_end,pdims_s%j_start:pdims_s%j_end, & 2:bl_levels), & - ! out dimensionless function for + ! OUT dimensionless function for ! non-gradient stresses rhof2(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels), & rhofsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & 2:bl_levels) - ! out f2 and fsc term shape profiles + ! OUT f2 and fsc term shape profiles ! multiplied by rho real(kind=r_bl), intent(out) :: & ft_nt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels+1), & - ! out Non-turbulent heat and moisture + ! OUT Non-turbulent heat and moisture fq_nt(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels+1) ! fluxes (rho*Km/s, rho*m/s) real(kind=r_bl), intent(out) :: & tothf_zh(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Total heat fluxes at inversions + ! OUT Total heat fluxes at inversions tothf_zhsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! (rho*Km/s) totqf_zh(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Total moisture fluxes at + ! OUT Total moisture fluxes at totqf_zhsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! inversions (rho*m/s) ft_nt_dscb(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Non-turbulent heat and moisture + ! OUT Non-turbulent heat and moisture fq_nt_dscb(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! flux at the base of the DSC layer grad_t_adj(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Temperature gradient adjustment + ! OUT Temperature gradient adjustment ! for non-local mixing in unstable ! turbulent boundary layer. grad_q_adj(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! out Humidity gradient adjustment + ! OUT Humidity gradient adjustment ! for non-local mixing in unstable ! turbulent boundary layer. - zhsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! out Cloud layer height (m). - ! The following are used in tracer mixing. - ! At 9B 3 elements were used - here just (i,j,2) is used. + zhsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & + ! OUT Cloud layer height (m). + zdsc_base(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) + ! OUT Height of base of K_top in DSC + +! The following are used in tracer mixing. real(kind=r_bl), intent(out) :: & we_lim(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3), & - ! out rho*entrainment rate implied by + ! OUT rho*entrainment rate implied by ! placing of subsidence zrzi_tr(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3), & - ! out (z-z_base)/(z_i-z_base) + ! OUT (z-z_base)/(z_i-z_base) t_frac_tr(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3), & - ! out a fraction of the timestep + ! OUT a fraction of the timestep we_lim_dsc(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3), & - ! out rho*entrainment rate implied by + ! OUT rho*entrainment rate implied by ! placing of subsidence zrzi_dsc_tr(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3), & - ! out (z-z_base)/(z_i-z_base) + ! OUT (z-z_base)/(z_i-z_base) t_frac_dsc_tr(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,3) - ! out a fraction of the timestep + ! OUT a fraction of the timestep !---------------------------------------------------------------------- ! Local and other symbolic constants :- @@ -655,9 +657,6 @@ subroutine kmkhz_9c ( & rhokh_dsct_ent(pdims%i_start:pdims%i_end, & pdims%j_start:pdims%j_end), & ! DSC top-driven entrainment KH - zdsc_base(pdims%i_start:pdims%i_end, & - pdims%j_start:pdims%j_end), & - ! Height of base of K_top in DSC we_parm(pdims%i_start:pdims%i_end, & pdims%j_start:pdims%j_end), & ! Parametrised entrainment rates (m/s)