Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27286) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27287) @@ -132,4 +132,6 @@ iomodel->ConstantToInput(inputs,elements,0.,CalvingratexEnum,P1Enum); iomodel->ConstantToInput(inputs,elements,0.,CalvingrateyEnum,P1Enum); + break; + case CalvingPollardEnum: break; @@ -228,4 +230,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_threshold",CalvingVelThresholdEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.calving.vel_upperbound",CalvingVelUpperboundEnum)); + break; + case CalvingPollardEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.rc",CalvingRcEnum)); break; default: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27286) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27287) @@ -297,5 +297,5 @@ this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,&values[0],P1Enum); }/*}}}*/ -void Element::CalvingFromRate(){/*{{{*/ +void Element::CalvingRateToVector(){/*{{{*/ /*We are provided a calving rate, figure out the x/y components*/ @@ -3694,4 +3694,7 @@ case CalvingLevermannEnum: this->CalvingRateLevermann(); + break; + case CalvingPollardEnum: + this->CalvingPollard(); break; case CalvingVonmisesEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27286) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27287) @@ -72,5 +72,5 @@ void BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation); void CalvingSetZeroRate(void); - void CalvingFromRate(void); + void CalvingRateToVector(void); void ComputeLambdaS(void); void ComputeNewDamage(); @@ -238,4 +238,5 @@ virtual void CalvingCrevasseDepth(void){_error_("not implemented yet");}; virtual void CalvingRateLevermann(void)=0; + virtual void CalvingPollard(void){_error_("not implemented yet");}; virtual void CalvingFluxLevelset(void){_error_("not implemented yet");}; virtual void CalvingMeltingFluxLevelset(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 27286) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 27287) @@ -324,5 +324,5 @@ /*Add input*/ this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingFromRate(); + this->CalvingRateToVector(); /*Extrude*/ @@ -365,5 +365,5 @@ /*Add input*/ this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingFromRate(); + this->CalvingRateToVector(); }/*}}}*/ void Penta::CalvingFluxLevelset(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27286) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27287) @@ -364,5 +364,5 @@ /*Add input*/ this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingFromRate(); + this->CalvingRateToVector(); } /*}}}*/ @@ -558,5 +558,93 @@ /*Add input*/ this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingFromRate(); + this->CalvingRateToVector(); +}/*}}}*/ +void Tria::CalvingPollard(){/*{{{*/ + + /*Intermediaries*/ + IssmDouble calvingrate[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dvx[2], dvy[2]; + IssmDouble B, n, H, bed, vx, vy, vel, smb; + IssmDouble ds, db, da, dt, dw, r, R; + + /*Retrieve all inputs and parameters we will need*/ + IssmDouble rc = FindParam(CalvingRcEnum); + IssmDouble rho_ice = FindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = FindParam(MaterialsRhoSeawaterEnum); + IssmDouble gravity = FindParam(ConstantsGEnum); + IssmDouble mig_max = FindParam(MigrationMaxEnum); + + /*Retrieve all inputs and parameters we will need */ + Input *bs_input = this->GetInput(BaseEnum); _assert_(bs_input); + Input *vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input *B_input = this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); + Input *n_input = this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input *H_input = this->GetInput(ThicknessEnum); _assert_(H_input); + Input *smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input); + + /* Start looping on the number of vertices: */ + GaussTria gauss; + for (int iv=0;ivGetInputValue(&bed,&gauss); + + /*Only calve if bed is below sea level, as always*/ + if(bed<=0.){ + + /*Get Triangle node coordinates*/ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Get strain rates*/ + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss); + + /*Get other inputs*/ + B_input->GetInputValue(&B,&gauss); + n_input->GetInputValue(&n,&gauss); + H_input->GetInputValue(&H,&gauss); + vx_input->GetInputValue(&vx,&gauss); + vy_input->GetInputValue(&vy,&gauss); + smb_input->GetInputValue(&smb,&gauss); + + /*1. with surface crevasses, ds*/ + ds = 2./(rho_ice*gravity) * B * pow( max(0.,dvx[0]) + max(0.,dvy[1]) , 1./n); + + /*2. basal crevasses*/ + db = (rho_ice)/(rho_water - rho_ice) * ds; + + /*3. "Additional" crevasse opening*/ + vel = sqrt(vx*vx + vy*vy)/(365.25*24*3600); + da = H* max(0., log(vel/1600.))/log(1.2); + + /*4. deal with shallow ice*/ + dt = H* max(0., min(1., (150. - H)/50.)); + + /*5. water induced opening*/ + dw = 0.; + R = smb*365.25*24*3600; //convert from m/s to m/yr + if(R>1.5 && R<=3.){ + dw = 4*1.5*(R - 1.5); + } + else if(R>3.){ + dw = R*R; + } + + /*Total calving rate*/ + r = (ds+db+da+dt+dw)/H; + calvingrate[iv]= mig_max * max(0., min(1., (r - rc)/(1 - rc))); //P&DC: mig_max = 3000 m/yr + _assert_(!xIsNan(calvingrate[iv])); + _assert_(!xIsInf(calvingrate[iv])); + } + else + calvingrate[iv]=0.; + } + + /*Add input*/ + this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); + this->CalvingRateToVector(); }/*}}}*/ void Tria::CalvingFluxLevelset(){/*{{{*/ @@ -951,5 +1039,5 @@ this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); - this->CalvingFromRate(); + this->CalvingRateToVector(); /*Clean up and return*/ @@ -4312,4 +4400,5 @@ case CalvingVonmisesEnum: case CalvingLevermannEnum: + case CalvingPollardEnum: case CalvingTestEnum: case CalvingParameterizationEnum: @@ -4357,4 +4446,5 @@ case CalvingParameterizationEnum: case CalvingLevermannEnum: + case CalvingPollardEnum: calvingratex_input->GetInputValue(&c[0],&gauss); calvingratey_input->GetInputValue(&c[1],&gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27286) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27287) @@ -58,4 +58,5 @@ void CalvingCrevasseDepth(); void CalvingRateLevermann(); + void CalvingPollard(); void CalvingFluxLevelset(); void CalvingMeltingFluxLevelset(); Index: /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 27286) +++ /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 27287) @@ -20,5 +20,5 @@ break; case DefaultCalvingEnum: - femmodel->ElementOperationx(&Element::CalvingFromRate); + femmodel->ElementOperationx(&Element::CalvingRateToVector); break; case CalvingCrevasseDepthEnum: @@ -31,5 +31,4 @@ break; case CalvingLevermannEnum: - if(VerboseModule()) _printf0_(" computing Levermann's calving rate\n"); femmodel->StrainRateparallelx(); femmodel->StrainRateperpendicularx(); @@ -46,4 +45,7 @@ femmodel->ElementOperationx(&Element::CalvingRateParameterization); break; + case CalvingPollardEnum: + femmodel->ElementOperationx(&Element::CalvingPollard); + break; default: _error_("Caving law "<