Merge branch 'dev/gfdl' into nonBous_dumbbell_initialization

src/core/MOM_boundary_update.F90

@@ -156,7 +156,7 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time)
       call shelfwave_set_OBC_data(OBC, CS%shelfwave_OBC_CSp, G, GV, US, h, Time)
   if (CS%use_dyed_channel) &
       call dyed_channel_update_flow(OBC, CS%dyed_channel_OBC_CSp, G, GV, US, Time)
-  if (OBC%needs_IO_for_data .or. OBC%add_tide_constituents)  &
+  if (OBC%any_needs_IO_for_data .or. OBC%add_tide_constituents)  &
       call update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
 
 end subroutine update_OBC_data

src/core/MOM_dynamics_split_RK2.F90

@@ -50,11 +50,11 @@ module MOM_dynamics_split_RK2
 use MOM_hor_index,             only : hor_index_type
 use MOM_hor_visc,              only : horizontal_viscosity, hor_visc_CS
 use MOM_hor_visc,              only : hor_visc_init, hor_visc_end
-use MOM_interface_heights,     only : find_eta, thickness_to_dz
+use MOM_interface_heights,     only : thickness_to_dz, find_col_avg_SpV
 use MOM_lateral_mixing_coeffs, only : VarMix_CS
 use MOM_MEKE_types,            only : MEKE_type
 use MOM_open_boundary,         only : ocean_OBC_type, radiation_open_bdry_conds
-use MOM_open_boundary,         only : open_boundary_zero_normal_flow
+use MOM_open_boundary,         only : open_boundary_zero_normal_flow, open_boundary_query
 use MOM_open_boundary,         only : open_boundary_test_extern_h, update_OBC_ramp
 use MOM_PressureForce,         only : PressureForce, PressureForce_CS
 use MOM_PressureForce,         only : PressureForce_init
@@ -344,6 +344,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   real, dimension(SZI_(G),SZJ_(G)) :: eta_pred ! The predictor value of the free surface height
                                                ! or column mass [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G)) :: SpV_avg  ! The column averaged specific volume [R-1 ~> m3 kg-1]
   real, dimension(SZI_(G),SZJ_(G)) :: deta_dt  ! A diagnostic of the time derivative of the free surface
                                                ! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1]
 
@@ -596,6 +597,14 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   if (.not.BT_cont_BT_thick) &
     call btcalc(h, G, GV, CS%barotropic_CSp, OBC=CS%OBC)
   call bt_mass_source(h, eta, .true., G, GV, CS%barotropic_CSp)
+
+  SpV_avg(:,:) = 0.0
+  if ((.not.GV%Boussinesq) .and. associated(CS%OBC)) then
+    ! Determine the column average specific volume if it is needed due to the
+    ! use of Flather open boundary conditions in non-Boussinesq mode.
+    if (open_boundary_query(CS%OBC, apply_Flather_OBC=.true.)) &
+      call find_col_avg_SpV(h, SpV_avg, tv, G, GV, US)
+  endif
   call cpu_clock_end(id_clock_btcalc)
 
   if (G%nonblocking_updates) &
@@ -625,7 +634,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   ! The CS%ADp argument here stores the weights for certain integrated diagnostics.
   call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
               CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
-              CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, CS%ADp, CS%OBC, CS%BT_cont, &
+              CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
               eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
   if (showCallTree) call callTree_leave("btstep()")
   call cpu_clock_end(id_clock_btstep)
@@ -847,7 +856,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   ! This is the corrector step call to btstep.
   call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
               CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
-              CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, CS%ADp, CS%OBC, CS%BT_cont, &
+              CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
               eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
   if (CS%id_deta_dt>0) then
     do j=js,je ; do i=is,ie ; deta_dt(i,j) = (eta_pred(i,j) - eta(i,j))*Idt_bc ; enddo ; enddo

src/core/MOM_interface_heights.F90

@@ -19,7 +19,7 @@ module MOM_interface_heights
 
 public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
 public calc_derived_thermo
-public find_rho_bottom
+public find_rho_bottom, find_col_avg_SpV
 
 !> Calculates the heights of the free surface or all interfaces from layer thicknesses.
 interface find_eta
@@ -323,6 +323,74 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
 end subroutine calc_derived_thermo
 
 
+!> Determine the column average specific volumes.
+subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size)
+  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, 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(inout) :: SpV_avg !< Column average specific volume [R-1 ~> m3 kg-1]
+                                                  ! SpV_avg is intent inout to retain excess halo values.
+  type(thermo_var_ptrs),    intent(in)    :: tv   !< Structure containing pointers to any available
+                                                  !! thermodynamic fields.
+  integer,        optional, intent(in)    :: halo_size !< width of halo points on which to work
+
+  ! Local variables
+  real :: h_tot(SZI_(G))        ! Sum of the layer thicknesses [H ~> m or kg m-3]
+  real :: SpV_x_h_tot(SZI_(G))  ! Vertical sum of the layer average specific volume times
+                                ! the layer thicknesses [H R-1 ~> m4 kg-1 or m]
+  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
+
+  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
+  nz = GV%ke
+
+  if (GV%Boussinesq) then
+    I_rho = 1.0 / GV%Rho0
+    do j=js,je ; do i=is,ie
+      SpV_avg(i,j) = I_rho
+    enddo ; enddo
+  elseif (.not.allocated(tv%SpV_avg)) then
+    do k=1,nz ; Spv_lay(k) = 1.0 / GV%Rlay(k) ; enddo
+    do j=js,je
+      do i=is,ie ; SpV_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo
+      do k=1,nz ; do i=is,ie
+        h_tot(i) = h_tot(i) + max(h(i,j,k), GV%H_subroundoff)
+        SpV_x_h_tot(i) = SpV_x_h_tot(i) + Spv_lay(k)*max(h(i,j,k), GV%H_subroundoff)
+      enddo ; enddo
+      do i=is,ie ; SpV_avg(i,j) = SpV_x_h_tot(i) / h_tot(i) ; enddo
+    enddo
+  else
+    ! Check that SpV_avg has been set.
+    if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < halo)) then
+      if (tv%valid_SpV_halo < 0) then
+        mesg = "invalid values of SpV_avg."
+      else
+        write(mesg, '("insufficiently large SpV_avg halos of width ", i2, " but ", i2," is needed.")') &
+                     tv%valid_SpV_halo, halo
+      endif
+      call MOM_error(FATAL, "find_col_avg_SpV called in fully non-Boussinesq mode with "//trim(mesg))
+    endif
+
+    do j=js,je
+      do i=is,ie ; SpV_x_h_tot(i) = 0.0 ; h_tot(i) = 0.0 ; enddo
+      do k=1,nz ; do i=is,ie
+        h_tot(i) = h_tot(i) + max(h(i,j,k), GV%H_subroundoff)
+        SpV_x_h_tot(i) = SpV_x_h_tot(i) + tv%SpV_avg(i,j,k)*max(h(i,j,k), GV%H_subroundoff)
+      enddo ; enddo
+      do i=is,ie ; SpV_avg(i,j) = SpV_x_h_tot(i) / h_tot(i) ; enddo
+    enddo
+  endif
+
+end subroutine find_col_avg_SpV
+
+
 !> 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.
 subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)