Index: ../trunk-jpl/src/m/classes/frictionjosh.m =================================================================== --- ../trunk-jpl/src/m/classes/frictionjosh.m (revision 21550) +++ ../trunk-jpl/src/m/classes/frictionjosh.m (revision 21551) @@ -6,12 +6,12 @@ classdef frictionjosh properties (SetAccess=public) coefficient = NaN; - temperature = NaN; + pressure_adjusted_temperature = NaN; end methods function self = extrude(self,md) % {{{ self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1); - self.temperature=project3d(md,'vector',self.temperature,'type','node','layer',1); + self.pressure_adjusted_temperature=project3d(md,'vector',self.pressure_adjusted_temperature,'type','node','layer',1); end % }}} function self = frictionjosh(varargin) % {{{ switch nargin @@ -31,8 +31,8 @@ %Early return if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','friction.temperature','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.coefficient','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','friction.pressure_adjusted_temperature','NaN',1,'Inf',1); %Check that temperature is provided md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); @@ -40,13 +40,13 @@ function disp(self) % {{{ disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)')); fielddisplay(self,'coefficient','friction coefficient [SI]'); - fielddisplay(self,'temperature','friction temperature [K]'); + fielddisplay(self,'pressure_adjusted_temperature','friction pressure_adjusted_temperature (T - Tpmp) [K]'); end % }}} function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer'); WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end % }}} end end Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21550) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21551) @@ -2032,13 +2032,17 @@ yearlytemperatures[2] = s[2]; this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); - /*Convert that to enthalpy for the enthalpy model*/ - IssmDouble enthalpy[6]; - GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); - ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.); - ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.); - ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.); - this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum)); + bool isenthalpy; + this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + if(isenthalpy){ + /*Convert that to enthalpy for the enthalpy model*/ + IssmDouble enthalpy[6]; + GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); + ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.); + ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.); + ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.); + this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum)); + } } this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21550) +++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 21551) @@ -246,6 +246,9 @@ case 8: GetAlpha2Sommers(palpha2,gauss); break; + case 9: + GetAlpha2Josh(palpha2,gauss); + break; default: _error_("Friction law "<< this->law <<" not supported"); } @@ -462,6 +465,47 @@ /*Assign output pointers:*/ *palpha2=alpha2; }/*}}}*/ +void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + /*Here, we want to parameterize the friction as a function of temperature + * + * alpha2 = alpha2_viscous * 1/f(T) + * + * where f(T) = exp((T-Tpmp)/gamma) + */ + + /*Intermediaries: */ + IssmDouble T,Tpmp,deltaT,deltaTref,pressure; + IssmDouble alpha2,time,gamma; + const IssmDouble yts = 365*24*3600.; + + /*Get viscous part*/ + this->GetAlpha2Viscous(&alpha2,gauss); + + /*Get delta Refs*/ + element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum); + + /*Compute delta T*/ + element->GetInputValue(&T,gauss,TemperatureEnum); + element->GetInputValue(&pressure,gauss,PressureEnum); + Tpmp = element->TMeltingPoint(pressure); + deltaT = T-Tpmp; + + /*Compute gamma*/ + element->parameters->FindParam(&time,TimeEnum); + if(time<25e3*yts){ + gamma = 10.; + } + else{ + gamma = 5.; + } + gamma = 5.; + + /*Compute scaling parameter*/ + alpha2 = alpha2 * exp((deltaTref - deltaT)/(2*gamma)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ /*This routine calculates the basal friction coefficient Index: ../trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 21550) +++ ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 21551) @@ -35,6 +35,7 @@ void GetAlpha2(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Sommers(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss); Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21550) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 21551) @@ -822,6 +822,10 @@ iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum); iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum); break; + case 9: + iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum); + iomodel->FetchDataToInput(elements,"md.friction.pressure_adjusted_temperature",FrictionPressureAdjustedTemperatureEnum); + break; default: _error_("not supported"); } Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21550) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21551) @@ -108,6 +108,7 @@ FrictionAsEnum, FrictionCoefficientEnum, FrictionCoefficientcoulombEnum, + FrictionPressureAdjustedTemperatureEnum, FrictionPEnum, FrictionQEnum, FrictionMEnum, @@ -273,6 +274,9 @@ MaterialsMantleDensityEnum, MaterialsEarthDensityEnum, MeshAverageVertexConnectivityEnum, + MeshXEnum, + MeshYEnum, + MeshZEnum, MeshElementsEnum, MeshNumberofelementsEnum, MeshNumberoflayersEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21550) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21551) @@ -114,6 +114,7 @@ case FrictionAsEnum : return "FrictionAs"; case FrictionCoefficientEnum : return "FrictionCoefficient"; case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb"; + case FrictionPressureAdjustedTemperatureEnum : return "FrictionPressureAdjustedTemperature"; case FrictionPEnum : return "FrictionP"; case FrictionQEnum : return "FrictionQ"; case FrictionMEnum : return "FrictionM"; @@ -279,6 +280,9 @@ case MaterialsMantleDensityEnum : return "MaterialsMantleDensity"; case MaterialsEarthDensityEnum : return "MaterialsEarthDensity"; case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity"; + case MeshXEnum : return "MeshX"; + case MeshYEnum : return "MeshY"; + case MeshZEnum : return "MeshZ"; case MeshElementsEnum : return "MeshElements"; case MeshNumberofelementsEnum : return "MeshNumberofelements"; case MeshNumberoflayersEnum : return "MeshNumberoflayers"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21550) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21551) @@ -114,6 +114,7 @@ else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum; else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum; else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum; + else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; else if (strcmp(name,"FrictionM")==0) return FrictionMEnum; @@ -135,11 +136,11 @@ else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; else if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum; - else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; + if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; + else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum; else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum; @@ -258,11 +259,11 @@ else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum; else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum; - else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; + if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; + else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum; else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; @@ -285,6 +286,9 @@ else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum; else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum; else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum; + else if (strcmp(name,"MeshX")==0) return MeshXEnum; + else if (strcmp(name,"MeshY")==0) return MeshYEnum; + else if (strcmp(name,"MeshZ")==0) return MeshZEnum; else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; else if (strcmp(name,"MeshNumberofelements")==0) return MeshNumberofelementsEnum; else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum; @@ -378,14 +382,14 @@ else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; - else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; + else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; @@ -501,14 +505,14 @@ else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; - else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; else if (strcmp(name,"Temperature")==0) return TemperatureEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; + else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; @@ -624,14 +628,14 @@ else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; - else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; + else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum; else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum; else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum; @@ -747,14 +751,14 @@ else if (strcmp(name,"Sset")==0) return SsetEnum; else if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; - else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"Gsl")==0) return GslEnum; + else if (strcmp(name,"Gsl")==0) return GslEnum; else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; @@ -870,14 +874,14 @@ else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else if (strcmp(name,"OptionCell")==0) return OptionCellEnum; else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; - else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; + else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum; else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; @@ -993,14 +997,14 @@ else if (strcmp(name,"Melange")==0) return MelangeEnum; else if (strcmp(name,"Water")==0) return WaterEnum; else if (strcmp(name,"DataSet")==0) return DataSetEnum; - else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Nodes")==0) return NodesEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Contours")==0) return ContoursEnum; + else if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"Parameters")==0) return ParametersEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"Results")==0) return ResultsEnum;