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.

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

Eap #19

Merged
merged 2 commits into from
Oct 3, 2021
Merged

Eap #19

Changes from all commits
Commits
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
94 changes: 38 additions & 56 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1939,58 +1939,38 @@ subroutine stepa (nx_block, ny_block, &
j = indxtj(ij)

! ne
call calc_ffrac(1, stressp_1(i,j), stressm_1(i,j), &
stress12_1(i,j), &
a11_1(i,j), &
mresult11)

call calc_ffrac(2, stressp_1(i,j), stressm_1(i,j), &
stress12_1(i,j), &
a12_1(i,j), &
mresult12)
call calc_ffrac(stressp_1(i,j), stressm_1(i,j), &
stress12_1(i,j), &
a11_1(i,j), a12_1(i,j), &
mresult11, mresult12)

a11_1(i,j) = (a11_1(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
a12_1(i,j) = (a12_1(i,j)*dtei - mresult12) * dteikth ! implicit


! nw
call calc_ffrac(1, stressp_2(i,j), stressm_2(i,j), &
stress12_2(i,j), &
a11_2(i,j), &
mresult11)

call calc_ffrac(2, stressp_2(i,j), stressm_2(i,j), &
stress12_2(i,j), &
a12_2(i,j), &
mresult12)
call calc_ffrac(stressp_2(i,j), stressm_2(i,j), &
stress12_2(i,j), &
a11_2(i,j), a12_2(i,j), &
mresult11, mresult12)

a11_2(i,j) = (a11_2(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
a12_2(i,j) = (a12_2(i,j)*dtei - mresult12) * dteikth ! implicit

! sw
call calc_ffrac(1, stressp_3(i,j), stressm_3(i,j), &
stress12_3(i,j), &
a11_3(i,j), &
mresult11)

call calc_ffrac(2, stressp_3(i,j), stressm_3(i,j), &
stress12_3(i,j), &
a12_3(i,j), &
mresult12)
call calc_ffrac(stressp_3(i,j), stressm_3(i,j), &
stress12_3(i,j), &
a11_3(i,j), a12_3(i,j), &
mresult11, mresult12)

a11_3(i,j) = (a11_3(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
a12_3(i,j) = (a12_3(i,j)*dtei - mresult12) * dteikth ! implicit

! se
call calc_ffrac(1, stressp_4(i,j), stressm_4(i,j), &
stress12_4(i,j), &
a11_4(i,j), &
mresult11)

call calc_ffrac(2, stressp_4(i,j), stressm_4(i,j), &
stress12_4(i,j), &
a12_4(i,j), &
mresult12)
call calc_ffrac(stressp_4(i,j), stressm_4(i,j), &
stress12_4(i,j), &
a11_4(i,j), a12_4(i,j), &
mresult11, mresult12)

a11_4(i,j) = (a11_4(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
a12_4(i,j) = (a12_4(i,j)*dtei - mresult12) * dteikth ! implicit
Expand All @@ -2010,19 +1990,17 @@ end subroutine stepa
! the ice floe re-orientation due to fracture
! Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2

subroutine calc_ffrac (blockno, stressp, stressm, &
stress12, &
a1x, &
mresult)

integer(kind=int_kind), intent(in) :: &
blockno
subroutine calc_ffrac (stressp, stressm, &
stress12, &
a1x, a2x, &
mresult1, mresult2)

real (kind=dbl_kind), intent(in) :: &
stressp, stressm, stress12, a1x
stressp, stressm, stress12, a1x, a2x

real (kind=dbl_kind), intent(out) :: &
mresult
mresult1, mresult2

! local variables

Expand All @@ -2042,11 +2020,12 @@ subroutine calc_ffrac (blockno, stressp, stressm, &
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

sigma11 = p5*(stressp+stressm)
sigma12 = stress12
sigma22 = p5*(stressp-stressm)
sigma11 = p5*(stressp+stressm)
sigma12 = stress12
sigma22 = p5*(stressp-stressm)

if ((sigma11-sigma22) == c0) then
! if ((sigma11-sigma22) == c0) then sigma11-sigma22 == 0 => stressn ==0
if (stressm == c0) then
gamma = p5*(pih)
else
gamma = p5*atan2((c2*sigma12),(sigma11-sigma22))
Expand All @@ -2068,24 +2047,27 @@ subroutine calc_ffrac (blockno, stressp, stressm, &

! Pure divergence
if ((sigma_1 >= c0).and.(sigma_2 >= c0)) then
mresult = c0
mresult1 = c0
mresult2 = c0

! Unconfined compression: cracking of blocks not along the axial splitting direction
! which leads to the loss of their shape, so we again model it through diffusion
elseif ((sigma_1 >= c0).and.(sigma_2 < c0)) then
if (blockno == 1) mresult = kfrac * (a1x - Q12Q12)
if (blockno == 2) mresult = kfrac * (a1x + Q11Q12)
mresult1 = kfrac * (a1x - Q12Q12)
mresult2 = kfrac * (a2x + Q11Q12)

! Shear faulting
elseif (sigma_2 == c0) then
mresult = c0
mresult1 = c0
mresult2 = c0
elseif ((sigma_1 <= c0).and.(sigma_1/sigma_2 <= threshold)) then
if (blockno == 1) mresult = kfrac * (a1x - Q12Q12)
if (blockno == 2) mresult = kfrac * (a1x + Q11Q12)
mresult1 = kfrac * (a1x - Q12Q12)
mresult2 = kfrac * (a2x + Q11Q12)

! Horizontal spalling
else
mresult = c0
else
mresult1 = c0
mresult2 = c0
endif

end subroutine calc_ffrac
Expand Down