Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16589) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16590) @@ -3815,5 +3815,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(2,2); + gauss=new GaussPenta(3,3); for(int ig=gauss->begin();igend();ig++){ @@ -3828,4 +3828,5 @@ pressure_input->GetInputValue(&pressure, gauss); kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.); + D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice; if(reCast(dt)) D_scalar_conduct=D_scalar_conduct*dt; @@ -3895,5 +3896,5 @@ else if(stabilization==2){ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); + tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice); for(i=0;iConstant(&this->heatcapacity,MaterialsHeatcapacityEnum); iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum); + iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum); iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum); iomodel->Constant(&this->beta,MaterialsBetaEnum); @@ -100,4 +101,5 @@ _printf_(" heatcapacity: " << heatcapacity << "\n"); _printf_(" thermalconductivity: " << thermalconductivity << "\n"); + _printf_(" temperateiceconductivity: " << temperateiceconductivity << "\n"); _printf_(" latentheat: " << latentheat << "\n"); _printf_(" beta: " << beta << "\n"); @@ -169,6 +171,9 @@ this->heatcapacity=constant; break; - case MaterialsThermalconductivityEnum: + case MaterialsThermalconductivityEnum: this->thermalconductivity=constant; + break; + case MaterialsTemperateiceconductivityEnum: + this->temperateiceconductivity=constant; break; case MaterialsLatentheatEnum: @@ -306,4 +311,9 @@ IssmDouble Matpar::GetThermalConductivity(){ return thermalconductivity; +} +/*}}}*/ +/*FUNCTION Matpar::GetTemperateIceConductivity {{{*/ +IssmDouble Matpar::GetTemperateIceConductivity(){ + return temperateiceconductivity; } /*}}}*/ @@ -383,9 +393,20 @@ */ - IssmDouble eps=0.05*heatcapacity; + IssmDouble eps=0.1*heatcapacity; IssmDouble hpmp=PureIceEnthalpy(pressure); IssmDouble kappa_c=thermalconductivity/heatcapacity; - IssmDouble kappa_t=thermalconductivity/heatcapacity*pow(10.,-1); - return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c; + IssmDouble kappa_t=temperateiceconductivity/heatcapacity; + + if(enthalpy<=hpmp-eps) + return kappa_c; + else if(enthalpy>=hpmp+eps) + return kappa_t; + else { + IssmDouble xi=enthalpy-hpmp; + IssmDouble pi=3.141592653589793238462643; + return kappa_c + (kappa_t-kappa_c)*((xi+eps)/(2*eps) + sin(pi*xi/eps)/(2*pi)); + } + + //return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16589) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16590) @@ -22,4 +22,5 @@ IssmDouble heatcapacity; IssmDouble thermalconductivity; + IssmDouble temperateiceconductivity; IssmDouble latentheat; IssmDouble beta; @@ -109,4 +110,5 @@ IssmDouble GetHeatCapacity(); IssmDouble GetThermalConductivity(); + IssmDouble GetTemperateIceConductivity(); IssmDouble GetLatentHeat(); IssmDouble GetBeta(); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16589) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16590) @@ -178,4 +178,5 @@ MaterialsThermalExchangeVelocityEnum, MaterialsThermalconductivityEnum, + MaterialsTemperateiceconductivityEnum, MaterialsLithosphereShearModulusEnum, MaterialsLithosphereDensityEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16589) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16590) @@ -186,4 +186,5 @@ case MaterialsThermalExchangeVelocityEnum : return "MaterialsThermalExchangeVelocity"; case MaterialsThermalconductivityEnum : return "MaterialsThermalconductivity"; + case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity"; case MaterialsLithosphereShearModulusEnum : return "MaterialsLithosphereShearModulus"; case MaterialsLithosphereDensityEnum : return "MaterialsLithosphereDensity"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16589) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16590) @@ -189,4 +189,5 @@ else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum; else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; + else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum; else if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum; else if (strcmp(name,"MaterialsLithosphereDensity")==0) return MaterialsLithosphereDensityEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; - else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum; + if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; + else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; - else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"IntParam")==0) return IntParamEnum; + if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; + else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; - else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; + else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; Index: /issm/trunk-jpl/src/m/classes/matice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.m (revision 16589) +++ /issm/trunk-jpl/src/m/classes/matice.m (revision 16590) @@ -13,4 +13,5 @@ latentheat = 0.; thermalconductivity = 0.; + temperateiceconductivity = 0.; meltingpoint = 0.; beta = 0.; @@ -69,4 +70,7 @@ %ice thermal conductivity (W/m/K) obj.thermalconductivity=2.4; + + %wet ice thermal conductivity (W/m/K) + obj.temperateiceconductivity=.24; %the melting point of ice at 1 atmosphere of pressure in K @@ -118,5 +122,6 @@ fielddisplay(obj,'mu_water','water viscosity [N s/m^2]'); fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]'); - fielddisplay(obj,'thermalconductivity','ice thermal conductivity [W/m/K]'); + fielddisplay(obj,'thermalconductivity',['ice thermal conductivity [W/m/K]']); + fielddisplay(obj,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]'); fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K'); fielddisplay(obj,'latentheat','latent heat of fusion [J/m^3]'); @@ -141,4 +146,5 @@ WriteData(fid,'object',obj,'class','materials','fieldname','latentheat','format','Double'); WriteData(fid,'object',obj,'class','materials','fieldname','thermalconductivity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','temperateiceconductivity','format','Double'); WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double'); WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double'); Index: /issm/trunk-jpl/src/m/classes/matice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.py (revision 16589) +++ /issm/trunk-jpl/src/m/classes/matice.py (revision 16590) @@ -21,4 +21,5 @@ self.latentheat = 0. self.thermalconductivity = 0. + self.temperateiceconductivity = 0. self.meltingpoint = 0. self.beta = 0. @@ -46,4 +47,5 @@ string="%s\n%s"%(string,fielddisplay(self,"heatcapacity","heat capacity [J/kg/K]")) string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]")) string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K")) string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]")) @@ -82,4 +84,7 @@ #ice thermal conductivity (W/m/K) self.thermalconductivity=2.4 + + #temperate ice thermal conductivity (W/m/K) + self.temperateiceconductivity=0.24 #the melting point of ice at 1 atmosphere of pressure in K @@ -130,4 +135,5 @@ WriteData(fid,'object',self,'class','materials','fieldname','latentheat','format','Double') WriteData(fid,'object',self,'class','materials','fieldname','thermalconductivity','format','Double') + WriteData(fid,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double') WriteData(fid,'object',self,'class','materials','fieldname','meltingpoint','format','Double') WriteData(fid,'object',self,'class','materials','fieldname','beta','format','Double') Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16589) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16590) @@ -178,4 +178,5 @@ def MaterialsThermalExchangeVelocityEnum(): return StringToEnum("MaterialsThermalExchangeVelocity")[0] def MaterialsThermalconductivityEnum(): return StringToEnum("MaterialsThermalconductivity")[0] +def MaterialsTemperateiceconductivityEnum(): return StringToEnum("MaterialsTemperateiceconductivity")[0] def MaterialsLithosphereShearModulusEnum(): return StringToEnum("MaterialsLithosphereShearModulus")[0] def MaterialsLithosphereDensityEnum(): return StringToEnum("MaterialsLithosphereDensity")[0] Index: /issm/trunk-jpl/src/m/enum/MaterialsTemperateiceconductivityEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaterialsTemperateiceconductivityEnum.m (revision 16590) +++ /issm/trunk-jpl/src/m/enum/MaterialsTemperateiceconductivityEnum.m (revision 16590) @@ -0,0 +1,11 @@ +function macro=MaterialsTemperateiceconductivityEnum() +%MATERIALSTEMPERATEICECONDUCTIVITYENUM - Enum of MaterialsTemperateiceconductivity +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=MaterialsTemperateiceconductivityEnum() + +macro=StringToEnum('MaterialsTemperateiceconductivity');