From bc38c073a0c058e8e79baf1e4df545464d570f45 Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Fri, 19 Nov 2021 11:35:39 -0500 Subject: [PATCH] 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 | 50 +++++++++++++------ cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 17 ++++--- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 35 ++++++++----- cicecore/cicedynB/general/ice_flux.F90 | 11 +++- .../infrastructure/ice_restart_driver.F90 | 37 ++++++++++++-- .../io/io_netcdf/ice_restart.F90 | 10 ++++ .../infrastructure/io/io_pio2/ice_restart.F90 | 10 ++++ 7 files changed, 132 insertions(+), 38 deletions(-) diff --git a/cicecore/cicedynB/analysis/ice_history.F90 b/cicecore/cicedynB/analysis/ice_history.F90 index 4ade53cfb..e0855aa1c 100644 --- a/cicecore/cicedynB/analysis/ice_history.F90 +++ b/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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 779e84452..a78cdd457 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 7ceb715af..fb0a65d68 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/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,13 +1352,13 @@ 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) @@ -1355,9 +1366,9 @@ subroutine principal_stress(nx_block, ny_block, & 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 diff --git a/cicecore/cicedynB/general/ice_flux.F90 b/cicecore/cicedynB/general/ice_flux.F90 index c593d91b1..00d9aac97 100644 --- a/cicecore/cicedynB/general/ice_flux.F90 +++ b/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') diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 index 1a5681b38..9a5a75bea 100644 --- a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 +++ b/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 diff --git a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 b/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 index f6002ff40..0bba6e36e 100644 --- a/cicecore/cicedynB/infrastructure/io/io_netcdf/ice_restart.F90 +++ b/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 diff --git a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 b/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 index 0ec6b7628..2e5338fc0 100644 --- a/cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90 +++ b/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