Skip to content

Commit

Permalink
Subroutines for CD-grid rheology (CICE-Consortium#10)
Browse files Browse the repository at this point in the history 
* In process of coding stress_T subroutine

* Almost done with stress_T subroutine...it compiles.

* Done with stress_T...it compiles

* Minor modif to deformations_T subroutine
  • Loading branch information
@JFLemieux73
JFLemieux73 authored Nov 16, 2021
1 parent 34a3ada commit e9aa132
Show file tree
Hide file tree
Showing 2 changed files with 241 additions and 3 deletions.
144 changes: 144 additions & 0 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -898,6 +898,150 @@ 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_T, deformations_T, &
viscous_coeffs_and_rep_pressure_T

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
!-----------------------------------------------------------------

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

enddo ! ij

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then

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

end subroutine stress_T

!=======================================================================

end module ice_dyn_evp
Expand Down
100 changes: 97 additions & 3 deletions cicecore/cicedynB/dynamics/ice_dyn_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,11 @@ 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, deformations_T, &
strain_rates, strain_rates_T, &
viscous_coeffs_and_rep_pressure, &
viscous_coeffs_and_rep_pressure_T, &
stack_velocity_field, unstack_velocity_field

! namelist parameters
Expand Down Expand Up @@ -1276,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: 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
Expand Down Expand Up @@ -1506,7 +1601,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
Expand Down

0 comments on commit e9aa132

Please sign in to comment.