Skip to content

Commit

Permalink
dynamics: refactor strain rates and deformation calculations
Browse files Browse the repository at this point in the history 
Create subroutines for strain rates and deformations computations, add
them to ice_dyn_shared and call them in ice_dyn_evp and ice_dyn_vp.

This reduces code duplication.
  • Loading branch information
@phil-blain
phil-blain committed Jul 13, 2020
1 parent d2871b0 commit 6af142e
Show file tree
Hide file tree
Showing 3 changed files with 297 additions and 194 deletions.
85 changes: 33 additions & 52 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -599,6 +599,8 @@ subroutine stress (nx_block, ny_block, &
rdg_conv, rdg_shear, &
str )

use ice_dyn_shared, only: strain_rates, deformations

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ksub , & ! subcycling step
Expand Down Expand Up @@ -676,58 +678,20 @@ subroutine stress (nx_block, ny_block, &
! strain rates
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------
! divergence = e_11 + e_22
divune = cyp(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxp(i,j)*vvel(i ,j ) - dxt(i,j)*vvel(i ,j-1)
divunw = cym(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxp(i,j)*vvel(i-1,j ) - dxt(i,j)*vvel(i-1,j-1)
divusw = cym(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxm(i,j)*vvel(i-1,j-1) + dxt(i,j)*vvel(i-1,j )
divuse = cyp(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxm(i,j)*vvel(i ,j-1) + dxt(i,j)*vvel(i ,j )

! tension strain rate = e_11 - e_22
tensionne = -cym(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxm(i,j)*vvel(i ,j ) + dxt(i,j)*vvel(i ,j-1)
tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxm(i,j)*vvel(i-1,j ) + dxt(i,j)*vvel(i-1,j-1)
tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxp(i,j)*vvel(i-1,j-1) - dxt(i,j)*vvel(i-1,j )
tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j )

! shearing strain rate = e_12
shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) &
- cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1)
shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) &
- cxm(i,j)*uvel(i-1,j ) - dxt(i,j)*uvel(i-1,j-1)
shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyt(i,j)*vvel(i ,j-1) &
- cxp(i,j)*uvel(i-1,j-1) + dxt(i,j)*uvel(i-1,j )
shearse = -cym(i,j)*vvel(i ,j-1) - dyt(i,j)*vvel(i-1,j-1) &
- cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j )

! Delta (in the denominator of zeta, eta)
Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2))
Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2))
Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2))
Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2))

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
tmp = p25*(Deltane + Deltanw + Deltase + Deltasw) * 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) = p25*tarear(i,j)*sqrt( &
(tensionne + tensionnw + tensionse + tensionsw)**2 &
+ (shearne + shearnw + shearse + shearsw)**2)

endif
call strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
divuse, divusw, &
tensionne, tensionnw, &
tensionse, tensionsw, &
shearne, shearnw, &
shearse, shearsw, &
Deltane, Deltanw, &
Deltase, Deltasw )

!-----------------------------------------------------------------
! strength/Delta ! kg/s
Expand Down Expand Up @@ -902,6 +866,23 @@ subroutine stress (nx_block, ny_block, &

enddo ! ij

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
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

!=======================================================================
Expand Down
197 changes: 196 additions & 1 deletion cicecore/cicedynB/dynamics/ice_dyn_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ module ice_dyn_shared
private
public :: init_evp, set_evp_parameters, stepu, principal_stress, &
dyn_prep1, dyn_prep2, dyn_finish, basal_stress_coeff, &
alloc_dyn_shared
alloc_dyn_shared, deformations, strain_rates

! namelist parameters

Expand Down Expand Up @@ -988,6 +988,201 @@ subroutine principal_stress(nx_block, ny_block, &

end subroutine principal_stress

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

! Compute deformations for mechanical redistribution
!
! author: Elizabeth C. Hunke, LANL
!
! 2019: subroutine created by Philippe Blain, ECCC

subroutine deformations (nx_block, ny_block, &
icellt, &
indxti, indxtj, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
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) :: &
uvel , & ! x-component of velocity (m/s)
vvel , & ! y-component of velocity (m/s)
dxt , & ! width of T-cell through the middle (m)
dyt , & ! height of T-cell through the middle (m)
cyp , & ! 1.5*HTE - 0.5*HTE
cxp , & ! 1.5*HTN - 0.5*HTN
cym , & ! 0.5*HTE - 1.5*HTE
cxm , & ! 0.5*HTN - 1.5*HTN
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) :: & ! at each corner :
divune, divunw, divuse, divusw , & ! divergence
tensionne, tensionnw, tensionse, tensionsw, & ! tension
shearne, shearnw, shearse, shearsw , & ! shearing
Deltane, Deltanw, Deltase, Deltasw , & ! Delta
tmp ! useful combination

character(len=*), parameter :: subname = '(deformations)'

!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu
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 (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
divuse, divusw, &
tensionne, tensionnw, &
tensionse, tensionsw, &
shearne, shearnw, &
shearse, shearsw, &
Deltane, Deltanw, &
Deltase, Deltasw )
!-----------------------------------------------------------------
! deformations for mechanical redistribution
!-----------------------------------------------------------------
divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
tmp = p25*(Deltane + Deltanw + Deltase + Deltasw) * 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) = p25*tarear(i,j)*sqrt( &
(tensionne + tensionnw + tensionse + tensionsw)**2 + &
(shearne + shearnw + shearse + shearsw )**2)

enddo ! ij

end subroutine deformations

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

! Compute strain rates
!
! author: Elizabeth C. Hunke, LANL
!
! 2019: subroutine created by Philippe Blain, ECCC

subroutine strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
divuse, divusw, &
tensionne, tensionnw, &
tensionse, tensionsw, &
shearne, shearnw, &
shearse, shearsw, &
Deltane, Deltanw, &
Deltase, Deltasw )

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block ! block dimensions

integer (kind=int_kind) :: &
i, j ! indices

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uvel , & ! x-component of velocity (m/s)
vvel , & ! y-component of velocity (m/s)
dxt , & ! width of T-cell through the middle (m)
dyt , & ! height of T-cell through the middle (m)
cyp , & ! 1.5*HTE - 0.5*HTE
cxp , & ! 1.5*HTN - 0.5*HTN
cym , & ! 0.5*HTE - 1.5*HTE
cxm ! 0.5*HTN - 1.5*HTN

real (kind=dbl_kind), intent(out):: & ! at each corner :
divune, divunw, divuse, divusw , & ! divergence
tensionne, tensionnw, tensionse, tensionsw, & ! tension
shearne, shearnw, shearse, shearsw , & ! shearing
Deltane, Deltanw, Deltase, Deltasw ! Delta

character(len=*), parameter :: subname = '(strain_rates)'

!-----------------------------------------------------------------
! strain rates
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------

! divergence = e_11 + e_22
divune = cyp(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxp(i,j)*vvel(i ,j ) - dxt(i,j)*vvel(i ,j-1)
divunw = cym(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxp(i,j)*vvel(i-1,j ) - dxt(i,j)*vvel(i-1,j-1)
divusw = cym(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxm(i,j)*vvel(i-1,j-1) + dxt(i,j)*vvel(i-1,j )
divuse = cyp(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxm(i,j)*vvel(i ,j-1) + dxt(i,j)*vvel(i ,j )

! tension strain rate = e_11 - e_22
tensionne = -cym(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxm(i,j)*vvel(i ,j ) + dxt(i,j)*vvel(i ,j-1)
tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxm(i,j)*vvel(i-1,j ) + dxt(i,j)*vvel(i-1,j-1)
tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxp(i,j)*vvel(i-1,j-1) - dxt(i,j)*vvel(i-1,j )
tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j )

! shearing strain rate = e_12
shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) &
- cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1)
shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) &
- cxm(i,j)*uvel(i-1,j ) - dxt(i,j)*uvel(i-1,j-1)
shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyt(i,j)*vvel(i ,j-1) &
- cxp(i,j)*uvel(i-1,j-1) + dxt(i,j)*uvel(i-1,j )
shearse = -cym(i,j)*vvel(i ,j-1) - dyt(i,j)*vvel(i-1,j-1) &
- cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j )

! Delta (in the denominator of zeta, eta)
Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2))
Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2))
Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2))
Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2))

end subroutine strain_rates

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

end module ice_dyn_shared
Expand Down
Loading

1 comment on commit 6af142e

@phil-blain
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@TillRasmussen we talked in the last teleconference of changes in my work for the VP solver that might affect bit4bit-ness of the standard EVP solver. I think this here commit is the main culprit.

Please sign in to comment.