source: issm/oecreview/Archive/21724-22754/ISSM-22469-22470.diff

Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 22469)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 22470)
@@ -811,7 +811,7 @@
 
                /*Update signed distance*/
                dmin = sqrt(dmin);
-               if(dmin>10000) dmin=10000;
+               // if(dmin>10000) dmin=10000;
                if(ls[j]>0){
                        ls[j] = dmin;
                }
@@ -2736,6 +2736,99 @@
        return TriaRef::NumberofNodes(this->VelocityInterpolation());
 }
 /*}}}*/
+void       Tria::PicoUpdateBoxid(int* pmax_boxid_basin){/*{{{*/
+
+       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;
+
+       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);
+
+       Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
+       dist_gl_input->GetInputValue(&dist_gl,gauss);
+       dist_cf_input->GetInputValue(&dist_cf,gauss);
+
+       rel_dist_gl=dist_gl/(dist_gl+dist_cf);
+
+       boxid=-1;
+       for(i=0;i<boxid_max;i++){
+               if(rel_dist_gl>=(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());}
+
+       this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));        
+
+       /*Cleanup & return: */
+       delete gauss;
+}
+/*}}}*/
+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
+       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(MaterialsThermalExchangeVelocityEnum);
+       lambda          = latentheat/c_p_ocean;
+       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
+
+       this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
+       this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
+       this->parameters->FindParam(&t_farocean,BasalforcingsPicoFarOceanTemperatureEnum);
+       this->parameters->FindParam(&s_farocean,BasalforcingsPicoFarOceanSalinityEnum);
+
+       this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+       this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+       this->inputs->GetInputValue(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+       this->inputs->GetInputValue(&thickness,ThicknessEnum);
+
+       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());
+
+}
+/*}}}*/
 void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
 
        IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
@@ -2967,7 +3060,7 @@
        /*Get out if this is not an element input*/
        if(!IsInput(control_enum)) return;
 
-       /*Prepare index list*/
+       /*hrepare index list*/
        GradientIndexing(&vertexpidlist[0],control_index);
 
        /*Get values on vertices*/