ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Typedefs | Functions
ERF_OrbCosZenith.H File Reference
#include <vector>
#include <cmath>
#include <ERF_Constants.H>
Include dependency graph for ERF_OrbCosZenith.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Typedefs

typedef double real
 

Functions

real orbital_cos_zenith (real &jday, real &lat, real &lon, real &declin, real dt_avg=-1., real uniform_angle=-1., real constant_zenith_angle_deg=-1.)
 
real orbital_avg_cos_zenith (real &jday, real &lat, real &lon, real &declin, real &dt_avg)
 
void orbital_params (int &iyear_AD, real &eccen, real &obliq, real &mvelp, real &obliqr, real &lambm0, real &mvelpp)
 
void orbital_decl (real &calday, real &eccen, real &mvelpp, real &lambm0, real &obliqr, real &delta, real &eccf)
 

Typedef Documentation

◆ real

typedef double real

Function Documentation

◆ orbital_avg_cos_zenith()

real orbital_avg_cos_zenith ( real &  jday,
real &  lat,
real &  lon,
real &  declin,
real &  dt_avg 
)
43 {
44  // adjust latitude so that its tangent will be defined
45  real del;
46  if (lat == PIoTwo) {
47  del = lat - real(1.0e-05);
48  } else if (lat == -PIoTwo) {
49  del = lat + real(1.0e-05);
50  } else {
51  del = lat;
52  }
53 
54  // adjust declination so that its tangent will be defined
55  real phi;
56  if (declin == PIoTwo) {
57  phi = declin - real(1.0e-05);
58  } else if (declin == -PIoTwo) {
59  phi = declin + real(1.0e-05);
60  } else {
61  phi = declin;
62  }
63 
64  // define the cosine of the half-day length
65  // adjust for cases of all daylight or all night
66  real h;
67  real cos_h = - std::tan(del) * std::tan(phi);
68  if (cos_h <= -1.0) {
69  h = PI;
70  } else if (cos_h >= 1.0) {
71  h = 0.0;
72  } else {
73  h = std::acos(cos_h);
74  }
75 
76  // Define Local Time t and t + dt
77  // adjust t to be between -pi and pi
78  real t1 = (jday - int(jday)) * 2.0*PI + lon - PI;
79  if (t1 >= PI) {
80  t1 = t1 - 2.0*PI;
81  } else if (t1 < -PI) {
82  t1 = t1 + 2.0*PI;
83  }
84 
85  real dt = dt_avg / real(86400.0) * 2.0*PI;
86  real t2 = t1 + dt;
87 
88  // Compute Cosine Solar Zenith angle
89  // define terms needed in the cosine zenith angle equation
90  real aa = std::sin(lat) * std::sin(declin);
91  real bb = std::cos(lat) * std::cos(declin);
92 
93  // define the hour angle
94  // force it to be between -h and h
95  // consider the situation when the night period is too short
96  real tt1,tt2,tt3,tt4;
97  if ( (t2 >= PI) && (t1 <= PI) && (PI - h <= dt) ) {
98  tt2 = h;
99  tt1 = std::min(std::max(t1, -h), h);
100  tt4 = std::min(std::max(t2, 2.0*PI - h), 2.0*PI + h);
101  tt3 = 2.0*PI - h;
102  } else if ( (t2 >= -PI) && (t1 <= -PI) && (PI - h <= dt) ) {
103  tt2 = - 2.0*PI + h;
104  tt1 = std::min(std::max(t1, -2.0*PI - h), -2.0*PI + h);
105  tt4 = std::min(std::max(t2, -h), h);
106  tt3 = -h;
107  } else {
108  if (t2 > PI) {
109  tt2 = std::min(std::max(t2 - 2.0*PI, -h), h);
110  } else if (t2 < - PI) {
111  tt2 = std::min(std::max(t2 + 2.0*PI, -h), h);
112  } else {
113  tt2 = std::min(std::max(t2 , -h), h);
114  }
115 
116  if (t1 > PI) {
117  tt1 = std::min(std::max(t1 - 2.0*PI, -h), h);
118  } else if (t1 < -PI) {
119  tt1 = std::min(std::max(t1 + 2.0*PI, -h), h);
120  } else {
121  tt1 = std::min(std::max(t1 , -h), h);
122  }
123  tt4 = 0.0;
124  tt3 = 0.0;
125  }
126 
127  // perform a time integration to obtain cosz if desired
128  // output is valid over the period from t to t + dt
129  if ( (tt2 > tt1) || (tt4 > tt3) ) {
130  return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
131  (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
132  } else {
133  return 0.0;
134  }
135 }
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
real
double real
Definition: ERF_OrbCosZenith.H:9

Referenced by orbital_cos_zenith().

Here is the caller graph for this function:

◆ orbital_cos_zenith()

real orbital_cos_zenith ( real &  jday,
real &  lat,
real &  lon,
real &  declin,
real  dt_avg = -1.,
real  uniform_angle = -1.,
real  constant_zenith_angle_deg = -1. 
)
11 {
12  // Constant zenith angle is true
13  if ( constant_zenith_angle_deg >= 0. ) {
14  return std::cos( constant_zenith_angle_deg * PI/180. );
15  }
16 
17  // Uniform angle is true
18  if ( uniform_angle >= 0.) {
19  return std::cos(uniform_angle);
20  }
21 
22  // Perform the calculation of Orbital_Cos_Zenith
23  bool use_dt_avg = false;
24  if (dt_avg > 0.) {
25  use_dt_avg = true;
26  }
27  // If dt for the average cosz is specified, then call the shr_orb_avg_cosz
28  if (use_dt_avg) {
29  return orbital_avg_cos_zenith(jday, lat, lon, declin, dt_avg);
30  } else {
31  return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
32  std::cos((jday-floor(jday))*real(2.0)*PI + lon);
33  }
34 }
orbital_avg_cos_zenith
real orbital_avg_cos_zenith(real &jday, real &lat, real &lon, real &declin, real &dt_avg)
Definition: ERF_OrbCosZenith.cpp:38

Referenced by Radiation::run_impl().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ orbital_decl()

void orbital_decl ( real &  calday,
real &  eccen,
real &  mvelpp,
real &  lambm0,
real &  obliqr,
real &  delta,
real &  eccf 
)
519 {
520 
521  real lambm; // Lambda m, mean long of perihelion (rad)
522  real lmm; // Intermediate argument involving lambm
523  real lamb; // Lambda, the earths long of perihelion
524  real invrho; // Inverse normalized sun/earth distance
525  real sinl; // Sine of lmm
526  static constexpr real dayspy = real(365.0); // Day per year
527  static constexpr real ve = real( 80.5); // Calday of vernal equinox
528 
529  /*
530  Compute eccentricity factor and solar declination using
531  day value where a round day (such as 213.0) refers to 0z at
532  Greenwich longitude.
533 
534  Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
535  Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
536  35:2362-2367.
537 
538  To get the earths true longitude (position in orbit; lambda in Berger
539  1978) which is necessary to find the eccentricity factor and declination,
540  must first calculate the mean longitude (lambda m in Berger 1978) at
541  the present day. This is done by adding to lambm0 (the mean longitude
542  at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
543  an increment (delta lambda m in Berger 1978) that is the number of
544  days past or before (a negative increment) the vernal equinox divided by
545  the days in a model year times the 2*pi radians in a complete orbit.
546  */
547 
548  lambm = lambm0 + (calday - ve)*2.0*PI/dayspy;
549  lmm = lambm - mvelpp;
550 
551  // The earths true longitude, in radians, is then found from
552  // the formula in Berger 1978:
553 
554  sinl = std::sin(lmm);
555  lamb = lambm + eccen*( 2.0*sinl + eccen*( 1.25*sin(2.0*lmm)
556  + eccen*((13.0/12.0)*std::sin(3.0*lmm) - 0.25*sinl)));
557 
558  /*
559  Using the obliquity, eccentricity, moving vernal equinox longitude of
560  perihelion (plus), and earths true longitude, the declination (delta)
561  and the normalized earth/sun distance (rho in Berger 1978; actually inverse
562  rho will be used), and thus the eccentricity factor (eccf), can be
563  calculated from formulas given in Berger 1978.
564  */
565 
566  invrho = (1.0 + eccen*std::cos(lamb - mvelpp)) / (1.0 - eccen*eccen);
567 
568  // Set solar declination and eccentricity factor
569 
570  delta = std::asin(std::sin(obliqr)*std::sin(lamb));
571  eccf = invrho*invrho;
572 }

Referenced by Radiation::run_impl().

Here is the caller graph for this function:

◆ orbital_params()

void orbital_params ( int &  iyear_AD,
real &  eccen,
real &  obliq,
real &  mvelp,
real &  obliqr,
real &  lambm0,
real &  mvelpp 
)
146 {
147  /*
148  !-------------------------------------------------------------------------------
149  !
150  ! Calculate earths orbital parameters using Dave Threshers formula which
151  ! came from Berger, Andre. 1978 "A Simple Algorithm to Compute Long-Term
152  ! Variations of Daily Insolation". Contribution 18, Institute of Astronomy
153  ! and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium
154  !
155  !-------------------------------------------------------------------------------
156  */
157 
158  int poblen = 47; // # of elements in series wrt obliquity
159  int pecclen = 19; // # of elements in series wrt eccentricity
160  int pmvelen = 78; // # of elements in series wrt vernal equinox
161  static constexpr real psecdeg = real(1.0)/real(3600.0); // arc sec to deg conversion
162  static constexpr real degrad = PI/real(180.); // degree to radian conversion factor
163 
164  // Cosine series data for computation of obliquity: amplitude (arc seconds),
165  // rate (arc seconds/year), phase (degrees).
166  // amplitudes for obliquity cos series
167  std::vector<real> obamp = {-2462.2214466, -857.3232075, -629.3231835,
168  -414.2804924, -311.7632587, 308.9408604,
169  -162.5533601, -116.1077911, 101.1189923,
170  -67.6856209, 24.9079067, 22.5811241,
171  -21.1648355, -15.6549876, 15.3936813,
172  14.6660938, -11.7273029, 10.2742696,
173  6.4914588, 5.8539148, -5.4872205,
174  -5.4290191, 5.1609570, 5.0786314,
175  -4.0735782, 3.7227167, 3.3971932,
176  -2.8347004, -2.6550721, -2.5717867,
177  -2.4712188, 2.4625410, 2.2464112,
178  -2.0755511, -1.9713669, -1.8813061,
179  -1.8468785, 1.8186742, 1.7601888,
180  -1.5428851, 1.4738838, -1.4593669,
181  1.4192259, -1.1818980, 1.1756474,
182  -1.1316126, 1.0896928};
183  // rates for obliquity cosine series
184  std::vector<real> obrate = {31.609974, 32.620504, 24.172203,
185  31.983787, 44.828336, 30.973257,
186  43.668246, 32.246691, 30.599444,
187  42.681324, 43.836462, 47.439436,
188  63.219948, 64.230478, 1.010530,
189  7.437771, 55.782177, 0.373813,
190  13.218362, 62.583231, 63.593761,
191  76.438310, 45.815258, 8.448301,
192  56.792707, 49.747842, 12.058272,
193  75.278220, 65.241008, 64.604291,
194  1.647247, 7.811584, 12.207832,
195  63.856665, 56.155990, 77.448840,
196  6.801054, 62.209418, 20.656133,
197  48.344406, 55.145460, 69.000539,
198  11.071350, 74.291298, 11.047742,
199  0.636717, 12.844549};
200 
201  // phases for obliquity cosine series
202  std::vector<real> obphas = {251.9025, 280.8325, 128.3057,
203  292.7252, 15.3747, 263.7951,
204  308.4258, 240.0099, 222.9725,
205  268.7809, 316.7998, 319.6024,
206  143.8050, 172.7351, 28.9300,
207  123.5968, 20.2082, 40.8226,
208  123.4722, 155.6977, 184.6277,
209  267.2772, 55.0196, 152.5268,
210  49.1382, 204.6609, 56.5233,
211  200.3284, 201.6651, 213.5577,
212  17.0374, 164.4194, 94.5422,
213  131.9124, 61.0309, 296.2073,
214  135.4894, 114.8750, 247.0691,
215  256.6114, 32.1008, 143.6804,
216  16.8784, 160.6835, 27.5932,
217  348.1074, 82.6496};
218 
219  // Cosine/sine series data for computation of eccentricity and fixed vernal
220  // equinox longitude of perihelion (fvelp): amplitude,
221  // rate (arc seconds/year), phase (degrees).
222 
223  // ampl for eccen/fvelp cos/sin series
224  std::vector<real> ecamp = { 0.01860798, 0.01627522, -0.01300660,
225  0.00988829, -0.00336700, 0.00333077,
226  -0.00235400, 0.00140015, 0.00100700,
227  0.00085700, 0.00064990, 0.00059900,
228  0.00037800, -0.00033700, 0.00027600,
229  0.00018200, -0.00017400, -0.00012400,
230  0.00001250};
231 
232  // rates for eccen/fvelp cos/sin series
233  std::vector<real> ecrate = { 4.2072050, 7.3460910, 17.8572630,
234  17.2205460, 16.8467330, 5.1990790,
235  18.2310760, 26.2167580, 6.3591690,
236  16.2100160, 3.0651810, 16.5838290,
237  18.4939800, 6.1909530, 18.8677930,
238  17.4255670, 6.1860010, 18.4174410,
239  0.6678630};
240 
241  // phases for eccen/fvelp cos/sin series
242  std::vector<real> ecphas = { 28.620089, 193.788772, 308.307024,
243  320.199637, 279.376984, 87.195000,
244  349.129677, 128.443387, 154.143880,
245  291.269597, 114.860583, 332.092251,
246  296.414411, 145.769910, 337.237063,
247  152.092288, 126.839891, 210.667199,
248  72.108838};
249 
250  // Sine series data for computation of moving vernal equinox longitude of
251  // perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees).
252 
253  // amplitudes for mvelp sine series
254  std::vector<real> mvamp = { 7391.0225890, 2555.1526947, 2022.7629188,
255  -1973.6517951, 1240.2321818, 953.8679112,
256  -931.7537108, 872.3795383, 606.3544732,
257  -496.0274038, 456.9608039, 346.9462320,
258  -305.8412902, 249.6173246, -199.1027200,
259  191.0560889, -175.2936572, 165.9068833,
260  161.1285917, 139.7878093, -133.5228399,
261  117.0673811, 104.6907281, 95.3227476,
262  86.7824524, 86.0857729, 70.5893698,
263  -69.9719343, -62.5817473, 61.5450059,
264  -57.9364011, 57.1899832, -57.0236109,
265  -54.2119253, 53.2834147, 52.1223575,
266  -49.0059908, -48.3118757, -45.4191685,
267  -42.2357920, -34.7971099, 34.4623613,
268  -33.8356643, 33.6689362, -31.2521586,
269  -30.8798701, 28.4640769, -27.1960802,
270  27.0860736, -26.3437456, 24.7253740,
271  24.6732126, 24.4272733, 24.0127327,
272  21.7150294, -21.5375347, 18.1148363,
273  -16.9603104, -16.1765215, 15.5567653,
274  15.4846529, 15.2150632, 14.5047426,
275  -14.3873316, 13.1351419, 12.8776311,
276  11.9867234, 11.9385578, 11.7030822,
277  11.6018181, -11.2617293, -10.4664199,
278  10.4333970, -10.2377466, 10.1934446,
279  -10.1280191, 10.0289441, -10.0034259};
280 
281  // rates for mvelp sine series
282  std::vector<real> mvrate = {31.609974, 32.620504, 24.172203,
283  0.636717, 31.983787, 3.138886,
284  30.973257, 44.828336, 0.991874,
285  0.373813, 43.668246, 32.246691,
286  30.599444, 2.147012, 10.511172,
287  42.681324, 13.650058, 0.986922,
288  9.874455, 13.013341, 0.262904,
289  0.004952, 1.142024, 63.219948,
290  0.205021, 2.151964, 64.230478,
291  43.836462, 47.439436, 1.384343,
292  7.437771, 18.829299, 9.500642,
293  0.431696, 1.160090, 55.782177,
294  12.639528, 1.155138, 0.168216,
295  1.647247, 10.884985, 5.610937,
296  12.658184, 1.010530, 1.983748,
297  14.023871, 0.560178, 1.273434,
298  12.021467, 62.583231, 63.593761,
299  76.438310, 4.280910, 13.218362,
300  17.818769, 8.359495, 56.792707,
301  8.448301, 1.978796, 8.863925,
302  0.186365, 8.996212, 6.771027,
303  45.815258, 12.002811, 75.278220,
304  65.241008, 18.870667, 22.009553,
305  64.604291, 11.498094, 0.578834,
306  9.237738, 49.747842, 2.147012,
307  1.196895, 2.133898, 0.173168};
308 
309  // phases for mvelp sine series
310  std::vector<real> mvphas = {251.9025, 280.8325, 128.3057,
311  348.1074, 292.7252, 165.1686,
312  263.7951, 15.3747, 58.5749,
313  40.8226, 308.4258, 240.0099,
314  222.9725, 106.5937, 114.5182,
315  268.7809, 279.6869, 39.6448,
316  126.4108, 291.5795, 307.2848,
317  18.9300, 273.7596, 143.8050,
318  191.8927, 125.5237, 172.7351,
319  316.7998, 319.6024, 69.7526,
320  123.5968, 217.6432, 85.5882,
321  156.2147, 66.9489, 20.2082,
322  250.7568, 48.0188, 8.3739,
323  17.0374, 155.3409, 94.1709,
324  221.1120, 28.9300, 117.1498,
325  320.5095, 262.3602, 336.2148,
326  233.0046, 155.6977, 184.6277,
327  267.2772, 78.9281, 123.4722,
328  188.7132, 180.1364, 49.1382,
329  152.5268, 98.2198, 97.4808,
330  221.5376, 168.2438, 161.1199,
331  55.0196, 262.6495, 200.3284,
332  201.6651, 294.6547, 99.8233,
333  213.5577, 154.1631, 232.7153,
334  138.3034, 204.6609, 106.5938,
335  250.4676, 332.3345, 27.3039};
336 
337  //---------------------------Local variables----------------------------------
338  real obsum; // Obliquity series summation
339  real cossum; // Cos series summation for eccentricity/fvelp
340  real sinsum; // Sin series summation for eccentricity/fvelp
341  real fvelp; // Fixed vernal equinox long of perihelion
342  real mvsum; // mvelp series summation
343  real beta; // Intermediate argument for lambm0
344  real years; // Years to time of interest ( pos <=> future)
345  real eccen2; // eccentricity squared
346  real eccen3; // eccentricity cubed
347  real yb4_1950AD; // number of years before 1950 AD
348 
349  // Check for flag to use input orbit parameters
350  if ( iyear_AD == ORB_UNDEF_INT ) {
351 
352  // The AMIP II settings (for a 1995 orbit) are adopted
353  obliq = 23.4441;
354  eccen = 0.016715;
355  mvelp = 102.7;
356  eccen2 = eccen*eccen;
357  eccen3 = eccen2*eccen;
358 
359  } else { // Otherwise calculate based on years before present
360 
361  /*
362  The following calculates the earths obliquity, orbital eccentricity
363  (and various powers of it) and vernal equinox mean longitude of
364  perihelion for years in the past (future = negative of years past),
365  using constants (see parameter section) given in the program of:
366 
367  Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations
368  of Daily Insolation. Contribution 18, Institute of Astronomy and
369  Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
370 
371  and formulas given in the paper (where less precise constants are also
372  given):
373 
374  Berger, Andre. 1978. Long-Term Variations of Daily Insolation and
375  Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367
376 
377  The algorithm is valid only to 1,000,000 years past or hence.
378  For a solution valid to 5-10 million years past see the above author.
379  Algorithm below is better for years closer to present than is the
380  5-10 million year solution.
381 
382  Years to time of interest must be negative of years before present
383  (1950) in formulas that follow.
384  */
385 
386  yb4_1950AD = real(1950.0) - real(iyear_AD);
387  years = - yb4_1950AD;
388 
389  /*
390  In the summations below, cosine or sine arguments, which end up in
391  degrees, must be converted to radians via multiplication by degrad.
392 
393  Summation of cosine series for obliquity (epsilon in Berger 1978) in
394  degrees. Convert the amplitudes and rates, which are in arc secs, into
395  degrees via multiplication by psecdeg (arc seconds to degrees conversion
396  factor). For obliq, first term is Berger 1978 epsilon star; second
397  term is series summation in degrees.
398  */
399 
400  obsum = 0.0;
401  for (int i(0); i<poblen; ++i) {
402  obsum = obsum + obamp[i]*psecdeg*std::cos( (obrate[i]*psecdeg*years + obphas[i]) * degrad );
403  }
404  obliq = real(23.320556) + obsum;
405 
406  /*
407  Summation of cosine and sine series for computation of eccentricity
408  (eccen; e in Berger 1978) and fixed vernal equinox longitude of
409  perihelion (fvelp; pi in Berger 1978), which is used for computation
410  of moving vernal equinox longitude of perihelion. Convert the rates,
411  which are in arc seconds, into degrees via multiplication by psecdeg.
412  */
413 
414  cossum = 0.0;
415  for (int i(0); i<pecclen; ++i) {
416  cossum = cossum + ecamp[i]*std::cos( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
417  }
418 
419  sinsum = 0.0;
420  for (int i(0); i<pecclen; ++i) {
421  sinsum = sinsum + ecamp[i]*std::sin( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
422  }
423 
424  // Use summations to calculate eccentricity
425 
426  eccen2 = cossum*cossum + sinsum*sinsum;
427  eccen = std::sqrt(eccen2);
428  eccen3 = eccen2*eccen;
429 
430  // A series of cases for fvelp, which is in radians.
431  if (std::fabs(cossum) <= 1.0e-8) {
432  if (sinsum == 0.0) {
433  fvelp = 0.0;
434  } else if (sinsum <= 0.0) {
435  fvelp = 1.5*PI;
436  } else if (sinsum > 0.0) {
437  fvelp = 0.5*PI;
438  }
439  } else if (cossum <= 0.0) {
440  fvelp = std::atan(sinsum/cossum) + PI;
441  } else {
442  if (sinsum < 0.0) {
443  fvelp = std::atan(sinsum/cossum) + 2.0*PI;
444  } else {
445  fvelp = std::atan(sinsum/cossum);
446  }
447  }
448 
449  /*
450  Summation of sin series for computation of moving vernal equinox long
451  of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp,
452  first term is fvelp in degrees; second term is Berger 1978 psi bar
453  times years and in degrees; third term is Berger 1978 zeta; fourth
454  term is series summation in degrees. Convert the amplitudes and rates,
455  which are in arc seconds, into degrees via multiplication by psecdeg.
456  Series summation plus second and third terms constitute Berger 1978
457  psi, which is the general precession.
458  */
459  mvsum = 0.0;
460  for (int i(0); i<pmvelen; ++i) {
461  mvsum = mvsum + mvamp[i]*psecdeg*std::sin( (mvrate[i]*psecdeg*years + mvphas[i]) * degrad);
462  }
463  mvelp = fvelp/degrad + real(50.439273)*psecdeg*years + real(3.392506) + mvsum;
464 
465  // Cases to make sure mvelp is between 0 and 360.
466  do {
467  mvelp += 360.0;
468  } while (mvelp < 0.0);
469 
470  do {
471  mvelp -= 360.0;
472  } while (mvelp >= 360.0);
473 
474  } // end of test on whether to calculate or use input orbital params
475 
476  // Orbit needs the obliquity in radians
477 
478  obliqr = obliq*degrad;
479 
480  /*
481  180 degrees must be added to mvelp since observations are made from the
482  earth and the sun is considered (wrongly for the algorithm) to go around
483  the earth. For a more graphic explanation see Appendix B in:
484 
485  A. Berger, M. Loutre and C. Tricot. 1993. Insolation and Earth Orbital
486  Periods. J. of Geophysical Research 98:10,341-10,362.
487 
488  Additionally, orbit will need this value in radians. So mvelp becomes
489  mvelpp (mvelp plus pi)
490  */
491 
492  mvelpp = (mvelp + 180.0)*degrad;
493 
494  // Set up an argument used several times in lambm0 calculation ahead.
495 
496  beta = std::sqrt(1.0 - eccen2);
497 
498  /*
499  The mean longitude at the vernal equinox (lambda m nought in Berger
500  1978; in radians) is calculated from the following formula given in
501  Berger 1978. At the vernal equinox the true longitude (lambda in Berger
502  1978) is 0.
503  */
504 
505  lambm0 = 2.0*( (.5*eccen + .125*eccen3)*(1.0 + beta)*std::sin(mvelpp)
506  - .250*eccen2*(.5 + beta)*std::sin(2.0*mvelpp)
507  + .125*eccen3*(1./3. + beta)*std::sin(3.0*mvelpp) );
508 }
ORB_UNDEF_INT
static constexpr int ORB_UNDEF_INT
Definition: ERF_Constants.H:96
module_model_constants::degrad
real(c_double), parameter degrad
Definition: ERF_module_model_constants.F90:75

Referenced by Radiation::run_impl().

Here is the caller graph for this function: