Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15869) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15870) @@ -5058,110 +5058,107 @@ /*FUNCTION Penta::ComputeBasalMeltrate{{{*/ void Penta::ComputeBasalMeltrate(void){ - /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/ - - /*Parameters*/ - int dim; - parameters->FindParam(&dim, MeshDimensionEnum); - - /* Intermediaries */ - bool isenthalpy; - int i, analysis_type; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble heatflux, geothermalflux_value; - IssmDouble vec_heatflux[dim]; - IssmDouble normal_base[dim], d1enthalpy[dim]; - IssmDouble kappa; - IssmDouble enthalpy, enthalpyup; - IssmDouble pressure, pressureup; - IssmDouble meltrate, watercolumn; - IssmDouble temperature, waterfraction; - IssmDouble latentheat; - IssmDouble basalfriction, alpha2; - IssmDouble vx,vy,vz; - IssmDouble dt; - Friction* friction=NULL; - GaussPenta* gauss=NULL; - GaussPenta* gaussup=NULL; - - /*check that dimension is 3*/ - if(dim!=3) return; - /* Only compute melt rates at the base of grounded ice*/ - if(!IsOnBed() || IsFloating()) return; - /*Check wether enthalpy is activated*/ - parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); - if(!isenthalpy) return; - - /*Fetch parameters: */ - latentheat=matpar->GetLatentHeat(); - - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - - /*Build friction element, needed later: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - friction=new Friction("3d",inputs,matpar,analysis_type); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(0,1,2,2); - gaussup=new GaussPenta(3,4,5,2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - gaussup->GaussPoint(ig); - - watercolumn_input->GetInputValue(&watercolumn, gauss); - enthalpy_input->GetInputValue(&enthalpy, gauss); - pressure_input->GetInputValue(&pressure, gauss); - - /*Calculate basal meltrate*/ - if((watercolumn>0) && (enthalpyPureIceEnthalpy(pressure))) - enthalpy=matpar->PureIceEnthalpy(pressure); - else if(enthalpyPureIceEnthalpy(pressure)){ - meltrate=0.; //TODO: set zero meltrate and watercolumn in model - watercolumn=0.; - return; - } - - /*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;iGetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); - for(i=0;iGetInputValue(&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=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; - - /*Update water column*/ - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - if(reCast(dt)) - watercolumn+=dt*meltrate; - else - watercolumn=meltrate; - // TODO: feed meltrate & watercolumn back to model - } + /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/ + + /* Intermediaries */ + bool isenthalpy; + int i,j,analysis_type; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble xyz_list_tria[NUMVERTICES2D][3]; + IssmDouble heatflux, geothermalflux_value; + IssmDouble vec_heatflux[3]; + IssmDouble normal_base[3], d1enthalpy[3]; + IssmDouble kappa; + IssmDouble enthalpy, enthalpyup; + IssmDouble pressure, pressureup; + IssmDouble meltrate, watercolumn; + IssmDouble temperature, waterfraction; + IssmDouble latentheat; + IssmDouble basalfriction, alpha2; + IssmDouble vx,vy,vz; + IssmDouble dt; + Friction *friction = NULL; + GaussPenta *gauss = NULL; + GaussPenta *gaussup = NULL; + + /* Only compute melt rates at the base of grounded ice*/ + if(!IsOnBed() || IsFloating()) return; + + /*Check wether enthalpy is activated*/ + parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + if(!isenthalpy) return; + + /*Fetch parameters and inputs */ + latentheat=matpar->GetLatentHeat(); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); + + for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); + friction=new Friction("3d",inputs,matpar,analysis_type); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(0,1,2,2); + gaussup=new GaussPenta(3,4,5,2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + gaussup->GaussPoint(ig); + + watercolumn_input->GetInputValue(&watercolumn, gauss); + enthalpy_input->GetInputValue(&enthalpy, gauss); + pressure_input->GetInputValue(&pressure, gauss); + + /*Calculate basal meltrate*/ + if((watercolumn>0.) && (enthalpyPureIceEnthalpy(pressure))){ + enthalpy=matpar->PureIceEnthalpy(pressure); + } + else if(enthalpyPureIceEnthalpy(pressure)){ + meltrate=0.; //TODO: set zero meltrate and watercolumn in model + watercolumn=0.; + return; + } + + /*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=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; + + /*Update water column*/ + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + if(reCast(dt)) + watercolumn+=dt*meltrate; + else + watercolumn=meltrate; + // TODO: feed meltrate & watercolumn back to model + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 15869) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 15870) @@ -398,8 +398,8 @@ IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){ if(enthalpy