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

Express rheology term as a function of viscous coefficients for EVP and rEVP #639

Merged
merged 10 commits into from
Oct 8, 2021
Merged

Express rheology term as a function of viscous coefficients for EVP and rEVP #639

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
d4313d9
added subroutine to calc zeta and eta
JFLemieux73 Sep 27, 2021
d0d98e3
new subroutine now also calculates the replacement pressure
JFLemieux73 Sep 28, 2021
38be1a4
stress calc in evp are now with zeta, eta and rep_prs
JFLemieux73 Sep 28, 2021
06b4d55
corrected minor compiling issues
JFLemieux73 Sep 28, 2021
8b866d8
improved comments and added references for papers
JFLemieux73 Oct 1, 2021
c709eb3
Modif to doc for new code with viscous coefficients
JFLemieux73 Oct 4, 2021
8e495e0
Minor modifs to doc
JFLemieux73 Oct 5, 2021
85d7b71
minor correction to QC testing doc
JFLemieux73 Oct 5, 2021
0fc0150
Small corrections following Philippes comments
JFLemieux73 Oct 6, 2021
ffdd705
Changed eta,zeta to etax2,zetax2
JFLemieux73 Oct 7, 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
89 changes: 44 additions & 45 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 @@ -6,23 +6,25 @@
! See:
!
! Hunke, E. C., and J. K. Dukowicz (1997). An elastic-viscous-plastic model
! for sea ice dynamics. {\em J. Phys. Oceanogr.}, {\bf 27}, 1849--1867.
! for sea ice dynamics. J. Phys. Oceanogr., 27, 1849-1867.
!
! Hunke, E. C. (2001). Viscous-Plastic Sea Ice Dynamics with the EVP Model:
! Linearization Issues. {\em Journal of Computational Physics}, {\bf 170},
! 18--38.
! Linearization Issues. J. Comput. Phys., 170, 18-38.
!
! Hunke, E. C., and J. K. Dukowicz (2002). The Elastic-Viscous-Plastic
! Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates
! on a Sphere---Incorporation of Metric Terms. {\em Monthly Weather Review},
! {\bf 130}, 1848--1865.
! on a Sphere - Incorporation of Metric Terms. Mon. Weather Rev.,
! 130, 1848-1865.
!
! Hunke, E. C., and J. K. Dukowicz (2003). The sea ice momentum
! equation in the free drift regime. Los Alamos Tech. Rep. LA-UR-03-2219.
!
! Bouillon, S., T. Fichefet, V. Legat and G. Madec (submitted 2013). The
! revised elastic-viscous-plastic method. Ocean Modelling.
! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys.
! Oceanogr., 9, 817-846.
!
! Bouillon, S., T. Fichefet, V. Legat and G. Madec (2013). The
! elastic-viscous-plastic method revisited. Ocean Model., 71, 2-12.
!
! author: Elizabeth C. Hunke, LANL
!
! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb (LANL)
Expand Down Expand Up @@ -590,8 +592,8 @@ subroutine stress (nx_block, ny_block, &
rdg_conv, rdg_shear, &
str )

use ice_dyn_shared, only: strain_rates, deformations

use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure
integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ksub , & ! subcycling step
Expand Down Expand Up @@ -640,9 +642,10 @@ subroutine stress (nx_block, ny_block, &
tensionne, tensionnw, tensionse, tensionsw, & ! tension
shearne, shearnw, shearse, shearsw , & ! shearing
Deltane, Deltanw, Deltase, Deltasw , & ! Delt
zetane, zetanw, zetase, zetasw , & ! zeta viscous coeff
etane, etanw, etase, etasw , & ! eta viscous coeff
rep_prsne, rep_prsnw, rep_prsse, rep_prssw, & ! replacement pressure
! puny , & ! puny
c0ne, c0nw, c0se, c0sw , & ! useful combinations
c1ne, c1nw, c1se, c1sw , &
ssigpn, ssigps, ssigpe, ssigpw , &
ssigmn, ssigms, ssigme, ssigmw , &
ssig12n, ssig12s, ssig12e, ssig12w , &
Expand All @@ -669,6 +672,7 @@ subroutine stress (nx_block, ny_block, &
! strain rates
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------

call strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
Expand All @@ -685,46 +689,41 @@ subroutine stress (nx_block, ny_block, &
Deltase, Deltasw )

!-----------------------------------------------------------------
! strength/Delta ! kg/s
! viscous coefficients and replacement pressure
!-----------------------------------------------------------------
c0ne = strength(i,j)/max(Deltane,tinyarea(i,j))
c0nw = strength(i,j)/max(Deltanw,tinyarea(i,j))
c0sw = strength(i,j)/max(Deltasw,tinyarea(i,j))
c0se = strength(i,j)/max(Deltase,tinyarea(i,j))

c1ne = c0ne*arlx1i
c1nw = c0nw*arlx1i
c1sw = c0sw*arlx1i
c1se = c0se*arlx1i

c0ne = c1ne*ecci
c0nw = c1nw*ecci
c0sw = c1sw*ecci
c0se = c1se*ecci


call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),&
Deltane, Deltanw, &
Deltase, Deltasw, &
zetane, zetanw, &
zetase, zetasw, &
etane, etanw, &
etase, etasw, &
rep_prsne, rep_prsnw, &
rep_prsse, rep_prssw )

!-----------------------------------------------------------------
! the stresses ! kg/s^2
! (1) northeast, (2) northwest, (3) southwest, (4) southeast
!-----------------------------------------------------------------

stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + c1ne*(divune*(c1+Ktens) - Deltane*(c1-Ktens))) &
* denom1
stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + c1nw*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens))) &
* denom1
stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + c1sw*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens))) &
* denom1
stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + c1se*(divuse*(c1+Ktens) - Deltase*(c1-Ktens))) &
* denom1

stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + c0ne*tensionne*(c1+Ktens)) * denom1
stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + c0nw*tensionnw*(c1+Ktens)) * denom1
stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + c0sw*tensionsw*(c1+Ktens)) * denom1
stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + c0se*tensionse*(c1+Ktens)) * denom1

stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + c0ne*shearne*p5*(c1+Ktens)) * denom1
stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + c0nw*shearnw*p5*(c1+Ktens)) * denom1
stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + c0sw*shearsw*p5*(c1+Ktens)) * denom1
stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + c0se*shearse*p5*(c1+Ktens)) * denom1
! NOTE: for comp. efficiency zeta and eta in this code are 2x zeta and eta
! as defined by Hibler 1979.

stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + arlx1i*(zetane*divune - rep_prsne))*denom1
stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + arlx1i*(zetanw*divunw - rep_prsnw))*denom1
stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + arlx1i*(zetasw*divusw - rep_prssw))*denom1
stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + arlx1i*(zetase*divuse - rep_prsse))*denom1

stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + arlx1i*etane*tensionne) * denom1
stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + arlx1i*etanw*tensionnw) * denom1
stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + arlx1i*etasw*tensionsw) * denom1
stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + arlx1i*etase*tensionse) * denom1

stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + arlx1i*p5*etane*shearne) * denom1
stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + arlx1i*p5*etanw*shearnw) * denom1
stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + arlx1i*p5*etasw*shearsw) * denom1
stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + arlx1i*p5*etase*shearse) * denom1

!-----------------------------------------------------------------
! Eliminate underflows.
Expand Down
73 changes: 72 additions & 1 deletion 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,6 +26,7 @@ module ice_dyn_shared
dyn_prep1, dyn_prep2, dyn_finish, &
seabed_stress_factor_LKD, seabed_stress_factor_prob, &
alloc_dyn_shared, deformations, strain_rates, &
viscous_coeffs_and_rep_pressure, &
stack_velocity_field, unstack_velocity_field

! namelist parameters
Expand Down Expand Up @@ -864,7 +865,7 @@ end subroutine dyn_finish
!
! Lemieux, J. F., F. Dupont, P. Blain, F. Roy, G.C. Smith, G.M. Flato (2016).
! Improving the simulation of landfast ice by combining tensile strength and a
! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121.
! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121, 7354-7368.
!
! author: JF Lemieux, Philippe Blain (ECCC)
!
Expand Down Expand Up @@ -1358,6 +1359,76 @@ subroutine strain_rates (nx_block, ny_block, &

end subroutine strain_rates

!=======================================================================
! Computes viscous coefficients and replacement pressure for stress
! calculations. Note that tensile strength is included here.
!
! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys.
! Oceanogr., 9, 817-846.
!
! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by
! adding tensile strength. J. Phys. Oceanogr. 40, 185-198.
!
! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice
! by combining tensile strength and a parameterization for grounded ridges.
! J. Geophys. Res. Oceans, 121, 7354-7368.

subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, &
Deltane, Deltanw, &
Deltase, Deltasw, &
zetane, zetanw, &
zetase, zetasw, &
etane, etanw, &
etase, etasw, &
rep_prsne, rep_prsnw, &
rep_prsse, rep_prssw )

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

real (kind=dbl_kind), intent(out):: &
zetane, zetanw, zetase, zetasw, & ! zeta at each corner
etane, etanw, etase, etasw, & ! eta at each corner
rep_prsne, rep_prsnw, rep_prsse, rep_prssw ! replacement pressure

! local variables
real (kind=dbl_kind) :: &
tmpzeta

! NOTE: for comp. efficiency zeta and eta in this code are
! 2x zeta and eta as defined by Hibler 1979.

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

tmpzeta = strength/max(Deltane,tinyarea) ! northeast
zetane = (c1+Ktens)*tmpzeta
rep_prsne = (c1-Ktens)*tmpzeta*Deltane
etane = ecci*zetane ! CHANGE FOR e_plasticpot

tmpzeta = strength/max(Deltanw,tinyarea) ! northwest
zetanw = (c1+Ktens)*tmpzeta
rep_prsnw = (c1-Ktens)*tmpzeta*Deltanw
etanw = ecci*zetanw ! CHANGE FOR e_plasticpot

tmpzeta = strength/max(Deltase,tinyarea) ! southeast
zetase = (c1+Ktens)*tmpzeta
rep_prsse = (c1-Ktens)*tmpzeta*Deltase
etase = ecci*zetase ! CHANGE FOR e_plasticpot

tmpzeta = strength/max(Deltasw,tinyarea) ! southwest
zetasw = (c1+Ktens)*tmpzeta
rep_prssw = (c1-Ktens)*tmpzeta*Deltasw
etasw = ecci*zetasw ! CHANGE FOR e_plasticpot

! else

! endif

end subroutine viscous_coeffs_and_rep_pressure

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

! Load velocity components into array for boundary updates
Expand Down
5 changes: 3 additions & 2 deletions configuration/scripts/machines/Macros.daley_intel
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,9 @@ FFLAGS := -fp-model source -convert big_endian -assume byterecl -ftz -traceb
#-xHost

ifeq ($(ICE_BLDDEBUG), true)
FFLAGS += -O0 -g -check -fpe0 -ftrapuv -fp-model except -check noarg_temp_created -init=snan,arrays
# -heap-arrays 1024
FFLAGS += -O0 -g -check -fpe0 -ftrapuv -fp-model except -check noarg_temp_created
#-heap-arrays 1024
#-init=snan,arrays
else
FFLAGS += -O2
endif
Expand Down
4 changes: 3 additions & 1 deletion doc/source/cice_index.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,7 @@ either Celsius or Kelvin units).
"eps13", "a small number", "10\ :math:`^{-13}`"
"eps16", "a small number", "10\ :math:`^{-16}`"
"esno(n)", "energy of melting of snow per unit area (in category n)", "J/m\ :math:`^2`"
"eta", "viscous coefficient (shear)", "kg/s"
Copy link
Member

Choose a reason for hiding this comment

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

Do we also want to point out here that this is actually 2*eta ? (and same for zeta below)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I will think about it a bit more and might correct it later in another PR.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think it's confusing to use eta and zeta, especially if they are different from the 'classic' definitions. Maybe etax2 and zetax2? Or eta2 and zeta2, although those could also be ^2 rather than x2. The code isn't consistent with this kind of thing.

"evap", "evaporative water flux", "kg/m\ :math:`^2`/s"
"ew_boundary_type", "type of east-west boundary condition", ""
"eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36"
Expand Down Expand Up @@ -515,7 +516,6 @@ either Celsius or Kelvin units).
"print_global", "if true, print global data", "F"
"print_points", "if true, print point data", "F"
"processor_shape", "descriptor for processor aspect ratio", ""
"prs_sig", "replacement pressure", "N/m"
"Pstar", "ice strength parameter", "2.75\ :math:`\times`\ 10\ :math:`^4`\ N/m\ :math:`^2`"
"puny", "a small positive number", "1\ :math:`\times`\ 10\ :math:`^{-11}`"
"**Q**", "", ""
Expand All @@ -540,6 +540,7 @@ either Celsius or Kelvin units).
"rdg_shear", "shear for ridging", "1/s"
"real_kind", "definition of single precision real", "selected_real_kind(6)"
"refindx", "refractive index of sea ice", "1.310"
"rep_prs", "replacement pressure", "N/m"
"revp", "real(revised_evp)", ""
"restart", "if true, initialize ice state from file", "T"
"restart_age", "if true, read age restart file", ""
Expand Down Expand Up @@ -734,6 +735,7 @@ either Celsius or Kelvin units).
"yieldstress11(12, 22)", "yield stress tensor components", ""
"year_init", "the initial year", ""
"**Z**", "", ""
"zeta", "viscous coefficient (bulk)", "kg/s"
"zlvl", "atmospheric level height (momentum)", "m"
"zlvs", "atmospheric level height (scalars)", "m"
"zref", "reference height for stability", "10. m"
Expand Down
30 changes: 14 additions & 16 deletions doc/source/science_guide/sg_dynamics.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -385,7 +385,7 @@ where :math:`k_t` should be set to a value between 0 and 1 (this can
be changed at runtime with the namelist parameter ``Ktens``). The ice
strength :math:`P` is a function of the ice thickness distribution as
described in the `Icepack
Documentation<https://cice-consortium-icepack.readthedocs.io/en/master/science_guide/index.html>`_.
Documentation <https://cice-consortium-icepack.readthedocs.io/en/master/science_guide/index.html>`_.

.. _stress-vp:

Expand All @@ -405,7 +405,7 @@ An elliptical yield curve is used, with the viscosities given by
\zeta = {P(1+k_t)\over 2\Delta},

.. math::
\eta = {P(1+k_t)\over {2\Delta e^2}},
\eta = e^{-2} \zeta,

where

Expand Down Expand Up @@ -455,29 +455,28 @@ parameter less than one. Including the modification proposed by :cite:`Bouillon1
.. math::
\begin{aligned}
{\partial\sigma_1\over\partial t} + {\sigma_1\over 2T}
+ {P_R(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta} D_D, \\
{\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {P(1+k_t)\over
2Te^2\Delta} D_T,\\
+ {P_R(1-k_t)\over 2T} &=& {\zeta \over T} D_D, \\
{\partial\sigma_2\over\partial t} + {\sigma_2\over 2T} &=& {\eta \over
T} D_T,\\
{\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 2T} &=&
{P(1+k_t)\over 4Te^2\Delta}D_S.\end{aligned}
{\eta \over 2T}D_S.\end{aligned}

Once discretized in time, these last three equations are written as

.. math::
\begin{aligned}
{(\sigma_1^{k+1}-\sigma_1^{k})\over\Delta t_e} + {\sigma_1^{k+1}\over 2T}
+ {P_R^k(1-k_t)\over 2T} &=& {P(1+k_t)\over 2T\Delta^k} D_D^k, \\
{(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {P(1+k_t)\over
2Te^2\Delta^k} D_T^k,\\
+ {P_R^k(1-k_t)\over 2T} &=& {\zeta^k\over T} D_D^k, \\
{(\sigma_2^{k+1}-\sigma_2^{k})\over\Delta t_e} + {\sigma_2^{k+1}\over 2T} &=& {\eta^k \over
T} D_T^k,\\
{(\sigma_{12}^{k+1}-\sigma_{12}^{k})\over\Delta t_e} + {\sigma_{12}^{k+1}\over 2T} &=&
{P(1+k_t)\over 4Te^2\Delta^k}D_S^k,\end{aligned}
{\eta^k \over 2T}D_S^k,\end{aligned}
:label: sigdisc


where :math:`k` denotes again the subcycling step. All coefficients on the left-hand side are constant except for
:math:`P_R`. This modification compensates for the decreased efficiency of including
the viscosity terms in the subcycling. (Note that the viscosities do not
appear explicitly.) Choices of the parameters used to define :math:`E`,
the viscosity terms in the subcycling. Choices of the parameters used to define :math:`E`,
:math:`T` and :math:`\Delta t_e` are discussed in
Sections :ref:`revp` and :ref:`parameters`.

Expand Down Expand Up @@ -721,11 +720,10 @@ Introducing another numerical parameter :math:`\alpha=2T \Delta t_e ^{-1}` :cite
.. math::
\begin{aligned}
{\alpha (\sigma_1^{k+1}-\sigma_1^{k})} + {\sigma_1^{k}}
+ {P_R^k(1-k_t)} &=& {P(1+k_t)\over \Delta^k} D_D^k, \\
{\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& {P(1+k_t)\over
e^2\Delta^k} D_T^k,\\
+ {P_R^k(1-k_t)} &=& 2 \zeta^k D_D^k, \\
{\alpha (\sigma_2^{k+1}-\sigma_2^{k})} + {\sigma_2^{k}} &=& 2 \eta^k D_T^k,\\
{\alpha (\sigma_{12}^{k+1}-\sigma_{12}^{k})} + {\sigma_{12}^{k}} &=&
{P(1+k_t)\over 2e^2\Delta^k}D_S^k,\end{aligned}
\eta^k D_S^k,\end{aligned}

where as opposed to the classic EVP, the second term in each equation is at iteration :math:`k` :cite:`Bouillon13`. Also, as opposed to the classic EVP,
:math:`\Delta t_e` times the number of subcycles (or iterations) does not need to be equal to the advective time step :math:`\Delta t`.
Expand Down
5 changes: 3 additions & 2 deletions doc/source/user_guide/ug_testing.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1121,7 +1121,7 @@ If the regression comparisons fail, then you may want to run the QC test,
# From the updated sandbox
# Generate the same test case(s) as the baseline using options or namelist changes to activate new code modifications

./cice.setup -m onyx -e intel --test smoke -g gx1 -p 44x1 -testid qc_test -s qc,medium
./cice.setup -m onyx -e intel --test smoke -g gx1 -p 44x1 --testid qc_test -s qc,medium
cd onyx_intel_smoke_gx1_44x1_medium_qc.qc_test
# modify ice_in to activate the namelist options that were determined above
./cice.build
Expand All @@ -1130,7 +1130,8 @@ If the regression comparisons fail, then you may want to run the QC test,
# Wait for runs to finish
# Perform the QC test

cp configuration/scripts/tests/QC/cice.t-test.py
# From the updated sandbox
cp configuration/scripts/tests/QC/cice.t-test.py .
./cice.t-test.py /p/work/turner/CICE_RUNS/onyx_intel_smoke_gx1_44x1_medium_qc.qc_base \
/p/work/turner/CICE_RUNS/onyx_intel_smoke_gx1_44x1_medium_qc.qc_test

Expand Down