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  // Any of obliq/eccen/mvelp that is non-negative on entry was supplied by the
292  // user and must survive this routine; the negative ones are the ones we fill
293  // in below, either from the Berger series or from the AMIP II defaults.
294  real obliq_in = obliq;
295  real eccen_in = eccen;
296  real mvelp_in = mvelp;
297 
298  // Check for flag to use input orbit parameters
299  if ( iyear_AD == ORB_UNDEF_INT ) {
300 
301  // No orbital year is available, so the AMIP II settings (for a 1995
302  // orbit) are adopted
303  obliq = amrex::Real(23.4441);
304  eccen = amrex::Real(0.016715);
305  mvelp = amrex::Real(102.7);
306  eccen2 = eccen*eccen;
307  eccen3 = eccen2*eccen;
308 
309  } else { // Otherwise calculate based on years before present
310 
311  /*
312  The following calculates the earths obliquity, orbital eccentricity
313  (and various powers of it) and vernal equinox mean longitude of
314  perihelion for years in the past (future = negative of years past),
315  using constants (see parameter section) given in the program of:
316 
317  Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations
318  of Daily Insolation. Contribution 18, Institute of Astronomy and
319  Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
320 
321  and formulas given in the paper (where less precise constants are also
322  given):
323 
324  Berger, Andre. amrex::Real(1978.) Long-Term Variations of Daily Insolation and
325  Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367
326 
327  The algorithm is valid only to 1,000,000 years past or hence.
328  For a solution valid to 5-10 million years past see the above author.
329  Algorithm below is better for years closer to present than is the
330  5-10 million year solution.
331 
332  Years to time of interest must be negative of years before present
333  (1950) in formulas that follow.
334  */
335 
336  yb4_1950AD = real(amrex::Real(1950.0)) - real(iyear_AD);
337  years = - yb4_1950AD;
338 
339  /*
340  In the summations below, cosine or sine arguments, which end up in
341  degrees, must be converted to radians via multiplication by degrad.
342 
343  Summation of cosine series for obliquity (epsilon in Berger 1978) in
344  degrees. Convert the amplitudes and rates, which are in arc secs, into
345  degrees via multiplication by psecdeg (arc seconds to degrees conversion
346  factor). For obliq, first term is Berger 1978 epsilon star; second
347  term is series summation in degrees.
348  */
349 
350  obsum = zero;
351  for (int i(0); i<poblen; ++i) {
352  obsum = obsum + obamp[i]*psecdeg*std::cos( (obrate[i]*psecdeg*years + obphas[i]) * degrad );
353  }
354  obliq = real(amrex::Real(23.320556)) + obsum;
355 
356  /*
357  Summation of cosine and sine series for computation of eccentricity
358  (eccen; e in Berger 1978) and fixed vernal equinox longitude of
359  perihelion (fvelp; pi in Berger 1978), which is used for computation
360  of moving vernal equinox longitude of perihelion. Convert the rates,
361  which are in arc seconds, into degrees via multiplication by psecdeg.
362  */
363 
364  cossum = zero;
365  for (int i(0); i<pecclen; ++i) {
366  cossum = cossum + ecamp[i]*std::cos( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
367  }
368 
369  sinsum = zero;
370  for (int i(0); i<pecclen; ++i) {
371  sinsum = sinsum + ecamp[i]*std::sin( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
372  }
373 
374  // Use summations to calculate eccentricity
375 
376  eccen2 = cossum*cossum + sinsum*sinsum;
377  eccen = std::sqrt(eccen2);
378  eccen3 = eccen2*eccen;
379 
380  // A series of cases for fvelp, which is in radians.
381  if (std::fabs(cossum) <= amrex::Real(1.0e-8)) {
382  if (sinsum == zero) {
383  fvelp = zero;
384  } else if (sinsum < zero) {
385  fvelp = amrex::Real(1.5)*PI;
386  } else if (sinsum > zero) {
387  fvelp = myhalf*PI;
388  }
389  } else if (cossum < zero) {
390  fvelp = std::atan(sinsum/cossum) + PI;
391  } else {
392  if (sinsum < zero) {
393  fvelp = std::atan(sinsum/cossum) + two*PI;
394  } else {
395  fvelp = std::atan(sinsum/cossum);
396  }
397  }
398 
399  /*
400  Summation of sin series for computation of moving vernal equinox long
401  of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp,
402  first term is fvelp in degrees; second term is Berger 1978 psi bar
403  times years and in degrees; third term is Berger 1978 zeta; fourth
404  term is series summation in degrees. Convert the amplitudes and rates,
405  which are in arc seconds, into degrees via multiplication by psecdeg.
406  Series summation plus second and third terms constitute Berger 1978
407  psi, which is the general precession.
408  */
409  mvsum = zero;
410  for (int i(0); i<pmvelen; ++i) {
411  mvsum = mvsum + mvamp[i]*psecdeg*std::sin( (mvrate[i]*psecdeg*years + mvphas[i]) * degrad);
412  }
413  mvelp = fvelp/degrad + real(amrex::Real(50.439273))*psecdeg*years + real(amrex::Real(3.392506)) + mvsum;
414 
415  } // end of test on whether to calculate or use input orbital params
416 
417  // Restore any parameter the user specified explicitly, overriding whatever
418  // value was just computed for it
419  if (obliq_in >= zero) { obliq = obliq_in; }
420  if (mvelp_in >= zero) { mvelp = mvelp_in; }
421  if (eccen_in >= zero) {
422  eccen = eccen_in;
423  eccen2 = eccen*eccen;
424  eccen3 = eccen2*eccen;
425  }
426 
427  // Cases to make sure mvelp is between 0 and amrex::Real(360.)
428  if (mvelp < zero) {
429  do {
430  mvelp += amrex::Real(360.0);
431  } while (mvelp < zero);
432  }
433  if (mvelp >= amrex::Real(360.0)) {
434  do {
435  mvelp -= amrex::Real(360.0);
436  } while (mvelp >= amrex::Real(360.0));
437  }
438 
439  // Orbit needs the obliquity in radians
440 
441  obliqr = obliq*degrad;
442 
443  /*
444  180 degrees must be added to mvelp since observations are made from the
445  earth and the sun is considered (wrongly for the algorithm) to go around
446  the earth. For a more graphic explanation see Appendix B in:
447 
448  A. Berger, M. Loutre and C. Tricot. amrex::Real(1993.) Insolation and Earth Orbital
449  Periods. J. of Geophysical Research 98:10,341-10,amrex::Real(362.)
450 
451  Additionally, orbit will need this value in radians. So mvelp becomes
452  mvelpp (mvelp plus pi)
453  */
454 
455  mvelpp = (mvelp + amrex::Real(180.0))*degrad;
456 
457  // Set up an argument used several times in lambm0 calculation ahead.
458 
459  beta = std::sqrt(one - eccen2);
460 
461  /*
462  The mean longitude at the vernal equinox (lambda m nought in Berger
463  1978; in radians) is calculated from the following formula given in
464  Berger amrex::Real(1978.) At the vernal equinox the true longitude (lambda in Berger
465  1978) is zero
466  */
467 
468  lambm0 = two*( (amrex::Real(.5)*eccen + amrex::Real(.125)*eccen3)*(one + beta)*std::sin(mvelpp)
469  - amrex::Real(.250)*eccen2*(amrex::Real(.5) + beta)*std::sin(two*mvelpp)
470  + amrex::Real(.125)*eccen3*(one/three + beta)*std::sin(three*mvelpp) );
471 }
472 
473 
474 AMREX_GPU_HOST
475 AMREX_FORCE_INLINE
476 real
477 orbital_avg_cos_zenith (real& jday,
478  real& lat,
479  real& lon,
480  real& declin,
481  real& dt_avg)
482 {
483  // adjust latitude so that its tangent will be defined
484  real del;
485  if (lat == PIoTwo) {
486  del = lat - real(amrex::Real(1.0e-05));
487  } else if (lat == -PIoTwo) {
488  del = lat + real(amrex::Real(1.0e-05));
489  } else {
490  del = lat;
491  }
492 
493  // adjust declination so that its tangent will be defined
494  real phi;
495  if (declin == PIoTwo) {
496  phi = declin - real(amrex::Real(1.0e-05));
497  } else if (declin == -PIoTwo) {
498  phi = declin + real(amrex::Real(1.0e-05));
499  } else {
500  phi = declin;
501  }
502 
503  // define the cosine of the myhalf-day length
504  // adjust for cases of all daylight or all night
505  real h;
506  real cos_h = - std::tan(del) * std::tan(phi);
507  if (cos_h <= -one) {
508  h = PI;
509  } else if (cos_h >= one) {
510  h = zero;
511  } else {
512  h = std::acos(cos_h);
513  }
514 
515  // Define Local Time t and t + dt
516  // adjust t to be between -pi and pi
517  real t1 = (jday - int(jday)) * two*PI + lon - PI;
518  if (t1 >= PI) {
519  t1 = t1 - two*PI;
520  } else if (t1 < -PI) {
521  t1 = t1 + two*PI;
522  }
523 
524  real dt = dt_avg / real(amrex::Real(86400.0)) * two*PI;
525  real t2 = t1 + dt;
526 
527  // Compute Cosine Solar Zenith angle
528  // define terms needed in the cosine zenith angle equation
529  real aa = std::sin(lat) * std::sin(declin);
530  real bb = std::cos(lat) * std::cos(declin);
531 
532  // define the hour angle
533  // force it to be between -h and h
534  // consider the situation when the night period is too short
535  real tt1,tt2,tt3,tt4;
536  if ( (t2 >= PI) && (t1 <= PI) && ((PI - h) <= dt) ) {
537  tt2 = h;
538  tt1 = std::min(std::max(t1, -h), h);
539  tt4 = std::min(std::max(t2, two*PI - h), two*PI + h);
540  tt3 = two*PI - h;
541  } else if ( (t2 >= -PI) && (t1 <= -PI) && ((PI - h) <= dt) ) {
542  tt2 = - two*PI + h;
543  tt1 = std::min(std::max(t1, -two*PI - h), -two*PI + h);
544  tt4 = std::min(std::max(t2, -h), h);
545  tt3 = -h;
546  } else {
547  if (t2 > PI) {
548  tt2 = std::min(std::max(t2 - two*PI, -h), h);
549  } else if (t2 < -PI) {
550  tt2 = std::min(std::max(t2 + two*PI, -h), h);
551  } else {
552  tt2 = std::min(std::max(t2 , -h), h);
553  }
554 
555  if (t1 > PI) {
556  tt1 = std::min(std::max(t1 - two*PI, -h), h);
557  } else if (t1 < -PI) {
558  tt1 = std::min(std::max(t1 + two*PI, -h), h);
559  } else {
560  tt1 = std::min(std::max(t1 , -h), h);
561  }
562  tt4 = zero;
563  tt3 = zero;
564  }
565 
566  // perform a time integration to obtain cosz if desired
567  // output is valid over the period from t to t + dt
568  if ( (tt2 > tt1) || (tt4 > tt3) ) {
569  return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
570  (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
571  } else {
572  return zero;
573  }
574 }
575 
576 
577 AMREX_GPU_HOST
578 AMREX_FORCE_INLINE
579 real
580 orbital_cos_zenith (real& jday,
581  real& lat,
582  real& lon,
583  real& declin,
584  real dt_avg = -one,
585  real uniform_angle = -one,
586  real constant_zenith_angle_deg = -one)
587 {
588  // Constant zenith angle is true
589  if ( constant_zenith_angle_deg >= zero ) {
590  return std::cos( constant_zenith_angle_deg * PI/amrex::Real(180.) );
591  }
592 
593  // Uniform angle is true
594  if ( uniform_angle >= zero) {
595  return std::cos(uniform_angle);
596  }
597 
598  // Perform the calculation of Orbital_Cos_Zenith
599  bool use_dt_avg = false;
600  if (dt_avg > zero) {
601  use_dt_avg = true;
602  }
603  // If dt for the average cosz is specified, then call the shr_orb_avg_cosz
604  if (use_dt_avg) {
605  return orbital_avg_cos_zenith(jday, lat, lon, declin, dt_avg);
606  } else {
607  return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
608  std::cos((jday-floor(jday))*real(two)*PI + lon);
609  }
610 }
611 #endif
ERF_Constants.H
three
constexpr amrex::Real three
Definition: ERF_Constants.H:11
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
ORB_UNDEF_INT
static constexpr int ORB_UNDEF_INT
Definition: ERF_Constants.H:138
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
PI
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:43
beta
amrex::Real beta
Definition: ERF_InitCustomPert_DataAssimilation_ISV.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:580
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:477
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