diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 77a2b4ea8b..5b19a7549c 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -880,21 +880,20 @@ subroutine check_remapping_grid( G, GV, h, dzInterface, msg ) !$OMP parallel do default(shared) do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, GV%Z_to_H*G%bathyT(i,j), h(i,j,:), dzInterface(i,j,:), msg ) + if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, h(i,j,:), dzInterface(i,j,:), msg ) enddo ; enddo end subroutine check_remapping_grid !> Check that the total thickness of new and old grids are consistent -subroutine check_grid_column( nk, depth, h, dzInterface, msg ) +subroutine check_grid_column( nk, h, dzInterface, msg ) integer, intent(in) :: nk !< Number of cells - real, intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h) character(len=*), intent(in) :: msg !< Message to append to errors ! Local variables integer :: k - real :: eps, total_h_old, total_h_new, h_new, z_old, z_new + real :: eps, total_h_old, total_h_new, h_new eps =1. ; eps = epsilon(eps) @@ -904,13 +903,8 @@ subroutine check_grid_column( nk, depth, h, dzInterface, msg ) total_h_old = total_h_old + h(k) enddo - ! Integrate upwards for the interfaces consistent with the rest of MOM6 - z_old = - depth - if (depth == 0.) z_old = - total_h_old total_h_new = 0. do k = nk,1,-1 - z_old = z_old + h(k) ! Old interface position above layer k - z_new = z_old + dzInterface(k) ! New interface position based on dzInterface h_new = h(k) + ( dzInterface(k) - dzInterface(k+1) ) ! New thickness if (h_new<0.) then write(0,*) 'k,h,hnew=',k,h(k),h_new @@ -1082,7 +1076,7 @@ subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) end subroutine filtered_grid_motion -!> Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004). +!> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004). !! z* is defined as !! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H . subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) @@ -1118,7 +1112,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) cycle endif - ! Local depth (G%bathyT is positive) + ! Local depth (G%bathyT is positive downward) nominalDepth = G%bathyT(i,j)*GV%Z_to_H ! Determine water column thickness @@ -1319,7 +1313,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel endif - ! Local depth (G%bathyT is positive) + ! Local depth (G%bathyT is positive downward) nominalDepth = G%bathyT(i,j)*GV%Z_to_H ! Determine total water column thickness @@ -1406,7 +1400,7 @@ end subroutine build_rho_grid !! density interpolated from the column profile and a clipping of depth for !! each interface to a fixed z* or p* grid. This should probably be (optionally?) !! changed to find the nearest location of the target density. -!! \remark { Based on Bleck, 2002: An oceanice general circulation model framed in +!! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 } subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) @@ -1575,7 +1569,9 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] - real :: depth + + ! Local variables + real :: depth ! Depth of the ocean relative to the mean sea surface height in thickness units [H ~> m or kg m-2] integer :: i, j, k, nz real :: h_neglect, h_neglect_edge @@ -1631,8 +1627,8 @@ end subroutine build_grid_SLight subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) type(regridding_CS), intent(in) :: CS !< Regridding control structure integer, intent(in) :: nk !< Number of layers in h_old - real, dimension(nk), intent(in) :: h_old !< Minium allowed thickness of h [H ~> m or kg m-2] - real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h [H ~> m or kg m-2] + real, dimension(nk), intent(in) :: h_old !< Minimum allowed thickness of h [H ~> m or kg m-2] + real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minimum allowed thickness of h [H ~> m or kg m-2] ! Local variables integer :: k real :: h_new, eps, h_total, h_err @@ -1710,8 +1706,8 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) real :: total_height real :: delta_h real :: max_depth - real :: eta ! local elevation - real :: local_depth + real :: eta ! local elevation [H ~> m or kg m-2] + real :: local_depth ! The local ocean depth relative to mean sea level in thickness units [H ~> m or kg m-2] real :: x1, y1, x2, y2 real :: x, t @@ -1769,7 +1765,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) endif enddo - ! Chnage in interface position + ! Change in interface position x = 0. ! Left boundary at x=0 dzInterface(i,j,1) = 0. do k = 2,nz @@ -1797,7 +1793,7 @@ subroutine inflate_vanished_layers_old( CS, G, GV, h ) ! objective is to make sure all layers are at least as thick as the minimum ! thickness allowed for regridding purposes (this parameter is set in the ! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they -! are inflated up to the minmum thickness. +! are inflated up to the minimum thickness. !------------------------------------------------------------------------------ ! Arguments @@ -1901,7 +1897,7 @@ function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) ! Arguments integer, intent(in) :: nk !< Number of cells in source grid character(len=*), intent(in) :: coordMode !< A string indicating the coordinate mode. - !! See the documenttion for regrid_consts + !! See the documentation for regrid_consts !! for the recognized values. real, intent(in) :: maxDepth !< The range of the grid values in some modes real, intent(in) :: rhoLight !< The minimum value of the grid in RHO mode diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 0c758fadaf..87019d46cf 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -373,7 +373,7 @@ function get_polynomial_coordinate( N, h, x_g, edge_values, ppoly_coefs, & real :: grad ! gradient during N-R iterations [A] integer :: i, k, iter ! loop indices integer :: k_found ! index of target cell - character(len=200) :: mesg + character(len=320) :: mesg logical :: use_2018_answers ! If true use older, less acccurate expressions. eps = NR_OFFSET diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 302ba0a714..04e151d5a7 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -38,7 +38,7 @@ module MOM_density_integrals !! required for calculating the finite-volume form pressure accelerations in a !! Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -78,13 +78,14 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] if (EOS_quadrature(EOS)) then call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) else call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif end subroutine int_density_dz @@ -93,8 +94,8 @@ end subroutine int_density_dz !> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, use_inaccurate_form) + EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dz_neglect, useMassWghtInterp, use_inaccurate_form, Z_0p) type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [degC] @@ -136,6 +137,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! interpolate T/S for top and bottom integrals. logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] + ! Local variables real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa] @@ -148,6 +151,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz ! The layer thickness [Z ~> m] + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A pressure-thickness below topography [Z ~> m] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m] real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2] @@ -173,6 +177,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p use_rho_ref = .true. if (present(use_inaccurate_form)) then if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form @@ -191,7 +196,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = z_t(i,j) - z_b(i,j) do n=1,5 T5(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz) enddo if (use_rho_ref) then if (rho_scale /= 1.0) then @@ -245,7 +250,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -300,7 +305,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) p5(n) = p5(n-1) + GxRho*0.25*dz @@ -336,7 +341,7 @@ end subroutine int_density_dz_generic_pcm subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, & intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, & - use_inaccurate_form) + use_inaccurate_form, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -379,6 +384,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & !! interpolate T/S for top and bottom integrals. logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -427,6 +433,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim] real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [degC] real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [ppt] + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] @@ -443,6 +450,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -473,7 +481,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & do i = Isq,Ieq+1 dz(i) = e(i,j,K) - e(i,j,K+1) do n=1,5 - p5(i*5+n) = -GxRho*(e(i,j,K) - 0.25*real(n-1)*dz(i)) + p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz(i)) ! Salinity and temperature points are linearly interpolated S5(i*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * S_b(i,j,k) T5(i*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * T_b(i,j,k) @@ -581,7 +589,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & S15(pos+1) = w_left*Stl + w_right*Str S15(pos+5) = w_left*Sbl + w_right*Sbr - p15(pos+1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i+1,j,K)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) ! Pressure do n=2,5 @@ -692,7 +700,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & S15(pos+1) = w_left*Stl + w_right*Str S15(pos+5) = w_left*Sbl + w_right*Sbr - p15(pos+1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i,j+1,K)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) ! Pressure do n=2,5 @@ -775,7 +783,7 @@ end subroutine int_density_dz_generic_plm !! are parabolic profiles subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, & - dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -816,6 +824,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & !! divided by the y grid spacing [R L2 T-2 ~> Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -854,6 +863,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: t6 ! PPM curvature coefficient for T [degC] real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] @@ -868,6 +878,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -900,7 +911,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & endif dz = e(i,j,K) - e(i,j,K+1) do n=1,5 - p5(n) = -GxRho*(e(i,j,K) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz) ! Salinity and temperature points are reconstructed with PPM S5(n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) ) T5(n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) ) @@ -978,7 +989,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i+1,j,K) - e(i+1,j,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i+1,j,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -1066,7 +1077,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i,j+1,K) - e(i,j+1,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i,j+1,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index ec7501c5f0..17729e586c 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -28,7 +28,7 @@ module MOM_interface_heights !! form for consistency with the calculation of the pressure gradient forces. !! Additionally, these height may be dilated for consistency with the !! corresponding time-average quantity from the barotropic calculation. -subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) +subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m, dZref) 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 @@ -37,14 +37,17 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !! thermodynamic variables. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta !< layer interface heights !! [Z ~> m] or [1/eta_to_m m]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total water - !! column mass per unit area (non-Boussinesq). This is used to dilate the layer. - !! thicknesses when calculating interfaceheights [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable + !! that gives the "correct" free surface height (Boussinesq) or total water + !! column mass per unit area (non-Boussinesq). This is used to dilate the layer + !! thicknesses when calculating interface heights [H ~> m or kg m-2]. + !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. integer, optional, intent(in) :: halo_size !< width of halo points on !! which to calculate eta. real, optional, intent(in) :: eta_to_m !< The conversion factor from - !! the units of eta to m; by default this is US%Z_to_m. + !! the units of eta to m; by default this is US%Z_to_m. + real, optional, intent(in) :: dZref !< The difference in the + !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. ! Local variables real :: p(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa] @@ -55,6 +58,8 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real :: I_gEarth ! The inverse of the gravitational acceleration times the ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names. + real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. + ! dZ_ref is 0 unless the optional argument dZref is present. integer i, j, k, isv, iev, jsv, jev, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -69,33 +74,35 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) H_to_eta = GV%H_to_Z * Z_to_eta H_to_rho_eta = GV%H_to_RZ * Z_to_eta I_gEarth = Z_to_eta / GV%g_Earth + dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref -!$OMP parallel default(shared) private(dilate,htot) -!$OMP do - do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -Z_to_eta*G%bathyT(i,j) ; enddo ; enddo + !$OMP parallel default(shared) private(dilate,htot) + !$OMP do + do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -Z_to_eta*(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo if (GV%Boussinesq) then -!$OMP do + !$OMP do do j=jsv,jev ; do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)*H_to_eta enddo ; enddo ; enddo if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! that is used for the dynamics. -!$OMP do + !$OMP do do j=jsv,jev do i=isv,iev dilate(i) = (eta_bt(i,j)*H_to_eta + Z_to_eta*G%bathyT(i,j)) / & - (eta(i,j,1) + Z_to_eta*G%bathyT(i,j)) + (eta(i,j,1) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) enddo do k=1,nz ; do i=isv,iev - eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*G%bathyT(i,j)) - Z_to_eta*G%bathyT(i,j) + eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo ; enddo enddo endif else if (associated(tv%eqn_of_state)) then -!$OMP do + !$OMP do do j=jsv,jev if (associated(tv%p_surf)) then do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo @@ -106,19 +113,19 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) p(i,j,K+1) = p(i,j,K) + GV%g_Earth*GV%H_to_RZ*h(i,j,k) enddo ; enddo enddo -!$OMP do + !$OMP do do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & 0.0, G%HI, tv%eqn_of_state, US, dz_geo(:,:,k), halo_size=halo) enddo -!$OMP do + !$OMP do do j=jsv,jev do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + I_gEarth * dz_geo(i,j,k) enddo ; enddo enddo else -!$OMP do + !$OMP do do j=jsv,jev ; do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + H_to_rho_eta*h(i,j,k) / GV%Rlay(k) enddo ; enddo ; enddo @@ -126,18 +133,19 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! from the time-averaged barotropic solution. -!$OMP do + !$OMP do do j=jsv,jev do i=isv,iev ; htot(i) = GV%H_subroundoff ; enddo do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo do k=1,nz ; do i=isv,iev - eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*G%bathyT(i,j)) - Z_to_eta*G%bathyT(i,j) + eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo ; enddo enddo endif endif -!$OMP end parallel + !$OMP end parallel end subroutine find_eta_3d @@ -145,7 +153,7 @@ end subroutine find_eta_3d !! with the calculation of the pressure gradient forces. Additionally, the sea !! surface height may be adjusted for consistency with the corresponding !! time-average quantity from the barotropic calculation. -subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) +subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m, dZref) 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 @@ -155,12 +163,16 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta !< free surface height relative to !! mean sea level (z=0) often [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total - !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2]. + !! variable that gives the "correct" free surface height (Boussinesq) or total + !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2]. + !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. integer, optional, intent(in) :: halo_size !< width of halo points on !! which to calculate eta. real, optional, intent(in) :: eta_to_m !< The conversion factor from - !! the units of eta to m; by default this is US%Z_to_m. + !! the units of eta to m; by default this is US%Z_to_m. + real, optional, intent(in) :: dZref !< The difference in the + !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. + ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & p ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa] @@ -170,6 +182,8 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real :: I_gEarth ! The inverse of the gravitational acceleration times the ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names. + real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. + ! dZ_ref is 0 unless the optional argument dZref is present. integer i, j, k, is, ie, js, je, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -180,26 +194,27 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) H_to_eta = GV%H_to_Z * Z_to_eta H_to_rho_eta = GV%H_to_RZ * Z_to_eta I_gEarth = Z_to_eta / GV%g_Earth + dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref -!$OMP parallel default(shared) private(htot) -!$OMP do - do j=js,je ; do i=is,ie ; eta(i,j) = -Z_to_eta*G%bathyT(i,j) ; enddo ; enddo + !$OMP parallel default(shared) private(htot) + !$OMP do + do j=js,je ; do i=is,ie ; eta(i,j) = -Z_to_eta*(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo if (GV%Boussinesq) then if (present(eta_bt)) then -!$OMP do + !$OMP do do j=js,je ; do i=is,ie - eta(i,j) = H_to_eta*eta_bt(i,j) + eta(i,j) = H_to_eta*eta_bt(i,j) - Z_to_eta*dZ_ref enddo ; enddo else -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + h(i,j,k)*H_to_eta enddo ; enddo ; enddo endif else if (associated(tv%eqn_of_state)) then -!$OMP do + !$OMP do do j=js,je if (associated(tv%p_surf)) then do i=is,ie ; p(i,j,1) = tv%p_surf(i,j) ; enddo @@ -211,17 +226,17 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) p(i,j,k+1) = p(i,j,k) + GV%g_Earth*GV%H_to_RZ*h(i,j,k) enddo ; enddo enddo -!$OMP do + !$OMP do do k = 1, nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, & G%HI, tv%eqn_of_state, US, dz_geo(:,:,k), halo_size=halo) enddo -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + I_gEarth * dz_geo(i,j,k) enddo ; enddo ; enddo else -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + H_to_rho_eta*h(i,j,k) / GV%Rlay(k) enddo ; enddo ; enddo @@ -229,18 +244,18 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) if (present(eta_bt)) then ! Dilate the water column to agree with the time-averaged column ! mass from the barotropic solution. -!$OMP do + !$OMP do do j=js,je do i=is,ie ; htot(i) = GV%H_subroundoff ; enddo do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo do i=is,ie - eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + Z_to_eta*G%bathyT(i,j)) - & - Z_to_eta*G%bathyT(i,j) + eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo enddo endif endif -!$OMP end parallel + !$OMP end parallel end subroutine find_eta_2d diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 61e20d14a6..318d10008c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3802,29 +3802,31 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ishift=0;jshift=0 if (segment%is_E_or_W) then allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed)) - normal_trans_bt(:,:)=0.0 + normal_trans_bt(:,:) = 0.0 if (segment%direction == OBC_DIRECTION_W) ishift=1 I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed - segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j)) - segment%Htot(I,j)=0.0 + segment%Htot(I,j) = 0.0 do k=1,GV%ke segment%h(I,j,k) = h(i+ishift,j,k) - segment%Htot(I,j)=segment%Htot(I,j)+segment%h(I,j,k) + segment%Htot(I,j) = segment%Htot(I,j) + segment%h(I,j,k) enddo + segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j)) + !### This should be: segment%Cg(I,j) = sqrt(GV%g_prime(1)*segment%Htot(I,j)*GV%H_to_Z) enddo else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S) allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB)) - normal_trans_bt(:,:)=0.0 + normal_trans_bt(:,:) = 0.0 if (segment%direction == OBC_DIRECTION_S) jshift=1 J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied - segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift)) - segment%Htot(i,J)=0.0 + segment%Htot(i,J) = 0.0 do k=1,GV%ke segment%h(i,J,k) = h(i,j+jshift,k) - segment%Htot(i,J)=segment%Htot(i,J)+segment%h(i,J,k) + segment%Htot(i,J) = segment%Htot(i,J) + segment%h(i,J,k) enddo + segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift)) + !### This should be: segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z) enddo endif @@ -4715,7 +4717,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth + real :: min_depth ! The minimum depth for ocean points [Z ~> m] integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list @@ -4730,7 +4732,6 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) allocate(color(G%isd:G%ied, G%jsd:G%jed)) ; color = 0 allocate(color2(G%isd:G%ied, G%jsd:G%jed)) ; color2 = 0 - ! Paint a frame around the outside. do j=G%jsd,G%jed color(G%isd,j) = cedge diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 035386f92d..d363b185f8 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -444,7 +444,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ hw = 0.5*(Hc(k-1)+Hc(k)) gp = gprime(K) if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then - if ( ((G%bathyT(i,j)-sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. & + !### Change to: if ( ((htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) .or. & ) ) + if ( ((G%bathyT(i,j) - sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. & ((l_mono_N2_depth >= 0.) .and. (sum_hc > l_mono_N2_depth))) .and. & (L2_to_Z2*gp > N2min*hw) ) then ! Filters out regions where N2 increases with depth but only in a lower fraction diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 6e74c3ffa3..23f22d8a24 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1253,7 +1253,7 @@ end subroutine analytic_int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -1292,6 +1292,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] + ! Local variables real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the ! desired units [R m3 kg-1 ~> 1] @@ -1322,11 +1324,11 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p=Z_0p) else call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp) + dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif case default call MOM_error(FATAL, "No analytic integration option is available with this EOS!") diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index e5cc9555b7..730687fbf6 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -408,7 +408,7 @@ end subroutine calculate_compress_wright !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -451,6 +451,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure !! into Pa [Pa T2 R-1 L-2 ~> 1]. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! Local variables real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d @@ -461,7 +462,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: g_Earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2] real :: I_Rho ! The inverse of the Boussinesq density [m3 kg-1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] - real :: p_ave, I_al0, I_Lzz + real :: p_ave ! The layer averaged pressure [Pa] + real :: I_al0, I_Lzz real :: dz ! The layer thickness [Z ~> m]. real :: hWght ! A pressure-thickness below topography [Z ~> m]. real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m]. @@ -470,10 +472,11 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. - real :: intz(5) ! The integrals of density with height at the - ! 5 sub-column locations [R L2 T-2 ~> Pa]. + real :: intz(5) ! The gravitational acceleration times the integrals of density + ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]. real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by ! pres_scale [R L2 T-2 Pa-1 ~> 1] or [1]. + real :: z0pres ! The height at which the pressure is zero [Z ~> m] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. @@ -499,6 +502,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & else rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0 endif + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -517,7 +521,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dz = z_t(i,j) - z_b(i,j) - p_ave = -0.5*GxRho*(z_t(i,j)+z_b(i,j)) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -561,8 +565,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j) dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) - p_ave = -0.5*GxRho*(wt_L*(z_t(i,j)+z_b(i,j)) + & - wt_R*(z_t(i+1,j)+z_b(i+1,j))) + p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i+1,j)+z_b(i+1,j))) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -603,8 +606,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i,j+1) dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) - p_ave = -0.5*GxRho*(wt_L*(z_t(i,j)+z_b(i,j)) + & - wt_R*(z_t(i,j+1)+z_b(i,j+1))) + p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i,j+1)+z_b(i,j+1))) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 2a9a381caa..89a59374a7 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -294,9 +294,9 @@ subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units) if (m_in_new_units == 1.0) return if (m_in_new_units < 0.0) & - call MOM_error(FATAL, "rescale_grid_bathymetry: Negative depth units are not permitted.") + call MOM_error(FATAL, "rescale_dyn_horgrid_bathymetry: Negative depth units are not permitted.") if (m_in_new_units == 0.0) & - call MOM_error(FATAL, "rescale_grid_bathymetry: Zero depth units are not permitted.") + call MOM_error(FATAL, "rescale_dyn_horgrid_bathymetry: Zero depth units are not permitted.") rescale = 1.0 / m_in_new_units do j=jsd,jed ; do i=isd,ied diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f8c772348..62b06a0b57 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -11,7 +11,7 @@ module MOM_state_initialization use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, NOTE, WARNING, is_root_pe +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_file_parser, only : log_version @@ -655,7 +655,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi !! only read parameters without changing h. ! Local variables - real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units. + real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units [Z ~> m]. integer :: inconsistent = 0 logical :: correct_thickness logical :: just_read ! If true, just read parameters but set nothing. @@ -696,7 +696,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi call MOM_read_data(filename, "eta", eta(:,:,:), G%Domain, scale=US%m_to_Z) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, eta, h) + call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=0.0) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then @@ -732,25 +732,30 @@ end subroutine initialize_thickness_from_file !! is dilated (expanded) to fill the void. !! @remark{There is a (hard-wired) "tolerance" parameter such that the !! criteria for adjustment must equal or exceed 10cm.} -subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h) +subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta) 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)+1), intent(inout) :: eta !< Interface heights [Z ~> m]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, optional, intent(in) :: dZ_ref_eta !< The difference between the + !! reference heights for bathyT and + !! eta [Z ~> m], 0 by default. ! Local variables integer :: i, j, k, is, ie, js, je, nz, contractions, dilations real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m] - real :: hTmp, eTmp, dilate + real :: dilate ! A factor by which the column is dilated [nondim] + real :: dZ_ref ! The difference in the reference heights for G%bathyT and eta [Z ~> m] character(len=100) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke hTolerance = 0.1*US%m_to_Z + dZ_ref = 0.0 ; if (present(dZ_ref_eta)) dZ_ref = dZ_ref_eta contractions = 0 do j=js,je ; do i=is,ie - if (-eta(i,j,nz+1) > G%bathyT(i,j) + hTolerance) then - eta(i,j,nz+1) = -G%bathyT(i,j) + if (-eta(i,j,nz+1) > (G%bathyT(i,j) + dZ_ref) + hTolerance) then + eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) contractions = contractions + 1 endif enddo ; enddo @@ -779,12 +784,12 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h) ! The whole column is dilated to accommodate deeper topography than ! the bathymetry would indicate. ! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. & - if (-eta(i,j,nz+1) < G%bathyT(i,j) - hTolerance) then + if (-eta(i,j,nz+1) < (G%bathyT(i,j) + dZ_ref) - hTolerance) then dilations = dilations + 1 if (eta(i,j,1) <= eta(i,j,nz+1)) then - do k=1,nz ; h(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo + do k=1,nz ; h(i,j,k) = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / real(nz) ; enddo else - dilate = (eta(i,j,1) + G%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1)) + dilate = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / (eta(i,j,1) - eta(i,j,nz+1)) do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo endif do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + h(i,j,k) ; enddo @@ -1746,11 +1751,11 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic !! variables. real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u !< The zonal velocity that is being - !! initialized [L T-1 ~> m s-1] + intent(in) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v !< The meridional velocity that is being - !! initialized [L T-1 ~> m s-1] + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control !! structure for this module (in layered mode). @@ -2182,9 +2187,9 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p ! calculate increments if input are full fields if (oda_inc) then ! input are increments - if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") + if (is_root_pe()) call MOM_mesg("incupd using increments fields ") else ! inputs are full fields - if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") + if (is_root_pe()) call MOM_mesg("incupd using full fields ") call calc_oda_increments(h, tv, u, v, G, GV, US, oda_incupd_CSp) if (save_inc) then call output_oda_incupd_inc(Time, G, GV, param_file, oda_incupd_CSp, US) @@ -2319,6 +2324,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor + ! relative to the surface [Z ~> m]. integer, dimension(SZI_(G),SZJ_(G)) :: nlevs real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa]. @@ -2513,6 +2520,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call pass_var(mask_z,G%Domain) call pass_var(rho_z,G%Domain) + do j=js,je ; do i=is,ie + Z_bottom(i,j) = -G%bathyT(i,j) + enddo ; enddo + ! Done with horizontal interpolation. ! Now remap to model coordinates if (useALEremapping) then @@ -2530,11 +2541,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param tmp_mask_in(i,j,1:kd) = mask_z(i,j,:) do k = 1, nkd if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then - zBottomOfCell = max( z_edges_in(k+1), -G%bathyT(i,j) ) + zBottomOfCell = max( z_edges_in(k+1), Z_bottom(i,j)) tmpT1dIn(i,j,k) = temp_z(i,j,k) tmpS1dIn(i,j,k) = salt_z(i,j,k) elseif (k>1) then - zBottomOfCell = -G%bathyT(i,j) + zBottomOfCell = Z_bottom(i,j) tmpT1dIn(i,j,k) = tmpT1dIn(i,j,k-1) tmpS1dIn(i,j,k) = tmpS1dIn(i,j,k-1) else ! This next block should only ever be reached over land @@ -2544,7 +2555,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param h1(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo - h1(i,j,kd) = h1(i,j,kd) + GV%Z_to_H * max(0., zTopOfCell + G%bathyT(i,j) ) + h1(i,j,kd) = h1(i,j,kd) + GV%Z_to_H * max(0., zTopOfCell - Z_bottom(i,j) ) ! The max here is in case the data data is shallower than model endif ! mask2dT enddo ; enddo @@ -2567,7 +2578,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! Build the target grid combining hTarget and topography zTopOfCell = 0. ; zBottomOfCell = 0. do k = 1, nz - zBottomOfCell = max( zTopOfCell - hTarget(k), -G%bathyT(i,j) ) + zBottomOfCell = max( zTopOfCell - hTarget(k), Z_bottom(i,j)) h(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo @@ -2624,11 +2635,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param nkml = 0 ; if (separate_mixed_layer) nkml = GV%nkml - call find_interfaces(rho_z, z_in, kd, Rb, G%bathyT, zi, G, GV, US, nlevs, nkml, & + call find_interfaces(rho_z, z_in, kd, Rb, Z_bottom, zi, G, GV, US, nlevs, nkml, & Hmix_depth, eps_z, eps_rho, density_extrap_bug) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, zi, h) + call adjustEtaToFitBathymetry(G, GV, US, zi, h, dZ_ref_eta=0.0) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (zi(i,j,K) < (zi(i,j,K+1) + GV%Angstrom_Z)) then @@ -2640,7 +2651,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param enddo ; enddo ; enddo inconsistent=0 do j=js,je ; do i=is,ie - if (abs(zi(i,j,nz+1) + G%bathyT(i,j)) > 1.0*US%m_to_Z) & + if (abs(zi(i,j,nz+1) - Z_bottom(i,j)) > 1.0*US%m_to_Z) & inconsistent = inconsistent + 1 enddo ; enddo call sum_across_PEs(inconsistent) @@ -2710,7 +2721,7 @@ end subroutine MOM_temp_salt_initialize_from_Z !> Find interface positions corresponding to interpolated depths in a density profile -subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, nkml, hml, & +subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, nkml, hml, & eps_z, eps_rho, density_extrap_bug) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure @@ -2720,7 +2731,8 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, n real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m]. real, dimension(SZK_(GV)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: depth !< ocean depth [Z ~> m]. + intent(in) :: Z_bot !< The (usually negative) height of the seafloor + !! relative to the surface [Z ~> m]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(out) :: zi !< The returned interface heights [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -2808,15 +2820,15 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, n ! Linearly interpolate to find the depth, zi_, where Rb would be found. slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho) zi_(K) = -1.0*(zin(k_int) + slope*(Rb(K)-rho_(k_int))) - zi_(K) = min(max(zi_(K), -depth(i,j)), -1.0*hml) + zi_(K) = min(max(zi_(K), Z_bot(i,j)), -1.0*hml) enddo - zi_(nz+1) = -depth(i,j) + zi_(nz+1) = Z_bot(i,j) if (nkml > 0) then ; do K=2,nkml+1 - zi_(K) = max(hml*((1.0-real(K))/real(nkml)), -depth(i,j)) + zi_(K) = max(hml*((1.0-real(K))/real(nkml)), Z_bot(i,j)) enddo ; endif do K=nz,max(nkml+2,2),-1 if (zi_(K) < zi_(K+1) + eps_Z) zi_(K) = zi_(K+1) + eps_Z - if (zi_(K) > -1.0*hml) zi_(K) = max(-1.0*hml, -depth(i,j)) + if (zi_(K) > -1.0*hml) zi_(K) = max(-1.0*hml, Z_bot(i,j)) enddo do K=1,nz+1 @@ -2864,7 +2876,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) S_t(k) = 35. - (0. * I_z_scale)*e(k) S(k) = 35. + (0. * I_z_scale)*z(k) S_b(k) = 35. - (0. * I_z_scale)*e(k+1) - call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*US%m_to_Z*z(k), & + call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*z(k), & rho(k), tv%eqn_of_state) P_tot = P_tot + GV%g_Earth * rho(k) * GV%H_to_Z*h(k) enddo diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 762b2edaea..9441ab7107 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -637,7 +637,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h end subroutine step_forward_MEKE -!> Calculates the equilibrium solutino where the source depends only on MEKE diffusivity +!> Calculates the equilibrium solution where the source depends only on MEKE diffusivity !! and there is no lateral diffusion of MEKE. !! Results is in MEKE%MEKE. subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass) @@ -667,6 +667,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1] real :: FatH ! Coriolis parameter at h points; to compute topographic beta [T-1 ~> s-1] real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1] + real :: dZ_neglect ! A negligible change in height [Z ~> m] integer :: i, j, is, ie, js, je, n1, n2 real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2]. logical :: useSecant, debugIteration @@ -680,6 +681,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m Ubg2 = CS%MEKE_Uscale**2 cd2 = CS%cdrag**2 tolerance = 1.0e-12*US%m_s_to_L_T**2 + dZ_neglect = GV%H_to_Z*GV%H_subroundoff !$OMP do do j=js,je ; do i=is,ie @@ -701,14 +703,14 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) @@ -853,9 +855,11 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & real :: SN ! The local Eady growth rate [T-1 ~> s-1] real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1] real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1] + real :: dZ_neglect ! A negligible change in height [Z ~> m] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + dZ_neglect = GV%H_to_Z*GV%H_subroundoff !$OMP do do j=js,je ; do i=is,ie @@ -878,14 +882,14 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c588a1faa4..b4f857dec4 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -496,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + & + grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + & (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) enddo ; enddo @@ -1389,7 +1389,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (grad_vel_mag_bt_h(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) else GME_coeff = 0.0 @@ -1405,8 +1405,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - if (grad_vel_mag_bt_q(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + if (grad_vel_mag_bt_q(I,J)>0) then + !### This expression is not rotationally invariant - bathyT is to the SW of q points, + ! and it needs parentheses in the sum of the 4 diffusivities. + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) else GME_coeff = 0.0 diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index da62ffc6b7..3b3d72576c 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -481,14 +481,19 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp do j=js,je ; do I=is-1,ie hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect) if (hu(I,j) /= 0.0) hu(I,j) = 1.0 + !### The same result would be accomplished with the following without a division: + ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1)) enddo ; enddo do J=js-1,je ; do i=is,ie hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) if (hv(i,J) /= 0.0) hv(i,J) = 1.0 + !### The same result would be accomplished with the following without a division: + ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1)) enddo ; enddo ! diagnose diffusivity at T-point + !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent. do j=js,je ; do i=is,ie Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j)) & +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) & @@ -505,7 +510,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index e122452368..419b012387 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -47,7 +47,7 @@ module MOM_ALE_sponge module procedure set_up_ALE_sponge_vel_field_varying end interface -!> Ddetermine the number of points which are within sponges in this computational domain. +!> Determine the number of points which are within sponges in this computational domain. !! !! Only points that have positive values of Iresttime and which mask2dT indicates are ocean !! points are included in the sponges. It also stores the target interface heights. @@ -90,31 +90,19 @@ module MOM_ALE_sponge !> ALE sponge control structure type, public :: ALE_sponge_CS ; private integer :: nz !< The total number of layers. - integer :: nz_data !< The total number of arbritary layers (used by older code). - integer :: isc !< The starting i-index of the computational domain at h. - integer :: iec !< The ending i-index of the computational domain at h. - integer :: jsc !< The starting j-index of the computational domain at h. - integer :: jec !< The ending j-index of the computational domain at h. - integer :: IscB !< The starting I-index of the computational domain at u/v. - integer :: IecB !< The ending I-index of the computational domain at u/v. - integer :: JscB !< The starting J-index of the computational domain at u/v. - integer :: JecB !< The ending J-index of the computational domain at h. - integer :: isd !< The starting i-index of the data domain at h. - integer :: ied !< The ending i-index of the data domain at h. - integer :: jsd !< The starting j-index of the data domain at h. - integer :: jed !< The ending j-index of the data domain at h. + integer :: nz_data !< The total number of arbitrary layers (used by older code). integer :: num_col !< The number of sponge tracer points within the computational domain. integer :: num_col_u !< The number of sponge u-points within the computational domain. integer :: num_col_v !< The number of sponge v-points within the computational domain. integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_sponge_field logical :: sponge_uv !< Control whether u and v are included in sponge - integer, pointer :: col_i(:) => NULL() !< Array of the i-indicies of each tracer columns being damped. - integer, pointer :: col_j(:) => NULL() !< Array of the j-indicies of each tracer columns being damped. - integer, pointer :: col_i_u(:) => NULL() !< Array of the i-indicies of each u-columns being damped. - integer, pointer :: col_j_u(:) => NULL() !< Array of the j-indicies of each u-columns being damped. - integer, pointer :: col_i_v(:) => NULL() !< Array of the i-indicies of each v-columns being damped. - integer, pointer :: col_j_v(:) => NULL() !< Array of the j-indicies of each v-columns being damped. + integer, pointer :: col_i(:) => NULL() !< Array of the i-indices of each tracer column being damped. + integer, pointer :: col_j(:) => NULL() !< Array of the j-indices of each tracer column being damped. + integer, pointer :: col_i_u(:) => NULL() !< Array of the i-indices of each u-column being damped. + integer, pointer :: col_j_u(:) => NULL() !< Array of the j-indices of each u-column being damped. + integer, pointer :: col_i_v(:) => NULL() !< Array of the i-indices of each v-column being damped. + integer, pointer :: col_j_v(:) => NULL() !< Array of the j-indices of each v-column being damped. real, pointer :: Iresttime_col(:) => NULL() !< The inverse restoring time of each tracer column [T-1 ~> s-1]. real, pointer :: Iresttime_col_u(:) => NULL() !< The inverse restoring time of each u-column [T-1 ~> s-1]. @@ -124,8 +112,8 @@ module MOM_ALE_sponge type(p2d) :: Ref_val(MAX_FIELDS_) !< The values to which the fields are damped. type(p2d) :: Ref_val_u !< The values to which the u-velocities are damped. type(p2d) :: Ref_val_v !< The values to which the v-velocities are damped. - type(p3d) :: var_u !< Pointer to the u velocities. that are being damped. - type(p3d) :: var_v !< Pointer to the v velocities. that are being damped. + type(p3d) :: var_u !< Pointer to the u velocities that are being damped. + type(p3d) :: var_v !< Pointer to the v velocities that are being damped. type(p2d) :: Ref_h !< Grid on which reference data is provided (older code). type(p2d) :: Ref_hu !< u-point grid on which reference data is provided (older code). type(p2d) :: Ref_hv !< v-point grid on which reference data is provided (older code). @@ -137,7 +125,7 @@ module MOM_ALE_sponge logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that !! recover the answers for remapping from the end of 2018. !! Otherwise, use more robust forms of the same expressions. - logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding + logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizontal regridding !! that recovers the answers from the end of 2018. Otherwise, use !! rotationally symmetric forms of the same expressions. @@ -241,9 +229,6 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, CS%time_varying_sponges = .false. CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%iscB = G%iscB ; CS%iecB = G%iecB; CS%jscB = G%jscB ; CS%jecB = G%jecB ! number of columns to be restored CS%num_col = 0 ; CS%fldno = 0 @@ -265,7 +250,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, col = col +1 endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data CS%nz_data = nz_data allocate(CS%Ref_h%p(CS%nz_data,CS%num_col)) do col=1,CS%num_col ; do K=1,CS%nz_data @@ -295,11 +280,11 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, if (present(Iresttime_u_in)) then Iresttime_u(:,:) = Iresttime_u_in(:,:) else - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) enddo ; enddo endif - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 enddo ; enddo @@ -312,7 +297,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! Store the column indices and restoring rates in the CS structure col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = I ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(I,j) @@ -320,7 +305,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) do col=1,CS%num_col_u I = CS%col_i_u(col) ; j = CS%col_j_u(col) @@ -339,11 +324,11 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, if (present(Iresttime_v_in)) then Iresttime_v(:,:) = Iresttime_v_in(:,:) else - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) enddo ; enddo endif - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 enddo ; enddo @@ -356,7 +341,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! pass indices, restoring time to the CS structure col = 1 - do J=CS%jscB,CS%jecB ; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB ; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) @@ -364,7 +349,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data allocate(CS%Ref_hv%p(CS%nz_data,CS%num_col_v)) do col=1,CS%num_col_v i = CS%col_i_v(col) ; J = CS%col_j_v(col) @@ -430,7 +415,7 @@ subroutine get_ALE_sponge_thicknesses(G, data_h, sponge_mask, CS) end subroutine get_ALE_sponge_thicknesses -!> This subroutine determines the number of points which are to be restoref in the computational +!> This subroutine determines the number of points which are to be restored in the computational !! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean !! points are included in the sponges. subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Iresttime_u_in, Iresttime_v_in) @@ -510,9 +495,6 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest CS%time_varying_sponges = .true. CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%iscB = G%iscB ; CS%iecB = G%iecB; CS%jscB = G%jscB ; CS%jecB = G%jecB ! number of columns to be restored CS%num_col = 0 ; CS%fldno = 0 @@ -551,12 +533,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest if (present(Iresttime_u_in)) then Iresttime_u(:,:) = Iresttime_u_in(:,:) else - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) enddo ; enddo endif CS%num_col_u = 0 ; - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 enddo ; enddo @@ -566,14 +548,14 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u = 0 ! pass indices, restoring time to the CS structure col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) col = col + 1 endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data endif total_sponge_cols_u = CS%num_col_u call sum_across_PEs(total_sponge_cols_u) @@ -583,12 +565,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest if (present(Iresttime_v_in)) then Iresttime_v(:,:) = Iresttime_v_in(:,:) else - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) enddo ; enddo endif CS%num_col_v = 0 ; - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 enddo ; enddo @@ -598,7 +580,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest allocate(CS%col_j_v(CS%num_col_v)) ; CS%col_j_v = 0 ! pass indices, restoring time to the CS structure col = 1 - do J=CS%jscB,CS%jecB ; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB ; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) @@ -630,16 +612,16 @@ subroutine init_ALE_sponge_diags(Time, G, diag, CS, US) CS%id_sp_tendency(1) = -1 CS%id_sp_tendency(1) = register_diag_field('ocean_model', 'sp_tendency_temp', diag%axesTL, Time, & - 'Time tendency due to temperature restoring', 'degC s-1',conversion=US%s_to_T) + 'Time tendency due to temperature restoring', 'degC s-1', conversion=US%s_to_T) CS%id_sp_tendency(2) = -1 CS%id_sp_tendency(2) = register_diag_field('ocean_model', 'sp_tendency_salt', diag%axesTL, Time, & - 'Time tendency due to salinity restoring', 'g kg-1 s-1',conversion=US%s_to_T) + 'Time tendency due to salinity restoring', 'g kg-1 s-1', conversion=US%s_to_T) CS%id_sp_u_tendency = -1 CS%id_sp_u_tendency = register_diag_field('ocean_model', 'sp_tendency_u', diag%axesCuL, Time, & - 'Zonal acceleration due to sponges', 'm s-2',conversion=US%L_T2_to_m_s2) + 'Zonal acceleration due to sponges', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_sp_v_tendency = -1 CS%id_sp_v_tendency = register_diag_field('ocean_model', 'sp_tendency_v', diag%axesCvL, Time, & - 'Meridional acceleration due to sponges', 'm s-2',conversion=US%L_T2_to_m_s2) + 'Meridional acceleration due to sponges', 'm s-2', conversion=US%L_T2_to_m_s2) end subroutine init_ALE_sponge_diags @@ -718,19 +700,18 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed CS%fldno = CS%fldno + 1 if (CS%fldno > MAX_FIELDS_) then - write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & - &the number of fields to be damped in the call to & - &initialize_ALE_sponge." )') CS%fldno + write(mesg, '("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease "//& + "the number of fields to be damped in the call to initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif - ! get a unique time interp id for this field. If sponge data is ongrid, then setup + ! get a unique time interp id for this field. If sponge data is on-grid, then setup ! to only read on the computational domain if (CS%spongeDataOngrid) then CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname, MOM_domain=G%Domain) else CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname) endif - fld_sz(1:4)=-1 + fld_sz(1:4) = -1 call get_external_field_info(CS%Ref_val(CS%fldno)%id, size=fld_sz) nz_data = fld_sz(3) CS%Ref_val(CS%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid @@ -748,17 +729,19 @@ end subroutine set_up_ALE_sponge_field_varying !> This subroutine stores the reference profile at u and v points for the variable !! whose address is given by u_ptr and v_ptr. subroutine set_up_ALE_sponge_vel_field_fixed(u_val, v_val, G, GV, u_ptr, v_ptr, CS) - type(ocean_grid_type), intent(in) :: G !< Grid structure (in). + type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). + type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u_val !< u field to be used in the sponge, it has arbritary number of layers but - !! not to exceed the total number of model layers + intent(in) :: u_val !< u field to be used in the sponge [L T-1 ~> m s-1], + !! it is provided on its own vertical grid that may + !! have fewer layers than the model itself, but not more. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v_val !< v field to be used in the sponge, it has arbritary number of layers but - !! not to exceed the number of model layers - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u pointer to the field to be damped - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v pointer to the field to be damped + intent(in) :: v_val !< v field to be used in the sponge [L T-1 ~> m s-1], + !! it is provided on its own vertical grid that may + !! have fewer layers than the model itself, but not more. + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] integer :: j, k, col, fld_sz(4) character(len=256) :: mesg ! String for error messages @@ -794,15 +777,16 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename character(len=*), intent(in) :: filename_v !< File name for v field character(len=*), intent(in) :: fieldname_v !< Name of v variable in file type(time_type), intent(in) :: Time !< Model time - type(ocean_grid_type), intent(in) :: G !< Ocean grid (in) - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u pointer to the field to be damped (in). - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v pointer to the field to be damped (in). + type(ocean_grid_type), intent(in) :: G !< Ocean grid (in) + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] + ! Local variables - real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge. - real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge. + real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge [L T-1 ~> m s-1]. + real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge [L T-1 ~> m s-1]. real, allocatable, dimension(:), target :: z_in, z_edges_in real :: missing_value @@ -892,10 +876,13 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real, dimension(:), allocatable :: tmpT1d integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz, nz_data integer :: col, total_sponge_cols - real, allocatable, dimension(:), target :: z_in, z_edges_in - real :: missing_value, Idt - real :: h_neglect, h_neglect_edge - real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m]. + real, allocatable, dimension(:), target :: z_in ! The depths (positive downward) in the input file [Z ~> m] + real, allocatable, dimension(:), target :: z_edges_in ! The depths (positive downward) of the + ! edges in the input file [Z ~> m] + real :: missing_value + real :: Idt ! The inverse of the timestep [T-1 ~> s-1] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] + real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m]. integer :: nPoints is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1020,7 +1007,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) call pass_var(sp_val, G%Domain) call pass_var(mask_z, G%Domain) - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec; do I=G%iscB,G%iecB sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data)) mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) enddo ; enddo @@ -1068,7 +1055,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) answers_2018=CS%hor_regrid_answers_2018) call pass_var(sp_val, G%Domain) call pass_var(mask_z, G%Domain) - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data)) mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) enddo ; enddo @@ -1107,7 +1094,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%Ref_val_u%nz_data allocate(tmp_val2(nz_data)) if (CS%id_sp_u_tendency > 0) then - allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz));tmp_u(:,:,:)=0.0 + allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz)) ; tmp_u(:,:,:)=0.0 endif ! u points do c=1,CS%num_col_u @@ -1137,7 +1124,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) endif ! v points if (CS%id_sp_v_tendency > 0) then - allocate(tmp_v(G%isd:G%ied,G%jsdB:G%jedB,nz));tmp_v(:,:,:)=0.0 + allocate(tmp_v(G%isd:G%ied,G%jsdB:G%jedB,nz)) ; tmp_v(:,:,:)=0.0 endif nz_data = CS%Ref_val_v%nz_data allocate(tmp_val2(nz_data)) @@ -1170,9 +1157,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) deallocate(tmp_val2) endif - - - end subroutine apply_ALE_sponge !> Rotate the ALE sponge fields from the input to the model index map.