Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22627) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22628) @@ -2777,16 +2777,16 @@ } /*}}}*/ -void Tria::PicoUpdateBoxid(int* pmax_boxid_basin){/*{{{*/ +void Tria::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/ if(!this->IsIceInElement() || !this->IsFloating()) return; - int i,boxid,basin_id,boxid_max; - IssmDouble dist_gl,dist_cf,rel_dist_gl,lowbound,highbound; + int basin_id; + IssmDouble dist_gl,dist_cf; inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum); - boxid_max=pmax_boxid_basin[basin_id]; + IssmDouble boxid_max=reCast(max_boxid_basin_list[basin_id]); Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input); - Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); + Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); /*Get dist_gl and dist_cf at center of element*/ @@ -2794,17 +2794,18 @@ dist_gl_input->GetInputValue(&dist_gl,gauss); dist_cf_input->GetInputValue(&dist_cf,gauss); + delete 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);} + dist_gl = fabs(dist_gl); + dist_cf = fabs(dist_cf); /*Compute relative distance to grounding line*/ - rel_dist_gl=dist_gl/(dist_gl+dist_cf); + IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf); /*Assign box numbers based on rel_dist_gl*/ - boxid=-1; - for(i=0;i=lowbound && rel_dist_gl<=highbound){ boxid=i; @@ -2815,9 +2816,5 @@ this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); - - /*Cleanup & return: */ - delete gauss; -} -/*}}}*/ +}/*}}}*/ void Tria::PicoUpdateFirstBox(){/*{{{*/ @@ -2829,26 +2826,22 @@ 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; + IssmDouble time, gamma_T, overturning_coeff,thickness; + IssmDouble pressure, T_star,p_coeff, q_coeff,toc_farocean,soc_farocean; + 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 + IssmDouble rhoi = this->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rhow = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble earth_grav = this->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_star = 1033.; // kg/m^3 + IssmDouble nu = rhoi/rhow; + IssmDouble latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble c_p_ocean = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + IssmDouble lambda = latentheat/c_p_ocean; + IssmDouble a = -0.0572; // K/psu + IssmDouble b = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum); //K + IssmDouble c = 7.77e-4; + IssmDouble alpha = 7.5e-5; // 1/K + IssmDouble Beta = 7.7e-4; // K /* Get parameters and inputs */ @@ -2857,19 +2850,15 @@ this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); - this->parameters->FindParam(&t_farocean, &num_basins, time, BasalforcingsPicoFarOceantemperatureEnum); - this->parameters->FindParam(&s_farocean, &num_basins, time, BasalforcingsPicoFarOceansalinityEnum); + this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); + this->parameters->FindParam(&soc_farocean, basinid, 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 area_box1 = boxareas[basinid*maxbox+boxid]; + IssmDouble g1 = area_box1*gamma_T; + IssmDouble s1 = soc_farocean/(nu*lambda); + + /*Define new inputs*/ IssmDouble basalmeltrates_shelf[NUMVERTICES]; IssmDouble potential_pressure_melting_point[NUMVERTICES]; @@ -2880,15 +2869,15 @@ /* 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))); + 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)); + if((0.25*pow(p_coeff,2)-q_coeff)<0){ + q_coeff = 0.25*p_coeff*p_coeff; } @@ -2907,6 +2896,4 @@ /*Cleanup and return */ delete gauss; - xDelete(t_farocean); - xDelete(s_farocean); xDelete(boxareas); } @@ -2920,29 +2907,27 @@ 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; + int basinid, maxbox, num_basins, numnodes,M; + IssmDouble gamma_T, overturning_coeff; + IssmDouble thickness, toc_farocean, soc_farocean; + IssmDouble pressure, T_star,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 + IssmDouble rhoi = this->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rhow = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble earth_grav = this->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_star = 1033.; // kg/m^3 + IssmDouble nu = rhoi/rhow; + IssmDouble latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble c_p_ocean = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + IssmDouble lambda = latentheat/c_p_ocean; + IssmDouble a = -0.0572; // K/psu + IssmDouble b = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum); //K + IssmDouble c = 7.77e-4; + IssmDouble alpha = 7.5e-5; // 1/K + IssmDouble Beta = 7.7e-4; // K /* Get parameters and inputs */ @@ -2954,16 +2939,15 @@ 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); + _assert_(basinid<=num_basins); + IssmDouble area_boxi = boxareas[basinid*maxbox+boxid]; + IssmDouble mean_toc = toc_weighted_avg[basinid]; + IssmDouble mean_soc = soc_weighted_avg[basinid]; + IssmDouble mean_overturning = overturning_weighted_avg[basinid]; + IssmDouble g1 = area_boxi*gamma_T; + IssmDouble g2 = g2/(nu*lambda); IssmDouble basalmeltrates_shelf[NUMVERTICES]; @@ -2974,11 +2958,11 @@ /* 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)); + pressure = (rhoi*earth_grav*1.e-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]); Index: /issm/trunk-jpl/src/c/classes/Params/BoolParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/BoolParam.h (revision 22627) +++ /issm/trunk-jpl/src/c/classes/Params/BoolParam.h (revision 22628) @@ -45,5 +45,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DataSetParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.h (revision 22628) @@ -48,5 +48,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.h (revision 22628) @@ -48,5 +48,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleParam.h (revision 22628) @@ -46,5 +46,4 @@ 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(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 22628) @@ -77,31 +77,23 @@ /*DoubleVecParam virtual functions definitions: */ -void DoubleVecParam::GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/ - IssmDouble* output=NULL; - int M; +void DoubleVecParam::GetParameterValue(IssmDouble** poutput,int* pM){/*{{{*/ - M=this->M; - output=xNew(M); + IssmDouble* output=xNew(M); xMemCpy(output,values,M); /*Assign output pointers:*/ if(pM) *pM=M; - *pIssmDoublearray=output; + *poutput=output; } /*}}}*/ -void DoubleVecParam::GetParameterValue(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/ - IssmDouble* output=NULL; - int M; - int N; +void DoubleVecParam::GetParameterValue(IssmDouble** poutput,int* pM,int* pN){/*{{{*/ - N=1; - M=this->M; - output=xNew(M); + IssmDouble* output=xNew(this->M); xMemCpy(output,values,M); /*Assign output pointers:*/ - if(pM) *pM=M; - if(pN) *pN=N; - *pIssmDoublearray=output; + if(pM) *pM=this->M; + if(pN) *pN=1; + *poutput=output; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 22627) +++ /issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/FileParam.h (revision 22628) @@ -45,5 +45,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/GenericParam.h (revision 22628) @@ -71,5 +71,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/IntMatParam.h (revision 22628) @@ -47,5 +47,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/IntParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/IntVecParam.h (revision 22628) @@ -47,5 +47,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/MatrixParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/Param.h (revision 22628) @@ -34,5 +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(IssmDouble* pdouble,int row, IssmDouble time){_error_("not implemented yet");}; 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/Parameters.cpp (revision 22628) @@ -300,12 +300,12 @@ } /*}}}*/ -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(IssmDouble* pscalar,int row,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(pscalar,row,time); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Params/Parameters.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/Parameters.h (revision 22627) +++ /issm/trunk-jpl/src/c/classes/Params/Parameters.h (revision 22628) @@ -40,5 +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(IssmDouble* pscalar, int row,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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/StringArrayParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/StringParam.h (revision 22628) @@ -46,5 +46,4 @@ 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 22627) +++ /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp (revision 22628) @@ -25,13 +25,13 @@ enum_type=in_enum_type; - N=in_N; //Number of time steps - M=in_M; //Number of rows + M=in_M; //Number of time steps + N=in_N; //Number of rows interpolation=interpolation_on; - values=xNew(N*M); - xMemCpy(values,in_values,N*M); + values=xNew(M*N); + xMemCpy(values,in_values,M*N); - timesteps=xNew(N); - xMemCpy(timesteps,in_time,N); + timesteps=xNew(M); + xMemCpy(timesteps,in_time,M); } /*}}}*/ @@ -39,11 +39,10 @@ xDelete(values); xDelete(timesteps); -} -/*}}}*/ +}/*}}}*/ /*Object virtual functions definitions:*/ Param* TransientArrayParam::copy() {/*{{{*/ - return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->N,this->M); + return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->M,this->N); } @@ -68,5 +67,5 @@ _printf_("TransientArrayParam:\n"); _printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n"); - _printf_(" size: " << this->N << " by " << this->M << "\n"); + _printf_(" size: " << this->M << " by " << this->N << "\n"); } @@ -83,9 +82,9 @@ MARSHALLING(N); if(marshall_direction==MARSHALLING_BACKWARD){ - values=xNew(N*M); - timesteps=xNew(N*M); + values = xNew(M*N); + timesteps = xNew(N); } - MARSHALLING_ARRAY(values,IssmDouble,N*M); - MARSHALLING_ARRAY(timesteps,IssmDouble,N*M); + MARSHALLING_ARRAY(values,IssmDouble,M*N); + MARSHALLING_ARRAY(timesteps,IssmDouble,N); }/*}}}*/ @@ -94,12 +93,11 @@ return TransientArrayParamEnum; -} -/*}}}*/ +}/*}}}*/ /*TransientArrayParam virtual functions definitions: */ -void TransientArrayParam::GetParameterValue(IssmDouble* pdouble,IssmDouble time, int row){/*{{{*/ +void TransientArrayParam::GetParameterValue(IssmDouble* pdouble,int row,IssmDouble time){/*{{{*/ IssmDouble output; - bool found; + bool found; /*Ok, we have the time, go through the timesteps, and figure out which interval we @@ -142,57 +140,2 @@ } /*}}}*/ -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]