Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22732) @@ -278,6 +278,5 @@ virtual int NumberofNodesVelocity(void)=0; virtual void PicoUpdateBoxid(int* pmax_boxid_basin)=0; - virtual void PicoUpdateFirstBox()=0; - virtual void PicoUpdateNextBox(int loopboxid)=0; + virtual void PicoUpdateBox(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 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 22732) @@ -140,6 +140,5 @@ 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 PicoUpdateBox(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 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 22732) @@ -129,6 +129,5 @@ 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 PicoUpdateBox(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 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 22732) @@ -137,6 +137,5 @@ 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 PicoUpdateBox(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 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22732) @@ -2817,102 +2817,16 @@ this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); }/*}}}*/ -void Tria::PicoUpdateFirstBox(){/*{{{*/ +void Tria::PicoUpdateBox(int loopboxid){/*{{{*/ 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, gamma_T, overturning_coeff,thickness; - IssmDouble pressure, T_star,p_coeff, q_coeff,toc_farocean,soc_farocean; - IssmDouble* boxareas = NULL; - - /*Get variables*/ - 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 */ - 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(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); - this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum); - this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); - this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); - this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - 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]; - 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*p_coeff*p_coeff; - } - - 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(boxareas); -} -/*}}}*/ -void Tria::PicoUpdateNextBox(int loopboxid){/*{{{*/ - - if(!this->IsIceInElement() || !this->IsFloating()) return; - int boxid; this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); if(loopboxid!=boxid) return; - int basinid, maxbox, num_basins, numnodes,M; - IssmDouble gamma_T, overturning_coeff; - IssmDouble thickness, toc_farocean, soc_farocean; + int basinid, maxbox, num_basins, numnodes, M; + IssmDouble gamma_T, overturning_coeff, thickness; 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; + IssmDouble* boxareas = NULL; /*Get variables*/ @@ -2931,52 +2845,103 @@ IssmDouble Beta = 7.7e-4; // K - /* Get parameters and inputs */ + /* Get non-box-specific 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); this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - - _assert_(basinid<=num_basins); + Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); + _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]; - 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*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]); - } - - this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); - - /*Cleanup and return */ - delete gauss; + IssmDouble potential_pressure_melting_point[NUMVERTICES]; + IssmDouble Tocs[NUMVERTICES]; + IssmDouble Socs[NUMVERTICES]; + + /* First box calculations */ + if(boxid==0){ + /* Get box1 parameters and inputs */ + IssmDouble time, toc_farocean, soc_farocean; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); + this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum); + IssmDouble s1 = soc_farocean/(nu*lambda); + IssmDouble overturnings[NUMVERTICES]; + + /* Start looping on the number of verticies 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*p_coeff*p_coeff; + } + + 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; + } + + /* Subsequent box calculations */ + else { + /* Get subsequent box parameters and inputs */ + IssmDouble* toc_weighted_avg = NULL; + IssmDouble* soc_weighted_avg = NULL; + IssmDouble* overturning_weighted_avg = NULL; + 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); + IssmDouble mean_toc = toc_weighted_avg[basinid]; + IssmDouble mean_soc = soc_weighted_avg[basinid]; + IssmDouble mean_overturning = overturning_weighted_avg[basinid]; + IssmDouble g2 = g1/(nu*lambda); + + /* Start looping on the number of verticies and calculate ocean vars */ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + thickness_input->GetInputValue(&thickness,gauss); + 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]); + } + + this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); + + /*Cleanup and return*/ + xDelete(toc_weighted_avg); + xDelete(soc_weighted_avg); + xDelete(overturning_weighted_avg); + delete gauss; + } + + /*Cleanup and return*/ xDelete(boxareas); - xDelete(toc_weighted_avg); - xDelete(soc_weighted_avg); - xDelete(overturning_weighted_avg); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22731) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22732) @@ -111,6 +111,5 @@ int NumberofNodesVelocity(void); void PicoUpdateBoxid(int* pmax_boxid_basin); - void PicoUpdateFirstBox(); - void PicoUpdateNextBox(int loopboxid); + void PicoUpdateBox(int loopboxid); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp (revision 22731) +++ /issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp (revision 22732) @@ -24,14 +24,14 @@ _assert_(in_values && in_time); - enum_type=in_enum_type; - M=in_M; //Number of time steps - N=in_N; //Number of rows - interpolation=interpolation_on; + this->enum_type=in_enum_type; + this->M=in_M; //Number of rows + this->N=in_N; //Number of timesteps + this->interpolation=interpolation_on; - values=xNew(M*N); + this->values=xNew(M*N); xMemCpy(values,in_values,M*N); - timesteps=xNew(M); - xMemCpy(timesteps,in_time,M); + this->timesteps=xNew(N); + xMemCpy(timesteps,in_time,N); } /*}}}*/ @@ -100,6 +100,7 @@ IssmDouble output; bool found; + _assert_(row>=0 && rowM); - /*Ok, we have the time, go through the timesteps, and figure out which interval we + /*Ok, we have the time and row, go through the timesteps, and figure out which interval we *fall within. Then interpolate the values on this interval: */ if(timetimesteps[0]){ @@ -118,5 +119,5 @@ if(time==this->timesteps[i]){ /*We are right on one step time: */ - output=this->values[row*this->N+i]; + output = this->values[row*this->N+i]; found=true; break; //we are done with the time interpolation. @@ -124,7 +125,7 @@ else{ if(this->timesteps[i]