Pass rescaled pressures to calculate_TFreeze

  Use the new pres_scale argument to TFreeze and pass rescaled pressures to
calculate_TFreeze and several instances of calculate_density.  All answers are
bitwise identical.

src/ice_shelf/MOM_ice_shelf.F90

@@ -158,7 +158,8 @@ module MOM_ice_shelf
   logical :: find_salt_root              !< If true, if true find Sbdry using a quadratic eq.
   real    :: TFr_0_0                     !< The freezing point at 0 pressure and 0 salinity [degC]
   real    :: dTFr_dS                     !< Partial derivative of freezing temperature with salinity [degC ppt-1]
-  real    :: dTFr_dp                     !< Partial derivative of freezing temperature with pressure [degC Pa-1]
+  real    :: dTFr_dp                     !< Partial derivative of freezing temperature with
+                                         !! pressure [degC T2 R-1 L-2 ~> degC Pa-1]
   !>@{ Diagnostic handles
   integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, &
              id_tfreeze = -1, id_tfl_shelf = -1, &
@@ -217,7 +218,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
                !< with temperature [kg m-3 degC-1].
     dR0_dS, &  !< Partial derivative of the mixed layer density
                !< with salinity [kg m-3 ppt-1].
-    p_int      !< The pressure at the ice-ocean interface [Pa].
+    p_int      !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa].
 
   real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
     exch_vel_t, &  !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1]
@@ -371,13 +372,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
   do j=js,je
     ! Find the pressure at the ice-ocean interface, averaged only over the
     ! part of the cell covered by ice shelf.
-    do i=is,ie ; p_int(i) = US%RL2_T2_to_Pa*CS%g_Earth * ISS%mass_shelf(i,j) ; enddo
+    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)
+             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)
+             dR0_dT, dR0_dS, is, ie-is+1, CS%eqn_of_state, pres_scale=US%RL2_T2_to_Pa)
 
     do i=is,ie
       if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. &
@@ -445,7 +446,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
 
           do it1 = 1,20
             ! Determine the potential temperature at the ice-ocean interface.
-            call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state)
+            call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, &
+                                   pres_scale=US%RL2_T2_to_Pa)
 
             dT_ustar = (ISS%tfreeze(i,j) - state%sst(i,j)) * ustar_h
             dS_ustar = (Sbdry(i,j) - state%sss(i,j)) * ustar_h
@@ -588,7 +590,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
           ! is specified and large enough that the ocean salinity at the interface
           ! is about the same as the boundary layer salinity.
 
-          call calculate_TFreeze(state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state)
+          call calculate_TFreeze(state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, &
+                                 pres_scale=US%RL2_T2_to_Pa)
 
           exch_vel_t(i,j) = CS%gamma_t
           ISS%tflux_ocn(i,j) = RhoCp * exch_vel_t(i,j) * (ISS%tfreeze(i,j) - state%sst(i,j))
@@ -1272,12 +1275,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl
                  "this is the freezing potential temperature at "//&
                  "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.)
     call get_param(param_file, mdl, "DTFREEZE_DS", CS%dTFr_dS, &
-                 "this is the derivative of the freezing potential "//&
-                 "temperature with salinity.", units="degC psu-1", default=-0.054, do_not_log=.true.)
+                 "this is the derivative of the freezing potential temperature with salinity.", &
+                 units="degC psu-1", default=-0.054, do_not_log=.true.)
     call get_param(param_file, mdl, "DTFREEZE_DP", CS%dTFr_dp, &
-                 "this is the derivative of the freezing potential "//&
-                 "temperature with pressure.", &
-                 units="degC Pa-1", default=0.0, do_not_log=.true.)
+                 "this is the derivative of the freezing potential temperature with pressure.", &
+                 units="degC Pa-1", default=0.0, scale=US%RL2_T2_to_Pa, do_not_log=.true.)
   endif
 
   call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, &

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -118,9 +118,10 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
   real, dimension(SZI_(G)) :: &
     fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
     T_freeze, & ! The freezing potential temperature at the current salinity [degC].
-    ps          ! pressure
+    ps          ! Surface pressure [R L2 T-2 ~> Pa]
   real, dimension(SZI_(G),SZK_(G)) :: &
-    pressure    ! The pressure at the middle of each layer [Pa].
+    pressure    ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
+  real :: H_to_RL2_T2  ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
   real :: hc    ! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1].
   logical :: T_fr_set  ! True if the freezing point has been calculated for a
                        ! row of points.
@@ -135,25 +136,27 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
 
   if (.not.CS%pressure_dependent_frazil) then
     do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
+  else
+    H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth
   endif
 !$OMP parallel do default(none) shared(is,ie,js,je,CS,G,GV,US,h,nz,tv,p_surf) &
 !$OMP                           private(fraz_col,T_fr_set,T_freeze,hc,ps)  &
 !$OMP                      firstprivate(pressure)    !pressure might be set above, so should be firstprivate
   do j=js,je
     ps(:) = 0.0
     if (PRESENT(p_surf)) then ; do i=is,ie
-      ps(i) = US%RL2_T2_to_Pa*p_surf(i,j)
+      ps(i) = p_surf(i,j)
     enddo ; endif
 
     do i=is,ie ; fraz_col(i) = 0.0 ; enddo
 
     if (CS%pressure_dependent_frazil) then
       do i=is,ie
-        pressure(i,1) = ps(i) + (0.5*GV%H_to_Pa)*h(i,j,1)
+        pressure(i,1) = ps(i) + (0.5*H_to_RL2_T2)*h(i,j,1)
       enddo
       do k=2,nz ; do i=is,ie
         pressure(i,k) = pressure(i,k-1) + &
-          (0.5*GV%H_to_Pa) * (h(i,j,k) + h(i,j,k-1))
+          (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1))
       enddo ; enddo
     endif
 
@@ -162,7 +165,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
       do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
         if (.not.T_fr_set) then
           call calculate_TFreeze(tv%S(i:,j,1), pressure(i:,1), T_freeze(i:), &
-                                 1, ie-i+1, tv%eqn_of_state)
+                                 1, ie-i+1, tv%eqn_of_state, pres_scale=US%RL2_T2_to_Pa)
           T_fr_set = .true.
         endif
 
@@ -188,7 +191,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
             ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
           if (.not.T_fr_set) then
             call calculate_TFreeze(tv%S(i:,j,k), pressure(i:,k), T_freeze(i:), &
-                                   1, ie-i+1, tv%eqn_of_state)
+                                   1, ie-i+1, tv%eqn_of_state, pres_scale=US%RL2_T2_to_Pa)
             T_fr_set = .true.
           endif
 

src/parameterizations/vertical/MOM_entrain_diffusive.F90

@@ -174,7 +174,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
   real :: g_2dt     ! 0.5 * G_Earth / dt, times unit conversion factors
                     ! [m3 H-2 s-2 T-1 ~> m s-3 or m7 kg-2 s-3].
   real, dimension(SZI_(G)) :: &
-    pressure, &      ! The pressure at an interface [Pa].
+    pressure, &      ! The pressure at an interface [R L2 T-2 ~> Pa].
     T_eos, S_eos, &  ! The potential temperature and salinity at which to
                      ! evaluate dRho_dT and dRho_dS [degC] and [ppt].
     dRho_dT, dRho_dS ! The partial derivatives of potential density with temperature and
@@ -836,12 +836,12 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
       do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo
       if (associated(tv%eqn_of_state)) then
         if (associated(fluxes%p_surf)) then
-          do i=is,ie ; pressure(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo
+          do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo
         else
           do i=is,ie ; pressure(i) = 0.0 ; enddo
         endif
         do K=2,nz
-          do i=is,ie ; pressure(i) = pressure(i) + GV%H_to_Pa*h(i,j,k-1) ; enddo
+          do i=is,ie ; pressure(i) = pressure(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) ; enddo
           do i=is,ie
             if (k==kb(i)) then
               T_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
@@ -851,8 +851,8 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
               S_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
             endif
           enddo
-          call calculate_density_derivs(T_eos, S_eos, pressure, &
-                  dRho_dT, dRho_dS, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+          call calculate_density_derivs(T_eos, S_eos, pressure, dRho_dT, dRho_dS, is, ie-is+1, &
+                             tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
           do i=is,ie
             if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,K) = 0.0
             else

src/parameterizations/vertical/MOM_full_convection.F90

@@ -345,7 +345,7 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h
   real :: c1(SZI_(G),SZK_(G))      ! tridiagonal solver.
   real :: T_f(SZI_(G),SZK_(G))     ! Filtered temperatures [degC]
   real :: S_f(SZI_(G),SZK_(G))     ! Filtered salinities [ppt]
-  real :: pres(SZI_(G))            ! Interface pressures [Pa].
+  real :: pres(SZI_(G))            ! Interface pressures [R L2 T-2 ~> Pa].
   real :: T_EOS(SZI_(G))           ! Filtered and vertically averaged temperatures [degC]
   real :: S_EOS(SZI_(G))           ! Filtered and vertically averaged salinities [ppt]
   real :: kap_dt_x2                ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4].
@@ -403,24 +403,24 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h
   endif
 
   if (associated(p_surf)) then
-    do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*p_surf(i,j) ; enddo
+    do i=is,ie ; pres(i) = p_surf(i,j) ; enddo
   else
     do i=is,ie ; pres(i) = 0.0 ; enddo
   endif
-  call calculate_density_derivs(T_f(:,1), S_f(:,1), pres, dR_dT(:,1), dR_dS(:,1), &
-                                is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
-  do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*GV%H_to_Pa ; enddo
+  call calculate_density_derivs(T_f(:,1), S_f(:,1), pres, dR_dT(:,1), dR_dS(:,1), is-G%isd+1, ie-is+1, &
+                                tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
+  do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*(GV%H_to_RZ*GV%g_Earth) ; enddo
   do K=2,nz
     do i=is,ie
       T_EOS(i) = 0.5*(T_f(i,k-1) + T_f(i,k))
       S_EOS(i) = 0.5*(S_f(i,k-1) + S_f(i,k))
     enddo
-    call calculate_density_derivs(T_EOS, S_EOS, pres, dR_dT(:,K), dR_dS(:,K), &
-                                  is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
-    do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*GV%H_to_Pa ; enddo
+    call calculate_density_derivs(T_EOS, S_EOS, pres, dR_dT(:,K), dR_dS(:,K), is-G%isd+1, ie-is+1, &
+                                  tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
+    do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*(GV%H_to_RZ*GV%g_Earth) ; enddo
   enddo
-  call calculate_density_derivs(T_f(:,nz), S_f(:,nz), pres, dR_dT(:,nz+1), dR_dS(:,nz+1), &
-                                is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+  call calculate_density_derivs(T_f(:,nz), S_f(:,nz), pres, dR_dT(:,nz+1), dR_dS(:,nz+1), is-G%isd+1, &
+                       ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
 
 
 end subroutine smoothed_dRdT_dRdS

src/parameterizations/vertical/MOM_internal_tide_input.F90

@@ -167,7 +167,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot)
   real, dimension(SZI_(G),SZK_(G)+1) :: &
     dRho_int      ! The unfiltered density differences across interfaces [R ~> kg m-3].
   real, dimension(SZI_(G)) :: &
-    pres, &       ! The pressure at each interface [Pa].
+    pres, &       ! The pressure at each interface [R L2 T-2 ~> Pa].
     Temp_int, &   ! The temperature at each interface [degC].
     Salin_int, &  ! The salinity at each interface [ppt].
     drho_bot, &   ! The density difference at the bottom of a layer [R ~> kg m-3]
@@ -199,18 +199,18 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot)
   do j=js,je
     if (associated(tv%eqn_of_state)) then
       if (associated(fluxes%p_surf)) then
-        do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo
+        do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
       else
         do i=is,ie ; pres(i) = 0.0 ; enddo
       endif
       do K=2,nz
         do i=is,ie
-          pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1)
+          pres(i) = pres(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1)
           Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1))
           Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1))
         enddo
-        call calculate_density_derivs(Temp_int, Salin_int, pres, &
-                 dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+        call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:), dRho_dS(:), &
+                     is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
         do i=is,ie
           dRho_int(i,K) = max(dRho_dT(i)*(T_f(i,j,k) - T_f(i,j,k-1)) + &
                               dRho_dS(i)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0)

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -868,7 +868,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
     dRho_dS            ! partial derivative of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
 
   real, dimension(SZI_(G)) :: &
-    pres,      &  ! pressure at each interface [Pa]
+    pres,      &  ! pressure at each interface [R L2 T-2 ~> Pa]
     Temp_int,  &  ! temperature at each interface [degC]
     Salin_int, &  ! salinity at each interface [ppt]
     drho_bot,  &  ! A density difference [R ~> kg m-3]
@@ -896,18 +896,18 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
   enddo
   if (associated(tv%eqn_of_state)) then
     if (associated(fluxes%p_surf)) then
-      do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo
+      do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
     else
       do i=is,ie ; pres(i) = 0.0 ; enddo
     endif
     do K=2,nz
       do i=is,ie
-        pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1)
+        pres(i) = pres(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1)
         Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1))
         Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1))
       enddo
-      call calculate_density_derivs(Temp_int, Salin_int, pres, &
-               dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R)
+      call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:,K), dRho_dS(:,K), &
+                   is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
       do i=is,ie
         dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + &
                             dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0)