From 9b86edb3d8f0acb13cfaa610cad17b4a740d7bd2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 12 Aug 2023 11:55:57 -0400 Subject: [PATCH] *Non-Boussinesq revision of set_diffusivity This commit revises the internal routines called by set_diffusivity to work in fully non-Boussinesq mode, eliminating all dependencies on the Boussinesq reference density when in non-Boussinesq mode. The publicly visible interfaces to this module and the external routines it calls have already been revised, so only this file needs to be updated. The specific changes include: - Refactored add_LOTW_BBL_diffusivity, add_MLrad_diffusivity, and set_BBL_TKE to work in units of layer vertical extent rather than layer thickness to give results in non-Boussinesq mode that avoid dependence on the Boussinesq reference density. - Work with internal variables in vertical distances in the denominator of diffusive flux calculations in find_TKE_to_Kd, while other expressions in this routine are recast in terms of thicknesses to avoid conversions. - Use tv%SpV_avg instead of 1/RHO_0 in find_TKE_to_Kd when in fully non-Boussinesq mode. - Use layer target density differences in place of g_prime in set_density_ratios in mathematically equivalent expressions when non-Boussinesq. - Use thickness_to_dz in 3 places to convert layer thicknesses into vertical distances. - The thickness argument to add_MLrad_diffusivity (in [H ~> m or kg m-2]) has been replaced with a vertical extent argument (in [Z ~> m]). - Use fluxes%tau_mag in place of fluxes%ustar in add_MLrad_diffusivity when in non-Boussinesq or semi-Boussinesq mode. There is a new thermo_var_ptrs type argument to the internal routine add_MLrad_diffusivity to permit this changes. - Use the in situ near-bottom density when adding certain contributions to non-Boussinesq diffusivities. This change includes the addition of a new bottom density (rho_bot) argument to add_int_tide_diffusivity, add_LOTW_BBL_diffusivity and add_drag_diffusivity. - Use GV%dZ_subroundoff in 4 spots in place of GV%H_to_Z*GV%H_subroundoff. - A long-standing comment questioning whether there is double-counting of tidal mixing has been addressed (there is not) and the comment has been revised accordingly. These changes involved changing the units of 21 internal variables and 1 element in the set_diffusivity_CS type. There are 11 new internal variables, while 9 internal variables were eliminated. A total of 44 thickness rescaling factors were eliminated, and there are 2 places where GV%Rho_0 was being used explicitly that were changed into equivalent rescaling factors. All answers are bitwise identical in Boussinesq mode, but some solutions will change in non-Boussinesq mode with this commit. --- .../vertical/MOM_set_diffusivity.F90 | 341 ++++++++++-------- 1 file changed, 200 insertions(+), 141 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 2aac478086..c792f5200e 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -80,8 +80,8 @@ module MOM_set_diffusivity real :: cdrag !< quadratic drag coefficient [nondim] real :: dz_BBL_avg_min !< A minimal distance over which to average to determine the average !! bottom boundary layer density [Z ~> m] - real :: IMax_decay !< inverse of a maximum decay scale for - !! bottom-drag driven turbulence [Z-1 ~> m-1]. + real :: IMax_decay !< Inverse of a maximum decay scale for + !! bottom-drag driven turbulence [H-1 ~> m-1 or m2 kg-1]. real :: Kv !< The interior vertical viscosity [H Z T-1 ~> m2 s-1 or Pa s] real :: Kd !< interior diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real :: Kd_min !< minimum diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -505,8 +505,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i endif ! Add the ML_Rad diffusivity. - if (CS%ML_radiation) & - call add_MLrad_diffusivity(h, fluxes, j, Kd_int_2d, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d) + if (CS%ML_radiation) then + call add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int_2d, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d) + endif ! Add the Nikurashin and / or tidal bottom-driven mixing if (CS%use_tidal_mixing) & @@ -517,11 +518,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i ! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag. if (CS%bottomdraglaw .and. (CS%BBL_effic > 0.0)) then if (CS%use_LOTW_BBL_diffusivity) then - call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int_2d, G, GV, US, CS, & - dd%Kd_BBL, Kd_lay_2d) + call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bot, Kd_int_2d, & + G, GV, US, CS, dd%Kd_BBL, Kd_lay_2d) else call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & - maxTKE, kb, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL) + maxTKE, kb, rho_bot, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL) endif endif @@ -567,7 +568,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (associated(dd%Kd_work)) then do k=1,nz ; do i=is,ie - dd%Kd_Work(i,j,k) = GV%H_to_RZ * Kd_lay_2d(i,k) * N2_lay(i,k) * GV%H_to_Z*h(i,j,k) ! Watt m-2 s = kg s-3 + dd%Kd_Work(i,j,k) = GV%H_to_RZ * Kd_lay_2d(i,k) * N2_lay(i,k) * dz(i,k) ! Watt m-2 = kg s-3 enddo ; enddo endif @@ -697,27 +698,30 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & ! across an interface times the difference across the ! interface above it [nondim] rho_0, & ! Layer potential densities relative to surface pressure [R ~> kg m-3] + dz, & ! Height change across layers [Z ~> m] maxEnt ! maxEnt is the maximum value of entrainment from below (with ! compensating entrainment from above to keep the layer ! density from changing) that will not deplete all of the - ! layers above or below a layer within a timestep [Z ~> m]. + ! layers above or below a layer within a timestep [H ~> m or kg m-2]. real, dimension(SZI_(G)) :: & htot, & ! total thickness above or below a layer, or the - ! integrated thickness in the BBL [Z ~> m]. - mFkb, & ! total thickness in the mixed and buffer layers times ds_dsp1 [Z ~> m]. + ! integrated thickness in the BBL [H ~> m or kg m-2]. + mFkb, & ! total thickness in the mixed and buffer layers times ds_dsp1 [H ~> m or kg m-2] p_ref, & ! array of tv%P_Ref pressures [R L2 T-2 ~> Pa] Rcv_kmb, & ! coordinate density in the lowest buffer layer [R ~> kg m-3] p_0 ! An array of 0 pressures [R L2 T-2 ~> Pa] real :: dh_max ! maximum amount of entrainment a layer could undergo before - ! entraining all fluid in the layers above or below [Z ~> m]. + ! entraining all fluid in the layers above or below [H ~> m or kg m-2] real :: dRho_lay ! density change across a layer [R ~> kg m-3] real :: Omega2 ! rotation rate squared [T-2 ~> s-2] - real :: G_Rho0 ! Gravitational acceleration divided by Boussinesq reference density [Z T-2 R-1 ~> m4 s-2 kg-1] - real :: G_IRho0 ! Alternate calculation of G_Rho0 for reproducibility [Z T-2 R-1 ~> m4 s-2 kg-1] - real :: I_Rho0 ! inverse of Boussinesq reference density [R-1 ~> m3 kg-1] + real :: grav ! Gravitational acceleration [Z T-1 ~> m s-2] + real :: G_Rho0 ! Gravitational acceleration divided by Boussinesq reference density + ! [Z R-1 T-2 ~> m4 s-2 kg-1] + real :: G_IRho0 ! Alternate calculation of G_Rho0 with thickness rescaling factors + ! [Z2 T-2 R-1 H-1 ~> m4 s-2 kg-1 or m7 kg-2 s-2] real :: I_dt ! 1/dt [T-1 ~> s-1] - real :: H_neglect ! negligibly small thickness [H ~> m or kg m-2] + real :: dz_neglect ! A negligibly small height change [Z ~> m] real :: hN2pO2 ! h (N^2 + Omega^2), in [Z T-2 ~> m s-2]. logical :: do_i(SZI_(G)) @@ -727,22 +731,25 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & I_dt = 1.0 / dt Omega2 = CS%omega**2 - H_neglect = GV%H_subroundoff - G_Rho0 = (US%L_to_Z**2 * GV%g_Earth) / (GV%Rho0) + dz_neglect = GV%dZ_subroundoff + grav = (US%L_to_Z**2 * GV%g_Earth) + G_Rho0 = grav / GV%Rho0 if (CS%answer_date < 20190101) then - I_Rho0 = 1.0 / (GV%Rho0) - G_IRho0 = (US%L_to_Z**2 * GV%g_Earth) * I_Rho0 + G_IRho0 = grav * GV%H_to_Z**2 * GV%RZ_to_H else - G_IRho0 = G_Rho0 + G_IRho0 = GV%H_to_Z*G_Rho0 endif + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, j, G, GV) + ! Simple but coordinate-independent estimate of Kd/TKE if (CS%simple_TKE_to_Kd) then do k=1,nz ; do i=is,ie - 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 T2 Z-1. - else; TKE_to_Kd(i,k) = 0.; endif + hN2pO2 = dz(i,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 T2 H-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 H Z2 T-3. @@ -793,18 +800,17 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & if (CS%bulkmixedlayer) then kmb = GV%nk_rho_varies do i=is,ie - htot(i) = GV%H_to_Z*h(i,j,kmb) + htot(i) = h(i,j,kmb) mFkb(i) = 0.0 - if (kb(i) < nz) & - mFkb(i) = ds_dsp1(i,kb(i)) * (GV%H_to_Z*(h(i,j,kmb) - GV%Angstrom_H)) + if (kb(i) < nz) mFkb(i) = ds_dsp1(i,kb(i)) * (h(i,j,kmb) - GV%Angstrom_H) enddo do k=1,kmb-1 ; do i=is,ie - htot(i) = htot(i) + GV%H_to_Z*h(i,j,k) - mFkb(i) = mFkb(i) + ds_dsp1(i,k+1)*(GV%H_to_Z*(h(i,j,k) - GV%Angstrom_H)) + htot(i) = htot(i) + h(i,j,k) + mFkb(i) = mFkb(i) + ds_dsp1(i,k+1)*(h(i,j,k) - GV%Angstrom_H) enddo ; enddo else do i=is,i - maxEnt(i,1) = 0.0 ; htot(i) = GV%H_to_Z*(h(i,j,1) - GV%Angstrom_H) + maxEnt(i,1) = 0.0 ; htot(i) = h(i,j,1) - GV%Angstrom_H enddo endif do k=kb_min,nz-1 ; do i=is,ie @@ -816,12 +822,12 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & else maxEnt(i,k) = ds_dsp1(i,k)*(maxEnt(i,k-1) + htot(i)) endif - htot(i) = htot(i) + GV%H_to_Z*(h(i,j,k) - GV%Angstrom_H) + htot(i) = htot(i) + (h(i,j,k) - GV%Angstrom_H) endif enddo ; enddo do i=is,ie - htot(i) = GV%H_to_Z*(h(i,j,nz) - GV%Angstrom_H) ; maxEnt(i,nz) = 0.0 + htot(i) = h(i,j,nz) - GV%Angstrom_H ; maxEnt(i,nz) = 0.0 do_i(i) = (G%mask2dT(i,j) > 0.0) enddo do k=nz-1,kb_min,-1 @@ -829,8 +835,8 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & do i=is,ie ; if (do_i(i)) then if (k This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, & - kb, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL) + kb, rho_bot, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1171,6 +1183,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, !! maximum-realizable thickness [H Z2 T-3 ~> m3 s-3 or W m-2] integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer !! layer, or -1 without a bulk mixed layer + real, dimension(SZI_(G)), intent(in) :: rho_bot !< In situ density averaged over a near-bottom + !! region [R ~> kg m-3] type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers, !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -1185,23 +1199,24 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, Rint ! coordinate density of an interface [R ~> kg m-3] real, dimension(SZI_(G)) :: & htot, & ! total thickness above or below a layer, or the - ! integrated thickness in the BBL [Z ~> m]. - rho_htot, & ! running integral with depth of density [R Z ~> kg m-2] + ! integrated thickness in the BBL [H ~> m or kg m-2]. + rho_htot, & ! running integral with depth of density [R H ~> kg m-2 or kg2 m-5] gh_sum_top, & ! BBL value of g'h that can be supported by - ! the local ustar, times R0_g [R Z ~> kg m-2] + ! the local ustar, times R0_g [R H ~> kg m-2 or kg2 m-5] Rho_top, & ! density at top of the BBL [R ~> kg m-3] TKE, & ! turbulent kinetic energy available to drive ! bottom-boundary layer mixing in a layer [H Z2 T-3 ~> m3 s-3 or W m-2] - I2decay ! inverse of twice the TKE decay scale [Z-1 ~> m-1]. + I2decay ! inverse of twice the TKE decay scale [H-1 ~> m-1 or m2 kg-1]. real :: TKE_to_layer ! TKE used to drive mixing in a layer [H Z2 T-3 ~> m3 s-3 or W m-2] real :: TKE_Ray ! TKE from layer Rayleigh drag used to drive mixing in layer [H Z2 T-3 ~> m3 s-3 or W m-2] real :: TKE_here ! TKE that goes into mixing in this layer [H Z2 T-3 ~> m3 s-3 or W m-2] real :: dRl, dRbot ! temporaries holding density differences [R ~> 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]. + real :: ustar_h ! Ustar at a thickness point rescaled into thickness + ! flux units [H T-1 ~> m s-1 or kg m-2 s-1]. real :: absf ! average absolute Coriolis parameter around a thickness point [T-1 ~> s-1] - real :: R0_g ! Rho0 / G_Earth [R T2 Z-1 ~> kg s2 m-4] + real :: R0_g ! Rho0 / G_Earth [R T2 H-1 ~> kg s2 m-4 or s2 m-1] real :: delta_Kd ! increment to Kd from the bottom boundary layer mixing [H Z T-1 ~> m2 s-1 or kg m-1 s-1] logical :: Rayleigh_drag ! Set to true if Rayleigh drag velocities ! defined in visc, on the assumption that this @@ -1221,7 +1236,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, TKE_Ray = 0.0 ; Rayleigh_drag = .false. if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true. - R0_g = GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth) + R0_g = GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth) do K=2,nz ; Rint(K) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo @@ -1231,9 +1246,14 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, ! Any turbulence that makes it into the mixed layers is assumed ! to be relatively small and is discarded. do i=is,ie - ustar_h = GV%H_to_Z*visc%ustar_BBL(i,j) - if (associated(fluxes%ustar_tidal)) & - ustar_h = ustar_h + fluxes%ustar_tidal(i,j) + ustar_h = visc%ustar_BBL(i,j) + if (associated(fluxes%ustar_tidal)) then + if (allocated(tv%SpV_avg)) then + ustar_h = ustar_h + GV%RZ_to_H*rho_bot(i) * fluxes%ustar_tidal(i,j) + else + ustar_h = ustar_h + GV%Z_to_H * fluxes%ustar_tidal(i,j) + endif + endif absf = 0.25 * ((abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + & (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1)))) if ((ustar_h > 0.0) .and. (absf > 0.5*CS%IMax_decay*ustar_h)) then @@ -1243,12 +1263,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)*(GV%H_to_Z*h(i,j,nz))) ) * & - visc%TKE_BBL(i,j) + TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)*h(i,j,nz)) ) * visc%TKE_BBL(i,j) if (associated(fluxes%TKE_tidal)) & TKE(i) = TKE(i) + fluxes%TKE_tidal(i,j) * GV%RZ_to_H * & - (CS%BBL_effic * exp(-I2decay(i)*(GV%H_to_Z*h(i,j,nz)))) + (CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz))) ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following ! Killworth & Edwards (1999) and Zilitikevich & Mironov (1996). @@ -1258,16 +1277,16 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, gh_sum_top(i) = R0_g * 400.0 * ustar_h**2 do_i(i) = (G%mask2dT(i,j) > 0.0) - htot(i) = GV%H_to_Z*h(i,j,nz) - rho_htot(i) = GV%Rlay(nz)*(GV%H_to_Z*h(i,j,nz)) + htot(i) = h(i,j,nz) + rho_htot(i) = GV%Rlay(nz)*(h(i,j,nz)) Rho_top(i) = GV%Rlay(1) if (CS%bulkmixedlayer .and. do_i(i)) Rho_top(i) = GV%Rlay(kb(i)-1) enddo do k=nz-1,2,-1 ; domore = .false. do i=is,ie ; if (do_i(i)) then - htot(i) = htot(i) + GV%H_to_Z*h(i,j,k) - rho_htot(i) = rho_htot(i) + GV%Rlay(k)*(GV%H_to_Z*h(i,j,k)) + htot(i) = htot(i) + h(i,j,k) + rho_htot(i) = rho_htot(i) + GV%Rlay(k)*(h(i,j,k)) if (htot(i)*GV%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i))) then ! The top of the mixing is in the interface atop the current layer. Rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i) @@ -1286,9 +1305,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, i_rem = i_rem + 1 ! Count the i-rows that are still being worked on. ! Apply vertical decay of the turbulent energy. This energy is ! simply lost. - TKE(i) = TKE(i) * exp(-I2decay(i) * (GV%H_to_Z*(h(i,j,k) + h(i,j,k+1)))) + TKE(i) = TKE(i) * exp(-I2decay(i) * (h(i,j,k) + h(i,j,k+1))) -! if (maxEnt(i,k) <= 0.0) cycle if (maxTKE(i,k) <= 0.0) cycle ! This is an analytic integral where diffusivity is a quadratic function of @@ -1377,7 +1395,7 @@ end subroutine add_drag_diffusivity !> Calculates a BBL diffusivity use a Prandtl number 1 diffusivity with a law of the !! wall turbulent viscosity, up to a BBL height where the energy used for mixing has !! consumed the mechanical TKE input. -subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int, & +subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bot, Kd_int, & G, GV, US, CS, Kd_BBL, Kd_lay) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -1396,6 +1414,8 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int integer, intent(in) :: j !< j-index of row to work on real, dimension(SZI_(G),SZK_(GV)+1), & intent(in) :: N2_int !< Square of Brunt-Vaisala at interfaces [T-2 ~> s-2] + real, dimension(SZI_(G)), intent(in) :: rho_bot !< In situ density averaged over a near-bottom + !! region [R ~> kg m-3] real, dimension(SZI_(G),SZK_(GV)+1), & intent(inout) :: Kd_int !< Interface net diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure @@ -1404,26 +1424,29 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int optional, intent(inout) :: Kd_lay !< Layer net diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] ! Local variables + real :: dz(SZI_(G),SZK_(GV)) ! Height change across layers [Z ~> m] real :: TKE_column ! net TKE input into the column [H Z2 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] 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] + real :: ustar ! value of ustar at a thickness point [H T-1 ~> m s-1 or kg m-2 s-1]. + real :: ustar2 ! The square of ustar [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2] real :: absf ! average absolute value of Coriolis parameter around a thickness point [T-1 ~> s-1] - real :: dh, dhm1 ! thickness of layers k and k-1, respectively [Z ~> m]. - real :: z_bot ! distance to interface k from bottom [Z ~> m]. - real :: D_minus_z ! distance to interface k from surface [Z ~> m]. - real :: total_thickness ! total thickness of water column [Z ~> m]. - real :: Idecay ! inverse of decay scale used for "Joule heating" loss of TKE with height [Z-1 ~> m-1]. + real :: dz_int ! Distance between the center of the layers around an interface [Z ~> m] + real :: z_bot ! Distance to interface K from bottom [Z ~> m] + real :: h_bot ! Total thickness between interface K and the bottom [H ~> m or kg m-2] + real :: D_minus_z ! Distance between interface k and the surface [Z ~> m] + real :: total_depth ! Total distance between the seafloor and the sea surface [Z ~> m] + real :: Idecay ! Inverse of decay scale used for "Joule heating" loss of TKE with + ! height [H-1 ~> m-1 or m2 kg-1]. real :: Kd_wall ! Law of the wall diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real :: Kd_lower ! diffusivity for lower interface [H Z T-1 ~> m2 s-1 or kg m-1 s-1] - real :: ustar_D ! u* x D [Z2 T-1 ~> m2 s-1]. + real :: ustar_D ! The extent of the water column times u* [H Z T-1 ~> m2 s-1 or Pa s]. real :: N2_min ! Minimum value of N2 to use in calculation of TKE_Kd_wall [T-2 ~> s-2] logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on ! the assumption that this extracted energy also drives diapycnal mixing. - integer :: i, k, km1 + integer :: i, k logical :: do_diag_Kd_BBL if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return @@ -1437,19 +1460,28 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true. cdrag_sqrt = sqrt(CS%cdrag) + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, j, G, GV) + do i=G%isc,G%iec ! Developed in single-column mode ! Column-wise parameters. absf = 0.25 * ((abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + & (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1)))) ! Non-zero on equator! - ! u* at the bottom [Z T-1 ~> m s-1]. - ustar = GV%H_to_Z*visc%ustar_BBL(i,j) + ! u* at the bottom [H T-1 ~> m s-1 or kg m-2 s-1]. + ustar = visc%ustar_BBL(i,j) ustar2 = ustar**2 - ! In add_drag_diffusivity(), fluxes%ustar_tidal is added in. This might be double counting - ! since ustar_BBL should already include all contributions to u*? -AJA - !### Examine the question of whether there is double counting of fluxes%ustar_tidal. - if (associated(fluxes%ustar_tidal)) ustar = ustar + fluxes%ustar_tidal(i,j) + ! In add_drag_diffusivity(), fluxes%ustar_tidal is also added in. There is no + ! double-counting because the logic surrounding the calls to add_drag_diffusivity() + ! and add_LOTW_BBL_diffusivity() only calls one of the two routines. + if (associated(fluxes%ustar_tidal)) then + if (allocated(tv%SpV_avg)) then + ustar = ustar + GV%RZ_to_H*rho_bot(i) * fluxes%ustar_tidal(i,j) + else + ustar = ustar + GV%Z_to_H * fluxes%ustar_tidal(i,j) + endif + endif ! The maximum decay scale should be something of order 200 m. We use the smaller of u*/f and ! (IMax_decay)^-1 as the decay scale. If ustar = 0, this is land so this value doesn't matter. @@ -1467,17 +1499,16 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int TKE_column = CS%BBL_effic * TKE_column ! Only use a fraction of the mechanical dissipation for mixing. TKE_remaining = TKE_column - total_thickness = ( sum(h(i,j,:)) + GV%H_subroundoff )* GV%H_to_Z ! Total column thickness [Z ~> m]. - ustar_D = ustar * total_thickness + total_depth = ( sum(dz(i,:)) + GV%dz_subroundoff ) ! Total column thickness [Z ~> m]. + ustar_D = ustar * total_depth + h_bot = 0. z_bot = 0. Kd_lower = 0. ! Diffusivity on bottom boundary. ! Work upwards from the bottom, accumulating work used until it exceeds the available TKE input ! at the bottom. - do k=GV%ke,2,-1 - dh = GV%H_to_Z * h(i,j,k) ! Thickness of this level [Z ~> m]. - km1 = max(k-1, 1) - dhm1 = GV%H_to_Z * h(i,j,km1) ! Thickness of level above [Z ~> m]. + do K=GV%ke,2,-1 + dz_int = 0.5 * (dz(i,k-1) + dz(i,k)) ! Add in additional energy input from bottom-drag against slopes (sides) if (Rayleigh_drag) TKE_remaining = TKE_remaining + & @@ -1489,23 +1520,24 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int ! Exponentially decay TKE across the thickness of the layer. ! This is energy loss in addition to work done as mixing, apparently to Joule heating. - TKE_remaining = exp(-Idecay*dh) * TKE_remaining + TKE_remaining = exp(-Idecay*h(i,j,k)) * TKE_remaining - z_bot = z_bot + h(i,j,k)*GV%H_to_Z ! Distance between upper interface of layer and the bottom [Z ~> m]. - D_minus_z = max(total_thickness - z_bot, 0.) ! Thickness above layer [Z ~> m]. + z_bot = z_bot + dz(i,k) ! Distance between upper interface of layer and the bottom [Z ~> m]. + h_bot = h_bot + h(i,j,k) ! Thickness between upper interface of layer and the bottom [H ~> m or kg m-2]. + D_minus_z = max(total_depth - z_bot, 0.) ! Thickness above layer [H ~> m or kg m-2]. - ! Diffusivity using law of the wall, limited by rotation, at height z [Z2 T-1 ~> m2 s-1]. + ! Diffusivity using law of the wall, limited by rotation, at height z [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. ! This calculation is at the upper interface of the layer - if ( ustar_D + absf * ( z_bot * D_minus_z ) == 0.) then + if ( ustar_D + absf * ( h_bot * D_minus_z ) == 0.) then Kd_wall = 0. else - Kd_wall = ((GV%Z_to_H*CS%von_karm * ustar2) * (z_bot * D_minus_z)) & - / (ustar_D + absf * (z_bot * D_minus_z)) + Kd_wall = ((CS%von_karm * ustar2) * (z_bot * D_minus_z)) & + / (ustar_D + absf * (h_bot * D_minus_z)) endif ! TKE associated with Kd_wall [H Z2 T-3 ~> m3 s-3 or W m-2]. - ! This calculation if for the volume spanning the interface. - TKE_Kd_wall = Kd_wall * 0.5 * (dh + dhm1) * max(N2_int(i,k), N2_min) + ! This calculation is for the volume spanning the interface. + TKE_Kd_wall = Kd_wall * dz_int * 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 @@ -1535,13 +1567,14 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int end subroutine add_LOTW_BBL_diffusivity !> This routine adds effects of mixed layer radiation to the layer diffusivities. -subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, Kd_lay) +subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_to_Kd, Kd_lay) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dz !< Height change across layers [Z ~> m] type(forcing), intent(in) :: fluxes !< Surface fluxes structure + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available + !! thermodynamic fields. integer, intent(in) :: j !< The j-index to work on real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -1557,24 +1590,26 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, ! This routine adds effects of mixed layer radiation to the layer diffusivities. - real, dimension(SZI_(G)) :: h_ml ! Mixed layer thickness [Z ~> m]. + real, dimension(SZI_(G)) :: h_ml ! Mixed layer thickness [Z ~> m] real, dimension(SZI_(G)) :: TKE_ml_flux ! Mixed layer TKE flux [H Z2 T-3 ~> m3 s-3 or W m-2] - real, dimension(SZI_(G)) :: I_decay ! A decay rate [Z-1 ~> m-1]. + real, dimension(SZI_(G)) :: I_decay ! A decay rate [Z-1 ~> m-1]. real, dimension(SZI_(G)) :: Kd_mlr_ml ! Diffusivities associated with mixed layer radiation ! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real :: f_sq ! The square of the local Coriolis parameter or a related variable [T-2 ~> s-2]. - real :: h_ml_sq ! The square of the mixed layer thickness [Z2 ~> m2]. + real :: h_ml_sq ! The square of the mixed layer thickness [Z2 ~> m2] + real :: u_star_H ! ustar converted to thickness based units [H T-1 ~> m s-1 or kg m-2 s-1] real :: ustar_sq ! ustar squared [Z2 T-2 ~> m2 s-2] real :: Kd_mlr ! A diffusivity associated with mixed layer turbulence radiation ! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] real :: C1_6 ! 1/6 [nondim] real :: Omega2 ! rotation rate squared [T-2 ~> s-2]. real :: z1 ! layer thickness times I_decay [nondim] - real :: dzL ! thickness converted to heights [Z ~> m]. - real :: I_decay_len2_TKE ! squared inverse decay lengthscale for - ! TKE, as used in the mixed layer code [Z-2 ~> m-2]. - real :: h_neglect ! negligibly small thickness [Z ~> m]. + real :: I_decay_len2_TKE ! Squared inverse decay lengthscale for TKE from the bulk mixed + ! layer code [Z-2 ~> m-2] + real :: dz_neglect ! A negligibly small height change [Z ~> m] logical :: do_any, do_i(SZI_(G)) integer :: i, k, is, ie, nz, kml @@ -1583,12 +1618,13 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, Omega2 = CS%omega**2 C1_6 = 1.0 / 6.0 kml = GV%nkml - h_neglect = GV%H_subroundoff*GV%H_to_Z + dz_neglect = GV%dz_subroundoff + 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. if (.not.CS%ML_radiation) return do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo - do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + GV%H_to_Z*h(i,j,k) ; enddo ; enddo + do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + dz(i,k) ; enddo ; enddo do i=is,ie ; if (do_i(i)) then if (CS%ML_omega_frac >= 1.0) then @@ -1600,21 +1636,31 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, f_sq = CS%ML_omega_frac * 4.0 * Omega2 + (1.0 - CS%ML_omega_frac) * f_sq endif - ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2 - - TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * (GV%Z_to_H*fluxes%ustar(i,j))) + ! Determine the energy flux out of the mixed layer and its vertical decay scale. + if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then + ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2 + u_star_H = GV%Z_to_H * fluxes%ustar(i,j) + elseif (allocated(tv%SpV_avg)) then + ustar_sq = max(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2) + u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + else ! This semi-Boussinesq form is mathematically equivalent to the Boussinesq version above. + ! Differs at roundoff: ustar_sq = max(fluxes%tau_mag(i,j) * I_rho, CS%ustar_min**2) + ustar_sq = max((sqrt(fluxes%tau_mag(i,j) * I_rho))**2, CS%ustar_min**2) + u_star_H = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * GV%Rho0) + endif + TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * u_star_H) I_decay_len2_TKE = CS%TKE_decay**2 * (f_sq / ustar_sq) if (CS%ML_rad_TKE_decay) & TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-h_ml(i) * sqrt(I_decay_len2_TKE)) ! Calculate the inverse decay scale - h_ml_sq = (CS%ML_rad_efold_coeff * (h_ml(i)+h_neglect))**2 + h_ml_sq = (CS%ML_rad_efold_coeff * (h_ml(i)+dz_neglect))**2 I_decay(i) = sqrt((I_decay_len2_TKE * h_ml_sq + 1.0) / h_ml_sq) ! Average the dissipation layer kml+1, using ! a more accurate Taylor series approximations for very thin layers. - z1 = (GV%H_to_Z*h(i,j,kml+1)) * I_decay(i) + z1 = dz(i,kml+1) * I_decay(i) if (z1 > 1e-5) then Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) * (1.0 - exp(-z1)) else @@ -1639,14 +1685,14 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, Kd_int, G, GV, US, CS, TKE_to_Kd, do k=kml+2,nz-1 do_any = .false. do i=is,ie ; if (do_i(i)) then - dzL = GV%H_to_Z*h(i,j,k) ; z1 = dzL*I_decay(i) + z1 = dz(i,k)*I_decay(i) if (CS%ML_Rad_bug) then ! These expressions are dimensionally inconsistent. -RWH ! This is supposed to be the integrated energy deposited in the layer, ! not the average over the layer as in these expressions. if (z1 > 1e-5) then Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of H Z T-1 - US%m_to_Z * ((1.0 - exp(-z1)) / dzL) ! Units of m-1 + US%m_to_Z * ((1.0 - exp(-z1)) / dz(i,k)) ! Units of m-1 else Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of H Z T-1 US%m_to_Z * (I_decay(i) * (1.0 - z1 * (0.5 - C1_6*z1))) ! Units of m-1 @@ -1698,23 +1744,23 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) ! boundary layer turbulence. real, dimension(SZI_(G)) :: & - htot ! total thickness above or below a layer, or the - ! integrated thickness in the BBL [Z ~> m]. + htot ! Running sum of the depth in the BBL [Z ~> m]. real, dimension(SZIB_(G)) :: & uhtot, & ! running integral of u in the BBL [Z L T-1 ~> m2 s-1] - ustar, & ! bottom boundary layer turbulence speed [Z T-1 ~> m s-1]. + ustar, & ! bottom boundary layer piston velocity [H T-1 ~> m s-1 or kg m-2 s-1]. u2_bbl ! square of the mean zonal velocity in the BBL [L2 T-2 ~> m2 s-2] real :: vhtot(SZI_(G)) ! running integral of v in the BBL [Z L T-1 ~> m2 s-1] real, dimension(SZI_(G),SZJB_(G)) :: & - vstar, & ! ustar at at v-points [Z T-1 ~> m s-1]. + vstar, & ! ustar at at v-points [H T-1 ~> m s-1 or kg m-2 s-1]. v2_bbl ! square of average meridional velocity in BBL [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & + dz ! The vertical distance between interfaces around a layer [Z ~> m] - real :: cdrag_sqrt ! square root of the drag coefficient [nondim] - real :: I_cdrag_sqrt ! The inverse of the square root of the drag coefficient [nondim] - real :: hvel ! thickness at velocity points [Z ~> m]. + real :: cdrag_sqrt ! Square root of the drag coefficient [nondim] + real :: hvel ! thickness at velocity points [Z ~> m] logical :: domore, do_i(SZI_(G)) integer :: i, j, k, is, ie, js, je, nz @@ -1748,7 +1794,9 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) endif cdrag_sqrt = sqrt(CS%cdrag) - I_cdrag_sqrt = 0.0 ; if (cdrag_sqrt > 0.0) I_cdrag_sqrt = 1.0 / cdrag_sqrt + + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) !$OMP parallel default(shared) private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl) !$OMP do @@ -1761,7 +1809,7 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (allocated(visc%Kv_bbl_v)) then do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then do_i(i) = .true. - vstar(i,J) = GV%H_to_Z*visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J)) + vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J)) endif ; enddo endif !### What about terms from visc%Ray? @@ -1781,12 +1829,12 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) ! Compute h based on OBC state if (has_obc) then if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then - hvel = GV%H_to_Z*h(i,j,k) + hvel = dz(i,j,k) else - hvel = GV%H_to_Z*h(i,j+1,k) + hvel = dz(i,j+1,k) endif else - hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i,j+1,k)) + hvel = 0.5*(dz(i,j,k) + dz(i,j+1,k)) endif if ((htot(i) + hvel) >= visc%bbl_thick_v(i,J)) then @@ -1802,7 +1850,7 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (.not.domore) exit enddo do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (htot(i) > 0.0)) then - v2_bbl(i,J) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i)) + v2_bbl(i,J) = (vhtot(i)*vhtot(i)) / (htot(i)*htot(i)) else v2_bbl(i,J) = 0.0 endif ; enddo @@ -1815,7 +1863,7 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (allocated(visc%bbl_thick_u)) then do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0)) then do_i(I) = .true. - ustar(I) = GV%H_to_Z*visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j)) + ustar(I) = visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j)) endif ; enddo endif @@ -1833,12 +1881,12 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) ! Compute h based on OBC state if (has_obc) then if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then - hvel = GV%H_to_Z*h(i,j,k) + hvel = dz(i,j,k) else ! OBC_DIRECTION_W - hvel = GV%H_to_Z*h(i+1,j,k) + hvel = dz(i+1,j,k) endif else - hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i+1,j,k)) + hvel = 0.5*(dz(i,j,k) + dz(i+1,j,k)) endif if ((htot(I) + hvel) >= visc%bbl_thick_u(I,j)) then @@ -1854,18 +1902,18 @@ subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC) if (.not.domore) exit enddo do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (htot(i) > 0.0)) then - u2_bbl(I) = (uhtot(I)*uhtot(I))/(htot(I)*htot(I)) + u2_bbl(I) = (uhtot(I)*uhtot(I)) / (htot(I)*htot(I)) else u2_bbl(I) = 0.0 endif ; enddo do i=is,ie - visc%ustar_BBL(i,j) = GV%Z_to_H*sqrt(0.5*G%IareaT(i,j) * & + visc%ustar_BBL(i,j) = sqrt(0.5*G%IareaT(i,j) * & ((G%areaCu(I-1,j)*(ustar(I-1)*ustar(I-1)) + & 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) = GV%Z_to_H*US%L_to_Z**2 * & + visc%TKE_BBL(i,j) = US%L_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)) + & @@ -1904,7 +1952,8 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) real :: a(SZK_(GV)), a_0(SZK_(GV)) ! nondimensional temporary variables [nondim] real :: p_ref(SZI_(G)) ! an array of tv%P_Ref pressures [R L2 T-2 ~> Pa] real :: Rcv(SZI_(G),SZK_(GV)) ! coordinate density in the mixed and buffer layers [R ~> kg m-3] - real :: I_Drho ! temporary variable [R-1 ~> m3 kg-1] + real :: I_Drho ! The inverse of the coordinate density difference between + ! layers [R-1 ~> m3 kg-1] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, k, k3, is, ie, nz, kmb @@ -1912,9 +1961,15 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) do k=2,nz-1 if (GV%g_prime(k+1) /= 0.0) then - do i=is,ie - ds_dsp1(i,k) = GV%g_prime(k) / GV%g_prime(k+1) - enddo + if (GV%Boussinesq .or. GV%Semi_Boussinesq) then + do i=is,ie + ds_dsp1(i,k) = GV%g_prime(k) / GV%g_prime(k+1) + enddo + else ! Use a mathematically equivalent form that avoids any dependency on RHO_0. + do i=is,ie + ds_dsp1(i,k) = (GV%Rlay(k) - GV%Rlay(k-1)) / (GV%Rlay(k+1) - GV%Rlay(k)) + enddo + endif else do i=is,ie ds_dsp1(i,k) = 1. @@ -1937,7 +1992,11 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) ! interfaces above and below the buffer layer and the next denser layer. k = kb(i) - I_Drho = g_R0 / GV%g_prime(k+1) + if (GV%Boussinesq .or. GV%Semi_Boussinesq) then + I_Drho = g_R0 / GV%g_prime(k+1) + else + I_Drho = 1.0 / (GV%Rlay(k+1) - GV%Rlay(k)) + endif ! The indexing convention for a is appropriate for the interfaces. do k3=1,kmb a(k3+1) = (GV%Rlay(k) - Rcv(i,k3)) * I_Drho @@ -2005,7 +2064,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ !! surface boundary layer. ! Local variables - real :: decay_length ! The maximum decay scale for the BBL diffusion [Z ~> m] + real :: decay_length ! The maximum decay scale for the BBL diffusion [H ~> m or kg m-2] logical :: ML_use_omega integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. @@ -2097,7 +2156,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "length scale.", default=.false.) if (CS%ML_radiation) then ! This give a minimum decay scale that is typically much less than Angstrom. - CS%ustar_min = 2e-4 * CS%omega * (GV%Angstrom_Z + GV%H_subroundoff*GV%H_to_Z) + CS%ustar_min = 2e-4 * CS%omega * (GV%Angstrom_Z + GV%dZ_subroundoff) call get_param(param_file, mdl, "ML_RAD_EFOLD_COEFF", CS%ML_rad_efold_coeff, & "A coefficient that is used to scale the penetration "//& @@ -2161,7 +2220,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//& "to penetrate as far as stratification and rotation permit. The default "//& "for now is 200 m. This is only used if BOTTOMDRAGLAW is true.", & - units="m", default=200.0, scale=US%m_to_Z) + units="m", default=200.0, scale=GV%m_to_H) CS%IMax_decay = 0.0 if (decay_length > 0.0) CS%IMax_decay = 1.0/decay_length