From 060c41290bcb4a3312bdc2bab419fb59dd473c3f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 31 Aug 2021 09:24:34 -0400 Subject: [PATCH] (*)Offset G%bathyT by -G%Z_ref Added code to adjust the value of G%bathyT in the ocean_grid_type by -G%Z_ref, with changes scattered throughout the code to compensate. This does not impact the value of bathyT in the dyn_hor_grid. Using a non-zero value of REFERENCE_HEIGHT leads to mathematically equivalent solutions that differ at round-off but are demonstrably similar for values of REFERENCE_HEIGHT ranging from -1e4 m to 1e4 m, and in cases that do not use the less exact ..._2018_ANSWERS algorithms the solutions are qualitatively equivalent for values ranging from -1e8 m to 1e8 m. (For these more extreme values, the 64-bit roundoff in the free surface height calculation starts to become physically significant at of order 0.1 mm.) This new capability is useful for developing wetting and drying and for identifying algorithmic shortcomings, but because answers are not bitwise identical for various reference heights, it is not such a good candidate for fully automated testing. By default, all answers are bitwise identical and all output files are the same. --- src/ALE/MOM_ALE.F90 | 4 +- src/ALE/MOM_regridding.F90 | 12 ++-- src/ALE/coord_adapt.F90 | 2 +- src/core/MOM.F90 | 18 ++--- src/core/MOM_PressureForce_FV.F90 | 42 ++++++------ src/core/MOM_PressureForce_Montgomery.F90 | 19 +++--- src/core/MOM_barotropic.F90 | 57 +++++++++------- src/core/MOM_grid.F90 | 11 +-- src/core/MOM_open_boundary.F90 | 1 + src/core/MOM_transcribe_grid.F90 | 26 +++---- src/diagnostics/MOM_PointAccel.F90 | 12 ++-- src/diagnostics/MOM_diagnostics.F90 | 12 ++-- src/diagnostics/MOM_sum_output.F90 | 12 ++-- src/diagnostics/MOM_wave_speed.F90 | 2 +- src/framework/MOM_diag_remap.F90 | 10 +-- src/framework/MOM_horizontal_regridding.F90 | 11 +-- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 68 +++++++++---------- .../MOM_state_initialization.F90 | 10 +-- .../MOM_tracer_initialization_from_Z.F90 | 2 +- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- .../lateral/MOM_internal_tides.F90 | 2 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 26 +++---- .../vertical/MOM_ALE_sponge.F90 | 6 +- .../vertical/MOM_diabatic_driver.F90 | 2 +- .../vertical/MOM_internal_tide_input.F90 | 4 +- .../vertical/MOM_set_viscosity.F90 | 18 ++--- src/parameterizations/vertical/MOM_sponge.F90 | 6 +- .../vertical/MOM_tidal_mixing.F90 | 12 ++-- src/tracer/MOM_tracer_Z_init.F90 | 8 +-- 30 files changed, 219 insertions(+), 200 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 63f8193b33..93696d3879 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -351,7 +351,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) if (CS%id_T_preale > 0) call post_data(CS%id_T_preale, tv%T, CS%diag) if (CS%id_S_preale > 0) call post_data(CS%id_S_preale, tv%S, CS%diag) if (CS%id_e_preale > 0) then - call find_eta(h, tv, G, GV, US, eta_preale) + call find_eta(h, tv, G, GV, US, eta_preale, dZref=G%Z_ref) call post_data(CS%id_e_preale, eta_preale, CS%diag) endif @@ -1304,7 +1304,7 @@ subroutine ALE_initThicknessToCoord( CS, G, GV, h ) integer :: i, j, k do j = G%jsd,G%jed ; do i = G%isd,G%ied - h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j) ) + h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref ) enddo ; enddo end subroutine ALE_initThicknessToCoord diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 5b19a7549c..e215fde06f 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1113,7 +1113,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) endif ! Local depth (G%bathyT is positive downward) - nominalDepth = G%bathyT(i,j)*GV%Z_to_H + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H ! Determine water column thickness totalThickness = 0.0 @@ -1203,7 +1203,7 @@ subroutine build_sigma_grid( CS, G, GV, h, dzInterface ) endif ! The rest of the model defines grids integrating up from the bottom - nominalDepth = G%bathyT(i,j)*GV%Z_to_H + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H ! Determine water column height totalThickness = 0.0 @@ -1314,7 +1314,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel ! Local depth (G%bathyT is positive downward) - nominalDepth = G%bathyT(i,j)*GV%Z_to_H + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H ! Determine total water column thickness totalThickness = 0.0 @@ -1444,7 +1444,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then - nominalDepth = G%bathyT(i,j) * GV%Z_to_H + nominalDepth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H if (ice_shelf) then totalThickness = 0.0 @@ -1592,7 +1592,7 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then - depth = G%bathyT(i,j) * GV%Z_to_H + depth = (G%bathyT(i,j)+G%Z_ref) * GV%Z_to_H z_col(1) = 0. ! Work downward rather than bottom up do K=1,nz z_col(K+1) = z_col(K) + h(i,j,k) @@ -1718,7 +1718,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) do i = G%isc-1,G%iec+1 ! Local depth - local_depth = G%bathyT(i,j)*GV%Z_to_H + local_depth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H ! Determine water column height total_height = 0.0 diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 8fa4b09fc5..fe3864fc7a 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -144,7 +144,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex zNext(nz+1) = zInt(i,j,nz+1) ! local depth for scaling diffusivity - depth = G%bathyT(i,j) * GV%Z_to_H + depth = (G%bathyT(i,j) + G%Z_ref) * GV%Z_to_H ! initialize del2sigma and the thickness change response to it zero del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 06850dca97..82533e567b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -850,7 +850,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! Determining the time-average sea surface height is part of the algorithm. ! This may be eta_av if Boussinesq, or need to be diagnosed if not. CS%time_in_cycle = CS%time_in_cycle + dt - call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0) + call find_eta(h, CS%tv, G, GV, US, ssh, CS%eta_av_bc, eta_to_m=1.0, dZref=G%Z_ref) do j=js,je ; do i=is,ie CS%ssh_rint(i,j) = CS%ssh_rint(i,j) + dt*ssh(i,j) enddo ; enddo @@ -2278,7 +2278,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ALLOC_(CS%ssh_rint(isd:ied,jsd:jed)) ; CS%ssh_rint(:,:) = 0.0 ALLOC_(CS%ave_ssh_ibc(isd:ied,jsd:jed)) ; CS%ave_ssh_ibc(:,:) = 0.0 - ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0 + ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0 ! -G%Z_ref CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0 ! Use the Wright equation of state by default, unless otherwise specified @@ -2789,9 +2789,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then if (CS%split) then - call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0) + call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0, dZref=G%Z_ref) else - call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0) + call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta_to_m=1.0, dZref=G%Z_ref) endif endif if (CS%split) deallocate(eta) @@ -2852,7 +2852,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) restart_CSp_tmp = restart_CSp call restart_registry_lock(restart_CSp_tmp, unlocked=.true.) allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1)) - call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0) + call find_eta(CS%h, CS%tv, G, GV, US, z_interface, eta_to_m=1.0, dZref=G%Z_ref) call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, & "Interface heights", "meter", z_grid='i') ! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface @@ -3401,10 +3401,10 @@ subroutine extract_surface_state(CS, sfc_state_in) numberOfErrors=0 ! count number of errors do j=js,je ; do i=is,ie if (G%mask2dT(i,j)>0.) then - localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) & + localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) - G%Z_ref & .or. sfc_state%sea_lev(i,j) >= CS%bad_val_ssh_max & .or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max & - .or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) < CS%bad_val_col_thick + .or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick if (use_temperature) localError = localError & .or. sfc_state%SSS(i,j)<0. & .or. sfc_state%SSS(i,j)>=CS%bad_val_sss_max & @@ -3420,7 +3420,7 @@ subroutine extract_surface_state(CS, sfc_state_in) 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), & - 'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & + 'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & 'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), & 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) @@ -3429,7 +3429,7 @@ subroutine extract_surface_state(CS, sfc_state_in) 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), & - 'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & + 'D=',CS%US%Z_to_m*(G%bathyT(i,j)+G%Z_ref), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) endif diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 89a7a1faff..23e58272ed 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -90,10 +90,9 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ !! or atmosphere-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure !! anomaly in each layer due to eta anomalies - !! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to - !! calculate PFu and PFv [H ~> m or kg m-2], with any tidal - !! contributions or compressibility compensation. + !! [L2 T-2 H-1 ~> m4 s-2 kg-1]. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to + !! calculate PFu and PFv [H ~> kg m-2]. ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: & @@ -301,7 +300,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! Find and add the tidal geopotential anomaly. !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth + SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref enddo ; enddo call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z) !$OMP parallel do default(shared) @@ -430,15 +429,16 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm !! or atmosphere-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure !! anomaly in each layer due to eta anomalies - !! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to - !! calculate PFu and PFv [H ~> m or kg m-2], with any - !! tidal contributions or compressibility compensation. + !! [L2 T-2 H-1 ~> m s-2]. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The sea-surface height used to + !! calculate PFu and PFv [H ~> m], with any + !! tidal contributions. ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)) :: & e_tidal, & ! The bottom geopotential anomaly due to tidal forces from ! astronomical sources and self-attraction and loading [Z ~> m]. + SSH, & ! The sea surface height anomaly, in depth units [Z ~> m]. dM ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G)) :: & @@ -451,7 +451,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm dpa, & ! The change in pressure anomaly between the top and bottom ! of a layer [R L2 T-2 ~> Pa]. intz_dpa ! The vertical integral in depth of the pressure anomaly less the - ! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa or kg m-2 Pa]. + ! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa]. real, dimension(SZIB_(G),SZJ_(G)) :: & intx_pa, & ! The zonal integral of the pressure anomaly along the interface ! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa]. @@ -485,7 +485,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. real :: Tl(5) ! copy and T in local stencil [degC] real :: mn_T ! mean of T in local stencil [degC] - real :: mn_T2 ! mean of T**2 in local stencil [degC] + real :: mn_T2 ! mean of T**2 in local stencil [degC2] real :: hl(5) ! Copy of local stencil of H [H ~> m] real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1] real, parameter :: C1_6 = 1.0/6.0 @@ -565,13 +565,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm !$OMP parallel do default(shared) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - e(i,j,1) = -G%bathyT(i,j) + SSH(i,j) = -G%bathyT(i,j) - G%Z_ref enddo do k=1,nz ; do i=Isq,Ieq+1 - e(i,j,1) = e(i,j,1) + h(i,j,k)*GV%H_to_Z + SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo - call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z) + call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z) endif ! Here layer interface heights, e, are calculated. @@ -637,13 +637,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm tv%eqn_of_state, EOSdom) endif do i=Isq,Ieq+1 - dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) + dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * (e(i,j,1) - G%Z_ref) enddo enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * GV%Rlay(1)) * e(i,j,1) + dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * GV%Rlay(1)) * (e(i,j,1) - G%Z_ref) enddo ; enddo endif endif @@ -667,12 +667,12 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (use_p_atm) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1) + p_atm(i,j) + pa(i,j) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref) + p_atm(i,j) enddo ; enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1) + pa(i,j) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref) enddo ; enddo endif !$OMP parallel do default(shared) @@ -700,17 +700,17 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, & useMassWghtInterp=CS%useMassWghtInterp, & - use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom) + use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, & - useMassWghtInterp=CS%useMassWghtInterp) + useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, Z_0p=G%Z_ref) endif !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index ac5cb6c84c..05e68aef12 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -77,8 +77,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb optional, intent(out) :: pbce !< The baroclinic pressure anomaly in !! each layer due to free surface height anomalies, !! [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< Free surface height [H ~> kg m-1]. - + real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to calculate + !! PFu and PFv [H ~> kg m-2]. ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & M, & ! The Montgomery potential, M = (p/rho + gz) [L2 T-2 ~> m2 s-2]. @@ -104,7 +104,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb SSH, & ! The sea surface height anomaly, in depth units [Z ~> m]. e_tidal, & ! Bottom geopotential anomaly due to tidal forces from ! astronomical sources and self-attraction and loading [Z ~> m]. - geopot_bot ! Bottom geopotential relative to time-mean sea level, + geopot_bot ! Bottom geopotential relative to a temporally fixed reference value, ! including any tidal contributions [L2 T-2 ~> m2 s-2]. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). @@ -183,7 +183,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb ! of self-attraction and loading. !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = -G%bathyT(i,j) + SSH(i,j) = -G%bathyT(i,j) - G%Z_ref enddo ; enddo if (use_EOS) then !$OMP parallel do default(shared) @@ -393,6 +393,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, ! the deepest variable density near-surface layer [R ~> kg m-3]. real :: h_star(SZI_(G),SZJ_(G)) ! Layer thickness after compensation ! for compressibility [Z ~> m]. + real :: SSH(SZI_(G),SZJ_(G)) ! The sea surface height anomaly, in depth units [Z ~> m]. real :: e_tidal(SZI_(G),SZJ_(G)) ! Bottom geopotential anomaly due to tidal ! forces from astronomical sources and self- ! attraction and loading, in depth units [Z ~> m]. @@ -444,12 +445,12 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, ! barotropic tides. !$OMP parallel do default(shared) do j=Jsq,Jeq+1 - do i=Isq,Ieq+1 ; e(i,j,1) = -G%bathyT(i,j) ; enddo + do i=Isq,Ieq+1 ; SSH(i,j) = -G%bathyT(i,j) - G%Z_ref ; enddo do k=1,nz ; do i=Isq,Ieq+1 - e(i,j,1) = e(i,j,1) + h(i,j,k)*GV%H_to_Z + SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo - call calc_tidal_forcing(CS%Time, e(:,:,1), e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z) + call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z) endif ! Here layer interface heights, e, are calculated. @@ -664,7 +665,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) - press(i) = -Rho0xG*e(i,j,1) + press(i) = -Rho0xG*(e(i,j,1) - G%Z_ref) enddo call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, & tv%eqn_of_state, EOSdom) @@ -673,7 +674,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) enddo do k=2,nz do i=Isq,Ieq+1 - press(i) = -Rho0xG*e(i,j,K) + press(i) = -Rho0xG*(e(i,j,K) - G%Z_ref) T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k)) S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k)) enddo diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 51f9a5cb85..7a4d1c6bf9 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -685,6 +685,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: isv, iev, jsv, jev ! The valid array size at the end of a step. integer :: stencil ! The stencil size of the algorithm, often 1 or 2. integer :: isvf, ievf, jsvf, jevf, num_cycles + integer :: err_count ! A counter to limit the volume of error messages written to stdout. integer :: i, j, k, n integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -700,6 +701,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB MS%isdw = CS%isdw ; MS%iedw = CS%iedw ; MS%jsdw = CS%jsdw ; MS%jedw = CS%jedw h_neglect = GV%H_subroundoff + err_count = 0 Idt = 1.0 / dt accel_underflow = CS%vel_underflow * Idt @@ -2356,13 +2358,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (GV%Boussinesq) then do j=js,je ; do i=is,ie - if (eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) & - call MOM_error(WARNING, "btstep: eta has dropped below bathyT.") + if ((eta(i,j) < -GV%Z_to_H*G%bathyT(i,j)) .and. (G%mask2dT(i,j) > 0.0)) then + write(mesg,'(ES24.16," vs. ",ES24.16)') GV%H_to_m*eta(i,j), -US%Z_to_m*G%bathyT(i,j) + if (err_count < 2) & + call MOM_error(WARNING, "btstep: eta has dropped below bathyT: "//trim(mesg), all_print=.true.) + err_count = err_count + 1 + endif enddo ; enddo else do j=js,je ; do i=is,ie - if (eta(i,j) < 0.0) & - call MOM_error(WARNING, "btstep: negative eta in a non-Boussinesq barotropic solver.") + if (eta(i,j) < 0.0) then + if (err_count < 2) & + call MOM_error(WARNING, "btstep: negative eta in a non-Boussinesq barotropic solver.", all_print=.true.) + err_count = err_count + 1 + endif enddo ; enddo endif @@ -2566,7 +2575,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%id_vaccel > 0) call post_data(CS%id_vaccel, v_accel_bt(isd:ied,JsdB:JedB), CS%diag) if (CS%id_eta_cor > 0) call post_data(CS%id_eta_cor, CS%eta_cor, CS%diag) - if (CS%id_eta_bt > 0) call post_data(CS%id_eta_bt, eta_out, CS%diag) + if (CS%id_eta_bt > 0) call post_data(CS%id_eta_bt, eta_out, CS%diag) ! - G%Z_ref? if (CS%id_gtotn > 0) call post_data(CS%id_gtotn, gtot_N(isd:ied,jsd:jed), CS%diag) if (CS%id_gtots > 0) call post_data(CS%id_gtots, gtot_S(isd:ied,jsd:jed), CS%diag) if (CS%id_gtote > 0) call post_data(CS%id_gtote, gtot_E(isd:ied,jsd:jed), CS%diag) @@ -3139,7 +3148,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_B if (segment%is_E_or_W .and. segment%Flather) then do j=segment%HI%jsd,segment%HI%jed ; do I=segment%HI%IsdB,segment%HI%IedB BT_OBC%ubt_outer(I,j) = segment%normal_vel_bt(I,j) - BT_OBC%eta_outer_u(I,j) = segment%eta(I,j) + BT_OBC%eta_outer_u(I,j) = segment%eta(I,j) + G%Z_ref*GV%Z_to_H enddo ; enddo endif enddo @@ -3193,7 +3202,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_B if (segment%is_N_or_S .and. segment%Flather) then do J=segment%HI%JsdB,segment%HI%JedB ; do i=segment%HI%isd,segment%HI%ied BT_OBC%vbt_outer(i,J) = segment%normal_vel_bt(i,J) - BT_OBC%eta_outer_v(i,J) = segment%eta(i,J) + BT_OBC%eta_outer_v(i,J) = segment%eta(i,J) + G%Z_ref*GV%Z_to_H enddo ; enddo endif enddo @@ -3268,8 +3277,10 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC) real :: Rh ! A ratio of summed thicknesses, nondim. real :: e_u(SZIB_(G),SZK_(GV)+1) ! The interface heights at u-velocity and real :: e_v(SZI_(G),SZK_(GV)+1) ! v-velocity points [H ~> m or kg m-2]. - real :: D_shallow_u(SZI_(G)) ! The shallower of the adjacent depths [H ~> m or kg m-2]. - real :: D_shallow_v(SZIB_(G))! The shallower of the adjacent depths [H ~> m or kg m-2]. + real :: D_shallow_u(SZI_(G)) ! The height of the shallower of the adjacent bathymetric depths + ! around a u-point (positive upward) [H ~> m or kg m-2] + real :: D_shallow_v(SZIB_(G))! The height of the shallower of the adjacent bathymetric depths + ! around a v-point (positive upward) [H ~> m or kg m-2] real :: htot ! The sum of the layer thicknesses [H ~> m or kg m-2]. real :: Ihtot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1]. @@ -4124,7 +4135,7 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec hs = 1 ; if (present(halo)) hs = max(halo,0) -!$OMP parallel default(none) shared(is,ie,js,je,hs,eta,GV,CS,Datu,Datv,add_max) & +!$OMP parallel default(none) shared(is,ie,js,je,hs,eta,GV,G,CS,Datu,Datv,add_max) & !$OMP private(H1,H2) if (present(eta)) then ! The use of harmonic mean thicknesses ensure positive definiteness. @@ -4163,31 +4174,27 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max) !$OMP do do j=js-hs,je+hs ; do I=is-1-hs,ie+hs Datu(I,j) = CS%dy_Cu(I,j) * GV%Z_to_H * & - (max(CS%bathyT(i+1,j), CS%bathyT(i,j)) + add_max) + max(max(CS%bathyT(i+1,j), CS%bathyT(i,j)) + (G%Z_ref + add_max), 0.0) enddo ; enddo !$OMP do do J=js-1-hs,je+hs ; do i=is-hs,ie+hs Datv(i,J) = CS%dx_Cv(i,J) * GV%Z_to_H * & - (max(CS%bathyT(i,j+1), CS%bathyT(i,j)) + add_max) + max(max(CS%bathyT(i,j+1), CS%bathyT(i,j)) + (G%Z_ref + add_max), 0.0) enddo ; enddo else !$OMP do do j=js-hs,je+hs ; do I=is-1-hs,ie+hs - Datu(I, j) = 0.0 - !Would be "if (G%mask2dCu(I,j)>0.) &" is G was valid on BT domain - if (CS%bathyT(i+1,j)+CS%bathyT(i,j)>0.) & - Datu(I,j) = 2.0*CS%dy_Cu(I,j) * GV%Z_to_H * & - (CS%bathyT(i+1,j) * CS%bathyT(i,j)) / & - (CS%bathyT(i+1,j) + CS%bathyT(i,j)) + H1 = (CS%bathyT(i,j) + G%Z_ref) * GV%Z_to_H ; H2 = (CS%bathyT(i+1,j) + G%Z_ref) * GV%Z_to_H + Datu(I,j) = 0.0 + if ((H1 > 0.0) .and. (H2 > 0.0)) & + Datu(I,j) = CS%dy_Cu(I,j) * (2.0 * H1 * H2) / (H1 + H2) enddo ; enddo !$OMP do do J=js-1-hs,je+hs ; do i=is-hs,ie+hs - Datv(i, J) = 0.0 - !Would be "if (G%mask2dCv(i,J)>0.) &" is G was valid on BT domain - if (CS%bathyT(i,j+1)+CS%bathyT(i,j)>0.) & - Datv(i,J) = 2.0*CS%dx_Cv(i,J) * GV%Z_to_H * & - (CS%bathyT(i,j+1) * CS%bathyT(i,j)) / & - (CS%bathyT(i,j+1) + CS%bathyT(i,j)) + H1 = (CS%bathyT(i,j) + G%Z_ref) * GV%Z_to_H ; H2 = (CS%bathyT(i,j+1) + G%Z_ref) * GV%Z_to_H + Datv(i,J) = 0.0 + if ((H1 > 0.0) .and. (H2 > 0.0)) & + Datv(i,J) = CS%dx_Cv(i,J) * (2.0 * H1 * H2) / (H1 + H2) enddo ; enddo endif !$OMP end parallel @@ -4660,7 +4667,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ! IareaT, IdxCu, and IdyCv need to be allocated with wide halos. ALLOC_(CS%IareaT(CS%isdw:CS%iedw,CS%jsdw:CS%jedw)) ; CS%IareaT(:,:) = 0.0 - ALLOC_(CS%bathyT(CS%isdw:CS%iedw,CS%jsdw:CS%jedw)) ; CS%bathyT(:,:) = GV%Angstrom_m !### Change to 0.0? + ALLOC_(CS%bathyT(CS%isdw:CS%iedw,CS%jsdw:CS%jedw)) ; CS%bathyT(:,:) = 0.0 ALLOC_(CS%IdxCu(CS%isdw-1:CS%iedw,CS%jsdw:CS%jedw)) ; CS%IdxCu(:,:) = 0.0 ALLOC_(CS%IdyCv(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) ; CS%IdyCv(:,:) = 0.0 ALLOC_(CS%dy_Cu(CS%isdw-1:CS%iedw,CS%jsdw:CS%jedw)) ; CS%dy_Cu(:,:) = 0.0 diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 1ac5e39dd5..e672252c24 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -283,10 +283,10 @@ subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_v G%bathymetry_at_vel = .false. if (present(bathymetry_at_vel)) G%bathymetry_at_vel = bathymetry_at_vel if (G%bathymetry_at_vel) then - ALLOC_(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = 0.0 - ALLOC_(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = 0.0 - ALLOC_(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = 0.0 - ALLOC_(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = 0.0 + ALLOC_(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = -G%Z_ref + ALLOC_(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = -G%Z_ref + ALLOC_(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = -G%Z_ref + ALLOC_(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = -G%Z_ref endif ! setup block indices. @@ -387,6 +387,7 @@ end subroutine MOM_grid_init subroutine rescale_grid_bathymetry(G, m_in_new_units) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid structure real, intent(in) :: m_in_new_units !< The new internal representation of 1 m depth. + ! It appears that this routine is never called. ! Local variables real :: rescale @@ -578,7 +579,7 @@ subroutine allocate_metrics(G) ALLOC_(G%IareaCu(IsdB:IedB,jsd:jed)) ; G%IareaCu(:,:) = 0.0 ALLOC_(G%IareaCv(isd:ied,JsdB:JedB)) ; G%IareaCv(:,:) = 0.0 - ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = 0.0 + ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = -G%Z_ref ALLOC_(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0 ALLOC_(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0 ALLOC_(G%dF_dy(isd:ied, jsd:jed)) ; G%dF_dy(:,:) = 0.0 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 059677c6f7..b83c4d1be8 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4726,6 +4726,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) + ! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref allocate(color(G%isd:G%ied, G%jsd:G%jed)) ; color = 0 allocate(color2(G%isd:G%ied, G%jsd:G%jed)) ; color2 = 0 diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index 51d44c1041..a9626a805c 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -20,7 +20,8 @@ module MOM_transcribe_grid contains !> Copies information from a dynamic (shared) horizontal grid type into an -!! ocean_grid_type. +!! ocean_grid_type. There may also be a change in the reference +!! height for topography between the two grids. subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) type(dyn_horgrid_type), intent(in) :: dG !< Common horizontal grid type type(ocean_grid_type), intent(inout) :: oG !< Ocean grid type @@ -54,7 +55,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%dxT(i,j) = dG%dxT(i+ido,j+jdo) oG%dyT(i,j) = dG%dyT(i+ido,j+jdo) oG%areaT(i,j) = dG%areaT(i+ido,j+jdo) - oG%bathyT(i,j) = dG%bathyT(i+ido,j+jdo) + oG%bathyT(i,j) = dG%bathyT(i+ido,j+jdo) - oG%Z_ref oG%dF_dx(i,j) = dG%dF_dx(i+ido,j+jdo) oG%dF_dy(i,j) = dG%dF_dy(i+ido,j+jdo) @@ -100,12 +101,12 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%bathymetry_at_vel = dG%bathymetry_at_vel if (oG%bathymetry_at_vel) then do I=IsdB,IedB ; do j=jsd,jed - oG%Dblock_u(I,j) = dG%Dblock_u(I+ido,j+jdo) - oG%Dopen_u(I,j) = dG%Dopen_u(I+ido,j+jdo) + oG%Dblock_u(I,j) = dG%Dblock_u(I+ido,j+jdo) - oG%Z_ref + oG%Dopen_u(I,j) = dG%Dopen_u(I+ido,j+jdo) - oG%Z_ref enddo ; enddo do i=isd,ied ; do J=JsdB,JedB - oG%Dblock_v(i,J) = dG%Dblock_v(i+ido,J+jdo) - oG%Dopen_v(i,J) = dG%Dopen_v(i+ido,J+jdo) + oG%Dblock_v(i,J) = dG%Dblock_v(i+ido,J+jdo) - oG%Z_ref + oG%Dopen_v(i,J) = dG%Dopen_v(i+ido,J+jdo) - oG%Z_ref enddo ; enddo endif @@ -164,7 +165,8 @@ end subroutine copy_dyngrid_to_MOM_grid !> Copies information from an ocean_grid_type into a dynamic (shared) -!! horizontal grid type. +!! horizontal grid type. There may also be a change in the reference +!! height for topography between the two grids. subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) type(ocean_grid_type), intent(in) :: oG !< Ocean grid type type(dyn_horgrid_type), intent(inout) :: dG !< Common horizontal grid type @@ -198,7 +200,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%dxT(i,j) = oG%dxT(i+ido,j+jdo) dG%dyT(i,j) = oG%dyT(i+ido,j+jdo) dG%areaT(i,j) = oG%areaT(i+ido,j+jdo) - dG%bathyT(i,j) = oG%bathyT(i+ido,j+jdo) + dG%bathyT(i,j) = oG%bathyT(i+ido,j+jdo) + oG%Z_ref dG%dF_dx(i,j) = oG%dF_dx(i+ido,j+jdo) dG%dF_dy(i,j) = oG%dF_dy(i+ido,j+jdo) @@ -244,12 +246,12 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%bathymetry_at_vel = oG%bathymetry_at_vel if (dG%bathymetry_at_vel) then do I=IsdB,IedB ; do j=jsd,jed - dG%Dblock_u(I,j) = oG%Dblock_u(I+ido,j+jdo) - dG%Dopen_u(I,j) = oG%Dopen_u(I+ido,j+jdo) + dG%Dblock_u(I,j) = oG%Dblock_u(I+ido,j+jdo) + oG%Z_ref + dG%Dopen_u(I,j) = oG%Dopen_u(I+ido,j+jdo) + oG%Z_ref enddo ; enddo do i=isd,ied ; do J=JsdB,JedB - dG%Dblock_v(i,J) = oG%Dblock_v(i+ido,J+jdo) - dG%Dopen_v(i,J) = oG%Dopen_v(i+ido,J+jdo) + dG%Dblock_v(i,J) = oG%Dblock_v(i+ido,J+jdo) + oG%Z_ref + dG%Dopen_v(i,J) = oG%Dopen_v(i+ido,J+jdo) + oG%Z_ref enddo ; enddo endif diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index 45b08cc799..b5a1a6bf0c 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -248,13 +248,13 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j+1,k)); enddo - e(nz+1) = -US%Z_to_m*G%bathyT(i,j) + e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo write(file,'(/,"e-: ",$)') write(file,'(ES10.3," ",$)') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K); enddo - e(nz+1) = -US%Z_to_m*G%bathyT(i+1,j) + e(nz+1) = -US%Z_to_m*(G%bathyT(i+1,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i+1,j,k) ; enddo write(file,'(/,"e+: ",$)') write(file,'(ES10.3," ",$)') e(ks) @@ -331,7 +331,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp (0.5*US%L_T_to_m_s*CS%v_av(i+1,J,k)*h_scale*(hin(i+1,j,k) + hin(i+1,j+1,k))); enddo endif - write(file,'(/,"D: ",2(ES10.3))') US%Z_to_m*G%bathyT(i,j),US%Z_to_m*G%bathyT(i+1,j) + write(file,'(/,"D: ",2(ES10.3))') US%Z_to_m*(G%bathyT(i,j) + G%Z_ref), US%Z_to_m*(G%bathyT(i+1,j) + G%Z_ref) ! From here on, the normalized accelerations are written. if (prev_avail) then @@ -584,13 +584,13 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp write(file,'(/,"h++: ",$)') do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i+1,j+1,k); enddo - e(nz+1) = -US%Z_to_m*G%bathyT(i,j) + e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k); enddo write(file,'(/,"e-: ",$)') write(file,'(ES10.3," ",$)') e(ks) do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(K); enddo - e(nz+1) = -US%Z_to_m*G%bathyT(i,j+1) + e(nz+1) = -US%Z_to_m*(G%bathyT(i,j+1) + G%Z_ref) do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j+1,k) ; enddo write(file,'(/,"e+: ",$)') write(file,'(ES10.3," ",$)') e(ks) @@ -667,7 +667,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp (CS%u_prev(I,j+1,k) * h_scale*0.5*(hin(i,j+1,k) + hin(i+1,j+1,k))); enddo endif - write(file,'(/,"D: ",2(ES10.3))') US%Z_to_m*G%bathyT(i,j),US%Z_to_m*G%bathyT(i,j+1) + write(file,'(/,"D: ",2(ES10.3))') US%Z_to_m*(G%bathyT(i,j) + G%Z_ref), US%Z_to_m*(G%bathyT(i,j+1) + G%Z_ref) ! From here on, the normalized accelerations are written. if (prev_avail) then diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 44b05cc081..5d0d800c1d 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -386,14 +386,14 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & endif if (associated(CS%e)) then - call find_eta(h, tv, G, GV, US, CS%e) + call find_eta(h, tv, G, GV, US, CS%e, dZref=G%Z_ref) if (CS%id_e > 0) call post_data(CS%id_e, CS%e, CS%diag) endif if (associated(CS%e_D)) then if (associated(CS%e)) then do k=1,nz+1 ; do j=js,je ; do i=is,ie - CS%e_D(i,j,k) = CS%e(i,j,k) + G%bathyT(i,j) + CS%e_D(i,j,k) = CS%e(i,j,k) + (G%bathyT(i,j) + G%Z_ref) enddo ; enddo ; enddo else call find_eta(h, tv, G, GV, US, CS%e_D) @@ -2132,7 +2132,8 @@ subroutine write_static_fields(G, GV, US, tv, diag) type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output ! Local variables - integer :: id + real :: work_2d(SZI_(G),SZJ_(G)) ! A 2-d temporary work array. + integer :: id, i, j logical :: use_temperature id = register_static_field('ocean_model', 'geolat', diag%axesT1, & @@ -2204,7 +2205,10 @@ subroutine write_static_fields(G, GV, US, tv, diag) cmor_field_name='deptho', cmor_long_name='Sea Floor Depth', & cmor_standard_name='sea_floor_depth_below_geoid', area=diag%axesT1%id_area, & x_cell_method='mean', y_cell_method='mean', area_cell_method='mean') - if (id > 0) call post_data(id, G%bathyT, diag, .true., mask=G%mask2dT) + if (id > 0) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; work_2d(i,j) = G%bathyT(i,j)+G%Z_ref ; enddo ; enddo + call post_data(id, work_2d, diag, .true., mask=G%mask2dT) + endif id = register_static_field('ocean_model', 'wet', diag%axesT1, & '0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 72523edfd3..d190cee7a3 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -542,7 +542,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ enddo ; enddo ; enddo mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP) - call find_eta(h, tv, G, GV, US, eta) + call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) do k=1,nz ; do j=js,je ; do i=is,ie tmp1(i,j,k) = US%Z_to_m*US%L_to_m**2*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo @@ -674,8 +674,8 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ hbelow = 0.0 do k=nz,1,-1 hbelow = hbelow + h(i,j,k) * GV%H_to_Z - hint = Z_0APE(K) + (hbelow - G%bathyT(i,j)) - hbot = Z_0APE(K) - G%bathyT(i,j) + hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref)) + hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref) hbot = (hbot + ABS(hbot)) * 0.5 PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j)) * (GV%Rho0*GV%g_prime(K)) * & (hint * hint - hbot * hbot) @@ -685,7 +685,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ do j=js,je ; do i=is,ie do k=nz,1,-1 hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. - hbot = max(Z_0APE(K) - G%bathyT(i,j), 0.0) + hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0) PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * & (hint * hint - hbot * hbot) enddo @@ -1166,7 +1166,7 @@ subroutine create_depth_list(G, DL, min_depth_inc) i_global = i + G%idg_offset - (G%isg-1) list_pos = (j_global-1)*G%Domain%niglobal + i_global - Dlist(list_pos) = G%bathyT(i,j) + Dlist(list_pos) = G%bathyT(i,j) + G%Z_ref Arealist(list_pos) = G%mask2dT(i,j) * G%areaT(i,j) enddo ; enddo @@ -1401,7 +1401,7 @@ subroutine get_depth_list_checksums(G, depth_chksum, area_chksum) ! Depth checksum do j=G%jsc,G%jec ; do i=G%isc,G%iec - field(i,j) = G%bathyT(i,j) + field(i,j) = G%bathyT(i,j) + G%Z_ref enddo ; enddo write(depth_chksum, '(Z16)') field_chksum(field(:,:)) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index d363b185f8..6a4d9660d7 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -445,7 +445,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ gp = gprime(K) if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then !### 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. & + if ( (((G%bathyT(i,j)+G%Z_ref) - sum_hc < l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) .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/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index d3eb21dcbe..bb11d92673 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -321,22 +321,22 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then call build_zstar_column(get_zlike_CS(remap_cs%regrid_cs), & - GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), & + GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), & zInterfaces, zScale=GV%Z_to_H) elseif (remap_cs%vertical_coord == coordinateMode('SIGMA')) then call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & - GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) + GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then call build_rho_column(get_rho_CS(remap_cs%regrid_cs), GV%ke, & - GV%Z_to_H*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & + GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then ! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & -! GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!") elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then ! call build_hycom1_column(remap_cs%regrid_cs, nz, & -! GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") endif do k = 1,nz diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index d1a4b7f45d..a18c8442b9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -450,10 +450,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, allocate(mask_in(id,jdp)) ; mask_in(:,:) = 0.0 endif - max_depth = maxval(G%bathyT) + max_depth = maxval(G%bathyT(:,:)) + G%Z_ref call max_across_PEs(max_depth) - if (z_edges_in(kd+1) CS%Grid ; CS%Grid_in => CS%Grid diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 1f8d45e88d..89c25f7525 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -89,9 +89,10 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: h_bdry_val => NULL() !< The ice thickness at inflowing boundaries [m]. real, pointer, dimension(:,:) :: t_bdry_val => NULL() !< The ice temperature at inflowing boundaries [degC]. - real, pointer, dimension(:,:) :: bed_elev => NULL() !< The bed elevation used for ice dynamics [Z ~> m]. - !! the same as bathyT, when below sea-level. - !!Sign convention: positive below sea-level, negative above. + real, pointer, dimension(:,:) :: bed_elev => NULL() !< The bed elevation used for ice dynamics [Z ~> m], + !! relative to mean sea-level. This is + !! the same as G%bathyT+Z_ref, when below sea-level. + !! Sign convention: positive below sea-level, negative above. real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area integrated nonlinear part of "linearized" !! basal stress [R Z L2 T-1 ~> kg s-1]. @@ -266,7 +267,7 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) allocate( CS%ground_frac(isd:ied,jsd:jed) ) ; CS%ground_frac(:,:) = 0.0 allocate( CS%taudx_shelf(IsdB:IedB,JsdB:JedB) ) ; CS%taudx_shelf(:,:) = 0.0 allocate( CS%taudy_shelf(IsdB:IedB,JsdB:JedB) ) ; CS%taudy_shelf(:,:) = 0.0 - allocate( CS%bed_elev(isd:ied,jsd:jed) ) ; CS%bed_elev(:,:)=G%bathyT(:,:)!CS%bed_elev(:,:) = 0.0 + allocate( CS%bed_elev(isd:ied,jsd:jed) ) ; CS%bed_elev(:,:) = G%bathyT(:,:) + G%Z_ref allocate( CS%u_bdry_val(IsdB:IedB,JsdB:JedB) ) ; CS%u_bdry_val(:,:) = 0.0 allocate( CS%v_bdry_val(IsdB:IedB,JsdB:JedB) ) ; CS%v_bdry_val(:,:) = 0.0 allocate( CS%u_face_mask_bdry(IsdB:IedB,JsdB:JedB) ) ; CS%u_face_mask_bdry(:,:) = -2.0 @@ -611,7 +612,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) do j=jsd,jed do i=isd,ied - OD = G%bathyT(i,j) - rhoi_rhow * ISS%h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * ISS%h_shelf(i,j) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -816,7 +817,7 @@ end subroutine ice_shelf_advect !>This subroutine computes u- and v-velocities of the ice shelf iterating on non-linear ice viscosity !subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) - subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, iters, Time) + subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, iters, Time) type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state @@ -869,13 +870,13 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0 - do j=G%jsc,G%jec - do i=G%isc,G%iec - if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then - float_cond(i,j) = 1.0 - CS%ground_frac(i,j) = 1.0 - endif - enddo + do j=G%jsc,G%jec + do i=G%isc,G%iec + if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then + float_cond(i,j) = 1.0 + CS%ground_frac(i,j) = 1.0 + endif + enddo enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) @@ -896,7 +897,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite do l=0,1 ; do k=0,1 if ((ISS%hmask(i,j) == 1) .and. & - (rhoi_rhow * H_node(i-1+k,j-1+l) - G%bathyT(i,j) <= 0)) then + (rhoi_rhow * H_node(i-1+k,j-1+l) - CS%bed_elev(i,j) <= 0)) then nodefloat = nodefloat + 1 endif enddo ; enddo @@ -936,7 +937,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite Au(:,:) = 0.0 ; Av(:,:) = 0.0 call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & + CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow) call pass_vector(Au,Av,G%domain,TO_ALL,BGRID_NE) @@ -992,7 +993,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite Au(:,:) = 0 ; Av(:,:) = 0 call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & + CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow) err_max = 0 @@ -1162,7 +1163,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H call pass_vector(DIAGu, DIAGv, G%domain, TO_ALL, BGRID_NE) call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & G, US, isc-1, iec+1, jsc-1, jec+1, rhoi_rhow) call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE) @@ -1216,7 +1217,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H Au(:,:) = 0 ; Av(:,:) = 0 call CG_action(Au, Av, Du, Dv, Phi, Phisub, CS%umask, CS%vmask, hmask, & - H_node, CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & + H_node, CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & G, US, is, ie, js, je, rhoi_rhow) ! Au, Av valid region moves in by 1 @@ -1820,21 +1821,16 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! prelim - go through and calculate S ! or is this faster? - !BASE(:,:) = -G%bathyT(:,:) + OD(:,:) BASE(:,:) = -CS%bed_elev(:,:) + OD(:,:) S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) ! check whether the ice is floating or grounded do j=jsc-G%domain%njhalo,jec+G%domain%njhalo - do i=isc-G%domain%nihalo,iec+G%domain%nihalo - -! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then - if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then - S(i,j)=(1 - rhoi_rhow)*ISS%h_shelf(i,j) - endif - - - enddo + do i=isc-G%domain%nihalo,iec+G%domain%nihalo + if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then + S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) + endif + enddo enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 @@ -1935,7 +1931,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) endif if (CS%ground_frac(i,j) == 1) then -! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) +! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2) neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 else neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 @@ -2086,8 +2082,9 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas intent(in) :: float_cond !< An array indicating where the ice !! shelf is floating: 0 if floating, 1 if not. real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points + !! relative to sea-level [Z ~> m]. + real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: basal_trac !< A field related to the nonlinear part of the !! "linearized" basal stress [R Z T-1 ~> kg m-2 s-1]. @@ -2206,7 +2203,8 @@ subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, real, dimension(2,2), intent(in) :: H !< The ice shelf thickness at nodal (corner) points [Z ~> m]. real, dimension(2,2), intent(in) :: U !< The zonal ice shelf velocity at vertices [L T-1 ~> m s-1] real, dimension(2,2), intent(in) :: V !< The meridional ice shelf velocity at vertices [L T-1 ~> m s-1] - real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m]. + real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points + !! relative to sea-level [Z ~> m]. real, intent(in) :: dens_ratio !< The density of ice divided by the density !! of seawater [nondim] real, dimension(2,2), intent(out) :: Ucontr !< The areal average of u-velocities where the ice shelf @@ -2338,7 +2336,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_diagonal_subgrid_basal(Phisub, Hcell, G%bathyT(i,j), dens_ratio, sub_ground) + call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground) do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi if (CS%umask(Itgt,Jtgt) == 1) then u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) @@ -2512,7 +2510,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, if (float_cond(i,j) == 1) then Ucell(:,:) = CS%u_bdry_val(i-1:i,j-1:j) ; Vcell(:,:) = CS%v_bdry_val(i-1:i,j-1:j) Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, G%bathyT(i,j), & + call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, CS%bed_elev(i,j), & dens_ratio, Usubcontr, Vsubcontr) if (CS%umask(I-1,J-1)==1) u_bdry_contr(I-1,J-1) = u_bdry_contr(I-1,J-1) + Usubcontr(1,1) * basal_trac(i,j) @@ -2693,7 +2691,7 @@ subroutine update_OD_ffrac_uncoupled(CS, G, h_shelf) do j=jsd,jed do i=isd,ied - OD = G%bathyT(i,j) - rhoi_rhow * h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * h_shelf(i,j) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1907a75c74..1bb03bdd85 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -232,7 +232,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! geopotential height, for use by the various initialization routines. G%bathyT has ! already been initialized in previous calls. do j=jsd,jed ; do i=isd,ied - depth_tot(i,j) = G%bathyT(i,j) + depth_tot(i,j) = G%bathyT(i,j) + G%Z_ref enddo ; enddo ! The remaining initialization calls are done, regardless of whether the @@ -706,7 +706,7 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f 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, dZ_ref_eta=0.0) + call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=G%Z_ref) 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 @@ -1086,7 +1086,7 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params) call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor) ! Convert thicknesses to interface heights. - call find_eta(h, tv, G, GV, US, eta) + call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then ! if (eta_sfc(i,j) < eta(i,j,nz+1)) then @@ -1199,7 +1199,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params) endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j), & + call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, & min_thickness, tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), & tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:), remap_CS, & z_tol=1.0e-5*US%m_to_Z, remap_answers_2018=remap_answers_2018) @@ -2658,7 +2658,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just Hmix_depth, eps_z, eps_rho, density_extrap_bug) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, zi, h, dZ_ref_eta=0.0) + call adjustEtaToFitBathymetry(G, GV, US, zi, h, dZ_ref_eta=G%Z_ref) 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 diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 48b67bf295..8a67d71fe2 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -152,7 +152,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ if (G%mask2dT(i,j)>0.) then ! Build the source grid zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 - z_bathy = G%bathyT(i,j) + z_bathy = G%bathyT(i,j) + G%Z_ref do k = 1, kd if (mask_z(i,j,k) > 0.) then zBottomOfCell = -min( z_edges_in(k+1), z_bathy ) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 8206fd9717..b4f850904d 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -288,7 +288,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%fixed_total_depth) then !$OMP parallel do default(shared) do j=js-1,je+1 ; do i=is-1,ie+1 - depth_tot(i,j) = G%bathyT(i,j) + depth_tot(i,j) = G%bathyT(i,j) + G%Z_ref enddo ; enddo else !$OMP parallel do default(shared) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 0e3ae5ab71..abad05f837 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -2359,7 +2359,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) do j=G%jsc,G%jec ; do i=G%isc,G%iec ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth. if (RMS_roughness_frac >= 0.0) then - h2(i,j) = max(min((RMS_roughness_frac*G%bathyT(i,j))**2, h2(i,j)), 0.0) + h2(i,j) = max(min((RMS_roughness_frac*(G%bathyT(i,j)+G%Z_ref))**2, h2(i,j)), 0.0) else h2(i,j) = max(h2(i,j), 0.0) endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index fb70f5d679..9af47b3cea 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -184,11 +184,11 @@ subroutine calc_depth_function(G, CS) expo = CS%depth_scaled_khth_exp !$OMP do do j=js,je ; do I=is-1,Ieq - CS%Depth_fn_u(I,j) = (MIN(1.0, 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j))/H0))**expo + CS%Depth_fn_u(I,j) = (MIN(1.0, (0.5*(G%bathyT(i,j) + G%bathyT(i+1,j)) + G%Z_ref)/H0))**expo enddo ; enddo !$OMP do do J=js-1,Jeq ; do i=is,ie - CS%Depth_fn_v(i,J) = (MIN(1.0, 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1))/H0))**expo + CS%Depth_fn_v(i,J) = (MIN(1.0, (0.5*(G%bathyT(i,j) + G%bathyT(i,j+1)) + G%Z_ref)/H0))**expo enddo ; enddo end subroutine calc_depth_function @@ -958,11 +958,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop enddo ; enddo ! SN above contains S^2*N^2*H, convert to vertical average of S*N do I=is-1,ie - !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom_Z ) )) + !### Replace G%bathT+G%Z_ref here with (e(i,j,1) - e(i,j,nz+1)). + !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(i,j), G%bathyT(i+1,j)) + (G%Z_ref + GV%Angstrom_Z) ) ) !The code below behaves better than the line above. Not sure why? AJA - if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff*GV%H_to_Z ) then + if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then CS%SN_u(I,j) = G%mask2dCu(I,j) * sqrt( CS%SN_u(I,j) / & - (max(G%bathyT(I,j), G%bathyT(I+1,j))) ) + (max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) ) else CS%SN_u(I,j) = 0.0 endif @@ -984,20 +985,21 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k) enddo ; enddo do i=is,ie - !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom_Z ) )) + !### Replace G%bathT+G%Z_ref here with (e(i,j,1) - e(i,j,nz+1)). + !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + (G%Z_ref + GV%Angstrom_Z) ) ) !The code below behaves better than the line above. Not sure why? AJA - if ( min(G%bathyT(I,j), G%bathyT(I+1,j)) > H_cutoff*GV%H_to_Z ) then + if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then CS%SN_v(i,J) = G%mask2dCv(i,J) * sqrt( CS%SN_v(i,J) / & - (max(G%bathyT(i,J), G%bathyT(i,J+1))) ) + (max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) ) else - CS%SN_v(I,j) = 0.0 + CS%SN_v(i,J) = 0.0 endif if (local_open_v_BC) then - l_seg = OBC%segnum_v(I,j) + l_seg = OBC%segnum_v(i,J) if (l_seg /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(I,j))%open) then - CS%SN_v(I,j) = 0. + if (OBC%segment(OBC%segnum_v(i,J))%open) then + CS%SN_v(i,J) = 0. endif endif endif diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index be1d265620..66610ba316 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -918,7 +918,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 do k=1,nz_data if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) + zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(CS%col_i(c),CS%col_j(c)) ) tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) elseif (k>1) then zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) @@ -1026,7 +1026,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_u(i,j,k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) + zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(i,j) ) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) else ! This next block should only ever be reached over land @@ -1074,7 +1074,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 do k=1,nz_data if (mask_v(i,j,k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) + zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(i,j) ) elseif (k>1) then zBottomOfCell = -G%bathyT(i,j) else ! This next block should only ever be reached over land diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 040f3c4ecb..f8e071e41d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -316,7 +316,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_T_predia > 0) call post_data(CS%id_T_predia, tv%T, CS%diag) if (CS%id_S_predia > 0) call post_data(CS%id_S_predia, tv%S, CS%diag) if (CS%id_e_predia > 0) then - call find_eta(h, tv, G, GV, US, eta, eta_to_m=1.0) + call find_eta(h, tv, G, GV, US, eta, eta_to_m=1.0, dZref=G%Z_ref) call post_data(CS%id_e_predia, eta, CS%diag) endif diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 68e2f39f0e..a6835d42ed 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -435,13 +435,13 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) do j=js,je ; do i=is,ie mask_itidal = 1.0 - if (G%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0 + if (G%bathyT(i,j) + G%Z_ref < min_zbot_itides) mask_itidal = 0.0 itide%tideamp(i,j) = itide%tideamp(i,j) * mask_itidal * G%mask2dT(i,j) ! Restrict rms topo to a fraction (often 10 percent) of the column depth. if (max_frac_rough >= 0.0) & - itide%h2(i,j) = min((max_frac_rough*G%bathyT(i,j))**2, itide%h2(i,j)) + itide%h2(i,j) = min((max_frac_rough*(G%bathyT(i,j)+G%Z_ref))**2, itide%h2(i,j)) ! Compute the fixed part of internal tidal forcing; units are [R Z3 T-2 ~> J m-2] here. CS%TKE_itidal_coef(i,j) = 0.5*US%L_to_Z*kappa_h2_factor*GV%Rho0*& diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 99bd91d8f8..1618dd442d 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -228,11 +228,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) real :: Rhtot ! Running sum of thicknesses times the layer potential ! densities [H R ~> kg m-2 or kg2 m-5]. real, dimension(SZIB_(G),SZJ_(G)) :: & - D_u, & ! Bottom depth interpolated to u points [Z ~> m]. + D_u, & ! Bottom depth linearly interpolated to u points [Z ~> m]. mask_u ! A mask that disables any contributions from u points that ! are land or past open boundary conditions [nondim], 0 or 1. real, dimension(SZI_(G),SZJB_(G)) :: & - D_v, & ! Bottom depth interpolated to v points [Z ~> m]. + D_v, & ! Bottom depth linearly interpolated to v points [Z ~> m]. mask_v ! A mask that disables any contributions from v points that ! are land or past open boundary conditions [nondim], 0 or 1. real, dimension(SZIB_(G),SZK_(GV)) :: & @@ -399,12 +399,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) !$OMP parallel do default(shared) do J=js-1,je ; do i=is-1,ie+1 - D_v(i,J) = 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1)) + D_v(i,J) = 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1)) + G%Z_ref mask_v(i,J) = G%mask2dCv(i,J) enddo ; enddo !$OMP parallel do default(shared) do j=js-1,je+1 ; do I=is-1,ie - D_u(I,j) = 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j)) + D_u(I,j) = 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j)) + G%Z_ref mask_u(I,j) = G%mask2dCu(I,j) enddo ; enddo @@ -414,13 +414,13 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB if (OBC%segment(n)%is_N_or_S .and. (J >= js-1) .and. (J <= je)) then do i = max(is-1,OBC%segment(n)%HI%isd), min(ie+1,OBC%segment(n)%HI%ied) - if (OBC%segment(n)%direction == OBC_DIRECTION_N) D_v(i,J) = G%bathyT(i,j) - if (OBC%segment(n)%direction == OBC_DIRECTION_S) D_v(i,J) = G%bathyT(i,j+1) + if (OBC%segment(n)%direction == OBC_DIRECTION_N) D_v(i,J) = G%bathyT(i,j) + G%Z_ref + if (OBC%segment(n)%direction == OBC_DIRECTION_S) D_v(i,J) = G%bathyT(i,j+1) + G%Z_ref enddo elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-1) .and. (I <= ie)) then do j = max(js-1,OBC%segment(n)%HI%jsd), min(je+1,OBC%segment(n)%HI%jed) - if (OBC%segment(n)%direction == OBC_DIRECTION_E) D_u(I,j) = G%bathyT(i,j) - if (OBC%segment(n)%direction == OBC_DIRECTION_W) D_u(I,j) = G%bathyT(i+1,j) + if (OBC%segment(n)%direction == OBC_DIRECTION_E) D_u(I,j) = G%bathyT(i,j) + G%Z_ref + if (OBC%segment(n)%direction == OBC_DIRECTION_W) D_u(I,j) = G%bathyT(i+1,j) + G%Z_ref enddo endif enddo ; endif @@ -809,6 +809,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) ! The drag within the bottommost bbl_thick 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. + + !### The harmonic mean edge depths here are not invariant to offsets! if (m==1) then D_vel = D_u(I,j) tmp = G%mask2dCu(I,j+1) * D_u(I,j+1) diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index 25b5406449..ebb9575974 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -409,17 +409,17 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml) enddo ; enddo ; enddo do j=js,je do i=is,ie - dilate(i) = G%bathyT(i,j) / (e_D(i,j,1) + G%bathyT(i,j)) + dilate(i) = (G%bathyT(i,j) + G%Z_ref) / (e_D(i,j,1) + G%bathyT(i,j)) enddo do k=1,nz+1 ; do i=is,ie - e_D(i,j,K) = dilate(i) * (e_D(i,j,K) + G%bathyT(i,j)) - G%bathyT(i,j) + e_D(i,j,K) = dilate(i) * (e_D(i,j,K) + G%bathyT(i,j)) - (G%bathyT(i,j) + G%Z_ref) enddo ; enddo enddo do k=2,nz do j=js,je ; do i=is,ie eta_anom(i,j) = e_D(i,j,k) - CS%Ref_eta_im(j,k) - if (CS%Ref_eta_im(j,K) < -G%bathyT(i,j)) eta_anom(i,j) = 0.0 + if (CS%Ref_eta_im(j,K) < -(G%bathyT(i,j) + G%Z_ref)) eta_anom(i,j) = 0.0 enddo ; enddo call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G, tmp_scale=US%Z_to_m) enddo diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 21562817c0..797ceb9a35 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -484,16 +484,16 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.1) do j=js,je ; do i=is,ie - if (G%bathyT(i,j) < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0 + if (G%bathyT(i,j)+G%Z_ref < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0 CS%tideamp(i,j) = CS%tideamp(i,j) * CS%mask_itidal(i,j) * G%mask2dT(i,j) ! Restrict rms topo to a fraction (often 10 percent) of the column depth. if (CS%answers_2018 .and. (max_frac_rough >= 0.0)) then - hamp = min(max_frac_rough*G%bathyT(i,j), sqrt(CS%h2(i,j))) + hamp = min(max_frac_rough*(G%bathyT(i,j)+G%Z_ref), sqrt(CS%h2(i,j))) CS%h2(i,j) = hamp*hamp else if (max_frac_rough >= 0.0) & - CS%h2(i,j) = min((max_frac_rough*G%bathyT(i,j))**2, CS%h2(i,j)) + CS%h2(i,j) = min((max_frac_rough*(G%bathyT(i,j)+G%Z_ref))**2, CS%h2(i,j)) endif utide = CS%tideamp(i,j) @@ -1678,7 +1678,7 @@ subroutine read_tidal_constituents(G, US, tidal_energy_file, CS) CS%h_src(k) = US%Z_to_m*(z_t(k)-z_w(k))*2.0 ! form tidal_qe_3d_in from weighted tidal constituents do j=js,je ; do i=is,ie - if ((z_t(k) <= G%bathyT(i,j)) .and. (z_w(k) > CS%tidal_diss_lim_tc)) & + if ((z_t(k) <= G%bathyT(i,j) + G%Z_ref) .and. (z_w(k) > CS%tidal_diss_lim_tc)) & CS%tidal_qe_3d_in(i,j,k) = C1_3*tc_m2(i,j,k) + C1_3*tc_s2(i,j,k) + & tidal_qk1(i,j)*tc_k1(i,j,k) + tidal_qo1(i,j)*tc_o1(i,j,k) enddo ; enddo @@ -1692,7 +1692,7 @@ subroutine read_tidal_constituents(G, US, tidal_energy_file, CS) ! do k=50,nz_in(1) ! write(1905,*) i,j,k ! write(1905,*) CS%tidal_qe_3d_in(i,j,k), tc_m2(i,j,k) - ! write(1905,*) z_t(k), G%bathyT(i,j), z_w(k),CS%tidal_diss_lim_tc + ! write(1905,*) z_t(k), G%bathyT(i,j)+G%Z_ref, z_w(k),CS%tidal_diss_lim_tc ! end do ! endif ! enddo @@ -1707,7 +1707,7 @@ subroutine read_tidal_constituents(G, US, tidal_energy_file, CS) !! collapse 3D q*E to 2D q*E !CS%tidal_qe_2d(:,:) = 0.0 !do k=1,nz_in(1) ; do j=js,je ; do i=is,ie - ! if (z_t(k) <= G%bathyT(i,j)) & + ! if (z_t(k) <= G%bathyT(i,j) + G%Z_ref) & ! CS%tidal_qe_2d(i,j) = CS%tidal_qe_2d(i,j) + CS%tidal_qe_3d_in(i,j,k) !enddo ; enddo ; enddo diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index ad0a997cc4..7fb71f9773 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -137,8 +137,8 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va do i=is,ie ; if (G%mask2dT(i,j)*htot(i) > 0.0) then ! Determine the z* heights of the model interfaces. - dilate = (G%bathyT(i,j) - 0.0) / htot(i) - e(nz+1) = -G%bathyT(i,j) + dilate = (G%bathyT(i,j) + G%Z_ref) / htot(i) + e(nz+1) = -G%bathyT(i,j) - G%Z_ref do k=nz,1,-1 ; e(K) = e(K+1) + dilate * h(i,j,k) ; enddo ! Create a single-column copy of tr_in. Efficiency is not an issue here. @@ -212,8 +212,8 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va do i=is,ie ; if (G%mask2dT(i,j)*htot(i) > 0.0) then ! Determine the z* heights of the model interfaces. - dilate = (G%bathyT(i,j) - 0.0) / htot(i) - e(nz+1) = -G%bathyT(i,j) + dilate = (G%bathyT(i,j) + G%Z_ref) / htot(i) + e(nz+1) = -G%bathyT(i,j) - G%Z_ref do k=nz,1,-1 ; e(K) = e(K+1) + dilate * h(i,j,k) ; enddo ! Create a single-column copy of tr_in. Efficiency is not an issue here.