Index: /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp (revision 24072) +++ /issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp (revision 24073) @@ -1,16 +1,6 @@ /* \file NyeCO2.cpp - * \ brief figure out B of CO2 ice for a certain temperature + * \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 */ @@ -20,105 +10,15 @@ IssmDouble NyeCO2(IssmDouble temperature){ - IssmDouble B,T; + /*Coefficients*/ + const IssmPDouble Rg = 8.3144598; /* J mol^-1 K^-1 */ + const IssmPDouble A_const = pow(10.,10.8); /*s^-1 MPa */ + const IssmPDouble Q = 63000.; /*J mol^-1 */ + const IssmPDouble n = 7.; /*Glen's exponent*/ - T=temperature; + /*Arrhenius Law*/ + IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa */ + IssmDouble B = 1e6*pow(A,-1/n); /*s^(1/n) Pa */ - 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 output*/ return B; } Index: /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp (revision 24072) +++ /issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp (revision 24073) @@ -1,4 +1,4 @@ /* \file NyeH2O.cpp - * \ brief figure out B of H2O ice for a certain temperature + * 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- ). @@ -20,76 +20,15 @@ IssmDouble NyeH2O(IssmDouble temperature){ - IssmDouble B,T; + /*Coefficients*/ + const IssmPDouble Rg = 8.3144598; /* J mol^-1 K^-1 */ + const IssmPDouble A_const = 9.e4; /*s^-1 MPa */ + const IssmPDouble Q = 60000.; /*J mol^-1 */ + const IssmPDouble n = 3.; /*Glen's exponent*/ - /*Switch to celsius from Kelvin: */ - T=temperature-273.15; + /*Arrhenius Law*/ + IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa */ + IssmDouble B = 1e6*pow(A,-1/n); /*s^(1/n) Pa */ - 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 output*/ return B; }