Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16025) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16026) @@ -306,4 +306,6 @@ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.h\ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp\ + ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\ + ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\ ./modules/Delta18oParameterizationx/Delta18oParameterizationx.h\ ./modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp\ Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 16025) +++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 16026) @@ -92,5 +92,7 @@ if(isenthalpy) InputToResultx(femmodel,WaterfractionEnum); if(isenthalpy) InputToResultx(femmodel,EnthalpyEnum); - if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum); + if(isenthalpy) InputToResultx(femmodel,WatercolumnEnum); + //if(!isenthalpy) InputToResultx(femmodel,BasalforcingsMeltingRateEnum); + InputToResultx(femmodel,BasalforcingsMeltingRateEnum); femmodel->RequestedOutputsx(requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16025) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16026) @@ -111,4 +111,5 @@ else{ enthalpy_core(femmodel); + PostprocessingEnthalpyx(femmodel); } #else Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16025) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16026) @@ -127,6 +127,6 @@ #ifdef _HAVE_THERMAL_ virtual void UpdateThermalBasalConstraints(void)=0; - //virtual void ComputeBasalMeltingrate(void)=0; - //virtual void DrainWaterfraction(void)=0; + virtual void ComputeBasalMeltingrate(void)=0; + virtual void DrainWaterfraction(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16025) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16026) @@ -4700,4 +4700,5 @@ IssmDouble scalar,enthalpy,enthalpyup; IssmDouble pressure,pressureup; + IssmDouble watercolumn; IssmDouble basis[NUMVERTICES]; Friction* friction=NULL; @@ -4724,6 +4725,7 @@ Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - - /*Build frictoin element, needed later: */ + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + + /*Build friction element, needed later: */ friction=new Friction("3d",inputs,matpar,analysis_type); @@ -4741,27 +4743,31 @@ enthalpy_input->GetInputValue(&enthalpy,gauss); pressure_input->GetInputValue(&pressure,gauss); -// if(enthalpy>matpar->PureIceEnthalpy(pressure)){ -// enthalpy_input->GetInputValue(&enthalpyup,gaussup); -// pressure_input->GetInputValue(&pressureup,gaussup); -// if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){ -// //do nothing, don't add heatflux -// } -// else{ -// //need to change spcenthalpy according to Aschwanden -// //NEED TO UPDATE -// } -// } -// else{ - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); - if(reCast(dt)) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; -// } + if(enthalpy>=matpar->PureIceEnthalpy(pressure)){ + enthalpy_input->GetInputValue(&enthalpyup,gaussup); + pressure_input->GetInputValue(&pressureup,gaussup); + if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ + // temperate ice has positive thickness: grad enthalpy*n=0. + } + else{ + // only base temperate, set Dirichlet BCs in Penta::UpdateThermalBasalConstraints() + } + } + else{ + watercolumn_input->GetInputValue(&watercolumn,gauss); + if(watercolumn==0.){ + /*add geothermal heat flux and basal friction*/ + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + if(reCast(dt)) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*basis[i]; + } + else { /*do nothing (water layer acts as insulation)*/ } + } } @@ -5202,4 +5208,9 @@ this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); break; + case LliboutryDuvalEnum: + B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0, material->GetN()); + for(i=0;imaterial->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); + break; default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet"); @@ -5220,12 +5231,13 @@ /*Intermediary*/ bool isenthalpy,isdynamicbasalspc,istemperatelayer; - int numindices; - IssmDouble h_pmp,pressure; - int *indices = NULL; + int numindices, numindicesup; + IssmDouble h_pmp,pressure, pressureup; + IssmDouble enthalpy, enthalpyup; + int *indices = NULL, *indicesup = NULL; /* Only update Constraints at the base of grounded ice*/ - if(!IsOnBed() || !IsFloating()) return; - - /*Check wether dynamic basal boudary conditions are activated -> TODO: Johannes :) */ + if(!IsOnBed() || IsFloating()) return; + + /*Check wether dynamic basal boundary conditions are activated */ parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); if(!isenthalpy) return; @@ -5234,42 +5246,60 @@ if(!isdynamicbasalspc) return; - - /*Fetch indices of basal nodes for this finite element*/ + /*Fetch indices of basal & surface nodes for this finite element*/ BasalNodeIndices(&numindices,&indices,this->element_type); + SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); + _assert_(numindices==numindicesup); /*Get parameters and inputs: */ Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + + /*if there is a temperate layer of positive thickness, set enthalpy=h_pmp at that node*/ GaussPenta* gauss=new GaussPenta(); + GaussPenta* gaussup=new GaussPenta(); for(int i=0;iGaussNode(this->element_type,indices[i]); + gaussup->GaussNode(this->element_type,indicesup[i]); // TODO: check: are the nodes corresponding? /*Check wether there is a temperate layer at the base or not -> TODO: Johannes:) */ - istemperatelayer = false; - - /*Add Dirichlet constraint to this node if there is a positive thickness of temperate ice*/ - if(istemperatelayer){ - + /*check if node is temperate, if not, return*/ + enthalpy_input->GetInputValue(&enthalpy, gauss); + pressure_input->GetInputValue(&pressure, gauss); + if (enthalpyPureIceEnthalpy(pressure)){ + // TODO: reset, if necessary, all spcs to non-valid + continue; + } + /*check if upper node is temperate. if yes, then we have a temperate layer of positive thickness. if not, continue.*/ + enthalpy_input->GetInputValue(&enthalpyup, gaussup); + pressure_input->GetInputValue(&pressureup, gaussup); + istemperatelayer = false; + if (enthalpyup>=matpar->PureIceEnthalpy(pressureup)) + istemperatelayer=true; + + /*Add Dirichlet constraint to this node if there is no layer of temperate ice with positive thickness*/ + if(!istemperatelayer){ /*Calculate enthalpy at pressure melting point */ - pressure_input->GetInputValue(&pressure,gauss); h_pmp=matpar->PureIceEnthalpy(pressure); - /*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*/ - nodes[indices[i]]->ApplyConstraint(0,h_pmp); - } + nodes[indices[i]]->ApplyConstraint(1,h_pmp); + } + else { + /*remove spc*/ + nodes[indices[i]]->DofInFSet(0); + } } /*Free ressources:*/ xDelete(indices); -} -/*}}}*/ -/*FUNCTION Penta::ComputeBasalMeltrate{{{*/ -void Penta::ComputeBasalMeltrate(void){ + xDelete(indicesup); +} +/*}}}*/ +/*FUNCTION Penta::ComputeBasalMeltingrate{{{*/ +void Penta::ComputeBasalMeltingrate(void){ /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/ /* Intermediaries */ - bool isenthalpy; + bool isenthalpy, checkpositivethickness, istemperatelayer; int i,j,analysis_type; IssmDouble xyz_list[NUMVERTICES][3]; @@ -5279,11 +5309,12 @@ IssmDouble normal_base[3], d1enthalpy[3]; IssmDouble kappa; - IssmDouble meltrate[3], watercolumn[3]; - IssmDouble enthalpy, enthalpyup; + IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy; + IssmDouble enthalpy[NUMVERTICES], enthalpyup; IssmDouble pressure, pressureup; IssmDouble temperature, waterfraction; - IssmDouble latentheat; + IssmDouble latentheat, rho_ice; IssmDouble basalfriction, alpha2; IssmDouble vx,vy,vz; + IssmDouble connectivity; IssmDouble dt; Friction *friction = NULL; @@ -5298,5 +5329,7 @@ /*Fetch parameters and inputs */ latentheat=matpar->GetLatentHeat(); + rho_ice=this->matpar->GetRhoIce(); Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basalmeltingrate_input); Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); @@ -5312,5 +5345,5 @@ friction=new Friction("3d",inputs,matpar,analysis_type); - /*Ok, get meltrates now from basal conditions*/ + /*Ok, get meltingrates now from basal conditions*/ GaussPenta* gauss=new GaussPenta(); GaussPenta* gaussup=new GaussPenta(); @@ -5318,63 +5351,76 @@ gauss->GaussVertex(iv); gaussup->GaussVertex(iv+3); - - // TODO: make sure that no node is computed twice (insert mask) - + + checkpositivethickness=true; watercolumn_input->GetInputValue(&watercolumn[iv], gauss); - enthalpy_input->GetInputValue(&enthalpy, gauss); + basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss); + enthalpy_input->GetInputValue(&enthalpy[iv], gauss); pressure_input->GetInputValue(&pressure, gauss); - /*Calculate basal meltrate*/ - if((watercolumn[iv]>0.) && (enthalpyPureIceEnthalpy(pressure))){ - enthalpy=matpar->PureIceEnthalpy(pressure); - } - else if(enthalpyPureIceEnthalpy(pressure)){ - meltrate[iv]=0.; - watercolumn[iv]=0.; - continue; - } - - /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice) */ - enthalpy_input->GetInputValue(&enthalpyup, gaussup); - pressure_input->GetInputValue(&pressureup, gaussup); - if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ - for(i=0;i<3;i++) vec_heatflux[i]=0.; - } - else{ - enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); - for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; - } - - /*Get normal vector to the bed */ - BedNormal(&normal_base[0],xyz_list_tria); - - heatflux=0.; - for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - - matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure); - - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); - meltrate[iv]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; - - /*Update water column*/ + /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/ + meltingrate_enthalpy=0.; + if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure))){ + /*ensure that no ice is at TPureIceEnthalpy(pressure); + //meltingrate_enthalpy=0.; + } + else if(enthalpy[iv]PureIceEnthalpy(pressure)){ + /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/ + meltingrate_enthalpy=0.; + checkpositivethickness=false; + } + else {/*do nothing, go to next check*/} + + if(checkpositivethickness){ + /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */ + istemperatelayer=false; + enthalpy_input->GetInputValue(&enthalpyup, gaussup); + pressure_input->GetInputValue(&pressureup, gaussup); + if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) + istemperatelayer=true; + if(istemperatelayer) + for(i=0;i<3;i++) vec_heatflux[i]=0.; + else{ + enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure); + for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; + } + + /*geothermal heatflux*/ + BedNormal(&normal_base[0],xyz_list_tria); + heatflux=0.; + for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; + geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); + + /*basal friction*/ + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); + + matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure); + meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent + } + + /*Update water column, basal meltingrate*/ + connectivity=IssmDouble(vertices[iv]->Connectivity()); + meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex + basalmeltingrate[iv]+=meltingrate_enthalpy; this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); if(reCast(dt)) - watercolumn[iv]+=dt*meltrate[iv]; + watercolumn[iv]+=dt*meltingrate_enthalpy; else - watercolumn[iv]=meltrate[iv]; + watercolumn[iv]+=meltingrate_enthalpy; } - /*update meltrate and watercolumn*/ + /*feed updated variables back into model*/ + this->inputs->AddInput(new PentaInput(EnthalpyEnum, enthalpy, P1Enum)); this->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumn, P1Enum)); - this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum, meltrate, P1Enum)); + this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum, basalmeltingrate, P1Enum)); /*Clean up and return*/ delete gauss; delete gaussup; + delete friction; } /*}}}*/ @@ -5382,6 +5428,4 @@ void Penta::DrainWaterfraction(void){ -// TODO: create vector to mark whether node has been drained already i.o. to not drain nodes multiple times - /*Intermediaries*/ bool isenthalpy; @@ -5389,4 +5433,5 @@ IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; IssmDouble latentheat, dt; + IssmDouble connectivity; GaussPenta* gauss; @@ -5401,16 +5446,16 @@ this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); latentheat=matpar->GetLatentHeat(); + gauss=new GaussPenta(); for(int iv=0;ivGaussVertex(iv); - /*TODO: Check whether Vertex has been drained already*/ enthalpy_input->GetInputValue(&enthalpy[iv], gauss); pressure_input->GetInputValue(&pressure[iv], gauss); matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); - /*drain water fraction*/ - waterfraction[iv]-=dt*DrainageFunctionWaterfraction(waterfraction[iv]); - if(waterfraction[iv]<0.) waterfraction[iv]=0.; - /*update enthalpy*/ + /*drain water fraction & update enthalpy*/ + // TODO: replace connectivity-wise draining by draining once per node per timestep + connectivity=IssmDouble(vertices[iv]->Connectivity()); + waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity; matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); } @@ -5418,5 +5463,5 @@ this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); - this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum)); + // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum)); // temperature should not change during drainage... } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16025) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16026) @@ -357,5 +357,5 @@ void InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong); void UpdateThermalBasalConstraints(void); - void ComputeBasalMeltrate(void); + void ComputeBasalMeltingrate(void); void DrainWaterfraction(void); #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16025) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16026) @@ -255,5 +255,6 @@ #ifdef _HAVE_THERMAL_ void UpdateThermalBasalConstraints(void){_error_("not implemented yet");}; - void ComputeBasalMeltrate(void){_error_("not implemented yet");}; + void ComputeBasalMeltingrate(void){_error_("not implemented yet");}; + void DrainWaterfraction(void){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 16025) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 16026) @@ -44,4 +44,6 @@ iomodel->FetchDataToInput(elements,WaterfractionEnum); iomodel->FetchDataToInput(elements,BasalforcingsGeothermalfluxEnum); + iomodel->FetchDataToInput(elements,WatercolumnEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 16025) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 16026) @@ -75,4 +75,5 @@ #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" #include "./PositiveDegreeDayx/PositiveDegreeDayx.h" +#include "./PostprocessingEnthalpyx/PostprocessingEnthalpyx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" #include "./Reduceloadx/Reduceloadx.h" Index: /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 16025) +++ /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 16026) @@ -8,20 +8,36 @@ /*FUNCTION IssmDouble DrainageFunctionWaterfraction()*/ -IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction){ +IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.){ /* DrainageFunctionWaterfraction returns how much of the waterfraction is drained per year */ - _assert_(waterfraction>0); + _assert_(waterfraction>=0.); + _assert_(dt>=0.); IssmDouble w0=0.01, w1=0.02, w2=0.03; - IssmDouble D0=0, D1=0.005, D2=0.05; + IssmDouble Dret, D0=0, D1=0.005, D2=0.05; + /*get drainage function value*/ if((w0==w1)||(w1==w2)||(w0==w2)) - perror ("error in DrainageFunctionWaterfraction. Abort"); + _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort"); if(waterfraction<=w0) - return D0; + Dret=D0; if((waterfraction>w0) && (waterfraction<=w1)) - return (D1-D0)/(w1-w0)*(waterfraction-w0)+D0; + Dret=(D1-D0)/(w1-w0)*(waterfraction-w0)+D0; if((waterfraction>w1) && (waterfraction<=w2)) - return (D2-D1)/(w2-w1)*(waterfraction-w1)+D1; + Dret=(D2-D1)/(w2-w1)*(waterfraction-w1)+D1; else - return D2; + Dret=D2; + + /*check if dt*Dret>waterfraction. If so, drain whole waterfraction*/ + if(dt==0.){ + if(Dret>waterfraction) + return waterfraction; + else + return Dret; + } + else{ + if(dt*Dret>waterfraction) + return waterfraction; + else + return dt*Dret; + } } Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 16025) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 16026) @@ -20,5 +20,5 @@ IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); -IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction); +IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.); /*Print arrays*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16025) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16026) @@ -236,4 +236,6 @@ SurfaceforcingsBPosEnum, SurfaceforcingsBNegEnum, + ThermalIsenthalpyEnum, + ThermalIsdynamicbasalspcEnum, ThermalMaxiterEnum, ThermalPenaltyFactorEnum, @@ -242,5 +244,4 @@ ThermalSpctemperatureEnum, ThermalStabilizationEnum, - ThermalIsenthalpyEnum, GiaMantleViscosityEnum, GiaLithosphereThicknessEnum, Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 16025) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 16026) @@ -671,4 +671,5 @@ if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end; if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end; + if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);end; %bedinfo and surface info Index: /issm/trunk-jpl/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.m (revision 16025) +++ /issm/trunk-jpl/src/m/classes/thermal.m (revision 16026) @@ -13,4 +13,5 @@ penalty_factor = 0; isenthalpy = 0; + isdynamicbasalspc = 0; end methods @@ -39,4 +40,7 @@ %Should we use cold ice (default) or enthalpy formulation obj.isenthalpy=0; + + %will basal boundary conditions be set dynamically + obj.isdynamicbasalspc=0; end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ @@ -52,4 +56,5 @@ md = checkfield(md,'thermal.spctemperature(find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN))','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos),'message','spctemperature should be below the adjusted melting point'); md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0 1]); + md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0 1]); end end % }}} @@ -63,5 +68,6 @@ fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)'); fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); - + fielddisplay(obj,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']); + end % }}} function marshall(obj,md,fid) % {{{ @@ -73,5 +79,6 @@ WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean'); - end % }}} + WriteData(fid,'object',obj,'fieldname','isdynamicbasalspc','format','Boolean'); + end % }}} end end