forked from VirtualPlanetaryLaboratory/vplanet
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathvplanet.h
2416 lines (2211 loc) · 116 KB
/
vplanet.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/**
@file vplanet.h
@brief The main entry point for the code. All the magic starts here.
@author Rory Barnes ([RoryBarnes](https://github.com/RoryBarnes/))
@date May 7 2014
*/
#include <assert.h>
#include <ctype.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <time.h>
#ifdef __x86_64__
#include <xmmintrin.h>
#endif
// Windows-specific
#ifdef VPLANET_ON_WINDOWS
#define unlink _unlink
#else
#include <unistd.h>
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/*! Top-level declarations */
/* Implemented Moduules
The number is a bit value that can be used to uniquely identify the modules
that have been applied to a specific body. The value is stored in
module.iModuleBitSum.
*/
#define EQTIDE 2
#define RADHEAT 4
#define ATMESC 8
#define DISTORB 16
#define DISTROT 32
#define STELLAR 64
#define THERMINT 128
#define POISE 256
#define FLARE 512
#define BINARY 1024
#define GALHABIT 2048
#define SPINBODY 4096
#define DISTRES 8192
#define MAGMOC 16384
/* Fundamental constants; Some of these are taken from the IAU working
group on Fundamental constants, as described in Prsa et al. 2016. */
/* Units: Calculations are done in SI */
#define BIGG \
6.67428e-11 // From Luzum et al., 2011; value recommended \
// by IAU NSFA in Prsa et al. 2016
#define PI M_PI // Quick fix
#define KGAUSS 0.01720209895 // Gauss' Gravitational Constamt
#define EPS 1e-10 // Precision for difference of doubles to be effectively 0
#define AUM 1.49597870700e11 // Exact m/AU per 31 AUG 2012 IAU resolution B2
#define AUPC 206265.0 // AU per parsec
#define LIGHTSPEED 299792458.0 // Speed of light
#define MEARTH 5.972186e24 // Earth's mass; Prsa et al. 2016
#define MSUN 1.988416e30 // Sun's mass; Prsa et al. 2016
#define RSUN 6.957e8 // Sun's radius; Prsa et al. 2016
#define YEARSEC 3.15576e7 // Seconds per year
#define DAYSEC 86400 // Seconds per day
#define REARTH 6.3781e6 // Earth's Equatorial Radius; Prsa et al. 2016
#define RJUP 7.1492e7 // Jupiter's Equatorial Radius; Prsa et al. 2016
#define RNEP 2.4764e7 // Neptune's Radius
#define MNEP 1.0244e26 // Neptune's Mass
#define RHOEARTH 5515 // Earth's Density
#define eEARTH 0.016710219 // Earth's Eccentricity
#define MJUP 1.898130e27 // Jupiter's mass; Prsa et al. 2016
#define YEARDAY 365.25 // Days per year
#define MSAT 5.6851e26 // Saturns' Mass
#define DEGRAD 0.017453292519444445 // Degrees per radian
#define TOMASS 1.39e21 // Mass of one terrestrial ocean in kg (TO)
#define ATOMMASS 1.660538921e-27 // Atomic Mass
#define OXYMASS 2.6568622736e-26 // Mass of Atomic Oxygen 16*ATOMMASS
#define PROTONMASS 1.6726219e-27 // Proton mass kg
#define SIGMA 5.670367e-8 // Stefan-Boltzmann Constant
#define RGAS 8.3144598 // gas constant in J K^-1 mol^-1
#define KBOLTZ 1.38064852e-23 // Boltzmann constant, J/K
#define ALPHA_STRUCT \
0.6 // Structural constant for spherical mass \
//distribution potential energy (E_pot = -ALPHA*BIGG*M^2/R)
// Angle unit IDs
#define U_RADIANS 0
#define U_DEGREES 1
// Length unit IDs
#define U_METER 0
#define U_CENTIMETER 1
#define U_KILOMETER 2
#define U_SOLARRADIUS 3
#define U_EARTHRADIUS 4
#define U_JUPRADIUS 5
#define U_AU 6
// Mass unit IDs
#define U_GRAM 0
#define U_KILOGRAM 1
#define U_SOLARMASS 2
#define U_EARTHMASS 3
#define U_JUPITERMASS 4
#define U_NEPTUNEMASS 5
// Temperature unit IDs
#define U_KELVIN 0
#define U_CELSIUS 1
#define U_FARENHEIT 2
// Time unit IDs
#define U_SECOND 0
#define U_DAY 1
#define U_YEAR 2
#define U_MYR 3
#define U_GYR 4
/* Do not change these declarations */
extern const double dHUGE;
extern const double dTINY;
/* Do not change these declarations */
// IO limits for files, lines, and names
#define OPTLEN 48 /* Maximum length of an option */
#define OPTDESCR 128 /* Number of characters in option description */
#define OPTLONDESCR 2048 /* Number of characters in option long description */
/* !!! Hack to get long description and module list formatting for Long Help.
It'd be better to keep LINE back at 256 for memory, and malloc the char's,
but I don't know how to do that. */
#define LINE 2048 /* Maximum number of characters in a line */
#define NAMELEN 100
#define MAXFILES 128 /* Maximum number of input files */
#define MAXARRAY \
128 /* Maximum number of options in \
* an option array */
#define NUMOPT \
1000 /* Number of options that could be \
* in MODULE */
#define MAXLINES \
256 /* Maximum Number of Lines in an \
* input file */
#define OUTLEN \
48 /* Maximum number of characters in an output column header \
*/
#define OUTDESCR 256 /* Number of characters in output description */
#define OUTLONDESCR 2048 /* Number of characters in output long description */
/* Forward declaration of structs.
This is necessary in order to add pointers to structs into typedef'd functions
*/
typedef struct BODY BODY;
typedef struct CONTROL CONTROL;
typedef struct EVOLVE EVOLVE;
typedef struct FILES FILES;
typedef struct HALT HALT;
typedef struct INFILE INFILE;
typedef struct IO IO;
typedef struct MODULE MODULE;
typedef struct OPTIONS OPTIONS;
typedef struct OUTFILE OUTFILE;
typedef struct OUTPUT OUTPUT;
typedef struct SYSTEM SYSTEM;
typedef struct UNITS UNITS;
typedef struct UPDATE UPDATE;
typedef struct VERIFY VERIFY;
/*! \brief BODY contains all the physical parameters for every object in the
* system.
*/
struct BODY {
/* Body Properties */
char cName[NAMELEN]; /**< Body's Name */
char sColor[OPTLEN]; /**< Body color (for plotting) */
int bMantle; /**< Is there a mantle? */
int bOcean; /**< Is there an ocean? */
int bEnv; /**< Is there an envelope? */
int iBodyType; /**< Type of object: 0=star, 1=rocky planet, 2 = giant */
double dAge; /**< Body's Age */
double dMass; /**< Body's Mass */
double dSolidMass; /**< Mass of a body's solid component */
double dRadius; /**< Radius of body */
double dDensity; /**< Bulk density of body*/
double dGravAccel; /**< Body's gravitational acceleration */
double dK2; /**< Body's Total Love number */
double dImK2; /**< Imaginary part of Love's k_2 (total) */
double dObliquity; /**< Body's Obliquity */
double dCosObl; /**< Cosine of body's obliquity */
double dRotRate; /**< Body's Rotation Rate */
double dRotPer; /**< Body's Rotation Period */
double dRotVel; /**< Body's Rotational Velocity */
double dRadGyra; /**< Body's Radius of Gyration */
/* Orbital Properties. By convention, these are stored in the
* second element in the BODY array and, if using binary
* in the secondary (1st (0, 1, ..)) body*/
double dSemi; /**< Body's Semi-major Axis */
double dEcc; /**< Body's Eccentricity */
double dMeanMotion; /**< Body's Mean Motion */
double dOrbPeriod; /**< Body's Orbital Period */
double dEccSq; /**< Eccentricity squared */
/* ATMESC Parameters */
int bAtmEsc; /**< Apply Module ATMESC? */
int bInstantO2Sink; /**< Is oxygen absorbed instantaneously at the surface? */
int bRunaway; /**< Is the planet experiencing a runaway greenhouse? */
int bCalcFXUV; /**< Does incidenx XUV flow need to be calculated every
time step? */
int bEnvelopeLostMessage; /**< Has the envelope lost message been printed? */
int bIgnoreRocheLobe; /**< Ignore Roche lobe overflow? */
int bUseEnergyLimited; /**< Use energy-limited escape */
int bUseBondiLimited; /**< Use Bondi-limited H mass loss */
int bUseRRLimited; /**< Use radiation/recombination-limited H mass loss */
int bAtmEscAuto; /**< Transition H escape regime depending on physics */
int bAutoThermTemp; /**< Calculate thermal temperature from environemnt? */
int bStopWaterLossInHZ; /**< Stop water loss once planet enters habitable
zone? */
int iWaterLossModel; /**< Water Loss and Oxygen Buildup Model */
int iAtmXAbsEffH2OModel; /**< Water X-ray/XUV absorption efficiency evolution
model */
int iPlanetRadiusModel; /**< Planet Radius model. */
int iWaterEscapeRegime; /**< Track water escape regime */
int iHEscapeRegime; /**< Tracks H escape regime */
double dSurfaceWaterMass; /**< Surface water mass */
double dMinSurfaceWaterMass; /**< Minimum surface water to avoid a halt */
double dOxygenMass; /**< Atmospheric oxygen mass */
double dOxygenMantleMass; /**< Mass of oxygen in the mantle */
double dEnvelopeMass; /**< Envelope mass */
double dMinEnvelopeMass; /**< Minimum envelope mass to avoid a halt */
double dXFrac; /**< XUV radius in units of planet's radius */
double dAtmXAbsEffH; /**< Effective XUV absorpation efficiency for hydrogen */
double dAtmXAbsEffH2O; /**< Effective XUV absorpation efficiency for water */
double dRGDuration; /**< Duration of runaway greenhouse phase */
double dKTide; /**< Tidal enhancement factor for mass loss */
double dMinKTide; /**< Minimum allowed value for KTide */
double dAtmEscXi; /**< Ratio of Roche radius to XUV radius */
double dMDotWater; /**< Water mass loss rate */
double dFHRef; /**< Reference hydrogen escape value */
double dOxygenEta; /**< Factor for drag of oxygen by hydrogen */
double dCrossoverMass; /**< Atomic mass at which drag begins */
double dFHDiffLim; /**< Diffusion-limited H escape rate */
double dRadXUV; /**< XUV Radius in Lehmer-Catling model */
double dRadSolid; /**< Solid planet radius in Lehmer-Catling model */
double dPresSurf; /**< Surface Pressure in Lehmer-Catling model */
double dPresXUV; /**< Pressure at XUV radius in Lehmer-Catling model */
double dScaleHeight; /**< Scale Height used in Lehmer-Catling model */
double dThermTemp; /**< Thermosphere's temperature in Lehmer-Catling model */
double dAtmGasConst; /**< Atmosphere's gas constant in Lehmer-Catling model */
double dFXUV; /**< XUV Flux at planet's atmosphere */
double dJeansTime; /**< Jeans timescale for atmospheric escape */
double dFlowTemp; /**< Temperature of the hydrodynamic flow */
double dRocheRadius; /**< Radius of the Roche lobe */
double dBondiRadius; /**< Bondi (Sonic) Radius */
double dEnvMassDt; /**< Time derivative of H envelope mass */
/* BINARY parameters */
int bBinary; /**< Apply BINARY module? */
double dR0; /**< Guiding Radius,initially equal to dSemi */
double dCBPR; /**< CBP orbital radius */
double dCBPZ; /**< CBP height above/below the orbital plane */
double dCBPPhi; /**< CBP azimuthal angle in orbital plane */
double dCBPRDot; /**< CBP radial orbital velocity */
double dCBPZDot; /**< CBP z orbital velocity */
double dCBPPhiDot; /**< CBP phi angular orbital velocity */
double dFreeEcc; /**< CBP's free eccentricity */
double dFreeInc; /**< CBP's free inclination, or binary's inclination */
double dInc; /**< Orbital inclication */
double dLL13N0; /**< CBP's Mean motion defined in LL13 eqn 12 */
double
dLL13K0; /**< CBP's radial epicyclic frequency defined in LL13 eqn 26 */
double dLL13V0; /**< CBP's vertical epicyclic frequency defined in LL13 eqn 36
*/
double dLL13PhiAB; /**< Binary's initial mean anomaly */
double dCBPM0; /**< CBP's initial mean anomaly */
double dCBPZeta; /**< CBP's z oscillation angle (see LL13 eqn 35) */
double dCBPPsi; /**< CBP's R, phi oscillation phase angle (see LL13 eqn 27) */
/* SPINBODY parameters */
int bSpiNBody; /**< Has module SPINBODY been implemented */
double dVelX; /**< x Component of the body's velocity */
double dVelY; /**< y Component of the body's velocity */
double dVelZ; /**< z Component of the body's velocity */
double dPositionX; /**< x Component of the body's position */
double dPositionY; /**< y Component of the body's position */
double dPositionZ; /**< z Component of the body's position */
double bUseOrbParams; /**< Boolean flag to use orbital parameters as inputs */
double *dDistance3; /**< Distance cubed to different perturbers */
double *dDistanceX; /**< X Distance between two bodies */
double *dDistanceY; /**< Y Distance between two bodies */
double *dDistanceZ; /**< Z Distance between two bodies */
double *dHCartPos; /**< Heliocentric Cartesian Position used for orbital
element calculations */
double *dHCartVel; /**< Heliocentric Cartesian Velocity used for orbital
element calculations */
double *dBCartPos; /**< Barycentric Cartesian Position used for orbital
element calculations */
double *dBCartVel; /**< Barycentric Cartesian Velocity used for orbital
element calculations */
double dGM; /**< GM for the star */
double dMu; /**< G(M+m) */
int iGravPertsSpiNBody; /**< Number of bodies that are orbitally relevent
(equal to for evolve->iNumBodies) */
/* DISTORB parameters */
int bDistOrb; /**< Has module DISTORB been implemented */
double dHecc; /**< Poincare H */
double dKecc; /**< Poincare K */
double dPinc; /**< Poincare P */
double dQinc; /**< Poincare Q */
double dSinc; /**< sin(0.5*Inclination) */
double dLongA; /**< Longitude of ascending node */
double dArgP; /**< Argument of pericenter */
double dLongP; /**< Longitude of pericenter */
double dMeanA; /**< Mean anomaly (only used for inv plane calculation) */
double dTrueL; /**< True longitude (only used for insolation calculation */
double dEccA; /**< Eccentric anomaly (only used for inv plane calculation) */
double *daCartPos; /**< Cartesian position of body (only used for inv plane
calc) */
double *daCartVel; /**< Cartesian velocity of body (only used for inv plane
calc) */
int iGravPerts; /**< Number of bodies which perturb the body */
int *iaGravPerts; /**< Which bodies are perturbers of the body */
int iEigFreqs; /**< Number of eigenfrequencies that control the body's motion
*/
int *iaEigFreqs; /**< Indices of eigenfrequencies */
int bGRCorr; /**< Use general relativistic correction in DistOrb+DistRot?*/
int iDistOrbModel; /**< Which orbital model to use (RD4 or LL2) */
double
dSemiPrev; /**< Semi-major axis at which LL2 eigensolution was calc'd */
double dEigenvalue; /**< User input eigenvalue (diagnostic only) */
double
dEigenvector; /**< User input eigenvector amplitude (diagnostic only) */
int bEigenSet; /**< Manually set an eigenvalue/frequency */
double *daLOrb; /**< Orbital angular momentum */
double *daLOrbTmp; /**< Temp copy of orbital angular momentum */
double dRPeri; /**< Pericenter distance */
double dRApo; /**< Apocenter distance */
/* DISTROT parameters */
int bDistRot;
double dPrecA; /**< Precession angle */
double dTrueApA; /**< True anomaly at equinox (used for invariable plane
conversion) */
double dDynEllip; /**< Dynamical ellipticity */
double dYobl; /**< sin(obliq)*sin(preca) */
double dXobl; /**< sin(obliq)*cos(preca) */
double dZobl; /**< cos(obliq) */
double *daLRot; /**< Spin angular momentum vector */
double *daLRotTmp; /**< Temp copy of spin angular momentum vector */
int bForcePrecRate; /**< Set precession rate to a fixed value */
double dPrecRate; /**< Value to set fixed precession rate to */
int bCalcDynEllip; /**< Calc dyn ellipticity from spin, radius, mass, inertia?
*/
int bRelaxDynEllip; /**< shape of planet relaxes when spun down */
int bReadOrbitData; /**< Use orbit data from file rather than distorb */
char cFileOrbitData[NAMELEN]; /**< read orbital data from this file
(distorb=0) */
double *daTimeSeries; /**< time series for orbital data */
double *daSemiSeries; /**< time series for orbital data */
double *daEccSeries; /**< time series for orbital data */
double *daIncSeries; /**< time series for orbital data */
double *daArgPSeries; /**< time series for orbital data */
double *daLongASeries; /**< time series for orbital data */
double *daMeanASeries; /**< time series for orbital data */
int iCurrentStep; /**< index for time series arrays */
double *daHeccSeries; /**< time series for orbital data */
double *daKeccSeries; /**< time series for orbital data */
double *daPincSeries; /**< time series for orbital data */
double *daQincSeries; /**< time series for orbital data */
double dPdot; /**< inclination derivative used for obliquity evol */
double dQdot; /**< inclination derivative used for obliquity evol */
int iNLines; /**< Number of lines of orbital data file */
double dSpecMomInertia; /**< C/M/R^2 used for dynamical ellipticity
calculation */
/* EQTIDE Parameters */
int bEqtide; /**< Apply Module EQTIDE? */
int bTideLock; /**< Is a body tidally locked? */
double dLockTime; /**< Time when body tidally-locked */
int bUseTidalRadius; /**< Set a fixed tidal radius? */
int bUseOuterTidalQ; /**< Set total Q to outer layer's value? */
double dTidalRadius; /**< Radius used by tidal evoltion equations (CPL only
currently) */
int iTidePerts; /**< Number of Tidal Perturbers */
int *iaTidePerts; /**< Body #'s of Tidal Perturbers */
char saTidePerts[MAXARRAY][NAMELEN]; /**< Names of Tidal Perturbers */
double dK2Man; /**< Mantle k2 love number */
double dK2Ocean; /**< Ocean's Love Number */
double dK2Env; /**< Envelope's Love Number */
double dTidalQMan; /**< Tidal Q of the Mantle */
double dTidalQOcean; /**< Body's Ocean Component to Tidal Q */
double dTidalQEnv; /**< Body's Envelope Component to Tidal Q */
double dImK2Man; /**< Mantle Im(k2) love number */
double dImK2ManOrbModel; /**< Mantle Im(k2) model for DB15 orbital eqns */
double dImK2Ocean; /**< Envelope Component to Imaginary part of Love's K_2 */
double dImK2Env; /**< Envelope Component to Imaginary part of Love's K_2 */
double dTidalQ; /**< Body's Tidal Q */
double dTidalTau; /**< Body's Tidal Time Lag */
double *dTidalZ; /**< As Defined in \cite HellerEtal2011 */
double *dTidalChi; /**< As Defined in \cite HellerEtal2011 */
double **dTidalF; /**< As Defined in \cite HellerEtal2011 */
double *dTidalBeta; /**< As Defined in \cite HellerEtal2011 */
int **iTidalEpsilon; /**< Signs of Phase Lags */
double dDeccDtEqtide; /**< Eccentricity time rate of change */
double *daDoblDtEqtide; /**< Obliquity time rate of change */
/* RADHEAT Parameters: H = Const*exp[-Time/HalfLife] */
int bRadheat; /**< Apply Module RADHEAT? */
double d26AlConstMan; /**< Body's Mantle 26Al Decay Constant */
double d26AlMassMan; /**< Body's Mantle Mass of 26Al */
double d26AlNumMan; /**< Body's Mantle Number of 26Al Atoms */
double d26AlPowerMan; /**< Body's Mantle Internal Power Due to 26Al Decay */
double d26AlConstCore; /**< Body's Core 26Al Decay Constant */
double d26AlMassCore; /**< Body's Core Mass in 26Al */
double d26AlNumCore; /**< Body's Core Number of 26Al Atoms */
double d26AlPowerCore; /**< Body's Core Power from 26Al */
double d40KConstMan; /**< Body's Mantle 40K Decay Constant */
double d40KMassMan; /**< Body's Mantle Mass of 40K */
double d40KNumMan; /**< Body's Mantle Number of 40K Atoms */
double d40KPowerMan; /**< Body's Mantle Internal Power Due to 40K Decay */
double d40KConstCore; /**< Body's Core 40K Decay Constant */
double d40KNumCore; /**< Body's Core Number of 40K Atoms */
double d40KPowerCore; /**< Body's Core Power due to 40K */
double d40KMassCore; /**< Body's Core Mass of 40K */
double d40KConstCrust; /**< Body's Crust 40K Decay Constant */
double d40KNumCrust; /**< Body's Crust Number of 40K Atoms */
double d40KPowerCrust; /**< Body's Crust Power due to 40K */
double d40KMassCrust; /**< Body's Crust Mass of 40K */
double d232ThConstMan; /**< Body's Thorium-232 Decay Constant */
double d232ThNumMan; /**< Body's Number of Thorium-232 Atoms */
double d232ThPowerMan; /**< Body's Internal Power Due to Thorium-232 Decay */
double d232ThMassMan; /**< Body's Total Mass of Thorium-232 Atoms */
double d232ThConstCore;
double d232ThNumCore;
double d232ThPowerCore;
double d232ThMassCore;
double d232ThConstCrust;
double d232ThNumCrust;
double d232ThPowerCrust;
double d232ThMassCrust;
double d238UConstMan; /**< Body's Uranium-238 Decay Constant */
double d238UNumMan; /**< Body's Number of Uranium-238 Atoms */
double d238UPowerMan; /**< Body's Internal Power Due to Uranium-238 Decay */
double d238UMassMan; /**< Body's Total Mass of Uranium-238 Atoms */
double d238UConstCore;
double d238UNumCore;
double d238UPowerCore;
double d238UMassCore;
double d238UConstCrust;
double d238UNumCrust;
double d238UPowerCrust;
double d238UMassCrust;
double d235UConstMan;
double d235UNumMan;
double d235UPowerMan;
double d235UMassMan;
double d235UConstCore;
double d235UNumCore;
double d235UPowerCore;
double d235UMassCore;
double d235UConstCrust;
double d235UNumCrust;
double d235UPowerCrust;
double d235UMassCrust;
double dRadPowerTotal; /**< Total planet Radiogenic Power */
double dRadPowerMan; /**< Total Mantle Radiogenic Power */
double dRadPowerCore; /**< Total Core Radiogenic Power */
double dRadPowerCrust; /**< Total Crust Radiogenic Power */
/* Thermint Parameters */
int bThermint; /**< Apply Module THERMINT? */
double dTSurf; /**< Surface Temperature */
double dTMan; /**< Temperature Mantle AVE */
double dTCore; /**< Temperature Core AVE */
double dTUMan; /**< Temperature UMTBL */
double dTLMan; /**< Temperature LMTBL */
double dTCMB; /**< Temperature CMB */
double dTICB; /**< Temperature ICB */
double dBLUMan; /**< UM TBL thickness */
double dBLLMan; /**< LM TBL thickness */
double dTJumpUMan; /**< Abs Temperature Jump across UMTBL */
double dTJumpLMan; /**< Abs Temperature Jump across LMTBL */
double dSignTJumpUMan; /**< Sign of Temperature Jump across UMTBL */
double dSignTJumpLMan; /**< Sign of Temperature Jump across LMTBL */
double dViscUManArr; /**< Viscosity UMTBL Arrhenius Law */
double dViscUMan; /**< Viscosity UMTBL */
double dViscLMan; /**< Viscosity LMTBL */
double dViscMMan; /**< Viscosity Mid (ave) mantle */
double dDynamViscos; /**< Dynamic viscosity of the mantle */
double dViscJumpMan; /**< Viscosity Jump UM to LM */
double dShmodUMan; /**< Shear modulus UMTBL */
double dShmodLMan; /**< Shear modulus LMTBL */
double dTsolUMan; /**< Solidus Temperature UMTBL */
double dTliqUMan; /**< Liquidus Temperature UMTBL */
double dTsolLMan; /**< Solidus Temperature LMTBL */
double dTliqLMan; /**< Liquidus Temperature LMTBL */
double dFMeltUMan; /**< Melt fraction UMTBL */
double dFMeltLMan; /**< Melt fraction LMTBL */
double dMeltfactorUMan; /**< Melt Phase Factor for Rheology */
double dMeltfactorLMan; /**< Melt Phase Factor for Rheology */
double dFixMeltfactorUMan; /**< Melt Phase Factor for Rheology */
double dViscMeltB; /**< Viscosity Melt Factor B */
double dViscMeltGamma; /**< Viscosity Melt Factor Gamma */
double dViscMeltDelta; /**< Viscosity Melt Factor Delta */
double dViscMeltXi; /**< Viscosity Melt Factor Xi */
double dViscMeltPhis; /**< Viscosity Melt Factor Phis */
double dDepthMeltMan; /**< Depth to base of UM Melt layer */
double dTDepthMeltMan; /**< Temp at base of UM Melt layer */
double dTJumpMeltMan; /**< Temp Jump to base of UM Melt layer */
double dMeltMassFluxMan; /**< Mantle upwelling melt mass flux */
double dRayleighMan; /**< Mantle Rayleigh Number */
/* Time Derivatives & Gradients */
double dTDotMan; /**< Time deriv of mean mantle temp */
double dTDotCore; /**< time deriv of mean core temp */
double dHfluxUMan; /**< hflux upper mantle thermal boundary layer (UMTBL) */
double dHflowUMan; /**< hflow UMTBL */
double dHfluxLMan; /**< hflux lower mantle thermal boundary layer (UMTBL) */
double dHflowLMan; /**< hflow LMTBL */
double dHfluxCMB; /**< hflux CMB */
double dHflowCMB; /**< hflow CMB */
double dHflowTidalMan; /**< hflow tidal dissipation in mantle */
double dHflowTidalCore; /**< hflow tidal dissipation in core */
double dHflowLatentMan; /**< latent hflow from solidification of mantle */
double dHflowMeltMan; /**< Eruptive Melt Hflow from mantle */
double dHflowSecMan; /**< Mantle Secular Heat flow */
double dMassICDot; /**< Mass Growth Rate of IC */
double dHflowLatentIC; /**< latent hflow from solidification of IC */
double dPowerGravIC; /**< latent hflow from solidification of IC */
double dHflowICB; /**< hflow across ICB */
double dHfluxSurf; /**< hflux surface of mantle */
double dHflowSurf; /**< hflow surface of mantle */
double dTidalPowMan; /**< Tidal Dissipation Power in Mantle */
/* Core Variables */
double dRIC; /**< IC radius */
double dDRICDTCMB; /**< d(R_ic)/d(T_cmb) */
double dDOC; /**< OC shell thickness */
double dThermConductOC; /**< Thermal conductivity OC */
double dThermConductIC; /**< Thermal conductivity IC */
double dChiOC; /**< OC light element concentration chi. */
double dChiIC; /**< IC light element concentration chi. */
double dMassOC; /**< OC Mass. */
double dMassIC; /**< IC Mass. */
double dMassChiOC; /**< OC Chi Mass. */
double dMassChiIC; /**< IC Chi Mass. */
double dDTChi; /**< Core Liquidus Depression */
double dHfluxCMBAd; /**< CMB Adiabatic Heat flux. */
double dHfluxCMBConv; /**< CMB Convective (super-adiabatic) Heat flux. */
double dCoreBuoyTherm; /**< Core Thermal buoyancy flux */
double dCoreBuoyCompo; /**< Core Compositional buoyancy flux */
double dCoreBuoyTotal; /**< Core total (therm+compo) buoyancy flux */
double dGravICB; /**< Gravity at ICB */
double dDensAnomICB; /**< Density anomaly across ICB (Delta rho_chi in DB14).
*/
double dRICDot; /**< Inner core growth rate */
/* Magnetic Field */
double dMagMom; /**< Core Dynamo Magnetic Moment scaling law. */
double dMagMomCoef; /**< Dynamo magnetic moment scaling law dipolarity
coefficient (gamma_d in DB14) */
double dPresSWind; /**< Stellar wind pressure at planets orbit. */
double dMagPauseRad; /**< Magnetopause stand-off radius from center of planet
*/
/* Constants */
double dViscRatioMan; /**< Viscosity Ratio Man */
double dEruptEff; /**< Mantle melt eruption efficiency */
double dViscRef; /**< Mantle Viscosity Reference (coefficient) */
double dTrefLind; /**< Core Liquidus Lindemann Reference (coefficient) */
double dDTChiRef; /**< Core Liquidus Depression Reference (E) */
double dStagLid; /**< Stagnant Lid heat flow switch (0 or 1)*/
double dManHFlowPref; /**< Mantle Hflow Prefix */
double dActViscMan; /**< Mantle viscosity activation energy */
double dShModRef; /**< reference kinematic mantle shear modulus */
double dStiffness; /**< effective stiffness of mantle */
double dDLind; /**< lindemann's law length scale for iron liquidus*/
double dDAdCore; /**< liq iron core adiabatic length scale */
double dAdJumpM2UM; /**< adiabatic temp jump from ave mantle to UM */
double dAdJumpM2LM; /**< adiabatic temp jump from ave mantle to LM */
double dAdJumpC2CMB; /**< adiabatic temp jump from ave core to CMB */
double dElecCondCore; /**< electrical conductivity of core */
/* STELLAR Parameters */
int bStellar;
double dLuminosity;
double dTemperature;
double dSatXUVFrac;
double dSatXUVTime;
double dXUVBeta;
int iStellarModel;
int iMagBrakingModel;
int iWindModel;
int iXUVModel;
double dLXUV; // Not really a STELLAR parameter
double iHZModel;
double dLostAngMom; /**< Angular momemntum lost to space via magnetic braking
*/
double dLostEng; /**< Energy lost to space, i.e. via stellar contraction */
int bRossbyCut; /**< Whether or not to shut off magnetic braking for
Ro>ROSSBYCRIT */
int bEvolveRG; /**< Whether or not to evolve radius of gyration? Defaults to 0
*/
double dLuminosityInitial;
double dLuminosityAmplitude;
double dLuminosityFrequency;
double dLuminosityPhase;
/* POISE parameters */
int bPoise; /**< Apply POISE module? */
double dAblateFF; /**< Scaling factor for ice ablation rate */
int bAccuracyMode; /**< This forces EBM to re-invert matrix every time step */
double dAlbedoGlobal; /**< Global average albedo (Bond albedo) */
double dAlbedoGlobalTmp; /**< A copy of global average albedo (sometimes
needed) */
double dAlbedoLand; /**< Sets base albedo of land (sea model) */
double dAlbedoWater; /**< Sets base albedo of water (sea model) */
int bAlbedoZA; /**< Use albedo based on zenith angle (ann model) */
double dAreaIceCov; /**< Tracks area of surface covered in permanent ice*/
double dAstroDist; /**< Distance between primary and planet */
int bCalcAB; /**< Calc A and B from Williams & Kasting 1997 */
int iClimateModel; /**< Which EBM to be used (ann or sea) */
int bColdStart; /**< Start from global glaciation (snowball) conditions */
double dCw_dt; /**< Heat capacity of water / EBM time step */
double dDiffCoeff; /**< Diffusion coefficient set by user */
int bDiffRot; /**< Adjust heat diffusion for rotation rate */
int bElevFB; /**< Apply elevation feedback to ice ablation */
double dFixIceLat; /**< Fixes ice line latitude to user set value */
double dFluxInGlobal; /**< Global mean of incoming flux */
double dFluxInGlobalTmp; /**< Copy of global mean incoming flux */
double dFluxOutGlobal; /**< Global mean of outgoing flux */
double dFluxOutGlobalTmp; /**< Copy of global mean outgoing flux */
int bForceObliq; /**< Force obliquity to evolve sinusoidally */
int bForceEcc; /**< Force eccentricity to evolve sinusoidally */
double dFrzTSeaIce; /**< Freezing temperature of sea water */
int iGeography; /**< Type of geography to use (uni3 or modn) */
int bHadley; /**< Use Hadley circ in tropics when calc'ing diffusion? */
double dHeatCapAnn; /**< Surface heat capacity in annual model */
double dHeatCapLand; /**< Heat capacity of land */
double dHeatCapWater; /**< Heat capacity of water */
double dIceAlbedo; /**< Base albedo of ice covered surfaces */
double dIceBalanceTot; /**< Total gain/loss in ice globally */
double dIceDepRate; /**< Snow deposition rate when below freezing */
double dIceFlowTot; /**< Total flow of ice (should be zero) */
double dIceMassTot; /**< Total ice mass over entire globe */
int bIceSheets; /**< Use ice sheet model? */
int iIceTimeStep; /**< Time step of ice sheet model (should be > iNumYears) */
double dInitIceHeight; /**< Initial height of ice sheet */
double dInitIceLat; /**< Initial latitude of ice line (ice cap only) */
double dLapseR; /**< Lapse rate used for elevation feedback of ice sheet */
double dLatentHeatIce; /**< Latent heat of fusion of ice over mixing depth*/
double dLatFHeatCp; /**< Latent heat of ice/heat capacity */
int bMEPDiff; /**< Compute diff from maximum entropy prod (D = B/4) */
double dMixingDepth; /**< Depth of mixing layer of ocean (for thermal
inertia)*/
int iNDays; /**< Number of days in planet's year/orbit */
int iNStepInYear; /**< Number of time steps in a year/orbit */
double dNuLandWater; /**< Land-ocean interaction term */
int iNumLats; /**< Number of latitude cells */
int iNumYears; /**< Number of orbits!!! to run seasonal model */
double dObliqAmp; /**< Amplitude of forced obliquity oscillation */
double dObliqPer; /**< Period of force obliquity oscillation */
double dObliq0; /**< Start obliquity for forced oscillation */
double dEccAmp; /**< Amplitude of forced eccentricity oscillation */
double dEccPer; /**< Period of force eccentricity oscillation */
double dEcc0; /**< Start eccentricity for forced oscillation */
int iOLRModel; /**< OLR fit (use with bCalcAB=1) from Kasting or Spiegel */
double dpCO2; /**< Partial pressure of CO2 (only if bCalcAB = 1) */
double dPlanckA; /**< Constant term in Blackbody linear approximation */
double dPlanckB; /**< Linear coeff in Blackbody linear approx (sensitivity) */
double dPrecA0; /**< Initial pA value used when distrot is not called */
int bReadOrbitOblData; /**< Use orbit and obliquity data from file rather
than distrot */
char cFileOrbitOblData[NAMELEN]; /**< read orbital and obliquity data from
this file (distorb=0) */
double *daOblSeries; /**< time series for obliquity data */
double *daPrecASeries; /**< time series for obliquity data */
double dRefHeight; /**< Ref height of "surface" in elevation feedback */
int iReRunSeas; /**< When to rerun EBM in ice sheet model */
double dSeaIceConduct; /**< Conductivity of sea ice */
int bSeaIceModel; /**< Use sea ice model? */
double dSeasDeltat; /**< Time step of seasonal model */
double dSeasDeltax; /**< Spacing of grid points in seasonal model */
double dSeasOutputTime; /**< When to output seasonal data */
double dSeasNextOutput; /**< Next time step to output seasonal data */
int bSkipSeas; /**< Ann model will be used if in snowball state */
int bSkipSeasEnabled; /**< Allow ann model to be used if in snowball state? */
int bSnowball; /**< Is planet in snowball state (oceans are frozen)? */
double
dSpinUpTol; /**< Tolerance for mean global temp change during spin up */
double dSurfAlbedo; /**< Base surface albedo used in ann model */
double dTGlobal; /**< Global mean temperature at surface */
double dTGlobalInit; /**< Initial estimate of global surface temperature */
double dTGlobalTmp; /**< Mean global surface temp */
int iWriteLat; /**< Stores index of latitude to be written in write fxn */
double dMinIceHeight; /**< Minimum ice thickness to count as icy */
/* Arrays used by seasonal and annual */
double *daAnnualInsol; /**< Annually averaged insolation at each latitude */
double *daDivFlux; /**< Divergence of surface flux */
double *daDMidPt; /**< Diffusion at edges of grid points */
double **daInsol; /**< Daily insolation at each latitude */
double *daFlux; /**< Meridional surface heat flux */
double *daFluxIn; /**< Incoming surface flux (insolation) */
double *daFluxOut; /**< Outgoing surface flux (longwave) */
double *daLats; /**< Latitude of each cell (centered); South Pole is 0 */
double *daPeakInsol; /**< Annually averaged insolation at each latitude */
double *daTGrad; /**< Gradient of temperature (meridional) */
/* Arrays for annual model */
double *daAlbedoAnn; /**< Albedo of each cell */
double
*daDiffusionAnn; /**< Diffusion coefficient of each latitude boundary */
double **daMEulerAnn; /**< Matrix used for Euler step in annual model */
double **daMEulerCopyAnn; /**< Temp copy of Euler matrix */
double **daInvMAnn; /**< Inverted matrix for annual model */
double *daLambdaAnn; /**< Diffusion terms for annual matrix */
double **daMClim; /**< Raw climate matrix for annual model */
double **daMDiffAnn; /**< Diffusion matrix for annual model */
double *daPlanckAAnn; /**< Array of Planck A values for ann model */
double *daPlanckBAnn; /**< Array of Planck B values for ann model */
int *iaRowswapAnn; /**< Array of interchanged rows in matrix inversion */
double *daScaleAnn; /**< Used in matrix inversion routine */
double *daSourceF; /**< Heating terms in EBM */
double *daTempAnn; /**< Surface temperature in each cell */
double *daTempTerms; /**< Temperature dependent terms in matrix */
double *daTmpTempAnn; /**< Temporary copy of temperature */
double *daTmpTempTerms; /**< Temporary copy of temp dependent terms */
double *daUnitVAnn; /**< Unit vector used in matrix inversion */
/* Arrays for seasonal model */
double *daAlbedoAvg; /**< Orbit average albedo by latitude */
double *daAlbedoAvgL; /**< Orbit average albedo by latitude on land */
double *daAlbedoAvgW; /**< Orbit average albedo by latitude on water */
double *daAlbedoLand; /**< Albedo of land by latitude */
double *daAlbedoLW; /**< Land-water averaged albedo */
double *daAlbedoWater; /**< Albedo of land by latitude */
double *daBasalFlow; /**< Basal flow of ice = d(u*h)/dy */
double *daBasalFlowMid; /**< Basal flow of ice d(u*h)/dy (midpoints) */
double *daBasalVel; /**< Basal velocity of ice */
double *daBedrockH; /**< Height of bedrock (can be negative) */
double *daBedrockHEq; /**< Equilibrium height of bedrock */
double *daDeclination; /**< Daily solar declination */
double *daDeltaTempL; /**< Keeps track of temp change on land for energy check
*/
double *daDeltaTempW; /**< Keeps track of temp change on water for energy
check */
double *daDIceHeightDy; /**< Gradient of ice height */
double *daDiffusionSea; /**< Diffusion coefficient for seasonal model */
double *daDivFluxAvg; /**< Divergence of flux averaged over orbit */
double **daDivFluxDaily; /**< Daily values of divergence of flux */
double *daEnergyResL; /**< Energy residuals on land */
double *daEnergyResW; /**< Energy residuals over water */
double *daEnerResLAnn; /**< Annually averaged energy residuals on land */
double *daEnerResWAnn; /**< Annually averaged energy residuals over water */
double *daFluxAvg; /**< Annually averaged meridional flux */
double *daFluxOutAvg; /**< Annually averaged outgoing flux */
double **daFluxDaily; /**< Daily meridional flux values */
double *daFluxInAvg; /**< Annually averaged incoming flux */
double **daFluxInDaily; /**< Daily incoming flux values */
double *daFluxInLand; /**< Annually averaged incoming flux on land */
double *daFluxInWater; /**< Annually averaged incoming flux on water */
double **daFluxOutDaily; /**< Daily outgoing flux values */
double *daFluxOutLand; /**< Annually averaged outgoing flux on land */
double *daFluxOutWater; /**< Annually averaged outgoing flux on water */
double *daFluxSeaIce; /**< Heat flux through sea ice */
double **daIceBalance; /**< Gain/loss of ice at each latitude and day */
double *daIceAblateTot; /**< Total ice loss per orbit */
double *daIceAccumTot; /**< Total ice gain per orbit */
double *daIceBalanceAnnual; /**< Net ice gain/loss over orbit */
double *daIceBalanceAvg; /**< Average ice gain/loss over orbit */
double *daIceBalanceTmp; /**< Temporary (current) ice gain/loss */
double *daIceFlow; /**< Flow of ice */
double *daIceFlowAvg; /**< Average flow of ice over orbit */
double *daIceFlowMid; /**< Flow of ice at boundaries of grid points */
double *daIceGamTmp; /**< Temporary variable used in ice sheet matrix */
double *daIceHeight; /**< Height of ice sheet */
double *daIceMass; /**< Ice mass per area */
double *daIceMassTmp; /**< Temporary copy of ice mass per area */
double *daIcePropsTmp; /**< Temporary array used in ice sheet matrix */
double *daIceSheetDiff; /**< Diffusion coefficient of ice sheet flow */
double **daIceSheetMat; /**< Matrix used in ice sheet flow */
double **daInvMSea; /**< Inverted matrix in seasonal EBM */
double *daLambdaSea; /**< Diffusion terms in seasonal EBM matrix */
double dLandFrac; /**< Land fraction input by user */
double *daLandFrac; /**< Fraction of cell which is land */
double **daMDiffSea; /**< Diffusion only matrix in seasonal EBM */
double **daMEulerCopySea; /**< Temporary copy of Euler time step matrix
(seasonal) */
double **daMEulerSea; /**< Euler time step matrix in seasonal EBM */
double **daMInit; /**< Temporary matrix used in constructing Euler matrix */
double **daMLand; /**< Land terms in seasonal matrix */
double **daMWater; /**< Water terms in seasonal matrix */
double *daPlanckASea; /**< Array of Planck A values in seasonal model */
double *daPlanckBSea; /**< Array of Planck B values in seasonal model */
double **daPlanckBDaily; /**< Array of Planck B values over seasonal cycle */
double *daPlanckBAvg; /**< Orbit averaged Planck B values in seasonal model */
int *iaRowswapSea; /**< Interchanged rows in seasonal matrix inversion */
double *daScaleSea; /**< Used in matrix inversion routine */
double *daSeaIceHeight; /**< Sea ice height by latitude */
double *daSeaIceK; /**< Heat conductivity of sea ice */
double *daSedShear; /**< sediment shear stress (for ice sheets) */
double *daSourceL; /**< Land heating terms: PlanckA - (1-albedo)*Insol */
double *
daSourceLW; /**< Combined heat terms what inverser matrix operates on */
double *daSourceW; /**< Water heating terms: PlanckA - (1-albedo)*Insol */
double *daTempAvg; /**< Temperature averaged over orbit and land/water */
double *daTempAvgL; /**< Land temp averaged over orbit */
double *daTempAvgW; /**< Water temp averaged over orbit */
double **daTempDaily; /**< Daily temp over seasonal cycle */
double *daTempLand; /**< Temperature over land (by latitude) */
double *daTempLW; /**< Surface temperature (avg over land & water) */
double *daTempMaxLW; /**< Maximum temperature over year */
double *daTempMaxLand; /**< Maximum temperature over year over land */
double *daTempMaxWater; /**< Maximum temperature over year over water */
double *daTempMinLW; /**< Minimum temperature over year */
double *daTempWater; /**< Temperature over ocean (by lat) */
double *daTmpTempSea; /**< Temporary copy of temp dependent terms (sea EBM)*/
double *daUnitVSea; /**< Unit vector used in matrix routines */
double *daWaterFrac; /**< Fraction of cell which is water */
double *daXBoundary; /**< Locations of grid boundaries in x = sin(lat) */
double *daYBoundary; /**< Locations of grid boundaries in y = R*lat */
// FLARE
int bFlare;
double dFlareYInt; /**< Flare function Y intercept /FFD linear coefficient*/
// double dFlareYIntErrorUpper; /**< Upper error of the Y intercept /FFD
// linear coefficient*/
// double dFlareYIntErrorLower; /**< Lower error of the Y intercept /FFD
// linear coefficient*/
double dFlareSlope; /**< Flare function slope /FFD angular coefficient*/
// double dFlareSlopeErrorUpper; /**< Upper error of slope /FFD angular
// coefficient*/
// double dFlareSlopeErrorLower; /**< Lower error of slope /FFD angular
// coefficient*/
double dFlareMinEnergy; /**< Flare minimum energy value to calculate the FFD*/
double dFlareMaxEnergy; /**< Flare maximum energy value to calculate the FFD*/
double dFlareFreq1; /**< First value of flare frequency range*/
double dFlareFreq2; /**< Second value of flare frequency range*/
double dFlareFreq3; /**< Third value of flare frequency range*/
double dFlareFreq4; /**< Fourth value of flare frequency range*/
double dFlareFreqMin; /**< Flare frequency of the flares with the lowest
energy*/
double dFlareFreqMid; /**< Flare frequency of the flares with the central
energy value in the energy range*/
double dFlareFreqMax; /**< Flare frequency of the flares with the highest
energy*/
double dFlareEnergy1; /**< First value of flare energy range*/
double dFlareEnergy2; /**< Second value of flare energy range*/
double dFlareEnergy3; /**< Third value of flare energy range*/
double dFlareEnergy4; /**< Fourth value of flare energy range*/
double dFlareEnergyMin; /**< Minimum value of flare energy in the energy
range*/
double dFlareEnergyMid; /**< Central value of flare energy in the energy
range*/
double dFlareEnergyMax; /**< Maximum value of flare energy in the energy
range*/
double dLXUVFlare; /**< XUV luminosity by flare*/
// double dLXUVFlareUpper; /**< Upper value of XUV luminosity by flare when
// the
// user include the slope and Y-intercept errors */
// double dLXUVFlareLower; /**< Lower value of XUV luminosity by flare when
// the
// user include the slope and Y-intercept errors */
double dLXUVTot; /**< XUV luminosity total, flare + stellar*/
double dLXUVFlareConst; /**< XUV luminosity given by the user*/
int iFlareFFD; /**< Flare mode*/
int iFlareBandPass; /**< Option to choose in which band pass the input energy
are*/
int iFlareSlopeUnits; /**< Mode to choose in which units the FFD slopes are*/
int iEnergyBin; /**< Number of energies consider between the minimum and
maximum energies to calculate the luminosity by flares*/
double *daEnergyERG;
double *daEnergyJOU;
double *daLogEner;
double *daEnerJOU;
double *daEnergyJOUXUV;
double *daEnergyERGXUV;
double *daLogEnerXUV;
double *daFFD;
double *daLXUVFlare;
// GALHABIT
int bGalHabit; /**< Use galhabit module */
double dPeriQ; /**< Pericenter distance */
int iDisrupt; /**< Secondary body has been disrupted */
int bGalacTides; /**< Enable galactic tides */
double dHostBinSemi; /**< Semi-major axis of host binary */
double dHostBinEcc; /**< Eccentricity of host binary */
double dHostBinInc; /**< Inclination of host binary */
double dHostBinArgP; /**< Arg pericenter of host binary */
double dHostBinLongA; /**< Long asc node of host binary */
double dHostBinMass1; /**< Mass of large host binary star */
int bHostBinary; /**< Model dynamics of inner host binary */
double
*daRelativeImpact; /**< Impact param of passing star relative to body */
double *daRelativeVel; /**< Velocity of passing star relative to body */
double dEccX; /**< X component of eccentricity vector */
double dEccY; /**< Y component of eccentricity vector */
double dEccZ; /**< Z component of eccentricity vector */
double dAngMX; /**< X component of orbital momentum vector */
double dAngMY; /**< Y component of orbital momentum vector */
double dAngMZ; /**< Z component of orbital momentum vector */
double dAngM; /**< Magnitude of orbital momentum vector */
double dEccXTmp; /**< Ecc X in the host binary reference plane */
double dEccYTmp; /**< Ecc Y in the host binary reference plane */
double dEccZTmp; /**< Ecc Z in the host binary reference plane */
double dAngMXTmp; /**< AngM X in the host binary reference plane */
double dAngMYTmp; /**< AngM Y in the host binary reference plane */
double dAngMZTmp; /**< AngM Z in the host binary reference plane */
double dArgPTmp; /**< Arg pericenter in the host binary reference plane */
double dLongATmp; /**< Long asc node in the host binary reference plane */
double dIncTmp; /**< Inclination in the host binary reference plane */
double dCosArgP; /**< Cosine of arg pericenter */
double dMinStellarApproach; /**< minimum allowed close approach of body to
host */
double dMassInterior; /**< Total mass of bodies interior to body */
int iBadImpulse; /**< Was there a bad impulse? */
double dMeanL; /**< Body's mean longitude */
// MAGMOC
/* HERE
* declare all variables used
*/
int bMagmOc; /**< Use magmoc model */
int bManSolid; /**< Mantle solidified */
int bAllFeOOxid; /**< All FeO in manlte oxidized to Fe2O3 */
int bLowPressSol; /**< Switch to low pressure treatment of solidus */
int bManStartSol; /**< Mantle starts to solidify */
int bCalcFugacity; /**< Need to calc oxygen fugacity */
int bPlanetDesiccated; /**< Atmosphere desiccated */
int bManQuasiSol; /**< Atmosphere desiccated & T_surf below 1000K */
int bMagmOcHaltSolid; /**< Mantle solidifed or atm desiccated */
int bMagmOcHaltDesicc; /**< Atm desiccated or escape stopped*/
int bEscapeStop; /**< Atmospheric escaped stopped */
int bCO2InAtmosphere; /**< Is CO2 present in the atmopshere? */
int iRadioHeatModel; /**< Which Radiogenic Heating model to use */
int iMagmOcAtmModel; /**< Which Atmopsheric Flux model to use */
int bOptManQuasiSol; /**< Solidify mantle inst. when melt frac = 0.4 at surf
*/
/* Primary variables */
double dPotTemp; /**< Potential Temp of the mantle [K] */
double dSurfTemp; /**< Surface Temp of the planet [K] */
double dSolidRadius; /**< Solidification radius of the mantle [m] */
double dWaterMassMOAtm; /**< Water mass in magma ocean and atmosphere [kg] */
double dWaterMassSol; /**< Water mass in the solidified mantle [kg] */
double dOxygenMassMOAtm; /**< Water mass in magma ocean and atmosphere [kg] */
double dOxygenMassSol; /**< Water mass in the solidified mantle [kg] */
double dHydrogenMassSpace; /**< Mass of hydrogen that is lost to space */
double dOxygenMassSpace; /**< Mass of oxygen that is lost to space */
double dCO2MassMOAtm; /**< Mass of CO2 in magma ocean and atmosphere [kg] */
double dCO2MassSol; /**< Mass of CO2 in solidified mantle [kg] */
/* Input variables */
double dCoreRadius; /**< Core radius of the planet [m] */
double dWaterMassAtm; /**< Water mass in the atmosphere [kg] */
double dManMeltDensity; /**< Density of the molten mantle [km/m^3] */
double dMassFracFeOIni; /**< Initial FeO mass fraction in the mantle */
double dWaterPartCoeff; /**< Water partition coefficient between melt and
solid */
double dDepthMO; /**< Initial depth of Magma Ocean [km] */
/* Other variables Thermal model */
double
dGravAccelSurf; /**< Graviational acceleration at the surface [m/s^2] */
double dSolidRadiusLocal; /**< Local variable for solidification radius of the
mantle [m] */
double dTransDepthSol; /**< Depth of transition from low to high pressure
solidus [Pa] */
double dPrefactorA; /**< Prefactor for linear solidus */
double dPrefactorB; /**< Prefactor for linear solidus */
double dMeltFraction; /**< Melt fraction of the mantle */
double dMeltFracSurf; /**< Melt fraction at the surface */
double dKinemViscos; /**< Kinematic viscosity of the mantle [m/s^2] */
double dFactorDerivative; /**< Factor to calculate the derivatives of Tpot and
Rsol */
double dManHeatFlux; /**< Mantle heat flux [W/m^2] */
double dRadioHeat; /**< Radiogenic heating rate GET FROM RADHEAT [W/kg] */
double dTidalHeat; /**< Tidal heating rate GET FROM EQTIDE [W/kg] */
double dNetFluxAtmo; /**< Net atmospheric flux OLR-ASR [W/m^2] */
double dAlbedo; /**< Albedo of the planet */
double dEffTempAtm; /**< Effective temperature of the planet's atmosphere */
/* Other variables Volatile model */
double dPressWaterAtm; /**< Water pressure in atmosphere [Pa] */