*Adjust Segment Data Thicknesses

  - previously, the code was treating the segment data thicknesses
    as given, without adjustments for model bathymetry.  This
    change starts with the reported dz and compresses to an Angstrom
    below the level of bathymetry in the case where the segment thicknesses
    are in excess. If the segments are deficient, then layers are uniformly
    expanded. This procedure follows model dz adjustments upon initilization.

src/core/MOM_open_boundary.F90

@@ -3229,6 +3229,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
             endif
           endif
 
+          call adjustSegmentEtaToFitBathymetry(G,GV,US,segment,m)
+
           if (segment%is_E_or_W) then
             ishift=1
             if (segment%direction == OBC_DIRECTION_E) ishift=0
@@ -4051,6 +4053,121 @@ subroutine open_boundary_register_restarts(HI, GV, OBC_CS,restart_CSp)
 
 end subroutine open_boundary_register_restarts
 
+!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
+!!
+!! If the bottom most interface is below the topography then the bottom-most
+!! layers are contracted to GV%Angstrom_m.
+!! If the bottom most interface is above the topography then the entire column
+!! is dilated (expanded) to fill the void.
+!!   @remark{There is a (hard-wired) "tolerance" parameter such that the
+!! criteria for adjustment must equal or exceed 10cm.}
+subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld)
+  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(OBC_segment_type),                       intent(inout) :: segment !< pointer to segment type
+  integer,                                    intent(in)    :: fld
+  ! Local variables
+  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
+  integer :: n, ishift,jshift
+  real, allocatable, dimension(:,:,:) :: eta ! Segment source data interface heights, [Z -> m]
+  real :: hTolerance = 0.1 !<  Tolerance to exceed adjustment criteria [Z ~> m]
+  real :: hTmp, eTmp, dilate
+  character(len=100) :: mesg
+
+  nz = G%ke
+  hTolerance = 0.1*US%m_to_Z
+
+  ishift=0;jshift=0
+  if (segment%is_E_or_W) then
+     if (segment%field(fld)%name == 'V' .or. segment%field(fld)%name == 'DVDX') then
+        is = segment%HI%isdB ; ie = segment%HI%iedB
+        js = segment%HI%jsdB ; je = segment%HI%jedB
+     else
+        is = segment%HI%isd ; ie = segment%HI%ied
+        js = segment%HI%jsd ; je = segment%HI%jed
+     endif
+     if (segment%direction == OBC_DIRECTION_W) ishift=1
+  else
+     if (segment%field(fld)%name == 'U' .or. segment%field(fld)%name == 'DUDY') then
+        is = segment%HI%isdB ; ie = segment%HI%iedB
+        js = segment%HI%jsdB ; je = segment%HI%jedB
+     else
+        is = segment%HI%isd ; ie = segment%HI%ied
+        js = segment%HI%jsd ; je = segment%HI%jed
+     endif
+     if (segment%direction == OBC_DIRECTION_S) jshift=1
+  endif
+  allocate(eta(is:ie,js:je,nz+1))
+  contractions=0; dilations=0
+  do j=js,je ; do i=is,ie
+    eta(i,j,1)=0.0   ! segment data are assumed to be located on a static grid
+    ! For remapping calls, the entire column will be dilated
+    ! by a factor equal to the ratio of the sum of the geopotential referenced
+    ! source data thicknesses, and the current model thicknesses. This could be
+    ! an issue to be addressed, for instance if we are placing open boundaries
+    ! under ice shelf cavities.
+    do k=2,nz+1
+      eta(i,j,k)=eta(i,j,k-1)-segment%field(fld)%dz_src(i,j,k)
+    enddo
+    ! The normal slope at the boundary is zero by a
+    ! previous call to open_boundary_impose_normal_slope
+    do k=nz+1,1,-1
+      if (-eta(i,j,k) > segment%htot(i,j) + hTolerance) then
+         eta(i,j,k) = -segment%htot(i,j)
+         contractions = contractions + 1
+      endif
+    enddo
+
+    do k=1,nz
+      ! Collapse layers to thinnest possible if the thickness less than
+      ! the thinnest possible (or negative).
+      if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then
+         eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
+         segment%field(fld)%dz_src(i,j,k) = GV%Angstrom_Z
+      else
+         segment%field(fld)%dz_src(i,j,k) = (eta(i,j,K) - eta(i,j,K+1))
+      endif
+    enddo
+
+    !   The whole column is dilated to accommodate deeper topography than
+    ! the bathymetry would indicate.
+    if (-eta(i,j,nz+1) < segment%htot(i,j) - hTolerance) then
+       dilations = dilations + 1
+       if (eta(i,j,1) <= eta(i,j,nz+1)) then
+         do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + segment%htot(i,j)) / real(nz) ; enddo
+       else
+         dilate = (eta(i,j,1) + segment%htot(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
+         do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k) * dilate ; enddo
+       endif
+       do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + segment%field(fld)%dz_src(i,j,k) ; enddo
+    endif
+    ! Now convert thicknesses to units of H.
+    do k=1,nz
+      segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*GV%Z_to_H
+    enddo
+  enddo; enddo
+
+  ! can not do communication call here since only PEs on the current segment are here
+
+  ! call sum_across_PEs(contractions)
+  ! if ((contractions > 0) .and. (is_root_pe())) then
+  !    write(mesg,'("Thickness OBCs were contracted ",'// &
+  !         '"to fit topography in ",I8," places.")') contractions
+  !    call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
+  ! endif
+  ! call sum_across_PEs(dilations)
+  ! if ((dilations > 0) .and. (is_root_pe())) then
+  !    write(mesg,'("Thickness OBCs were dilated ",'// &
+  !         '"to fit topography in ",I8," places.")') dilations
+  !    call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
+  ! endif
+  deallocate(eta)
+
+
+
+end subroutine adjustSegmentEtaToFitBathymetry
+
 !> \namespace mom_open_boundary
 !! This module implements some aspects of internal open boundary
 !! conditions in MOM.