Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 24059) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 24060) @@ -174,4 +174,6 @@ ./shared/Elements/Paterson.cpp\ ./shared/Elements/Arrhenius.cpp\ + ./shared/Elements/NyeCO2.cpp\ + ./shared/Elements/NyeH2O.cpp\ ./shared/Elements/LliboutryDuval.cpp \ ./shared/Elements/PrintArrays.cpp\ Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24059) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24060) @@ -1606,4 +1606,12 @@ element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); break; + case NyeH2OEnum: + for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); + break; + case NyeCO2Enum: + for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); + break; case ArrheniusEnum:{ element->GetVerticesCoordinates(&xyz_list); Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 24059) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 24060) @@ -874,6 +874,14 @@ element->AddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); break; + case NyeH2OEnum: + for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); + break; + case NyeCO2Enum: + for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); + break; case ArrheniusEnum:{ - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); for(i=0;iAddInput(MaterialsRheologyBEnum,&B[0],element->GetElementType()); Index: /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp (revision 24060) +++ /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp (revision 24060) @@ -0,0 +1,124 @@ +/* \file NyeCO2.cpp + * \ brief figure out B of CO2 ice for a certain temperature + * INPUT function B=NyeCO2(temperature) + * where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Nye, p -COMPLETE- ). + + VALUES CALCULATED USING THE FOLLOWING MATLAB CODE + + Rg = 8.3144598; % J mol^-1 K^-1 + A_const = 10^(10.8); % s^-1 MPa + Q = 63000; % J mol^-1 + n = 7; % Glen's exponent + T = 97.5:5:252.5; % K + A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa + B = A.^(-1/n)*1e6; % s^(1/n) Pa + */ + +#include +#include "../Numerics/types.h" + +IssmDouble NyeCO2(IssmDouble temperature){ + + IssmDouble B,T; + + T=temperature; + + if(T<=100.0){ + B=1899806746.55828; + } + else if((T>=100.0) && (T<=105.0)){ + B=1105378610.10266; + } + else if((T>=105.0) && (T<=110.0)){ + B=676381320.528217; + } + else if((T>=110.0) && (T<=115.0)){ + B=432348426.904519; + } + else if((T>=115.0) && (T<=120.0)){ + B=287089456.910568; + } + else if((T>=120.0) && (T<=125.0)){ + B=197113445.038555; + } + else if((T>=125.0) && (T<=130.0)){ + B=139387289.896988; + } + else if((T>=130.0) && (T<=135.0)){ + B=101178460.700527; + } + else if((T>=135.0) && (T<=140.0)){ + B=75174730.2531849; + } + else if((T>=140.0) && (T<=145.0)){ + B=57030799.0697205; + } + else if((T>=145.0) && (T<=150.0)){ + B=44083907.2876590; + } + else if((T>=150.0) && (T<=155.0)){ + B=34656426.3665295; + } + else if((T>=155.0) && (T<=160.0)){ + B=27664457.7106089; + } + else if((T>=160.0) && (T<=165.0)){ + B=22391479.9130613; + } + else if((T>=165.0) && (T<=170.0)){ + B=18353816.3612526; + } + else if((T>=170.0) && (T<=175.0)){ + B=15218650.2018789; + } + else if((T>=175.0) && (T<=180.0)){ + B=12752901.3547654; + } + else if((T>=180.0) && (T<=185.0)){ + B=10790666.8437103; + } + else if((T>=185.0) && (T<=190.0)){ + B=9212075.16184763; + } + else if((T>=190.0) && (T<=195.0)){ + B=7929303.02671413; + } + else if((T>=195.0) && (T<=200.0)){ + B=6877172.17359016; + } + else if((T>=200.0) && (T<=205.0)){ + B=6006726.68822080; + } + else if((T>=205.0) && (T<=210.0)){ + B=5280781.91812612; + } + else if((T>=210.0) && (T<=215.0)){ + B=4670797.58917498; + } + else if((T>=215.0) && (T<=220.0)){ + B=4154653.09694929; + } + else if((T>=220.0) && (T<=225.0)){ + B=3715045.71619004; + } + else if((T>=225.0) && (T<=230.0)){ + B=3338324.34308717; + } + else if((T>=230.0) && (T<=235.0)){ + B=3013631.36885997; + } + else if((T>=235.0) && (T<=240.0)){ + B=2732264.98697388; + } + else if((T>=240.0) && (T<=245.0)){ + B=2487200.85823651; + } + else if((T>=245.0) && (T<=250.0)){ + B=2272730.12712433; + } + else{ + B=2084183.18884150; + } + + return B; +} Index: /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp (revision 24060) +++ /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp (revision 24060) @@ -0,0 +1,95 @@ +/* \file NyeH2O.cpp + * \ brief figure out B of H2O ice for a certain temperature + * INPUT function B=NyeH2O(temperature) + * where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Nye, p -COMPLETE- ). + + VALUES CALCULATED USING THE FOLLOWING MATLAB CODE + + Rg = 8.3144598; % J mol^-1 K^-1 + A_const = 9e4; % s^-1 MPa + Q = 60000; % J mol^-1 + n = 3; % Glen's exponent + T = -102.5:5:2.5; T = T + 273; % K + A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa + B = A.^(-1/n)*1e6; % s^(1/n) Pa + */ + +#include +#include "../Numerics/types.h" + +IssmDouble NyeH2O(IssmDouble temperature){ + + IssmDouble B,T; + + /*Switch to celsius from Kelvin: */ + T=temperature-273.15; + + if(T<=-100){ + B=29901899089.3257; + } + else if((T>=-100.0) && (T<=-95.0)){ + B=20004935214.6693; + } + else if((T>=-95.0) && (T<=-90.0)){ + B=13685054490.7813; + } + else if((T>=-90.0) && (T<=-85.0)){ + B=9555312498.71763; + } + else if((T>=-85.0) && (T<=-80.0)){ + B=6798800864.99471; + } + else if((T>=-80.0) && (T<=-75.0)){ + B=4922440695.83365; + } + else if((T>=-75.0) && (T<=-70.0)){ + B=3621794191.06010; + } + else if((T>=-70.0) && (T<=-65.0)){ + B=2704902125.98878; + } + else if((T>=-65.0) && (T<=-60.0)){ + B=2048338724.47737; + } + else if((T>=-60.0) && (T<=-55.0)){ + B=1571285974.28101; + } + else if((T>=-55.0) && (T<=-50.0)){ + B=1219918332.04930; + } + else if((T>=-50.0) && (T<=-45.0)){ + B=957813559.547873; + } + else if((T>=-45.0) && (T<=-40.0)){ + B=759956402.731125; + } + else if((T>=-40.0) && (T<=-35.0)){ + B=608924663.484756; + } + else if((T>=-35.0) && (T<=-30.0)){ + B=492424321.822968; + } + else if((T>=-30.0) && (T<=-25.0)){ + B=401672414.952195; + } + else if((T>=-25.0) && (T<=-20.0)){ + B=330320977.865296; + } + else if((T>=-20.0) && (T<=-15.0)){ + B=273731390.267536; + } + else if((T>=-15.0) && (T<=-10.0)){ + B=228478811.983490; + } + else if((T>=-10.0) && (T<=-5.0)){ + B=192009688.386368; + } + else if((T>=-5.0) && (T<=0.0)){ + B=162402355.392225; + } + else{ + B=138197905.740112; + } + + return B; +} Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 24059) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 24060) @@ -13,4 +13,6 @@ IssmDouble Paterson(IssmDouble temperature); IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); +IssmDouble NyeH2O(IssmDouble temperature); +IssmDouble NyeCO2(IssmDouble temperature); IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat); // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24059) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24060) @@ -1124,4 +1124,6 @@ syn keyword cConstant NoneEnum syn keyword cConstant NumberedcostfunctionEnum +syn keyword cConstant NyeCO2Enum +syn keyword cConstant NyeH2OEnum syn keyword cConstant NumericalfluxEnum syn keyword cConstant OldGradientEnum @@ -1274,5 +1276,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1281,6 +1282,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1295,6 +1296,6 @@ syn keyword cType ElementHook syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1303,10 +1304,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1323,5 +1323,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1334,5 +1333,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1347,6 +1345,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1361,10 +1359,10 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1386,4 +1384,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1394,5 +1396,4 @@ syn keyword cType BalancethicknessSoftAnalysis syn keyword cType BalancevelocityAnalysis -syn keyword cType DamageCalvingAnalysis syn keyword cType DamageEvolutionAnalysis syn keyword cType DepthAverageAnalysis @@ -1405,6 +1406,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GLheightadvectionAnalysis syn keyword cType GiaIvinsAnalysis -syn keyword cType GLheightadvectionAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24059) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24060) @@ -1122,4 +1122,6 @@ NoneEnum, NumberedcostfunctionEnum, + NyeCO2Enum, + NyeH2OEnum, NumericalfluxEnum, OldGradientEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24059) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24060) @@ -1126,4 +1126,6 @@ case NoneEnum : return "None"; case NumberedcostfunctionEnum : return "Numberedcostfunction"; + case NyeCO2Enum : return "NyeCO2"; + case NyeH2OEnum : return "NyeH2O"; case NumericalfluxEnum : return "Numericalflux"; case OldGradientEnum : return "OldGradient"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24059) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24060) @@ -1153,4 +1153,6 @@ else if (strcmp(name,"None")==0) return NoneEnum; else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; + else if (strcmp(name,"NyeCO2")==0) return NyeCO2Enum; + else if (strcmp(name,"NyeH2O")==0) return NyeH2OEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; - else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; - else if (strcmp(name,"Tetra")==0) return TetraEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"TetraInput")==0) return TetraInputEnum; + if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; + else if (strcmp(name,"Tetra")==0) return TetraEnum; + else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 24059) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 24060) @@ -153,5 +153,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'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); + fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'', ''LliboutryDuval'', ''NyeCO2'', or ''NyeH2O''']); case 'litho' disp(sprintf(' \nLitho:')); @@ -183,5 +183,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' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); case 'litho' if ~ismember('LoveAnalysis',analyses), return; end Index: /issm/trunk-jpl/src/m/classes/matice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.m (revision 24059) +++ /issm/trunk-jpl/src/m/classes/matice.m (revision 24060) @@ -118,5 +118,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' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]); @@ -150,5 +150,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'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); + fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'', ''LliboutryDuval'', ''NyeH2O'', or ''NyeCO2''']); fielddisplay(self,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]'); fielddisplay(self,'lithosphere_density','Lithosphere density [g/cm^-3]'); Index: /issm/trunk-jpl/src/m/materials/nye.m =================================================================== --- /issm/trunk-jpl/src/m/materials/nye.m (revision 24060) +++ /issm/trunk-jpl/src/m/materials/nye.m (revision 24060) @@ -0,0 +1,37 @@ +function rigidity = nye(temperature,ice_type) + %NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature + % rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) (Nye, p2000). + % temperature is in Kelvin degrees + % Usage: + % rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice + + if (any(temperature < 130) || any(temperature > 273)) + error('Input temperature out of bounds (130