= SI3-10_EAP_rheology Last edition: '''[[Wikinfo(changed_ts)]]''' by '''[[Wikinfo(changed_by)]]''' The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected. [[PageOutline(2, , inline)]] == Summary ||=Action || SI3-10_EAP_rheology || ||=PI(S) || Stefanie Rynders || ||=Digest ||Add new module icedyn_rhg_eap.F90 with similar structure to the existing icedyn_rhg_evp.F90 module to include EAP rheology as an option in SI3. || ||=Dependencies || || ||=Branch || source:/NEMO/branches/{YEAR}/dev_r{REV}_{ACTION_NAME} || ||=Previewer(s) || Andrew Coward || ||=Reviewer(s) || Clement Rousset? || ||=Ticket || #2325 || == Description Current EVP rheology assumes ice is isotropic, but observations show linear kinematic features (LKF) which are sign of anisotropic behaviour. Elasitc anisotropic plastic (EAP) rheology has been developed by UCL to address this shortcoming (refs 1-3 theorical development, refs 4-5 for implemation in CICE). It is routinely used in CICE simulations and shows more realistic LKF patterns.[[BR]] The IMMERSE project has a milestone and deliverable to implement EAP rheology in SI3. == Implementation === Development cautions === - EAP used in CICE with explicit solver, hasn't been tested with implicit solver (check CICE manual) - EAP hasn't been used with adaptive solver (Kimmritz) ' stability unknown ' not used for now - EAP hasn't been used with Lemieux landfast ice changes ' stability unknown ' not used for now - need 8 functions, though coding rules say it is to be avoided - all recycled CICE code will need to be rewritten to match NEMO coding conventions. === Module list === modules to be changed: {{{ trunk/src/ICE/icedyn_adv_pra.F90 trunk/src/ICE/icedyn_adv_umx.F90 trunk/src/ICE/icedyn_rhg.F90 }}} new module {{{ trunk/src/ICE/icedyn_rhg_eap.F90 }}} === Code modifications === EAP code based on Heorton et al. 2018 implementation in CICE '''icedyn_rhg.F90''' This is the main rheology module, contains two subroutines: ice_dyn_rhg and ice_dyn_rhg_init. Some placeholders for EAP rheology already defined The subroutine ice_dyn_rhg does timing and conservation checks, calls the rheologies (select case pattern) and calls the writing of restart files The subroutine ice_dyn_rhg_init reads the rheology parameters from the namelist and writes to the log no namelist parameters needed for eap specifically, need to use nn_evp and rn_relast for the subcycling ''actions'' - select case parameter np_rhg_EAP to choose EAP rheology, defined on line 35 to be uncommented - namelist parameter to be added for EAP after line 38 {{{ LOGICAL :: ln_rhg_EAP ! EAP rheology }}} - np_rhg_EAP used to switch on eap rheology on line 139 to be uncommented - namelist parameter to be added in line 112 {{{ NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_EAP }}} - line of log output to be added for eap, proposed after line 133 ? {{{ WRITE(numout,*) ' rheology EAP (icedyn_rhg_eap) ln_rhg_EAP = ', ln_rhg_EAP }}} - to be added USE statement, after line 18 {{{ USE icedyn_rhg_eap ! sea-ice: EAP rheology }}} - to be added CASE option, after line 82 {{{ CASE ( np_rhgEAP) CALL ice_dyn_rhg_eap( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) }}} - to be added call to restart file writing in ice_dyn_rhg, after line 86 {{{ IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'WRITE', kt ) }}} - to be added call to restart file reading in ice_dyn_rhg_init, after line 142 {{{ IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'READ' ) !* read or initialize all required files }}} '''icedyn_rhg_eap.F90''' new module to be added The structure will be the same as icedyn_rhg_evp.F90, with some code and extra functions added The subroutine ice_dyn_rhg_eap contains the stress calculation and momentum balance and subcycling * compute mask at F point: no changes needed * compute ice, snow mass + ice strength: initialisation of variables to be added: uses 8 functions for the variables w1, w2, s11kr, s12kr, s22r, s11ks, s12ks, s22ks * keep landfast ice option, but test later unsure whether extra tensile strength will be stable hasn't been used with EAP in other models * compute terms of momentum equation: wind/ocean stress, mass terms, coriolis term: no changes needed * solve momentum equation iteratively: use Bouillon option (more stable), possibly test Kimmritz option later * insert stress update: insert structure tension evolution (stepa subroutine), line 688 unsure whether this will be stable with implicit solver instead of explicit one ->ask Harry (note in stepa says it is implicit, clarify) * calls calc_ffrac 4 times, 4 corners of cell and average: ( warning: s11, s22 not defined on same point as s12 in SI3 need 4 corners, check?) * calculate delta, shear, divergence for redistribution: no changes needed * diagnostics: add output of anisotropy with iom_put calls need to define fields in xml files * The subroutine rhg_eap_rst: add structure tensor components (anistropy) restart writing/reading variables a11_1, a11_2, a11_3, a11_4 and a12_1, a12_2, a12_3, a12_4; pay attention to order of reads 8 functions for yield surface lookup table inside functions everything scalar * calc_ffrac subroutine? uses stresses and anisotropy tensor (both defined on corners), but inside function everything scalar - separate subroutine will not be needed if only one grid point is enough '''icedyn_adv_pra.F90''' The first advection module, containing tracers: add anisotropy tracer '''icedyn_adv_umx.F90''' The second advection module, containing tracers: add anisotropy tracer === References === 1) Wilchinsky, A.V. and Feltham, D.L. A continuum anisotropic model of sea-ice dynamics. Proc. R. Soc. Lond. A, 460(2047): 2105-2140, July 2004. doi: 10.1098/rspa.2004.1282 . 2) Wilchinsky, A.V. and Feltham, D.L. Dependence of Sea Ice Yield-Curve Shape on Ice Thickness. J. Phys. Oceanogr., 34(12): 2852-2856, 2004. doi: 10.1175/JPO2667.1 3) Wilchinsky, A.V. and Feltham, D.L. Modelling the rheology of sea ice as a collection of diamond-shaped floes. J. Nonnewton. Fluid Mech., 138(1): 22-32, September 2006. doi: 10.1016/j.jnnfm.2006.05.001. 4) Tsamados, M.; Feltham, D. L., and Wilchinsky, A. V. Impact of a new anisotropic rheology on simulations of Arctic sea ice. J. Geophys. Res., 118 (1):91-107, January 2013. doi: 10.1029/2012JC007990. 5) Heorton, H. D. B. S. and Feltham, D. L. and Tsamados, M. Stress and deformation characteristics of sea ice in a high-resolution, anisotropic sea ice model. Phil. Trans. R. Soc. A, 376(2129), 2018. doi: 10.1098/rsta.2017.0349 === Test case === We wil use the ice_adv2d testcase modified to examine the effect of the rheology by imposing diffent controlled boudary conditions. The approach is inspired by the idealised simulations in Heorton et al. 2018. === Documentation updates Section on EAP rheology to be included SI3 reference manual section 3 : ice dynamics. == Preview {{{#!box width=50em info [[Include(wiki:Developers/DevProcess#preview_)]] }}} update: EAP rheology has an impact on the closing rate for ridging/rafting. Currently this is hardcoded appropriate for EVP in icedyn_rdgrft.F90. An option will be added in icedyn_rdgrft.F90 to use EAP rate instead. update 2: Current implementation in CICE ignores the impact of advection on the anisotropy. Initial implementation in SI3 will do the same, so no extra tracer is needed for now. This can be reassed later. == Tests {{{#!box width=50em info [[Include(wiki:Developers/DevProcess#tests)]] }}} ''...'' == Review {{{#!box width=50em info [[Include(wiki:Developers/DevProcess#review)]] }}} ''...''