From 63529538662f8f382b134edacc1c4745be862d42 Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Tue, 27 Jan 2026 09:43:41 +0000 Subject: [PATCH 1/3] Full debug holds,next fast then vernier timings --- .../site/meto/groups/groups_lfric_atm.cylc | 5 + ...l9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 + ...l9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 + .../source/boundary_layer/excf_nl_9c.F90 | 3246 +++++++++-------- 4 files changed, 1648 insertions(+), 1621 deletions(-) create mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt create mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index 174b5698f..be3d7267a 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -222,6 +222,11 @@ "lfric_atm_clim_gal9_chem_1T-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_clim_gal9_chem_2T-C12_ex1a_cce_fast-debug-32bit", ], + "ex1a_omp_C48_cce_full": [ + "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_2T-C48_MG_ex1a_cce_full-debug-32bit", + "lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit", + ], "ex1a_omp_C48_cce": [ "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_gal9_2T-C48_MG_ex1a_cce_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt new file mode 100644 index 000000000..7d0bbf9c4 --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 48D7FD20 +Inner product checksum theta = 5392A6D8 +Inner product checksum u = 6A97B848 +Inner product checksum mr1 = 41CD0DCE +Inner product checksum mr2 = 39CF631F +Inner product checksum mr3 = 37AD770F +Inner product checksum mr4 = 39604773 +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_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt new file mode 100644 index 000000000..41334f2a4 --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 48D7FD32 +Inner product checksum theta = 5392A6F0 +Inner product checksum u = 6A97B6F0 +Inner product checksum mr1 = 41CD0826 +Inner product checksum mr2 = 39C92715 +Inner product checksum mr3 = 37B19C28 +Inner product checksum mr4 = 396070D2 +Inner product checksum mr5 = 0 +Inner product checksum mr6 = 0 diff --git a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 index 4525be05d..23c8cbc70 100644 --- a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 @@ -700,7 +700,7 @@ subroutine excf_nl_9c ( & integer :: i1, j1, l -integer :: jj !Cache blocking - loop index +integer :: jj, ii !Cache blocking - loop index(s) integer(kind=jpim), parameter :: zhook_in = 0 integer(kind=jpim), parameter :: zhook_out = 1 @@ -765,9 +765,9 @@ subroutine excf_nl_9c ( & !----------------------------------------------------------------------- ! 0. Calculate top-of-b.l. velocity scales and Prandtl number. !----------------------------------------------------------------------- - +j = 1 !$OMP PARALLEL DEFAULT(SHARED) & -!$OMP private(i, j, k, jj, i_wt, c_ws, wstar3, pr_neut, pr_conv, w_m_neut, & +!$OMP private(i, k, ii, jj, i_wt, c_ws, wstar3, pr_neut, pr_conv, w_m_neut, & !$OMP zeta_s_fac, sf_term, sf_shear_term, zeta_r_sq, ir_term, zr, & !$OMP evap_term, dz_inv, l_rad, alpha_t, dr_term, zil_corr, & !$OMP rhokh_ent, frac_top, zh_pr, factor, rhokm_dsct, & @@ -780,102 +780,102 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - rhokh_surf_ent(i,j) = zero - rhokh_top_ent(i,j) = zero - rhokh_dsct_ent(i,j) = zero - rhokm_surf_ent(i,j) = zero - rhokm_top_ent(i,j) = zero - rhokm_dsct_ent(i,j) = zero - v_top(i,j) = zero - v_sum(i,j) = zero - v_top_dsc(i,j) = zero - v_sum_dsc(i,j) = zero - rho_we(i,j) = zero - rho_we_sml(i,j) = zero - rho_we_dsc(i,j) = zero - - if (fb_surf(i,j) >= zero) then - - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) - - ! Free-convective velocity scale cubed - - if (coupled(i,j)) then - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) - else - wstar3 = zhsc(i,j) * fb_surf(i,j) - end if +! do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + + rhokh_surf_ent(i,j) = zero + rhokh_top_ent(i,j) = zero + rhokh_dsct_ent(i,j) = zero + rhokm_surf_ent(i,j) = zero + rhokm_top_ent(i,j) = zero + rhokm_dsct_ent(i,j) = zero + v_top(i,j) = zero + v_sum(i,j) = zero + v_top_dsc(i,j) = zero + v_sum_dsc(i,j) = zero + rho_we(i,j) = zero + rho_we_sml(i,j) = zero + rho_we_dsc(i,j) = zero + + if (fb_surf(i,j) >= zero) then + + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) + + ! Free-convective velocity scale cubed + + if (coupled(i,j)) then + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) else - wstar3 = zh(i,j) * fb_surf(i,j) + wstar3 = zhsc(i,j) * fb_surf(i,j) end if + else + wstar3 = zh(i,j) * fb_surf(i,j) + end if - if (flux_grad == Locketal2000) then + if (flux_grad == Locketal2000) then - ! Turbulent velocity scale for momentum + ! Turbulent velocity scale for momentum - c_ws = 0.25_r_bl - w_m_top(i,j) = (v_s(i,j)*v_s(i,j)*v_s(i,j) + & - c_ws*wstar3)**one_third + c_ws = 0.25_r_bl + w_m_top(i,j) = (v_s(i,j)*v_s(i,j)*v_s(i,j) + & + c_ws*wstar3)**one_third - ! Turbulent Prandtl number and velocity scale for scalars - ! gives 0.375= zero ) then - !----------------------------------------------------------- - ! Calculate the surface buoyancy flux term - !----------------------------------------------------------- - zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) - sf_term = a_ent_1 * max ( zero , & - ( (one - zeta_s_fac) * bflux_surf(i,j) & - + zeta_s_fac * bflux_surf_sat(i,j) ) ) - !----------------------------------------------------------- - ! Calculate the surface shear term - !----------------------------------------------------------- - sf_shear_term = a_ent_shr * v_s(i,j) * v_s(i,j) * v_s(i,j) & - * rho_mix(i,j,k) / zh(i,j) - !----------------------------------------------------------- - ! Calculate the indirect radiative term - !----------------------------------------------------------- - zeta_r_sq = zeta_r(i,j)*zeta_r(i,j) - ir_term = ( bt_top(i,j)*zeta_r_sq + & - btt_top(i,j)*(one-zeta_r_sq) ) & - * a_ent_1 * df_top_over_cp(i,j) - !----------------------------------------------------------- - ! Calculate the evaporative term - !----------------------------------------------------------- - if ( db_top(i,j) > zero) then - zr = sqrt( zc(i,j) / zh(i,j) ) - evap_term = a_ent_2 * rho_mix(i,j,k) & - * chi_s_top(i,j) * chi_s_top(i,j) & - * zr * zr * zr * db_top_cld(i,j) & - * sqrt( zh(i,j) * db_top(i,j) ) - else - evap_term = zero - end if - !----------------------------------------------------------- - ! Combine forcing terms to calculate the representative - ! velocity scales - !----------------------------------------------------------- - v_sum(i,j) = ( (sf_term + sf_shear_term + & - ir_term + evap_term) & - * zh(i,j) /(a_ent_1*rho_mix(i,j,k)) )**one_third - v_top(i,j) = ( (ir_term+evap_term) * zh(i,j) & - / (a_ent_1*rho_mix(i,j,k)) )**one_third - !----------------------------------------------------------- - ! Calculate the direct radiative term - ! can only calculate for DB_TOP > 0 - !----------------------------------------------------------- - if ( db_top(i,j) > zero) then - dz_inv = min( v_sum(i,j)*v_sum(i,j) / db_top(i,j) ,100.0_r_bl ) - l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc(i,j)+rbl_eps) ) - alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback - alpha_t = alpha_t + br_fback(i,j)*(one-alpha_t) - dr_term = btc_top(i,j) * alpha_t * df_top_over_cp(i,j) - !---------------------------------------------------------- - ! Combine terms to calculate the entrainment - ! mixing coefficients - !---------------------------------------------------------- - zil_corr = c_t * ( (sf_term + sf_shear_term + & - ir_term + evap_term) / & - (rho_mix(i,j,k) * sqrt(zh(i,j))) )**two_thirds - - rho_we_sml(i,j) = (sf_term + sf_shear_term & - + ir_term + evap_term + dr_term) & - / ( db_top(i,j) + zil_corr ) - - rhokh_ent = rho_we_sml(i,j)/ rdz(i,j,k) - - frac_top = v_top(i,j) / ( v_top(i,j)+w_h_top(i,j)+rbl_eps ) - if ( l_use_sml_dsc_fixes ) then - ! If bug-fix is on, leave surface-driven entrainment flux set to - ! zero in cumulus layers if we are using the buoyancy-flux - ! integration method to calculate the non-local mixing profile - ! (code originally converged on a consistent negative buoyancy - ! flux at the top and then added the entrainment flux onto - ! this, which is double-counting). - l_apply_surf_ent = ( .not. ( ( kprof_cu==buoy_integ .or. & - kprof_cu==buoy_integ_low ) .and. & - cumulus(i,j) ) ) - else - ! If bug-fix is off, always add on surface-driven entrainment flux - l_apply_surf_ent = .true. - end if - if ( l_apply_surf_ent ) then - ! max to avoid rounding errors giving small negative numbers - rhokh_surf_ent(i,j) = max(rhokh_ent * ( one - frac_top ), zero) - rhokm_surf_ent(i,j) = prandtl_top(i,j) * rhokh_surf_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - end if - rhokh_top_ent(i,j) = rhokh_ent * frac_top - rhokm_top_ent(i,j) = 0.75_r_bl * rhokh_top_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - if (.not. l_use_var_fixes) then - rhokm(i,j,k) = rhokm_surf_ent(i,j) - rhokm_top(i,j,k) = rhokm_top_ent(i,j) - end if - - end if ! test on DB_TOP gt 0 - end if - !---------------------------------------------------------------- - ! then the top of the DSC (if coupled use ZHSC length-scale) - !---------------------------------------------------------------- - if ( ntdsc(i,j) > 0 ) then - k = ntdsc(i,j)+1 - if (coupled(i,j)) then - !-------------------------------------------------------- - ! Calculate the surface buoyancy flux term - !-------------------------------------------------------- - zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) - sf_term = a_ent_1 * max ( zero , & - ( (one - zeta_s_fac) * bflux_surf(i,j) & - + zeta_s_fac * bflux_surf_sat(i,j) ) ) - !-------------------------------------------------------- - ! Calculate the surface shear term - !-------------------------------------------------------- - if (fb_surf(i,j) >= zero) then - sf_shear_term = a_ent_shr * v_s(i,j)*v_s(i,j)*v_s(i,j) & - * rho_mix(i,j,k) / zhsc(i,j) - else - sf_shear_term = zero - end if - if ( entr_smooth_dec == on .or. entr_smooth_dec == entr_taper_zh ) then - ! taper surface terms to zero depending on svl_diff_frac - sf_term = sf_term * svl_diff_frac(i,j) - sf_shear_term = sf_shear_term * svl_diff_frac(i,j) - end if - else - sf_term = zero - sf_shear_term = zero - end if +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + !----------------------------------------------------------------------- + ! 1.2 Calculate top-of-b.l. entrainment mixing coefficients + ! and store b.l. top quantities for later use. + !----------------------------------------------------------------------- + ! FIRST the top of the SML (if not coupled) + !----------------------------------------------- + k = ntml(i,j)+1 + if ( .not. coupled(i,j) .and. fb_surf(i,j) >= zero ) then + !----------------------------------------------------------- + ! Calculate the surface buoyancy flux term + !----------------------------------------------------------- + zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) + sf_term = a_ent_1 * max ( zero , & + ( (one - zeta_s_fac) * bflux_surf(i,j) & + + zeta_s_fac * bflux_surf_sat(i,j) ) ) + !----------------------------------------------------------- + ! Calculate the surface shear term + !----------------------------------------------------------- + sf_shear_term = a_ent_shr * v_s(i,j) * v_s(i,j) * v_s(i,j) & + * rho_mix(i,j,k) / zh(i,j) !----------------------------------------------------------- ! Calculate the indirect radiative term !----------------------------------------------------------- - zeta_r_sq = zeta_r_dsc(i,j)*zeta_r_dsc(i,j) - ir_term = ( bt_dsct(i,j)*zeta_r_sq + & - btt_dsct(i,j)*(one-zeta_r_sq) ) & - * a_ent_1 * df_dsct_over_cp(i,j) + zeta_r_sq = zeta_r(i,j)*zeta_r(i,j) + ir_term = ( bt_top(i,j)*zeta_r_sq + & + btt_top(i,j)*(one-zeta_r_sq) ) & + * a_ent_1 * df_top_over_cp(i,j) !----------------------------------------------------------- ! Calculate the evaporative term !----------------------------------------------------------- - if (db_dsct(i,j) > zero) then - zr = sqrt( zc_dsc(i,j) / dscdepth(i,j) ) - evap_term = a_ent_2 * rho_mix(i,j,k) & - * chi_s_dsct(i,j) * chi_s_dsct(i,j) & - * zr * zr * zr * db_dsct_cld(i,j) & - * sqrt( dscdepth(i,j) * db_dsct(i,j) ) - else - evap_term = zero - end if + if ( db_top(i,j) > zero) then + zr = sqrt( zc(i,j) / zh(i,j) ) + evap_term = a_ent_2 * rho_mix(i,j,k) & + * chi_s_top(i,j) * chi_s_top(i,j) & + * zr * zr * zr * db_top_cld(i,j) & + * sqrt( zh(i,j) * db_top(i,j) ) + else + evap_term = zero + end if !----------------------------------------------------------- ! Combine forcing terms to calculate the representative ! velocity scales !----------------------------------------------------------- - v_sum_dsc(i,j) = ( (sf_term + sf_shear_term + & - ir_term + evap_term) & - * dscdepth(i,j) / (a_ent_1*rho_mix(i,j,k)) )**one_third - v_top_dsc(i,j) =( (ir_term + evap_term) * dscdepth(i,j) / & - (a_ent_1*rho_mix(i,j,k)) )**one_third + v_sum(i,j) = ( (sf_term + sf_shear_term + & + ir_term + evap_term) & + * zh(i,j) /(a_ent_1*rho_mix(i,j,k)) )**one_third + v_top(i,j) = ( (ir_term+evap_term) * zh(i,j) & + / (a_ent_1*rho_mix(i,j,k)) )**one_third !----------------------------------------------------------- ! Calculate the direct radiative term + ! can only calculate for DB_TOP > 0 !----------------------------------------------------------- - if (db_dsct(i,j) > zero) then - dz_inv = min( v_sum_dsc(i,j)*v_sum_dsc(i,j)/db_dsct(i,j), & - 100.0_r_bl ) - l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc_dsc(i,j)+one) ) - alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback - alpha_t = alpha_t + br_fback_dsc(i,j)*(one-alpha_t) - dr_term = btc_top(i,j) * alpha_t * df_dsct_over_cp(i,j) - !---------------------------------------------------------- - ! Finally combine terms to calculate the entrainment - ! rate and mixing coefficients - !---------------------------------------------------------- - zil_corr = c_t * ( (sf_term + sf_shear_term + & - ir_term + evap_term) / & - (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds - rho_we_dsc(i,j) = ( sf_term + sf_shear_term & - + ir_term + evap_term + dr_term ) & - / ( db_dsct(i,j) + zil_corr ) - rhokh_dsct_ent(i,j) = rho_we_dsc(i,j)/ rdz(i,j,k) - rhokm_dsct_ent(i,j) = 0.75_r_bl * rhokh_dsct_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - if (.not. l_use_var_fixes) rhokm_top(i,j,k) = rhokm_dsct_ent(i,j) - end if ! test on DB_DSCT gt 0 + if ( db_top(i,j) > zero) then + dz_inv = min( v_sum(i,j)*v_sum(i,j) / db_top(i,j) ,100.0_r_bl ) + l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc(i,j)+rbl_eps) ) + alpha_t = one - exp(-one_half*dz_inv/l_rad) + ! Make enhancement due to buoyancy reversal feedback + alpha_t = alpha_t + br_fback(i,j)*(one-alpha_t) + dr_term = btc_top(i,j) * alpha_t * df_top_over_cp(i,j) + !---------------------------------------------------------- + ! Combine terms to calculate the entrainment + ! mixing coefficients + !---------------------------------------------------------- + zil_corr = c_t * ( (sf_term + sf_shear_term + & + ir_term + evap_term) / & + (rho_mix(i,j,k) * sqrt(zh(i,j))) )**two_thirds + + rho_we_sml(i,j) = (sf_term + sf_shear_term & + + ir_term + evap_term + dr_term) & + / ( db_top(i,j) + zil_corr ) + + rhokh_ent = rho_we_sml(i,j)/ rdz(i,j,k) + + frac_top = v_top(i,j) / ( v_top(i,j)+w_h_top(i,j)+rbl_eps ) + if ( l_use_sml_dsc_fixes ) then + ! If bug-fix is on, leave surface-driven entrainment flux set to + ! zero in cumulus layers if we are using the buoyancy-flux + ! integration method to calculate the non-local mixing profile + ! (code originally converged on a consistent negative buoyancy + ! flux at the top and then added the entrainment flux onto + ! this, which is double-counting). + l_apply_surf_ent = ( .not. ( ( kprof_cu==buoy_integ .or. & + kprof_cu==buoy_integ_low ) .and. & + cumulus(i,j) ) ) + else + ! If bug-fix is off, always add on surface-driven entrainment flux + l_apply_surf_ent = .true. + end if + if ( l_apply_surf_ent ) then + ! max to avoid rounding errors giving small negative numbers + rhokh_surf_ent(i,j) = max(rhokh_ent * ( one - frac_top ), zero) + rhokm_surf_ent(i,j) = prandtl_top(i,j) * rhokh_surf_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + end if + rhokh_top_ent(i,j) = rhokh_ent * frac_top + rhokm_top_ent(i,j) = 0.75_r_bl * rhokh_top_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + if (.not. l_use_var_fixes) then + rhokm(i,j,k) = rhokm_surf_ent(i,j) + rhokm_top(i,j,k) = rhokm_top_ent(i,j) + end if + + end if ! test on DB_TOP gt 0 + end if + !---------------------------------------------------------------- + ! then the top of the DSC (if coupled use ZHSC length-scale) + !---------------------------------------------------------------- + if ( ntdsc(i,j) > 0 ) then + k = ntdsc(i,j)+1 + if (coupled(i,j)) then + !-------------------------------------------------------- + ! Calculate the surface buoyancy flux term + !-------------------------------------------------------- + zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) + sf_term = a_ent_1 * max ( zero , & + ( (one - zeta_s_fac) * bflux_surf(i,j) & + + zeta_s_fac * bflux_surf_sat(i,j) ) ) + !-------------------------------------------------------- + ! Calculate the surface shear term + !-------------------------------------------------------- + if (fb_surf(i,j) >= zero) then + sf_shear_term = a_ent_shr * v_s(i,j)*v_s(i,j)*v_s(i,j) & + * rho_mix(i,j,k) / zhsc(i,j) + else + sf_shear_term = zero + end if + if ( entr_smooth_dec == on .or. entr_smooth_dec == entr_taper_zh ) then + ! taper surface terms to zero depending on svl_diff_frac + sf_term = sf_term * svl_diff_frac(i,j) + sf_shear_term = sf_shear_term * svl_diff_frac(i,j) + end if + else + sf_term = zero + sf_shear_term = zero end if - end do + !----------------------------------------------------------- + ! Calculate the indirect radiative term + !----------------------------------------------------------- + zeta_r_sq = zeta_r_dsc(i,j)*zeta_r_dsc(i,j) + ir_term = ( bt_dsct(i,j)*zeta_r_sq + & + btt_dsct(i,j)*(one-zeta_r_sq) ) & + * a_ent_1 * df_dsct_over_cp(i,j) + !----------------------------------------------------------- + ! Calculate the evaporative term + !----------------------------------------------------------- + if (db_dsct(i,j) > zero) then + zr = sqrt( zc_dsc(i,j) / dscdepth(i,j) ) + evap_term = a_ent_2 * rho_mix(i,j,k) & + * chi_s_dsct(i,j) * chi_s_dsct(i,j) & + * zr * zr * zr * db_dsct_cld(i,j) & + * sqrt( dscdepth(i,j) * db_dsct(i,j) ) + else + evap_term = zero + end if + !----------------------------------------------------------- + ! Combine forcing terms to calculate the representative + ! velocity scales + !----------------------------------------------------------- + v_sum_dsc(i,j) = ( (sf_term + sf_shear_term + & + ir_term + evap_term) & + * dscdepth(i,j) / (a_ent_1*rho_mix(i,j,k)) )**one_third + v_top_dsc(i,j) =( (ir_term + evap_term) * dscdepth(i,j) / & + (a_ent_1*rho_mix(i,j,k)) )**one_third + !----------------------------------------------------------- + ! Calculate the direct radiative term + !----------------------------------------------------------- + if (db_dsct(i,j) > zero) then + dz_inv = min( v_sum_dsc(i,j)*v_sum_dsc(i,j)/db_dsct(i,j), & + 100.0_r_bl ) + l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc_dsc(i,j)+one) ) + alpha_t = one - exp(-one_half*dz_inv/l_rad) + ! Make enhancement due to buoyancy reversal feedback + alpha_t = alpha_t + br_fback_dsc(i,j)*(one-alpha_t) + dr_term = btc_top(i,j) * alpha_t * df_dsct_over_cp(i,j) + !---------------------------------------------------------- + ! Finally combine terms to calculate the entrainment + ! rate and mixing coefficients + !---------------------------------------------------------- + zil_corr = c_t * ( (sf_term + sf_shear_term + & + ir_term + evap_term) / & + (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds + rho_we_dsc(i,j) = ( sf_term + sf_shear_term & + + ir_term + evap_term + dr_term ) & + / ( db_dsct(i,j) + zil_corr ) + rhokh_dsct_ent(i,j) = rho_we_dsc(i,j)/ rdz(i,j,k) + rhokm_dsct_ent(i,j) = 0.75_r_bl * rhokh_dsct_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + if (.not. l_use_var_fixes) rhokm_top(i,j,k) = rhokm_dsct_ent(i,j) + end if ! test on DB_DSCT gt 0 + end if end do +!end do !$OMP end do ! If there is no turbulence generation in DSC layer, ignore it. !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( dsc(i,j) .and. v_top_dsc(i,j) <= zero ) then - dsc_removed(i,j) = 1 - dsc(i,j) = .false. - ntdsc(i,j) = 0 - zhsc(i,j) = zero - zc_dsc(i,j) = zero - dscdepth(i,j) = zero - coupled(i,j) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + if ( dsc(i,j) .and. v_top_dsc(i,j) <= zero ) then + dsc_removed(i,j) = 1 + dsc(i,j) = .false. + ntdsc(i,j) = 0 + zhsc(i,j) = zero + zc_dsc(i,j) = zero + dscdepth(i,j) = zero + coupled(i,j) = .false. + end if end do +!end do !$OMP end do ! ---------------------------------------------------------------------- @@ -1144,13 +1144,13 @@ subroutine excf_nl_9c ( & ! Need to use pdims%k_end(=bl_levels-1) here because bl_levels ! in this routine is actually nl_bl_levels (< bl_levels) do k = pdims%k_start, pdims%k_end - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - wbmix(i,j,k) = zero ! WB if were diag as well-mixed - wbend(i,j,k) = zero ! WB after dec diag - wbend_sml(i,j,k) = zero ! WB after ksurf_iterate - end do + ! do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + wbmix(i,j,k) = zero ! WB if were diag as well-mixed + wbend(i,j,k) = zero ! WB after dec diag + wbend_sml(i,j,k) = zero ! WB after ksurf_iterate end do + ! end do end do !$OMP end do end if ! model_type @@ -1160,40 +1160,40 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - zsml_top(i,j) = zh(i,j) - zsml_base(i,j) = 0.1_r_bl * zh(i,j) - zdsc_base(i,j) = 0.1_r_bl * zhsc(i,j) - z_top_lim(i,j) = zero ! initialising - z_inv(i,j) = zero ! inversion height (top of K profiles) - zwb0(i,j) = zero ! height at which WB goes to zero - wbp_int(i,j) = zero - wbn_int(i,j) = zero - wb_surf_int(i,j) = zero - wb_dzrad_int(i,j) = zero - dzrad(i,j) = 100.0_r_bl - tothf_zi(i,j) = zero - totqf_zi(i,j) = zero - kstatus(i,j)= .true. - kwb0(i,j) = 2 - ntop(i,j) = -1 - ksurf(i,j) = 1 - dec_thres(i,j) = dec_thres_cloud ! use cloudy by default - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + zsml_top(i,j) = zh(i,j) + zsml_base(i,j) = 0.1_r_bl * zh(i,j) + zdsc_base(i,j) = 0.1_r_bl * zhsc(i,j) + z_top_lim(i,j) = zero ! initialising + z_inv(i,j) = zero ! inversion height (top of K profiles) + zwb0(i,j) = zero ! height at which WB goes to zero + wbp_int(i,j) = zero + wbn_int(i,j) = zero + wb_surf_int(i,j) = zero + wb_dzrad_int(i,j) = zero + dzrad(i,j) = 100.0_r_bl + tothf_zi(i,j) = zero + totqf_zi(i,j) = zero + kstatus(i,j)= .true. + kwb0(i,j) = 2 + ntop(i,j) = -1 + ksurf(i,j) = 1 + dec_thres(i,j) = dec_thres_cloud ! use cloudy by default end do +!end do !$OMP end do ! Find KSURF, the first theta-level above the surface layer !$OMP do SCHEDULE(STATIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( z_tq(i,j,k-1) < 0.1_r_bl*zh(i,j) ) ksurf(i,j) = k - end do + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if ( z_tq(i,j,k-1) < 0.1_r_bl*zh(i,j) ) ksurf(i,j) = k + !end do end do end do end do @@ -1205,117 +1205,117 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - test_well_mixed(i,j) = .false. - ksurf_iterate(i,j)= .false. - ktop_iterate(i,j) = .false. - if ( ntdsc(i,j) > 2 ) then - ktop_iterate(i,j) = .true. - end if - end do ! I -end do ! J +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + test_well_mixed(i,j) = .false. + ksurf_iterate(i,j)= .false. + ktop_iterate(i,j) = .false. + if ( ntdsc(i,j) > 2 ) then + ktop_iterate(i,j) = .true. + end if +end do ! I +!end do ! J !$OMP end do !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( bflux_surf(i,j) > zero) then - ! can only be coupled to an unstable SML - if ( ntdsc(i,j) > 2 ) then - ! well-resolved DSC layer - ! ...test for recoupling with SML - test_well_mixed(i,j) = .true. - z_inv(i,j) = zhsc(i,j) - z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) - cf_ml(i,j) = cf_dsc(i,j) - v_ktop(i,j) = v_top_dsc(i,j) - v_ksum(i,j) = v_sum_dsc(i,j) - df_ctop(i,j)= df_dsct_over_cp(i,j) - rho_we(i,j) = rho_we_dsc(i,j) - dsl(i,j) = dsl_dsc(i,j) - dqw(i,j) = dqw_dsc(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) - ntop(i,j) = ntdsc(i,j) - ! assuming wb goes to zero by the lowest of ZH - ! or cloud-base, but above surface layer - else if ( .not. dsc(i,j) .and. .not. cumulus(i,j) .and. & - ntml(i,j) > 2) then - ! well-resolved SML - ! ...test for decoupling - ! Note: code can only deal with one DSC layer at a time so - ! can't decouple SML if a DSC layer already exists. - ! --------------------------------------------------------- - ! If the BL layer is cloud-free then use a less restrictive - ! threshold - ideally, the parcel ascent would have - ! found the correct BL top in this case but this test is - ! kept to keep negative buoyancy fluxes under control - ! (ie. DEC_THRES_CLEAR=1 ensures wbn_int < |wbp_int|) - if (abs(zc(i,j)) < rbl_eps) dec_thres(i,j) = dec_thres_clear - ! --------------------------------------------------------- - test_well_mixed(i,j) = .true. - z_inv(i,j) = zh(i,j) - z_cbase(i,j)= z_inv(i,j) - zc(i,j) - cf_ml(i,j) = cf_sml(i,j) - v_ktop(i,j) = v_top(i,j) - v_ksum(i,j) = v_sum(i,j) - df_ctop(i,j)= df_top_over_cp(i,j) - rho_we(i,j) = rho_we_sml(i,j) - dsl(i,j) = dsl_sml(i,j) - dqw(i,j) = dqw_sml(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zh(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zh(i,j) - ntop(i,j) = ntml(i,j) - end if +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + if ( bflux_surf(i,j) > zero) then + ! can only be coupled to an unstable SML + if ( ntdsc(i,j) > 2 ) then + ! well-resolved DSC layer + ! ...test for recoupling with SML + test_well_mixed(i,j) = .true. + z_inv(i,j) = zhsc(i,j) + z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) + cf_ml(i,j) = cf_dsc(i,j) + v_ktop(i,j) = v_top_dsc(i,j) + v_ksum(i,j) = v_sum_dsc(i,j) + df_ctop(i,j)= df_dsct_over_cp(i,j) + rho_we(i,j) = rho_we_dsc(i,j) + dsl(i,j) = dsl_dsc(i,j) + dqw(i,j) = dqw_dsc(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) + ntop(i,j) = ntdsc(i,j) + ! assuming wb goes to zero by the lowest of ZH + ! or cloud-base, but above surface layer + else if ( .not. dsc(i,j) .and. .not. cumulus(i,j) .and. & + ntml(i,j) > 2) then + ! well-resolved SML + ! ...test for decoupling + ! Note: code can only deal with one DSC layer at a time so + ! can't decouple SML if a DSC layer already exists. + ! --------------------------------------------------------- + ! If the BL layer is cloud-free then use a less restrictive + ! threshold - ideally, the parcel ascent would have + ! found the correct BL top in this case but this test is + ! kept to keep negative buoyancy fluxes under control + ! (ie. DEC_THRES_CLEAR=1 ensures wbn_int < |wbp_int|) + if (abs(zc(i,j)) < rbl_eps) dec_thres(i,j) = dec_thres_clear + ! --------------------------------------------------------- + test_well_mixed(i,j) = .true. + z_inv(i,j) = zh(i,j) + z_cbase(i,j)= z_inv(i,j) - zc(i,j) + cf_ml(i,j) = cf_sml(i,j) + v_ktop(i,j) = v_top(i,j) + v_ksum(i,j) = v_sum(i,j) + df_ctop(i,j)= df_top_over_cp(i,j) + rho_we(i,j) = rho_we_sml(i,j) + dsl(i,j) = dsl_sml(i,j) + dqw(i,j) = dqw_sml(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zh(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zh(i,j) + ntop(i,j) = ntml(i,j) end if - end do ! I -end do ! J + end if +end do ! I +!end do ! J !$OMP end do if ( l_converge_ga ) then ! Set gradient adjustment terms consistent with zh = z_inv !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) ) then - ! The gradient adjustment scales with 1/zh. The terms were calculated - ! in kmkhz_9c using zh_prev (mixed-layer height from previous timestep), - ! but here we are setting the mixed-layer height to z_inv, - ! so need to scale the gradient adjustment by zh_prev/z_inv - ga_fac(i,j) = zh_prev(i,j) / z_inv(i,j) - else - ! Set to 1 where not used. - ga_fac(i,j) = one - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) ) then + ! The gradient adjustment scales with 1/zh. The terms were calculated + ! in kmkhz_9c using zh_prev (mixed-layer height from previous timestep), + ! but here we are setting the mixed-layer height to z_inv, + ! so need to scale the gradient adjustment by zh_prev/z_inv + ga_fac(i,j) = zh_prev(i,j) / z_inv(i,j) + else + ! Set to 1 where not used. + ga_fac(i,j) = one + end if end do + !end do !$OMP end do !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Updated gradient-adjusted db is non-adjusted value + ga_fac times - ! the original adjustment. - db_ga_dry_n(i,j,k) = db_noga_dry(i,j,k) & - + ( db_ga_dry(i,j,k) - db_noga_dry(i,j,k) ) * ga_fac(i,j) - db_ga_cld_n(i,j,k) = db_noga_cld(i,j,k) & - + ( db_ga_cld(i,j,k) - db_noga_cld(i,j,k) ) * ga_fac(i,j) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! Updated gradient-adjusted db is non-adjusted value + ga_fac times + ! the original adjustment. + db_ga_dry_n(i,j,k) = db_noga_dry(i,j,k) & + + ( db_ga_dry(i,j,k) - db_noga_dry(i,j,k) ) * ga_fac(i,j) + db_ga_cld_n(i,j,k) = db_noga_cld(i,j,k) & + + ( db_ga_cld(i,j,k) - db_noga_cld(i,j,k) ) * ga_fac(i,j) end do + ! end do end do !$OMP end do else ! ( .not. l_converge_ga ) ! Use original gradient adjustment set consistent with zh_prev !$OMP do SCHEDULE(STATIC) do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - db_ga_dry_n(i,j,k) = db_ga_dry(i,j,k) - db_ga_cld_n(i,j,k) = db_ga_cld(i,j,k) - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + db_ga_dry_n(i,j,k) = db_ga_dry(i,j,k) + db_ga_cld_n(i,j,k) = db_ga_cld(i,j,k) end do + !end do end do !$OMP end do end if ! ( l_converge_ga ) @@ -1323,18 +1323,18 @@ subroutine excf_nl_9c ( & ! Find kwb0, level with lowest positive cloud-free buoyancy gradient !$OMP do SCHEDULE(STATIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%j_start, pdims%i_end, bl_segment_size do k = 2, bl_levels - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if (kstatus(i,j)) then - if ( (db_ga_dry_n(i,j,k) <= zero) .or. & - (k >= ntml(i,j)) ) then - kstatus(i,j)=.false. - kwb0(i,j)=k - end if + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if (kstatus(i,j)) then + if ( (db_ga_dry_n(i,j,k) <= zero) .or. & + (k >= ntml(i,j)) ) then + kstatus(i,j)=.false. + kwb0(i,j)=k end if - end do + end if + !end do end do end do end do @@ -1353,99 +1353,99 @@ subroutine excf_nl_9c ( & ! ---------------------------------------------------------------------- !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) ) then - dzrad(i,j) = 100.0_r_bl - - if ( dzrad_disc_opt == dzrad_ntm1 ) then - ! Original discretisation of cloud-top radiatively-cooled layer; - ! set the base of the layer at theta-level ntdsc-1 - ! (which is between 1.5 and 2.5 model-levels below cloud-top, - ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), - ! then lower it by extra whole model-levels if needed to increase the - ! depth to 100m. +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) ) then + dzrad(i,j) = 100.0_r_bl - ntop(i,j) = ntop(i,j) - 1 - if ( ntop(i,j) > 2 ) then - inner_loop1: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) - ntop(i,j) = ntop(i,j) - 1 - if (ntop(i,j) == 2 ) then - exit inner_loop1 - end if - end do inner_loop1 - end if - dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) - - else if ( dzrad_disc_opt == dzrad_1p5dz ) then - ! Alternative method; - ! set the base of the layer 1.5 model-levels below cloud-top so that - ! it always moves smoothly with zhsc instead of jumping suddenly - ! (1.5 levels is the minimum distance below zhsc which will always - ! keep the base of the layer below theta-level ntdsc, which is a - ! requirement of the discretisation). - ! The depth is also forced to be at least 100m. - - ! Find smoothly-varying height 1.5 model-levels below z_inv - k_inv = ntop(i,j) - if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then - interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) - z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) - else - interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) - z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) - end if - ! Go down to 100m below z_inv (or cloud-base) if that is lower, - ! but don't allow lower than theta-level 1 - z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) - z_pr = max( z_pr, z_tq(i,j,1) ) - ! Set ntop to the rho-level straddling the base of the layer - do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) + if ( dzrad_disc_opt == dzrad_ntm1 ) then + ! Original discretisation of cloud-top radiatively-cooled layer; + ! set the base of the layer at theta-level ntdsc-1 + ! (which is between 1.5 and 2.5 model-levels below cloud-top, + ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), + ! then lower it by extra whole model-levels if needed to increase the + ! depth to 100m. + + ntop(i,j) = ntop(i,j) - 1 + if ( ntop(i,j) > 2 ) then + inner_loop1: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) ntop(i,j) = ntop(i,j) - 1 - end do - ! Set new cloud-top radiatively-cooled layer-depth - dzrad(i,j) = z_inv(i,j) - z_pr - - ! Fraction of rho-level ntop below the dzrad layer - interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) - ! Add on SL and qw differences between the discontinuous inversion base - ! and the base of the dzrad layer - dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & - - interp *sl(i,j,ntop(i,j)) - dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & - - interp *qw(i,j,ntop(i,j)) - - end if ! ( dzrad_disc_opt ) - - wsl_dzrad_int = dzrad(i,j) * & - ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) - wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) - - wb_dzrad_int(i,j) = wsl_dzrad_int * ( & - (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & - + wqw_dzrad_int * ( & - (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) - wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) - wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) - - ! Include WB_DZRAD_INT in WBP_INT as it set to be >0 - wbp_int(i,j) = wbp_int(i,j) + wb_dzrad_int(i,j) + if (ntop(i,j) == 2 ) then + exit inner_loop1 + end if + end do inner_loop1 + end if + dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) + + else if ( dzrad_disc_opt == dzrad_1p5dz ) then + ! Alternative method; + ! set the base of the layer 1.5 model-levels below cloud-top so that + ! it always moves smoothly with zhsc instead of jumping suddenly + ! (1.5 levels is the minimum distance below zhsc which will always + ! keep the base of the layer below theta-level ntdsc, which is a + ! requirement of the discretisation). + ! The depth is also forced to be at least 100m. + + ! Find smoothly-varying height 1.5 model-levels below z_inv + k_inv = ntop(i,j) + if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then + interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + z_pr = (one-interp) * z_tq(i,j,k_inv-1) & + + interp * z_uv(i,j,k_inv) + else + interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + z_pr = (one-interp) * z_uv(i,j,k_inv) & + + interp * z_tq(i,j,k_inv) + end if + ! Go down to 100m below z_inv (or cloud-base) if that is lower, + ! but don't allow lower than theta-level 1 + z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) + z_pr = max( z_pr, z_tq(i,j,1) ) + ! Set ntop to the rho-level straddling the base of the layer + do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) + ntop(i,j) = ntop(i,j) - 1 + end do + ! Set new cloud-top radiatively-cooled layer-depth + dzrad(i,j) = z_inv(i,j) - z_pr + + ! Fraction of rho-level ntop below the dzrad layer + interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + ! Add on SL and qw differences between the discontinuous inversion base + ! and the base of the dzrad layer + dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & + - interp *sl(i,j,ntop(i,j)) + dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & + - interp *qw(i,j,ntop(i,j)) + + end if ! ( dzrad_disc_opt ) + + wsl_dzrad_int = dzrad(i,j) * & + ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) + wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) + + wb_dzrad_int(i,j) = wsl_dzrad_int * ( & + (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & + + wqw_dzrad_int * ( & + (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) + wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) + wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) + + ! Include WB_DZRAD_INT in WBP_INT as it set to be >0 + wbp_int(i,j) = wbp_int(i,j) + wb_dzrad_int(i,j) - else + else - wb_dzrad_int(i,j) = -one ! To identify not calculated + wb_dzrad_int(i,j) = -one ! To identify not calculated - end if - end do ! I -end do ! J + end if +end do ! I +!end do ! J !$OMP end do ! For WB diagnostics, convert integrated WB to uniform profile @@ -1453,19 +1453,19 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(STATIC) do k=2, bl_levels - do j=pdims%j_start, pdims%j_end - do i=pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) .and. k >= ntop(i,j)+1 ) then - if ( k <= ntdsc(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - else if ( k <= ntml(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if + !do j=pdims%j_start, pdims%j_end + do i=pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) .and. k >= ntop(i,j)+1 ) then + if ( k <= ntdsc(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + else if ( k <= ntml(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) end if - end do ! i - end do ! j + end if + end do ! i + !end do ! j end do ! k !$OMP end do @@ -1480,48 +1480,48 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - if ( test_well_mixed(i,j) ) then - if ( kwb0(i,j) == ntml(i,j) ) then - zwb0(i,j) = zh(i,j) - else if ( kwb0(i,j) == 2 ) then - zwb0(i,j) = z_uv(i,j,2) - else - k=kwb0(i,j) - ! now DB_GA_DRY(K) le 0 and DB_GA_DRY(K-1) gt 0 - ! so interpolate: - db_ratio = db_ga_dry_n(i,j,k-1) & - / ( db_ga_dry_n(i,j,k-1) - db_ga_dry_n(i,j,k) ) - db_ratio = max( zero, db_ratio ) ! trap for rounding error - zwb0(i,j)=z_uv(i,j,k-1) + & - db_ratio * (z_uv(i,j,k)-z_uv(i,j,k-1)) - end if - wb_surf_int(i,j) = bflux_surf(i,j) * z_tq(i,j,ksurf(i,j)) * & - ( one - z_tq(i,j,ksurf(i,j))/(2.0_r_bl*zwb0(i,j))) - wb_surf_int(i,j) = max( rbl_eps, wb_surf_int(i,j) ) +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + + if ( test_well_mixed(i,j) ) then + if ( kwb0(i,j) == ntml(i,j) ) then + zwb0(i,j) = zh(i,j) + else if ( kwb0(i,j) == 2 ) then + zwb0(i,j) = z_uv(i,j,2) else - ! only include surface layer contribution for unstable mixing - wb_surf_int(i,j) = rbl_eps + k=kwb0(i,j) + ! now DB_GA_DRY(K) le 0 and DB_GA_DRY(K-1) gt 0 + ! so interpolate: + db_ratio = db_ga_dry_n(i,j,k-1) & + / ( db_ga_dry_n(i,j,k-1) - db_ga_dry_n(i,j,k) ) + db_ratio = max( zero, db_ratio ) ! trap for rounding error + zwb0(i,j)=z_uv(i,j,k-1) + & + db_ratio * (z_uv(i,j,k)-z_uv(i,j,k-1)) end if + wb_surf_int(i,j) = bflux_surf(i,j) * z_tq(i,j,ksurf(i,j)) * & + ( one - z_tq(i,j,ksurf(i,j))/(2.0_r_bl*zwb0(i,j))) + wb_surf_int(i,j) = max( rbl_eps, wb_surf_int(i,j) ) + else + ! only include surface layer contribution for unstable mixing + wb_surf_int(i,j) = rbl_eps + end if - wbp_int(i,j) = wbp_int(i,j) + wb_surf_int(i,j) ! must be >0 + wbp_int(i,j) = wbp_int(i,j) + wb_surf_int(i,j) ! must be >0 - end do ! I -end do ! J +end do ! I +!end do ! J !$OMP end do if (model_type == mt_single_column) then !$OMP do SCHEDULE(STATIC) ! Save surface and bl-top layer integral for diagnostics - do j=pdims%j_start, pdims%j_end - do i=pdims%i_start, pdims%i_end - wbmix(i,j,ksurf(i,j)) = wb_surf_int(i,j) - wbend(i,j,ksurf(i,j)) = wb_surf_int(i,j) - end do ! i - end do ! j + !do j=pdims%j_start, pdims%j_end + do i=pdims%i_start, pdims%i_end + wbmix(i,j,ksurf(i,j)) = wb_surf_int(i,j) + wbend(i,j,ksurf(i,j)) = wb_surf_int(i,j) + end do ! i + !end do ! j !$OMP end do end if ! model_type @@ -1531,135 +1531,135 @@ subroutine excf_nl_9c ( & !----------------------------------------------------------------------- !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels-1 - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - - if ( test_well_mixed(i,j) ) then - ! ---------------------------------------------- - ! worth testing layer as well-mixed to cloud-top - ! ---------------------------------------------- - zb_ktop = 0.1_r_bl*z_inv(i,j) - zinv_pr(i,j) = z_inv(i,j) - zb_ktop - ! DB(K)is the K to K-1 difference and already - ! integrated up to KSURF, so start this loop at KSURF+1 - if ( (k >= ksurf(i,j)+1) .and. (k <= ntop(i,j)) ) then - - khtop(i,j) = zero - khsurf(i,j)= zero - f2 = zero - fsc = zero - - z_pr = z_uv(i,j,k) - zb_ktop - if (z_pr > zero .and. z_pr < zinv_pr(i,j)) then - z_ratio = z_pr/zinv_pr(i,j) - - if (flux_grad == LockWhelan2006) then - khtop(i,j) = 3.6_r_bl * vkman * rho_mix(i,j,k) *v_ktop(i,j) & - * zinv_pr(i,j) * (z_ratio**3) & - * ( (one-z_ratio)*(one-z_ratio) ) - f2 = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_ksum(i,j) > zero ) then - fsc = rho_mix(i,j,k) * 3.5_r_bl*(v_ktop(i,j)/v_ksum(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if - else - ! max to avoid rounding errors giving small negative numbers - khtop(i,j) = g1 * vkman * rho_mix(i,j,k) * v_ktop(i,j) & - * ((max(one - z_ratio,zero))**0.8_r_bl) & - * z_pr * z_ratio - end if - end if + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end - z_pr = z_uv(i,j,k) - if ( z_pr < z_inv(i,j)) then - !-------------------------------- - ! include surface-driven profile - !-------------------------------- - khsurf(i,j) = vkman * rho_mix(i,j,k) * & - w_h_top(i,j)*z_pr*( one - z_pr/z_inv(i,j) ) & - *( one - z_pr/z_inv(i,j) ) - end if + if ( test_well_mixed(i,j) ) then + ! ---------------------------------------------- + ! worth testing layer as well-mixed to cloud-top + ! ---------------------------------------------- + zb_ktop = 0.1_r_bl*z_inv(i,j) + zinv_pr(i,j) = z_inv(i,j) - zb_ktop + ! DB(K)is the K to K-1 difference and already + ! integrated up to KSURF, so start this loop at KSURF+1 + if ( (k >= ksurf(i,j)+1) .and. (k <= ntop(i,j)) ) then + + khtop(i,j) = zero + khsurf(i,j)= zero + f2 = zero + fsc = zero + + z_pr = z_uv(i,j,k) - zb_ktop + if (z_pr > zero .and. z_pr < zinv_pr(i,j)) then + z_ratio = z_pr/zinv_pr(i,j) if (flux_grad == LockWhelan2006) then - wslng = (f2+fsc)*tothf_zi(i,j) - ft_nt(i,j,k) - wqwng = (f2+fsc)*totqf_zi(i,j) - fq_nt(i,j,k) - end if - - if ( z_tq(i,j,k) <= z_cbase(i,j) ) then - ! Completely below cloud-base so use cloud-free formula - wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & - khtop(i,j) * db_noga_dry(i,j,k) - if (flux_grad == LockWhelan2006) then - wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & - ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) + khtop(i,j) = 3.6_r_bl * vkman * rho_mix(i,j,k) *v_ktop(i,j) & + * zinv_pr(i,j) * (z_ratio**3) & + * ( (one-z_ratio)*(one-z_ratio) ) + f2 = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_ksum(i,j) > zero ) then + fsc = rho_mix(i,j,k) * 3.5_r_bl*(v_ktop(i,j)/v_ksum(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - wb_cld = zero - else if (z_tq(i,j,k-1) >= z_cbase(i,j)) then - ! Completely above cloud-base so use cloudy formula - wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & - khtop(i,j) * db_noga_cld(i,j,k) ) - if (flux_grad == LockWhelan2006) then - wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & - ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & - btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & - ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & - bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) - end if - wb_scld = zero else - ! cloud-base within this integration range - ! so treat cloud and sub-cloud layer wb separately - wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & - khtop(i,j) * db_noga_dry(i,j,k) - wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & - khtop(i,j) * db_noga_cld(i,j,k) ) - if (flux_grad == LockWhelan2006) then - wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & - ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) - wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & - ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & - btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & - ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & - bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) - end if - cld_frac = (z_tq(i,j,k)-z_cbase(i,j)) & - /(z_tq(i,j,k)-z_tq(i,j,k-1)) - wb_cld = cld_frac * wb_cld - wb_scld = (one-cld_frac) * wb_scld + ! max to avoid rounding errors giving small negative numbers + khtop(i,j) = g1 * vkman * rho_mix(i,j,k) * v_ktop(i,j) & + * ((max(one - z_ratio,zero))**0.8_r_bl) & + * z_pr * z_ratio end if + end if - if ( dzrad_disc_opt == dzrad_1p5dz .and. k == ntop(i,j) ) then - ! At top of Sc layer, only include the part of the integral - ! that is below the radiatively-cooled cloud-top layer. - interp = ( z_inv(i,j) - dzrad(i,j) - z_tq(i,j,k-1) ) & - / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - wb_cld = wb_cld * interp - wb_scld = wb_scld * interp - end if + z_pr = z_uv(i,j,k) + if ( z_pr < z_inv(i,j)) then + !-------------------------------- + ! include surface-driven profile + !-------------------------------- + khsurf(i,j) = vkman * rho_mix(i,j,k) * & + w_h_top(i,j)*z_pr*( one - z_pr/z_inv(i,j) ) & + *( one - z_pr/z_inv(i,j) ) + end if - if (wb_cld >= zero) then - wbp_int(i,j) = wbp_int(i,j) + wb_cld - else - wbn_int(i,j) = wbn_int(i,j) - wb_cld + if (flux_grad == LockWhelan2006) then + wslng = (f2+fsc)*tothf_zi(i,j) - ft_nt(i,j,k) + wqwng = (f2+fsc)*totqf_zi(i,j) - fq_nt(i,j,k) + end if + + if ( z_tq(i,j,k) <= z_cbase(i,j) ) then + ! Completely below cloud-base so use cloud-free formula + wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & + khtop(i,j) * db_noga_dry(i,j,k) + if (flux_grad == LockWhelan2006) then + wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & + ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) end if - if (wb_scld >= zero) then - wbp_int(i,j) = wbp_int(i,j) + wb_scld - else - wbn_int(i,j) = wbn_int(i,j) - wb_scld + wb_cld = zero + else if (z_tq(i,j,k-1) >= z_cbase(i,j)) then + ! Completely above cloud-base so use cloudy formula + wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & + khtop(i,j) * db_noga_cld(i,j,k) ) + if (flux_grad == LockWhelan2006) then + wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & + ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & + btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & + ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & + bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) end if - - if (model_type == mt_single_column) then - wbmix(i,j,k) = wb_cld+wb_scld - wbend(i,j,k) = wb_cld+wb_scld + wb_scld = zero + else + ! cloud-base within this integration range + ! so treat cloud and sub-cloud layer wb separately + wb_scld = khsurf(i,j) * db_ga_dry_n(i,j,k) + & + khtop(i,j) * db_noga_dry(i,j,k) + wb_cld = ( khsurf(i,j) * db_ga_cld_n(i,j,k) + & + khtop(i,j) * db_noga_cld(i,j,k) ) + if (flux_grad == LockWhelan2006) then + wb_scld = wb_scld + ( g/rdz(i,j,k) ) * & + ( btm(i,j,k-1)*wslng + bqm(i,j,k-1)*wqwng ) + wb_cld = wb_cld + ( g/rdz(i,j,k) ) * ( & + ( btm(i,j,k-1)*(one-cf_ml(i,j)) + & + btm_cld(i,j,k-1)*cf_ml(i,j) )*wslng + & + ( bqm(i,j,k-1)*(one-cf_ml(i,j)) + & + bqm_cld(i,j,k-1)*cf_ml(i,j) )*wqwng ) end if + cld_frac = (z_tq(i,j,k)-z_cbase(i,j)) & + /(z_tq(i,j,k)-z_tq(i,j,k-1)) + wb_cld = cld_frac * wb_cld + wb_scld = (one-cld_frac) * wb_scld + end if - end if ! K + if ( dzrad_disc_opt == dzrad_1p5dz .and. k == ntop(i,j) ) then + ! At top of Sc layer, only include the part of the integral + ! that is below the radiatively-cooled cloud-top layer. + interp = ( z_inv(i,j) - dzrad(i,j) - z_tq(i,j,k-1) ) & + / ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + wb_cld = wb_cld * interp + wb_scld = wb_scld * interp + end if - end if ! TEST_WELL_MIXED - end do ! I + if (wb_cld >= zero) then + wbp_int(i,j) = wbp_int(i,j) + wb_cld + else + wbn_int(i,j) = wbn_int(i,j) - wb_cld + end if + if (wb_scld >= zero) then + wbp_int(i,j) = wbp_int(i,j) + wb_scld + else + wbn_int(i,j) = wbn_int(i,j) - wb_scld + end if + + if (model_type == mt_single_column) then + wbmix(i,j,k) = wb_cld+wb_scld + wbend(i,j,k) = wb_cld+wb_scld + end if + + end if ! K + + end if ! TEST_WELL_MIXED + !end do ! I end do ! J end do ! K end do @@ -1671,182 +1671,182 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (l_wtrac) iset_wtrac(i,j) = 0 - if ( test_well_mixed(i,j) ) then +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + if (l_wtrac) iset_wtrac(i,j) = 0 + if ( test_well_mixed(i,j) ) then + + wb_ratio(i,j) = wbn_int(i,j)/wbp_int(i,j) + + if ( wb_ratio(i,j) <= dec_thres(i,j) ) then + ! No need to test depth of mixing any further as within + ! well-mixed layer buoyancy flux integral criteria. + ! SML will simply stay well-mixed (and so use defaults) + ksurf_iterate(i,j)= .false. + ktop_iterate(i,j) = .false. + if ( dsc(i,j) ) then + ! Recouple DSC with SML: + dsc_removed(i,j) = 2 ! set indicator flag + ! move surface driven entrainment + ! RHOKH(z_i) = rho * w_e * DZL and w_e ~ 1/DB_TOP, so: + if ( db_top(i,j) > zero .and. db_dsct(i,j) > 0.01_r_bl ) then + ! can't calc Zil. term + rhokh_surf_ent(i,j) = rhokh_surf_ent(i,j) * & + ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & + rdz(i,j,ntdsc(i,j)+1) ) + rhokm_surf_ent(i,j) = rhokm_surf_ent(i,j) * & + ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & + rdz(i,j,ntdsc(i,j)+1) ) + if (.not. l_use_var_fixes) & + rhokm(i,j,ntdsc(i,j)+1) = rhokm_surf_ent(i,j) + end if + ! redesignate top-driven entrainment at ZHSC + ! (ignore that calculated at ZH) + rhokh_top_ent(i,j) = rhokh_dsct_ent(i,j) + rhokm_top_ent(i,j) = rhokm_dsct_ent(i,j) + zh(i,j) = zhsc(i,j) + ntml(i,j) = ntdsc(i,j) + v_top(i,j) = v_top_dsc(i,j) + zsml_base(i,j) = 0.1_r_bl * zh(i,j) + zc(i,j) = zc_dsc(i,j) + zhsc(i,j) = zero + ntdsc(i,j) = 0 + v_top_dsc(i,j) = zero + zdsc_base(i,j) = zero + zc_dsc(i,j) = zero + ft_nt_zh(i,j) = ft_nt_zhsc(i,j) + ft_nt_zhsc(i,j) = zero + fq_nt_zh(i,j) = fq_nt_zhsc(i,j) + fq_nt_zhsc(i,j) = zero + if (l_wtrac) iset_wtrac(i,j) = 1 + dsc(i,j) = .false. + dzh(i,j) = rmdi ! SML inversion subsumed into new mixed layer + cumulus(i,j) = .false. + coupled(i,j) = .false. + end if ! recoupled DSC layer + else ! buoyancy flux threshold violated - wb_ratio(i,j) = wbn_int(i,j)/wbp_int(i,j) + if ( l_use_sml_dsc_fixes ) then + ! Under fix l_use_sml_dsc_fixes: + ! If there is already a DSC layer diagnosed above the existing + ! ntml, we have so far just tested for recoupling of that DSC + ! with the SML; we have not tested whether the existing SML + ! should split off another decoupled layer. + ! In this case, we don't need to turn on ktop_iterate + ! (already setup earlier), and we should only turn on ksurf_iterate + ! if the SML is cumulus capped and using buoyancy flux integration + ! at cumulus points in general. - if ( wb_ratio(i,j) <= dec_thres(i,j) ) then - ! No need to test depth of mixing any further as within - ! well-mixed layer buoyancy flux integral criteria. - ! SML will simply stay well-mixed (and so use defaults) - ksurf_iterate(i,j)= .false. - ktop_iterate(i,j) = .false. if ( dsc(i,j) ) then - ! Recouple DSC with SML: - dsc_removed(i,j) = 2 ! set indicator flag - ! move surface driven entrainment - ! RHOKH(z_i) = rho * w_e * DZL and w_e ~ 1/DB_TOP, so: - if ( db_top(i,j) > zero .and. db_dsct(i,j) > 0.01_r_bl ) then - ! can't calc Zil. term - rhokh_surf_ent(i,j) = rhokh_surf_ent(i,j) * & - ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & - rdz(i,j,ntdsc(i,j)+1) ) - rhokm_surf_ent(i,j) = rhokm_surf_ent(i,j) * & - ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & - rdz(i,j,ntdsc(i,j)+1) ) - if (.not. l_use_var_fixes) & - rhokm(i,j,ntdsc(i,j)+1) = rhokm_surf_ent(i,j) + ! Existing DSC layer confirmed not to be merged with the SML. + if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low ) ) then + ! Cumulus layer under the DSC; only turn on ksurf_iterate if + ! using appropriate kprof_cu option. + ksurf_iterate(i,j) = .true. + dec_thres(i,j) = dec_thres_cu end if - ! redesignate top-driven entrainment at ZHSC - ! (ignore that calculated at ZH) - rhokh_top_ent(i,j) = rhokh_dsct_ent(i,j) - rhokm_top_ent(i,j) = rhokm_dsct_ent(i,j) - zh(i,j) = zhsc(i,j) - ntml(i,j) = ntdsc(i,j) - v_top(i,j) = v_top_dsc(i,j) - zsml_base(i,j) = 0.1_r_bl * zh(i,j) - zc(i,j) = zc_dsc(i,j) - zhsc(i,j) = zero - ntdsc(i,j) = 0 - v_top_dsc(i,j) = zero - zdsc_base(i,j) = zero - zc_dsc(i,j) = zero - ft_nt_zh(i,j) = ft_nt_zhsc(i,j) - ft_nt_zhsc(i,j) = zero - fq_nt_zh(i,j) = fq_nt_zhsc(i,j) - fq_nt_zhsc(i,j) = zero - if (l_wtrac) iset_wtrac(i,j) = 1 - dsc(i,j) = .false. - dzh(i,j) = rmdi ! SML inversion subsumed into new mixed layer - cumulus(i,j) = .false. - coupled(i,j) = .false. - end if ! recoupled DSC layer - else ! buoyancy flux threshold violated - - if ( l_use_sml_dsc_fixes ) then - ! Under fix l_use_sml_dsc_fixes: - ! If there is already a DSC layer diagnosed above the existing - ! ntml, we have so far just tested for recoupling of that DSC - ! with the SML; we have not tested whether the existing SML - ! should split off another decoupled layer. - ! In this case, we don't need to turn on ktop_iterate - ! (already setup earlier), and we should only turn on ksurf_iterate - ! if the SML is cumulus capped and using buoyancy flux integration - ! at cumulus points in general. - - if ( dsc(i,j) ) then - ! Existing DSC layer confirmed not to be merged with the SML. - if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low ) ) then - ! Cumulus layer under the DSC; only turn on ksurf_iterate if - ! using appropriate kprof_cu option. - ksurf_iterate(i,j) = .true. - dec_thres(i,j) = dec_thres_cu - end if - else - ! No existing DSC so we're splitting a new DSC off the SML; - ! turn on ktop_iterate to find the new DSC's base - ktop_iterate(i,j) = .true. - if ( (.not. cumulus(i,j)) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low ) then - ! Switch to using b-flux integration to find new lower SML-top, - ! but not in Cu layers unless using appropriate kprof_cu option. - ksurf_iterate(i,j) = .true. - if ( cumulus(i,j) ) dec_thres(i,j) = dec_thres_cu - end if + else + ! No existing DSC so we're splitting a new DSC off the SML; + ! turn on ktop_iterate to find the new DSC's base + ktop_iterate(i,j) = .true. + if ( (.not. cumulus(i,j)) .or. kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low ) then + ! Switch to using b-flux integration to find new lower SML-top, + ! but not in Cu layers unless using appropriate kprof_cu option. + ksurf_iterate(i,j) = .true. + if ( cumulus(i,j) ) dec_thres(i,j) = dec_thres_cu end if + end if - else ! ( .not. l_use_sml_dsc_fixes ) - ! If the bug-fix is off, we turn on ksurf_iterate at non-cumulus - ! points even if dsc was already true - ! (when all we've done is confirm the DSC shouldn't merge with the - ! SML, so we had no reason to be meddling with the SML). - ! This creates some inconsistencies down the line... - - !--------------------------------- - ! Extent of mixing must be reduced - !--------------------------------- - if ( .not. cumulus(i,j) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low) & - ksurf_iterate(i,j) = .true. - if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low) ) & - dec_thres(i,j) = dec_thres_cu - ktop_iterate(i,j) = .true. + else ! ( .not. l_use_sml_dsc_fixes ) + ! If the bug-fix is off, we turn on ksurf_iterate at non-cumulus + ! points even if dsc was already true + ! (when all we've done is confirm the DSC shouldn't merge with the + ! SML, so we had no reason to be meddling with the SML). + ! This creates some inconsistencies down the line... + + !--------------------------------- + ! Extent of mixing must be reduced + !--------------------------------- + if ( .not. cumulus(i,j) .or. kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low) & + ksurf_iterate(i,j) = .true. + if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low) ) & + dec_thres(i,j) = dec_thres_cu + ktop_iterate(i,j) = .true. - end if ! ( l_use_sml_dsc_fixes ) + end if ! ( l_use_sml_dsc_fixes ) - if ( .not. dsc(i,j) ) then - if (l_use_var_fixes .and. v_top(i,j) <= zero) then - ! If v_top=0 we have no mechanism to generate turbulence in a DSC - ktop_iterate(i,j) = .false. - z_top_lim(i,j) = zh(i,j) ! needs to be set here as zhsc not used - else - ! Set up a `COUPLED' decoupled layer, - ! implies no explicit `entrainment' at ZH. - ! Note a new ZH (and thence NTML) will be calculated by - ! wb integral iteration. - dsc(i,j) = .true. - coupled(i,j) = .true. - ntdsc(i,j) = ntml(i,j) - zhsc(i,j) = zh(i,j) - zc_dsc(i,j) = zc(i,j) - v_top_dsc(i,j) = v_top(i,j) - if ( l_use_sml_dsc_fixes ) then - ! Need to ensure SML-top entrainment flux is turned off when - ! using buoyancy-flux integration to find a new zh - v_top(i,j) = zero - ! Set coupling strength to maximal when just diagnosed decoupled - svl_diff_frac(i,j) = one - end if - v_sum_dsc(i,j) = v_sum(i,j) - ft_nt_zhsc(i,j) = ft_nt_zh(i,j) - fq_nt_zhsc(i,j) = fq_nt_zh(i,j) - if (l_wtrac) iset_wtrac(i,j) = 2 - ! put all entrainment into RHOKH_TOP - rhokh_dsct_ent(i,j) = rhokh_top_ent(i,j) & - + rhokh_surf_ent(i,j) - rhokh_top_ent(i,j) = zero - rhokh_surf_ent(i,j) = zero - rhokm_dsct_ent(i,j) = rhokm_top_ent(i,j) & - + rhokm_surf_ent(i,j) - rhokm_top_ent(i,j) = zero - rhokm_surf_ent(i,j) = zero - if (.not. l_use_var_fixes) then - rhokm_top(i,j,ntml(i,j)+1) = rhokm_top(i,j,ntml(i,j)+1) & - + rhokm(i,j,ntml(i,j)+1) - rhokm(i,j,ntml(i,j)+1) = zero - end if - end if ! test on l_use_var_fixes and v_top <= zero - end if ! test no not.dsc - end if ! test on WB_RATIO le DEC_THRES - end if ! testing for well-mixed layer (TEST_WELL_MIXED) + if ( .not. dsc(i,j) ) then + if (l_use_var_fixes .and. v_top(i,j) <= zero) then + ! If v_top=0 we have no mechanism to generate turbulence in a DSC + ktop_iterate(i,j) = .false. + z_top_lim(i,j) = zh(i,j) ! needs to be set here as zhsc not used + else + ! Set up a `COUPLED' decoupled layer, + ! implies no explicit `entrainment' at ZH. + ! Note a new ZH (and thence NTML) will be calculated by + ! wb integral iteration. + dsc(i,j) = .true. + coupled(i,j) = .true. + ntdsc(i,j) = ntml(i,j) + zhsc(i,j) = zh(i,j) + zc_dsc(i,j) = zc(i,j) + v_top_dsc(i,j) = v_top(i,j) + if ( l_use_sml_dsc_fixes ) then + ! Need to ensure SML-top entrainment flux is turned off when + ! using buoyancy-flux integration to find a new zh + v_top(i,j) = zero + ! Set coupling strength to maximal when just diagnosed decoupled + svl_diff_frac(i,j) = one + end if + v_sum_dsc(i,j) = v_sum(i,j) + ft_nt_zhsc(i,j) = ft_nt_zh(i,j) + fq_nt_zhsc(i,j) = fq_nt_zh(i,j) + if (l_wtrac) iset_wtrac(i,j) = 2 + ! put all entrainment into RHOKH_TOP + rhokh_dsct_ent(i,j) = rhokh_top_ent(i,j) & + + rhokh_surf_ent(i,j) + rhokh_top_ent(i,j) = zero + rhokh_surf_ent(i,j) = zero + rhokm_dsct_ent(i,j) = rhokm_top_ent(i,j) & + + rhokm_surf_ent(i,j) + rhokm_top_ent(i,j) = zero + rhokm_surf_ent(i,j) = zero + if (.not. l_use_var_fixes) then + rhokm_top(i,j,ntml(i,j)+1) = rhokm_top(i,j,ntml(i,j)+1) & + + rhokm(i,j,ntml(i,j)+1) + rhokm(i,j,ntml(i,j)+1) = zero + end if + end if ! test on l_use_var_fixes and v_top <= zero + end if ! test no not.dsc + end if ! test on WB_RATIO le DEC_THRES + end if ! testing for well-mixed layer (TEST_WELL_MIXED) - end do ! I -end do ! J +end do ! I +!end do ! J !$OMP end do ! Set water tracer fields in the same manner as water in last loop if (l_wtrac) then do i_wt = 1, n_wtrac !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (iset_wtrac(i,j) == 1) then - fq_nt_zh_wtrac(i,j,i_wt) = fq_nt_zhsc_wtrac(i,j,i_wt) - fq_nt_zhsc_wtrac(i,j,i_wt) = zero - else if (iset_wtrac(i,j) == 2) then - fq_nt_zhsc_wtrac(i,j,i_wt) = fq_nt_zh_wtrac(i,j,i_wt) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (iset_wtrac(i,j) == 1) then + fq_nt_zh_wtrac(i,j,i_wt) = fq_nt_zhsc_wtrac(i,j,i_wt) + fq_nt_zhsc_wtrac(i,j,i_wt) = zero + else if (iset_wtrac(i,j) == 2) then + fq_nt_zhsc_wtrac(i,j,i_wt) = fq_nt_zh_wtrac(i,j,i_wt) + end if end do + !end do !$OMP end do end do @@ -1865,57 +1865,57 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - if ( ksurf_iterate(i,j) ) then - !----------------------------------------------------------- - ! Mixing must extend just above surface layer - ! (not clear precisely how to define this here: for now use - ! KSURF calculated from ZH) - !----------------------------------------------------------- - z_bot_lim(i,j)=z_uv(i,j,ksurf(i,j)+1) & - + 0.1_r_bl * (z_uv(i,j,ksurf(i,j)+2)-z_uv(i,j,ksurf(i,j)+1)) - if ( l_converge_ga ) then - ! Allow K-surf to go up to the Sc-top (increasing zsml_top within - ! the cloud-top layer can still make the buoyancy fluxes below more - ! negative because this reduces the gradient adjustment) - z_top_lim(i,j)=max( z_top_lim(i,j), & - max( z_bot_lim(i,j), zhsc(i,j) ) ) - ! ksurf_iterate currently only true if test_well_mixed was true and - ! decoupling was diagnosed - if (cumulus(i,j)) then - ! ignore cloudy component of wb for cumulus sml iteration as the - ! cloudy flux is carried by the convection scheme - z_cbase(i,j) = z_top_lim(i,j) - else - ! sml has decoupled from the cloud layer above. For consistency with - ! the test_well_mixed calculation, continue to use saturated wb above - ! the now decoupled layer's cloud base. - z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - end if - else ! ( .not. l_converge_ga ) - ! limit K-surf to below cloud-top radiatively cooled layer - ! or original zh (held in z_top_lim) if no top-driven turbulence - z_top_lim(i,j)=max( z_top_lim(i,j), & - max( z_bot_lim(i,j), zhsc(i,j) - dzrad(i,j) ) ) - ! Always use saturated wb above decoupled-layer cloud-base +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + + if ( ksurf_iterate(i,j) ) then + !----------------------------------------------------------- + ! Mixing must extend just above surface layer + ! (not clear precisely how to define this here: for now use + ! KSURF calculated from ZH) + !----------------------------------------------------------- + z_bot_lim(i,j)=z_uv(i,j,ksurf(i,j)+1) & + + 0.1_r_bl * (z_uv(i,j,ksurf(i,j)+2)-z_uv(i,j,ksurf(i,j)+1)) + if ( l_converge_ga ) then + ! Allow K-surf to go up to the Sc-top (increasing zsml_top within + ! the cloud-top layer can still make the buoyancy fluxes below more + ! negative because this reduces the gradient adjustment) + z_top_lim(i,j)=max( z_top_lim(i,j), & + max( z_bot_lim(i,j), zhsc(i,j) ) ) + ! ksurf_iterate currently only true if test_well_mixed was true and + ! decoupling was diagnosed + if (cumulus(i,j)) then + ! ignore cloudy component of wb for cumulus sml iteration as the + ! cloudy flux is carried by the convection scheme + z_cbase(i,j) = z_top_lim(i,j) + else + ! sml has decoupled from the cloud layer above. For consistency with + ! the test_well_mixed calculation, continue to use saturated wb above + ! the now decoupled layer's cloud base. z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - end if ! ( l_converge_ga ) - !----------------------------------------------------- - ! Initial increment to ZSML_TOP found by dividing - ! up depth of layer within which it is allowed: - ! Start with ZSML_TOP at lower limit and work upwards - !----------------------------------------------------- - z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) - zsml_top(i,j) = z_bot_lim(i,j) + end if + else ! ( .not. l_converge_ga ) + ! limit K-surf to below cloud-top radiatively cooled layer + ! or original zh (held in z_top_lim) if no top-driven turbulence + z_top_lim(i,j)=max( z_top_lim(i,j), & + max( z_bot_lim(i,j), zhsc(i,j) - dzrad(i,j) ) ) + ! Always use saturated wb above decoupled-layer cloud-base + z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) + end if ! ( l_converge_ga ) + !----------------------------------------------------- + ! Initial increment to ZSML_TOP found by dividing + ! up depth of layer within which it is allowed: + ! Start with ZSML_TOP at lower limit and work upwards + !----------------------------------------------------- + z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) + zsml_top(i,j) = z_bot_lim(i,j) - wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES + wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES - end if ! KSURF_ITERATE + end if ! KSURF_ITERATE - end do end do +!end do !$OMP end do ! ---------------------------------------------------------------------- @@ -1927,53 +1927,53 @@ subroutine excf_nl_9c ( & if (kprof_cu == buoy_integ_low) lcl_fac = one_half !$OMP do SCHEDULE(DYNAMIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .and. .not. dsc(i,j) .and. & - .not. ksurf_iterate(i,j)) then - ! pure cumulus layer and necessary parameters not already set up - ksurf_iterate(i,j) = .true. - dec_thres(i,j) = dec_thres_cu ! use cu threshold - ! limit top of K_surf to be between... - z_bot_lim(i,j) = lcl_fac*z_lcl(i,j) ! ...a fraction of the LCL and... - z_top_lim(i,j) = z_lcl(i,j) + 1000.0_r_bl ! ...1km above the LCL - ntop(i,j)=ntml(i,j) - do while ( z_uv(i,j,ntop(i,j)+1) <= z_top_lim(i,j) .and. & - ntop(i,j)+1 < bl_levels-2 ) - ntop(i,j) = ntop(i,j) + 1 ! z_uv(ntop+1) > z_top_lim - end do - z_top_lim(i,j)=z_uv(i,j,ntop(i,j)) ! highest z_uv below z_top_lim - z_cbase(i,j) = z_top_lim(i,j) - !----------------------------------------------------- - ! Initial increment to ZSML_TOP found by dividing - ! up depth of layer within which it is allowed: - ! Start with ZSML_TOP at lower limit and work upwards - !----------------------------------------------------- - z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) & - / real(n_steps, r_bl) - zsml_top(i,j) = z_bot_lim(i,j) - wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES - end if - end do ! I - end do ! J + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .and. .not. dsc(i,j) .and. & + .not. ksurf_iterate(i,j)) then + ! pure cumulus layer and necessary parameters not already set up + ksurf_iterate(i,j) = .true. + dec_thres(i,j) = dec_thres_cu ! use cu threshold + ! limit top of K_surf to be between... + z_bot_lim(i,j) = lcl_fac*z_lcl(i,j) ! ...a fraction of the LCL and... + z_top_lim(i,j) = z_lcl(i,j) + 1000.0_r_bl ! ...1km above the LCL + ntop(i,j)=ntml(i,j) + do while ( z_uv(i,j,ntop(i,j)+1) <= z_top_lim(i,j) .and. & + ntop(i,j)+1 < bl_levels-2 ) + ntop(i,j) = ntop(i,j) + 1 ! z_uv(ntop+1) > z_top_lim + end do + z_top_lim(i,j)=z_uv(i,j,ntop(i,j)) ! highest z_uv below z_top_lim + z_cbase(i,j) = z_top_lim(i,j) + !----------------------------------------------------- + ! Initial increment to ZSML_TOP found by dividing + ! up depth of layer within which it is allowed: + ! Start with ZSML_TOP at lower limit and work upwards + !----------------------------------------------------- + z_inc(i,j)=(z_top_lim(i,j)-z_bot_lim(i,j)) & + / real(n_steps, r_bl) + zsml_top(i,j) = z_bot_lim(i,j) + wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES + end if + end do ! I + !end do ! J !$OMP end do end if ! test in kprof_cu ! ---------------------------------------------------------------------- !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) - ind_todo(l) = l - up(l) = 1 - if (ksurf_iterate(i,j)) then - to_do(l) = .true. - else - to_do(l) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) + ind_todo(l) = l + up(l) = 1 + if (ksurf_iterate(i,j)) then + to_do(l) = .true. + else + to_do(l) = .false. + end if end do +!end do !$OMP end do !$OMP MASTER @@ -2180,56 +2180,56 @@ subroutine excf_nl_9c ( & end do ! n_sweep !$OMP do SCHEDULE(DYNAMIC) -do jj = pdims%j_start, pdims%j_end,bl_segment_size - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - ntml_new(i,j) = 2 - status_ntml(i,j)=.true. - end do +do ii = pdims%i_start, pdims%i_end,bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + ntml_new(i,j) = 2 + status_ntml(i,j)=.true. + !end do end do do k = 2, bl_levels-2 !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - if ( ksurf_iterate(i,j) .and. status_ntml(i,j) ) then - ! ------------- - ! find new NTML - ! ------------- - if (z_uv(i,j,k+1) < zsml_top(i,j)) then - ntml_new(i,j) = k+1 + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + if ( ksurf_iterate(i,j) .and. status_ntml(i,j) ) then + ! ------------- + ! find new NTML + ! ------------- + if (z_uv(i,j,k+1) < zsml_top(i,j)) then + ntml_new(i,j) = k+1 + else + status_ntml(i,j)=.false. + end if + ! -------------------------------------------------------- + ! Rounding error previously found to give + ! ZSML_TOP > Z_TOP_LIM = ZHSC + ! Test on ZSML_TOP hitting thresholds consequently changed + ! but also include the following failsafe tests here. + ! -------------------------------------------------------- + if (dsc(i,j)) then + ntml(i,j) = min( ntdsc(i,j), ntml_new(i,j)-1 ) + zh(i,j) = min( zhsc(i,j), zsml_top(i,j) ) + else if (.not. status_ntml(i,j)) then + if (kprof_cu == buoy_integ_low) then + ! just use zh from buoyancy flux integration + ntml(i,j) = ntml_new(i,j)-1 + zh(i,j) = zsml_top(i,j) else - status_ntml(i,j)=.false. - end if - ! -------------------------------------------------------- - ! Rounding error previously found to give - ! ZSML_TOP > Z_TOP_LIM = ZHSC - ! Test on ZSML_TOP hitting thresholds consequently changed - ! but also include the following failsafe tests here. - ! -------------------------------------------------------- - if (dsc(i,j)) then - ntml(i,j) = min( ntdsc(i,j), ntml_new(i,j)-1 ) - zh(i,j) = min( zhsc(i,j), zsml_top(i,j) ) - else if (.not. status_ntml(i,j)) then - if (kprof_cu == buoy_integ_low) then - ! just use zh from buoyancy flux integration - ntml(i,j) = ntml_new(i,j)-1 - zh(i,j) = zsml_top(i,j) + ! for kprof_cu eq buoy_integ use max of ksurf top calculated + ! from buoyancy flux integration and original LCL based values + ntml(i,j) = max( ntml(i,j), ntml_new(i,j)-1 ) + if (l_use_var_fixes) then + zh(i,j) = max( zh(i,j), zsml_top(i,j) ) else - ! for kprof_cu eq buoy_integ use max of ksurf top calculated - ! from buoyancy flux integration and original LCL based values - ntml(i,j) = max( ntml(i,j), ntml_new(i,j)-1 ) - if (l_use_var_fixes) then - zh(i,j) = max( zh(i,j), zsml_top(i,j) ) - else - zh(i,j) = zsml_top(i,j) - end if + zh(i,j) = zsml_top(i,j) end if end if + end if - end if ! KSURF_ITERATE true + end if ! KSURF_ITERATE true - end do + !end do end do end do end do @@ -2241,164 +2241,166 @@ subroutine excf_nl_9c ( & ! ---------------------------------------------------------------------- !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - if ( ktop_iterate(i,j) ) then - - ! Lower limit on base of DSC layer - z_bot_lim(i,j) = 0.1_r_bl * zh(i,j) - ! Upper limit on base of DSC layer - z_top_lim(i,j) = zhsc(i,j) - dzrad(i,j) - ! If cumulus limit base of top-driven mixing to above ZH - ! (not applying this limit if converging the gradient adjustment, - ! as in this case we need to allow kh_surf to go as high as it wants - ! so the DSC kh needs to be able to overlap with it). - if ( cumulus(i,j) .and. (.not. l_converge_ga) ) then - z_bot_lim(i,j) = min( zh(i,j), z_top_lim(i,j) ) - end if - - z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + + if ( ktop_iterate(i,j) ) then + + ! Lower limit on base of DSC layer + z_bot_lim(i,j) = 0.1_r_bl * zh(i,j) + ! Upper limit on base of DSC layer + z_top_lim(i,j) = zhsc(i,j) - dzrad(i,j) + ! If cumulus limit base of top-driven mixing to above ZH + ! (not applying this limit if converging the gradient adjustment, + ! as in this case we need to allow kh_surf to go as high as it wants + ! so the DSC kh needs to be able to overlap with it). + if ( cumulus(i,j) .and. (.not. l_converge_ga) ) then + z_bot_lim(i,j) = min( zh(i,j), z_top_lim(i,j) ) + end if - !..Divide up depth of layer within which ZDSC_BASE is allowed - z_inc(i,j) = (z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) - zdsc_base(i,j) = z_bot_lim(i,j) - ! will start at Z_BOT_LIM+Z_INC + z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - if (l_reset_dec_thres) then - dec_thres(i,j) = dec_thres_cloud ! reset after potential use - ! of dec_thres_cu or - ! dec_thres_clear above - end if - wb_ratio(i,j) = dec_thres(i,j) + one ! to be > DEC_THRES + !..Divide up depth of layer within which ZDSC_BASE is allowed + z_inc(i,j) = (z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) + zdsc_base(i,j) = z_bot_lim(i,j) + ! will start at Z_BOT_LIM+Z_INC - end if ! KTOP_ITERATE + if (l_reset_dec_thres) then + dec_thres(i,j) = dec_thres_cloud ! reset after potential use + ! of dec_thres_cu or + ! dec_thres_clear above + end if + wb_ratio(i,j) = dec_thres(i,j) + one ! to be > DEC_THRES - end do + end if ! KTOP_ITERATE - do i = pdims%i_start, pdims%i_end - if ( ktop_iterate(i,j) .and. wb_dzrad_int(i,j) < zero ) then - !------------------------------------------------------------- - ! Estimation of wb integral over radiatively cooled cloud-top - ! region not yet performed (ie. DSC over stable surface) so - ! do it now. - !------------------------------------------------------------- - z_inv(i,j) = zhsc(i,j) - z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) - cf_ml(i,j) = cf_dsc(i,j) - df_ctop(i,j)= df_dsct_over_cp(i,j) - rho_we(i,j) = rho_we_dsc(i,j) - dsl(i,j) = dsl_dsc(i,j) - dqw(i,j) = dqw_dsc(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) - ntop(i,j) = ntdsc(i,j) - dzrad(i,j) = 100.0_r_bl - if ( dzrad_disc_opt == dzrad_ntm1 ) then - ! Original discretisation of cloud-top radiatively-cooled layer; - ! set the base of the layer at theta-level ntdsc-1 - ! (which is between 1.5 and 2.5 model-levels below cloud-top, - ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), - ! then lower it by extra whole model-levels if needed to increase the - ! depth to 100m. +end do +!$OMP end do - ntop(i,j) = ntop(i,j) - 1 - if ( ntop(i,j) > 1 ) then - inner_loop2: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) - ntop(i,j) = ntop(i,j) - 1 - if (ntop(i,j) == 1 ) then - exit inner_loop2 - end if - end do inner_loop2 - end if - dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) - - else if ( dzrad_disc_opt == dzrad_1p5dz ) then - ! Alternative method; - ! set the base of the layer 1.5 model-levels below cloud-top so that - ! it always moves smoothly with zhsc instead of jumping suddenly - ! (1.5 levels is the minimum distance below zhsc which will always - ! keep the base of the layer below theta-level ntdsc, which is a - ! requirement of the discretisation). - ! The depth is also forced to be at least 100m. - - ! Find smoothly-varying height 1.5 model-levels below z_inv - k_inv = ntop(i,j) - if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then - interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) - z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) - else - interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) - z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) - end if - ! Go down to 100m below z_inv (or cloud-base) if that is lower, - ! but don't allow lower than theta-level 1 - z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) - z_pr = max( z_pr, z_tq(i,j,1) ) - ! Set ntop to the rho-level straddling the base of the layer - do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) +!$OMP do SCHEDULE(DYNAMIC) +do i = pdims%i_start, pdims%i_end + if ( ktop_iterate(i,j) .and. wb_dzrad_int(i,j) < zero ) then + !------------------------------------------------------------- + ! Estimation of wb integral over radiatively cooled cloud-top + ! region not yet performed (ie. DSC over stable surface) so + ! do it now. + !------------------------------------------------------------- + z_inv(i,j) = zhsc(i,j) + z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) + cf_ml(i,j) = cf_dsc(i,j) + df_ctop(i,j)= df_dsct_over_cp(i,j) + rho_we(i,j) = rho_we_dsc(i,j) + dsl(i,j) = dsl_dsc(i,j) + dqw(i,j) = dqw_dsc(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) + ntop(i,j) = ntdsc(i,j) + dzrad(i,j) = 100.0_r_bl + if ( dzrad_disc_opt == dzrad_ntm1 ) then + ! Original discretisation of cloud-top radiatively-cooled layer; + ! set the base of the layer at theta-level ntdsc-1 + ! (which is between 1.5 and 2.5 model-levels below cloud-top, + ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), + ! then lower it by extra whole model-levels if needed to increase the + ! depth to 100m. + + ntop(i,j) = ntop(i,j) - 1 + if ( ntop(i,j) > 1 ) then + inner_loop2: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) ntop(i,j) = ntop(i,j) - 1 - end do - ! Set new cloud-top radiatively-cooled layer-depth - dzrad(i,j) = z_inv(i,j) - z_pr - - ! Fraction of rho-level ntop below the dzrad layer - interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) - ! Add on SL and qw differences between the discontinuous inversion base - ! and the base of the dzrad layer - dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & - - interp *sl(i,j,ntop(i,j)) - dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & - - interp *qw(i,j,ntop(i,j)) - - end if ! ( dzrad_disc_opt ) - - wsl_dzrad_int = dzrad(i,j) * & - ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) - wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) - - wb_dzrad_int(i,j) = wsl_dzrad_int * ( & - (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & - + wqw_dzrad_int * ( & - (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) - wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) - wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) - - if (model_type == mt_single_column) then - do k = ntop(i,j)+1, ntdsc(i,j)+1 - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end do + if (ntop(i,j) == 1 ) then + exit inner_loop2 + end if + end do inner_loop2 end if - + dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) + + else if ( dzrad_disc_opt == dzrad_1p5dz ) then + ! Alternative method; + ! set the base of the layer 1.5 model-levels below cloud-top so that + ! it always moves smoothly with zhsc instead of jumping suddenly + ! (1.5 levels is the minimum distance below zhsc which will always + ! keep the base of the layer below theta-level ntdsc, which is a + ! requirement of the discretisation). + ! The depth is also forced to be at least 100m. + + ! Find smoothly-varying height 1.5 model-levels below z_inv + k_inv = ntop(i,j) + if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then + interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + z_pr = (one-interp) * z_tq(i,j,k_inv-1) & + + interp * z_uv(i,j,k_inv) + else + interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + z_pr = (one-interp) * z_uv(i,j,k_inv) & + + interp * z_tq(i,j,k_inv) + end if + ! Go down to 100m below z_inv (or cloud-base) if that is lower, + ! but don't allow lower than theta-level 1 + z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) + z_pr = max( z_pr, z_tq(i,j,1) ) + ! Set ntop to the rho-level straddling the base of the layer + do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) + ntop(i,j) = ntop(i,j) - 1 + end do + ! Set new cloud-top radiatively-cooled layer-depth + dzrad(i,j) = z_inv(i,j) - z_pr + + ! Fraction of rho-level ntop below the dzrad layer + interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + ! Add on SL and qw differences between the discontinuous inversion base + ! and the base of the dzrad layer + dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & + - interp *sl(i,j,ntop(i,j)) + dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & + - interp *qw(i,j,ntop(i,j)) + + end if ! ( dzrad_disc_opt ) + + wsl_dzrad_int = dzrad(i,j) * & + ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) + wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) + + wb_dzrad_int(i,j) = wsl_dzrad_int * ( & + (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & + + wqw_dzrad_int * ( & + (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) + wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) + wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) + + if (model_type == mt_single_column) then + do k = ntop(i,j)+1, ntdsc(i,j)+1 + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end do end if - wb_dzrad_int(i,j) = max( rbl_eps, wb_dzrad_int(i,j) ) + end if + + wb_dzrad_int(i,j) = max( rbl_eps, wb_dzrad_int(i,j) ) - end do ! I -end do ! J +end do ! I +!end do ! J !$OMP end do !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) - ind_todo(l) = l - up(l) = 1 - if (ktop_iterate(i,j)) then - to_do(l) = .true. - else - to_do(l) = .false. - end if - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + l = i - pdims%i_start + 1 + pdims%i_len * (j - pdims%j_start) + ind_todo(l) = l + up(l) = 1 + if (ktop_iterate(i,j)) then + to_do(l) = .true. + else + to_do(l) = .false. + end if end do +!end do !$OMP end do !$OMP MASTER @@ -2679,60 +2681,62 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels-1 - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! Include sml - if ( .not. zdsc_base(i,j) < zsml_top(i,j) ) then - ! Note: wbend already includes the SML buoyancy flux - ! in the case where the SML and DSC mixing profiles overlap - ! (since the DSC b-flux integration sums both the SML and DSC - ! mixing contributions and stores this in wbend). - ! So only add the SML contribution back on if they do not overlap. - wbend(i,j,k) = wbend(i,j,k) + wbend_sml(i,j,k) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! Include sml + if ( .not. zdsc_base(i,j) < zsml_top(i,j) ) then + ! Note: wbend already includes the SML buoyancy flux + ! in the case where the SML and DSC mixing profiles overlap + ! (since the DSC b-flux integration sums both the SML and DSC + ! mixing contributions and stores this in wbend). + ! So only add the SML contribution back on if they do not overlap. + wbend(i,j,k) = wbend(i,j,k) + wbend_sml(i,j,k) + end if end do + !end do end do !$OMP end do -!$OMP do SCHEDULE(STATIC) + ! Note parallelised over j as k isn't independent - do j = pdims%j_start, pdims%j_end - do k = 1, bl_levels-1 - do i = pdims%i_start, pdims%i_end - ! convert to m2/s-3 - if ( k <= ntop(i,j) ) then - if ( k <= ksurf(i,j) ) then - gamma_wbs = ( (wbmix(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & - - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) - wbmix(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) - - gamma_wbs = ( (wbend(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & - - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) - wbend(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) - else - wbmix(i,j,k)=wbmix(i,j,k)*rdz(i,j,k) - wbend(i,j,k)=wbend(i,j,k)*rdz(i,j,k) - end if +!do j = pdims%j_start, pdims%j_end + do k = 1, bl_levels-1 + !$OMP do SCHEDULE(STATIC) + do i = pdims%i_start, pdims%i_end + ! convert to m2/s-3 + if ( k <= ntop(i,j) ) then + if ( k <= ksurf(i,j) ) then + gamma_wbs = ( (wbmix(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & + - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) + wbmix(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) + + gamma_wbs = ( (wbend(i,j,ksurf(i,j))/z_tq(i,j,ksurf(i,j))) & + - bflux_surf(i,j) )*2.0_r_bl/z_tq(i,j,ksurf(i,j)) + wbend(i,j,k) = bflux_surf(i,j) + gamma_wbs*z_uv(i,j,k) + else + wbmix(i,j,k)=wbmix(i,j,k)*rdz(i,j,k) + wbend(i,j,k)=wbend(i,j,k)*rdz(i,j,k) end if - end do + end if end do + !$OMP end do end do -!$OMP end do + !end do + !$OMP do SCHEDULE(STATIC) do k = 1, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - ! convert to m2/s-3 - if ( k >= ntop(i,j)+1 .and. k <= ntdsc(i,j)+1 ) then - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if - if ( k >= ntop(i,j)+1 .and. k <= ntml(i,j)+1 ) then - wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end if - end do + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + ! convert to m2/s-3 + if ( k >= ntop(i,j)+1 .and. k <= ntdsc(i,j)+1 ) then + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end if + if ( k >= ntop(i,j)+1 .and. k <= ntml(i,j)+1 ) then + wbmix(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end if end do + !end do end do !$OMP end do end if ! model_type @@ -2748,74 +2752,74 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .or. & - ( dsc(i,j) .and. zdsc_base(i,j) < zh(i,j) ) ) then - ! ignore SML `cloud-top' driven mixing - zsml_base(i,j) = zh(i,j) - v_top(i,j) = zero - else - zsml_base(i,j) = 0.1_r_bl*zh(i,j) - end if - end do ! loop over j - end do ! loop over I + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .or. & + ( dsc(i,j) .and. zdsc_base(i,j) < zh(i,j) ) ) then + ! ignore SML `cloud-top' driven mixing + zsml_base(i,j) = zh(i,j) + v_top(i,j) = zero + else + zsml_base(i,j) = 0.1_r_bl*zh(i,j) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do else ! entr_smooth_dec off !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) .or. coupled(i,j) ) then - ! ignore SML `cloud-top' driven mixing - zsml_base(i,j) = zh(i,j) - v_top(i,j) = zero - else - zsml_base(i,j) = 0.1_r_bl*zh(i,j) - end if - end do ! loop over j - end do ! loop over I + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) .or. coupled(i,j) ) then + ! ignore SML `cloud-top' driven mixing + zsml_base(i,j) = zh(i,j) + v_top(i,j) = zero + else + zsml_base(i,j) = 0.1_r_bl*zh(i,j) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do end if ! test on entr_smooth_dec !cdir collapse !$OMP do SCHEDULE(STATIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - zsml_top(i,j) = zh(i,j) - if (bl_res_inv /= off .and. dzh(i,j) > zero) & - zsml_top(i,j) = zh(i,j)+dzh(i,j) - end do ! loop over j -end do ! loop over I +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + zsml_top(i,j) = zh(i,j) + if (bl_res_inv /= off .and. dzh(i,j) > zero) & + zsml_top(i,j) = zh(i,j)+dzh(i,j) +end do ! loop over j +!end do ! loop over I !$OMP end do if ( kprof_cu == klcl_entr ) then !cdir collapse !$OMP do SCHEDULE(STATIC) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if ( cumulus(i,j) ) then - zsml_top(i,j) = zh(i,j)+max_cu_depth - ! Keep top of Kprof for cu below DSC base: - if (dsc(i,j)) then - zsml_top(i,j) = min( zsml_top(i,j), zdsc_base(i,j) ) - end if - ! ...but at least 1.1*z_lcl: - zsml_top(i,j) = max( 1.1_r_bl*z_lcl(i,j), zsml_top(i,j) ) - ! ...but no higher than parcel top: - zsml_top(i,j) = min( zhpar(i,j), zsml_top(i,j) ) + !do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if ( cumulus(i,j) ) then + zsml_top(i,j) = zh(i,j)+max_cu_depth + ! Keep top of Kprof for cu below DSC base: + if (dsc(i,j)) then + zsml_top(i,j) = min( zsml_top(i,j), zdsc_base(i,j) ) + end if + ! ...but at least 1.1*z_lcl: + zsml_top(i,j) = max( 1.1_r_bl*z_lcl(i,j), zsml_top(i,j) ) + ! ...but no higher than parcel top: + zsml_top(i,j) = min( zhpar(i,j), zsml_top(i,j) ) - cu_depth_scale(i,j) = (zsml_top(i,j)-zh(i,j))/3.0_r_bl + cu_depth_scale(i,j) = (zsml_top(i,j)-zh(i,j))/3.0_r_bl - rhokh_lcl(i,j) = zero - ! Use the BL entrainment parametrization as calculated above - rhokh_lcl(i,j) = min( rhokh_surf_ent(i,j), 5.0_r_bl) - end if - end do ! loop over j - end do ! loop over I + rhokh_lcl(i,j) = zero + ! Use the BL entrainment parametrization as calculated above + rhokh_lcl(i,j) = min( rhokh_surf_ent(i,j), 5.0_r_bl) + end if + end do ! loop over j + !end do ! loop over I !$OMP end do end if ! test on kprof_cu !----------------------------------------------------------------------- @@ -2826,78 +2830,78 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - k=ntml(i,j)+1 - kh_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( & - rhokh_surf_ent(i,j) / & - ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) +! do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + k=ntml(i,j)+1 + kh_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( & + rhokh_surf_ent(i,j) / & + ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) + if (l_use_var_fixes) then + km_top_factor(i,j) = max( 0.7_r_bl, one-sqrt( rhokm_surf_ent(i,j) / & + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + else + km_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( rhokm(i,j,k) / & + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + end if + scdepth(i,j) = zh(i,j) - zsml_base(i,j) + factor = g1 * rho_mix(i,j,k) * v_top(i,j) *vkman *scdepth(i,j) + if ( factor > zero) then + kh_sct_factor(i,j) = one - ( rhokh_top_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 + else + kh_sct_factor(i,j) = one + end if + factor = g1 * rho_wet_tq(i,j,k-1) * v_top(i,j) * & + vkman * scdepth(i,j) * 0.75_r_bl + if ( factor > zero) then if (l_use_var_fixes) then - km_top_factor(i,j) = max( 0.7_r_bl, one-sqrt( rhokm_surf_ent(i,j) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + km_sct_factor(i,j) = one - & + ( rhokm_top_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 else - km_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( rhokm(i,j,k) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) + km_sct_factor(i,j) = one - & + ( rhokm_top(i,j,k) / factor )**1.25_r_bl end if - scdepth(i,j) = zh(i,j) - zsml_base(i,j) - factor = g1 * rho_mix(i,j,k) * v_top(i,j) *vkman *scdepth(i,j) + else + km_sct_factor(i,j) = one + end if + + if (ntdsc(i,j) > 0) then + !------------------------------------------------------------- + ! Set up factors to ensure K profile continuity at ZHSC; + ! no need to limit size of factor as precise shape of top-down + ! mixing profile not important. + ! Only calculate _DSCT_FACTORs when a decoupled stratocumulus + ! layer exists, i.e. NTDSC > 0. + !------------------------------------------------------------- + k=ntdsc(i,j)+1 + dscdepth(i,j) = zhsc(i,j) - zdsc_base(i,j) + factor = g1*rho_mix(i,j,k)*v_top_dsc(i,j)*vkman*dscdepth(i,j) if ( factor > zero) then - kh_sct_factor(i,j) = one - ( rhokh_top_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 + kh_dsct_factor(i,j) = one - & + ( rhokh_dsct_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 else - kh_sct_factor(i,j) = one + kh_dsct_factor(i,j) = one end if - factor = g1 * rho_wet_tq(i,j,k-1) * v_top(i,j) * & - vkman * scdepth(i,j) * 0.75_r_bl + + factor = 0.75_r_bl * g1 * rho_wet_tq(i,j,k-1) * v_top_dsc(i,j) * & + vkman * dscdepth(i,j) if ( factor > zero) then if (l_use_var_fixes) then - km_sct_factor(i,j) = one - & - ( rhokm_top_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 + km_dsct_factor(i,j) = one - & + ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 else - km_sct_factor(i,j) = one - & - ( rhokm_top(i,j,k) / factor )**1.25_r_bl + km_dsct_factor(i,j) = one - & + ( rhokm_top(i,j,k) / factor )**1.25_r_bl end if else - km_sct_factor(i,j) = one - end if - - if (ntdsc(i,j) > 0) then - !------------------------------------------------------------- - ! Set up factors to ensure K profile continuity at ZHSC; - ! no need to limit size of factor as precise shape of top-down - ! mixing profile not important. - ! Only calculate _DSCT_FACTORs when a decoupled stratocumulus - ! layer exists, i.e. NTDSC > 0. - !------------------------------------------------------------- - k=ntdsc(i,j)+1 - dscdepth(i,j) = zhsc(i,j) - zdsc_base(i,j) - factor = g1*rho_mix(i,j,k)*v_top_dsc(i,j)*vkman*dscdepth(i,j) - if ( factor > zero) then - kh_dsct_factor(i,j) = one - & - ( rhokh_dsct_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 - else - kh_dsct_factor(i,j) = one - end if - - factor = 0.75_r_bl * g1 * rho_wet_tq(i,j,k-1) * v_top_dsc(i,j) * & - vkman * dscdepth(i,j) - if ( factor > zero) then - if (l_use_var_fixes) then - km_dsct_factor(i,j) = one - & - ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 - else - km_dsct_factor(i,j) = one - & - ( rhokm_top(i,j,k) / factor )**1.25_r_bl - end if - else - km_dsct_factor(i,j) = one - end if + km_dsct_factor(i,j) = one end if - end do + end if end do +! end do !$OMP end do !----------------------------------------------------------------------- @@ -2910,12 +2914,12 @@ subroutine excf_nl_9c ( & !$OMP MASTER !cdir collapse -do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - scbase(i,j) = .false. - nbdsc(i,j) = 0 - end do +!do j = pdims%j_start, pdims%j_end +do i = pdims%i_start, pdims%i_end + scbase(i,j) = .false. + nbdsc(i,j) = 0 end do +!end do !$OMP end MASTER !$OMP BARRIER @@ -2932,150 +2936,150 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*g1) !$OMP do SCHEDULE(STATIC) -do jj = pdims%j_start, pdims%j_end, bl_segment_size +do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - - ! Calculate the height of u,v-level above the surface - ! *APL: z0m removed from z in K(z) - zk_uv = z_uv(i,j,k) - - ! Calculate the height of T,q-level above the surface - - zk_tq = z_tq(i,j,k-1) - - if ( zk_uv < zh(i,j) .and. & - zk_uv > zsml_base(i,j) ) then - z_pr = zk_uv - zsml_base(i,j) - zh_pr = zh(i,j) - zsml_base(i,j) - z_ratio = z_pr/zh_pr - if (flux_grad == LockWhelan2006) then - - rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) * v_top(i,j) & - * zh_pr * (z_ratio**3) & - * (( one - z_ratio )**2) - - if ( .not. coupled(i,j) ) then - rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_sum(i,j) > zero ) then - rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & - * (v_top(i,j)/v_sum(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + + ! Calculate the height of u,v-level above the surface + ! *APL: z0m removed from z in K(z) + zk_uv = z_uv(i,j,k) + + ! Calculate the height of T,q-level above the surface + + zk_tq = z_tq(i,j,k-1) + + if ( zk_uv < zh(i,j) .and. & + zk_uv > zsml_base(i,j) ) then + z_pr = zk_uv - zsml_base(i,j) + zh_pr = zh(i,j) - zsml_base(i,j) + z_ratio = z_pr/zh_pr + if (flux_grad == LockWhelan2006) then + + rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) * v_top(i,j) & + * zh_pr * (z_ratio**3) & + * (( one - z_ratio )**2) + + if ( .not. coupled(i,j) ) then + rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_sum(i,j) > zero ) then + rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & + * (v_top(i,j)/v_sum(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - - else ! Not LockWhelan2006 - ! max to avoid rounding errors giving small negative numbers - rhokh_top(i,j,k) = rho_mix(i,j,k) * v_top(i,j) * g1 * & - vkman * ((max(one - kh_sct_factor(i,j)*z_ratio,zero))**0.8_r_bl )& - * z_pr * z_ratio end if + else ! Not LockWhelan2006 + ! max to avoid rounding errors giving small negative numbers + rhokh_top(i,j,k) = rho_mix(i,j,k) * v_top(i,j) * g1 * & + vkman * ((max(one - kh_sct_factor(i,j)*z_ratio,zero))**0.8_r_bl )& + * z_pr * z_ratio end if - !------------------------------------------------------- - ! For LockWhelan2006, KM_TOP could be changed to match - ! the shape of KH_TOP. This has not been done on the - ! grounds that the change in shape arises with the - ! inclusion of the other non-gradient terms. - !------------------------------------------------------- - if ( zk_tq < zh(i,j) .and. & - zk_tq > zsml_base(i,j) ) then - z_pr = zk_tq - zsml_base(i,j) - zh_pr = zh(i,j) - zsml_base(i,j) - rhokm_top(i,j,k) = 0.75_r_bl * rho_wet_tq(i,j,k-1) * v_top(i,j) * & - g1 * vkman * & - ( (max(one - km_sct_factor(i,j)*z_pr/zh_pr, zero))**0.8_r_bl ) & - * z_pr * z_pr / zh_pr - ! PRANDTL=0.75 - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,k)=rhokm_top(i,j,k)*c_tke*v_top(i,j)/zh(i,j) - end if + + end if + !------------------------------------------------------- + ! For LockWhelan2006, KM_TOP could be changed to match + ! the shape of KH_TOP. This has not been done on the + ! grounds that the change in shape arises with the + ! inclusion of the other non-gradient terms. + !------------------------------------------------------- + if ( zk_tq < zh(i,j) .and. & + zk_tq > zsml_base(i,j) ) then + z_pr = zk_tq - zsml_base(i,j) + zh_pr = zh(i,j) - zsml_base(i,j) + rhokm_top(i,j,k) = 0.75_r_bl * rho_wet_tq(i,j,k-1) * v_top(i,j) * & + g1 * vkman * & + ( (max(one - km_sct_factor(i,j)*z_pr/zh_pr, zero))**0.8_r_bl ) & + * z_pr * z_pr / zh_pr + ! PRANDTL=0.75 + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,k)=rhokm_top(i,j,k)*c_tke*v_top(i,j)/zh(i,j) + end if + end if + !------------------------------------------------------------- + ! Add contribution to top-down mixing coefficient + ! profiles for decoupled stratocumulus layers when + ! one exists + !------------------------------------------------------------- + if ( zk_uv < zhsc(i,j) .and. & + zk_uv > zdsc_base(i,j) ) then + if (.not. scbase(i,j) ) then + scbase(i,j) = .true. + ! identifies lowest layer below which there is mixing + nbdsc(i,j) = k end if - !------------------------------------------------------------- - ! Add contribution to top-down mixing coefficient - ! profiles for decoupled stratocumulus layers when - ! one exists - !------------------------------------------------------------- - if ( zk_uv < zhsc(i,j) .and. & - zk_uv > zdsc_base(i,j) ) then - if (.not. scbase(i,j) ) then - scbase(i,j) = .true. - ! identifies lowest layer below which there is mixing - nbdsc(i,j) = k + !----------------------------------------------------------- + ! Calculate RHOK(H/M)_TOP, top-down turbulent mixing + ! profiles and add to any generated in the surface mixing + ! layer. + ! This is a variation on an up-side-down version of the + ! cubic surface-forced profiles above. Implement between + ! at least the top of the `surface layer' (at Z=0.1*ZH) and + ! ZHSC. + !----------------------------------------------------------- + z_pr = zk_uv - zdsc_base(i,j) + zh_pr = zhsc(i,j) - zdsc_base(i,j) + z_ratio = z_pr/zh_pr + + if (flux_grad == LockWhelan2006) then + + rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) & + * v_top_dsc(i,j) * zh_pr * (z_ratio**3) & + * (( one - z_ratio )**2) + + rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & + * 2.0_r_bl**( z_ratio**4 ) + if ( v_sum_dsc(i,j) > zero ) then + rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & + * (v_top_dsc(i,j)/v_sum_dsc(i,j)) & + * (z_ratio**3) * (one-z_ratio) end if - !----------------------------------------------------------- - ! Calculate RHOK(H/M)_TOP, top-down turbulent mixing - ! profiles and add to any generated in the surface mixing - ! layer. - ! This is a variation on an up-side-down version of the - ! cubic surface-forced profiles above. Implement between - ! at least the top of the `surface layer' (at Z=0.1*ZH) and - ! ZHSC. - !----------------------------------------------------------- - z_pr = zk_uv - zdsc_base(i,j) - zh_pr = zhsc(i,j) - zdsc_base(i,j) - z_ratio = z_pr/zh_pr - if (flux_grad == LockWhelan2006) then - - rhokh_top(i,j,k) = 3.6_r_bl*vkman * rho_mix(i,j,k) & - * v_top_dsc(i,j) * zh_pr * (z_ratio**3) & - * (( one - z_ratio )**2) - - rhof2(i,j,k) = rho_mix(i,j,k) * one_half * z_ratio & - * 2.0_r_bl**( z_ratio**4 ) - if ( v_sum_dsc(i,j) > zero ) then - rhofsc(i,j,k) = 3.5_r_bl * rho_mix(i,j,k) & - * (v_top_dsc(i,j)/v_sum_dsc(i,j)) & - * (z_ratio**3) * (one-z_ratio) - end if - - else ! Not LockWhelan2006 - ! max to avoid rounding errors giving small negative numbers - rhokh_top(i,j,k) = rhokh_top(i,j,k) + & - rho_mix(i,j,k)*v_top_dsc(i,j)*g1*vkman* & - ( (max(one - kh_dsct_factor(i,j)*z_ratio,zero))**0.8_r_bl ) & - * z_pr * z_ratio - end if - end if - !------------------------------------------------------------- - ! Now momentum - !------------------------------------------------------------- - if ( zk_tq < zhsc(i,j) .and. & - zk_tq > zdsc_base(i,j) ) then - z_pr = zk_tq - zdsc_base(i,j) - zh_pr = zhsc(i,j) - zdsc_base(i,j) + else ! Not LockWhelan2006 ! max to avoid rounding errors giving small negative numbers - rhokm_dsct = & - 0.75_r_bl*rho_wet_tq(i,j,k-1)*v_top_dsc(i,j)*g1*vkman* & - ( (max(one - km_dsct_factor(i,j)*z_pr/zh_pr,zero))**0.8_r_bl ) & - * z_pr * z_pr / zh_pr - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,k) = tke_nl(i,j,k) + & - rhokm_dsct*c_tke*v_top_dsc(i,j)/dscdepth(i,j) - end if - rhokm_top(i,j,k) = rhokm_top(i,j,k) + rhokm_dsct + rhokh_top(i,j,k) = rhokh_top(i,j,k) + & + rho_mix(i,j,k)*v_top_dsc(i,j)*g1*vkman* & + ( (max(one - kh_dsct_factor(i,j)*z_ratio,zero))**0.8_r_bl ) & + * z_pr * z_ratio end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + end if + !------------------------------------------------------------- + ! Now momentum + !------------------------------------------------------------- + if ( zk_tq < zhsc(i,j) .and. & + zk_tq > zdsc_base(i,j) ) then + z_pr = zk_tq - zdsc_base(i,j) + zh_pr = zhsc(i,j) - zdsc_base(i,j) + ! max to avoid rounding errors giving small negative numbers + rhokm_dsct = & + 0.75_r_bl*rho_wet_tq(i,j,k-1)*v_top_dsc(i,j)*g1*vkman* & + ( (max(one - km_dsct_factor(i,j)*z_pr/zh_pr,zero))**0.8_r_bl ) & + * z_pr * z_pr / zh_pr + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,k) = tke_nl(i,j,k) + & + rhokm_dsct*c_tke*v_top_dsc(i,j)/dscdepth(i,j) + end if + rhokm_top(i,j,k) = rhokm_top(i,j,k) + rhokm_dsct + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) .and. & + zk_tq > zsml_base(i,j) ) then + BL_diag%tke(i,j,k) = c_tke*v_top(i,j)/zh(i,j) + end if + if ( zk_tq < zhsc(i,j) .and. & + zk_tq > zdsc_base(i,j) ) then ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) .and. & - zk_tq > zsml_base(i,j) ) then - BL_diag%tke(i,j,k) = c_tke*v_top(i,j)/zh(i,j) - end if - if ( zk_tq < zhsc(i,j) .and. & - zk_tq > zdsc_base(i,j) ) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*v_top_dsc(i,j)/dscdepth(i,j) ) - end if + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*v_top_dsc(i,j)/dscdepth(i,j) ) end if + end if - end do + !end do end do end do end do @@ -3095,109 +3099,109 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*c_ws**two_thirds) !$OMP do SCHEDULE(STATIC) - do jj = pdims%j_start, pdims%j_end, bl_segment_size + do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end - zk_uv = z_uv(i,j,k) - zk_tq = z_tq(i,j,k-1) + zk_uv = z_uv(i,j,k) + zk_tq = z_tq(i,j,k-1) - if (fb_surf(i,j) >= zero) then + if (fb_surf(i,j) >= zero) then - ! Calculate the free-convective scaling velocity at z(k) + ! Calculate the free-convective scaling velocity at z(k) - if (coupled(i,j)) then ! coupled - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) - else - wstar3 = zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) else - wstar3 = zh(i,j) * fb_surf(i,j) + wstar3 = zhsc(i,j) * fb_surf(i,j) end if + else + wstar3 = zh(i,j) * fb_surf(i,j) + end if - if (zk_uv <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation - w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) - else - ! Outer layer calculation - w_s_cubed_uv = c_ws * wstar3 - end if + if (zk_uv <= 0.1_r_bl*zh(i,j)) then + ! Surface layer calculation + w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) + else + ! Outer layer calculation + w_s_cubed_uv = c_ws * wstar3 + end if - if (zk_tq <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation - w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) - else - ! Outer layer calculation - w_s_cubed_tq = c_ws * wstar3 - end if + if (zk_tq <= 0.1_r_bl*zh(i,j)) then + ! Surface layer calculation + w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) + else + ! Outer layer calculation + w_s_cubed_tq = c_ws * wstar3 + end if - ! Turbulent velocity scale for scalars + ! Turbulent velocity scale for scalars - w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv ) & - **one_third - w_h_uv = w_m_neut/pr_neut + w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv ) & + **one_third + w_h_uv = w_m_neut/pr_neut - ! Also calc on TQ levels for W_M_TQ - w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq ) & - **one_third - w_h_tq = w_m_neut/pr_neut + ! Also calc on TQ levels for W_M_TQ + w_m_neut = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq ) & + **one_third + w_h_tq = w_m_neut/pr_neut - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) - ! Turbulent Prandtl number and velocity scale for scalars + ! Turbulent Prandtl number and velocity scale for scalars - Prandtl = pr_neut* & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*w_s_cubed_tq*w_m_neut ) / & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_tq* & - w_m_neut ) + Prandtl = pr_neut* & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*w_s_cubed_tq*w_m_neut ) / & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_tq* & + w_m_neut ) - w_m_tq = Prandtl * w_h_tq + w_m_tq = Prandtl * w_h_tq - if ( zk_uv < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKH(w_h,z/z_h) - !--------------------------------------------------------- + if ( zk_uv < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKH(w_h,z/z_h) + !--------------------------------------------------------- - rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & - ( one - ( zk_uv / zh(i,j) ) ) * & - ( one - ( zk_uv / zh(i,j) ) ) + rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & + ( one - ( zk_uv / zh(i,j) ) ) * & + ( one - ( zk_uv / zh(i,j) ) ) - end if - if ( zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKM(w_m,z/z_h) - !--------------------------------------------------------- + end if + if ( zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKM(w_m,z/z_h) + !--------------------------------------------------------- - rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman*zk_tq * & - ( one - ( zk_tq / zh(i,j) ) ) * & - ( one - ( zk_tq / zh(i,j) ) ) + rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman*zk_tq * & + ( one - ( zk_tq / zh(i,j) ) ) * & + ( one - ( zk_tq / zh(i,j) ) ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) ) then - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*w_m_tq/zh(i,j) ) - end if + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) ) then + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*w_m_tq/zh(i,j) ) end if end if - end do + end if + !end do end do end do end do @@ -3220,148 +3224,148 @@ subroutine excf_nl_9c ( & c_tke = 1.33_r_bl/(vkman*c_ws**two_thirds) !$OMP do SCHEDULE(STATIC) - do jj = pdims%j_start, pdims%j_end, bl_segment_size + do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end - ! Calculate the height of u,v-level above the surface - ! *APL: z0m removed from z in K(z) - zk_uv = z_uv(i,j,k) + ! Calculate the height of u,v-level above the surface + ! *APL: z0m removed from z in K(z) + zk_uv = z_uv(i,j,k) - ! Calculate the height of T,q-level above the surface + ! Calculate the height of T,q-level above the surface - zk_tq = z_tq(i,j,k-1) + zk_tq = z_tq(i,j,k-1) - if (fb_surf(i,j) >= zero) then + if (fb_surf(i,j) >= zero) then - ! Calculate the free-convective scaling velocity at z(k) + ! Calculate the free-convective scaling velocity at z(k) - if (zk_uv <= 0.1_r_bl*zh(i,j)) then + if (zk_uv <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation + ! Surface layer calculation - w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) - else + w_s_cubed_uv = 10.0_r_bl*c_ws * zk_uv * fb_surf(i,j) + else - ! Outer layer calculation + ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - w_s_cubed_uv = c_ws & - * ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) - else - w_s_cubed_uv = c_ws * zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + w_s_cubed_uv = c_ws & + * ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) else - w_s_cubed_uv = c_ws * zh(i,j) * fb_surf(i,j) + w_s_cubed_uv = c_ws * zhsc(i,j) * fb_surf(i,j) end if + else + w_s_cubed_uv = c_ws * zh(i,j) * fb_surf(i,j) end if + end if - if (zk_tq <= 0.1_r_bl*zh(i,j)) then + if (zk_tq <= 0.1_r_bl*zh(i,j)) then - ! Surface layer calculation + ! Surface layer calculation - w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) - else + w_s_cubed_tq = 10.0_r_bl*c_ws * zk_tq * fb_surf(i,j) + else - ! Outer layer calculation + ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - w_s_cubed_tq = c_ws & - * ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) - else - w_s_cubed_tq = c_ws * zhsc(i,j) * fb_surf(i,j) - end if + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + w_s_cubed_tq = c_ws & + * ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) * fb_surf(i,j) else - w_s_cubed_tq = c_ws * zh(i,j) * fb_surf(i,j) + w_s_cubed_tq = c_ws * zhsc(i,j) * fb_surf(i,j) end if + else + w_s_cubed_tq = c_ws * zh(i,j) * fb_surf(i,j) end if + end if - ! Turbulent velocity scale for momentum - - w_m_uv = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv) & - **one_third - - w_m_tq = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq) & - **one_third - - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) - - ! Turbulent Prandtl number and velocity scale for scalars - - Prandtl = pr_neut* & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*w_s_cubed_uv*w_m_uv ) / & - ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & - (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_uv* & - w_m_uv ) - w_h_uv = w_m_uv / Prandtl - - if ( zk_uv < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKH(w_h,z/z_h) - !--------------------------------------------------------- - - rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & - ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) * & - ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) - else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then - if ( zk_uv < zsml_top(i,j) ) then - ! Exponential decay from ZH but tends to zero - ! at zsml_top - rhokh(i,j,k) = rhokh_lcl(i,j) * & - exp(-(zk_uv-zh(i,j))/cu_depth_scale(i,j)) * & - (one-(zk_uv-zh(i,j))/(zsml_top(i,j)-zh(i,j))) - end if + ! Turbulent velocity scale for momentum + + w_m_uv = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_uv) & + **one_third + + w_m_tq = (v_s(i,j)*v_s(i,j)*v_s(i,j) + w_s_cubed_tq) & + **one_third + + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) + + ! Turbulent Prandtl number and velocity scale for scalars + + Prandtl = pr_neut* & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*w_s_cubed_uv*w_m_uv ) / & + ( v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j)*v_s_dbl(i,j) + & + (one/(c_ws*25.0_r_bl))*(pr_neut/pr_conv)*w_s_cubed_uv* & + w_m_uv ) + w_h_uv = w_m_uv / Prandtl + + if ( zk_uv < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKH(w_h,z/z_h) + !--------------------------------------------------------- + + rhokh(i,j,k) = rho_mix(i,j,k) * w_h_uv * vkman * zk_uv * & + ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) * & + ( one - kh_top_factor(i,j) * ( zk_uv / zh(i,j) ) ) + else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then + if ( zk_uv < zsml_top(i,j) ) then + ! Exponential decay from ZH but tends to zero + ! at zsml_top + rhokh(i,j,k) = rhokh_lcl(i,j) * & + exp(-(zk_uv-zh(i,j))/cu_depth_scale(i,j)) * & + (one-(zk_uv-zh(i,j))/(zsml_top(i,j)-zh(i,j))) end if - if ( zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate RHOKM(w_m,z/z_h) - !--------------------------------------------------------- + end if + if ( zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate RHOKM(w_m,z/z_h) + !--------------------------------------------------------- - rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman* zk_tq * & - ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) * & - ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) + rhokm(i,j,k) = rho_wet_tq(i,j,k-1)*w_m_tq*vkman* zk_tq * & + ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) * & + ( one - km_top_factor(i,j) * ( zk_tq / zh(i,j) ) ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + ! save Km/timescale for TKE diag, completed in bdy_expl2 + tke_nl(i,j,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) + end if + else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then + if ( zk_tq < zsml_top(i,j) ) then + ! Exponential decay from ZH but tends to zero + ! at zsml_top + rhokm(i,j,k) = prandtl_top(i,j) * rhokh_lcl(i,j) * & + exp(-(zk_tq-zh(i,j))/cu_depth_scale(i,j)) * & + (one-(zk_tq-zh(i,j))/(zsml_top(i,j)-zh(i,j))) + if (BL_diag%l_tke .and. var_diags_opt==split_tke_and_inv) then ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if - else if ( kprof_cu == klcl_entr .and. cumulus(i,j) ) then - if ( zk_tq < zsml_top(i,j) ) then - ! Exponential decay from ZH but tends to zero - ! at zsml_top - rhokm(i,j,k) = prandtl_top(i,j) * rhokh_lcl(i,j) * & - exp(-(zk_tq-zh(i,j))/cu_depth_scale(i,j)) * & - (one-(zk_tq-zh(i,j))/(zsml_top(i,j)-zh(i,j))) - if (BL_diag%l_tke .and. var_diags_opt==split_tke_and_inv) then - ! save Km/timescale for TKE diag, completed in bdy_expl2 - tke_nl(i,j,k) = tke_nl(i,j,k) + & - rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) - end if + tke_nl(i,j,k) = tke_nl(i,j,k) + & + rhokm(i,j,k)*c_tke*w_m_tq/zh(i,j) end if end if - if (BL_diag%l_tke .and. var_diags_opt == original_vars) then - ! save 1/timescale for TKE diag, completed in bdy_expl2 - if ( zk_tq < zsml_top(i,j) ) then - BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & - c_tke*w_m_tq/zh(i,j) ) - end if + end if + if (BL_diag%l_tke .and. var_diags_opt == original_vars) then + ! save 1/timescale for TKE diag, completed in bdy_expl2 + if ( zk_tq < zsml_top(i,j) ) then + BL_diag%tke(i,j,k) = max( BL_diag%tke(i,j,k), & + c_tke*w_m_tq/zh(i,j) ) end if - end if - end do + + end if + !end do end do end do end do @@ -3374,59 +3378,59 @@ subroutine excf_nl_9c ( & ng_stress == BrownGrant97_original ) then !$OMP do SCHEDULE(STATIC) - do jj = pdims%j_start, pdims%j_end, bl_segment_size + do ii = pdims%i_start, pdims%i_end, bl_segment_size do k = 2, bl_levels !cdir collapse - do j = jj, min(jj+bl_segment_size-1, pdims%j_end) - do i = pdims%i_start, pdims%i_end - zk_tq = z_tq(i,j,k-1) ! stresses are calc on theta-levs - if ( fb_surf(i,j) > zero .and. zk_tq < zh(i,j) ) then - !--------------------------------------------------------- - ! Calculate non-gradient stress function - ! (Brown and Grant 1997) - ! Shape function chosen such that non-gradient stress - ! either goes to zero at 0.1*ZH and ZH (ng_stress==BrownGrant97 - ! and ng_stress==BrownGrant97_limited), or goes to zero at - ! at the ground surface and ZH (ng_stress==BrownGrant97_original) - !--------------------------------------------------------- - ng_stress_calculate = .false. - - if (ng_stress == BrownGrant97 .or. & - ng_stress == BrownGrant97_limited) then - if ( zk_tq > 0.1_r_bl*zh(i,j) ) then - z_pr = zk_tq - 0.1_r_bl*zh(i,j) - zh_pr = 0.9_r_bl*zh(i,j) - ng_stress_calculate = .true. - end if - else if (ng_stress==BrownGrant97_original) then - z_pr = zk_tq - zh_pr = zh(i,j) + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end + zk_tq = z_tq(i,j,k-1) ! stresses are calc on theta-levs + if ( fb_surf(i,j) > zero .and. zk_tq < zh(i,j) ) then + !--------------------------------------------------------- + ! Calculate non-gradient stress function + ! (Brown and Grant 1997) + ! Shape function chosen such that non-gradient stress + ! either goes to zero at 0.1*ZH and ZH (ng_stress==BrownGrant97 + ! and ng_stress==BrownGrant97_limited), or goes to zero at + ! at the ground surface and ZH (ng_stress==BrownGrant97_original) + !--------------------------------------------------------- + ng_stress_calculate = .false. + + if (ng_stress == BrownGrant97 .or. & + ng_stress == BrownGrant97_limited) then + if ( zk_tq > 0.1_r_bl*zh(i,j) ) then + z_pr = zk_tq - 0.1_r_bl*zh(i,j) + zh_pr = 0.9_r_bl*zh(i,j) ng_stress_calculate = .true. end if + else if (ng_stress==BrownGrant97_original) then + z_pr = zk_tq + zh_pr = zh(i,j) + ng_stress_calculate = .true. + end if - if ( ng_stress_calculate ) then - ! Outer layer calculation - if (coupled(i,j)) then ! coupled and cloudy - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) - else - wstar3 = zhsc(i,j) * fb_surf(i,j) - end if + if ( ng_stress_calculate ) then + ! Outer layer calculation + if (coupled(i,j)) then ! coupled and cloudy + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) ) *fb_surf(i,j) else - wstar3 = zh(i,j) * fb_surf(i,j) + wstar3 = zhsc(i,j) * fb_surf(i,j) end if + else + wstar3 = zh(i,j) * fb_surf(i,j) + end if - ! Use the Holtslag and Boville velocity scale for - ! non-gradient stress stability dependence, as in BG97 - w_m_hb_3 = v_s(i,j)*v_s(i,j)*v_s(i,j) + 0.6_r_bl*wstar3 - f_ngstress(i,j,k) =(rho_wet_tq(i,j,k-1)/rhostar_gb(i,j)) & - * s_m * ( a_ngs * wstar3 / w_m_hb_3 ) & - * ( z_pr / zh_pr ) * ( one - ( z_pr / zh_pr ) ) * & - ( one - ( z_pr / zh_pr ) ) - end if ! ng_stress_calculate - end if ! fb_surf>0 and z0 and z Date: Wed, 4 Feb 2026 10:17:15 +0000 Subject: [PATCH 2/3] Commit up blocking behaviour which grants really good performance --- .../source/boundary_layer/excf_nl_9c.F90 | 1559 +++++++++-------- 1 file changed, 795 insertions(+), 764 deletions(-) diff --git a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 index 23c8cbc70..1246882cf 100644 --- a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 @@ -781,100 +781,103 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) ! do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - - rhokh_surf_ent(i,j) = zero - rhokh_top_ent(i,j) = zero - rhokh_dsct_ent(i,j) = zero - rhokm_surf_ent(i,j) = zero - rhokm_top_ent(i,j) = zero - rhokm_dsct_ent(i,j) = zero - v_top(i,j) = zero - v_sum(i,j) = zero - v_top_dsc(i,j) = zero - v_sum_dsc(i,j) = zero - rho_we(i,j) = zero - rho_we_sml(i,j) = zero - rho_we_dsc(i,j) = zero - - if (fb_surf(i,j) >= zero) then - - ! The calculations below involve v_s**4. In very rare circumstances - ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) - ! which is outside the range of single precision calculations - ! Hence we create a double precision version to enforce the correct - ! calculation - v_s_dbl(i,j) = v_s(i,j) - - ! Free-convective velocity scale cubed - - if (coupled(i,j)) then - if ( entr_smooth_dec == entr_taper_zh ) then - wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & - + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + + rhokh_surf_ent(i,j) = zero + rhokh_top_ent(i,j) = zero + rhokh_dsct_ent(i,j) = zero + rhokm_surf_ent(i,j) = zero + rhokm_top_ent(i,j) = zero + rhokm_dsct_ent(i,j) = zero + v_top(i,j) = zero + v_sum(i,j) = zero + v_top_dsc(i,j) = zero + v_sum_dsc(i,j) = zero + rho_we(i,j) = zero + rho_we_sml(i,j) = zero + rho_we_dsc(i,j) = zero + + if (fb_surf(i,j) >= zero) then + + ! The calculations below involve v_s**4. In very rare circumstances + ! v_s can be order(10^-11), therefore v_s**4 is order(10^-44) + ! which is outside the range of single precision calculations + ! Hence we create a double precision version to enforce the correct + ! calculation + v_s_dbl(i,j) = v_s(i,j) + + ! Free-convective velocity scale cubed + + if (coupled(i,j)) then + if ( entr_smooth_dec == entr_taper_zh ) then + wstar3 = ( svl_diff_frac(i,j) * zhsc(i,j) & + + (one-svl_diff_frac(i,j)) * zh(i,j) )* fb_surf(i,j) + else + wstar3 = zhsc(i,j) * fb_surf(i,j) + end if else - wstar3 = zhsc(i,j) * fb_surf(i,j) + wstar3 = zh(i,j) * fb_surf(i,j) end if - else - wstar3 = zh(i,j) * fb_surf(i,j) - end if - if (flux_grad == Locketal2000) then + if (flux_grad == Locketal2000) then - ! Turbulent velocity scale for momentum + ! Turbulent velocity scale for momentum - c_ws = 0.25_r_bl - w_m_top(i,j) = (v_s(i,j)*v_s(i,j)*v_s(i,j) + & - c_ws*wstar3)**one_third + c_ws = 0.25_r_bl + w_m_top(i,j) = (v_s(i,j)*v_s(i,j)*v_s(i,j) + & + c_ws*wstar3)**one_third - ! Turbulent Prandtl number and velocity scale for scalars - ! gives 0.375= zero ) then - !----------------------------------------------------------- - ! Calculate the surface buoyancy flux term - !----------------------------------------------------------- - zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) - sf_term = a_ent_1 * max ( zero , & - ( (one - zeta_s_fac) * bflux_surf(i,j) & - + zeta_s_fac * bflux_surf_sat(i,j) ) ) - !----------------------------------------------------------- - ! Calculate the surface shear term - !----------------------------------------------------------- - sf_shear_term = a_ent_shr * v_s(i,j) * v_s(i,j) * v_s(i,j) & - * rho_mix(i,j,k) / zh(i,j) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + !----------------------------------------------------------------------- + ! 1.2 Calculate top-of-b.l. entrainment mixing coefficients + ! and store b.l. top quantities for later use. + !----------------------------------------------------------------------- + ! FIRST the top of the SML (if not coupled) + !----------------------------------------------- + k = ntml(i,j)+1 + if ( .not. coupled(i,j) .and. fb_surf(i,j) >= zero ) then + !----------------------------------------------------------- + ! Calculate the surface buoyancy flux term + !----------------------------------------------------------- + zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) + sf_term = a_ent_1 * max ( zero , & + ( (one - zeta_s_fac) * bflux_surf(i,j) & + + zeta_s_fac * bflux_surf_sat(i,j) ) ) + !----------------------------------------------------------- + ! Calculate the surface shear term + !----------------------------------------------------------- + sf_shear_term = a_ent_shr * v_s(i,j) * v_s(i,j) * v_s(i,j) & + * rho_mix(i,j,k) / zh(i,j) + !----------------------------------------------------------- + ! Calculate the indirect radiative term + !----------------------------------------------------------- + zeta_r_sq = zeta_r(i,j)*zeta_r(i,j) + ir_term = ( bt_top(i,j)*zeta_r_sq + & + btt_top(i,j)*(one-zeta_r_sq) ) & + * a_ent_1 * df_top_over_cp(i,j) + !----------------------------------------------------------- + ! Calculate the evaporative term + !----------------------------------------------------------- + if ( db_top(i,j) > zero) then + zr = sqrt( zc(i,j) / zh(i,j) ) + evap_term = a_ent_2 * rho_mix(i,j,k) & + * chi_s_top(i,j) * chi_s_top(i,j) & + * zr * zr * zr * db_top_cld(i,j) & + * sqrt( zh(i,j) * db_top(i,j) ) + else + evap_term = zero + end if + !----------------------------------------------------------- + ! Combine forcing terms to calculate the representative + ! velocity scales + !----------------------------------------------------------- + v_sum(i,j) = ( (sf_term + sf_shear_term + & + ir_term + evap_term) & + * zh(i,j) /(a_ent_1*rho_mix(i,j,k)) )**one_third + v_top(i,j) = ( (ir_term+evap_term) * zh(i,j) & + / (a_ent_1*rho_mix(i,j,k)) )**one_third + !----------------------------------------------------------- + ! Calculate the direct radiative term + ! can only calculate for DB_TOP > 0 + !----------------------------------------------------------- + if ( db_top(i,j) > zero) then + dz_inv = min( v_sum(i,j)*v_sum(i,j) / db_top(i,j) ,100.0_r_bl ) + l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc(i,j)+rbl_eps) ) + alpha_t = one - exp(-one_half*dz_inv/l_rad) + ! Make enhancement due to buoyancy reversal feedback + alpha_t = alpha_t + br_fback(i,j)*(one-alpha_t) + dr_term = btc_top(i,j) * alpha_t * df_top_over_cp(i,j) + !---------------------------------------------------------- + ! Combine terms to calculate the entrainment + ! mixing coefficients + !---------------------------------------------------------- + zil_corr = c_t * ( (sf_term + sf_shear_term + & + ir_term + evap_term) / & + (rho_mix(i,j,k) * sqrt(zh(i,j))) )**two_thirds + + rho_we_sml(i,j) = (sf_term + sf_shear_term & + + ir_term + evap_term + dr_term) & + / ( db_top(i,j) + zil_corr ) + + rhokh_ent = rho_we_sml(i,j)/ rdz(i,j,k) + + frac_top = v_top(i,j) / ( v_top(i,j)+w_h_top(i,j)+rbl_eps ) + if ( l_use_sml_dsc_fixes ) then + ! If bug-fix is on, leave surface-driven entrainment flux set to + ! zero in cumulus layers if we are using the buoyancy-flux + ! integration method to calculate the non-local mixing profile + ! (code originally converged on a consistent negative buoyancy + ! flux at the top and then added the entrainment flux onto + ! this, which is double-counting). + l_apply_surf_ent = ( .not. ( ( kprof_cu==buoy_integ .or. & + kprof_cu==buoy_integ_low ) .and. & + cumulus(i,j) ) ) + else + ! If bug-fix is off, always add on surface-driven entrainment flux + l_apply_surf_ent = .true. + end if + if ( l_apply_surf_ent ) then + ! max to avoid rounding errors giving small negative numbers + rhokh_surf_ent(i,j) = max(rhokh_ent * ( one - frac_top ), zero) + rhokm_surf_ent(i,j) = prandtl_top(i,j) * rhokh_surf_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + end if + rhokh_top_ent(i,j) = rhokh_ent * frac_top + rhokm_top_ent(i,j) = 0.75_r_bl * rhokh_top_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + if (.not. l_use_var_fixes) then + rhokm(i,j,k) = rhokm_surf_ent(i,j) + rhokm_top(i,j,k) = rhokm_top_ent(i,j) + end if + + end if ! test on DB_TOP gt 0 + end if + !---------------------------------------------------------------- + ! then the top of the DSC (if coupled use ZHSC length-scale) + !---------------------------------------------------------------- + if ( ntdsc(i,j) > 0 ) then + k = ntdsc(i,j)+1 + if (coupled(i,j)) then + !-------------------------------------------------------- + ! Calculate the surface buoyancy flux term + !-------------------------------------------------------- + zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) + sf_term = a_ent_1 * max ( zero , & + ( (one - zeta_s_fac) * bflux_surf(i,j) & + + zeta_s_fac * bflux_surf_sat(i,j) ) ) + !-------------------------------------------------------- + ! Calculate the surface shear term + !-------------------------------------------------------- + if (fb_surf(i,j) >= zero) then + sf_shear_term = a_ent_shr * v_s(i,j)*v_s(i,j)*v_s(i,j) & + * rho_mix(i,j,k) / zhsc(i,j) + else + sf_shear_term = zero + end if + if ( entr_smooth_dec == on .or. entr_smooth_dec == entr_taper_zh ) then + ! taper surface terms to zero depending on svl_diff_frac + sf_term = sf_term * svl_diff_frac(i,j) + sf_shear_term = sf_shear_term * svl_diff_frac(i,j) + end if + else + sf_term = zero + sf_shear_term = zero + end if !----------------------------------------------------------- ! Calculate the indirect radiative term !----------------------------------------------------------- - zeta_r_sq = zeta_r(i,j)*zeta_r(i,j) - ir_term = ( bt_top(i,j)*zeta_r_sq + & - btt_top(i,j)*(one-zeta_r_sq) ) & - * a_ent_1 * df_top_over_cp(i,j) + zeta_r_sq = zeta_r_dsc(i,j)*zeta_r_dsc(i,j) + ir_term = ( bt_dsct(i,j)*zeta_r_sq + & + btt_dsct(i,j)*(one-zeta_r_sq) ) & + * a_ent_1 * df_dsct_over_cp(i,j) !----------------------------------------------------------- ! Calculate the evaporative term !----------------------------------------------------------- - if ( db_top(i,j) > zero) then - zr = sqrt( zc(i,j) / zh(i,j) ) - evap_term = a_ent_2 * rho_mix(i,j,k) & - * chi_s_top(i,j) * chi_s_top(i,j) & - * zr * zr * zr * db_top_cld(i,j) & - * sqrt( zh(i,j) * db_top(i,j) ) - else - evap_term = zero - end if + if (db_dsct(i,j) > zero) then + zr = sqrt( zc_dsc(i,j) / dscdepth(i,j) ) + evap_term = a_ent_2 * rho_mix(i,j,k) & + * chi_s_dsct(i,j) * chi_s_dsct(i,j) & + * zr * zr * zr * db_dsct_cld(i,j) & + * sqrt( dscdepth(i,j) * db_dsct(i,j) ) + else + evap_term = zero + end if !----------------------------------------------------------- ! Combine forcing terms to calculate the representative ! velocity scales !----------------------------------------------------------- - v_sum(i,j) = ( (sf_term + sf_shear_term + & - ir_term + evap_term) & - * zh(i,j) /(a_ent_1*rho_mix(i,j,k)) )**one_third - v_top(i,j) = ( (ir_term+evap_term) * zh(i,j) & - / (a_ent_1*rho_mix(i,j,k)) )**one_third + v_sum_dsc(i,j) = ( (sf_term + sf_shear_term + & + ir_term + evap_term) & + * dscdepth(i,j) / (a_ent_1*rho_mix(i,j,k)) )**one_third + v_top_dsc(i,j) =( (ir_term + evap_term) * dscdepth(i,j) / & + (a_ent_1*rho_mix(i,j,k)) )**one_third !----------------------------------------------------------- ! Calculate the direct radiative term - ! can only calculate for DB_TOP > 0 !----------------------------------------------------------- - if ( db_top(i,j) > zero) then - dz_inv = min( v_sum(i,j)*v_sum(i,j) / db_top(i,j) ,100.0_r_bl ) - l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc(i,j)+rbl_eps) ) - alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback - alpha_t = alpha_t + br_fback(i,j)*(one-alpha_t) - dr_term = btc_top(i,j) * alpha_t * df_top_over_cp(i,j) - !---------------------------------------------------------- - ! Combine terms to calculate the entrainment - ! mixing coefficients - !---------------------------------------------------------- - zil_corr = c_t * ( (sf_term + sf_shear_term + & - ir_term + evap_term) / & - (rho_mix(i,j,k) * sqrt(zh(i,j))) )**two_thirds - - rho_we_sml(i,j) = (sf_term + sf_shear_term & - + ir_term + evap_term + dr_term) & - / ( db_top(i,j) + zil_corr ) - - rhokh_ent = rho_we_sml(i,j)/ rdz(i,j,k) - - frac_top = v_top(i,j) / ( v_top(i,j)+w_h_top(i,j)+rbl_eps ) - if ( l_use_sml_dsc_fixes ) then - ! If bug-fix is on, leave surface-driven entrainment flux set to - ! zero in cumulus layers if we are using the buoyancy-flux - ! integration method to calculate the non-local mixing profile - ! (code originally converged on a consistent negative buoyancy - ! flux at the top and then added the entrainment flux onto - ! this, which is double-counting). - l_apply_surf_ent = ( .not. ( ( kprof_cu==buoy_integ .or. & - kprof_cu==buoy_integ_low ) .and. & - cumulus(i,j) ) ) - else - ! If bug-fix is off, always add on surface-driven entrainment flux - l_apply_surf_ent = .true. - end if - if ( l_apply_surf_ent ) then - ! max to avoid rounding errors giving small negative numbers - rhokh_surf_ent(i,j) = max(rhokh_ent * ( one - frac_top ), zero) - rhokm_surf_ent(i,j) = prandtl_top(i,j) * rhokh_surf_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - end if - rhokh_top_ent(i,j) = rhokh_ent * frac_top - rhokm_top_ent(i,j) = 0.75_r_bl * rhokh_top_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - if (.not. l_use_var_fixes) then - rhokm(i,j,k) = rhokm_surf_ent(i,j) - rhokm_top(i,j,k) = rhokm_top_ent(i,j) - end if - - end if ! test on DB_TOP gt 0 - end if - !---------------------------------------------------------------- - ! then the top of the DSC (if coupled use ZHSC length-scale) - !---------------------------------------------------------------- - if ( ntdsc(i,j) > 0 ) then - k = ntdsc(i,j)+1 - if (coupled(i,j)) then - !-------------------------------------------------------- - ! Calculate the surface buoyancy flux term - !-------------------------------------------------------- - zeta_s_fac = (one - zeta_s(i,j)) * (one - zeta_s(i,j)) - sf_term = a_ent_1 * max ( zero , & - ( (one - zeta_s_fac) * bflux_surf(i,j) & - + zeta_s_fac * bflux_surf_sat(i,j) ) ) - !-------------------------------------------------------- - ! Calculate the surface shear term - !-------------------------------------------------------- - if (fb_surf(i,j) >= zero) then - sf_shear_term = a_ent_shr * v_s(i,j)*v_s(i,j)*v_s(i,j) & - * rho_mix(i,j,k) / zhsc(i,j) - else - sf_shear_term = zero - end if - if ( entr_smooth_dec == on .or. entr_smooth_dec == entr_taper_zh ) then - ! taper surface terms to zero depending on svl_diff_frac - sf_term = sf_term * svl_diff_frac(i,j) - sf_shear_term = sf_shear_term * svl_diff_frac(i,j) - end if - else - sf_term = zero - sf_shear_term = zero - end if - !----------------------------------------------------------- - ! Calculate the indirect radiative term - !----------------------------------------------------------- - zeta_r_sq = zeta_r_dsc(i,j)*zeta_r_dsc(i,j) - ir_term = ( bt_dsct(i,j)*zeta_r_sq + & - btt_dsct(i,j)*(one-zeta_r_sq) ) & - * a_ent_1 * df_dsct_over_cp(i,j) - !----------------------------------------------------------- - ! Calculate the evaporative term - !----------------------------------------------------------- - if (db_dsct(i,j) > zero) then - zr = sqrt( zc_dsc(i,j) / dscdepth(i,j) ) - evap_term = a_ent_2 * rho_mix(i,j,k) & - * chi_s_dsct(i,j) * chi_s_dsct(i,j) & - * zr * zr * zr * db_dsct_cld(i,j) & - * sqrt( dscdepth(i,j) * db_dsct(i,j) ) - else - evap_term = zero + if (db_dsct(i,j) > zero) then + dz_inv = min( v_sum_dsc(i,j)*v_sum_dsc(i,j)/db_dsct(i,j), & + 100.0_r_bl ) + l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc_dsc(i,j)+one) ) + alpha_t = one - exp(-one_half*dz_inv/l_rad) + ! Make enhancement due to buoyancy reversal feedback + alpha_t = alpha_t + br_fback_dsc(i,j)*(one-alpha_t) + dr_term = btc_top(i,j) * alpha_t * df_dsct_over_cp(i,j) + !---------------------------------------------------------- + ! Finally combine terms to calculate the entrainment + ! rate and mixing coefficients + !---------------------------------------------------------- + zil_corr = c_t * ( (sf_term + sf_shear_term + & + ir_term + evap_term) / & + (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds + rho_we_dsc(i,j) = ( sf_term + sf_shear_term & + + ir_term + evap_term + dr_term ) & + / ( db_dsct(i,j) + zil_corr ) + rhokh_dsct_ent(i,j) = rho_we_dsc(i,j)/ rdz(i,j,k) + rhokm_dsct_ent(i,j) = 0.75_r_bl * rhokh_dsct_ent(i,j) & + * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & + * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) + if (.not. l_use_var_fixes) rhokm_top(i,j,k) = rhokm_dsct_ent(i,j) + end if ! test on DB_DSCT gt 0 end if - !----------------------------------------------------------- - ! Combine forcing terms to calculate the representative - ! velocity scales - !----------------------------------------------------------- - v_sum_dsc(i,j) = ( (sf_term + sf_shear_term + & - ir_term + evap_term) & - * dscdepth(i,j) / (a_ent_1*rho_mix(i,j,k)) )**one_third - v_top_dsc(i,j) =( (ir_term + evap_term) * dscdepth(i,j) / & - (a_ent_1*rho_mix(i,j,k)) )**one_third - !----------------------------------------------------------- - ! Calculate the direct radiative term - !----------------------------------------------------------- - if (db_dsct(i,j) > zero) then - dz_inv = min( v_sum_dsc(i,j)*v_sum_dsc(i,j)/db_dsct(i,j), & - 100.0_r_bl ) - l_rad = 15.0_r_bl * max( one, 200.0_r_bl/(zc_dsc(i,j)+one) ) - alpha_t = one - exp(-one_half*dz_inv/l_rad) - ! Make enhancement due to buoyancy reversal feedback - alpha_t = alpha_t + br_fback_dsc(i,j)*(one-alpha_t) - dr_term = btc_top(i,j) * alpha_t * df_dsct_over_cp(i,j) - !---------------------------------------------------------- - ! Finally combine terms to calculate the entrainment - ! rate and mixing coefficients - !---------------------------------------------------------- - zil_corr = c_t * ( (sf_term + sf_shear_term + & - ir_term + evap_term) / & - (rho_mix(i,j,k) * sqrt(dscdepth(i,j))) )**two_thirds - rho_we_dsc(i,j) = ( sf_term + sf_shear_term & - + ir_term + evap_term + dr_term ) & - / ( db_dsct(i,j) + zil_corr ) - rhokh_dsct_ent(i,j) = rho_we_dsc(i,j)/ rdz(i,j,k) - rhokm_dsct_ent(i,j) = 0.75_r_bl * rhokh_dsct_ent(i,j) & - * rdz(i,j,k) * (z_uv(i,j,k)-z_uv(i,j,k-1)) & - * rho_wet_tq(i,j,k-1) / rho_mix(i,j,k) - if (.not. l_use_var_fixes) rhokm_top(i,j,k) = rhokm_dsct_ent(i,j) - end if ! test on DB_DSCT gt 0 - end if -end do + end do ! I +end do ! II !end do !$OMP end do @@ -1220,57 +1226,60 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - if ( bflux_surf(i,j) > zero) then - ! can only be coupled to an unstable SML - if ( ntdsc(i,j) > 2 ) then - ! well-resolved DSC layer - ! ...test for recoupling with SML - test_well_mixed(i,j) = .true. - z_inv(i,j) = zhsc(i,j) - z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) - cf_ml(i,j) = cf_dsc(i,j) - v_ktop(i,j) = v_top_dsc(i,j) - v_ksum(i,j) = v_sum_dsc(i,j) - df_ctop(i,j)= df_dsct_over_cp(i,j) - rho_we(i,j) = rho_we_dsc(i,j) - dsl(i,j) = dsl_dsc(i,j) - dqw(i,j) = dqw_dsc(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) - ntop(i,j) = ntdsc(i,j) - ! assuming wb goes to zero by the lowest of ZH - ! or cloud-base, but above surface layer - else if ( .not. dsc(i,j) .and. .not. cumulus(i,j) .and. & - ntml(i,j) > 2) then - ! well-resolved SML - ! ...test for decoupling - ! Note: code can only deal with one DSC layer at a time so - ! can't decouple SML if a DSC layer already exists. - ! --------------------------------------------------------- - ! If the BL layer is cloud-free then use a less restrictive - ! threshold - ideally, the parcel ascent would have - ! found the correct BL top in this case but this test is - ! kept to keep negative buoyancy fluxes under control - ! (ie. DEC_THRES_CLEAR=1 ensures wbn_int < |wbp_int|) - if (abs(zc(i,j)) < rbl_eps) dec_thres(i,j) = dec_thres_clear - ! --------------------------------------------------------- - test_well_mixed(i,j) = .true. - z_inv(i,j) = zh(i,j) - z_cbase(i,j)= z_inv(i,j) - zc(i,j) - cf_ml(i,j) = cf_sml(i,j) - v_ktop(i,j) = v_top(i,j) - v_ksum(i,j) = v_sum(i,j) - df_ctop(i,j)= df_top_over_cp(i,j) - rho_we(i,j) = rho_we_sml(i,j) - dsl(i,j) = dsl_sml(i,j) - dqw(i,j) = dqw_sml(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zh(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zh(i,j) - ntop(i,j) = ntml(i,j) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + if ( bflux_surf(i,j) > zero) then + ! can only be coupled to an unstable SML + if ( ntdsc(i,j) > 2 ) then + ! well-resolved DSC layer + ! ...test for recoupling with SML + test_well_mixed(i,j) = .true. + z_inv(i,j) = zhsc(i,j) + z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) + cf_ml(i,j) = cf_dsc(i,j) + v_ktop(i,j) = v_top_dsc(i,j) + v_ksum(i,j) = v_sum_dsc(i,j) + df_ctop(i,j)= df_dsct_over_cp(i,j) + rho_we(i,j) = rho_we_dsc(i,j) + dsl(i,j) = dsl_dsc(i,j) + dqw(i,j) = dqw_dsc(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) + ntop(i,j) = ntdsc(i,j) + ! assuming wb goes to zero by the lowest of ZH + ! or cloud-base, but above surface layer + else if ( .not. dsc(i,j) .and. .not. cumulus(i,j) .and. & + ntml(i,j) > 2) then + ! well-resolved SML + ! ...test for decoupling + ! Note: code can only deal with one DSC layer at a time so + ! can't decouple SML if a DSC layer already exists. + ! --------------------------------------------------------- + ! If the BL layer is cloud-free then use a less restrictive + ! threshold - ideally, the parcel ascent would have + ! found the correct BL top in this case but this test is + ! kept to keep negative buoyancy fluxes under control + ! (ie. DEC_THRES_CLEAR=1 ensures wbn_int < |wbp_int|) + if (abs(zc(i,j)) < rbl_eps) dec_thres(i,j) = dec_thres_clear + ! --------------------------------------------------------- + test_well_mixed(i,j) = .true. + z_inv(i,j) = zh(i,j) + z_cbase(i,j)= z_inv(i,j) - zc(i,j) + cf_ml(i,j) = cf_sml(i,j) + v_ktop(i,j) = v_top(i,j) + v_ksum(i,j) = v_sum(i,j) + df_ctop(i,j)= df_top_over_cp(i,j) + rho_we(i,j) = rho_we_sml(i,j) + dsl(i,j) = dsl_sml(i,j) + dqw(i,j) = dqw_sml(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zh(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zh(i,j) + ntop(i,j) = ntml(i,j) + end if end if - end if -end do ! I + end do ! I +end do ! II !end do ! J !$OMP end do @@ -1354,97 +1363,100 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) ) then - dzrad(i,j) = 100.0_r_bl - - if ( dzrad_disc_opt == dzrad_ntm1 ) then - ! Original discretisation of cloud-top radiatively-cooled layer; - ! set the base of the layer at theta-level ntdsc-1 - ! (which is between 1.5 and 2.5 model-levels below cloud-top, - ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), - ! then lower it by extra whole model-levels if needed to increase the - ! depth to 100m. - - ntop(i,j) = ntop(i,j) - 1 - if ( ntop(i,j) > 2 ) then - inner_loop1: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) - ntop(i,j) = ntop(i,j) - 1 - if (ntop(i,j) == 2 ) then - exit inner_loop1 - end if - end do inner_loop1 - end if - dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) - - else if ( dzrad_disc_opt == dzrad_1p5dz ) then - ! Alternative method; - ! set the base of the layer 1.5 model-levels below cloud-top so that - ! it always moves smoothly with zhsc instead of jumping suddenly - ! (1.5 levels is the minimum distance below zhsc which will always - ! keep the base of the layer below theta-level ntdsc, which is a - ! requirement of the discretisation). - ! The depth is also forced to be at least 100m. - - ! Find smoothly-varying height 1.5 model-levels below z_inv - k_inv = ntop(i,j) - if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then - interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) - z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) - else - interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) - z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) - end if - ! Go down to 100m below z_inv (or cloud-base) if that is lower, - ! but don't allow lower than theta-level 1 - z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) - z_pr = max( z_pr, z_tq(i,j,1) ) - ! Set ntop to the rho-level straddling the base of the layer - do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + if ( test_well_mixed(i,j) ) then + dzrad(i,j) = 100.0_r_bl + + if ( dzrad_disc_opt == dzrad_ntm1 ) then + ! Original discretisation of cloud-top radiatively-cooled layer; + ! set the base of the layer at theta-level ntdsc-1 + ! (which is between 1.5 and 2.5 model-levels below cloud-top, + ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), + ! then lower it by extra whole model-levels if needed to increase the + ! depth to 100m. + ntop(i,j) = ntop(i,j) - 1 - end do - ! Set new cloud-top radiatively-cooled layer-depth - dzrad(i,j) = z_inv(i,j) - z_pr - - ! Fraction of rho-level ntop below the dzrad layer - interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) - ! Add on SL and qw differences between the discontinuous inversion base - ! and the base of the dzrad layer - dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & - - interp *sl(i,j,ntop(i,j)) - dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & - - interp *qw(i,j,ntop(i,j)) - - end if ! ( dzrad_disc_opt ) - - wsl_dzrad_int = dzrad(i,j) * & - ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) - wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) - - wb_dzrad_int(i,j) = wsl_dzrad_int * ( & - (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & - + wqw_dzrad_int * ( & - (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) - wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) - wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) - - ! Include WB_DZRAD_INT in WBP_INT as it set to be >0 - wbp_int(i,j) = wbp_int(i,j) + wb_dzrad_int(i,j) + if ( ntop(i,j) > 2 ) then + inner_loop1: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) + ntop(i,j) = ntop(i,j) - 1 + if (ntop(i,j) == 2 ) then + exit inner_loop1 + end if + end do inner_loop1 + end if + dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) + + else if ( dzrad_disc_opt == dzrad_1p5dz ) then + ! Alternative method; + ! set the base of the layer 1.5 model-levels below cloud-top so that + ! it always moves smoothly with zhsc instead of jumping suddenly + ! (1.5 levels is the minimum distance below zhsc which will always + ! keep the base of the layer below theta-level ntdsc, which is a + ! requirement of the discretisation). + ! The depth is also forced to be at least 100m. + + ! Find smoothly-varying height 1.5 model-levels below z_inv + k_inv = ntop(i,j) + if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then + interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + z_pr = (one-interp) * z_tq(i,j,k_inv-1) & + + interp * z_uv(i,j,k_inv) + else + interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + z_pr = (one-interp) * z_uv(i,j,k_inv) & + + interp * z_tq(i,j,k_inv) + end if + ! Go down to 100m below z_inv (or cloud-base) if that is lower, + ! but don't allow lower than theta-level 1 + z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) + z_pr = max( z_pr, z_tq(i,j,1) ) + ! Set ntop to the rho-level straddling the base of the layer + do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) + ntop(i,j) = ntop(i,j) - 1 + end do + ! Set new cloud-top radiatively-cooled layer-depth + dzrad(i,j) = z_inv(i,j) - z_pr + + ! Fraction of rho-level ntop below the dzrad layer + interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + ! Add on SL and qw differences between the discontinuous inversion base + ! and the base of the dzrad layer + dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & + - interp *sl(i,j,ntop(i,j)) + dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & + - interp *qw(i,j,ntop(i,j)) + + end if ! ( dzrad_disc_opt ) + + wsl_dzrad_int = dzrad(i,j) * & + ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) + wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) + + wb_dzrad_int(i,j) = wsl_dzrad_int * ( & + (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & + + wqw_dzrad_int * ( & + (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) + wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) + wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) + + ! Include WB_DZRAD_INT in WBP_INT as it set to be >0 + wbp_int(i,j) = wbp_int(i,j) + wb_dzrad_int(i,j) - else + else - wb_dzrad_int(i,j) = -one ! To identify not calculated + wb_dzrad_int(i,j) = -one ! To identify not calculated - end if -end do ! I + end if + end do ! I +end do ! II !end do ! J !$OMP end do @@ -1481,33 +1493,36 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end - if ( test_well_mixed(i,j) ) then - if ( kwb0(i,j) == ntml(i,j) ) then - zwb0(i,j) = zh(i,j) - else if ( kwb0(i,j) == 2 ) then - zwb0(i,j) = z_uv(i,j,2) + if ( test_well_mixed(i,j) ) then + if ( kwb0(i,j) == ntml(i,j) ) then + zwb0(i,j) = zh(i,j) + else if ( kwb0(i,j) == 2 ) then + zwb0(i,j) = z_uv(i,j,2) + else + k=kwb0(i,j) + ! now DB_GA_DRY(K) le 0 and DB_GA_DRY(K-1) gt 0 + ! so interpolate: + db_ratio = db_ga_dry_n(i,j,k-1) & + / ( db_ga_dry_n(i,j,k-1) - db_ga_dry_n(i,j,k) ) + db_ratio = max( zero, db_ratio ) ! trap for rounding error + zwb0(i,j)=z_uv(i,j,k-1) + & + db_ratio * (z_uv(i,j,k)-z_uv(i,j,k-1)) + end if + wb_surf_int(i,j) = bflux_surf(i,j) * z_tq(i,j,ksurf(i,j)) * & + ( one - z_tq(i,j,ksurf(i,j))/(2.0_r_bl*zwb0(i,j))) + wb_surf_int(i,j) = max( rbl_eps, wb_surf_int(i,j) ) else - k=kwb0(i,j) - ! now DB_GA_DRY(K) le 0 and DB_GA_DRY(K-1) gt 0 - ! so interpolate: - db_ratio = db_ga_dry_n(i,j,k-1) & - / ( db_ga_dry_n(i,j,k-1) - db_ga_dry_n(i,j,k) ) - db_ratio = max( zero, db_ratio ) ! trap for rounding error - zwb0(i,j)=z_uv(i,j,k-1) + & - db_ratio * (z_uv(i,j,k)-z_uv(i,j,k-1)) + ! only include surface layer contribution for unstable mixing + wb_surf_int(i,j) = rbl_eps end if - wb_surf_int(i,j) = bflux_surf(i,j) * z_tq(i,j,ksurf(i,j)) * & - ( one - z_tq(i,j,ksurf(i,j))/(2.0_r_bl*zwb0(i,j))) - wb_surf_int(i,j) = max( rbl_eps, wb_surf_int(i,j) ) - else - ! only include surface layer contribution for unstable mixing - wb_surf_int(i,j) = rbl_eps - end if - wbp_int(i,j) = wbp_int(i,j) + wb_surf_int(i,j) ! must be >0 + wbp_int(i,j) = wbp_int(i,j) + wb_surf_int(i,j) ! must be >0 + end do ! I end do ! I !end do ! J !$OMP end do @@ -1672,163 +1687,166 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - if (l_wtrac) iset_wtrac(i,j) = 0 - if ( test_well_mixed(i,j) ) then - - wb_ratio(i,j) = wbn_int(i,j)/wbp_int(i,j) - - if ( wb_ratio(i,j) <= dec_thres(i,j) ) then - ! No need to test depth of mixing any further as within - ! well-mixed layer buoyancy flux integral criteria. - ! SML will simply stay well-mixed (and so use defaults) - ksurf_iterate(i,j)= .false. - ktop_iterate(i,j) = .false. - if ( dsc(i,j) ) then - ! Recouple DSC with SML: - dsc_removed(i,j) = 2 ! set indicator flag - ! move surface driven entrainment - ! RHOKH(z_i) = rho * w_e * DZL and w_e ~ 1/DB_TOP, so: - if ( db_top(i,j) > zero .and. db_dsct(i,j) > 0.01_r_bl ) then - ! can't calc Zil. term - rhokh_surf_ent(i,j) = rhokh_surf_ent(i,j) * & - ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & - rdz(i,j,ntdsc(i,j)+1) ) - rhokm_surf_ent(i,j) = rhokm_surf_ent(i,j) * & - ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & - rdz(i,j,ntml(i,j)+1) ) / & - ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & - rdz(i,j,ntdsc(i,j)+1) ) - if (.not. l_use_var_fixes) & - rhokm(i,j,ntdsc(i,j)+1) = rhokm_surf_ent(i,j) - end if - ! redesignate top-driven entrainment at ZHSC - ! (ignore that calculated at ZH) - rhokh_top_ent(i,j) = rhokh_dsct_ent(i,j) - rhokm_top_ent(i,j) = rhokm_dsct_ent(i,j) - zh(i,j) = zhsc(i,j) - ntml(i,j) = ntdsc(i,j) - v_top(i,j) = v_top_dsc(i,j) - zsml_base(i,j) = 0.1_r_bl * zh(i,j) - zc(i,j) = zc_dsc(i,j) - zhsc(i,j) = zero - ntdsc(i,j) = 0 - v_top_dsc(i,j) = zero - zdsc_base(i,j) = zero - zc_dsc(i,j) = zero - ft_nt_zh(i,j) = ft_nt_zhsc(i,j) - ft_nt_zhsc(i,j) = zero - fq_nt_zh(i,j) = fq_nt_zhsc(i,j) - fq_nt_zhsc(i,j) = zero - if (l_wtrac) iset_wtrac(i,j) = 1 - dsc(i,j) = .false. - dzh(i,j) = rmdi ! SML inversion subsumed into new mixed layer - cumulus(i,j) = .false. - coupled(i,j) = .false. - end if ! recoupled DSC layer - else ! buoyancy flux threshold violated +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + if (l_wtrac) iset_wtrac(i,j) = 0 + if ( test_well_mixed(i,j) ) then - if ( l_use_sml_dsc_fixes ) then - ! Under fix l_use_sml_dsc_fixes: - ! If there is already a DSC layer diagnosed above the existing - ! ntml, we have so far just tested for recoupling of that DSC - ! with the SML; we have not tested whether the existing SML - ! should split off another decoupled layer. - ! In this case, we don't need to turn on ktop_iterate - ! (already setup earlier), and we should only turn on ksurf_iterate - ! if the SML is cumulus capped and using buoyancy flux integration - ! at cumulus points in general. + wb_ratio(i,j) = wbn_int(i,j)/wbp_int(i,j) + if ( wb_ratio(i,j) <= dec_thres(i,j) ) then + ! No need to test depth of mixing any further as within + ! well-mixed layer buoyancy flux integral criteria. + ! SML will simply stay well-mixed (and so use defaults) + ksurf_iterate(i,j)= .false. + ktop_iterate(i,j) = .false. if ( dsc(i,j) ) then - ! Existing DSC layer confirmed not to be merged with the SML. - if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low ) ) then - ! Cumulus layer under the DSC; only turn on ksurf_iterate if - ! using appropriate kprof_cu option. - ksurf_iterate(i,j) = .true. - dec_thres(i,j) = dec_thres_cu + ! Recouple DSC with SML: + dsc_removed(i,j) = 2 ! set indicator flag + ! move surface driven entrainment + ! RHOKH(z_i) = rho * w_e * DZL and w_e ~ 1/DB_TOP, so: + if ( db_top(i,j) > zero .and. db_dsct(i,j) > 0.01_r_bl ) then + ! can't calc Zil. term + rhokh_surf_ent(i,j) = rhokh_surf_ent(i,j) * & + ( rho_mix(i,j,ntdsc(i,j)+1) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_mix(i,j,ntml(i,j)+1) * db_dsct(i,j) * & + rdz(i,j,ntdsc(i,j)+1) ) + rhokm_surf_ent(i,j) = rhokm_surf_ent(i,j) * & + ( rho_wet_tq(i,j,ntdsc(i,j)) * db_top(i,j) * & + rdz(i,j,ntml(i,j)+1) ) / & + ( rho_wet_tq(i,j,ntml(i,j)) * db_dsct(i,j) * & + rdz(i,j,ntdsc(i,j)+1) ) + if (.not. l_use_var_fixes) & + rhokm(i,j,ntdsc(i,j)+1) = rhokm_surf_ent(i,j) end if - else - ! No existing DSC so we're splitting a new DSC off the SML; - ! turn on ktop_iterate to find the new DSC's base - ktop_iterate(i,j) = .true. - if ( (.not. cumulus(i,j)) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low ) then - ! Switch to using b-flux integration to find new lower SML-top, - ! but not in Cu layers unless using appropriate kprof_cu option. - ksurf_iterate(i,j) = .true. - if ( cumulus(i,j) ) dec_thres(i,j) = dec_thres_cu + ! redesignate top-driven entrainment at ZHSC + ! (ignore that calculated at ZH) + rhokh_top_ent(i,j) = rhokh_dsct_ent(i,j) + rhokm_top_ent(i,j) = rhokm_dsct_ent(i,j) + zh(i,j) = zhsc(i,j) + ntml(i,j) = ntdsc(i,j) + v_top(i,j) = v_top_dsc(i,j) + zsml_base(i,j) = 0.1_r_bl * zh(i,j) + zc(i,j) = zc_dsc(i,j) + zhsc(i,j) = zero + ntdsc(i,j) = 0 + v_top_dsc(i,j) = zero + zdsc_base(i,j) = zero + zc_dsc(i,j) = zero + ft_nt_zh(i,j) = ft_nt_zhsc(i,j) + ft_nt_zhsc(i,j) = zero + fq_nt_zh(i,j) = fq_nt_zhsc(i,j) + fq_nt_zhsc(i,j) = zero + if (l_wtrac) iset_wtrac(i,j) = 1 + dsc(i,j) = .false. + dzh(i,j) = rmdi ! SML inversion subsumed into new mixed layer + cumulus(i,j) = .false. + coupled(i,j) = .false. + end if ! recoupled DSC layer + else ! buoyancy flux threshold violated + + if ( l_use_sml_dsc_fixes ) then + ! Under fix l_use_sml_dsc_fixes: + ! If there is already a DSC layer diagnosed above the existing + ! ntml, we have so far just tested for recoupling of that DSC + ! with the SML; we have not tested whether the existing SML + ! should split off another decoupled layer. + ! In this case, we don't need to turn on ktop_iterate + ! (already setup earlier), and we should only turn on ksurf_iterate + ! if the SML is cumulus capped and using buoyancy flux integration + ! at cumulus points in general. + + if ( dsc(i,j) ) then + ! Existing DSC layer confirmed not to be merged with the SML. + if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low ) ) then + ! Cumulus layer under the DSC; only turn on ksurf_iterate if + ! using appropriate kprof_cu option. + ksurf_iterate(i,j) = .true. + dec_thres(i,j) = dec_thres_cu + end if + else + ! No existing DSC so we're splitting a new DSC off the SML; + ! turn on ktop_iterate to find the new DSC's base + ktop_iterate(i,j) = .true. + if ( (.not. cumulus(i,j)) .or. kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low ) then + ! Switch to using b-flux integration to find new lower SML-top, + ! but not in Cu layers unless using appropriate kprof_cu option. + ksurf_iterate(i,j) = .true. + if ( cumulus(i,j) ) dec_thres(i,j) = dec_thres_cu + end if end if - end if - else ! ( .not. l_use_sml_dsc_fixes ) - ! If the bug-fix is off, we turn on ksurf_iterate at non-cumulus - ! points even if dsc was already true - ! (when all we've done is confirm the DSC shouldn't merge with the - ! SML, so we had no reason to be meddling with the SML). - ! This creates some inconsistencies down the line... - - !--------------------------------- - ! Extent of mixing must be reduced - !--------------------------------- - if ( .not. cumulus(i,j) .or. kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low) & - ksurf_iterate(i,j) = .true. - if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & - kprof_cu == buoy_integ_low) ) & - dec_thres(i,j) = dec_thres_cu - ktop_iterate(i,j) = .true. + else ! ( .not. l_use_sml_dsc_fixes ) + ! If the bug-fix is off, we turn on ksurf_iterate at non-cumulus + ! points even if dsc was already true + ! (when all we've done is confirm the DSC shouldn't merge with the + ! SML, so we had no reason to be meddling with the SML). + ! This creates some inconsistencies down the line... + + !--------------------------------- + ! Extent of mixing must be reduced + !--------------------------------- + if ( .not. cumulus(i,j) .or. kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low) & + ksurf_iterate(i,j) = .true. + if ( cumulus(i,j) .and. ( kprof_cu == buoy_integ .or. & + kprof_cu == buoy_integ_low) ) & + dec_thres(i,j) = dec_thres_cu + ktop_iterate(i,j) = .true. - end if ! ( l_use_sml_dsc_fixes ) + end if ! ( l_use_sml_dsc_fixes ) - if ( .not. dsc(i,j) ) then - if (l_use_var_fixes .and. v_top(i,j) <= zero) then - ! If v_top=0 we have no mechanism to generate turbulence in a DSC - ktop_iterate(i,j) = .false. - z_top_lim(i,j) = zh(i,j) ! needs to be set here as zhsc not used - else - ! Set up a `COUPLED' decoupled layer, - ! implies no explicit `entrainment' at ZH. - ! Note a new ZH (and thence NTML) will be calculated by - ! wb integral iteration. - dsc(i,j) = .true. - coupled(i,j) = .true. - ntdsc(i,j) = ntml(i,j) - zhsc(i,j) = zh(i,j) - zc_dsc(i,j) = zc(i,j) - v_top_dsc(i,j) = v_top(i,j) - if ( l_use_sml_dsc_fixes ) then - ! Need to ensure SML-top entrainment flux is turned off when - ! using buoyancy-flux integration to find a new zh - v_top(i,j) = zero - ! Set coupling strength to maximal when just diagnosed decoupled - svl_diff_frac(i,j) = one - end if - v_sum_dsc(i,j) = v_sum(i,j) - ft_nt_zhsc(i,j) = ft_nt_zh(i,j) - fq_nt_zhsc(i,j) = fq_nt_zh(i,j) - if (l_wtrac) iset_wtrac(i,j) = 2 - ! put all entrainment into RHOKH_TOP - rhokh_dsct_ent(i,j) = rhokh_top_ent(i,j) & - + rhokh_surf_ent(i,j) - rhokh_top_ent(i,j) = zero - rhokh_surf_ent(i,j) = zero - rhokm_dsct_ent(i,j) = rhokm_top_ent(i,j) & - + rhokm_surf_ent(i,j) - rhokm_top_ent(i,j) = zero - rhokm_surf_ent(i,j) = zero - if (.not. l_use_var_fixes) then - rhokm_top(i,j,ntml(i,j)+1) = rhokm_top(i,j,ntml(i,j)+1) & - + rhokm(i,j,ntml(i,j)+1) - rhokm(i,j,ntml(i,j)+1) = zero - end if - end if ! test on l_use_var_fixes and v_top <= zero - end if ! test no not.dsc - end if ! test on WB_RATIO le DEC_THRES - end if ! testing for well-mixed layer (TEST_WELL_MIXED) + if ( .not. dsc(i,j) ) then + if (l_use_var_fixes .and. v_top(i,j) <= zero) then + ! If v_top=0 we have no mechanism to generate turbulence in a DSC + ktop_iterate(i,j) = .false. + z_top_lim(i,j) = zh(i,j) ! needs to be set here as zhsc not used + else + ! Set up a `COUPLED' decoupled layer, + ! implies no explicit `entrainment' at ZH. + ! Note a new ZH (and thence NTML) will be calculated by + ! wb integral iteration. + dsc(i,j) = .true. + coupled(i,j) = .true. + ntdsc(i,j) = ntml(i,j) + zhsc(i,j) = zh(i,j) + zc_dsc(i,j) = zc(i,j) + v_top_dsc(i,j) = v_top(i,j) + if ( l_use_sml_dsc_fixes ) then + ! Need to ensure SML-top entrainment flux is turned off when + ! using buoyancy-flux integration to find a new zh + v_top(i,j) = zero + ! Set coupling strength to maximal when just diagnosed decoupled + svl_diff_frac(i,j) = one + end if + v_sum_dsc(i,j) = v_sum(i,j) + ft_nt_zhsc(i,j) = ft_nt_zh(i,j) + fq_nt_zhsc(i,j) = fq_nt_zh(i,j) + if (l_wtrac) iset_wtrac(i,j) = 2 + ! put all entrainment into RHOKH_TOP + rhokh_dsct_ent(i,j) = rhokh_top_ent(i,j) & + + rhokh_surf_ent(i,j) + rhokh_top_ent(i,j) = zero + rhokh_surf_ent(i,j) = zero + rhokm_dsct_ent(i,j) = rhokm_top_ent(i,j) & + + rhokm_surf_ent(i,j) + rhokm_top_ent(i,j) = zero + rhokm_surf_ent(i,j) = zero + if (.not. l_use_var_fixes) then + rhokm_top(i,j,ntml(i,j)+1) = rhokm_top(i,j,ntml(i,j)+1) & + + rhokm(i,j,ntml(i,j)+1) + rhokm(i,j,ntml(i,j)+1) = zero + end if + end if ! test on l_use_var_fixes and v_top <= zero + end if ! test no not.dsc + end if ! test on WB_RATIO le DEC_THRES + end if ! testing for well-mixed layer (TEST_WELL_MIXED) + end do ! I end do ! I !end do ! J !$OMP end do @@ -1865,8 +1883,10 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end +!do i = pdims%i_start, pdims%i_end if ( ksurf_iterate(i,j) ) then !----------------------------------------------------------- @@ -1913,7 +1933,7 @@ subroutine excf_nl_9c ( & wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES end if ! KSURF_ITERATE - + end do end do !end do !$OMP end do @@ -1928,7 +1948,9 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) + !do i = pdims%i_start, pdims%i_end if ( cumulus(i,j) .and. .not. dsc(i,j) .and. & .not. ksurf_iterate(i,j)) then ! pure cumulus layer and necessary parameters not already set up @@ -1955,6 +1977,7 @@ subroutine excf_nl_9c ( & wb_ratio(i,j) = dec_thres(i,j) - one ! to be < DEC_THRES end if end do ! I +end do ! II !end do ! J !$OMP end do @@ -2242,149 +2265,154 @@ subroutine excf_nl_9c ( & !cdir collapse !$OMP do SCHEDULE(DYNAMIC) !do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - - if ( ktop_iterate(i,j) ) then - - ! Lower limit on base of DSC layer - z_bot_lim(i,j) = 0.1_r_bl * zh(i,j) - ! Upper limit on base of DSC layer - z_top_lim(i,j) = zhsc(i,j) - dzrad(i,j) - ! If cumulus limit base of top-driven mixing to above ZH - ! (not applying this limit if converging the gradient adjustment, - ! as in this case we need to allow kh_surf to go as high as it wants - ! so the DSC kh needs to be able to overlap with it). - if ( cumulus(i,j) .and. (.not. l_converge_ga) ) then - z_bot_lim(i,j) = min( zh(i,j), z_top_lim(i,j) ) - end if +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + + if ( ktop_iterate(i,j) ) then + + ! Lower limit on base of DSC layer + z_bot_lim(i,j) = 0.1_r_bl * zh(i,j) + ! Upper limit on base of DSC layer + z_top_lim(i,j) = zhsc(i,j) - dzrad(i,j) + ! If cumulus limit base of top-driven mixing to above ZH + ! (not applying this limit if converging the gradient adjustment, + ! as in this case we need to allow kh_surf to go as high as it wants + ! so the DSC kh needs to be able to overlap with it). + if ( cumulus(i,j) .and. (.not. l_converge_ga) ) then + z_bot_lim(i,j) = min( zh(i,j), z_top_lim(i,j) ) + end if - z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) + z_cbase(i,j) = zhsc(i,j) - zc_dsc(i,j) - !..Divide up depth of layer within which ZDSC_BASE is allowed - z_inc(i,j) = (z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) - zdsc_base(i,j) = z_bot_lim(i,j) - ! will start at Z_BOT_LIM+Z_INC + !..Divide up depth of layer within which ZDSC_BASE is allowed + z_inc(i,j) = (z_top_lim(i,j)-z_bot_lim(i,j)) / real(n_steps, r_bl) + zdsc_base(i,j) = z_bot_lim(i,j) + ! will start at Z_BOT_LIM+Z_INC - if (l_reset_dec_thres) then - dec_thres(i,j) = dec_thres_cloud ! reset after potential use - ! of dec_thres_cu or - ! dec_thres_clear above - end if - wb_ratio(i,j) = dec_thres(i,j) + one ! to be > DEC_THRES + if (l_reset_dec_thres) then + dec_thres(i,j) = dec_thres_cloud ! reset after potential use + ! of dec_thres_cu or + ! dec_thres_clear above + end if + wb_ratio(i,j) = dec_thres(i,j) + one ! to be > DEC_THRES - end if ! KTOP_ITERATE + end if ! KTOP_ITERATE + end do end do !$OMP end do !$OMP do SCHEDULE(DYNAMIC) -do i = pdims%i_start, pdims%i_end - if ( ktop_iterate(i,j) .and. wb_dzrad_int(i,j) < zero ) then - !------------------------------------------------------------- - ! Estimation of wb integral over radiatively cooled cloud-top - ! region not yet performed (ie. DSC over stable surface) so - ! do it now. - !------------------------------------------------------------- - z_inv(i,j) = zhsc(i,j) - z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) - cf_ml(i,j) = cf_dsc(i,j) - df_ctop(i,j)= df_dsct_over_cp(i,j) - rho_we(i,j) = rho_we_dsc(i,j) - dsl(i,j) = dsl_dsc(i,j) - dqw(i,j) = dqw_dsc(i,j) - tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) - totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) - ntop(i,j) = ntdsc(i,j) - dzrad(i,j) = 100.0_r_bl - if ( dzrad_disc_opt == dzrad_ntm1 ) then - ! Original discretisation of cloud-top radiatively-cooled layer; - ! set the base of the layer at theta-level ntdsc-1 - ! (which is between 1.5 and 2.5 model-levels below cloud-top, - ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), - ! then lower it by extra whole model-levels if needed to increase the - ! depth to 100m. - - ntop(i,j) = ntop(i,j) - 1 - if ( ntop(i,j) > 1 ) then - inner_loop2: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& - .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) - ntop(i,j) = ntop(i,j) - 1 - if (ntop(i,j) == 1 ) then - exit inner_loop2 - end if - end do inner_loop2 - end if - dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) - - else if ( dzrad_disc_opt == dzrad_1p5dz ) then - ! Alternative method; - ! set the base of the layer 1.5 model-levels below cloud-top so that - ! it always moves smoothly with zhsc instead of jumping suddenly - ! (1.5 levels is the minimum distance below zhsc which will always - ! keep the base of the layer below theta-level ntdsc, which is a - ! requirement of the discretisation). - ! The depth is also forced to be at least 100m. - - ! Find smoothly-varying height 1.5 model-levels below z_inv - k_inv = ntop(i,j) - if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then - interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & - / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) - z_pr = (one-interp) * z_tq(i,j,k_inv-1) & - + interp * z_uv(i,j,k_inv) - else - interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & - / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) - z_pr = (one-interp) * z_uv(i,j,k_inv) & - + interp * z_tq(i,j,k_inv) - end if - ! Go down to 100m below z_inv (or cloud-base) if that is lower, - ! but don't allow lower than theta-level 1 - z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) - z_pr = max( z_pr, z_tq(i,j,1) ) - ! Set ntop to the rho-level straddling the base of the layer - do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +! do i = pdims%i_start, pdims%i_end + if ( ktop_iterate(i,j) .and. wb_dzrad_int(i,j) < zero ) then + !------------------------------------------------------------- + ! Estimation of wb integral over radiatively cooled cloud-top + ! region not yet performed (ie. DSC over stable surface) so + ! do it now. + !------------------------------------------------------------- + z_inv(i,j) = zhsc(i,j) + z_cbase(i,j)= z_inv(i,j) - zc_dsc(i,j) + cf_ml(i,j) = cf_dsc(i,j) + df_ctop(i,j)= df_dsct_over_cp(i,j) + rho_we(i,j) = rho_we_dsc(i,j) + dsl(i,j) = dsl_dsc(i,j) + dqw(i,j) = dqw_dsc(i,j) + tothf_zi(i,j)= - rho_we(i,j)*dsl(i,j) + ft_nt_zhsc(i,j) + totqf_zi(i,j)= - rho_we(i,j)*dqw(i,j) + fq_nt_zhsc(i,j) + ntop(i,j) = ntdsc(i,j) + dzrad(i,j) = 100.0_r_bl + if ( dzrad_disc_opt == dzrad_ntm1 ) then + ! Original discretisation of cloud-top radiatively-cooled layer; + ! set the base of the layer at theta-level ntdsc-1 + ! (which is between 1.5 and 2.5 model-levels below cloud-top, + ! depending on where zhsc is between z_uv(ntdsc+1) and z_uv(ntdsc+2)), + ! then lower it by extra whole model-levels if needed to increase the + ! depth to 100m. + ntop(i,j) = ntop(i,j) - 1 - end do - ! Set new cloud-top radiatively-cooled layer-depth - dzrad(i,j) = z_inv(i,j) - z_pr - - ! Fraction of rho-level ntop below the dzrad layer - interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & - / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) - ! Add on SL and qw differences between the discontinuous inversion base - ! and the base of the dzrad layer - dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & - - interp *sl(i,j,ntop(i,j)) - dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & - - interp *qw(i,j,ntop(i,j)) - - end if ! ( dzrad_disc_opt ) - - wsl_dzrad_int = dzrad(i,j) * & - ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) - wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) - - wb_dzrad_int(i,j) = wsl_dzrad_int * ( & - (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & - + wqw_dzrad_int * ( & - (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & - cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) - wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) - wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) - - if (model_type == mt_single_column) then - do k = ntop(i,j)+1, ntdsc(i,j)+1 - wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) - end do - end if + if ( ntop(i,j) > 1 ) then + inner_loop2: do while ( z_tq(i,j,ntop(i,j)) > z_inv(i,j)-dzrad(i,j)& + .and. z_tq(i,j,ntop(i,j)-1) > z_inv(i,j)-z_cbase(i,j)) + ntop(i,j) = ntop(i,j) - 1 + if (ntop(i,j) == 1 ) then + exit inner_loop2 + end if + end do inner_loop2 + end if + dzrad(i,j) = z_inv(i,j) - z_tq(i,j,ntop(i,j)) + + else if ( dzrad_disc_opt == dzrad_1p5dz ) then + ! Alternative method; + ! set the base of the layer 1.5 model-levels below cloud-top so that + ! it always moves smoothly with zhsc instead of jumping suddenly + ! (1.5 levels is the minimum distance below zhsc which will always + ! keep the base of the layer below theta-level ntdsc, which is a + ! requirement of the discretisation). + ! The depth is also forced to be at least 100m. + + ! Find smoothly-varying height 1.5 model-levels below z_inv + k_inv = ntop(i,j) + if ( z_inv(i,j) < z_tq(i,j,k_inv+1) ) then + interp = ( z_inv(i,j) - z_uv(i,j,k_inv+1) ) & + / ( z_tq(i,j,k_inv+1) - z_uv(i,j,k_inv+1) ) + z_pr = (one-interp) * z_tq(i,j,k_inv-1) & + + interp * z_uv(i,j,k_inv) + else + interp = ( z_inv(i,j) - z_tq(i,j,k_inv+1) ) & + / ( z_uv(i,j,k_inv+2) - z_tq(i,j,k_inv+1) ) + z_pr = (one-interp) * z_uv(i,j,k_inv) & + + interp * z_tq(i,j,k_inv) + end if + ! Go down to 100m below z_inv (or cloud-base) if that is lower, + ! but don't allow lower than theta-level 1 + z_pr = min( z_pr, max( z_inv(i,j) - dzrad(i,j), z_cbase(i,j) ) ) + z_pr = max( z_pr, z_tq(i,j,1) ) + ! Set ntop to the rho-level straddling the base of the layer + do while ( z_tq(i,j,ntop(i,j)-1) > z_pr ) + ntop(i,j) = ntop(i,j) - 1 + end do + ! Set new cloud-top radiatively-cooled layer-depth + dzrad(i,j) = z_inv(i,j) - z_pr + + ! Fraction of rho-level ntop below the dzrad layer + interp = ( z_pr - z_tq(i,j,ntop(i,j)-1) ) & + / ( z_tq(i,j,ntop(i,j)) - z_tq(i,j,ntop(i,j)-1) ) + ! Add on SL and qw differences between the discontinuous inversion base + ! and the base of the dzrad layer + dsl(i,j) = dsl(i,j) + sl(i,j,k_inv) - (one-interp)*sl(i,j,ntop(i,j)-1) & + - interp *sl(i,j,ntop(i,j)) + dqw(i,j) = dqw(i,j) + qw(i,j,k_inv) - (one-interp)*qw(i,j,ntop(i,j)-1) & + - interp *qw(i,j,ntop(i,j)) + + end if ! ( dzrad_disc_opt ) + + wsl_dzrad_int = dzrad(i,j) * & + ( 0.66_r_bl*df_ctop(i,j) - rho_we(i,j)*dsl(i,j) ) + wqw_dzrad_int = - dzrad(i,j) * rho_we(i,j) * dqw(i,j) + + wb_dzrad_int(i,j) = wsl_dzrad_int * ( & + (one-cf_ml(i,j))*btm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*btm_cld(i,j,ntop(i,j)+1) ) & + + wqw_dzrad_int * ( & + (one-cf_ml(i,j))*bqm(i,j,ntop(i,j)+1) + & + cf_ml(i,j)*bqm_cld(i,j,ntop(i,j)+1) ) + wb_dzrad_int(i,j) = g * wb_dzrad_int(i,j) + wb_dzrad_int(i,j) = max( zero, wb_dzrad_int(i,j) ) - end if + if (model_type == mt_single_column) then + do k = ntop(i,j)+1, ntdsc(i,j)+1 + wbend(i,j,k) = wb_dzrad_int(i,j)/dzrad(i,j) + end do + end if - wb_dzrad_int(i,j) = max( rbl_eps, wb_dzrad_int(i,j) ) + end if -end do ! I + wb_dzrad_int(i,j) = max( rbl_eps, wb_dzrad_int(i,j) ) + end do ! I +end do ! II !end do ! J !$OMP end do @@ -2831,76 +2859,79 @@ subroutine excf_nl_9c ( & !$OMP do SCHEDULE(DYNAMIC) ! do j = pdims%j_start, pdims%j_end -do i = pdims%i_start, pdims%i_end - k=ntml(i,j)+1 - kh_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( & - rhokh_surf_ent(i,j) / & - ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) - if (l_use_var_fixes) then - km_top_factor(i,j) = max( 0.7_r_bl, one-sqrt( rhokm_surf_ent(i,j) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) - else - km_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( rhokm(i,j,k) / & - ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) - end if - scdepth(i,j) = zh(i,j) - zsml_base(i,j) - factor = g1 * rho_mix(i,j,k) * v_top(i,j) *vkman *scdepth(i,j) - if ( factor > zero) then - kh_sct_factor(i,j) = one - ( rhokh_top_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 - else - kh_sct_factor(i,j) = one - end if - factor = g1 * rho_wet_tq(i,j,k-1) * v_top(i,j) * & - vkman * scdepth(i,j) * 0.75_r_bl - if ( factor > zero) then +do ii = pdims%i_start, pdims%i_end, bl_segment_size + do i = ii, min(ii+bl_segment_size-1, pdims%i_end) +!do i = pdims%i_start, pdims%i_end + k=ntml(i,j)+1 + kh_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( & + rhokh_surf_ent(i,j) / & + ( rho_mix(i,j,k)*w_h_top(i,j)*vkman*zh(i,j) ) ) ) if (l_use_var_fixes) then - km_sct_factor(i,j) = one - & - ( rhokm_top_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 + km_top_factor(i,j) = max( 0.7_r_bl, one-sqrt( rhokm_surf_ent(i,j) / & + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) else - km_sct_factor(i,j) = one - & - ( rhokm_top(i,j,k) / factor )**1.25_r_bl + km_top_factor(i,j) = max( 0.7_r_bl , one - sqrt( rhokm(i,j,k) / & + ( rho_wet_tq(i,j,k-1)*w_m_top(i,j)*vkman*zh(i,j) ) ) ) end if - else - km_sct_factor(i,j) = one - end if - - if (ntdsc(i,j) > 0) then - !------------------------------------------------------------- - ! Set up factors to ensure K profile continuity at ZHSC; - ! no need to limit size of factor as precise shape of top-down - ! mixing profile not important. - ! Only calculate _DSCT_FACTORs when a decoupled stratocumulus - ! layer exists, i.e. NTDSC > 0. - !------------------------------------------------------------- - k=ntdsc(i,j)+1 - dscdepth(i,j) = zhsc(i,j) - zdsc_base(i,j) - factor = g1*rho_mix(i,j,k)*v_top_dsc(i,j)*vkman*dscdepth(i,j) + scdepth(i,j) = zh(i,j) - zsml_base(i,j) + factor = g1 * rho_mix(i,j,k) * v_top(i,j) *vkman *scdepth(i,j) if ( factor > zero) then - kh_dsct_factor(i,j) = one - & - ( rhokh_dsct_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 + kh_sct_factor(i,j) = one - ( rhokh_top_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 else - kh_dsct_factor(i,j) = one + kh_sct_factor(i,j) = one end if - - factor = 0.75_r_bl * g1 * rho_wet_tq(i,j,k-1) * v_top_dsc(i,j) * & - vkman * dscdepth(i,j) + factor = g1 * rho_wet_tq(i,j,k-1) * v_top(i,j) * & + vkman * scdepth(i,j) * 0.75_r_bl if ( factor > zero) then if (l_use_var_fixes) then - km_dsct_factor(i,j) = one - & - ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl - ! 1.25=1/0.8 + km_sct_factor(i,j) = one - & + ( rhokm_top_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 else - km_dsct_factor(i,j) = one - & + km_sct_factor(i,j) = one - & ( rhokm_top(i,j,k) / factor )**1.25_r_bl end if else - km_dsct_factor(i,j) = one + km_sct_factor(i,j) = one end if - end if -end do + + if (ntdsc(i,j) > 0) then + !------------------------------------------------------------- + ! Set up factors to ensure K profile continuity at ZHSC; + ! no need to limit size of factor as precise shape of top-down + ! mixing profile not important. + ! Only calculate _DSCT_FACTORs when a decoupled stratocumulus + ! layer exists, i.e. NTDSC > 0. + !------------------------------------------------------------- + k=ntdsc(i,j)+1 + dscdepth(i,j) = zhsc(i,j) - zdsc_base(i,j) + factor = g1*rho_mix(i,j,k)*v_top_dsc(i,j)*vkman*dscdepth(i,j) + if ( factor > zero) then + kh_dsct_factor(i,j) = one - & + ( rhokh_dsct_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 + else + kh_dsct_factor(i,j) = one + end if + + factor = 0.75_r_bl * g1 * rho_wet_tq(i,j,k-1) * v_top_dsc(i,j) * & + vkman * dscdepth(i,j) + if ( factor > zero) then + if (l_use_var_fixes) then + km_dsct_factor(i,j) = one - & + ( rhokm_dsct_ent(i,j) / factor )**1.25_r_bl + ! 1.25=1/0.8 + else + km_dsct_factor(i,j) = one - & + ( rhokm_top(i,j,k) / factor )**1.25_r_bl + end if + else + km_dsct_factor(i,j) = one + end if + end if + end do !I +end do !II ! end do !$OMP end do From 88deed071aadce24a74588798375a16c567bb1dd Mon Sep 17 00:00:00 2001 From: MetBenjaminWent Date: Wed, 4 Feb 2026 10:22:59 +0000 Subject: [PATCH 3/3] Remove groups to bring in line with trunk --- rose-stem/site/meto/groups/groups_lfric_atm.cylc | 5 ----- ..._atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 --------- ..._atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt | 9 --------- 3 files changed, 23 deletions(-) delete mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt delete mode 100644 rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index be3d7267a..174b5698f 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -222,11 +222,6 @@ "lfric_atm_clim_gal9_chem_1T-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_clim_gal9_chem_2T-C12_ex1a_cce_fast-debug-32bit", ], - "ex1a_omp_C48_cce_full": [ - "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit", - "lfric_atm_nwp_gal9_2T-C48_MG_ex1a_cce_full-debug-32bit", - "lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit", - ], "ex1a_omp_C48_cce": [ "lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_gal9_2T-C48_MG_ex1a_cce_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt deleted file mode 100644 index 7d0bbf9c4..000000000 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_1T-C48_MG_ex1a_cce_full-debug-32bit.txt +++ /dev/null @@ -1,9 +0,0 @@ -Inner product checksum rho = 48D7FD20 -Inner product checksum theta = 5392A6D8 -Inner product checksum u = 6A97B848 -Inner product checksum mr1 = 41CD0DCE -Inner product checksum mr2 = 39CF631F -Inner product checksum mr3 = 37AD770F -Inner product checksum mr4 = 39604773 -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_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt deleted file mode 100644 index 41334f2a4..000000000 --- a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_nwp_gal9_4T-C48_MG_ex1a_cce_full-debug-32bit.txt +++ /dev/null @@ -1,9 +0,0 @@ -Inner product checksum rho = 48D7FD32 -Inner product checksum theta = 5392A6F0 -Inner product checksum u = 6A97B6F0 -Inner product checksum mr1 = 41CD0826 -Inner product checksum mr2 = 39C92715 -Inner product checksum mr3 = 37B19C28 -Inner product checksum mr4 = 396070D2 -Inner product checksum mr5 = 0 -Inner product checksum mr6 = 0