ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_OrbCosZenith.H
Go to the documentation of this file.
1 #ifndef ERF_ORB_COS_ZENITH_H
2 #define ERF_ORB_COS_ZENITH_H
3 
4 #include <vector>
5 #include <cmath>
6 
7 #include <ERF_Constants.H>
8 
9 typedef double real;
10 
11 
12 AMREX_GPU_HOST
13 AMREX_FORCE_INLINE
14 void
15 orbital_decl (real& calday,
16  real& eccen,
17  real& mvelpp,
18  real& lambm0,
19  real& obliqr,
20  real& delta,
21  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(amrex::Real(365.0)); // Day per year
30  static constexpr real ve = real( amrex::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 amrex::Real(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-amrex::Real(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)*two*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*( two*sinl + eccen*( amrex::Real(1.25)*sin(two*lmm)
59  + eccen*((amrex::Real(13.0)/amrex::Real(12.0))*std::sin(three*lmm) - fourth*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 amrex::Real(1978.)
67  */
68 
69  invrho = (one + eccen*std::cos(lamb - mvelpp)) / (one - 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 }
76 
77 
78 AMREX_GPU_HOST
79 AMREX_FORCE_INLINE
80 void
81 orbital_params (int& iyear_AD,
82  real& eccen,
83  real& obliq,
84  real& mvelp,
85  real& obliqr,
86  real& lambm0,
87  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(one)/real(amrex::Real(3600.0)); // arc sec to deg conversion
104  static constexpr real degrad = PI/real(amrex::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 = {-amrex::Real(2462.2214466), -amrex::Real(857.3232075), -amrex::Real(629.3231835),
110  -amrex::Real(414.2804924), -amrex::Real(311.7632587), amrex::Real(308.9408604),
111  -amrex::Real(162.5533601), -amrex::Real(116.1077911), amrex::Real(101.1189923),
112  -amrex::Real(67.6856209), amrex::Real(24.9079067), amrex::Real(22.5811241),
113  -amrex::Real(21.1648355), -amrex::Real(15.6549876), amrex::Real(15.3936813),
114  amrex::Real(14.6660938), -amrex::Real(11.7273029), amrex::Real(10.2742696),
115  amrex::Real(6.4914588), amrex::Real(5.8539148), -amrex::Real(5.4872205),
116  -amrex::Real(5.4290191), amrex::Real(5.1609570), amrex::Real(5.0786314),
117  -amrex::Real(4.0735782), amrex::Real(3.7227167), amrex::Real(3.3971932),
118  -amrex::Real(2.8347004), -amrex::Real(2.6550721), -amrex::Real(2.5717867),
119  -amrex::Real(2.4712188), amrex::Real(2.4625410), amrex::Real(2.2464112),
120  -amrex::Real(2.0755511), -amrex::Real(1.9713669), -amrex::Real(1.8813061),
121  -amrex::Real(1.8468785), amrex::Real(1.8186742), amrex::Real(1.7601888),
122  -amrex::Real(1.5428851), amrex::Real(1.4738838), -amrex::Real(1.4593669),
123  amrex::Real(1.4192259), -amrex::Real(1.1818980), amrex::Real(1.1756474),
124  -amrex::Real(1.1316126), amrex::Real(1.0896928)};
125  // rates for obliquity cosine series
126  std::vector<real> obrate = {amrex::Real(31.609974), amrex::Real(32.620504), amrex::Real(24.172203),
127  amrex::Real(31.983787), amrex::Real(44.828336), amrex::Real(30.973257),
128  amrex::Real(43.668246), amrex::Real(32.246691), amrex::Real(30.599444),
129  amrex::Real(42.681324), amrex::Real(43.836462), amrex::Real(47.439436),
130  amrex::Real(63.219948), amrex::Real(64.230478), amrex::Real(1.010530),
131  amrex::Real(7.437771), amrex::Real(55.782177), amrex::Real(0.373813),
132  amrex::Real(13.218362), amrex::Real(62.583231), amrex::Real(63.593761),
133  amrex::Real(76.438310), amrex::Real(45.815258), amrex::Real(8.448301),
134  amrex::Real(56.792707), amrex::Real(49.747842), amrex::Real(12.058272),
135  amrex::Real(75.278220), amrex::Real(65.241008), amrex::Real(64.604291),
136  amrex::Real(1.647247), amrex::Real(7.811584), amrex::Real(12.207832),
137  amrex::Real(63.856665), amrex::Real(56.155990), amrex::Real(77.448840),
138  amrex::Real(6.801054), amrex::Real(62.209418), amrex::Real(20.656133),
139  amrex::Real(48.344406), amrex::Real(55.145460), amrex::Real(69.000539),
140  amrex::Real(11.071350), amrex::Real(74.291298), amrex::Real(11.047742),
141  amrex::Real(0.636717), amrex::Real(12.844549)};
142 
143  // phases for obliquity cosine series
144  std::vector<real> obphas = {amrex::Real(251.9025), amrex::Real(280.8325), amrex::Real(128.3057),
145  amrex::Real(292.7252), amrex::Real(15.3747), amrex::Real(263.7951),
146  amrex::Real(308.4258), amrex::Real(240.0099), amrex::Real(222.9725),
147  amrex::Real(268.7809), amrex::Real(316.7998), amrex::Real(319.6024),
148  amrex::Real(143.8050), amrex::Real(172.7351), amrex::Real(28.9300),
149  amrex::Real(123.5968), amrex::Real(20.2082), amrex::Real(40.8226),
150  amrex::Real(123.4722), amrex::Real(155.6977), amrex::Real(184.6277),
151  amrex::Real(267.2772), amrex::Real(55.0196), amrex::Real(152.5268),
152  amrex::Real(49.1382), amrex::Real(204.6609), amrex::Real(56.5233),
153  amrex::Real(200.3284), amrex::Real(201.6651), amrex::Real(213.5577),
154  amrex::Real(17.0374), amrex::Real(164.4194), amrex::Real(94.5422),
155  amrex::Real(131.9124), amrex::Real(61.0309), amrex::Real(296.2073),
156  amrex::Real(135.4894), amrex::Real(114.8750), amrex::Real(247.0691),
157  amrex::Real(256.6114), amrex::Real(32.1008), amrex::Real(143.6804),
158  amrex::Real(16.8784), amrex::Real(160.6835), amrex::Real(27.5932),
159  amrex::Real(348.1074), amrex::Real(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 = { amrex::Real(0.01860798), amrex::Real(0.01627522), -amrex::Real(0.01300660),
167  amrex::Real(0.00988829), -amrex::Real(0.00336700), amrex::Real(0.00333077),
168  -amrex::Real(0.00235400), amrex::Real(0.00140015), amrex::Real(0.00100700),
169  amrex::Real(0.00085700), amrex::Real(0.00064990), amrex::Real(0.00059900),
170  amrex::Real(0.00037800), -amrex::Real(0.00033700), amrex::Real(0.00027600),
171  amrex::Real(0.00018200), -amrex::Real(0.00017400), -amrex::Real(0.00012400),
172  amrex::Real(0.00001250)};
173 
174  // rates for eccen/fvelp cos/sin series
175  std::vector<real> ecrate = { amrex::Real(4.2072050), amrex::Real(7.3460910), amrex::Real(17.8572630),
176  amrex::Real(17.2205460), amrex::Real(16.8467330), amrex::Real(5.1990790),
177  amrex::Real(18.2310760), amrex::Real(26.2167580), amrex::Real(6.3591690),
178  amrex::Real(16.2100160), amrex::Real(3.0651810), amrex::Real(16.5838290),
179  amrex::Real(18.4939800), amrex::Real(6.1909530), amrex::Real(18.8677930),
180  amrex::Real(17.4255670), amrex::Real(6.1860010), amrex::Real(18.4174410),
181  amrex::Real(0.6678630)};
182 
183  // phases for eccen/fvelp cos/sin series
184  std::vector<real> ecphas = { amrex::Real(28.620089), amrex::Real(193.788772), amrex::Real(308.307024),
185  amrex::Real(320.199637), amrex::Real(279.376984), amrex::Real(87.195000),
186  amrex::Real(349.129677), amrex::Real(128.443387), amrex::Real(154.143880),
187  amrex::Real(291.269597), amrex::Real(114.860583), amrex::Real(332.092251),
188  amrex::Real(296.414411), amrex::Real(145.769910), amrex::Real(337.237063),
189  amrex::Real(152.092288), amrex::Real(126.839891), amrex::Real(210.667199),
190  amrex::Real(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 = { amrex::Real(7391.0225890), amrex::Real(2555.1526947), amrex::Real(2022.7629188),
197  -amrex::Real(1973.6517951), amrex::Real(1240.2321818), amrex::Real(953.8679112),
198  -amrex::Real(931.7537108), amrex::Real(872.3795383), amrex::Real(606.3544732),
199  -amrex::Real(496.0274038), amrex::Real(456.9608039), amrex::Real(346.9462320),
200  -amrex::Real(305.8412902), amrex::Real(249.6173246), -amrex::Real(199.1027200),
201  amrex::Real(191.0560889), -amrex::Real(175.2936572), amrex::Real(165.9068833),
202  amrex::Real(161.1285917), amrex::Real(139.7878093), -amrex::Real(133.5228399),
203  amrex::Real(117.0673811), amrex::Real(104.6907281), amrex::Real(95.3227476),
204  amrex::Real(86.7824524), amrex::Real(86.0857729), amrex::Real(70.5893698),
205  -amrex::Real(69.9719343), -amrex::Real(62.5817473), amrex::Real(61.5450059),
206  -amrex::Real(57.9364011), amrex::Real(57.1899832), -amrex::Real(57.0236109),
207  -amrex::Real(54.2119253), amrex::Real(53.2834147), amrex::Real(52.1223575),
208  -amrex::Real(49.0059908), -amrex::Real(48.3118757), -amrex::Real(45.4191685),
209  -amrex::Real(42.2357920), -amrex::Real(34.7971099), amrex::Real(34.4623613),
210  -amrex::Real(33.8356643), amrex::Real(33.6689362), -amrex::Real(31.2521586),
211  -amrex::Real(30.8798701), amrex::Real(28.4640769), -amrex::Real(27.1960802),
212  amrex::Real(27.0860736), -amrex::Real(26.3437456), amrex::Real(24.7253740),
213  amrex::Real(24.6732126), amrex::Real(24.4272733), amrex::Real(24.0127327),
214  amrex::Real(21.7150294), -amrex::Real(21.5375347), amrex::Real(18.1148363),
215  -amrex::Real(16.9603104), -amrex::Real(16.1765215), amrex::Real(15.5567653),
216  amrex::Real(15.4846529), amrex::Real(15.2150632), amrex::Real(14.5047426),
217  -amrex::Real(14.3873316), amrex::Real(13.1351419), amrex::Real(12.8776311),
218  amrex::Real(11.9867234), amrex::Real(11.9385578), amrex::Real(11.7030822),
219  amrex::Real(11.6018181), -amrex::Real(11.2617293), -amrex::Real(10.4664199),
220  amrex::Real(10.4333970), -amrex::Real(10.2377466), amrex::Real(10.1934446),
221  -amrex::Real(10.1280191), amrex::Real(10.0289441), -amrex::Real(10.0034259)};
222 
223  // rates for mvelp sine series
224  std::vector<real> mvrate = {amrex::Real(31.609974), amrex::Real(32.620504), amrex::Real(24.172203),
225  amrex::Real(0.636717), amrex::Real(31.983787), amrex::Real(3.138886),
226  amrex::Real(30.973257), amrex::Real(44.828336), amrex::Real(0.991874),
227  amrex::Real(0.373813), amrex::Real(43.668246), amrex::Real(32.246691),
228  amrex::Real(30.599444), amrex::Real(2.147012), amrex::Real(10.511172),
229  amrex::Real(42.681324), amrex::Real(13.650058), amrex::Real(0.986922),
230  amrex::Real(9.874455), amrex::Real(13.013341), amrex::Real(0.262904),
231  amrex::Real(0.004952), amrex::Real(1.142024), amrex::Real(63.219948),
232  amrex::Real(0.205021), amrex::Real(2.151964), amrex::Real(64.230478),
233  amrex::Real(43.836462), amrex::Real(47.439436), amrex::Real(1.384343),
234  amrex::Real(7.437771), amrex::Real(18.829299), amrex::Real(9.500642),
235  amrex::Real(0.431696), amrex::Real(1.160090), amrex::Real(55.782177),
236  amrex::Real(12.639528), amrex::Real(1.155138), amrex::Real(0.168216),
237  amrex::Real(1.647247), amrex::Real(10.884985), amrex::Real(5.610937),
238  amrex::Real(12.658184), amrex::Real(1.010530), amrex::Real(1.983748),
239  amrex::Real(14.023871), amrex::Real(0.560178), amrex::Real(1.273434),
240  amrex::Real(12.021467), amrex::Real(62.583231), amrex::Real(63.593761),
241  amrex::Real(76.438310), amrex::Real(4.280910), amrex::Real(13.218362),
242  amrex::Real(17.818769), amrex::Real(8.359495), amrex::Real(56.792707),
243  amrex::Real(8.448301), amrex::Real(1.978796), amrex::Real(8.863925),
244  amrex::Real(0.186365), amrex::Real(8.996212), amrex::Real(6.771027),
245  amrex::Real(45.815258), amrex::Real(12.002811), amrex::Real(75.278220),
246  amrex::Real(65.241008), amrex::Real(18.870667), amrex::Real(22.009553),
247  amrex::Real(64.604291), amrex::Real(11.498094), amrex::Real(0.578834),
248  amrex::Real(9.237738), amrex::Real(49.747842), amrex::Real(2.147012),
249  amrex::Real(1.196895), amrex::Real(2.133898), amrex::Real(0.173168)};
250 
251  // phases for mvelp sine series
252  std::vector<real> mvphas = {amrex::Real(251.9025), amrex::Real(280.8325), amrex::Real(128.3057),
253  amrex::Real(348.1074), amrex::Real(292.7252), amrex::Real(165.1686),
254  amrex::Real(263.7951), amrex::Real(15.3747), amrex::Real(58.5749),
255  amrex::Real(40.8226), amrex::Real(308.4258), amrex::Real(240.0099),
256  amrex::Real(222.9725), amrex::Real(106.5937), amrex::Real(114.5182),
257  amrex::Real(268.7809), amrex::Real(279.6869), amrex::Real(39.6448),
258  amrex::Real(126.4108), amrex::Real(291.5795), amrex::Real(307.2848),
259  amrex::Real(18.9300), amrex::Real(273.7596), amrex::Real(143.8050),
260  amrex::Real(191.8927), amrex::Real(125.5237), amrex::Real(172.7351),
261  amrex::Real(316.7998), amrex::Real(319.6024), amrex::Real(69.7526),
262  amrex::Real(123.5968), amrex::Real(217.6432), amrex::Real(85.5882),
263  amrex::Real(156.2147), amrex::Real(66.9489), amrex::Real(20.2082),
264  amrex::Real(250.7568), amrex::Real(48.0188), amrex::Real(8.3739),
265  amrex::Real(17.0374), amrex::Real(155.3409), amrex::Real(94.1709),
266  amrex::Real(221.1120), amrex::Real(28.9300), amrex::Real(117.1498),
267  amrex::Real(320.5095), amrex::Real(262.3602), amrex::Real(336.2148),
268  amrex::Real(233.0046), amrex::Real(155.6977), amrex::Real(184.6277),
269  amrex::Real(267.2772), amrex::Real(78.9281), amrex::Real(123.4722),
270  amrex::Real(188.7132), amrex::Real(180.1364), amrex::Real(49.1382),
271  amrex::Real(152.5268), amrex::Real(98.2198), amrex::Real(97.4808),
272  amrex::Real(221.5376), amrex::Real(168.2438), amrex::Real(161.1199),
273  amrex::Real(55.0196), amrex::Real(262.6495), amrex::Real(200.3284),
274  amrex::Real(201.6651), amrex::Real(294.6547), amrex::Real(99.8233),
275  amrex::Real(213.5577), amrex::Real(154.1631), amrex::Real(232.7153),
276  amrex::Real(138.3034), amrex::Real(204.6609), amrex::Real(106.5938),
277  amrex::Real(250.4676), amrex::Real(332.3345), amrex::Real(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 = amrex::Real(23.4441);
296  eccen = amrex::Real(0.016715);
297  mvelp = amrex::Real(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. amrex::Real(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(amrex::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 = zero;
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(amrex::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 = zero;
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 = zero;
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) <= amrex::Real(1.0e-8)) {
374  if (sinsum == zero) {
375  fvelp = zero;
376  } else if (sinsum < zero) {
377  fvelp = amrex::Real(1.5)*PI;
378  } else if (sinsum > zero) {
379  fvelp = myhalf*PI;
380  }
381  } else if (cossum < zero) {
382  fvelp = std::atan(sinsum/cossum) + PI;
383  } else {
384  if (sinsum < zero) {
385  fvelp = std::atan(sinsum/cossum) + two*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 = zero;
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(amrex::Real(50.439273))*psecdeg*years + real(amrex::Real(3.392506)) + mvsum;
406 
407  // Cases to make sure mvelp is between 0 and amrex::Real(360.)
408  if (mvelp < zero) {
409  do {
410  mvelp += amrex::Real(360.0);
411  } while (mvelp < zero);
412  }
413  if (mvelp >= amrex::Real(360.0)) {
414  do {
415  mvelp -= amrex::Real(360.0);
416  } while (mvelp >= amrex::Real(360.0));
417  }
418 
419  } // end of test on whether to calculate or use input orbital params
420 
421  // Orbit needs the obliquity in radians
422 
423  obliqr = obliq*degrad;
424 
425  /*
426  180 degrees must be added to mvelp since observations are made from the
427  earth and the sun is considered (wrongly for the algorithm) to go around
428  the earth. For a more graphic explanation see Appendix B in:
429 
430  A. Berger, M. Loutre and C. Tricot. amrex::Real(1993.) Insolation and Earth Orbital
431  Periods. J. of Geophysical Research 98:10,341-10,amrex::Real(362.)
432 
433  Additionally, orbit will need this value in radians. So mvelp becomes
434  mvelpp (mvelp plus pi)
435  */
436 
437  mvelpp = (mvelp + amrex::Real(180.0))*degrad;
438 
439  // Set up an argument used several times in lambm0 calculation ahead.
440 
441  beta = std::sqrt(one - eccen2);
442 
443  /*
444  The mean longitude at the vernal equinox (lambda m nought in Berger
445  1978; in radians) is calculated from the following formula given in
446  Berger amrex::Real(1978.) At the vernal equinox the true longitude (lambda in Berger
447  1978) is zero
448  */
449 
450  lambm0 = two*( (amrex::Real(.5)*eccen + amrex::Real(.125)*eccen3)*(one + beta)*std::sin(mvelpp)
451  - amrex::Real(.250)*eccen2*(amrex::Real(.5) + beta)*std::sin(two*mvelpp)
452  + amrex::Real(.125)*eccen3*(one/three + beta)*std::sin(three*mvelpp) );
453 }
454 
455 
456 AMREX_GPU_HOST
457 AMREX_FORCE_INLINE
458 real
459 orbital_avg_cos_zenith (real& jday,
460  real& lat,
461  real& lon,
462  real& declin,
463  real& dt_avg)
464 {
465  // adjust latitude so that its tangent will be defined
466  real del;
467  if (lat == PIoTwo) {
468  del = lat - real(amrex::Real(1.0e-05));
469  } else if (lat == -PIoTwo) {
470  del = lat + real(amrex::Real(1.0e-05));
471  } else {
472  del = lat;
473  }
474 
475  // adjust declination so that its tangent will be defined
476  real phi;
477  if (declin == PIoTwo) {
478  phi = declin - real(amrex::Real(1.0e-05));
479  } else if (declin == -PIoTwo) {
480  phi = declin + real(amrex::Real(1.0e-05));
481  } else {
482  phi = declin;
483  }
484 
485  // define the cosine of the myhalf-day length
486  // adjust for cases of all daylight or all night
487  real h;
488  real cos_h = - std::tan(del) * std::tan(phi);
489  if (cos_h <= -one) {
490  h = PI;
491  } else if (cos_h >= one) {
492  h = zero;
493  } else {
494  h = std::acos(cos_h);
495  }
496 
497  // Define Local Time t and t + dt
498  // adjust t to be between -pi and pi
499  real t1 = (jday - int(jday)) * two*PI + lon - PI;
500  if (t1 >= PI) {
501  t1 = t1 - two*PI;
502  } else if (t1 < -PI) {
503  t1 = t1 + two*PI;
504  }
505 
506  real dt = dt_avg / real(amrex::Real(86400.0)) * two*PI;
507  real t2 = t1 + dt;
508 
509  // Compute Cosine Solar Zenith angle
510  // define terms needed in the cosine zenith angle equation
511  real aa = std::sin(lat) * std::sin(declin);
512  real bb = std::cos(lat) * std::cos(declin);
513 
514  // define the hour angle
515  // force it to be between -h and h
516  // consider the situation when the night period is too short
517  real tt1,tt2,tt3,tt4;
518  if ( (t2 >= PI) && (t1 <= PI) && ((PI - h) <= dt) ) {
519  tt2 = h;
520  tt1 = std::min(std::max(t1, -h), h);
521  tt4 = std::min(std::max(t2, two*PI - h), two*PI + h);
522  tt3 = two*PI - h;
523  } else if ( (t2 >= -PI) && (t1 <= -PI) && ((PI - h) <= dt) ) {
524  tt2 = - two*PI + h;
525  tt1 = std::min(std::max(t1, -two*PI - h), -two*PI + h);
526  tt4 = std::min(std::max(t2, -h), h);
527  tt3 = -h;
528  } else {
529  if (t2 > PI) {
530  tt2 = std::min(std::max(t2 - two*PI, -h), h);
531  } else if (t2 < -PI) {
532  tt2 = std::min(std::max(t2 + two*PI, -h), h);
533  } else {
534  tt2 = std::min(std::max(t2 , -h), h);
535  }
536 
537  if (t1 > PI) {
538  tt1 = std::min(std::max(t1 - two*PI, -h), h);
539  } else if (t1 < -PI) {
540  tt1 = std::min(std::max(t1 + two*PI, -h), h);
541  } else {
542  tt1 = std::min(std::max(t1 , -h), h);
543  }
544  tt4 = zero;
545  tt3 = zero;
546  }
547 
548  // perform a time integration to obtain cosz if desired
549  // output is valid over the period from t to t + dt
550  if ( (tt2 > tt1) || (tt4 > tt3) ) {
551  return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
552  (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
553  } else {
554  return zero;
555  }
556 }
557 
558 
559 AMREX_GPU_HOST
560 AMREX_FORCE_INLINE
561 real
562 orbital_cos_zenith (real& jday,
563  real& lat,
564  real& lon,
565  real& declin,
566  real dt_avg = -one,
567  real uniform_angle = -one,
568  real constant_zenith_angle_deg = -one)
569 {
570  // Constant zenith angle is true
571  if ( constant_zenith_angle_deg >= zero ) {
572  return std::cos( constant_zenith_angle_deg * PI/amrex::Real(180.) );
573  }
574 
575  // Uniform angle is true
576  if ( uniform_angle >= zero) {
577  return std::cos(uniform_angle);
578  }
579 
580  // Perform the calculation of Orbital_Cos_Zenith
581  bool use_dt_avg = false;
582  if (dt_avg > zero) {
583  use_dt_avg = true;
584  }
585  // If dt for the average cosz is specified, then call the shr_orb_avg_cosz
586  if (use_dt_avg) {
587  return orbital_avg_cos_zenith(jday, lat, lon, declin, dt_avg);
588  } else {
589  return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
590  std::cos((jday-floor(jday))*real(two)*PI + lon);
591  }
592 }
593 #endif
ERF_Constants.H
three
constexpr amrex::Real three
Definition: ERF_Constants.H:9
two
constexpr amrex::Real two
Definition: ERF_Constants.H:8
ORB_UNDEF_INT
static constexpr int ORB_UNDEF_INT
Definition: ERF_Constants.H:106
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:16
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:17
beta
amrex::Real beta
Definition: ERF_InitCustomPert_IsentropicVortex.H:10
real
double real
Definition: ERF_OrbCosZenith.H:9
orbital_cos_zenith
AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_cos_zenith(real &jday, real &lat, real &lon, real &declin, real dt_avg=-one, real uniform_angle=-one, real constant_zenith_angle_deg=-one)
Definition: ERF_OrbCosZenith.H:562
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:459
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)
Definition: ERF_OrbCosZenith.H:15
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)
Definition: ERF_OrbCosZenith.H:81
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
module_model_constants::degrad
real(c_double), parameter degrad
Definition: ERF_module_model_constants.F90:75