diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index b207b1ff1c..bf0db51a44 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -20,6 +20,7 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_interface_heights, only : thickness_to_dz +! use MOM_intrinsic_functions, only : cuberoot use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE, OBC_DIRECTION_E @@ -258,40 +259,15 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: crv ! crv is the curvature of the bottom depth across a ! cell, times the cell width squared [Z ~> m]. real :: crv_3 ! crv/3 [Z ~> m]. - real :: C24_crv ! 24/crv [Z-1 ~> m-1]. real :: slope ! The absolute value of the bottom depth slope across ! a cell times the cell width [Z ~> m]. - real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim]. - real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] - ! All of the following "volumes" have units of vertical heights because they are normalized - ! by the full horizontal area of a velocity cell. real :: Vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel ! drag parameterization, normalized by the full horizontal area ! of the velocity cell [Z ~> m]. - real :: Vol_open ! The cell volume above which it is open [Z ~> m]. - real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct - ! solution of a cubic equation for L. - real :: Vol_2_reg ! The cell volume above which there are two separate - ! open areas that must be integrated [Z ~> m]. - real :: vol ! The volume below the interface whose normalized - ! width is being sought [Z ~> m]. - real :: vol_below ! The volume below the interface below the one that - ! is currently under consideration [Z ~> m]. - real :: Vol_err ! The error in the volume with the latest estimate of - ! L, or the error for the interface below [Z ~> m]. - real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m]. - real :: Vol_tol ! A volume error tolerance [Z ~> m]. + real :: vol_below(SZK_(GV)+1) ! The volume below each interface, normalized by the full + ! horizontal area of a velocity cell [Z ~> m]. real :: L(SZK_(GV)+1) ! The fraction of the full cell width that is open at ! the depth of each interface [nondim]. - real :: L_direct ! The value of L above volume Vol_direct [nondim]. - real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim]. - real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m] - real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m] - real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m]. - real :: L0 ! The value of L above volume Vol_0 [nondim]. - real :: dVol ! vol - Vol_0 [Z ~> m]. - real :: dV_dL2 ! The partial derivative of volume with L squared - ! evaluated at L=L0 [Z ~> m]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -315,14 +291,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] real :: C2pi_3 ! An irrational constant, 2/3 pi. [nondim] real :: tmp ! A temporary variable, sometimes in [Z ~> m] - real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] - real :: curv_tol ! Numerator of curvature cubed, used to estimate - ! accuracy of a single L(:) Newton iteration [Z5 ~> m5] - logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration logical :: use_BBL_EOS, do_i(SZIB_(G)) integer, dimension(2) :: EOSdom ! The computational domain for the equation of state integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml - integer :: itt, maxitt=20 type(ocean_OBC_type), pointer :: OBC => NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -332,7 +303,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) dz_neglect = GV%dZ_subroundoff Rho0x400_G = 400.0*(GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth)) - Vol_quit = (0.9*GV%Angstrom_Z + dz_neglect) C2pi_3 = 8.0*atan(1.0)/3.0 if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//& @@ -442,8 +412,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G,C2pi_3, & !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, & !$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, & - !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v,pbv) & - !$OMP firstprivate(Vol_quit) + !$OMP OBC,D_u,D_v,mask_u,mask_v,pbv) do j=Jsq,Jeq ; do m=1,2 if (m==1) then @@ -865,12 +834,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (CS%body_force_drag) bbl_thick = dz_bbl_drag(i) if (CS%Channel_drag) then - ! The drag within the bottommost Vol_bbl_chan is applied as a part of - ! an enhanced bottom viscosity, while above this the drag is applied - ! directly to the layers in question as a Rayleigh drag term. - ! Restrict the volume over which the channel drag is applied. - if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) + vol_below(nz+1) = 0.0 + do K=nz,1,-1 + vol_below(K) = vol_below(K+1) + dz_vel(i,k) + enddo !### The harmonic mean edge depths here are not invariant to offsets! if (m==1) then @@ -887,162 +855,33 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) Dm = 2.0 * D_vel * tmp / (D_vel + tmp) endif if (Dm > Dp) then ; tmp = Dp ; Dp = Dm ; Dm = tmp ; endif - crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 slope = Dp - Dm + ! If the curvature is small enough, there is no reason not to assume ! a uniformly sloping or flat bottom. if (abs(crv) < 1e-2*(slope + CS%BBL_thick_min)) crv = 0.0 - ! Each cell extends from x=-1/2 to 1/2, and has a topography - ! given by D(x) = crv*x^2 + slope*x + D - crv/12. - - ! Calculate the volume above which the entire cell is open and the - ! other volumes at which the equation that is solved for L changes. - if (crv > 0.0) then - if (slope >= crv) then - Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open - else - tmp = slope/crv - Vol_open = 0.25*slope*tmp + C1_12*crv - Vol_2_reg = 0.5*tmp**2 * (crv - C1_3*slope) - endif - ! Define some combinations of crv & slope for later use. - C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) - apb_4a = (slope+crv)/(4.0*crv) ; a2x48_apb3 = (48.0*(crv*crv))*(Iapb**3) - ax2_3apb = 2.0*C1_3*crv*Iapb - elseif (crv == 0.0) then - Vol_open = 0.5*slope - if (slope > 0) Iapb = 1.0/slope - else ! crv < 0.0 - Vol_open = D_vel - Dm - if (slope >= -crv) then - Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope) - Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0 - else - C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) - L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0 - Vol_direct = -C1_6*crv*L_direct**3 - endif - Ibma_2 = 2.0 / (slope - crv) - endif - L(nz+1) = 0.0 ; vol = 0.0 ; Vol_err = 0.0 ; BBL_visc_frac = 0.0 - ! Determine the normalized open length at each interface. - do K=nz,1,-1 - vol_below = vol - - vol = vol + dz_vel(i,k) - h_vel_pos = h_vel(i,k) + h_neglect - - if (vol >= Vol_open) then ; L(K) = 1.0 - elseif (crv == 0) then ! The bottom has no curvature. - L(K) = sqrt(2.0*vol*Iapb) - elseif (crv > 0) then - ! There may be a minimum depth, and there are - ! analytic expressions for L for all cases. - if (vol < Vol_2_reg) then - ! In this case, there is a contiguous open region and - ! vol = 0.5*L^2*(slope + crv/3*(3-4L)). - if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7? - ! There is a very good approximation here for massless layers. - L0 = sqrt(2.0*vol*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) - else - L(K) = apb_4a * (1.0 - & - 2.0 * cos(C1_3*acos(a2x48_apb3*vol - 1.0) - C2pi_3)) - endif - ! To check the answers. - ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - else ! There are two separate open regions. - ! vol = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L) - ! At the deepest volume, L = slope/crv, at the top L = 1. - !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol)) - C2pi_3) - tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol) - tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1)) - L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3) - ! To check the answers. - ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol - endif - else ! a < 0. - if (vol <= Vol_direct) then - ! Both edges of the cell are bounded by walls. - L(K) = (-0.25*C24_crv*vol)**C1_3 - else - ! x_R is at 1/2 but x_L is in the interior & L is found by solving - ! vol = 0.5*L^2*(slope + crv/3*(3-4L)) - - ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + crv_3*(3.0-4.0*L(K+1))) - vol_below - ! Change to ... - ! if (min(vol_below + Vol_err, vol) <= Vol_direct) then ? - if (vol_below + Vol_err <= Vol_direct) then - L0 = L_direct ; Vol_0 = Vol_direct - else - L0 = L(K+1) ; Vol_0 = vol_below + Vol_err - ! Change to Vol_0 = min(vol_below + Vol_err, vol) ? - endif - - ! Try a relatively simple solution that usually works well - ! for massless layers. - dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = (vol-Vol_0) - ! dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = max(vol-Vol_0, 0.0) - - use_L0 = .false. - do_one_L_iter = .false. - if (CS%answer_date < 20190101) then - curv_tol = GV%Angstrom_Z*dV_dL2**2 & - * (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol) - do_one_L_iter = (crv * crv * dVol**3) < curv_tol - else - ! The following code is more robust when GV%Angstrom_H=0, but - ! it changes answers. - use_L0 = (dVol <= 0.) - - Vol_tol = max(0.5 * GV%Angstrom_Z + dz_neglect, 1e-14 * vol) - Vol_quit = max(0.9 * GV%Angstrom_Z + dz_neglect, 1e-14 * vol) + ! Determine the normalized open length (L) at each interface. + if (crv == 0.0) then + call find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) + elseif (crv > 0.0) then + call find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + else ! crv < 0.0 + call find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) + endif ! end of crv<0 cases. - curv_tol = Vol_tol * dV_dL2**2 & - * (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol) - do_one_L_iter = (crv * crv * dVol**3) < curv_tol - endif + ! Determine the Rayliegh drag contributions. - if (use_L0) then - L(K) = L0 - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - elseif (do_one_L_iter) then - ! One iteration of Newton's method should give an estimate - ! that is accurate to within Vol_tol. - L(K) = sqrt(L0*L0 + dVol / dV_dL2) - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - else - if (dV_dL2*(1.0-L0*L0) < dVol + & - dV_dL2 * (Vol_open - Vol)*Ibma_2) then - L_max = sqrt(1.0 - (Vol_open - Vol)*Ibma_2) - else - L_max = sqrt(L0*L0 + dVol / dV_dL2) - endif - L_min = sqrt(L0*L0 + dVol / (0.5*(slope+crv) - crv*L_max)) + ! The drag within the bottommost Vol_bbl_chan is applied as a part of an enhanced bottom + ! viscosity, while above this the drag is applied directly to the layers in question as a + ! Rayleigh drag term. - Vol_err_min = 0.5*(L_min**2)*(slope + crv_3*(3.0-4.0*L_min)) - vol - Vol_err_max = 0.5*(L_max**2)*(slope + crv_3*(3.0-4.0*L_max)) - vol - ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then - if (abs(Vol_err_min) <= Vol_quit) then - L(K) = L_min ; Vol_err = Vol_err_min - else - L(K) = sqrt((L_min**2*Vol_err_max - L_max**2*Vol_err_min) / & - (Vol_err_max - Vol_err_min)) - do itt=1,maxitt - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - if (abs(Vol_err) <= Vol_quit) exit - ! Take a Newton's method iteration. This equation has proven - ! robust enough not to need bracketing. - L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0*crv*L(K))) - ! This would be a Newton's method iteration for L^2: - ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+crv) - crv*L(K))) - enddo - endif ! end of iterative solver - endif ! end of 1-boundary alternatives. - endif ! end of a<0 cases. - endif + ! Restrict the volume over which the channel drag is applied from the previously determined valu.e + if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) + BBL_visc_frac = 0.0 + do K=nz,1,-1 !modify L(K) for porous barrier parameterization if (m==1) then ; L(K) = L(K)*pbv%por_layer_widthU(I,j,K) else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K); endif @@ -1050,8 +889,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! Determine the drag contributing to the bottom boundary layer ! and the Rayleigh drag that acts on each layer. if (L(K) > L(K+1)) then - if (vol_below < Vol_bbl_chan) then - BBL_frac = (1.0-vol_below/Vol_bbl_chan)**2 + if (vol_below(K+1) < Vol_bbl_chan) then + BBL_frac = (1.0-vol_below(K+1)/Vol_bbl_chan)**2 BBL_visc_frac = BBL_visc_frac + BBL_frac*(L(K) - L(K+1)) else BBL_frac = 0.0 @@ -1063,6 +902,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) cdrag_conv = cdrag_L_to_H endif + h_vel_pos = h_vel(i,k) + h_neglect if (m==1) then ; Cell_width = G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k) else ; Cell_width = G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k) ; endif gam = 1.0 - L(K+1)/L(K) @@ -1085,7 +925,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) else ; visc%Ray_v(i,J,k) = 0.0 ; endif endif - enddo ! k loop to determine L(K). + enddo ! k loop to determine visc%Ray_[uv]. ! Set the near-bottom viscosity to a value which will give ! the correct stress when the shear occurs over bbl_thick. @@ -1202,6 +1042,299 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) end subroutine set_viscous_BBL +!> Determine the normalized open length of each interface, given the edge depths and normalized +!! volumes below each interface. +subroutine find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + + ! Local variables + real :: slope ! The absolute value of the bottom depth slope across a cell times the cell width [Z ~> m]. + real :: I_slope ! The inverse of the normalized slope [Z-1 ~> m-1] + real :: Vol_open ! The cell volume above which it is open [Z ~> m]. + integer :: K, nz + + nz = GV%ke + + slope = abs(Dp - Dm) + if (slope == 0.0) then + L(1:nz) = 1.0 ; L(nz+1) = 0.0 + else + Vol_open = 0.5*slope + I_slope = 1.0 / slope + + L(nz+1) = 0.0 + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then ; L(K) = 1.0 + else + ! With a uniformly sloping bottom, the calculation of L(K) is the solution of a simple quadratic equation. + L(K) = sqrt(2.0*vol_below(K)*I_slope) + endif + enddo + endif + +end subroutine find_L_open_uniform_slope + +!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using +!! analytic expressions. In this case there can be two separate open regions. +subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: C2pi_3 !< An irrational constant, 2/3 pi [nondim] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + ! The following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m]. + real :: Vol_2_reg ! The cell volume above which there are two separate + ! open areas that must be integrated [Z ~> m]. + real :: C24_crv ! 24/crv [Z-1 ~> m-1]. + real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim]. + real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] + real :: L0 ! A linear estimate of L approprate for tiny volumes [nondim]. + real :: slope_crv ! The slope divided by the curvature [nondim] + real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] + real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] + integer :: K, nz + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv < 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there are two separate open regions. + if (slope >= crv) then + Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open + else + slope_crv = slope / crv + Vol_open = 0.25*slope*slope_crv + C1_12*crv + Vol_2_reg = 0.5*slope_crv**2 * (crv - C1_3*slope) + endif + ! Define some combinations of crv & slope for later use. + C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) + apb_4a = (slope+crv)/(4.0*crv) ; a2x48_apb3 = (48.0*(crv*crv))*(Iapb**3) + ax2_3apb = 2.0*C1_3*crv*Iapb + + L(nz+1) = 0.0 + ! Determine the normalized open length (L) at each interface. + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then ! The whole cell is open. + L(K) = 1.0 + elseif (vol_below(K) < Vol_2_reg) then + ! In this case, there is a contiguous open region and + ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)). + if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7? + ! There is a very good approximation here for massless layers. + L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) + else + L(K) = apb_4a * (1.0 - & + 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)) + endif + ! To check the answers. + ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + else ! There are two separate open regions. + ! vol_below(K) = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L) + ! At the deepest volume, L = slope/crv, at the top L = 1. + ! L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol_below(K))) - C2pi_3) + tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol_below(K)) + tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1)) + L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3) + ! To check the answers. + ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol_below(K) + endif + enddo ! k loop to determine L(K) in the concave case + +end subroutine find_L_open_concave_analytic + + +!> Determine the normalized open length of each interface for convex bathymetry (from the ocean perspective) using +!! Newton's method iterations. In this case there is a single open region with the minimum depth +!! at one edge of the cell. +subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(set_visc_CS), intent(in) :: CS !< The control structure returned by a previous + !! call to set_visc_init. + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + ! All of the following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_err ! The error in the volume with the latest estimate of + ! L, or the error for the interface below [Z ~> m]. + real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m]. + real :: Vol_tol ! A volume error tolerance [Z ~> m]. + real :: Vol_open ! The cell volume above which the face is fully open [Z ~> m]. + real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct + ! solution of a cubic equation for L. + real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m] + real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m] + real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m]. + real :: dVol ! vol - Vol_0 [Z ~> m]. + real :: dV_dL2 ! The partial derivative of volume with L squared + ! evaluated at L=L0 [Z ~> m]. + real :: L_direct ! The value of L above volume Vol_direct [nondim]. + real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim]. + real :: L0 ! The value of L above volume Vol_0 [nondim]. + real :: Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] + real :: C24_crv ! 24/crv [Z-1 ~> m-1]. + real :: curv_tol ! Numerator of curvature cubed, used to estimate + ! accuracy of a single L(:) Newton iteration [Z5 ~> m5] + real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0 ! Rational constants [nondim] + logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration + integer :: K, nz, itt, maxitt=20 + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv > 0.0) call MOM_error(FATAL, "find_L_open_convex should only be called with a negative curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there is a direct solution. + Vol_open = D_vel - Dm + if (slope >= -crv) then + Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope) + Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0 + else + C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) + L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0 + Vol_direct = -C1_6*crv*L_direct**3 + endif + Ibma_2 = 2.0 / (slope - crv) + + if (CS%answer_date < 20190101) Vol_quit = (0.9*GV%Angstrom_Z + GV%dZ_subroundoff) + + L(nz+1) = 0.0 ; Vol_err = 0.0 + ! Determine the normalized open length (L) at each interface. + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then + L(K) = 1.0 + elseif (vol_below(K) <= Vol_direct) then + ! Both edges of the cell are bounded by walls. + ! if (CS%answer_date < 20240101)) then + L(K) = (-0.25*C24_crv*vol_below(K))**C1_3 + ! else + ! L(K) = cuberoot(-0.25*C24_crv*vol_below(K)) + ! endif + else + ! x_R is at 1/2 but x_L is in the interior & L is found by iteratively solving + ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)) + + ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + crv_3*(3.0-4.0*L(K+1))) - vol_below(K+1) + ! Change to ... + ! if (min(vol_below(K+1) + Vol_err, vol_below(K)) <= Vol_direct) then ? + if (vol_below(K+1) + Vol_err <= Vol_direct) then + L0 = L_direct ; Vol_0 = Vol_direct + else + L0 = L(K+1) ; Vol_0 = vol_below(K+1) + Vol_err + ! Change to Vol_0 = min(vol_below(K+1) + Vol_err, vol_below(K)) ? + endif + + ! Try a relatively simple solution that usually works well + ! for massless layers. + dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = (vol_below(K)-Vol_0) + ! dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = max(vol_below(K)-Vol_0, 0.0) + + use_L0 = .false. + do_one_L_iter = .false. + if (CS%answer_date < 20190101) then + curv_tol = GV%Angstrom_Z*dV_dL2**2 & + * (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol) + do_one_L_iter = (crv * crv * dVol**3) < curv_tol + else + ! The following code is more robust when GV%Angstrom_H=0, but + ! it changes answers. + use_L0 = (dVol <= 0.) + + Vol_tol = max(0.5 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K)) + Vol_quit = max(0.9 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K)) + + curv_tol = Vol_tol * dV_dL2**2 & + * (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol) + do_one_L_iter = (crv * crv * dVol**3) < curv_tol + endif + + if (use_L0) then + L(K) = L0 + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + elseif (do_one_L_iter) then + ! One iteration of Newton's method should give an estimate + ! that is accurate to within Vol_tol. + L(K) = sqrt(L0*L0 + dVol / dV_dL2) + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + else + if (dV_dL2*(1.0-L0*L0) < dVol + & + dV_dL2 * (Vol_open - vol_below(K))*Ibma_2) then + L_max = sqrt(1.0 - (Vol_open - vol_below(K))*Ibma_2) + else + L_max = sqrt(L0*L0 + dVol / dV_dL2) + endif + L_min = sqrt(L0*L0 + dVol / (0.5*(slope+crv) - crv*L_max)) + + Vol_err_min = 0.5*(L_min**2)*(slope + crv_3*(3.0-4.0*L_min)) - vol_below(K) + Vol_err_max = 0.5*(L_max**2)*(slope + crv_3*(3.0-4.0*L_max)) - vol_below(K) + ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then + if (abs(Vol_err_min) <= Vol_quit) then + L(K) = L_min ; Vol_err = Vol_err_min + else + L(K) = sqrt((L_min**2*Vol_err_max - L_max**2*Vol_err_min) / & + (Vol_err_max - Vol_err_min)) + do itt=1,maxitt + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + if (abs(Vol_err) <= Vol_quit) exit + ! Take a Newton's method iteration. This equation has proven + ! robust enough not to need bracketing. + L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0*crv*L(K))) + ! This would be a Newton's method iteration for L^2: + ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+crv) - crv*L(K))) + enddo + endif ! end of iterative solver + endif ! end of 1-boundary alternatives. + endif ! end of 0, 1- and 2- boundary cases. + enddo ! k loop to determine L(K) in the convex case + +end subroutine find_L_open_convex + !> This subroutine finds a thickness-weighted value of v at the u-points. function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure