From 180e0ea12db99a6512aae70efd2a07b36178b4a1 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Wed, 26 Sep 2018 14:56:37 -0400
Subject: [PATCH] (*)Correct the grid used for allocate_forcing_type

  Use the correct grid type in allocate_forcing_type and allocate_mech_forcing
calls inside initialize_ice_shelf, and added error checks to verify that
compatible grids are being used by the ice shelf and ocean.  All answers in the
existing MOM6 test cases are bitwise identical, but this corrects a problem in
which symmetric memory ocean models fail for the ISOMIP test case unless the
ice shelf grid also uses symmetric memory.
---
 src/ice_shelf/MOM_ice_shelf.F90 | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index e6989caa54..24389af17f 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/src/ice_shelf/MOM_ice_shelf.F90
@@ -754,6 +754,10 @@ subroutine add_shelf_forces(G, CS, forces, do_shelf_area)
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed
 
+  if ((CS%grid%isc /= G%isc) .or. (CS%grid%iec /= G%iec) .or. &
+      (CS%grid%jsc /= G%jsc) .or. (CS%grid%jec /= G%jec)) &
+    call MOM_error(FATAL,"add_shelf_forces: Incompatible ocean and ice shelf grids.")
+
   ISS => CS%ISS
 
   find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area
@@ -823,6 +827,10 @@ subroutine add_shelf_pressure(G, CS, fluxes)
   integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
+  if ((CS%grid%isc /= G%isc) .or. (CS%grid%iec /= G%iec) .or. &
+      (CS%grid%jsc /= G%jsc) .or. (CS%grid%jec /= G%jec)) &
+    call MOM_error(FATAL,"add_shelf_pressure: Incompatible ocean and ice shelf grids.")
+
   do j=js,je ; do i=is,ie
     press_ice = (CS%ISS%area_shelf_h(i,j) * G%IareaT(i,j)) * (CS%g_Earth * CS%ISS%mass_shelf(i,j))
     if (associated(fluxes%p_surf)) then
@@ -877,6 +885,10 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed
 
+  if ((CS%grid%isc /= G%isc) .or. (CS%grid%iec /= G%iec) .or. &
+      (CS%grid%jsc /= G%jsc) .or. (CS%grid%jec /= G%jec)) &
+    call MOM_error(FATAL,"add_shelf_flux: Incompatible ocean and ice shelf grids.")
+
   ISS => CS%ISS
 
   call add_shelf_pressure(G, CS, fluxes)
@@ -1098,16 +1110,17 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl
   call set_grid_metrics(dG, param_file)
   ! call set_diag_mediator_grid(CS%grid, CS%diag)
 
-  ! The ocean grid is possibly different
-  if (associated(ocn_grid)) CS%ocn_grid => ocn_grid
+  ! The ocean grid possibly uses different symmetry.
+  if (associated(ocn_grid)) then ; CS%ocn_grid => ocn_grid
+  else ; CS%ocn_grid => CS%grid ; endif
 
   ! Convenience pointers
   G => CS%grid
   OG => CS%ocn_grid
 
   if (is_root_pe()) then
-   write(0,*) 'OG: ', OG%isd, OG%isc, OG%iec, OG%ied, OG%jsd, OG%jsc, OG%jsd, OG%jed
-   write(0,*) 'IG: ', G%isd, G%isc, G%iec, G%ied, G%jsd, G%jsc, G%jsd, G%jed
+    write(0,*) 'OG: ', OG%isd, OG%isc, OG%iec, OG%ied, OG%jsd, OG%jsc, OG%jsd, OG%jed
+    write(0,*) 'IG: ', G%isd, G%isc, G%iec, G%ied, G%jsd, G%jsc, G%jsd, G%jed
   endif
 
   CS%Time = Time ! ### This might not be in the right place?
@@ -1344,10 +1357,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl
      ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
      ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
     if (present(fluxes)) &
-      call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., &
+      call allocate_forcing_type(CS%ocn_grid, fluxes, ustar=.true., shelf=.true., &
                                  press=.true., water=CS%isthermo, heat=CS%isthermo)
     if (present(forces)) &
-      call allocate_mech_forcing(G, forces, ustar=.true., shelf=.true., press=.true.)
+      call allocate_mech_forcing(CS%ocn_grid, forces, ustar=.true., shelf=.true., press=.true.)
   else
     call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
     if (present(fluxes)) &