From e4e8b8336603699e61cd5333fd8ff919af0dadc9 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 16 Nov 2021 16:15:04 +0000 Subject: [PATCH 1/4] In process of coding stress_T subroutine --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 153 +++++++++++++++++++++ 1 file changed, 153 insertions(+) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index afb8eadb7..36ed385d9 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -898,6 +898,159 @@ subroutine stress (nx_block, ny_block, & end subroutine stress +!======================================================================= + +! Computes the strain rates and internal stress components for T points +! Computes stress terms for the momentum equation + +! author: JF Lemieux, ECCC +! Nov 2021 + + subroutine stress_T (nx_block, ny_block, & + ksub, icellt, & + indxti, indxtj, & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + tarear, tinyarea, & + strength, & + stresspT, stressmT, & + stress12T, & + shear, divu, & + rdg_conv, rdg_shear ) + +! use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure + + integer (kind=int_kind), intent(in) :: & + nx_block, ny_block, & ! block dimensions + ksub , & ! subcycling step + 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) + strength , & ! ice strength (N/m) + tarear , & ! 1/tarea + tinyarea ! puny*tarea + + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & + stresspT , & ! sigma11+sigma22 + stressmT , & ! sigma11-sigma22 + stress12T ! sigma12 + + 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) :: & + divT, tensionT, shearT, DeltaT, & ! strain rates at T point + zetax2T , & ! 2 x zeta (visc coeff) at T point + etax2T , & ! 2 x eta (visc coeff) at T point + rep_prsT ! replacement pressure at T point + + logical :: capping ! of the viscous coef + + character(len=*), parameter :: subname = '(stress_T)' + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + + capping = .true. ! could be later included in ice_in + + do ij = 1, icellt + i = indxti(ij) + j = indxtj(ij) + + !----------------------------------------------------------------- + ! strain rates at T point + ! NOTE these are actually strain rates * area (m^2/s) + !----------------------------------------------------------------- + +! call strain_rates_T (nx_block, ny_block, & +! i, j, & +! uvelE, vvelE, & +! uvelN, vvelN, & +! dxN, dyE, & +! dxT, dyT, & +! divT, tensionT, & +! shearT, DeltaT ) + + !----------------------------------------------------------------- + ! viscous coefficients and replacement pressure at T point + !----------------------------------------------------------------- + ! TODO JFL +! call viscous_coeffs_and_rep_pressure_T (strength(i,j), & +! tinyarea(i,j), & +! DeltaT, zetax2T, etax2T, & +! rep_prsT, capping ) + + !----------------------------------------------------------------- + ! the stresses ! kg/s^2 + !----------------------------------------------------------------- + + ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code + + stresspT(i,j) = (stresspT(i,j)*(c1-arlx1i*revp) + & + arlx1i*(zetax2T*divT - rep_prsT)) * denom1 + + stressmT(i,j) = (stressmT(i,j)*(c1-arlx1i*revp) + & + arlx1i*etax2T*tensionT) * denom1 + + stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) + & + arlx1i*p5*etax2T*shearT) * denom1 + + !----------------------------------------------------------------- + ! for dF/dx (u momentum) + !----------------------------------------------------------------- + + + !----------------------------------------------------------------- + ! for dF/dy (v momentum) + !----------------------------------------------------------------- + + + enddo ! ij + + !----------------------------------------------------------------- + ! on last subcycle, save quantities for mechanical redistribution + !----------------------------------------------------------------- + if (ksub == ndte) then + ! TODO JFL +! call deformations (nx_block , ny_block , & +! icellt , & +! indxti , indxtj , & +! uvel , vvel , & +! dxt , dyt , & +! cxp , cyp , & +! cxm , cym , & +! tarear , & +! shear , divu , & +! rdg_conv , rdg_shear ) + + endif + + end subroutine stress_T + !======================================================================= end module ice_dyn_evp From e5e7d4b1f6e170750cb3a01a0ab43d76fd4bc89b Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 16 Nov 2021 16:52:44 +0000 Subject: [PATCH 2/4] Almost done with stress_T subroutine...it compiles. --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 29 +++++++++---------- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 5 ++-- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 36ed385d9..42c7a48e9 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -920,7 +920,7 @@ subroutine stress_T (nx_block, ny_block, & shear, divu, & rdg_conv, rdg_shear ) -! use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure + use ice_dyn_shared, only: strain_rates_T, viscous_coeffs_and_rep_pressure_T integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -986,23 +986,23 @@ subroutine stress_T (nx_block, ny_block, & ! NOTE these are actually strain rates * area (m^2/s) !----------------------------------------------------------------- -! call strain_rates_T (nx_block, ny_block, & -! i, j, & -! uvelE, vvelE, & -! uvelN, vvelN, & -! dxN, dyE, & -! dxT, dyT, & -! divT, tensionT, & -! shearT, DeltaT ) + call strain_rates_T (nx_block, ny_block, & + i, j, & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + divT, tensionT, & + shearT, DeltaT ) !----------------------------------------------------------------- ! viscous coefficients and replacement pressure at T point !----------------------------------------------------------------- - ! TODO JFL -! call viscous_coeffs_and_rep_pressure_T (strength(i,j), & -! tinyarea(i,j), & -! DeltaT, zetax2T, etax2T, & -! rep_prsT, capping ) + + call viscous_coeffs_and_rep_pressure_T (strength(i,j), & + tinyarea(i,j), & + DeltaT, zetax2T, etax2T, & + rep_prsT, capping ) !----------------------------------------------------------------- ! the stresses ! kg/s^2 @@ -1028,7 +1028,6 @@ subroutine stress_T (nx_block, ny_block, & ! for dF/dy (v momentum) !----------------------------------------------------------------- - enddo ! ij !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index fb86f520e..e70e443d7 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -26,8 +26,10 @@ module ice_dyn_shared public :: init_dyn, set_evp_parameters, stepu, principal_stress, & dyn_prep1, dyn_prep2, dyn_finish, & seabed_stress_factor_LKD, seabed_stress_factor_prob, & - alloc_dyn_shared, deformations, strain_rates, & + alloc_dyn_shared, deformations, & + strain_rates, strain_rates_T, & viscous_coeffs_and_rep_pressure, & + viscous_coeffs_and_rep_pressure_T, & stack_velocity_field, unstack_velocity_field ! namelist parameters @@ -1506,7 +1508,6 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & ! endif end subroutine viscous_coeffs_and_rep_pressure - !======================================================================= ! Computes viscous coefficients and replacement pressure for stress From 7506ce5d5314a1e4ec33429d56441537fa67e0c5 Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 16 Nov 2021 17:12:57 +0000 Subject: [PATCH 3/4] Done with stress_T...it compiles --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 34 +++---- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 97 ++++++++++++++++++- 2 files changed, 108 insertions(+), 23 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 42c7a48e9..1e308b65a 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -920,7 +920,8 @@ subroutine stress_T (nx_block, ny_block, & shear, divu, & rdg_conv, rdg_shear ) - use ice_dyn_shared, only: strain_rates_T, viscous_coeffs_and_rep_pressure_T + use ice_dyn_shared, only: strain_rates_T, deformations_T, & + viscous_coeffs_and_rep_pressure_T integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1019,32 +1020,23 @@ subroutine stress_T (nx_block, ny_block, & stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) + & arlx1i*p5*etax2T*shearT) * denom1 - !----------------------------------------------------------------- - ! for dF/dx (u momentum) - !----------------------------------------------------------------- - - - !----------------------------------------------------------------- - ! for dF/dy (v momentum) - !----------------------------------------------------------------- - enddo ! ij !----------------------------------------------------------------- ! on last subcycle, save quantities for mechanical redistribution !----------------------------------------------------------------- if (ksub == ndte) then - ! TODO JFL -! call deformations (nx_block , ny_block , & -! icellt , & -! indxti , indxtj , & -! uvel , vvel , & -! dxt , dyt , & -! cxp , cyp , & -! cxm , cym , & -! tarear , & -! shear , divu , & -! rdg_conv , rdg_shear ) + + call deformations_T (nx_block , ny_block , & + icellt , & + indxti , indxtj , & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + tarear , & + shear , divu , & + rdg_conv , rdg_shear ) endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index e70e443d7..db79c3568 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -26,7 +26,8 @@ module ice_dyn_shared public :: init_dyn, set_evp_parameters, stepu, principal_stress, & dyn_prep1, dyn_prep2, dyn_finish, & seabed_stress_factor_LKD, seabed_stress_factor_prob, & - alloc_dyn_shared, deformations, & + alloc_dyn_shared, & + deformations, deformations_T, & strain_rates, strain_rates_T, & viscous_coeffs_and_rep_pressure, & viscous_coeffs_and_rep_pressure_T, & @@ -1278,8 +1279,100 @@ subroutine deformations (nx_block, ny_block, & enddo ! ij - end subroutine deformations + end subroutine deformations + +!======================================================================= + +! Compute deformations for mechanical redistribution at T point +! +! 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 ) + + use ice_constants, only: p25, 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 + + 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) :: & + divT, tensionT, shearT, DeltaT, & ! strain rates at T point + tmp ! useful combination + + character(len=*), parameter :: subname = '(deformations_T)' + + do ij = 1, icellt + i = indxti(ij) + j = indxtj(ij) + + !----------------------------------------------------------------- + ! strain rates + ! NOTE these are actually strain rates * area (m^2/s) + !----------------------------------------------------------------- + + call strain_rates_T (nx_block, ny_block, & + i, j, & + uvelE, vvelE, & + uvelN, vvelN, & + dxN, dyE, & + dxT, dyT, & + divT, tensionT, & + shearT, DeltaT ) + + !----------------------------------------------------------------- + ! deformations for mechanical redistribution + !----------------------------------------------------------------- + divu(i,j) = divT * tarear(i,j) + tmp = Deltat * tarear(i,j) + rdg_conv(i,j) = -min(divu(i,j),c0) + rdg_shear(i,j) = p5*(tmp-abs(divu(i,j))) + + ! diagnostic only + ! shear = sqrt(tension**2 + shearing**2) + shear(i,j) = tarear(i,j)*sqrt( tensionT**2 + shearT**2 ) + + enddo ! ij + end subroutine deformations_T + !======================================================================= ! Compute strain rates From 0f75f068d63725062da4ea944443fbdd6dbed2fb Mon Sep 17 00:00:00 2001 From: Jean-Francois Lemieux Date: Tue, 16 Nov 2021 19:03:28 +0000 Subject: [PATCH 4/4] Minor modif to deformations_T subroutine --- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index db79c3568..1d27ff792 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -1299,7 +1299,7 @@ subroutine deformations_T (nx_block, ny_block, & shear, divu, & rdg_conv, rdg_shear ) - use ice_constants, only: p25, p5 + use ice_constants, only: p5 integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions