From 96923336ef754917c0c37f2d632e71ff1d65a502 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 8 Jan 2025 07:40:50 -0500 Subject: [PATCH 1/4] +Rename visc%TKE_BBL to visc%BBL_meanKE_loss Revised the name of visc%TKE_BBL to visc%BBL_meanKE_loss to better reflect what these variables contain, and the fact that the efficiency of the conversion from mean kinetic energy loss to turbulent kinetic energy has not been applied yet. The units of visc%BBL_meanKE_loss are [H L2 T-3 ~> m3 s-3 or W m-2], whereas those of visc%TKE_BBL were [H Z2 T-3 ~> m3 s-3 or W m-2]. The factor rescaling between the units of mean kinetic energy and those of turbulent kinetic energy have been incorporated into set_diffusivity_CS%BBL_effic and energetic_PBL_CS%ePBL_BBL_effic. All answers are bitwise identical. --- src/core/MOM_variables.F90 | 4 +- .../vertical/MOM_energetic_PBL.F90 | 13 ++++--- .../vertical/MOM_set_diffusivity.F90 | 38 ++++++++++--------- .../vertical/MOM_set_viscosity.F90 | 4 +- 4 files changed, 33 insertions(+), 26 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 65e915705a..efd558cdb3 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -237,8 +237,8 @@ module MOM_variables kv_bbl_v, & !< The bottom boundary layer viscosity at the v-points [H Z T-1 ~> m2 s-1 or Pa s] ustar_BBL, & !< The turbulence velocity in the bottom boundary layer at !! h points [H T-1 ~> m s-1 or kg m-2 s-1]. - TKE_BBL, & !< A term related to the bottom boundary layer source of turbulent kinetic - !! energy, currently in [H Z2 T-3 ~> m3 s-3 or W m-2]. + BBL_meanKE_loss, & !< The viscous loss of mean kinetic energy in the bottom boundary layer + !! [H L2 T-3 ~> m3 s-3 or W m-2]. taux_shelf, & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa]. tauy_shelf !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa]. real, allocatable, dimension(:,:) :: tbl_thick_shelf_u diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index e32bc8de2d..fdbd9c13f0 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -162,7 +162,9 @@ module MOM_energetic_PBL real :: Max_Enhance_M = 5. !< The maximum allowed LT enhancement to the mixing [nondim]. !/ Bottom boundary layer mixing related options - real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL [nondim] + real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL, times + !! conversion factors between the natural units of mean kinetic energy + !! and those those used for TKE [Z2 L-2 ~> nondim]. logical :: Use_BBLD_iteration !< If true, use the proximity to the top of the actively turbulent !! bottom boundary layer to constrain the mixing lengths. real :: TKE_decay_BBL !< The ratio of the natural Ekman depth to the TKE decay scale for @@ -654,7 +656,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (CS%ePBL_BBL_effic > 0.0) then if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j) BBLD_in = BBLD_io - BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) @@ -730,7 +732,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, mixvel, mixlen, GV, US, CS_tmp2, eCD_tmp, Waves, G, i, j) else BLD_1 = BBLD_in ; BLD_2 = BBLD_in - BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & u_star_BBL, Kd_1, BLD_1, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp1, eCD_tmp) @@ -760,7 +762,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, enddo ! j-loop if (CS%debug .and. (CS%ePBL_BBL_effic > 0.0)) then - call hchksum(visc%TKE_BBL, "ePBL visc%TKE_BBL", G%HI, unscale=GV%H_to_MKS*US%Z_to_m**2*US%s_to_T**3) + call hchksum(visc%BBL_meanKE_loss, "ePBL visc%BBL_meanKE_loss", G%HI, & + unscale=GV%H_to_MKS*US%L_T_to_m_s**2*US%s_to_T) call hchksum(visc%ustar_BBL, "ePBL visc%ustar_BBL", G%HI, unscale=GV%H_to_MKS*US%s_to_T) call hchksum(Kd_int, "End of ePBL Kd_int", G%HI, unscale=GV%H_to_MKS*US%Z_to_m*US%s_to_T) call hchksum(diag_Velocity_Scale, "ePBL Velocity_Scale", G%HI, unscale=US%Z_to_m*US%s_to_T) @@ -3670,7 +3673,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, & "The efficiency of bottom boundary layer mixing via ePBL. Setting this to a "//& "value that is greater than 0 to enable bottom boundary layer mixing from EPBL.", & - units="nondim", default=0.0) + units="nondim", default=0.0, scale=US%L_to_Z**2) no_BBL = (CS%ePBL_BBL_effic<=0.0) call get_param(param_file, mdl, "USE_BBLD_ITERATION", CS%Use_BBLD_iteration, & "A logical that specifies whether or not to use the distance to the top of the "//& diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 2f54eb2b86..e4de7ee081 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -75,8 +75,10 @@ module MOM_set_diffusivity logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity. real :: Von_Karm !< The von Karman constant as used in the BBL diffusivity calculation !! [nondim]. See (http://en.wikipedia.org/wiki/Von_Karman_constant) - real :: BBL_effic !< efficiency with which the energy extracted - !! by bottom drag drives BBL diffusion in the original BBL scheme [nondim] + real :: BBL_effic !< Efficiency with which the energy extracted + !! by bottom drag drives BBL diffusion in the original BBL scheme, times + !! conversion factors between the natural units of mean kinetic energy + !! and those those used for TKE [Z2 L-2 ~> nondim]. real :: ePBL_BBL_effic !< efficiency with which the energy extracted !! by bottom drag drives BBL diffusion in the ePBL BBL scheme [nondim] real :: cdrag !< quadratic drag coefficient [nondim] @@ -1412,10 +1414,10 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, ! If ustar_h = 0, this is land so this value doesn't matter. I2decay(i) = 0.5*CS%IMax_decay endif - TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%TKE_BBL(i,j) + TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j) if (associated(fluxes%TKE_tidal)) & - TKE(i) = TKE(i) + fluxes%TKE_tidal(i,j) * GV%RZ_to_H * & + TKE(i) = TKE(i) + US%Z_to_L**2*fluxes%TKE_tidal(i,j) * GV%RZ_to_H * & (CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz))) ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following @@ -1471,7 +1473,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, else ; TKE_to_layer = 0.0 ; endif ! TKE_Ray has been initialized to 0 above. - if (Rayleigh_drag) TKE_Ray = 0.5*CS%BBL_effic * US%L_to_Z**2 * G%IareaT(i,j) * & + if (Rayleigh_drag) TKE_Ray = 0.5*CS%BBL_effic * G%IareaT(i,j) * & (((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + & (G%areaCu(I,j) * visc%Ray_u(I,j,k) * u(I,j,k)**2)) + & ((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + & @@ -1576,6 +1578,8 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo real :: dz(SZI_(G),SZK_(GV)) ! Height change across layers [Z ~> m] real :: dz_above(SZK_(GV)+1) ! Distance from each interface to the surface [Z ~> m] real :: TKE_column ! net TKE input into the column [H Z2 T-3 ~> m3 s-3 or W m-2] + real :: BBL_meanKE_dis ! Sum of tidal and mean kinetic energy dissipation in the bottom boundary layer, which + ! can act as a source of TKE [H L2 T-3 ~> m3 s-3 or W m-2] real :: TKE_remaining ! remaining TKE available for mixing in this layer and above [H Z2 T-3 ~> m3 s-3 or W m-2] real :: TKE_consumed ! TKE used for mixing in this layer [H Z2 T-3 ~> m3 s-3 or W m-2] real :: TKE_Kd_wall ! TKE associated with unlimited law of the wall mixing [H Z2 T-3 ~> m3 s-3 or W m-2] @@ -1639,14 +1643,14 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo if ((ustar > 0.0) .and. (absf > CS%IMax_decay * ustar)) Idecay = absf / ustar ! Energy input at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2]. - ! (Note that visc%TKE_BBL is in [H Z2 T-3 ~> m3 s-3 or W m-2], set in set_BBL_TKE().) - ! I am still unsure about sqrt(cdrag) in this expressions - AJA - TKE_column = cdrag_sqrt * visc%TKE_BBL(i,j) + ! (Note that visc%BBL_meanKE_loss is in [H L2 T-3 ~> m3 s-3 or W m-2], set in set_BBL_TKE().) + !### I am still unsure about sqrt(cdrag) in this expressions - AJA + BBL_meanKE_dis = cdrag_sqrt * visc%BBL_meanKE_loss(i,j) ! Add in tidal dissipation energy at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2]. ! Note that TKE_tidal is in [R Z3 T-3 ~> W m-2]. if (associated(fluxes%TKE_tidal)) & - TKE_column = TKE_column + fluxes%TKE_tidal(i,j) * GV%RZ_to_H - TKE_column = CS%BBL_effic * TKE_column ! Only use a fraction of the mechanical dissipation for mixing. + BBL_meanKE_dis = BBL_meanKE_dis + US%Z_to_L**2*fluxes%TKE_tidal(i,j) * GV%RZ_to_H + TKE_column = CS%BBL_effic * BBL_meanKE_dis ! Only use a fraction of the mechanical dissipation for mixing. TKE_remaining = TKE_column if (CS%LOTW_BBL_answer_date > 20240630) then @@ -1670,7 +1674,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo ! Add in additional energy input from bottom-drag against slopes (sides) if (Rayleigh_drag) TKE_remaining = TKE_remaining + & - 0.5*CS%BBL_effic * US%L_to_Z**2 * G%IareaT(i,j) * & + 0.5*CS%BBL_effic * G%IareaT(i,j) * & (((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + & (G%areaCu(I,j) * visc%Ray_u(I,j,k) * u(I,j,k)**2)) + & ((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + & @@ -1949,8 +1953,8 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (allocated(visc%ustar_BBL)) then do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo endif - if (allocated(visc%TKE_BBL)) then - do j=js,je ; do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ; enddo ; enddo + if (allocated(visc%BBL_meanKE_loss)) then + do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss(i,j) = 0.0 ; enddo ; enddo endif return endif @@ -2075,7 +2079,7 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) (G%areaCu(I,j)*(ustar(I)*ustar(I)))) + & ((G%areaCv(i,J-1)*(vstar(i,J-1)*vstar(i,J-1))) + & (G%areaCv(i,J)*(vstar(i,J)*vstar(i,J)))) ) ) - visc%TKE_BBL(i,j) = US%L_to_Z**2 * & + visc%BBL_meanKE_loss(i,j) = & ((((G%areaCu(I-1,j)*(ustar(I-1)*u2_bbl(I-1))) + & (G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + & ((G%areaCv(i,J-1)*(vstar(i,J-1)*v2_bbl(i,J-1))) + & @@ -2358,9 +2362,9 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "velocity field to the bottom stress. CDRAG is only used "//& "if BOTTOMDRAGLAW is true.", units="nondim", default=0.003) call get_param(param_file, mdl, "BBL_EFFIC", CS%BBL_effic, & - "The efficiency with which the energy extracted by "//& - "bottom drag drives BBL diffusion. This is only "//& - "used if BOTTOMDRAGLAW is true.", units="nondim", default=0.20) + "The efficiency with which the energy extracted by bottom drag drives BBL "//& + "diffusion. This is only used if BOTTOMDRAGLAW is true.", & + units="nondim", default=0.20, scale=US%L_to_Z**2) call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, & units="nondim", default=0.0,do_not_log=.true.) call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index b77483456c..bebdc4e7d0 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -3139,7 +3139,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS allocate(visc%kv_bbl_u(IsdB:IedB,jsd:jed), source=0.0) allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0) allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0) - allocate(visc%TKE_bbl(isd:ied,jsd:jed), source=0.0) + allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0) CS%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', & diag%axesCu1, Time, 'BBL thickness at u points', 'm', conversion=US%Z_to_m) @@ -3214,7 +3214,7 @@ subroutine set_visc_end(visc, CS) if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear) if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu) if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) - if (allocated(visc%TKE_bbl)) deallocate(visc%TKE_bbl) + if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss) if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf) if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf) if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u) From 2181e38f113c779e20c798825868742715a90d36 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 8 Jan 2025 12:50:22 -0500 Subject: [PATCH 2/4] +Rename fluxes%TKE_tidal to fluxes%BBL_tidal_dis Revised the name of fluxes%TKE_tidal to fluxes%BBL_tidal_dis to better reflect what this field holds,and the fact that the efficiency of the conversion from mean kinetic energy loss to turbulent kinetic energy has not been applied yet. The units of fluxes%BBL_tidal_dis are [R Z L2 T-3 ~> W m-2], whereas those of fluxes%TKE_tidal were [R Z3 T-3 ~> W m-2]. The factor rescaling between the units of mean kinetic energy and those of turbulent kinetic energy were already in set_diffusivity_CS%BBL_effic, and these have been cancelled out by this change, but this is offset by the addition of rescaling factors in the term setting this array in the NUOPC and mct convert_IOB_to_fluxes routines. In the FMS_cap version of convert_IOB_to_fluxes, the extra rescaling factor is rolled into the scaling factor used in the get_param call for rho_TKE_tidal. All answers are bitwise identical, but there is a change in the name and rescaled units of an element of a transparent type. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 29 ++++++++++-------- .../STALE_mct_cap/mom_surface_forcing_mct.F90 | 4 +-- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 4 +-- src/core/MOM_forcing_type.F90 | 30 ++++++++++--------- .../vertical/MOM_set_diffusivity.F90 | 10 +++---- 5 files changed, 42 insertions(+), 35 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index b686c59a1f..f59a7c7fe1 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -88,8 +88,8 @@ module MOM_surface_forcing_gfdl real :: gust_const !< Constant unresolved background gustiness for ustar [R L Z T-2 ~> Pa] logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied from an input file. real, pointer, dimension(:,:) :: & - TKE_tidal => NULL() !< Turbulent kinetic energy introduced to the bottom boundary layer - !! by drag on the tidal flows [R Z3 T-3 ~> W m-2]. + BBL_tidal_dis => NULL() !< Tidal energy dissipation in the bottom boundary layer that can act as a + !! source of energy for bottom boundary layer mixing [R Z L2 T-3 ~> W m-2] real, pointer, dimension(:,:) :: & gust => NULL() !< A spatially varying unresolved background gustiness that !! contributes to ustar [R L Z T-2 ~> Pa]. gust is used when read_gust_2d is true. @@ -302,7 +302,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) - call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%BBL_tidal_dis,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) @@ -311,7 +311,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed) do j=js-2,je+2 ; do i=is-2,ie+2 - fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) + fluxes%BBL_tidal_dis(i,j) = CS%BBL_tidal_dis(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) enddo ; enddo @@ -1304,11 +1304,15 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) ! Local variables real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJ_(G)) :: & + utide_2d ! A 2d array of RMS tidal velocities [Z T-1 ~> m s-1]. real :: Flux_const_dflt ! A default piston velocity for restoring surface properties [m day-1] logical :: Boussinesq ! If true, this run is fully Boussinesq logical :: semi_Boussinesq ! If true, this run is partially non-Boussinesq - real :: rho_TKE_tidal ! The constant bottom density used to translate tidal amplitudes into the - ! tidal bottom TKE input used with INT_TIDE_DISSIPATION [R ~> kg m-3] + real :: rho_TKE_tidal ! The constant bottom density used to translate tidal amplitudes into + ! the tidal bottom TKE input used with INT_TIDE_DISSIPATION, times the + ! factor rescaling from the units of TKE to those of mean kinetic + ! energy [R L2 Z-2 ~> kg m-3] logical :: new_sim ! False if this simulation was started from a restart file ! or other equivalent files. logical :: iceberg_flux_diags ! If true, diagnostics of fluxes from icebergs are available. @@ -1573,27 +1577,28 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) call get_param(param_file, mdl, "TKE_TIDAL_RHO", rho_TKE_tidal, & "The constant bottom density used to translate tidal amplitudes into the tidal "//& "bottom TKE input used with INT_TIDE_DISSIPATION.", & - units="kg m-3", default=CS%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, & + units="kg m-3", default=CS%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R*US%Z_to_L**2, & do_not_log=.not.(CS%read_TIDEAMP.or.(CS%utide>0.0))) - call safe_alloc_ptr(CS%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%BBL_tidal_dis,isd,ied,jsd,jed) call safe_alloc_ptr(CS%ustar_tidal,isd,ied,jsd,jed) if (CS%read_TIDEAMP) then TideAmp_file = trim(CS%inputdir) // trim(TideAmp_file) ! NOTE: There are certain cases where FMS is unable to read this file, so ! we use read_netCDF_data in place of MOM_read_data. - call read_netCDF_data(TideAmp_file, 'tideamp', CS%TKE_tidal, G%Domain, & + utide_2d(:,:) = 0.0 + call read_netCDF_data(TideAmp_file, 'tideamp', utide_2d, G%Domain, & rescale=US%m_to_Z*US%T_to_s) do j=jsd, jed; do i=isd, ied - utide = CS%TKE_tidal(i,j) - CS%TKE_tidal(i,j) = G%mask2dT(i,j)*rho_TKE_tidal*CS%cd_tides*(utide*utide*utide) + utide = utide_2d(i,j) + CS%BBL_tidal_dis(i,j) = G%mask2dT(i,j)*rho_TKE_tidal*CS%cd_tides*(utide*utide*utide) CS%ustar_tidal(i,j) = sqrt(CS%cd_tides)*utide enddo ; enddo else do j=jsd,jed; do i=isd,ied utide = CS%utide - CS%TKE_tidal(i,j) = rho_TKE_tidal*CS%cd_tides*(utide*utide*utide) + CS%BBL_tidal_dis(i,j) = rho_TKE_tidal*CS%cd_tides*(utide*utide*utide) CS%ustar_tidal(i,j) = sqrt(CS%cd_tides)*utide enddo ; enddo endif diff --git a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 index 720046d517..ab3964e626 100644 --- a/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/STALE_mct_cap/mom_surface_forcing_mct.F90 @@ -295,7 +295,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) - call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%BBL_tidal_dis,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed) if (CS%allow_flux_adjustments) then @@ -304,7 +304,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, endif do j=js-2,je+2 ; do i=is-2,ie+2 - fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) + fluxes%BBL_tidal_dis(i,j) = US%Z_to_L**2*CS%TKE_tidal(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) enddo ; enddo diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 5d09c58917..61afee89ba 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -317,7 +317,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) - call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%BBL_tidal_dis,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed) if (CS%allow_flux_adjustments) then @@ -326,7 +326,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, endif do j=js-2,je+2 ; do i=is-2,ie+2 - fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) + fluxes%BBL_tidal_dis(i,j) = US%Z_to_L**2*CS%TKE_tidal(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) enddo ; enddo diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 494e36e966..2433c3c5c3 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -168,7 +168,8 @@ module MOM_forcing_type ! tide related inputs real, pointer, dimension(:,:) :: & - TKE_tidal => NULL(), & !< tidal energy source driving mixing in bottom boundary layer [R Z3 T-3 ~> W m-2] + BBL_tidal_dis => NULL(), & !< Tidal energy dissipation in the bottom boundary layer that can act + !! as a source of energy for bottom boundary layer mixing [R Z L2 T-3 ~> W m-2] ustar_tidal => NULL() !< tidal contribution to bottom ustar [Z T-1 ~> m s-1] ! iceberg related inputs @@ -1313,8 +1314,9 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) - if (associated(fluxes%TKE_tidal)) & - call hchksum(fluxes%TKE_tidal, mesg//" fluxes%TKE_tidal", G%HI, haloshift=hshift, unscale=US%RZ3_T3_to_W_m2) + if (associated(fluxes%BBL_tidal_dis)) & + call hchksum(fluxes%BBL_tidal_dis, mesg//" fluxes%BBL_tidal_dis", G%HI, haloshift=hshift, & + unscale=US%L_to_Z**2*US%RZ3_T3_to_W_m2) if (associated(fluxes%ustar_tidal)) & call hchksum(fluxes%ustar_tidal, mesg//" fluxes%ustar_tidal", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(fluxes%lrunoff)) & @@ -1438,7 +1440,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg) call locMsg(fluxes%seaice_melt_heat,'seaice_melt_heat') call locMsg(fluxes%p_surf,'p_surf') call locMsg(fluxes%salt_flux,'salt_flux') - call locMsg(fluxes%TKE_tidal,'TKE_tidal') + call locMsg(fluxes%BBL_tidal_dis,'BBL_tidal_dis') call locMsg(fluxes%ustar_tidal,'ustar_tidal') call locMsg(fluxes%lrunoff,'lrunoff') call locMsg(fluxes%frunoff,'frunoff') @@ -1546,7 +1548,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_standard_name='sea_water_pressure_at_sea_water_surface') handles%id_TKE_tidal = register_diag_field('ocean_model', 'TKE_tidal', diag%axesT1, Time, & - 'Tidal source of BBL mixing', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + 'Tidal source of BBL mixing', 'W m-2', conversion=US%L_to_Z**2*US%RZ3_T3_to_W_m2) if (.not. use_temperature) then handles%id_buoy = register_diag_field('ocean_model', 'buoy', diag%axesT1, Time, & @@ -3174,8 +3176,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_psurf > 0) .and. associated(fluxes%p_surf)) & call post_data(handles%id_psurf, fluxes%p_surf, diag) - if ((handles%id_TKE_tidal > 0) .and. associated(fluxes%TKE_tidal)) & - call post_data(handles%id_TKE_tidal, fluxes%TKE_tidal, diag) + if ((handles%id_TKE_tidal > 0) .and. associated(fluxes%BBL_tidal_dis)) & + call post_data(handles%id_TKE_tidal, fluxes%BBL_tidal_dis, diag) if ((handles%id_buoy > 0) .and. associated(fluxes%buoy)) & call post_data(handles%id_buoy, fluxes%buoy, diag) @@ -3362,8 +3364,8 @@ subroutine allocate_forcing_by_ref(fluxes_ref, G, fluxes, turns) call myAlloc(fluxes%buoy, G%isd, G%ied, G%jsd, G%jed, & associated(fluxes_ref%buoy)) - call myAlloc(fluxes%TKE_tidal, G%isd, G%ied, G%jsd, G%jed, & - associated(fluxes_ref%TKE_tidal)) + call myAlloc(fluxes%BBL_tidal_dis, G%isd, G%ied, G%jsd, G%jed, & + associated(fluxes_ref%BBL_tidal_dis)) call myAlloc(fluxes%ustar_tidal, G%isd, G%ied, G%jsd, G%jed, & associated(fluxes_ref%ustar_tidal)) @@ -3580,7 +3582,7 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%salt_flux)) deallocate(fluxes%salt_flux) if (associated(fluxes%p_surf_full)) deallocate(fluxes%p_surf_full) if (associated(fluxes%p_surf)) deallocate(fluxes%p_surf) - if (associated(fluxes%TKE_tidal)) deallocate(fluxes%TKE_tidal) + if (associated(fluxes%BBL_tidal_dis)) deallocate(fluxes%BBL_tidal_dis) if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal) if (associated(fluxes%ustar_shelf)) deallocate(fluxes%ustar_shelf) if (associated(fluxes%iceshelf_melt)) deallocate(fluxes%iceshelf_melt) @@ -3727,8 +3729,8 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns) if (associated(fluxes_in%buoy)) & call rotate_array(fluxes_in%buoy, turns, fluxes%buoy) - if (associated(fluxes_in%TKE_tidal)) & - call rotate_array(fluxes_in%TKE_tidal, turns, fluxes%TKE_tidal) + if (associated(fluxes_in%BBL_tidal_dis)) & + call rotate_array(fluxes_in%BBL_tidal_dis, turns, fluxes%BBL_tidal_dis) if (associated(fluxes_in%ustar_tidal)) & call rotate_array(fluxes_in%ustar_tidal, turns, fluxes%ustar_tidal) @@ -3999,8 +4001,8 @@ subroutine homogenize_forcing(fluxes, G, GV, US) if (associated(fluxes%buoy)) & call homogenize_field_t(fluxes%buoy, G, tmp_scale=US%L_to_m**2*US%s_to_T**3) - if (associated(fluxes%TKE_tidal)) & - call homogenize_field_t(fluxes%TKE_tidal, G, tmp_scale=US%RZ3_T3_to_W_m2) + if (associated(fluxes%BBL_tidal_dis)) & + call homogenize_field_t(fluxes%BBL_tidal_dis, G, tmp_scale=US%L_to_Z**2*US%RZ3_T3_to_W_m2) if (associated(fluxes%ustar_tidal)) & call homogenize_field_t(fluxes%ustar_tidal, G, tmp_scale=US%Z_to_m*US%s_to_T) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index e4de7ee081..3e5b2d44c8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -1416,8 +1416,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, endif TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j) - if (associated(fluxes%TKE_tidal)) & - TKE(i) = TKE(i) + US%Z_to_L**2*fluxes%TKE_tidal(i,j) * GV%RZ_to_H * & + if (associated(fluxes%BBL_tidal_dis)) & + TKE(i) = TKE(i) + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H * & (CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz))) ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following @@ -1647,9 +1647,9 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo !### I am still unsure about sqrt(cdrag) in this expressions - AJA BBL_meanKE_dis = cdrag_sqrt * visc%BBL_meanKE_loss(i,j) ! Add in tidal dissipation energy at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2]. - ! Note that TKE_tidal is in [R Z3 T-3 ~> W m-2]. - if (associated(fluxes%TKE_tidal)) & - BBL_meanKE_dis = BBL_meanKE_dis + US%Z_to_L**2*fluxes%TKE_tidal(i,j) * GV%RZ_to_H + ! Note that BBL_tidal_dis is in [R Z L2 T-3 ~> W m-2]. + if (associated(fluxes%BBL_tidal_dis)) & + BBL_meanKE_dis = BBL_meanKE_dis + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H TKE_column = CS%BBL_effic * BBL_meanKE_dis ! Only use a fraction of the mechanical dissipation for mixing. TKE_remaining = TKE_column From d55fa2b92710de2b0be8bf89eb1e6f3c5296b48e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 19 Jan 2025 11:50:52 -0500 Subject: [PATCH 3/4] +Add EPBL_BBL_TIDAL_EFFIC Added the option to use the energy source in fluxes%BBL_tidal_dis to drive mixing in ePBL_BBL_column(), similarly to what is done in add_drag_diffusivity() and add_LOTW_BBL_diffusivity(). Because this was omitted when ePBL_BBL_column() was first added, a separate runtime parameter, EPBL_BBL_TIDAL_EFFIC, is used to determine the extent to which this tidal energy source is used. The default for EPBL_BBL_TIDAL_EFFIC is currently set to 0 to avoid changing answers, but perhaps it should be changed follow the value of EPBL_BBL_EFFIC. By default all answers are bitwise identical, but there is a new runtime parameter in the MOM_parameter_doc files for cases that use ePBL. --- .../vertical/MOM_energetic_PBL.F90 | 54 +++++++++++++------ 1 file changed, 37 insertions(+), 17 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index fdbd9c13f0..6e3cdfa1a5 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -162,9 +162,14 @@ module MOM_energetic_PBL real :: Max_Enhance_M = 5. !< The maximum allowed LT enhancement to the mixing [nondim]. !/ Bottom boundary layer mixing related options - real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL, times + real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL driven by + !! the bottom drag dissipation of mean kinetic energy, times !! conversion factors between the natural units of mean kinetic energy - !! and those those used for TKE [Z2 L-2 ~> nondim]. + !! and those used for TKE [Z2 L-2 ~> nondim]. + real :: ePBL_tidal_effic !< The efficiency of bottom boundary layer mixing via ePBL driven by + !! the bottom drag dissipation of tides, times conversion factors + !! between the natural units of mean kinetic energy and those used for + !! TKE [Z2 L-2 ~> nondim]. logical :: Use_BBLD_iteration !< If true, use the proximity to the top of the actively turbulent !! bottom boundary layer to constrain the mixing lengths. real :: TKE_decay_BBL !< The ratio of the natural Ekman depth to the TKE decay scale for @@ -464,7 +469,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, type(ePBL_column_diags) :: eCD_tmp ! A container for not passing around diagnostics. type(energetic_PBL_CS) :: CS_tmp1, CS_tmp2 ! Copies of the energetic PBL control structure that ! can be modified to test for sensitivities - + logical :: BBL_mixing ! If true, there is bottom boundary layer mixing. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -487,6 +492,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 ! This is not used when fully non-Boussinesq. I_dt = 0.0 ; if (dt > 0.0) I_dt = 1.0 / dt I_rho0dt = 1.0 / (GV%Rho0 * dt) ! This is not used when fully non-Boussinesq. + BBL_mixing = ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0)) ! Zero out diagnostics before accumulation. if (CS%TKE_diagnostics) then @@ -497,7 +503,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, diag_TKE_mixing(i,j) = 0.0 ; diag_TKE_mech_decay(i,j) = 0.0 diag_TKE_conv_decay(i,j) = 0.0 !; diag_TKE_unbalanced(i,j) = 0.0 enddo ; enddo - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then !!OMP parallel do default(shared) do j=js,je ; do i=is,ie diag_TKE_BBL(i,j) = 0.0 ; diag_TKE_BBL_mixing(i,j) = 0.0 @@ -507,7 +513,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, endif if (CS%debug .or. (CS%id_Mixing_Length>0)) diag_Mixing_Length(:,:,:) = 0.0 if (CS%debug .or. (CS%id_Velocity_Scale>0)) diag_Velocity_Scale(:,:,:) = 0.0 - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) BBL_Mix_Length(:,:,:) = 0.0 if (CS%debug .or. (CS%id_BBL_Vel_Scale>0)) BBL_Vel_Scale(:,:,:) = 0.0 if (CS%id_Kd_BBL > 0) Kd_BBL_3d(:,:,:) = 0.0 @@ -545,7 +551,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (CS%id_opt_diff_hML_depth > 0) diff_hML_depth(:,:) = 0.0 endif - !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,I_dt, & + !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,I_dt,BBL_mixing, & !!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int) do j=js,je ! Copy the thicknesses and other fields to 2-d arrays. @@ -653,11 +659,17 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, endif ! Add the diffusivity due to bottom boundary layer mixing, if there is energy to drive this mixing. - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j) BBLD_in = BBLD_io BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) + + ! Add in tidal dissipation energy at the bottom, noting that fluxes%BBL_tidal_dis is + ! in [R Z L2 T-3 ~> W m-2], unlike visc%BBL_meanKE_loss. + if ((CS%ePBL_tidal_effic > 0.0) .and. associated(fluxes%BBL_tidal_dis)) & + BBL_TKE = BBL_TKE + CS%ePBL_tidal_effic * dt * fluxes%BBL_tidal_dis(i,j) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) @@ -694,7 +706,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (CS%debug .or. (CS%id_Velocity_Scale > 0)) then ; do K=1,nz+1 diag_Velocity_Scale(i,j,K) = mixvel(K) enddo ; endif - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) then ; do k=1,nz BBL_Mix_Length(i,j,k) = mixlen_BBL(k) enddo ; endif @@ -711,7 +723,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (report_avg_its) then CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(eCD%OBL_its)) CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then CS%sum_its_BBL(1) = CS%sum_its_BBL(1) + real_to_EFP(real(eCD%BBL_its)) CS%sum_its_BBL(2) = CS%sum_its_BBL(2) + real_to_EFP(1.0) endif @@ -733,6 +745,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, else BLD_1 = BBLD_in ; BLD_2 = BBLD_in BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) + if ((CS%ePBL_tidal_effic > 0.0) .and. associated(fluxes%BBL_tidal_dis)) & + BBL_TKE = BBL_TKE + CS%ePBL_tidal_effic * dt * fluxes%BBL_tidal_dis(i,j) u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & u_star_BBL, Kd_1, BLD_1, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp1, eCD_tmp) @@ -761,7 +775,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, enddo ! j-loop - if (CS%debug .and. (CS%ePBL_BBL_effic > 0.0)) then + if (CS%debug .and. BBL_mixing) then call hchksum(visc%BBL_meanKE_loss, "ePBL visc%BBL_meanKE_loss", G%HI, & unscale=GV%H_to_MKS*US%L_T_to_m_s**2*US%s_to_T) call hchksum(visc%ustar_BBL, "ePBL visc%ustar_BBL", G%HI, unscale=GV%H_to_MKS*US%s_to_T) @@ -786,7 +800,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, diag_Mixing_Length, CS%diag) if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, diag_Velocity_Scale, CS%diag) if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, diag_mStar_MIX, CS%diag) - if (CS%ePBL_BBL_effic > 0.0) then + if (BBL_mixing) then if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, Kd_BBL_3d, CS%diag) if (CS%id_BBL_Mix_Length > 0) call post_data(CS%id_BBL_Mix_Length, BBL_Mix_Length, CS%diag) if (CS%id_BBL_Vel_Scale > 0) call post_data(CS%id_BBL_Vel_Scale, BBL_Vel_Scale, CS%diag) @@ -2095,7 +2109,7 @@ subroutine ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, absf, & no_MKE_conversion = ((CS%direct_calc) ) ! .and. (CS%MKE_to_TKE_effic == 0.0)) ! Add bottom boundary layer mixing if there is energy to support it. - if ((CS%ePBL_BBL_effic <= 0.0) .or. (BBL_TKE_in <= 0.0)) then + if (((CS%ePBL_BBL_effic <= 0.0) .and. (CS%ePBL_tidal_effic <= 0.0)) .or. (BBL_TKE_in <= 0.0)) then ! There is no added bottom boundary layer mixing. BBLD_io = 0.0 Kd_BBL(:) = 0.0 @@ -3393,7 +3407,8 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) integer :: mstar_mode, LT_enhance, wT_mode integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: use_omega - logical :: no_BBL ! If true, EPBL_BBL_EFFIC < 0 and bottom boundary layer mixing is not enabled. + logical :: no_BBL ! If true, EPBL_BBL_EFFIC < 0 and EPBL_BBL_TIDAL_EFFIC < 0, so + ! bottom boundary layer mixing is not enabled. logical :: use_la_windsea isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -3674,7 +3689,12 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "The efficiency of bottom boundary layer mixing via ePBL. Setting this to a "//& "value that is greater than 0 to enable bottom boundary layer mixing from EPBL.", & units="nondim", default=0.0, scale=US%L_to_Z**2) - no_BBL = (CS%ePBL_BBL_effic<=0.0) + call get_param(param_file, mdl, "EPBL_BBL_TIDAL_EFFIC", CS%ePBL_tidal_effic, & + "The efficiency of bottom boundary layer mixing via ePBL driven by the "//& + "bottom drag dissipation of tides, as provided in fluxes%BBL_tidal_dis.", & + units="nondim", default=0.0, scale=US%L_to_Z**2) !### Change the default to follow EPBL_BBL_EFFIC? + no_BBL = ((CS%ePBL_BBL_effic <= 0.0) .and. (CS%ePBL_tidal_effic <= 0.0)) + call get_param(param_file, mdl, "USE_BBLD_ITERATION", CS%Use_BBLD_iteration, & "A logical that specifies whether or not to use the distance to the top of the "//& "actively turbulent bottom boundary layer to help set the EPBL length scale.", & @@ -3885,7 +3905,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) Time, 'Velocity Scale that is used.', units='m s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, & Time, 'Total mstar that is used.', 'nondim') - if (CS%ePBL_BBL_effic > 0.0) then + if ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0)) then CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_ePBL_BBL', diag%axesTi, & Time, 'ePBL bottom boundary layer diffusivity', units='m2 s-1', conversion=GV%HZ_T_to_m2_s) CS%id_BBL_Mix_Length = register_diag_field('ocean_model', 'BBL_Mixing_Length', diag%axesTi, & @@ -3936,7 +3956,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) CS%TKE_diagnostics = (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & CS%id_TKE_conv_decay) > 0) - if (CS%ePBL_BBL_effic > 0.0) then + if ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0)) then CS%TKE_diagnostics = CS%TKE_diagnostics .or. & (max(CS%id_TKE_BBL, CS%id_TKE_BBL_mixing, CS%id_TKE_BBL_decay) > 0) endif @@ -3962,7 +3982,7 @@ subroutine energetic_PBL_end(CS) write (mesg,*) "Average ePBL iterations = ", avg_its call MOM_mesg(mesg) - if (CS%ePBL_BBL_effic > 0.0) then + if ((CS%ePBL_BBL_effic > 0.0) .or. (CS%ePBL_tidal_effic > 0.0)) then call EFP_sum_across_PEs(CS%sum_its_BBL, 2) avg_its = EFP_to_real(CS%sum_its_BBL(1)) / EFP_to_real(CS%sum_its_BBL(2)) write (mesg,*) "Average ePBL BBL iterations = ", avg_its From dbfdfa81c3073d4e40c964bc0463616314c1a7ee Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 20 Jan 2025 09:57:36 -0500 Subject: [PATCH 4/4] +(*)EPBL_BBL_EFFIC_BUG & DRAG_DIFFUSIVITY_ANSWER_DATE Previously, visc%BBL_meanKE_loss (which was called visc%TKE_BBL before it was renamed) was missing a factor of the square root of the drag coefficient compared with the net loss of kinetic energy. This was compensated for by multiplication by factors of the square root of the drag coefficient in add_drag_diffusivity and add_LOTW_BBL_diffusivity, but when an equivalent expression was added in the ePBL BBL code this was erroneously omitted. Moreover, Alistair has had a comment questioning this added factor in add_LOTW_BBL_diffusivity for a decade without adequate resolution. To correct this confusing situation, visc%BBL_meanKE_loss has been changed to include the missing factor of the square root of the drag coefficient, while the new variable visc%BBL_meanKE_loss_sqrtCd was added to allow for the previous answers to be recovered when the new runtime parameter EPBL_BBL_EFFIC_BUG is set to true and DRAG_DIFFUSIVITY_ANSWER_DATE is set below 20250302. Because the ePBL bottom boundary layer was only added to dev/gfdl a month ago and has not yet been merged into main, we can be confident that it has only received very limited use as yet, so the default for EPBL_BBL_EFFIC_BUG is false but this default will change answers when EPBL_BBL_EFFIC > 0. The default for DRAG_DIFFUSIVITY_ANSWER_DATE is 20250101, which will preserve the previous answers, but the default should later be taken from DEFAULT_ANSWER_DATE. By default, answers are unchanged in any cases that are more than a month old, but answers can change by default in a few very recent experiments. There are two new runtime parameters in some MOM_parameter_doc files. --- src/core/MOM_variables.F90 | 3 ++ .../vertical/MOM_energetic_PBL.F90 | 13 ++++++- .../vertical/MOM_set_diffusivity.F90 | 38 ++++++++++++++++--- .../vertical/MOM_set_viscosity.F90 | 2 + 4 files changed, 49 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index efd558cdb3..cabda1cf58 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -239,6 +239,9 @@ module MOM_variables !! h points [H T-1 ~> m s-1 or kg m-2 s-1]. BBL_meanKE_loss, & !< The viscous loss of mean kinetic energy in the bottom boundary layer !! [H L2 T-3 ~> m3 s-3 or W m-2]. + BBL_meanKE_loss_sqrtCd, & !< The viscous loss of mean kinetic energy in the bottom boundary layer + !! divided by the square root of the drag coefficient [H L2 T-3 ~> m3 s-3 or W m-2]. + !! This is being set only to retain old answers, and should be phased out. taux_shelf, & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa]. tauy_shelf !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa]. real, allocatable, dimension(:,:) :: tbl_thick_shelf_u diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 6e3cdfa1a5..6caaace9e9 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -198,6 +198,9 @@ module MOM_energetic_PBL !! energetics that accounts for an exponential decay of TKE from a !! near-bottom source and an assumed piecewise linear linear profile !! of the buoyancy flux response to a change in a diffusivity. + logical :: BBL_effic_bug !< If true, overestimate the efficiency of the non-tidal ePBL bottom boundary + !! layer diffusivity by a factor of 1/sqrt(CDRAG), which is often a factor of + !! about 18.3. !/ Options for documenting differences from parameter choices integer :: options_diff !< If positive, this is a coded integer indicating a pair of @@ -662,7 +665,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, if (BBL_mixing) then if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j) BBLD_in = BBLD_io - BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) + if (CS%BBL_effic_bug) then + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss_sqrtCd(i,j) + else + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%BBL_meanKE_loss(i,j) + endif u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) ! Add in tidal dissipation energy at the bottom, noting that fluxes%BBL_tidal_dis is @@ -3752,6 +3759,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "length scale by rotation in the bottom boundary layer. Making this larger "//& "decreases the bottom boundary layer diffusivity.", & units="nondim", default=CS%Ekman_scale_coef, do_not_log=no_BBL) + call get_param(param_file, mdl, "EPBL_BBL_EFFIC_BUG", CS%BBL_effic_bug, & + "If true, overestimate the efficiency of the non-tidal ePBL bottom boundary "//& + "layer diffusivity by a factor of 1/sqrt(CDRAG), which is often a factor of "//& + "about 18.3.", default=.false., do_not_log=(CS%ePBL_BBL_effic<=0.0)) call get_param(param_file, mdl, "DECAY_ADJUSTED_BBL_TKE", CS%decay_adjusted_BBL_TKE, & "If true, include an adjustment factor in the bottom boundary layer energetics "//& diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 3e5b2d44c8..1b313aaca3 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -163,11 +163,17 @@ module MOM_set_diffusivity !! calculations. Values below 20190101 recover the answers from the !! end of 2018, while higher values use updated and more robust forms !! of the same expressions. Values above 20240630 use more accurate - !! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true. + !! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true. Values + !! above 20250301 use less confusing expressions to set the bottom-drag + !! generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false. integer :: LOTW_BBL_answer_date !< The vintage of the order of arithmetic and expressions !! in the LOTW_BBL calculations. Values below 20240630 recover the !! original answers, while higher values use more accurate expressions. !! This only applies when USE_LOTW_BBL_DIFFUSIVITY is true. + integer :: drag_diff_answer_date !< The vintage of the order of arithmetic in the drag diffusivity + !! calculations. Values above 20250301 use less confusing expressions + !! to set the bottom-drag generated diffusivity when + !! USE_LOTW_BBL_DIFFUSIVITY is false. character(len=200) :: inputdir !< The directory in which input files are found type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() !< Control structure for a child module @@ -1414,7 +1420,11 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, ! If ustar_h = 0, this is land so this value doesn't matter. I2decay(i) = 0.5*CS%IMax_decay endif - TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j) + if (CS%drag_diff_answer_date <= 20250301) then + TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss_sqrtCd(i,j) + else + TKE(i) = (CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz)) ) * visc%BBL_meanKE_loss(i,j) + endif if (associated(fluxes%BBL_tidal_dis)) & TKE(i) = TKE(i) + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H * & @@ -1644,8 +1654,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bo ! Energy input at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2]. ! (Note that visc%BBL_meanKE_loss is in [H L2 T-3 ~> m3 s-3 or W m-2], set in set_BBL_TKE().) - !### I am still unsure about sqrt(cdrag) in this expressions - AJA - BBL_meanKE_dis = cdrag_sqrt * visc%BBL_meanKE_loss(i,j) + BBL_meanKE_dis = visc%BBL_meanKE_loss(i,j) ! Add in tidal dissipation energy at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2]. ! Note that BBL_tidal_dis is in [R Z L2 T-3 ~> W m-2]. if (associated(fluxes%BBL_tidal_dis)) & @@ -1956,6 +1965,9 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (allocated(visc%BBL_meanKE_loss)) then do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss(i,j) = 0.0 ; enddo ; enddo endif + if (allocated(visc%BBL_meanKE_loss_sqrtCd)) then + do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss_sqrtCd(i,j) = 0.0 ; enddo ; enddo + endif return endif @@ -2079,7 +2091,13 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) (G%areaCu(I,j)*(ustar(I)*ustar(I)))) + & ((G%areaCv(i,J-1)*(vstar(i,J-1)*vstar(i,J-1))) + & (G%areaCv(i,J)*(vstar(i,J)*vstar(i,J)))) ) ) - visc%BBL_meanKE_loss(i,j) = & + visc%BBL_meanKE_loss(i,j) = cdrag_sqrt * & + ((((G%areaCu(I-1,j)*(ustar(I-1)*u2_bbl(I-1))) + & + (G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + & + ((G%areaCv(i,J-1)*(vstar(i,J-1)*v2_bbl(i,J-1))) + & + (G%areaCv(i,J) * (vstar(i,J)*v2_bbl(i,J)))) )*G%IareaT(i,j)) + ! The following line could be omitted if SET_DIFF_ANSWER_DATE > 20250301 and EPBL_BBL_EFFIC_BUG is false. + visc%BBL_meanKE_loss_sqrtCd(i,j) = & ((((G%areaCu(I-1,j)*(ustar(I-1)*u2_bbl(I-1))) + & (G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + & ((G%areaCv(i,J-1)*(vstar(i,J-1)*v2_bbl(i,J-1))) + & @@ -2289,7 +2307,9 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ call get_param(param_file, mdl, "SET_DIFF_ANSWER_DATE", CS%answer_date, & "The vintage of the order of arithmetic and expressions in the set diffusivity "//& "calculations. Values below 20190101 recover the answers from the end of 2018, "//& - "while higher values use updated and more robust forms of the same expressions.", & + "while higher values use updated and more robust forms of the same expressions. "//& + "Values above 20250301 also use less confusing expressions to set the bottom-drag "//& + "generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) @@ -2405,6 +2425,12 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "USE_LOTW_BBL_DIFFUSIVITY is true.", & default=20190101, do_not_log=.not.CS%use_LOTW_BBL_diffusivity) !### Set default as default=default_answer_date, or use SET_DIFF_ANSWER_DATE. + call get_param(param_file, mdl, "DRAG_DIFFUSIVITY_ANSWER_DATE", CS%drag_diff_answer_date, & + "The vintage of the order of arithmetic in the drag diffusivity calculations. "//& + "Values above 20250301 use less confusing expressions to set the bottom-drag "//& + "generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false. ", & + default=20250101, do_not_log=CS%use_LOTW_BBL_diffusivity.or.(CS%BBL_effic<=0.0)) + !### Set default as default=default_answer_date, or use SET_DIFF_ANSWER_DATE. CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, Time, & 'Bottom Boundary Layer Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index bebdc4e7d0..e701837ae1 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -3140,6 +3140,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0) allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0) allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0) + allocate(visc%BBL_meanKE_loss_sqrtCd(isd:ied,jsd:jed), source=0.0) CS%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', & diag%axesCu1, Time, 'BBL thickness at u points', 'm', conversion=US%Z_to_m) @@ -3215,6 +3216,7 @@ subroutine set_visc_end(visc, CS) if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu) if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl) if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss) + if (allocated(visc%BBL_meanKE_loss_sqrtCd)) deallocate(visc%BBL_meanKE_loss_sqrtCd) if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf) if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf) if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)