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

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -811 +811 @@
 
                 /*Update signed distance*/
                 dmin = sqrt(dmin);
-                if(dmin>10000) dmin=10000;
+                // if(dmin>10000) dmin=10000;
                 if(ls[j]>0){
                         ls[j] = dmin;
                 }
@@ -2736 +2736 @@
         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 +3060 @@
         /*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*/