Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 18967) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 18968) @@ -220,4 +220,5 @@ ./shared/Elements/PddSurfaceMassBalance.cpp\ ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\ + ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\ ./shared/Elements/DrainageFunctionWaterfraction.cpp\ ./shared/String/sharedstring.h\ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18968) @@ -120,5 +120,5 @@ int stabilization,finiteelement,smb_model; bool dakota_analysis; - bool isdelta18o; + bool isdelta18o,ismungsm; bool isgroundingline; bool islevelset; @@ -179,6 +179,7 @@ case SMBpddEnum: iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); - if(isdelta18o){ + if(isdelta18o || ismungsm){ iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum); @@ -187,7 +188,8 @@ } else{ - iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); } + break; case SMBgradientsEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18968) @@ -180,4 +180,5 @@ virtual void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum)=0; virtual void Delta18oParameterization(void)=0; + virtual void MungsmtpParameterization(void)=0; virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; @@ -251,5 +252,5 @@ virtual int NumberofNodesPressure(void)=0; virtual int NumberofNodesVelocity(void)=0; - virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0; + virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18968) @@ -608,5 +608,4 @@ input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec } } @@ -647,4 +646,70 @@ this->inputs->AddInput(NewPrecipitationInput); + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + + /*clean-up*/ + delete gauss; +} +void Penta::MungsmtpParameterization(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],TemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble TdiffTime,PfacTime; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4); + 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(&TemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,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); @@ -2082,30 +2147,28 @@ } /*}}}*/ -void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/ - +void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/ + + int i; IssmDouble agd[NUMVERTICES]; // surface mass balance IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; IssmDouble h[NUMVERTICES],s[NUMVERTICES]; IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; IssmDouble PfacTime,TdiffTime,sealevTime; - IssmDouble sconv; //rhow_rain/rhoi / 12 months /*Get material parameters :*/ rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - //rho_ice=matpar->GetRhoIce(); - //rho_water=matpar->GetRhoFreshwater(); - - sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month + + /*Get some pdd parameters*/ + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*Recover monthly temperatures and precipitation*/ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); - Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3); - Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4); - Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5); - Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6); GaussPenta* gauss=new GaussPenta(); IssmDouble time,yts; @@ -2119,20 +2182,18 @@ monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - monthlyprec[iv][month]=monthlyprec[iv][month]*sconv; // convertion to m/yr - input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius - input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); - TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius - input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; - input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); - PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month]; } } - /*Recover Pfac, Tdiff and sealev at time t:*/ - this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); - this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); - this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + /*Recover Pfac, Tdiff and sealev at time t: + This parameters are used to interpolate the temperature + and precipitaton between PD and LGM when ismungsm==1 */ + if (ismungsm==1){ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + } + else { + TdiffTime=0; + sealevTime=0; + } /*Recover info at the vertices: */ @@ -2140,27 +2201,40 @@ GetInputListOnVertices(&s[0],SurfaceEnum); - /*Get other pdd parameters*/ - desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); - s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); - s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); - rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); - rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ for (int iv = 0; iv < NUMVERTICES; iv++){ - agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], - &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], - &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], + agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds,pds, signorm, yts, h[iv], s[iv], - desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime); + desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime, + rho_water,rho_ice); } /*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(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); - //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); - this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); - - /*clean-up*/ - delete gauss; + // this->inputs->AddInput(NewTemperatureInput); + // this->inputs->AddInput(NewPrecipitationInput); + // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); + this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); + // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + + /*clean-up*/ + delete gauss; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18968) @@ -59,4 +59,5 @@ ElementMatrix* CreateBasalMassMatrix(void); void Delta18oParameterization(void); + void MungsmtpParameterization(void); void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); @@ -137,5 +138,5 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18968) @@ -54,4 +54,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; + void MungsmtpParameterization(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");}; @@ -129,5 +130,5 @@ int NumberofNodesPressure(void){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18968) @@ -54,4 +54,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; + void MungsmtpParameterization(void){_error_("not implemented yet");}; IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; @@ -135,5 +136,5 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18968) @@ -609,5 +609,4 @@ input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec } } @@ -626,4 +625,67 @@ Delta18oPresent, Delta18oLgm, Delta18oTime, &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::MungsmtpParameterization(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],TemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble TdiffTime,PfacTime; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4); + 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(&TemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivGetMaterialParameter(MaterialsRhoSeawaterEnum); rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - //rho_ice=matpar->GetRhoIce(); - //rho_water=matpar->GetRhoFreshwater(); - - sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month + + /*Get some pdd parameters*/ + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*Recover monthly temperatures and precipitation and Present day and LGm ones*/ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); - Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3); - Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4); - Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5); - Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6); GaussTria* gauss=new GaussTria(); IssmDouble time,yts; this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); + + for(int month=0;month<12;month++) { for(int iv=0;ivGetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - monthlyprec[iv][month]=monthlyprec[iv][month]*sconv*yts; // convertion in m/y - input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius - input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); - TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius - input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; - input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); - PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month]; } } - /*Recover Pfac, Tdiff and sealev at time t:*/ - this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); - this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); - this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + /*Recover Pfac, Tdiff and sealev at time t: + This parameters are used to interpolate the temperature + and precipitaton between PD and LGM when ismungsm==1 */ + if (ismungsm==1){ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + } + else { + TdiffTime=0; + sealevTime=0; + } /*Recover info at the vertices: */ GetInputListOnVertices(&h[0],ThicknessEnum); GetInputListOnVertices(&s[0],SurfaceEnum); - - /*Get other pdd parameters*/ - desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); - s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); - s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); - rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); - rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ - 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(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); + // this->inputs->AddInput(NewTemperatureInput); + // this->inputs->AddInput(NewPrecipitationInput); // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); + + //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); + // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); /*clean-up*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18967) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18968) @@ -63,4 +63,5 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); void Delta18oParameterization(void); + void MungsmtpParameterization(void); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -109,5 +110,5 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18968) @@ -22,5 +22,5 @@ char** requestedoutputs = NULL; IssmDouble time; - bool isdelta18o; + bool isdelta18o,ismungsm; /*parameters for mass flux:*/ @@ -108,4 +108,5 @@ case SMBpddEnum: parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum)); @@ -114,19 +115,22 @@ parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapslgmEnum)); iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); - - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsPfacEnum); + iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); + + if(ismungsm){ + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsPfacEnum); - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum); - - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsSealevEnum); + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum); + + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsSealevEnum); + } if(isdelta18o){ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18968) @@ -11,5 +11,5 @@ /*Intermediaties*/ int smb_model; - bool isdelta18o; + bool isdelta18o,ismungsm; /*First, get SMB model from parameters*/ @@ -23,7 +23,12 @@ case SMBpddEnum: femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum); if(isdelta18o){ - if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n"); + if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n"); Delta18oParameterizationx(femmodel); + } + if(ismungsm){ + if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n"); + MungsmtpParameterizationx(femmodel); } if(VerboseSolution()) _printf0_(" call positive degree day module\n"); @@ -115,4 +120,12 @@ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); element->Delta18oParameterization(); + } + +}/*}}}*/ +void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->MungsmtpParameterization(); } @@ -144,4 +157,6 @@ IssmDouble tstar; // monthly mean surface temp + bool ismungsm; + IssmDouble *pdds = NULL; IssmDouble *pds = NULL; @@ -150,4 +165,7 @@ pdds=xNew(NPDMAX+1); pds=xNew(NPDCMAX+1); + + // Get ismungsm parameter + femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum); /* initialize PDD (creation of a lookup table)*/ @@ -205,7 +223,6 @@ for(i=0;ielements->Size();i++){ element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->PositiveDegreeDay(pdds,pds,signorm); - } - + element->PositiveDegreeDay(pdds,pds,signorm,ismungsm); + } /*free ressouces: */ xDelete(pdds); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18967) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18968) @@ -12,4 +12,5 @@ void SmbGradientsx(FemModel* femmodel); void Delta18oParameterizationx(FemModel* femmodel); +void MungsmtpParameterizationx(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 18967) +++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 18968) @@ -1,4 +1,5 @@ /* file: PddSurfaceMassBlance.cpp Calculating the surface mass balance using the positive degree day method. + Updating the precipitation and temperature to the new elevation */ @@ -6,15 +7,10 @@ #include "../Numerics/numerics.h" -IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday, - IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, - IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, - IssmDouble s0t, - IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, - IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ - - // CURENTLY monthlytemperatures and monthlyprec ARE NOT USED HERE. - // THEY ARE REPLACE BY tdiffh No usage of deltO18 anymore +IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, + IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, + IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, + IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm, + IssmDouble TdiffTime,IssmDouble sealevTime, + IssmDouble rho_water,IssmDouble rho_ice){ // output: @@ -28,13 +24,13 @@ IssmDouble water; //water=rain + snowmelt IssmDouble runoff; //meltwater only, does not include rain + IssmDouble sconv; //rhow_rain/rhoi / 12 months //IssmDouble sealev=0.; // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper //IssmDouble Pfac=0.5,Tdiff=0.5; - IssmDouble rtlaps; + IssmDouble rtlaps; // IssmDouble lapser=6.5 // lapse rate // IssmDouble desfac = 0.3; // desert elevation factor // IssmDouble s0p=0.; // should be set to elevation from precip source // IssmDouble s0t=0.; // should be set to elevation from temperature source - IssmDouble tdiffh; IssmDouble st; // elevation between altitude of the temp record and current altitude IssmDouble sp; // elevation between altitude of the prec record and current altitude @@ -49,5 +45,4 @@ IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - IssmDouble premt; // monthly precipitation (m of ice equivalent) IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. IssmDouble qm = 0.; // snow part of the precipitation @@ -73,4 +68,5 @@ IssmDouble pddtj, hmx2; + sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months /*PDD constant*/ @@ -90,6 +86,6 @@ /*********compute Surface temperature *******/ - tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] ); - tstar = tdiffh + TemperaturesPresentday[imonth] - rtlaps *max(st,sealevTime*0.001); + monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001); + tstar = monthlytemperatures[imonth]; Tsurf = tstar*deltm+Tsurf; @@ -106,8 +102,6 @@ else {q = 1.0;} - premt =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime)); // [m/month] - - qmt= qmt + premt; - qmpt= q*premt; + 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; @@ -124,4 +118,8 @@ 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 //****************************************************************** Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 18967) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 18968) @@ -13,16 +13,19 @@ IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat); // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); -IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday, - IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts, - IssmDouble h, IssmDouble s, IssmDouble desfac, - IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, - IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime); +IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, + IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts, + IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t, + IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, + IssmDouble TdiffTime,IssmDouble sealevTime, + IssmDouble rho_water, IssmDouble rho_ice); void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, - IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, - IssmDouble* PrecipitationsPresentday, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); + IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, + IssmDouble* PrecipitationsPresentday, + IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); +void ComputeMungsmTemperaturePrecipitation(IssmDouble TdiffTime, IssmDouble PfacTime, + IssmDouble* PrecipitationsLgm,IssmDouble* PrecipitationsPresentday, + IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, + 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 18967) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18968) @@ -344,4 +344,5 @@ SurfaceforcingsDelta18oSurfaceEnum, SurfaceforcingsIsdelta18oEnum, + SurfaceforcingsIsmungsmEnum, SurfaceforcingsPrecipitationsPresentdayEnum, SurfaceforcingsPrecipitationsLgmEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18968) @@ -350,4 +350,5 @@ case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface"; case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o"; + case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm"; case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday"; case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18967) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18968) @@ -356,4 +356,5 @@ else if (strcmp(name,"SurfaceforcingsDelta18oSurface")==0) return SurfaceforcingsDelta18oSurfaceEnum; else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum; + else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum; - else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; + if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; + else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; @@ -505,9 +506,9 @@ 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 stage=5; } if(stage==5){ - if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; + if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; + else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; @@ -628,9 +629,9 @@ 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 stage=6; } if(stage==6){ - if (strcmp(name,"StrainRate")==0) return StrainRateEnum; + if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; + else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Gset")==0) return GsetEnum; else if (strcmp(name,"Index")==0) return IndexEnum; - else if (strcmp(name,"Indexed")==0) return IndexedEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"Intersect")==0) return IntersectEnum; + if (strcmp(name,"Indexed")==0) return IndexedEnum; + else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; Index: /issm/trunk-jpl/src/m/classes/SMBpdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 18967) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 18968) @@ -17,4 +17,5 @@ sealev = NaN; isdelta18o = 0; + ismungsm = 0; delta18o = NaN; delta18o_surface = NaN; @@ -34,20 +35,22 @@ end % }}} function self = extrude(self,md) % {{{ - - self.precipitation=project3d(md,'vector',self.precipitation,'type','node'); - self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node'); + 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),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end if(self.isdelta18o),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end if(self.isdelta18o),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end if(self.isdelta18o),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end - + if(self.ismungsm),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end + if(self.ismungsm),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end + 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 end % }}} function self = initialize(self,md) % {{{ - - if isnan(self.precipitation), - self.precipitation=zeros(md.mesh.numberofvertices,1); - disp(' no SMBpdd.precipitation specified: values set as zero'); - end + + % if isnan(self.precipitation), + % self.precipitation=zeros(md.mesh.numberofvertices,1); + % disp(' no SMBpdd.precipitation specified: values set as zero'); + % end end % }}} @@ -55,4 +58,5 @@ obj.isdelta18o = 0; + obj.ismungsm = 0; obj.desfac = 0.5; obj.s0p = 0; @@ -70,19 +74,26 @@ md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1); md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1); - md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]); - md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,NaN]); - md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,NaN]); - if(obj.isdelta18o==0) - md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1); - else - md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + if(obj.isdelta18o==0 & obj.ismungsm==0) + md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1); + elseif(obj.isdelta18o==1) + md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1); + elseif(obj.ismungsm==1) + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]); + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1); end - end + end end % }}} function disp(obj) % {{{ @@ -91,4 +102,5 @@ disp(sprintf('\n PDD and deltaO18 parameters:')); fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'); + fielddisplay(obj,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)'); fielddisplay(obj,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]'); fielddisplay(obj,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'); @@ -96,16 +108,25 @@ fielddisplay(obj,'rlaps','present day lapse rate [degree/km]'); fielddisplay(obj,'rlapslgm','LGM lapse rate [degree/km]'); - fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year)'); - fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year)'); - fielddisplay(obj,'sealev','sea level [m], 1D(year)'); - fielddisplay(obj,'monthlytemperatures',['CURRENTLY NOT USED monthly surface temperatures [K], required if pdd is activated and delta18o not activated']); - fielddisplay(obj,'precipitation',['CURRENTLY NOT USED monthly surface precipitation [m/yr water eq]']); - fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated'); - fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'); - + if(obj.isdelta18o==0 & obj.ismungsm==0) + fielddisplay(obj,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']); + fielddisplay(obj,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']); + elseif(obj.isdelta18o==1) + fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated'); + fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'); + fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); + fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated'); + elseif(obj.ismungsm==1) + fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); + fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated'); + end end % }}} function marshall(obj,md,fid) % {{{ @@ -115,5 +136,6 @@ WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer'); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean'); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','ismungsm','format','Boolean'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','desfac','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0p','format','Double'); @@ -121,9 +143,10 @@ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlaps','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlapslgm','format','Double'); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean'); - if obj.isdelta18o + + if(obj.isdelta18o==0 & obj.ismungsm==0) + %WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1); + elseif obj.isdelta18o WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1); @@ -132,6 +155,14 @@ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1); - else - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); + elseif obj.ismungsm + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); end end % }}} Index: /issm/trunk-jpl/src/m/classes/SMBpdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 18967) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 18968) @@ -26,4 +26,5 @@ self.sealev = float('NaN') self.isdelta18o = 0 + self.ismungsm = 0 self.delta18o = float('NaN') self.delta18o_surface = float('NaN') @@ -31,4 +32,5 @@ self.temperatures_lgm = float('NaN') self.precipitations_presentday = float('NaN') + self.precipitations_lgm = float('NaN') #set defaults @@ -38,39 +40,55 @@ string=" surface forcings parameters:" - string="%s\n%s"%(string,fielddisplay(self,'precipitation','surface precipitation [m/yr water eq]')) + string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)')) + string="%s\n%s"%(string,fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]')) - string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]')) string="%s\n%s"%(string,fielddisplay(self,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]')) string="%s\n%s"%(string,fielddisplay(self,'rlaps','present day lapse rate [degree/km]')) string="%s\n%s"%(string,fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]')) - string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D (year)')) - string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D (year)')) - string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year)')) - string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated')) - string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated')) - + if not (self.isdelta18o and self.ismungsm): + string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated'])) + string="%s\n%s"%(string,fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated'])) + if self.isdelta18o: + string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated')) + string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated')) + if self.ismungsm: + string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated')) return string #}}} def extrude(self,md): # {{{ - self.precipitation=project3d(md,'vector',self.precipitation,'type','node'); - self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node'); + if not (self.isdelta18o and self.ismungsm): + self.precipitation=project3d(md,'vector',self.precipitation,'type','node') + self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node') if self.isdelta18o: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node') if self.isdelta18o: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node') if self.isdelta18o: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node') + if self.isdelta18o: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node') + if self.ismungsm: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node') + if self.ismungsm: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node') + if self.ismungsm: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node') + if self.ismungsm: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node') return self #}}} def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.precipitation)): - self.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) - print " no SMBpdd.precipitation specified: values set as zero" - - return self + # if numpy.all(numpy.isnan(self.precipitation)): + # self.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) + # print " no SMBpdd.precipitation specified: values set as zero" + # + # return self #}}} def setdefaultparameters(self): # {{{ @@ -78,4 +96,5 @@ #pdd method not used in default mode self.isdelta18o = 0 + self.ismungsm = 0 self.desfac = 0.5 self.s0p = 0. @@ -89,16 +108,14 @@ if MasstransportAnalysisEnum() in analyses: - md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]); - md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]); - md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]); - if not self.isdelta18o: - md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1) - md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1) - else: + md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]) + + if not (self.isdelta18o and self.ismungsm): + md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitation','NaN',1) + elif self.isdelta18o: md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1) md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1) @@ -106,4 +123,15 @@ md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1) md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]) + elif self.ismungsm: + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]) return md @@ -113,23 +141,33 @@ yts=365.0*24.0*3600.0 - WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer'); + WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer') - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean') - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); - if self.isdelta18o: + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double') + + if not (self.isdelta18o and self.ismungsm): + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1) + elif self.isdelta18o: WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1) - else: - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1) + elif self.ismungsm: + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1) # }}} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18967) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18968) @@ -342,4 +342,5 @@ def SurfaceforcingsDelta18oSurfaceEnum(): return StringToEnum("SurfaceforcingsDelta18oSurface")[0] def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0] +def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0] def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0] def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0] Index: /issm/trunk-jpl/test/NightlyRun/test236.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test236.m (revision 18967) +++ /issm/trunk-jpl/test/NightlyRun/test236.m (revision 18968) @@ -3,10 +3,16 @@ md=parameterize(md,'../Par/SquareShelf.par'); +%md.verbose=verbose('all'); + % Use of ispdd and isdelta18o methods md.surfaceforcings = SMBpdd(); md.surfaceforcings.isdelta18o=1; +md.surfaceforcings.ismungsm=0; + +%md.surfaceforcings.precipitation(1:md.mesh.numberofvertices,1:12)=0; +%md.surfaceforcings.monthlytemperatures(1:md.mesh.numberofvertices,1:12)=273; % Add temperature, precipitation and delta18o needed to measure the surface mass balance -% creating delta18o +% creating delta18o load '../Data/delta18o.data' md.surfaceforcings.delta18o=delta18o; @@ -26,5 +32,6 @@ end -% creating initialization and spc temperatures initialization and spc +% creating initialization and spc temperatures initialization and +% spc md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1); md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step); @@ -37,19 +44,26 @@ for imonth=0:11 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + % Time for the last line: md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); end + +% Interpolation factors +md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5; +% Year of each data point +md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; % time steps and resolution md.timestepping.time_step=20; +md.settings.output_frequency=1; md.timestepping.final_time=60; - -md.surfaceforcings.Pfac=[1;1]; -md.surfaceforcings.Tdiff=[1;1]; -md.surfaceforcings.sealev=[1;1]; % md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'}; md=setflowequation(md,'SSA','all'); -md.cluster=generic('name',oshostname(),'np',3); +md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster md=solve(md,TransientSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test236.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test236.py (revision 18967) +++ /issm/trunk-jpl/test/NightlyRun/test236.py (revision 18968) @@ -16,4 +16,5 @@ md.surfaceforcings = SMBpdd(); md.surfaceforcings.isdelta18o=1 +md.surfaceforcings.ismungsm=0 # Add temperature, precipitation and delta18o needed to measure the surface mass balance @@ -50,4 +51,13 @@ md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 + md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + +# Interpolation factors +md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5; +# Year of each data point +md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; # time steps and resolution @@ -55,7 +65,4 @@ md.timestepping.final_time=60. -md.surfaceforcings.Pfac=[1,1] -md.surfaceforcings.Tdiff=[1,1] -md.surfaceforcings.sealev=[1,1] # Index: /issm/trunk-jpl/test/NightlyRun/test237.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test237.m (revision 18967) +++ /issm/trunk-jpl/test/NightlyRun/test237.m (revision 18968) @@ -1,18 +1,26 @@ md=triangle(model(),'../Exp/Square.exp',600000.);%180000 md=setmask(md,'all',''); - -% Use of ispdd and isdelta18o methods -md.surfaceforcings = SMBpdd(); -md.surfaceforcings.isdelta18o=1; - md=parameterize(md,'../Par/SquareShelf.par'); +%md.verbose=verbose('all'); + +% Use of ispdd methods +md.surfaceforcings = SMBpdd(); +md.surfaceforcings.isdelta18o=0; +md.surfaceforcings.ismungsm=1; + +if md.surfaceforcings.isdelta18o==0 & md.surfaceforcings.ismungsm==0 + md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); + md.surfaceforcings.monthlytemperatures=273*ones(md.mesh.numberofvertices,1); +end + + % Add temperature, precipitation and delta18o needed to measure the surface mass balance -% creating delta18o -load '../Data/delta18o.data' -md.surfaceforcings.delta18o=delta18o; -% creating delta18oSurface -md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; -md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); +% % creating delta18o +% load '../Data/delta18o.data' +% md.surfaceforcings.delta18o=delta18o; +% % creating delta18oSurface +% md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; +% md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); % creating Present day and lgm temperatures @@ -38,10 +46,17 @@ for imonth=0:11 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + % Time for the last line: md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); end -md.surfaceforcings.Pfac=[1;1]; -md.surfaceforcings.Tdiff=[1;1]; -md.surfaceforcings.sealev=[1;1]; +md.surfaceforcings.Pfac(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5; +% Year of each data point +md.surfaceforcings.Pfac(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; % time steps and resolution @@ -52,8 +67,9 @@ % md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'}; -md=extrude(md,3,1.); +md=extrude(md,3,1); + md=setflowequation(md,'SSA','all'); -md.cluster=generic('name',oshostname(),'np',1); -md=solve(md,TransientSolutionEnum()); +md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster +md=solve(md,TransientSolutionEnum); %Fields and tolerances to track changes Index: /issm/trunk-jpl/test/NightlyRun/test237.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test237.py (revision 18967) +++ /issm/trunk-jpl/test/NightlyRun/test237.py (revision 18968) @@ -12,18 +12,19 @@ md=triangle(model(),'../Exp/Square.exp',600000) #180000 md=setmask(md,'all','') +md=parameterize(md,'../Par/SquareShelf.py') # Use of ispdd and isdelta18o methods md.surfaceforcings = SMBpdd(); -md.surfaceforcings.isdelta18o=1 +md.surfaceforcings.isdelta18o=0 +md.surfaceforcings.ismungsm=1 -md=parameterize(md,'../Par/SquareShelf.py') -# Add temperature, precipitation and delta18o needed to measure the surface mass balance -# creating delta18o -delta18o=numpy.loadtxt('../Data/delta18o.data') -md.surfaceforcings.delta18o=delta18o -# creating delta18oSurface -md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1))) -md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:] +# # Add temperature, precipitation and delta18o needed to measure the surface mass balance +# # creating delta18o +# delta18o=numpy.loadtxt('../Data/delta18o.data') +# md.surfaceforcings.delta18o=delta18o +# # creating delta18oSurface +# md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1))) +# md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:] # creating Present day and lgm temperatures @@ -52,4 +53,15 @@ md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 + md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + +# Interpolation factors +md.surfaceforcings.Pfac[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5; +# Year of each data point +md.surfaceforcings.Pfac[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; # time steps and resolution @@ -57,8 +69,4 @@ md.settings.output_frequency=1 md.timestepping.final_time=60. - -md.surfaceforcings.Pfac=[1,1] -md.surfaceforcings.Tdiff=[1,1] -md.surfaceforcings.sealev=[1,1] #