Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 20621) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 20622) @@ -1493,4 +1493,8 @@ element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); break; + case CuffeyTemperateEnum: + for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); + break; case PatersonEnum: for(i=0;i + +#include "../Numerics/types.h" + +IssmDouble CuffeyTemperate(IssmDouble temperature, IssmDouble waterfraction, IssmDouble stressexp){ + + /*output: */ + IssmDouble B; + B=Cuffey(temperature)*pow(1+181.25*max(0., min(0.01, waterfraction)), -1./stressexp); + + return B; +} Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 20621) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 20622) @@ -9,4 +9,5 @@ IssmDouble Cuffey(IssmDouble temperature); +IssmDouble CuffeyTemperate(IssmDouble temperature, IssmDouble waterfraction); IssmDouble Paterson(IssmDouble temperature); IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20621) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 20622) @@ -1037,4 +1037,5 @@ /*Rheology law (move to Material) {{{*/ CuffeyEnum, + CuffeyTemperateEnum, PatersonEnum, ArrheniusEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20621) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 20622) @@ -993,4 +993,5 @@ case OptionStructEnum : return "OptionStruct"; case CuffeyEnum : return "Cuffey"; + case CuffeyTemperateEnum : return "CuffeyTemperate"; case PatersonEnum : return "Paterson"; case ArrheniusEnum : return "Arrhenius"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20621) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 20622) @@ -1017,4 +1017,5 @@ else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; + else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; Index: /issm/trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 20621) +++ /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 20622) @@ -114,5 +114,5 @@ md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); if ismember(GiaAnalysisEnum(),analyses), @@ -145,5 +145,5 @@ fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/n)]'); fielddisplay(self,'rheology_n','Glen''s flow law exponent'); - fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); + fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]'); Index: /issm/trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 20621) +++ /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 20622) @@ -60,5 +60,5 @@ string="%s\n%s"%(string,fielddisplay(self,"rheology_B","flow law parameter [Pa/s^(1/n)]")) string="%s\n%s"%(string,fielddisplay(self,"rheology_n","Glen's flow law exponent")) - string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'Cuffey', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) string="%s\n%s"%(string,fielddisplay(self,"lithosphere_shear_modulus","Lithosphere shear modulus [Pa]")) string="%s\n%s"%(string,fielddisplay(self,"lithosphere_density","Lithosphere density [g/cm^-3]")) @@ -135,5 +135,5 @@ md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements]) - md = checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey','Paterson','Arrhenius','LliboutryDuval']) + md = checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey', 'CuffeyTemperate', 'Paterson','Arrhenius','LliboutryDuval']) md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',[1]); md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',[1]); Index: /issm/trunk-jpl/src/m/classes/matice.js =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.js (revision 20621) +++ /issm/trunk-jpl/src/m/classes/matice.js (revision 20622) @@ -76,5 +76,5 @@ fielddisplay(this,'rheology_B','flow law parameter [Pa/s^(1/n)]'); fielddisplay(this,'rheology_n',"Glen's flow law exponent"); - fielddisplay(this,'rheology_law',"law for the temperature dependance of the rheology: 'None', 'Cuffey', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"); + fielddisplay(this,'rheology_law',"law for the temperature dependance of the rheology: 'None', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"); fielddisplay(this,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); fielddisplay(this,'lithosphere_density','Lithosphere density [g/cm^-3]'); @@ -94,5 +94,5 @@ checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements,1]); - checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey','Paterson','Arrhenius','LliboutryDuval']); + checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval']); if(ArrayAnyEqual(ArrayIsMember(GiaAnalysisEnum(),analyses),1)){ Index: /issm/trunk-jpl/src/m/classes/matice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.m (revision 20621) +++ /issm/trunk-jpl/src/m/classes/matice.m (revision 20622) @@ -114,5 +114,5 @@ md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); if ismember(GiaAnalysisEnum(),analyses), @@ -144,5 +144,5 @@ fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/n)]'); fielddisplay(self,'rheology_n','Glen''s flow law exponent'); - fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); + fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]'); Index: /issm/trunk-jpl/src/m/classes/matice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.py (revision 20621) +++ /issm/trunk-jpl/src/m/classes/matice.py (revision 20622) @@ -61,5 +61,5 @@ string="%s\n%s"%(string,fielddisplay(self,"rheology_B","flow law parameter [Pa/s^(1/n)]")) string="%s\n%s"%(string,fielddisplay(self,"rheology_n","Glen's flow law exponent")) - string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'Cuffey', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) + string="%s\n%s"%(string,fielddisplay(self,"rheology_law","law for the temperature dependance of the rheology: 'None', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) string="%s\n%s"%(string,fielddisplay(self,"lithosphere_shear_modulus","Lithosphere shear modulus [Pa]")) string="%s\n%s"%(string,fielddisplay(self,"lithosphere_density","Lithosphere density [g/cm^-3]")) @@ -136,5 +136,5 @@ md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1) md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements]) - md = checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey','Paterson','Arrhenius','LliboutryDuval']) + md = checkfield(md,'fieldname','materials.rheology_law','values',['None','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval']) md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',[1]); md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',[1]); Index: /issm/trunk-jpl/src/m/enum/CuffeyTemperateEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/CuffeyTemperateEnum.m (revision 20622) +++ /issm/trunk-jpl/src/m/enum/CuffeyTemperateEnum.m (revision 20622) @@ -0,0 +1,11 @@ +function macro=CuffeyTemperateEnum() +%CUFFEYTEMPERATEENUM - Enum of CuffeyTemperate +% +% 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=CuffeyTemperateEnum() + +macro=StringToEnum('CuffeyTemperate'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.js =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.js (revision 20621) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.js (revision 20622) @@ -978,33 +978,34 @@ function OptionStructEnum(){ return 974;} function CuffeyEnum(){ return 975;} -function PatersonEnum(){ return 976;} -function ArrheniusEnum(){ return 977;} -function LliboutryDuvalEnum(){ return 978;} -function SpclevelsetEnum(){ return 979;} -function ExtrapolationVariableEnum(){ return 980;} -function IceMaskNodeActivationEnum(){ return 981;} -function LevelsetfunctionSlopeXEnum(){ return 982;} -function LevelsetfunctionSlopeYEnum(){ return 983;} -function LevelsetfunctionPicardEnum(){ return 984;} -function LevelsetReinitFrequencyEnum(){ return 985;} -function SealevelriseSolutionEnum(){ return 986;} -function SealevelriseAnalysisEnum(){ return 987;} -function SealevelEnum(){ return 988;} -function SealevelEustaticEnum(){ return 989;} -function SealevelriseDeltathicknessEnum(){ return 990;} -function SealevelriseMaxiterEnum(){ return 991;} -function SealevelriseReltolEnum(){ return 992;} -function SealevelriseAbstolEnum(){ return 993;} -function SealevelriseLoveHEnum(){ return 994;} -function SealevelriseLoveKEnum(){ return 995;} -function SealevelriseTideLoveHEnum(){ return 996;} -function SealevelriseTideLoveKEnum(){ return 997;} -function SealevelriseRigidEnum(){ return 998;} -function SealevelriseElasticEnum(){ return 999;} -function SealevelriseRotationEnum(){ return 1000;} -function SealevelriseGElasticEnum(){ return 1001;} -function SealevelriseDegaccEnum(){ return 1002;} -function SealevelriseTransitionsEnum(){ return 1003;} -function SealevelriseRequestedOutputsEnum(){ return 1004;} -function SealevelriseNumRequestedOutputsEnum(){ return 1005;} -function MaximumNumberOfDefinitionsEnum(){ return 1006;} +function CuffeyTemperateEnum(){ return 976;} +function PatersonEnum(){ return 977;} +function ArrheniusEnum(){ return 978;} +function LliboutryDuvalEnum(){ return 979;} +function SpclevelsetEnum(){ return 980;} +function ExtrapolationVariableEnum(){ return 981;} +function IceMaskNodeActivationEnum(){ return 982;} +function LevelsetfunctionSlopeXEnum(){ return 983;} +function LevelsetfunctionSlopeYEnum(){ return 984;} +function LevelsetfunctionPicardEnum(){ return 985;} +function LevelsetReinitFrequencyEnum(){ return 986;} +function SealevelriseSolutionEnum(){ return 987;} +function SealevelriseAnalysisEnum(){ return 988;} +function SealevelEnum(){ return 989;} +function SealevelEustaticEnum(){ return 990;} +function SealevelriseDeltathicknessEnum(){ return 991;} +function SealevelriseMaxiterEnum(){ return 992;} +function SealevelriseReltolEnum(){ return 993;} +function SealevelriseAbstolEnum(){ return 994;} +function SealevelriseLoveHEnum(){ return 995;} +function SealevelriseLoveKEnum(){ return 996;} +function SealevelriseTideLoveHEnum(){ return 997;} +function SealevelriseTideLoveKEnum(){ return 998;} +function SealevelriseRigidEnum(){ return 999;} +function SealevelriseElasticEnum(){ return 1000;} +function SealevelriseRotationEnum(){ return 1001;} +function SealevelriseGElasticEnum(){ return 1002;} +function SealevelriseDegaccEnum(){ return 1003;} +function SealevelriseTransitionsEnum(){ return 1004;} +function SealevelriseRequestedOutputsEnum(){ return 1005;} +function SealevelriseNumRequestedOutputsEnum(){ return 1006;} +function MaximumNumberOfDefinitionsEnum(){ return 1007;} Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 20621) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 20622) @@ -985,4 +985,5 @@ def OptionStructEnum(): return StringToEnum("OptionStruct")[0] def CuffeyEnum(): return StringToEnum("Cuffey")[0] +def CuffeyTemperateEnum(): return StringToEnum("CuffeyTemperate")[0] def PatersonEnum(): return StringToEnum("Paterson")[0] def ArrheniusEnum(): return StringToEnum("Arrhenius")[0] Index: /issm/trunk-jpl/src/m/materials/cuffeytemperate.m =================================================================== --- /issm/trunk-jpl/src/m/materials/cuffeytemperate.m (revision 20622) +++ /issm/trunk-jpl/src/m/materials/cuffeytemperate.m (revision 20622) @@ -0,0 +1,21 @@ +function rigidity=cuffeytemperate(temperature, waterfraction, stressexp) +%CUFFEYTEMPERATE - calculates ice rigidity as a function of temperature and waterfraction +% +% rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) +% (Cuffey and Paterson, p75). +% temperature is in Kelvin degrees +% +% Usage: +% rigidity=cuffeytemperate(temperature, waterfraction, stressexp) + +if (any(size(temperature)~=size(waterfraction))), + error('input temperature and waterfraction should have same size!'); +end +if any(temperature<0) + error('input temperature should be in Kelvin (positive)'); +end +if any(waterfraction<0 | waterfraction>1) + error('input waterfraction should be between 0 and 1'); +end + +rigidity=reshape(cuffey(temperature), size(temperature)).*(1+181.25*max(0., min(0.01, waterfraction))).^(-1/stressexp); Index: /issm/trunk-jpl/src/m/materials/cuffeytemperate.py =================================================================== --- /issm/trunk-jpl/src/m/materials/cuffeytemperate.py (revision 20622) +++ /issm/trunk-jpl/src/m/materials/cuffeytemperate.py (revision 20622) @@ -0,0 +1,29 @@ +import numpy +import cuffey + +def cuffeytemperate(temperature, waterfraction, stressexp) + + """ + CUFFEYTEMPERATE - calculates ice rigidity as a function of temperature and waterfraction + + rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) + (Cuffey and Paterson, p75). + temperature is in Kelvin degrees + + Usage: + rigidity=cuffeytemperate(temperature, waterfraction, stressexp) + """ + + if numpy.any(temperature<0.): + raise RuntimeError("input temperature should be in Kelvin (positive)") + + if (numpy.any(temperature.shape~=waterfraction.shape)), + error('input temperature and waterfraction should have same size!'); + end + if numpy.any(waterfraction<0 | waterfraction>1) + error('input waterfraction should be between 0 and 1'); + end + + rigidity=numpy.multiply(cuffey(temperature).reshape(temperature.shape), (1*numpy.ones(waterfraction.shape)+181.25*numpy.maximum(numpy.zeros(waterfraction.shape), numpy.minimum(0.01*numpy.ones(waterfraction.shape), waterfraction)))**(-1/stressexp)); + + return rigidity