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.