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 <AMReX_REAL.H>
8 #include <AMReX_GpuQualifiers.H>
9 #include <ERF_Constants.H>
10 
11 // The orbital routines work in double regardless of amrex::Real: a calendar
12 // day at ~1e2 with a fraction and the trig series need the precision.
13 
14 
15 AMREX_GPU_HOST
16 AMREX_FORCE_INLINE
17 void
18 orbital_decl (double& calday,
19  double& eccen,
20  double& mvelpp,
21  double& lambm0,
22  double& obliqr,
23  double& delta,
24  double& eccf)
25 {
26 
27  double lambm; // Lambda m, mean long of perihelion (rad)
28  double lmm; // Intermediate argument involving lambm
29  double lamb; // Lambda, the earths long of perihelion
30  double invrho; // Inverse normalized sun/earth distance
31  double sinl; // Sine of lmm
32  static constexpr double dayspy = double(amrex::Real(365.0)); // Day per year
33  static constexpr double ve = double( amrex::Real(80.5)); // Calday of vernal equinox
34 
35  /*
36  Compute eccentricity factor and solar declination using
37  day value where a round day (such as amrex::Real(213.0)) refers to 0z at
38  Greenwich longitude.
39 
40  Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
41  Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
42  35:2362-amrex::Real(2367.)
43 
44  To get the earths true longitude (position in orbit; lambda in Berger
45  1978) which is necessary to find the eccentricity factor and declination,
46  must first calculate the mean longitude (lambda m in Berger 1978) at
47  the present day. This is done by adding to lambm0 (the mean longitude
48  at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
49  an increment (delta lambda m in Berger 1978) that is the number of
50  days past or before (a negative increment) the vernal equinox divided by
51  the days in a model year times the 2*pi radians in a complete orbit.
52  */
53 
54  lambm = lambm0 + (calday - ve)*two*PI/dayspy;
55  lmm = lambm - mvelpp;
56 
57  // The earths true longitude, in radians, is then found from
58  // the formula in Berger 1978:
59 
60  sinl = std::sin(lmm);
61  lamb = lambm + eccen*( two*sinl + eccen*( amrex::Real(1.25)*sin(two*lmm)
62  + eccen*((amrex::Real(13.0)/amrex::Real(12.0))*std::sin(three*lmm) - fourth*sinl)));
63 
64  /*
65  Using the obliquity, eccentricity, moving vernal equinox longitude of
66  perihelion (plus), and earths true longitude, the declination (delta)
67  and the normalized earth/sun distance (rho in Berger 1978; actually inverse
68  rho will be used), and thus the eccentricity factor (eccf), can be
69  calculated from formulas given in Berger amrex::Real(1978.)
70  */
71 
72  invrho = (one + eccen*std::cos(lamb - mvelpp)) / (one - eccen*eccen);
73 
74  // Set solar declination and eccentricity factor
75 
76  delta = std::asin(std::sin(obliqr)*std::sin(lamb));
77  eccf = invrho*invrho;
78 }
79 
80 
81 AMREX_GPU_HOST
82 AMREX_FORCE_INLINE
83 void
84 orbital_params (int& iyear_AD,
85  double& eccen,
86  double& obliq,
87  double& mvelp,
88  double& obliqr,
89  double& lambm0,
90  double& mvelpp)
91 {
92  /*
93  !-------------------------------------------------------------------------------
94  !
95  ! Calculate earths orbital parameters using Dave Threshers formula which
96  ! came from Berger, Andre. 1978 "A Simple Algorithm to Compute Long-Term
97  ! Variations of Daily Insolation". Contribution 18, Institute of Astronomy
98  ! and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium
99  !
100  !-------------------------------------------------------------------------------
101  */
102 
103  int poblen = 47; // # of elements in series wrt obliquity
104  int pecclen = 19; // # of elements in series wrt eccentricity
105  int pmvelen = 78; // # of elements in series wrt vernal equinox
106  static constexpr double psecdeg = double(one)/double(amrex::Real(3600.0)); // arc sec to deg conversion
107  static constexpr double degrad = PI/double(amrex::Real(180.)); // degree to radian conversion factor
108 
109  // Cosine series data for computation of obliquity: amplitude (arc seconds),
110  // rate (arc seconds/year), phase (degrees).
111  // amplitudes for obliquity cos series
112  std::vector<double> obamp = {-amrex::Real(2462.2214466), -amrex::Real(857.3232075), -amrex::Real(629.3231835),
113  -amrex::Real(414.2804924), -amrex::Real(311.7632587), amrex::Real(308.9408604),
114  -amrex::Real(162.5533601), -amrex::Real(116.1077911), amrex::Real(101.1189923),
115  -amrex::Real(67.6856209), amrex::Real(24.9079067), amrex::Real(22.5811241),
116  -amrex::Real(21.1648355), -amrex::Real(15.6549876), amrex::Real(15.3936813),
117  amrex::Real(14.6660938), -amrex::Real(11.7273029), amrex::Real(10.2742696),
118  amrex::Real(6.4914588), amrex::Real(5.8539148), -amrex::Real(5.4872205),
119  -amrex::Real(5.4290191), amrex::Real(5.1609570), amrex::Real(5.0786314),
120  -amrex::Real(4.0735782), amrex::Real(3.7227167), amrex::Real(3.3971932),
121  -amrex::Real(2.8347004), -amrex::Real(2.6550721), -amrex::Real(2.5717867),
122  -amrex::Real(2.4712188), amrex::Real(2.4625410), amrex::Real(2.2464112),
123  -amrex::Real(2.0755511), -amrex::Real(1.9713669), -amrex::Real(1.8813061),
124  -amrex::Real(1.8468785), amrex::Real(1.8186742), amrex::Real(1.7601888),
125  -amrex::Real(1.5428851), amrex::Real(1.4738838), -amrex::Real(1.4593669),
126  amrex::Real(1.4192259), -amrex::Real(1.1818980), amrex::Real(1.1756474),
127  -amrex::Real(1.1316126), amrex::Real(1.0896928)};
128  // rates for obliquity cosine series
129  std::vector<double> obrate = {amrex::Real(31.609974), amrex::Real(32.620504), amrex::Real(24.172203),
130  amrex::Real(31.983787), amrex::Real(44.828336), amrex::Real(30.973257),
131  amrex::Real(43.668246), amrex::Real(32.246691), amrex::Real(30.599444),
132  amrex::Real(42.681324), amrex::Real(43.836462), amrex::Real(47.439436),
133  amrex::Real(63.219948), amrex::Real(64.230478), amrex::Real(1.010530),
134  amrex::Real(7.437771), amrex::Real(55.782177), amrex::Real(0.373813),
135  amrex::Real(13.218362), amrex::Real(62.583231), amrex::Real(63.593761),
136  amrex::Real(76.438310), amrex::Real(45.815258), amrex::Real(8.448301),
137  amrex::Real(56.792707), amrex::Real(49.747842), amrex::Real(12.058272),
138  amrex::Real(75.278220), amrex::Real(65.241008), amrex::Real(64.604291),
139  amrex::Real(1.647247), amrex::Real(7.811584), amrex::Real(12.207832),
140  amrex::Real(63.856665), amrex::Real(56.155990), amrex::Real(77.448840),
141  amrex::Real(6.801054), amrex::Real(62.209418), amrex::Real(20.656133),
142  amrex::Real(48.344406), amrex::Real(55.145460), amrex::Real(69.000539),
143  amrex::Real(11.071350), amrex::Real(74.291298), amrex::Real(11.047742),
144  amrex::Real(0.636717), amrex::Real(12.844549)};
145 
146  // phases for obliquity cosine series
147  std::vector<double> obphas = {amrex::Real(251.9025), amrex::Real(280.8325), amrex::Real(128.3057),
148  amrex::Real(292.7252), amrex::Real(15.3747), amrex::Real(263.7951),
149  amrex::Real(308.4258), amrex::Real(240.0099), amrex::Real(222.9725),
150  amrex::Real(268.7809), amrex::Real(316.7998), amrex::Real(319.6024),
151  amrex::Real(143.8050), amrex::Real(172.7351), amrex::Real(28.9300),
152  amrex::Real(123.5968), amrex::Real(20.2082), amrex::Real(40.8226),
153  amrex::Real(123.4722), amrex::Real(155.6977), amrex::Real(184.6277),
154  amrex::Real(267.2772), amrex::Real(55.0196), amrex::Real(152.5268),
155  amrex::Real(49.1382), amrex::Real(204.6609), amrex::Real(56.5233),
156  amrex::Real(200.3284), amrex::Real(201.6651), amrex::Real(213.5577),
157  amrex::Real(17.0374), amrex::Real(164.4194), amrex::Real(94.5422),
158  amrex::Real(131.9124), amrex::Real(61.0309), amrex::Real(296.2073),
159  amrex::Real(135.4894), amrex::Real(114.8750), amrex::Real(247.0691),
160  amrex::Real(256.6114), amrex::Real(32.1008), amrex::Real(143.6804),
161  amrex::Real(16.8784), amrex::Real(160.6835), amrex::Real(27.5932),
162  amrex::Real(348.1074), amrex::Real(82.6496)};
163 
164  // Cosine/sine series data for computation of eccentricity and fixed vernal
165  // equinox longitude of perihelion (fvelp): amplitude,
166  // rate (arc seconds/year), phase (degrees).
167 
168  // ampl for eccen/fvelp cos/sin series
169  std::vector<double> ecamp = { amrex::Real(0.01860798), amrex::Real(0.01627522), -amrex::Real(0.01300660),
170  amrex::Real(0.00988829), -amrex::Real(0.00336700), amrex::Real(0.00333077),
171  -amrex::Real(0.00235400), amrex::Real(0.00140015), amrex::Real(0.00100700),
172  amrex::Real(0.00085700), amrex::Real(0.00064990), amrex::Real(0.00059900),
173  amrex::Real(0.00037800), -amrex::Real(0.00033700), amrex::Real(0.00027600),
174  amrex::Real(0.00018200), -amrex::Real(0.00017400), -amrex::Real(0.00012400),
175  amrex::Real(0.00001250)};
176 
177  // rates for eccen/fvelp cos/sin series
178  std::vector<double> ecrate = { amrex::Real(4.2072050), amrex::Real(7.3460910), amrex::Real(17.8572630),
179  amrex::Real(17.2205460), amrex::Real(16.8467330), amrex::Real(5.1990790),
180  amrex::Real(18.2310760), amrex::Real(26.2167580), amrex::Real(6.3591690),
181  amrex::Real(16.2100160), amrex::Real(3.0651810), amrex::Real(16.5838290),
182  amrex::Real(18.4939800), amrex::Real(6.1909530), amrex::Real(18.8677930),
183  amrex::Real(17.4255670), amrex::Real(6.1860010), amrex::Real(18.4174410),
184  amrex::Real(0.6678630)};
185 
186  // phases for eccen/fvelp cos/sin series
187  std::vector<double> ecphas = { amrex::Real(28.620089), amrex::Real(193.788772), amrex::Real(308.307024),
188  amrex::Real(320.199637), amrex::Real(279.376984), amrex::Real(87.195000),
189  amrex::Real(349.129677), amrex::Real(128.443387), amrex::Real(154.143880),
190  amrex::Real(291.269597), amrex::Real(114.860583), amrex::Real(332.092251),
191  amrex::Real(296.414411), amrex::Real(145.769910), amrex::Real(337.237063),
192  amrex::Real(152.092288), amrex::Real(126.839891), amrex::Real(210.667199),
193  amrex::Real(72.108838)};
194 
195  // Sine series data for computation of moving vernal equinox longitude of
196  // perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees).
197 
198  // amplitudes for mvelp sine series
199  std::vector<double> mvamp = { amrex::Real(7391.0225890), amrex::Real(2555.1526947), amrex::Real(2022.7629188),
200  -amrex::Real(1973.6517951), amrex::Real(1240.2321818), amrex::Real(953.8679112),
201  -amrex::Real(931.7537108), amrex::Real(872.3795383), amrex::Real(606.3544732),
202  -amrex::Real(496.0274038), amrex::Real(456.9608039), amrex::Real(346.9462320),
203  -amrex::Real(305.8412902), amrex::Real(249.6173246), -amrex::Real(199.1027200),
204  amrex::Real(191.0560889), -amrex::Real(175.2936572), amrex::Real(165.9068833),
205  amrex::Real(161.1285917), amrex::Real(139.7878093), -amrex::Real(133.5228399),
206  amrex::Real(117.0673811), amrex::Real(104.6907281), amrex::Real(95.3227476),
207  amrex::Real(86.7824524), amrex::Real(86.0857729), amrex::Real(70.5893698),
208  -amrex::Real(69.9719343), -amrex::Real(62.5817473), amrex::Real(61.5450059),
209  -amrex::Real(57.9364011), amrex::Real(57.1899832), -amrex::Real(57.0236109),
210  -amrex::Real(54.2119253), amrex::Real(53.2834147), amrex::Real(52.1223575),
211  -amrex::Real(49.0059908), -amrex::Real(48.3118757), -amrex::Real(45.4191685),
212  -amrex::Real(42.2357920), -amrex::Real(34.7971099), amrex::Real(34.4623613),
213  -amrex::Real(33.8356643), amrex::Real(33.6689362), -amrex::Real(31.2521586),
214  -amrex::Real(30.8798701), amrex::Real(28.4640769), -amrex::Real(27.1960802),
215  amrex::Real(27.0860736), -amrex::Real(26.3437456), amrex::Real(24.7253740),
216  amrex::Real(24.6732126), amrex::Real(24.4272733), amrex::Real(24.0127327),
217  amrex::Real(21.7150294), -amrex::Real(21.5375347), amrex::Real(18.1148363),
218  -amrex::Real(16.9603104), -amrex::Real(16.1765215), amrex::Real(15.5567653),
219  amrex::Real(15.4846529), amrex::Real(15.2150632), amrex::Real(14.5047426),
220  -amrex::Real(14.3873316), amrex::Real(13.1351419), amrex::Real(12.8776311),
221  amrex::Real(11.9867234), amrex::Real(11.9385578), amrex::Real(11.7030822),
222  amrex::Real(11.6018181), -amrex::Real(11.2617293), -amrex::Real(10.4664199),
223  amrex::Real(10.4333970), -amrex::Real(10.2377466), amrex::Real(10.1934446),
224  -amrex::Real(10.1280191), amrex::Real(10.0289441), -amrex::Real(10.0034259)};
225 
226  // rates for mvelp sine series
227  std::vector<double> mvrate = {amrex::Real(31.609974), amrex::Real(32.620504), amrex::Real(24.172203),
228  amrex::Real(0.636717), amrex::Real(31.983787), amrex::Real(3.138886),
229  amrex::Real(30.973257), amrex::Real(44.828336), amrex::Real(0.991874),
230  amrex::Real(0.373813), amrex::Real(43.668246), amrex::Real(32.246691),
231  amrex::Real(30.599444), amrex::Real(2.147012), amrex::Real(10.511172),
232  amrex::Real(42.681324), amrex::Real(13.650058), amrex::Real(0.986922),
233  amrex::Real(9.874455), amrex::Real(13.013341), amrex::Real(0.262904),
234  amrex::Real(0.004952), amrex::Real(1.142024), amrex::Real(63.219948),
235  amrex::Real(0.205021), amrex::Real(2.151964), amrex::Real(64.230478),
236  amrex::Real(43.836462), amrex::Real(47.439436), amrex::Real(1.384343),
237  amrex::Real(7.437771), amrex::Real(18.829299), amrex::Real(9.500642),
238  amrex::Real(0.431696), amrex::Real(1.160090), amrex::Real(55.782177),
239  amrex::Real(12.639528), amrex::Real(1.155138), amrex::Real(0.168216),
240  amrex::Real(1.647247), amrex::Real(10.884985), amrex::Real(5.610937),
241  amrex::Real(12.658184), amrex::Real(1.010530), amrex::Real(1.983748),
242  amrex::Real(14.023871), amrex::Real(0.560178), amrex::Real(1.273434),
243  amrex::Real(12.021467), amrex::Real(62.583231), amrex::Real(63.593761),
244  amrex::Real(76.438310), amrex::Real(4.280910), amrex::Real(13.218362),
245  amrex::Real(17.818769), amrex::Real(8.359495), amrex::Real(56.792707),
246  amrex::Real(8.448301), amrex::Real(1.978796), amrex::Real(8.863925),
247  amrex::Real(0.186365), amrex::Real(8.996212), amrex::Real(6.771027),
248  amrex::Real(45.815258), amrex::Real(12.002811), amrex::Real(75.278220),
249  amrex::Real(65.241008), amrex::Real(18.870667), amrex::Real(22.009553),
250  amrex::Real(64.604291), amrex::Real(11.498094), amrex::Real(0.578834),
251  amrex::Real(9.237738), amrex::Real(49.747842), amrex::Real(2.147012),
252  amrex::Real(1.196895), amrex::Real(2.133898), amrex::Real(0.173168)};
253 
254  // phases for mvelp sine series
255  std::vector<double> mvphas = {amrex::Real(251.9025), amrex::Real(280.8325), amrex::Real(128.3057),
256  amrex::Real(348.1074), amrex::Real(292.7252), amrex::Real(165.1686),
257  amrex::Real(263.7951), amrex::Real(15.3747), amrex::Real(58.5749),
258  amrex::Real(40.8226), amrex::Real(308.4258), amrex::Real(240.0099),
259  amrex::Real(222.9725), amrex::Real(106.5937), amrex::Real(114.5182),
260  amrex::Real(268.7809), amrex::Real(279.6869), amrex::Real(39.6448),
261  amrex::Real(126.4108), amrex::Real(291.5795), amrex::Real(307.2848),
262  amrex::Real(18.9300), amrex::Real(273.7596), amrex::Real(143.8050),
263  amrex::Real(191.8927), amrex::Real(125.5237), amrex::Real(172.7351),
264  amrex::Real(316.7998), amrex::Real(319.6024), amrex::Real(69.7526),
265  amrex::Real(123.5968), amrex::Real(217.6432), amrex::Real(85.5882),
266  amrex::Real(156.2147), amrex::Real(66.9489), amrex::Real(20.2082),
267  amrex::Real(250.7568), amrex::Real(48.0188), amrex::Real(8.3739),
268  amrex::Real(17.0374), amrex::Real(155.3409), amrex::Real(94.1709),
269  amrex::Real(221.1120), amrex::Real(28.9300), amrex::Real(117.1498),
270  amrex::Real(320.5095), amrex::Real(262.3602), amrex::Real(336.2148),
271  amrex::Real(233.0046), amrex::Real(155.6977), amrex::Real(184.6277),
272  amrex::Real(267.2772), amrex::Real(78.9281), amrex::Real(123.4722),
273  amrex::Real(188.7132), amrex::Real(180.1364), amrex::Real(49.1382),
274  amrex::Real(152.5268), amrex::Real(98.2198), amrex::Real(97.4808),
275  amrex::Real(221.5376), amrex::Real(168.2438), amrex::Real(161.1199),
276  amrex::Real(55.0196), amrex::Real(262.6495), amrex::Real(200.3284),
277  amrex::Real(201.6651), amrex::Real(294.6547), amrex::Real(99.8233),
278  amrex::Real(213.5577), amrex::Real(154.1631), amrex::Real(232.7153),
279  amrex::Real(138.3034), amrex::Real(204.6609), amrex::Real(106.5938),
280  amrex::Real(250.4676), amrex::Real(332.3345), amrex::Real(27.3039)};
281 
282  //---------------------------Local variables----------------------------------
283  double obsum; // Obliquity series summation
284  double cossum; // Cos series summation for eccentricity/fvelp
285  double sinsum; // Sin series summation for eccentricity/fvelp
286  double fvelp; // Fixed vernal equinox long of perihelion
287  double mvsum; // mvelp series summation
288  double beta; // Intermediate argument for lambm0
289  double years; // Years to time of interest ( pos <=> future)
290  double eccen2; // eccentricity squared
291  double eccen3; // eccentricity cubed
292  double yb4_1950AD; // number of years before 1950 AD
293 
294  // Any of obliq/eccen/mvelp that is non-negative on entry was supplied by the
295  // user and must survive this routine; the negative ones are the ones we fill
296  // in below, either from the Berger series or from the AMIP II defaults.
297  double obliq_in = obliq;
298  double eccen_in = eccen;
299  double mvelp_in = mvelp;
300 
301  // Check for flag to use input orbit parameters
302  if ( iyear_AD == ORB_UNDEF_INT ) {
303 
304  // No orbital year is available, so the AMIP II settings (for a 1995
305  // orbit) are adopted
306  obliq = amrex::Real(23.4441);
307  eccen = amrex::Real(0.016715);
308  mvelp = amrex::Real(102.7);
309  eccen2 = eccen*eccen;
310  eccen3 = eccen2*eccen;
311 
312  } else { // Otherwise calculate based on years before present
313 
314  /*
315  The following calculates the earths obliquity, orbital eccentricity
316  (and various powers of it) and vernal equinox mean longitude of
317  perihelion for years in the past (future = negative of years past),
318  using constants (see parameter section) given in the program of:
319 
320  Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations
321  of Daily Insolation. Contribution 18, Institute of Astronomy and
322  Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
323 
324  and formulas given in the paper (where less precise constants are also
325  given):
326 
327  Berger, Andre. amrex::Real(1978.) Long-Term Variations of Daily Insolation and
328  Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367
329 
330  The algorithm is valid only to 1,000,000 years past or hence.
331  For a solution valid to 5-10 million years past see the above author.
332  Algorithm below is better for years closer to present than is the
333  5-10 million year solution.
334 
335  Years to time of interest must be negative of years before present
336  (1950) in formulas that follow.
337  */
338 
339  yb4_1950AD = double(amrex::Real(1950.0)) - double(iyear_AD);
340  years = - yb4_1950AD;
341 
342  /*
343  In the summations below, cosine or sine arguments, which end up in
344  degrees, must be converted to radians via multiplication by degrad.
345 
346  Summation of cosine series for obliquity (epsilon in Berger 1978) in
347  degrees. Convert the amplitudes and rates, which are in arc secs, into
348  degrees via multiplication by psecdeg (arc seconds to degrees conversion
349  factor). For obliq, first term is Berger 1978 epsilon star; second
350  term is series summation in degrees.
351  */
352 
353  obsum = zero;
354  for (int i(0); i<poblen; ++i) {
355  obsum = obsum + obamp[i]*psecdeg*std::cos( (obrate[i]*psecdeg*years + obphas[i]) * degrad );
356  }
357  obliq = double(amrex::Real(23.320556)) + obsum;
358 
359  /*
360  Summation of cosine and sine series for computation of eccentricity
361  (eccen; e in Berger 1978) and fixed vernal equinox longitude of
362  perihelion (fvelp; pi in Berger 1978), which is used for computation
363  of moving vernal equinox longitude of perihelion. Convert the rates,
364  which are in arc seconds, into degrees via multiplication by psecdeg.
365  */
366 
367  cossum = zero;
368  for (int i(0); i<pecclen; ++i) {
369  cossum = cossum + ecamp[i]*std::cos( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
370  }
371 
372  sinsum = zero;
373  for (int i(0); i<pecclen; ++i) {
374  sinsum = sinsum + ecamp[i]*std::sin( (ecrate[i]*psecdeg*years+ecphas[i]) * degrad );
375  }
376 
377  // Use summations to calculate eccentricity
378 
379  eccen2 = cossum*cossum + sinsum*sinsum;
380  eccen = std::sqrt(eccen2);
381  eccen3 = eccen2*eccen;
382 
383  // A series of cases for fvelp, which is in radians.
384  if (std::fabs(cossum) <= amrex::Real(1.0e-8)) {
385  if (sinsum == zero) {
386  fvelp = zero;
387  } else if (sinsum < zero) {
388  fvelp = amrex::Real(1.5)*PI;
389  } else if (sinsum > zero) {
390  fvelp = myhalf*PI;
391  }
392  } else if (cossum < zero) {
393  fvelp = std::atan(sinsum/cossum) + PI;
394  } else {
395  if (sinsum < zero) {
396  fvelp = std::atan(sinsum/cossum) + two*PI;
397  } else {
398  fvelp = std::atan(sinsum/cossum);
399  }
400  }
401 
402  /*
403  Summation of sin series for computation of moving vernal equinox long
404  of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp,
405  first term is fvelp in degrees; second term is Berger 1978 psi bar
406  times years and in degrees; third term is Berger 1978 zeta; fourth
407  term is series summation in degrees. Convert the amplitudes and rates,
408  which are in arc seconds, into degrees via multiplication by psecdeg.
409  Series summation plus second and third terms constitute Berger 1978
410  psi, which is the general precession.
411  */
412  mvsum = zero;
413  for (int i(0); i<pmvelen; ++i) {
414  mvsum = mvsum + mvamp[i]*psecdeg*std::sin( (mvrate[i]*psecdeg*years + mvphas[i]) * degrad);
415  }
416  mvelp = fvelp/degrad + double(amrex::Real(50.439273))*psecdeg*years + double(amrex::Real(3.392506)) + mvsum;
417 
418  } // end of test on whether to calculate or use input orbital params
419 
420  // Restore any parameter the user specified explicitly, overriding whatever
421  // value was just computed for it
422  if (obliq_in >= zero) { obliq = obliq_in; }
423  if (mvelp_in >= zero) { mvelp = mvelp_in; }
424  if (eccen_in >= zero) {
425  eccen = eccen_in;
426  eccen2 = eccen*eccen;
427  eccen3 = eccen2*eccen;
428  }
429 
430  // Cases to make sure mvelp is between 0 and amrex::Real(360.)
431  if (mvelp < zero) {
432  do {
433  mvelp += amrex::Real(360.0);
434  } while (mvelp < zero);
435  }
436  if (mvelp >= amrex::Real(360.0)) {
437  do {
438  mvelp -= amrex::Real(360.0);
439  } while (mvelp >= amrex::Real(360.0));
440  }
441 
442  // Orbit needs the obliquity in radians
443 
444  obliqr = obliq*degrad;
445 
446  /*
447  180 degrees must be added to mvelp since observations are made from the
448  earth and the sun is considered (wrongly for the algorithm) to go around
449  the earth. For a more graphic explanation see Appendix B in:
450 
451  A. Berger, M. Loutre and C. Tricot. amrex::Real(1993.) Insolation and Earth Orbital
452  Periods. J. of Geophysical Research 98:10,341-10,amrex::Real(362.)
453 
454  Additionally, orbit will need this value in radians. So mvelp becomes
455  mvelpp (mvelp plus pi)
456  */
457 
458  mvelpp = (mvelp + amrex::Real(180.0))*degrad;
459 
460  // Set up an argument used several times in lambm0 calculation ahead.
461 
462  beta = std::sqrt(one - eccen2);
463 
464  /*
465  The mean longitude at the vernal equinox (lambda m nought in Berger
466  1978; in radians) is calculated from the following formula given in
467  Berger amrex::Real(1978.) At the vernal equinox the true longitude (lambda in Berger
468  1978) is zero
469  */
470 
471  lambm0 = two*( (amrex::Real(.5)*eccen + amrex::Real(.125)*eccen3)*(one + beta)*std::sin(mvelpp)
472  - amrex::Real(.250)*eccen2*(amrex::Real(.5) + beta)*std::sin(two*mvelpp)
473  + amrex::Real(.125)*eccen3*(one/three + beta)*std::sin(three*mvelpp) );
474 }
475 
476 
477 /**
478  * Day of the year plus fraction of the day, with calday 1.0 at 00:00 UTC on
479  * 1 January, from a UTC calendar date: the form orbital_decl takes. Leap
480  * years add a day after February, as in the RRTMGP interface.
481  */
482 AMREX_GPU_HOST
483 AMREX_FORCE_INLINE
484 double
485 orbital_calday (int year, int mon, int day, int sec)
486 {
487  static constexpr double dpy[] = {zero, double(amrex::Real(31.0)), double(amrex::Real(59.0)),
488  double(amrex::Real(90.0)), double(amrex::Real(120.0)), double(amrex::Real(151.0)),
489  double(amrex::Real(181.0)), double(amrex::Real(212.0)), double(amrex::Real(243.0)),
490  double(amrex::Real(273.0)), double(amrex::Real(304.0)), double(amrex::Real(334.0))};
491  const bool leap = (year % 4 == 0) && ((year % 100 != 0) || (year % 400 == 0));
492  double calday = one + dpy[mon-1] + double(day - 1) + double(sec) / double(amrex::Real(86400.0));
493  if (leap && mon > 2) { calday += one; }
494  return calday;
495 }
496 
497 /**
498  * Instantaneous cosine of the solar zenith angle at latitude lat and
499  * longitude lon (radians, east positive) at calday for declination declin
500  * (radians). This is the formula orbital_cos_zenith uses when it is not
501  * asked to average over an interval; it is callable from device code so a
502  * per-column kernel can evaluate it.
503  */
504 AMREX_GPU_HOST_DEVICE
505 AMREX_FORCE_INLINE
506 double
507 orbital_cos_zenith_instant (double jday, double lat, double lon, double declin)
508 {
509  return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
510  std::cos((jday - std::floor(jday))*double(two)*PI + lon);
511 }
512 
513 
514 AMREX_GPU_HOST
515 AMREX_FORCE_INLINE
516 double
517 orbital_avg_cos_zenith (double& jday,
518  double& lat,
519  double& lon,
520  double& declin,
521  double& dt_avg)
522 {
523  // adjust latitude so that its tangent will be defined
524  double del;
525  if (lat == PIoTwo) {
526  del = lat - double(amrex::Real(1.0e-05));
527  } else if (lat == -PIoTwo) {
528  del = lat + double(amrex::Real(1.0e-05));
529  } else {
530  del = lat;
531  }
532 
533  // adjust declination so that its tangent will be defined
534  double phi;
535  if (declin == PIoTwo) {
536  phi = declin - double(amrex::Real(1.0e-05));
537  } else if (declin == -PIoTwo) {
538  phi = declin + double(amrex::Real(1.0e-05));
539  } else {
540  phi = declin;
541  }
542 
543  // define the cosine of the myhalf-day length
544  // adjust for cases of all daylight or all night
545  double h;
546  double cos_h = - std::tan(del) * std::tan(phi);
547  if (cos_h <= -one) {
548  h = PI;
549  } else if (cos_h >= one) {
550  h = zero;
551  } else {
552  h = std::acos(cos_h);
553  }
554 
555  // Define Local Time t and t + dt
556  // adjust t to be between -pi and pi
557  double t1 = (jday - int(jday)) * two*PI + lon - PI;
558  if (t1 >= PI) {
559  t1 = t1 - two*PI;
560  } else if (t1 < -PI) {
561  t1 = t1 + two*PI;
562  }
563 
564  double dt = dt_avg / double(amrex::Real(86400.0)) * two*PI;
565  double t2 = t1 + dt;
566 
567  // Compute Cosine Solar Zenith angle
568  // define terms needed in the cosine zenith angle equation
569  double aa = std::sin(lat) * std::sin(declin);
570  double bb = std::cos(lat) * std::cos(declin);
571 
572  // define the hour angle
573  // force it to be between -h and h
574  // consider the situation when the night period is too short
575  double tt1,tt2,tt3,tt4;
576  if ( (t2 >= PI) && (t1 <= PI) && ((PI - h) <= dt) ) {
577  tt2 = h;
578  tt1 = std::min(std::max(t1, -h), h);
579  tt4 = std::min(std::max(t2, two*PI - h), two*PI + h);
580  tt3 = two*PI - h;
581  } else if ( (t2 >= -PI) && (t1 <= -PI) && ((PI - h) <= dt) ) {
582  tt2 = - two*PI + h;
583  tt1 = std::min(std::max(t1, -two*PI - h), -two*PI + h);
584  tt4 = std::min(std::max(t2, -h), h);
585  tt3 = -h;
586  } else {
587  if (t2 > PI) {
588  tt2 = std::min(std::max(t2 - two*PI, -h), h);
589  } else if (t2 < -PI) {
590  tt2 = std::min(std::max(t2 + two*PI, -h), h);
591  } else {
592  tt2 = std::min(std::max(t2 , -h), h);
593  }
594 
595  if (t1 > PI) {
596  tt1 = std::min(std::max(t1 - two*PI, -h), h);
597  } else if (t1 < -PI) {
598  tt1 = std::min(std::max(t1 + two*PI, -h), h);
599  } else {
600  tt1 = std::min(std::max(t1 , -h), h);
601  }
602  tt4 = zero;
603  tt3 = zero;
604  }
605 
606  // perform a time integration to obtain cosz if desired
607  // output is valid over the period from t to t + dt
608  if ( (tt2 > tt1) || (tt4 > tt3) ) {
609  return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
610  (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
611  } else {
612  return zero;
613  }
614 }
615 
616 
617 AMREX_GPU_HOST
618 AMREX_FORCE_INLINE
619 double
620 orbital_cos_zenith (double& jday,
621  double& lat,
622  double& lon,
623  double& declin,
624  double dt_avg = -one,
625  double uniform_angle = -one,
626  double constant_zenith_angle_deg = -one)
627 {
628  // Constant zenith angle is true
629  if ( constant_zenith_angle_deg >= zero ) {
630  return std::cos( constant_zenith_angle_deg * PI/amrex::Real(180.) );
631  }
632 
633  // Uniform angle is true
634  if ( uniform_angle >= zero) {
635  return std::cos(uniform_angle);
636  }
637 
638  // Perform the calculation of Orbital_Cos_Zenith
639  bool use_dt_avg = false;
640  if (dt_avg > zero) {
641  use_dt_avg = true;
642  }
643  // If dt for the average cosz is specified, then call the shr_orb_avg_cosz
644  if (use_dt_avg) {
645  return orbital_avg_cos_zenith(jday, lat, lon, declin, dt_avg);
646  } else {
647  return orbital_cos_zenith_instant(jday, lat, lon, declin);
648  }
649 }
650 #endif
ERF_Constants.H
ORB_UNDEF_INT
static constexpr int ORB_UNDEF_INT
Definition: ERF_Constants.H:68
beta
amrex::Real beta
Definition: ERF_InitCustomPert_DataAssimilation_ISV.H:10
three
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:32
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
one
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:30
fourth
constexpr amrex::Real fourth
Definition: ERF_NumericalConstants.H:35
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
PI
constexpr amrex::Real PI
Definition: ERF_NumericalConstants.H:39
PIoTwo
constexpr amrex::Real PIoTwo
Definition: ERF_NumericalConstants.H:40
orbital_avg_cos_zenith
AMREX_GPU_HOST AMREX_FORCE_INLINE double orbital_avg_cos_zenith(double &jday, double &lat, double &lon, double &declin, double &dt_avg)
Definition: ERF_OrbCosZenith.H:517
orbital_calday
AMREX_GPU_HOST AMREX_FORCE_INLINE double orbital_calday(int year, int mon, int day, int sec)
Definition: ERF_OrbCosZenith.H:485
orbital_cos_zenith_instant
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double orbital_cos_zenith_instant(double jday, double lat, double lon, double declin)
Definition: ERF_OrbCosZenith.H:507
orbital_cos_zenith
AMREX_GPU_HOST AMREX_FORCE_INLINE double orbital_cos_zenith(double &jday, double &lat, double &lon, double &declin, double dt_avg=-one, double uniform_angle=-one, double constant_zenith_angle_deg=-one)
Definition: ERF_OrbCosZenith.H:620
orbital_decl
AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_decl(double &calday, double &eccen, double &mvelpp, double &lambm0, double &obliqr, double &delta, double &eccf)
Definition: ERF_OrbCosZenith.H:18
orbital_params
AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_params(int &iyear_AD, double &eccen, double &obliq, double &mvelp, double &obliqr, double &lambm0, double &mvelpp)
Definition: ERF_OrbCosZenith.H:84
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