Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26157) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26158) @@ -5534,14 +5534,83 @@ dI_list[2] = +4*M_PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; } - else if(masks->isiceonly[this->lid]){ - IssmDouble rho_ice, I; - - /*recover material parameters: */ + if(masks->isiceonly[this->lid]){ + + IssmDouble rho_ice,I; + bool computerigid= false; + + bool notfullygrounded=false; + + bool scaleoceanarea= false; + int glfraction=1; + IssmDouble phi=1.0; + /*{{{*/ + + + + if (masks->isfullyfloating[this->lid]){ + I=0; + this->AddInput(EffectivePressureEnum,&I,P0Enum); + return; + } + + + + if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. + + /*recover some parameters:*/ rho_ice=FindParam(MaterialsRhoIceEnum); - - /*Compute ice thickness change: */ + this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ + this->Element::GetInputValue(&area,AreaEnum); + + /*Compute fraction of the element that is grounded: */ + if(notfullygrounded){ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + if(glfraction==0)phi=1; + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); + #endif + } + else phi=1.0; + + /*Retrieve surface load for ice: */ Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); - if (!deltathickness_input)_error_("DeltaIceThicknessEnum delta ice thickness input needed to compute sea level change!"); - deltathickness_input->GetInputAverage(&I); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); + + /*/Average ice thickness over grounded area of the element only: {{{*/ + if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); + else{ + IssmDouble total_weight=0; + bool mainlyfloating = true; + int point1; + IssmDouble fraction1,fraction2; + + /*Recover portion of element that is grounded*/ + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + + /* Start looping on the number of gaussian points and average over these gaussian points: */ + total_weight=0; + I=0; + while(gauss->next()){ + IssmDouble Ig=0; + deltathickness_input->GetInputValue(&Ig,gauss); + I+=Ig*gauss->weight; + total_weight+=gauss->weight; + } + I=I/total_weight; + delete gauss; + } + /*}}}*/ + + /*}}}*/ + /*convert to kg/m^2*/ + I=I*phi; dI_list[0] = -4*M_PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; @@ -6053,5 +6122,7 @@ /*diverse:*/ int gsize; - IssmDouble I, S, BP; //change in relative ice thickness and sea level + IssmDouble I=0; + IssmDouble S=0; + IssmDouble BP=0; IssmDouble rho_ice,rho_water; int horiz; @@ -6065,4 +6136,10 @@ /*computational flags:*/ bool computeelastic= false; + bool computerigid= false; + bool notfullygrounded=false; + bool scaleoceanarea= false; + int glfraction=1; + IssmDouble area; + IssmDouble phi=1.0; /*early return if we are not on the ocean or on an ice cap:*/ @@ -6116,14 +6193,67 @@ } } - else if (masks->isiceonly[this->lid]){ - - /*Compute ice thickness change: */ + if (masks->isiceonly[this->lid]){ + + if (masks->isfullyfloating[this->lid]) return; + + + if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. + + /*recover some parameters:*/ + rho_ice=FindParam(MaterialsRhoIceEnum); + rho_water=FindParam(MaterialsRhoSeawaterEnum); + this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); + + /*Get area of element: precomputed in the sealevelchange_geometry:*/ + this->Element::GetInputValue(&area,AreaEnum); + + /*Compute fraction of the element that is grounded: */ + if(notfullygrounded){ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + if(glfraction==0)phi=1; + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticMaskEnum,&phi,P0Enum); + #endif + } + else phi=1.0; + + /*Retrieve surface load for ice: */ Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum); - if (!deltathickness_input)_error_("DeltaIceThicknessEnum delta thickness input needed to compute sea level change!"); - deltathickness_input->GetInputAverage(&I); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!"); + + /*/Average ice thickness over grounded area of the element only: {{{*/ + if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); + else{ + IssmDouble total_weight=0; + bool mainlyfloating = true; + int point1; + IssmDouble fraction1,fraction2; + + /*Recover portion of element that is grounded*/ + this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2); + + /* Start looping on the number of gaussian points and average over these gaussian points: */ + total_weight=0; + I=0; + while(gauss->next()){ + IssmDouble Ig=0; + deltathickness_input->GetInputValue(&Ig,gauss); + I+=Ig*gauss->weight; + total_weight+=gauss->weight; + } + I=I/total_weight; + delete gauss; + } + /*}}}*/ /*convert to kg/m^2*/ - I=I*rho_ice; - + I=I*rho_ice*phi; + for(int i=0;i