Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 19171) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 19172) @@ -19,5 +19,5 @@ if test "x$ADJOINTBALANCETHICKNESS" = "xyes"; then HAVE_ADJOINTBALANCETHICKNESS=yes - AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethicknesscapability]) + AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethickness capability]) fi AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes]) @@ -33,5 +33,5 @@ if test "x$ADJOINTBALANCETHICKNESS2" = "xyes"; then HAVE_ADJOINTBALANCETHICKNESS2=yes - AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2capability]) + AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2 capability]) fi AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes]) @@ -47,5 +47,5 @@ if test "x$ADJOINTHORIZ" = "xyes"; then HAVE_ADJOINTHORIZ=yes - AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHorizcapability]) + AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHoriz capability]) fi AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes]) @@ -61,5 +61,5 @@ if test "x$BALANCETHICKNESS" = "xyes"; then HAVE_BALANCETHICKNESS=yes - AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethicknesscapability]) + AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethickness capability]) fi AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes]) @@ -75,5 +75,5 @@ if test "x$BALANCETHICKNESS2" = "xyes"; then HAVE_BALANCETHICKNESS2=yes - AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2capability]) + AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2 capability]) fi AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes]) @@ -89,5 +89,5 @@ if test "x$BALANCETHICKNESSSOFT" = "xyes"; then HAVE_BALANCETHICKNESSSOFT=yes - AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoftcapability]) + AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoft capability]) fi AM_CONDITIONAL([BALANCETHICKNESSSOFT], [test x$HAVE_BALANCETHICKNESSSOFT = xyes]) @@ -103,5 +103,5 @@ if test "x$BALANCEVELOCITY" = "xyes"; then HAVE_BALANCEVELOCITY=yes - AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocitycapability]) + AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocity capability]) fi AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes]) @@ -117,5 +117,5 @@ if test "x$L2PROJECTIONEPL" = "xyes"; then HAVE_L2PROJECTIONEPL=yes - AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPLcapability]) + AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPL capability]) fi AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes]) @@ -131,5 +131,5 @@ if test "x$L2PROJECTIONBASE" = "xyes"; then HAVE_L2PROJECTIONBASE=yes - AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBasecapability]) + AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBase capability]) fi AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes]) @@ -145,5 +145,5 @@ if test "x$DAMAGEEVOLUTION" = "xyes"; then HAVE_DAMAGEEVOLUTION=yes - AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolutioncapability]) + AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolution capability]) fi AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes]) @@ -159,5 +159,5 @@ if test "x$STRESSBALANCE" = "xyes"; then HAVE_STRESSBALANCE=yes - AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalancecapability]) + AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalance capability]) fi AM_CONDITIONAL([STRESSBALANCE], [test x$HAVE_STRESSBALANCE = xyes]) @@ -173,5 +173,5 @@ if test "x$STRESSBALANCESIA" = "xyes"; then HAVE_STRESSBALANCESIA=yes - AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIAcapability]) + AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIA capability]) fi AM_CONDITIONAL([STRESSBALANCESIA], [test x$HAVE_STRESSBALANCESIA = xyes]) @@ -187,5 +187,5 @@ if test "x$STRESSBALANCEVERTICAL" = "xyes"; then HAVE_STRESSBALANCEVERTICAL=yes - AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVerticalcapability]) + AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVertical capability]) fi AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes]) @@ -201,5 +201,5 @@ if test "x$ENTHALPY" = "xyes"; then HAVE_ENTHALPY=yes - AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpycapability]) + AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpy capability]) fi AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes]) @@ -215,5 +215,5 @@ if test "x$HYDROLOGYSHREVE" = "xyes"; then HAVE_HYDROLOGYSHREVE=yes - AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShrevecapability]) + AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShreve capability]) fi AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes]) @@ -229,5 +229,5 @@ if test "x$HYDROLOGYDCINEFFICIENT" = "xyes"; then HAVE_HYDROLOGYDCINEFFICIENT=yes - AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficientcapability]) + AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficient capability]) fi AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes]) @@ -243,5 +243,5 @@ if test "x$HYDROLOGYDCEFFICIENT" = "xyes"; then HAVE_HYDROLOGYDCEFFICIENT=yes - AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficientcapability]) + AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficient capability]) fi AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes]) @@ -257,5 +257,5 @@ if test "x$MELTING" = "xyes"; then HAVE_MELTING=yes - AC_DEFINE([_HAVE_MELTING_],[1],[with Meltingcapability]) + AC_DEFINE([_HAVE_MELTING_],[1],[with Melting capability]) fi AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes]) @@ -271,5 +271,5 @@ if test "x$MASSTRANSPORT" = "xyes"; then HAVE_MASSTRANSPORT=yes - AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransportcapability]) + AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransport capability]) fi AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes]) @@ -285,5 +285,5 @@ if test "x$FREESURFACEBASE" = "xyes"; then HAVE_FREESURFACEBASE=yes - AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBasecapability]) + AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBase capability]) fi AM_CONDITIONAL([FREESURFACEBASE], [test x$HAVE_FREESURFACEBASE = xyes]) @@ -299,5 +299,5 @@ if test "x$FREESURFACETOP" = "xyes"; then HAVE_FREESURFACETOP=yes - AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTopcapability]) + AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTop capability]) fi AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes]) @@ -313,5 +313,5 @@ if test "x$EXTRUDEFROMBASE" = "xyes"; then HAVE_EXTRUDEFROMBASE=yes - AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBasecapability]) + AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBase capability]) fi AM_CONDITIONAL([EXTRUDEFROMBASE], [test x$HAVE_EXTRUDEFROMBASE = xyes]) @@ -327,5 +327,5 @@ if test "x$EXTRUDEFROMTOP" = "xyes"; then HAVE_EXTRUDEFROMTOP=yes - AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTopcapability]) + AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTop capability]) fi AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes]) @@ -341,5 +341,5 @@ if test "x$DEPTHAVERAGE" = "xyes"; then HAVE_DEPTHAVERAGE=yes - AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAveragecapability]) + AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAverage capability]) fi AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes]) @@ -355,5 +355,5 @@ if test "x$SMOOTH" = "xyes"; then HAVE_SMOOTH=yes - AC_DEFINE([_HAVE_SMOOTH_],[1],[with Smoothcapability]) + AC_DEFINE([_HAVE_SMOOTH_],[1],[with Smooth capability]) fi AM_CONDITIONAL([SMOOTH], [test x$HAVE_SMOOTH = xyes]) @@ -369,5 +369,5 @@ if test "x$THERMAL" = "xyes"; then HAVE_THERMAL=yes - AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermalcapability]) + AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermal capability]) fi AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes]) @@ -383,5 +383,5 @@ if test "x$UZAWAPRESSURE" = "xyes"; then HAVE_UZAWAPRESSURE=yes - AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressurecapability]) + AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressure capability]) fi AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes]) @@ -397,22 +397,8 @@ if test "x$GIA" = "xyes"; then HAVE_GIA=yes - AC_DEFINE([_HAVE_GIA_],[1],[with Giacapability]) + AC_DEFINE([_HAVE_GIA_],[1],[with Gia capability]) 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{{{ @@ -425,5 +411,5 @@ if test "x$MESHDEFORMATION" = "xyes"; then HAVE_MESHDEFORMATION=yes - AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformationcapability]) + AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformation capability]) fi AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes]) @@ -439,5 +425,5 @@ if test "x$LEVELSET" = "xyes"; then HAVE_LEVELSET=yes - AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelsetcapability]) + AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelset capability]) fi AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes]) @@ -453,5 +439,5 @@ if test "x$EXTRAPOLATION" = "xyes"; then HAVE_EXTRAPOLATION=yes - AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolationcapability]) + AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolation capability]) fi AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes]) @@ -467,5 +453,5 @@ if test "x$LSFREINITIALIZATION" = "xyes"; then HAVE_LSFREINITIALIZATION=yes - AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitializationcapability]) + AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitialization capability]) fi AM_CONDITIONAL([LSFREINITIALIZATION], [test x$HAVE_LSFREINITIALIZATION = xyes]) Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 19171) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 19172) @@ -114,4 +114,5 @@ ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\ ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\ + ./shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp\ ./shared/Elements/DrainageFunctionWaterfraction.cpp\ ./shared/String/DescriptorIndex.cpp\ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 19172) @@ -120,5 +120,5 @@ int stabilization,finiteelement,smb_model; bool dakota_analysis; - bool isdelta18o,ismungsm; + bool isdelta18o,ismungsm,isd18opd; bool isgroundingline; bool islevelset; @@ -180,4 +180,5 @@ iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); + iomodel->Constant(&isd18opd,SurfaceforcingsIsd18opdEnum); iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); if(isdelta18o || ismungsm){ @@ -187,9 +188,12 @@ iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum); } + else if (isd18opd){ + iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum); + } else{ iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); } - break; case SMBgradientsEnum: Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 19172) @@ -35,5 +35,5 @@ } else{ - IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement); + IoModelToDynamicConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,finiteelement); } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19172) @@ -183,4 +183,5 @@ virtual void Delta18oParameterization(void)=0; virtual void MungsmtpParameterization(void)=0; + virtual void Delta18opdParameterization(void)=0; virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 19172) @@ -713,4 +713,69 @@ this->inputs->AddInput(NewPrecipitationInput); + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Penta::Delta18opdParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int i; + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble TemperaturesPresentday[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble Delta18oTime; + IssmDouble dpermil; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Get some pdd parameters*/ + dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2); + GaussPenta* gauss=new GaussPenta(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivAddTimeInput(newmonthinput1,time+imonth/12.*yts); + + for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 19172) @@ -60,4 +60,5 @@ void Delta18oParameterization(void); void MungsmtpParameterization(void); + void Delta18opdParameterization(void); void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 19172) @@ -55,4 +55,5 @@ void Delta18oParameterization(void){_error_("not implemented yet");}; void MungsmtpParameterization(void){_error_("not implemented yet");}; + void Delta18opdParameterization(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 19172) @@ -55,4 +55,5 @@ void Delta18oParameterization(void){_error_("not implemented yet");}; void MungsmtpParameterization(void){_error_("not implemented yet");}; + void Delta18opdParameterization(void){_error_("not implemented yet");}; IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19172) @@ -750,4 +750,66 @@ &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0], &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], + &monthlytemperatures[iv][0], &monthlyprec[iv][0]); + } + + /*Update inputs*/ + TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum); + TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum); + for (int imonth=0;imonth<12;imonth++) { + for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); + + for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Tria::Delta18opdParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int i; + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble TemperaturesPresentday[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble Delta18oTime; + IssmDouble dpermil; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Get some pdd parameters*/ + dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2); + GaussTria* gauss=new GaussTria(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivConstant(&this->rlaps,SurfaceforcingsRlapsEnum); iomodel->Constant(&this->rlapslgm,SurfaceforcingsRlapslgmEnum); + iomodel->Constant(&this->dpermil,SurfaceforcingsDpermilEnum); break; case SMBgradientsEnum: @@ -134,4 +135,5 @@ _printf_(" rlaps: " << rlaps << "\n"); _printf_(" rlapslgm: " << rlapslgm << "\n"); + _printf_(" dpermil: " << dpermil << "\n"); return; } @@ -179,4 +181,5 @@ matpar->rlaps=this->rlaps; matpar->rlapslgm=this->rlapslgm; + matpar->dpermil=this->dpermil; matpar->sediment_compressibility=this->sediment_compressibility; @@ -273,4 +276,7 @@ case SurfaceforcingsRlapslgmEnum: this->rlapslgm=constant; + break; + case SurfaceforcingsDpermilEnum: + this->dpermil=constant; break; default: @@ -400,4 +406,5 @@ case SurfaceforcingsRlapsEnum: return this->rlaps; case SurfaceforcingsRlapslgmEnum: return this->rlapslgm; + case SurfaceforcingsDpermilEnum: return this->dpermil; case MaterialsLithosphereShearModulusEnum: return this->lithosphere_shear_modulus; case MaterialsLithosphereDensityEnum: return this->lithosphere_density; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 19171) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 19172) @@ -36,4 +36,5 @@ IssmDouble rlaps; IssmDouble rlapslgm; + IssmDouble dpermil; /*hydrology Dual Porous Continuum: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 19172) @@ -22,5 +22,5 @@ char** requestedoutputs = NULL; IssmDouble time; - bool isdelta18o,ismungsm; + bool isdelta18o,ismungsm,isd18opd; /*parameters for mass flux:*/ @@ -107,4 +107,5 @@ parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsd18opdEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum)); @@ -114,4 +115,5 @@ iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); + iomodel->Constant(&isd18opd,SurfaceforcingsIsd18opdEnum); if(ismungsm){ @@ -131,5 +133,4 @@ iomodel->DeleteData(temp,SurfaceforcingsSealevEnum); } - if(isdelta18o){ iomodel->Constant(&yts,ConstantsYtsEnum); @@ -144,4 +145,14 @@ parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oSurfaceEnum,&temp[0],&temp[M],M)); iomodel->DeleteData(temp,SurfaceforcingsDelta18oSurfaceEnum); + } + if(isd18opd){ + iomodel->Constant(&yts,ConstantsYtsEnum); + + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsDelta18oEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsDelta18oEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsDelta18oEnum); + + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDpermilEnum)); } break; Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 19172) @@ -11,5 +11,5 @@ /*Intermediaties*/ int smb_model; - bool isdelta18o,ismungsm; + bool isdelta18o,ismungsm,isd18opd; /*First, get SMB model from parameters*/ @@ -24,4 +24,5 @@ femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum); + femmodel->parameters->FindParam(&isd18opd,SurfaceforcingsIsd18opdEnum); if(isdelta18o){ if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n"); @@ -31,4 +32,8 @@ if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n"); MungsmtpParameterizationx(femmodel); + } + if(isd18opd){ + if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n"); + Delta18opdParameterizationx(femmodel); } if(VerboseSolution()) _printf0_(" call positive degree day module\n"); @@ -128,4 +133,12 @@ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); element->MungsmtpParameterization(); + } + +}/*}}}*/ +void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->Delta18opdParameterization(); } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 19171) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 19172) @@ -13,4 +13,5 @@ void Delta18oParameterizationx(FemModel* femmodel); void MungsmtpParameterizationx(FemModel* femmodel); +void Delta18opdParameterizationx(FemModel* femmodel); void PositiveDegreeDayx(FemModel* femmodel); void SmbHenningx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 19172) @@ -90,133 +90,133 @@ Tsurf = tstar*deltm+Tsurf; - /*********compute PD ****************/ - if (tstar < PDup){ - pd = 1.; - if (tstar >= -siglimc){ pd = pds[reCast(tstar/DT + siglim0c)];}} - else { - pd = 0.;} - - /******exp des/elev precip reduction*******/ - 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; - 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 (tstar >= siglim) {pdd = pdd + tstar*deltm;} - else if (tstar> -siglim){ - pddsig=pdds[reCast(tstar/DT + siglim0)]; - pdd = pdd + pddsig*deltm; - frzndd = frzndd - (tstar-pddsig)*deltm;} - else{frzndd = frzndd - tstar*deltm; } - - /*Assign output pointer*/ - *(monthlytemperatures+imonth) = monthlytemperatures[imonth]; - *(monthlyprec+imonth) = monthlyprec[imonth]; + /*********compute PD ****************/ + if (tstar < PDup){ + pd = 1.; + if (tstar >= -siglimc){ pd = pds[reCast(tstar/DT + siglim0c)];}} + else { + pd = 0.;} + + /******exp des/elev precip reduction*******/ + 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; + 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 (tstar >= siglim) {pdd = pdd + tstar*deltm;} + else if (tstar> -siglim){ + pddsig=pdds[reCast(tstar/DT + siglim0)]; + pdd = pdd + pddsig*deltm; + frzndd = frzndd - (tstar-pddsig)*deltm;} + else{frzndd = frzndd - tstar*deltm; } + + /*Assign output pointer*/ + *(monthlytemperatures+imonth) = monthlytemperatures[imonth]; + *(monthlyprec+imonth) = monthlyprec[imonth]; } // end of seasonal loop //****************************************************************** - saccu = qm; - prect = qmp; // total precipitation during 1 year taking into account des. ef. - Tsum=Tsum/3; - - /***** determine PDD factors *****/ - if(Tsum< -1.) { - snwmf=2.65*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd - smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002) - } - else if(Tsum< 10){ - snwmf = (0.15*Tsum + 2.8)*0.001; - smf = (0.0067*pow((10.-Tsum),3) + 8.3)*0.001; - } - else{ - snwmf=4.3*0.001; - smf=8.3*0.001; - } - snwmf=0.95*snwmf; - smf=0.95*smf; - - /***** compute PDD ablation and refreezing *****/ - pddt = pdd *365; - snwm = snwmf*pddt; // snow that could have been melted in a year - hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz - - if(snwm < saccu) { - water=prect-saccu + snwm; //water=rain + snowmelt - // l 2.2= capillary factor - // Should refreezing be controlled by frzndd or by mean annual Tsurf? - // dCovrLm concept is of warming of active layer (thickness =d@=1- - // >2m) - // problem with water seepage into ice: should be sealed after - // refreezing - // so everything needs to be predicated on 1 year scale, except for - // thermal - // conductivity through ice - // also, need to account that melt season has low accum, so what's - // going to - // hold the meltwater around for refreezing? And melt-time will have - // low seasonal frzndd - - // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002 - - supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice - supcap=min(2.2*(saccu-snwm),water); - runoff=snwm - supice; //meltwater only, does not include rain - } - else { //all snow melted - supice= min(hmx2*CovrLm*frzndd, prect ); - runoff= saccu + smf*(pddt-saccu/snwmf) - supice; - supcap=0; - } - // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it - // except pdd melt heat source is atmosphere, while refreeze is - // ground/ice stored interim - // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0 - // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be - // <0 - // assume ndd>pdd, little melt => little supice - // bottom line: compare for Tsurf<0 : supice and no supice case, - // expect Tsurf difference - // except some of cooling flux comes from atmosphere// - // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C - // does supice make sense when H< 0.1m? then d=thermoactive ice layer //// - // < 0.1 - - // make more sense to just use residual pdd-ndd except that pdd - // residual not clear yet - // frzndd should not be used up by refreezing in snow, so stick in - // supcap. - diffndd=0; - if (frzndd>0) { - diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd); - frzndd=frzndd-diffndd; - } - if(runoff<0){ - saccu= saccu -runoff; - smelt = 0; - precrunoff=prect-saccu; - //here assume pdd residual is 0, => - Tsurf= max(Tsurf,-frzndd); - } - else { - smelt = runoff; - precrunoff=prect-max(0.,supice)-saccu;} - //here really need pdd balance, try 0.5 fudge factor? - //at least runoff>0 => it's fairly warm, so Tsurf is !<<0, - //yet from site plots, can be ice free with Tsurf=-5.5C - if(Tsurf<0) { - Tsurf= min(Tsurf+fsupT*diffndd , 0.);} - - B = -smelt+saccu; - B = B/yts; - pddtj=pddt; + saccu = qm; + prect = qmp; // total precipitation during 1 year taking into account des. ef. + Tsum=Tsum/3; + + /***** determine PDD factors *****/ + if(Tsum< -1.) { + snwmf=2.65*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd + smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002) + } + else if(Tsum< 10){ + snwmf = (0.15*Tsum + 2.8)*0.001; + smf = (0.0067*pow((10.-Tsum),3) + 8.3)*0.001; + } + else{ + snwmf=4.3*0.001; + smf=8.3*0.001; + } + snwmf=0.95*snwmf; + smf=0.95*smf; + + /***** compute PDD ablation and refreezing *****/ + pddt = pdd *365; + snwm = snwmf*pddt; // snow that could have been melted in a year + hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz + + if(snwm < saccu) { + water=prect-saccu + snwm; //water=rain + snowmelt + // l 2.2= capillary factor + // Should refreezing be controlled by frzndd or by mean annual Tsurf? + // dCovrLm concept is of warming of active layer (thickness =d@=1- + // >2m) + // problem with water seepage into ice: should be sealed after + // refreezing + // so everything needs to be predicated on 1 year scale, except for + // thermal + // conductivity through ice + // also, need to account that melt season has low accum, so what's + // going to + // hold the meltwater around for refreezing? And melt-time will have + // low seasonal frzndd + + // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002 + + supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice + supcap=min(2.2*(saccu-snwm),water); + runoff=snwm - supice; //meltwater only, does not include rain + } + else { //all snow melted + supice= min(hmx2*CovrLm*frzndd, prect ); + runoff= saccu + smf*(pddt-saccu/snwmf) - supice; + supcap=0; + } + // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it + // except pdd melt heat source is atmosphere, while refreeze is + // ground/ice stored interim + // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0 + // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be + // <0 + // assume ndd>pdd, little melt => little supice + // bottom line: compare for Tsurf<0 : supice and no supice case, + // expect Tsurf difference + // except some of cooling flux comes from atmosphere// + // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C + // does supice make sense when H< 0.1m? then d=thermoactive ice layer //// + // < 0.1 + + // make more sense to just use residual pdd-ndd except that pdd + // residual not clear yet + // frzndd should not be used up by refreezing in snow, so stick in + // supcap. + diffndd=0; + if (frzndd>0) { + diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd); + frzndd=frzndd-diffndd; + } + if(runoff<0){ + saccu= saccu -runoff; + smelt = 0; + precrunoff=prect-saccu; + //here assume pdd residual is 0, => + Tsurf= max(Tsurf,-frzndd); + } + else { + smelt = runoff; + precrunoff=prect-max(0.,supice)-saccu;} + //here really need pdd balance, try 0.5 fudge factor? + //at least runoff>0 => it's fairly warm, so Tsurf is !<<0, + //yet from site plots, can be ice free with Tsurf=-5.5C + if(Tsurf<0) { + Tsurf= min(Tsurf+fsupT*diffndd , 0.);} + + B = -smelt+saccu; + B = B/yts; + pddtj=pddt; return B; Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 19171) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 19172) @@ -28,4 +28,7 @@ IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); +void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil, + IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.); IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19171) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19172) @@ -348,4 +348,5 @@ SurfaceforcingsIsdelta18oEnum, SurfaceforcingsIsmungsmEnum, + SurfaceforcingsIsd18opdEnum, SurfaceforcingsPrecipitationsPresentdayEnum, SurfaceforcingsPrecipitationsLgmEnum, @@ -361,4 +362,5 @@ SurfaceforcingsTdiffEnum, SurfaceforcingsSealevEnum, + SurfaceforcingsDpermilEnum, SMBgradientsEnum, SurfaceforcingsMonthlytemperaturesEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19172) @@ -354,4 +354,5 @@ case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o"; case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm"; + case SurfaceforcingsIsd18opdEnum : return "SurfaceforcingsIsd18opd"; case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday"; case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm"; @@ -367,4 +368,5 @@ case SurfaceforcingsTdiffEnum : return "SurfaceforcingsTdiff"; case SurfaceforcingsSealevEnum : return "SurfaceforcingsSealev"; + case SurfaceforcingsDpermilEnum : return "SurfaceforcingsDpermil"; case SMBgradientsEnum : return "SMBgradients"; case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19171) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19172) @@ -360,4 +360,5 @@ else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum; else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum; + else if (strcmp(name,"SurfaceforcingsIsd18opd")==0) return SurfaceforcingsIsd18opdEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum; @@ -373,4 +374,5 @@ else if (strcmp(name,"SurfaceforcingsTdiff")==0) return SurfaceforcingsTdiffEnum; else if (strcmp(name,"SurfaceforcingsSealev")==0) return SurfaceforcingsSealevEnum; + else if (strcmp(name,"SurfaceforcingsDpermil")==0) return SurfaceforcingsDpermilEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; - else if (strcmp(name,"SurfaceforcingsAccumulation")==0) return SurfaceforcingsAccumulationEnum; - else if (strcmp(name,"SurfaceforcingsEvaporation")==0) return SurfaceforcingsEvaporationEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum; + if (strcmp(name,"SurfaceforcingsAccumulation")==0) return SurfaceforcingsAccumulationEnum; + else if (strcmp(name,"SurfaceforcingsEvaporation")==0) return SurfaceforcingsEvaporationEnum; + else if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"Penpair")==0) return PenpairEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; - else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; - else if (strcmp(name,"Masscon")==0) return MassconEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MassconName")==0) return MassconNameEnum; + if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; + else if (strcmp(name,"Masscon")==0) return MassconEnum; + else if (strcmp(name,"MassconName")==0) return MassconNameEnum; else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; - else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; - else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; + if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; + else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; + else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; - else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; - else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum; else stage=7; } if(stage==7){ - if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum; + if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; + else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum; + else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum; else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; @@ -873,10 +875,10 @@ else if (strcmp(name,"XY")==0) return XYEnum; else if (strcmp(name,"XYZ")==0) return XYZEnum; - else if (strcmp(name,"Dense")==0) return DenseEnum; - else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + if (strcmp(name,"Dense")==0) return DenseEnum; + else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; + else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; Index: /issm/trunk-jpl/src/m/classes/SMBpdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 19171) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 19172) @@ -12,5 +12,6 @@ s0t = 0; rlaps = 0; - rlapslgm = 0; + rlapslgm = 0; + dpermil = 0; Pfac = NaN; Tdiff = NaN; @@ -18,4 +19,5 @@ isdelta18o = 0; ismungsm = 0; + isd18opd = 0; delta18o = NaN; delta18o_surface = NaN; @@ -35,6 +37,6 @@ end % }}} function self = extrude(self,md) % {{{ - if(self.isdelta18o==0 & self.ismungsm==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end - if(self.isdelta18o==0 & self.ismungsm==0),self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');end + if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end + if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0),self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');end 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 @@ -45,4 +47,6 @@ if(self.ismungsm),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end if(self.ismungsm),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end + if(self.isd18opd),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end + if(self.isd18opd),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end end % }}} @@ -59,4 +63,5 @@ self.isdelta18o = 0; self.ismungsm = 0; + self.isd18opd = 0; self.desfac = 0.5; self.s0p = 0; @@ -64,4 +69,5 @@ self.rlaps = 6.5; self.rlapslgm = 6.5; + self.dpermil = 2.4; end % }}} @@ -74,5 +80,5 @@ md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1); md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1); - if(self.isdelta18o==0 & self.ismungsm==0) + if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0) md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','timeseries',1,'NaN',1); md = checkfield(md,'fieldname','surfaceforcings.precipitation','timeseries',1,'NaN',1); @@ -94,4 +100,9 @@ md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1); md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1); + elseif(self.isd18opd==1) + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','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.delta18o','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.dpermil','>=',0,'numel',1); end end @@ -103,4 +114,5 @@ fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'); fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)'); + fielddisplay(self,'isd18opd','is delta18o parametrisation from present day temperature and precipitation activated (0 or 1, default is 0)'); fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]'); fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'); @@ -108,5 +120,5 @@ fielddisplay(self,'rlaps','present day lapse rate [degree/km]'); fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]'); - if(self.isdelta18o==0 & self.ismungsm==0) + if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0) fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']); fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']); @@ -116,16 +128,21 @@ fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); - fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); - fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated'); + fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm is activated'); fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'); elseif(self.ismungsm==1) - fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); + fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated'); fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); - fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); - fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated'); + fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm is activated'); fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated'); fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'); + elseif(self.isd18opd==1) + fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated'); + fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated'); + fielddisplay(self,'delta18o','delta18o, required if pdd is activated and d18opd activated'); + fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated'); end end % }}} @@ -138,4 +155,5 @@ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean'); WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean'); + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isd18opd','format','Boolean'); WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double'); WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double'); @@ -144,5 +162,5 @@ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double'); - if(self.isdelta18o==0 & self.ismungsm==0) + if(self.isdelta18o==0 & self.ismungsm==0 & self.isd18opd==0) %WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1); WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1); @@ -165,4 +183,9 @@ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); + elseif self.isd18opd + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1); + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','dpermil','format','Double'); end end % }}} Index: /issm/trunk-jpl/src/m/classes/clusters/acenet.m =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/acenet.m (revision 19171) +++ /issm/trunk-jpl/src/m/classes/clusters/acenet.m (revision 19172) @@ -68,16 +68,38 @@ fprintf(fid,'#$ -cwd\n'); fprintf(fid,'#$ -N issm\n'); - %fprintf(fid,'#$ -l h_rt=25:00:0\n'); - %fprintf(fid,'#$ -l h_rt=47:59:00\n'); - %fprintf(fid,'#$ -l h_rt=72:00:0\n'); - fprintf(fid,'#$ -l h_rt=96:00:0\n'); - fprintf(fid,'#$ -l h_vmem=4G\n'); - fprintf(fid,'#$ -pe ompi* %i\n',cluster.np); - fprintf(fid,'#$ -j y\n'); - fprintf(fid,'#$ -l h=cl27*|cl28*|cl29*|cl30*|cl31*|cl320|cl267|cl268|cl269|cl338 \n'); + % fprintf(fid,'#$ -l h_rt=25:00:0\n'); + fprintf(fid,'#$ -l h_rt=47:59:00\n'); + % fprintf(fid,'#$ -l h_rt=72:00:0\n'); + % fprintf(fid,'#$ -l h_rt=96:00:0\n'); + % fprintf(fid,'#$ -l h_rt=336:00:0\n'); + + fprintf(fid,'#$ -l h_vmem=3G\n'); + % if cluster.np>10 + % fprintf(fid,'#$ -l h_vmem=3G\n'); + % else + % fprintf(fid,'#$ -l h_vmem=3G\n'); + % end + + % ---- Which queue to use ---- + %fprintf(fid,'#$ -q !tarasov.q\n'); % + %fprintf(fid,'#$ -q medium.q@*,short.q@*\n'); + fprintf(fid,'#$ -q short.q@*\n'); + + % ---- Which node are selected ---- + % fprintf(fid,'#$ -l h=cl27*|cl28*|cl29*|cl30*|cl31*|cl320|cl267|cl268|cl269|cl338 \n'); + %fprintf(fid,'#$ -l h=cl0* \n'); %fprintf(fid,'#$ -l h=cl338 \n'); - %fprintf(fid,'#$ -pe openmp 20 \n'); - %fprintf(fid,'#$ -q !tarasov.q\n'); % - fprintf(fid,'#$ -pe openmp 8\n'); + % Acenet nodes with 16cpus and more than 60G mem + % fprintf(fid,'#$ -l h=cl001|cl002|cl003|cl004|cl005|cl006|cl007|cl008|cl009|cl010|cl011|cl012|cl021|cl022|cl023|cl024 \n'); + + % ---- cpus on different nodes ---- + fprintf(fid,'#$ -pe ompi %i\n',cluster.np); % To avoid green acenet that does not have InfiniBand + % fprintf(fid,'#$ -pe ompi* %i\n',cluster.np); + % -------- All cpus in the same node -------- + % fprintf(fid,'#$ -pe openmp %i\n',cluster.np); + + % ---- misc ---- + fprintf(fid,'#$ -j y\n'); + fprintf(fid,'module purge\n'); %fprintf(fid,'module load gcc openmpi/gcc\n'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19171) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19172) @@ -346,4 +346,5 @@ def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0] def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0] +def SurfaceforcingsIsd18opdEnum(): return StringToEnum("SurfaceforcingsIsd18opd")[0] def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0] def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0] @@ -359,4 +360,5 @@ def SurfaceforcingsTdiffEnum(): return StringToEnum("SurfaceforcingsTdiff")[0] def SurfaceforcingsSealevEnum(): return StringToEnum("SurfaceforcingsSealev")[0] +def SurfaceforcingsDpermilEnum(): return StringToEnum("SurfaceforcingsDpermil")[0] def SMBgradientsEnum(): return StringToEnum("SMBgradients")[0] def SurfaceforcingsMonthlytemperaturesEnum(): return StringToEnum("SurfaceforcingsMonthlytemperatures")[0]