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.

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

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

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

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

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