Skip to content

Implementation of plastic potential for VP and EVP. #649

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 8 commits into from
Nov 9, 2021
Merged

Implementation of plastic potential for VP and EVP. #649

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
f7b9f7b
Implementation of plastic potential
JFLemieux73 Oct 27, 2021
46c0d7c
Modified ice_in
JFLemieux73 Oct 27, 2021
07d6989
Corrected minor issue for variable declaration
JFLemieux73 Oct 27, 2021
58ee4cd
Modifs to the doc for plastic potential
JFLemieux73 Oct 27, 2021
8fcc93e
Minor modifs to the doc
JFLemieux73 Oct 27, 2021
91643ae
More details in doc about the plastic potential
JFLemieux73 Oct 28, 2021
c1b0297
Added journal TC=the cryosphere in master_list.bib
JFLemieux73 Oct 28, 2021
6b2abeb
Modified the doc following comments from Elizabeth and Philippe
JFLemieux73 Oct 31, 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
32 changes: 17 additions & 15 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 @@ -63,8 +63,11 @@ module ice_dyn_shared

real (kind=dbl_kind), public :: &
revp , & ! 0 for classic EVP, 1 for revised EVP
e_ratio , & ! e = EVP ellipse aspect ratio
ecci , & ! 1/e^2
e_yieldcurve, & ! VP aspect ratio of elliptical yield curve
e_plasticpot, & ! VP aspect ratio of elliptical plastic potential
epp2i , & ! 1/(e_plasticpot)^2
e_factor , & ! (e_yieldcurve)^2/(e_plasticpot)^4
ecci , & ! temporary for 1d evp
dtei , & ! 1/dte, where dte is subcycling timestep (1/s)
! dte2T , & ! dte/2T
denom1 ! constants for stress equation
Expand Down Expand Up @@ -220,7 +223,6 @@ subroutine set_evp_parameters (dt)

!real (kind=dbl_kind) :: &
!dte , & ! subcycling timestep for EVP dynamics, s
!ecc , & ! (ratio of major to minor ellipse axes)^2
!tdamp2 ! 2*(wave damping time scale T)

character(len=*), parameter :: subname = '(set_evp_parameters)'
Expand All @@ -230,10 +232,10 @@ subroutine set_evp_parameters (dt)
!dtei = c1/dte ! 1/s
dtei = real(ndte,kind=dbl_kind)/dt

! major/minor axis length ratio, squared
!ecc = e_ratio**2
!ecci = c1/ecc ! 1/ecc
ecci = c1/e_ratio**2 ! 1/ecc
! variables for elliptical yield curve and plastic potential
epp2i = c1/e_plasticpot**2
Copy link
Copy Markdown
Member

@phil-blain phil-blain Oct 29, 2021

Choose a reason for hiding this comment

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

not sure I like that name.. but I don't have anything much better to propose :P

Copy link
Copy Markdown
Contributor Author

@JFLemieux73 JFLemieux73 Oct 30, 2021

Choose a reason for hiding this comment

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

At least it looks like the previous variable ecci...

e_factor = e_yieldcurve**2 / e_plasticpot**4
ecci = c1/e_yieldcurve**2 ! temporary for 1d evp

! constants for stress equation
!tdamp2 = c2*eyc*dt ! s
Expand Down Expand Up @@ -1352,10 +1354,10 @@ subroutine strain_rates (nx_block, ny_block, &
- cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j )

! Delta (in the denominator of zeta, eta)
Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2))
Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2))
Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2))
Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2))
Deltane = sqrt(divune**2 + e_factor*(tensionne**2 + shearne**2))
Deltanw = sqrt(divunw**2 + e_factor*(tensionnw**2 + shearnw**2))
Deltasw = sqrt(divusw**2 + e_factor*(tensionsw**2 + shearsw**2))
Deltase = sqrt(divuse**2 + e_factor*(tensionse**2 + shearse**2))

end subroutine strain_rates

Expand Down Expand Up @@ -1419,19 +1421,19 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, &

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

zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest
rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw
etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot
etax2nw = epp2i*zetax2nw

zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest
rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw
etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot
etax2sw = epp2i*zetax2sw

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

! else

Expand Down
5 changes: 4 additions & 1 deletion cicecore/cicedynB/dynamics/ice_dyn_vp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ module ice_dyn_vp
use ice_domain, only: nblocks, distrb_info
use ice_domain_size, only: max_blocks
use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, &
ecci, cosw, sinw, fcor_blk, uvel_init, vvel_init, &
cosw, sinw, fcor_blk, uvel_init, vvel_init, &
seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, &
seabed_stress, Ktens, stack_velocity_field, unstack_velocity_field
use ice_fileunits, only: nu_diag
Expand Down Expand Up @@ -1320,6 +1320,9 @@ end subroutine calc_zeta_dPr

! Computes the VP stresses (as diagnostic)

! Lemieux, J.-F., and Dupont, F. (2020), On the calculation of normalized
! viscous-plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769,

subroutine stress_vp (nx_block , ny_block , &
icellt , &
indxti , indxtj , &
Expand Down
32 changes: 18 additions & 14 deletions cicecore/cicedynB/general/ice_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,9 +105,10 @@ subroutine input_data
use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
evp_algorithm, &
seabed_stress, seabed_stress_method, &
k1, k2, alphab, threshold_hw, &
Ktens, e_ratio, coriolis, ssh_stress, &
kridge, brlx, arlx
k1, k2, alphab, threshold_hw, Ktens, &
e_yieldcurve, e_plasticpot, coriolis, &
ssh_stress, kridge, brlx, arlx

use ice_dyn_vp, only: maxits_nonlin, precond, dim_fgmres, dim_pgmres, maxits_fgmres, &
maxits_pgmres, monitor_nonlin, monitor_fgmres, &
monitor_pgmres, reltol_nonlin, reltol_fgmres, reltol_pgmres, &
Expand Down Expand Up @@ -208,20 +209,20 @@ subroutine input_data

namelist /dynamics_nml/ &
kdyn, ndte, revised_evp, yield_curve, &
evp_algorithm, &
evp_algorithm, &
brlx, arlx, ssh_stress, &
advection, coriolis, kridge, ktransport, &
kstrength, krdg_partic, krdg_redist, mu_rdg, &
e_ratio, Ktens, Cf, seabed_stress, &
k1, maxits_nonlin, precond, dim_fgmres, &
e_yieldcurve, e_plasticpot, Ktens, &
maxits_nonlin, precond, dim_fgmres, &
dim_pgmres, maxits_fgmres, maxits_pgmres, monitor_nonlin, &
monitor_fgmres, monitor_pgmres, reltol_nonlin, reltol_fgmres, &
reltol_pgmres, algo_nonlin, dim_andacc, reltol_andacc, &
damping_andacc, start_andacc, fpfunc_andacc, use_mean_vrel, &
ortho_type, &
k2, alphab, threshold_hw, &
seabed_stress_method, Pstar, Cstar

ortho_type, seabed_stress, seabed_stress_method, &
k1, k2, alphab, threshold_hw, &
Cf, Pstar, Cstar
namelist /shortwave_nml/ &
shortwave, albedo_type, &
albicev, albicei, albsnowv, albsnowi, &
Expand Down Expand Up @@ -367,7 +368,8 @@ subroutine input_data
alphab = 20.0_dbl_kind ! alphab=Cb factor in Lemieux et al 2015
threshold_hw = 30.0_dbl_kind ! max water depth for grounding
Ktens = 0.0_dbl_kind ! T=Ktens*P (tensile strength: see Konig and Holland, 2010)
e_ratio = 2.0_dbl_kind ! VP ellipse aspect ratio
e_yieldcurve = 2.0_dbl_kind ! VP aspect ratio of elliptical yield curve
e_plasticpot = 2.0_dbl_kind ! VP aspect ratio of elliptical plastic potential
maxits_nonlin = 4 ! max nb of iteration for nonlinear solver
precond = 'pgmres' ! preconditioner for fgmres: 'ident' (identity), 'diag' (diagonal), 'pgmres' (Jacobi-preconditioned GMRES)
dim_fgmres = 50 ! size of fgmres Krylov subspace
Expand Down Expand Up @@ -729,7 +731,8 @@ subroutine input_data
call broadcast_scalar(alphab, master_task)
call broadcast_scalar(threshold_hw, master_task)
call broadcast_scalar(Ktens, master_task)
call broadcast_scalar(e_ratio, master_task)
call broadcast_scalar(e_yieldcurve, master_task)
call broadcast_scalar(e_plasticpot, master_task)
call broadcast_scalar(advection, master_task)
call broadcast_scalar(conserv_check, master_task)
call broadcast_scalar(shortwave, master_task)
Expand Down Expand Up @@ -1440,9 +1443,10 @@ subroutine input_data
if (kdyn == 1 .or. kdyn == 3) then
write(nu_diag,1030) ' yield_curve = ', trim(yield_curve)
if (trim(yield_curve) == 'ellipse') &
write(nu_diag,1002) ' e_ratio = ', e_ratio, ' : aspect ratio of ellipse'
write(nu_diag,1002) ' e_yieldcurve = ', e_yieldcurve, ' : aspect ratio of yield curve'
write(nu_diag,1002) ' e_plasticpot = ', e_plasticpot, ' : aspect ratio of plastic potential'
endif

if (trim(coriolis) == 'latitude') then
tmpstr2 = ' : latitude-dependent Coriolis parameter'
elseif (trim(coriolis) == 'contant') then
Expand Down
3 changes: 2 additions & 1 deletion configuration/scripts/ice_in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,8 @@
Cstar = 20
Cf = 17.
Ktens = 0.
e_ratio = 2.
e_yieldcurve = 2.
e_plasticpot = 2.
seabed_stress = .false.
seabed_stress_method = 'LKD'
k1 = 7.5
Expand Down
2 changes: 1 addition & 1 deletion configuration/scripts/options/set_nml.alt03
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ sw_dtemp = 0.02d0
tfrz_option = 'linear_salt'
revised_evp = .false.
Ktens = 0.
e_ratio = 2.
e_yieldcurve = 2.
seabed_stress = .true.
use_bathymetry = .true.
l_mpond_fresh = .true.
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 @@ -205,7 +205,9 @@ either Celsius or Kelvin units).
"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"
"eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36"
"e_yieldcurve", "yield curve minor/major axis ratio", "2"
"e_plasticpot", "plastic potential minor/major axis ratio", "2"
"**F**", "", ""
"faero_atm", "aerosol deposition rate", "kg/m\ :math:`^2`/s"
"faero_ocn", "aerosol flux to the ocean", "kg/m\ :math:`^2`/s"
Expand Down
13 changes: 12 additions & 1 deletion doc/source/master_list.bib
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ @string{CRST
@string{IJHPCA={Int. J High Perform. Comput. Appl}}
@string{PTRSA={Philos. Trans. Royal Soc. A}}
@string{SIAMJCP={SIAM J. Sci. Comput.}}

@string{TC={The Cryosphere}}


% **********************************************
Expand Down Expand Up @@ -979,6 +979,17 @@ @article{Roach18
volume = {123},
year = {2018}
}

@Article{Ringeisen21
author = "D. Ringeisen and L.B. Tremblay and M. Losch",
title = "{Non-normal flow rules affect fracture angles in sea ice viscous-plastic rheologies}",
journal = TC,
year = {2021},
volume = {15},
pages = {2873-2888},
url = {https://doi.org/10.5194/tc-15-2873-2021}
}

@Article{Tsujino18,
author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.",
title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}",
Expand Down
Loading