Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16815) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16816) @@ -192,5 +192,5 @@ virtual void UpdateBasalConstraintsEnthalpy(void)=0; virtual void ComputeBasalMeltingrate(void)=0; - virtual void DrainWaterfraction(void)=0; + virtual void DrainWaterfraction(IssmDouble* drainrate_element)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16815) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16816) @@ -4750,5 +4750,5 @@ /* Intermediaries */ bool isenthalpy, checkpositivethickness, istemperatelayer; - int i,j,analysis_type; + int i,j,iv,analysis_type; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble xyz_list_tria[NUMVERTICES2D][3]; @@ -4764,7 +4764,9 @@ IssmDouble geothermalflux[NUMVERTICES]; IssmDouble dt, yts; - IssmDouble melting_overshoot,meltingrate_enthalpy; - IssmDouble lambda,heating; + IssmDouble melting_overshoot,meltingrate_enthalpy[NUMVERTICES2D]; + IssmDouble drainrate_element[NUMVERTICES2D],drainrate_column[NUMVERTICES2D]; + IssmDouble lambda,heating[NUMVERTICES2D]; Friction *friction = NULL; + Penta *penta = NULL; /* Only compute melt rates at the base of grounded ice*/ @@ -4795,7 +4797,7 @@ friction=new Friction("3d",inputs,matpar,analysis_type); - /*Ok, get meltingrates now from basal conditions*/ + /******** MELTING RATES ************************************/ GaussPenta* gauss=new GaussPenta(); - for(int iv=0;iv<3;iv++){ + for(iv=0;ivGaussVertex(iv); @@ -4805,6 +4807,6 @@ /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/ - meltingrate_enthalpy=0.; - heating=0.; + meltingrate_enthalpy[iv]=0.; + heating[iv]=0.; if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure[iv]))){ /*ensure that no ice is at T=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true; + if(enthalpy[iv+NUMVERTICES2D]>=matpar->PureIceEnthalpy(pressure[iv+NUMVERTICES2D])) istemperatelayer=true; if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp else{ @@ -4839,29 +4841,53 @@ matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]); // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66 - heating=(heatflux+basalfriction+geothermalflux[iv]); - meltingrate_enthalpy=heating/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent - } - - /*Update water column, basal meltingrate*/ - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + heating[iv]=(heatflux+basalfriction+geothermalflux[iv]); + meltingrate_enthalpy[iv]=heating[iv]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent + } + } + // enthalpy might have been changed, update + this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); + + /******** DRAINAGE *****************************************/ + for(iv=0; ivDrainWaterfraction(&drainrate_element[0]); // TODO: make sure every vertex is only drained once + for(iv=0; ivIsOnSurface()) break; + penta=penta->GetUpperElement(); + } + // add drained water to melting rate + for(iv=0; ivparameters->FindParam(&dt,TimesteppingTimeStepEnum); + for(iv=0;iv(dt)){ - if(watercolumn[iv]+meltingrate_enthalpy*dt<0.){ - melting_overshoot=watercolumn[iv]+meltingrate_enthalpy*dt; - lambda=melting_overshoot/(meltingrate_enthalpy*dt); _assert_(lambda>0); _assert_(lambda<1); - basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy; + if(watercolumn[iv]+meltingrate_enthalpy[iv]*dt<0.){ + melting_overshoot=watercolumn[iv]+meltingrate_enthalpy[iv]*dt; + lambda=melting_overshoot/(meltingrate_enthalpy[iv]*dt); _assert_(lambda>0); _assert_(lambda<1); + basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy[iv]; watercolumn[iv]=0.; yts=365*24*60*60; - enthalpy[iv]+=dt/yts*lambda*heating; + enthalpy[iv]+=dt/yts*lambda*heating[iv]; } else{ - basalmeltingrate[iv]=meltingrate_enthalpy; - watercolumn[iv]+=dt*meltingrate_enthalpy; + basalmeltingrate[iv]=meltingrate_enthalpy[iv]; + watercolumn[iv]+=dt*meltingrate_enthalpy[iv]; } } else{ - basalmeltingrate[iv]=meltingrate_enthalpy; - watercolumn[iv]+=meltingrate_enthalpy; + basalmeltingrate[iv]=meltingrate_enthalpy[iv]; + watercolumn[iv]+=meltingrate_enthalpy[iv]; } - } + + _assert_(watercolumn[iv]>=0.); + } + /*feed updated variables back into model*/ this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); @@ -4875,14 +4901,15 @@ /*}}}*/ /*FUNCTION Penta::DrainWaterfraction{{{*/ -void Penta::DrainWaterfraction(void){ +void Penta::DrainWaterfraction(IssmDouble* drainrate_element){ /*Intermediaries*/ bool isenthalpy; + int iv, index0; IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES]; - IssmDouble watercolumnbase[NUMVERTICES]; + IssmDouble dw[NUMVERTICES]; IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; - IssmDouble latentheat, dt; - IssmDouble dw, dwc; - Penta *pentabase = NULL; + IssmDouble dt, height_element; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble rho_water, rho_ice; /*Check wether enthalpy is activated*/ @@ -4890,30 +4917,34 @@ if(!isenthalpy) return; - /*get basal element, needed for basal watercolumn*/ - pentabase=(Penta*)this->GetBasalElement(); - + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); this->GetInputListOnVertices(&pressure[0],PressureEnum); - pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - latentheat=matpar->GetLatentHeat(); - - for(int iv=0;ivEnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); - dw=DrainageFunctionWaterfraction(waterfraction[iv], dt); - /*drain water fraction & update enthalpy*/ - waterfraction[iv]-=dw; - matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); - /*add drained water to watercolumn at base*/ - dwc=dw*this->IceVolume(); - watercolumnbase[iv%NUMVERTICES2D]+=dwc; - } - - /*feed updated results back into model*/ + dw[iv]=DrainageFunctionWaterfraction(waterfraction[iv], dt); + } + + /*drain waterfraction, feed updated variables back into model*/ + for(iv=0;iv(dt)) + waterfraction[iv]-=dw[iv]*dt; + else + waterfraction[iv]-=dw[iv]; + matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); + } this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); - pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum)); - + + /*return meltwater column equivalent to drained water*/ + for(iv=0;iv division by zero. Abort"); + if(waterfraction<=w0) Dret=D0; - if((waterfraction>w0) && (waterfraction<=w1)) + else if((waterfraction>w0) && (waterfraction<=w1)) Dret=(D1-D0)/(w1-w0)*(waterfraction-w0)+D0; - if((waterfraction>w1) && (waterfraction<=w2)) + else if((waterfraction>w1) && (waterfraction<=w2)) Dret=(D2-D1)/(w2-w1)*(waterfraction-w1)+D1; else Dret=D2; - /*check if dt*Dret>waterfraction. If so, drain whole waterfraction*/ + /*drain only up to w0*/ if(dt==0.){ - if(Dret>waterfraction) - return waterfraction; + if((waterfraction>w0) && (waterfraction-Dret*ytsw0) && (waterfraction-dt*Dretwaterfraction) - return waterfraction; - else - return dt*Dret; - } } Index: /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 16815) +++ /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 16816) @@ -36,8 +36,8 @@ /*check feasibility*/ - _assert_(pressure>0); + _assert_(pressure>=0); _assert_(n>0); - _assert_(betaCC>0); - _assert_(referencetemperature>0); + _assert_(betaCC>=0); + _assert_(referencetemperature>=0); _assert_(heatcapacity>0); _assert_(latentheat>0); @@ -58,4 +58,6 @@ Tstar=Tpmp; waterfraction=(enthalpy - H_sp)/latentheat; + if (waterfraction > 0.01) + waterfraction = 0.01; } @@ -75,37 +77,2 @@ return B; } - -// /*Get stiffness from temperature, waterfraction and depth*/ -// IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){ -// /*Use parameterization for the rheology: Greve and Blatter 2009 -// * get enthalpy from temperature and water fraction, -// * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */ -// -// /*TODO: update params from model*/ -// IssmDouble rho_ice=910; // kg/m^3 -// IssmDouble g=9.81; //kg*m/s^2 -// IssmDouble heatcapacity=2009; // J/kg/K -// IssmDouble referencetemperature=253.15; -// IssmDouble betaCC=7.9*pow(10.,-8.); -// IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012 -// -// IssmDouble Tstar, enthalpy, pressure, B; -// _assert_(temperature>0); -// _assert_(waterfraction>0); -// _assert_(depth>0); -// -// /*get pressure*/ -// pressure= rho_ice*g*depth; -// Tstar=temperature-betaCC*pressure; // TODO: check whether plus or minus -// /*get enthalpy*/ -// if (Tstar < 273.15){ -// enthalpy=heatcapacity*(Tstar - referencetemperature); -// } -// else{ -// enthalpy=heatcapacity*(273.15 - referencetemperature) + waterfraction*latentheat; -// } -// -// B=LliboutryDuval(enthalpy, pressure, n, betaCC, referencetemperature, heatcapacity, latentheat); -// -// return B; -// }