sea ice comments by ML

[ifenty]

On 3. Dec 2019, at 23:54, Fenty, Ian G (US 329C) <[email protected]> wrote:
Martin,
I tried increasing SEAICEnonLinIterMax to 10 per your suggestion but seaice_check throws an error:
Need to increase MPSEUDOTIMESTEPS IN SEAICE_PARAMS.h
Is MPSEUDOTIMESTEPS some kind of upper bound on SEAICEnonLinIterMax?
Ian


Hi Ian,

didn’t think of that. For the adjoint (and only for the adjoint), one needs to hard-code a upper limit for the number of nonlinear steps (formerly called NPSEUDOTIMESTEPS) because it needs to be know when tapes are defined (in the_main_loop.F). That’s the parameter MPSEUDOTIMESTEPS. Just set it to 10. This will increase the memory footprint, but not too much as long as SEAICE_LSR_ADJOINT_ITER is undefined (which you want to do because you don’t even want to use the adjoint of the dynamics).

There are a couple of these upper bounds because of the adjoint (see SEAICE_SIZE.h)

Martin

# SEAICE parameters
# taking /MITgcm_contrib/ecco_utils/ecco_v4_release3_devel
 &SEAICE_PARM01

      SEAICEpresH0=2.,
      SEAICEpresPow0=1,
      SEAICEpresPow1=1,

      SEAICE_strength = 2.25e4,
# I would try to get rid of this in favor of SEAICE_multDim > 1
      SEAICE_area_max = 0.97,      

      SEAICE_no_slip     = .TRUE.,

# this is already the default (with recent code)
      SEAICE_drag=0.001,
      OCEAN_drag=0.001,

#####
      SEAICEuseTILT=.FALSE.,
# I suspect that this makes it necessary to have
# this weird SEAICE_area_max = 0.97
      SEAICE_multDim=1,
      SEAICErestoreUnderIce=.TRUE.,

      SEAICE_salt0=4.,

# I recommend 1.e-5 (only slightly more expensive)
      LSR_ERROR          = 2.e-4,
# I recomment (makes the dynamics solver 5 times more expensive than default 2)
# SEAICEnonLinIterMax=10
# this is already the default:
     SEAICEuseDYNAMICS  = .TRUE.,
#default = -50.
      MIN_ATEMP          = -40.,
#default = -50.
      MIN_TICE           = -40.,
      SEAICEadvScheme    = 30,
# this is already the default:
      SEAICEuseFluxForm = .TRUE.,
# this is already the default:
      SEAICEadvSnow      = .TRUE.,
# this is truely disturbing. Why do you do/need that?
      SEAICEdiffKhHeff   = 400.,
      SEAICEdiffKhArea   = 400.,
      SEAICEdiffKhSnow   = 400.,
# this is already the default:
      SEAICEuseFlooding  = .TRUE.,
# I don't know what these 5 parameter do exactly
      SEAICE_mcPheePiston= 3.858024691358025E-05,
      SEAICE_frazilFrac  = 1.,
      SEAICE_mcPheeTaper = 0.,
      SEAICE_areaLossFormula=2,
      SEAICEheatConsFix  = .TRUE.,
# why not use the default with variable freezing point?
      SEAICE_tempFrz0    = -1.96,
      SEAICE_dTempFrz_dS = 0.,
# this is already the default
      SEAICEuseMetricTerms = .TRUE.,
# why do you need that?
      SEAICE_clipVelocities = .TRUE.,
# there are new defaults:
# will lead to zero velocities where there is no ice (no need for clipping
# velocities anymore)
# SEAICEscaleSurfStress = .TRUE.,
# not important just more consistent
# SEAICEaddSnowMass = .TRUE.,
# this makes more sense than 30, but may not work with adjoint. If you
# don't use the adjoint, then I would use this default or 33
# SEAICEadvScheme = 77,
# only makes sense with SEAICE_multDim > 1, but then it makes a lot of sense
# SEAICE_useMultDimSnow = .TRUE.
# the following are for better stability of the solver
# SEAICE_OLx/y = OLx/y - 2
# SEAICEetaZmethod = 3
#end of new defaults

#take 33% out of (1-albedo)
      SEAICE_dryIceAlb   = 0.84,
      SEAICE_wetIceAlb   = 0.78,
      SEAICE_drySnowAlb  = 0.90,
      SEAICE_wetSnowAlb  = 0.8 ,
#default albedos
      SEAICE_dryIceAlb_south   = 0.75
      SEAICE_wetIceAlb_south   = 0.66
      SEAICE_drySnowAlb_south  = 0.84
      SEAICE_wetSnowAlb_south  = 0.7 
 /
#
 &SEAICE_PARM02
 /

CPP -flags:

C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/ecco_utils/ecco_v4_release3_devel/code/SEAICE_OPTIONS.h,v 1.1 2017/05/04 17:46:37 ou.wang Exp $
C $Name:  $

C     *==========================================================*
C     | SEAICE_OPTIONS.h
C     | o CPP options file for sea ice package.
C     *==========================================================*
C     | Use this file for selecting options within the sea ice
C     | package.
C     *==========================================================*

#ifndef SEAICE_OPTIONS_H
#define SEAICE_OPTIONS_H
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

#ifdef ALLOW_SEAICE
C     Package-specific Options & Macros go here

C--   Write "text-plots" of certain fields in STDOUT for debugging.
#undef SEAICE_DEBUG

C--   Allow sea-ice dynamic code.
C     This option is provided to allow use of TAMC
C     on the thermodynamics component of the code only.
C     Sea-ice dynamics can also be turned off at runtime
C     using variable SEAICEuseDYNAMICS.
#define SEAICE_ALLOW_DYNAMICS

C--   By default, the sea-ice package uses its own integrated bulk
C     formulae to compute fluxes (fu, fv, EmPmR, Qnet, and Qsw) over
C     open-ocean.  When this flag is set, these variables are computed
C     in a separate external package, for example, pkg/exf, and then
C     modified for sea-ice effects by pkg/seaice.
#define SEAICE_EXTERNAL_FLUXES

C--   This CPP flag has been retired.  The number of ice categories
C     used to solve for seaice flux is now specified by run-time
C     parameter SEAICE_multDim.
C     Note: be aware of pickup_seaice.* compatibility issues when
C     restarting a simulation with a different number of categories.
c#define SEAICE_MULTICATEGORY

C--   run with sea Ice Thickness Distribution (ITD);
C     set number of categories (nITD) in SEAICE_SIZE.h
#undef SEAICE_ITD

C--   Since the missing sublimation term is now included
C     this flag is needed for backward compatibility
#undef SEAICE_DISABLE_SUBLIM

C--   Suspected missing term in coupled ocn-ice heat budget (to be confirmed)
#undef SEAICE_DISABLE_HEATCONSFIX

C--   Default is constant seaice salinity (SEAICE_salt0); Define the following
C     flag to consider (space & time) variable salinity: advected and forming
C     seaice with a fraction (=SEAICE_saltFrac) of freezing seawater salinity.
C- Note: SItracer also offers an alternative way to handle variable salinity.
#undef SEAICE_VARIABLE_SALINITY

C--   Tracers of ice and/or ice cover.
#undef ALLOW_SITRACER
#ifdef ALLOW_SITRACER
C--   To try avoid 'spontaneous generation' of tracer maxima by advdiff.
# define ALLOW_SITRACER_ADVCAP
#endif

C--   Enable grease ice parameterization
C     The grease ice parameterization delays formation of solid
C     sea ice from frazil ice by a time constant and provides a
C     dynamic calculation of the initial solid sea ice thickness
C     HO as a function of winds, currents and available grease ice
C     volume. Grease ice does not significantly reduce heat loss
C     from the ocean in winter and area covered by grease is thus
C     handled like open water.
C     (For details see Smedsrud and Martin, 2014, Ann.Glac.)
C     Set SItrName(1) = 'grease' in namelist SEAICE_PARM03 in data.seaice
C     then output SItr01 is SItrNameLong(1) = 'grease ice volume fraction',
C     with SItrUnit(1) = '[0-1]', which needs to be multiplied by SIheff
C     to yield grease ice volume. Additionally, the actual grease ice
C     layer thickness (diagnostic SIgrsLT) can be saved.
#undef SEAICE_GREASE
C--   grease ice uses SItracer:
#ifdef SEAICE_GREASE
# define ALLOW_SITRACER
# define ALLOW_SITRACER_ADVCAP
#endif

C--   By default the seaice model is discretized on a B-Grid (for
C     historical reasons). Define the following flag to use a new
C     (not thoroughly) test version on a C-grid
#define SEAICE_CGRID

C--   Only for the C-grid version it is possible to
#ifdef SEAICE_CGRID
C     enable JFNK code by defining the following flag
# undef  SEAICE_ALLOW_JFNK
C     enable LSR to use global (multi-tile) tri-diagonal solver
# undef SEAICE_GLOBAL_3DIAG_SOLVER
C     enable EVP code by defining the following flag
CML this can be undefined since you dont use it
# define SEAICE_ALLOW_EVP
# ifdef SEAICE_ALLOW_EVP
C--   When set use SEAICE_zetaMin and SEAICE_evpDampC to limit viscosities
C     from below and above in seaice_evp: not necessary, and not recommended
#  undef SEAICE_ALLOW_CLIPZETA
# endif /* SEAICE_ALLOW_EVP */
C     regularize zeta to zmax with a smooth tanh-function instead
C     of a min(zeta,zmax). This improves convergence of iterative
C     solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP
# undef SEAICE_ZETA_SMOOTHREG
C     allow the truncated ellipse rheology (runtime flag SEAICEuseTEM)
# undef SEAICE_ALLOW_TEM
C     Use LSR vector code; not useful on non-vector machines, because it
C     slows down convergence considerably, but the extra iterations are
C     more than made up by the much faster code on vector machines. For
C     the only regularly test vector machine these flags a specified
C     in the build options file SUPER-UX_SX-8_sxf90_awi, so that we comment
C     them out here.
C# define SEAICE_VECTORIZE_LSR
C# ifdef SEAICE_VECTORIZE_LSR
C     Use modified LSR vector code that splits vector loop into two with
C     step size 2. This modification improves the convergence of the vector
C     code dramatically, so that is may actually be useful in general, but
C     that needs to be tested.
C#  define SEAICE_VECTORIZE_LSR_ZEBRA
CML replacing this with new flag:
# undef SEAICE_VECTORIZE_LSR
C     Use zebra-method (alternate lines) for line-successive-relaxation
C     This modification improves the convergence of the vector code
C     dramatically, so that is may actually be useful in general, but
C     that needs to be tested. Can be used without vectorization options.
# undef SEAICE_LSR_ZEBRA
C     Use parameterisation of grounding ice for a better representation
C     of fastice in shallow seas
C# endif
#else /* not SEAICE_CGRID, but old B-grid */
C--   By default for B-grid dynamics solver wind stress under sea-ice is
C     set to the same value as it would be if there was no sea-ice.
C     Define following CPP flag for B-grid ice-ocean stress coupling.
# define SEAICE_BICE_STRESS

C--   By default for B-grid dynamics solver surface tilt is obtained
C     indirectly via geostrophic velocities. Define following CPP
C     in order to use ETAN instead.
# define EXPLICIT_SSH_SLOPE
C--   Defining this flag turns on FV-discretization of the B-grid LSOR solver.
C     It is smoother and includes all metric terms, similar to C-grid solvers.
C     It is here for completeness, but its usefulness is unclear.
# undef SEAICE_LSRBNEW
#endif /* SEAICE_CGRID */

C--   When set limit the Ice-Loading to mass of 1/5 of Surface ocean grid-box
#undef SEAICE_CAP_ICELOAD
C--   When set use SEAICE_clipVelocties = .true., to clip U/VICE at 40cm/s,
C     not recommended
#define SEAICE_ALLOW_CLIPVELS
C--   When set cap the sublimation latent heat flux in solve4temp according
C     to the available amount of ice+snow. Otherwise this term is treated
C     like all of the others -- residuals heat and fw stocks are passed to
C     the ocean at the end of seaice_growth in a conservative manner.
C     SEAICE_CAP_SUBLIM is not needed as of now, but kept just in case.
#undef SEAICE_CAP_SUBLIM

C--   Enable free drift code
#define SEAICE_ALLOW_FREEDRIFT

C--   pkg/seaice cost functions compile flags
c       >>> Sea-ice volume (requires pkg/cost)
#undef ALLOW_COST_ICE
c       >>> Sea-ice misfit to obs (requires pkg/cost and ecco)
#undef ALLOW_SEAICE_COST_SMR_AREA

C--   enforce cfl condition without cuting sensitivity flow
c#define ALLOW_CFL_FIX

C--   cut the adjoint dependency to hactual, etc.
c# undef SEAICE_SIMPLIFY_GROWTH_ADJ

C--   go through heff and open ocean
c#define SEAICE_MODIFY_GROWTH_ADJ

#endif /* ALLOW_SEAICE */
#endif /* SEAICE_OPTIONS_H */

CEH3 ;;; Local Variables: ***
CEH3 ;;; mode:fortran ***
CEH3 ;;; End: ***

[owang01]

In MITgcm/MITgcm@dc26f15, the default SEAICEnonLinIterMax is still 2 since we use SEAICEuseLSR=.TRUE. See lines 839-844 of seaice_readparms.F

https://github.com/MITgcm/MITgcm/blob/f397060a0550404d6ba24560e5b359909fc19356/pkg/seaice/seaice_readparms.F#L839-L844