Add CD-grid stress arrays (#28)

* ice_flux: add stress arrays at CD-grid locations (T,U)

* ice_dyn_shared: initialize CD-grid stress arrays

* infrastructure: add CD-grid stress arrays to restarts

* ice_dyn_shared: generalize 'principal_stress' arguments names

In a subsequent commit we will call 'principal_stress' with the CD-grid
stress arrays 'stress{p,m,12}T' to compute the principal stresses at the
tracer point when using the CD grid.

In that light, remove the '_1' suffix from the stress arguments since
they won't always be located at the NE corner anymore.

* ice_history: compute 'sig[12P]' at T-point for CD-grid

For the CD-grid, compute the principal stresses and the ice pressure at
the T point by passing the appropriate arrays to 'principal_stress'.

These three history variables are computed using the NE-corner values on
the B-grid, but this is not reflected in the description of the fields
in the history output.

Add the location at which the stresses are computed to the comment
argument in the call to 'define_hist_field', for both the B and CD grid.

* ice_dyn_evp: pass 'stress{p,m,12}[TU]' to dyn_prep2

For the CD-grid, the 'stress{p,m,12}_[1-4]' arrays are not used. Pass
the CD-grid location stress arrays 'stress{p,m,12}[TU]' to 'dyn_prep2',
which zero-initializes them anywhere icetmask is zero.

Also, add a TODO in the `if (grid_type) == 'tripole'` block, since it is
not yet clear how the halo updates should be done for the new stress
arrays on the tripole grid.

cicecore/cicedynB/analysis/ice_history.F90

@@ -94,6 +94,7 @@ subroutine init_hist (dt)
       integer (kind=int_kind), dimension(max_nstrm) :: &
          ntmp
       integer (kind=int_kind) :: nml_error ! namelist i/o error flag
+      character(len=char_len) :: description
       character(len=*), parameter :: subname = '(init_hist)'
 
       !-----------------------------------------------------------------
@@ -1199,19 +1200,26 @@ subroutine init_hist (dt)
              "none", secday*c100, c0,                                &
              ns1, f_shear)
 
+         select case (grid_system)
+         case('B')
+            description = ", on U grid  (NE corner values)"
+         case ('CD')
+            description = ", on T grid"
+         end select
+
          call define_hist_field(n_sig1,"sig1","1",ustr2D, ucstr, &
              "norm. principal stress 1",                       &
-             "sig1 is instantaneous", c1, c0,                  &
+             "sig1 is instantaneous" // trim(description), c1, c0, &
              ns1, f_sig1)
 
          call define_hist_field(n_sig2,"sig2","1",ustr2D, ucstr, &
              "norm. principal stress 2",                       &
-             "sig2 is instantaneous", c1, c0,                  &
+             "sig2 is instantaneous" // trim(description), c1, c0, &
              ns1, f_sig2)
 
          call define_hist_field(n_sigP,"sigP","1",ustr2D, ucstr, &
              "ice pressure",                       &
-             "sigP is instantaneous", c1, c0,                  &
+             "sigP is instantaneous" // trim(description), c1, c0, &
              ns1, f_sigP)
 
          call define_hist_field(n_dvidtt,"dvidtt","cm/day",tstr2D, tcstr, &
@@ -1979,7 +1987,7 @@ subroutine accum_hist (dt)
       use ice_blocks, only: block, get_block, nx_block, ny_block
       use ice_domain, only: blocks_ice, nblocks
       use ice_domain_size, only: nfsd
-      use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu
+      use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu, grid_system
       use ice_calendar, only: new_year, write_history, &
                               write_ic, timesecs, histfreq, nstreams, mmonth, &
                               new_month
@@ -2001,6 +2009,7 @@ subroutine accum_hist (dt)
           fm, fmN, fmE, daidtt, dvidtt, daidtd, dvidtd, fsurf, &
           fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, snwcnt, &
           stressp_1, stressm_1, stress12_1, &
+          stresspT, stressmT, stress12T, &
           stressp_2, &
           stressp_3, &
           stressp_4, sig1, sig2, sigP, &
@@ -4287,17 +4296,28 @@ subroutine accum_hist (dt)
       ! snapshots
       !---------------------------------------------------------------
 
-          ! compute sig1 and sig2
-
-           call principal_stress (nx_block,  ny_block,  &
-                                  stressp_1 (:,:,iblk), &
-                                  stressm_1 (:,:,iblk), &
-                                  stress12_1(:,:,iblk), &
-                                  strength  (:,:,iblk), &
-                                  sig1      (:,:,iblk), &
-                                  sig2      (:,:,iblk), &
-                                  sigP      (:,:,iblk))
-
+            ! compute sig1 and sig2
+            select case (grid_system)
+            case('B')
+               call principal_stress (nx_block,  ny_block,  &
+                                      stressp_1 (:,:,iblk), &
+                                      stressm_1 (:,:,iblk), &
+                                      stress12_1(:,:,iblk), &
+                                      strength  (:,:,iblk), &
+                                      sig1      (:,:,iblk), &
+                                      sig2      (:,:,iblk), &
+                                      sigP      (:,:,iblk))
+            case('CD')
+               call principal_stress (nx_block,  ny_block,  &
+                                      stresspT  (:,:,iblk), &
+                                      stressmT  (:,:,iblk), &
+                                      stress12T (:,:,iblk), &
+                                      strength  (:,:,iblk), &
+                                      sig1      (:,:,iblk), &
+                                      sig2      (:,:,iblk), &
+                                      sigP      (:,:,iblk))
+            end select
+
            do j = jlo, jhi
            do i = ilo, ihi
               if (.not. tmask(i,j,iblk)) then ! mask out land points

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -95,7 +95,9 @@ subroutine evp (dt)
           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
+          stress12_1, stress12_2, stress12_3, stress12_4, &
+          stresspT, stressmT, stress12T, &
+          stresspU, stressmU, stress12U
       use ice_grid, only: tmask, umask, nmask, emask, dxt, dyt, &
           dxhy, dyhx, cxp, cyp, cxm, cym, &
           tarear, uarear, tinyarea, grid_average_X2Y, &
@@ -410,11 +412,11 @@ subroutine evp (dt)
                          taubxN    (:,:,iblk), taubyN    (:,:,iblk), & 
                          waterxN   (:,:,iblk), wateryN   (:,:,iblk), & 
                          forcexN   (:,:,iblk), forceyN   (:,:,iblk), & 
-                         stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & 
+                         stresspT  (:,:,iblk), stressp_2 (:,:,iblk), & 
                          stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & 
-                         stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & 
+                         stressmT  (:,:,iblk), stressm_2 (:,:,iblk), & 
                          stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & 
-                         stress12_1(:,:,iblk), stress12_2(:,:,iblk), & 
+                         stress12T (:,:,iblk), stress12_2(:,:,iblk), & 
                          stress12_3(:,:,iblk), stress12_4(:,:,iblk), & 
                          uvelN_init (:,:,iblk), vvelN_init (:,:,iblk), &
                          uvelN      (:,:,iblk), vvelN      (:,:,iblk), &
@@ -443,11 +445,11 @@ subroutine evp (dt)
                          taubxE    (:,:,iblk), taubyE    (:,:,iblk), & 
                          waterxE   (:,:,iblk), wateryE   (:,:,iblk), & 
                          forcexE   (:,:,iblk), forceyE   (:,:,iblk), & 
-                         stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & 
+                         stresspU  (:,:,iblk), stressp_2 (:,:,iblk), & 
                          stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & 
-                         stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & 
+                         stressmU  (:,:,iblk), stressm_2 (:,:,iblk), & 
                          stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & 
-                         stress12_1(:,:,iblk), stress12_2(:,:,iblk), & 
+                         stress12U (:,:,iblk), stress12_2(:,:,iblk), & 
                          stress12_3(:,:,iblk), stress12_4(:,:,iblk), & 
                          uvelE_init (:,:,iblk), vvelE_init (:,:,iblk), &
                          uvelE      (:,:,iblk), vvelE      (:,:,iblk), &
@@ -726,6 +728,7 @@ subroutine evp (dt)
 
       ! Force symmetry across the tripole seam
       if (trim(grid_type) == 'tripole') then
+         ! TODO: CD-grid
       if (maskhalo_dyn) then
          !-------------------------------------------------------
          ! set halomask to zero because ice_HaloMask always keeps

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -158,7 +158,9 @@ subroutine init_dyn (dt)
       use ice_flux, only: rdg_conv, rdg_shear, iceumask, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
-          stress12_1, stress12_2, stress12_3, stress12_4
+          stress12_1, stress12_2, stress12_3, stress12_4, &
+          stresspT, stressmT, stress12T, &
+          stresspU, stressmU, stress12U
       use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN, divu, shear
       use ice_grid, only: ULAT, NLAT, ELAT
 
@@ -247,6 +249,15 @@ subroutine init_dyn (dt)
          stress12_3(i,j,iblk) = c0
          stress12_4(i,j,iblk) = c0
 
+         if (grid_system == 'CD') then
+            stresspT  (i,j,iblk) = c0
+            stressmT  (i,j,iblk) = c0
+            stress12T (i,j,iblk) = c0
+            stresspU  (i,j,iblk) = c0
+            stressmU  (i,j,iblk) = c0
+            stress12U (i,j,iblk) = c0
+         endif
+
          ! ice extent mask on velocity points
          iceumask(i,j,iblk) = .false.
 
@@ -1341,23 +1352,23 @@ end subroutine seabed_stress_factor_prob
 !=======================================================================
 
 ! Computes principal stresses for comparison with the theoretical
-! yield curve; northeast values
+! yield curve
 !
 ! author: Elizabeth C. Hunke, LANL
 
       subroutine principal_stress(nx_block,   ny_block,  &
-                                  stressp_1,  stressm_1, &
-                                  stress12_1, strength,  &
+                                  stressp,    stressm,   &
+                                  stress12,   strength,  &
                                   sig1,       sig2,      &
                                   sigP)
 
       integer (kind=int_kind), intent(in) :: &
          nx_block, ny_block  ! block dimensions
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
-         stressp_1 , & ! sigma11 + sigma22
-         stressm_1 , & ! sigma11 - sigma22
-         stress12_1, & ! sigma12
+         stressp   , & ! sigma11 + sigma22
+         stressm   , & ! sigma11 - sigma22
+         stress12  , & ! sigma12
          strength      ! for normalization of sig1 and sig2
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
@@ -1382,14 +1393,14 @@ subroutine principal_stress(nx_block,   ny_block,  &
       do i = 1, nx_block
          if (strength(i,j) > puny) then
             ! ice internal pressure          
-            sigP(i,j) = -p5*stressp_1(i,j) 
+            sigP(i,j) = -p5*stressp(i,j) 
 
             ! normalized principal stresses
-            sig1(i,j) = (p5*(stressp_1(i,j) &
-                      + sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
+            sig1(i,j) = (p5*(stressp(i,j) &
+                      + sqrt(stressm(i,j)**2+c4*stress12(i,j)**2))) &
                       / strength(i,j)
-            sig2(i,j) = (p5*(stressp_1(i,j) &
-                      - sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
+            sig2(i,j) = (p5*(stressp(i,j) &
+                      - sqrt(stressm(i,j)**2+c4*stress12(i,j)**2))) &
                       / strength(i,j)         
          else
             sig1(i,j) = spval_dbl

cicecore/cicedynB/general/ice_flux.F90

@@ -134,7 +134,10 @@ module ice_flux
        ! ice stress tensor in each corner of T cell (kg/s^2)
          stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22
          stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22
-         stress12_1,stress12_2,stress12_3,stress12_4    ! sigma12
+         stress12_1,stress12_2,stress12_3,stress12_4, & ! sigma12
+       ! ice stress tensor at U and T locations (grid_system = 'CD') (kg/s^2)
+         stresspT, stressmT, stress12T, & ! sigma11+sigma22, sigma11-sigma22, sigma12
+         stresspU, stressmU, stress12U    ! "
 
       logical (kind=log_kind), &
          dimension (:,:,:), allocatable, public :: &
@@ -623,6 +626,12 @@ subroutine alloc_flux
          iceemask   (nx_block,ny_block,max_blocks), & ! ice extent mask (E-cell)
          fmE        (nx_block,ny_block,max_blocks), & ! Coriolis param. * mass in E-cell (kg/s)
          TbE        (nx_block,ny_block,max_blocks), & ! factor for seabed stress (landfast ice)
+         stresspT   (nx_block,ny_block,max_blocks), & ! sigma11+sigma22
+         stressmT   (nx_block,ny_block,max_blocks), & ! sigma11-sigma22
+         stress12T  (nx_block,ny_block,max_blocks), & ! sigma12
+         stresspU   (nx_block,ny_block,max_blocks), & ! sigma11+sigma22
+         stressmU   (nx_block,ny_block,max_blocks), & ! sigma11-sigma22
+         stress12U  (nx_block,ny_block,max_blocks), & ! sigma12
          stat=ierr)
       if (ierr/=0) call abort_ice('(alloc_flux): Out of memory')
 

cicecore/cicedynB/infrastructure/ice_restart_driver.F90

@@ -57,8 +57,11 @@ subroutine dumpfile(filename_spec)
           strocnxT, strocnyT, sst, frzmlt, iceumask, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
-          stress12_1, stress12_2, stress12_3, stress12_4
+          stress12_1, stress12_2, stress12_3, stress12_4, &
+          stresspT, stressmT, stress12T, &
+          stresspU, stressmU, stress12U
       use ice_flux, only: coszen
+      use ice_grid, only: grid_system
       use ice_state, only: aicen, vicen, vsnon, trcrn, uvel, vvel
 
       character(len=char_len_long), intent(in), optional :: filename_spec
@@ -164,6 +167,15 @@ subroutine dumpfile(filename_spec)
       call write_restart_field(nu_dump,0,stress12_2,'ruf8','stress12_2',1,diag)
       call write_restart_field(nu_dump,0,stress12_4,'ruf8','stress12_4',1,diag)
 
+      if (grid_system == 'CD') then
+         call write_restart_field(nu_dump,0,stresspT ,'ruf8','stresspT' ,1,diag)
+         call write_restart_field(nu_dump,0,stressmT ,'ruf8','stressmT' ,1,diag)
+         call write_restart_field(nu_dump,0,stress12T,'ruf8','stress12T',1,diag)
+         call write_restart_field(nu_dump,0,stresspU ,'ruf8','stresspU' ,1,diag)
+         call write_restart_field(nu_dump,0,stressmU ,'ruf8','stressmU' ,1,diag)
+         call write_restart_field(nu_dump,0,stress12U,'ruf8','stress12U',1,diag)
+      endif
+
       !-----------------------------------------------------------------
       ! ice mask for dynamics
       !-----------------------------------------------------------------
@@ -206,9 +218,11 @@ subroutine restartfile (ice_ic)
           strocnxT, strocnyT, sst, frzmlt, iceumask, &
           stressp_1, stressp_2, stressp_3, stressp_4, &
           stressm_1, stressm_2, stressm_3, stressm_4, &
-          stress12_1, stress12_2, stress12_3, stress12_4
+          stress12_1, stress12_2, stress12_3, stress12_4, &
+          stresspT, stressmT, stress12T, &
+          stresspU, stressmU, stress12U
       use ice_flux, only: coszen
-      use ice_grid, only: tmask, grid_type
+      use ice_grid, only: tmask, grid_type, grid_system
       use ice_state, only: trcr_depend, aice, vice, vsno, trcr, &
           aice0, aicen, vicen, vsnon, trcrn, aice_init, uvel, vvel, &
           trcr_base, nt_strata, n_trcr_strata
@@ -367,6 +381,21 @@ subroutine restartfile (ice_ic)
       call read_restart_field(nu_restart,0,stress12_4,'ruf8', &
            'stress12_4',1,diag,field_loc_center,field_type_scalar) ! stress12_4
 
+      if (grid_system == 'CD') then
+         call read_restart_field(nu_restart,0,stresspT,'ruf8', &
+            'stresspT' ,1,diag,field_loc_center,field_type_scalar) ! stresspT
+         call read_restart_field(nu_restart,0,stressmT,'ruf8', &
+            'stressmT' ,1,diag,field_loc_center,field_type_scalar) ! stressmT
+         call read_restart_field(nu_restart,0,stress12T,'ruf8', &
+            'stress12T',1,diag,field_loc_center,field_type_scalar) ! stress12T
+         call read_restart_field(nu_restart,0,stresspU,'ruf8', &
+            'stresspU' ,1,diag,field_loc_center,field_type_scalar) ! stresspU
+         call read_restart_field(nu_restart,0,stressmU,'ruf8', &
+            'stressmU' ,1,diag,field_loc_center,field_type_scalar) ! stressmU
+         call read_restart_field(nu_restart,0,stress12U,'ruf8', &
+            'stress12U',1,diag,field_loc_center,field_type_scalar) ! stress12U
+      endif
+
       if (trim(grid_type) == 'tripole') then
          call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info, &
                                     field_loc_center,  field_type_scalar)
@@ -394,6 +423,7 @@ subroutine restartfile (ice_ic)
                                     field_loc_center,  field_type_scalar)
          call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info, &
                                     field_loc_center,  field_type_scalar)
+         ! TODO: CD-grid
       endif
 
       !-----------------------------------------------------------------
@@ -465,6 +495,7 @@ subroutine restartfile (ice_ic)
             stress12_4(i,j,iblk) = c0
          enddo
          enddo
+         ! TODO: CD-grid ?
       enddo
       !$OMP END PARALLEL DO
 

cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90

@@ -137,6 +137,7 @@ subroutine init_restart_write(filename_spec)
                                  n_dic, n_don, n_fed, n_fep, nfsd
       use ice_arrays_column, only: oceanmixed_ice
       use ice_dyn_shared, only: kdyn
+      use ice_grid, only: grid_system
 
       character(len=char_len_long), intent(in), optional :: filename_spec
 
@@ -273,6 +274,15 @@ subroutine init_restart_write(filename_spec)
          call define_rest_field(ncid,'stress12_3',dims)
          call define_rest_field(ncid,'stress12_4',dims)
 
+         if (grid_system == 'CD') then
+            call define_rest_field(ncid,'stresspT' ,dims)
+            call define_rest_field(ncid,'stressmT' ,dims)
+            call define_rest_field(ncid,'stress12T',dims)
+            call define_rest_field(ncid,'stresspU' ,dims)
+            call define_rest_field(ncid,'stressmU' ,dims)
+            call define_rest_field(ncid,'stress12U',dims)
+         endif
+
          call define_rest_field(ncid,'iceumask',dims)
 
          if (oceanmixed_ice) then

cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90

@@ -145,6 +145,7 @@ subroutine init_restart_write(filename_spec)
                                  n_dic, n_don, n_fed, n_fep, nfsd
       use ice_dyn_shared, only: kdyn
       use ice_arrays_column, only: oceanmixed_ice
+      use ice_grid, only: grid_system
 
       logical (kind=log_kind) :: &
           solve_zsal, skl_bgc, z_tracers
@@ -276,6 +277,15 @@ subroutine init_restart_write(filename_spec)
          call define_rest_field(File,'stress12_3',dims)
          call define_rest_field(File,'stress12_4',dims)
 
+         if (grid_system == 'CD') then
+            call define_rest_field(File,'stresspT' ,dims)
+            call define_rest_field(File,'stressmT' ,dims)
+            call define_rest_field(File,'stress12T',dims)
+            call define_rest_field(File,'stresspU' ,dims)
+            call define_rest_field(File,'stressmU' ,dims)
+            call define_rest_field(File,'stress12U',dims)
+         endif
+
          call define_rest_field(File,'iceumask',dims)
 
          if (oceanmixed_ice) then