diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index f668f24508..3748684fd4 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -208,7 +208,8 @@ module MOM_variables ustar_BBL => NULL() !< The turbulence velocity in the bottom boundary layer at h points [Z s-1 ~> m s-1]. real, pointer, dimension(:,:) :: TKE_BBL => NULL() !< A term related to the bottom boundary layer source of turbulent kinetic - !! energy, currently in [m3 s-3], but will later be changed to [W m-2]. + !! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed + !! to [kg Z3 m-3 T-3 ~> W m-2]. real, pointer, dimension(:,:) :: & taux_shelf => NULL(), & !< The zonal stresses on the ocean under shelves [Pa]. tauy_shelf => NULL() !< The meridional stresses on the ocean under shelves [Pa]. diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 93f1c33b06..c5fe83a9e7 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -173,14 +173,14 @@ module MOM_set_diffusivity N2_3d => NULL(), & !< squared buoyancy frequency at interfaces [T-2 ~> s-2] Kd_user => NULL(), & !< user-added diffusivity at interfaces [Z2 T-1 ~> m2 s-1] Kd_BBL => NULL(), & !< BBL diffusivity at interfaces [Z2 T-1 ~> m2 s-1] - Kd_work => NULL(), & !< layer integrated work by diapycnal mixing [kg T-3 ~> W m-2] - maxTKE => NULL(), & !< energy required to entrain to h_max [m3 T-3 ~> m3 s-3] + Kd_work => NULL(), & !< layer integrated work by diapycnal mixing [kg Z3 m-3 T-3 ~> W m-2] + maxTKE => NULL(), & !< energy required to entrain to h_max [Z3 T-3 ~> m3 s-3] KT_extra => NULL(), & !< double diffusion diffusivity for temp [Z2 T-1 ~> m2 s-1]. KS_extra => NULL() !< double diffusion diffusivity for saln [Z2 T-1 ~> m2 s-1]. real, pointer, dimension(:,:,:) :: TKE_to_Kd => NULL() !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between TKE !! dissipated within a layer and Kd in that layer - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> T2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] end type diffusivity_diags @@ -247,7 +247,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & maxTKE, & !< energy required to entrain to h_max [m3 T-3] TKE_to_Kd !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between !< TKE dissipated within a layer and Kd in that layer - !< [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !< [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(G)+1) :: & N2_int, & !< squared buoyancy frequency associated at interfaces [T-2 ~> s-2] @@ -395,7 +395,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call sfc_bkgnd_mixing(G, US, CS%bkgnd_mixing_csp) !$OMP parallel do default(shared) private(dRho_int, N2_lay, N2_int, N2_bot, KT_extra, & - !$OMP KS_extra, TKE_to_Kd,maxTKE, dissip, kb) + !$OMP KS_extra, TKE_to_Kd, maxTKE, dissip, kb) do j=js,je ! Set up variables related to the stratification. @@ -523,7 +523,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_work)) then do k=1,nz ; do i=is,ie - dd%Kd_Work(i,j,k) = GV%Rho0 * (US%Z_to_m**3 * Kd_lay(i,j,k)) * N2_lay(i,k) * & + dd%Kd_Work(i,j,k) = GV%Rho0 * Kd_lay(i,j,k) * N2_lay(i,k) * & GV%H_to_Z*h(i,j,k) ! Watt m-2 s = kg s-3 enddo ; enddo endif @@ -680,9 +680,9 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & !! TKE dissipated within a layer and the !! diapycnal diffusivity witin that layer, !! usually (~Rho_0 / (G_Earth * dRho_lay)) - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(G)), intent(out) :: maxTKE !< The energy required to for a layer to entrain - !! to its maximum realizable thickness [m3 T-3] + !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3] integer, dimension(SZI_(G)), intent(out) :: kb !< Index of lightest layer denser than the buffer !! layer, or -1 without a bulk mixed layer. ! Local variables @@ -736,13 +736,13 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & ! Simple but coordinate-independent estimate of Kd/TKE if (CS%simple_TKE_to_Kd) then do k=1,nz ; do i=is,ie - hN2pO2 = US%Z_to_m**3 * (GV%H_to_Z * h(i,j,k)) * (N2_lay(i,k) + Omega2) ! Units of m3 Z-2 T-2. + hN2pO2 = (GV%H_to_Z * h(i,j,k)) * (N2_lay(i,k) + Omega2) ! Units of Z T-2. if (hN2pO2>0.) then - TKE_to_Kd(i,k) = 1.0 / hN2pO2 ! Units of Z2 s2 m-3. + TKE_to_Kd(i,k) = 1.0 / hN2pO2 ! Units of T2 Z-1. else; TKE_to_Kd(i,k) = 0.; endif ! The maximum TKE conversion we allow is really a statement ! about the upper diffusivity we allow. Kd_max must be set. - maxTKE(i,k) = hN2pO2 * CS%Kd_max ! Units of m3 T-3. + maxTKE(i,k) = hN2pO2 * CS%Kd_max ! Units of Z3 T-3. enddo ; enddo kb(is:ie) = -1 ! kb should not be used by any code in non-layered mode -AJA return @@ -837,7 +837,7 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & enddo do k=2,kmb ; do i=is,ie maxTKE(i,k) = 0.0 - TKE_to_Kd(i,k) = US%m_to_Z**3 / ((N2_lay(i,k) + Omega2) * & + TKE_to_Kd(i,k) = 1.0 / ((N2_lay(i,k) + Omega2) * & (GV%H_to_Z*(h(i,j,k) + H_neglect))) enddo ; enddo do k=kmb+1,kb_min-1 ; do i=is,ie @@ -858,10 +858,10 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & ! dRho_int should already be non-negative, so the max is redundant? dh_max = maxEnt(i,k) * (1.0 + dsp1_ds(i,k)) dRho_lay = 0.5 * max(dRho_int(i,K) + dRho_int(i,K+1), 0.0) - maxTKE(i,k) = US%Z_to_m**3 * I_dt * (G_IRho0 * & + maxTKE(i,k) = I_dt * (G_IRho0 * & (0.5*max(dRho_int(i,K+1) + dsp1_ds(i,k)*dRho_int(i,K), 0.0))) * & - ((GV%H_to_Z*h(i,j,k) + dh_max) * maxEnt(i,k)) - TKE_to_Kd(i,k) = US%m_to_Z**3 / (G_Rho0 * dRho_lay + & + ((GV%H_to_Z*h(i,j,k) + dh_max) * maxEnt(i,k)) + TKE_to_Kd(i,k) = 1.0 / (G_Rho0 * dRho_lay + & CS%omega**2 * GV%H_to_Z*(h(i,j,k) + H_neglect)) endif enddo ; enddo @@ -1150,7 +1150,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & !! TKE dissipated within a layer and the !! diapycnal diffusivity witin that layer, !! usually (~Rho_0 / (G_Earth * dRho_lay)) - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(G)), intent(in) :: maxTKE !< The energy required to for a layer to entrain !! to its maximum-realizable thickness [m3 T-3 ~> m3 s-3] integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer @@ -1176,12 +1176,12 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & ! the local ustar, times R0_g [kg m-2] Rho_top, & ! density at top of the BBL [kg m-3] TKE, & ! turbulent kinetic energy available to drive - ! bottom-boundary layer mixing in a layer [m3 T-3 ~> m3 s-3] + ! bottom-boundary layer mixing in a layer [Z3 T-3 ~> m3 s-3] I2decay ! inverse of twice the TKE decay scale [Z-1 ~> m-1]. - real :: TKE_to_layer ! TKE used to drive mixing in a layer [m3 T-3 ~> m3 s-3] - real :: TKE_Ray ! TKE from layer Rayleigh drag used to drive mixing in layer [m3 T-3 ~> m3 s-3] - real :: TKE_here ! TKE that goes into mixing in this layer [m3 T-3 ~> m3 s-3] + real :: TKE_to_layer ! TKE used to drive mixing in a layer [Z3 T-3 ~> m3 s-3] + real :: TKE_Ray ! TKE from layer Rayleigh drag used to drive mixing in layer [Z3 T-3 ~> m3 s-3] + real :: TKE_here ! TKE that goes into mixing in this layer [Z3 T-3 ~> m3 s-3] real :: dRl, dRbot ! temporaries holding density differences [kg m-3] real :: cdrag_sqrt ! square root of the drag coefficient [nondim] real :: ustar_h ! value of ustar at a thickness point [Z T-1 ~> m s-1]. @@ -1230,12 +1230,11 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & ! 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)*(GV%H_to_Z*h(i,j,nz))) ) * & - (US%T_to_s**3 * visc%TKE_BBL(i,j)) + TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*(GV%H_to_Z*h(i,j,nz))) ) * & + visc%TKE_BBL(i,j) if (associated(fluxes%TKE_tidal)) & - TKE(i) = TKE(i) + (US%T_to_s**3 * fluxes%TKE_tidal(i,j)) * I_Rho0 * & + TKE(i) = TKE(i) + (US%T_to_s**3 * US%m_to_Z**3 * fluxes%TKE_tidal(i,j)) * I_Rho0 * & (CS%BBL_effic * exp(-I2decay(i)*(GV%H_to_Z*h(i,j,nz)))) ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following @@ -1286,14 +1285,14 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & TKE_to_layer = TKE(i) else dRl = Rint(K+1) - Rint(K) ; dRbot = Rint(K+1) - Rho_top(i) - TKE_to_layer = TKE(i) * dRl * (3.0*dRbot*(Rint(K) - Rho_top(i)) +& - dRl**2) / dRbot**3 + TKE_to_layer = TKE(i) * dRl * & + (3.0*dRbot*(Rint(K) - Rho_top(i)) + dRl**2) / dRbot**3 endif 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%Z_to_m * G%IareaT(i,j) * & - US%T_to_s**3 * & + if (Rayleigh_drag) TKE_Ray = 0.5*CS%BBL_effic * G%IareaT(i,j) * & + US%m_to_Z**2 * US%T_to_s**3 * & ((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 + & @@ -1327,13 +1326,12 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & TKE(i) = TKE(i) + TKE_Ray elseif (Kd_lay(i,j,k) + (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k) > & maxTKE(i,k) * TKE_to_Kd(i,k)) then - TKE_here = ((TKE_to_layer + TKE_Ray) + Kd_lay(i,j,k) / TKE_to_Kd(i,k)) & - - maxTKE(i,k) + TKE_here = ((TKE_to_layer + TKE_Ray) + Kd_lay(i,j,k) / TKE_to_Kd(i,k)) - maxTKE(i,k) ! ### Non-bracketed ternary sum TKE(i) = TKE(i) - TKE_here + TKE_Ray else TKE_here = TKE_to_layer + TKE_Ray - TKE(i) = TKE(i) - TKE_to_Layer + TKE(i) = TKE(i) - TKE_to_layer endif if (TKE(i) < 0.0) TKE(i) = 0.0 ! This should be unnecessary? @@ -1395,11 +1393,10 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1] ! Local variables - real :: TKE_column ! net TKE input into the column [m3 T-3 ~> m3 s-3] - real :: TKE_to_layer ! TKE used to drive mixing in a layer [m3 T-3 ~> m3 s-3] - real :: TKE_remaining ! remaining TKE available for mixing in this layer and above [m3 T-3 ~> m3 s-3] - real :: TKE_consumed ! TKE used for mixing in this layer [m3 T-3 ~> m3 s-3] - real :: TKE_Kd_wall ! TKE associated with unlimited law of the wall mixing [m3 T-3 ~> m3 s-3] + real :: TKE_column ! net TKE input into the column [Z3 T-3 ~> m3 s-3] + real :: TKE_remaining ! remaining TKE available for mixing in this layer and above [Z3 T-3 ~> m3 s-3] + real :: TKE_consumed ! TKE used for mixing in this layer [Z3 T-3 ~> m3 s-3] + real :: TKE_Kd_wall ! TKE associated with unlimited law of the wall mixing [Z3 T-3 ~> m3 s-3] real :: cdrag_sqrt ! square root of the drag coefficient [nondim] real :: ustar ! value of ustar at a thickness point [Z T-1 ~> m s-1]. real :: ustar2 ! square of ustar, for convenience [Z2 T-2 ~> m2 s-2] @@ -1450,13 +1447,14 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & Idecay = CS%IMax_decay if ((ustar > 0.0) .and. (absf > CS%IMax_decay * ustar)) Idecay = absf / ustar - ! Energy input at the bottom [m3 s-3]. - ! (Note that visc%TKE_BBL is in m3 s-3, set in set_BBL_TKE().) + ! Energy input at the bottom [Z3 T-3 ~> m3 s-3]. + ! (Note that visc%TKE_BBL is in [Z3 T-3 ~> m3 s-3], set in set_BBL_TKE().) ! I am still unsure about sqrt(cdrag) in this expressions - AJA - TKE_column = cdrag_sqrt * (US%T_to_s**3 * visc%TKE_BBL(i,j)) + TKE_column = cdrag_sqrt * visc%TKE_BBL(i,j) ! Add in tidal dissipation energy at the bottom [m3 s-3]. ! Note that TKE_tidal is in [W m-2]. - if (associated(fluxes%TKE_tidal)) TKE_column = TKE_column + (US%T_to_s**3 * fluxes%TKE_tidal(i,j)) * I_Rho0 + if (associated(fluxes%TKE_tidal)) & + TKE_column = TKE_column + US%m_to_Z**3*US%T_to_s**3 * fluxes%TKE_tidal(i,j) * I_Rho0 TKE_column = CS%BBL_effic * TKE_column ! Only use a fraction of the mechanical dissipation for mixing. TKE_remaining = TKE_column @@ -1474,8 +1472,8 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & ! Add in additional energy input from bottom-drag against slopes (sides) if (Rayleigh_drag) TKE_remaining = TKE_remaining + & - US%T_to_s**3 * & - 0.5*CS%BBL_effic * US%Z_to_m * G%IareaT(i,j) * & + US%m_to_Z**2 * US%T_to_s**3 * & + 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 + & @@ -1499,7 +1497,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & ! TKE associated with Kd_wall [m3 s-2]. ! This calculation if for the volume spanning the interface. - TKE_Kd_wall = US%Z_to_m**3 * Kd_wall * 0.5 * (dh + dhm1) * max(N2_int(i,k), N2_min) + TKE_Kd_wall = Kd_wall * 0.5 * (dh + dhm1) * max(N2_int(i,k), N2_min) ! Now bound Kd such that the associated TKE is no greater than available TKE for mixing. if (TKE_Kd_wall > 0.) then @@ -1544,7 +1542,7 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, !! TKE dissipated within a layer and the !! diapycnal diffusivity witin that layer, !! usually (~Rho_0 / (G_Earth * dRho_lay)) - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), & optional, intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces !! [Z2 T-1 ~> m2 s-1]. @@ -1594,8 +1592,7 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, ustar_sq = max(US%T_to_s * fluxes%ustar(i,j), CS%ustar_min)**2 - TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) & - * (US%Z_to_m**3 * ustar_sq * (US%T_to_s * fluxes%ustar(i,j))) + TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * (US%T_to_s * fluxes%ustar(i,j))) I_decay_len2_TKE = CS%TKE_decay**2 * (f_sq / ustar_sq) if (CS%ML_rad_TKE_decay) & @@ -1611,8 +1608,7 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, if (z1 > 1e-5) then Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) * (1.0 - exp(-z1)) else - Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) & - * (z1 * (1.0 - z1 * (0.5 - C1_6 * z1))) + Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) * (z1 * (1.0 - z1 * (0.5 - C1_6 * z1))) endif Kd_mlr_ml(i) = min(Kd_mlr, CS%ML_rad_kd_max) TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-z1) @@ -1635,11 +1631,13 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, do i=is,ie ; if (do_i(i)) then dzL = GV%H_to_Z*h(i,j,k) ; z1 = dzL*I_decay(i) if (z1 > 1e-5) then - Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) & - * US%m_to_Z * ((1.0 - exp(-z1)) / dzL) + !### I think that this might be dimensionally inconsistent, but untested. -RWH + Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of Z2 T-1 ? + US%m_to_Z * ((1.0 - exp(-z1)) / dzL) ! Units of m-1 ? else - Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) & - * US%m_to_Z * (I_decay(i) * (1.0 - z1 * (0.5 - C1_6*z1))) + !### I think that this might be dimensionally inconsistent, but untested. -RWH + Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of Z2 T-1 ? + US%m_to_Z * (I_decay(i) * (1.0 - z1 * (0.5 - C1_6*z1))) ! Units of m-1 ? endif Kd_mlr = min(Kd_mlr, CS%ML_rad_kd_max) Kd_lay(i,j,k) = Kd_lay(i,j,k) + Kd_mlr @@ -1649,8 +1647,7 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, endif TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-z1) - if (TKE_ml_flux(i) * I_decay(i) & - < 0.1 * CS%Kd_min * US%Z_to_m**3 * Omega2) then + if (TKE_ml_flux(i) * I_decay(i) < 0.1 * CS%Kd_min * Omega2) then do_i(i) = .false. else ; do_any = .true. ; endif endif ; enddo @@ -1787,7 +1784,7 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS) 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%Z_to_m * & + visc%TKE_BBL(i,j) = US%T_to_s**3 * US%m_to_Z**2 * & (((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)) + & @@ -2169,12 +2166,12 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, diag_to_Z CS%tm_csp%Lowmode_itidal_dissipation) then CS%id_Kd_Work = register_diag_field('ocean_model', 'Kd_Work', diag%axesTL, Time, & - 'Work done by Diapycnal Mixing', 'W m-2', conversion=US%s_to_T**3) + 'Work done by Diapycnal Mixing', 'W m-2', conversion=US%Z_to_m**3*US%s_to_T**3) CS%id_maxTKE = register_diag_field('ocean_model', 'maxTKE', diag%axesTL, Time, & - 'Maximum layer TKE', 'm3 s-3', conversion=US%s_to_T**3) + 'Maximum layer TKE', 'm3 s-3', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_TKE_to_Kd = register_diag_field('ocean_model', 'TKE_to_Kd', diag%axesTL, Time, & 'Convert TKE to Kd', 's2 m', & - conversion=US%Z2_T_to_m2_s*(US%T_to_s**3)) + conversion=US%Z2_T_to_m2_s*(US%m_to_Z**3*US%T_to_s**3)) CS%id_N2 = register_diag_field('ocean_model', 'N2', diag%axesTi, Time, & 'Buoyancy frequency squared', 's-2', cmor_field_name='obvfsq', & cmor_long_name='Square of seawater buoyancy frequency', & @@ -2190,8 +2187,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, diag_to_Z "Buoyancy frequency, interpolated to z", z_grid='z') CS%id_N2_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time, conversion=US%s_to_T**2) if (CS%user_change_diff) & - CS%id_Kd_user_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time, & - conversion=US%Z2_T_to_m2_s) + CS%id_Kd_user_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time, conversion=US%Z2_T_to_m2_s) endif endif diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 2af17d734b..b82313dc6c 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -45,11 +45,11 @@ module MOM_tidal_mixing type, public :: tidal_mixing_diags ; private real, pointer, dimension(:,:,:) :: & Kd_itidal => NULL(),& !< internal tide diffusivity at interfaces [Z2 s-1 ~> m2 s-1]. - Fl_itidal => NULL(),& !< vertical flux of tidal turbulent dissipation [m3 s-3] + Fl_itidal => NULL(),& !< vertical flux of tidal turbulent dissipation [Z3 T-3 ~> m3 s-3] Kd_Niku => NULL(),& !< lee-wave diffusivity at interfaces [Z2 s-1 ~> m2 s-1]. - Kd_Niku_work => NULL(),& !< layer integrated work by lee-wave driven mixing [W m-2] - Kd_Itidal_Work => NULL(),& !< layer integrated work by int tide driven mixing [W m-2] - Kd_Lowmode_Work => NULL(),& !< layer integrated work by low mode driven mixing [W m-2] + Kd_Niku_work => NULL(),& !< layer integrated work by lee-wave driven mixing [kg Z3 m-3 T-3 ~> W m-2] + Kd_Itidal_Work => NULL(),& !< layer integrated work by int tide driven mixing [kg Z3 m-3 T-3 ~> W m-2] + Kd_Lowmode_Work => NULL(),& !< layer integrated work by low mode driven mixing [kg Z3 m-3 T-3 ~> W m-2] N2_int => NULL(),& !< Bouyancy frequency squared at interfaces [s-2] vert_dep_3d => NULL(),& !< The 3-d mixing energy deposition [W m-3] Schmittner_coeff_3d => NULL() !< The coefficient in the Schmittner et al mixing scheme, in UNITS? @@ -58,9 +58,9 @@ module MOM_tidal_mixing real, pointer, dimension(:,:,:) :: Kd_lowmode => NULL() !< internal tide diffusivity at interfaces !! due to propagating low modes [Z2 s-1 ~> m2 s-1]. real, pointer, dimension(:,:,:) :: Fl_lowmode => NULL() !< vertical flux of tidal turbulent - !! dissipation due to propagating low modes [m3 s-3] + !! dissipation due to propagating low modes [Z3 T-3 ~> m3 s-3] real, pointer, dimension(:,:) :: & - TKE_itidal_used => NULL(),& !< internal tide TKE input at ocean bottom [W m-2] + TKE_itidal_used => NULL(),& !< internal tide TKE input at ocean bottom [kg Z3 m-3 T-3 ~> W m-2] N2_bot => NULL(),& !< bottom squared buoyancy frequency [s-2] N2_meanz => NULL(),& !< vertically averaged buoyancy frequency [s-2] Polzin_decay_scale_scaled => NULL(),& !< vertical scale of decay for tidal dissipation @@ -123,10 +123,10 @@ module MOM_tidal_mixing real :: Polzin_min_decay_scale !< minimum decay scale of the tidal dissipation !! profile in Polzin formulation [Z ~> m]. - real :: TKE_itide_max !< maximum internal tide conversion [W m-2] + real :: TKE_itide_max !< maximum internal tide conversion [kg Z3 m-3 T-3 ~> W m-2] !! available to mix above the BBL - real :: utide !< constant tidal amplitude [m s-1] used if + real :: utide !< constant tidal amplitude [Z T-1 ~> m s-1] if READ_TIDEAMP is false. real :: kappa_itides !< topographic wavenumber and non-dimensional scaling [Z-1 ~> m-1]. real :: kappa_h2_factor !< factor for the product of wavenumber * rms sgs height character(len=200) :: inputdir !< The directory in which to find input files @@ -146,9 +146,10 @@ module MOM_tidal_mixing type(remapping_CS) :: remap_CS !< The control structure for remapping ! Data containers - real, pointer, dimension(:,:) :: TKE_Niku => NULL() !< Lee wave driven Turbulent Kinetic Energy input [W m-2] + real, pointer, dimension(:,:) :: TKE_Niku => NULL() !< Lee wave driven Turbulent Kinetic Energy input + !! [kg Z3 m-3 T-3 ~> W m-2] real, pointer, dimension(:,:) :: TKE_itidal => NULL() !< The internal Turbulent Kinetic Energy input divided - !! by the bottom stratfication [J m-2]. + !! by the bottom stratfication [kg Z3 m-3 T-2 ~> J m-2]. real, pointer, dimension(:,:) :: Nb => NULL() !< The near bottom buoyancy frequency [s-1]. real, pointer, dimension(:,:) :: mask_itidal => NULL() !< A mask of where internal tide energy is input real, pointer, dimension(:,:) :: h2 => NULL() !< Squared bottom depth variance [m2]. @@ -421,7 +422,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ call get_param(param_file, mdl, "UTIDE", CS%utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & - units="m s-1", default=0.0) + units="m s-1", default=0.0, scale=US%m_to_Z*US%T_to_s) call safe_alloc_ptr(CS%tideamp,is,ie,js,je) ; CS%tideamp(:,:) = CS%utide call get_param(param_file, mdl, "KAPPA_H2_FACTOR", CS%kappa_h2_factor, & @@ -430,7 +431,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_max, & "The maximum internal tide energy source available to mix \n"//& "above the bottom boundary layer with INT_TIDE_DISSIPATION.", & - units="W m-2", default=1.0e3) + units="W m-2", default=1.0e3, scale=US%m_to_Z**3*US%T_to_s**3) call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, & "If true, read a file (given by TIDEAMP_FILE) containing \n"//& @@ -445,7 +446,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc") filename = trim(CS%inputdir) // trim(tideamp_file) call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename) - call MOM_read_data(filename, 'tideamp', CS%tideamp, G%domain, timelevel=1) + call MOM_read_data(filename, 'tideamp', CS%tideamp, G%domain, timelevel=1, scale=US%m_to_Z*US%T_to_s) endif call get_param(param_file, mdl, "H2_FILE", h2_file, & @@ -466,8 +467,8 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ CS%h2(i,j) = hamp*hamp utide = CS%tideamp(i,j) - ! Compute the fixed part of internal tidal forcing; units are [J m-2 = kg s-2] here. - CS%TKE_itidal(i,j) = 0.5*US%Z_to_m * CS%kappa_h2_factor*GV%Rho0*& + ! Compute the fixed part of internal tidal forcing; units are [kg Z3 m-3 T-2 ~> J m-2 = kg s-2] here. + CS%TKE_itidal(i,j) = 0.5 * CS%kappa_h2_factor * GV%Rho0 * & CS%kappa_itides * CS%h2(i,j) * utide*utide enddo ; enddo @@ -488,7 +489,8 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", & filename) call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je) ; CS%TKE_Niku(:,:) = 0.0 - call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1 ) ! ??? timelevel -aja + call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1, & ! ??? timelevel -aja + scale=US%m_to_Z**3*US%T_to_s**3) CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & @@ -582,21 +584,25 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ else CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, & - 'Internal Tide Driven Turbulent Kinetic Energy', 'W m-2') + 'Internal Tide Driven Turbulent Kinetic Energy', 'W m-2', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, & 'Bottom Buoyancy Frequency', 's-1') CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, & - 'Internal Tide Driven Diffusivity (from propagating low modes)', 'm2 s-1', conversion=US%Z_to_m**2) + 'Internal Tide Driven Diffusivity (from propagating low modes)', & + 'm2 s-1', conversion=US%Z_to_m**2) CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation', 'm3 s-3') + 'Vertical flux of tidal turbulent dissipation', & + 'm3 s-3', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'm3 s-3') + 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', & + 'm3 s-3', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,Time, & - 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'm', conversion=US%Z_to_m) + 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', & + 'm', conversion=US%Z_to_m) CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model', & 'Polzin_decay_scale_scaled', diag%axesT1, Time, & @@ -610,17 +616,18 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ 'Buoyancy frequency squared averaged over the water column', 's-2') CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing', 'W m-2') + 'Work done by Internal Tide Diapycnal Mixing', 'W m-2', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, & - 'Work done by Nikurashin Lee Wave Drag Scheme', 'W m-2') + 'Work done by Nikurashin Lee Wave Drag Scheme', 'W m-2', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing (low modes)', 'W m-2') + 'Work done by Internal Tide Diapycnal Mixing (low modes)', & + 'W m-2', conversion=(US%Z_to_m**3*US%s_to_T**3)) if (CS%Lee_wave_dissipation) then CS%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,Time, & - 'Lee wave Driven Turbulent Kinetic Energy', 'W m-2') + 'Lee wave Driven Turbulent Kinetic Energy', 'W m-2', conversion=(US%Z_to_m**3*US%s_to_T**3)) CS%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,Time, & 'Lee Wave Driven Diffusivity', 'm2 s-1', conversion=US%Z_to_m**2) endif @@ -666,12 +673,12 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, C !! interfaces [T-2 ~> s-2]. integer, intent(in) :: j !< The j-index to work on real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE - !! TKE dissipated within a layer and the - !! diapycnal diffusivity witin that layer, + !! dissipated within a layer and the + !! diapycnal diffusivity within that layer, !! usually (~Rho_0 / (G_Earth * dRho_lay)) - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(G)), intent(in) :: max_TKE !< The energy required to for a layer to entrain - !! to its maximum realizable thickness [m3 T-3 ~> m3 s-3] + !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3] type(tidal_mixing_cs), pointer :: CS !< The control structure for this module real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(inout) :: Kd_lay !< The diapycnal diffusvity in layers [Z2 T-1 ~> m2 s-1]. @@ -946,12 +953,12 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, !! layers [T-2 ~> s-2]. integer, intent(in) :: j !< The j-index to work on real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE - !! TKE dissipated within a layer and the - !! diapycnal diffusivity witin that layer, + !! dissipated within a layer and the + !! diapycnal diffusivity within that layer, !! usually (~Rho_0 / (G_Earth * dRho_lay)) - !! [Z2 T-1 / m3 T-3 = Z2 T2 m-3 ~> s2 m-1] + !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1] real, dimension(SZI_(G),SZK_(G)), intent(in) :: max_TKE !< The energy required to for a layer to entrain - !! to its maximum realizable thickness [m3 T-3 ~> m3 s-3] + !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3] type(tidal_mixing_cs), pointer :: CS !< The control structure for this module real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(inout) :: Kd_lay !< The diapycnal diffusvity in layers [Z2 T-1 ~> m2 s-1]. @@ -969,9 +976,9 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, htot, & ! total thickness above or below a layer, or the ! integrated thickness in the BBL [Z ~> m]. htot_WKB, & ! WKB scaled distance from top to bottom [Z ~> m]. - TKE_itidal_bot, & ! internal tide TKE at ocean bottom [m3 s-3] - TKE_Niku_bot, & ! lee-wave TKE at ocean bottom [m3 s-3] - TKE_lowmode_bot, & ! internal tide TKE at ocean bottom lost from all remote low modes [m3 s-3] (BDM) + TKE_itidal_bot, & ! internal tide TKE at ocean bottom [Z3 T-3 ~> m3 s-3] + TKE_Niku_bot, & ! lee-wave TKE at ocean bottom [Z3 T-3 ~> m3 s-3] + TKE_lowmode_bot, & ! internal tide TKE at ocean bottom lost from all remote low modes [Z3 T-3 ~> m3 s-3] (BDM) Inv_int, & ! inverse of TKE decay for int tide over the depth of the ocean [nondim] Inv_int_lee, & ! inverse of TKE decay for lee waves over the depth of the ocean [nondim] Inv_int_low, & ! inverse of TKE decay for low modes over the depth of the ocean [nondim] (BDM) @@ -981,9 +988,9 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz) ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz N2_meanz, & ! vertically averaged squared buoyancy frequency [s-2] for WKB scaling - TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) - TKE_Niku_rem, & ! remaining lee-wave TKE - TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) (BDM) + TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) [Z3 T-3 ~> m3 s-3] + TKE_Niku_rem, & ! remaining lee-wave TKE [Z3 T-3 ~> m3 s-3] + TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) [Z3 T-3 ~> m3 s-3] (BDM) TKE_frac_top, & ! fraction of bottom TKE that should appear at top of a layer [nondim] TKE_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer [nondim] TKE_frac_top_lowmode, & @@ -993,14 +1000,14 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, real :: I_rho0 ! 1 / RHO0 [m3 kg-1] real :: Kd_add ! diffusivity to add in a layer [Z2 s-1 ~> m2 s-1]. - real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) [m3 s-3] - real :: TKE_Niku_lay ! lee-wave TKE imparted to a layer [m3 s-3] - real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) [m3 s-3] (BDM) + real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) [Z3 T-3 ~> m3 s-3] + real :: TKE_Niku_lay ! lee-wave TKE imparted to a layer [Z3 T-3 ~> m3 s-3] + real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) [Z3 T-3 ~> m3 s-3] (BDM) real :: frac_used ! fraction of TKE that can be used in a layer [nondim] real :: Izeta ! inverse of TKE decay scale [Z-1 ~> m-1]. real :: Izeta_lee ! inverse of TKE decay scale for lee waves [Z-1 ~> m-1]. real :: z0_psl ! temporary variable [Z ~> m]. - real :: TKE_lowmode_tot ! TKE from all low modes [W m-2] (BDM) + real :: TKE_lowmode_tot ! TKE from all low modes [kg Z3 m-3 T-3 ~> W m-2] (BDM) logical :: use_Polzin, use_Simmons character(len=160) :: mesg ! The text of an error message @@ -1080,9 +1087,9 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, !### In the code below 1.0e-14 is a dimensional constant in [s-3] if ((CS%tideamp(i,j) > 0.0) .and. & (CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 > 1.0e-14) ) then - z0_polzin(i) = US%m_to_Z * CS%Polzin_decay_scale_factor * CS%Nu_Polzin * & + z0_polzin(i) = CS%Polzin_decay_scale_factor * CS%Nu_Polzin * & CS%Nbotref_Polzin**2 * CS%tideamp(i,j) / & - ( CS%kappa_itides**2 * CS%h2(i,j) * CS%Nb(i,j)**3 ) + ( CS%kappa_itides**2 * CS%h2(i,j) * US%T_to_s * CS%Nb(i,j)**3 ) if (z0_polzin(i) < CS%Polzin_min_decay_scale) & z0_polzin(i) = CS%Polzin_min_decay_scale if (N2_meanz(i) > 1.0e-14 ) then !### Here 1.0e-14 has dimensions of s-2. @@ -1137,7 +1144,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, ! Both Polzin and Simmons: do i=is,ie ! Dissipation of locally trapped internal tide (non-propagating high modes) - TKE_itidal_bot(i) = min(CS%TKE_itidal(i,j)*CS%Nb(i,j),CS%TKE_itide_max) + TKE_itidal_bot(i) = min(CS%TKE_itidal(i,j)*US%T_to_s*CS%Nb(i,j), CS%TKE_itide_max) if (associated(dd%TKE_itidal_used)) & dd%TKE_itidal_used(i,j) = TKE_itidal_bot(i) TKE_itidal_bot(i) = (I_rho0 * CS%Mu_itides * CS%Gamma_itides) * TKE_itidal_bot(i) @@ -1186,8 +1193,8 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)* TKE_frac_top_lowmode(i) ! Actual power expended may be less than predicted if stratification is weak; adjust - if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > (US%s_to_T**3 * max_TKE(i,k))) then - frac_used = (US%s_to_T**3 * max_TKE(i,k)) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > (max_TKE(i,k))) then + frac_used = (max_TKE(i,k)) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) TKE_itide_lay = frac_used * TKE_itide_lay TKE_Niku_lay = frac_used * TKE_Niku_lay TKE_lowmode_lay = frac_used * TKE_lowmode_lay @@ -1199,7 +1206,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay ! Convert power to diffusivity - Kd_add = (US%T_to_s**2 * TKE_to_Kd(i,k)) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + Kd_add = US%s_to_T * TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) Kd_lay(i,j,k) = Kd_lay(i,j,k) + (US%T_to_s * Kd_add) @@ -1213,7 +1220,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, if (associated(dd%Kd_itidal)) then ! If at layers, dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay ! The following sets the interface diagnostics. - Kd_add = (US%T_to_s**2 * TKE_to_Kd(i,k)) * TKE_itide_lay + Kd_add = US%s_to_T * TKE_to_Kd(i,k) * TKE_itide_lay if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add if (k= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add if (k= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add if (k (US%s_to_T**3 * max_TKE(i,k))) then - frac_used = (US%s_to_T**3 * max_TKE(i,k)) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > (max_TKE(i,k))) then + frac_used = (max_TKE(i,k)) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) TKE_itide_lay = frac_used * TKE_itide_lay TKE_Niku_lay = frac_used * TKE_Niku_lay TKE_lowmode_lay = frac_used * TKE_lowmode_lay @@ -1287,7 +1294,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay ! Convert power to diffusivity - Kd_add = (US%T_to_s**2 * TKE_to_Kd(i,k)) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) + Kd_add = US%s_to_T * TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay) if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) Kd_lay(i,j,k) = Kd_lay(i,j,k) + (US%T_to_s * Kd_add) @@ -1301,7 +1308,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, if (associated(dd%Kd_itidal)) then ! If at layers, this is just dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay ! The following sets the interface diagnostics. - Kd_add = (US%T_to_s**2 * TKE_to_Kd(i,k)) * TKE_itide_lay + Kd_add = US%s_to_T * TKE_to_Kd(i,k) * TKE_itide_lay if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_itidal(i,j,K) = dd%Kd_itidal(i,j,K) + 0.5*Kd_add if (k= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_Niku(i,j,K) = dd%Kd_Niku(i,j,K) + 0.5*Kd_add if (k= 0.0) Kd_add = min(Kd_add, Kd_max) if (k>1) dd%Kd_lowmode(i,j,K) = dd%Kd_lowmode(i,j,K) + 0.5*Kd_add if (k