Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 22624) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 22625) @@ -223,4 +223,5 @@ ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp\ ./modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp\ + ./modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp\ ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\ ./modules/SpcNodesx/SpcNodesx.cpp\ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22625) @@ -2341,29 +2341,4 @@ } /*}}}*/ -void Element::PicoFloatingiceMeltingRate(){/*{{{*/ - - int numvertices = this->GetNumberOfVertices(); - IssmDouble deepwaterel,upperwaterel,deepwatermelt; - IssmDouble* base = xNew(numvertices); - IssmDouble* values = xNew(numvertices); - IssmDouble time; - - parameters->FindParam(&time,TimeEnum); - parameters->FindParam(&deepwaterel,BasalforcingsDeepwaterElevationEnum,time); - parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum,time); - parameters->FindParam(&deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum,time); - - this->GetInputListOnVertices(base,BaseEnum); - for(int i=0;iupperwaterel) values[i]=0; - else if (base[i]AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); - xDelete(base); - xDelete(values); - -}/*}}}*/ void Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22625) @@ -142,5 +142,4 @@ ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); - void PicoFloatingiceMeltingRate(); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); IssmDouble PureIceEnthalpy(IssmDouble pressure); @@ -213,4 +212,5 @@ virtual Element* GetBasalElement(void)=0; virtual int GetElementType(void)=0; + virtual IssmDouble GetHorizontalSurfaceArea(void){_error_("not implemented");}; virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0; virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0; @@ -278,5 +278,6 @@ virtual int NumberofNodesVelocity(void)=0; virtual void PicoUpdateBoxid(int* pmax_boxid_basin)=0; - virtual void PicoUpdateFirstBox(void){_error_("not implemented");}; + virtual void PicoUpdateFirstBox()=0; + virtual void PicoUpdateNextBox(int loopboxid)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22625) @@ -140,4 +140,6 @@ int NumberofNodesVelocity(void); void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; + void PicoUpdateFirstBox(){_error_("not implemented yet");}; + void PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");}; 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 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22625) @@ -129,4 +129,6 @@ int NumberofNodesVelocity(void){_error_("not implemented yet");}; void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; + void PicoUpdateFirstBox(){_error_("not implemented yet");}; + void PicoUpdateNextBox(int loopboxid){_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 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22625) @@ -137,4 +137,6 @@ int NumberofNodesVelocity(void); void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; + void PicoUpdateFirstBox(){_error_("not implemented yet");}; + void PicoUpdateNextBox(int loopboxid){_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 22624) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22625) @@ -1107,4 +1107,9 @@ _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0); return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2; +} +/*}}}*/ +IssmDouble Tria::GetHorizontalSurfaceArea(void){/*{{{*/ + + return this->GetArea(); } /*}}}*/ @@ -2776,32 +2781,36 @@ if(!this->IsIceInElement() || !this->IsFloating()) return; - //Load inputs and params - int i, boxid; - int basin_id, boxid_max; - IssmDouble dist_gl; - IssmDouble dist_cf; - IssmDouble rel_dist_gl; + int i,boxid,basin_id,boxid_max; + IssmDouble dist_gl,dist_cf,rel_dist_gl,lowbound,highbound; inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum); - boxid_max=pmax_boxid_basin[basin_id]; //0-based - - Input* dist_gl_input=NULL; - Input* dist_cf_input=NULL; - dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input); - dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); - + boxid_max=pmax_boxid_basin[basin_id]; + + Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input); + Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); + + /*Get dist_gl and dist_cf at center of element*/ Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); dist_gl_input->GetInputValue(&dist_gl,gauss); dist_cf_input->GetInputValue(&dist_cf,gauss); - + + /*Ensure values are positive for floating ice*/ + if(dist_gl<0){dist_gl=dist_gl*(-1);} + if(dist_cf<0){dist_cf=dist_cf*(-1);} + + /*Compute relative distance to grounding line*/ rel_dist_gl=dist_gl/(dist_gl+dist_cf); + /*Assign box numbers based on rel_dist_gl*/ boxid=-1; for(i=0;i=(1-sqrt((boxid_max-i)/boxid_max)) && rel_dist_gl<=(1-sqrt((boxid_max-i-1)/boxid_max))){ - boxid=i; break; - } - } - if(boxid==-1){_error_("no boxid found for element" << this->Sid());} + lowbound = 1-sqrt((boxid_max-i-0.0)/boxid_max); + highbound = 1-sqrt((boxid_max-i-1.0)/boxid_max); + if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){ + boxid=i; + break; + } + } + if(boxid==-1){_error_("No boxid found for element " << this->Sid() << "!");} this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); @@ -2813,58 +2822,177 @@ void Tria::PicoUpdateFirstBox(){/*{{{*/ - IssmDouble rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c; - IssmDouble alpha, Beta, gamma_T, overturning_coeff, t_farocean, s_farocean; - IssmDouble basinid, boxid, thickness, toc_farocean, soc_farocean, area_box1; - IssmDouble pressure, T_star, g1, s1, p_coeff, q_coeff, Toc, Soc, potential_pressure_melting_point; - IssmDouble basalmeltrate_shelf, overturning, maxbox; - - //Get variables + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int boxid; + this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); + if(!boxid==0) return; + + int basinid, maxbox, num_basins, numnodes, M; + IssmDouble time, rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c; + IssmDouble alpha, Beta, gamma_T, overturning_coeff; + IssmDouble thickness, toc_farocean, soc_farocean, area_box1; + IssmDouble pressure, T_star, g1, s1, p_coeff, q_coeff; + IssmDouble* boxareas=NULL; + IssmDouble* t_farocean=NULL; + IssmDouble* s_farocean=NULL; + + /*Get variables*/ rhoi = this->GetMaterialParameter(MaterialsRhoIceEnum); rhow = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); earth_grav = this->GetMaterialParameter(ConstantsGEnum); - rho_star = 1033.; //kg/m^3 + rho_star = 1033.; // kg/m^3 nu = rhoi/rhow; latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum); - c_p_ocean = this->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum); + c_p_ocean = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); lambda = latentheat/c_p_ocean; - a = -0.0572; //K/psu + a = -0.0572; // K/psu b = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum); //K - c = this->GetMaterialParameter(MaterialsBetaEnum); - alpha = 0.000075; //1/K - Beta = 0.00077; //K - + c = 7.77e-4; + alpha = 7.5e-5; // 1/K + Beta = 7.7e-4; // K + + /* Get parameters and inputs */ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum); this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); - this->parameters->FindParam(&t_farocean,BasalforcingsPicoFarOceansalinityEnum); - this->parameters->FindParam(&s_farocean,BasalforcingsPicoFarOceansalinityEnum); - + this->parameters->FindParam(&t_farocean, &num_basins, time, BasalforcingsPicoFarOceantemperatureEnum); + this->parameters->FindParam(&s_farocean, &num_basins, time, BasalforcingsPicoFarOceansalinityEnum); + this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); + M = num_basins*maxbox; + this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + + Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); + + toc_farocean = t_farocean[basinid]; + soc_farocean = s_farocean[basinid]; + area_box1 = boxareas[basinid*maxbox+boxid]; + g1 = area_box1*gamma_T; + s1 = soc_farocean/(nu*lambda); + + IssmDouble basalmeltrates_shelf[NUMVERTICES]; + IssmDouble potential_pressure_melting_point[NUMVERTICES]; + IssmDouble Tocs[NUMVERTICES]; + IssmDouble Socs[NUMVERTICES]; + IssmDouble overturnings[NUMVERTICES]; + + /* Start looping on the number of nodes and calculate ocean vars */ + Gauss* gauss=this->NewGauss(); + for (int i=0;iGaussVertex(i); + thickness_input->GetInputValue(&thickness,gauss); + pressure = (rhoi*earth_grav*1e-4)*thickness; + T_star = a*soc_farocean+b-c*pressure-toc_farocean; + p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha)); + q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha))); + + /* To avoid negatives under the square root */ + if ((0.25*pow(p_coeff,2)-q_coeff)<0){ + q_coeff = (0.25*pow(p_coeff,2)); + } + + Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff)); + Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]); + potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; + basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); + overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i])); + } + + this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum); + + /*Cleanup and return */ + delete gauss; + xDelete(t_farocean); + xDelete(s_farocean); + xDelete(boxareas); +} +/*}}}*/ +void Tria::PicoUpdateNextBox(int loopboxid){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int boxid; this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); - this->inputs->GetInputValue(&maxbox,BasalforcingsPicoMaxboxcountEnum); - this->inputs->GetInputValue(&thickness,ThicknessEnum); - - _error_("to be continued"); - /* - - toc_farocean = t_farocean[basinid]; - soc_farocean = s_farocean[basinid]; - area_box1 = areas[basinid*maxbox+boxid]; - pressure = (rhoi*earth_grav)*thickness; - T_star = a*soc_farocean+b-c*pressure-toc_farocean; - g1 = area_box1*gamma_T; - s1 = soc_farocean/(nu*lambda); - p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha)); - q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha))); - if ((0.25*(p_coeff)^(2) - q_coeff) < 0) - {q_coeff = (0.25*(p_coeff)^(2))} - Toc = toc_farocean-(-0.5*p_coeff+sqrt(0.25*p_coeff^(2)-q_coeff)); - Soc = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Toc); - potential_pressure_melting_point = a*Soc+b-c*pressure; - basalmeltrate_shelf = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point-Toc); - overturning = overturning_coeff*rho_star*(Beta*(soc_farocean-Soc)-alpha*(toc_farocean-Toc)); - - this->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,this->GetElementType()); - */ - + if(loopboxid!=boxid) return; + + int basinid, maxbox, num_basins, numnodes,M; + IssmDouble rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c; + IssmDouble alpha, Beta, gamma_T, overturning_coeff; + IssmDouble thickness, toc_farocean, soc_farocean, area_boxi; + IssmDouble mean_toc,mean_soc,mean_overturning; + IssmDouble pressure, T_star, g1, g2, p_coeff, q_coeff; + IssmDouble* toc_weighted_avg=NULL; + IssmDouble* soc_weighted_avg=NULL; + IssmDouble* overturning_weighted_avg=NULL; + IssmDouble* boxareas=NULL; + + /*Get variables*/ + rhoi = this->GetMaterialParameter(MaterialsRhoIceEnum); + rhow = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); + earth_grav = this->GetMaterialParameter(ConstantsGEnum); + rho_star = 1033.; // kg/m^3 + nu = rhoi/rhow; + latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum); + c_p_ocean = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + lambda = latentheat/c_p_ocean; + a = -0.0572; // K/psu + b = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum); //K + c = 7.77e-4; + alpha = 7.5e-5; // 1/K + Beta = 7.7e-4; // K + + /* Get parameters and inputs */ + this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum); + this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); + this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); + this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); + this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum); + this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum); + this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum); + M=num_basins*maxbox; + this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); + this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + + Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); + + area_boxi = boxareas[basinid*maxbox+boxid]; + mean_toc = toc_weighted_avg[basinid]; + mean_soc = soc_weighted_avg[basinid]; + mean_overturning = overturning_weighted_avg[basinid]; + g1 = area_boxi*gamma_T; + g2 = g2/(nu*lambda); + + IssmDouble basalmeltrates_shelf[NUMVERTICES]; + IssmDouble potential_pressure_melting_point[NUMVERTICES]; + IssmDouble Tocs[NUMVERTICES]; + IssmDouble Socs[NUMVERTICES]; + + /* Start looping on the number of nodes and calculate ocean vars */ + Gauss* gauss=this->NewGauss(); + for (int i=0;iGaussVertex(i); + thickness_input->GetInputValue(&thickness,gauss); + pressure = (rhoi*earth_grav*1e-4)*thickness; + T_star = a*mean_soc+b-c*pressure-mean_toc; + Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc)); + Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); + potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; + basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); + } + + this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); + + /*Cleanup and return */ + delete gauss; + xDelete(boxareas); + xDelete(toc_weighted_avg); + xDelete(soc_weighted_avg); + xDelete(overturning_weighted_avg); } /*}}}*/ @@ -3720,4 +3848,10 @@ if(minvalue>fieldvalue || maxvalueGetInputListOnVertices(&lsf[0],fieldenum); + int nrice=0; + for(int i=0;i* potential_sheet_ungrounding); int PressureInterpolation(); @@ -165,4 +166,5 @@ void DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum); IssmDouble GetArea(void); + IssmDouble GetHorizontalSurfaceArea(void); IssmDouble GetArea3D(void); IssmDouble GetAreaIce(void); Index: /issm/trunk-jpl/src/c/classes/Params/BoolParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/BoolParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/BoolParam.h (revision 22625) @@ -45,4 +45,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/DataSetParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DataSetParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/DataSetParam.h (revision 22625) @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(FILE** pfile){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a file pointer");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} Index: /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.h (revision 22625) @@ -48,4 +48,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.h (revision 22625) @@ -48,4 +48,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/DoubleParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleParam.h (revision 22625) @@ -45,6 +45,7 @@ void GetParameterValue(int** pintarray,int* pM,int* pN); void GetParameterValue(IssmDouble* pIssmDouble){*pIssmDouble=value;} + void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} - void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM); Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 22625) @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/FileParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/FileParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/FileParam.h (revision 22625) @@ -45,4 +45,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/GenericParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/GenericParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/GenericParam.h (revision 22625) @@ -71,4 +71,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/IntMatParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/IntMatParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/IntMatParam.h (revision 22625) @@ -47,4 +47,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/IntParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/IntParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/IntParam.h (revision 22625) @@ -39,5 +39,5 @@ int ObjectEnum(); /*}}}*/ - /*Param vritual function definitions: {{{*/ + /*Param virtual function definitions: {{{*/ void GetParameterValue(bool* pbool){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a bool");} void GetParameterValue(int* pinteger){*pinteger=value;} @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/IntVecParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/IntVecParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/IntVecParam.h (revision 22625) @@ -47,4 +47,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/MatrixParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/MatrixParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/MatrixParam.h (revision 22625) @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/Param.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/Param.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/Param.h (revision 22625) @@ -34,4 +34,5 @@ virtual void GetParameterValue(IssmDouble* pIssmDouble)=0; virtual void GetParameterValue(IssmDouble* pdouble,IssmDouble time)=0; + virtual void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time)=0; virtual void GetParameterValue(char** pstring)=0; virtual void GetParameterValue(char*** pstringarray,int* pM)=0; Index: /issm/trunk-jpl/src/c/classes/Params/Parameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/Parameters.cpp (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/Parameters.cpp (revision 22625) @@ -300,4 +300,14 @@ } /*}}}*/ +void Parameters::FindParam(IssmDouble** pIssmDoublearray, int* pM, IssmDouble time, int param_enum){ _assert_(this);/*{{{*/ + + _assert_(param_enum>ParametersSTARTEnum); + _assert_(param_enumparams[index]) _error_("Parameter " << EnumToStringx(param_enum) <<" not set"); + this->params[index]->GetParameterValue(pIssmDoublearray, pM, time); +} +/*}}}*/ void Parameters::FindParam(char** pstring,int param_enum){ _assert_(this);/*{{{*/ @@ -349,6 +359,8 @@ int index = param_enum - ParametersSTARTEnum -1; + if(!this->params[index]) _error_("Parameter " << EnumToStringx(param_enum) <<" not set"); this->params[index]->GetParameterValue(pIssmDoublearray,pM); + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Params/Parameters.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/Parameters.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/Parameters.h (revision 22625) @@ -40,4 +40,5 @@ void FindParam(IssmDouble* pscalar, int enum_type); void FindParam(IssmDouble* pscalar, int enum_type,IssmDouble time); + void FindParam(IssmDouble** pIssmDoublearray, int* pM, IssmDouble time, int enum_type); void FindParam(char** pstring,int enum_type); void FindParam(char*** pstringarray,int* pM,int enum_type); Index: /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.h (revision 22625) @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} void GetParameterValue(char*** pstringarray,int* pM); Index: /issm/trunk-jpl/src/c/classes/Params/StringParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/StringParam.h (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/StringParam.h (revision 22625) @@ -46,4 +46,5 @@ void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");} void GetParameterValue(char** pstring); void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp (revision 22624) +++ /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp (revision 22625) @@ -142,2 +142,57 @@ } /*}}}*/ +void TransientArrayParam::GetParameterValue(IssmDouble** pIssmDoubleArray, int* pM, IssmDouble time){/*{{{*/ + + IssmDouble* output=NULL; + bool found; + int M=*pM; + + output=xNew(M); + + /*Ok, we have the time, go through the timesteps, and figure out which interval we + *fall within. Then interpolate the values on this interval: */ + if(timetimesteps[0]){ + /*get values for the first time: */ + for(int k=0; kvalues[k*this->N]; + found=true; + } + } + else if(time>this->timesteps[this->N-1]){ + /*get values for the last time: */ + for(int k=0; kvalues[(k+1)*this->N-1]; + found=true; + } + } + else{ + /*Find which interval we fall within: */ + for(int i=0;iN;i++){ + if(time==this->timesteps[i]){ + /*We are right on one step time: */ + for(int k=0; kvalues[k*this->N+i]; + } + found=true; + break; //we are done with the time interpolation. + } + else{ + if(this->timesteps[i]