Add an option to calculate bottom pressure from initialization for SAL

src/core/MOM.F90

@@ -3553,14 +3553,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
       endif
     endif
   elseif (CS%use_RK2) then
-    call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, Time, G, GV, US,     &
-            param_file, diag, CS%dyn_unsplit_RK2_CSp,                      &
+    call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, CS%tv, Time, G, GV,  &
+            US, param_file, diag, CS%dyn_unsplit_RK2_CSp,                  &
             CS%ADp, CS%CDp, MOM_internal_state, CS%OBC,                    &
             CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, &
             CS%ntrunc, cont_stencil=CS%cont_stencil)
   else
-    call initialize_dyn_unsplit(CS%u, CS%v, CS%h, Time, G, GV, US,         &
-            param_file, diag, CS%dyn_unsplit_CSp,                          &
+    call initialize_dyn_unsplit(CS%u, CS%v, CS%h, CS%tv, Time, G, GV,      &
+            US, param_file, diag, CS%dyn_unsplit_CSp,                      &
             CS%ADp, CS%CDp, MOM_internal_state, CS%OBC,                    &
             CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, &
             CS%ntrunc, cont_stencil=CS%cont_stencil)

src/core/MOM_dynamics_split_RK2.F90

@@ -1552,7 +1552,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
+  if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
   if (CS%use_tides) then
     call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp)
     HA_CSp => CS%HA_CSp

src/core/MOM_dynamics_split_RK2b.F90

@@ -1435,7 +1435,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
+  if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
   if (CS%use_tides) then
     call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp)
     HA_CSp => CS%HA_CSp

src/core/MOM_dynamics_unsplit.F90

@@ -576,7 +576,7 @@ subroutine register_restarts_dyn_unsplit(HI, GV, param_file, CS)
 end subroutine register_restarts_dyn_unsplit
 
 !> Initialize parameters and allocate memory associated with the unsplit dynamics module.
-subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS, &
+subroutine initialize_dyn_unsplit(u, v, h, tv, Time, G, GV, US, param_file, diag, CS, &
                                   Accel_diag, Cont_diag, MIS, &
                                   OBC, update_OBC_CSp, ALE_CSp, set_visc, &
                                   visc, dirs, ntrunc, cont_stencil)
@@ -589,6 +589,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS
                                   intent(inout) :: v          !< The meridional velocity [L T-1 ~> m s-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                   intent(inout) :: h          !< Layer thicknesses [H ~> m or kg m-2]
+  type(thermo_var_ptrs),          intent(in)    :: tv         !< Thermodynamic type
   type(time_type),        target, intent(in)    :: Time       !< The current model time.
   type(param_file_type),          intent(in)    :: param_file !< A structure to parse
                                                               !! for run-time parameters.
@@ -710,7 +711,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
+  if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp)
   if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, CS%ADp, &
                           CS%SAL_CSp, CS%tides_CSp)

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -526,7 +526,7 @@ subroutine register_restarts_dyn_unsplit_RK2(HI, GV, param_file, CS)
 end subroutine register_restarts_dyn_unsplit_RK2
 
 !> Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module.
-subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag, CS, &
+subroutine initialize_dyn_unsplit_RK2(u, v, h, tv, Time, G, GV, US, param_file, diag, CS, &
                                       Accel_diag, Cont_diag, MIS, &
                                       OBC, update_OBC_CSp, ALE_CSp, set_visc, &
                                       visc, dirs, ntrunc, cont_stencil)
@@ -536,6 +536,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u   !< The zonal velocity [L T-1 ~> m s-1].
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v   !< The meridional velocity [L T-1 ~> m s-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h    !< Layer thicknesses [H ~> m or kg m-2]
+  type(thermo_var_ptrs),                     intent(in)    :: tv   !< Thermodynamic type
   type(time_type),                   target, intent(in)    :: Time !< The current model time.
   type(param_file_type),                     intent(in)    :: param_file !< A structure to parse
                                                                          !! for run-time parameters.
@@ -673,7 +674,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp)
+  if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp)
   if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, CS%ADp, &
                           CS%SAL_CSp, CS%tides_CSp)

src/core/MOM_interface_heights.F90

@@ -3,7 +3,7 @@ module MOM_interface_heights
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
-use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol
+use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol, int_density_dz
 use MOM_debugging,     only : hchksum
 use MOM_error_handler, only : MOM_error, FATAL
 use MOM_EOS,           only : calculate_density, average_specific_vol, EOS_type, EOS_domain
@@ -20,7 +20,7 @@ module MOM_interface_heights
 public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
 public calc_derived_thermo
 public convert_MLD_to_ML_thickness
-public find_rho_bottom, find_col_avg_SpV
+public find_rho_bottom, find_col_avg_SpV, find_col_mass
 
 !> Calculates the heights of the free surface or all interfaces from layer thicknesses.
 interface find_eta
@@ -73,7 +73,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
                             ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
   real :: dZ_ref    ! The difference in the reference height between G%bathyT and eta [Z ~> m].
                     ! dZ_ref is 0 unless the optional argument dZref is present.
-  integer i, j, k, isv, iev, jsv, jev, nz, halo
+  integer :: i, j, k, isv, iev, jsv, jev, nz, halo
 
   halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
 
@@ -191,7 +191,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
                             ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1]
   real :: dZ_ref    ! The difference in the reference height between G%bathyT and eta [Z ~> m].
                     ! dZ_ref is 0 unless the optional argument dZref is present.
-  integer i, j, k, is, ie, js, je, nz, halo
+  integer :: i, j, k, is, ie, js, je, nz, halo
 
   halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
   is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo
@@ -345,7 +345,7 @@ subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size)
   real :: I_rho                 ! The inverse of the Boussiensq reference density [R-1 ~> m3 kg-1]
   real :: SpV_lay(SZK_(GV))     ! The inverse of the layer target potential densities [R-1 ~> m3 kg-1]
   character(len=128) :: mesg    ! A string for error messages
-  integer i, j, k, is, ie, js, je, nz, halo
+  integer :: i, j, k, is, ie, js, je, nz, halo
 
   halo = 0 ; if (present(halo_size)) halo = max(0,halo_size)
 
@@ -391,6 +391,71 @@ subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size)
 
 end subroutine find_col_avg_SpV
 
+!> Calculate the integrated mass of the water column.
+subroutine find_col_mass(h, tv, G, GV, US, mass, p_bot, p_surf)
+  type(ocean_grid_type),                      intent(in)  :: G    !< The ocean's grid structure.
+  type(verticalGrid_type),                    intent(in)  :: GV   !< The ocean's vertical grid structure.
+  type(unit_scale_type),                      intent(in)  :: US   !< A dimensional unit scaling type
+  type(thermo_var_ptrs),                      intent(in)  :: tv   !< A structure pointing to various
+                                                                  !! thermodynamic variables.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: mass !< Integrated mass of the water column
+                                                                  !! [R Z ~> kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: p_bot  !< Bottom pressure = g * mass + psurf
+                                                                    !! [R L2 T-2 ~> Pa]
+  real, dimension(:,:),             optional, pointer     :: p_surf !< A pointer to surface pressure
+                                                                    !! [R L2 T-2 ~> Pa]
+
+  ! Local variables
+  real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1]
+  real, dimension(SZI_(G),SZJ_(G)) :: &
+    z_top, & ! Height of the top of a layer [Z ~> m].
+    z_bot, & ! Height of the bottom of a layer [Z ~> m].
+    dp       ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa].
+  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+  isq = G%iscB ; ieq = G%iecB ; jsq = G%jscB ; jeq = G%jecB
+  nz = GV%ke
+
+  do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
+  if (GV%Boussinesq) then
+    if (associated(tv%eqn_of_state)) then
+      I_gEarth = 1.0 / GV%g_Earth
+      do j=jsq,jeq+1 ; do i=isq,ieq+1 ; z_bot(i,j) = 0.0 ; enddo ; enddo
+      do k=1,nz
+        ! NOTE: int_density_z expects z_top and z_bot values from [ij]sq to [ij]eq+1
+        do j=jsq,jeq+1 ; do i=isq,ieq+1
+          z_top(i,j) = z_bot(i,j)
+          z_bot(i,j) = z_top(i,j) - GV%H_to_Z * h(i,j,k)
+        enddo ; enddo
+        call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, &
+                            G%HI, tv%eqn_of_state, US, dp)
+        do j=js,je ; do i=is,ie
+          mass(i,j) = mass(i,j) + dp(i,j) * I_gEarth
+        enddo ; enddo
+      enddo
+    else
+      do k=1,nz ; do j=js,je ; do i=is,ie
+        mass(i,j) = mass(i,j) + (GV%H_to_Z * GV%Rlay(k)) * h(i,j,k)
+      enddo ; enddo ; enddo
+    endif
+  else
+    do k=1,nz ; do j=js,je ; do i=is,ie
+      mass(i,j) = mass(i,j) + GV%H_to_RZ * h(i,j,k)
+    enddo ; enddo ; enddo
+  endif
+
+  if (present(p_bot)) then
+    do j=js,je ; do i=is,ie
+      p_bot(i,j) = GV%g_Earth * mass(i,j)
+    enddo ; enddo
+    if (present(p_surf) .and. associated(p_surf)) then ; do j=js,je ; do i=is,ie
+      p_bot(i,j) = p_bot(i,j) + p_surf(i,j)
+    enddo ; enddo ; endif
+  endif
+
+end subroutine find_col_mass
 
 !> Determine the in situ density averaged over a specified distance from the bottom,
 !! calculating it as the inverse of the mass-weighted average specific volume.

src/diagnostics/MOM_diagnostics.F90

@@ -24,7 +24,7 @@ module MOM_diagnostics
 use MOM_error_handler,     only : MOM_error, FATAL, WARNING
 use MOM_file_parser,       only : get_param, log_version, param_file_type
 use MOM_grid,              only : ocean_grid_type
-use MOM_interface_heights, only : find_eta
+use MOM_interface_heights, only : find_eta, find_col_mass
 use MOM_spatial_means,     only : global_area_mean, global_layer_mean
 use MOM_spatial_means,     only : global_volume_mean, global_area_integral
 use MOM_tracer_registry,   only : tracer_registry_type, post_tracer_transport_diagnostics
@@ -39,7 +39,6 @@ module MOM_diagnostics
 #include <MOM_memory.h>
 
 public calculate_diagnostic_fields, register_time_deriv, write_static_fields
-public find_eta
 public register_surface_diags, post_surface_dyn_diags, post_surface_thermo_diags
 public register_transport_diags, post_transport_diagnostics
 public MOM_diagnostics_init, MOM_diagnostics_end
@@ -903,7 +902,6 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
               ! This is the column mass multiplied by gravity plus the pressure
               ! at the ocean surface [R L2 T-2 ~> Pa].
-    dpress, & ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa].
     tr_int    ! vertical integral of a tracer times density,
               ! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2].
   real    :: IG_Earth  ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1].
@@ -943,52 +941,14 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     call post_data(CS%id_col_ht, z_bot, CS%diag)
   endif
 
-  ! NOTE: int_density_z expects z_top and z_btm values from [ij]sq to [ij]eq+1
   if (CS%id_col_mass > 0 .or. CS%id_pbo > 0) then
-    do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
-    if (GV%Boussinesq) then
-      if (associated(tv%eqn_of_state)) then
-        IG_Earth = 1.0 / GV%g_Earth
-        do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1
-          z_bot(i,j) = 0.0
-        enddo ; enddo
-        do k=1,nz
-          do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1
-            z_top(i,j) = z_bot(i,j)
-            z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k)
-          enddo ; enddo
-          call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, &
-                              G%HI, tv%eqn_of_state, US, dpress)
-          do j=js,je ; do i=is,ie
-            mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth
-          enddo ; enddo
-        enddo
-      else
-        do k=1,nz ; do j=js,je ; do i=is,ie
-          mass(i,j) = mass(i,j) + (GV%H_to_Z*GV%Rlay(k))*h(i,j,k)
-        enddo ; enddo ; enddo
-      endif
-    else
-      do k=1,nz ; do j=js,je ; do i=is,ie
-        mass(i,j) = mass(i,j) + GV%H_to_RZ*h(i,j,k)
-      enddo ; enddo ; enddo
-    endif
-    if (CS%id_col_mass > 0) then
-      call post_data(CS%id_col_mass, mass, CS%diag)
-    endif
     if (CS%id_pbo > 0) then
-      do j=js,je ; do i=is,ie ; btm_pres(i,j) = 0.0 ; enddo ; enddo
-      ! 'pbo' is defined as the sea water pressure at the sea floor
-      !     pbo = (mass * g) + p_surf
-      ! where p_surf is the sea water pressure at sea water surface.
-      do j=js,je ; do i=is,ie
-        btm_pres(i,j) = GV%g_Earth * mass(i,j)
-        if (associated(p_surf)) then
-          btm_pres(i,j) = btm_pres(i,j) + p_surf(i,j)
-        endif
-      enddo ; enddo
+      call find_col_mass(h, tv, G, GV, US, mass, btm_pres, p_surf)
       call post_data(CS%id_pbo, btm_pres, CS%diag)
+    else
+      call find_col_mass(h, tv, G, GV, US, mass)
     endif
+    if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag)
   endif
 
 end subroutine calculate_vertical_integrals