Skip to content

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

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 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 all 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
100 changes: 55 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
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
! 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,52 @@ 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, &
zetax2ne, zetax2nw, &
zetax2se, zetax2sw, &
etax2ne, etax2nw, &
etax2se, etax2sw, &
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 2 x zeta and 2 x eta are used in the code

stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + &
arlx1i*(zetax2ne*divune - rep_prsne)) * denom1
stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + &
arlx1i*(zetax2nw*divunw - rep_prsnw)) * denom1
stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + &
arlx1i*(zetax2sw*divusw - rep_prssw)) * denom1
stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + &
arlx1i*(zetax2se*divuse - rep_prsse)) * denom1

stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + &
arlx1i*etax2ne*tensionne) * denom1
stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + &
arlx1i*etax2nw*tensionnw) * denom1
stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + &
arlx1i*etax2sw*tensionsw) * denom1
stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + &
arlx1i*etax2se*tensionse) * denom1

stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + &
arlx1i*p5*etax2ne*shearne) * denom1
stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + &
arlx1i*p5*etax2nw*shearnw) * denom1
stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + &
arlx1i*p5*etax2sw*shearsw) * denom1
stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + &
arlx1i*p5*etax2se*shearse) * denom1

!-----------------------------------------------------------------
! Eliminate underflows.
Expand Down
72 changes: 71 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,75 @@ 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, &
zetax2ne, zetax2nw, &
zetax2se, zetax2sw, &
etax2ne, etax2nw, &
etax2se, etax2sw, &
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):: &
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

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

! 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
etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot

tmpcalc = strength/max(Deltanw,tinyarea) ! northwest
zetax2nw = (c1+Ktens)*tmpcalc
rep_prsnw = (c1-Ktens)*tmpcalc*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
etax2sw = ecci*zetax2sw ! 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`"
"etax2", "2 x eta (shear viscous coefficient)", "kg/s"
"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**", "", ""
"zetax2", "2 x zeta (bulk viscous coefficient)", "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