Merge branch 'Hallberg-NOAA-rescale_pressure' into rescale_pressure

- Resolved conflicts

src/ALE/MOM_regridding.F90

@@ -1889,8 +1889,7 @@ subroutine convective_adjustment(G, GV, h, tv)
   do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
 
     ! Compute densities within current water column
-    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, &
-                            densities, 1, GV%ke, tv%eqn_of_state )
+    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state )
 
     ! Repeat restratification until complete
     do
@@ -1909,8 +1908,7 @@ subroutine convective_adjustment(G, GV, h, tv)
           tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0
           h(i,j,k)    = h1 ; h(i,j,k+1)    = h0
           ! Recompute densities at levels k and k+1
-          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
-                                  densities(k), tv%eqn_of_state )
+          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state)
           call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
                                   densities(k+1), tv%eqn_of_state )
           stratified = .false.

src/ALE/coord_adapt.F90

@@ -156,7 +156,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
          0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * GV%H_to_Pa, &
-         alpha, beta, 2, nz - 1, tv%eqn_of_state)
+         alpha, beta, tv%eqn_of_state, dom=(/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + &
@@ -168,7 +168,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
          0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * GV%H_to_Pa, &
-         alpha, beta, 2, nz - 1, tv%eqn_of_state)
+         alpha, beta, tv%eqn_of_state, dom=(/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + &
@@ -180,7 +180,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
          0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * GV%H_to_Pa, &
-         alpha, beta, 2, nz - 1, tv%eqn_of_state)
+         alpha, beta, tv%eqn_of_state, dom=(/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + &
@@ -192,7 +192,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
          0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), &
          0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), &
          0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * GV%H_to_Pa, &
-         alpha, beta, 2, nz - 1, tv%eqn_of_state)
+         alpha, beta, tv%eqn_of_state, dom=(/2,nz/))
 
     del2sigma(2:nz) = del2sigma(2:nz) + &
          (alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + &
@@ -206,7 +206,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
   ! a positive curvature means we're too light relative to adjacent columns,
   ! so del2sigma needs to be positive too (push the interface deeper)
   call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * GV%H_to_Pa, &
-       alpha, beta, 1, nz + 1, tv%eqn_of_state)
+       alpha, beta, tv%eqn_of_state, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below.
   do K = 2, nz
     ! TODO make lower bound here configurable
     del2sigma(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &

src/ALE/coord_hycom.F90

@@ -131,7 +131,7 @@ subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
   z_scale = 1.0 ; if (present(zScale)) z_scale = zScale
 
   ! Work bottom recording potential density
-  call calculate_density(T, S, p_col, rho_col, 1, nz, eqn_of_state)
+  call calculate_density(T, S, p_col, rho_col, eqn_of_state)
   ! This ensures the potential density profile is monotonic
   ! although not necessarily single valued.
   do k = nz-1, 1, -1

src/ALE/coord_rho.F90

@@ -124,7 +124,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
 
     ! Compute densities on source column
     pres(:) = CS%ref_pressure
-    call calculate_density(T, S, pres, densities, 1, nz, eqn_of_state)
+    call calculate_density(T, S, pres, densities, eqn_of_state)
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))
     enddo

src/ALE/coord_slight.F90

@@ -251,7 +251,7 @@ subroutine build_slight_column(CS, eqn_of_state, H_to_pres, H_subroundoff, &
     dz = (z_col(nz+1) - z_col(1)) / real(nz)
     do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo
   else
-    call calculate_density(T_col, S_col, p_col, rho_col, 1, nz, eqn_of_state)
+    call calculate_density(T_col, S_col, p_col, rho_col, eqn_of_state)
 
     ! Find the locations of the target potential densities, flagging
     ! locations in apparently unstable regions as not reliable.
@@ -372,10 +372,10 @@ subroutine build_slight_column(CS, eqn_of_state, H_to_pres, H_subroundoff, &
       enddo
       T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz)
       p_IS(nz+1) = z_col(nz+1) * H_to_pres
-      call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, 2, nz-1, &
-                                    eqn_of_state)
-      call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, &
-                                    eqn_of_state)
+      call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, &
+                                    eqn_of_state, dom=(/2,nz/))
+      call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, &
+                                    eqn_of_state, dom=(/2,nz/))
       if (CS%compressibility_fraction > 0.0) then
         call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state)
       else

src/core/MOM.F90

@@ -74,7 +74,7 @@ module MOM
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
 use MOM_dyn_horgrid,           only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
-use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze
+use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze, EOS_domain
 use MOM_fixed_initialization,  only : MOM_initialize_fixed
 use MOM_forcing_type,          only : allocate_forcing_type, allocate_mech_forcing
 use MOM_forcing_type,          only : deallocate_mech_forcing, deallocate_forcing_type
@@ -2904,8 +2904,8 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
     ! 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(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, G%HI, &
-                               tv%eqn_of_state)
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, &
+                               tv%eqn_of_state, dom=EOS_domain(G%HI))
         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

src/core/MOM_PressureForce_Montgomery.F90

@@ -124,12 +124,13 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
   real :: alpha_Lay(SZK_(G)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
   real :: dalpha_int(SZK_(G)+1) ! The change in specific volume across each
                              ! interface [R-1 ~> m3 kg-1].
-  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
+  integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
+  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer :: i, j, k
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
   nkmb=GV%nk_rho_varies
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
-  start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2
+  EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
   use_p_atm = .false.
   if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
@@ -228,8 +229,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
         do k=1,nkmb ; do i=Isq,Ieq+1
           tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
         enddo ; enddo
-        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
-                               tv%eqn_of_state)
+        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
+                               tv%eqn_of_state, dom=EOSdom)
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
           if (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)
@@ -245,8 +246,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
     endif
     !$OMP parallel do default(shared) private(rho_in_situ)
     do k=1,nz ; do j=Jsq,Jeq+1
-      call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, start, npts, &
-                             tv%eqn_of_state)
+      call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, &
+                             tv%eqn_of_state, dom=EOSdom)
       do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo
     enddo ; enddo
   endif                                               ! use_EOS
@@ -409,13 +410,14 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
                              ! gradient terms are to be split into
                              ! barotropic and baroclinic pieces.
   type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
-  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
+  integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
+  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer :: i, j, k
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
   nkmb=GV%nk_rho_varies
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
-  start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2
+  EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
   use_p_atm = .false.
   if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
@@ -484,8 +486,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
         do k=1,nkmb ; do i=Isq,Ieq+1
           tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
         enddo ; enddo
-        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
-                               tv%eqn_of_state)
+        call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
+                               tv%eqn_of_state, dom=EOSdom)
 
         do k=nkmb+1,nz ; do i=Isq,Ieq+1
           if (GV%Rlay(k) < Rho_cv_BL(i)) then
@@ -505,8 +507,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
     ! will come down with a fatal error if there is any compressibility.
     !$OMP parallel do default(shared)
     do k=1,nz+1 ; do j=Jsq,Jeq+1
-      call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), start, npts, &
-                             tv%eqn_of_state)
+      call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
+                             tv%eqn_of_state, dom=EOSdom)
       do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo
     enddo ; enddo
   endif                                               ! use_EOS
@@ -632,10 +634,11 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
                              ! an equation of state.
   real :: z_neglect          ! A thickness that is so small it is usually lost
                              ! in roundoff and can be neglected [Z ~> m].
-  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts
+  integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
+  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
 
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
-  start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2
+  EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
   Rho0xG = Rho0 * GV%g_Earth
   G_Rho0 = GV%g_Earth / GV%Rho0
@@ -662,8 +665,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
           Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
           press(i) = -Rho0xG*e(i,j,1)
         enddo
-        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, start, npts, &
-                               tv%eqn_of_state)
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
+                               tv%eqn_of_state, dom=EOSdom)
         do i=Isq,Ieq+1
           pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
         enddo
@@ -673,8 +676,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
             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
-          call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, start, npts, &
-                                        tv%eqn_of_state)
+          call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, &
+                                        tv%eqn_of_state, dom=EOSdom)
           do i=Isq,Ieq+1
             pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * &
                ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * &
@@ -733,10 +736,11 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
   real :: dp_neglect         ! A thickness that is so small it is usually lost
                              ! in roundoff and can be neglected [R L2 T-2 ~> Pa].
   logical :: use_EOS         ! If true, density is calculated from T & S using an equation of state.
-  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts
+  integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
+  integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
 
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
-  start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2
+  EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
   use_EOS = associated(tv%eqn_of_state)
 
@@ -759,8 +763,8 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
     else
       !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
       do j=Jsq,Jeq+1
-        call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, start, npts, &
-                               tv%eqn_of_state)
+        call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, &
+                               tv%eqn_of_state, dom=EOSdom)
         do i=Isq,Ieq+1
           C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
           pbce(i,j,nz) = dP_dH / (rho_in_situ(i))
@@ -770,9 +774,9 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
             T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
             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, start, npts, tv%eqn_of_state)
-          call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, start, npts, &
-                                        tv%eqn_of_state)
+          call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, dom=EOSdom)
+          call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
+                                        tv%eqn_of_state, dom=EOSdom)
           do i=Isq,Ieq+1
             pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) *  &
                 ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &