Baroclinic sea level calculation

src/core/MOM.F90

@@ -101,7 +101,7 @@ module MOM
 use MOM_harmonic_analysis,     only : HA_accum, harmonic_analysis_CS
 use MOM_hor_index,             only : hor_index_type, hor_index_init
 use MOM_hor_index,             only : rotate_hor_index
-use MOM_interface_heights,     only : find_eta, calc_derived_thermo, thickness_to_dz
+use MOM_interface_heights,     only : find_eta, find_bsl, calc_derived_thermo, thickness_to_dz
 use MOM_interface_filter,      only : interface_filter, interface_filter_init, interface_filter_end
 use MOM_interface_filter,      only : interface_filter_CS
 use MOM_internal_tides,        only : int_tide_CS
@@ -332,6 +332,7 @@ module MOM
                                          !! segment data with every call to update_OBC_segment_data.
   type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period.
   type(time_type) :: dt_obc_seg_time     !< The next time OBC segment update is applied to OBGC tracers.
+  real :: rho_bsl                      !< Surface density for calculating the baroclinic sea level [R ~> kg m-3]
 
   real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied
   !! by ice shelf [nondim]
@@ -409,6 +410,8 @@ module MOM
     !< Pointer to the control structure used for an alternate version of the mode-split RK2 dynamics
   type(harmonic_analysis_CS),    pointer :: HA_CSp => NULL()
     !< Pointer to the control structure for harmonic analysis
+  logical                                :: HA_bsl
+    !< If true, perform harmonic analysis of baroclinic sea level
   type(thickness_diffuse_CS) :: thickness_diffuse_CSp
     !< Pointer to the control structure used for the isopycnal height diffusive transport.
     !! This is also common referred to as Gent-McWilliams diffusion
@@ -597,6 +600,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
                 ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
   real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
     ssh         ! sea surface height, which may be based on eta_av [Z ~> m]
+  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
+    bsl         ! Baroclinic sea level [Z ~> m]
   real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: &
     dz          ! Vertical distance across layers [Z ~> m]
 
@@ -1057,6 +1062,10 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
                           G, GV, US, CS%diagnostics_CSp)
       call post_tracer_diagnostics_at_sync(CS%Tracer_reg, h, CS%diag_pre_sync, CS%diag, G, GV, CS%t_dyn_rel_diag)
       call diag_copy_diag_to_storage(CS%diag_pre_sync, h, CS%diag)
+      if (associated(CS%HA_CSp) .and. CS%HA_bsl) then
+        call find_bsl(h, CS%tv, G, GV, US, CS%rho_bsl, bsl, dZref=G%Z_ref)
+        call HA_accum('bsl', bsl, Time_local, G, CS%HA_CSp)
+      endif
       if (showCallTree) call callTree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
       call disable_averaging(CS%diag)
       CS%t_dyn_rel_diag = 0.0
@@ -3576,6 +3585,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     CS%dt_obc_seg_time = Time + CS%dt_obc_seg_interval
   endif
 
+  if (associated(CS%HA_CSp)) then
+    call get_param(param_file, "MOM", "HA_BSL",CS%HA_bsl, &
+                   "If true, perform harmonic analysis of baroclinic sea level.", &
+                   default=.false., do_not_log=.true.)
+    call get_param(param_file, "MOM", "RHO_BSL", CS%rho_bsl, &
+                   "Surface density for baroclinic sea level calculation.", &
+                   units='kg m-3', default=1025.0, scale=US%kg_m3_to_R, do_not_log=.not.CS%HA_bsl)
+  else
+    CS%HA_bsl = .false.
+  endif
+
   call callTree_waypoint("dynamics initialized (initialize_MOM)")
 
   CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, &

src/core/MOM_interface_heights.F90

@@ -19,7 +19,7 @@ module MOM_interface_heights
 
 #include <MOM_memory.h>
 
-public find_eta, find_dz_for_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
+public find_eta, find_bsl, find_dz_for_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, find_col_mass
@@ -290,6 +290,105 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref)
 
 end subroutine find_eta_2d
 
+!> Calculates the baroclinic sea level, following Xu et al., to be submitted to JPO
+subroutine find_bsl(h, tv, G, GV, US, rho_s, bsl, dZref)
+  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
+  real,                                       intent(in)  :: rho_s !< Surface density [R ~> kg m-3]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(in)  :: h   !< Layer thicknesses [H ~> m or kg m-2]
+  type(thermo_var_ptrs),                      intent(in)  :: tv  !< A structure pointing to various
+                                                                 !! thermodynamic variables
+  real, dimension(SZI_(G),SZJ_(G)),           intent(out) :: bsl !< Baroclinic sea level [Z ~> m]
+  real,                             optional, intent(in)  :: dZref !< The difference in the
+                    !! reference height between G%bathyT and eta [Z ~> m]. The default is 0
+
+  ! Local variables
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta             ! layer interface heights [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G)) :: &
+    bathyT, &       ! Bathymetry at T points plus dZ_ref [Z ~> m]
+    pt, &           ! Pressure at the top of a layer [R L2 T-2 ~> Pa]
+    pb, &           ! Pressure at the bottom of a layer [R L2 T-2 ~> Pa]
+    gz, &           ! Geopotential at the bottom of a layer [L2 T-2 ~> m2 s-2]
+    dp, &           ! Pressure change across a layer in Boussinesq mode [R L2 T-2 ~> Pa]
+    dg, &           ! Geopotential change across a layer in non-Boussinesq mode [L2 T-2 ~> m2 s-2]
+    dp_int, &       ! Layer-integrated pressure change in Boussinesq mode [R L2 Z T-2 ~> Pa m]
+    dg_int, &       ! Layer-integrated geopotential change in non-Boussinesq mode [R L4 T-4 ~> Pa m2 s-2]
+    p_int           ! Vertical integral of pressure at the bottom of a layer [R L2 Z T-2 ~> Pa m]
+                    ! or that scaled by GV%g_Earth in non-Boussinesq and EOS mode [R L4 T-4 ~> Pa m2 s-2]
+                    ! or that normalized by GV%g_Earth in non-EOS mode [R Z2 ~> Pa s2]
+  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
+  real :: I_gEarth  ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
+  real :: SpV_s     ! Specific volume of the surface layer [R-1 ~> m3 kg-1]
+  logical, dimension(SZI_(G),SZJ_(G)) :: maskT ! Mask at T points for skipping land points in calculations
+  integer :: i, j, k, is, ie, js, je, nz
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+  dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref
+
+  I_gEarth = 1.0 / GV%g_Earth; SpV_s = 1.0 / rho_s
+
+  call find_eta(h, tv, G, GV, US, eta, halo_size=1, dZref=dZ_ref)
+
+  !$OMP parallel default(shared)
+  !$OMP do
+  do j=js,je ; do i=is,ie
+    pt(i,j) = 0.0 ; pb(i,j) = 0.0 ; gz(i,j) = 0.0 ; p_int(i,j) = 0.0 ; bsl(i,j) = 0.0
+    bathyT(i,j) = G%bathyT(i,j) + dZ_ref
+    maskT(i,j) = G%mask2dT(i,j) > 0.0 .and. bathyT(i,j) > 0.0
+  enddo ; enddo
+
+  if (associated(tv%eqn_of_state)) then
+    if (GV%Boussinesq) then
+      do k=1,nz
+        call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), eta(:,:,k), eta(:,:,k+1), rho_s, &
+                            GV%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dp, dp_int)
+        !$OMP do
+        do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+          p_int(i,j) = p_int(i,j) + (pt(i,j) * (GV%H_to_Z * h(i,j,k)) + dp_int(i,j))
+          pt(i,j) = pt(i,j) + dp(i,j)
+        endif ; enddo ; enddo
+      enddo
+      !$OMP do
+      do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+        bsl(i,j) = - (p_int(i,j) * I_gEarth) / (rho_s * bathyT(i,j))
+      endif ; enddo ; enddo
+    else ! (.not. GV%Boussinesq)
+      do k=1,nz
+        !$OMP do
+        do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+          dp(i,j) = GV%g_Earth * (GV%H_to_RZ * h(i,j,k))
+          pb(i,j) = pt(i,j) + dp(i,j)
+        endif ; enddo ; enddo
+        call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), pt, pb, SpV_s, G%HI, &
+                                 tv%eqn_of_state, US, dg, dg_int)
+        !$OMP do
+        do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+          gz(i,j) = gz(i,j) + dg(i,j)
+          p_int(i,j) = p_int(i,j) + (gz(i,j) * dp(i,j) + dg_int(i,j))
+          pt(i,j) = pb(i,j)
+        endif ; enddo ; enddo
+      enddo
+      !$OMP do
+      do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+        bsl(i,j) = (p_int(i,j) * (I_gEarth * I_gEarth)) / (rho_s * bathyT(i,j))
+      endif ;enddo ; enddo
+    endif ! (GV%Boussinesq)
+  else ! (.not. associated(tv%eqn_of_state))
+    !$OMP do
+    do j=js,je ; do i=is,ie ; if (maskT(i,j)) then
+      do k=2,nz
+        p_int(i,j) = p_int(i,j) + (GV%Rlay(k) - GV%Rlay(k-1)) * &
+                                  ((eta(i,j,k) + bathyT(i,j)) * (eta(i,j,k) + bathyT(i,j)))
+      enddo
+      bsl(i,j) = - 0.5 * (p_int(i,j) / (rho_s * bathyT(i,j)))
+    endif ; enddo ; enddo
+  endif ! (associated(tv%eqn_of_state))
+  !$OMP end parallel
+
+end subroutine find_bsl
 
 !> Calculate derived thermodynamic quantities for re-use later.
 subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)

src/diagnostics/MOM_diagnostics.F90

@@ -27,7 +27,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, find_dz_for_eta, find_col_mass
+use MOM_interface_heights, only : find_eta, find_bsl, find_dz_for_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
@@ -64,6 +64,7 @@ module MOM_diagnostics
                                        !! non-Boussinesq layer thicknesses as are used to find the free
                                        !! surface height, instead of using an approximate thickness
                                        !! based on division by the mid-layer density.
+  real :: rho_bsl                      !< Surface density for calculating the baroclinic sea level [R ~> kg m-3]
 
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                        !! regulate the timing of diagnostic output.
@@ -126,6 +127,7 @@ module MOM_diagnostics
   integer :: id_h_pre_sync     = -1
   integer :: id_tosq           = -1, id_sosq           = -1
   integer :: id_t20d           = -1, id_t17d           = -1
+  integer :: id_bsl            = -1
 
   !>@}
   type(wave_speed_CS) :: wave_speed  !< Wave speed control struct
@@ -920,6 +922,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     z_top, &  ! Height of the top of a layer or the ocean [Z ~> m].
     z_bot, &  ! Height of the bottom of a layer (for id_mass) or the
               ! (positive) depth of the ocean (for id_col_ht) [Z ~> m].
+    bsl, &    ! Baroclinic sea level [Z ~> m]
     mass, &   ! integrated mass of the water column [R Z ~> kg m-2].  For
               ! non-Boussinesq models this is rho*dz. For Boussinesq
               ! models, this is either the integral of in-situ density
@@ -983,6 +986,12 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
     endif
     if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag)
   endif
+
+  if (CS%id_bsl > 0) then
+    call find_bsl(h, tv, G, GV, US, CS%rho_bsl, bsl, dZref=G%Z_ref)
+    call post_data(CS%id_bsl, bsl, CS%diag)
+  endif
+
   if (CS%id_t20d > 0 .or. CS%id_t17d > 0) then
     call PPM%init(GV%ke, h_neglect=0.)
     do j=js,je ; do i=is,ie
@@ -2294,6 +2303,12 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   CS%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, Time, &
       long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
       units='Pa', conversion=US%RL2_T2_to_Pa)
+  CS%id_bsl = register_diag_field('ocean_model', 'bsl', diag%axesT1, Time, &
+      long_name='Baroclinic Sea Level', standard_name='baroclinic_sea_level', &
+      units='m', conversion=US%Z_to_m)
+  call get_param(param_file, mdl, "RHO_BSL", CS%rho_bsl, &
+                 "Surface density for baroclinic sea level calculation.", &
+                 units='kg m-3', default=1025.0, scale=US%kg_m3_to_R, do_not_log=CS%id_bsl<=0)
 
   ! Register time derivatives and allocate memory for diagnostics that need
   ! access from across several modules.

src/diagnostics/MOM_harmonic_analysis.F90

@@ -93,7 +93,7 @@ subroutine HA_init(Time, US, param_file, nc, CS)
   type(HA_type) :: ha1                              !< A temporary, null field used for initializing CS%list
   real :: HA_start_time                             !< Start time of harmonic analysis [T ~> s]
   real :: HA_end_time                               !< End time of harmonic analysis [T ~> s]
-  logical :: HA_ssh, HA_ubt, HA_vbt
+  logical :: HA_ssh, HA_bsl, HA_ubt, HA_vbt
   character(len=40)  :: mdl="MOM_harmonic_analysis" !< This module's name
   character(len=255) :: mesg
   integer :: year, month, day, hour, minute, second
@@ -268,8 +268,14 @@ subroutine HA_init(Time, US, param_file, nc, CS)
 
   ! Register variables/fields to be analyzed
   call get_param(param_file, mdl, "HA_SSH", HA_ssh, &
-                 "If true, perform harmonic analysis of sea serface height.", default=.false.)
+                 "If true, perform harmonic analysis of sea surface height.", default=.false.)
   if (HA_ssh) call HA_register('ssh', 'h', CS)
+  call get_param(param_file, mdl, "HA_BSL", HA_bsl, &
+                 "If true, perform harmonic analysis of baroclinic sea level.", default=.false.)
+  if (HA_bsl) then
+    call HA_register('bsl', 'h', CS)
+  endif
+
   call get_param(param_file, mdl, "HA_UBT", HA_ubt, &
                  "If true, perform harmonic analysis of zonal barotropic velocity.", default=.false.)
   if (HA_ubt) call HA_register('ubt', 'u', CS)