Changeset 24073


Ignore:
Timestamp:
07/09/19 08:28:35 (6 years ago)
Author:
Mathieu Morlighem
Message:

BUG: improved poorly coded law

Location:
issm/trunk-jpl/src/c/shared/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/shared/Elements/NyeCO2.cpp

    r24060 r24073  
    11/* \file NyeCO2.cpp
    2  * \    brief figure out B of CO2 ice for a certain temperature
     2 * \brief figure out B of CO2 ice for a certain temperature
    33 *              INPUT function B=NyeCO2(temperature)
    44 *      where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Nye, p -COMPLETE- ).
    5 
    6  VALUES CALCULATED USING THE FOLLOWING MATLAB CODE
    7 
    8  Rg = 8.3144598; % J mol^-1 K^-1
    9  A_const = 10^(10.8); % s^-1 MPa
    10  Q = 63000; % J mol^-1
    11  n = 7; % Glen's exponent
    12  T = 97.5:5:252.5; % K
    13  A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
    14  B = A.^(-1/n)*1e6; % s^(1/n) Pa
    155 */
    166
    … …  
    2010IssmDouble NyeCO2(IssmDouble temperature){
    2111
    22         IssmDouble B,T;
     12        /*Coefficients*/
     13        const IssmPDouble Rg      = 8.3144598;     /* J mol^-1 K^-1 */
     14        const IssmPDouble A_const = pow(10.,10.8); /*s^-1 MPa       */
     15        const IssmPDouble Q       = 63000.;        /*J mol^-1       */
     16        const IssmPDouble n       = 7.;            /*Glen's exponent*/
    2317
    24         T=temperature;
     18        /*Arrhenius Law*/
     19        IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa   */
     20        IssmDouble B = 1e6*pow(A,-1/n);                    /*s^(1/n) Pa */
    2521
    26         if(T<=100.0){
    27                 B=1899806746.55828;
    28         }
    29         else if((T>=100.0) && (T<=105.0)){
    30                 B=1105378610.10266;
    31         }
    32         else if((T>=105.0) && (T<=110.0)){
    33                 B=676381320.528217;
    34         }
    35         else if((T>=110.0) && (T<=115.0)){
    36                 B=432348426.904519;
    37         }
    38         else if((T>=115.0) && (T<=120.0)){
    39                 B=287089456.910568;
    40         }
    41         else if((T>=120.0) && (T<=125.0)){
    42                 B=197113445.038555;
    43         }
    44         else if((T>=125.0) && (T<=130.0)){
    45                 B=139387289.896988;
    46         }
    47         else if((T>=130.0) && (T<=135.0)){
    48                 B=101178460.700527;
    49         }
    50         else if((T>=135.0) && (T<=140.0)){
    51                 B=75174730.2531849;
    52         }
    53         else if((T>=140.0) && (T<=145.0)){
    54                 B=57030799.0697205;
    55         }
    56         else if((T>=145.0) && (T<=150.0)){
    57                 B=44083907.2876590;
    58         }
    59         else if((T>=150.0) && (T<=155.0)){
    60                 B=34656426.3665295;
    61         }
    62         else if((T>=155.0) && (T<=160.0)){
    63                 B=27664457.7106089;
    64         }
    65         else if((T>=160.0) && (T<=165.0)){
    66                 B=22391479.9130613;
    67         }
    68         else if((T>=165.0) && (T<=170.0)){
    69                 B=18353816.3612526;
    70         }
    71         else if((T>=170.0) && (T<=175.0)){
    72                 B=15218650.2018789;
    73         }
    74         else if((T>=175.0) && (T<=180.0)){
    75                 B=12752901.3547654;
    76         }
    77         else if((T>=180.0) && (T<=185.0)){
    78                 B=10790666.8437103;
    79         }
    80         else if((T>=185.0) && (T<=190.0)){
    81                 B=9212075.16184763;
    82         }
    83         else if((T>=190.0) && (T<=195.0)){
    84                 B=7929303.02671413;
    85         }
    86         else if((T>=195.0) && (T<=200.0)){
    87                 B=6877172.17359016;
    88         }
    89         else if((T>=200.0) && (T<=205.0)){
    90                 B=6006726.68822080;
    91         }
    92         else if((T>=205.0) && (T<=210.0)){
    93                 B=5280781.91812612;
    94         }
    95         else if((T>=210.0) && (T<=215.0)){
    96                 B=4670797.58917498;
    97         }
    98         else if((T>=215.0) && (T<=220.0)){
    99                 B=4154653.09694929;
    100         }
    101         else if((T>=220.0) && (T<=225.0)){
    102                 B=3715045.71619004;
    103         }
    104         else if((T>=225.0) && (T<=230.0)){
    105                 B=3338324.34308717;
    106         }
    107         else if((T>=230.0) && (T<=235.0)){
    108                 B=3013631.36885997;
    109         }
    110         else if((T>=235.0) && (T<=240.0)){
    111                 B=2732264.98697388;
    112         }
    113         else if((T>=240.0) && (T<=245.0)){
    114                 B=2487200.85823651;
    115         }
    116         else if((T>=245.0) && (T<=250.0)){
    117                 B=2272730.12712433;
    118         }
    119         else{
    120                 B=2084183.18884150;
    121         }
    122 
     22        /*Return output*/
    12323        return B;
    12424}

  • issm/trunk-jpl/src/c/shared/Elements/NyeH2O.cpp

    r24060 r24073  
    11/* \file NyeH2O.cpp
    2  * \    brief figure out B of H2O ice for a certain temperature
     2 * brief figure out B of H2O ice for a certain temperature
    33 *              INPUT function B=NyeH2O(temperature)
    44 *      where rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Nye, p -COMPLETE- ).
    … …  
    2020IssmDouble NyeH2O(IssmDouble temperature){
    2121
    22         IssmDouble B,T;
     22        /*Coefficients*/
     23        const IssmPDouble Rg      = 8.3144598; /* J mol^-1 K^-1 */
     24        const IssmPDouble A_const = 9.e4;      /*s^-1 MPa       */
     25        const IssmPDouble Q       = 60000.;    /*J mol^-1       */
     26        const IssmPDouble n       = 3.;        /*Glen's exponent*/
    2327
    24         /*Switch to celsius from Kelvin: */
    25         T=temperature-273.15;
     28        /*Arrhenius Law*/
     29        IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa   */
     30        IssmDouble B = 1e6*pow(A,-1/n);                    /*s^(1/n) Pa */
    2631
    27         if(T<=-100){
    28                 B=29901899089.3257;
    29         }
    30         else if((T>=-100.0) && (T<=-95.0)){
    31                 B=20004935214.6693;
    32         }
    33         else if((T>=-95.0) && (T<=-90.0)){
    34                 B=13685054490.7813;
    35         }
    36         else if((T>=-90.0) && (T<=-85.0)){
    37                 B=9555312498.71763;
    38         }
    39         else if((T>=-85.0) && (T<=-80.0)){
    40                 B=6798800864.99471;
    41         }
    42         else if((T>=-80.0) && (T<=-75.0)){
    43                 B=4922440695.83365;
    44         }
    45         else if((T>=-75.0) && (T<=-70.0)){
    46                 B=3621794191.06010;
    47         }
    48         else if((T>=-70.0) && (T<=-65.0)){
    49                 B=2704902125.98878;
    50         }
    51         else if((T>=-65.0) && (T<=-60.0)){
    52                 B=2048338724.47737;
    53         }
    54         else if((T>=-60.0) && (T<=-55.0)){
    55                 B=1571285974.28101;
    56         }
    57         else if((T>=-55.0) && (T<=-50.0)){
    58                 B=1219918332.04930;
    59         }
    60         else if((T>=-50.0) && (T<=-45.0)){
    61                 B=957813559.547873;
    62         }
    63         else if((T>=-45.0) && (T<=-40.0)){
    64                 B=759956402.731125;
    65         }
    66         else if((T>=-40.0) && (T<=-35.0)){
    67                 B=608924663.484756;
    68         }
    69         else if((T>=-35.0) && (T<=-30.0)){
    70                 B=492424321.822968;
    71         }
    72         else if((T>=-30.0) && (T<=-25.0)){
    73                 B=401672414.952195;
    74         }
    75         else if((T>=-25.0) && (T<=-20.0)){
    76                 B=330320977.865296;
    77         }
    78         else if((T>=-20.0) && (T<=-15.0)){
    79                 B=273731390.267536;
    80         }
    81         else if((T>=-15.0) && (T<=-10.0)){
    82                 B=228478811.983490;
    83         }
    84         else if((T>=-10.0) && (T<=-5.0)){
    85                 B=192009688.386368;
    86         }
    87         else if((T>=-5.0) && (T<=0.0)){
    88                 B=162402355.392225;
    89         }
    90         else{
    91                 B=138197905.740112;
    92         }
    93 
     32        /*Return output*/
    9433        return B;
    9534}
Note: See TracChangeset for help on using the changeset viewer.