From 6ff1a713a55889713077b5125245756f49e74b85 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Sun, 10 Apr 2022 00:49:31 +0000 Subject: [PATCH] New deformationsC_T subroutine for more consistent calculation (#89) * Added new subroutine to calc deformations at T using shearTsqr * Renamed deformations_T suroutines and clean up --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 30 ++-- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 135 ++++++++++++++++-- 2 files changed, 140 insertions(+), 25 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index a3e916619..f18e60802 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -113,7 +113,8 @@ subroutine evp (dt) use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, & ice_dyn_evp_1d_copyout use ice_dyn_shared, only: evp_algorithm, stack_fields, unstack_fields, & - DminTarea, visc_method, deformations, deformations_T, strain_rates_U, & + DminTarea, visc_method, deformations, deformationsC_T, deformationsCD_T, & + strain_rates_U, & dyn_haloUpdate real (kind=dbl_kind), intent(in) :: & @@ -853,16 +854,19 @@ subroutine evp (dt) ! on last subcycle, save quantities for mechanical redistribution !----------------------------------------------------------------- if (ksub == ndte) then - call deformations_T (nx_block , ny_block , & + + call deformationsC_T (nx_block , ny_block , & icellt (iblk), & indxti (:,iblk), indxtj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & dxN (:,:,iblk), dyE (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & - tarear (:,:,iblk), & + tarear (:,:,iblk), uarea (:,:,iblk), & + shearU (:,:,iblk), & shear (:,:,iblk), divu (:,:,iblk), & rdg_conv(:,:,iblk), rdg_shear(:,:,iblk)) + endif enddo !$OMP END PARALLEL DO @@ -996,16 +1000,16 @@ subroutine evp (dt) ! on last subcycle, save quantities for mechanical redistribution !----------------------------------------------------------------- if (ksub == ndte) then - call deformations_T (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - uvelE (:,:,iblk), vvelE (:,:,iblk), & - uvelN (:,:,iblk), vvelN (:,:,iblk), & - dxN (:,:,iblk), dyE (:,:,iblk), & - dxT (:,:,iblk), dyT (:,:,iblk), & - tarear (:,:,iblk), & - shear (:,:,iblk), divu (:,:,iblk), & - rdg_conv(:,:,iblk), rdg_shear(:,:,iblk)) + call deformationsCD_T (nx_block , ny_block , & + icellt (iblk), & + indxti (:,iblk), indxtj (:,iblk), & + uvelE (:,:,iblk), vvelE (:,:,iblk), & + uvelN (:,:,iblk), vvelN (:,:,iblk), & + dxN (:,:,iblk), dyE (:,:,iblk), & + dxT (:,:,iblk), dyT (:,:,iblk), & + tarear (:,:,iblk), & + shear (:,:,iblk), divu (:,:,iblk), & + rdg_conv(:,:,iblk), rdg_shear(:,:,iblk)) endif enddo !$OMP END PARALLEL DO diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index ee33eb09e..bc7f3abb1 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -27,7 +27,7 @@ module ice_dyn_shared principal_stress, init_dyn, dyn_prep1, dyn_prep2, dyn_finish, & seabed_stress_factor_LKD, seabed_stress_factor_prob, & alloc_dyn_shared, & - deformations, deformations_T, & + deformations, deformationsC_T, deformationsCD_T, & strain_rates, strain_rates_T, strain_rates_U, & visc_replpress, & dyn_haloUpdate, & @@ -1733,16 +1733,16 @@ end subroutine deformations ! author: JF Lemieux, ECCC ! Nov 2021 - subroutine deformations_T (nx_block, ny_block, & - icellt, & - indxti, indxtj, & - uvelE, vvelE, & - uvelN, vvelN, & - dxN, dyE, & - dxT, dyT, & - tarear, & - shear, divu, & - rdg_conv, rdg_shear ) + subroutine deformationsCD_T (nx_block, ny_block, & + icellt, & + indxti, indxtj, & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + tarear, & + shear, divu, & + rdg_conv, rdg_shear ) use ice_constants, only: p5 @@ -1820,7 +1820,118 @@ subroutine deformations_T (nx_block, ny_block, & enddo ! ij - end subroutine deformations_T + end subroutine deformationsCD_T + + +!======================================================================= +! Compute deformations for mechanical redistribution at T point +! +! author: JF Lemieux, ECCC +! Nov 2021 + + subroutine deformationsC_T (nx_block, ny_block, & + icellt, & + indxti, indxtj, & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + tarear, uarea, & + shearU, & + shear, divu, & + rdg_conv, rdg_shear ) + + use ice_constants, only: p5 + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + icellt ! no. of cells where icetmask = 1 + + integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & + indxti , & ! compressed index in i-direction + indxtj ! compressed index in j-direction + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & + uvelE , & ! x-component of velocity (m/s) at the E point + vvelE , & ! y-component of velocity (m/s) at the N point + uvelN , & ! x-component of velocity (m/s) at the E point + vvelN , & ! y-component of velocity (m/s) at the N point + dxN , & ! width of N-cell through the middle (m) + dyE , & ! height of E-cell through the middle (m) + dxT , & ! width of T-cell through the middle (m) + dyT , & ! height of T-cell through the middle (m) + tarear , & ! 1/tarea + uarea , & ! area of u cell + shearU ! shearU + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & + shear , & ! strain rate II component (1/s) + divu , & ! strain rate I component, velocity divergence (1/s) + rdg_conv , & ! convergence term for ridging (1/s) + rdg_shear ! shear term for ridging (1/s) + + ! local variables + + integer (kind=int_kind) :: & + i, j, ij + + real (kind=dbl_kind), dimension (nx_block,ny_block) :: & + divT , & ! divergence at T point + tensionT , & ! tension at T point + shearT , & ! shear at T point + DeltaT ! delt at T point + + real (kind=dbl_kind) :: & + tmp , & ! useful combination + shearTsqr ! strain rates squared at T point + + character(len=*), parameter :: subname = '(deformations_T2)' + + !----------------------------------------------------------------- + ! strain rates + ! NOTE these are actually strain rates * area (m^2/s) + !----------------------------------------------------------------- + + call strain_rates_T (nx_block , ny_block , & + icellt , & + indxti(:) , indxtj (:) , & + uvelE (:,:), vvelE (:,:), & + uvelN (:,:), vvelN (:,:), & + dxN (:,:), dyE (:,:), & + dxT (:,:), dyT (:,:), & + divT (:,:), tensionT(:,:), & + shearT(:,:), DeltaT (:,:) ) + + ! DeltaT is calc by strain_rates_T but replaced by calculation below. + + do ij = 1, icellt + i = indxti(ij) + j = indxtj(ij) + + !----------------------------------------------------------------- + ! deformations for mechanical redistribution + !----------------------------------------------------------------- + + shearTsqr = (shearU(i ,j )**2 * uarea(i ,j ) & + + shearU(i ,j-1)**2 * uarea(i ,j-1) & + + shearU(i-1,j-1)**2 * uarea(i-1,j-1) & + + shearU(i-1,j )**2 * uarea(i-1,j )) & + / (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j)) + + DeltaT(i,j) = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr)) + + divu(i,j) = divT(i,j) * tarear(i,j) + tmp = DeltaT(i,j) * tarear(i,j) + rdg_conv(i,j) = -min(divu(i,j),c0) + rdg_shear(i,j) = p5*(tmp-abs(divu(i,j))) + + ! diagnostic only...maybe we dont want to use shearTsqr here???? + ! shear = sqrt(tension**2 + shearing**2) + shear(i,j) = tarear(i,j)*sqrt( tensionT(i,j)**2 + shearT(i,j)**2 ) + + enddo ! ij + + end subroutine deformationsC_T !======================================================================= ! Compute strain rates