Skip to content

ice_dyn_vp: express rheology in terms of bulk and shear viscosities #647

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

Merged
apcraig merged 11 commits into from
Oct 21, 2021
Merged

ice_dyn_vp: express rheology in terms of bulk and shear viscosities #647

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
a613ccc
Added array for eta and changed name of zetaD...BFB
JFLemieux73 Oct 14, 2021
d0648e9
Added array for replacement pressure...BFB
JFLemieux73 Oct 14, 2021
2c32426
Replacement pressure now used on calc...BFB
JFLemieux73 Oct 14, 2021
cfcec17
eta is now used in the calc...BFB
JFLemieux73 Oct 14, 2021
c3f8e29
Ktens included in zeta,eta and rep_prs...roundoff errors...BFB only i…
JFLemieux73 Oct 14, 2021
e24cf6b
Small modif to viscous coeff subroutine for passing logical capping
JFLemieux73 Oct 14, 2021
932197e
Modifs to viscous_coeffs suroutine to be used for VP and EVP
JFLemieux73 Oct 14, 2021
9a87e52
Further modifs to viscous_coeff subroutine
JFLemieux73 Oct 14, 2021
af1b43b
cosmetic change: order of calc is ne-nw-sw-se
JFLemieux73 Oct 15, 2021
0aa8cc9
vp solver also uses viscous_coeffs_and_rep_pressure subroutine
JFLemieux73 Oct 15, 2021
1eb3096
Minor modifications following PR review
JFLemieux73 Oct 20, 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
26 changes: 15 additions & 11 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 @@ -638,13 +638,13 @@ subroutine stress (nx_block, ny_block, &
i, j, ij

real (kind=dbl_kind) :: &
divune, divunw, divuse, divusw , & ! divergence
tensionne, tensionnw, tensionse, tensionsw, & ! tension
shearne, shearnw, shearse, shearsw , & ! shearing
Deltane, Deltanw, Deltase, Deltasw , & ! Delt
zetax2ne, zetax2nw, zetax2se, zetax2sw , & ! 2 x zeta (visc coeff)
etax2ne, etax2nw, etax2se, etax2sw , & ! 2 x eta (visc coeff)
rep_prsne, rep_prsnw, rep_prsse, rep_prssw, & ! replacement pressure
divune, divunw, divusw, divuse , & ! divergence
tensionne, tensionnw, tensionsw, tensionse, & ! tension
shearne, shearnw, shearsw, shearse , & ! shearing
Deltane, Deltanw, Deltasw, Deltase , & ! Delt
zetax2ne, zetax2nw, zetax2sw, zetax2se , & ! 2 x zeta (visc coeff)
etax2ne, etax2nw, etax2sw, etax2se , & ! 2 x eta (visc coeff)
rep_prsne, rep_prsnw, rep_prssw, rep_prsse, & ! replacement pressure

@phil-blain phil-blain Oct 19, 2021

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

These changes to the argument declaration are not strictly needed and add noise to the diff, personally I would drop them.

@phil-blain phil-blain Oct 19, 2021

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Same thing for the arguments to viscous_coeffs_and_rep_pressure below

@JFLemieux73 JFLemieux73 Oct 20, 2021

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

ok corrected.

! puny , & ! puny
ssigpn, ssigps, ssigpe, ssigpw , &
ssigmn, ssigms, ssigme, ssigmw , &
Expand All @@ -656,13 +656,16 @@ subroutine stress (nx_block, ny_block, &
str12ew, str12we, str12ns, str12sn , &
strp_tmp, strm_tmp, tmp

logical :: capping ! of the viscous coef

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

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

str(:,:,:) = c0
capping = .true. ! to be improved

@phil-blain phil-blain Oct 19, 2021

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Do we want to keep that comment ? as it stands I find it a little confusing, i.e. it could mean either "the fact that capping is hardcoded here could be improved" or "it would be an improvement to allow capping = .false. here"...

@JFLemieux73 JFLemieux73 Oct 20, 2021

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

replaced by: ! could be later included in ice_in


do ij = 1, icellt
i = indxti(ij)
Expand Down Expand Up @@ -694,13 +697,14 @@ subroutine stress (nx_block, ny_block, &

call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),&
Deltane, Deltanw, &
Deltase, Deltasw, &
Deltasw, Deltase, &
zetax2ne, zetax2nw, &
zetax2se, zetax2sw, &
zetax2sw, zetax2se, &
etax2ne, etax2nw, &
etax2se, etax2sw, &
etax2sw, etax2se, &
rep_prsne, rep_prsnw, &
rep_prsse, rep_prssw )
rep_prssw, rep_prsse, &
capping)

!-----------------------------------------------------------------
! the stresses ! kg/s^2
Expand Down
57 changes: 34 additions & 23 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 @@ -1375,53 +1375,64 @@ end subroutine strain_rates

subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, &
Deltane, Deltanw, &
Deltase, Deltasw, &
Deltasw, Deltase, &
zetax2ne, zetax2nw, &
zetax2se, zetax2sw, &
zetax2sw, zetax2se, &
etax2ne, etax2nw, &
etax2se, etax2sw, &
etax2sw, etax2se, &
rep_prsne, rep_prsnw,&
rep_prsse, rep_prssw )
rep_prssw, rep_prsse,&
capping)

real (kind=dbl_kind), intent(in):: &
strength, tinyarea ! at the t-point

real (kind=dbl_kind), intent(in):: &
Deltane, Deltanw, Deltase, Deltasw ! Delta at each corner
Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner

logical , intent(in):: capping

real (kind=dbl_kind), intent(out):: &
zetax2ne, zetax2nw, zetax2se, zetax2sw, & ! zetax2 at each corner
etax2ne, etax2nw, etax2se, etax2sw, & ! etax2 at each corner
rep_prsne, rep_prsnw, rep_prsse, rep_prssw ! replacement pressure
zetax2ne, zetax2nw, zetax2sw, zetax2se, & ! zetax2 at each corner
etax2ne, etax2nw, etax2sw, etax2se, & ! etax2 at each corner
rep_prsne, rep_prsnw, rep_prssw, rep_prsse ! replacement pressure

! local variables
real (kind=dbl_kind) :: &
tmpcalc
tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse

! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code

! if (trim(yield_curve) == 'ellipse') then

tmpcalc = strength/max(Deltane,tinyarea) ! northeast
zetax2ne = (c1+Ktens)*tmpcalc
rep_prsne = (c1-Ktens)*tmpcalc*Deltane
if (capping) then
tmpcalcne = strength/max(Deltane,tinyarea)
tmpcalcnw = strength/max(Deltanw,tinyarea)
tmpcalcsw = strength/max(Deltasw,tinyarea)
tmpcalcse = strength/max(Deltase,tinyarea)
else
tmpcalcne = strength/(Deltane + tinyarea)
tmpcalcnw = strength/(Deltanw + tinyarea)
tmpcalcsw = strength/(Deltasw + tinyarea)
tmpcalcse = strength/(Deltase + tinyarea)
endif

zetax2ne = (c1+Ktens)*tmpcalcne ! northeast
rep_prsne = (c1-Ktens)*tmpcalcne*Deltane
etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot

tmpcalc = strength/max(Deltanw,tinyarea) ! northwest
zetax2nw = (c1+Ktens)*tmpcalc
rep_prsnw = (c1-Ktens)*tmpcalc*Deltanw
zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest
rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw
etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot

tmpcalc = strength/max(Deltase,tinyarea) ! southeast
zetax2se = (c1+Ktens)*tmpcalc
rep_prsse = (c1-Ktens)*tmpcalc*Deltase
etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot

tmpcalc = strength/max(Deltasw,tinyarea) ! southwest
zetax2sw = (c1+Ktens)*tmpcalc
rep_prssw = (c1-Ktens)*tmpcalc*Deltasw
zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest
rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw
etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot

zetax2se = (c1+Ktens)*tmpcalcse ! southeast
rep_prsse = (c1-Ktens)*tmpcalcse*Deltase
etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot

! else

! endif
Expand Down
Loading