Use new HI_1d forms of calculate_density calls

  Revised numerous calls to calculate_density and calculate_density_derivs to
use the new form with domain extents indicated by a hor_index_type argument.
Internal density variables were also rescaled in a few cases.  All answers are
bitwise identical.

src/core/MOM.F90

@@ -2685,11 +2685,11 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
   type(verticalGrid_type),           intent(in)    :: GV  !< ocean vertical grid structure
   type(unit_scale_type),             intent(in)    :: US  !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: ssh !< time mean surface height [m]
-  real, dimension(:,:),    optional, pointer       :: p_atm !< atmospheric pressure [R L2 T-2 ~> Pa]
+  real, dimension(:,:),    optional, pointer       :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa]
   logical,                 optional, intent(in)    :: use_EOS !< If true, calculate the density for
                                                        !! the SSH correction using the equation of state.
 
-  real :: Rho_conv    ! The density used to convert surface pressure to
+  real :: Rho_conv(SZI_(G))  ! The density used to convert surface pressure to
                       ! a corrected effective SSH [R ~> kg m-3].
   real :: IgR0        ! The SSH conversion factor from R L2 T-2 to m [m T2 R-1 L-2 ~> m Pa-1].
   logical :: calc_rho
@@ -2699,18 +2699,21 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
   if (present(p_atm)) then ; if (associated(p_atm)) then
     calc_rho = associated(tv%eqn_of_state)
     if (present(use_EOS) .and. calc_rho) calc_rho = use_EOS
-    ! Correct the output sea surface height for the contribution from the
-    ! atmospheric pressure
-    do j=js,je ; do i=is,ie
+    ! Correct the output sea surface height for the contribution from the ice pressure.
+    do j=js,je
       if (calc_rho) then
-        call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, Rho_conv, &
-                               tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, G%HI, &
+                               tv%eqn_of_state, US)
+        do i=is,ie
+          IgR0 = US%Z_to_m / (Rho_conv(i) * GV%g_Earth)
+          ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
+        enddo
       else
-        Rho_conv = GV%Rho0
+        do i=is,ie
+          ssh(i,j) = ssh(i,j) + p_atm(i,j) * (US%Z_to_m / (GV%Rho0 * GV%g_Earth))
+        enddo
       endif
-      IgR0 = US%Z_to_m / (Rho_conv * GV%g_Earth)
-      ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
-    enddo ; enddo
+    enddo
   endif ; endif
 
 end subroutine adjust_ssh_for_p_atm

src/core/MOM_PressureForce_Montgomery.F90

@@ -766,9 +766,6 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
             S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
           enddo
           call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, G%HI, tv%eqn_of_state, US, halo=1)
-!          call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
-!                                 Isq, Ieq-Isq+2, tv%eqn_of_state, &
-!                                 scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
           call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, G%HI, &
                                         tv%eqn_of_state, US, halo=1)
           do i=Isq,Ieq+1

src/core/MOM_PressureForce_blocked_AFV.F90

@@ -301,8 +301,8 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
     if (use_EOS) then
       !$OMP parallel do default(shared) private(rho_in_situ)
       do j=Jsq,Jeq+1
-        call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, Isq, Ieq-Isq+2, &
-                               tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
+        call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, G%HI, &
+                               tv%eqn_of_state, US, halo=1)
 
         do i=Isq,Ieq+1
           dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
@@ -586,11 +586,11 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1
         if (use_p_atm) then
-          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, Isq, &
-                       Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
+          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, G%HI, &
+                                 tv%eqn_of_state, US, halo=1)
         else
-          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &
-                                 Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
+          call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, G%HI, &
+                                 tv%eqn_of_state, US, halo=1)
         endif
         do i=Isq,Ieq+1
           dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)

src/core/MOM_isopycnal_slopes.F90

@@ -175,9 +175,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
         T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1)))
         S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
       enddo
-      call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, &
-                   drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, scale=US%kg_m3_to_R, &
-                   pres_scale=US%RL2_T2_to_Pa)
+      call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, (is-IsdB+1)-1, ie-is+2, &
+                   tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
     endif
 
     do I=is-1,ie
@@ -262,9 +261,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
         T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
         S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
       enddo
-      call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, &
-                   drho_dS_v, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, &
-                   pres_scale=US%RL2_T2_to_Pa)
+      call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, is, ie-is+1, &
+                   tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
     endif
     do i=is,ie
       if (use_EOS) then

src/ice_shelf/MOM_ice_shelf.F90

@@ -197,27 +197,27 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
                                                 !! describe the surface state of the ocean.  The
                                                 !! intent is only inout to allow for halo updates.
   type(forcing),         intent(inout) :: fluxes !< structure containing pointers to any possible
-                                                 !! thermodynamic or mass-flux forcing fields.
+                                                !! thermodynamic or mass-flux forcing fields.
   type(time_type),       intent(in)    :: Time  !< Start time of the fluxes.
-  real,                  intent(in)    :: time_step !< Length of time over which
-                                                    !! these fluxes will be applied [s].
-  type(ice_shelf_CS),    pointer       :: CS !< A pointer to the control structure
-                                             !! returned by a previous call to
-                                             !! initialize_ice_shelf.
+  real,                  intent(in)    :: time_step !< Length of time over which these fluxes
+                                                !! will be applied [s].
+  type(ice_shelf_CS),    pointer       :: CS    !< A pointer to the control structure returned
+                                                !! by a previous call to initialize_ice_shelf.
   type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
 
-  type(ocean_grid_type), pointer :: G => NULL() ! The grid structure used by the ice shelf.
-  type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing
-                                                 ! various unit conversion factors
+  ! Local variables
+  type(ocean_grid_type), pointer :: G => NULL()  !< The grid structure used by the ice shelf.
+  type(unit_scale_type), pointer :: US => NULL() !< Pointer to a structure containing
+                                                 !! various unit conversion factors
   type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe
-                                          !! the ice-shelf state
+                                                 !! the ice-shelf state
 
   real, dimension(SZI_(CS%grid)) :: &
-    Rhoml, &   !< Ocean mixed layer density [kg m-3].
+    Rhoml, &   !< Ocean mixed layer density [R ~> kg m-3].
     dR0_dT, &  !< Partial derivative of the mixed layer density
-               !< with temperature [kg m-3 degC-1].
+               !< with temperature [R degC-1 ~> kg m-3 degC-1].
     dR0_dS, &  !< Partial derivative of the mixed layer density
-               !< with salinity [kg m-3 ppt-1].
+               !< with salinity [R ppt-1 ~> kg m-3 ppt-1].
     p_int      !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa].
 
   real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
@@ -235,8 +235,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
                !! viscosity is linearly increasing [nondim]. (Was 1/8. Why?)
   real, parameter :: RC    = 0.20     ! critical flux Richardson number.
   real :: I_ZETA_N !< The inverse of ZETA_N [nondim].
-  real :: I_LF     !< The inverse of the latent Heat of fusion [Q-1 ~> kg J-1].
-  real :: I_VK     !< The inverse of VK.
+  real :: I_LF     !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1].
+  real :: I_VK     !< The inverse of the Von Karman constant [nondim].
   real :: PR, SC   !< The Prandtl number and Schmidt number [nondim].
 
   ! 3 equations formulation variables
@@ -263,7 +263,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
                    ! boundary layer salinity times the friction velocity [ppt Z T-1 ~> ppt m s-1]
   real :: ustar_h  ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1]
   real :: Gam_turb ! [nondim]
-  real :: Gam_mol_t, Gam_mol_s
+  real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivites [nondim]
   real :: RhoCp     ! A typical ocean density times the heat capacity of water [Q R ~> J m-3]
   real :: ln_neut
   real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1]
@@ -306,8 +306,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
   RhoCp = CS%Rho_ocn * CS%Cp
 
   !first calculate molecular component
-  Gam_mol_t = 12.5 * (PR**c2_3) - 6
-  Gam_mol_s = 12.5 * (SC**c2_3) - 6
+  Gam_mol_t = 12.5 * (PR**c2_3) - 6.0
+  Gam_mol_s = 12.5 * (SC**c2_3) - 6.0
 
   ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS)
   ! these fields are already set to zero during initialization
@@ -375,10 +375,10 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
     do i=is,ie ; p_int(i) = CS%g_Earth * ISS%mass_shelf(i,j) ; enddo
 
     ! Calculate insitu densities and expansion coefficients
-    call calculate_density(state%sst(:,j), state%sss(:,j), p_int, &
-             Rhoml(:), is, ie-is+1, CS%eqn_of_state, pres_scale=US%RL2_T2_to_Pa)
-    call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, &
-             dR0_dT, dR0_dS, is, ie-is+1, CS%eqn_of_state, pres_scale=US%RL2_T2_to_Pa)
+    call calculate_density(state%sst(:,j), state%sss(:,j), p_int, Rhoml(:), G%HI, &
+                           CS%eqn_of_state, US)
+    call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, dR0_dT, dR0_dS, G%HI, &
+                                  CS%eqn_of_state, US)
 
     do i=is,ie
       if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. &

src/initialization/MOM_state_initialization.F90

@@ -2185,7 +2185,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
   allocate(area_shelf_h(isd:ied,jsd:jed))
   allocate(frac_shelf_h(isd:ied,jsd:jed))
 
-  press(:) = tv%p_ref
+  press(:) = tv%P_Ref
 
   ! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
   call convert_temp_salt_for_TEOS10(temp_z, salt_z, press, G, kd, mask_z, eos)
@@ -2399,8 +2399,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
 
   if (adjust_temperature .and. .not. useALEremapping) then
     call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
-            GV%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, US, eos)
-
+            GV%Rlay(1:nz), tv%P_Ref, niter, missing_value, h(is:ie,js:je,:), ks, US, eos)
   endif
 
   deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -291,9 +291,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
     Net_salt, & ! The surface salt flux into the ocean over a time step, ppt H.
     Idecay_len_TKE, &  ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
     p_ref, &    !   Reference pressure for the potential density governing mixed
-                ! layer dynamics, almost always 0 (or 1e5) Pa.
+                ! layer dynamics, almost always 0 (or 1e5) [R L2 T-2 ~> Pa].
     p_ref_cv, & !   Reference pressure for the potential density which defines
-                ! the coordinate variable, set to P_Ref [Pa].
+                ! the coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
     dR0_dT, &   !   Partial derivative of the mixed layer potential density with
                 ! temperature [R degC-1 ~> kg m-3 degC-1].
     dRcv_dT, &  !   Partial derivative of the coordinate variable potential
@@ -376,7 +376,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
   Idt_diag = 1.0 / (dt__diag)
   write_diags = .true. ; if (present(last_call)) write_diags = last_call
 
-  p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref
+  p_ref(:) = 0.0 ; p_ref_cv(:) = US%kg_m3_to_R*US%m_s_to_L_T**2*tv%P_Ref
 
   nsw = CS%nsw
 
@@ -464,17 +464,14 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
     ! Calculate an estimate of the mid-mixed layer pressure [Pa]
     do i=is,ie ; p_ref(i) = 0.0 ; enddo
     do k=1,CS%nkml ; do i=is,ie
-      p_ref(i) = p_ref(i) + 0.5*GV%H_to_Pa*h(i,k)
+      p_ref(i) = p_ref(i) + 0.5*(GV%H_to_RZ*GV%g_Earth)*h(i,k)
     enddo ; enddo
-    call calculate_density_derivs(T(:,1), S(:,1), p_ref, dR0_dT, dR0_dS, &
-                                  is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
-    call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, &
-                                  is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+    call calculate_density_derivs(T(:,1), S(:,1), p_ref, dR0_dT, dR0_dS, G%HI, tv%eqn_of_state, US)
+    call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, G%HI, &
+                                  tv%eqn_of_state, US)
     do k=1,nz
-      call calculate_density(T(:,k), S(:,k), p_ref, R0(:,k), is, ie-is+1, &
-                             tv%eqn_of_state, scale=US%kg_m3_to_R)
-      call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), is, &
-                             ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+      call calculate_density(T(:,k), S(:,k), p_ref, R0(:,k), G%HI, tv%eqn_of_state, US)
+      call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state, US)
     enddo
     if (id_clock_EOS>0) call cpu_clock_end(id_clock_EOS)
 

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -408,11 +408,11 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
   real :: dzbr(SZI_(G)) ! Cumulative depth over which brine is distributed [H ~> m to kg m-2]
   real :: inject_layer(SZI_(G),SZJ_(G)) ! diagnostic
 
-  real :: p_ref_cv(SZI_(G))
+  real :: p_ref_cv(SZI_(G))       ! The pressure used to calculate the coordinate density [R L2 T-2 ~> Pa]
   real :: T(SZI_(G),SZK_(G))
   real :: S(SZI_(G),SZK_(G))
   real :: h_2d(SZI_(G),SZK_(G))   ! A 2-d slice of h with a minimum thickness [H ~> m to kg m-2]
-  real :: Rcv(SZI_(G),SZK_(G))
+  real :: Rcv(SZI_(G),SZK_(G))    ! The coordinate density [R ~> kg m-3]
   real :: s_new,R_new,t0,scale, cdz
   integer :: i, j, k, is, ie, js, je, nz, ks
 
@@ -427,7 +427,7 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
   ! because it is not convergent when resolution becomes very fine. I think that this whole
   ! subroutine needs to be revisited.- RWH
 
-  p_ref_cv(:)  = tv%P_ref
+  p_ref_cv(:) = US%kg_m3_to_R*US%m_s_to_L_T**2*tv%P_Ref
   brine_dz = 1.0*GV%m_to_H
 
   inject_layer(:,:) = nz
@@ -447,8 +447,7 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
         h_2d(i,k) = MAX(h(i,j,k), GV%Angstrom_H)
       enddo
 
-      call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), is, &
-                             ie-is+1, tv%eqn_of_state)
+      call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state, US)
     enddo
 
     ! First, try to find an interior layer where inserting all the salt
@@ -459,7 +458,7 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay)
       if ((G%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
         s_new = S(i,k) + salt(i) / (GV%H_to_RZ * h_2d(i,k))
         t0 = T(i,k)
-        call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state)
+        call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state, scale=US%kg_m3_to_R)
         if (R_new < 0.5*(Rcv(i,k)+Rcv(i,k+1)) .and. s_new<s_max) then
           dzbr(i) = dzbr(i)+h_2d(i,k)
           inject_layer(i,j) = min(inject_layer(i,j),real(k))

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -1979,9 +1979,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   integer :: kb(SZI_(G),SZJ_(G)) ! index of the lightest layer denser
                                  ! than the buffer layer [nondim]
 
-  real :: p_ref_cv(SZI_(G))      ! Reference pressure for the potential
-                                 ! density which defines the coordinate
-                                 ! variable, set to P_Ref [Pa].
+  real :: p_ref_cv(SZI_(G))      ! Reference pressure for the potential density that defines the
+                                 ! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
 
   logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below,
                                   ! where massive is defined as sufficiently thick that
@@ -2681,11 +2680,11 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
     call cpu_clock_begin(id_clock_sponge)
     ! Layer mode sponge
     if (CS%bulkmixedlayer .and. associated(tv%eqn_of_state)) then
-      do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo
+      do i=is,ie ; p_ref_cv(i) = US%kg_m3_to_R*US%m_s_to_L_T**2*tv%P_Ref ; enddo
       !$OMP parallel do default(shared)
       do j=js,je
-         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), &
-                             is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), G%HI, &
+                                tv%eqn_of_state, US)
       enddo
       call apply_sponge(h, dt, G, GV, US, ea, eb, CS%sponge_CSp, Rcv_ml)
     else