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