Skip to content

Commit

Permalink
New deformationsC_T subroutine for more consistent calculation (CICE-…
Browse files Browse the repository at this point in the history 
…Consortium#89)

* Added new subroutine to calc deformations at T using shearTsqr

* Renamed deformations_T suroutines and clean up
  • Loading branch information
@JFLemieux73
JFLemieux73 committed Apr 10, 2022
1 parent 014852a commit 6ff1a71
Show file tree
Hide file tree
Showing 2 changed files with 140 additions and 25 deletions.
30 changes: 17 additions & 13 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) :: &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
135 changes: 123 additions & 12 deletions cicecore/cicedynB/dynamics/ice_dyn_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 6ff1a71

Please sign in to comment.