Index: /issm/trunk-jpl/externalpackages/matlab/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/matlab/install.sh (revision 12747) +++ /issm/trunk-jpl/externalpackages/matlab/install.sh (revision 12748) @@ -5,9 +5,11 @@ #Select or create a new simlink -ln -s /usr/local/pkgs/matlab-7.6/ install +#ln -s /usr/local/pkgs/matlab-7.6/ install #ln -s /discover/vis/mathworks/matlab_r2011b/ install #ln -s /usr/local/matlab704/ install #ln -s /usr/local/matlab711/ install #ln -s /usr/local/matlab712/ install + +ln -s /home/kevin/perso.dir/matlab10 install # Macintosh (OSX) simlink Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12747) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12748) @@ -149,4 +149,10 @@ SettingsResultsAsPatchesEnum, SettingsWaitonlockEnum, + SurfaceforcingsDelta18oEnum, + SurfaceforcingsDelta18oTemperaturesPresentdayEnum, + SurfaceforcingsDelta18oTemperaturesLgmEnum, + SurfaceforcingsDelta18oSurfaceEnum, + SurfaceforcingsIsdelta18oEnum, + SurfaceforcingsPrecipitationsPresentdayEnum, DebugProfilingEnum, ProfilingCurrentMemEnum, Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12747) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12748) @@ -212,4 +212,5 @@ ./shared/Elements/GetNumberOfDofs.cpp\ ./shared/Elements/PddSurfaceMassBalance.cpp\ + ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\ ./shared/String/sharedstring.h\ ./shared/Wrapper/wrappershared.h\ @@ -294,4 +295,6 @@ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.h\ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp\ + ./modules/Delta18oParameterizationx/Delta18oParameterizationx.h\ + ./modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp\ ./modules/SmbGradientsx/SmbGradientsx.h\ ./modules/SmbGradientsx/SmbGradientsx.cpp\ @@ -818,5 +821,5 @@ matlabpetsc_sources= ./matlab/io/MatlabMatrixToPetscMatrix.cpp\ ./matlab/io/MatlabVectorToPetscVector.cpp - + #}}} #Modules sources{{{ Index: /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp (revision 12748) +++ /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp (revision 12748) @@ -0,0 +1,23 @@ +/*!\file Delta18oParameterizationx + * \brief: Create past temperature field from present day fields + */ + +#include "./Delta18oParameterizationx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../io/io.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void Delta18oParameterizationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + + /*Intermediary*/ + int i; + Element* element = NULL; + + /*Compute temperature and precipitation fields: */ + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->Delta18oParameterization(); + } +} Index: /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.h (revision 12748) +++ /issm/trunk-jpl/src/c/modules/Delta18oParameterizationx/Delta18oParameterizationx.h (revision 12748) @@ -0,0 +1,14 @@ +/*!\file: Delta18oParameterizationx.cpp + * \brief header file for + */ + +#ifndef _DELTA18OPARAMETERIZATIONX_H +#define _DELTA18OPARAMETERIZATIONX_H + +#include "../../Container/Container.h" +#include "../../objects/objects.h" + +/* local prototypes: */ +void Delta18oParameterizationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); + +#endif /* _DELTA18OPARAMETERIZATIONX_H*/ Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12748) @@ -154,4 +154,10 @@ case SettingsResultsAsPatchesEnum : return "SettingsResultsAsPatches"; case SettingsWaitonlockEnum : return "SettingsWaitonlock"; + case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o"; + case SurfaceforcingsDelta18oTemperaturesPresentdayEnum : return "SurfaceforcingsDelta18oTemperaturesPresentday"; + case SurfaceforcingsDelta18oTemperaturesLgmEnum : return "SurfaceforcingsDelta18oTemperaturesLgm"; + case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface"; + case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o"; + case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday"; case DebugProfilingEnum : return "DebugProfiling"; case ProfilingCurrentMemEnum : return "ProfilingCurrentMem"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12748) @@ -25,5 +25,8 @@ Parameters *parameters = NULL; IssmDouble *requestedoutputs = NULL; - + bool isdelta18o; + IssmDouble *Delta18oT, *Delta18oSurfaceT; + int nlinesD18o, ncolsD18o, nlinesD18osurf, ncolsD18osurf; + if(*pparameters)return; //do not create parameters twice! @@ -90,5 +93,17 @@ parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIssmbgradientsEnum)); + + iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + if(isdelta18o){ + iomodel->FetchData(&Delta18oT,&nlinesD18o,&ncolsD18o,SurfaceforcingsDelta18oEnum); + _assert_(nlinesD18o==2); + parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oEnum,&Delta18oT[0],&Delta18oT[ncolsD18o],ncolsD18o)); + + iomodel->FetchData(&Delta18oSurfaceT,&nlinesD18osurf,&ncolsD18osurf,SurfaceforcingsDelta18oSurfaceEnum); + _assert_(nlinesD18osurf==2); + parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oSurfaceEnum,&Delta18oSurfaceT[0],&Delta18oSurfaceT[ncolsD18osurf],ncolsD18osurf)); + } /*some parameters that did not come with the iomodel: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12748) @@ -20,5 +20,4 @@ int stabilization; bool dakota_analysis; - bool ispdd; bool issmbgradients; @@ -29,5 +28,4 @@ iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); iomodel->FetchData(1,MeshElementsEnum); - iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); iomodel->Constant(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); @@ -69,22 +67,16 @@ iomodel->FetchDataToInput(elements,TemperatureEnum); } - if(ispdd){ - iomodel->FetchDataToInput(elements,VyEnum); - iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); + if(issmbgradients){ + iomodel->FetchDataToInput(elements,SurfaceforcingsHcEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMaxEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMinEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsAPosEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBPosEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsANegEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum); } - if(issmbgradients){ - iomodel->FetchDataToInput(elements,SurfaceforcingsHcEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMaxEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMinEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsAPosEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsBPosEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsANegEnum); - iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum); - } - else{ - iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); - } + //else{ + // iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + //} /*Free data: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp (revision 12748) @@ -16,5 +16,29 @@ void UpdateElementsTransient(Elements* elements, Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type){ - /*nothing for now: */ + bool ispdd; + bool isdelta18o; + IssmDouble *size, Delta18oTimeSerie,Delta18oSurfaceTimeSerie ; + /*Fetch data needed: */ + iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); + iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + + if(ispdd){ + iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); + if(isdelta18o){ + iomodel->FetchDataToInput(elements,SurfaceforcingsDelta18oEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsDelta18oSurfaceEnum); + + iomodel->FetchDataToInput(elements,SurfaceforcingsDelta18oTemperaturesLgmEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsDelta18oTemperaturesPresentdayEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum); + } + else{ + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); + } + } + else{ + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + } } Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12748) @@ -158,4 +158,10 @@ else if (strcmp(name,"SettingsResultsAsPatches")==0) return SettingsResultsAsPatchesEnum; else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; + else if (strcmp(name,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum; + else if (strcmp(name,"SurfaceforcingsDelta18oTemperaturesPresentday")==0) return SurfaceforcingsDelta18oTemperaturesPresentdayEnum; + else if (strcmp(name,"SurfaceforcingsDelta18oTemperaturesLgm")==0) return SurfaceforcingsDelta18oTemperaturesLgmEnum; + 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,"DebugProfiling")==0) return DebugProfilingEnum; else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; @@ -255,5 +261,8 @@ else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; - else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; @@ -261,8 +270,5 @@ else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"Element")==0) return ElementEnum; + else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ElementResult")==0) return ElementResultEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; @@ -378,5 +384,8 @@ else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; - else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum; else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; @@ -384,8 +393,5 @@ else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; + else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 12747) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 12748) @@ -25,4 +25,5 @@ #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h" #include "./DakotaResponsesx/DakotaResponsesx.h" +#include "./Delta18oParameterizationx/Delta18oParameterizationx.h" #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" #include "./ElementConnectivityx/ElementConnectivityx.h" Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12747) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12748) @@ -72,4 +72,5 @@ virtual void PotentialSheetUngrounding(Vector* potential_sheet_ungrounding)=0; virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0; + virtual void Delta18oParameterization(void)=0; virtual void SmbGradients()=0; virtual int UpdatePotentialSheetUngrounding(IssmDouble* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12748) @@ -683,4 +683,64 @@ this->results=new Results(); +} +/*}}}*/ +/*FUNCTION Penta::Delta18oParameterization{{{*/ +void Penta::Delta18oParameterization(void){ + + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble Delta18oTemperaturesPresentday[NUMVERTICES][12],Delta18oTemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; + IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; + IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; + IssmDouble time,yts,finaltime; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsDelta18oTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsDelta18oTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + GaussPenta* gauss=new GaussPenta(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&Delta18oTemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&Delta18oTemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr + } + } + + /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ + this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime*yts); + this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-21000)*yts); + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time*yts); + this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime*yts); + this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-21000)*yts); + this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time*yts); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivAddTimeInput(newmonthinput1,imonth/12.*yts); + PentaP1Input* newmonthinput2 = new PentaP1Input(SurfaceforcingsPrecipitationEnum,&monthlyprec[iv][imonth]); + NewPrecipitationInput->AddTimeInput(newmonthinput2,imonth/12.*yts); + } + } + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12747) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12748) @@ -84,4 +84,5 @@ void CreatePVector(Vector* pf); void CreateJacobianMatrix(Matrix* Jff); + void Delta18oParameterization(void); void DeleteResults(void); int GetNodeIndex(Node* node); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12748) @@ -906,4 +906,64 @@ this->results=new Results(); +} +/*}}}*/ +/*FUNCTION Tria::Delta18oParameterization{{{*/ +void Tria::Delta18oParameterization(void){ + + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble Delta18oTemperaturesPresentday[NUMVERTICES][12],Delta18oTemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; + IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; + IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; + IssmDouble time,yts,finaltime; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsDelta18oTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsDelta18oTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + GaussTria* gauss=new GaussTria(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&Delta18oTemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&Delta18oTemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr + } + } + + /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ + this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime*yts); + this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-21000)*yts); + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time*yts); + this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime*yts); + this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-21000)*yts); + this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time*yts); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivAddTimeInput(newmonthinput1,imonth/12.*yts); + TriaP1Input* newmonthinput2 = new TriaP1Input(SurfaceforcingsPrecipitationEnum,&monthlyprec[iv][imonth]); + NewPrecipitationInput->AddTimeInput(newmonthinput2,imonth/12.*yts); + } + } + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); } /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12747) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12748) @@ -80,4 +80,5 @@ void CreatePVector(Vector* pf); void CreateJacobianMatrix(Matrix* Jff); + void Delta18oParameterization(void); int GetNodeIndex(Node* node); int Sid(); Index: /issm/trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp (revision 12748) +++ /issm/trunk-jpl/src/c/shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp (revision 12748) @@ -0,0 +1,30 @@ +/* file: ComputeTemperaturePrecipitation.cpp + Compute the temperature and precipitation at time t from + the data at present day and delta18O + */ + +#include "./elements.h" + +void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, + IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, + IssmDouble* PrecipitationsPresentday, + IssmDouble* Delta18oTemperaturesLgm, IssmDouble* Delta18oTemperaturesPresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){ + + IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12]; + IssmDouble delta18oLapseRate=-6.2*pow(10.,-3); + IssmDouble glacialindex; // used to vary present day temperature + + glacialindex = (Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent)) + /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); + + for (int imonth = 0; imonth<12; imonth++){ + monthlytemperaturestmp[imonth] = glacialindex*Delta18oTemperaturesLgm[imonth] + (1-glacialindex)*Delta18oTemperaturesPresentday[imonth]; + monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-0)/10)); //equation from rob's paper + /*Assign output pointer*/ + *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth]; + *(monthlyprecout+imonth) = monthlyprectmp[imonth]; + } + + +} Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 12748) @@ -43,5 +43,5 @@ IssmDouble frzndd = 0.; - IssmDouble tstar; // monthly mean surface temp + IssmDouble tstar; // monthly mean surface temp IssmDouble Tsum= 0.; // average summer (JJA) temperature IssmDouble Tsurf = 0.; // average annual temperature Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12747) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 12748) @@ -13,5 +13,11 @@ IssmDouble Paterson(IssmDouble temperature); IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); -IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water); +IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, + IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water); +void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, + IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, + IssmDouble* PrecipitationsPresentday, + IssmDouble* Delta18oTemperaturesLgm, IssmDouble* Delta18oTemperaturesPresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); void GetVerticesCoordinates(IssmDouble* xyz, Node** nodes, int numvertices); int GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum); Index: /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 12748) @@ -16,5 +16,4 @@ /*parameters: */ bool save_results; - bool ispdd; bool issmbgradients; @@ -24,11 +23,6 @@ /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); - if(ispdd){ - if(VerboseSolution()) _pprintLine_(" call positive degree day module"); - PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - } if(issmbgradients){ Index: /issm/trunk-jpl/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 12747) +++ /issm/trunk-jpl/src/c/solutions/transient_core.cpp (revision 12748) @@ -24,5 +24,5 @@ /*parameters: */ IssmDouble starttime,finaltime,dt,yts; - bool isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy; + bool isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy,ispdd,isdelta18o; bool save_results,dakota_analysis; bool time_adapt=false; @@ -53,4 +53,6 @@ femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum); + femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); + femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); /*initialize: */ @@ -106,4 +108,13 @@ #endif } + + if(ispdd){ + if(isdelta18o){ + if(VerboseSolution()) _pprintLine_(" call Delta18oParametrization module"); + Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(VerboseSolution()) _pprintLine_(" call positive degree day module"); + PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } if(isdiagnostic){ Index: /issm/trunk-jpl/src/m/classes/surfaceforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 12747) +++ /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 12748) @@ -10,4 +10,5 @@ ispdd = 0; issmbgradients = 0; + isdelta18o = 0; hc = NaN; smb_pos_max = NaN; @@ -18,4 +19,9 @@ b_neg = NaN; monthlytemperatures = NaN; + delta18o = NaN; + delta18o_surface = NaN; + temperatures_presentday = NaN; + temperatures_lgm = NaN; + precipitations_presentday = NaN; end methods @@ -33,4 +39,5 @@ obj.ispdd=0; obj.issmbgradients=0; + obj.isdelta18o=0; end % }}} @@ -41,5 +48,14 @@ checkfield(md,'surfaceforcings.issmbgradients','numel',1,'values',[0 1]); if(obj.ispdd) - md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); + if(obj.isdelta18o==0) + md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); + md = checkfield(md,'surfaceforcings.precipitation','forcing',1,'NaN',1); + else + md = checkfield(md,'surfaceforcings.delta18o','NaN',1); + md = checkfield(md,'surfaceforcings.delta18o_surface','NaN',1); + md = checkfield(md,'surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + end elseif(obj.issmbgradients) checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1); @@ -61,8 +77,13 @@ disp(sprintf(' surface forcings parameters:')); - fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]'); fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]'); fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)'); - fielddisplay(obj,'monthlytemperatures','monthly surface temperatures required if pdd is activated'); + fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'); + fielddisplay(obj,'monthlytemperatures','monthly surface temperatures [Kelvin], required if pdd is activated and delta18o not activated'); + fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]'); + fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [Kelvin], required if pdd is activated and delta18o activated'); + fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [Kelvin], 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,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'); fielddisplay(obj,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)'); fielddisplay(obj,'hc',' elevation of intersection between accumulation and ablation regime required if smb gradients is activated'); @@ -79,6 +100,15 @@ WriteData(fid,'object',obj,'fieldname','mass_balance','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','ispdd','format','Boolean'); + WriteData(fid,'object',obj,'fieldname','isdelta18o','format','Boolean'); if obj.ispdd, + if obj.isdelta18o + WriteData(fid,'object',obj,'fieldname','temperatures_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','temperatures_lgm','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','precipitations_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','delta18o_surface','format','DoubleMat','mattype',1); + else WriteData(fid,'object',obj,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1); + end end WriteData(fid,'object',obj,'fieldname','issmbgradients','format','Boolean'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 12747) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 12748) @@ -1377,4 +1377,64 @@ return StringToEnum('SettingsWaitonlock') +def SurfaceforcingsDelta18oEnum(): + """ + SURFACEFORCINGSDELTA18OENUM - Enum of SurfaceforcingsDelta18o + + Usage: + macro=SurfaceforcingsDelta18oEnum() + """ + + return StringToEnum('SurfaceforcingsDelta18o') + +def SurfaceforcingsDelta18oTemperaturesPresentdayEnum(): + """ + SURFACEFORCINGSDELTA18OTEMPERATURESPRESENTDAYENUM - Enum of SurfaceforcingsDelta18oTemperaturesPresentday + + Usage: + macro=SurfaceforcingsDelta18oTemperaturesPresentdayEnum() + """ + + return StringToEnum('SurfaceforcingsDelta18oTemperaturesPresentday') + +def SurfaceforcingsDelta18oTemperaturesLgmEnum(): + """ + SURFACEFORCINGSDELTA18OTEMPERATURESLGMENUM - Enum of SurfaceforcingsDelta18oTemperaturesLgm + + Usage: + macro=SurfaceforcingsDelta18oTemperaturesLgmEnum() + """ + + return StringToEnum('SurfaceforcingsDelta18oTemperaturesLgm') + +def SurfaceforcingsDelta18oSurfaceEnum(): + """ + SURFACEFORCINGSDELTA18OSURFACEENUM - Enum of SurfaceforcingsDelta18oSurface + + Usage: + macro=SurfaceforcingsDelta18oSurfaceEnum() + """ + + return StringToEnum('SurfaceforcingsDelta18oSurface') + +def SurfaceforcingsIsdelta18oEnum(): + """ + SURFACEFORCINGSISDELTA18OENUM - Enum of SurfaceforcingsIsdelta18o + + Usage: + macro=SurfaceforcingsIsdelta18oEnum() + """ + + return StringToEnum('SurfaceforcingsIsdelta18o') + +def SurfaceforcingsPrecipitationsPresentdayEnum(): + """ + SURFACEFORCINGSPRECIPITATIONSPRESENTDAYENUM - Enum of SurfaceforcingsPrecipitationsPresentday + + Usage: + macro=SurfaceforcingsPrecipitationsPresentdayEnum() + """ + + return StringToEnum('SurfaceforcingsPrecipitationsPresentday') + def DebugProfilingEnum(): """ @@ -4505,4 +4565,4 @@ """ - return 449 - + return 455 + Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 12747) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 12748) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=449; +macro=455; Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsDelta18oEnum() +%SURFACEFORCINGSDELTA18OENUM - Enum of SurfaceforcingsDelta18o +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsDelta18oEnum() + +macro=StringToEnum('SurfaceforcingsDelta18o'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oSurfaceEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oSurfaceEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oSurfaceEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsDelta18oSurfaceEnum() +%SURFACEFORCINGSDELTA18OSURFACEENUM - Enum of SurfaceforcingsDelta18oSurface +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsDelta18oSurfaceEnum() + +macro=StringToEnum('SurfaceforcingsDelta18oSurface'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesLgmEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesLgmEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesLgmEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsDelta18oTemperaturesLgmEnum() +%SURFACEFORCINGSDELTA18OTEMPERATURESLGMENUM - Enum of SurfaceforcingsDelta18oTemperaturesLgm +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsDelta18oTemperaturesLgmEnum() + +macro=StringToEnum('SurfaceforcingsDelta18oTemperaturesLgm'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesPresentdayEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesPresentdayEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsDelta18oTemperaturesPresentdayEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsDelta18oTemperaturesPresentdayEnum() +%SURFACEFORCINGSDELTA18OTEMPERATURESPRESENTDAYENUM - Enum of SurfaceforcingsDelta18oTemperaturesPresentday +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsDelta18oTemperaturesPresentdayEnum() + +macro=StringToEnum('SurfaceforcingsDelta18oTemperaturesPresentday'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsdelta18oEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsdelta18oEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsdelta18oEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsIsdelta18oEnum() +%SURFACEFORCINGSISDELTA18OENUM - Enum of SurfaceforcingsIsdelta18o +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsIsdelta18oEnum() + +macro=StringToEnum('SurfaceforcingsIsdelta18o'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsPrecipitationsPresentdayEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsPrecipitationsPresentdayEnum.m (revision 12748) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsPrecipitationsPresentdayEnum.m (revision 12748) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsPrecipitationsPresentdayEnum() +%SURFACEFORCINGSPRECIPITATIONSPRESENTDAYENUM - Enum of SurfaceforcingsPrecipitationsPresentday +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=SurfaceforcingsPrecipitationsPresentdayEnum() + +macro=StringToEnum('SurfaceforcingsPrecipitationsPresentday');