From 104a1c97b6843ca79f977f3a3e4bce9bd68f9989 Mon Sep 17 00:00:00 2001 From: apcraig Date: Thu, 15 Jul 2021 19:00:59 -0600 Subject: [PATCH] Fix 1x1 and 2x2 blocksize error - Fix 1x1 and 2x2 non-bit-for-bit results. It turns out the 1x1 and 2x2 blocks were bit-for-bit physicallly. The problem was that a set of blocks were "land block eliminated" near the equator incorrectly. These blocks never have sea ice form, so it's not a problem scientifically, but it changes the masking and restart files. The fix is a one line change in ice_domain.F90 which that sets the minimum "flat" value to 1. This error could only happen in cases where the blocks were within 0.5 degrees of the equator, and a set of conditions were required where an active gridcell was set to inactive which then resulted in the entire block being eliminated. - Update the initialization of ARRAY_G in the gather_global_ext interfaces. Get rid of older, more complex approach which left some blocks uninitialized. - Set psalt(n) to zero by default. Was uninitialized in cases where vice=0. - Add a test to check bit-for-bit in log files for 1x1 blocks with bfbflag set to reprosum. --- .../cicedynB/analysis/ice_diagnostics.F90 | 1 + .../comm/mpi/ice_gather_scatter.F90 | 136 ++---------------- .../comm/serial/ice_gather_scatter.F90 | 38 +---- .../cicedynB/infrastructure/ice_domain.F90 | 2 +- configuration/scripts/tests/reprosum_suite.ts | 1 + 5 files changed, 23 insertions(+), 155 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_diagnostics.F90 b/cicecore/cicedynB/analysis/ice_diagnostics.F90 index 6b9b32301..760952fc5 100644 --- a/cicecore/cicedynB/analysis/ice_diagnostics.F90 +++ b/cicecore/cicedynB/analysis/ice_diagnostics.F90 @@ -825,6 +825,7 @@ subroutine runtime_diags (dt) enddo endif endif + psalt(n) = c0 if (vice(i,j,iblk) /= c0) psalt(n) = work2(i,j,iblk)/vice(i,j,iblk) pTsfc(n) = trcr(i,j,nt_Tsfc,iblk) ! ice/snow sfc temperature pevap(n) = evap(i,j,iblk)*dt/rhoi ! sublimation/condensation diff --git a/cicecore/cicedynB/infrastructure/comm/mpi/ice_gather_scatter.F90 b/cicecore/cicedynB/infrastructure/comm/mpi/ice_gather_scatter.F90 index 010a5c8c4..5f0d0af89 100644 --- a/cicecore/cicedynB/infrastructure/comm/mpi/ice_gather_scatter.F90 +++ b/cicecore/cicedynB/infrastructure/comm/mpi/ice_gather_scatter.F90 @@ -636,6 +636,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) else special_value = spval_dbl endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -744,92 +745,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - !*** fill land blocks with special values - - else if (src_dist%blockLocation(n) == 0) then - - this_block = get_block(n,n) - - ! interior - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - -#ifdef CICE_IN_NEMO -!echmod: this code is temporarily wrapped for nemo pending further testing elsewhere - ! fill ghost cells - if (this_block%jblock == 1) then - ! south block - do j=1, nghost - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost,j) = special_value - end do - end do - if (this_block%iblock == 1) then - ! southwest corner - do j=1, nghost - do i=1, nghost - ARRAY_G(i,j) = special_value - end do - end do - endif - endif - if (this_block%jblock == nblocks_y) then - ! north block - do j=1, nghost - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - ny_global + nghost + j) = special_value - end do - end do - if (this_block%iblock == nblocks_x) then - ! northeast corner - do j=1, nghost - do i=1, nghost - ARRAY_G(nx-i+1, ny-j+1) = special_value - end do - end do - endif - endif - if (this_block%iblock == 1) then - ! west block - do j=this_block%jlo,this_block%jhi - do i=1, nghost - ARRAY_G(i,this_block%j_glob(j)+nghost) = special_value - end do - end do - if (this_block%jblock == nblocks_y) then - ! northwest corner - do j=1, nghost - do i=1, nghost - ARRAY_G(i, ny-j+1) = special_value - end do - end do - endif - endif - if (this_block%iblock == nblocks_x) then - ! east block - do j=this_block%jlo,this_block%jhi - do i=1, nghost - ARRAY_G(nx_global + nghost + i, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - if (this_block%jblock == 1) then - ! southeast corner - do j=1, nghost - do i=1, nghost - ARRAY_G( nx-i+1,j) = special_value - end do - end do - endif - endif -#endif - - endif + endif ! src_dist%blockLocation end do @@ -939,7 +855,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) ! !----------------------------------------------------------------------- - else + else ! master task allocate(snd_request(nblocks_tot), & snd_status (MPI_STATUS_SIZE, nblocks_tot)) @@ -960,7 +876,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) call MPI_WAITALL(nsends, snd_request, snd_status, ierr) deallocate(snd_request, snd_status) - endif + endif ! master task if (add_mpi_barriers) then call ice_barrier() @@ -1028,8 +944,9 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) if (present(spc_val)) then special_value = spc_val else - special_value = 0 !MHRI NOTE: 0,1,-999,??? + special_value = -9999 endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -1138,21 +1055,7 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - !*** fill land blocks with special values - - else if (src_dist%blockLocation(n) == 0) then - - this_block = get_block(n,n) - - ! interior - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - - endif + endif ! src_dist%blockLocation end do @@ -1262,7 +1165,7 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) ! !----------------------------------------------------------------------- - else + else ! master task allocate(snd_request(nblocks_tot), & snd_status (MPI_STATUS_SIZE, nblocks_tot)) @@ -1283,7 +1186,7 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) call MPI_WAITALL(nsends, snd_request, snd_status, ierr) deallocate(snd_request, snd_status) - endif + endif ! master task if (add_mpi_barriers) then call ice_barrier() @@ -1353,6 +1256,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) else special_value = .false. !MHRI NOTE: .true./.false. ??? endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -1461,21 +1365,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - !*** fill land blocks with special values - - else if (src_dist%blockLocation(n) == 0) then - - this_block = get_block(n,n) - - ! interior - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - - endif + endif ! src_dist%blockLocation end do @@ -1585,7 +1475,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) ! !----------------------------------------------------------------------- - else + else ! master task allocate(snd_request(nblocks_tot), & snd_status (MPI_STATUS_SIZE, nblocks_tot)) @@ -1606,7 +1496,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) call MPI_WAITALL(nsends, snd_request, snd_status, ierr) deallocate(snd_request, snd_status) - endif + endif ! master task if (add_mpi_barriers) then call ice_barrier() diff --git a/cicecore/cicedynB/infrastructure/comm/serial/ice_gather_scatter.F90 b/cicecore/cicedynB/infrastructure/comm/serial/ice_gather_scatter.F90 index 418c80f61..515cdfa65 100644 --- a/cicecore/cicedynB/infrastructure/comm/serial/ice_gather_scatter.F90 +++ b/cicecore/cicedynB/infrastructure/comm/serial/ice_gather_scatter.F90 @@ -373,6 +373,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) else special_value = spval_dbl endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -477,16 +478,7 @@ subroutine gather_global_ext_dbl(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - else !*** fill land blocks with special values - - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - - endif + endif ! src_dist%blockLocation end do @@ -537,8 +529,9 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) if (present(spc_val)) then special_value = spc_val else - special_value = 0 !MHRI: 0,1,999,-999 ?? + special_value = -9999 endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -643,16 +636,7 @@ subroutine gather_global_ext_int(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - else !*** fill land blocks with special values - - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - - endif + endif ! src_dist%blockLocation end do @@ -705,6 +689,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) else special_value = .false. !MHRI: true/false endif + ARRAY_G = special_value nx = nx_global + 2*nghost ny = ny_global + 2*nghost @@ -809,16 +794,7 @@ subroutine gather_global_ext_log(ARRAY_G, ARRAY, dst_task, src_dist, spc_val) endif endif - else !*** fill land blocks with special values - - do j=this_block%jlo,this_block%jhi - do i=this_block%ilo,this_block%ihi - ARRAY_G(this_block%i_glob(i)+nghost, & - this_block%j_glob(j)+nghost) = special_value - end do - end do - - endif + endif ! src_dist%blockLocation end do diff --git a/cicecore/cicedynB/infrastructure/ice_domain.F90 b/cicecore/cicedynB/infrastructure/ice_domain.F90 index 52f0da850..1dfdd0428 100644 --- a/cicecore/cicedynB/infrastructure/ice_domain.F90 +++ b/cicecore/cicedynB/infrastructure/ice_domain.F90 @@ -441,7 +441,7 @@ subroutine init_domain_distribution(KMTG,ULATG) !---------------------------------------------------------------------- if (distribution_wght == 'latitude') then - flat = NINT(abs(ULATG*rad_to_deg), int_kind) ! linear function + flat = max(NINT(abs(ULATG*rad_to_deg), int_kind),1) ! linear function else flat = 1 endif diff --git a/configuration/scripts/tests/reprosum_suite.ts b/configuration/scripts/tests/reprosum_suite.ts index d65370e0a..dd6a6d56b 100644 --- a/configuration/scripts/tests/reprosum_suite.ts +++ b/configuration/scripts/tests/reprosum_suite.ts @@ -8,4 +8,5 @@ logbfb gx3 1x20x5x29x80 dsectrobin,diag1,short,reprosum l logbfb gx3 8x2x8x10x20 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum logbfb gx3 6x2x50x58x1 droundrobin,diag1,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum logbfb gx3 6x2x4x29x18 dspacecurve,diag1,maskhalo,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum +logbfb gx3 17x2x1x1x800 droundrobin,diag1,maskhalo,reprosum logbfb_gx3_4x2x25x29x4_diag1_dslenderX2_reprosum #logbfb gx3 8x2x8x10x20 droundrobin,diag1 logbfb_gx3_4x2x25x29x4_diag1_dslenderX2