Index: ../trunk-jpl/test/Archives/Archive327.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive120.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive432.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive325.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive121.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive326.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive122.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/test/Archives/Archive431.nc =================================================================== Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16589) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16590) @@ -177,6 +177,7 @@ MaterialsMuWaterEnum, MaterialsThermalExchangeVelocityEnum, MaterialsThermalconductivityEnum, + MaterialsTemperateiceconductivityEnum, MaterialsLithosphereShearModulusEnum, MaterialsLithosphereDensityEnum, MaterialsMantleShearModulusEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16589) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16590) @@ -185,6 +185,7 @@ case MaterialsMuWaterEnum : return "MaterialsMuWater"; case MaterialsThermalExchangeVelocityEnum : return "MaterialsThermalExchangeVelocity"; case MaterialsThermalconductivityEnum : return "MaterialsThermalconductivity"; + case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity"; case MaterialsLithosphereShearModulusEnum : return "MaterialsLithosphereShearModulus"; case MaterialsLithosphereDensityEnum : return "MaterialsLithosphereDensity"; case MaterialsMantleShearModulusEnum : return "MaterialsMantleShearModulus"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16589) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16590) @@ -188,6 +188,7 @@ else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; 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; else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum; @@ -258,11 +259,11 @@ else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; 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; else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"Input")==0) return InputEnum; 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; else if (strcmp(name,"Matice")==0) return MaticeEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; 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; else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16589) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16590) @@ -33,6 +33,7 @@ iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum); iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum); iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum); + iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum); iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum); iomodel->Constant(&this->beta,MaterialsBetaEnum); iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum); @@ -99,6 +100,7 @@ _printf_(" mu_water: " << mu_water << "\n"); _printf_(" heatcapacity: " << heatcapacity << "\n"); _printf_(" thermalconductivity: " << thermalconductivity << "\n"); + _printf_(" temperateiceconductivity: " << temperateiceconductivity << "\n"); _printf_(" latentheat: " << latentheat << "\n"); _printf_(" beta: " << beta << "\n"); _printf_(" meltingpoint: " << meltingpoint << "\n"); @@ -168,9 +170,12 @@ case MaterialsHeatcapacityEnum: this->heatcapacity=constant; break; - case MaterialsThermalconductivityEnum: + case MaterialsThermalconductivityEnum: this->thermalconductivity=constant; break; + case MaterialsTemperateiceconductivityEnum: + this->temperateiceconductivity=constant; + break; case MaterialsLatentheatEnum: this->latentheat=constant; break; @@ -307,6 +312,11 @@ return thermalconductivity; } /*}}}*/ +/*FUNCTION Matpar::GetTemperateIceConductivity {{{*/ +IssmDouble Matpar::GetTemperateIceConductivity(){ + return temperateiceconductivity; +} +/*}}}*/ /*FUNCTION Matpar::GetThermalExchangeVelocity {{{*/ IssmDouble Matpar::GetThermalExchangeVelocity(){ return thermal_exchange_velocity; @@ -382,11 +392,22 @@ return thermalconductivity/heatcapacity*pow(10.,-5); */ - 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; } /*}}}*/ /*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/ Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16589) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16590) @@ -21,6 +21,7 @@ IssmDouble mu_water; IssmDouble heatcapacity; IssmDouble thermalconductivity; + IssmDouble temperateiceconductivity; IssmDouble latentheat; IssmDouble beta; IssmDouble meltingpoint; @@ -108,6 +109,7 @@ IssmDouble GetThermalExchangeVelocity(); IssmDouble GetHeatCapacity(); IssmDouble GetThermalConductivity(); + IssmDouble GetTemperateIceConductivity(); IssmDouble GetLatentHeat(); IssmDouble GetBeta(); IssmDouble GetMeltingPoint(); Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16589) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16590) @@ -3814,7 +3814,7 @@ if (stabilization==2) diameter=MinEdgeLength(xyz_list); /* 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++){ gauss->GaussPoint(ig); @@ -3827,6 +3827,7 @@ enthalpy_input->GetInputValue(&enthalpy, gauss); 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; @@ -3894,7 +3895,7 @@ } 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;i