Index: ../trunk-jpl/m4/analyses.m4
===================================================================
--- ../trunk-jpl/m4/analyses.m4	(revision 18949)
+++ ../trunk-jpl/m4/analyses.m4	(revision 18950)
@@ -10,862 +10,466 @@
 [
 
 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 AdjointBalancethicknesscapability])
-
 fi
-
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes])
-
 AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS)
-
 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 AdjointBalancethickness2capability])
-
 fi
-
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes])
-
 AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS2)
-
 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 AdjointHorizcapability])
-
 fi
-
 AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes])
-
 AC_MSG_RESULT($HAVE_ADJOINTHORIZ)
-
 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 Balancethicknesscapability])
-
 fi
-
 AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes])
-
 AC_MSG_RESULT($HAVE_BALANCETHICKNESS)
-
 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 Balancethickness2capability])
-
 fi
-
 AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes])
-
 AC_MSG_RESULT($HAVE_BALANCETHICKNESS2)
-
 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 BalancethicknessSoftcapability])
-
 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 Balancevelocitycapability])
-
 fi
-
 AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes])
-
 AC_MSG_RESULT($HAVE_BALANCEVELOCITY)
-
 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 L2ProjectionEPLcapability])
-
 fi
-
 AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes])
-
 AC_MSG_RESULT($HAVE_L2PROJECTIONEPL)
-
 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 L2ProjectionBasecapability])
-
 fi
-
 AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes])
-
 AC_MSG_RESULT($HAVE_L2PROJECTIONBASE)
-
 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 DamageEvolutioncapability])
-
 fi
-
 AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes])
-
 AC_MSG_RESULT($HAVE_DAMAGEEVOLUTION)
-
 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 Stressbalancecapability])
-
 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 StressbalanceSIAcapability])
-
 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 StressbalanceVerticalcapability])
-
 fi
-
 AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes])
-
 AC_MSG_RESULT($HAVE_STRESSBALANCEVERTICAL)
-
 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 Enthalpycapability])
-
 fi
-
 AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes])
-
 AC_MSG_RESULT($HAVE_ENTHALPY)
-
 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 HydrologyShrevecapability])
-
 fi
-
 AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes])
-
 AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE)
-
 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 HydrologyDCInefficientcapability])
-
 fi
-
 AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes])
-
 AC_MSG_RESULT($HAVE_HYDROLOGYDCINEFFICIENT)
-
 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 HydrologyDCEfficientcapability])
-
 fi
-
 AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes])
-
 AC_MSG_RESULT($HAVE_HYDROLOGYDCEFFICIENT)
-
 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 Meltingcapability])
-
 fi
-
 AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes])
-
 AC_MSG_RESULT($HAVE_MELTING)
-
 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 Masstransportcapability])
-
 fi
-
 AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
-
 AC_MSG_RESULT($HAVE_MASSTRANSPORT)
-
 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 FreeSurfaceBasecapability])
-
 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 FreeSurfaceTopcapability])
-
 fi
-
 AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes])
-
 AC_MSG_RESULT($HAVE_FREESURFACETOP)
-
 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 ExtrudeFromBasecapability])
-
 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 ExtrudeFromTopcapability])
-
 fi
-
 AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes])
-
 AC_MSG_RESULT($HAVE_EXTRUDEFROMTOP)
-
 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 DepthAveragecapability])
-
 fi
-
 AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes])
-
 AC_MSG_RESULT($HAVE_DEPTHAVERAGE)
-
 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 Smoothcapability])
-
 fi
-
 AM_CONDITIONAL([SMOOTH], [test x$HAVE_SMOOTH = xyes])
-
 AC_MSG_RESULT($HAVE_SMOOTH)
-
 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 Thermalcapability])
-
 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 UzawaPressurecapability])
-
 fi
-
 AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes])
-
 AC_MSG_RESULT($HAVE_UZAWAPRESSURE)
-
 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-Seaice{{{
-
 AC_ARG_WITH([Seaice],
-
 	AS_HELP_STRING([--with-Seaice = YES], [compile with Seaice capabilities (default is yes)]),
-
 	[SEAICE=$withval],[SEAICE=yes])
-
 AC_MSG_CHECKING(for Seaice capability compilation)
 
-
 HAVE_SEAICE=no 
-
 if test "x$SEAICE" = "xyes"; then
-
 	HAVE_SEAICE=yes
-
 	AC_DEFINE([_HAVE_SEAICE_],[1],[with Seaicecapability])
-
 fi
-
 AM_CONDITIONAL([SEAICE], [test x$HAVE_SEAICE = xyes])
-
 AC_MSG_RESULT($HAVE_SEAICE)
-
 dnl }}}
 dnl with-Meshdeformation{{{
-
 AC_ARG_WITH([Meshdeformation],
-
 	AS_HELP_STRING([--with-Meshdeformation = YES], [compile with Meshdeformation capabilities (default is yes)]),
-
 	[MESHDEFORMATION=$withval],[MESHDEFORMATION=yes])
-
 AC_MSG_CHECKING(for Meshdeformation capability compilation)
 
-
 HAVE_MESHDEFORMATION=no 
-
 if test "x$MESHDEFORMATION" = "xyes"; then
-
 	HAVE_MESHDEFORMATION=yes
-
 	AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformationcapability])
-
 fi
-
 AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes])
-
 AC_MSG_RESULT($HAVE_MESHDEFORMATION)
-
 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 Levelsetcapability])
-
 fi
-
 AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes])
-
 AC_MSG_RESULT($HAVE_LEVELSET)
-
 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 Extrapolationcapability])
-
 fi
-
 AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes])
-
 AC_MSG_RESULT($HAVE_EXTRAPOLATION)
-
 dnl }}}
 dnl with-LsfReinitialization{{{
-
 AC_ARG_WITH([LsfReinitialization],
-
 	AS_HELP_STRING([--with-LsfReinitialization = YES], [compile with LsfReinitialization capabilities (default is yes)]),
-
 	[LSFREINITIALIZATION=$withval],[LSFREINITIALIZATION=yes])
-
 AC_MSG_CHECKING(for LsfReinitialization capability compilation)
 
-
 HAVE_LSFREINITIALIZATION=no 
-
 if test "x$LSFREINITIALIZATION" = "xyes"; then
-
 	HAVE_LSFREINITIALIZATION=yes
-
 	AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitializationcapability])
-
 fi
-
 AM_CONDITIONAL([LSFREINITIALIZATION], [test x$HAVE_LSFREINITIALIZATION = xyes])
-
 AC_MSG_RESULT($HAVE_LSFREINITIALIZATION)
-
 dnl }}}
 
 ])
Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18950)
@@ -183,6 +183,7 @@
 				iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum);
 				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum);
+				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum);
 			}
 			else{
 				iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18950)
@@ -341,6 +341,7 @@
 	SurfaceforcingsDelta18oSurfaceEnum,
 	SurfaceforcingsIsdelta18oEnum,
 	SurfaceforcingsPrecipitationsPresentdayEnum,
+	SurfaceforcingsPrecipitationsLgmEnum,
 	SurfaceforcingsTemperaturesPresentdayEnum,
 	SurfaceforcingsTemperaturesLgmEnum,
 	SurfaceforcingsPrecipitationEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18950)
@@ -347,6 +347,7 @@
 		case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface";
 		case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o";
 		case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday";
+		case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm";
 		case SurfaceforcingsTemperaturesPresentdayEnum : return "SurfaceforcingsTemperaturesPresentday";
 		case SurfaceforcingsTemperaturesLgmEnum : return "SurfaceforcingsTemperaturesLgm";
 		case SurfaceforcingsPrecipitationEnum : return "SurfaceforcingsPrecipitation";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18950)
@@ -353,6 +353,7 @@
 	      else if (strcmp(name,"SurfaceforcingsDelta18oSurface")==0) return SurfaceforcingsDelta18oSurfaceEnum;
 	      else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum;
 	      else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum;
+	      else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum;
 	      else if (strcmp(name,"SurfaceforcingsTemperaturesPresentday")==0) return SurfaceforcingsTemperaturesPresentdayEnum;
 	      else if (strcmp(name,"SurfaceforcingsTemperaturesLgm")==0) return SurfaceforcingsTemperaturesLgmEnum;
 	      else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
@@ -381,11 +382,11 @@
 	      else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
 	      else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
 	      else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
-	      else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
+	      if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
+	      else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
 	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
@@ -504,11 +505,11 @@
 	      else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
 	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
 	      else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
-	      else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
+	      if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
+	      else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
 	      else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum;
 	      else if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
 	      else if (strcmp(name,"MassconaxpbyBeta")==0) return MassconaxpbyBetaEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
 	      else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
 	      else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
-	      else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
+	      if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
+	      else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
 	      else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
 	      else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
 	      else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
@@ -750,11 +751,11 @@
 	      else if (strcmp(name,"Intersect")==0) return IntersectEnum;
 	      else if (strcmp(name,"Nodal")==0) return NodalEnum;
 	      else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
-	      else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
+	      if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
+	      else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
 	      else if (strcmp(name,"RootPath")==0) return RootPathEnum;
 	      else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
 	      else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum;
Index: ../trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Elements/elements.h	(revision 18950)
@@ -14,8 +14,9 @@
 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
 IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
 				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
+				IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
 				IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
-				IssmDouble h, IssmDouble s, IssmDouble rho_ice,IssmDouble rho_water, IssmDouble desfac,
+				IssmDouble h, IssmDouble s, IssmDouble desfac,
 				IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
 				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime);
 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
Index: ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp	(revision 18950)
@@ -16,9 +16,12 @@
   IssmDouble delta18oLapseRate=-6.2*pow(10.,-3);
   IssmDouble glacialindex; // used to vary present day temperature
 
-  glacialindex = 0.;//(Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))
-  //  /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent));
+  // CURENTLY THE TEMPERATURE COMPUTED HERE ARE NOT USED.
+  // Tdiff is used isnstead in PddSurfaceMassBalance.cpp to get the interpolated temperature
 
+  glacialindex = (Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))
+    /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); // Tarasov 2004 paper
+
   for (int imonth = 0; imonth<12; imonth++){
     monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1.-glacialindex)*TemperaturesPresentday[imonth];
     //monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-273.15-0)/10)); //equation from rob's paper
Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 18950)
@@ -7,12 +7,16 @@
 
 IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, 
 				IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
+				IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
 				IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, 
-				IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, 
-				IssmDouble rho_water, IssmDouble desfac, IssmDouble s0t,
+				IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, 
+				IssmDouble s0t,
 				IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, 
 				IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ 
 
+  // CURENTLY monthlytemperatures and monthlyprec ARE NOT USED HERE.
+  // THEY ARE REPLACE BY tdiffh No usage of deltO18 anymore
+
   // output:
   IssmDouble B;    // surface mass balance, melt+accumulation
   int    iqj,imonth;
@@ -23,18 +27,18 @@
   IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
   IssmDouble water; //water=rain + snowmelt 
   IssmDouble runoff; //meltwater only, does not include rain 
-  IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
   //IssmDouble  sealev=0.;         // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
   //IssmDouble  Pfac=0.5,Tdiff=0.5;
   IssmDouble  rtlaps;
   // IssmDouble lapser=6.5         // lapse rate
-  // IssmDouble desfac = 0.5;      // desert elevation factor
+  // IssmDouble desfac = 0.3;      // desert elevation factor
   // IssmDouble s0p=0.;            // should be set to elevation from precip source
   // IssmDouble s0t=0.;         // should be set to elevation from temperature source
   IssmDouble tdiffh;  
   IssmDouble st;             // elevation between altitude of the temp record and current altitude
   IssmDouble sp;             // elevation between altitude of the prec record and current altitude
+  IssmDouble deselcut=1.0;
 
   // PDD and PD constants and variables
   IssmDouble siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
@@ -44,10 +48,11 @@
   IssmDouble DT = 0.02;
   IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
 
-  IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
-  IssmDouble qm = 0.;        // snow part of the precipitation 
-  IssmDouble qmt = 0.;       // precipitation without desertification effect adjustment
-  IssmDouble qmp = 0.;       // desertification taken into account
+  IssmDouble premt;     // monthly precipitation (m of ice equivalent)
+  IssmDouble q, qmpt;   // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
+  IssmDouble qm = 0.;   // snow part of the precipitation 
+  IssmDouble qmt = 0.;  // precipitation without desertification effect adjustment
+  IssmDouble qmp = 0.;  // desertification taken into account
   IssmDouble pdd = 0.;     
   IssmDouble frzndd = 0.;  
 
@@ -67,7 +72,6 @@
   IssmDouble fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
   IssmDouble pddtj, hmx2;
 
-  sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
 
   /*PDD constant*/
   siglim = 2.5*signorm; 
@@ -97,19 +101,21 @@
 	pd = 0.;}
 
       /******exp des/elev precip reduction*******/
-      sp=(s-s0p)/1000.; // deselev effect is wrt chng in topo
+      sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
       if (sp>0.0){q = exp(-desfac*sp);}
       else {q = 1.0;}
 
-      qmt= qmt + monthlyprec[imonth]*sconv;  //*sconv to convert in m of ice equivalent per month
-      qmpt= q*monthlyprec[imonth]*sconv;           
+      premt =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime));   // [m/month]
+
+      qmt= qmt + premt;  
+      qmpt= q*premt;           
       qmp= qmp + qmpt;
       qm= qm + qmpt*pd;
 
       /*********compute PDD************/
       // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
       // gaussian=T_m, so ndd=-(Tsurf-pdd)
-      if (iqj>5 &&  iqj<9){ Tsum=Tsum+tstar;} 
+      if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;} 
 
       if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
       else if (tstar> -siglim){
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18950)
@@ -2214,15 +2214,25 @@
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+   IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
+   IssmDouble sconv; //rhow_rain/rhoi / 12 months
+   
+   /*Get material parameters :*/
+   rho_ice=matpar->GetRhoIce();
+   rho_water=matpar->GetRhoFreshwater();
 
+   sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
+
    /*Recover monthly temperatures and precipitation and Present day and LGm ones*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
    Input*     input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input3);
    Input*     input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input4);
+   Input*     input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum);  _assert_(input5);
+   Input*     input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);   _assert_(input6);
    GaussTria* gauss=new GaussTria();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
@@ -2233,11 +2243,15 @@
        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
-       monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion in m/y
+       monthlyprec[iv][month]=monthlyprec[iv][month]*sconv*yts; // convertion in m/y
        input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
        TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
        input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
        TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
+       input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+       PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; 
+       input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
+       PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];  
      }
    }
 
@@ -2250,10 +2264,6 @@
   GetInputListOnVertices(&h[0],ThicknessEnum);
   GetInputListOnVertices(&s[0],SurfaceEnum);
 
-  /*Get material parameters :*/
-  rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-  rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
-
   /*Get other pdd parameters*/
   desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
   s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
@@ -2265,7 +2275,8 @@
    for (int iv = 0; iv<NUMVERTICES; iv++){
      agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
 				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
-				  pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water, 
+				  &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], 
+				  pdds, pds, signorm, yts, h[iv], s[iv],
 				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
    }
 
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18949)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18950)
@@ -2011,30 +2011,45 @@
    IssmDouble agd[NUMVERTICES];             // surface mass balance
    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
+   IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    IssmDouble PfacTime,TdiffTime,sealevTime;
+   IssmDouble sconv; //rhow_rain/rhoi / 12 months
 
+   /*Get material parameters :*/
+   rho_ice=matpar->GetRhoIce();
+   rho_water=matpar->GetRhoFreshwater();
+
+   sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
+
    /*Recover monthly temperatures and precipitation*/
    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
    Input*     input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);  _assert_(input3);
    Input*     input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);   _assert_(input4);
+   Input*     input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum);  _assert_(input5);
+   Input*     input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);   _assert_(input6);
    GaussPenta* gauss=new GaussPenta();
    IssmDouble time,yts;
    this->parameters->FindParam(&time,TimeEnum);
    this->parameters->FindParam(&yts,ConstantsYtsEnum);
+   
    for(int month=0;month<12;month++) {
      for(int iv=0;iv<NUMVERTICES;iv++) {
        gauss->GaussVertex(iv);
        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
-       monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
+       monthlyprec[iv][month]=monthlyprec[iv][month]*sconv; // convertion to m/yr
        input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
        TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
        input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
        TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
+       input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
+       PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; 
+       input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
+       PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];  
      }
    } 
 
@@ -2045,12 +2060,8 @@
 
   /*Recover info at the vertices: */
   GetInputListOnVertices(&h[0],ThicknessEnum);
-  GetInputListOnVertices(&s[0],SurfaceEnum);
+  GetInputListOnVertices(&s[0],SurfaceEnum); 
 
-  /*Get material parameters :*/
-  rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-  rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
-
   /*Get other pdd parameters*/
   desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
   s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
@@ -2062,7 +2073,8 @@
    for (int iv = 0; iv < NUMVERTICES; iv++){
      agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],  
 				  &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], 
-				  pdds,pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water,
+				  &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], 
+				  pdds,pds, signorm, yts, h[iv], s[iv],
 				  desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
    }
 
Index: ../trunk-jpl/src/m/plot/processmesh.m
===================================================================
--- ../trunk-jpl/src/m/plot/processmesh.m	(revision 18949)
+++ ../trunk-jpl/src/m/plot/processmesh.m	(revision 18950)
@@ -20,16 +20,15 @@
 	return;
 end
 
-%first load x,y, etc ... to speed up plot
-if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon'),
-	x=md.mesh.x;
-	if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end
-	y=md.mesh.y;
-	if isprop(md.mesh,'y2d'), y2d=md.mesh.y2d; end
-else
-	x=md.mesh.long;
-	y=md.mesh.lat;
-end
+	if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon'),
+		x=md.mesh.x;
+		x2d=md.mesh.x2d; %if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end
+		y=md.mesh.y;
+		y2d=md.mesh.y2d; %if isprop(md.mesh,'y2d'), y2d=md.mesh.y2d; end
+	else
+		x=md.mesh.long;
+		y=md.mesh.lat;
+	end
 
 if isprop(md.mesh,'z'),
 	z=md.mesh.z;
@@ -41,8 +40,8 @@
 	z=md.(z);
 end
 
-if isprop(md.mesh,'elements2d'), elements2d=md.mesh.elements2d; end
-elements=md.mesh.elements;
+	elements2d=md.mesh.elements2d;%if isprop(md.mesh,'elements2d'), elements2d=md.mesh.elements2d; end
+	elements=md.mesh.elements;
 
 %is it a 2d plot?
 if md.mesh.dimension()==2,
Index: ../trunk-jpl/src/m/classes/SMBpdd.m
===================================================================
--- ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18949)
+++ ../trunk-jpl/src/m/classes/SMBpdd.m	(revision 18950)
@@ -10,17 +10,18 @@
 		desfac                    = 0;
 		s0p                       = 0;
 		s0t                       = 0;
-      rlaps                     = 0;
-      rlapslgm                  = 0;                
-      Pfac                      = NaN;
-      Tdiff                     = NaN;
-      sealev                    = NaN;
+                rlaps                     = 0;
+                rlapslgm                  = 0;                
+                Pfac                      = NaN;
+                Tdiff                     = NaN;
+                sealev                    = NaN;
 		isdelta18o                = 0;
 		delta18o                  = NaN;
 		delta18o_surface          = NaN;
 		temperatures_presentday   = NaN;
 		temperatures_lgm          = NaN;
 		precipitations_presentday = NaN;
+		precipitations_lgm        = NaN;
 	end
 	methods
 		function obj = SMBpdd(varargin) % {{{
@@ -38,6 +39,7 @@
 			if(self.isdelta18o),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end
 			if(self.isdelta18o),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
 			if(self.isdelta18o),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
+			if(self.isdelta18o),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end
 
 
 		end % }}}
@@ -79,6 +81,7 @@
 					md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
 					md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
 					md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
+					md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
 				end
 			end
 		end % }}}
@@ -95,11 +98,12 @@
 			fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year)');
 			fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year)');
 			fielddisplay(obj,'sealev','sea level [m], 1D(year)');
-			fielddisplay(obj,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated');
-			fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]');
+			fielddisplay(obj,'monthlytemperatures',['CURRENTLY NOT USED monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
+			fielddisplay(obj,'precipitation',['CURRENTLY NOT USED monthly surface precipitation [m/yr water eq]']);
 			fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated');
 			fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated');
 			fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
+			fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
 			fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated');
 			fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated');
 
@@ -124,6 +128,7 @@
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
+				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1);
 				WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1);
 			else
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18949)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18950)
@@ -339,6 +339,7 @@
 def SurfaceforcingsDelta18oSurfaceEnum(): return StringToEnum("SurfaceforcingsDelta18oSurface")[0]
 def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0]
 def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0]
+def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0]
 def SurfaceforcingsTemperaturesPresentdayEnum(): return StringToEnum("SurfaceforcingsTemperaturesPresentday")[0]
 def SurfaceforcingsTemperaturesLgmEnum(): return StringToEnum("SurfaceforcingsTemperaturesLgm")[0]
 def SurfaceforcingsPrecipitationEnum(): return StringToEnum("SurfaceforcingsPrecipitation")[0]
