Update the automated max_blocks calculation (CICE-Consortium#954)

Update support for max_blocks=-1. This update computes the blocks required on
each MPI task and then sets that as max_blocks if max_blocks=-1 in namelist.
This is done in ice_distribution and is a function of the decomposition among
other things. Refactor the decomposition computation to defer usage of max_blocks
and eliminate the blockIndex array. Update some indentation formatting in
ice_distribution.F90.

Modify cice.setup and cice_decomp.csh to set max_blocks=-1 unless it's explicitly
defined by the cice.setup -p setting.

Fix a bug in ice_gather_scatter related to zero-ing out of the halo with the
field_loc_noupdate setting. This was zero-ing out the blocks extra times and
there were no problems as long as max_blocks was the same value on all MPI tasks.
With the new implementation of max_blocks=-1, max_blocks can be different values
on different MPI tasks. An error was generated and then the implementation
was fixed so each block on each task is now zeroed out exactly once.

Update diagnostics related to max_block information. Write out the min and max
max_blocks values across MPI tasks.

Add extra allocation/deallocation checks in ice_distribution.F90 and add
a function, ice_memusage_allocErr, to ice_memusage.F90 that checks the
alloc/dealloc return code, writes an error message, and aborts. This
function could be used in other parts of the code as well.

Fix a bug in the io_binary restart output where each task was writing some
output when it should have just been the master task.

Update test cases

Update documentation

cice.setup

@@ -684,7 +684,7 @@ EOF
         set thrd  = `echo ${pesx} | cut -d x -f 2`
         set blckx = `echo ${pesx} | cut -d x -f 3`
         set blcky = `echo ${pesx} | cut -d x -f 4`
-        set mblck = 0
+        set mblck = -1
         if (${task} == 0 || ${thrd} == 0 || ${blckx} == 0 || ${blcky} == 0) then
           echo "${0}: ERROR in -p argument, cannot have zeros"
           exit -1
@@ -696,7 +696,7 @@ EOF
           set thrd = `echo ${pesx} | cut -d x -f 2`
           set blckx = 0
           set blcky = 0
-          set mblck = 0
+          set mblck = -1
           if (${task} == 0 || ${thrd} == 0) then
             echo "${0}: ERROR in -p argument, cannot have zeros"
             exit -1
@@ -708,7 +708,7 @@ EOF
             set thrd = 1
             set blckx = 0
             set blcky = 0
-            set mblck = 0
+            set mblck = -1
             if (${task} == 0) then
               echo "${0}: ERROR in -p argument, cannot have zeros"
               exit -1
@@ -757,7 +757,7 @@ EOF
 
     # update pesx based on use defined settings and machine limits to reflect actual value
     set pesx = ${task}x${thrd}x${blckx}x${blcky}x${mblck}
-    if (${mblck} == 0) then
+    if (${mblck} <= 0) then
       set pesx = ${task}x${thrd}x${blckx}x${blcky}
     endif
     if (${blckx} == 0 || ${blcky} == 0) then

cicecore/cicedyn/infrastructure/comm/mpi/ice_gather_scatter.F90

@@ -1836,12 +1836,12 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) == my_task+1 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block
@@ -1867,8 +1867,8 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, &
          enddo
          enddo
 
-         endif
-      enddo
+       endif
+     enddo
    endif
 
    if (add_mpi_barriers) then
@@ -2222,12 +2222,12 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) == my_task+1 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block
@@ -2253,8 +2253,8 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, &
          enddo
          enddo
 
-         endif
-      enddo
+       endif
+     enddo
    endif
 
    if (add_mpi_barriers) then
@@ -2608,12 +2608,12 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) == my_task+1 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block
@@ -2639,8 +2639,8 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, &
          enddo
          enddo
 
-         endif
-      enddo
+       endif
+     enddo
    endif
 
    if (add_mpi_barriers) then

cicecore/cicedyn/infrastructure/comm/serial/ice_gather_scatter.F90

@@ -1002,12 +1002,12 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) /= 0 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block
@@ -1033,8 +1033,8 @@ subroutine scatter_global_dbl(ARRAY, ARRAY_G, src_task, dst_dist, &
          enddo
          enddo
 
-         endif
-      enddo
+       endif
+     enddo
    endif
 
 !-----------------------------------------------------------------------
@@ -1250,12 +1250,12 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) /= 0 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block
@@ -1281,8 +1281,8 @@ subroutine scatter_global_real(ARRAY, ARRAY_G, src_task, dst_dist, &
          enddo
          enddo
 
-         endif
-      enddo
+       endif
+     enddo
    endif
 
 !-----------------------------------------------------------------------
@@ -1498,12 +1498,12 @@ subroutine scatter_global_int(ARRAY, ARRAY_G, src_task, dst_dist, &
    !-----------------------------------------------------------------
 
    if (field_loc == field_loc_noupdate) then
-      do n=1,nblocks_tot
+     do n=1,nblocks_tot
+       if (dst_dist%blockLocation(n) /= 0 .and. &
+           dst_dist%blockLocalID(n) > 0) then
+
          dst_block = dst_dist%blockLocalID(n)
          this_block = get_block(n,n)
-
-         if (dst_block > 0) then
-
          ! north edge
          do j = this_block%jhi+1,ny_block
          do i = 1, nx_block

cicecore/cicedyn/infrastructure/ice_domain.F90

@@ -223,7 +223,12 @@ subroutine init_domain_blocks
    call broadcast_scalar(nx_global,         master_task)
    call broadcast_scalar(ny_global,         master_task)
 
-   ! Set nprocs if not set in namelist
+!----------------------------------------------------------------------
+!
+! Set nprocs if not explicitly set to valid value in namelist
+!
+!----------------------------------------------------------------------
+
 #ifdef CESMCOUPLED
    nprocs = get_num_procs()
 #else
@@ -235,18 +240,6 @@ subroutine init_domain_blocks
    endif
 #endif
 
-   ! Determine max_blocks if not set
-   if (max_blocks < 1) then
-      call proc_decomposition(nprocs, nprocs_x, nprocs_y)
-      max_blocks=((nx_global-1)/block_size_x/nprocs_x+1) * &
-                  ((ny_global-1)/block_size_y/nprocs_y+1)
-      max_blocks=max(1,max_blocks)
-      if (my_task == master_task) then
-         write(nu_diag,'(/,a52,i6,/)') &
-            '(ice_domain): max_block < 1: max_block estimated to ',max_blocks
-      endif
-   endif
-
 !----------------------------------------------------------------------
 !
 !  perform some basic checks on domain
@@ -321,6 +314,7 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice)
    use ice_boundary, only: ice_HaloCreate
    use ice_distribution, only: create_distribution, create_local_block_ids, ice_distributionGet
    use ice_domain_size, only: max_blocks, nx_global, ny_global
+   use ice_global_reductions, only: global_sum, global_maxval
 
    real (dbl_kind), dimension(nx_global,ny_global), intent(in) :: &
       KMTG           ,&! global topography
@@ -608,9 +602,9 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice)
       work_per_block = 0
    end where
    if (my_task == master_task) then
-      write(nu_diag,*) 'ice_domain work_unit, max_work_unit = ',work_unit, max_work_unit
-      write(nu_diag,*) 'ice_domain nocn = ',minval(nocn),maxval(nocn),sum(nocn)
-      write(nu_diag,*) 'ice_domain work_per_block = ',minval(work_per_block),maxval(work_per_block),sum(work_per_block)
+      write(nu_diag,'(2a,4i9)') subname,' work_unit      = ',work_unit, max_work_unit
+      write(nu_diag,'(2a,4i9)') subname,' nocn           = ',minval(nocn),maxval(nocn),sum(nocn)
+      write(nu_diag,'(2a,4i9)') subname,' work_per_block = ',minval(work_per_block),maxval(work_per_block),sum(work_per_block)
    endif
    deallocate(nocn)
 
@@ -634,8 +628,42 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice)
 
    call create_local_block_ids(blocks_ice, distrb_info)
 
-   ! write out block distribution
-   ! internal check of icedistributionGet as part of verification process
+!----------------------------------------------------------------------
+!
+!  check block sizes and max_blocks
+!
+!----------------------------------------------------------------------
+
+   if (associated(blocks_ice)) then
+      nblocks = size(blocks_ice)
+   else
+      nblocks = 0
+   endif
+
+   tblocks_tmp = global_sum(nblocks, distrb_info)
+   nblocks_max = global_maxval(nblocks, distrb_info)
+
+   if (my_task == master_task) then
+      write(nu_diag,'(2a,i8)') subname,' total number of blocks is', tblocks_tmp
+   endif
+
+   if (nblocks > max_blocks) then
+      write(nu_diag,'(2a,2i6)') subname,' ERROR: nblocks, max_blocks = ',nblocks,max_blocks
+      write(nu_diag,'(2a,2i6)') subname,' ERROR: max_blocks too small: increase to', nblocks_max
+      call abort_ice(subname//' ERROR max_blocks too small', file=__FILE__, line=__LINE__)
+   else if (nblocks_max < max_blocks) then
+      if (my_task == master_task) then
+          write(nu_diag,'(2a,2i6)') subname,' NOTE: max_blocks too large: decrease to', nblocks_max
+      endif
+   endif
+
+!----------------------------------------------------------------------
+!
+! write out block distribution
+! internal check of icedistributionGet as part of verification process
+!
+!----------------------------------------------------------------------
+
    if (debug_blocks) then
 
       call flush_fileunit(nu_diag)
@@ -713,38 +741,6 @@ subroutine init_domain_distribution(KMTG,ULATG,grid_ice)
       endif
    endif
 
-   if (associated(blocks_ice)) then
-      nblocks = size(blocks_ice)
-   else
-      nblocks = 0
-   endif
-   nblocks_max = 0
-   tblocks_tmp = 0
-   do n=0,distrb_info%nprocs - 1
-      nblocks_tmp = nblocks
-      call broadcast_scalar(nblocks_tmp, n)
-      nblocks_max = max(nblocks_max,nblocks_tmp)
-      tblocks_tmp = tblocks_tmp + nblocks_tmp
-   end do
-
-   if (my_task == master_task) then
-      write(nu_diag,*) &
-          'ice: total number of blocks is', tblocks_tmp
-   endif
-
-   if (nblocks_max > max_blocks) then
-      write(outstring,*) ' ERROR: num blocks exceed max: increase max to', nblocks_max
-      call abort_ice(subname//trim(outstring), file=__FILE__, line=__LINE__)
-   else if (nblocks_max < max_blocks) then
-      write(outstring,*) 'WARNING: ice no. blocks too large: decrease max to', nblocks_max
-      if (my_task == master_task) then
-          write(nu_diag,*) ' ********WARNING***********'
-          write(nu_diag,*) subname,trim(outstring)
-          write(nu_diag,*) ' **************************'
-          write(nu_diag,*) ' '
-      endif
-   endif
-
 !----------------------------------------------------------------------
 !
 !  Set up ghost cell updates for each distribution.

cicecore/cicedyn/infrastructure/ice_grid.F90

@@ -301,6 +301,10 @@ subroutine init_grid1
       real (kind=dbl_kind), dimension(:,:), allocatable :: &
          work_g1, work_g2
 
+      integer (kind=int_kind) :: &
+         max_blocks_min, & ! min value of max_blocks across procs
+         max_blocks_max    ! max value of max_blocks across procs
+
       real (kind=dbl_kind) :: &
          rad_to_deg
 
@@ -390,9 +394,15 @@ subroutine init_grid1
       ! write additional domain information
       !-----------------------------------------------------------------
 
+      max_blocks_min = global_minval(max_blocks, distrb_info)
+      max_blocks_max = global_maxval(max_blocks, distrb_info)
       if (my_task == master_task) then
-        write(nu_diag,'(a26,i6)') '  Block size:  nx_block = ',nx_block
-        write(nu_diag,'(a26,i6)') '               ny_block = ',ny_block
+        write(nu_diag,*        ) ''
+        write(nu_diag,'(2a)'   ) subname,' Block size:'
+        write(nu_diag,'(2a,i8)') subname,'   nx_block        = ',nx_block
+        write(nu_diag,'(2a,i8)') subname,'   ny_block        = ',ny_block
+        write(nu_diag,'(2a,i8)') subname,'   min(max_blocks) = ',max_blocks_min
+        write(nu_diag,'(2a,i8)') subname,'   max(max_blocks) = ',max_blocks_max
       endif
 
       end subroutine init_grid1