+Changed units of GV%Rlay to [R]

  Changed the units of GV%Rlay from [kg m-3] to [R] for dimensional consistency
testing.  This required the addition of unit_scale_type arguments to several
interfaces.  All answers are bitwise identical, but new arguments have been
added to several public interfaces.

src/ALE/MOM_regridding.F90

@@ -311,11 +311,11 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
     endif
     allocate(dz(ke))
     if (ke==1) then
-      dz(:) = uniformResolution(ke, coord_mode, tmpReal, GV%Rlay(1), GV%Rlay(1))
+      dz(:) = uniformResolution(ke, coord_mode, tmpReal, US%R_to_kg_m3*GV%Rlay(1), US%R_to_kg_m3*GV%Rlay(1))
     else
       dz(:) = uniformResolution(ke, coord_mode, tmpReal, &
-                   GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2)), &
-                   GV%Rlay(ke)+0.5*(GV%Rlay(ke)-GV%Rlay(ke-1)) )
+                   US%R_to_kg_m3*(GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2))), &
+                   US%R_to_kg_m3*(GV%Rlay(ke)+0.5*(GV%Rlay(ke)-GV%Rlay(ke-1))) )
     endif
     if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
                    trim(message), units=trim(coord_units))
@@ -491,7 +491,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
 
   ! \todo This line looks like it would overwrite the target densities set just above?
   elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
-    call set_target_densities_from_GV(GV, CS)
+    call set_target_densities_from_GV(GV, US, CS)
     call log_param(param_file, mdl, "!TARGET_DENSITIES", CS%target_density, &
              'RHO target densities for interfaces', units=coordinateUnits(coord_mode))
   endif
@@ -1991,15 +1991,16 @@ subroutine setCoordinateResolution( dz, CS, scale )
 end subroutine setCoordinateResolution
 
 !> Set target densities based on the old Rlay variable
-subroutine set_target_densities_from_GV( GV, CS )
+subroutine set_target_densities_from_GV( GV, US, CS )
   type(verticalGrid_type), intent(in)    :: GV !< Ocean vertical grid structure
+  type(unit_scale_type),   intent(in)    :: US  !< A dimensional unit scaling type
   type(regridding_CS),     intent(inout) :: CS !< Regridding control structure
   ! Local variables
   integer :: k, nz
 
   nz = CS%nk
-  CS%target_density(1)    = GV%Rlay(1)+0.5*(GV%Rlay(1)-GV%Rlay(2))
-  CS%target_density(nz+1) = GV%Rlay(nz)+0.5*(GV%Rlay(nz)-GV%Rlay(nz-1))
+  CS%target_density(1)    = US%R_to_kg_m3*(GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(2)))
+  CS%target_density(nz+1) = US%R_to_kg_m3*(GV%Rlay(nz) + 0.5*(GV%Rlay(nz)-GV%Rlay(nz-1)))
   do k = 2,nz
     CS%target_density(k) = CS%target_density(k-1) + CS%coordinateResolution(k)
   enddo

src/core/MOM.F90

@@ -2795,7 +2795,7 @@ subroutine extract_surface_state(CS, sfc_state)
           sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * CS%tv%T(i,j,k)
           sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * CS%tv%S(i,j,k)
         else
-          sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * GV%Rlay(k)
+          sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * US%R_to_kg_m3*GV%Rlay(k)
         endif
         depth(i) = depth(i) + dh
       enddo ; enddo

src/core/MOM_PressureForce_Montgomery.F90

@@ -150,7 +150,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
 
   I_gEarth = 1.0 / (US%L_T_to_m_s**2 * GV%g_Earth)
   dp_neglect = GV%H_to_Pa * GV%H_subroundoff
-  do k=1,nz ; alpha_Lay(k) = 1.0 / GV%Rlay(k) ; enddo
+  do k=1,nz ; alpha_Lay(k) = 1.0 / (US%R_to_kg_m3*GV%Rlay(k)) ; enddo
   do k=2,nz ; dalpha_int(K) = alpha_Lay(k-1) - alpha_Lay(k) ; enddo
 
   if (use_p_atm) then
@@ -235,7 +235,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
         call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -491,7 +491,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
 
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -745,7 +745,7 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
   dP_dH = US%m_s_to_L_T**2*GV%H_to_Pa
   dp_neglect = GV%H_to_Pa * GV%H_subroundoff
 
-  do k=1,nz ; alpha_Lay(k) = 1.0 / GV%Rlay(k) ; enddo
+  do k=1,nz ; alpha_Lay(k) = 1.0 / (US%R_to_kg_m3*GV%Rlay(k)) ; enddo
   do k=2,nz ; dalpha_int(K) = alpha_Lay(k-1) - alpha_Lay(k) ; enddo
 
   if (use_EOS) then

src/core/MOM_PressureForce_analytic_FV.F90

@@ -231,7 +231,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
         call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -286,7 +286,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
                                useMassWghtInterp = CS%useMassWghtInterp)
       endif
     else
-      alpha_anom = 1.0/GV%Rlay(k) - alpha_ref
+      alpha_anom = 1.0/(US%R_to_kg_m3*GV%Rlay(k)) - alpha_ref
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dp(i,j) = GV%H_to_Pa * h(i,j,k)
         dza(i,j,k) = alpha_anom * dp(i,j)
@@ -349,7 +349,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dM(i,j) = (CS%GFS_scale - 1.0) * US%m_s_to_L_T**2 * &
-          (p(i,j,1)*(1.0/GV%Rlay(1) - alpha_ref) + za(i,j))
+          (p(i,j,1)*(1.0/(US%R_to_kg_m3*GV%Rlay(1)) - alpha_ref) + za(i,j))
       enddo ; enddo
     endif
 !  else
@@ -590,7 +590,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
 
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -622,7 +622,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
     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 * US%R_to_kg_m3*GV%Rlay(1)) * e(i,j,1)
       enddo ; enddo
     endif
   endif
@@ -702,16 +702,16 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dz(i,j) = g_Earth_z * GV%H_to_Z*h(i,j,k)
-        dpa(i,j) = (GV%Rlay(k) - rho_ref)*dz(i,j)
-        intz_dpa(i,j) = 0.5*(GV%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k)
+        dpa(i,j) = (US%R_to_kg_m3*GV%Rlay(k) - rho_ref)*dz(i,j)
+        intz_dpa(i,j) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref)*dz(i,j)*h(i,j,k)
       enddo ; enddo
       !$OMP parallel do default(shared)
       do j=js,je ; do I=Isq,Ieq
-        intx_dpa(I,j) = 0.5*(GV%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j))
+        intx_dpa(I,j) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref) * (dz(i,j)+dz(i+1,j))
       enddo ; enddo
       !$OMP parallel do default(shared)
       do J=Jsq,Jeq ; do i=is,ie
-        inty_dpa(i,J) = 0.5*(GV%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1))
+        inty_dpa(i,J) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref) * (dz(i,j)+dz(i,j+1))
       enddo ; enddo
     endif
 

src/core/MOM_PressureForce_blocked_AFV.F90

@@ -227,7 +227,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
         call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, &
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -251,7 +251,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
                                inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, &
                                useMassWghtInterp = CS%useMassWghtInterp)
     else
-      alpha_anom = 1.0/GV%Rlay(k) - alpha_ref
+      alpha_anom = 1.0/(US%R_to_kg_m3*GV%Rlay(k)) - alpha_ref
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dp(i,j) = GV%H_to_Pa * h(i,j,k)
         dza(i,j,k) = alpha_anom * dp(i,j)
@@ -314,7 +314,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dM(i,j) = (CS%GFS_scale - 1.0) * &
-          US%m_s_to_L_T**2*(p(i,j,1)*(1.0/GV%Rlay(1) - alpha_ref) + za(i,j))
+          US%m_s_to_L_T**2*(p(i,j,1)*(1.0/(US%R_to_kg_m3*GV%Rlay(1)) - alpha_ref) + za(i,j))
       enddo ; enddo
     endif
 !  else
@@ -574,7 +574,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
                         Rho_cv_BL(:), Isq, Ieq-Isq+2, tv%eqn_of_state)
 
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
-          if (GV%Rlay(k) < Rho_cv_BL(i)) then
+          if (US%R_to_kg_m3*GV%Rlay(k) < Rho_cv_BL(i)) then
             tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
           else
             tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
@@ -606,7 +606,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
     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 * US%R_to_kg_m3*GV%Rlay(1)) * e(i,j,1)
       enddo ; enddo
     endif
   endif
@@ -698,14 +698,14 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
         do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1
           i = ib+ioff_bk ; j = jb+joff_bk
           dz_bk(ib,jb) = g_Earth_z*GV%H_to_Z*h(i,j,k)
-          dpa_bk(ib,jb) = (GV%Rlay(k) - rho_ref)*dz_bk(ib,jb)
-          intz_dpa_bk(ib,jb) = 0.5*(GV%Rlay(k) - rho_ref)*dz_bk(ib,jb)*h(i,j,k)
+          dpa_bk(ib,jb) = (US%R_to_kg_m3*GV%Rlay(k) - rho_ref)*dz_bk(ib,jb)
+          intz_dpa_bk(ib,jb) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref)*dz_bk(ib,jb)*h(i,j,k)
         enddo ; enddo
         do jb=js_bk,je_bk ; do Ib=Isq_bk,Ieq_bk
-          intx_dpa_bk(Ib,jb) = 0.5*(GV%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib+1,jb))
+          intx_dpa_bk(Ib,jb) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib+1,jb))
         enddo ; enddo
         do Jb=Jsq_bk,Jeq_bk ; do ib=is_bk,ie_bk
-          inty_dpa_bk(ib,Jb) = 0.5*(GV%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib,jb+1))
+          inty_dpa_bk(ib,Jb) = 0.5*(US%R_to_kg_m3*GV%Rlay(k) - rho_ref) * (dz_bk(ib,jb)+dz_bk(ib,jb+1))
         enddo ; enddo
       endif
 

src/core/MOM_interface_heights.F90

@@ -116,7 +116,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
     else
 !$OMP do
       do j=jsv,jev ;  do k=nz,1,-1; do i=isv,iev
-        eta(i,j,K) = eta(i,j,K+1) + H_to_rho_eta*h(i,j,k)/GV%Rlay(k)
+        eta(i,j,K) = eta(i,j,K+1) + H_to_rho_eta*h(i,j,k) / (US%R_to_kg_m3*GV%Rlay(k))
       enddo ; enddo ; enddo
     endif
     if (present(eta_bt)) then
@@ -214,7 +214,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
     else
 !$OMP do
       do j=js,je ; do k=1,nz ; do i=is,ie
-        eta(i,j) = eta(i,j) + H_to_rho_eta*h(i,j,k)/GV%Rlay(k)
+        eta(i,j) = eta(i,j) + H_to_rho_eta*h(i,j,k) / (US%R_to_kg_m3*GV%Rlay(k))
       enddo ; enddo ; enddo
     endif
     if (present(eta_bt)) then

src/core/MOM_isopycnal_slopes.F90

@@ -166,7 +166,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
   do j=js,je ; do K=nz,2,-1
     if (.not.(use_EOS)) then
       drdiA = 0.0 ; drdiB = 0.0
-      drdkL = GV%Rlay(k)-GV%Rlay(k-1) ; drdkR = GV%Rlay(k)-GV%Rlay(k-1)
+      drdkL = US%R_to_kg_m3*(GV%Rlay(k)-GV%Rlay(k-1)) ; drdkR = US%R_to_kg_m3*(GV%Rlay(k)-GV%Rlay(k-1))
     endif
 
     ! Calculate the zonal isopycnal slope.
@@ -253,7 +253,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
   do j=js-1,je ; do K=nz,2,-1
     if (.not.(use_EOS)) then
       drdjA = 0.0 ; drdjB = 0.0
-      drdkL = GV%Rlay(k)-GV%Rlay(k-1) ; drdkR = GV%Rlay(k)-GV%Rlay(k-1)
+      drdkL = US%R_to_kg_m3*(GV%Rlay(k)-GV%Rlay(k-1)) ; drdkR = US%R_to_kg_m3*(GV%Rlay(k)-GV%Rlay(k-1))
     endif
 
     if (use_EOS) then

src/core/MOM_verticalGrid.F90

@@ -49,7 +49,7 @@ module MOM_verticalGrid
                         !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m.
   real, allocatable, dimension(:) :: &
     g_prime, &          !< The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2].
-    Rlay                !< The target coordinate value (potential density) in each layer [kg m-3].
+    Rlay                !< The target coordinate value (potential density) in each layer [R ~> kg m-3].
   integer :: nkml = 0   !< The number of layers at the top that should be treated
                         !! as parts of a homogeneous region.
   integer :: nk_rho_varies = 0 !< The number of layers at the top where the
@@ -272,23 +272,24 @@ function get_tr_flux_units(GV, tr_units, tr_vol_conc_units, tr_mass_conc_units)
 end function get_tr_flux_units
 
 !> This sets the coordinate data for the "layer mode" of the isopycnal model.
-subroutine setVerticalGridAxes( Rlay, GV )
+subroutine setVerticalGridAxes( Rlay, GV, scale )
   type(verticalGrid_type), intent(inout) :: GV   !< The container for vertical grid data
-  real, dimension(GV%ke),  intent(in)    :: Rlay !< The layer target density
+  real, dimension(GV%ke),  intent(in)    :: Rlay !< The layer target density [R ~> kg m-3]
+  real,                    intent(in)    :: scale !< A unit scaling factor for Rlay
   ! Local variables
   integer :: k, nk
 
   nk = GV%ke
 
   GV%zAxisLongName = 'Target Potential Density'
   GV%zAxisUnits = 'kg m-3'
-  do k=1,nk ; GV%sLayer(k) = Rlay(k) ; enddo
+  do k=1,nk ; GV%sLayer(k) = scale*Rlay(k) ; enddo
   if (nk > 1) then
-    GV%sInterface(1) = 1.5*Rlay(1) - 0.5*Rlay(2)
-    do K=2,nk ; GV%sInterface(K) = 0.5*( Rlay(k-1) + Rlay(k) ) ; enddo
-    GV%sInterface(nk+1) = 1.5*Rlay(nk) - 0.5*Rlay(nk-1)
+    GV%sInterface(1) = scale * (1.5*Rlay(1) - 0.5*Rlay(2))
+    do K=2,nk ; GV%sInterface(K) = scale * 0.5*( Rlay(k-1) + Rlay(k) ) ; enddo
+    GV%sInterface(nk+1) = scale * (1.5*Rlay(nk) - 0.5*Rlay(nk-1))
   else
-    GV%sInterface(1) = 0.0 ; GV%sInterface(nk+1) = 2.0*Rlay(nk)
+    GV%sInterface(1) = 0.0 ; GV%sInterface(nk+1) = 2.0*scale*Rlay(nk)
   endif
 
 end subroutine setVerticalGridAxes

src/diagnostics/MOM_diagnostics.F90

@@ -471,7 +471,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       enddo ; enddo
     else ! Rcv should not be used much in this case, so fill in sensible values.
       do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
-        Rcv(i,j,k) = GV%Rlay(k)
+        Rcv(i,j,k) = US%R_to_kg_m3*GV%Rlay(k)
       enddo ; enddo ; enddo
     endif
     if (CS%id_Rml > 0) call post_data(CS%id_Rml, Rcv, CS%diag)
@@ -489,7 +489,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
           CS%h_Rlay(i,j,k) = h(i,j,k)
         enddo ; enddo
         do k=1,nkmb ; do i=is,ie
-          call find_weights(GV%Rlay, Rcv(i,j,k), k_list, nz, wt, wt_p)
+          call find_weights(US%R_to_kg_m3*GV%Rlay, Rcv(i,j,k), k_list, nz, wt, wt_p)
           CS%h_Rlay(i,j,k_list)   = CS%h_Rlay(i,j,k_list)   + h(i,j,k)*wt
           CS%h_Rlay(i,j,k_list+1) = CS%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
         enddo ; enddo
@@ -511,7 +511,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
         enddo ; enddo
         k_list = nz/2
         do k=1,nkmb ; do I=Isq,Ieq
-          call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
+          call find_weights(US%R_to_kg_m3*GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
           CS%uh_Rlay(I,j,k_list)   = CS%uh_Rlay(I,j,k_list)   + uh(I,j,k)*wt
           CS%uh_Rlay(I,j,k_list+1) = CS%uh_Rlay(I,j,k_list+1) + uh(I,j,k)*wt_p
         enddo ; enddo
@@ -532,7 +532,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
           CS%vh_Rlay(i,J,k) = vh(i,J,k)
         enddo ; enddo
         do k=1,nkmb ; do i=is,ie
-          call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
+          call find_weights(US%R_to_kg_m3*GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
           CS%vh_Rlay(i,J,k_list)   = CS%vh_Rlay(i,J,k_list)   + vh(i,J,k)*wt
           CS%vh_Rlay(i,J,k_list+1) = CS%vh_Rlay(i,J,k_list+1) + vh(i,J,k)*wt_p
         enddo ; enddo
@@ -553,7 +553,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
           CS%uhGM_Rlay(I,j,k) = CDp%uhGM(I,j,k)
         enddo ; enddo
         do k=1,nkmb ; do I=Isq,Ieq
-          call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
+          call find_weights(US%R_to_kg_m3*GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i+1,j,k)), k_list, nz, wt, wt_p)
           CS%uhGM_Rlay(I,j,k_list)   = CS%uhGM_Rlay(I,j,k_list)   + CDp%uhGM(I,j,k)*wt
           CS%uhGM_Rlay(I,j,k_list+1) = CS%uhGM_Rlay(I,j,k_list+1) + CDp%uhGM(I,j,k)*wt_p
         enddo ; enddo
@@ -574,7 +574,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
           CS%vhGM_Rlay(i,J,k) = CDp%vhGM(i,J,k)
         enddo ; enddo
         do k=1,nkmb ; do i=is,ie
-          call find_weights(GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
+          call find_weights(US%R_to_kg_m3*GV%Rlay, 0.5*(Rcv(i,j,k)+Rcv(i,j+1,k)), k_list, nz, wt, wt_p)
           CS%vhGM_Rlay(i,J,k_list)   = CS%vhGM_Rlay(i,J,k_list)   + CDp%vhGM(i,J,k)*wt
           CS%vhGM_Rlay(i,J,k_list+1) = CS%vhGM_Rlay(i,J,k_list+1) + CDp%vhGM(i,J,k)*wt_p
         enddo ; enddo
@@ -850,7 +850,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
         enddo
       else
         do k=1,nz ; do j=js,je ; do i=is,ie
-          mass(i,j) = mass(i,j) + (GV%H_to_m*GV%Rlay(k))*h(i,j,k)
+          mass(i,j) = mass(i,j) + (GV%H_to_m*US%R_to_kg_m3*GV%Rlay(k))*h(i,j,k)
         enddo ; enddo ; enddo
       endif
     else

src/diagnostics/MOM_wave_speed.F90

@@ -192,10 +192,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
 
           ! Start a new layer
           H_here(i) = h(i,j,k)*GV%H_to_Z
-          HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k)
+          HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*US%R_to_kg_m3*GV%Rlay(k)
         else
           H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z
-          HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k)
+          HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*US%R_to_kg_m3*GV%Rlay(k)
         endif
       enddo ; enddo
       do i=is,ie ; if (H_here(i) > 0.0) then
@@ -649,10 +649,10 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
 
           ! Start a new layer
           H_here(i) = h(i,j,k)*GV%H_to_Z
-          HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k)
+          HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*US%R_to_kg_m3*GV%Rlay(k)
         else
           H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z
-          HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k)
+          HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*US%R_to_kg_m3*GV%Rlay(k)
         endif
       enddo ; enddo
       do i=is,ie ; if (H_here(i) > 0.0) then