Fixes for uninitialized variables (snan testing) (#579)

* CICE_RunMod: query 'Lsub' from Icepack

When CICE is compiled with the 'CICE_IN_NEMO' CPP macro, and 'calc_Tsfc'
is false, ice_step::coupling_prep calls CICE_RunMod::srcflux_to_ocn to
transfer heat fluxes on ice-free grid points to the ocean.

The calculation in srcflux_to_ocn uses the Icepack parameter 'Lsub', but
without declaring this variable nor querying it from Icepack, which
results in a compile-time failure when using the standalone driver with
'CICE_IN_NEMO' defined.

Fix that by correctly declaring 'Lsub' and querying it using the Icepack
interface.

* ice_flux: remove 'CICE_IN_NEMO' CPP around sst initialization

In subroutine 'init_coupler_flux', the initialization of the 'sst' array
to 'Tf' is protected by a '#ifndef CICE_IN_NEMO' preprocessor directive.
This has been the case since CICE-Consortium/CICE-svn-trunk@151b9af
(cice-4.0~17, 2008-04-28), though that commit does not explain why.
This is probably because in some conditions (depending on active CPPs at
NEMO compilation time, as well as NEMO namelist parameters), the SST
needs to be passed early to CICE, before calling CICE_Initialize (see
the CICE coupling interface in NEMO [1]), and initializing it to 'Tf' in
'init_coupler_flux' would erase the values already passed from the
ocean.

If however the SST is *not* passed to CICE before CICE_Initialize is
called, and 'CICE_IN_NEMO' is in use, then 'ice_init::set_state_var'
reads from the uninitialized 'sst' array when placing ice where the
ocean surface is cold.

In ECCC's in-house version of the coupling interface (sbc_ice_cice),
extensive changes have been made to work around bugs in the original
version and limitations due to the sequence in which CICE and NEMO fields
are initialized, such that we manually call 'ice_init::init_state' a
second time at the end of subroutine sbc_ice_cice::cice_sbc_init.

To avoid using a uninitialized array, remove the '#ifdef
CICE_IN_NEMO' around the initialization of the 'sst' array to 'Tf', so
that it is always initialized. These values will anyway be replaced by
the "correct" ones when init_state is called a second time in our
in-house version.

[1] http://forge.ipsl.jussieu.fr/nemo/browser/NEMO/releases/release-3.6/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90#L176

* Fix uninitialized values and an array bounds issue in spacecurve

Picked up by testing with initialized snan and debug options in decomps
Update test suite to add debug tests for multiple decomps

* Fix several decomp bugs and add new decomp tests

- Fix bug in ice_boundary fill in halo, was not taking into account padding correctly
- Fix bug in ice_blocks for cases where the active size is 1 gridcell wide
- Update ice_grid initialization to clean up prior workarounds.
- Update some intel compilers to add -init=snan,arrays for debug builds
- Add new tests to the decomp test suite including tests with debugging on,
  tests with the active size 1 gridcell wide, and tests with a single block
  that's bigger than the whole grid.

* Update some uninitialized variables.
Back off on some of the full debug compiler flags.

* revert cheyenne_gnu compiler flags

* add 1x1, 2x2, 3x3 block tests

* add abort to default bathymetry if kmt gt nlevel

* update bathymetry k/kmt calculation

Co-authored-by: Philippe Blain <[email protected]>

cicecore/cicedynB/general/ice_flux.F90

@@ -654,9 +654,7 @@ subroutine init_coupler_flux
       enddo
       enddo
 
-#ifndef CICE_IN_NEMO
       sst   (:,:,:) = Tf(:,:,:)       ! sea surface temp (C)
-#endif
       qdp   (:,:,:) = c0              ! deep ocean heat flux (W/m^2)
       hmix  (:,:,:) = c20             ! ocean mixed layer depth (m)
       hwater(:,:,:) = bathymetry(:,:,:) ! ocean water depth (m)

cicecore/cicedynB/infrastructure/ice_blocks.F90

@@ -252,8 +252,8 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, &
             !*** set last physical point if padded domain
 
             else if (j_global(j,n) == ny_global .and. &
-                     j > all_blocks(n)%jlo      .and. &
-                     j < all_blocks(n)%jhi) then
+                     j >= all_blocks(n)%jlo      .and. &
+                     j <  all_blocks(n)%jhi) then
                all_blocks(n)%jhi = j   ! last physical point in padded domain
             endif
          end do
@@ -300,8 +300,8 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, &
             !*** last physical point in padded domain
 
             else if (i_global(i,n) == nx_global .and. &
-                     i > all_blocks(n)%ilo      .and. &
-                     i < all_blocks(n)%ihi) then
+                     i >= all_blocks(n)%ilo      .and. &
+                     i <  all_blocks(n)%ihi) then
                all_blocks(n)%ihi = i
             endif
          end do

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -26,7 +26,7 @@ module ice_grid
       use ice_domain, only: blocks_ice, nblocks, halo_info, distrb_info, &
           ew_boundary_type, ns_boundary_type, init_domain_distribution
       use ice_fileunits, only: nu_diag, nu_grid, nu_kmt, &
-          get_fileunit, release_fileunit
+          get_fileunit, release_fileunit, flush_fileunit
       use ice_gather_scatter, only: gather_global, scatter_global
       use ice_read_write, only: ice_read, ice_read_nc, ice_read_global, &
           ice_read_global_nc, ice_open, ice_open_nc, ice_close_nc
@@ -384,11 +384,9 @@ subroutine init_grid2
       ! T-grid cell and U-grid cell quantities
       !-----------------------------------------------------------------
 
-!     tarea(:,:,:) = c0
-
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -486,7 +484,7 @@ subroutine init_grid2
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
       !$OMP                     angle_0,angle_w,angle_s,angle_sw)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -642,7 +640,7 @@ subroutine popgrid
       kmt(:,:,:) = c0
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -785,7 +783,7 @@ subroutine popgrid_nc
       kmt(:,:,:) = c0
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -1104,7 +1102,7 @@ subroutine latlongrid
 
      !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -1198,15 +1196,9 @@ subroutine rectgrid
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
 
-      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
-      do iblk = 1, nblocks
-         do j = 1, ny_block
-         do i = 1, nx_block
-            ANGLE(i,j,iblk) = c0              ! "square with the world"
-         enddo
-         enddo
-      enddo
-      !$OMP END PARALLEL DO
+      hm (:,:,:) = c0
+      kmt(:,:,:) = c0
+      angle(:,:,:) = c0   ! "square with the world"
 
       allocate(work_g1(nx_global,ny_global))
 
@@ -1396,7 +1388,7 @@ subroutine cpomgrid
       kmt(:,:,:) = c0
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -1636,11 +1628,10 @@ subroutine makemask
       !-----------------------------------------------------------------
 
       bm = c0
-!     uvm = c0
 
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -1663,12 +1654,19 @@ subroutine makemask
                            field_loc_center, field_type_scalar)
       call ice_timer_stop(timer_bound)
 
-      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+      !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         do j = 1, ny_block
-         do i = 1, nx_block
-            tmask(i,j,iblk) = .false.
-            umask(i,j,iblk) = .false.
+         this_block = get_block(blocks_ice(iblk),iblk)
+         ilo = this_block%ilo
+         ihi = this_block%ihi
+         jlo = this_block%jlo
+         jhi = this_block%jhi
+
+         ! needs to cover halo (no halo update for logicals)
+         tmask(:,:,iblk) = .false.
+         umask(:,:,iblk) = .false.
+         do j = jlo-nghost, jhi+nghost
+         do i = ilo-nghost, ihi+nghost
             if ( hm(i,j,iblk) > p5) tmask(i,j,iblk) = .true.
             if (uvm(i,j,iblk) > p5) umask(i,j,iblk) = .true.
          enddo
@@ -1684,11 +1682,14 @@ subroutine makemask
          tarean(:,:,iblk) = c0
          tareas(:,:,iblk) = c0
 
-         do j = 1, ny_block
-         do i = 1, nx_block
+         do j = jlo,jhi
+         do i = ilo,ihi
 
-            if (ULAT(i,j,iblk) >= -puny) lmask_n(i,j,iblk) = .true. ! N. Hem.
-            if (ULAT(i,j,iblk) <  -puny) lmask_s(i,j,iblk) = .true. ! S. Hem.
+            if (ULAT(i,j,iblk) >= -puny) then
+               lmask_n(i,j,iblk) = .true. ! N. Hem.
+            else
+               lmask_s(i,j,iblk) = .true. ! S. Hem.
+            endif
 
             ! N hemisphere area mask (m^2)
             if (lmask_n(i,j,iblk)) tarean(i,j,iblk) = tarea(i,j,iblk) &
@@ -1743,7 +1744,7 @@ subroutine Tlatlon
       !$OMP                     x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4, &
       !$OMP                     tx,ty,tz,da)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -1915,7 +1916,7 @@ subroutine to_ugrid(work1,work2)
 
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -2000,7 +2001,7 @@ subroutine to_tgrid(work1, work2)
 
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -2073,7 +2074,7 @@ subroutine gridbox_corners
 
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
       do iblk = 1, nblocks
-         this_block = get_block(blocks_ice(iblk),iblk)         
+         this_block = get_block(blocks_ice(iblk),iblk)
          ilo = this_block%ilo
          ihi = this_block%ihi
          jlo = this_block%jlo
@@ -2400,8 +2401,9 @@ subroutine get_bathymetry
          do iblk = 1, nblocks
             do j = 1, ny_block
             do i = 1, nx_block
-               k = min(nint(kmt(i,j,iblk)),nlevel)
-               if (k > puny) bathymetry(i,j,iblk) = depth(k)
+               k = nint(kmt(i,j,iblk))
+               if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error')
+               if (k > 0) bathymetry(i,j,iblk) = depth(k)
             enddo
             enddo
          enddo
@@ -2491,8 +2493,8 @@ subroutine get_bathymetry_popfile
       do iblk = 1, nblocks
          do j = 1, ny_block
          do i = 1, nx_block
-            k = kmt(i,j,iblk)
-            if (k > nlevel) call abort_ice(subname//' kmt/nlevel error')
+            k = nint(kmt(i,j,iblk))
+            if (k > nlevel) call abort_ice(subname//' kmt gt nlevel error')
             if (k > 0) bathymetry(i,j,iblk) = depth(k)
          enddo
          enddo

cicecore/cicedynB/infrastructure/ice_restoring.F90

@@ -98,6 +98,11 @@ subroutine ice_HaloRestore_init
              vsnon_rest(nx_block,ny_block,ncat,max_blocks), &
              trcrn_rest(nx_block,ny_block,ntrcr,ncat,max_blocks))
 
+   aicen_rest(:,:,:,:) = c0
+   vicen_rest(:,:,:,:) = c0
+   vsnon_rest(:,:,:,:) = c0
+   trcrn_rest(:,:,:,:,:) = c0
+
 !-----------------------------------------------------------------------
 ! initialize
 ! halo cells have to be filled manually at this stage

cicecore/drivers/standalone/cice/CICE_RunMod.F90

@@ -628,11 +628,12 @@ subroutine sfcflux_to_ocn(nx_block,   ny_block,     &
 
       real (kind=dbl_kind)    :: &
           puny, &          !
+          Lsub, &          !
           rLsub            ! 1/Lsub
 
       character(len=*), parameter :: subname = '(sfcflux_to_ocn)'
 
-      call icepack_query_parameters(puny_out=puny)
+      call icepack_query_parameters(puny_out=puny, Lsub_out=Lsub)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)

cicecore/shared/ice_init_column.F90

@@ -408,7 +408,7 @@ subroutine init_shortwave
          do i = ilo, ihi
 
                if (aicen(i,j,n,iblk) > puny) then
-                  
+
                   alvdf(i,j,iblk) = alvdf(i,j,iblk) &
                        + alvdfn(i,j,n,iblk)*aicen(i,j,n,iblk)
                   alidf(i,j,iblk) = alidf(i,j,iblk) &
@@ -417,7 +417,7 @@ subroutine init_shortwave
                        + alvdrn(i,j,n,iblk)*aicen(i,j,n,iblk)
                   alidr(i,j,iblk) = alidr(i,j,iblk) &
                        + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk)
-                  
+
                   netsw = swvdr(i,j,iblk) + swidr(i,j,iblk) &
                         + swvdf(i,j,iblk) + swidf(i,j,iblk)
                   if (netsw > puny) then ! sun above horizon
@@ -428,12 +428,12 @@ subroutine init_shortwave
                      albpnd(i,j,iblk) = albpnd(i,j,iblk) &
                           + albpndn(i,j,n,iblk)*aicen(i,j,n,iblk)
                   endif
-                  
+
                   apeff_ai(i,j,iblk) = apeff_ai(i,j,iblk) &
                        + apeffn(i,j,n,iblk)*aicen(i,j,n,iblk)
                   snowfrac(i,j,iblk) = snowfrac(i,j,iblk) &
                        + snowfracn(i,j,n,iblk)*aicen(i,j,n,iblk)
-               
+
                endif ! aicen > puny
 
          enddo ! i