Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17027) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17028) @@ -391,7 +391,5 @@ #}}} #Thermal sources {{{ -thermal_sources = ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\ - ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\ - ./analyses/ThermalAnalysis.h\ +thermal_sources = ./analyses/ThermalAnalysis.h\ ./analyses/ThermalAnalysis.cpp\ ./analyses/EnthalpyAnalysis.h\ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17027) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17028) @@ -271,10 +271,4 @@ virtual void UpdateConstraintsExtrudeFromTop(void)=0; - #ifdef _HAVE_THERMAL_ - virtual void UpdateBasalConstraintsEnthalpy(void)=0; - virtual void ComputeBasalMeltingrate(void)=0; - virtual void DrainWaterfraction(IssmDouble* drainrate_element)=0; - #endif - #ifdef _HAVE_HYDROLOGY_ virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17027) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17028) @@ -3282,287 +3282,4 @@ #endif -#ifdef _HAVE_THERMAL_ -/*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/ -void Penta::UpdateBasalConstraintsEnthalpy(void){ - - /*Intermediary*/ - bool isenthalpy,isdynamicbasalspc,setspc; - int numindices, numindicesup; - IssmDouble pressure, pressureup; - IssmDouble h_pmp, enthalpy, enthalpyup; - IssmDouble watercolumn; - int *indices = NULL, *indicesup = NULL; - - /* Only update Constraints at the base of grounded ice*/ - if(!IsOnBed() || IsFloating()) return; - - /*Check wether dynamic basal boundary conditions are activated */ - parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); - if(!isenthalpy) return; - parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); - if(!isdynamicbasalspc) return; - - /*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); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); //_assert_(watercolumn_input); - - /*if there is a temperate layer of zero thickness, set spc 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]); - - /*Check wether there is a temperate layer at the base or not */ - /*check if node is temperate, if not, continue*/ - enthalpy_input->GetInputValue(&enthalpy, gauss); - pressure_input->GetInputValue(&pressure, gauss); - watercolumn_input->GetInputValue(&watercolumn,gauss); - setspc = false; - if (enthalpy>=matpar->PureIceEnthalpy(pressure)){ - /*check if upper node is temperate, too. - if yes, then we have a temperate layer of positive thickness and reset the spc. - if not, apply dirichlet BC.*/ - enthalpy_input->GetInputValue(&enthalpyup, gaussup); - pressure_input->GetInputValue(&pressureup, gaussup); - setspc=((enthalpyupPureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; - } - - if (setspc) { - /*Calculate enthalpy at pressure melting point */ - 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(1,h_pmp); - } - else { - /*remove spc*/ - nodes[indices[i]]->DofInFSet(0); - } - } - - /*Free ressources:*/ - xDelete(indices); - xDelete(indicesup); - delete gauss; - delete gaussup; -} -/*}}}*/ -/*FUNCTION Penta::ComputeBasalMeltingrate{{{*/ -void Penta::ComputeBasalMeltingrate(void){ - /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/ - /* melting rate is positive when melting, negative when refreezing*/ - - /* Intermediaries */ - bool isenthalpy, checkpositivethickness, istemperatelayer; - int i,j,iv,analysis_type; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble heatflux,kappa; - IssmDouble vec_heatflux[3]; - IssmDouble normal_base[3], d1enthalpy[3]; - IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES]; - IssmDouble enthalpy[NUMVERTICES],pressure[NUMVERTICES]; - IssmDouble temperature, waterfraction; - IssmDouble latentheat, rho_ice; - IssmDouble basalfriction, alpha2; - IssmDouble vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES]; - IssmDouble geothermalflux[NUMVERTICES]; - IssmDouble dt, yts; - 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*/ - if(!IsOnBed() || IsFloating()) return; - - /*Check wether enthalpy is activated*/ - parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); - if(!isenthalpy) return; - - /*Fetch parameters and inputs */ - latentheat=matpar->GetLatentHeat(); - rho_ice=matpar->GetRhoIce(); - GetInputListOnVertices(&vx[0],VxEnum); - GetInputListOnVertices(&vy[0],VyEnum); - GetInputListOnVertices(&vz[0],VzEnum); - GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); - GetInputListOnVertices(&pressure[0],PressureEnum); - GetInputListOnVertices(&watercolumn[0],WatercolumnEnum); - GetInputListOnVertices(&basalmeltingrate[0],BasalforcingsMeltingRateEnum); - GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_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); - - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - friction=new Friction(this,3); - - /******** MELTING RATES ************************************/ - GaussPenta* gauss=new GaussPenta(); - for(iv=0;ivGaussVertex(iv); - checkpositivethickness=true; - - _assert_(watercolumn[iv]>=0.); - - /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/ - meltingrate_enthalpy[iv]=0.; - heating[iv]=0.; - if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure[iv]))){ - /*ensure that no ice is at TPureIceEnthalpy(pressure[iv]); - } - else if(enthalpy[iv]PureIceEnthalpy(pressure[iv])){ - /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/ - checkpositivethickness=false; // cold base, skip next test - } - else {/*we have a temperate base, go to next test*/} - - if(checkpositivethickness){ - /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */ - istemperatelayer=false; - 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{ - enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,&enthalpy[0],&pressure[0]); _assert_(kappa>0.); - for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; - } - - /*geothermal heatflux*/ - NormalBase(&normal_base[0],&xyz_list_tria[0][0]); - heatflux=0.; - for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; - - /*basal friction*/ - friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); - basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0)); - - matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]); - // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66 - 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->GetUpperPenta(); - } - // add drained water to melting rate - for(iv=0; ivparameters->FindParam(&dt,TimesteppingTimeStepEnum); - for(iv=0;iv(dt)){ - 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[iv]; - } - else{ - basalmeltingrate[iv]=meltingrate_enthalpy[iv]; - watercolumn[iv]+=dt*meltingrate_enthalpy[iv]; - } - } - else{ - 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)); - this->inputs->AddInput(new PentaInput(WatercolumnEnum,watercolumn,P1Enum)); - this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum)); - - /*Clean up and return*/ - delete gauss; - delete friction; -} -/*}}}*/ -/*FUNCTION Penta::DrainWaterfraction{{{*/ -void Penta::DrainWaterfraction(IssmDouble* drainrate_element){ - - /*Intermediaries*/ - bool isenthalpy; - int iv, index0; - IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES]; - IssmDouble dw[NUMVERTICES]; - IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; - IssmDouble dt, height_element; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble rho_water, rho_ice; - - /*Check wether enthalpy is activated*/ - parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); - if(!isenthalpy) return; - - 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); - - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - for(iv=0;ivEnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); - 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)); - - /*return meltwater column equivalent to drained water*/ - for(iv=0;iv