Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 26098) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 26099) @@ -10,451 +10,835 @@ [ +dnl with-HydrologyTws{{{ + +AC_ARG_WITH([HydrologyTws], + + AS_HELP_STRING([--with-HydrologyTws = YES], [compile with HydrologyTws capabilities (default is yes)]), + + [HYDROLOGYTWS=$withval],[HYDROLOGYTWS=yes]) + +AC_MSG_CHECKING(for HydrologyTws capability compilation) + + +HAVE_HYDROLOGYTWS=no + +if test "x$HYDROLOGYTWS" = "xyes"; then + + HAVE_HYDROLOGYTWS=yes + + AC_DEFINE([_HAVE_HYDROLOGYTWS_],[1],[with HydrologyTws capability]) + +fi + +AM_CONDITIONAL([HYDROLOGYTWS], [test x$HAVE_HYDROLOGYTWS = xyes]) + +AC_MSG_RESULT($HAVE_HYDROLOGYTWS) + +dnl }}} dnl with-AdjointBalancethickness2{{{ + AC_ARG_WITH([AdjointBalancethickness2], + AS_HELP_STRING([--with-AdjointBalancethickness2 = YES], [compile with AdjointBalancethickness2 capabilities (default is yes)]), + [ADJOINTBALANCETHICKNESS2=$withval],[ADJOINTBALANCETHICKNESS2=yes]) + AC_MSG_CHECKING(for AdjointBalancethickness2 capability compilation) + HAVE_ADJOINTBALANCETHICKNESS2=no + if test "x$ADJOINTBALANCETHICKNESS2" = "xyes"; then + HAVE_ADJOINTBALANCETHICKNESS2=yes + AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2 capability]) -fi + +fi + AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes]) + AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS2) + dnl }}} dnl with-AdjointBalancethickness{{{ + AC_ARG_WITH([AdjointBalancethickness], + AS_HELP_STRING([--with-AdjointBalancethickness = YES], [compile with AdjointBalancethickness capabilities (default is yes)]), + [ADJOINTBALANCETHICKNESS=$withval],[ADJOINTBALANCETHICKNESS=yes]) + AC_MSG_CHECKING(for AdjointBalancethickness capability compilation) + HAVE_ADJOINTBALANCETHICKNESS=no + if test "x$ADJOINTBALANCETHICKNESS" = "xyes"; then + HAVE_ADJOINTBALANCETHICKNESS=yes + AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethickness capability]) -fi + +fi + AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes]) + AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS) + dnl }}} dnl with-AdjointHoriz{{{ + AC_ARG_WITH([AdjointHoriz], + AS_HELP_STRING([--with-AdjointHoriz = YES], [compile with AdjointHoriz capabilities (default is yes)]), + [ADJOINTHORIZ=$withval],[ADJOINTHORIZ=yes]) + AC_MSG_CHECKING(for AdjointHoriz capability compilation) + HAVE_ADJOINTHORIZ=no + if test "x$ADJOINTHORIZ" = "xyes"; then + HAVE_ADJOINTHORIZ=yes + AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHoriz capability]) -fi + +fi + AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes]) + AC_MSG_RESULT($HAVE_ADJOINTHORIZ) + dnl }}} dnl with-Balancethickness2{{{ + AC_ARG_WITH([Balancethickness2], + AS_HELP_STRING([--with-Balancethickness2 = YES], [compile with Balancethickness2 capabilities (default is yes)]), + [BALANCETHICKNESS2=$withval],[BALANCETHICKNESS2=yes]) + AC_MSG_CHECKING(for Balancethickness2 capability compilation) + HAVE_BALANCETHICKNESS2=no + if test "x$BALANCETHICKNESS2" = "xyes"; then + HAVE_BALANCETHICKNESS2=yes + AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2 capability]) -fi + +fi + AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes]) + AC_MSG_RESULT($HAVE_BALANCETHICKNESS2) + dnl }}} dnl with-Balancethickness{{{ + AC_ARG_WITH([Balancethickness], + AS_HELP_STRING([--with-Balancethickness = YES], [compile with Balancethickness capabilities (default is yes)]), + [BALANCETHICKNESS=$withval],[BALANCETHICKNESS=yes]) + AC_MSG_CHECKING(for Balancethickness capability compilation) + HAVE_BALANCETHICKNESS=no + if test "x$BALANCETHICKNESS" = "xyes"; then + HAVE_BALANCETHICKNESS=yes + AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethickness capability]) -fi + +fi + AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes]) + AC_MSG_RESULT($HAVE_BALANCETHICKNESS) + dnl }}} dnl with-BalancethicknessSoft{{{ + AC_ARG_WITH([BalancethicknessSoft], + AS_HELP_STRING([--with-BalancethicknessSoft = YES], [compile with BalancethicknessSoft capabilities (default is yes)]), + [BALANCETHICKNESSSOFT=$withval],[BALANCETHICKNESSSOFT=yes]) + AC_MSG_CHECKING(for BalancethicknessSoft capability compilation) + HAVE_BALANCETHICKNESSSOFT=no + if test "x$BALANCETHICKNESSSOFT" = "xyes"; then + HAVE_BALANCETHICKNESSSOFT=yes + AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoft capability]) -fi + +fi + AM_CONDITIONAL([BALANCETHICKNESSSOFT], [test x$HAVE_BALANCETHICKNESSSOFT = xyes]) + AC_MSG_RESULT($HAVE_BALANCETHICKNESSSOFT) + dnl }}} dnl with-Balancevelocity{{{ + AC_ARG_WITH([Balancevelocity], + AS_HELP_STRING([--with-Balancevelocity = YES], [compile with Balancevelocity capabilities (default is yes)]), + [BALANCEVELOCITY=$withval],[BALANCEVELOCITY=yes]) + AC_MSG_CHECKING(for Balancevelocity capability compilation) + HAVE_BALANCEVELOCITY=no + if test "x$BALANCEVELOCITY" = "xyes"; then + HAVE_BALANCEVELOCITY=yes + AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocity capability]) -fi + +fi + AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes]) + AC_MSG_RESULT($HAVE_BALANCEVELOCITY) + dnl }}} dnl with-DamageEvolution{{{ + AC_ARG_WITH([DamageEvolution], + AS_HELP_STRING([--with-DamageEvolution = YES], [compile with DamageEvolution capabilities (default is yes)]), + [DAMAGEEVOLUTION=$withval],[DAMAGEEVOLUTION=yes]) + AC_MSG_CHECKING(for DamageEvolution capability compilation) + HAVE_DAMAGEEVOLUTION=no + if test "x$DAMAGEEVOLUTION" = "xyes"; then + HAVE_DAMAGEEVOLUTION=yes + AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolution capability]) -fi + +fi + AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes]) + AC_MSG_RESULT($HAVE_DAMAGEEVOLUTION) + dnl }}} dnl with-DepthAverage{{{ + AC_ARG_WITH([DepthAverage], + AS_HELP_STRING([--with-DepthAverage = YES], [compile with DepthAverage capabilities (default is yes)]), + [DEPTHAVERAGE=$withval],[DEPTHAVERAGE=yes]) + AC_MSG_CHECKING(for DepthAverage capability compilation) + HAVE_DEPTHAVERAGE=no + if test "x$DEPTHAVERAGE" = "xyes"; then + HAVE_DEPTHAVERAGE=yes + AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAverage capability]) -fi + +fi + AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes]) + AC_MSG_RESULT($HAVE_DEPTHAVERAGE) + dnl }}} dnl with-Enthalpy{{{ + AC_ARG_WITH([Enthalpy], + AS_HELP_STRING([--with-Enthalpy = YES], [compile with Enthalpy capabilities (default is yes)]), + [ENTHALPY=$withval],[ENTHALPY=yes]) + AC_MSG_CHECKING(for Enthalpy capability compilation) + HAVE_ENTHALPY=no + if test "x$ENTHALPY" = "xyes"; then + HAVE_ENTHALPY=yes + AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpy capability]) -fi + +fi + AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes]) + AC_MSG_RESULT($HAVE_ENTHALPY) + dnl }}} dnl with-Esa{{{ + AC_ARG_WITH([Esa], + AS_HELP_STRING([--with-Esa = YES], [compile with Esa capabilities (default is yes)]), + [ESA=$withval],[ESA=yes]) + AC_MSG_CHECKING(for Esa capability compilation) + HAVE_ESA=no + if test "x$ESA" = "xyes"; then + HAVE_ESA=yes + AC_DEFINE([_HAVE_ESA_],[1],[with Esa capability]) -fi + +fi + AM_CONDITIONAL([ESA], [test x$HAVE_ESA = xyes]) + AC_MSG_RESULT($HAVE_ESA) + dnl }}} dnl with-Extrapolation{{{ + AC_ARG_WITH([Extrapolation], + AS_HELP_STRING([--with-Extrapolation = YES], [compile with Extrapolation capabilities (default is yes)]), + [EXTRAPOLATION=$withval],[EXTRAPOLATION=yes]) + AC_MSG_CHECKING(for Extrapolation capability compilation) + HAVE_EXTRAPOLATION=no + if test "x$EXTRAPOLATION" = "xyes"; then + HAVE_EXTRAPOLATION=yes + AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolation capability]) -fi + +fi + AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes]) + AC_MSG_RESULT($HAVE_EXTRAPOLATION) + dnl }}} dnl with-ExtrudeFromBase{{{ + AC_ARG_WITH([ExtrudeFromBase], + AS_HELP_STRING([--with-ExtrudeFromBase = YES], [compile with ExtrudeFromBase capabilities (default is yes)]), + [EXTRUDEFROMBASE=$withval],[EXTRUDEFROMBASE=yes]) + AC_MSG_CHECKING(for ExtrudeFromBase capability compilation) + HAVE_EXTRUDEFROMBASE=no + if test "x$EXTRUDEFROMBASE" = "xyes"; then + HAVE_EXTRUDEFROMBASE=yes + AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBase capability]) -fi + +fi + AM_CONDITIONAL([EXTRUDEFROMBASE], [test x$HAVE_EXTRUDEFROMBASE = xyes]) + AC_MSG_RESULT($HAVE_EXTRUDEFROMBASE) + dnl }}} dnl with-ExtrudeFromTop{{{ + AC_ARG_WITH([ExtrudeFromTop], + AS_HELP_STRING([--with-ExtrudeFromTop = YES], [compile with ExtrudeFromTop capabilities (default is yes)]), + [EXTRUDEFROMTOP=$withval],[EXTRUDEFROMTOP=yes]) + AC_MSG_CHECKING(for ExtrudeFromTop capability compilation) + HAVE_EXTRUDEFROMTOP=no + if test "x$EXTRUDEFROMTOP" = "xyes"; then + HAVE_EXTRUDEFROMTOP=yes + AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTop capability]) -fi + +fi + AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes]) + AC_MSG_RESULT($HAVE_EXTRUDEFROMTOP) + dnl }}} dnl with-FreeSurfaceBase{{{ + AC_ARG_WITH([FreeSurfaceBase], + AS_HELP_STRING([--with-FreeSurfaceBase = YES], [compile with FreeSurfaceBase capabilities (default is yes)]), + [FREESURFACEBASE=$withval],[FREESURFACEBASE=yes]) + AC_MSG_CHECKING(for FreeSurfaceBase capability compilation) + HAVE_FREESURFACEBASE=no + if test "x$FREESURFACEBASE" = "xyes"; then + HAVE_FREESURFACEBASE=yes + AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBase capability]) -fi + +fi + AM_CONDITIONAL([FREESURFACEBASE], [test x$HAVE_FREESURFACEBASE = xyes]) + AC_MSG_RESULT($HAVE_FREESURFACEBASE) + dnl }}} dnl with-FreeSurfaceTop{{{ + AC_ARG_WITH([FreeSurfaceTop], + AS_HELP_STRING([--with-FreeSurfaceTop = YES], [compile with FreeSurfaceTop capabilities (default is yes)]), + [FREESURFACETOP=$withval],[FREESURFACETOP=yes]) + AC_MSG_CHECKING(for FreeSurfaceTop capability compilation) + HAVE_FREESURFACETOP=no + if test "x$FREESURFACETOP" = "xyes"; then + HAVE_FREESURFACETOP=yes + AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTop capability]) -fi + +fi + AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes]) + AC_MSG_RESULT($HAVE_FREESURFACETOP) + dnl }}} dnl with-GLheightadvection{{{ + AC_ARG_WITH([GLheightadvection], + AS_HELP_STRING([--with-GLheightadvection = YES], [compile with GLheightadvection capabilities (default is yes)]), + [GLHEIGHTADVECTION=$withval],[GLHEIGHTADVECTION=yes]) + AC_MSG_CHECKING(for GLheightadvection capability compilation) + HAVE_GLHEIGHTADVECTION=no + if test "x$GLHEIGHTADVECTION" = "xyes"; then + HAVE_GLHEIGHTADVECTION=yes + AC_DEFINE([_HAVE_GLHEIGHTADVECTION_],[1],[with GLheightadvection capability]) -fi + +fi + AM_CONDITIONAL([GLHEIGHTADVECTION], [test x$HAVE_GLHEIGHTADVECTION = xyes]) + AC_MSG_RESULT($HAVE_GLHEIGHTADVECTION) -dnl }}} -dnl with-Gia{{{ -AC_ARG_WITH([Gia], - AS_HELP_STRING([--with-Gia= YES], [compile with Gia capabilities (default is yes)]), - [GIA=$withval],[GIA=yes]) -AC_MSG_CHECKING(for Gia capability compilation) - -HAVE_GIA=no -if test "x$GIA" = "xyes"; then - HAVE_GIA=yes - AC_DEFINE([_HAVE_GIA_],[1],[with Giacapability]) -fi -AM_CONDITIONAL([GIA], [test x$HAVE_GIA = xyes]) -AC_MSG_RESULT($HAVE_GIA) + dnl }}} dnl with-HydrologyDCEfficient{{{ + AC_ARG_WITH([HydrologyDCEfficient], + AS_HELP_STRING([--with-HydrologyDCEfficient = YES], [compile with HydrologyDCEfficient capabilities (default is yes)]), + [HYDROLOGYDCEFFICIENT=$withval],[HYDROLOGYDCEFFICIENT=yes]) + AC_MSG_CHECKING(for HydrologyDCEfficient capability compilation) + HAVE_HYDROLOGYDCEFFICIENT=no + if test "x$HYDROLOGYDCEFFICIENT" = "xyes"; then + HAVE_HYDROLOGYDCEFFICIENT=yes + AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficient capability]) -fi + +fi + AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes]) + AC_MSG_RESULT($HAVE_HYDROLOGYDCEFFICIENT) + dnl }}} dnl with-HydrologyDCInefficient{{{ + AC_ARG_WITH([HydrologyDCInefficient], + AS_HELP_STRING([--with-HydrologyDCInefficient = YES], [compile with HydrologyDCInefficient capabilities (default is yes)]), + [HYDROLOGYDCINEFFICIENT=$withval],[HYDROLOGYDCINEFFICIENT=yes]) + AC_MSG_CHECKING(for HydrologyDCInefficient capability compilation) + HAVE_HYDROLOGYDCINEFFICIENT=no + if test "x$HYDROLOGYDCINEFFICIENT" = "xyes"; then + HAVE_HYDROLOGYDCINEFFICIENT=yes + AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficient capability]) -fi + +fi + AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes]) + AC_MSG_RESULT($HAVE_HYDROLOGYDCINEFFICIENT) + +dnl }}} +dnl with-HydrologyGlaDS{{{ + +AC_ARG_WITH([HydrologyGlaDS], + + AS_HELP_STRING([--with-HydrologyGlaDS = YES], [compile with HydrologyGlaDS capabilities (default is yes)]), + + [HYDROLOGYGLADS=$withval],[HYDROLOGYGLADS=yes]) + +AC_MSG_CHECKING(for HydrologyGlaDS capability compilation) + + +HAVE_HYDROLOGYGLADS=no + +if test "x$HYDROLOGYGLADS" = "xyes"; then + + HAVE_HYDROLOGYGLADS=yes + + AC_DEFINE([_HAVE_HYDROLOGYGLADS_],[1],[with HydrologyGlaDS capability]) + +fi + +AM_CONDITIONAL([HYDROLOGYGLADS], [test x$HAVE_HYDROLOGYGLADS = xyes]) + +AC_MSG_RESULT($HAVE_HYDROLOGYGLADS) + +dnl }}} +dnl with-HydrologyPism{{{ + +AC_ARG_WITH([HydrologyPism], + + AS_HELP_STRING([--with-HydrologyPism = YES], [compile with HydrologyPism capabilities (default is yes)]), + + [HYDROLOGYPISM=$withval],[HYDROLOGYPISM=yes]) + +AC_MSG_CHECKING(for HydrologyPism capability compilation) + + +HAVE_HYDROLOGYPISM=no + +if test "x$HYDROLOGYPISM" = "xyes"; then + + HAVE_HYDROLOGYPISM=yes + + AC_DEFINE([_HAVE_HYDROLOGYPISM_],[1],[with HydrologyPism capability]) + +fi + +AM_CONDITIONAL([HYDROLOGYPISM], [test x$HAVE_HYDROLOGYPISM = xyes]) + +AC_MSG_RESULT($HAVE_HYDROLOGYPISM) + +dnl }}} +dnl with-HydrologyShakti{{{ + +AC_ARG_WITH([HydrologyShakti], + + AS_HELP_STRING([--with-HydrologyShakti = YES], [compile with HydrologyShakti capabilities (default is yes)]), + + [HYDROLOGYSHAKTI=$withval],[HYDROLOGYSHAKTI=yes]) + +AC_MSG_CHECKING(for HydrologyShakti capability compilation) + + +HAVE_HYDROLOGYSHAKTI=no + +if test "x$HYDROLOGYSHAKTI" = "xyes"; then + + HAVE_HYDROLOGYSHAKTI=yes + + AC_DEFINE([_HAVE_HYDROLOGYSHAKTI_],[1],[with HydrologyShakti capability]) + +fi + +AM_CONDITIONAL([HYDROLOGYSHAKTI], [test x$HAVE_HYDROLOGYSHAKTI = xyes]) + +AC_MSG_RESULT($HAVE_HYDROLOGYSHAKTI) + dnl }}} dnl with-HydrologyShreve{{{ + AC_ARG_WITH([HydrologyShreve], + AS_HELP_STRING([--with-HydrologyShreve = YES], [compile with HydrologyShreve capabilities (default is yes)]), + [HYDROLOGYSHREVE=$withval],[HYDROLOGYSHREVE=yes]) + AC_MSG_CHECKING(for HydrologyShreve capability compilation) + HAVE_HYDROLOGYSHREVE=no + if test "x$HYDROLOGYSHREVE" = "xyes"; then + HAVE_HYDROLOGYSHREVE=yes + AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShreve capability]) -fi + +fi + AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes]) + AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE) -dnl }}} -dnl with-HydrologyTws{{{ -AC_ARG_WITH([HydrologyTws], - AS_HELP_STRING([--with-HydrologyTws = YES], [compile with HydrologyTws capabilities (default is yes)]), - [HYDROLOGYTWS=$withval],[HYDROLOGYTWS=yes]) -AC_MSG_CHECKING(for HydrologyTws capability compilation) - -HAVE_HYDROLOGYTWS=no -if test "x$HYDROLOGYTWS" = "xyes"; then - HAVE_HYDROLOGYTWS=yes - AC_DEFINE([_HAVE_HYDROLOGYTWS_],[1],[with HydrologyTws capability]) -fi -AM_CONDITIONAL([HYDROLOGYTWS], [test x$HAVE_HYDROLOGYTWS = xyes]) -AC_MSG_RESULT($HAVE_HYDROLOGYTWS) -dnl }}} -dnl with-HydrologyGlaDS{{{ -AC_ARG_WITH([HydrologyGlaDS], - AS_HELP_STRING([--with-HydrologyGlaDS = YES], [compile with HydrologyGlaDS capabilities (default is yes)]), - [HYDROLOGYGLADS=$withval],[HYDROLOGYGLADS=yes]) -AC_MSG_CHECKING(for HydrologyGlaDS capability compilation) - -HAVE_HYDROLOGYGLADS=no -if test "x$HYDROLOGYGLADS" = "xyes"; then - HAVE_HYDROLOGYGLADS=yes - AC_DEFINE([_HAVE_HYDROLOGYGLADS_],[1],[with HydrologyGlaDS capability]) -fi -AM_CONDITIONAL([HYDROLOGYGLADS], [test x$HAVE_HYDROLOGYGLADS = xyes]) -AC_MSG_RESULT($HAVE_HYDROLOGYGLADS) -dnl }}} -dnl with-HydrologyShakti{{{ -AC_ARG_WITH([HydrologyShakti], - AS_HELP_STRING([--with-HydrologyShakti = YES], [compile with HydrologyShakti capabilities (default is yes)]), - [HYDROLOGYSHAKTI=$withval],[HYDROLOGYSHAKTI=yes]) -AC_MSG_CHECKING(for HydrologyShakti capability compilation) - -HAVE_HYDROLOGYSHAKTI=no -if test "x$HYDROLOGYSHAKTI" = "xyes"; then - HAVE_HYDROLOGYSHAKTI=yes - AC_DEFINE([_HAVE_HYDROLOGYSHAKTI_],[1],[with HydrologyShakti capability]) -fi -AM_CONDITIONAL([HYDROLOGYSHAKTI], [test x$HAVE_HYDROLOGYSHAKTI = xyes]) -AC_MSG_RESULT($HAVE_HYDROLOGYSHAKTI) -dnl }}} -dnl with-HydrologyPism{{{ -AC_ARG_WITH([HydrologyPism], - AS_HELP_STRING([--with-HydrologyPism = YES], [compile with HydrologyPism capabilities (default is yes)]), - [HYDROLOGYPISM=$withval],[HYDROLOGYPISM=yes]) -AC_MSG_CHECKING(for HydrologyPism capability compilation) - -HAVE_HYDROLOGYPISM=no -if test "x$HYDROLOGYPISM" = "xyes"; then - HAVE_HYDROLOGYPISM=yes - AC_DEFINE([_HAVE_HYDROLOGYPISM_],[1],[with HydrologyPism capability]) -fi -AM_CONDITIONAL([HYDROLOGYPISM], [test x$HAVE_HYDROLOGYPISM = xyes]) -AC_MSG_RESULT($HAVE_HYDROLOGYPISM) + dnl }}} dnl with-L2ProjectionBase{{{ + AC_ARG_WITH([L2ProjectionBase], + AS_HELP_STRING([--with-L2ProjectionBase = YES], [compile with L2ProjectionBase capabilities (default is yes)]), + [L2PROJECTIONBASE=$withval],[L2PROJECTIONBASE=yes]) + AC_MSG_CHECKING(for L2ProjectionBase capability compilation) + HAVE_L2PROJECTIONBASE=no + if test "x$L2PROJECTIONBASE" = "xyes"; then + HAVE_L2PROJECTIONBASE=yes + AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBase capability]) -fi + +fi + AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes]) + AC_MSG_RESULT($HAVE_L2PROJECTIONBASE) + dnl }}} dnl with-L2ProjectionEPL{{{ + AC_ARG_WITH([L2ProjectionEPL], + AS_HELP_STRING([--with-L2ProjectionEPL = YES], [compile with L2ProjectionEPL capabilities (default is yes)]), + [L2PROJECTIONEPL=$withval],[L2PROJECTIONEPL=yes]) + AC_MSG_CHECKING(for L2ProjectionEPL capability compilation) + HAVE_L2PROJECTIONEPL=no + if test "x$L2PROJECTIONEPL" = "xyes"; then + HAVE_L2PROJECTIONEPL=yes + AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPL capability]) -fi + +fi + AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes]) + AC_MSG_RESULT($HAVE_L2PROJECTIONEPL) + dnl }}} dnl with-Levelset{{{ + AC_ARG_WITH([Levelset], + AS_HELP_STRING([--with-Levelset = YES], [compile with Levelset capabilities (default is yes)]), + [LEVELSET=$withval],[LEVELSET=yes]) + AC_MSG_CHECKING(for Levelset capability compilation) + HAVE_LEVELSET=no + if test "x$LEVELSET" = "xyes"; then + HAVE_LEVELSET=yes + AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelset capability]) -fi + +fi + AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes]) + AC_MSG_RESULT($HAVE_LEVELSET) + dnl }}} dnl with-Love{{{ + AC_ARG_WITH([Love], + AS_HELP_STRING([--with-Love = YES], [compile with Love capabilities (default is yes)]), + [LOVE=$withval],[LOVE=yes]) + AC_MSG_CHECKING(for Love capability compilation) + HAVE_LOVE=no + if test "x$LOVE" = "xyes"; then + HAVE_LOVE=yes + AC_DEFINE([_HAVE_LOVE_],[1],[with Love capability]) -fi + +fi + AM_CONDITIONAL([LOVE], [test x$HAVE_LOVE = xyes]) + AC_MSG_RESULT($HAVE_LOVE) + dnl }}} dnl with-Masstransport{{{ + AC_ARG_WITH([Masstransport], + AS_HELP_STRING([--with-Masstransport = YES], [compile with Masstransport capabilities (default is yes)]), + [MASSTRANSPORT=$withval],[MASSTRANSPORT=yes]) + AC_MSG_CHECKING(for Masstransport capability compilation) + HAVE_MASSTRANSPORT=no + if test "x$MASSTRANSPORT" = "xyes"; then + HAVE_MASSTRANSPORT=yes + AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransport capability]) -fi + +fi + AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes]) + AC_MSG_RESULT($HAVE_MASSTRANSPORT) -dnl }}} -dnl with-Oceanmasstransport{{{ -AC_ARG_WITH([Oceanmasstransport], - AS_HELP_STRING([--with-Oceanmasstransport = YES], [compile with Oceanmasstransport capabilities (default is yes)]), - [OCEANMASSTRANSPORT=$withval],[OCEANMASSTRANSPORT=yes]) -AC_MSG_CHECKING(for Oceanmasstransport capability compilation) - -HAVE_OCEANMASSTRANSPORT=no -if test "x$OCEANMASSTRANSPORT" = "xyes"; then - HAVE_OCEANMASSTRANSPORT=yes - AC_DEFINE([_HAVE_OCEANMASSTRANSPORT_],[1],[with Oceanmasstransport capability]) -fi -AM_CONDITIONAL([OCEANMASSTRANSPORT], [test x$HAVE_OCEANMASSTRANSPORT = xyes]) -AC_MSG_RESULT($HAVE_OCEANMASSTRANSPORT) + dnl }}} dnl with-Melting{{{ + AC_ARG_WITH([Melting], + AS_HELP_STRING([--with-Melting = YES], [compile with Melting capabilities (default is yes)]), + [MELTING=$withval],[MELTING=yes]) + AC_MSG_CHECKING(for Melting capability compilation) + HAVE_MELTING=no + if test "x$MELTING" = "xyes"; then + HAVE_MELTING=yes + AC_DEFINE([_HAVE_MELTING_],[1],[with Melting capability]) -fi + +fi + AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes]) + AC_MSG_RESULT($HAVE_MELTING) + +dnl }}} +dnl with-Oceantransport{{{ + +AC_ARG_WITH([Oceantransport], + + AS_HELP_STRING([--with-Oceantransport = YES], [compile with Oceantransport capabilities (default is yes)]), + + [OCEANTRANSPORT=$withval],[OCEANTRANSPORT=yes]) + +AC_MSG_CHECKING(for Oceantransport capability compilation) + + +HAVE_OCEANTRANSPORT=no + +if test "x$OCEANTRANSPORT" = "xyes"; then + + HAVE_OCEANTRANSPORT=yes + + AC_DEFINE([_HAVE_OCEANTRANSPORT_],[1],[with Oceantransport capability]) + +fi + +AM_CONDITIONAL([OCEANTRANSPORT], [test x$HAVE_OCEANTRANSPORT = xyes]) + +AC_MSG_RESULT($HAVE_OCEANTRANSPORT) + +dnl }}} +dnl with-Recovery{{{ + +AC_ARG_WITH([Recovery], + + AS_HELP_STRING([--with-Recovery = YES], [compile with Recovery capabilities (default is yes)]), + + [RECOVERY=$withval],[RECOVERY=yes]) + +AC_MSG_CHECKING(for Recovery capability compilation) + + +HAVE_RECOVERY=no + +if test "x$RECOVERY" = "xyes"; then + + HAVE_RECOVERY=yes + + AC_DEFINE([_HAVE_RECOVERY_],[1],[with Recovery capability]) + +fi + +AM_CONDITIONAL([RECOVERY], [test x$HAVE_RECOVERY = xyes]) + +AC_MSG_RESULT($HAVE_RECOVERY) + dnl }}} dnl with-Sampling{{{ @@ -484,115 +868,211 @@ dnl }}} -dnl with-Sealevelrise{{{ -AC_ARG_WITH([Sealevelrise], - AS_HELP_STRING([--with-Sealevelrise = YES], [compile with Sealevelrise capabilities (default is yes)]), - [SEALEVELRISE=$withval],[SEALEVELRISE=yes]) -AC_MSG_CHECKING(for Sealevelrise capability compilation) - -HAVE_SEALEVELRISE=no -if test "x$SEALEVELRISE" = "xyes"; then - HAVE_SEALEVELRISE=yes - AC_DEFINE([_HAVE_SEALEVELRISE_],[1],[with Sealevelrise capability]) -fi -AM_CONDITIONAL([SEALEVELRISE], [test x$HAVE_SEALEVELRISE = xyes]) -AC_MSG_RESULT($HAVE_SEALEVELRISE) +dnl with-Sealevelchange{{{ + +AC_ARG_WITH([Sealevelchange], + + AS_HELP_STRING([--with-Sealevelchange = YES], [compile with Sealevelchange capabilities (default is yes)]), + + [SEALEVELCHANGE=$withval],[SEALEVELCHANGE=yes]) + +AC_MSG_CHECKING(for Sealevelchange capability compilation) + + +HAVE_SEALEVELCHANGE=no + +if test "x$SEALEVELCHANGE" = "xyes"; then + + HAVE_SEALEVELCHANGE=yes + + AC_DEFINE([_HAVE_SEALEVELCHANGE_],[1],[with Sealevelchange capability]) + +fi + +AM_CONDITIONAL([SEALEVELCHANGE], [test x$HAVE_SEALEVELCHANGE = xyes]) + +AC_MSG_RESULT($HAVE_SEALEVELCHANGE) + dnl }}} dnl with-Smb{{{ + AC_ARG_WITH([Smb], + AS_HELP_STRING([--with-Smb = YES], [compile with Smb capabilities (default is yes)]), + [SMB=$withval],[SMB=yes]) + AC_MSG_CHECKING(for Smb capability compilation) + HAVE_SMB=no + if test "x$SMB" = "xyes"; then + HAVE_SMB=yes + AC_DEFINE([_HAVE_SMB_],[1],[with Smb capability]) -fi + +fi + AM_CONDITIONAL([SMB], [test x$HAVE_SMB = xyes]) + AC_MSG_RESULT($HAVE_SMB) + dnl }}} dnl with-Smooth{{{ + AC_ARG_WITH([Smooth], + AS_HELP_STRING([--with-Smooth = YES], [compile with Smooth capabilities (default is yes)]), + [SMOOTH=$withval],[SMOOTH=yes]) + AC_MSG_CHECKING(for Smooth capability compilation) + HAVE_SMOOTH=no + if test "x$SMOOTH" = "xyes"; then + HAVE_SMOOTH=yes + AC_DEFINE([_HAVE_SMOOTH_],[1],[with Smooth capability]) -fi + +fi + AM_CONDITIONAL([SMOOTH], [test x$HAVE_SMOOTH = xyes]) + AC_MSG_RESULT($HAVE_SMOOTH) + dnl }}} dnl with-Stressbalance{{{ + AC_ARG_WITH([Stressbalance], + AS_HELP_STRING([--with-Stressbalance = YES], [compile with Stressbalance capabilities (default is yes)]), + [STRESSBALANCE=$withval],[STRESSBALANCE=yes]) + AC_MSG_CHECKING(for Stressbalance capability compilation) + HAVE_STRESSBALANCE=no + if test "x$STRESSBALANCE" = "xyes"; then + HAVE_STRESSBALANCE=yes + AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalance capability]) -fi + +fi + AM_CONDITIONAL([STRESSBALANCE], [test x$HAVE_STRESSBALANCE = xyes]) + AC_MSG_RESULT($HAVE_STRESSBALANCE) + dnl }}} dnl with-StressbalanceSIA{{{ + AC_ARG_WITH([StressbalanceSIA], + AS_HELP_STRING([--with-StressbalanceSIA = YES], [compile with StressbalanceSIA capabilities (default is yes)]), + [STRESSBALANCESIA=$withval],[STRESSBALANCESIA=yes]) + AC_MSG_CHECKING(for StressbalanceSIA capability compilation) + HAVE_STRESSBALANCESIA=no + if test "x$STRESSBALANCESIA" = "xyes"; then + HAVE_STRESSBALANCESIA=yes + AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIA capability]) -fi + +fi + AM_CONDITIONAL([STRESSBALANCESIA], [test x$HAVE_STRESSBALANCESIA = xyes]) + AC_MSG_RESULT($HAVE_STRESSBALANCESIA) + dnl }}} dnl with-StressbalanceVertical{{{ + AC_ARG_WITH([StressbalanceVertical], + AS_HELP_STRING([--with-StressbalanceVertical = YES], [compile with StressbalanceVertical capabilities (default is yes)]), + [STRESSBALANCEVERTICAL=$withval],[STRESSBALANCEVERTICAL=yes]) + AC_MSG_CHECKING(for StressbalanceVertical capability compilation) + HAVE_STRESSBALANCEVERTICAL=no + if test "x$STRESSBALANCEVERTICAL" = "xyes"; then + HAVE_STRESSBALANCEVERTICAL=yes + AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVertical capability]) -fi + +fi + AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes]) + AC_MSG_RESULT($HAVE_STRESSBALANCEVERTICAL) + dnl }}} dnl with-Thermal{{{ + AC_ARG_WITH([Thermal], + AS_HELP_STRING([--with-Thermal = YES], [compile with Thermal capabilities (default is yes)]), + [THERMAL=$withval],[THERMAL=yes]) + AC_MSG_CHECKING(for Thermal capability compilation) + HAVE_THERMAL=no + if test "x$THERMAL" = "xyes"; then + HAVE_THERMAL=yes + AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermal capability]) -fi + +fi + AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes]) + AC_MSG_RESULT($HAVE_THERMAL) + dnl }}} dnl with-UzawaPressure{{{ + AC_ARG_WITH([UzawaPressure], + AS_HELP_STRING([--with-UzawaPressure = YES], [compile with UzawaPressure capabilities (default is yes)]), + [UZAWAPRESSURE=$withval],[UZAWAPRESSURE=yes]) + AC_MSG_CHECKING(for UzawaPressure capability compilation) + HAVE_UZAWAPRESSURE=no + if test "x$UZAWAPRESSURE" = "xyes"; then + HAVE_UZAWAPRESSURE=yes + AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressure capability]) -fi + +fi + AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes]) + AC_MSG_RESULT($HAVE_UZAWAPRESSURE) + dnl }}} Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 26098) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 26099) @@ -212,4 +212,5 @@ ./modules/NodalValuex/NodalValuex.cpp \ ./modules/VertexCoordinatesx/VertexCoordinatesx.cpp \ + ./modules/ElementCoordinatesx/ElementCoordinatesx.cpp \ ./modules/OutputResultsx/OutputResultsx.cpp \ ./modules/InputDepthAverageAtBasex/InputDepthAverageAtBasex.cpp \ @@ -480,5 +481,5 @@ issm_sources += ./analyses/MasstransportAnalysis.cpp endif -if OCEANMASSTRANSPORT +if OCEANTRANSPORT issm_sources += ./analyses/OceantransportAnalysis.cpp endif @@ -517,6 +518,45 @@ endif #}}} -#Gia sources (only if have fortran){{{ -if GIA +#Love sources (only if have fortran){{{ +if LOVE +if FORTRAN +issm_sources += \ + ./cores/love_core.cpp \ + ./analyses/LoveAnalysis.cpp \ + ./modules/FourierLoveCorex/FourierLoveCorex.cpp \ + ./modules/FourierLoveCorex/lnb_param.f90 \ + ./modules/FourierLoveCorex/model.f90 \ + ./modules/FourierLoveCorex/util.f90 \ + ./modules/FourierLoveCorex/lovenb_sub.f90 \ + ./modules/FourierLoveCorex/love_numbers.f90 +endif +endif +#}}} +#Esa sources {{{ +if ESA +issm_sources += \ + ./cores/esa_core.cpp \ + ./analyses/EsaAnalysis.cpp +endif +#}}} +#Oceansources {{{ +if OCEAN +issm_sources += ./modules/OceanExchangeDatax/OceanExchangeDatax.cpp +endif +#}}} +#Sampling sources {{{ +if SAMPLING +issm_sources += \ + ./cores/sampling_core.cpp \ + ./analyses/SamplingAnalysis.cpp +endif +# }}} +#Slc sources {{{ +if SEALEVELCHANGE +issm_sources += \ + ./cores/sealevelchange_core.cpp \ + ./analyses/SealevelchangeAnalysis.cpp + +#gia ivins physics (only if have fortran) if FORTRAN issm_sources += \ @@ -530,45 +570,5 @@ ./modules/GiaDeflectionCorex/what0.f endif -endif -#}}} -#Love sources (only if have fortran){{{ -if LOVE -if FORTRAN -issm_sources += \ - ./cores/love_core.cpp \ - ./analyses/LoveAnalysis.cpp \ - ./modules/FourierLoveCorex/FourierLoveCorex.cpp \ - ./modules/FourierLoveCorex/lnb_param.f90 \ - ./modules/FourierLoveCorex/model.f90 \ - ./modules/FourierLoveCorex/util.f90 \ - ./modules/FourierLoveCorex/lovenb_sub.f90 \ - ./modules/FourierLoveCorex/love_numbers.f90 -endif -endif -#}}} -#Esa sources {{{ -if ESA -issm_sources += \ - ./cores/esa_core.cpp \ - ./analyses/EsaAnalysis.cpp -endif -#}}} -#Oceansources {{{ -if OCEAN -issm_sources += ./modules/OceanExchangeDatax/OceanExchangeDatax.cpp -endif -#}}} -#Sampling sources {{{ -if SAMPLING -issm_sources += \ - ./cores/sampling_core.cpp \ - ./analyses/SamplingAnalysis.cpp -endif -# }}} -#Slc sources {{{ -if SEALEVELRISE -issm_sources += \ - ./cores/sealevelchange_core.cpp \ - ./analyses/SealevelchangeAnalysis.cpp + endif #}}} Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 26099) @@ -14,4 +14,7 @@ switch(analysis_enum){ + #ifdef _HAVE_HYDROLOGYTWS_ + case HydrologyTwsAnalysisEnum : return new HydrologyTwsAnalysis(); + #endif #ifdef _HAVE_ADJOINTBALANCETHICKNESS2_ case AdjointBalancethickness2AnalysisEnum : return new AdjointBalancethickness2Analysis(); @@ -71,15 +74,15 @@ case HydrologyDCInefficientAnalysisEnum : return new HydrologyDCInefficientAnalysis(); #endif - #ifdef _HAVE_HYDROLOGYSHREVE_ - case HydrologyShreveAnalysisEnum : return new HydrologyShreveAnalysis(); - #endif #ifdef _HAVE_HYDROLOGYGLADS_ case HydrologyGlaDSAnalysisEnum : return new HydrologyGlaDSAnalysis(); + #endif + #ifdef _HAVE_HYDROLOGYPISM_ + case HydrologyPismAnalysisEnum : return new HydrologyPismAnalysis(); #endif #ifdef _HAVE_HYDROLOGYSHAKTI_ case HydrologyShaktiAnalysisEnum : return new HydrologyShaktiAnalysis(); #endif - #ifdef _HAVE_HYDROLOGYPISM_ - case HydrologyPismAnalysisEnum : return new HydrologyPismAnalysis(); + #ifdef _HAVE_HYDROLOGYSHREVE_ + case HydrologyShreveAnalysisEnum : return new HydrologyShreveAnalysis(); #endif #ifdef _HAVE_L2PROJECTIONBASE_ @@ -101,8 +104,11 @@ case MeltingAnalysisEnum : return new MeltingAnalysis(); #endif + #ifdef _HAVE_OCEANTRANSPORT_ + case OceantransportAnalysisEnum : return new OceantransportAnalysis(); + #endif #ifdef _HAVE_SAMPLING_ case SamplingAnalysisEnum : return new SamplingAnalysis(); #endif - #ifdef _HAVE_SEALEVELRISE_ + #ifdef _HAVE_SEALEVELCHANGE_ case SealevelchangeAnalysisEnum : return new SealevelchangeAnalysis(); #endif Index: /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/analyses/OceantransportAnalysis.cpp (revision 26099) @@ -18,5 +18,5 @@ }/*}}}*/ int OceantransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ - return 1; + return 3; }/*}}}*/ void OceantransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ @@ -59,4 +59,6 @@ iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaBottomPressureEnum,P1Enum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.bottompressure",BottomPressureEnum); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.dsl",DslEnum); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.str",StrEnum); Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26099) @@ -110,7 +110,9 @@ IssmDouble* partitionice=NULL; IssmDouble* partitionhydro=NULL; + IssmDouble* partitionocean=NULL; IssmDouble* bslcice_partition=NULL; IssmDouble* bslchydro_partition=NULL; - int npartice,nparthydro,nel; + IssmDouble* bslcocean_partition=NULL; + int npartice,nparthydro,npartocean,nel; @@ -159,4 +161,5 @@ xDelete(partitionice); } + iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); @@ -168,4 +171,15 @@ parameters->AddObject(new DoubleMatParam(CumBslcHydroPartitionEnum,bslchydro_partition,nparthydro,1)); xDelete(partitionhydro); + } + + iomodel->FetchData(&npartocean,"md.solidearth.npartocean"); + parameters->AddObject(new IntParam(SolidearthNpartOceanEnum,npartocean)); + if(npartocean){ + iomodel->FetchData(&partitionocean,&nel,NULL,"md.solidearth.partitionocean"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionOceanEnum,partitionocean,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartocean",SolidearthNpartOceanEnum)); + bslcocean_partition=xNewZeroInit(npartocean); + parameters->AddObject(new DoubleMatParam(CumBslcOceanPartitionEnum,bslcocean_partition,npartocean,1)); + xDelete(partitionocean); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26099) @@ -238,4 +238,5 @@ virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; + virtual void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false)=0; virtual int FiniteElement(void)=0; virtual IssmDouble FloatingArea(bool scaled)=0; @@ -374,5 +375,5 @@ virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0; #endif - #ifdef _HAVE_SEALEVELRISE_ + #ifdef _HAVE_SEALEVELCHANGE_ virtual void SetSealevelMasks(SealevelMasks* masks)=0; virtual IssmDouble GetArea3D(void)=0; @@ -383,5 +384,6 @@ virtual IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0; virtual void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0; - virtual void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0; + virtual void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, + IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0; virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26099) @@ -70,4 +70,5 @@ void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method); void ElementResponse(IssmDouble* presponse,int response_enum); + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); int FiniteElement(void); @@ -214,9 +215,10 @@ void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif - #ifdef _HAVE_SEALEVELRISE_ + #ifdef _HAVE_SEALEVELCHANGE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, + IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26099) @@ -54,4 +54,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum,int control_interp){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; int FiniteElement(void); @@ -169,8 +170,9 @@ void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif -#ifdef _HAVE_SEALEVELRISE_ +#ifdef _HAVE_SEALEVELCHANGE_ void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, + IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26099) @@ -53,4 +53,5 @@ IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); void FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3); @@ -176,8 +177,9 @@ void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif -#ifdef _HAVE_SEALEVELRISE_ +#ifdef _HAVE_SEALEVELCHANGE_ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, + IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea){_error_("not implemented yet!");}; IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26099) @@ -5411,5 +5411,5 @@ /*}}}*/ #endif -#ifdef _HAVE_SEALEVELRISE_ +#ifdef _HAVE_SEALEVELCHANGE_ IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ @@ -5541,5 +5541,6 @@ } /*}}}*/ -void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){ /*{{{*/ +void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ + /*diverse:*/ int gsize; @@ -5684,5 +5685,4 @@ /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelchangeIndicesEnum,this->lid,indices,gsize); this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,gsize); if(computeelastic){ @@ -6219,4 +6219,446 @@ } /*}}}*/ +void Tria::ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze, bool spherical){/*{{{*/ + + /*Look for x,y,z coordinates:*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xe,ye,ze; + + /*Element XYZ list*/ + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + + /*Find centroid:*/ + xe=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + ye=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + ze=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; + + if(!spherical){ + vxe->SetValue(this->sid,xe,INS_VAL); + vye->SetValue(this->sid,ye,INS_VAL); + vze->SetValue(this->sid,ze,INS_VAL); + } + else _error_("spherical coordinates not supported yet!"); + return; +} +/*}}}*/ + +//new logic +void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ + + /*diverse:*/ + int nel; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble area,planetarea,planetradius; + IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) + IssmDouble rho_earth; + IssmDouble lati,longi; + IssmDouble constant; + IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; + int sidlist[NUMVERTICES]; + + #ifdef _HAVE_RESTRICT_ + IssmDouble* __restrict__ G=NULL; + IssmDouble* __restrict__ GU=NULL; + IssmDouble* __restrict__ GN=NULL; + IssmDouble* __restrict__ GE=NULL; + IssmDouble* __restrict__ G_elastic_precomputed=NULL; + IssmDouble* __restrict__ G_rigid_precomputed=NULL; + IssmDouble* __restrict__ U_elastic_precomputed=NULL; + IssmDouble* __restrict__ H_elastic_precomputed=NULL; + IssmDouble* __restrict__ lates=NULL; + IssmDouble* __restrict__ longes=NULL; + #else + IssmDouble* G=NULL; + IssmDouble* GU=NULL; + IssmDouble* GN=NULL; + IssmDouble* GE=NULL; + IssmDouble* G_elastic_precomputed=NULL; + IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* U_elastic_precomputed=NULL; + IssmDouble* H_elastic_precomputed=NULL; + IssmDouble* lates=NULL; + IssmDouble* longes=NULL; + #endif + + /*elastic green function:*/ + int index; + int M; + + /*Computational flags:*/ + bool computerigid = false; + bool computeelastic = false; + int horiz; + + /*recover parameters: */ + rho_earth=FindParam(MaterialsEarthDensityEnum); + this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); + this->parameters->FindParam(&nel,MeshNumberofelementsEnum); + this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); + this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); + this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + + /*compute area and add to inputs:*/ + area=GetAreaSpherical(); + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); + if(!computerigid)return; + + /*recover precomputed green function kernels:*/ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + + /*allocate arrays:*/ + G=xNewZeroInit(3*nel); + longes=xNewZeroInit(nel); + lates=xNewZeroInit(nel); + + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + + GU=xNewZeroInit(3*nel); + if(horiz){ + GN=xNewZeroInit(3*nel); + GE=xNewZeroInit(3*nel); + } + } + + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*compute centroids of all elements:*/ + for(int i=0;iGetVerticesSidList(&sidlist[0]); + + for(int e=0;eM_PI)delLambda=2*M_PI-delLambda; + alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); + index=reCast( alpha/M_PI*reCast(M-1) ); + + /*Rigid earth gravitational perturbation: */ + G[i*nel+e]+=G_rigid_precomputed[index]; + + if(computeelastic){ + G[i*nel+e]+=G_elastic_precomputed[index]; + } + G[i*nel+e]=G[i*nel+e]*constant; + + /*Elastic components:*/ + if(computeelastic){ + GU[i*nel+e] = constant * U_elastic_precomputed[index]; + if(horiz){ + /*Compute azimuths, both north and east components: */ + x = xx[sid]; y = yy[sid]; z = zz[sid]; + if(latitude[sid]==90){ + x=1e-12; y=1e-12; + } + if(latitude[sid]==-90){ + x=1e-12; y=1e-12; + } + dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + + GN[i*nel+e] = constant*H_elastic_precomputed[index]*N_azim; + GE[i*nel+e] = constant*H_elastic_precomputed[index]*E_azim; + } + } + } + } + + /*Add in inputs:*/ + this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3); + if(computeelastic){ + this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3); + if(horiz){ + this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3); + this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3); + } + } + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); + + /*Free allocations:*/ + #ifdef _HAVE_RESTRICT_ + delete G; + if(computeelastic){ + delete GU; + if(horiz){ + delete GN; + delete GE; + } + } + delete lates; + delete longes; + #else + xDelete(G); + if(computeelastic){ + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); + } + } + delete lates; + delete longes; + #endif + + return; +} +/*}}}*/ +void Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ + + /*diverse:*/ + IssmDouble area; + IssmDouble phi_ice=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic + IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic + IssmDouble I,W,BP; //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4) + bool notfullygrounded=false; + bool scaleoceanarea= false; + bool computerigid= false; + int glfraction=1; + int npartice,nparthydro,npartocean; + + /*flags:*/ + bool isocean=false; + bool ishydro=false; + bool isice=false; + + /*output: */ + IssmDouble bslcice=0; + IssmDouble bslchydro=0; + IssmDouble bslcbp=0; + + /*ice properties: */ + IssmDouble rho_ice,rho_water,rho_freshwater; + + /*constants:*/ + IssmDouble constant=0; + + /*recover parameters:*/ + this->parameters->FindParam(&isice,TransientIsmasstransportEnum); + this->parameters->FindParam(&isocean,TransientIsoceantransportEnum); + this->parameters->FindParam(&ishydro,TransientIshydrologyEnum); + + /*early return if we are not on an ice cap, and we are not requesting + *hydrology or bottom pressure loads :*/ + if(!isocean && !ishydro){ + if(!masks->isiceonly[this->lid]){ + #ifdef _ISSM_DEBUG_ + constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif + /*barystatic_contribution[0]+=0; + barystatic_contribution[1]+=0; + barystatic_contribution[2]+=0;*/ + return; + } + } + + /*early return if we are fully floating and we are not doing bottom pressure loads:*/ + if(!isocean){ + if (masks->isfullyfloating[this->lid]){ + constant=0; + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif + /*barystatic_contribution[0]+=0; + barystatic_contribution[1]+=0; + barystatic_contribution[2]+=0;*/ + return; + } + } + + /*early return if we are not on the ocean and we are not doing ice mass transport of + * hydrology:*/ + if(!isice && !ishydro){ + if(!masks->isoceanin[this->lid]){ + /*barystatic_contribution[0]+=0; + barystatic_contribution[1]+=0; + barystatic_contribution[2]+=0;*/ + return; + } + + } + + /*if we are an ice shelf, are we fully grounded or not? (used later):*/ + if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. + + /*recover some parameters:*/ + rho_ice=FindParam(MaterialsRhoIceEnum); + rho_water=FindParam(MaterialsRhoSeawaterEnum); + this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); + this->parameters->FindParam(&npartice,SolidearthNpartIceEnum); + this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ + this->Element::GetInputValue(&area,AreaEnum); + + /*Deal with ice loads if we are on grounded ice:*/ + if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/ + + /*Compute fraction of the element that is grounded: */ + if(notfullygrounded){ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + phi_ice=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + if(glfraction==0)phi_ice=1; +#ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&phi_ice,P0Enum); +#endif + } + else phi_ice=1.0; + + /*Inform mask: */ + constant+=100; //1 for ice. + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); + #endif + + /*Retrieve surface load for ice: */ + Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); + + /*/Average ice thickness over grounded area of the element only: {{{*/ + if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); + else{ + IssmDouble total_weight=0; + bool mainlyfloating = true; + int point1; + IssmDouble fraction1,fraction2; + + /*Recover portion of element that is grounded*/ + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + + /* Start looping on the number of gaussian points and average over these gaussian points: */ + total_weight=0; + I=0; + while(gauss->next()){ + IssmDouble Ig=0; + deltathickness_input->GetInputValue(&Ig,gauss); + I+=Ig*gauss->weight; + total_weight+=gauss->weight; + } + I=I/total_weight; + delete gauss; + } + /*}}}*/ + + /*Compute barystatic contribution:*/ + _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + bslcice = rho_ice*area*I/(oceanarea*rho_water); + _assert_(!xIsNan(bslcice)); + + /*Transfer thickness change into kg/m^2:*/ + I=I*rho_ice*phi_ice; + } /*}}}*/ + + /*Deal with water loads if we are on ground:*/ + if(!masks->isfullyfloating[this->lid]){ + + /*Retrieve water height at vertices: */ + Input* deltathickness_input=this->GetInput(DeltaTwsEnum); + if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!"); + deltathickness_input->GetInputAverage(&W); + + /*Compute barystatic component:*/ + _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water); + _assert_(!xIsNan(bslchydro)); + + /*convert from m to kg/m^2:*/ + W=W*rho_freshwater*phi_water; + } + + /*Deal with ocean bottom pressures:*/ + if(masks->isoceanin[this->lid]){ + + /*Retrieve bottom pressure change and average over the element: */ + Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum); + if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!"); + bottompressure_change_input->GetInputAverage(&BP); + + /*Compute barystatic component:*/ + bslcbp = rho_water*area*BP/(oceanarea*rho_water); + + /*convert from m to kg/m^2:*/ + BP=BP*rho_water; + } + + /*Plug all loads into total load vector:*/ + localloads[this->lid]+=I+W+BP; + + /*Plug bslcice into barystatic contribution vector:*/ + if(barystatic_contribution_onpartition){ + int idi=reCast(partition[this->Sid()])+1; + int idj=0; + idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL); + idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL); + idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL); + } + + barystatic_contribution[0]+=bslcice; + barystatic_contribution[1]+=bslchydro; + barystatic_contribution[2]+=bslcbp; + +} +/*}}}*/ +IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/ + + /*sal green function:*/ + IssmDouble* G=NULL; + IssmDouble Sealevel[NUMVERTICES]={0,0,0}; + + bool sal = false; + int size; + int nel; + + this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&nel,MeshNumberofelementsEnum); + if(sal){ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); + + for(int i=0;iAddInput(SealevelRSLEnum,Sealevel,P1Enum); + } + +} /*}}}*/ + #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26099) @@ -70,4 +70,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum,int control_interp); void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum); + void ElementCoordinates(Vector* vxe,Vector* vye,Vector* vze,bool spherical=false); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -160,9 +161,11 @@ void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz); #endif - #ifdef _HAVE_SEALEVELRISE_ + #ifdef _HAVE_SEALEVELCHANGE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); + void SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); @@ -170,5 +173,7 @@ void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); + void SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea); + IssmDouble SealevelchangeConvolution(IssmDouble* loads); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26099) @@ -844,8 +844,9 @@ case TransientSolutionEnum:{ /*We have multiple analyses here, process one by one*/ - bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslc,isesa,issampling; + bool isSIA,isFS,isthermal,isenthalpy,ismasstransport,isoceantransport,isgroundingline,isstressbalance,ismovingfront,ishydrology,isdamage,issmb,isslc,isesa,issampling; iomodel->FindConstant(&isthermal,"md.transient.isthermal"); iomodel->FindConstant(&ismovingfront,"md.transient.ismovingfront"); iomodel->FindConstant(&ismasstransport,"md.transient.ismasstransport"); + iomodel->FindConstant(&isoceantransport,"md.transient.isoceantransport"); iomodel->FindConstant(&isstressbalance,"md.transient.isstressbalance"); iomodel->FindConstant(&isgroundingline,"md.transient.isgroundingline"); @@ -892,4 +893,7 @@ if(isdamage){ analyses_temp[numanalyses++]=DamageEvolutionAnalysisEnum; + } + if(isoceantransport){ + analyses_temp[numanalyses++]=OceantransportAnalysisEnum; } if(isslc){ @@ -4744,5 +4748,5 @@ /*}}}*/ #endif -#ifdef _HAVE_SEALEVELRISE_ +#ifdef _HAVE_SEALEVELCHANGE_ void FemModel::SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/ @@ -4756,4 +4760,5 @@ IssmDouble oceanarea_cpu = 0.; int bp_compute_fingerprints= 0; + bool isoceantransport=false; Vector* bslcice_partition=NULL; @@ -4765,7 +4770,11 @@ IssmDouble* bslchydro_partition_serial=NULL; IssmDouble* partitionhydro=NULL; + bool istws=0; int nparthydro; - bool istws; - + + int npartocean; + Vector* bslcocean_partition=NULL; + IssmDouble* bslcocean_partition_serial=NULL; + IssmDouble* partitionocean=NULL; /*Initialize temporary vector that will be used to sum barystatic components @@ -4800,4 +4809,31 @@ bslchydro_partition= new Vector(nparthydro); } + + this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum); + if(npartocean){ + this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum); + bslchydro_partition= new Vector(npartocean); + } + /*For later: + npartbarystatic=npartice; + if(nparthydro>npartbarystatic)npartbarystatic=nparthydro; + if(npartocean>npartbarystatic)npartbarystatic=npartocean; + bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3); + + bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0; + for(Object* & object : this->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea); + } + MPI Bcast localloads -> loads + + for(Object* & object : this->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeConvolution(loads); + } + */ + + + @@ -4821,5 +4857,6 @@ /*Call the barystatic sea level change core for bottom pressures: */ this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum); - if(bp_compute_fingerprints){ + this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum); + if(bp_compute_fingerprints && isoceantransport){ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); @@ -5108,4 +5145,8 @@ element->GiaDeflection(wg,dwgdt, matlitho, x,y); } + + /*Assemble parallel vector:*/ + dwgdt->Assemble(); + wg->Assemble(); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 26099) @@ -163,5 +163,5 @@ void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif - #ifdef _HAVE_SEALEVELRISE_ + #ifdef _HAVE_SEALEVELCHANGE_ void SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition, SealevelMasks* masks); void SealevelchangeSal(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); Index: /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp (revision 26099) @@ -82,4 +82,11 @@ int index= EnumToIndex(enum_in); if(this->inputs[index])this->inputs[index]->DeepEcho(); + return; +} +/*}}}*/ +void Inputs::Echo(int enum_in){/*{{{*/ + _printf_("Inputs Echo:\n"); + int index= EnumToIndex(enum_in); + if(this->inputs[index])this->inputs[index]->Echo(); return; } Index: /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h (revision 26098) +++ /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h (revision 26099) @@ -50,4 +50,5 @@ void DeepEcho(int enum_in); void Echo(void); + void Echo(int enum_in); bool Exist(int enum_type); void GetInputsInterpolations(int* pnuminputs,int** pinterpolations,int** penum); Index: /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26099) @@ -21,4 +21,5 @@ bool dakota_analysis; int solution_type; + Vector* ug = NULL; /*activate configuration*/ @@ -41,9 +42,8 @@ * and OceantransportSpcstrEnum inputs in each element, assemble into a vector and feed to * InputUpdateFromSolutionx which will deal with accumulating such inputs:*/ - Vector* ug = NULL; GetSolutionFromInputsx(&ug,femmodel); InputUpdateFromSolutionx(femmodel,ug); - delete ug; + SolidEarthOceanUpdates(femmodel); if(solution_type==OceantransportSolutionEnum)femmodel->RequestedDependentsx(); @@ -51,4 +51,8 @@ /*profiler*/ femmodel->profiler->Stop(OCEANTRANSPORTCORE); + + /*free ressources:*/ + delete ug; + } /*}}}*/ void SolidEarthOceanUpdates(FemModel* femmodel){ /*{{{*/ Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26099) @@ -267,4 +267,5 @@ /*get some of the updates into elements:*/ InputUpdateFromVectorx(femmodel,U_grd,SealevelUGrdEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,N_grd,SealevelEnum,VertexSIdEnum); InputUpdateFromVectorx(femmodel,N_grd,SealevelNGrdEnum,VertexSIdEnum); if(RSLg)InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); @@ -306,4 +307,5 @@ int step; int computesealevel=0; + bool isocean=false; IssmDouble time; @@ -316,4 +318,8 @@ femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); if (!computesealevel)return; + + /*early return if we have no ocean transport:*/ + femmodel->parameters->FindParam(&isocean,TransientIsoceantransportEnum); + if(!isocean)return; /*Verbose: */ @@ -406,6 +412,11 @@ IssmDouble *yy = NULL; IssmDouble *zz = NULL; + IssmDouble *xxe = NULL; + IssmDouble *yye = NULL; + IssmDouble *zze = NULL; + int horiz; bool geometrydone = false; + bool optim=false; @@ -424,4 +435,5 @@ /*first, recover lat,long and radius vectors from vertices: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + ElementCoordinatesx(&xxe,&yye,&zze,femmodel->elements); if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); @@ -430,5 +442,6 @@ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); - element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz); + element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze); + //if(optim) element->SealevelchangeGeometryOptim(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze); } Index: /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26099) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp (revision 26099) @@ -0,0 +1,48 @@ +/*!\file ElementCoordinatesx + * \brief: compute a vector xe,ye and ze of element centroid coordinates by + * marching through all our elements. + */ + +#include "./ElementCoordinatesx.h" + +#include "../../shared/shared.h" +#include "../../toolkits/toolkits.h" + +void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,Elements* elements,bool spherical) { + + /*figure out how many vertices we have: */ + int numberofelements=elements->NumberOfElements(); + + Vector* vxe=new Vector(numberofelements); + Vector* vye=new Vector(numberofelements); + Vector* vze=new Vector(numberofelements); + + /*march through our elements: */ + for(Object* & object : elements->objects){ + Element* element=(Element*)object; + element->ElementCoordinates(vxe,vye,vze,spherical); + } + + /*Assemble*/ + vxe->Assemble(); + vye->Assemble(); + vze->Assemble(); + + /*serialize: */ + IssmDouble* xe=vxe->ToMPISerial(); + IssmDouble* ye=vye->ToMPISerial(); + IssmDouble* ze=vze->ToMPISerial(); + + /*Free ressources: */ + delete vxe; + delete vye; + delete vze; + + /*output: */ + if(pxe) *pxe=xe; + else xDelete(xe); + if(pye) *pye=ye; + else xDelete(ye); + if(pze) *pze=ze; + else xDelete(ze); +} Index: /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h (revision 26099) +++ /issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h (revision 26099) @@ -0,0 +1,12 @@ +/*!\file: ElementCoordinatesx.h + */ + +#ifndef _ELEMENT_COORDINATESX_H +#define _ELEMENT_COORDINATESX_H + +#include "../../classes/classes.h" + +/* local prototypes: */ +void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,Elements* elements,bool spherical=false); + +#endif /* _ELEMENT_COORDINATESX_H */ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 26098) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 26099) @@ -103,4 +103,5 @@ #include "./UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h" #include "./VertexCoordinatesx/VertexCoordinatesx.h" +#include "./ElementCoordinatesx/ElementCoordinatesx.h" #ifdef _HAVE_OCEAN_ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26098) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26099) @@ -113,6 +113,8 @@ syn keyword cConstant CumBslcIceEnum syn keyword cConstant CumBslcHydroEnum +syn keyword cConstant CumBslcOceanEnum syn keyword cConstant CumBslcIcePartitionEnum syn keyword cConstant CumBslcHydroPartitionEnum +syn keyword cConstant CumBslcOceanPartitionEnum syn keyword cConstant CumGmtslcEnum syn keyword cConstant CumGmslcEnum @@ -290,4 +292,5 @@ syn keyword cConstant MeshNumberoflayersEnum syn keyword cConstant MeshNumberofverticesEnum +syn keyword cConstant MeshNumberofelementsEnum syn keyword cConstant MigrationMaxEnum syn keyword cConstant ModelIdEnum @@ -338,5 +341,7 @@ syn keyword cConstant SolidearthPartitionIceEnum syn keyword cConstant SolidearthPartitionHydroEnum +syn keyword cConstant SolidearthPartitionOceanEnum syn keyword cConstant SolidearthNpartIceEnum +syn keyword cConstant SolidearthNpartOceanEnum syn keyword cConstant SolidearthNpartHydroEnum syn keyword cConstant SolidearthPlanetRadiusEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26098) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26099) @@ -107,6 +107,8 @@ CumBslcIceEnum, CumBslcHydroEnum, + CumBslcOceanEnum, CumBslcIcePartitionEnum, CumBslcHydroPartitionEnum, + CumBslcOceanPartitionEnum, CumGmtslcEnum, CumGmslcEnum, @@ -284,4 +286,5 @@ MeshNumberoflayersEnum, MeshNumberofverticesEnum, + MeshNumberofelementsEnum, MigrationMaxEnum, ModelIdEnum, @@ -332,5 +335,7 @@ SolidearthPartitionIceEnum, SolidearthPartitionHydroEnum, + SolidearthPartitionOceanEnum, SolidearthNpartIceEnum, + SolidearthNpartOceanEnum, SolidearthNpartHydroEnum, SolidearthPlanetRadiusEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26099) @@ -115,6 +115,8 @@ case CumBslcIceEnum : return "CumBslcIce"; case CumBslcHydroEnum : return "CumBslcHydro"; + case CumBslcOceanEnum : return "CumBslcOcean"; case CumBslcIcePartitionEnum : return "CumBslcIcePartition"; case CumBslcHydroPartitionEnum : return "CumBslcHydroPartition"; + case CumBslcOceanPartitionEnum : return "CumBslcOceanPartition"; case CumGmtslcEnum : return "CumGmtslc"; case CumGmslcEnum : return "CumGmslc"; @@ -292,4 +294,5 @@ case MeshNumberoflayersEnum : return "MeshNumberoflayers"; case MeshNumberofverticesEnum : return "MeshNumberofvertices"; + case MeshNumberofelementsEnum : return "MeshNumberofelements"; case MigrationMaxEnum : return "MigrationMax"; case ModelIdEnum : return "ModelId"; @@ -340,5 +343,7 @@ case SolidearthPartitionIceEnum : return "SolidearthPartitionIce"; case SolidearthPartitionHydroEnum : return "SolidearthPartitionHydro"; + case SolidearthPartitionOceanEnum : return "SolidearthPartitionOcean"; case SolidearthNpartIceEnum : return "SolidearthNpartIce"; + case SolidearthNpartOceanEnum : return "SolidearthNpartOcean"; case SolidearthNpartHydroEnum : return "SolidearthNpartHydro"; case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26098) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26099) @@ -115,6 +115,8 @@ else if (strcmp(name,"CumBslcIce")==0) return CumBslcIceEnum; else if (strcmp(name,"CumBslcHydro")==0) return CumBslcHydroEnum; + else if (strcmp(name,"CumBslcOcean")==0) return CumBslcOceanEnum; else if (strcmp(name,"CumBslcIcePartition")==0) return CumBslcIcePartitionEnum; else if (strcmp(name,"CumBslcHydroPartition")==0) return CumBslcHydroPartitionEnum; + else if (strcmp(name,"CumBslcOceanPartition")==0) return CumBslcOceanPartitionEnum; else if (strcmp(name,"CumGmtslc")==0) return CumGmtslcEnum; else if (strcmp(name,"CumGmslc")==0) return CumGmslcEnum; @@ -135,10 +137,10 @@ else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum; else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; - else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; - else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"DslModel")==0) return DslModelEnum; + if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; + else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; + else if (strcmp(name,"DslModel")==0) return DslModelEnum; else if (strcmp(name,"DslModelid")==0) return DslModelidEnum; else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; else if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum; - else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; - else if (strcmp(name,"LoveG0")==0) return LoveG0Enum; else stage=3; } if(stage==3){ - if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; + if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; + else if (strcmp(name,"LoveG0")==0) return LoveG0Enum; + else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; @@ -298,4 +300,5 @@ else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum; else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum; + else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum; else if (strcmp(name,"MigrationMax")==0) return MigrationMaxEnum; else if (strcmp(name,"ModelId")==0) return ModelIdEnum; @@ -346,5 +349,7 @@ else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum; else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum; + else if (strcmp(name,"SolidearthPartitionOcean")==0) return SolidearthPartitionOceanEnum; else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum; + else if (strcmp(name,"SolidearthNpartOcean")==0) return SolidearthNpartOceanEnum; else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum; else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum; @@ -378,13 +383,13 @@ else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; - else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum; else if (strcmp(name,"SealevelchangeUElastic")==0) return SealevelchangeUElasticEnum; else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; + else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; else if (strcmp(name,"SettingsCheckpointFrequency")==0) return SettingsCheckpointFrequencyEnum; else if (strcmp(name,"SettingsResultsOnNodes")==0) return SettingsResultsOnNodesEnum; @@ -501,13 +506,13 @@ else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; - else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; else if (strcmp(name,"TransientIssampling")==0) return TransientIssamplingEnum; else if (strcmp(name,"TransientIsslc")==0) return TransientIsslcEnum; else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; + else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; @@ -624,13 +629,13 @@ else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum; else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; - else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; else if (strcmp(name,"EsaRotationrate")==0) return EsaRotationrateEnum; else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum; else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum; else if (strcmp(name,"EsaStrainrateyy")==0) return EsaStrainrateyyEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum; + else if (strcmp(name,"EsaUmotion")==0) return EsaUmotionEnum; else if (strcmp(name,"EsaXmotion")==0) return EsaXmotionEnum; else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum; @@ -747,13 +752,13 @@ else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; - else if (strcmp(name,"Sample")==0) return SampleEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Sample")==0) return SampleEnum; else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum; else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum; + else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum; else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum; else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum; @@ -870,13 +875,13 @@ else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; - else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SmbTa")==0) return SmbTaEnum; else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; + else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; @@ -993,13 +998,13 @@ else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; - else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; + else if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum; @@ -1116,13 +1121,13 @@ else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Contact")==0) return ContactEnum; - else if (strcmp(name,"Contour")==0) return ContourEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Contour")==0) return ContourEnum; else if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum; else if (strcmp(name,"ControlInputMaxs")==0) return ControlInputMaxsEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"ControlInputMins")==0) return ControlInputMinsEnum; + else if (strcmp(name,"ControlInputMins")==0) return ControlInputMinsEnum; else if (strcmp(name,"ControlInputValues")==0) return ControlInputValuesEnum; else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum; @@ -1239,13 +1244,13 @@ else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; - else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; + else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; else if (strcmp(name,"LoveHi")==0) return LoveHiEnum; else if (strcmp(name,"LoveHr")==0) return LoveHrEnum; @@ -1362,13 +1367,13 @@ else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; - else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; else if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"SealevelInertiaTensorYZ")==0) return SealevelInertiaTensorYZEnum; + else if (strcmp(name,"SealevelInertiaTensorYZ")==0) return SealevelInertiaTensorYZEnum; else if (strcmp(name,"SealevelInertiaTensorZZ")==0) return SealevelInertiaTensorZZEnum; else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum; Index: /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh (revision 26098) +++ /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh (revision 26099) @@ -17,5 +17,5 @@ #Deal with Analyses -if false ; then +if true ; then #Build EnumToAnalysis.cpp {{{ #Header @@ -39,5 +39,5 @@ #core cat temp | grep [a-zA-Z0-9]Analysis | \ - grep -v DefaultAnalysis | grep -v FlaimAnalysis | grep -v SurfaceSlopeAnalysis | grep -v BedSlopeAnalysis | \ + grep -v DefaultAnalysis | grep -v RecoveryAnalysis | grep -v FlaimAnalysis | grep -v SurfaceSlopeAnalysis | grep -v BedSlopeAnalysis | \ awk '{print "\t\t#ifdef _HAVE_"toupper(substr($1,1,length($1)-12))"_\n\t\t" "case " $1" : return new " substr($1,1,length($1)-4) "();\n\t\t#endif"}' \ >> $ISSM_DIR/src/c/analyses/EnumToAnalysis.cpp Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 26098) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 26099) @@ -22,4 +22,6 @@ sealevel = NaN; bottompressure = NaN; + dsl = NaN; + str = NaN; sample = NaN; end @@ -39,4 +41,6 @@ self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1); self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1); + self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1); + self.str=project3d(md,'vector',self.str,'type','node','layer',1); %Lithostatic pressure by default @@ -66,5 +70,7 @@ if ismember('OceantransportAnalysis',analyses) if strcmp(solution,'TransientSolution') & md.transient.isslc & md.transient.isoceantransport, - md = checkfield(md,'fieldname','initialization.bottompressure','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','initialization.bottompressure','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','initialization.dsl','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','initialization.str','NaN',1,'Inf',1,'size',[1 1]); end end @@ -156,4 +162,7 @@ fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]'); fielddisplay(self,'sample','Realization of a Gaussian random field'); + fielddisplay(self,'bottompressure','Bottom pressures'); + fielddisplay(self,'dsl','Dynamic sea level.'); + fielddisplay(self,'str','Steric sea level.'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -167,4 +176,6 @@ WriteData(fid,prefix,'object',self,'fieldname','sealevel','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','bottompressure','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','str','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','dsl','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','temperature','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','waterfraction','format','DoubleMat','mattype',1); Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 26098) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 26099) @@ -15,4 +15,5 @@ partitionice = []; partitionhydro = []; + partitionocean = []; end methods (Static) @@ -55,4 +56,5 @@ self.partitionice=[]; self.partitionhydro=[]; + self.partitionocean=[]; %no external solutions by default: @@ -94,4 +96,5 @@ fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); + fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution'); if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end self.settings.disp(); @@ -118,4 +121,10 @@ nparthydro=0; end + if ~isempty(self.partitionocean), + npartocean=max(self.partitionocean)+2; + else + npartocean=0; + end + @@ -124,4 +133,6 @@ WriteData(fid,prefix,'object',self,'fieldname','partitionhydro','mattype',1,'format','DoubleMat'); WriteData(fid,prefix,'data',nparthydro,'format','Integer','name','md.solidearth.nparthydro'); + WriteData(fid,prefix,'object',self,'fieldname','partitionocean','mattype',1,'format','DoubleMat'); + WriteData(fid,prefix,'data',npartocean,'format','Integer','name','md.solidearth.npartocean'); self.settings.marshall(prefix,md,fid); Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26098) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26099) @@ -5,19 +5,26 @@ md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh +%Geometry for the bed, arbitrary thickness of 1000: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); +md.geometry.base=md.geometry.bed; +md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1); +md.geometry.surface=md.geometry.bed+md.geometry.thickness; + + %parameterize solidearth solution: %solidearth loading: {{{ -md.masstransport.spcthickness=zeros(md.mesh.numberofelements+1,1); -md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.masstransport.spcthickness=[md.geometry.thickness;0]; md.dsl.global_average_thermosteric_sea_level=[0;0]; md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1); md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %antarctica late=sum(md.mesh.lat(md.mesh.elements),2)/3; longe=sum(md.mesh.long(md.mesh.elements),2)/3; pos=find(late < -80); -md.masstransport.spcthickness(pos)=-100; +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; %greenland pos=find(late>70 & late<80 & longe>-60 & longe<-30); -md.masstransport.spcthickness(pos)=-100; +md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; %elastic loading from love numbers: @@ -35,18 +42,14 @@ md.mask.ocean_levelset=-icemask; -%make sure that the elements that have loads are fully grounded: -pos=find(md.masstransport.spcthickness); -md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; %make sure wherever there is an ice load, that the mask is set to ice: -pos=find(md.masstransport.spcthickness); +pos=find(md.masstransport.spcthickness(1:end-1)); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % }}} -%Geometry for the bed, arbitrary: -md.geometry.bed=-ones(md.mesh.numberofvertices,1); -md.geometry.thickness=zeros(md.mesh.numberofvertices,1); -md.geometry.surface=zeros(md.mesh.numberofvertices,1); -md.geometry.base=zeros(md.mesh.numberofvertices,1); +%time stepping: +md.timestepping.start_time=0; +md.timestepping.time_step=1; +md.timestepping.final_time=1; @@ -55,4 +58,8 @@ md.initialization.vx=zeros(md.mesh.numberofvertices,1); md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1); +md.initialization.dsl=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; %Materials: @@ -76,5 +83,7 @@ md.transient.isthermal=0; md.transient.ismasstransport=1; +md.transient.isoceantransport=1; md.transient.isslc=1; +md.solidearth.requested_outputs={'Sealevel'}; @@ -100,4 +109,5 @@ md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; +md.solidearth.requested_outputs={'Sealevel'}; md=solve(md,'tr'); Selastic=md.results.TransientSolution.Sealevel;