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

AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_decl (real &calday, real &eccen, real &mvelpp, real &lambm0, real &obliqr, real &delta, real &eccf)
 
AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_params (int &iyear_AD, real &eccen, real &obliq, real &mvelp, real &obliqr, real &lambm0, real &mvelpp)
 
AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_avg_cos_zenith (real &jday, real &lat, real &lon, real &declin, real &dt_avg)
 
AMREX_GPU_HOST AMREX_FORCE_INLINE 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.)
 

Typedef Documentation

◆ real

typedef double real

Function Documentation

◆ orbital_avg_cos_zenith()

AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_avg_cos_zenith ( real &  jday,
real &  lat,
real &  lon,
real &  declin,
real &  dt_avg 
)
461 {
462  // adjust latitude so that its tangent will be defined
463  real del;
464  if (lat == PIoTwo) {
465  del = lat - real(1.0e-05);
466  } else if (lat == -PIoTwo) {
467  del = lat + real(1.0e-05);
468  } else {
469  del = lat;
470  }
471 
472  // adjust declination so that its tangent will be defined
473  real phi;
474  if (declin == PIoTwo) {
475  phi = declin - real(1.0e-05);
476  } else if (declin == -PIoTwo) {
477  phi = declin + real(1.0e-05);
478  } else {
479  phi = declin;
480  }
481 
482  // define the cosine of the half-day length
483  // adjust for cases of all daylight or all night
484  real h;
485  real cos_h = - std::tan(del) * std::tan(phi);
486  if (cos_h <= -1.0) {
487  h = PI;
488  } else if (cos_h >= 1.0) {
489  h = 0.0;
490  } else {
491  h = std::acos(cos_h);
492  }
493 
494  // Define Local Time t and t + dt
495  // adjust t to be between -pi and pi
496  real t1 = (jday - int(jday)) * 2.0*PI + lon - PI;
497  if (t1 >= PI) {
498  t1 = t1 - 2.0*PI;
499  } else if (t1 < -PI) {
500  t1 = t1 + 2.0*PI;
501  }
502 
503  real dt = dt_avg / real(86400.0) * 2.0*PI;
504  real t2 = t1 + dt;
505 
506  // Compute Cosine Solar Zenith angle
507  // define terms needed in the cosine zenith angle equation
508  real aa = std::sin(lat) * std::sin(declin);
509  real bb = std::cos(lat) * std::cos(declin);
510 
511  // define the hour angle
512  // force it to be between -h and h
513  // consider the situation when the night period is too short
514  real tt1,tt2,tt3,tt4;
515  if ( (t2 >= PI) && (t1 <= PI) && (PI - h <= dt) ) {
516  tt2 = h;
517  tt1 = std::min(std::max(t1, -h), h);
518  tt4 = std::min(std::max(t2, 2.0*PI - h), 2.0*PI + h);
519  tt3 = 2.0*PI - h;
520  } else if ( (t2 >= -PI) && (t1 <= -PI) && (PI - h <= dt) ) {
521  tt2 = - 2.0*PI + h;
522  tt1 = std::min(std::max(t1, -2.0*PI - h), -2.0*PI + h);
523  tt4 = std::min(std::max(t2, -h), h);
524  tt3 = -h;
525  } else {
526  if (t2 > PI) {
527  tt2 = std::min(std::max(t2 - 2.0*PI, -h), h);
528  } else if (t2 < - PI) {
529  tt2 = std::min(std::max(t2 + 2.0*PI, -h), h);
530  } else {
531  tt2 = std::min(std::max(t2 , -h), h);
532  }
533 
534  if (t1 > PI) {
535  tt1 = std::min(std::max(t1 - 2.0*PI, -h), h);
536  } else if (t1 < -PI) {
537  tt1 = std::min(std::max(t1 + 2.0*PI, -h), h);
538  } else {
539  tt1 = std::min(std::max(t1 , -h), h);
540  }
541  tt4 = 0.0;
542  tt3 = 0.0;
543  }
544 
545  // perform a time integration to obtain cosz if desired
546  // output is valid over the period from t to t + dt
547  if ( (tt2 > tt1) || (tt4 > tt3) ) {
548  return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
549  (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
550  } else {
551  return 0.0;
552  }
553 }
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()

AMREX_GPU_HOST AMREX_FORCE_INLINE 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. 
)
566 {
567  // Constant zenith angle is true
568  if ( constant_zenith_angle_deg >= 0. ) {
569  return std::cos( constant_zenith_angle_deg * PI/180. );
570  }
571 
572  // Uniform angle is true
573  if ( uniform_angle >= 0.) {
574  return std::cos(uniform_angle);
575  }
576 
577  // Perform the calculation of Orbital_Cos_Zenith
578  bool use_dt_avg = false;
579  if (dt_avg > 0.) {
580  use_dt_avg = true;
581  }
582  // If dt for the average cosz is specified, then call the shr_orb_avg_cosz
583  if (use_dt_avg) {
584  return orbital_avg_cos_zenith(jday, lat, lon, declin, dt_avg);
585  } else {
586  return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
587  std::cos((jday-floor(jday))*real(2.0)*PI + lon);
588  }
589 }
orbital_avg_cos_zenith
AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_avg_cos_zenith(real &jday, real &lat, real &lon, real &declin, real &dt_avg)
Definition: ERF_OrbCosZenith.H:456

Referenced by Radiation::run_impl().

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

◆ orbital_decl()

AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_decl ( real &  calday,
real &  eccen,
real &  mvelpp,
real &  lambm0,
real &  obliqr,
real &  delta,
real &  eccf 
)
22 {
23 
24  real lambm; // Lambda m, mean long of perihelion (rad)
25  real lmm; // Intermediate argument involving lambm
26  real lamb; // Lambda, the earths long of perihelion
27  real invrho; // Inverse normalized sun/earth distance
28  real sinl; // Sine of lmm
29  static constexpr real dayspy = real(365.0); // Day per year
30  static constexpr real ve = real( 80.5); // Calday of vernal equinox
31 
32  /*
33  Compute eccentricity factor and solar declination using
34  day value where a round day (such as 213.0) refers to 0z at
35  Greenwich longitude.
36 
37  Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
38  Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
39  35:2362-2367.
40 
41  To get the earths true longitude (position in orbit; lambda in Berger
42  1978) which is necessary to find the eccentricity factor and declination,
43  must first calculate the mean longitude (lambda m in Berger 1978) at
44  the present day. This is done by adding to lambm0 (the mean longitude
45  at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
46  an increment (delta lambda m in Berger 1978) that is the number of
47  days past or before (a negative increment) the vernal equinox divided by
48  the days in a model year times the 2*pi radians in a complete orbit.
49  */
50 
51  lambm = lambm0 + (calday - ve)*2.0*PI/dayspy;
52  lmm = lambm - mvelpp;
53 
54  // The earths true longitude, in radians, is then found from
55  // the formula in Berger 1978:
56 
57  sinl = std::sin(lmm);
58  lamb = lambm + eccen*( 2.0*sinl + eccen*( 1.25*sin(2.0*lmm)
59  + eccen*((13.0/12.0)*std::sin(3.0*lmm) - 0.25*sinl)));
60 
61  /*
62  Using the obliquity, eccentricity, moving vernal equinox longitude of
63  perihelion (plus), and earths true longitude, the declination (delta)
64  and the normalized earth/sun distance (rho in Berger 1978; actually inverse
65  rho will be used), and thus the eccentricity factor (eccf), can be
66  calculated from formulas given in Berger 1978.
67  */
68 
69  invrho = (1.0 + eccen*std::cos(lamb - mvelpp)) / (1.0 - eccen*eccen);
70 
71  // Set solar declination and eccentricity factor
72 
73  delta = std::asin(std::sin(obliqr)*std::sin(lamb));
74  eccf = invrho*invrho;
75 }

Referenced by Radiation::run_impl().

Here is the caller graph for this function:

◆ orbital_params()

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