Refactor strocnxT, strocnyT implementation (CICE-Consortium#764)

* Refactor strocnxT, strocnyT implementation - add aiU to ice_state - migrate computation of strocnxT and strocnyT to ice_step, needed for thermodynamics, better code reuse. - add strocnxT_sf, strocnyT_sf as coupling field, could be computed differently than the thermodynanics version. The _sf field computation should be in scale fluxes, but because scale_fluxes is called on a block and the _sf fields require a halo update among other things, the computation can't be done in scale_fluxes. - Update the coupling layers to use the _sf version of the fields. - CICE-Consortium#761 suggests the values of strocnxT, strocnyT should not be scaled for use in thermodynamics. This commit does not make that change yet, but allows for that change to be made easily. - These changes are bit-for-bit for a full test suite on cheyenne with 3 compilers. * Update computation of strocnxT, strocnyT passed into icepack_step_therm1 - No longer divided by aice - strocnxT_sf, strocnyT_sf are still computed in the same way as before * Rename strocn[x,y]T_sf to strocn[x,y]T_iavg Revert strocn[x,y]T passed into thermodynamics to be the version divided by aice, specifically strocn[x,y]T_iavg. This is identical to earlier implementations.
NOAA-EMC · DeniseWorthen · Mar 29, 2023 · Jul 15, 2022 · Jul 15, 2022 · Jul 15, 2022
commit 8c6ba0437b00a53fcd935fd717f931e9a5c5a6cf
diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90
@@ -138,15 +138,15 @@ subroutine eap (dt)
       use ice_flux, only: rdg_conv, strairxT, strairyT, &
           strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
           strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, &
-          strocnxT, strocnyT, strax, stray, &
+          strax, stray, &
           TbU, hwater, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4
       use ice_grid, only: tmask, umask, dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
           tarear, uarear, grid_average_X2Y, iceumask, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
-      use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, &
+      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
           ice_timer_start, ice_timer_stop
@@ -182,7 +182,6 @@ subroutine eap (dt)
          wateryU    , & ! for ocean stress calculation, y (m/s)
          forcexU    , & ! work array: combined atm stress and ocn tilt, x
          forceyU    , & ! work array: combined atm stress and ocn tilt, y
-         aiU        , & ! ice fraction on u-grid
          umass      , & ! total mass of ice and snow (u grid)
          umassdti       ! mass of U-cell/dte (kg/m^2 s)
 
@@ -205,10 +204,6 @@ subroutine eap (dt)
       type (block) :: &
          this_block     ! block information for current block
 
-      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
-         work1      , & ! temporary
-         work2          ! temporary
-
       character(len=*), parameter :: subname = '(eap)'
 
       call ice_timer_start(timer_dynamics) ! dynamics
@@ -567,27 +562,6 @@ subroutine eap (dt)
       enddo
       !$OMP END PARALLEL DO
 
-      ! strocn computed on U, N, E as needed. Map strocn U divided by aiU to T
-      ! TODO: This should be done elsewhere as part of generalization?
-      ! conservation requires aiU be divided before averaging
-      work1 = c0
-      work2 = c0
-      !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
-      do iblk = 1, nblocks
-      do ij = 1, icellu(iblk)
-         i = indxui(ij,iblk)
-         j = indxuj(ij,iblk)
-         work1(i,j,iblk) = strocnxU(i,j,iblk)/aiU(i,j,iblk)
-         work2(i,j,iblk) = strocnyU(i,j,iblk)/aiU(i,j,iblk)
-      enddo
-      enddo
-      call ice_HaloUpdate (work1,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call ice_HaloUpdate (work2,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call grid_average_X2Y('F', work1, 'U', strocnxT, 'T')    ! shift
-      call grid_average_X2Y('F', work2, 'U', strocnyT, 'T')
-
       call ice_timer_stop(timer_dynamics)    ! dynamics
 
       end subroutine eap

diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90
@@ -85,7 +85,7 @@ subroutine evp (dt)
       use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
           strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, &
-          strocnxT, strocnyT, strax, stray, &
+          strax, stray, &
           TbU, hwater, &
           strairxN, strairyN, fmN, &
           strtltxN, strtltyN, strocnxN, strocnyN, strintxN, strintyN, taubxN, taubyN, &
@@ -106,7 +106,7 @@ subroutine evp (dt)
           tarear, uarear, earear, narear, grid_average_X2Y, uarea, &
           grid_type, grid_ice, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
-      use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, &
+      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, uvelN, vvelN, &
           uvelE, vvelE, divu, shear, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
@@ -155,7 +155,6 @@ subroutine evp (dt)
          wateryU  , & ! for ocean stress calculation, y (m/s)
          forcexU  , & ! work array: combined atm stress and ocn tilt, x
          forceyU  , & ! work array: combined atm stress and ocn tilt, y
-         aiU      , & ! ice fraction on u-grid
          umass    , & ! total mass of ice and snow (u grid)
          umassdti     ! mass of U-cell/dte (kg/m^2 s)
 
@@ -217,10 +216,6 @@ subroutine evp (dt)
       type (block) :: &
          this_block   ! block information for current block
 
-      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
-         work1, &     ! temporary
-         work2        ! temporary
-
       logical (kind=log_kind), save :: &
          first_time = .true. ! first time logical
 
@@ -1326,29 +1321,6 @@ subroutine evp (dt)
 
       endif
 
-      ! strocn computed on U, N, E as needed. Map strocn U divided by aiU to T
-      ! TODO: This should be done elsewhere as part of generalization?
-      ! TODO: Rename strocn[x,y]T since it's different than strocn[x,y][U,N,E]
-      ! conservation requires aiU be divided before averaging
-      work1 = c0
-      work2 = c0
-      !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) SCHEDULE(runtime)
-      do iblk = 1, nblocks
-      do ij = 1, icellu(iblk)
-         i = indxui(ij,iblk)
-         j = indxuj(ij,iblk)
-         work1(i,j,iblk) = strocnxU(i,j,iblk)/aiU(i,j,iblk)
-         work2(i,j,iblk) = strocnyU(i,j,iblk)/aiU(i,j,iblk)
-      enddo
-      enddo
-      !$OMP END PARALLEL DO
-      call ice_HaloUpdate (work1,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call ice_HaloUpdate (work2,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call grid_average_X2Y('F', work1, 'U', strocnxT, 'T')    ! shift
-      call grid_average_X2Y('F', work2, 'U', strocnyT, 'T')
-
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
          call grid_average_X2Y('S', strintxE, 'E', strintxU, 'U')    ! diagnostic
          call grid_average_X2Y('S', strintyN, 'N', strintyU, 'U')    ! diagnostic

diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90
@@ -171,15 +171,15 @@ subroutine implicit_solver (dt)
       use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
           strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, &
-          strocnxT, strocnyT, strax, stray, &
+          strax, stray, &
           TbU, hwater, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4
       use ice_grid, only: tmask, umask, dxT, dyT, cxp, cyp, cxm, cym, &
           tarear, grid_type, grid_average_X2Y, iceumask, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
-      use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, &
+      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
           ice_timer_start, ice_timer_stop
@@ -210,7 +210,6 @@ subroutine implicit_solver (dt)
          Cb       , & ! seabed stress coefficient
          fpresx   , & ! fixed point residual vector, x components: fx = uvel - uprev_k
          fpresy   , & ! fixed point residual vector, y components: fy = vvel - vprev_k
-         aiU      , & ! ice fraction on u-grid
          umass    , & ! total mass of ice and snow (u grid)
          umassdti     ! mass of U-cell/dte (kg/m^2 s)
 
@@ -234,10 +233,6 @@ subroutine implicit_solver (dt)
       real (kind=dbl_kind), allocatable :: &
          sol(:)          ! solution vector
 
-      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
-         work1, &      ! temporary
-         work2         ! temporary
-
       character(len=*), parameter :: subname = '(implicit_solver)'
 
       call ice_timer_start(timer_dynamics) ! dynamics
@@ -640,27 +635,6 @@ subroutine implicit_solver (dt)
       enddo
       !$OMP END PARALLEL DO
 
-      ! strocn computed on U, N, E as needed. Map strocn U divided by aiU to T
-      ! TODO: This should be done elsewhere as part of generalization?
-      ! conservation requires aiU be divided before averaging
-      work1 = c0
-      work2 = c0
-      !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
-      do iblk = 1, nblocks
-      do ij = 1, icellu(iblk)
-         i = indxui(ij,iblk)
-         j = indxuj(ij,iblk)
-         work1(i,j,iblk) = strocnxU(i,j,iblk)/aiU(i,j,iblk)
-         work2(i,j,iblk) = strocnyU(i,j,iblk)/aiU(i,j,iblk)
-      enddo
-      enddo
-      call ice_HaloUpdate (work1,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call ice_HaloUpdate (work2,              halo_info, &
-                           field_loc_NEcorner, field_type_vector)
-      call grid_average_X2Y('F',work1,'U',strocnxT,'T')    ! shift
-      call grid_average_X2Y('F',work2,'U',strocnyT,'T')
-
 ! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport
 ! commented out in order to focus on EVP for now within the cdgrid
 ! should be used when routine is ready

diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90
@@ -35,6 +35,16 @@ module ice_flux
       ! Dynamics component
       ! All variables are assumed to be on the atm or ocn thermodynamic
       ! grid except as noted
+      !
+      ! scale_fluxes divides several of these by aice "in place", so
+      ! the state of some of these variables is not well defined.  In the
+      ! future, we need to refactor and add "_iavg" versions of the
+      ! fields to clearly differentiate fields that have been divided
+      ! by aice and others that are not.  The challenge is that we need
+      ! to go thru each field carefully to see which version is used.
+      ! For instance, in diagnostics, there are places where these
+      ! fields are multiplied by aice to compute things properly.
+      ! strocn[x,y]T_iavg is the first field defined using _iavg.
       !-----------------------------------------------------------------
 
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
@@ -56,8 +66,8 @@ module ice_flux
 
        ! out to ocean          T-cell (kg/m s^2)
        ! Note, CICE_IN_NEMO uses strocnx and strocny for coupling
-         strocnxT, & ! ice-ocean stress, x-direction at T points, per ice fraction
-         strocnyT    ! ice-ocean stress, y-direction at T points, per ice fraction
+         strocnxT_iavg, & ! ice-ocean stress, x-direction at T points, per ice fraction (scaled flux)
+         strocnyT_iavg    ! ice-ocean stress, y-direction at T points, per ice fraction (scaled flux)
 
        ! diagnostic
 
@@ -389,8 +399,8 @@ subroutine alloc_flux
          hwater     (nx_block,ny_block,max_blocks), & ! water depth for seabed stress calc (landfast ice)
          strairxT   (nx_block,ny_block,max_blocks), & ! stress on ice by air, x-direction
          strairyT   (nx_block,ny_block,max_blocks), & ! stress on ice by air, y-direction
-         strocnxT   (nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction
-         strocnyT   (nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction
+         strocnxT_iavg(nx_block,ny_block,max_blocks), & ! ice-ocean stress, x-direction, per ice area
+         strocnyT_iavg(nx_block,ny_block,max_blocks), & ! ice-ocean stress, y-direction, per ice area
          sig1       (nx_block,ny_block,max_blocks), & ! normalized principal stress component
          sig2       (nx_block,ny_block,max_blocks), & ! normalized principal stress component
          sigP       (nx_block,ny_block,max_blocks), & ! internal ice pressure (N/m)
@@ -765,8 +775,8 @@ subroutine init_coupler_flux
       ! fluxes sent to ocean
       !-----------------------------------------------------------------
 
-      strocnxT(:,:,:) = c0    ! ice-ocean stress, x-direction (T-cell)
-      strocnyT(:,:,:) = c0    ! ice-ocean stress, y-direction (T-cell)
+      strocnxT_iavg (:,:,:) = c0 ! ice-ocean stress, x-direction (T-cell)
+      strocnyT_iavg (:,:,:) = c0 ! ice-ocean stress, y-direction (T-cell)
       fresh   (:,:,:) = c0
       fsalt   (:,:,:) = c0
       fpond   (:,:,:) = c0

diff --git a/cicecore/cicedynB/general/ice_state.F90 b/cicecore/cicedynB/general/ice_state.F90
@@ -54,7 +54,8 @@ module ice_state
 
       real (kind=dbl_kind), dimension(:,:,:), allocatable, &
          public :: &
-         aice  , & ! concentration of ice
+         aice  , & ! concentration of ice on T grid
+         aiU   , & ! concentration of ice on U grid
          vice  , & ! volume per unit area of ice          (m)
          vsno      ! volume per unit area of snow         (m)
 
@@ -151,7 +152,8 @@ subroutine alloc_state
           file=__FILE__, line=__LINE__)
 
       allocate ( &
-         aice      (nx_block,ny_block,max_blocks) , & ! concentration of ice
+         aice      (nx_block,ny_block,max_blocks) , & ! concentration of ice T grid
+         aiU       (nx_block,ny_block,max_blocks) , & ! concentration of ice U grid
          vice      (nx_block,ny_block,max_blocks) , & ! volume per unit area of ice (m)
          vsno      (nx_block,ny_block,max_blocks) , & ! volume per unit area of snow (m)
          aice0     (nx_block,ny_block,max_blocks) , & ! concentration of open water

diff --git a/cicecore/cicedynB/general/ice_step_mod.F90 b/cicecore/cicedynB/general/ice_step_mod.F90
@@ -12,7 +12,14 @@
       module ice_step_mod
 
       use ice_kinds_mod
+      use ice_blocks, only: block, get_block
+      use ice_blocks, only: nx_block, ny_block
       use ice_constants, only: c0, c1, c1000, c4, p25
+      use ice_constants, only: field_loc_center, field_loc_NEcorner, &
+          field_loc_Nface, field_loc_Eface, &
+          field_type_scalar, field_type_vector
+      use ice_domain, only: halo_info, nblocks, blocks_ice
+      use ice_domain_size, only: max_blocks
       use ice_exit, only: abort_ice
       use ice_fileunits, only: nu_diag
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
@@ -95,8 +102,6 @@ end subroutine step_prep
 
       subroutine prep_radiation (iblk)
 
-      use ice_blocks, only: block, get_block
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr
       use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
           alvdr_ai, alvdf_ai, alidr_ai, alidf_ai, &
@@ -188,14 +193,9 @@ subroutine step_therm1 (dt, iblk)
           hfreebd, hdraft, hridge, distrdg, hkeel, dkeel, lfloe, dfloe, &
           fswsfcn, fswintn, Sswabsn, Iswabsn, meltsliqn, meltsliq, &
           fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
-      use ice_blocks, only: block, get_block
-#ifdef CICE_IN_NEMO
-      use ice_blocks, only: nx_block, ny_block
-#endif
       use ice_calendar, only: yday
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero
-      use ice_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fbot, Tbot, Tsnice, &
+      use ice_flux, only: frzmlt, sst, Tf, strocnxT_iavg, strocnyT_iavg, rside, fbot, Tbot, Tsnice, &
           meltsn, melttn, meltbn, congeln, snoicen, uatmT, vatmT, fside, &
           wind, rhoa, potT, Qa, zlvl, zlvs, strax, stray, flatn, fsensn, fsurfn, fcondtopn, &
           flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, fsloss, &
@@ -448,8 +448,8 @@ subroutine step_therm1 (dt, iblk)
                       sst          = sst         (i,j,  iblk), &
                       sss          = sss         (i,j,  iblk), &
                       Tf           = Tf          (i,j,  iblk), &
-                      strocnxT     = strocnxT    (i,j,  iblk), &
-                      strocnyT     = strocnyT    (i,j,  iblk), &
+                      strocnxT    = strocnxT_iavg(i,j,  iblk), &
+                      strocnyT    = strocnyT_iavg(i,j,  iblk), &
                       fbot         = fbot        (i,j,  iblk), &
                       Tbot         = Tbot        (i,j,  iblk), &
                       Tsnice       = Tsnice      (i,j,  iblk), &
@@ -602,9 +602,7 @@ subroutine step_therm2 (dt, iblk)
           wave_spectrum, wavefreq, dwavefreq, &
           first_ice, bgrid, cgrid, igrid, floe_rad_c, floe_binwidth, &
           d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld
-      use ice_blocks, only: block, get_block
       use ice_calendar, only: yday
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr, nblyr, nfsd
       use ice_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset, &
           update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, &
@@ -744,8 +742,6 @@ end subroutine step_therm2
 
       subroutine update_state (dt, daidt, dvidt, dagedt, offset)
 
-      use ice_blocks, only: nx_block, ny_block
-      use ice_domain, only: nblocks
       use ice_domain_size, only: ncat
 !     use ice_grid, only: tmask
       use ice_state, only: aicen, trcrn, vicen, vsnon, &
@@ -862,8 +858,6 @@ subroutine step_dyn_wave (dt)
 
       use ice_arrays_column, only: wave_spectrum, &
           d_afsd_wave, floe_rad_l, floe_rad_c, wavefreq, dwavefreq
-      use ice_blocks, only: block, get_block
-      use ice_domain, only: blocks_ice, nblocks
       use ice_domain_size, only: ncat, nfsd, nfreq
       use ice_state, only: trcrn, aicen, aice, vice
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_column, &
@@ -937,18 +931,34 @@ end subroutine step_dyn_wave
 
       subroutine step_dyn_horiz (dt)
 
+      use ice_boundary, only: ice_HaloUpdate
       use ice_dyn_evp, only: evp
       use ice_dyn_eap, only: eap
       use ice_dyn_vp, only: implicit_solver
       use ice_dyn_shared, only: kdyn
+      use ice_flux, only: strocnxU, strocnyU, strocnxT_iavg, strocnyT_iavg
       use ice_flux, only: init_history_dyn
       use ice_grid, only: grid_average_X2Y
-      use ice_state, only: uvel, vvel, uvelT, vvelT
+      use ice_state, only: aiU, uvel, vvel, uvelT, vvelT
       use ice_transport_driver, only: advection, transport_upwind, transport_remap
 
       real (kind=dbl_kind), intent(in) :: &
          dt      ! dynamics time step
 
+      ! local variables
+
+      type (block) :: &
+         this_block      ! block information for current block
+
+      integer (kind=int_kind) :: &
+         ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
+         iblk,            & ! block index
+         i, j               ! horizontal indices
+
+      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
+         work1, &     ! temporary
+         work2        ! temporary
+
       character(len=*), parameter :: subname = '(step_dyn_horiz)'
 
       call init_history_dyn     ! initialize dynamic history variables
@@ -969,6 +979,38 @@ subroutine step_dyn_horiz (dt)
       call grid_average_X2Y('A', uvel, 'U', uvelT, 'T')
       call grid_average_X2Y('A', vvel, 'U', vvelT, 'T')
 
+      !-----------------------------------------------------------------
+      ! Compute strocnxT_iavg, strocnyT_iavg for thermo and coupling
+      !-----------------------------------------------------------------
+
+      ! strocn computed on U, N, E as needed. Map strocn U divided by aiU to T
+      ! conservation requires aiU be divided before averaging
+      work1 = c0
+      work2 = 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)
+         ilo = this_block%ilo
+         ihi = this_block%ihi
+         jlo = this_block%jlo
+         jhi = this_block%jhi
+         do j = jlo, jhi
+         do i = ilo, ihi
+            if (aiU(i,j,iblk) /= c0) then
+               work1(i,j,iblk) = strocnxU(i,j,iblk)/aiU(i,j,iblk)
+               work2(i,j,iblk) = strocnyU(i,j,iblk)/aiU(i,j,iblk)
+            endif
+        enddo
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+      call ice_HaloUpdate (work1,              halo_info, &
+                           field_loc_NEcorner, field_type_vector)
+      call ice_HaloUpdate (work2,              halo_info, &
+                           field_loc_NEcorner, field_type_vector)
+      call grid_average_X2Y('F', work1, 'U', strocnxT_iavg, 'T')    ! shift
+      call grid_average_X2Y('F', work2, 'U', strocnyT_iavg, 'T')
+
       !-----------------------------------------------------------------
       ! Horizontal ice transport
       !-----------------------------------------------------------------
@@ -991,8 +1033,6 @@ end subroutine step_dyn_horiz
       subroutine step_dyn_ridge (dt, ndtd, iblk)
 
       use ice_arrays_column, only: hin_max, fzsal, first_ice
-      use ice_blocks, only: block, get_block
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr
       use ice_flux, only: &
           rdg_conv, rdg_shear, dardg1dt, dardg2dt, &
@@ -1115,9 +1155,7 @@ end subroutine step_dyn_ridge
 
       subroutine step_snow (dt, iblk)
 
-      use ice_blocks, only: block, get_block
       use ice_calendar, only: nstreams
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nslyr, nilyr
       use ice_flux, only: snwcnt, wind, fresh, fhocn, fsloss, fsnow
       use ice_state, only: trcrn, vsno, vsnon, vicen, aicen, aice
@@ -1233,9 +1271,7 @@ subroutine step_radiation (dt, iblk)
           alvdrn, alidrn, alvdfn, alidfn, apeffn, trcrn_sw, snowfracn, &
           kaer_tab, waer_tab, gaer_tab, kaer_bc_tab, waer_bc_tab, &
           gaer_bc_tab, bcenh, swgrid, igrid
-      use ice_blocks, only: block, get_block
       use ice_calendar, only: calendar_type, days_per_year, nextsw_cday, yday, msec
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, n_aero, nilyr, nslyr, n_zaero, n_algae, nblyr
       use ice_flux, only: swvdr, swvdf, swidr, swidf, coszen, fsnow
       use ice_grid, only: TLAT, TLON, tmask
@@ -1423,7 +1459,6 @@ end subroutine step_radiation
       subroutine ocean_mixed_layer (dt, iblk)
 
       use ice_arrays_column, only: Cdn_atm, Cdn_atm_ratio
-      use ice_blocks, only: nx_block, ny_block
       use ice_flux, only: sst, Tf, Qa, uatmT, vatmT, wind, potT, rhoa, zlvl, &
            frzmlt, fhocn, fswthru, flw, flwout_ocn, fsens_ocn, flat_ocn, evap_ocn, &
            alvdr_ocn, alidr_ocn, alvdf_ocn, alidf_ocn, swidf, swvdf, swidr, swvdr, &
@@ -1580,8 +1615,6 @@ subroutine biogeochemistry (dt, iblk)
                            snow_bio_net, fswthrun, Rayleigh_criteria, &
                            ocean_bio_all, sice_rho, fzsal, fzsal_g, &
                            bgrid, igrid, icgrid, cgrid
-      use ice_blocks, only: block, get_block
-      use ice_domain, only: blocks_ice
       use ice_domain_size, only: nblyr, nilyr, nslyr, n_algae, n_zaero, ncat, &
                                  n_doc, n_dic,  n_don, n_fed, n_fep
       use ice_flux, only: meltbn, melttn, congeln, snoicen, &

diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90
@@ -56,7 +56,7 @@ subroutine dumpfile(filename_spec)
       use ice_domain, only: nblocks
       use ice_domain_size, only: nilyr, nslyr, ncat, max_blocks
       use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
-          strocnxT, strocnyT, sst, frzmlt, &
+          strocnxT_iavg, strocnyT_iavg, sst, frzmlt, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4, &
@@ -175,8 +175,8 @@ subroutine dumpfile(filename_spec)
       !-----------------------------------------------------------------
       ! ocean stress (for bottom heat flux in thermo)
       !-----------------------------------------------------------------
-      call write_restart_field(nu_dump,0,strocnxT,'ruf8','strocnxT',1,diag)
-      call write_restart_field(nu_dump,0,strocnyT,'ruf8','strocnyT',1,diag)
+      call write_restart_field(nu_dump,0,strocnxT_iavg,'ruf8','strocnxT',1,diag)
+      call write_restart_field(nu_dump,0,strocnyT_iavg,'ruf8','strocnyT',1,diag)
 
       !-----------------------------------------------------------------
       ! internal stress
@@ -277,7 +277,7 @@ subroutine restartfile (ice_ic)
       use ice_domain_size, only: nilyr, nslyr, ncat, &
           max_blocks
       use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
-          strocnxT, strocnyT, sst, frzmlt, &
+          strocnxT_iavg, strocnyT_iavg, sst, frzmlt, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4, &
@@ -431,9 +431,9 @@ subroutine restartfile (ice_ic)
       if (my_task == master_task) &
            write(nu_diag,*) 'min/max ocean stress components'
 
-      call read_restart_field(nu_restart,0,strocnxT,'ruf8', &
+      call read_restart_field(nu_restart,0,strocnxT_iavg,'ruf8', &
            'strocnxT',1,diag,field_loc_center, field_type_vector)
-      call read_restart_field(nu_restart,0,strocnyT,'ruf8', &
+      call read_restart_field(nu_restart,0,strocnyT_iavg,'ruf8', &
            'strocnyT',1,diag,field_loc_center, field_type_vector)
 
       !-----------------------------------------------------------------
@@ -711,7 +711,7 @@ subroutine restartfile_v4 (ice_ic)
       use ice_domain_size, only: nilyr, nslyr, ncat, nx_global, ny_global, &
           max_blocks
       use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
-          strocnxT, strocnyT, sst, frzmlt, &
+          strocnxT_iavg, strocnyT_iavg, sst, frzmlt, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4
@@ -876,9 +876,9 @@ subroutine restartfile_v4 (ice_ic)
       if (my_task == master_task) &
            write(nu_diag,*) 'min/max ocean stress components'
 
-      call ice_read(nu_restart,0,strocnxT,'ruf8',diag, &
+      call ice_read(nu_restart,0,strocnxT_iavg,'ruf8',diag, &
                        field_loc_center, field_type_vector)
-      call ice_read(nu_restart,0,strocnyT,'ruf8',diag, &
+      call ice_read(nu_restart,0,strocnyT_iavg,'ruf8',diag, &
                        field_loc_center, field_type_vector)
 
       !-----------------------------------------------------------------

diff --git a/cicecore/drivers/mct/cesm1/ice_import_export.F90 b/cicecore/drivers/mct/cesm1/ice_import_export.F90
@@ -9,7 +9,7 @@ module ice_import_export
   use ice_constants     , only: field_type_vector, c100
   use ice_constants     , only: p001, p5
   use ice_blocks        , only: block, get_block, nx_block, ny_block
-  use ice_flux          , only: strairxT, strairyT, strocnxT, strocnyT
+  use ice_flux          , only: strairxT, strairyT, strocnxT_iavg, strocnyT_iavg
   use ice_flux          , only: alvdr, alidr, alvdf, alidf, Tref, Qref, Uref
   use ice_flux          , only: flat, fsens, flwout, evap, fswabs, fhocn, fswthru
   use ice_flux          , only: fresh, fsalt, zlvl, uatm, vatm, potT, Tair, Qa
@@ -571,8 +571,8 @@ subroutine ice_export( i2x )
                              + workx*sin(ANGLET(i,j,iblk))
 
              ! ice/ocean stress (on POP T-grid:  convert to lat-lon)
-             workx = -strocnxT(i,j,iblk)                            ! N/m^2
-             worky = -strocnyT(i,j,iblk)                            ! N/m^2
+             workx = -strocnxT_iavg(i,j,iblk)                       ! N/m^2
+             worky = -strocnyT_iavg(i,j,iblk)                       ! N/m^2
              tauxo(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) &
                              - worky*sin(ANGLET(i,j,iblk))
              tauyo(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &

diff --git a/cicecore/drivers/mct/cesm1/ice_prescribed_mod.F90 b/cicecore/drivers/mct/cesm1/ice_prescribed_mod.F90
@@ -585,11 +585,11 @@ subroutine ice_prescribed_phys
    ! set non-computed fluxes, ice velocities, ice-ocn stresses to zero
    !--------------------------------------------------------------------
 
-   frzmlt    (:,:,:) = c0
-   uvel      (:,:,:) = c0
-   vvel      (:,:,:) = c0
-   strocnxT  (:,:,:) = c0
-   strocnyT  (:,:,:) = c0
+   frzmlt       (:,:,:) = c0
+   uvel         (:,:,:) = c0
+   vvel         (:,:,:) = c0
+   strocnxT_iavg(:,:,:) = c0
+   strocnyT_iavg(:,:,:) = c0
 
    !-----------------------------------------------------------------
    ! other atm and ocn fluxes

diff --git a/cicecore/drivers/nuopc/cmeps/ice_import_export.F90 b/cicecore/drivers/nuopc/cmeps/ice_import_export.F90
@@ -10,7 +10,7 @@ module ice_import_export
   use ice_domain         , only : nblocks, blocks_ice, halo_info, distrb_info
   use ice_domain_size    , only : nx_global, ny_global, block_size_x, block_size_y, max_blocks, ncat
   use ice_exit           , only : abort_ice
-  use ice_flux           , only : strairxT, strairyT, strocnxT, strocnyT
+  use ice_flux           , only : strairxT, strairyT, strocnxT_iavg, strocnyT_iavg
   use ice_flux           , only : alvdr, alidr, alvdf, alidf, Tref, Qref, Uref
   use ice_flux           , only : flat, fsens, flwout, evap, fswabs, fhocn, fswthru
   use ice_flux           , only : fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf
@@ -914,8 +914,8 @@ subroutine ice_export( exportState, rc )
              tauya(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) + workx*sin(ANGLET(i,j,iblk))
 
              ! ice/ocean stress (on POP T-grid:  convert to lat-lon)
-             workx = -strocnxT(i,j,iblk)                            ! N/m^2
-             worky = -strocnyT(i,j,iblk)                            ! N/m^2
+             workx = -strocnxT_iavg(i,j,iblk)                       ! N/m^2
+             worky = -strocnyT_iavg(i,j,iblk)                       ! N/m^2
              tauxo(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) - worky*sin(ANGLET(i,j,iblk))
              tauyo(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) + workx*sin(ANGLET(i,j,iblk))
           enddo

diff --git a/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90 b/cicecore/drivers/nuopc/cmeps/ice_prescribed_mod.F90
@@ -472,11 +472,11 @@ subroutine ice_prescribed_phys()
     ! set non-computed fluxes, ice velocities, ice-ocn stresses to zero
     !--------------------------------------------------------------------
 
-    frzmlt    (:,:,:) = c0
-    uvel      (:,:,:) = c0
-    vvel      (:,:,:) = c0
-    strocnxT  (:,:,:) = c0
-    strocnyT  (:,:,:) = c0
+    frzmlt       (:,:,:) = c0
+    uvel         (:,:,:) = c0
+    vvel         (:,:,:) = c0
+    strocnxT_iavg(:,:,:) = c0
+    strocnyT_iavg(:,:,:) = c0
 
     !-----------------------------------------------------------------
     ! other atm and ocn fluxes

diff --git a/cicecore/drivers/nuopc/dmi/cice_cap.info b/cicecore/drivers/nuopc/dmi/cice_cap.info
@@ -1023,8 +1023,8 @@ module cice_cap
           dataPtr_vice    (i1,j1,iblk) = vice(i,j,iblk)   ! sea ice volume
           dataPtr_vsno    (i1,j1,iblk) = vsno(i,j,iblk)   ! snow volume
           dataPtr_fswthru (i1,j1,iblk) = fswthru(i,j,iblk) ! short wave penetration through ice
-          ui = strocnxT(i,j,iblk)
-          vj = strocnyT(i,j,iblk)
+          ui = strocnxT_iavg(i,j,iblk)
+          vj = strocnyT_iavg(i,j,iblk)
           angT = ANGLET(i,j,iblk)
           dataPtr_strocnxT(i1,j1,iblk) =  ui*cos(-angT) + vj*sin(angT)  ! ice ocean stress
           dataPtr_strocnyT(i1,j1,iblk) = -ui*sin(angT)  + vj*cos(-angT)  ! ice ocean stress