Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19237) @@ -503,4 +503,265 @@ *pdmudB=dmudB; +} +/*}}}*/ +void Element::Delta18oParameterization(void){/*{{{*/ + + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12]; + IssmDouble TemperaturesPresentday[numvertices][12],TemperaturesLgm[numvertices][12]; + IssmDouble PrecipitationsPresentday[numvertices][12]; + IssmDouble tmp[numvertices]; + IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; + IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; + IssmDouble time,yts,finaltime,time_yr; + + /*Recover parameters*/ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + time_yr=floor(time/yts)*yts; + + /*Recover present day temperature and precipitation*/ + Input* input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + 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); + } + } + + /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ + this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime); + this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-(21000*yts))); + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); + this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime); + this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-(21000*yts))); + this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivObjectEnum()){ + case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + for(i=0;iObjectEnum()){ + case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + switch(this->ObjectEnum()){ + case TriaEnum: break; + case PentaEnum: + case TetraEnum: + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Element::MungsmtpParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int numvertices = this->GetNumberOfVertices(); + + 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,time_yr; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*Recover present day temperature and precipitation*/ + Input* input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=this->inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4); + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + 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;ivObjectEnum()){ + case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + for(i=0;iObjectEnum()){ + case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + switch(this->ObjectEnum()){ + case TriaEnum: break; + case PentaEnum: + case TetraEnum: + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +void Element::Delta18opdParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12]; + IssmDouble TemperaturesPresentday[numvertices][12]; + IssmDouble PrecipitationsPresentday[numvertices][12]; + IssmDouble tmp[numvertices]; + IssmDouble Delta18oTime; + IssmDouble dpermil; + IssmDouble time,yts,time_yr,month; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*Get some pdd parameters*/ + dpermil=this->matpar->GetMaterialParameter(SurfaceforcingsDpermilEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2); + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivObjectEnum()){ + case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + for(i=0;iObjectEnum()){ + case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; + default: _error_("Not implemented yet"); + } + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + switch(this->ObjectEnum()){ + case TriaEnum: break; + case PentaEnum: + case TetraEnum: + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + break; + default: _error_("Not implemented yet"); + } + + /*clean-up*/ + delete gauss; } /*}}}*/ @@ -1378,4 +1639,119 @@ ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/ return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +void Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/ + + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble agd[numvertices]; // surface mass balance + IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12]; + IssmDouble yearlytemperatures[numvertices]; memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble)); + IssmDouble tmp[numvertices]; + IssmDouble h[numvertices],s[numvertices]; + IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; + IssmDouble PfacTime,TdiffTime,sealevTime; + IssmDouble mavg=1./12.; //factor for monthly average + + /*Get material parameters :*/ + rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + /*Get some pdd parameters*/ + desfac=this->matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=this->matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=this->matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=this->matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=this->matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); + + /*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/ + Input* input=this->inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); + Input* input2=this->inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); + IssmDouble time,yts,time_yr; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + time_yr=floor(time/yts)*yts; + + /*loop over vertices: */ + Gauss* gauss=this->NewGauss(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&monthlytemperatures[iv][month],gauss,time_yr+month/12.*yts); + yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin + monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module + input2->GetInputValue(&monthlyprec[iv][month],gauss,time_yr+month/12.*yts); + } + } + + /*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); + + /*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); + + switch(this->ObjectEnum()){ + case TriaEnum: + this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); + break; + case PentaEnum: + this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + break; + case TetraEnum: + this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); + this->inputs->AddInput(new TetraInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + break; + default: _error_("Not implemented yet"); + } + // 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*/ + delete gauss; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19237) @@ -66,4 +66,7 @@ void DeleteInput(int input_enum); void DeleteMaterials(void); + void Delta18oParameterization(void); + void MungsmtpParameterization(void); + void Delta18opdParameterization(void); IssmDouble Divergence(void); void dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); @@ -119,9 +122,9 @@ void LinearFloatingiceMeltingRate(); void MigrateGroundingLine(IssmDouble* sheet_ungrounding); - void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses); - + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses); ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); IssmDouble PureIceEnthalpy(IssmDouble pressure); void ResetIceLevelset(void); @@ -183,7 +186,4 @@ virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; 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 Delta18opdParameterization(void)=0; virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; @@ -258,5 +258,4 @@ virtual int NumberofNodesPressure(void)=0; virtual int NumberofNodesVelocity(void)=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/ElementHook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 19237) @@ -77,5 +77,6 @@ int i; - bool* hnodes_null=NULL; /*intermediary needed*/ + bool* hnodesi_null=NULL; /*intermediary needed*/ + bool hnodes_null=true; /*this could be NULL on empty constructor*/ bool hneighbors_null=true; /*don't deal with hneighbors, unless explicitely asked to*/ @@ -85,8 +86,11 @@ if(marshall_direction==MARSHALLING_FORWARD || marshall_direction==MARSHALLING_SIZE){ if(this->hneighbors)hneighbors_null=false; - hnodes_null=xNew(numanalyses); - for(i=0;ihnodes[i])hnodes_null[i]=false; - else hnodes_null[i]=true; + if(this->hnodes){ + hnodes_null=false; + hnodesi_null=xNew(numanalyses); + for(i=0;ihnodes[i])hnodesi_null[i]=false; + else hnodesi_null[i]=true; + } } } @@ -96,22 +100,29 @@ MARSHALLING(numanalyses); MARSHALLING(hneighbors_null); - MARSHALLING_DYNAMIC(hnodes_null,bool,numanalyses); + MARSHALLING(hnodes_null); + MARSHALLING_DYNAMIC(hnodesi_null,bool,numanalyses); if(marshall_direction==MARSHALLING_BACKWARD){ - this->hnodes = new Hook*[numanalyses]; + if (!hnodes_null)this->hnodes = new Hook*[numanalyses]; + else this->hnodes=NULL; this->hvertices = new Hook(); this->hmaterial = new Hook(); this->hmatpar = new Hook(); if(!hneighbors_null)this->hneighbors = new Hook(); + else this->hneighbors=NULL; /*Initialize hnodes: */ - for(int i=0;inumanalyses;i++){ - if(!hnodes_null[i])this->hnodes[i]=new Hook(); - else this->hnodes[i]=NULL; + if (this->hnodes){ + for(int i=0;inumanalyses;i++){ + if(!hnodesi_null[i])this->hnodes[i]=new Hook(); + else this->hnodes[i]=NULL; + } } } - for (i=0;ihnodes[i])this->hnodes[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + if (this->hnodes){ + for (i=0;ihnodes[i])this->hnodes[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + } this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); this->hmatpar->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); @@ -120,5 +131,5 @@ /*Free ressources: */ - xDelete(hnodes_null); + xDelete(hnodesi_null); } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 19237) @@ -580,209 +580,4 @@ }/*}}}*/ -void Penta::Delta18oParameterization(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]; - IssmDouble tmp[NUMVERTICES]; - 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(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _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(&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); - } - } - - /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ - this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime); - this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,finaltime-(21000*yts)); - this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); - this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime); - this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,finaltime-(21000*yts)); - this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,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); - - /*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); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -void Penta::Delta18opdParameterization(void){/*{{{*/ - /*Are we on the base? If not, return*/ - if(!IsOnBase()) return; - - int i; - IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; - IssmDouble tmp[NUMVERTICES]; - IssmDouble Delta18oTime; - IssmDouble dpermil; - IssmDouble time,yts; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - - /*Get some pdd parameters*/ - dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum); - - /*Recover present day temperature and precipitation*/ - Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2); - GaussPenta* gauss=new GaussPenta(); - for(int month=0;month<12;month++) { - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - } - } - - /*Recover interpolation parameters at time t*/ - this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); - - /*Compute the temperature and precipitation*/ - for(int iv=0;ivAddTimeInput(newmonthinput1,time+imonth/12.*yts); - - for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); - } - NewTemperatureInput->Configure(this->parameters); - NewPrecipitationInput->Configure(this->parameters); - - this->inputs->AddInput(NewTemperatureInput); - this->inputs->AddInput(NewPrecipitationInput); - - this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); - this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ void Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/ @@ -2224,101 +2019,4 @@ return PentaEnum; -} -/*}}}*/ -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 yearlytemperatures[NUMVERTICES]={0.}; - IssmDouble tmp[NUMVERTICES]; - IssmDouble h[NUMVERTICES],s[NUMVERTICES]; - IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; - IssmDouble PfacTime,TdiffTime,sealevTime; - IssmDouble mavg=1./12.; //factor for monthly average - - /*Get material parameters :*/ - rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); - rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - /*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); - GaussPenta* gauss=new GaussPenta(); - IssmDouble time,yts; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - - for(int month=0;month<12;month++) { - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts); - yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin - monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module - input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - } - } - - /*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); - - - /*measure the surface mass balance*/ - for (int iv = 0; iv < NUMVERTICES; iv++){ - agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], - pdds,pds, signorm, yts, h[iv], s[iv], - 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(TemperatureEnum,&yearlytemperatures[0],P1Enum)); - this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); - // this->inputs->AddInput(NewTemperatureInput); - // this->inputs->AddInput(NewPrecipitationInput); - // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); - this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); - this->InputExtrude(TemperatureEnum,-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 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 19237) @@ -59,7 +59,4 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); ElementMatrix* CreateBasalMassMatrix(void); - void Delta18oParameterization(void); - void MungsmtpParameterization(void); - void Delta18opdParameterization(void); void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); @@ -141,5 +138,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - 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 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 19237) @@ -54,7 +54,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; 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 Delta18opdParameterization(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; @@ -133,5 +130,4 @@ int NumberofNodesPressure(void){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_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 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 19237) @@ -54,7 +54,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; 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 Delta18opdParameterization(void){_error_("not implemented yet");}; IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; @@ -139,5 +136,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - 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 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19237) @@ -123,4 +123,8 @@ Element::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses); TriaRef::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + + vertices = (Vertex**)this->hvertices->deliverp(); + material = (Material*)this->hmaterial->delivers(); + matpar = (Matpar*)this->hmatpar->delivers(); } @@ -654,199 +658,4 @@ }/*}}}*/ -void Tria::Delta18oParameterization(void){/*{{{*/ - - int i; - IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; - IssmDouble tmp[NUMVERTICES]; - IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; - IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; - IssmDouble time,yts,finaltime; - - /*Recover parameters*/ - 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(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _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(&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); - } - } - - /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ - this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime); - this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-(21000*yts))); - this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); - this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime); - this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-(21000*yts))); - this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,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); - - /*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;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); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ -void Tria::Delta18opdParameterization(void){/*{{{*/ - /*Are we on the base? If not, return*/ - if(!IsOnBase()) return; - - int i; - IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12]; - IssmDouble tmp[NUMVERTICES]; - IssmDouble Delta18oTime; - IssmDouble dpermil; - IssmDouble time,yts; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - - /*Get some pdd parameters*/ - dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum); - - /*Recover present day temperature and precipitation*/ - Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2); - GaussTria* gauss=new GaussTria(); - for(int month=0;month<12;month++) { - for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - } - } - - /*Recover interpolation parameters at time t*/ - this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time); - - /*Compute the temperature and precipitation*/ - for(int iv=0;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); - - /*clean-up*/ - delete gauss; -} -/*}}}*/ int Tria::EdgeOnBaseIndex(void){/*{{{*/ @@ -2638,8 +2447,4 @@ } /*}}}*/ -int Tria::PressureInterpolation(void){/*{{{*/ - return TriaRef::PressureInterpolation(this->element_type); -} -/*}}}*/ int Tria::NumberofNodesPressure(void){/*{{{*/ return TriaRef::NumberofNodes(this->PressureInterpolation()); @@ -2648,99 +2453,4 @@ int Tria::NumberofNodesVelocity(void){/*{{{*/ return TriaRef::NumberofNodes(this->VelocityInterpolation()); -} -/*}}}*/ -void Tria::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 yearlytemperatures[NUMVERTICES]={0.}; - IssmDouble tmp[NUMVERTICES]; - IssmDouble h[NUMVERTICES],s[NUMVERTICES]; - IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; - IssmDouble PfacTime,TdiffTime,sealevTime; - IssmDouble mavg=1./12.; //factor for monthly average - - /*Get material parameters :*/ - rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); - rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - - /*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 compute the yearly mean temperatures*/ - Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); - Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); - 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;ivGaussVertex(iv); - input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts); - yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin - monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module - input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - } - } - - /*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); - - /*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(TemperatureEnum,&yearlytemperatures[0],P1Enum)); - 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*/ - delete gauss; } /*}}}*/ @@ -2770,4 +2480,8 @@ } } +} +/*}}}*/ +int Tria::PressureInterpolation(void){/*{{{*/ + return TriaRef::PressureInterpolation(this->element_type); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19236) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19237) @@ -64,7 +64,4 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); - void Delta18oParameterization(void); - void MungsmtpParameterization(void); - void Delta18opdParameterization(void); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); @@ -114,5 +111,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation();