Update outer boundary implementation to support open boundary conditions (CICE-Consortium#1062)

Refactor CICE to support open boundary conditions where values may be imposed on the outer boundary. Prior to this update, the HaloUpdate was zeroing out the halo before updating the halo. This meant any values set on the halo were lost anytime a HaloUpdate was executed. Several other issues upstream and downstream were found as this was implemented. In particular, the global index set on the outer boundary for open and closed boundary conditions was "special", a zero global index was used internally as a flag, and both the scatter_global and haloUpdate was covering up some issues with uninitialized data.

Update the i_glob and j_glob indexing for open and closed boundary conditions on the outer boundary. For open, closed, and tripole (south) boundaries, now extrapolate index values on the outer boundary. This might mean indices are less than 0 or greater than n[x,y]_global+2*nghost. In the original implementation, these outer boundary indices were either set to some internal values so they would be set with the scatter_global method or to zero which was treated as a special value. The zero special value feature has been removed and the scatter_global implementation has been refactored to fill all gridpoints with valid global indices and leave others unset. The only special case remaining in the global indexing is for tripole grids on the north boundary where the j global index is set to a negative value.

- Refactor HaloUpdate in ice_boundary in serial and mpi to avoid zeroing out the halo. The ew and ns directions are handled separately. The new implementation will always zero the internal halo points. It will also zero the outer halo if cyclic or tripole is specified or if a fillvalue is passed into the haloUpdate. Otherwise, the outer halo is left unchanged.
   - Added ewBoundaryType and nsBoundaryType to halo datatype to keep track of the boundary conditions for each halo datatype.
   - Zero out the outer boundary during HaloUpdate only under certain conditions including cyclic or tripole bcs or if a fillValue is explicitly passed
- Add a few fillvalue=c0 arguments to HaloUpdate calls during initialization to force the outer boundary values to zero, mostly during initialization of grid data and similar cases.
- Add several HaloExtrapolate calls to initialize boundary values for grid fields. This replaces those fields being set by the incorrect global indices and the scatter.
- Explicitly zero out all allocated variables at initialization
- Explicitly zero out restart fields before reading to be extra careful
- Explicitly zero out some arrays in the evp1d before use
- Replace the G_HTN and G_HTE initialization with a call to gather_global_ext on HTN and HTE respectively after HTN and HTE are set and the halos updated.
- Move the closed boundary condition abort to initialization in ice_domain.F90
- Refactor halochk unit test, add explicit fill tests
- Update documentation

Additional Features

- Update ice_HaloExtrapolate to support more than 1 ghost cell
- Add support for kmt_type = none with in subroutine rectgrid for rectangular grids. This allows a user to setup a rectangular grid with no land.
- Update the restart_ext implementation in ice_read_write. There was a bug if restart_ext=.false. was passed, nothing would happen. The update also allows some extra restart_ext logic in io_netcdf/ice_restart.F90 to be cleaned up.
- Carry out some code cleanup in init_domain_distribution and init_grid2
- Rename iglb/jglb to isrc/jsrc consistent in scatter_global_ext to be consistent with other scatter methods.

All tests are bit-for-bit except the halochk unittest which was refactored and passes. There is an issue in the evp1d implementation where on some compilers when the debug flag is on, the code will trap on a floating point error. If the floating point error isn't trapped, the cases run to completion and are bit-for-bit with the current trunk. This suggests there is something like a 0/0 calculation being trapped in the evp1d which has no impact on the solution, is only caught by a subset of compilers, and probably requires an if test to avoid. I spent some time trying to track down the error without any success.

Author: apcraig

cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90

@@ -94,7 +94,7 @@ subroutine dyn_evp1d_init
     endif
 
     ! gather from blks to global
-    call gather_static(G_uarear,  G_dxT, G_dyT, G_tmask)
+    call gather_static(G_uarear, G_dxT, G_dyT, G_tmask)
 
     ! calculate number of water points (T and U). Only needed for the static version
     ! tmask in ocean/ice
@@ -349,7 +349,6 @@ subroutine evp1d_alloc_static_na(na0)
        call abort_ice(subname//' ERROR: allocating', file=__FILE__, line=__LINE__)
     endif
 
-
     allocate(indxTi(1:na0), &
              indxTj(1:na0), &
              stat=ierr)
@@ -628,6 +627,11 @@ subroutine gather_static(G_uarear, G_dxT, G_dyT, G_Tmask)
 
      character(len=*), parameter :: subname = '(gather_static)'
 
+     G_uarear = c0
+     G_dyT = c0
+     G_dxT = c0
+     G_tmask = .false.
+
      ! copy from distributed I_* to G_*
      call gather_global_ext(G_uarear, uarear, master_task, distrb_info)
      call gather_global_ext(G_dxT   , dxT   , master_task, distrb_info)
@@ -977,6 +981,37 @@ subroutine convert_1d_2d_dyn(na0         , navel0      ,
      integer(kind=int_kind) ::  lo, up, iw, i, j
      character(len=*), parameter :: subname = '(convert_1d_2d_dyn)'
 
+     G_stressp_1  = c0
+     G_stressp_2  = c0
+     G_stressp_3  = c0
+     G_stressp_4  = c0
+     G_stressm_1  = c0
+     G_stressm_2  = c0
+     G_stressm_3  = c0
+     G_stressm_4  = c0
+     G_stress12_1 = c0
+     G_stress12_2 = c0
+     G_stress12_3 = c0
+     G_stress12_4 = c0
+     G_strength   = c0
+     G_cdn_ocn    = c0
+     G_aiu        = c0
+     G_uocn       = c0
+     G_vocn       = c0
+     G_waterxU    = c0
+     G_wateryU    = c0
+     G_forcexU    = c0
+     G_forceyU    = c0
+     G_umassdti   = c0
+     G_fmU        = c0
+     G_strintxU   = c0
+     G_strintyU   = c0
+     G_Tbu        = c0
+     G_uvel       = c0
+     G_vvel       = c0
+     G_taubxU     = c0
+     G_taubyU     = c0
+
      lo=1
      up=na0
      do iw = lo, up

cicecore/cicedyn/dynamics/ice_dyn_shared.F90

@@ -193,13 +193,24 @@ subroutine alloc_dyn_shared
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of memory')
 
+      uvel_init  = c0
+      vvel_init  = c0
+      iceTmask   = .false.
+      iceUmask   = .false.
+      fcor_blk   = c0
+      DminTarea  = c0
+
       allocate( &
          fld2(nx_block,ny_block,2,max_blocks), &
          fld3(nx_block,ny_block,3,max_blocks), &
          fld4(nx_block,ny_block,4,max_blocks), &
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of memory')
 
+      fld2 = c0
+      fld3 = c0
+      fld4 = c0
+
       allocate( &
          cyp(nx_block,ny_block,max_blocks), & ! 1.5*HTE - 0.5*HTW
          cxp(nx_block,ny_block,max_blocks), & ! 1.5*HTN - 0.5*HTS
@@ -208,12 +219,19 @@ subroutine alloc_dyn_shared
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of memory')
 
+      cyp = c0
+      cxp = c0
+      cym = c0
+      cxm = c0
+
       if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
          allocate( &
             dxhy(nx_block,ny_block,max_blocks), & ! 0.5*(HTE - HTW)
             dyhx(nx_block,ny_block,max_blocks), & ! 0.5*(HTN - HTS)
             stat=ierr)
          if (ierr/=0) call abort_ice(subname//': Out of memory')
+         dxhy = c0
+         dyhx = c0
       endif
 
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
@@ -228,6 +246,14 @@ subroutine alloc_dyn_shared
             fcorN_blk  (nx_block,ny_block,max_blocks), &   ! Coriolis
             stat=ierr)
          if (ierr/=0) call abort_ice(subname//': Out of memory')
+         uvelE_init  = c0
+         vvelE_init  = c0
+         uvelN_init  = c0
+         vvelN_init  = c0
+         iceEmask    = .false.
+         iceNmask    = .false.
+         fcorE_blk   = c0
+         fcorN_blk   = c0
       endif
 
       end subroutine alloc_dyn_shared