Ice Sheet System Model  4.18
Code documentation
Functions
elements.h File Reference

prototypes for elements.h More...

#include "../Numerics/types.h"

Go to the source code of this file.

Functions

IssmDouble Cuffey (IssmDouble temperature)
 
IssmDouble BuddJacka (IssmDouble temperature)
 
IssmDouble CuffeyTemperate (IssmDouble temperature, IssmDouble waterfraction, IssmDouble stressexp)
 
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 EstarLambdaS (IssmDouble epseff, IssmDouble epsprime_norm)
 
void EstarOmega (IssmDouble *omega, IssmDouble vx, IssmDouble vy, IssmDouble vz, IssmDouble vmag, IssmDouble *dvx, IssmDouble *dvy, IssmDouble *dvz, IssmDouble *dvmag)
 
void EstarStrainrateQuantities (IssmDouble *pepsprime_norm, IssmDouble vx, IssmDouble vy, IssmDouble vz, IssmDouble vmag, IssmDouble *dvx, IssmDouble *dvy, IssmDouble *dvz, IssmDouble *dvmag)
 
IssmDouble PddSurfaceMassBalance (IssmDouble *monthlytemperatures, IssmDouble *monthlyprec, IssmDouble *pdds, IssmDouble *pds, IssmDouble *melt, IssmDouble *accu, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, IssmDouble s0t, IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, IssmDouble TdiffTime, IssmDouble sealevTime, IssmDouble pddsnowfac, IssmDouble pddicefac, IssmDouble rho_water, IssmDouble rho_ice)
 
IssmDouble PddSurfaceMassBalanceSicopolis (IssmDouble *monthlytemperatures, IssmDouble *monthlyprec, IssmDouble *melt, IssmDouble *accu, IssmDouble *melt_star, IssmDouble *t_ampl, IssmDouble *p_ampl, IssmDouble yts, IssmDouble s, IssmDouble desfac, IssmDouble s0t, IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice)
 
void ComputeDelta18oTemperaturePrecipitation (IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, IssmDouble *PrecipitationsPresentday, IssmDouble *TemperaturesLgm, IssmDouble *TemperaturesPresentday, IssmDouble *monthlytemperaturesout, IssmDouble *monthlyprecout)
 
void ComputeMungsmTemperaturePrecipitation (IssmDouble TdiffTime, IssmDouble PfacTime, IssmDouble *PrecipitationsLgm, IssmDouble *PrecipitationsPresentday, IssmDouble *TemperaturesLgm, IssmDouble *TemperaturesPresentday, IssmDouble *monthlytemperaturesout, IssmDouble *monthlyprecout)
 
void ComputeD18OTemperaturePrecipitationFromPD (IssmDouble d018, IssmDouble dpermil, bool isTemperatureScaled, bool isPrecipScaled, IssmDouble f, IssmDouble *PrecipitationPresentday, IssmDouble *TemperaturePresentday, IssmDouble *PrecipitationReconstructed, IssmDouble *TemperatureReconstructed, IssmDouble *monthlytemperaturesout, IssmDouble *monthlyprecout)
 
IssmDouble DrainageFunctionWaterfraction (IssmDouble waterfraction, IssmDouble dt=0.)
 
IssmDouble StressIntensityIntegralWeight (IssmDouble depth, IssmDouble water_depth, IssmDouble thickness)
 
void EnthalpyToThermal (IssmDouble *ptemperature, IssmDouble *pwaterfraction, IssmDouble enthalpy, IssmDouble pressure)
 
void printarray (IssmPDouble *array, int lines, int cols=1)
 
void printarray_matlab (const char *filename, int *array, int lines, int cols=1)
 
void printarray (int *array, int lines, int cols=1)
 
void printarray (bool *array, int lines, int cols=1)
 
void printsparsity (IssmPDouble *array, int lines, int cols=1)
 
void printbinary (int n)
 

Detailed Description

prototypes for elements.h

Definition in file elements.h.

Function Documentation

◆ Cuffey()

IssmDouble Cuffey ( IssmDouble  temperature)

Definition at line 11 of file Cuffey.cpp.

11  {
12 
13  /*output: */
14  IssmDouble B,T;
15 
16  /*Switch to celsius from Kelvin: */
17  T=temperature-273.15;
18 
19  if(T<=-45.0){
20  B=1.e+8*(-0.000396645116301*pow(T+50.,3)+ 0.013345579471334*pow(T+50.,2) -0.356868703259105*(T+50.)+7.272363035371383);
21  }
22  else if((T>=-45.0) && (T<=-40.0)){
23  B=1.e+8*(-0.000396645116301*pow(T+45.,3)+ 0.007395902726819*pow(T+45.,2) -0.253161292268336*(T+45.)+5.772078366321591);
24  }
25  else if((T>=-40.0) && (T<=-35.0)){
26  B=1.e+8*(0.000408322072669*pow(T+40.,3)+ 0.001446225982305*pow(T+40.,2) -0.208950648722716*(T+40.)+4.641588833612773);
27  }
28  else if((T>=-35.0) && (T<=-30.0)){
29  B=1.e+8*(-0.000423888728124*pow(T+35.,3)+ 0.007571057072334*pow(T+35.,2) -0.163864233449525*(T+35.)+3.684031498640382);
30  }
31  else if((T>=-30.0) && (T<=-25.0)){
32  B=1.e+8*(0.000147154327025*pow(T+30.,3)+ 0.001212726150476*pow(T+30.,2) -0.119945317335478*(T+30.)+3.001000667185614);
33  }
34  else if((T>=-25.0) && (T<=-20.0)){
35  B=1.e+8*(-0.000193435838672*pow(T+25.,3)+ 0.003420041055847*pow(T+25.,2) -0.096781481303861*(T+25.)+2.449986525148220);
36  }
37  else if((T>=-20.0) && (T<=-15.0)){
38  B=1.e+8*(0.000219771255067*pow(T+20.,3)+ 0.000518503475772*pow(T+20.,2) -0.077088758645767*(T+20.)+2.027400665191131);
39  }
40  else if((T>=-15.0) && (T<=-10.0)){
41  B=1.e+8*(-0.000653438900191*pow(T+15.,3)+ 0.003815072301777*pow(T+15.,2) -0.055420879758021*(T+15.)+1.682390865739973);
42  }
43  else if((T>=-10.0) && (T<=-5.0)){
44  B=1.e+8*(0.000692439419762*pow(T+10.,3) -0.005986511201093 *pow(T+10.,2) -0.066278074254598*(T+10.)+1.418983411970382);
45  }
46  else if((T>=-5.0) && (T<=-2.0)){
47  B=1.e+8*(-0.000132282004110*pow(T+5.,3) +0.004400080095332*pow(T+5.,2) -0.074210229783403*(T+5.)+ 1.024485188140279);
48  }
49  else{
50  B=1.e+8*(-0.000132282004110*pow(T+2.,3) +0.003209542058346*pow(T+2.,2) -0.051381363322371*(T+2.)+ 0.837883605537096);
51  }
52 
53  /*B cannot be negative!*/
54  if(B<0) B=1.e+6;
55 
56  return B;
57 }

◆ BuddJacka()

IssmDouble BuddJacka ( IssmDouble  temperature)

Definition at line 10 of file BuddJacka.cpp.

10  {
11 
12  /*output: */
13  IssmDouble B,T;
14 
15  /*Switch to celsius from Kelvin: */
16  T=temperature-273.15;
17 
18  if(T<=-40.){
19  B=1e9*(-0.000031098521204*pow(T+50.,3)+ 0.002234792114381*pow(T+50.,2)-0.065051516643164*(T+50.)+1.005181071430026);
20  }
21  else if((-40.<T) && (T<=-35.)){
22  B=1e9*(-0.000031098521204*pow(T+40.,3)+ 0.001301836478264*pow(T+40.,2)-0.029685230716715*(T+40.)+0.547046595232583);
23  }
24  else if((-35.<T) && (T<=-30.)){
25  B=1e9*(-0.000038394040864*pow(T+35.,3)+ 0.000835358660205*pow(T+35.,2)-0.018999255024368*(T+35.)+0.427279038455119);
26  }
27  else if((-30.<T) && (T<=-25.)){
28  B=1e9*(-0.000007037062330*pow(T+30.,3)+ 0.000259448047242*pow(T+30.,2)-0.013525221487131*(T+30.)+0.348367474730384);
29  }
30  else if((-25.<T) && (T<=-20.)){
31  B=1e9*( 0.000000905055684*pow(T+25.,3)+ 0.000153892112291*pow(T+25.,2)-0.011458520689465*(T+25.)+0.286347935684521);
32  }
33  else if((-20.<T) && (T<=-15.)){
34  B=1e9*(-0.000002025865930*pow(T+20.,3)+ 0.000167467947546*pow(T+20.,2)-0.009851720390281*(T+20.)+0.233015767004928);
35  }
36  else if((-15.<T) && (T<=-10.)){
37  B=1e9*(-0.000014464671112*pow(T+15.,3)+ 0.000137079958603*pow(T+15.,2)-0.008328980859537*(T+15.)+0.187690630500981);
38  }
39  else if((-10.<T) && (T<=-5.)){
40  B=1e9*(-0.000014230086582*pow(T+10.,3)+-0.000079890108083*pow(T+10.,2)-0.008043031606935*(T+10.)+0.147664641279324);
41  }
42  else if((-5.<T) && (T<=-2.)){
43  B=1e9*( 0.000022694046251*pow(T+5.,3)+-0.000293341406806*pow(T+5.,2)-0.009909189181377*(T+5. )+0.103673469719891);
44  }
45  else if((-2.<T) && (T<=-1.)){
46  B=1e9*( 0.000056280347425*pow(T+2.,3)+-0.000089094990549*pow(T+2.,2)-0.011056498373441*(T+2. )+0.071918568763277);
47  }
48  else if((-1.<T)){
49  B=1e9*( 0.000056280347425*pow(T+1.,3)+ 0.000079746051725*pow(T+1.,2)-0.011065847312265*(T+1. )+0.060829255746712);
50  }
51  else{
52  /*FIXME: just copying previous case for now.... Felicity?*/
53  B=1e9*( 0.000056280347425*pow(T+1.,3)+ 0.000079746051725*pow(T+1.,2)-0.011065847312265*(T+1. )+0.060829255746712);
54  }
55 
56  /*B cannot be negative!*/
57  if(B<0) B=1.e+6;
58 
59  return B;
60 }

◆ CuffeyTemperate()

IssmDouble CuffeyTemperate ( IssmDouble  temperature,
IssmDouble  waterfraction,
IssmDouble  stressexp 
)

Definition at line 11 of file CuffeyTemperate.cpp.

11  {
12 
13  return Cuffey(temperature)*pow(1+181.25*max(0., min(0.01, waterfraction)), -1./stressexp);
14 
15 }

◆ Paterson()

IssmDouble Paterson ( IssmDouble  temperature)

Definition at line 11 of file Paterson.cpp.

11  {
12 
13  IssmDouble B,T;
14 
15  /*Switch to celsius from Kelvin: */
16  T=temperature-273.15;
17 
18  if(T<=-45.0){
19  B=1.e+8*(-0.000292866376675*pow(T+50.,3)+ 0.011672640664130*pow(T+50.,2) -0.325004442485481*(T+50.)+ 6.524779401948101);
20  }
21  else if((T>=-45.0) && (T<=-40.0)){
22  B=1.e+8*(-0.000292866376675*pow(T+45.,3)+ 0.007279645014004*pow(T+45.,2) -0.230243014094813*(T+45.)+ 5.154964909039554);
23  }
24  else if((T>=-40.0) && (T<=-35.0)){
25  B=1.e+8*(0.000072737147457*pow(T+40.,3)+ 0.002886649363879*pow(T+40.,2) -0.179411542205399*(T+40.)+ 4.149132666831214);
26  }
27  else if((T>=-35.0) && (T<=-30.0)){
28  B=1.e+8*(-0.000086144770023*pow(T+35.,3)+ 0.003977706575736*pow(T+35.,2) -0.145089762507325*(T+35.)+ 3.333333333333331);
29  }
30  else if((T>=-30.0) && (T<=-25.0)){
31  B=1.e+8*(-0.000043984685769*pow(T+30.,3)+ 0.002685535025386*pow(T+30.,2) -0.111773554501713*(T+30.)+ 2.696559088937191);
32  }
33  else if((T>=-25.0) && (T<=-20.0)){
34  B=1.e+8*(-0.000029799523463*pow(T+25.,3)+ 0.002025764738854*pow(T+25.,2) -0.088217055680511*(T+25.)+ 2.199331606342181);
35  }
36  else if((T>=-20.0) && (T<=-15.0)){
37  B=1.e+8*(0.000136920904777*pow(T+20.,3)+ 0.001578771886910*pow(T+20.,2) -0.070194372551690*(T+20.)+ 1.805165505978111);
38  }
39  else if((T>=-15.0) && (T<=-10.0)){
40  B=1.e+8*(-0.000899763781026*pow(T+15.,3)+ 0.003632585458564*pow(T+15.,2) -0.044137585824322*(T+15.)+ 1.510778053489523);
41  }
42  else if((T>=-10.0) && (T<=-5.0)){
43  B=1.e+8*(0.001676964325070*pow(T+10.,3)- 0.009863871256831*pow(T+10.,2) -0.075294014815659*(T+10.)+ 1.268434288203714);
44  }
45  else if((T>=-5.0) && (T<=-2.0)){
46  B=1.e+8*(-0.003748937622487*pow(T+5.,3)+0.015290593619213*pow(T+5.,2) -0.048160403003748*(T+5.)+ 0.854987973338348);
47  }
48  else{
49  B=1.e+8*(-0.003748937622488*pow(T+2.,3)-0.018449844983174*pow(T+2.,2) -0.057638157095631*(T+2.)+ 0.746900791092860);
50  }
51 
52  /*B cannot be negative!*/
53  if(B<0) B=1.e+6;
54 
55  return B;
56 }

◆ Arrhenius()

IssmDouble Arrhenius ( IssmDouble  temperature,
IssmDouble  depth,
IssmDouble  n 
)

Definition at line 9 of file Arrhenius.cpp.

9  {
10  /*Use EISMINT Parameterization for the rheology: Payne2000
11  *
12  * A(T*) = A0 exp(-Q/RT*)
13  *
14  * A0 constant of proportionality
15  * = 3.61 * 10^-13 if T*<263.15K
16  * = 1.73 * 10^3 if T*>263.15K
17  * Q Activation energy for creep
18  * = 6.0 * 10^4 if T*<263.15K
19  * = 13.9 * 10^4 if T*>263.15K
20  * R Universal gas constant
21  * = 8.314
22  * T* Absolute temperature corrected for the dependence of Tpmp on P
23  * = T - beta (s-z)
24  *
25  * Convert A to B : B = A^(-1/n) */
26 
27  /*Some physical constants (Payne2000)*/
28  IssmDouble beta=8.66*1.e-4;
29  IssmDouble R=8.314;
30 
31  /*Intermediaries*/
32  IssmDouble A,B,Tstar;
33 
34  /*convert temperature to absolute temperature*/
35  _assert_(depth>0);
36  Tstar=temperature-beta*depth;
37  _assert_(Tstar>0);
38 
39  /*Get A*/
40  if(Tstar<263.15){
41  A=3.61e-13*exp( -6.e+4/(R*Tstar));
42  }
43  else{
44  A=1.73e+3 *exp(-13.9e+4/(R*Tstar));
45  }
46 
47  /*Convert to B*/
48  B=pow(A,-1./n);
49 
50  return B;
51 }

◆ NyeH2O()

IssmDouble NyeH2O ( IssmDouble  temperature)

Definition at line 11 of file NyeH2O.cpp.

11  {
12 
13  /*Coefficients*/
14  const IssmPDouble Rg = 8.3144598; /* J mol^-1 K^-1 */
15  const IssmPDouble A_const = 9.e4; /*s^-1 MPa */
16  const IssmPDouble Q = 60000.; /*J mol^-1 */
17  const IssmPDouble n = 3.; /*Glen's exponent*/
18 
19  /*Arrhenius Law*/
20  IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /*s^-1 MPa */
21  IssmDouble B = 1e6*pow(A,-1/n); /*s^(1/n) Pa */
22 
23  /*Beyond-melting-point case*/
24  if(temperature>=273.15) _printf0_("H2O ICE - GUARANTEED MELTING. Some temperature values are beyond 273.15K.\n");
25 
26  /*Return output*/
27  return B;
28 }

◆ NyeCO2()

IssmDouble NyeCO2 ( IssmDouble  temperature)

Definition at line 11 of file NyeCO2.cpp.

11  {
12 
13  /*Coefficients*/
14  const IssmPDouble Rg = 8.3144598; /* J mol^-1 K^-1 */
15  const IssmPDouble A_const = pow(10.,13.0); /* s^-1 MPa */
16  const IssmPDouble Q = 66900.; /* J mol^-1 */
17  const IssmPDouble n = 8.; /* Glen's exponent */
18 
19  /*Arrhenius Law*/
20  IssmDouble A = A_const *exp(-Q/(temperature*Rg)); /* s^-1 MPa */
21  IssmDouble B = 1e6*pow(A,-1/n); /* s^(1/n) Pa */
22 
23  /*Beyond-melting-point cases*/
24  if((temperature>200.)&&(temperature<220.)) _printf0_("CO2 ICE - POSSIBLE MELTING. Some temperature values are between 200K and 220K.\n");
25  else if(temperature>=220.) _printf0_("CO2 ICE - GUARANTEED MELTING. Some temperature values are beyond 220K.\n");
26 
27  /*Return output*/
28  return B;
29 }

◆ LliboutryDuval()

IssmDouble LliboutryDuval ( IssmDouble  enthalpy,
IssmDouble  pressure,
IssmDouble  n,
IssmDouble  betaCC,
IssmDouble  referencetemperature,
IssmDouble  heatcapacity,
IssmDouble  latentheat 
)

Definition at line 10 of file LliboutryDuval.cpp.

10  {
11  /*Use Lliboutry & Duval's 1985 parameterization for the rheology:
12  * see also: Grewe/Blatter 2009, Aschwanden et al. 2012
13  *
14  * ISSM uses enthalpy/temperature values that are not corrected for pressure.
15  *
16  * A(H,p) = A0 exp(-Q/RT(H,p)), if H < H_s(p)
17  * = A0 exp(-Q/RTpmp) (1+181.25w(H,p)), if H_s(p) \le H < H_l(p)
18  *
19  * T(H,p) = Tref + H/c_i, if H < H_s(p)
20  * = Tpmp , if H_s(p) \le H \le H_l(p)
21  *
22  * w(H,p) = 0, if H < H_s(p)
23  * = (H - H_s(p))/L
24  *
25  * H_s(p) = c_i (Tpmp - Tref)
26  *
27  * Tpmp = T - betaCC p;
28  *
29  * A0 constant of proportionality
30  * = 3.61 * 10^-13 if T*<263.15K
31  * = 1.73 * 10^3 if T*>263.15K
32  * Q Activation energy for creep
33  * = 6.0 * 10^4 if T*<263.15K
34  * = 13.9 * 10^4 if T*>263.15K
35  * R Universal gas constant
36  * = 8.314
37  *
38  * Convert A to B : B = A^(-1/n) */
39  /*check feasibility*/
40  _assert_(pressure+1.e-4>=0); // deal with pressure instability at ice surface
41  _assert_(n>0);
42  _assert_(betaCC>=0);
43  _assert_(referencetemperature>=0);
44  _assert_(heatcapacity>0);
45  _assert_(latentheat>0);
46 
47  /*Some physical constants*/
48  IssmDouble R=8.314;
49 
50  /*Intermediaries*/
51  IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction;
52 
53  Tpmp=273.15-betaCC*pressure; //pressure melting point temperature
54  H_sp=heatcapacity*(Tpmp-referencetemperature); //pressure melting point enthalpy
55  if (enthalpy<H_sp){ //compute homologous temperature and water fraction
56  Tstar=referencetemperature+enthalpy/heatcapacity+betaCC*pressure;
57  waterfraction=0.;
58  }
59  else{
60  Tstar=273.15;
61  waterfraction=(enthalpy-H_sp)/latentheat;
62  if (waterfraction>0.01) waterfraction=0.01; // limit softness of ice
63  }
64 
65  /*Get A*/
66  if(Tstar<=263.15){A=3.61e-13*exp(-6.e+4/(R*Tstar));}
67  else{A=1.73e3*exp(-13.9e+4/(R*Tstar));}
68  A*=(1.+181.25*waterfraction);
69 
70  /*Convert to B*/
71  B=pow(A,-1./n);
72  return B;
73 }

◆ EstarLambdaS()

IssmDouble EstarLambdaS ( IssmDouble  epseff,
IssmDouble  epsprime_norm 
)

Definition at line 118 of file EstarComponents.cpp.

118  {/*{{{*/
119  IssmDouble lambdas;
120 
121  _assert_(epsprime_norm>=0.);
122  if(epseff==0.){
123  lambdas=0.;
124  }
125  else{
126  lambdas=sqrt(epsprime_norm*epsprime_norm/(epseff*epseff));
127  }
128  return lambdas;
129 }/*}}}*/

◆ EstarOmega()

void EstarOmega ( IssmDouble *  omega,
IssmDouble  vx,
IssmDouble  vy,
IssmDouble  vz,
IssmDouble  vmag,
IssmDouble *  dvx,
IssmDouble *  dvy,
IssmDouble *  dvz,
IssmDouble *  dvmag 
)

Definition at line 72 of file EstarComponents.cpp.

72  {/*{{{*/
73 
74  /*Intermediaries*/
75  IssmDouble omega_norm;
76  IssmDouble omega_rigid[3];
77 
78  /*Create vorticity vector*/
79  _assert_(dvx && dvy && dvz && dvmag);
80  if(vmag<1e-12)vmag=1e-12;
81 
82  /*Create vorticity vector, corrected for rigid body rotation
83  * \overline{\omega} =\omega - \omega_rigid
84  * =\nabla\times{\bf v} - 2*U*\kappa*\hat{b};
85  * =\nabla\times{\bf v} - (2*{\bf v}\times(({\bf v}\cdot\nabla)*{\bf v}))/U^2
86  * check the magnitude of the second term -- if it is small, then the two
87  * vorticities (omega and first term in omega) are approx. equal
88  *
89  * */
90  omega_rigid[0] = 2*(vy*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]) - vz*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]))/(vmag*vmag);
91  omega_rigid[1] = 2*(vz*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]) - vx*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]))/(vmag*vmag);
92  omega_rigid[2] = 2*(vx*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]) - vy*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]))/(vmag*vmag);
93 
94  omega[0] = (dvz[1] - dvy[2]) - omega_rigid[0];
95  omega[1] = (dvx[2] - dvz[0]) - omega_rigid[1];
96  omega[2] = (dvy[0] - dvx[1]) - omega_rigid[2];
97 
98  /*Take out vorticity component aligned with the velocity*/
99  IssmDouble wdotv = vx/vmag*omega[0] + vy/vmag*omega[1] + vz/vmag*omega[2];
100  omega[0] = omega[0] - wdotv*vx/vmag;
101  omega[1] = omega[1] - wdotv*vy/vmag;
102  omega[2] = omega[2] - wdotv*vz/vmag;
103 
104  /*Normalize*/
105  omega_norm = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
106  if(omega_norm==0){
107  omega[0] = 0.;
108  omega[1] = 0.;
109  omega[2] = 0.;
110  }
111  else{
112  omega[0] =omega[0]/omega_norm;
113  omega[1] =omega[1]/omega_norm;
114  omega[2] =omega[2]/omega_norm;
115  }
116 
117 }/*}}}*/

◆ EstarStrainrateQuantities()

void EstarStrainrateQuantities ( IssmDouble *  pepsprime_norm,
IssmDouble  vx,
IssmDouble  vy,
IssmDouble  vz,
IssmDouble  vmag,
IssmDouble *  dvx,
IssmDouble *  dvy,
IssmDouble *  dvz,
IssmDouble *  dvmag 
)

Definition at line 5 of file EstarComponents.cpp.

5  {/*{{{*/
6 
7  /*Intermediaries*/
8  IssmDouble omega[3];
9  IssmDouble nrsp[3],nrsp_norm;
10  IssmDouble eps[3][3];
11  IssmDouble epsprime[3],epsprime_norm;
12 
13  /*Get omega, correction for rigid body rotation*/
14  EstarOmega(&omega[0],vx,vy,vz,vmag,dvx,dvy,dvz,dvmag);
15 
16  /*Non-rotating shear plane*/
17  nrsp[0] = vy*omega[2] - vz*omega[1];
18  nrsp[1] = vz*omega[0] - vx*omega[2];
19  nrsp[2] = vx*omega[1] - vy*omega[0];
20 
21  /*Normalize*/
22  nrsp_norm = sqrt(nrsp[0]*nrsp[0] + nrsp[1]*nrsp[1] + nrsp[2]*nrsp[2]);
23  if(nrsp_norm==0){
24  nrsp[0] = 0.;
25  nrsp[1] = 0.;
26  nrsp[2] = 0.;
27  }
28  else{
29  nrsp[0] =nrsp[0]/nrsp_norm;
30  nrsp[1] =nrsp[1]/nrsp_norm;
31  nrsp[2] =nrsp[2]/nrsp_norm;
32  }
33 
34  /*Build strain rate tensor*/
35  eps[0][0] = dvx[0]; eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]);
36  eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]);
37  eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
38 
39  /*Compute the shear strain rate on the non rotating shear plane*/
40  epsprime[0]=0.;
41  epsprime[1]=0.;
42  epsprime[2]=0.;
43  /*term #1: eps'.n */
44  for(int i=0;i<3;i++){
45  for(int j=0;j<3;j++){
46  epsprime[i] += eps[i][j]*nrsp[j];
47  }
48  }
49  /*term #2: ((eps'.n).n)n */
50  for(int i=0;i<3;i++){
51  for(int j=0;j<3;j++){
52  for(int k=0;k<3;k++){
53  epsprime[j] += -nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
54  }
55  }
56  }
57  /*term #3: ((eps'.n).omega)omega */
58  for(int i=0;i<3;i++){
59  for(int j=0;j<3;j++){
60  for(int k=0;k<3;k++){
61  epsprime[j] += -nrsp[i]*eps[i][k]*omega[k]*omega[j];
62  }
63  }
64  }
65 
66  /*Get norm of epsprime*/
67  epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
68 
69  /*Assign output pointers*/
70  *pepsprime_norm=epsprime_norm;
71 }/*}}}*/

◆ PddSurfaceMassBalance()

IssmDouble PddSurfaceMassBalance ( IssmDouble *  monthlytemperatures,
IssmDouble *  monthlyprec,
IssmDouble *  pdds,
IssmDouble *  pds,
IssmDouble *  melt,
IssmDouble *  accu,
IssmDouble  signorm,
IssmDouble  yts,
IssmDouble  h,
IssmDouble  s,
IssmDouble  desfac,
IssmDouble  s0t,
IssmDouble  s0p,
IssmDouble  rlaps,
IssmDouble  rlapslgm,
IssmDouble  TdiffTime,
IssmDouble  sealevTime,
IssmDouble  pddsnowfac,
IssmDouble  pddicefac,
IssmDouble  rho_water,
IssmDouble  rho_ice 
)

Definition at line 11 of file PddSurfaceMassBalance.cpp.

16  {
17 
18  // output:
19  IssmDouble B; // surface mass balance, melt+accumulation
20  int iqj,imonth;
21 
22  IssmDouble saccu; // yearly surface accumulation
23  IssmDouble smelt; // yearly melt
24  IssmDouble precrunoff; // yearly runoff
25  IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
26  IssmDouble water; //water=rain + snowmelt
27  IssmDouble runoff; //meltwater only, does not include rain
28  IssmDouble sconv; //rhow_rain/rhoi / 12 months
29 
30  //IssmDouble sealev=0.; // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
31  //IssmDouble Pfac=0.5,Tdiff=0.5;
32  IssmDouble rtlaps;
33  // IssmDouble lapser=6.5 // lapse rate
34  // IssmDouble desfac = 0.3; // desert elevation factor
35  // IssmDouble s0p=0.; // should be set to elevation from precip source
36  // IssmDouble s0t=0.; // should be set to elevation from temperature source
37  IssmDouble st; // elevation between altitude of the temp record and current altitude
38  IssmDouble sp; // elevation between altitude of the prec record and current altitude
39  IssmDouble deselcut=1.0;
40 
41  // PDD and PD constants and variables
42  IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
43  IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
44  IssmDouble siglimc, siglim0, siglim0c;
45  IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
46  IssmDouble DT = 0.02;
47  IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
48 
49  IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
50  IssmDouble qm = 0.; // snow part of the precipitation
51  IssmDouble qmt = 0.; // precipitation without desertification effect adjustment
52  IssmDouble qmp = 0.; // desertification taken into account
53  IssmDouble pdd = 0.;
54  IssmDouble frzndd = 0.;
55 
56  IssmDouble tstar; // monthly mean surface temp
57  IssmDouble Tsum= 0.; // average summer (JJA) temperature
58  IssmDouble Tsurf = 0.; // average annual temperature
59 
60  IssmDouble deltm=1./12.;
61  int ismon[12]={11,0,1,2,3,4,5,6,7,8,9,10};
62 
63  IssmDouble snwm; // snow that could have been melted in a year.
64  IssmDouble snwmf; // ablation factor for snow per positive degree day.
65  IssmDouble smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
66 
67  IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
68  IssmDouble supice,supcap,diffndd;
69  IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
70  IssmDouble pddtj, hmx2;
71  IssmDouble pddsnowfac0=4.3, pddicefac0=8.3;
72  IssmDouble snowfac, icefac;
73 
74  sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
75 
76  /*PDD constant*/
77  siglim = 2.5*signorm;
78  siglimc = 2.5*signormc;
79  siglim0 = siglim/DT + 0.5;
80  siglim0c = siglimc/DT + 0.5;
81  PDup = siglimc+PDCUT;
82 
83  // seasonal loop
84  for (iqj = 0; iqj < 12; iqj++){
85  imonth = ismon[iqj];
86 
87  /*********compute lapse rate ****************/
88  st=(s-s0t)/1000.;
89  rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
90 
91  /*********compute Surface temperature *******/
92  monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
93  tstar = monthlytemperatures[imonth];
94  Tsurf = tstar*deltm+Tsurf;
95 
96  /*********compute PD ****************/
97  if (tstar < PDup){
98  pd = 1.;
99  if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
100  else {
101  pd = 0.;}
102 
103  /******exp des/elev precip reduction*******/
104  sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
105  if (sp>0.0){q = exp(-desfac*sp);}
106  else {q = 1.0;}
107 
108  qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month
109  qmpt= q*monthlyprec[imonth]*sconv;
110  qmp= qmp + qmpt;
111  qm= qm + qmpt*pd;
112 
113  /*********compute PDD************/
114  // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
115  // gaussian=T_m, so ndd=-(Tsurf-pdd)
116  if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;}
117 
118  if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
119  else if (tstar> -siglim){
120  pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
121  pdd = pdd + pddsig*deltm;
122  frzndd = frzndd - (tstar-pddsig)*deltm;}
123  else{frzndd = frzndd - tstar*deltm; }
124 
125  /*Assign output pointer*/
126  *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
127  *(monthlyprec+imonth) = monthlyprec[imonth];
128  } // end of seasonal loop
129  //******************************************************************
130 
131  saccu = qm;
132  prect = qmp; // total precipitation during 1 year taking into account des. ef.
133  Tsum=Tsum/3;
134 
135  snowfac=pddsnowfac0;
136  icefac=pddicefac0;
137  if (pddsnowfac>0) {
138  if (pddsnowfac<1.65) {
139  _printf0_("WARNING: Pdd snow factor input, " << pddsnowfac << ", results in a negative value. It will be ignored. \n");
140  }
141  else{
142  snowfac=pddsnowfac;
143  }
144  }
145  if (pddicefac>0) {
146  if (pddicefac>17.22) {
147  _printf0_("WARNING: Pdd ice factor input, " << pddicefac << ", results in a negative value. It will be ignored. \n");
148  }
149  else{
150  icefac=pddicefac;
151  }
152  }
153 
154  /***** determine PDD factors *****/
155  if(Tsum< -1.) {
156  snwmf=(2.65+snowfac-pddsnowfac0)*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
157  smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
158  }
159  else if(Tsum< 10){
160  snwmf = (0.15*(Tsum+1) + (2.65+snowfac-pddsnowfac0))*0.001;
161  smf = (((17.22-icefac)/(pow(11,3)))*pow((10.-Tsum),3) + pddicefac0)*0.001;
162  //JC,smf = (((icefac-pddicefac0)/(Tsum+1))*pow((10.-Tsum),3) + pddicefac0)*0.001;
163  }
164  else{
165  snwmf=snowfac*0.001;
166  smf=icefac*0.001;
167  }
168  snwmf=0.95*snwmf;
169  smf=0.95*smf;
170 
171  /***** compute PDD ablation and refreezing *****/
172  pddt = pdd *365;
173  snwm = snwmf*pddt; // snow that could have been melted in a year
174  hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz
175 
176  if(snwm < saccu) {
177  water=prect-saccu + snwm; //water=rain + snowmelt
178  // l 2.2= capillary factor
179  // Should refreezing be controlled by frzndd or by mean annual Tsurf?
180  // dCovrLm concept is of warming of active layer (thickness =d@=1-
181  // >2m)
182  // problem with water seepage into ice: should be sealed after
183  // refreezing
184  // so everything needs to be predicated on 1 year scale, except for
185  // thermal
186  // conductivity through ice
187  // also, need to account that melt season has low accum, so what's
188  // going to
189  // hold the meltwater around for refreezing? And melt-time will have
190  // low seasonal frzndd
191 
192  // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002
193 
194  supice= min(hmx2*CovrLm*frzndd+2.2*(saccu-snwm), water); // superimposed ice
195  supcap=min(2.2*(saccu-snwm),water);
196  runoff=snwm - supice; //meltwater only, does not include rain
197  }
198  else { //all snow melted
199  supice= min(hmx2*CovrLm*frzndd, prect );
200  runoff= saccu + smf*(pddt-saccu/snwmf) - supice;
201  supcap=0;
202  }
203  // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
204  // except pdd melt heat source is atmosphere, while refreeze is
205  // ground/ice stored interim
206  // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
207  // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
208  // <0
209  // assume ndd>pdd, little melt => little supice
210  // bottom line: compare for Tsurf<0 : supice and no supice case,
211  // expect Tsurf difference
212  // except some of cooling flux comes from atmosphere//
213  // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
214  // does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
215  // < 0.1
216 
217  // make more sense to just use residual pdd-ndd except that pdd
218  // residual not clear yet
219  // frzndd should not be used up by refreezing in snow, so stick in
220  // supcap.
221  diffndd=0;
222  if (frzndd>0) {
223  diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd);
224  frzndd=frzndd-diffndd;
225  }
226  if(runoff<0){
227  saccu= saccu -runoff;
228  smelt = 0;
229  precrunoff=prect-saccu;
230  //here assume pdd residual is 0, =>
231  Tsurf= max(Tsurf,-frzndd);
232  }
233  else {
234  smelt = runoff;
235  precrunoff=prect-max(0.,supice)-saccu;}
236  //here really need pdd balance, try 0.5 fudge factor?
237  //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
238  //yet from site plots, can be ice free with Tsurf=-5.5C
239  if(Tsurf<0) {
240  Tsurf= min(Tsurf+fsupT*diffndd , 0.);}
241 
242  melt[0]=smelt/yts;
243  accu[0]=saccu/yts;
244  B = saccu - smelt;
245  B = B/yts;
246  pddtj=pddt;
247 
248  return B;
249 }

◆ PddSurfaceMassBalanceSicopolis()

IssmDouble PddSurfaceMassBalanceSicopolis ( IssmDouble *  monthlytemperatures,
IssmDouble *  monthlyprec,
IssmDouble *  melt,
IssmDouble *  accu,
IssmDouble *  melt_star,
IssmDouble *  t_ampl,
IssmDouble *  p_ampl,
IssmDouble  yts,
IssmDouble  s,
IssmDouble  desfac,
IssmDouble  s0t,
IssmDouble  s0p,
IssmDouble  rlaps,
IssmDouble  rho_water,
IssmDouble  rho_ice 
)

Definition at line 10 of file PddSurfaceMassBalanceSicopolis.cpp.

14  {
15 
16  int imonth; // month counter
17  IssmDouble B; // output: surface mass balance (m/a IE), melt+accumulation
18  IssmDouble frac_solid, snowfall, rainfall, runoff;
19  IssmDouble saccu; // yearly surface accumulation (m/a IE)
20  IssmDouble smelt; // yearly melt (m/a IE)
21  IssmDouble smelt_star; // yearly ...
22  IssmDouble precip; // total precipitation during 1 year
23  IssmDouble sconv; //rhow_rain/rhoi / 12 months
24  IssmDouble st; // elevation between altitude of the temp record and current altitude
25  IssmDouble sp; // elevation between altitude of the prec record and current altitude
26  IssmDouble q; // q is desert/elev. fact
27  IssmDouble pdd; // pdd factor (a * degC)
28  IssmDouble tstar; // monthly temp. after lapse rate correction (degC)
29  IssmDouble precip_star; // monthly precip after correction (m/a IE)
30  IssmDouble beta1 = 2.73; // 3 mm IE/(d*deg C), ablation factor for snow per positive degree day.
31  IssmDouble beta2 = 7.28; // 8 mm IE/(d*deg C), ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
32  IssmDouble Pmax = 0.6;
33  IssmDouble inv_twelve=1./12.;
34 
35  sconv=(rho_water/rho_ice); //rhow_rain/rhoi
36 
37  /* FIXME betas shoud be user input */
38  beta1=beta1*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
39  beta2=beta2*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
40 
41  /* initalize fields */
42  precip=0.0;
43  tstar=0.0;
44  snowfall=0.0;
45  pdd=0.0;
46  /* seasonal loop */
47  for(imonth=0;imonth<12;imonth++){
48 
49  /********* Surface temperature correction *******/
50  st=(s-s0t)/1000.;
51 
52  // FIXME rlaps ??
53  rlaps=-6.309e-03+(-5.426e-03-(-6.309e-03))*sin((imonth+1-4)*PI/6.0)*1000.0;
54  monthlytemperatures[imonth]=monthlytemperatures[imonth]-rlaps*st;//*max(st,1e-3);
55  tstar=monthlytemperatures[imonth]+t_ampl[0];
56 
57  /********* Precipitation correction *************/
58  /* Ref: Vizcaino et al 2010; DOI 10.1007/s00382-009-0591-y */
59  if(s0p<2000.0)
60  q=exp(desfac*(max(s,2000.0)-2000.0));
61  else
62  q=exp(desfac*(max(s,2000.0)-s0p));
63 
64  precip_star=q*monthlyprec[imonth]*sconv*p_ampl[0]*yts; // convert precip from m/s -> m/a
65  precip=precip+precip_star*inv_twelve;
66 
67  /********* compute PDD **************************/
68  /* Ref: Calov & Greve 2005 Journal of Glaciology, Vol. 51, No. 172, 2005, Correspondence */
69  IssmDouble s_stat=5.0;
70  IssmDouble inv_sqrt2pi =1.0/sqrt(2.0*PI);
71  IssmDouble inv_s_stat =1.0/s_stat;
72  IssmDouble inv_sqrt2 =1.0/sqrt(2.0);
73 
74  #if !defined(_HAVE_ADOLC_)
75  pdd=pdd+(s_stat*inv_sqrt2pi*exp(-0.5*pow(tstar*inv_s_stat,2))
76  +0.5*tstar*erfc(-tstar*inv_s_stat*inv_sqrt2))*inv_twelve;
77  #else
78  _error_("Cannot differentiate erfc, talk to ADOLC folks (http://functions.wolfram.com/GammaBetaErf/Erfc/20/01/)");
79  #endif
80 
81  /*Partition of precip in solid and liquid parts, Bales et al. (2009) */
82  IssmDouble temp_rain=7.2; // Threshold monthly mean temperature for
83  // precipitation = 101% rain, in deg C
84  IssmDouble temp_snow=-11.6; // Threshold monthly mean temperature for
85  // precipitation = 100% snow, in deg C
86 
87  IssmDouble coeff1=5.4714e-01; // Coefficients
88  IssmDouble coeff2=-9.1603e-02; // of
89  IssmDouble coeff3=-3.314e-03; // the
90  IssmDouble coeff4= 4.66e-04; // fifth-order
91  IssmDouble coeff5=3.8e-05; // polynomial
92  IssmDouble coeff6=6.0e-07; // fit
93 
94  if(tstar>=temp_rain)
95  frac_solid = 0.0;
96  else if(tstar<=temp_snow)
97  frac_solid = 1.0;
98  else{
99  frac_solid=coeff1+tstar*(coeff2
100  +tstar*(coeff3+tstar*(coeff4+tstar*(coeff5+tstar*coeff6))));
101  }
102 
103  snowfall=snowfall+precip_star*frac_solid*inv_twelve;
104  }
105  /* end of seasonal loop */
106 
107  rainfall=precip-snowfall;
108  if(snowfall<0.0) snowfall=0.0; // correction of
109  if(rainfall<0.0) rainfall=0.0; // negative values
110 
111  if(rainfall<=(Pmax*snowfall)){
112  if((rainfall+beta1*pdd)<=(Pmax*snowfall)) {
113  smelt_star = rainfall+beta1*pdd;
114  smelt = 0.0;
115  runoff = smelt;
116  }
117  else{
118  smelt_star = Pmax*snowfall;
119  smelt = beta2*(pdd-(smelt_star-rainfall)/beta1);
120  runoff = smelt;
121  }
122  }
123  else{
124  smelt_star = Pmax*snowfall;
125  smelt = beta2*pdd;
126  runoff = smelt+rainfall-Pmax*snowfall;
127  }
128 
129  saccu = precip;
130 
131  /* asign output*/
132  melt[0]=runoff/yts;
133  accu[0]=saccu/yts;
134  melt_star[0]=smelt_star/yts;
135  B=(saccu-runoff)/yts;
136 
137  return B;
138 }

◆ ComputeDelta18oTemperaturePrecipitation()

void ComputeDelta18oTemperaturePrecipitation ( IssmDouble  Delta18oSurfacePresent,
IssmDouble  Delta18oSurfaceLgm,
IssmDouble  Delta18oSurfaceTime,
IssmDouble  Delta18oPresent,
IssmDouble  Delta18oLgm,
IssmDouble  Delta18oTime,
IssmDouble *  PrecipitationsPresentday,
IssmDouble *  TemperaturesLgm,
IssmDouble *  TemperaturesPresentday,
IssmDouble *  monthlytemperaturesout,
IssmDouble *  monthlyprecout 
)

Definition at line 10 of file ComputeDelta18oTemperaturePrecipitation.cpp.

14  {
15 
16  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
17  IssmDouble delta18oLapseRate=-6.2*pow(10.,-3);
18  IssmDouble glacialindex; // used to vary present day temperature
19 
20  glacialindex = (Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))
21  /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent)); // Tarasov 2004 paper
22 
23  for (int imonth = 0; imonth<12; imonth++){
24  monthlytemperaturestmp[imonth] = glacialindex*TemperaturesLgm[imonth] + (1.-glacialindex)*TemperaturesPresentday[imonth];
25  monthlyprectmp[imonth] = PrecipitationsPresentday[imonth];
26 
27  /*Assign output pointer*/
28  *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
29  *(monthlyprecout+imonth) = monthlyprectmp[imonth];
30  }
31 }

◆ ComputeMungsmTemperaturePrecipitation()

void ComputeMungsmTemperaturePrecipitation ( IssmDouble  TdiffTime,
IssmDouble  PfacTime,
IssmDouble *  PrecipitationsLgm,
IssmDouble *  PrecipitationsPresentday,
IssmDouble *  TemperaturesLgm,
IssmDouble *  TemperaturesPresentday,
IssmDouble *  monthlytemperaturesout,
IssmDouble *  monthlyprecout 
)

Definition at line 11 of file ComputeMungsmTemperaturePrecipitation.cpp.

14  {
15 
16  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
17  IssmDouble tdiffh;
18 
19  for (int imonth = 0; imonth<12; imonth++){
20  tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
21  monthlytemperaturestmp[imonth] = tdiffh + TemperaturesPresentday[imonth] ;
22 
23  monthlyprectmp[imonth] =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime)); // [m/yr]
24 
25  /*Assign output pointer*/
26  *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
27  *(monthlyprecout+imonth) = monthlyprectmp[imonth];
28  }
29  // printf(" tempera %f\n",monthlytemperaturestmp[1]);
30 }

◆ ComputeD18OTemperaturePrecipitationFromPD()

void ComputeD18OTemperaturePrecipitationFromPD ( IssmDouble  d018,
IssmDouble  dpermil,
bool  isTemperatureScaled,
bool  isPrecipScaled,
IssmDouble  f,
IssmDouble *  PrecipitationPresentday,
IssmDouble *  TemperaturePresentday,
IssmDouble *  PrecipitationReconstructed,
IssmDouble *  TemperatureReconstructed,
IssmDouble *  monthlytemperaturesout,
IssmDouble *  monthlyprecout 
)

Definition at line 10 of file ComputeD18OTemperaturePrecipitationFromPD.cpp.

13  {
14 
15  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
16  IssmDouble deltaTemp;
17 
18  /* Constants */
19  // dpermil = 2.4;/*degrees C per mil*/
20 
21  /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/
22  deltaTemp = dpermil * (d018+34.83);
23 
24  for(int imonth = 0; imonth<12; imonth++){
25 
26  if(isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
27  else{
28  monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth];
29  deltaTemp=TemperatureReconstructed[imonth]-TemperaturePresentday[imonth];
30  }
31 
32  if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
33  else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
34 
35  /*Assign output pointer*/
36  *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
37  *(monthlyprecout+imonth) = monthlyprectmp[imonth];
38  }
39 }

◆ DrainageFunctionWaterfraction()

IssmDouble DrainageFunctionWaterfraction ( IssmDouble  waterfraction,
IssmDouble  dt = 0. 
)

Definition at line 9 of file DrainageFunctionWaterfraction.cpp.

9  {
10  /* DrainageFunctionWaterfraction returns how much of the waterfraction is drained per year */
11  _assert_(waterfraction>=0.);
12  _assert_(dt>=0.);
13 
14  IssmDouble w0=0.01, w1=0.02, w2=0.03;
15  IssmDouble yts=365.*24.*60.*60.;
16  IssmDouble Dret, D0=0., D1=0.005/yts, D2=0.05/yts;
17 
18  /*get drainage function value*/
19  if((w0==w1)||(w1==w2)||(w0==w2))
20  _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
21 
22  if(waterfraction<=w0)
23  Dret=D0;
24  else if((waterfraction>w0) && (waterfraction<=w1))
25  Dret=(D1-D0)/(w1-w0)*(waterfraction-w0)+D0;
26  else if((waterfraction>w1) && (waterfraction<=w2))
27  Dret=(D2-D1)/(w2-w1)*(waterfraction-w1)+D1;
28  else
29  Dret=D2;
30 
31  /*drain only up to w0*/
32  if(dt==0.){
33  if(waterfraction>w0)
34  return waterfraction-w0;
35  else
36  return Dret;
37  }
38  else{
39  if((waterfraction>w0) && (waterfraction-dt*Dret<w0))
40  return (waterfraction-w0)/dt;
41  else
42  return Dret;
43  }
44 }

◆ StressIntensityIntegralWeight()

IssmDouble StressIntensityIntegralWeight ( IssmDouble  depth,
IssmDouble  water_depth,
IssmDouble  thickness 
)

Definition at line 9 of file StressIntensityIntegralWeight.cpp.

9  {
10 
11  /*output: */
12  IssmDouble beta;
13 
14  /*intermediaries: */
15  IssmDouble M1,M2,M3,x,y,d;
16  const double pi = 3.141592653589793;
17  x = water_depth/thickness;
18  y = depth;
19  d = water_depth;
20 
21  M1 = 0.0719768-1.513476*x-61.1001*pow(x,2)+1554.95*pow(x,3)-14583.8*pow(x,4)+71590.7*pow(x,5)-205384*pow(x,6)+356469*pow(x,7)-368270*pow(x,8)+208233*pow(x,9)-49544*pow(x,10);
22  //printf("M1 : %g",M1);
23  M2 = 0.246984+6.47583*x+176.456*pow(x,2)-4058.76*pow(x,3)+37303.8*pow(x,4)-181755*pow(x,5)+520551*pow(x,6)-904370*pow(x,7)+936863*pow(x,8)-531940*pow(x,9)+127291*pow(x,10);
24  //printf("M2 : %g",M2);
25  M3 = 0.529659-22.3235*x+532.074*pow(x,2)-5479.53*pow(x,3)+28592.2*pow(x,4)-81388.6*pow(x,5)+128746*pow(x,6)-106246*pow(x,7)+35780.7*pow(x,8);
26  //printf("M3 : %g",M3);
27 
28  beta = 2/sqrt(2*pi*(d-y))*(1+M1*sqrt(1-y/d)+M2*(1-y/d)+M3*pow((1-y/d),1.5));
29 
30  return beta;
31 }

◆ EnthalpyToThermal()

void EnthalpyToThermal ( IssmDouble *  ptemperature,
IssmDouble *  pwaterfraction,
IssmDouble  enthalpy,
IssmDouble  pressure 
)

◆ printarray() [1/3]

void printarray ( IssmPDouble *  array,
int  lines,
int  cols = 1 
)

Definition at line 6 of file PrintArrays.cpp.

6  {
7  _printf_("\n");
8  for(int i=0;i<lines;i++){
9  _printf_(" [ ");
10  for(int j=0;j<cols;j++) _printf_( " " << setw(11) << setprecision (5) << array[i*cols+j]);
11  _printf_(" ]\n");
12  }
13  _printf_("\n");
14 }

◆ printarray_matlab()

void printarray_matlab ( const char *  filename,
int *  array,
int  lines,
int  cols = 1 
)

Definition at line 52 of file PrintArrays.cpp.

52  {
53  FILE *f = fopen(filename,"w");
54  fprintf(f,"%% Matrix of size %ix%i\n",lines,cols);
55  fprintf(f,"\n");
56  fprintf(f,"A=[...\n");
57  for(int i=0;i<lines;i++){
58  for(int j=0;j<cols;j++) fprintf(f," %i",array[i*cols+j]);
59  fprintf(f,"\n");
60  }
61  fprintf(f,"];\n");
62  fclose(f);
63 }

◆ printarray() [2/3]

void printarray ( int *  array,
int  lines,
int  cols = 1 
)

Definition at line 43 of file PrintArrays.cpp.

43  {
44  _printf_("\n");
45  for(int i=0;i<lines;i++){
46  _printf_(" [ ");
47  for(int j=0;j<cols;j++) _printf_( " " << setw(11) << setprecision (5) << array[i*cols+j]);
48  _printf_(" ]\n");
49  }
50  _printf_("\n");
51 }

◆ printarray() [3/3]

void printarray ( bool *  array,
int  lines,
int  cols = 1 
)

Definition at line 64 of file PrintArrays.cpp.

64  {
65  _printf_("\n");
66  for(int i=0;i<lines;i++){
67  _printf_(" [ ");
68  for(int j=0;j<cols;j++){
69  if(array[i*cols+j]) _printf_( " 1");
70  else _printf_( " 0");
71  }
72  _printf_(" ]\n");
73  }
74  _printf_("\n");
75 }

◆ printsparsity()

void printsparsity ( IssmPDouble *  array,
int  lines,
int  cols = 1 
)

Definition at line 26 of file PrintArrays.cpp.

26  {
27  int count;
28  _printf_("\n");
29  for(int i=0;i<lines;i++){
30  _printf_(" [ ");
31  count = 0;
32  for(int j=0;j<cols;j++){
33  if(array[i*cols+j]!=0.0){
34  _printf_( " X"); count++;
35  }
36  else
37  _printf_( " .");
38  }
39  _printf_(" ] "<<i<<" => "<<count << "\n");
40  }
41  _printf_("\n");
42 }

◆ printbinary()

void printbinary ( int  n)

Definition at line 76 of file PrintArrays.cpp.

76  {
77  unsigned int i=1L<<(sizeof(n)*8-1);
78  while (i>0) {
79  if (n&i)
80  _printf_("1");
81  else
82  _printf_("0");
83  i>>=1;
84  }
85 }
melt
void melt(IssmDouble *pM, IssmDouble *pMs, IssmDouble *pR, IssmDouble *pmAdd, IssmDouble *pdz_add, IssmDouble **pT, IssmDouble **pd, IssmDouble **pdz, IssmDouble **pW, IssmDouble **pa, IssmDouble **pre, IssmDouble **pgdn, IssmDouble **pgsp, int *pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid)
Definition: Gembx.cpp:1374
Cuffey
IssmDouble Cuffey(IssmDouble temperature)
Definition: Cuffey.cpp:11
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
EstarOmega
void EstarOmega(IssmDouble *omega, IssmDouble vx, IssmDouble vy, IssmDouble vz, IssmDouble vmag, IssmDouble *dvx, IssmDouble *dvy, IssmDouble *dvz, IssmDouble *dvmag)
Definition: EstarComponents.cpp:72
PI
const double PI
Definition: constants.h:11
R
const double R
Definition: Gembx.cpp:30
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
bamg::D2
P2< double, double > D2
Definition: Metric.h:11
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
Generated on Thu Jul 2 2020 08:09:24 for Ice Sheet System Model by   doxygen 1.8.19