ice_dyn_vp: do not skip halo updates in 'pgmres' under 'bfbflag'

The 'pgmres' subroutine implements a separate GMRES solver and is used
as a preconditioner for the FGMRES linear solver. Since it is only a
preconditioner, it was decided to skip the halo updates after computing
the matrix-vector product (in 'matvec'), for efficiency.

This leads to non-reproducibility since the content of the non-updated
halos depend on the block / MPI distribution.

Add the required halo updates, but only perform them when we are
explicitely asking for bit-for-bit global sums, i.e. when 'bfbflag' is
set to something else than 'not'.

Adjust the interfaces of 'pgmres' and 'precondition' (from which
'pgmres' is called) to accept 'halo_info_mask', since it is needed for
masked updates.

Closes CICE-Consortium#518

cicecore/cicedynB/dynamics/ice_dyn_vp.F90

@@ -2914,6 +2914,7 @@ subroutine fgmres (zetax2   , etax2   , &
             call precondition(zetax2      , etax2          , &
                               Cb          , vrel           , &
                               umassdti    ,                  &
+                              halo_info_mask,                &
                               arnoldi_basis_x(:,:,:,initer), &
                               arnoldi_basis_y(:,:,:,initer), &
                               diagx       , diagy          , &
@@ -3109,13 +3110,19 @@ end subroutine fgmres
       subroutine pgmres (zetax2   , etax2   , &
                          Cb       , vrel    , &
                          umassdti ,           &
+                         halo_info_mask,      &
                          bx       , by      , &
                          diagx    , diagy   , &
                          tolerance, maxinner, &
                          maxouter ,           &
                          solx     , soly    , &
                          nbiter)
 
+      use ice_boundary, only: ice_HaloUpdate
+      use ice_domain, only: maskhalo_dyn, halo_info
+      use ice_fileunits, only: bfbflag
+      use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
+
       real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
          zetax2   , & ! zetax2 = 2*zeta (bulk viscosity)
          etax2        ! etax2  = 2*eta  (shear viscosity)
@@ -3125,6 +3132,9 @@ subroutine pgmres (zetax2   , etax2   , &
          Cb    , & ! seabed stress coefficient
          umassdti  ! mass of U-cell/dte (kg/m^2 s)
 
+      type (ice_halo), intent(in) :: &
+         halo_info_mask !  ghost cell update info for masked halo
+
       real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
          bx       , & ! Right hand side of the linear system (x components)
          by           ! Right hand side of the linear system (y components)
@@ -3293,14 +3303,29 @@ subroutine pgmres (zetax2   , etax2   , &
             call precondition(zetax2      , etax2          , &
                               Cb          , vrel           , &
                               umassdti    ,                  &
+                              halo_info_mask,                &
                               arnoldi_basis_x(:,:,:,initer), &
                               arnoldi_basis_y(:,:,:,initer), &
                               diagx       , diagy          , &
                               precond_type,                  &
                               workspace_x , workspace_y)
 
-            ! NOTE: halo updates for (workspace_x, workspace_y)
-            ! are skipped here for efficiency since this is just a preconditioner
+            ! Update workspace with boundary values
+            ! NOTE: skipped for efficiency since this is just a preconditioner
+            ! unless bfbflag is active            
+            if (bfbflag /= 'not') then
+               call stack_fields(workspace_x, workspace_y, fld2)
+               call ice_timer_start(timer_bound)
+               if (maskhalo_dyn) then
+                  call ice_HaloUpdate (fld2,               halo_info_mask, &
+                                       field_loc_NEcorner, field_type_vector)
+               else
+                  call ice_HaloUpdate (fld2,               halo_info, &
+                                       field_loc_NEcorner, field_type_vector)
+               endif
+               call ice_timer_stop(timer_bound)
+               call unstack_fields(fld2, workspace_x, workspace_y)
+            endif
 
             !$OMP PARALLEL DO PRIVATE(iblk)
             do iblk = 1, nblocks
@@ -3423,6 +3448,7 @@ subroutine pgmres (zetax2   , etax2   , &
          call precondition(zetax2      , etax2      , &
                            Cb          , vrel       , &
                            umassdti    ,              &
+                           halo_info_mask,            &
                            workspace_x , workspace_y, &
                            diagx       , diagy      , &
                            precond_type,              &
@@ -3489,6 +3515,7 @@ end subroutine pgmres
       subroutine precondition(zetax2      , etax2, &
                               Cb          , vrel , &
                               umassdti    ,        &
+                              halo_info_mask,      &
                               vx          , vy   , &
                               diagx       , diagy, &
                               precond_type,        &
@@ -3503,6 +3530,9 @@ subroutine precondition(zetax2      , etax2, &
          Cb    , & ! seabed stress coefficient
          umassdti  ! mass of U-cell/dte (kg/m^2 s)
 
+      type (ice_halo), intent(in) :: &
+         halo_info_mask !  ghost cell update info for masked halo
+
       real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
          vx       , & ! input vector (x components)
          vy           ! input vector (y components)
@@ -3564,6 +3594,7 @@ subroutine precondition(zetax2      , etax2, &
          call pgmres (zetax2,    etax2   , &
                       Cb       , vrel    , &
                       umassdti ,           &
+                      halo_info_mask     , &
                       vx       , vy      , &
                       diagx    , diagy   , &
                       tolerance, maxinner, &