From 5bb179ddc27b1bbb84884776e2523b71c59b08b2 Mon Sep 17 00:00:00 2001 From: Martin Losch <30285667+mjlosch@users.noreply.github.com> Date: Thu, 17 Oct 2024 20:00:27 +0200 Subject: [PATCH] landfast ice parameterisation with lateral drag (#741) * cleaned up version of Yuqing Liu's lateral drag parameterisation - Liu et al 2021, JGR, A new parameterization of coastal drag to simulate landfast ice in deep marginal seas in the Arctic * update SEAICE_OPTIONS.h and compile SEAICE_ALLOW_SIDEDRAG code * add new file required for new parameterization * add diagnostics of lateral drag terms * add minimal documenation for basal and lateral drag by listing the runtime parameters and add references * make lateral drag work with EVP, JFNK, and Krylov solvers * fix a typo * add minimal documentation for some new fields and parameters * check for more erroneously set parameters * make it work with TAF, more store directives and include code -compile new code in offline_exf_seaice.ad without using it * more consistent names and some extra explanations * redefine the position of u/vCoastLineFile to be at U/V-points - add some more explanation in SEAICE_PARAMS.h - remove confusing averaging but keep normalisation * add missing variable declaration * remove unnecessary extra error message * update changed variable names in taf-store directive spotted by jm-c (thanks!) * more consistent upper/lower case * tweak some comments and description * document adding landfast ice parameterisation --------- Co-authored-by: Jean-Michel Campin Co-authored-by: Jean-Michel Campin --- doc/manual_references.bib | 27 +++ doc/phys_pkgs/seaice.rst | 15 ++ doc/tag-index | 2 + pkg/seaice/SEAICE.h | 14 ++ pkg/seaice/SEAICE_OPTIONS.h | 4 + pkg/seaice/SEAICE_PARAMS.h | 37 ++++- pkg/seaice/seaice_ad_check_lev1_dir.h | 4 + pkg/seaice/seaice_ad_check_lev2_dir.h | 4 + pkg/seaice/seaice_ad_check_lev3_dir.h | 4 + pkg/seaice/seaice_ad_check_lev4_dir.h | 4 + pkg/seaice/seaice_ad_diff.list | 1 + pkg/seaice/seaice_calc_lhs.F | 8 + pkg/seaice/seaice_calc_residual.F | 12 +- pkg/seaice/seaice_check.F | 30 +++- pkg/seaice/seaice_diagnostics_init.F | 18 ++ pkg/seaice/seaice_dynsolver.F | 30 ++++ pkg/seaice/seaice_evp.F | 17 ++ pkg/seaice/seaice_init_fixed.F | 73 +++++++- pkg/seaice/seaice_init_varia.F | 4 + pkg/seaice/seaice_lsr.F | 28 ++++ pkg/seaice/seaice_readparms.F | 15 ++ pkg/seaice/seaice_sidedrag_stress.F | 156 ++++++++++++++++++ pkg/seaice/seaice_summary.F | 8 + .../code_ad/MDSIO_BUFF_WH.h | 2 +- .../code_ad/SEAICE_OPTIONS.h | 4 + .../seaice_obcs/code/SEAICE_OPTIONS.h | 4 + 26 files changed, 508 insertions(+), 17 deletions(-) create mode 100644 pkg/seaice/seaice_sidedrag_stress.F diff --git a/doc/manual_references.bib b/doc/manual_references.bib index 92470ca3cf..b79fdee9d0 100644 --- a/doc/manual_references.bib +++ b/doc/manual_references.bib @@ -1298,6 +1298,19 @@ @Article{lemieux:12 pages = {5926--5944}, doi = {10.1016/j.jcp.2012.05.024}} +@Article{lemieux:15, + author = {Jean-Fran{\c{c}}ois Lemieux and L. Bruno Tremblay, + Fr{\'e}d{\'e}ric Dupont and Mathieu Plante and + Gregory C. Smith and Dany Dumont}, + title = {A Basal Stress Parameterization for Modeling + Land-Fast Ice}, + journal = jgr, + year = 2015, + volume = 120, + pages = {3157--3173}, + doi = {10.1002/2014JC010678} +} + @Article{leppaeranta:83, author = {Lepp{\"a}ranta, M.}, title = {A Growth Model for Black Ice, Snow Ican and Snow Thickness in Subarctic Basins}, @@ -1338,6 +1351,20 @@ @Article{lhans:74 pages = {118-133} } +@Article{liu:22, + author = {Yuqing Liu and Martin Losch and Nils Hutter and + Longjiang Mu}, + title = {A New Parameterization of Coastal Drag to Simulate + Landfast Ice in Deep Marginal Seas in the {A}rctic.}, + journal = jgr, + year = 2022, + volume = 127, + number = 6, + pages = {e2022JC018413}, + url = {https://doi.org/10.1029/2022JC018413}, + doi = {10.1029/2022JC018413} +} + @Article{lipscomb:01, author = {Lipscomb, W. H.}, title = {Remapping the thickness distribution in sea ice diff --git a/doc/phys_pkgs/seaice.rst b/doc/phys_pkgs/seaice.rst index b733b788d0..6f0eb547bd 100644 --- a/doc/phys_pkgs/seaice.rst +++ b/doc/phys_pkgs/seaice.rst @@ -73,6 +73,7 @@ automatically undefines more recent features, see :filelink:`SEAICE_OPTIONS.h :varlink:`SEAICE_ZETA_SMOOTHREG`, #define, use differentiable regularization for viscosities :varlink:`SEAICE_DELTA_SMOOTHREG`, #undef, use differentiable regularization :math:`\Delta_{\mathrm{reg}}=\sqrt{\Delta^2+\Delta_{\min}}` instead of :math:`\max`-function for :math:`1/\Delta_{\mathrm{reg}}` :varlink:`SEAICE_ALLOW_BOTTOMDRAG`, #undef, enable grounding parameterization for improved fastice in shallow seas + :varlink:`SEAICE_ALLOW_SIDEDRAG`, #undef, enable lateral drag parameterization for improved fastice along coastlines and islands :varlink:`SEAICE_BGRID_DYNAMICS`, #undef, use sea ice dynamics code on legacy B-grid; most of the previous flags are not available with B-grid :varlink:`SEAICE_BICE_STRESS`, #undef, B-grid only for backward compatiblity: turn on ice-stress on ocean; defined by default if :varlink:`SEAICE_BGRID_DYNAMICS` is defined :varlink:`EXPLICIT_SSH_SLOPE`, #undef, B-grid only for backward compatiblity: use ETAN for tilt computations rather than geostrophic velocities; defined by default if :varlink:`SEAICE_BGRID_DYNAMICS` is defined @@ -331,6 +332,20 @@ General flags and parameters +------------------------------------+------------------------------+-------------------------------------------------------------------------+ | :varlink:`SEAICE_useMultDimSnow` | TRUE | use same fixed pdf for snow as for multi-thickness-category ice | +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`SEAICEbasalDragK1` | 8.0 | basal drag parameter K\ :sub:`1` :cite:`lemieux:15` | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`SEAICEbasalDragK2` | 0.0 | basal drag parameter K\ :sub:`2` | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`SEAICE_cBasalStar` | :varlink:`SEAICE_cStar` value| basal drag parameter (no units) | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`SEAICEbasalDragU0` | 5.E-5 | basal drag parameter (m/s) | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`SEAICESideDrag` | 0.0 | lateral drag coefficient :cite:`liu:22` | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`uCoastLineFile` | unset | filename for coastline length for u-equation | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ + | :varlink:`vCoastLineFile` | unset | filename for coastline length for v-equation | + +------------------------------------+------------------------------+-------------------------------------------------------------------------+ The following dynamical ice thickness distribution and ridging parameters in diff --git a/doc/tag-index b/doc/tag-index index 57fce49858..d555a29a2b 100644 --- a/doc/tag-index +++ b/doc/tag-index @@ -1,6 +1,8 @@ Notes on tags used in MITgcmUV ============================== +o pkg/seaice: + - add new lateral drag parameterisation for landfast ice (Liu etal 2022, JGR) o verification/*/results: - update 5 adm + 3 tlm ref. output affected by TAF new version 6.4.5 & 6.5.0 o pkg/ecco: diff --git a/pkg/seaice/SEAICE.h b/pkg/seaice/SEAICE.h index 1fa0ae6fe2..f758dde194 100644 --- a/pkg/seaice/SEAICE.h +++ b/pkg/seaice/SEAICE.h @@ -188,6 +188,20 @@ C CbobC :: (linear) bottom drag coefficient for basals stress param. _RL CbotC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) # endif /* SEAICE_ALLOW_BOTTOMDRAG */ +# ifdef SEAICE_ALLOW_SIDEDRAG +C coastRoughU/V :: coast line roughness (w/out units) in U and V direction, +C computed from the coastline length at corner points, +C interpolated to U/V points, and scaled by the grid cell +C width +C sideDragU/V :: drag coefficients for lateral drag a parameterisation + COMMON/SEAICE_SIDEDRAG/ sideDragU, sideDragV, + & coastRoughU, coastRoughV + _RL sideDragU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL sideDragV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL coastRoughU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL coastRoughV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +# endif /* SEAICE_ALLOW_SIDEDRAG */ + # if ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) C diagnostics for the JFNK and Krylov solver INTEGER totalNewtonIters diff --git a/pkg/seaice/SEAICE_OPTIONS.h b/pkg/seaice/SEAICE_OPTIONS.h index 316661c94d..10720573ec 100644 --- a/pkg/seaice/SEAICE_OPTIONS.h +++ b/pkg/seaice/SEAICE_OPTIONS.h @@ -188,6 +188,10 @@ C iteration depends on the residual. Good for when a non-linear C convergence criterion must be satified # undef SEAICE_ALLOW_LSR_FLEX +C Use parameterisation of explicit lateral drag for a better +C representation of fastice along coast lines and islands +# undef SEAICE_ALLOW_SIDEDRAG + #endif /* SEAICE_CGRID */ #ifdef SEAICE_BGRID_DYNAMICS diff --git a/pkg/seaice/SEAICE_PARAMS.h b/pkg/seaice/SEAICE_PARAMS.h index 8c72ef1b9b..396b3d5901 100644 --- a/pkg/seaice/SEAICE_PARAMS.h +++ b/pkg/seaice/SEAICE_PARAMS.h @@ -318,12 +318,26 @@ C & SEAICE_debugPointJ C-- COMMON /SEAICE_PARM_C/ Character valued sea ice model parameters. -C AreaFile :: File containing initial sea-ice concentration -C HsnowFile :: File containing initial snow thickness -C HsaltFile :: File containing initial sea ice salt content -C HeffFile :: File containing initial sea-ice thickness -C uIceFile :: File containing initial sea-ice U comp. velocity -C vIceFile :: File containing initial sea-ice V comp. velocity +C AreaFile :: File containing initial sea-ice concentration +C HsnowFile :: File containing initial snow thickness +C HsaltFile :: File containing initial sea ice salt content +C HeffFile :: File containing initial sea-ice thickness +C uIceFile :: File containing initial sea-ice U comp. velocity +C vIceFile :: File containing initial sea-ice V comp. velocity +C uCoastLineFile :: File containing the some measure of coastline +C roughness length (in m) at the U-points in the +C X-direction (i.e. for the U-equation). +C vCoastLineFile :: Files containing the some measure of coastline +C roughness length (in m) at the V-points in the +C Y-direction (i.e. for the V-equation). +C +C This roughness length can be the subgrid +C scale length of the coastline in a grid cell +C projected in the direction normal to the u/v- +C direction as in Liu et al. (2022), but it can +C also be anything that is a good proxy of coast +C line roughness. +C C !!! NOTE !!! Initial sea-ice thickness can also be set using C SEAICE_initialHEFF below. But a constant initial condition C can mean large artificial fluxes of heat and freshwater in @@ -335,9 +349,11 @@ C CHARACTER*(MAX_LEN_FNAM) HeffFile CHARACTER*(MAX_LEN_FNAM) uIceFile CHARACTER*(MAX_LEN_FNAM) vIceFile + CHARACTER*(MAX_LEN_FNAM) uCoastLineFile + CHARACTER*(MAX_LEN_FNAM) vCoastLineFile COMMON /SEAICE_PARM_C/ & AreaFile, HsnowFile, HsaltFile, HeffFile, - & uIceFile, vIceFile + & uIceFile, vIceFile, uCoastLineFile, vCoastLineFile C-- COMMON /SEAICE_PARM_RL/ Real valued parameters of sea ice model. C SEAICE_deltaTtherm :: Seaice timestep for thermodynamic equations (s) @@ -415,6 +431,11 @@ C SEAICEbasalDragU0 (default = 5e-5) C SEAICEbasalDragK1 (default = 8) C SEAICEbasalDragK2 :: if > 0, turns on basal drag C (default = 0, Lemieux suggests 15) +C SEAICEsideDrag :: if > 0, turns on lateral static drag +C if < 0, turns on lateral quadratic drag +C both are different landfast ice parameterisations +C (Liu et al 2022 use 2e-4, +C the default = 0 turns off the parameterisations) C C SEAICE_wetAlbTemp :: Temp (deg.C) above which wet-albedo values are used C SEAICE_waterAlbedo :: water albedo @@ -508,6 +529,7 @@ C _RL SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south _RL SEAICE_cBasalStar, SEAICEbasalDragU0 _RL SEAICEbasalDragK1, SEAICEbasalDragK2 + _RL SEAICEsideDrag _RL SEAICE_wetAlbTemp, SEAICE_waterAlbedo _RL SEAICE_strength, SEAICE_cStar, SEAICEpressReplFac _RL SEAICE_tensilFac, SEAICE_tensilDepth @@ -566,6 +588,7 @@ C & SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south, & SEAICE_cBasalStar, SEAICEbasalDragU0, & SEAICEbasalDragK1, SEAICEbasalDragK2, + & SEAICEsideDrag, & SEAICE_wetAlbTemp, SEAICE_waterAlbedo, & SEAICE_strength, SEAICE_cStar, SEAICE_eccen, SEAICE_eccfr, & SEAICEtdMU, SEAICEmcMu, diff --git a/pkg/seaice/seaice_ad_check_lev1_dir.h b/pkg/seaice/seaice_ad_check_lev1_dir.h index ce280d1192..c96039a79d 100644 --- a/pkg/seaice/seaice_ad_check_lev1_dir.h +++ b/pkg/seaice/seaice_ad_check_lev1_dir.h @@ -30,6 +30,10 @@ CADJ STORE dwatn = comlev1, key=ikey_dynamics, kind=isbyte # ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE cbotc = comlev1, key=ikey_dynamics, kind=isbyte # endif +# ifdef SEAICE_ALLOW_SIDEDRAG +CADJ STORE sideDragU = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE sideDragV = comlev1, key=ikey_dynamics, kind=isbyte +# endif CADJ STORE seaicemasku = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE seaicemaskv = comlev1, key=ikey_dynamics, kind=isbyte # ifdef SEAICE_ALLOW_EVP diff --git a/pkg/seaice/seaice_ad_check_lev2_dir.h b/pkg/seaice/seaice_ad_check_lev2_dir.h index 107de2e625..aef185efa6 100644 --- a/pkg/seaice/seaice_ad_check_lev2_dir.h +++ b/pkg/seaice/seaice_ad_check_lev2_dir.h @@ -15,6 +15,10 @@ CADJ STORE VICE = tapelev2, key = ilev_2 CADJ STORE TICES = tapelev2, key = ilev_2 # ifdef SEAICE_CGRID CADJ STORE dwatn = tapelev2, key = ilev_2 +# ifdef SEAICE_ALLOW_SIDEDRAG +CADJ STORE sideDragU = tapelev2, key = ilev_2 +CADJ STORE sideDragV = tapelev2, key = ilev_2 +# endif CADJ STORE stressDivergenceX = tapelev2, key = ilev_2 CADJ STORE stressDivergenceY = tapelev2, key = ilev_2 # ifdef SEAICE_ALLOW_EVP diff --git a/pkg/seaice/seaice_ad_check_lev3_dir.h b/pkg/seaice/seaice_ad_check_lev3_dir.h index aed2dcc9d0..a200472f8b 100644 --- a/pkg/seaice/seaice_ad_check_lev3_dir.h +++ b/pkg/seaice/seaice_ad_check_lev3_dir.h @@ -15,6 +15,10 @@ CADJ STORE VICE = tapelev3, key = ilev_3 CADJ STORE TICES = tapelev3, key = ilev_3 # ifdef SEAICE_CGRID CADJ STORE dwatn = tapelev3, key = ilev_3 +# ifdef SEAICE_ALLOW_SIDEDRAG +CADJ STORE sideDragU = tapelev3, key = ilev_3 +CADJ STORE sideDragV = tapelev3, key = ilev_3 +# endif CADJ STORE stressDivergenceX = tapelev3, key = ilev_3 CADJ STORE stressDivergenceY = tapelev3, key = ilev_3 # ifdef SEAICE_ALLOW_EVP diff --git a/pkg/seaice/seaice_ad_check_lev4_dir.h b/pkg/seaice/seaice_ad_check_lev4_dir.h index e603ef5528..d23c68ebad 100644 --- a/pkg/seaice/seaice_ad_check_lev4_dir.h +++ b/pkg/seaice/seaice_ad_check_lev4_dir.h @@ -15,6 +15,10 @@ CADJ STORE VICE = tapelev4, key = ilev_4 CADJ STORE TICES = tapelev4, key = ilev_4 # ifdef SEAICE_CGRID CADJ STORE dwatn = tapelev4, key = ilev_4 +# ifdef SEAICE_ALLOW_SIDEDRAG +CADJ STORE sideDragU = tapelev4, key = ilev_4 +CADJ STORE sideDragV = tapelev4, key = ilev_4 +# endif CADJ STORE stressDivergenceX = tapelev4, key = ilev_4 CADJ STORE stressDivergenceY = tapelev4, key = ilev_4 # ifdef SEAICE_ALLOW_EVP diff --git a/pkg/seaice/seaice_ad_diff.list b/pkg/seaice/seaice_ad_diff.list index c99968e899..3568524288 100644 --- a/pkg/seaice/seaice_ad_diff.list +++ b/pkg/seaice/seaice_ad_diff.list @@ -35,3 +35,4 @@ seaice_calc_ice_strength.f seaice_reg_ridge.f seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f +seaice_sidedrag_stress.f diff --git a/pkg/seaice/seaice_calc_lhs.F b/pkg/seaice/seaice_calc_lhs.F index 9eea1c16a9..725e43bd7b 100644 --- a/pkg/seaice/seaice_calc_lhs.F +++ b/pkg/seaice/seaice_calc_lhs.F @@ -176,6 +176,14 @@ SUBROUTINE SEAICE_CALC_LHS( & + DWATN(I,J-1,bi,bj) * 0.5 _d 0 * & (uIceLoc(I,J-1,bi,bj)+uIceLoc(I+1,J-1,bi,bj)) & ) ) * areaS(I,J) +#ifdef SEAICE_ALLOW_SIDEDRAG +C lateral drag terms: + sideDragU*uIce + uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) + & + sideDragU(I,J,bi,bj)*uIceLoc(I,J,bi,bj) +C + sideDragV*vIce + vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) + & + sideDragV(I,J,bi,bj)*vIceLoc(I,J,bi,bj) +#endif /* SEAICE_ALLOW_SIDEDRAG */ C apply masks for interior (important when we have open boundaries) uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj)*maskinW(I,J,bi,bj) vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj)*maskinS(I,J,bi,bj) diff --git a/pkg/seaice/seaice_calc_residual.F b/pkg/seaice/seaice_calc_residual.F index 77558f21ed..56f9647cf8 100644 --- a/pkg/seaice/seaice_calc_residual.F +++ b/pkg/seaice/seaice_calc_residual.F @@ -25,9 +25,9 @@ SUBROUTINE SEAICE_CALC_RESIDUAL( #include "SIZE.h" #include "EEPARAMS.h" #include "SEAICE_SIZE.h" -#ifdef SEAICE_ALLOW_BOTTOMDRAG -#include "SEAICE_PARAMS.h" -#endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#if ( defined SEAICE_ALLOW_BOTTOMDRAG || defined SEAICE_ALLOW_SIDEDRAG ) +# include "SEAICE_PARAMS.h" +#endif #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: @@ -91,6 +91,12 @@ SUBROUTINE SEAICE_CALC_RESIDUAL( O CbotC, I krylovIter, myTime, myIter, myThid ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#ifdef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS( + I uIceLoc, vIceLoc, coastRoughU, coastRoughV, AREA, + O sideDragU, sideDragV, + I krylovIter, myTime, myIter, myThid ) +#endif /* SEAICE_ALLOW_SIDEDRAG */ CALL SEAICE_CALC_STRAINRATES( I uIceLoc, vIceLoc, O e11, e22, e12, diff --git a/pkg/seaice/seaice_check.F b/pkg/seaice/seaice_check.F index 4bff67dbb2..1d4941d13e 100644 --- a/pkg/seaice/seaice_check.F +++ b/pkg/seaice/seaice_check.F @@ -1169,7 +1169,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) IF ( SEAICEpressReplFac .LT. 0. _d 0 .OR. & SEAICEpressReplFac .GT. 1. _d 0 ) THEN - WRITE(msgBuf,'(A,I2)') + WRITE(msgBuf,'(A,F5.2)') & 'SEAICE_CHECK: SEAICEpressReplFac = ', SEAICEpressReplFac CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') 'SEAICE_CHECK: cannot be < 0 or > 1' @@ -1212,7 +1212,7 @@ SUBROUTINE SEAICE_CHECK( myThid ) #ifndef SEAICE_ALLOW_BOTTOMDRAG IF ( SEAICEbasalDragK2 .GT. 0. _d 0 ) THEN - WRITE(msgBuf,'(A,I2)') + WRITE(msgBuf,'(A,E12.5)') & 'SEAICE_CHECK: SEAICEbasalDragK2 = ', SEAICEbasalDragK2 CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(2A)') 'SEAICE_CHECK: is greater than 0, ', @@ -1222,6 +1222,32 @@ SUBROUTINE SEAICE_CHECK( myThid ) ENDIF #endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#ifndef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0. _d 0 ) THEN + WRITE(msgBuf,'(A,E12.5,A)') + & 'SEAICE_CHECK: SEAICEsideDrag = ', SEAICEsideDrag, + & ' is set, but SEAICE_ALLOW_SIDEDRAG is not defined' + CALL PRINT_ERROR( msgBuf, myThid ) + errCount = errCount + 1 + ENDIF + IF ( uCoastLineFile .NE. ' ' ) THEN + i = ILNBLNK(uCoastLineFile) + WRITE(msgBuf,'(3A)') + & 'SEAICE_CHECK: uCoastLineFile = ', uCoastLineFile(1:i), + & ' is set, but SEAICE_ALLOW_SIDEDRAG is not defined' + CALL PRINT_ERROR( msgBuf, myThid ) + errCount = errCount + 1 + ENDIF + IF( vCoastLineFile .NE. ' ' ) THEN + i = ILNBLNK(vCoastLineFile) + WRITE(msgBuf,'(3A)') + & 'SEAICE_CHECK: vCoastLineFile = ', vCoastLineFile(1:i), + & ' is set, but SEAICE_ALLOW_SIDEDRAG is not defined' + CALL PRINT_ERROR( msgBuf, myThid ) + errCount = errCount + 1 + ENDIF +#endif + #ifdef SEAICE_ITD C The ice thickness distribution (ITD) module can only be used with C the zero-layer thermodynamics of S/R SEAICE_GROWTH and the diff --git a/pkg/seaice/seaice_diagnostics_init.F b/pkg/seaice/seaice_diagnostics_init.F index 297b8b9028..5c14575f90 100644 --- a/pkg/seaice/seaice_diagnostics_init.F +++ b/pkg/seaice/seaice_diagnostics_init.F @@ -343,6 +343,24 @@ SUBROUTINE SEAICE_DIAGNOSTICS_INIT( myThid ) c CALL DIAGNOSTICS_ADDTOLIST( diagNum, c I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) +# ifdef SEAICE_ALLOW_SIDEDRAG + diagName = 'SIlatDgU' + diagTitle = 'SEAICE lateral coastal drag on U momentum' + diagUnits = 'N/m^2 ' + diagCode = 'UU M1 ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) + + diagName = 'SIlatDgV' + diagTitle = 'SEAICE lateral coastal drag on V momentum' + diagUnits = 'N/m^2 ' + diagCode = 'VV M1 ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) +#endif + C SIqnet, Qnet, and QNETtave are identical. C With #undef NONLIN_FRSURF SIqnet is identical to -(TFLUX-TRELAX). C Except over land and under sea ice, SIqnet is also identical to diff --git a/pkg/seaice/seaice_dynsolver.F b/pkg/seaice/seaice_dynsolver.F index 851bcdeffe..3fcd4d2d27 100644 --- a/pkg/seaice/seaice_dynsolver.F +++ b/pkg/seaice/seaice_dynsolver.F @@ -415,6 +415,36 @@ SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) CALL DIAGNOSTICS_FILL(eta ,'SIeta ',0,1,0,1,1,myThid) CALL DIAGNOSTICS_FILL(press ,'SIpress ',0,1,0,1,1,myThid) CALL DIAGNOSTICS_FILL(deltaC ,'SIdelta ',0,1,0,1,1,myThid) +# ifdef SEAICE_ALLOW_SIDEDRAG +C recompute lateral coast drag terms + IF ( DIAGNOSTICS_IS_ON('SIlatDgU',myThid) ) THEN + DO bj = myByLo(myThid), myByHi(myThid) + DO bi = myBxLo(myThid), myBxHi(myThid) +C use sig1 as a temporary field + DO j=1,sNy + DO i=1,sNx+1 + sig1(i,j) = sideDragU(i,j,bi,bj)*uIce(i,j,bi,bj) + ENDDO + ENDDO + CALL DIAGNOSTICS_FILL(sig1,'SIlatDgU',0,1,2,bi,bj,myThid) + ENDDO + ENDDO + ENDIF + IF ( DIAGNOSTICS_IS_ON('SIlatDgV',myThid) ) THEN + DO bj = myByLo(myThid), myByHi(myThid) + DO bi = myBxLo(myThid), myBxHi(myThid) +C use sig1 as a temporary field + DO j=1,sNy+1 + DO i=1,sNx + sig1(i,j) = sideDragV(i,j,bi,bj)*vIce(i,j,bi,bj) + ENDDO + ENDDO + CALL DIAGNOSTICS_FILL(sig1,'SIlatDgV',0,1,2,bi,bj,myThid) + ENDDO + ENDDO + ENDIF +# endif /* SEAICE_ALLOW_SIDEDRAG */ + IF ( DIAGNOSTICS_IS_ON('SItensil',myThid) ) THEN DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) diff --git a/pkg/seaice/seaice_evp.F b/pkg/seaice/seaice_evp.F index 95d9a88352..f9d050e2eb 100644 --- a/pkg/seaice/seaice_evp.F +++ b/pkg/seaice/seaice_evp.F @@ -323,7 +323,9 @@ SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) I iEVPstep, myTime, myIter, myThid ) #ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE e11 = comlev1_evp,key = evpkey, byte = isbyte CADJ STORE e22 = comlev1_evp,key = evpkey, byte = isbyte +CADJ STORE e12 = comlev1_evp,key = evpkey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ DO bj=myByLo(myThid),myByHi(myThid) @@ -748,11 +750,22 @@ SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) O CbotC, I iEVPstep, myTime, myIter, myThid ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#ifdef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS( + I uIce, vIce, coastRoughU, coastRoughV, AREA, + O sideDragU, sideDragV, + I iEVPstep, myTime, myIter, myThid ) +#endif /* SEAICE_ALLOW_SIDEDRAG */ + #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE DWATN = comlev1_evp, key = evpkey, byte=isbyte # ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE CbotC = comlev1_evp, key = evpkey, byte=isbyte # endif +# ifdef SEAICE_ALLOW_SIDEDRAG +CADJ STORE sideDragU = comlev1_evp, key = evpkey, byte=isbyte +CADJ STORE sideDragV = comlev1_evp, key = evpkey, byte=isbyte +# endif #endif DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) @@ -882,6 +895,10 @@ SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) denomV = denomV + areaS(i,j,bi,bj) & * 0.5 _d 0*( CbotC(i,j,bi,bj) + CbotC(i,j-1,bi,bj) ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#ifdef SEAICE_ALLOW_SIDEDRAG + denomU = denomU + sideDragU(i,j,bi,bj) + denomV = denomV + sideDragV(i,j,bi,bj) +#endif IF ( denomU .EQ. 0. _d 0 ) denomU = 1. _d 0 IF ( denomV .EQ. 0. _d 0 ) denomV = 1. _d 0 uIce(i,j,bi,bj) = seaiceMaskU(i,j,bi,bj) * diff --git a/pkg/seaice/seaice_init_fixed.F b/pkg/seaice/seaice_init_fixed.F index 246f930c9a..a3e5756444 100644 --- a/pkg/seaice/seaice_init_fixed.F +++ b/pkg/seaice/seaice_init_fixed.F @@ -1,15 +1,20 @@ #include "SEAICE_OPTIONS.h" -CStartOfInterface +CBOP +C !ROUTINE: SEAICE_INIT_FIXED + +C !INTERFACE: ========================================================== SUBROUTINE SEAICE_INIT_FIXED( myThid ) +C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE SEAICE_INIT_FIXED C | o Initialization of sea ice model. C *==========================================================* C *==========================================================* IMPLICIT NONE +C \ev -C === Global variables === +C !USES: =============================================================== #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" @@ -20,12 +25,12 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) #include "SEAICE.h" #include "SEAICE_TRACER.h" -C === Routine arguments === +C !INPUT PARAMETERS: =================================================== C myThid - Thread no. that called this routine. INTEGER myThid CEndOfInterface -C === Local variables === +C !LOCAL VARIABLES: ==================================================== #ifdef SEAICE_ITD C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf @@ -50,6 +55,10 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) #ifdef SEAICE_BGRID_DYNAMICS _RL mask_uice #endif +#ifdef SEAICE_ALLOW_SIDEDRAG + LOGICAL readCoastLineFields +#endif +CEOP IF ( usingPCoords ) THEN kSrf = Nr @@ -435,6 +444,62 @@ SUBROUTINE SEAICE_INIT_FIXED( myThid ) ENDDO ENDDO +#ifdef SEAICE_ALLOW_SIDEDRAG + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy+1,sNy+OLy-1 + DO i=1-OLx+1,sNx+OLx-1 + coastRoughU(i,j,bi,bj) = ( 3. _d 0 + & - SIMaskU(i,j-1,bi,bj) + & - SIMaskU(i,j,bi,bj) + & - SIMaskU(i,j+1,bi,bj)) * SIMaskU(i,j,bi,bj) + coastRoughV(i,j,bi,bj) = ( 3. _d 0 + & - SIMaskV(i-1,j,bi,bj) + & - SIMaskV(i,j,bi,bj) + & - SIMaskV(i+1,j,bi,bj) ) * SIMaskV(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + + readCoastLineFields = + & uCoastLineFile .NE. ' ' .OR. vCoastLineFile .NE. ' ' + IF ( readCoastLineFields ) THEN + IF ( uCoastLineFile .NE. ' ' ) THEN + CALL READ_FLD_XY_RL( uCoastLineFile,' ',coastRoughU,0,myThid ) + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx + coastRoughU(i,j,bi,bj) = coastRoughU(i,j,bi,bj) + & * SIMaskU(i,j,bi,bj) * recip_dxC(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF + IF ( vCoastLineFile .NE. ' ' ) THEN + CALL READ_FLD_XY_RL( vCoastLineFile,' ',coastRoughV,0,myThid ) + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + coastRoughV(i,j,bi,bj) = coastRoughV(i,j,bi,bj) + & * SIMaskV(i,j,bi,bj) * recip_dyC(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF + ENDIF + + CALL EXCH_UV_XY_RL( coastRoughU, coastRoughV, .FALSE., myThid ) + CALL WRITE_FLD_XY_RL( 'coastRoughUnormalized', + & ' ', coastRoughU, -1, myThid ) + CALL WRITE_FLD_XY_RL( 'coastRoughVnormalized', + & ' ', coastRoughV, -1, myThid ) +#endif /* SEAICE_ALLOW_SIDEDRAG */ + #ifdef ALLOW_SHELFICE IF ( useShelfIce ) THEN DO bj=myByLo(myThid),myByHi(myThid) diff --git a/pkg/seaice/seaice_init_varia.F b/pkg/seaice/seaice_init_varia.F index 692808de92..d4b69b61c7 100644 --- a/pkg/seaice/seaice_init_varia.F +++ b/pkg/seaice/seaice_init_varia.F @@ -115,6 +115,10 @@ SUBROUTINE SEAICE_INIT_VARIA( myThid ) # ifdef SEAICE_ALLOW_BOTTOMDRAG CbotC(i,j,bi,bj) = 0. _d 0 # endif /* SEAICE_ALLOW_BOTTOMDRAG */ +# ifdef SEAICE_ALLOW_SIDEDRAG + sideDragU(i,j,bi,bj) = 0. _d 0 + sideDragV(i,j,bi,bj) = 0. _d 0 +# endif /* SEAICE_ALLOW_SIDEDRAG */ #endif /* SEAICE_CGRID */ #ifdef SEAICE_BGRID_DYNAMICS uIceB(i,j,bi,bj) = 0. _d 0 diff --git a/pkg/seaice/seaice_lsr.F b/pkg/seaice/seaice_lsr.F index fa31bf26a4..4d33e3b93f 100644 --- a/pkg/seaice/seaice_lsr.F +++ b/pkg/seaice/seaice_lsr.F @@ -391,6 +391,17 @@ SUBROUTINE SEAICE_LSR( myTime, myIter, myThid ) O CbotC, I ipass, myTime, myIter, myThid ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ + +#ifdef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS( + I uIceC, vIceC, coastRoughU, coastRoughV, AREA, + O sideDragU, sideDragV, + I ipass, myTime, myIter, myThid ) +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE sideDragU = comlev1_dynsol, kind=isbyte, key = nlkey +CADJ STORE sideDragV = comlev1_dynsol, kind=isbyte, key = nlkey +#endif /* ALLOW_AUTODIFF_TAMC */ +#endif /* SEAICE_ALLOW_SIDEDRAG */ C C some abbreviations C @@ -532,6 +543,23 @@ SUBROUTINE SEAICE_LSR( myTime, myIter, myThid ) O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2, I iMin, iMax, jMin, jMax, myTime, myIter, myThid ) +#ifdef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0. _d 0 ) THEN + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=jMin,jMax + DO i=iMin,iMax + BU(i,j,bi,bj) = BU(i,j,bi,bj) + & + seaiceMaskU(i,j,bi,bj) * sideDragU(i,j,bi,bj) + BV(i,j,bi,bj) = BV(i,j,bi,bj) + & + seaiceMaskV(i,j,bi,bj) * sideDragV(i,j,bi,bj) + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF +#endif /* SEAICE_ALLOW_SIDEDRAG */ + #ifndef OBCS_UVICE_OLD C-- prevent tri-diagonal solver from modifying OB values: DO bj=myByLo(myThid),myByHi(myThid) diff --git a/pkg/seaice/seaice_readparms.F b/pkg/seaice/seaice_readparms.F index 1ad1433595..8120dff973 100644 --- a/pkg/seaice/seaice_readparms.F +++ b/pkg/seaice/seaice_readparms.F @@ -163,6 +163,7 @@ SUBROUTINE SEAICE_READPARMS( myThid ) & SEAICE_drag_south, SEAICE_waterDrag_south, & SEAICE_dryIceAlb_south, SEAICE_wetIceAlb_south, & SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south, + & SEAICEsideDrag, uCoastLineFile, vCoastLineFile, & SEAICE_cBasalStar, SEAICEbasalDragU0, SEAICEbasalDragK1, & SEAICEbasalDragK2, SEAICE_wetAlbTemp, SEAICE_waterAlbedo, & SEAICE_strength, SEAICE_cStar, SEAICE_eccen, @@ -394,6 +395,7 @@ SUBROUTINE SEAICE_READPARMS( myThid ) SEAICE_drag = 0.001 _d 0 OCEAN_drag = 0.001 _d 0 SEAICE_waterDrag = 0.0055 _d 0 + SEAICEsideDrag = 0.0 _d 0 SEAICEdWatMin = 0.25 _d 0 SEAICE_dryIceAlb = 0.75 _d 0 SEAICE_wetIceAlb = 0.66 _d 0 @@ -509,6 +511,8 @@ SUBROUTINE SEAICE_READPARMS( myThid ) HeffFile = ' ' uIceFile = ' ' vIceFile = ' ' + uCoastLineFile = ' ' + vCoastLineFile = ' ' IMAX_TICE = 10 postSolvTempIter = 2 C LSR parameters @@ -1438,6 +1442,17 @@ SUBROUTINE SEAICE_READPARMS( myThid ) & SQUEEZE_RIGHT, myThid ) ENDIF +#ifdef SEAICE_ALLOW_SIDEDRAG + IF ( SEAICEsideDrag .NE. 0.0 _d 0 .AND. SEAICE_no_slip ) THEN + WRITE(msgBuf,'(3A)') '** WARNING ** SEAICE_READPARMS: ', + & 'SEAICEsideDrag .NE. 0.,', + & ' resetting SEAICE_no_slip = .FALSE.' + CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, + & SQUEEZE_RIGHT, myThid ) + SEAICE_no_slip = .FALSE. + ENDIF +#endif + IF ( DIFF1 .EQ. UNSET_RL ) THEN DIFF1 = 0. _d 0 chkFlag = .FALSE. diff --git a/pkg/seaice/seaice_sidedrag_stress.F b/pkg/seaice/seaice_sidedrag_stress.F new file mode 100644 index 0000000000..845ae9a099 --- /dev/null +++ b/pkg/seaice/seaice_sidedrag_stress.F @@ -0,0 +1,156 @@ +#include "SEAICE_OPTIONS.h" +#ifdef ALLOW_OBCS +# include "OBCS_OPTIONS.h" +#endif +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif +#undef SEAICE_ALTERNATIVE_MODULUS_U_AND_V + +CBOP +C !ROUTINE: SEAICE_SIDEDRAG +C !INTERFACE: + SUBROUTINE SEAICE_SIDEDRAG_STRESS( + I uIceArg, vIceArg, + I coastRoughUarg, coastRoughVarg, AREAarg, + O sideDragUarg,sideDragVarg, + I iStep, myTime, myIter, myThid ) + +C !DESCRIPTION: \bv +C *==========================================================* +C | SUBROUTINE SEAICE_BOTTOMDRAG_COEFFS +C | o Compute the side drag stress for ice-side drag, +C | as a parameterization for island-drag fastice +C | (Liu et al 2022) +C *==========================================================* +C | written by Yuqing Liu, Jul 2020 +C *==========================================================* +C \ev + +C !USES: + IMPLICIT NONE + +C === Global variables === +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "DYNVARS.h" +#include "SEAICE_SIZE.h" +#include "SEAICE_PARAMS.h" +#include "SEAICE.h" + +C !INPUT/OUTPUT PARAMETERS: +C === Routine arguments === +C u/vIceArg :: local copies of the current ice velocity +C coastRoughUarg :: local copy of normalized coast length +C coastRoughVarg :: local copy of normalized coast length +C SideDraXarg :: side drag in x direction +C SideDraYarg :: side drag in y direction +C iStep :: current sub-time step iterate +C myTime :: Simulation time +C myIter :: Simulation timestep number +C myThid :: my Thread Id. number + _RL uIceArg (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL vIceArg (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL AREAarg (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL sideDragUarg (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL sideDragVarg (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL coastRoughUarg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL coastRoughVarg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + + INTEGER iStep + _RL myTime + INTEGER myIter + INTEGER myThid + +#ifdef SEAICE_ALLOW_SIDEDRAG +C === local variables === +C i,j,bi,bj,ksrf :: loop indices + INTEGER i,j,bi,bj + _RL tmp, tmpx, tmpy + _RL uSpd(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL u_0 + + u_0 = 0.0005 _d 0 + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) +C initialize fields + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + sideDragUarg(i,j,bi,bj) = 0. _d 0 + sideDragVarg(i,j,bi,bj) = 0. _d 0 + uSpd(i,j) = 0. _d 0 + ENDDO + ENDDO + + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 + uSpd(i,j) = 0.5 _d 0 * + & SQRT((uIceArg(i, j,bi,bj)+uIceArg(i+1,j, bi,bj))**2 + & +(vIceArg(i, j,bi,bj)+vIceArg(i, j+1,bi,bj))**2) + ENDDO + ENDDO +C calculate the sidedrag coeff for landfast ice parameterisations +C (Liu et al 2022 use SEAICEsideDrag 2e-4) + IF (SEAICEsideDrag .GT. 0. _d 0 ) THEN + DO j=1-OLy+1,sNy+OLy-1 + DO i=1-OLx+1,sNx+OLx-1 + IF ( AREAarg(i,j,bi,bj) .GT. 0.01 _d 0 ) THEN +# ifndef SEAICE_ALTERNATIVE_MODULUS_U_AND_V + tmpx = 0.5 _d 0 * ( uSpd(i,j) + uSpd(i-1,j) ) + tmpy = 0.5 _d 0 * ( uSpd(i,j) + uSpd(i,j-1) ) +# else /* ifdef SEAICE_ALTERNATIVE_MODULUS_U_AND_V */ +C A different way to compute |u|, |v|, deprecated + tmpx = sqrt(uIceArg(i,j,bi,bj)**2+ (0.25 _d 0* + & (vIceArg(i, j, bi,bj)+vIceArg(i-1,j,bi,bj)+ + & vIceArg(i-1,j+1,bi,bj)+vIceArg(i, j+1,bi,bj)))**2) + tmpy = sqrt(vIceArg(i,j,bi,bj)**2+ (0.25 _d 0* + & (uIceArg(i, j, bi,bj)+uIceArg(i, j-1,bi,bj)+ + & uIceArg(i+1,j-1,bi,bj)+uIceArg(i+1,j,bi,bj)))**2) +#endif + + sideDragUarg(i,j,bi,bj) = + & SEAICEsideDrag/(tmpx + u_0) * coastRoughUarg(i,j,bi,bj) + & * seaiceMassU(i,j,bi,bj) + sideDragVarg(i,j,bi,bj) = + & SEAICEsideDrag/(tmpy + u_0) * coastRoughVarg(i,j,bi,bj) + & * seaiceMassV(i,j,bi,bj) + ENDIF + ENDDO + ENDDO + ELSEIF ( SEAICEsideDrag .LT. 0. _d 0 ) THEN +C this is the quadratic drag law version, deprecated + DO j=1-OLy+1,sNy+OLy-1 + DO i=1-OLx+1,sNx+OLx-1 + IF ( AREAarg(i,j,bi,bj) .GT. 0.01 _d 0 ) THEN +# ifndef SEAICE_ALTERNATIVE_MODULUS_U_AND_V + tmpx = 0.5 _d 0 * ( uSpd(i,j) + uSpd(i-1,j) ) + tmpy = 0.5 _d 0 * ( uSpd(i,j) + uSpd(i,j-1) ) +# else /* ifdef SEAICE_ALTERNATIVE_MODULUS_U_AND_V */ +C A different way to compute |u|, |v|, deprecated + tmpx = sqrt(uIceArg(i,j,bi,bj)**2+ (0.25 _d 0* + & (vIceArg(i, j, bi,bj)+vIceArg(i-1,j,bi,bj)+ + & vIceArg(i-1,j+1,bi,bj)+vIceArg(i, j+1,bi,bj)))**2) + tmpy = sqrt(vIceArg(i,j,bi,bj)**2+ (0.25 _d 0* + & (uIceArg(i, j, bi,bj)+uIceArg(i, j-1,bi,bj)+ + & uIceArg(i+1,j-1,bi,bj)+uIceArg(i+1,j,bi,bj)))**2) +# endif + + sideDragUarg(i,j,bi,bj) = + & - SEAICEsideDrag * tmpx * coastRoughUarg(i,j,bi,bj) + & * seaiceMassU(i,j,bi,bj) + sideDragVarg(i,j,bi,bj) = + & - SEAICEsideDrag * tmpy * coastRoughVarg(i,j,bi,bj) + & * seaiceMassV(i,j,bi,bj) + ENDIF + ENDDO + ENDDO + ENDIF + ENDDO + ENDDO + +#endif /* SEAICE_ALLOW_SIDEDRAG */ + + RETURN + END diff --git a/pkg/seaice/seaice_summary.F b/pkg/seaice/seaice_summary.F index ece2412b1b..850d1baba5 100644 --- a/pkg/seaice/seaice_summary.F +++ b/pkg/seaice/seaice_summary.F @@ -161,6 +161,14 @@ SUBROUTINE SEAICE_SUMMARY( myThid ) CALL WRITE_0D_RL( SEAICEbasalDragK2 ,INDEX_NONE, & 'SEAICEbasalDragK2 =', ' /* Basal drag parameter */') #endif /* SEAICE_ALLOW_BOTTOMDRAG */ +#ifdef SEAICE_ALLOW_SIDEDRAG + CALL WRITE_0D_RL( SEAICEsideDrag ,INDEX_NONE, + & 'SEAICEsideDrag =', ' /* lateral drag parameter */') + CALL WRITE_0D_C( uCoastLineFile, -1, INDEX_NONE, + & 'uCoastLineFile =', ' /* u-coastline length file */') + CALL WRITE_0D_C( vCoastLineFile, -1, INDEX_NONE, + & 'vCoastLineFile =', ' /* v-coastline length file */') +#endif /* SEAICE_ALLOW_SIDEDRAG */ CALL WRITE_0D_L ( SEAICEuseTilt, INDEX_NONE, & 'SEAICEuseTilt =', ' /* include surface tilt in dyna. */') CALL WRITE_0D_L ( SEAICEuseTEM, INDEX_NONE, diff --git a/verification/offline_exf_seaice/code_ad/MDSIO_BUFF_WH.h b/verification/offline_exf_seaice/code_ad/MDSIO_BUFF_WH.h index fff5e500ae..2bec7f3df7 100644 --- a/verification/offline_exf_seaice/code_ad/MDSIO_BUFF_WH.h +++ b/verification/offline_exf_seaice/code_ad/MDSIO_BUFF_WH.h @@ -36,7 +36,7 @@ C (during read) or copy to (during write). LOGICAL writeWh COMMON /MDS_WH_BUFFERS_3D_I/ iWh, jWh INTEGER nWh, iWh, jWh - PARAMETER (nWh=93) + PARAMETER (nWh=95) COMMON /MDS_WH_BUFFERS_3D_R8/ fld3d_procbuff_r8 # ifdef INCLUDE_WHIO_GLOBUFF_3D & , fld3d_globuff_r8 diff --git a/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h b/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h index be18d11b1d..74bc1c6aa7 100644 --- a/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h +++ b/verification/offline_exf_seaice/code_ad/SEAICE_OPTIONS.h @@ -189,6 +189,10 @@ C Use parameterisation of grounding ice for a better representation C of fastice in shallow seas # undef SEAICE_ALLOW_BOTTOMDRAG +C Use parameterisation of explicit lateral drag for a better +C representation of fastice along coast lines and islands +# define SEAICE_ALLOW_SIDEDRAG + #endif /* SEAICE_CGRID */ #ifdef SEAICE_BGRID_DYNAMICS diff --git a/verification/seaice_obcs/code/SEAICE_OPTIONS.h b/verification/seaice_obcs/code/SEAICE_OPTIONS.h index 3d2ae1779d..8e39f3b4f4 100644 --- a/verification/seaice_obcs/code/SEAICE_OPTIONS.h +++ b/verification/seaice_obcs/code/SEAICE_OPTIONS.h @@ -183,6 +183,10 @@ C Use parameterisation of grounding ice for a better representation C of fastice in shallow seas # undef SEAICE_ALLOW_BOTTOMDRAG +C Use parameterisation of explicit lateral drag for a better +C representation of fastice along coast lines and islands +# define SEAICE_ALLOW_SIDEDRAG + #endif /* SEAICE_CGRID */ #ifdef SEAICE_BGRID_DYNAMICS