Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Subroutines for CD-grid rheology #10

Merged
merged 4 commits into from
Nov 16, 2021
Merged

Subroutines for CD-grid rheology #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e4e8b83
In process of coding stress_T subroutine
JFLemieux73 Nov 16, 2021
e5e7d4b
Almost done with stress_T subroutine...it compiles.
JFLemieux73 Nov 16, 2021
7506ce5
Done with stress_T...it compiles
JFLemieux73 Nov 16, 2021
0f75f06
Minor modif to deformations_T subroutine
JFLemieux73 Nov 16, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 144 additions & 0 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Open in desktop
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
Open in desktop
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 )
Comment on lines +1291 to +1300
Copy link
Collaborator

@phil-blain phil-blain Nov 16, 2021 •
edited
Loading

Choose a reason for hiding this comment

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

This whole subroutine looks like it's the same as deformations, only calling strain_ratesT instead of strain_rates.

Should we instead keep the existing deformations subroutine, and add an if (grid_system == 'B') .. else .. block for calling strain_rates vs strain_ratesT ?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Good eyes, @phil-blain .
I would prefer to reuse as much code as possible. Initially it's fine to have duplicative code, while developing, but if it makes sense and is easy to merge/clean up enroute, let's do it...

Copy link
Collaborator

Choose a reason for hiding this comment

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

Yeah, upon further inspection with @JFLemieux73, it's not exactly the same, as the B grid version does the average of the 4 corners:
https://github.com/CICE-Consortium/CICE/blob/0f75f068d63725062da4ea944443fbdd6dbed2fb/cicecore/cicedynB/dynamics/ice_dyn_shared.F90#L1269-L1278

whereas the CD-grid version already has the quantities at the tracer point, so it just computes the deformations without averaging:
https://github.com/CICE-Consortium/CICE/blob/0f75f068d63725062da4ea944443fbdd6dbed2fb/cicecore/cicedynB/dynamics/ice_dyn_shared.F90#L1363-L1370


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