diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 50ca66db0..4f91fde35 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -7,4 +7,5 @@ | mo-marqh | mark Hedley | Met Office | 2025-12-11 | | yaswant | Yaswant Pradhan | Met Office | 2025-12-16 | | oakleybrunt | Oakley Brunt | Met Office | 2025-12-19 | -| harry-shepherd | Harry Shepherd | Met Office | 2026-01-08 | \ No newline at end of file +| harry-shepherd | Harry Shepherd | Met Office | 2026-01-08 | +| Adrian-Lock | Adrian Lock | Met Office | 2026-01-09 | diff --git a/applications/lfric_atm/optimisation/meto-ex1a/transmute/kernel/jules_exp_kernel_mod.py b/applications/lfric_atm/optimisation/meto-ex1a/transmute/kernel/jules_exp_kernel_mod.py index 7c28ef7f4..055c52370 100644 --- a/applications/lfric_atm/optimisation/meto-ex1a/transmute/kernel/jules_exp_kernel_mod.py +++ b/applications/lfric_atm/optimisation/meto-ex1a/transmute/kernel/jules_exp_kernel_mod.py @@ -36,7 +36,6 @@ def trans(psyir): "sea_ice_pensolar", "rhostar_2d", "recip_l_mo_sea_2d", - "h_blend_orog_2d", "t1_sd_2d", "q1_sd_2d", "surf_interp", diff --git a/applications/ngarch/source/algorithm/boundary_layer_timestep_mod.x90 b/applications/ngarch/source/algorithm/boundary_layer_timestep_mod.x90 index 27cac265b..e2ba03076 100644 --- a/applications/ngarch/source/algorithm/boundary_layer_timestep_mod.x90 +++ b/applications/ngarch/source/algorithm/boundary_layer_timestep_mod.x90 @@ -57,7 +57,7 @@ contains water_extraction, lake_evap, & theta_star_surf, qv_star_surf, & recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d + t1_sd_2d, q1_sd_2d type( field_array_type ), pointer :: mr type(field_collection_type), pointer :: prognostic_fields, & @@ -121,7 +121,7 @@ contains turbulence_fields, convection_fields, cloud_fields, & surface_fields, soil_fields, snow_fields, & aerosol_fields, recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d ) + t1_sd_2d, q1_sd_2d ) ! Call the algorithm call log_event( "Running BOUNDARY LAYER EXP", LOG_LEVEL_INFO ) @@ -136,7 +136,7 @@ contains cloud_fields, & surface_fields, & recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d, & + t1_sd_2d, q1_sd_2d, & modeldb%clock ) call log_event( "Running BOUNDARY LAYER IMP", LOG_LEVEL_INFO ) diff --git a/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 b/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 index 18e785d09..a1d51d235 100644 --- a/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 +++ b/interfaces/jules_interface/source/algorithm/jules_exp_alg_mod.x90 @@ -64,7 +64,6 @@ contains !>@param[in] aerosol_fields Fields for aerosol scheme !>@param[in,out] recip_l_mo_sea Inverse Obukhov length over sea only !>@param[in,out] rhostar Surface density - !>@param[in,out] h_blend_orog Orographic blending height !>@param[in,out] t1_sd_2d StDev of level 1 temperature !>@param[in,out] q1_sd_2d StDev of level 1 humidity subroutine jules_exp_alg(modeldb, theta, mr_n, & @@ -73,7 +72,7 @@ contains turbulence_fields, convection_fields, & cloud_fields, surface_fields, soil_fields, & snow_fields, aerosol_fields, recip_l_mo_sea, & - rhostar, h_blend_orog, t1_sd_2d, q1_sd_2d) + rhostar, t1_sd_2d, q1_sd_2d) use psykal_lite_phys_mod, only: invoke_jules_exp_kernel_type use xios, only: xios_date, xios_get_current_date, & @@ -98,7 +97,6 @@ contains type( field_type ), intent( inout ) :: recip_l_mo_sea type( field_type ), intent( inout ) :: rhostar - type( field_type ), intent( inout ) :: h_blend_orog type( field_type ), intent( inout ) :: t1_sd_2d type( field_type ), intent( inout ) :: q1_sd_2d @@ -376,7 +374,6 @@ contains call zh%copy_field_properties(recip_l_mo_sea) call zh%copy_field_properties(rhostar) - call zh%copy_field_properties(h_blend_orog) call zh%copy_field_properties(t1_sd_2d) call zh%copy_field_properties(q1_sd_2d) @@ -448,7 +445,7 @@ contains flux_e, flux_h, urbwrr, & urbhwr, urbhgt, urbztm, urbdisp, & rhostar, recip_l_mo_sea, & - h_blend_orog, t1_sd_2d, q1_sd_2d, & + t1_sd_2d, q1_sd_2d, & gross_prim_prod, z0h_eff, ocn_cpl_point, & stencil_depth) diff --git a/interfaces/jules_interface/source/algorithm/jules_timestep_alg_mod.x90 b/interfaces/jules_interface/source/algorithm/jules_timestep_alg_mod.x90 index f529642ac..50a9db08c 100644 --- a/interfaces/jules_interface/source/algorithm/jules_timestep_alg_mod.x90 +++ b/interfaces/jules_interface/source/algorithm/jules_timestep_alg_mod.x90 @@ -218,7 +218,7 @@ contains ! For jules_exp_alg ! not required for standalone, but required for boundary layer - type( field_type ) :: recip_l_mo_sea, rhostar, h_blend_orog + type( field_type ) :: recip_l_mo_sea, rhostar type( field_type ) :: t1_sd_2d, q1_sd_2d ! For jules_imp_alg @@ -288,7 +288,7 @@ contains turbulence_fields, convection_fields, cloud_fields, & surface_fields, soil_fields, snow_fields, & aerosol_fields, recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d) + t1_sd_2d, q1_sd_2d) ! Fields set up for jules_imp ! ncells is set to 1 after jules_exp so needs to be reset diff --git a/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 b/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 index 27f10473f..24e2aceb5 100644 --- a/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 +++ b/interfaces/jules_interface/source/kernel/jules_exp_kernel_mod.F90 @@ -46,7 +46,7 @@ module jules_exp_kernel_mod !> type, public, extends(kernel_type) :: jules_exp_kernel_type private - type(arg_type) :: meta_args(109) = (/ & + type(arg_type) :: meta_args(108) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! theta_in_wth arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! exner_in_wth arg_type(GH_FIELD, GH_REAL, GH_READ, W3, STENCIL(REGION)), &! u_in_w3 @@ -150,7 +150,6 @@ module jules_exp_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), &! urbdisp arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! rhostar arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! recip_l_mo_sea - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! h_blend_orog arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! t1_sd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! q1_sd arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1),&! diag__gross_prim_prod @@ -273,7 +272,6 @@ module jules_exp_kernel_mod !> @param[in] urbdisp Urban displacement height !> @param[in,out] rhostar_2d Surface density !> @param[in,out] recip_l_mo_sea_2d Inverse Obukhov length over sea only - !> @param[in,out] h_blend_orog_2d Orographic blending height !> @param[in,out] t1_sd_2d StDev of level 1 temperature !> @param[in,out] q1_sd_2d StDev of level 1 humidity !> @param[in,out] gross_prim_prod Diagnostic: Gross Primary Productivity @@ -424,7 +422,6 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & urbdisp, & rhostar_2d, & recip_l_mo_sea_2d, & - h_blend_orog_2d, & t1_sd_2d, & q1_sd_2d, & gross_prim_prod, & @@ -603,7 +600,6 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & z0m_eff, & ustar, & soil_moist_avail, & - h_blend_orog_2d, & recip_l_mo_sea_2d, & rhostar_2d, & t1_sd_2d, q1_sd_2d @@ -1678,7 +1674,6 @@ subroutine jules_exp_code(nlayers, seg_len, seg_len_halo, & ! variables passed to explicit BL rhostar_2d(map_2d(1,i)) = rhostar(i,1) recip_l_mo_sea_2d(map_2d(1,i)) = recip_l_mo_sea(i,1) - h_blend_orog_2d(map_2d(1,i)) = jules_vars%h_blend_orog_ij(i,1) t1_sd_2d(map_2d(1,i)) = t1_sd(i,1) q1_sd_2d(map_2d(1,i)) = q1_sd(i,1) diff --git a/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 index a514beb49..642a0b8dd 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90 @@ -67,7 +67,6 @@ contains !>@param[in,out] surface_fields Fields for surface scheme !>@param[in] recip_l_mo_sea Inverse Obukhov length over sea only !>@param[in] rhostar Surface density - !>@param[in h_blend_orog Orographic blending height !>@param[in] t1_sd_2d StDev of level 1 temperature !>@param[in] q1_sd_2d StDev of level 1 humidity !>@param[in] model_clock Time in the model @@ -76,8 +75,7 @@ contains microphysics_fields, dmr_mphys, orography_fields, & turbulence_fields, convection_fields, & cloud_fields, surface_fields, & - recip_l_mo_sea, & - rhostar, h_blend_orog, t1_sd_2d, q1_sd_2d, & + recip_l_mo_sea, rhostar, t1_sd_2d, q1_sd_2d, & model_clock) use bl_exp_kernel_mod, only: bl_exp_kernel_type @@ -101,7 +99,6 @@ contains type( field_type ), intent( in ) :: recip_l_mo_sea type( field_type ), intent( in ) :: rhostar - type( field_type ), intent( in ) :: h_blend_orog type( field_type ), intent( in ) :: t1_sd_2d type( field_type ), intent( in ) :: q1_sd_2d @@ -380,7 +377,7 @@ contains cumulus, tile_fraction, sd_orog, & peak_to_trough_orog, silhouette_area_orog, & tile_temperature, rhostar, recip_l_mo_sea, & - h_blend_orog, t1_sd_2d, q1_sd_2d, & + t1_sd_2d, q1_sd_2d, & dtl_mphys, dmt_mphys, sw_heating_rate, & lw_heating_rate, cf_bulk, cf_liquid, & rh_crit, tnuc, tnuc_nlcl, & diff --git a/interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90 index f5ab6605d..c63bc632f 100644 --- a/interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90 @@ -40,7 +40,7 @@ module bl_exp_kernel_mod !> type, public, extends(kernel_type) :: bl_exp_kernel_type private - type(arg_type) :: meta_args(94) = (/ & + type(arg_type) :: meta_args(93) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! theta_in_wth arg_type(GH_FIELD, GH_REAL, GH_READ, W3), &! rho_in_w3 arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! rho_in_wth @@ -72,7 +72,6 @@ module bl_exp_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2),&! tile_temperature arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! rhostar arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! recip_l_mo_sea - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! h_blend_orog arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! t1_sd arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1),&! q1_sd arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &! dtl_mphys @@ -181,7 +180,6 @@ module bl_exp_kernel_mod !> @param[in] tile_temperature Surface tile temperatures !> @param[in] rhostar_2d Surface density !> @param[in] recip_l_mo_sea_2d Inverse Obukhov length over sea only - !> @param[in] h_blend_orog_2d Orographic blending height !> @param[in] t1_sd_2d StDev of level 1 temperature !> @param[in] q1_sd_2d StDev of level 1 humidity !> @param[in] dtl_mphys Microphysics liq temperature increment @@ -297,7 +295,6 @@ subroutine bl_exp_code(nlayers, seg_len, & tile_temperature, & rhostar_2d, & recip_l_mo_sea_2d, & - h_blend_orog_2d, & t1_sd_2d, & q1_sd_2d, & dtl_mphys, & @@ -479,8 +476,7 @@ subroutine bl_exp_code(nlayers, seg_len, & cumulus_2d, & shallow_flag, & level_parcel_top - real(kind=r_def), dimension(undf_2d), intent(in) :: h_blend_orog_2d, & - recip_l_mo_sea_2d, & + real(kind=r_def), dimension(undf_2d), intent(in) :: recip_l_mo_sea_2d, & rhostar_2d, & t1_sd_2d, q1_sd_2d, & sea_u_current, & @@ -549,7 +545,7 @@ subroutine bl_exp_code(nlayers, seg_len, & zh, dzh, wstar, wthvs, u_0_p, v_0_p, zlcl_uv, qsat_lcl, delthvu, & bl_type_1, bl_type_2, bl_type_3, bl_type_4, bl_type_5, bl_type_6, & bl_type_7, uw0, vw0, zhnl, rhostar, & - h_blend_orog, recip_l_mo_sea, flandg, t1_sd, q1_sd, qcl_inv_top, & + recip_l_mo_sea, flandg, t1_sd, q1_sd, qcl_inv_top, & fb_surf, rib_gb, z0m_eff_gb, zhsc, ustargbm, cos_theta_latitude, & max_diff, delta_smag, tnuc_nlcl_um real(r_um), dimension(seg_len,1) :: surf_dep_flux, zeroes @@ -652,7 +648,6 @@ subroutine bl_exp_code(nlayers, seg_len, & ustargbm(i,1) = ustar(map_2d(1,i)) rhostar(i,1) = rhostar_2d(map_2d(1,i)) recip_l_mo_sea(i,1) = recip_l_mo_sea_2d(map_2d(1,i)) - h_blend_orog(i,1) = h_blend_orog_2d(map_2d(1,i)) rib_gb(i,1) = gradrinr(map_wth(1,i)) z0m_eff_gb(i,1) = z0m_eff(map_2d(1,i)) ftl(i,1,1) = heat_flux_bl(map_w3(1,i)) @@ -895,8 +890,8 @@ subroutine bl_exp_code(nlayers, seg_len, & ! IN cloud/moisture data : bulk_cloud_fraction,q,qcf,qcl,temperature,qw,tl, & ! IN everything not covered so far : - rad_hr,micro_tends,fb_surf,ustargbm,p_star,tstar,h_blend_orog, & - zh_prev, zhpar,zlcl,ho2r2_orog_gb,sd_orog,wtrac_as, & + rad_hr,micro_tends,fb_surf,ustargbm,p_star,tstar, & + zh_prev, zhpar,zlcl,ho2r2_orog_gb,sd_orog,wtrac_as, & ! 2 IN 3 INOUT for Smagorinsky delta_smag, max_diff, rneutml_sq, visc_m, visc_h, & ! SCM Diagnostics (dummy values in full UM) & stash diag diff --git a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 index 71e723b3a..a064523c8 100644 --- a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 +++ b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 @@ -338,7 +338,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &soil_moist_avail, snow_unload_rate, albedo_obs_scaling, soil_clay, soil_sand, dust_mrel, dust_flux, day_of_year, second_of_day, & flux_e, flux_h, urbwrr, urbhwr, urbhgt, urbztm, urbdisp, & &rhostar, recip_l_mo_sea, & -&h_blend_orog, t1_sd_2d, q1_sd_2d, gross_prim_prod, z0h_eff, ocn_cpl_point, stencil_depth) +&t1_sd_2d, q1_sd_2d, gross_prim_prod, z0h_eff, ocn_cpl_point, stencil_depth) USE jules_exp_kernel_mod, ONLY: jules_exp_code USE mesh_mod, ONLY: mesh_type USE stencil_dofmap_mod, ONLY: STENCIL_REGION @@ -358,7 +358,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &chr1p5m_tile, resfs_tile, gc_tile, canhc_tile, tile_water_extract, z0m_eff, ustar, soil_moist_avail, snow_unload_rate, & &albedo_obs_scaling, soil_clay, soil_sand, dust_mrel, dust_flux, & urbwrr, urbhwr, urbhgt, urbztm, urbdisp, & -rhostar, recip_l_mo_sea, h_blend_orog, t1_sd_2d, q1_sd_2d, & +rhostar, recip_l_mo_sea, t1_sd_2d, q1_sd_2d, & &gross_prim_prod, z0h_eff TYPE(integer_field_type), intent(in) :: n_snow_layers, blend_height_tq, ocn_cpl_point INTEGER(KIND=i_def), intent(in) :: stencil_depth, ncells, ncells_halo, day_of_year, second_of_day @@ -385,7 +385,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth &tile_water_extract_proxy, z0m_eff_proxy, ustar_proxy, soil_moist_avail_proxy, snow_unload_rate_proxy, albedo_obs_scaling_proxy, & &soil_clay_proxy, soil_sand_proxy, dust_mrel_proxy, dust_flux_proxy, & urbwrr_proxy, urbhwr_proxy, urbhgt_proxy, urbztm_proxy, urbdisp_proxy, & -rhostar_proxy, recip_l_mo_sea_proxy, h_blend_orog_proxy, & +rhostar_proxy, recip_l_mo_sea_proxy, & &t1_sd_2d_proxy, q1_sd_2d_proxy, gross_prim_prod_proxy, z0h_eff_proxy INTEGER(KIND=i_def), pointer :: map_adspc10_dust_mrel(:,:) => null(), map_adspc1_zh(:,:) => null(), & &map_adspc2_tile_fraction(:,:) => null(), map_adspc3_leaf_area_index(:,:) => null(), & @@ -515,7 +515,6 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth urbdisp_proxy = urbdisp%get_proxy() rhostar_proxy = rhostar%get_proxy() recip_l_mo_sea_proxy = recip_l_mo_sea%get_proxy() - h_blend_orog_proxy = h_blend_orog%get_proxy() t1_sd_2d_proxy = t1_sd_2d%get_proxy() q1_sd_2d_proxy = q1_sd_2d%get_proxy() gross_prim_prod_proxy = gross_prim_prod%get_proxy() @@ -673,7 +672,7 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth urbwrr_proxy%data, urbhwr_proxy%data, urbhgt_proxy%data, urbztm_proxy%data, & urbdisp_proxy%data, & rhostar_proxy%data, recip_l_mo_sea_proxy%data, & -&h_blend_orog_proxy%data, t1_sd_2d_proxy%data, q1_sd_2d_proxy%data, gross_prim_prod_proxy%data, & +&t1_sd_2d_proxy%data, q1_sd_2d_proxy%data, gross_prim_prod_proxy%data, & z0h_eff_proxy%data, ocn_cpl_point_proxy%data, ndf_wtheta, & &undf_wtheta, map_wtheta, ndf_w3, undf_w3, map_w3, ndf_adspc1_zh, undf_adspc1_zh, map_adspc1_zh, & &ndf_adspc2_tile_fraction, undf_adspc2_tile_fraction, map_adspc2_tile_fraction, ndf_adspc3_leaf_area_index, & @@ -723,7 +722,6 @@ SUBROUTINE invoke_jules_exp_kernel_type(ncells, ncells_halo, theta, exner_in_wth CALL dust_flux_proxy%set_dirty() CALL rhostar_proxy%set_dirty() CALL recip_l_mo_sea_proxy%set_dirty() - CALL h_blend_orog_proxy%set_dirty() CALL t1_sd_2d_proxy%set_dirty() CALL q1_sd_2d_proxy%set_dirty() CALL gross_prim_prod_proxy%set_dirty() diff --git a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 index 191a42bce..b1ca9c1fe 100644 --- a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 +++ b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 @@ -338,9 +338,8 @@ subroutine um_physics_init() lem_adjust, interactive_fluxes, specified_fluxes_only, & except_disc_inv, ntml_level_corrn, free_trop_layers, sharpest, & lem_stability, sg_shear_enh_lambda, l_new_kcloudtop, buoy_integ, & - l_reset_dec_thres, DynDiag_ZL_CuOnly, var_diags_opt, & - i_interp_local, i_interp_local_gradients, & - split_tke_and_inv, l_noice_in_turb, l_use_var_fixes, & + l_reset_dec_thres, DynDiag_ZL_CuOnly, i_interp_local, & + i_interp_local_gradients, l_noice_in_turb, l_use_var_fixes, & i_interp_local_cf_dbdz, tke_diag_fac, a_ent_2, dec_thres_cloud, & dec_thres_cu, near_neut_z_on_l, blend_gridindep_fa, & specified_fluxes_tstar, buoy_integ_low, num_sweeps_bflux, & @@ -766,10 +765,9 @@ subroutine um_physics_init() sg_orog_mixing = sg_shear_enh_lambda end select - ! Switch for alternative TKE and variance diagnostics - var_diags_opt = split_tke_and_inv - tke_diag_fac = 1.0_r_bl + ! Switch for corrections to variance diagnostics l_use_var_fixes = .true. + tke_diag_fac = 1.0_r_bl zhloc_depth_fac = real(zhloc_depth_fac_in, r_bl) if (topography == topography_horizon) then diff --git a/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 b/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 index f2bbe547f..bb629c95c 100644 --- a/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 +++ b/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 @@ -274,7 +274,7 @@ contains type( field_type ) :: conv_liquid_fraction, conv_frozen_fraction type( field_type ) :: conv_liquid_mmr, conv_frozen_mmr, conv_frozen_number ! ...for JULES & Boundary Layer - type( field_type ) :: recip_l_mo_sea, rhostar, h_blend_orog + type( field_type ) :: recip_l_mo_sea, rhostar type( field_type ) :: t1_sd_2d, q1_sd_2d ! ...for spectral GWD type( field_type ) :: du_spectral_gwd, dv_spectral_gwd @@ -848,14 +848,14 @@ contains turbulence_fields, convection_fields, cloud_fields, & surface_fields, soil_fields, snow_fields, & aerosol_fields, recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d) + t1_sd_2d, q1_sd_2d) call bl_exp_alg(theta, rho, exner, mr_n, & derived_fields, radiation_fields, & microphysics_fields, dmr_mphys, orography_fields, & turbulence_fields, convection_fields, cloud_fields, & surface_fields, & recip_l_mo_sea, rhostar, & - h_blend_orog, t1_sd_2d, q1_sd_2d, clock) + t1_sd_2d, q1_sd_2d, clock) end if !-------------------------------------------------------------------- diff --git a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 index af0136c06..557125c8a 100644 --- a/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 +++ b/science/physics_schemes/source/boundary_layer/bdy_expl2.F90 @@ -46,7 +46,7 @@ subroutine bdy_expl2 ( & ! in cloud/moisture data : cf_bulk,q,qcf,qcl,t,qw,tl, & ! in everything not covered so far : - rad_hr,micro_tends,fb_surf,u_s,pstar,tstar,h_blend_orog, & + rad_hr,micro_tends,fb_surf,u_s,pstar,tstar, & zh_prev,zhpar,z_lcl,ho2r2_orog,sd_orog,wtrac_as, & ! 2 in 3 INOUT for Smagorinsky delta_smag, max_diff, rneutml_sq, visc_m, visc_h, & @@ -77,11 +77,11 @@ subroutine bdy_expl2 ( & off, max_t_grad, a_grad_adj, sg_orog_mixing, l_use_surf_in_ri, & h_scale, t_drain, idyndiag, DynDiag_ZL, l_noice_in_turb, & DynDiag_ZL_corrn, DynDiag_ZL_CuOnly, DynDiag_Ribased, & - RiCrit_sharp, zhloc_depth_fac, non_local_bl, on, l_full_lambdas, & + RiCrit_sharp, zhloc_depth_fac, non_local_bl, on, & nl_bl_levels, local_fa, free_trop_layers, to_sharp_across_1km, & sbl_op, equilibrium_sbl, one_third, two_thirds, blending_option, & blend_except_cu, blend_cth_shcu_only, sg_shear, sg_shear_enh_lambda, & - max_tke, var_diags_opt, original_vars, split_tke_and_inv, tke_diag_fac, & + max_tke, tke_diag_fac, & i_interp_local, i_interp_local_gradients, i_interp_local_cf_dbdz, & shallow_cu_maxtop, sc_cftol, near_neut_z_on_l, zero, one, one_half use cloud_inputs_mod, only: i_rhcpt, forced_cu, i_cld_vn, i_pc2_init_method, & @@ -257,9 +257,6 @@ subroutine bdy_expl2 ( & ! in Surface pressure (Pascals). tstar(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end), & ! in Surface temperature (K). - h_blend_orog(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end), & - ! in Blending height used as part - ! of effective roughness scheme zh_prev(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end), & ! in boundary layer height from ! previous timestep @@ -808,12 +805,9 @@ subroutine bdy_expl2 ( & real(kind=r_bl), parameter :: max_abs_obkhov = 1.0e6_r_bl ! Maximum permitted magnitude of the Obukhov ! length (m). -real(kind=r_bl), parameter :: small_sh = 0.01_r_bl - ! Minimum value of sh used in - ! original_vars variance diagnostics real(kind=r_bl), parameter :: small_tke = 1.0e-6_r_bl ! Minimum required value of TKE before - ! split_tke_and_inv variance diagnostics are calculated + ! variance diagnostics are calculated real(kind=r_bl), parameter :: max_ri = 0.01_r_bl*sqrt(huge(one)) ! Maximum (absolute) Richardson number which ensures that ! the stability functions (~ri^2) remain real-valued at @@ -957,11 +951,10 @@ subroutine bdy_expl2 ( & ! heterogeneous land surface this is poorly defined and we can't use Rib ! from the surface scheme as vertically averaging Ri is numerically ! unstable. So, over land, only the average temperature gradient is used -if ( .not. l_use_surf_in_ri .and. (var_diags_opt > original_vars .or. & - i_interp_local == i_interp_local_cf_dbdz ) ) then +if (.not. l_use_surf_in_ri) then ! if not using surface variables in Ri (l_use_surf_in_ri=false) we - ! extrapolate dbdz itself from level 2, but the sl and qw gradients are - ! used in the newer variance calculations and with i_interp_local_cf_dbdz + ! extrapolate dbdz itself from level 2, with the sl and qw gradients being + ! used in the variance calculations and with i_interp_local_cf_dbdz ! so extrapolate them here !$OMP do SCHEDULE(STATIC) do j = pdims%j_start, pdims%j_end @@ -2002,11 +1995,10 @@ subroutine bdy_expl2 ( & !----------------------------------------------------------------------- call ex_coef ( & ! in levels/logicals - bl_levels,k_log_layr,nSCMDpkgs,L_SCMDiags,BL_diag, & + bl_levels,k_log_layr,BL_diag, & ! in fields - sigma_h,flandg,dbdz,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m_eff_gb, & - h_blend_orog,zhpar,ntpar,ntml_nl,ntdsc,nbdsc,u_p,v_p,u_s,fb_surf, & - qw,tl,l_shallow_cth,rmlmax2, rneutml_sq, delta_smag, & + sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m_eff_gb,zhpar,ntpar, & + ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq,delta_smag, & ! in/out fields cumulus,weight_1dbl, & ! out fields @@ -2067,27 +2059,6 @@ subroutine bdy_expl2 ( & ! Include mixing length, ELH, in RHOKH. ! Code moved from EX_COEF to avoid interpolation !------------------------------------------------------------ - if ( k >= ntml_local(i,j)+2 .and. l_full_lambdas .and. & - local_fa == to_sharp_across_1km ) then - ! Assuming only LOCAL_FA = "to_sharp_across_1km" option - ! will have L_FULL_LAMBDAS. - ! If other LOCAL_FA options are coded here then - ! changes must be included in section 2.1 of ex_coef - if (l_rp2) then - lambdah = max ( lambda_min , & - par_mezcla_rp(rp_idx)*zh_local(i,j) ) - else - lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) - end if - z_scale = 1000.0_r_bl - weight1 = one_half*( one - & - tanh(3.0_r_bl*((z_uv(i,j,k)/z_scale )-one) ) ) - lambdah = lambdah * weight1 & - + lambda_min*( one - weight1) - ! no need to do log profile correction as klog_layr eq 2 - vkz = vkman * ( z_uv(i,j,k) + z0m_eff_gb(i,j) ) - elh_rho(i,j,k) = vkz / (one + vkz/lambdah ) - end if ! Reinstate UKV drainage flow bug here, where lambdah was not ! enhanced as intended (and as was done in ex_coef)! if (sg_orog_mixing == sg_shear_enh_lambda) then @@ -2097,7 +2068,7 @@ subroutine bdy_expl2 ( & else lambdah = max ( lambda_min , 0.15_r_bl*zh_local(i,j) ) end if - if (k >= ntml_local(i,j)+2 .and. .not. l_full_lambdas) then + if (k >= ntml_local(i,j)+2) then lambdah = lambda_min end if if (k <= k_log_layr) then @@ -2226,65 +2197,37 @@ subroutine bdy_expl2 ( & ! function of ustar, wstar) but is currently just set to zero if (BL_diag%l_tke) then - if (var_diags_opt == original_vars) then - ! BL_diag%tke currently contains (the reciprocal of) the non-local - ! (SML and DSC) mixed layer timescale, calculated in excf_nl - -!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED(BL_diag, rho_wet_tq, rhokm, dbdz, bl_levels, pdims) & -!$OMP private(i, j, k, recip_time_sbl, recip_time_cbl) - do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - - ! stable timescale - recip_time_sbl = sqrt( max(dbdz(i,j,k),1.0e-5_r_bl) )/0.7_r_bl - recip_time_cbl = BL_diag%tke(i,j,k) - - ! TKE diagnostic - taking 5 m2/s2 as a suitable maximum - BL_diag%tke(i,j,k)= min( max_tke, & - ( rhokm(i,j,k) / rho_wet_tq(i,j,k-1) )*( & - recip_time_cbl + recip_time_sbl ) ) - - end do - end do - end do -!$OMP end PARALLEL do - - else if (var_diags_opt == split_tke_and_inv) then - ! Combine the, separately calculated, local and non-local TKE diagnostics + ! Combine the, separately calculated, local and non-local TKE diagnostics !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED(BL_diag, tke_nl, tke_loc, rho_wet_tq, weight_1dbl, & !$OMP tke_diag_fac, bl_levels, pdims) & !$OMP private(i, j, k) - do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end + do k = 2, bl_levels + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end - ! TKE diagnostic - ! Currently tke_nl is strictly rho*sigma_w^2 while tke_loc is TKE - ! Assume isotropic turb so TKE_nl = 3/2 sigma_w^2 - tke_nl(i,j,k) = weight_1dbl(i,j,k) * 1.5_r_bl * & - tke_nl(i,j,k)/rho_wet_tq(i,j,k-1) - BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,k) ) + ! TKE diagnostic + ! Currently tke_nl is strictly rho*sigma_w^2 while tke_loc is TKE + ! Assume isotropic turb so TKE_nl = 3/2 sigma_w^2 + tke_nl(i,j,k) = weight_1dbl(i,j,k) * 1.5_r_bl * & + tke_nl(i,j,k)/rho_wet_tq(i,j,k-1) + BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,k) ) - ! Multiply by tuning factor - BL_diag%tke(i,j,k) = tke_diag_fac * BL_diag%tke(i,j,k) + ! Multiply by tuning factor + BL_diag%tke(i,j,k) = tke_diag_fac * BL_diag%tke(i,j,k) - ! Keep TKE below a sensible max value of max_tke - BL_diag%tke(i,j,k) = min( max_tke, BL_diag%tke(i,j,k) ) - ! Applying this limit can occasionally cause the length-scale - ! Km / sqrt(w_var) to become unrealistically large, since no - ! equivalent limiting is done on Km. + ! Keep TKE below a sensible max value of max_tke + BL_diag%tke(i,j,k) = min( max_tke, BL_diag%tke(i,j,k) ) + ! Applying this limit can occasionally cause the length-scale + ! Km / sqrt(w_var) to become unrealistically large, since no + ! equivalent limiting is done on Km. - end do end do end do + end do !$OMP end PARALLEL do - end if ! var_diags_opt - if ( i_bm_ez_opt == i_bm_ez_entpar ) then ! Calculate mixing-length to pass to bimodal cloud scheme, ! using Km and TKE. @@ -2362,7 +2305,7 @@ subroutine bdy_expl2 ( & ! to the microphysics turbulence call !$OMP PARALLEL DEFAULT(none) private(k,j,i) & -!$OMP SHARED(tdims, bl_levels,pdims,BL_diag,bl_w_var,var_diags_opt) +!$OMP SHARED(tdims, bl_levels,pdims,BL_diag,bl_w_var) !$OMP do SCHEDULE(STATIC) do k = 2, tdims%k_end+1 @@ -2374,73 +2317,41 @@ subroutine bdy_expl2 ( & end do !$OMP end do - if (var_diags_opt == original_vars) then - ! Original version: w_var = TKE + ! w_var = 2/3 TKE !$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 (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then - bl_w_var(i,j,k) = BL_diag%tke(i,j,k) - end if - end do - end do - end do -!$OMP end do - else if (var_diags_opt == split_tke_and_inv) then - ! Improved version: w_var = 2/3 TKE -!$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 (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then - bl_w_var(i,j,k) = two_thirds * BL_diag%tke(i,j,k) - end if - end do + do k = 2, bl_levels + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (BL_diag%tke(i,j,k) > 1.0e-12_r_bl) then + bl_w_var(i,j,k) = two_thirds * BL_diag%tke(i,j,k) + end if end do end do + end do !$OMP end do - end if ! test on var_diags_opt !$OMP end PARALLEL end if ! l_subgrid_qcl_mp .or. l_wvar_for_conv - if (var_diags_opt == original_vars) then - ! at this point, BL_diag%tke really contains 1.5*sigma_w^2. To make it look - ! a bit more like TKE near the surface, we will keep it constant below - ! the max of rhokm_surf -!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP SHARED( pdims, rhokmz, BL_diag ) private(i, j, k, kmax ) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - kmax = maxloc(rhokmz(i,j,:), dim=1) - do k = 2, kmax - BL_diag%tke(i,j,k) = BL_diag%tke(i,j,kmax) - end do - end do - end do -!$OMP end PARALLEL do - else if (var_diags_opt == split_tke_and_inv) then - ! at this point, tke_nl really contained 1.5*sigma_w^2. To make it look - ! a bit more like TKE near the surface, we will keep it constant below - ! the max of rhokm_surf (here we find the first local maximum in case there - ! is a larger rhokmz in a resolved inversion + ! At this point, tke_nl really contained 1.5*sigma_w^2. To make it look + ! a bit more like TKE near the surface, we will keep it constant below + ! the max of rhokm_surf (here we find the first local maximum in case there + ! is a larger rhokmz in a resolved inversion !$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & !$OMP SHARED( pdims, bl_levels, rhokmz, BL_diag, tke_loc, tke_nl ) & !$OMP private(i, j, k, kp) - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - kp=2 - do while ( rhokmz(i,j,kp+1) > rhokmz(i,j,kp) .and. kp < bl_levels ) - kp = kp +1 - end do - do k = 2, kp-1 - BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,kp) ) - end do + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + kp=2 + do while ( rhokmz(i,j,kp+1) > rhokmz(i,j,kp) .and. kp < bl_levels ) + kp = kp +1 + end do + do k = 2, kp-1 + BL_diag%tke(i,j,k) = max( tke_loc(i,j,k), tke_nl(i,j,kp) ) end do end do + end do !$OMP end PARALLEL do - end if ! test on var_diags_opt end if ! BL_diag%L_tke @@ -2915,14 +2826,6 @@ subroutine bdy_expl2 ( & ! slvar,qwvar,slqw are on theta-level K ! tke is on theta-level K-1 ! exner is on theta-level K-1 - - ! For var_diags_opt = original_vars we use the following variables: - ! dsldzm is on theta-level K-1 - ! dqwdzm is on theta-level K-1 - ! elm is on theta-level K-1 - ! rhokh is on rho-level K - - ! For var_diags_opt = split_tke_and_inv ! dsldz is on rho-level K ! dqwdz is on rho-level K ! ftl is on rho-level k @@ -2943,85 +2846,47 @@ subroutine bdy_expl2 ( & call qsat_wat(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1),tdims%i_len) end if - if (var_diags_opt == original_vars) then - - do i = tdims%i_start, tdims%i_end - ! calculate sh, don't interpolate with the surface flux or divide by 0 - if ( BL_diag%tke(i,j,k) > 1.0e-10_r_bl ) then - sh = rhokh(i,j,k) & - / ( rho_wet_tq(i,j,k-1) * elm(i,j,k) & - * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) - else - sh = small_sh - end if - if ( sh < small_sh ) sh = small_sh - ! calculate exner - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - ! calculate the variance, use gradient interpolated between levs 1 and 2 - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & - * elm(i,j,k)**2 * dsldz(i,j,k)**2 & - + a_qs(i,j,k-1)**2 * b2 * sh & - * elm(i,j,k)**2 * dqwdz(i,j,k)**2 & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & - * sh * elm(i,j,k)**2 * dsldz(i,j,k) * dqwdz(i,j,k) + do i = tdims%i_start, tdims%i_end + ! calculate the variance + sl_var = zero + qw_var = zero + sl_qw = zero + if ( BL_diag%tke(i,j,k) > small_tke ) then + ! vertical interpolation weights + weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) + weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) + ! var_fac=b2*L/sqrt(TKE) + var_fac = b2 * rhokm(i,j,k) / ( weight1 * BL_diag%tke(i,j,k) * & + rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1) ) + ! Note that flux*gradient can be negative so the absolute values + ! are used + qw_var= abs( -var_fac*( weight2 * fqw(i,j,k) * dqwdz(i,j,k) + & + weight3 * fqw(i,j,k-1)* dqwdz(i,j,k-1) ) ) + sl_var= abs( -var_fac*( weight2 * ftl(i,j,k) * dsldz(i,j,k) + & + weight3 * ftl(i,j,k-1)* dsldz(i,j,k-1) ) ) + sl_qw = - one_half * var_fac*( & + weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & + weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& + ) if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k)**2 + BL_diag%slvar(i,j,k-1) = sl_var end if if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dqwdz(i,j,k)**2 + BL_diag%qwvar(i,j,k-1) = qw_var end if if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k) & - * dqwdz(i,j,k) + BL_diag%slqw(i,j,k-1) = sl_qw end if - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - end do - - else ! var_diags_opt - - do i = tdims%i_start, tdims%i_end - ! calculate the variance - sl_var = zero - qw_var = zero - sl_qw = zero - if ( BL_diag%tke(i,j,k) > small_tke ) then - ! vertical interpolation weights - weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - ! var_fac=b2*L/sqrt(TKE) - var_fac = b2 * rhokm(i,j,k) / ( weight1 * BL_diag%tke(i,j,k) * & - rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1) ) - ! Note that flux*gradient can be negative so the absolute values - ! are used - qw_var= abs( -var_fac*( weight2 * fqw(i,j,k) * dqwdz(i,j,k) + & - weight3 * fqw(i,j,k-1)* dqwdz(i,j,k-1) ) ) - sl_var= abs( -var_fac*( weight2 * ftl(i,j,k) * dsldz(i,j,k) + & - weight3 * ftl(i,j,k-1)* dsldz(i,j,k-1) ) ) - sl_qw = - one_half * var_fac*( & - weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & - weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& - ) - if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = sl_var - end if - if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = qw_var - end if - if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = sl_qw - end if - end if - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & - + a_qs(i,j,k-1)**2 * qw_var & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) + end if + exner = ( p_theta_levels(i,j,k-1) / pref )**kappa + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & + + a_qs(i,j,k-1)**2 * qw_var & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) - end do !i - end if ! var_diags_opt + end do !i if (i_rhcpt == rhcpt_tke_based) then do i = tdims%i_start, tdims%i_end @@ -3057,110 +2922,61 @@ subroutine bdy_expl2 ( & call qsat_wat(qsw_arr,tl(:,j,k-1),p_theta_levels(:,j,k-1),tdims%i_len) end if - if (var_diags_opt == original_vars) then - - do i = tdims%i_start, tdims%i_end - ! calculate sh, don't divide by 0 - if ( BL_diag%tke(i,j,k) > 1.0e-10_r_bl ) then - weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1)-r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - sh = ( weight2 * rhokh(i,j,k) + weight3 * rhokh(i,j,k-1) ) & - / ( weight1 * rho_wet_tq(i,j,k-1) * elm(i,j,k) & - * sqrt( 2.0_r_bl * BL_diag%tke(i,j,k) ) ) - else - sh = small_sh - end if - if ( sh < small_sh ) sh = small_sh - ! calculate exner - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - ! calculate the variance - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * b2 * sh & - * elm(i,j,k)**2 * dsldzm(i,j,k)**2 & - + a_qs(i,j,k-1)**2 * b2 * sh & - * elm(i,j,k)**2 * dqwdzm(i,j,k)**2 & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * b2 & - * sh * elm(i,j,k)**2 * dsldzm(i,j,k) * dqwdzm(i,j,k) + do i = tdims%i_start, tdims%i_end + ! calculate the variance + sl_var = zero + qw_var = zero + sl_qw = zero + if ( BL_diag%tke(i,j,k) > small_tke ) then + ! vertical interpolation weights + weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) + weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) + weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) + var_fac = b2 * rhokm(i,j,k) / ( weight1*BL_diag%tke(i,j,k) * & + rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1)) + kp=k + km=k-1 + ! Don't use level ntml (or ntdsc) if disc_inv=2 as this indicates + ! the inversion has just risen and the gradients between ntml and + ! ntml-1 are likely to give excessive variances + if ( (kp == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & + (kp == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) kp = km + if ( (km == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & + (km == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) km = kp + ! Note that flux*gradient can be negative so the absolute values + ! are used + qw_var= abs( -var_fac*( weight2 * fqw(i,j,kp) * dqwdz(i,j,kp) + & + weight3 * fqw(i,j,km) * dqwdz(i,j,km) ) ) + sl_var= abs( -var_fac*( weight2 * ftl(i,j,kp) * dsldz(i,j,kp) + & + weight3 * ftl(i,j,km) * dsldz(i,j,km) ) ) + sl_qw = - one_half * var_fac*( & + weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & + weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& + ) if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldz(i,j,k)**2 + BL_diag%slvar(i,j,k-1) = sl_var end if if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dqwdz(i,j,k)**2 + BL_diag%qwvar(i,j,k-1) = qw_var end if if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = b2 * sh * elm(i,j,k)**2 * dsldzm(i,j,k) & - * dqwdz(i,j,k) + BL_diag%slqw(i,j,k-1) = sl_qw end if - end do !i - if (i_rhcpt == rhcpt_tke_based) then - do i = tdims%i_start, tdims%i_end - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - ! calculate rhcrit, with appropriate limits - rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) - end do !i end if - - else ! var_diags_opt - + exner = ( p_theta_levels(i,j,k-1) / pref )**kappa + sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & + + a_qs(i,j,k-1)**2 * qw_var & + - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw + end do !i + if (i_rhcpt == rhcpt_tke_based) then do i = tdims%i_start, tdims%i_end - ! calculate the variance - sl_var = zero - qw_var = zero - sl_qw = zero - if ( BL_diag%tke(i,j,k) > small_tke ) then - ! vertical interpolation weights - weight1 = r_rho_levels(i,j,k) - r_rho_levels(i,j,k-1) - weight2 = r_theta_levels(i,j,k-1) - r_rho_levels(i,j,k-1) - weight3 = r_rho_levels(i,j,k) - r_theta_levels(i,j,k-1) - var_fac = b2 * rhokm(i,j,k) / ( weight1*BL_diag%tke(i,j,k) * & - rho_wet_tq(i,j,k-1) * rho_mix_tq(i,j,k-1)) - kp=k - km=k-1 - ! Don't use level ntml (or ntdsc) if disc_inv=2 as this indicates - ! the inversion has just risen and the gradients between ntml and - ! ntml-1 are likely to give excessive variances - if ( (kp == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & - (kp == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) kp = km - if ( (km == ntml(i,j) .and. sml_disc_inv(i,j) == 2) .or. & - (km == ntdsc(i,j) .and. dsc_disc_inv(i,j) == 2) ) km = kp - ! Note that flux*gradient can be negative so the absolute values - ! are used - qw_var= abs( -var_fac*( weight2 * fqw(i,j,kp) * dqwdz(i,j,kp) + & - weight3 * fqw(i,j,km) * dqwdz(i,j,km) ) ) - sl_var= abs( -var_fac*( weight2 * ftl(i,j,kp) * dsldz(i,j,kp) + & - weight3 * ftl(i,j,km) * dsldz(i,j,km) ) ) - sl_qw = - one_half * var_fac*( & - weight2*( ftl(i,j,k)*dqwdz(i,j,k) + fqw(i,j,k)*dsldz(i,j,k) ) + & - weight3*( ftl(i,j,k-1)*dqwdz(i,j,k-1) + fqw(i,j,k-1)*dsldz(i,j,k-1))& - ) - if (BL_diag%l_slvar) then - BL_diag%slvar(i,j,k-1) = sl_var - end if - if (BL_diag%l_qwvar) then - BL_diag%qwvar(i,j,k-1) = qw_var - end if - if (BL_diag%l_slqw) then - BL_diag%slqw(i,j,k-1) = sl_qw - end if - end if - exner = ( p_theta_levels(i,j,k-1) / pref )**kappa - sgm(i) = a_dqsdt(i,j,k-1)**2 * exner**2 * sl_var & - + a_qs(i,j,k-1)**2 * qw_var & - - 2.0_r_bl * a_qs(i,j,k-1) * a_dqsdt(i,j,k-1) * exner * sl_qw + ! do this for safety, not sure if it's really needed + sgm(i) = sqrt ( max( sgm(i), zero ) ) + ! calculate rhcrit, with appropriate limits + rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & + one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) end do !i - if (i_rhcpt == rhcpt_tke_based) then - do i = tdims%i_start, tdims%i_end - ! do this for safety, not sure if it's really needed - sgm(i) = sqrt ( max( sgm(i), zero ) ) - ! calculate rhcrit, with appropriate limits - rhcpt(i,j,k-1) = min( max_rhcpt(i,j), max( min_rhcpt(i,j), & - one - root6 * sgm(i) / (a_qs(i,j,k-1) * qsw_arr(i)))) - end do !i - end if - - end if ! var_diags_opt + end if end do !j !$OMP end do NOWAIT diff --git a/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 b/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 index d6fccca90..faaa5a96a 100644 --- a/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 +++ b/science/physics_schemes/source/boundary_layer/bl_option_mod.F90 @@ -64,7 +64,7 @@ module bl_option_mod integer, parameter :: sharp_sea_long_land = 2 integer, parameter :: mes_tails = 3 integer, parameter :: louis_tails = 4 -integer, parameter :: depth_based = 5 +integer, parameter :: depth_based = 5 ! not available in LFRic integer, parameter :: sharp_sea_mes_land = 6 integer, parameter :: lem_stability = 7 integer, parameter :: sharp_sea_louis_land = 8 @@ -291,11 +291,6 @@ module bl_option_mod ! 13 Improved method to calculate cloud-top radiative flux jump logical :: l_new_kcloudtop = .false. -! 14 Switch for alternative TKE and variance diagnostics -integer :: var_diags_opt = imdi -integer, parameter :: original_vars = 0 -integer, parameter :: split_tke_and_inv = 1 - ! 15 Multiplicative tuning factor in TKE diagnosis, usually 1.0 real(kind=r_bl) :: tke_diag_fac = rmdi @@ -503,12 +498,6 @@ module bl_option_mod ! Switch for coupled gradient method in Equilibrium SBL model logical, parameter :: L_SBLco = .true. -! LambdaM=2*LambdaH (operational setting) -logical, parameter :: l_lambdam2 = .false. - -! Lambdas not reduced above NTML_LOCAL+1 -logical, parameter :: l_full_lambdas = .false. - ! logical for whether to skip calculations based on start of timestep ! quantities when using semi-lagrangian cycling with Endgame ! This is set to true in dynamics_input_mod as Endgame always uses it, @@ -531,7 +520,7 @@ module bl_option_mod l_use_surf_in_ri, lambda_min_nml, ritrans, c_gust, & dzrad_disc_opt, num_sweeps_bflux, l_converge_ga, & local_fa, Keep_Ri_FA, l_bl_mix_qcf, l_conv_tke, l_use_var_fixes, & - l_reset_neg_q, var_diags_opt, tke_diag_fac, i_interp_local, & + l_reset_neg_q, tke_diag_fac, i_interp_local, & sg_orog_mixing, fric_heating, calc_prob_of_vis, z_nl_bl_levels, & idyndiag, zhloc_depth_fac, flux_grad, entr_smooth_dec, & sc_diag_opt, kprof_cu, bl_res_inv, blending_option, & @@ -632,8 +621,6 @@ subroutine print_nlist_run_bl() call umprint(linebuffer,src='bl_option_mod') write(linebuffer,'(A,L1)') 'l_converge_ga = ',l_converge_ga call umprint(linebuffer,src='bl_option_mod') -write(linebuffer,'(A,I4)') 'var_diags_opt = ',var_diags_opt -call umprint(linebuffer,src='bl_option_mod') write(linebuffer,'(A,L1)') 'l_use_var_fixes = ',l_use_var_fixes call umprint(linebuffer,src='bl_option_mod') write(linebuffer,'(A,ES12.4)') 'tke_diag_fac = ',tke_diag_fac diff --git a/science/physics_schemes/source/boundary_layer/ex_coef.F90 b/science/physics_schemes/source/boundary_layer/ex_coef.F90 index fd9ffe3e1..6127e8dfa 100644 --- a/science/physics_schemes/source/boundary_layer/ex_coef.F90 +++ b/science/physics_schemes/source/boundary_layer/ex_coef.F90 @@ -24,11 +24,10 @@ module ex_coef_mod subroutine ex_coef ( & ! in levels/logicals - bl_levels,k_log_layr,nSCMDpkgs,L_SCMDiags, BL_diag, & + bl_levels, k_log_layr, BL_diag, & ! in fields - sigma_h,flandg,dbdz,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m,h_blend,zhpar, & - ntpar,ntml_nl,ntdsc,nbdsc,u_p,v_p,v_s,fb_surf,qw,tl,l_shallow_cth,rmlmax2, & - rneutml_sq, delta_smag, & + sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m,zhpar,ntpar, & + ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq, delta_smag, & ! in/out fields cumulus,weight_1dbl, & ! out fields @@ -39,9 +38,9 @@ subroutine ex_coef ( & use atm_fields_bounds_mod, only: pdims, tdims, pdims_s use bl_diags_mod, only: strnewbldiag use bl_option_mod, only: WeightLouisToLong, Variable_RiC, cbl_op, & - sg_orog_mixing, ricrit_sharp, pr_max, l_lambdam2, l_full_lambdas, & + sg_orog_mixing, ricrit_sharp, pr_max, & local_fa,Prandtl,ishear_bl,L_SBLco,Muw_SBL,Mwt_SBL,sbl_op, & - LockMailhot2004, depth_based, lem_stability, lem_std, lem_conven, & + LockMailhot2004, lem_stability, lem_std, lem_conven, & lem_adjust, cbl_mix_fac_nml, & off, on, sharpest, sharp_sea_long_land, sharp_sea_mes_land, & louis_tails, sharp_sea_louis_land, long_tails, mes_tails, ritrans, & @@ -49,8 +48,7 @@ subroutine ex_coef ( & lambda_fac, beta_bl, beta_fa, rlinfac, linear0, & to_sharp_across_1km, ntml_level_corrn, free_trop_layers, two_thirds, & blending_option, blend_except_cu, blend_gridindep_fa, blend_cth_shcu_only, & - extended_tail, sg_shear_enh_lambda, l_use_var_fixes, var_diags_opt, & - split_tke_and_inv, zero, one, one_half + extended_tail, zero, one, one_half use conversions_mod, only: pi => pi_bl use gen_phys_inputs_mod, only: l_mr_physics @@ -63,8 +61,6 @@ subroutine ex_coef ( & use turb_diff_mod, only: l_subfilter_vert, l_subfilter_horiz use model_domain_mod, only: model_type, mt_single_column -use s_scmop_mod, only: default_streams, & - t_avg, d_sl, scmdiag_bl use yomhook, only: lhook, dr_hook use parkind1, only: jprb, jpim @@ -98,20 +94,6 @@ subroutine ex_coef ( & bl_levels), & ! in density on theta levels; ! used in RHOKM so wet density - dbdz(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & - 2:bl_levels), & - ! in Buoyancy gradient across lower - ! interface of layer. - u_p(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels), & - ! in Westerly wind component horizontally - ! interpolated to P-grid. (m/s) - v_p(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels), & - ! in Southerly wind component horizontally - ! interpolated to P-grid. (m/s) - qw(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Total water content (kg per kg air). - tl(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & - ! in Liquid/frozen water temperature (K). rmlmax2(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels), & ! in Square of asymptotic mixing length for Smagorinsky scheme z_uv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels+1), & @@ -123,18 +105,11 @@ subroutine ex_coef ( & ! in Height of top of initial parcel ascent z0m(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! in Roughness length for momentum (m). - h_blend(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Blending height for effective - ! roughness length scheme dvdzm(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & 2:bl_levels), & ! in Modulus of wind shear. ri(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,2:bl_levels), & ! in Local Richardson number. - v_s(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Surface friction velocity (m/s) - fb_surf(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! in Surface buoyancy flux over density (m^2/s^3). flandg(pdims_s%i_start:pdims_s%i_end,pdims_s%j_start:pdims_s%j_end), & ! in Land fraction on all tiles. rneutml_sq(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels), & @@ -146,13 +121,6 @@ subroutine ex_coef ( & l_shallow_cth(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) ! in Flag to indicate shallow convection based on cl-top -! Additional variables for SCM diagnostics which are dummy in full UM -integer, intent(in) :: & - nSCMDpkgs ! in No of SCM diagnostics packages - -logical, intent(in) :: & - L_SCMDiags(nSCMDpkgs) ! Logicals for SCM diagnostics packages - ! Declaration of new BL diagnostics. type (strnewbldiag), intent(in out) :: BL_diag @@ -216,7 +184,7 @@ subroutine ex_coef ( & character(len=*), parameter :: RoutineName = 'EX_COEF' -real(kind=r_bl) :: eh,em,g0,dh,dm,a_lambda,r_c_tke +real(kind=r_bl) :: eh,em,g0,dh,dm,r_c_tke real(kind=r_bl) :: subbmin,subbmax,subcmin,subcmax real(kind=r_bl) :: a_ri,b_ri @@ -225,8 +193,6 @@ subroutine ex_coef ( & ! Used in calc of stability function FH. em=4.0_r_bl, & ! Used in calc of stability function FM. - a_lambda=2.0_r_bl, & - ! used in calc of LAMBDA_EFF r_c_tke=one/0.41_r_bl & ! used in calc of TKE (1/stress-energy ratio, see UMDP25) ) @@ -246,20 +212,6 @@ subroutine ex_coef ( & ! Used in calc of unstability function FM. ) -! Equilibrium SBL model constants -real(kind=r_bl) :: RtestMin -integer :: gn,NGz,kMINh -parameter ( & - RtestMin=zero, & - ! Threshold for Rtest - gn=19, & - ! Size of "G"-tables (No. HonL values) - NGz=90, & - ! No. z/h steps in each "G" integration - kMINh=2 & - ! Level of minimum SBL height (>=2) -) - ! Define local storage. real(kind=r_bl) :: & ricrit(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & @@ -268,21 +220,6 @@ subroutine ex_coef ( & ! 2D variable for SBL stabiliy function options sharp(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! 2D variable for SHARP stabiliy function - invMOsurf(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! Inverse of sfce M-O length - ! Note: Inverse is used so that neutral conditions - ! can be handled (M-O length --> infinity) - zh_esbl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! Ht of equilib SBL (sub-grid - HLtab(gn), & - ! Lookup tables (Gx calcs - GHsav(gn,NGz),gmsav(gn,NGz), & - ! in equilib SBL scheme) - THv_TQ(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end, & - bl_levels), & - ! Virtual potential temperature on theta levels - THv(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end,bl_levels), & - ! THv_TQ interpolated to U,V levels prandtl_number(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & ! = KM/KH (currently only calculated for stable) BL_weight(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end, & @@ -296,10 +233,6 @@ subroutine ex_coef ( & weight_bltop(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) ! weight_1dbl at the top of the PBL -integer :: & - ntml_esbl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) - ! No. UV-levels inside equilibrium SBL - logical :: & topbl(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) ! Flag for having reached the @@ -319,25 +252,12 @@ subroutine ex_coef ( & ! z/sigma_h rpr ! !/Pr -! Variables for boundary layer depth based formulation -real(kind=r_bl) :: & - h_tkeb(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! TKE budget based BL depth - MOsurf(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end), & - ! surface Obukhov length - diff_min(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end) real(kind=r_bl) :: & - h_est, & - rifb, & - ! Bulk flux Richardson number pr_n, & ! neutral Prandtl number - r_pr_n, & + r_pr_n ! 1 / neutral Prandtl number - m_tau,m_buoy, & - ! Indices for implied stress and buoyancy flux profs - ind,diff real(kind=r_bl) :: & subb, subc, subg, ric, ricinv, rifac @@ -373,12 +293,7 @@ subroutine ex_coef ( & ! top of boundary layer mixing zfa, & ! height to use beta_fa in blendin - lambda_eff ! Effective mixing length used with effective - ! roughness length scheme. - -!Equilibrium SBL model temporary real scalar variables -real(kind=r_bl) :: & - km, u1, zz + zz integer :: & i,j, & @@ -388,38 +303,15 @@ subroutine ex_coef ( & kb, kt ! Base and top level of unstable Ri layers -!Equilibrium SBL model temporary integer scalar variables -integer :: & - kZtop,kZbot,gk,kG0, & - ! Temporary loop counters - iERRSBL ! SBLequil error status - -!Equilibrium SBL model logical variables logical :: & - GcalcDone, & - ! Calculation of Gx values has been performed - subcrit, & - ! flag for being in a subcritical ri layer - subgrid ! Will perform subgrid SBL depth calculation - -!Switch to enable subgrid SBL depth diagnosis -logical :: sg_enabled -parameter (sg_enabled=.true.) - -!Equilibrium SBL model SAVED variables -save HLtab,GHsav,gmsav,GcalcDone - -!Equilibrium SBL model data statements -data HLtab /0.0001_r_bl,0.001_r_bl,0.002_r_bl,0.005_r_bl,0.01_r_bl,0.02_r_bl, & - 0.05_r_bl,0.1_r_bl,0.2_r_bl,one_half,one,2.0_r_bl,5.0_r_bl, & - 10.0_r_bl,20.0_r_bl,50.0_r_bl,100.0_r_bl,200.0_r_bl,500.0_r_bl/ -data GcalcDone /.false./ + subcrit ! flag for being in a subcritical ri layer integer(kind=jpim), parameter :: zhook_in = 0 integer(kind=jpim), parameter :: zhook_out = 1 real(kind=jprb) :: zhook_handle if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + !----------------------------------------------------------------------- ! if stochastic physics random parameters is used set the parameter ! used to vary the stability function to a perturbed value, if not @@ -443,7 +335,6 @@ subroutine ex_coef ( & !--------------------------------------------------------------- pr_n = one if (Prandtl == LockMailhot2004) pr_n = 0.7_r_bl -if (sbl_op == depth_based) pr_n = 0.7_r_bl ! Use pr_n=0.7 if any LEM stability function selected if (sbl_op == lem_stability .or. cbl_op == lem_std & .or. cbl_op == lem_conven & @@ -504,7 +395,6 @@ subroutine ex_coef ( & ! 0. Initialise flag for having reached top of turbulently mixed layer !----------------------------------------------------------------------- topbl(i,j) = .false. - prandtl_number(i,j) = pr_n ! initialise blending weight at top of BL to one weight_bltop(i,j) = one @@ -512,97 +402,6 @@ subroutine ex_coef ( & end do !$OMP end do NOWAIT !----------------------------------------------------------------------- -! Set-up a BL weighting function, =1 near the ground (ie in the BL) -! =0 in the free troposphere -! Rate and height at which transition occurs varys depending on choices: -!----------------------------------------------------------------------- -!$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 - BL_weight(i,j,k) = one - end do - end do -end do -!$OMP end do - -if (local_fa == to_sharp_across_1km) then - !--------------------------------------------------------- - ! Additional code to allow the local Ri scheme to use - ! SHARPEST in the free atmosphere, ie above the BL top, - ! regardless of the tail option selected above. - ! Set Z_SCALE to 1km to mimic old value of BL_LEVELS, - ! gives BL_weight~0 by 2km, ~0.95 at 500m - !--------------------------------------------------------- - z_scale = 1000.0_r_bl -!$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 - zpr = z_tq(i,j,k-1)/z_scale - BL_weight(i,j,k) = one_half*(one - tanh(3.0_r_bl*(zpr-one) ) ) - end do - end do - end do -!$OMP end do NOWAIT -end if - -if (sg_orog_mixing /= off) then - !----------------------------------------------------------------- - ! Subgrid orographic height dependence for SBL tail (option 1) - ! or orographic dependence of mixing lengths, lambdam,h (opt 2) - ! Gives BL_weight~[1,0.95,0.5,0] at ZPR=[0,0.6,1,1.7] - !---------------------------------------------------------------- -!$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 (sigma_h(i,j) > one ) then - zpr = z_tq(i,j,k-1)/sigma_h(i,j) - BL_weight(i,j,k) = one_half*( one - & - tanh(4.0_r_bl*(zpr-one) ) ) - end if - end do - end do - end do -!$OMP end do NOWAIT -end if - -if (l_use_var_fixes) then -!$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 - turb_length(i,j,k) = lambda_min*rlambda_fac - end do - end do - end do -!$OMP end do NOWAIT - if (blending_option == blend_cth_shcu_only) then - ! use Smag mixing length as background length scale if smaller - ! than lambda_min (ie ignore lambda_min for high res simulations) -!$OMP do SCHEDULE(STATIC) - do k = 2, bl_levels - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - turb_length(i,j,k) = min( turb_length(i,j,k), sqrt(rmlmax2(i,j,k)) ) - end do - end do - end do -!$OMP end do NOWAIT - end if -else -!$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 - turb_length(i,j,k)=lambda_min - end do - end do - end do -!$OMP end do NOWAIT -end if -!----------------------------------------------------------------------- ! Set critical Richardson number !----------------------------------------------------------------------- if (l_rp2) then @@ -699,20 +498,18 @@ subroutine ex_coef ( & end if !----------------------------------------------------------------------- -! 1.1 Loop over levels calculating Richardson numbers. +! 1.1 Use Richardson number profile to calculate BL depth, zh !----------------------------------------------------------------------- - !$OMP PARALLEL DEFAULT(none) private(k,j,i) & !$OMP SHARED(bl_levels,pdims,topbl,ri,ricrit,local_fa,ntml_local,zh_local,z_uv) do k = 2, bl_levels !$OMP do SCHEDULE(STATIC) do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end - !------------------------------------------------------------------ - ! 1.2 If either a stable layer (Ri>RiCrit) or the maximum BL - ! height has been reached, set boundary layer height (ZH_LOCAL) to - ! the height of the lower boundary of the current layer + ! If either a stable layer (Ri>RiCrit) or the maximum BL + ! height has been reached, set boundary layer height (ZH_LOCAL) to + ! the height of the lower boundary of the current layer !------------------------------------------------------------------ if ( .not. topbl(i,j) .and. & (ri(i,j,k) > ricrit(i,j) .or. k == bl_levels) ) then @@ -741,18 +538,17 @@ subroutine ex_coef ( & end if !----------------------------------------------------------------------- -! In CUMULUS layers the local scheme is capped at the LCL (given in -! this case by NTML_NL). Save local BL depth as SCM diagnostic. -!----------------------------------------------------------------------- -!----------------------------------------------------------------------- -! If NTML_LOCAL is greater than the top of the parcel ascent (NTPAR) -! for a cumulus-capped layer, shear driven mixing is allowed to -! dominate (if ISHEAR_BL=1 selected) +! 1.2 In CUMULUS layers the local scheme is capped at the LCL (given in +! this case by NTML_NL). +! If NTML_LOCAL is greater than the top of the parcel ascent (NTPAR) +! for a cumulus-capped layer, shear driven mixing is allowed to +! dominate (if ISHEAR_BL=1 selected) !----------------------------------------------------------------------- !$OMP PARALLEL DEFAULT(none) & -!$OMP SHARED( pdims, ishear_bl, l_use_var_fixes, ntml_local, ntpar, cumulus, & -!$OMP ntml_nl, zh_local, z_uv ) & -!$OMP private( i, j ) +!$OMP SHARED( pdims, ishear_bl, ntml_local, ntpar, cumulus, & +!$OMP ntml_nl, zh_local, z_uv, bl_levels, lambda_min, rlambda_fac, & +!$OMP turb_length, blending_option, rmlmax2) & +!$OMP private( i, j, k ) !$OMP do SCHEDULE(STATIC) do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end @@ -762,30 +558,37 @@ subroutine ex_coef ( & end do end do !$OMP end do - -if ( .not. l_use_var_fixes ) then - ! Use sub-cloud layer as local PBL depth (as for non-local). - ! Found to give large TKE and variances with new diagnostics - ! and isn't particularly justifiable +!----------------------------------------------------------------------- +! 1.3 Search for sub-critical layers above the PBL and set the +! mixing length to scale with these layer depths +!----------------------------------------------------------------------- !$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 ( cumulus(i,j) ) then - ntml_local(i,j) = ntml_nl(i,j) - zh_local(i,j) = z_uv(i,j,ntml_local(i,j)+1) - end if + turb_length(i,j,k) = lambda_min*rlambda_fac end do end do -!$OMP end do +end do +!$OMP end do NOWAIT +if (blending_option == blend_cth_shcu_only) then + ! use Smag mixing length as background length scale if smaller + ! than lambda_min (ie ignore lambda_min for high res simulations) +!$OMP do SCHEDULE(STATIC) + do k = 2, bl_levels + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + turb_length(i,j,k) = min( turb_length(i,j,k), sqrt(rmlmax2(i,j,k)) ) + end do + end do + end do +!$OMP end do NOWAIT end if !$OMP end PARALLEL -!----------------------------------------------------------------------- -! 1.3 Search for sub-critical layers above the PBL and set the -! mixing length to scale with these layer depths -!----------------------------------------------------------------------- + if (local_fa == free_trop_layers) then !$OMP PARALLEL do DEFAULT(none) SCHEDULE(STATIC) & -!$OMP SHARED( pdims, bl_levels, ntml_local, ri, ricrit, z_uv, l_use_var_fixes, & +!$OMP SHARED( pdims, bl_levels, ntml_local, ri, ricrit, z_uv, & !$OMP turb_length, rlambda_fac, lambda_min ) & !$OMP private( i, j, k, subcrit, kb, kt, kl, turb_length_layer ) do j = pdims%j_start, pdims%j_end @@ -806,17 +609,10 @@ subroutine ex_coef ( & ! turb_length(k) is held, with Ri(k), on th-level(k-1) !--------------------------------------------------------- turb_length_layer = z_uv(i,j,kt) - z_uv(i,j,kb-1) - if (l_use_var_fixes) then - do kl = kb, kt - turb_length(i,j,kl) = max( turb_length(i,j,kl), & - min(turb_length_layer,lambda_max_nml*rlambda_fac) ) - end do - else - do kl = kb, kt - turb_length(i,j,kl) = max( lambda_min*rlambda_fac, & - min(turb_length_layer,lambda_max_nml*rlambda_fac) ) - end do - end if + do kl = kb, kt + turb_length(i,j,kl) = max( turb_length(i,j,kl), & + min(turb_length_layer,lambda_max_nml*rlambda_fac) ) + end do end if end do @@ -851,94 +647,256 @@ subroutine ex_coef ( & !$OMP end PARALLEL do end if !----------------------------------------------------------------------- -! 2. Richardson Number based local mixing scheme +! 2.0 Loop over levels; calculate the mixing lengths !----------------------------------------------------------------------- -! 2.0 Loop round "boundary" levels; calculate the stability- -! dependent turbulent mixing coefficients. -!----------------------------------------------------------------------- -! TKE budget based depth diagnosis +do k = 2, bl_levels +!$OMP PARALLEL DEFAULT(none) & +!$OMP PRIVATE(z_scale,j,i,lambdam,lambdah, & +!$OMP lambdah_rho,vkz,f_log,zz,zht,zfa,beta) & +!$OMP SHARED(k,pdims,ri,ricrit,flandg,ntml_local,ntml_nl,subb,dh,z_tq, & +!$OMP l_rp2,lambda_min,par_mezcla_rp,zh_local,turb_length,k_log_layr, & +!$OMP z_uv,z0m,elm,elh,elh_rho,blending_option,cumulus,l_shallow_cth,zhpar, & +!$OMP ntdsc,weight_1dbl,weight_bltop,delta_smag,rneutml_sq,BL_diag, & +!$OMP cbl_op,rho_wet_tq,l_subfilter_vert,l_subfilter_horiz,local_fa) + !----------------------------------------------------------------- + ! 2.1 Calculate asymptotic mixing lengths LAMBDAM and LAMBDAH + !----------------------------------------------------------------- +!$OMP do SCHEDULE(STATIC) + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end + if (l_rp2) then + lambdam = max ( lambda_min , par_mezcla_rp(rp_idx)*zh_local(i,j) ) + else + lambdam = max ( lambda_min , lambda_fac*zh_local(i,j) ) + end if + !----------------------------------------------------------------- + ! Reduce mixing lengths above BL + !----------------------------------------------------------------- + if (k >= ntml_local(i,j)+2) then + lambdam = lambda_min + end if -! Starting with the definition of the flux Richardson -! number, assuming similarity profiles for -! stress and buoyancy flux, and vertically integrating -! gives an expression for the stable boundary layer -! depth which is based just on surface fluxes and -! the wind speed change across the boundary layer. -! ---------------------------------------------------- + lambdah = lambdam + lambdah_rho = lambdah -if (sbl_op == depth_based) then + if ( local_fa == free_trop_layers ) then + lambdam = max( lambdam, lambda_fac*turb_length(i,j,k) ) + lambdah = max( lambdah, lambda_fac*turb_length(i,j,k) ) + ! lambdah_rho does not need to be recalculated under + ! local_fa option "free_trop_layers" as the full KH profile + ! will be interpolated in bdy_expl2 + end if + !----------------------------------------------------------------------- + ! 2.2 Calculate mixing lengths ELH, ELM coincident with RI(K) and so + ! at Z_TQ(K-1) + !----------------------------------------------------------------------- + ! Incorporate log profile corrections to the vertical finite + ! differences into the definitions of ELM and ELH. + ! Note that ELH_RHO is calculated (on rho levels) for direct inclusion + ! in RHOKH and also (as elh) on theta levels for the unstable + ! stability functions and inclusion in RHOKH before interpolation + ! (under local_fa option "free_trop_layers"). + ! To save computing logarithms for all K, the values of ELM and ELH + ! are unchanged for K > K_LOG_LAYR. - ! Index for assumed buoyancy profile - m_buoy=one + if (k <= k_log_layr) then + vkz = vkman * ( z_uv(i,j,k) - z_uv(i,j,k-1) ) + f_log = log( ( z_uv(i,j,k) + z0m(i,j) ) / & + ( z_uv(i,j,k-1) + z0m(i,j) ) ) + elm(i,j,k) = vkz / ( f_log + vkz/lambdam ) + elh(i,j,k) = vkz / ( f_log + vkz/lambdah ) + vkz = vkman * ( z_tq(i,j,k) - z_tq(i,j,k-1) ) + f_log = log( ( z_tq(i,j,k) + z0m(i,j) ) / & + ( z_tq(i,j,k-1) + z0m(i,j) ) ) + elh_rho(i,j,k) = vkz / ( f_log + vkz/lambdah_rho ) + else + vkz = vkman * ( z_tq(i,j,k-1) + z0m(i,j) ) + elm(i,j,k) = vkz / (one + vkz/lambdam ) + elh(i,j,k) = vkz / (one + vkz/lambdah ) + vkz = vkman * ( z_uv(i,j,k) + z0m(i,j) ) + elh_rho(i,j,k) = vkz / (one + vkz/lambdah_rho ) + end if + end do + end do +!$OMP end do +!---------------------------------------------------------------- +! 2.3 Blend mixing lengths between 1D and 3D Smagorinsky +!---------------------------------------------------------------- + if (blending_option /= off) then +!$OMP do SCHEDULE(STATIC) + do j = pdims%j_start, pdims%j_end + do i = pdims%i_start, pdims%i_end - ! Index for assumed stress profile - m_tau=one + zz = z_tq(i,j,k-1) ! height of rhokm(k) + ! turb_length is the greater of the local and non-local + ! BL depths up to that bl top + z_scale = max( zz, turb_length(i,j,k) ) + ! zht = interface between BL and FA + zht = max( z_uv(i,j,ntml_nl(i,j)+1) , zh_local(i,j) ) + ! Relevant scale in cumulus layers can be cloud top height, zhpar + if ( cumulus(i,j) .and. ( blending_option /= blend_cth_shcu_only .or. & + l_shallow_cth(i,j) ) ) then + z_scale = max( z_scale, zhpar(i,j) ) + zht = max( zht, zhpar(i,j) ) + end if + ! BL top includes decoupled stratocu layer, if it exists + if (ntdsc(i,j) > 0) zht = max( zht, z_uv(i,j,ntdsc(i,j)+1) ) + ! Need to restrict z_scale to dsc depth within a dsc layer + ! (given by turb_length) and to distance from dsc top below the + ! dsc layer + if ( k-1 <= ntdsc(i,j) ) then + z_scale = min( z_scale, & + max( turb_length(i,j,k), z_uv(i,j,ntdsc(i,j)+1)-zz ) ) + end if - ! Effective bulk flux Richardson number - rifb=0.3_r_bl + ! Finally calculate 1D BL weighting factor + if ( blending_option == blend_except_cu .and. & + cumulus(i,j) .and. ntdsc(i,j) == 0) then + ! pure cumulus layer so revert to 1D BL scheme + weight_1dbl(i,j,k) = one + else - ind=m_buoy-m_tau+one + if ( blending_option == blend_gridindep_fa .or. & + blending_option == blend_cth_shcu_only ) then + if (zz <= zht) then + weight_1dbl(i,j,k) = & + one - tanh( beta_bl*z_scale/delta_smag(i,j)) * & + max( zero, & + min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) + weight_bltop(i,j) = weight_1dbl(i,j,k) + else ! above PBL + ! Above the PBL top (at zht) increase weight to one smoothly + ! between zht and zfa in order to default to 1D BL when not + ! turbulent. There is some arbitrariness here but: + ! a) we want to use a physical height, to avoid grid dependence + ! b) for shallow PBLs at high resolution it seems sensible to + ! get well (a PBL depth) above the resolved PBL before + ! reverting to 1D + ! c) for deep PBLs we still want to revert to 1D reasonably + ! quickly, hence within at most 1km of zht + zfa=min( 2.0_r_bl*zht, zht+1000.0_r_bl ) + if (zz <= zfa ) then + weight_1dbl(i,j,k) = one + one_half * & + (weight_bltop(i,j) - one) * & + ( one + cos(pi*(zz-zht)/(zfa-zht)) ) + else + weight_1dbl(i,j,k) = one + end if + if ( local_fa == free_trop_layers .and. & + ri(i,j,k) < ricrit(i,j) ) then + ! Except in an elevated turbulent layer where we still use + ! the standard blending weight + z_scale = turb_length(i,j,k) + weight_1dbl(i,j,k) = & + one - tanh( beta_bl*z_scale/delta_smag(i,j)) * & + max( zero, & + min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) + end if + end if ! test on zz < zht + else + zfa=zht+1000.0_r_bl + if (zz <= zht) then + beta=beta_bl + else if (zz <= zfa) then + beta = beta_bl*(zfa-zz)/(zfa-zht) + & + beta_fa*(zz-zht)/(zfa-zht) + else + beta=beta_fa + end if + weight_1dbl(i,j,k) = & + one - tanh( beta*z_scale/delta_smag(i,j)) * max( zero, & + min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) + end if + end if + + elm(i,j,k) = elm(i,j,k)*weight_1dbl(i,j,k) + & + sqrt(rneutml_sq(i,j,k-1))*(one-weight_1dbl(i,j,k)) + elh(i,j,k) = elh(i,j,k)*weight_1dbl(i,j,k) + & + sqrt(rneutml_sq(i,j,k-1))*(one-weight_1dbl(i,j,k)) + end do + end do +!$OMP end do + end if ! test on blending_option +!$OMP end PARALLEL +end do ! loop over levels +!---------------------------------------------------------------- +! 3.0 Calculate stability functions +!---------------------------------------------------------------- +! 3.1 Set-up a BL weighting function, =1 near the ground (ie in the BL) +! =0 in the free troposphere +! Rate and height at which transition occurs varys depending on choices +!----------------------------------------------------------------------- +!$OMP PARALLEL DEFAULT(none) private( i, j, k, z_scale, zpr) & +!$OMP SHARED( pdims, bl_levels, BL_weight, local_fa, z_tq, sg_orog_mixing, & +!$OMP sigma_h) +!$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 - - ! Set diff_min to a large initial value - diff_min(i,j)=1000.0_r_bl - - ! Surface Obukhov length - MOsurf(i,j)= -v_s(i,j)*v_s(i,j)*v_s(i,j) & - /(vkman*fb_surf(i,j)) + BL_weight(i,j,k) = one end do end do +end do +!$OMP end do +if (local_fa == to_sharp_across_1km) then + !--------------------------------------------------------- + ! Additional code to allow the local Ri scheme to use + ! SHARPEST in the free atmosphere, ie above the BL top, + ! regardless of the tail option selected above. + ! Set Z_SCALE to 1km to mimic old value of BL_LEVELS, + ! gives BL_weight~0 by 2km, ~0.95 at 500m + !--------------------------------------------------------- + z_scale = 1000.0_r_bl +!$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 + zpr = z_tq(i,j,k-1)/z_scale + BL_weight(i,j,k) = one_half*(one - tanh(3.0_r_bl*(zpr-one) ) ) + end do + end do + end do +!$OMP end do NOWAIT +end if - ! The wind speed change from level k to the surface - u1=sqrt(u_p(i,j,k)*u_p(i,j,k)+v_p(i,j,k)*v_p(i,j,k)) - - ! h_est is the estimate of the stable boundary layer - ! depth using the TKE based formula - h_est=vkman*MOsurf(i,j)*ind*rifb*u1/v_s(i,j) - - ! Absolute difference between height and estimate - diff=abs(z_uv(i,j,k)-h_est) - - ! If h_est is closer than the previous closest value - ! (diff_min) reset the h_tkeb to h_est - - if (diff < diff_min(i,j)) then - diff_min(i,j)=diff - h_tkeb(i,j)=h_est +if (sg_orog_mixing /= off) then + !----------------------------------------------------------------- + ! Subgrid orographic height dependence for SBL tail (option 1) + ! or orographic dependence of mixing lengths, lambdam,h (opt 2) + ! Gives BL_weight~[1,0.95,0.5,0] at ZPR=[0,0.6,1,1.7] + !---------------------------------------------------------------- +!$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 (sigma_h(i,j) > one ) then + zpr = z_tq(i,j,k-1)/sigma_h(i,j) + BL_weight(i,j,k) = one_half*( one - tanh(4.0_r_bl*(zpr-one) ) ) end if - end do end do end do - -end if ! SBL_OP = Depth_based - +!$OMP end do NOWAIT +end if +!$OMP end PARALLEL ! ---------------------------------------------------------------- -! Main loop over levels +! 3.2 calculating stable stability function +! ---------------------------------------------------------------- +! Load up 2D array FUNC with selected stability function for Ri>=0 +! +! SBL_OP Option +! +! Long_tails Long tails +! Sharpest SHARPEST function +! Sharp_sea_long_land SHARPEST over sea ; Long tails over land +! Mes_tails MESOSCALE model: Louis/SHARPEST blend +! Louis_tails Louis function +! Sharp_sea_mes_land SHARP over sea; Mes over land +! Sharp_sea_Louis_land SHARP over sea; Louis over land ! ---------------------------------------------------------------- - do k = 2, bl_levels - ! ---------------------------------------------------------------- - ! Load up 2D array FUNC with selected stability function for Ri>=0 - - ! SBL_OP Option - - ! Long_tails Long tails - ! Sharpest SHARPEST function - ! Sharp_sea_long_land SHARPEST over sea ; Long tails over land - ! Mes_tails MESOSCALE model: Louis/SHARPEST blend - ! Louis_tails Louis function - ! Depth_based Boundary layer depth based formulation - ! Sharp_sea_mes_land SHARP over sea; Mes over land - ! Sharp_sea_Louis_land SHARP over sea; Louis over land - ! ---------------------------------------------------------------- - select case (sbl_op) !-------------------------------------------- @@ -1061,19 +1019,6 @@ subroutine ex_coef ( & end do end do - !-------------------------------------------- - ! long TAILS FOR use WITH DEPTH BASED SCHEME - !-------------------------------------------- - case (depth_based) - ! long TAILS - do j = pdims%j_start, pdims%j_end - do i = pdims%i_start, pdims%i_end - if (ri(i,j,k) >= zero) then - func(i,j)=one / ( one + g0 * ri(i,j,k) ) - end if - end do - end do - !-------------------------------------------- ! SHARP TAILS OVER SEA; MES TAILS OVER LAND !-------------------------------------------- @@ -1152,9 +1097,8 @@ subroutine ex_coef ( & end select ! SBL_OP !------------------------------------------------------------------ - ! Additional code to allow the local Ri scheme to use - ! SHARPEST in the free atmosphere, ie above the BL top, - ! regardless of the tail option selected above. + ! Additional option to use SHARPEST in the free atmosphere, ie above + ! the BL top, regardless of the tail option selected above. !------------------------------------------------------------------ if (local_fa == to_sharp_across_1km) then @@ -1214,7 +1158,7 @@ subroutine ex_coef ( & end if !--------------------------------------------------------------- - ! Set stable Prandtl number (=KM/KH) + ! 3.3 Set stable Prandtl number (=KM/KH) !--------------------------------------------------------------- if (sbl_op == lem_stability) then !$OMP PARALLEL do DEFAULT(none) SCHEDULE(STATIC) & @@ -1244,223 +1188,20 @@ subroutine ex_coef ( & end do !$OMP end PARALLEL do end if - -!$OMP PARALLEL do SCHEDULE(STATIC) DEFAULT(none) & -!$OMP private(z_scale,fm,j,i,lambdam,lambdah,lambda_eff, & -!$OMP lambdah_rho,vkz,f_log,zz,zht,zfa,beta,fh,rtmri,km,rpr) & -!$OMP SHARED(k,sbl_op,var_diags_opt,bl_levels,pdims,g0,ri,ricrit, & -!$OMP flandg,ntml_local,ntml_nl,subb,dh,z_tq, & -!$OMP BL_weight,sg_orog_mixing,sigma_h,pr_n, & -!$OMP l_rp2,lambda_min,par_mezcla_rp,zh_local, & -!$OMP h_blend,turb_length,k_log_layr, & -!$OMP z_uv,z0m,elm,elh,elh_rho,blending_option,cumulus,l_shallow_cth,zhpar, & -!$OMP ntdsc,weight_1dbl,weight_bltop,delta_smag,rneutml_sq,BL_diag, & -!$OMP r_pr_n,cbl_op,subc,dm,h_tkeb,v_s,MOsurf,rho_wet_tq, & -!$OMP l_subfilter_vert,l_subfilter_horiz,fm_3d,fh_3d,rhokm,tke_loc, & -!$OMP dvdzm,l_mr_physics,rhokh,local_fa,fb_surf,func,prandtl_number) - do j = pdims%j_start, pdims%j_end +!---------------------------------------------------------------- +! 3.4 Calculate (values of) stability functions FH, FM. +!---------------------------------------------------------------- +!$OMP PARALLEL DEFAULT(none) & +!$OMP SHARED( k, pdims,BL_diag,elm,elh,ri,func,prandtl_number,cbl_op,r_pr_n, & +!$OMP l_subfilter_vert,l_subfilter_horiz,fm_3d,fh_3d,rhokm,rhokh, & +!$OMP rho_wet_tq,dvdzm,l_mr_physics,local_fa,tke_loc,subb,subc,g0,dm, & +!$OMP dh ) & +!$OMP PRIVATE( i, j, fm, fh, rtmri, rpr ) +!$OMP do SCHEDULE(STATIC) + do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end - !----------------------------------------------------------------- - ! 2.1 Calculate asymptotic mixing lengths LAMBDAM and LAMBDAH - ! (may be equal or LambdaM=2*LambdaH (operational setting)). - !----------------------------------------------------------------- - if (l_lambdam2) then - if (l_rp2) then - lambdam = max (lambda_min , 2*par_mezcla_rp(rp_idx)*zh_local(i,j)) - lambdah = max (lambda_min , par_mezcla_rp(rp_idx)*zh_local(i,j)) - else - lambdam = max (lambda_min , 0.30_r_bl*zh_local(i,j)) - lambdah = max (lambda_min , 0.15_r_bl*zh_local(i,j)) - end if - else - if (l_rp2) then - lambdam = max ( lambda_min , par_mezcla_rp(rp_idx)*zh_local(i,j) ) - else - lambdam = max ( lambda_min , lambda_fac*zh_local(i,j) ) - end if - lambdah = lambdam - end if - - if (sg_orog_mixing == sg_shear_enh_lambda) then - !-------------------------------------------------------------- - ! Use orographic mixing length for heat too, and reduce both - ! above sigma_h smoothly - ! NOTE: THIS CODE WILL not ENHANCE LAMBDAH because it only - ! gets used in bdy_expl2 where the calculation is redone as - ! standard - this was a mistake but is now operational in - ! the UKV! - !-------------------------------------------------------------- - if (k >= ntml_local(i,j)+2) then - lambdam = lambda_min - lambdah = lambda_min - end if - lambdah = max (lambdah, & - BL_weight(i,j,k)*a_lambda*h_blend(i,j) ) - lambda_eff = max (lambdam, & - BL_weight(i,j,k)*a_lambda*h_blend(i,j) ) - else - lambda_eff = max (lambdam, a_lambda*h_blend(i,j) ) - !------------------------------------------------------------ - ! Optionally reduce mixing length above local BL top - !------------------------------------------------------------ - if (k >= ntml_local(i,j)+2 .and. .not. l_full_lambdas) then - lambdam = lambda_min - lambdah = lambda_min - if (z_tq(i,j,k-1) > a_lambda*h_blend(i,j)) & - lambda_eff=lambda_min - end if - if ( k >= ntml_local(i,j)+2 .and. l_full_lambdas .and. & - local_fa == to_sharp_across_1km ) then - ! Weight lambda to lambda_min with height - ! Assuming only local_fa == to_sharp_across_1km will have - ! L_FULL_LAMBDAS. If other LOCAL_FA options are coded here - ! then changes must be included in section 5.3 of bdy_expl2 - - lambda_eff = lambda_eff * BL_weight(i,j,k) & - + lambda_min*( one - BL_weight(i,j,k) ) - lambdah = lambdah * BL_weight(i,j,k) & - + lambda_min*( one - BL_weight(i,j,k) ) - end if - end if - - lambdah_rho = lambdah - - if ( local_fa == free_trop_layers ) then - lambda_eff = max( lambda_eff, lambda_fac*turb_length(i,j,k) ) - lambdah = max( lambdah, lambda_fac*turb_length(i,j,k) ) - ! lambdah_rho does not need to be recalculated under - ! local_fa option "free_trop_layers" as the full KH profile - ! will be interpolated in bdy_expl2 - end if - !----------------------------------------------------------------------- - ! 2.2 Calculate mixing lengths ELH, ELM coincident with RI(K) and so - ! at Z_TQ(K-1) - !----------------------------------------------------------------------- - ! Incorporate log profile corrections to the vertical finite - ! differences into the definitions of ELM and ELH. - ! Note that ELH_RHO is calculated (on rho levels) for direct inclusion - ! in RHOKH and also (as elh) on theta levels for the unstable - ! stability functions and inclusion in RHOKH before interpolation - ! (under local_fa option "free_trop_layers"). - ! To save computing logarithms for all K, the values of ELM and ELH - ! are unchanged for K > K_LOG_LAYR. - - if (k <= k_log_layr) then - vkz = vkman * ( z_uv(i,j,k) - z_uv(i,j,k-1) ) - f_log = log( ( z_uv(i,j,k) + z0m(i,j) ) / & - ( z_uv(i,j,k-1) + z0m(i,j) ) ) - elm(i,j,k) = vkz / ( f_log + vkz/lambda_eff ) - elh(i,j,k) = vkz / ( f_log + vkz/lambdah ) - vkz = vkman * ( z_tq(i,j,k) - z_tq(i,j,k-1) ) - f_log = log( ( z_tq(i,j,k) + z0m(i,j) ) / & - ( z_tq(i,j,k-1) + z0m(i,j) ) ) - elh_rho(i,j,k) = vkz / ( f_log + vkz/lambdah_rho ) - else - vkz = vkman * ( z_tq(i,j,k-1) + z0m(i,j) ) - elm(i,j,k) = vkz / (one + vkz/lambda_eff ) - elh(i,j,k) = vkz / (one + vkz/lambdah ) - vkz = vkman * ( z_uv(i,j,k) + z0m(i,j) ) - elh_rho(i,j,k) = vkz / (one + vkz/lambdah_rho ) - end if - - if (blending_option /= off) then - - zz = z_tq(i,j,k-1) ! height of rhokm(k) - ! turb_length is the greater of the local and non-local - ! BL depths up to that bl top - z_scale = max( zz, turb_length(i,j,k) ) - ! zht = interface between BL and FA - zht = max( z_uv(i,j,ntml_nl(i,j)+1) , zh_local(i,j) ) - ! Relevant scale in cumulus layers can be cloud top height, zhpar - if ( cumulus(i,j) .and. ( blending_option /= blend_cth_shcu_only .or. & - l_shallow_cth(i,j) ) ) then - z_scale = max( z_scale, zhpar(i,j) ) - zht = max( zht, zhpar(i,j) ) - end if - ! BL top includes decoupled stratocu layer, if it exists - if (ntdsc(i,j) > 0) zht = max( zht, z_uv(i,j,ntdsc(i,j)+1) ) - ! Need to restrict z_scale to dsc depth within a dsc layer - ! (given by turb_length) and to distance from dsc top below the - ! dsc layer - if ( k-1 <= ntdsc(i,j) ) then - z_scale = min( z_scale, & - max( turb_length(i,j,k), z_uv(i,j,ntdsc(i,j)+1)-zz ) ) - end if - - ! Finally calculate 1D BL weighting factor - if ( blending_option == blend_except_cu .and. & - cumulus(i,j) .and. ntdsc(i,j) == 0) then - ! pure cumulus layer so revert to 1D BL scheme - weight_1dbl(i,j,k) = one - else - - if ( blending_option == blend_gridindep_fa .or. & - blending_option == blend_cth_shcu_only ) then - if (zz <= zht) then - weight_1dbl(i,j,k) = & - one - tanh( beta_bl*z_scale/delta_smag(i,j)) * & - max( zero, & - min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) - weight_bltop(i,j) = weight_1dbl(i,j,k) - else ! above PBL - ! Above the PBL top (at zht) increase weight to one smoothly - ! between zht and zfa in order to default to 1D BL when not - ! turbulent. There is some arbitrariness here but: - ! a) we want to use a physical height, to avoid grid dependence - ! b) for shallow PBLs at high resolution it seems sensible to - ! get well (a PBL depth) above the resolved PBL before - ! reverting to 1D - ! c) for deep PBLs we still want to revert to 1D reasonably - ! quickly, hence within at most 1km of zht - zfa=min( 2.0_r_bl*zht, zht+1000.0_r_bl ) - if (zz <= zfa ) then - weight_1dbl(i,j,k) = one + one_half * & - (weight_bltop(i,j) - one) * & - ( one + cos(pi*(zz-zht)/(zfa-zht)) ) - else - weight_1dbl(i,j,k) = one - end if - if ( local_fa == free_trop_layers .and. & - ri(i,j,k) < ricrit(i,j) ) then - ! Except in an elevated turbulent layer where we still use - ! the standard blending weight - z_scale = turb_length(i,j,k) - weight_1dbl(i,j,k) = & - one - tanh( beta_bl*z_scale/delta_smag(i,j)) * & - max( zero, & - min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) - - end if - end if ! test on zz < zht - else - zfa=zht+1000.0_r_bl - if (zz <= zht) then - beta=beta_bl - else if (zz <= zfa) then - beta = beta_bl*(zfa-zz)/(zfa-zht) + & - beta_fa*(zz-zht)/(zfa-zht) - else - beta=beta_fa - end if - weight_1dbl(i,j,k) = & - one - tanh( beta*z_scale/delta_smag(i,j)) * max( zero, & - min( one, (linear0-delta_smag(i,j)/z_scale)*rlinfac) ) - end if - end if - - elm(i,j,k) = elm(i,j,k)*weight_1dbl(i,j,k) + & - sqrt(rneutml_sq(i,j,k-1))*(one-weight_1dbl(i,j,k)) - elh(i,j,k) = elh(i,j,k)*weight_1dbl(i,j,k) + & - sqrt(rneutml_sq(i,j,k-1))*(one-weight_1dbl(i,j,k)) - - end if ! test on blending_option if (BL_diag%l_elm3d) BL_diag%elm3d(i,j,k)=elm(i,j,k) - - !---------------------------------------------------------------- - ! 2.4 Calculate (values of) stability functions FH, FM. - !---------------------------------------------------------------- - if (ri(i,j,k) >= zero) then !----------------------------------------------------------- ! Note that we choose to include the Pr dependence such that @@ -1489,7 +1230,7 @@ subroutine ex_coef ( & end if !------------------------------------------------------------------ - ! 2.5_r_bl Calculate exchange coefficients RHO*KM(K), RHO*KH(K) + ! 4.0 Calculate exchange coefficients RHO*KM(K), RHO*KH(K) ! both on TH-level K-1 at this stage (RHOKH will be interpolated ! onto uv-levels and then be multiplied by ELH in BDY_EXPL2 if ! local_fa is not "free_trop_layers") @@ -1517,7 +1258,7 @@ subroutine ex_coef ( & if (local_fa == free_trop_layers) & rhokh(i,j,k) = rhokh(i,j,k) * elh(i,j,k) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + if (BL_diag%l_tke) then rpr = fh / max(fm, tiny(one) ) tke_loc(i,j,k) = ( r_c_tke*elm(i,j,k)*dvdzm(i,j,k)*dvdzm(i,j,k) & *(rhokm(i,j,k)/rho_wet_tq(i,j,k-1)) & @@ -1525,47 +1266,13 @@ subroutine ex_coef ( & one - ri(i,j,k)*rpr ) ) & )**two_thirds end if - ! ------------------------------------------- - ! Boundary layer depth based formulation - ! ------------------------------------------- - - if (sbl_op == depth_based .and. & - fb_surf(i,j) <= zero) then - if (z_tq(i,j,k-1) < h_tkeb(i,j)) then - - ! Formula for diffusion coefficient - ! see Beare et al 2006, Boundary layer Met. - - km = v_s(i,j) * vkman * z_tq(i,j,k-1) * & - ( (one-z_tq(i,j,k-1)/h_tkeb(i,j))**(1.5_r_bl) ) & - / (one + 4.7_r_bl*z_tq(i,j,k-1)/MOsurf(i,j)) - rhokm(i,j,k)= rho_wet_tq(i,j,k-1) * km - if (l_mr_physics) then - ! Note "RHO" here is always wet density (RHO_WET_TQ) so - ! save multiplication of RHOKH to after interpolation - rhokh(i,j,k)= km*r_pr_n / elm(i,j,k) - else - rhokh(i,j,k) = rho_wet_tq(i,j,k-1)*km*r_pr_n /elm(i,j,k) - end if - else - rhokm(i,j,k)=zero - rhokh(i,j,k)=zero - end if - end if !SBL_OP == Depth_based - - end do !I + end do !i end do !j -!$OMP end PARALLEL do +!$OMP end do +!$OMP end PARALLEL end do ! bl_levels -!----------------------------------------------------------------------- -! 3. Equilibrium Stable Boundary Layer (SBL) model. -!----------------------------------------------------------------------- -!----------------------------------------------------------------------- -! Finish up -!----------------------------------------------------------------------- - if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) return end subroutine ex_coef 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..a708083c2 100644 --- a/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/excf_nl_9c.F90 @@ -56,7 +56,7 @@ subroutine excf_nl_9c ( & flux_grad, Locketal2000, HoltBov1993, LockWhelan2006, entr_smooth_dec, & entr_taper_zh, kprof_cu, klcl_entr, buoy_integ, buoy_integ_low, & max_cu_depth, bl_res_inv, a_ent_shr_nml, a_ent_2, one_third, two_thirds, & - l_reset_dec_thres, var_diags_opt, original_vars, split_tke_and_inv, & + l_reset_dec_thres, & l_use_var_fixes, dzrad_disc_opt, dzrad_ntm1, dzrad_1p5dz, & num_sweeps_bflux, l_converge_ga, l_use_sml_dsc_fixes, zero, one, one_half use model_domain_mod, only: model_type, mt_single_column @@ -2990,7 +2990,7 @@ subroutine excf_nl_9c ( & ( (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 + if (BL_diag%l_tke) 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 @@ -3054,26 +3054,13 @@ subroutine excf_nl_9c ( & 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 + if (BL_diag%l_tke) 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 - 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 @@ -3183,19 +3170,12 @@ subroutine excf_nl_9c ( & ( 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 + if (BL_diag%l_tke) 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 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 if end do end do @@ -3333,7 +3313,7 @@ subroutine excf_nl_9c ( & ( 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) 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) @@ -3345,20 +3325,13 @@ subroutine excf_nl_9c ( & 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 + if (BL_diag%l_tke) 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 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 if end do diff --git a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 index 24008fac1..32fb1559c 100644 --- a/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 +++ b/science/physics_schemes/source/boundary_layer/kmkhz_9c.F90 @@ -57,8 +57,7 @@ subroutine kmkhz_9c ( & a_grad_adj, max_t_grad, flux_grad, Locketal2000, & HoltBov1993, LockWhelan2006, entr_smooth_dec, entr_taper_zh, & kprof_cu, l_use_sml_dsc_fixes, l_converge_ga, & - bl_res_inv, cosine_inv_flux, target_inv_profile, pr_max, var_diags_opt, & - split_tke_and_inv, l_noice_in_turb, & + bl_res_inv, cosine_inv_flux, target_inv_profile, pr_max, l_noice_in_turb, & one_third, two_thirds, zero, one, one_half use conversions_mod, only: pi => pi_bl use cv_run_mod, only: l_param_conv @@ -4165,7 +4164,7 @@ subroutine kmkhz_9c ( & if (res_inv(i,j) == 1) then Prandtl = min( rhokm(i,j,k)/(rbl_eps+rhokh_surf_ent(i,j)), & pr_max ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + if (BL_diag%l_tke) then ! need velocity scale for TKE diagnostic w_m = ( v_s(i,j)*v_s(i,j)*v_s(i,j) + & c_ws * zh(i,j) * fb_surf(i,j) ) ** one_third @@ -4189,7 +4188,7 @@ subroutine kmkhz_9c ( & * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + if (BL_diag%l_tke) then ! save Km/timescale for TKE diag, completed in bdy_expl2 tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) end if @@ -4237,7 +4236,7 @@ subroutine kmkhz_9c ( & * rdz(i,j,kl) * (z_uv(i,j,kl)-z_uv(i,j,kl-1)) & * rho_wet_tq(i,j,kl-1) / rho_mix(i,j,kl) rhokm(i,j,kl) = max( rhokm(i,j,kl), rhok_inv ) - if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then + if (BL_diag%l_tke) then ! save Km/timescale for TKE diag, completed in bdy_expl2 tke_nl(i,j,kl) = max( tke_nl(i,j,kl), rhok_inv*c_tke*w_m/zh(i,j)) end if @@ -4702,7 +4701,7 @@ subroutine kmkhz_9c ( & ! Estimate turbulent w-variance scale at discontinuous inversions !----------------------------------------------------------------------- -if (BL_diag%l_tke .and. var_diags_opt == split_tke_and_inv) then +if (BL_diag%l_tke) then !$OMP do SCHEDULE(STATIC) do j = pdims%j_start, pdims%j_end do i = pdims%i_start, pdims%i_end