dynamics: remove 'ice_HaloUpdate_vel' and introduce '(un)stack_veloci…

…ty_field'

The 'ice_HaloUpdate_vel' subroutine feels out of place in
'ice_dyn_shared', and moving it to 'ice_boundary' introduces other
problems, including circular dependencies.

Remove it, and introduce two simple subroutines, 'stack_velocity_field' and
'unstack_velocity_field', that are responsible to load the 'uvel' and
'vvel' arrays into the 'fld2' array used for the halo updates.

Use these new subroutines in ice_dyn_{evp,eap,vp}.

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -122,7 +122,8 @@ subroutine eap (dt)
       use ice_dyn_shared, only: fcor_blk, ndte, dtei, &
           denom1, uvel_init, vvel_init, arlx1i, &
           dyn_prep1, dyn_prep2, stepu, dyn_finish, &
-          basal_stress_coeff, basalstress, ice_HaloUpdate_vel
+          basal_stress_coeff, basalstress, &
+          stack_velocity_field, unstack_velocity_field
       use ice_flux, only: rdg_conv, strairxT, strairyT, &
           strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
           strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, &
@@ -172,6 +173,8 @@ subroutine eap (dt)
          umass    , & ! total mass of ice and snow (u grid)
          umassdti     ! mass of U-cell/dte (kg/m^2 s)
 
+      real (kind=dbl_kind), allocatable :: fld2(:,:,:,:)
+
       real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
          strtmp       ! stress combinations for momentum equation
 
@@ -195,6 +198,8 @@ subroutine eap (dt)
       ! Initialize
       !-----------------------------------------------------------------
 
+      allocate(fld2(nx_block,ny_block,2,max_blocks))
+
        ! This call is needed only if dt changes during runtime.
 !      call set_evp_parameters (dt)
 
@@ -373,7 +378,17 @@ subroutine eap (dt)
       endif
 
       ! velocities may have changed in dyn_prep2
-      call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+      call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
       !-----------------------------------------------------------------
       ! basal stress coefficients (landfast ice)
@@ -480,10 +495,21 @@ subroutine eap (dt)
          enddo
          !$TCXOMP END PARALLEL DO
 
-         call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+         call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
       enddo                     ! subcycling
 
+      deallocate(fld2)
       if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)
 
       !-----------------------------------------------------------------

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -94,7 +94,7 @@ subroutine evp (dt)
           ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d
       use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, &
           ice_dyn_evp_1d_copyout
-      use ice_dyn_shared, only: kevp_kernel, ice_HaloUpdate_vel
+      use ice_dyn_shared, only: kevp_kernel, stack_velocity_field, unstack_velocity_field
 
       real (kind=dbl_kind), intent(in) :: &
          dt      ! time step
@@ -127,6 +127,8 @@ subroutine evp (dt)
          umass    , & ! total mass of ice and snow (u grid)
          umassdti     ! mass of U-cell/dte (kg/m^2 s)
 
+      real (kind=dbl_kind), allocatable :: fld2(:,:,:,:)
+
       real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
          strtmp       ! stress combinations for momentum equation
 
@@ -152,6 +154,8 @@ subroutine evp (dt)
       ! Initialize
       !-----------------------------------------------------------------
 
+      allocate(fld2(nx_block,ny_block,2,max_blocks))
+
        ! This call is needed only if dt changes during runtime.
 !      call set_evp_parameters (dt)
 
@@ -316,7 +320,17 @@ subroutine evp (dt)
       endif
 
       ! velocities may have changed in dyn_prep2
-      call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+      call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
       !-----------------------------------------------------------------
       ! basal stress coefficients (landfast ice)
@@ -429,12 +443,23 @@ subroutine evp (dt)
          enddo
          !$TCXOMP END PARALLEL DO
 
-         call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+         call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
       enddo                     ! subcycling
       endif  ! kevp_kernel
       call ice_timer_stop(timer_evp_2d)
 
+      deallocate(fld2)
       if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)
 
       ! Force symmetry across the tripole seam

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -24,7 +24,8 @@ module ice_dyn_shared
       private
       public :: init_evp, set_evp_parameters, stepu, principal_stress, &
                 dyn_prep1, dyn_prep2, dyn_finish, basal_stress_coeff,  &
-                alloc_dyn_shared, deformations, strain_rates, ice_HaloUpdate_vel
+                alloc_dyn_shared, deformations, strain_rates, &
+                stack_velocity_field, unstack_velocity_field
 
       ! namelist parameters
 
@@ -94,9 +95,6 @@ module ice_dyn_shared
                          ! see keel data from Amundrud et al. 2004 (JGR)
          u0 = 5e-5_dbl_kind ! residual velocity for basal stress (m/s)
 
-      real (kind=dbl_kind), allocatable :: &
-         fld2(:,:,:,:)    ! work array for boundary updates
-
 !=======================================================================
 
       contains
@@ -112,7 +110,6 @@ subroutine alloc_dyn_shared
       allocate( &
          uvel_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep
          vvel_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep
-         fld2      (nx_block,ny_block,2,max_blocks), & ! work array for boundary updates
          stat=ierr)
       if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory')
 
@@ -1183,29 +1180,25 @@ end subroutine strain_rates
 
 !=======================================================================
 
-! Perform a halo update for the velocity field
-! author: Philippe Blain, ECCC
-
-      subroutine ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+! Load velocity components into array for boundary updates
 
-      use ice_boundary, only: ice_HaloUpdate, ice_halo
-      use ice_constants, only: field_loc_NEcorner, field_type_vector
-      use ice_domain, only: halo_info, maskhalo_dyn, nblocks
-      use ice_timers, only: timer_bound, ice_timer_start, ice_timer_stop
+   subroutine stack_velocity_field(uvel, vvel, fld2)
 
-      type (ice_halo), intent(in) :: &
-         halo_info_mask !  ghost cell update info for masked halo
+      use ice_domain, only: nblocks
 
-      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
+      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
          uvel    , & ! u components of velocity vector
          vvel        ! v components of velocity vector
 
+      real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks), intent(out) :: &
+         fld2        ! work array for boundary updates
+
       ! local variables
 
       integer (kind=int_kind) :: &
          iblk        ! block index
 
-      character(len=*), parameter :: subname = '(ice_HaloUpdate_vel)'
+      character(len=*), parameter :: subname = '(stack_velocity_field)'
 
       ! load velocity into array for boundary updates
       !$OMP PARALLEL DO PRIVATE(iblk)
@@ -1215,28 +1208,42 @@ subroutine ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
       enddo
       !$OMP END PARALLEL DO
 
-      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)
+      end subroutine stack_velocity_field
+
+!=======================================================================
+
+! Unload velocity components from array after boundary updates
+
+   subroutine unstack_velocity_field(fld2, uvel, vvel)
+
+      use ice_domain, only: nblocks
+
+      real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks), intent(in) :: &
+         fld2        ! work array for boundary updates
 
-      ! Unload
+      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(out) :: &
+         uvel    , & ! u components of velocity vector
+         vvel        ! v components of velocity vector
+
+      ! local variables
+
+      integer (kind=int_kind) :: &
+         iblk        ! block index
+
+      character(len=*), parameter :: subname = '(unstack_velocity_field)'
+
+      ! Unload velocity from array after boundary updates
       !$OMP PARALLEL DO PRIVATE(iblk)
       do iblk = 1, nblocks
          uvel(:,:,iblk) = fld2(:,:,1,iblk)
          vvel(:,:,iblk) = fld2(:,:,2,iblk)
       enddo
       !$OMP END PARALLEL DO
 
-      end subroutine ice_HaloUpdate_vel
+      end subroutine unstack_velocity_field
 
 !=======================================================================
-
+      
       end module ice_dyn_shared
 
 !=======================================================================

cicecore/cicedynB/dynamics/ice_dyn_vp.F90

@@ -47,7 +47,8 @@ module ice_dyn_vp
       use ice_domain_size, only: max_blocks
       use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, &
           ecci, cosw, sinw, fcor_blk, uvel_init,  &
-          vvel_init, basal_stress_coeff, basalstress, Ktens, ice_HaloUpdate_vel
+          vvel_init, basal_stress_coeff, basalstress, Ktens, &
+          stack_velocity_field,  unstack_velocity_field
       use ice_fileunits, only: nu_diag
       use ice_flux, only: fm
       use ice_global_reductions, only: global_sum, global_allreduce_sum
@@ -102,6 +103,9 @@ module ice_dyn_vp
          indxui(:,:)  , & ! compressed index in i-direction
          indxuj(:,:)      ! compressed index in j-direction
 
+      real (kind=dbl_kind), allocatable :: & 
+         fld2(:,:,:,:)    ! work array for boundary updates
+
 !=======================================================================
 
       contains
@@ -145,6 +149,7 @@ subroutine init_vp (dt)
                indxtj(nx_block*ny_block, max_blocks), &
                indxui(nx_block*ny_block, max_blocks), &
                indxuj(nx_block*ny_block, max_blocks))
+      allocate(fld2(nx_block,ny_block,2,max_blocks))
 
       ! Redefine tinyarea using min_strain_rate
 
@@ -433,7 +438,17 @@ subroutine implicit_solver (dt)
       endif
 
       ! velocities may have changed in dyn_prep2
-      call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+      call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
       !-----------------------------------------------------------------
       ! basal stress coefficients (landfast ice)
@@ -664,13 +679,13 @@ subroutine anderson_solver (icellt  , icellu, &
       use ice_blocks, only: nx_block, ny_block
       use ice_boundary, only: ice_HaloUpdate
       use ice_constants, only: c1
-      use ice_domain, only: maskhalo_dyn
+      use ice_domain, only: maskhalo_dyn, halo_info
       use ice_domain_size, only: max_blocks
       use ice_flux, only:   uocn, vocn, fm, Tbu
       use ice_grid, only: dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, &
           uarear, tinyarea
       use ice_state, only: uvel, vvel, strength
-      use ice_timers, only: ice_timer_start, ice_timer_stop
+      use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
 
       integer (kind=int_kind), intent(in) :: &
          ntot         ! size of problem for Anderson
@@ -1069,7 +1084,17 @@ subroutine anderson_solver (icellt  , icellu, &
                              uvel (:,:,:), vvel (:,:,:))
 
          ! Do halo update so that halo cells contain up to date info for advection
-         call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask)
+         call stack_velocity_field(uvel, vvel, 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_velocity_field(fld2, uvel, vvel)
 
          ! Compute "progress" residual norm
          !$OMP PARALLEL DO PRIVATE(iblk)
@@ -2644,6 +2669,10 @@ subroutine fgmres (zetaD    ,           &
                          solx     , soly    , &
                          nbiter)
 
+      use ice_boundary, only: ice_HaloUpdate
+      use ice_domain, only: maskhalo_dyn, halo_info
+      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) :: &
          zetaD   ! zetaD = 2*zeta (viscous coefficient)
 
@@ -2839,8 +2868,18 @@ subroutine fgmres (zetaD    ,           &
             orig_basis_y(:,:,:,initer) = workspace_y
 
             ! Update workspace with boundary values
-            call ice_HaloUpdate_vel(workspace_x, workspace_y,  &
-                                    halo_info_mask)
+            call stack_velocity_field(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_velocity_field(fld2, workspace_x, workspace_y)
+
             !$OMP PARALLEL DO PRIVATE(iblk)
             do iblk = 1, nblocks
                call matvec (nx_block               , ny_block             , &