Line data Source code
1 : !=======================================================================
2 :
3 : ! Orbital parameters computed from date
4 : ! author: Bruce P. Briegleb, NCAR
5 : !
6 : ! 2006 ECH: Converted to free source form (F90)
7 : ! 2014 ECH: Moved routines from csm_share/shr_orb_mod.F90
8 :
9 : module icepack_orbital
10 :
11 : use icepack_kinds
12 : use icepack_parameters, only: c2, p5, pi, secday
13 : use icepack_warnings, only: warnstr, icepack_warnings_add
14 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
15 :
16 : implicit none
17 : private
18 : public :: compute_coszen, &
19 : icepack_init_orbit, &
20 : icepack_query_orbit
21 :
22 : ! orbital parameters
23 : integer (kind=int_kind) :: iyear_AD ! Year to calculate orbit for
24 : real(kind=dbl_kind) :: eccen ! Earth's orbital eccentricity
25 : real(kind=dbl_kind) :: obliqr ! Earth's obliquity in radians
26 : real(kind=dbl_kind) :: lambm0 ! Mean longitude of perihelion at the
27 : ! vernal equinox (radians)
28 : real(kind=dbl_kind) :: mvelpp ! Earth's moving vernal equinox longitude
29 : ! of perihelion + pi (radians)
30 : real(kind=dbl_kind) :: obliq ! obliquity in degrees
31 : real(kind=dbl_kind) :: mvelp ! moving vernal equinox long
32 : real(kind=dbl_kind) :: decln ! solar declination angle in radians
33 : real(kind=dbl_kind) :: eccf ! earth orbit eccentricity factor
34 : logical(kind=log_kind) :: log_print ! Flags print of status/error
35 :
36 : !=======================================================================
37 :
38 : contains
39 :
40 : !=======================================================================
41 :
42 : !autodocument_start icepack_init_orbit
43 : ! Compute orbital parameters for the specified date.
44 :
45 160 : subroutine icepack_init_orbit(iyear_AD_in, eccen_in, obliqr_in, &
46 : lambm0_in, mvelpp_in, obliq_in, mvelp_in, decln_in, eccf_in, &
47 : log_print_in)
48 :
49 : integer(kind=int_kind), optional, intent(in) :: iyear_AD_in ! Year to calculate orbit for
50 : real(kind=dbl_kind), optional, intent(in) :: eccen_in ! Earth's orbital eccentricity
51 : real(kind=dbl_kind), optional, intent(in) :: obliqr_in ! Earth's obliquity in radians
52 : real(kind=dbl_kind), optional, intent(in) :: lambm0_in ! Mean longitude of perihelion at the
53 : ! vernal equinox (radians)
54 : real(kind=dbl_kind), optional, intent(in) :: mvelpp_in ! Earth's moving vernal equinox longitude
55 : ! of perihelion + pi (radians)
56 : real(kind=dbl_kind), optional, intent(in) :: obliq_in ! obliquity in degrees
57 : real(kind=dbl_kind), optional, intent(in) :: mvelp_in ! moving vernal equinox long
58 : real(kind=dbl_kind), optional, intent(in) :: decln_in ! solar declination angle in radians
59 : real(kind=dbl_kind), optional, intent(in) :: eccf_in ! earth orbit eccentricity factor
60 : logical(kind=log_kind), optional, intent(in) :: log_print_in ! Flags print of status/error
61 :
62 : !autodocument_end
63 :
64 : character(len=*),parameter :: subname='(icepack_init_orbit)'
65 :
66 : !call icepack_warnings_add(subname//' ')
67 160 : iyear_AD = 1950
68 160 : log_print = .false. ! if true, write out orbital parameters
69 :
70 : #ifndef CESMCOUPLED
71 : call shr_orb_params( iyear_AD, eccen , obliq , mvelp , &
72 160 : obliqr , lambm0, mvelpp, log_print)
73 160 : if (icepack_warnings_aborted(subname)) return
74 : #endif
75 :
76 160 : if (present(iyear_AD_in)) iyear_AD = iyear_AD_in
77 160 : if (present(eccen_in)) eccen = eccen_in
78 160 : if (present(obliqr_in)) obliqr = obliqr_in
79 160 : if (present(lambm0_in)) lambm0 = lambm0_in
80 160 : if (present(mvelpp_in)) mvelpp = mvelpp_in
81 160 : if (present(obliq_in)) obliq = obliq_in
82 160 : if (present(mvelp_in)) mvelp = mvelp_in
83 160 : if (present(decln_in)) decln = decln_in
84 160 : if (present(eccf_in)) eccf = eccf_in
85 160 : if (present(log_print_in)) log_print = log_print_in
86 :
87 : end subroutine icepack_init_orbit
88 :
89 : !=======================================================================
90 :
91 : !autodocument_start icepack_query_orbit
92 : ! Compute orbital parameters for the specified date.
93 :
94 0 : subroutine icepack_query_orbit(iyear_AD_out, eccen_out, obliqr_out, &
95 : lambm0_out, mvelpp_out, obliq_out, mvelp_out, decln_out, eccf_out, &
96 : log_print_out)
97 :
98 : integer(kind=int_kind), optional, intent(out) :: iyear_AD_out ! Year to calculate orbit for
99 : real(kind=dbl_kind), optional, intent(out) :: eccen_out ! Earth's orbital eccentricity
100 : real(kind=dbl_kind), optional, intent(out) :: obliqr_out ! Earth's obliquity in radians
101 : real(kind=dbl_kind), optional, intent(out) :: lambm0_out ! Mean longitude of perihelion at the
102 : ! vernal equinox (radians)
103 : real(kind=dbl_kind), optional, intent(out) :: mvelpp_out ! Earth's moving vernal equinox longitude
104 : ! of perihelion + pi (radians)
105 : real(kind=dbl_kind), optional, intent(out) :: obliq_out ! obliquity in degrees
106 : real(kind=dbl_kind), optional, intent(out) :: mvelp_out ! moving vernal equinox long
107 : real(kind=dbl_kind), optional, intent(out) :: decln_out ! solar declination angle in radians
108 : real(kind=dbl_kind), optional, intent(out) :: eccf_out ! earth orbit eccentricity factor
109 : logical(kind=log_kind), optional, intent(out) :: log_print_out ! Flags print of status/error
110 :
111 : !autodocument_end
112 :
113 : character(len=*),parameter :: subname='(icepack_query_orbit)'
114 :
115 : !call icepack_warnings_add(subname//' ')
116 0 : iyear_AD = 1950
117 0 : log_print = .false. ! if true, write out orbital parameters
118 :
119 : #ifndef CESMCOUPLED
120 : call shr_orb_params( iyear_AD, eccen , obliq , mvelp , &
121 0 : obliqr , lambm0, mvelpp, log_print)
122 0 : if (icepack_warnings_aborted(subname)) return
123 : #endif
124 :
125 0 : if (present(iyear_AD_out)) iyear_AD_out = iyear_AD
126 0 : if (present(eccen_out)) eccen_out = eccen
127 0 : if (present(obliqr_out)) obliqr_out = obliqr
128 0 : if (present(lambm0_out)) lambm0_out = lambm0
129 0 : if (present(mvelpp_out)) mvelpp_out = mvelpp
130 0 : if (present(obliq_out)) obliq_out = obliq
131 0 : if (present(mvelp_out)) mvelp_out = mvelp
132 0 : if (present(decln_out)) decln_out = decln
133 0 : if (present(eccf_out)) eccf_out = eccf
134 0 : if (present(log_print_out)) log_print_out = log_print
135 :
136 : end subroutine icepack_query_orbit
137 :
138 : !=======================================================================
139 :
140 : ! Uses orbital and lat/lon info to compute cosine solar zenith angle
141 : ! for the specified date.
142 : !
143 : ! author: Bruce P. Briegleb, NCAR
144 :
145 909228 : subroutine compute_coszen (tlat, tlon, &
146 : yday, sec, coszen, &
147 : days_per_year, nextsw_cday, &
148 0 : calendar_type)
149 :
150 : #ifdef CESMCOUPLED
151 : use shr_orb_mod, only: shr_orb_decl
152 : #endif
153 :
154 : real (kind=dbl_kind), intent(in) :: &
155 : tlat, tlon ! latitude and longitude (radians)
156 :
157 : integer (kind=int_kind), intent(in) :: &
158 : sec ! elapsed seconds into date
159 :
160 : real (kind=dbl_kind), intent(in) :: &
161 : yday ! day of the year
162 :
163 : real (kind=dbl_kind), intent(inout) :: &
164 : coszen ! cosine solar zenith angle
165 : ! negative for sun below horizon
166 :
167 : integer (kind=int_kind), intent(in), optional :: &
168 : days_per_year ! number of days in one year
169 :
170 : real (kind=dbl_kind), intent(in), optional :: &
171 : nextsw_cday ! julian day of next shortwave calculation
172 :
173 : character (len=char_len), intent(in), optional :: &
174 : calendar_type ! differentiates Gregorian from other calendars
175 :
176 : ! local variables
177 :
178 329481 : real (kind=dbl_kind) :: ydayp1 ! day of year plus one time step
179 :
180 : character(len=*),parameter :: subname='(compute_coszen)'
181 :
182 : ! Solar declination for next time step
183 :
184 : #ifdef CESMCOUPLED
185 : if (calendar_type == "GREGORIAN") then
186 : ydayp1 = min(nextsw_cday, real(days_per_year,kind=dbl_kind))
187 : else
188 : ydayp1 = nextsw_cday
189 : endif
190 :
191 : !--- update coszen when nextsw_cday valid
192 : if (ydayp1 > -0.5_dbl_kind) then
193 : #else
194 909228 : ydayp1 = yday + sec/secday
195 : #endif
196 :
197 : call shr_orb_decl(ydayp1, eccen, mvelpp, lambm0, &
198 909228 : obliqr, decln, eccf)
199 909228 : if (icepack_warnings_aborted(subname)) return
200 :
201 : coszen = sin(tlat)*sin(decln) &
202 : + cos(tlat)*cos(decln) &
203 909228 : *cos((sec/secday-p5)*c2*pi + tlon) !cos(hour angle)
204 :
205 : #ifdef CESMCOUPLED
206 : endif
207 : #endif
208 :
209 : end subroutine compute_coszen
210 :
211 : !===============================================================================
212 :
213 : #ifndef CESMCOUPLED
214 160 : SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , &
215 : & obliqr , lambm0, mvelpp, log_print)
216 :
217 : !-------------------------------------------------------------------------------
218 : !
219 : ! Calculate earths orbital parameters using Dave Threshers formula which
220 : ! came from Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term
221 : ! Variations of Daily Insolation. Contribution 18, Institute of Astronomy
222 : ! and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium
223 : !
224 : !------------------------------Code history-------------------------------------
225 : !
226 : ! Original Author: Erik Kluzek
227 : ! Date: Oct/97
228 : !
229 : !-------------------------------------------------------------------------------
230 :
231 : !----------------------------- Arguments ------------------------------------
232 : integer(int_kind),intent(in) :: iyear_AD ! Year to calculate orbit for
233 : real (dbl_kind),intent(inout) :: eccen ! orbital eccentricity
234 : real (dbl_kind),intent(inout) :: obliq ! obliquity in degrees
235 : real (dbl_kind),intent(inout) :: mvelp ! moving vernal equinox long
236 : real (dbl_kind),intent(out) :: obliqr ! Earths obliquity in rad
237 : real (dbl_kind),intent(out) :: lambm0 ! Mean long of perihelion at
238 : ! vernal equinox (radians)
239 : real (dbl_kind),intent(out) :: mvelpp ! moving vernal equinox long
240 : ! of perihelion plus pi (rad)
241 : logical(log_kind),intent(in) :: log_print ! Flags print of status/error
242 :
243 : !------------------------------ Parameters ----------------------------------
244 : real (dbl_kind),parameter :: SHR_ORB_UNDEF_REAL = 1.e36_dbl_kind ! undefined real
245 : integer(int_kind),parameter :: SHR_ORB_UNDEF_INT = 2000000000 ! undefined int
246 : integer(int_kind),parameter :: poblen =47 ! # of elements in series wrt obliquity
247 : integer(int_kind),parameter :: pecclen=19 ! # of elements in series wrt eccentricity
248 : integer(int_kind),parameter :: pmvelen=78 ! # of elements in series wrt vernal equinox
249 : real (dbl_kind),parameter :: psecdeg = 1.0_dbl_kind/3600.0_dbl_kind ! arc sec to deg conversion
250 :
251 60 : real (dbl_kind) :: yb4_1950AD ! number of years before 1950 AD
252 :
253 : real (dbl_kind),parameter :: SHR_ORB_ECCEN_MIN = 0.0_dbl_kind ! min value for eccen
254 : real (dbl_kind),parameter :: SHR_ORB_ECCEN_MAX = 0.1_dbl_kind ! max value for eccen
255 : real (dbl_kind),parameter :: SHR_ORB_OBLIQ_MIN = -90.0_dbl_kind ! min value for obliq
256 : real (dbl_kind),parameter :: SHR_ORB_OBLIQ_MAX = +90.0_dbl_kind ! max value for obliq
257 : real (dbl_kind),parameter :: SHR_ORB_MVELP_MIN = 0.0_dbl_kind ! min value for mvelp
258 : real (dbl_kind),parameter :: SHR_ORB_MVELP_MAX = 360.0_dbl_kind ! max value for mvelp
259 :
260 : ! Cosine series data for computation of obliquity: amplitude (arc seconds),
261 : ! rate (arc seconds/year), phase (degrees).
262 :
263 : real (dbl_kind), parameter :: obamp(poblen) = & ! amplitudes for obliquity cos series
264 : & (/ -2462.2214466_dbl_kind, -857.3232075_dbl_kind, -629.3231835_dbl_kind, &
265 : & -414.2804924_dbl_kind, -311.7632587_dbl_kind, 308.9408604_dbl_kind, &
266 : & -162.5533601_dbl_kind, -116.1077911_dbl_kind, 101.1189923_dbl_kind, &
267 : & -67.6856209_dbl_kind, 24.9079067_dbl_kind, 22.5811241_dbl_kind, &
268 : & -21.1648355_dbl_kind, -15.6549876_dbl_kind, 15.3936813_dbl_kind, &
269 : & 14.6660938_dbl_kind, -11.7273029_dbl_kind, 10.2742696_dbl_kind, &
270 : & 6.4914588_dbl_kind, 5.8539148_dbl_kind, -5.4872205_dbl_kind, &
271 : & -5.4290191_dbl_kind, 5.1609570_dbl_kind, 5.0786314_dbl_kind, &
272 : & -4.0735782_dbl_kind, 3.7227167_dbl_kind, 3.3971932_dbl_kind, &
273 : & -2.8347004_dbl_kind, -2.6550721_dbl_kind, -2.5717867_dbl_kind, &
274 : & -2.4712188_dbl_kind, 2.4625410_dbl_kind, 2.2464112_dbl_kind, &
275 : & -2.0755511_dbl_kind, -1.9713669_dbl_kind, -1.8813061_dbl_kind, &
276 : & -1.8468785_dbl_kind, 1.8186742_dbl_kind, 1.7601888_dbl_kind, &
277 : & -1.5428851_dbl_kind, 1.4738838_dbl_kind, -1.4593669_dbl_kind, &
278 : & 1.4192259_dbl_kind, -1.1818980_dbl_kind, 1.1756474_dbl_kind, &
279 : & -1.1316126_dbl_kind, 1.0896928_dbl_kind/)
280 :
281 : real (dbl_kind), parameter :: obrate(poblen) = & ! rates for obliquity cosine series
282 : & (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, &
283 : & 31.983787_dbl_kind, 44.828336_dbl_kind, 30.973257_dbl_kind, &
284 : & 43.668246_dbl_kind, 32.246691_dbl_kind, 30.599444_dbl_kind, &
285 : & 42.681324_dbl_kind, 43.836462_dbl_kind, 47.439436_dbl_kind, &
286 : & 63.219948_dbl_kind, 64.230478_dbl_kind, 1.010530_dbl_kind, &
287 : & 7.437771_dbl_kind, 55.782177_dbl_kind, 0.373813_dbl_kind, &
288 : & 13.218362_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, &
289 : & 76.438310_dbl_kind, 45.815258_dbl_kind, 8.448301_dbl_kind, &
290 : & 56.792707_dbl_kind, 49.747842_dbl_kind, 12.058272_dbl_kind, &
291 : & 75.278220_dbl_kind, 65.241008_dbl_kind, 64.604291_dbl_kind, &
292 : & 1.647247_dbl_kind, 7.811584_dbl_kind, 12.207832_dbl_kind, &
293 : & 63.856665_dbl_kind, 56.155990_dbl_kind, 77.448840_dbl_kind, &
294 : & 6.801054_dbl_kind, 62.209418_dbl_kind, 20.656133_dbl_kind, &
295 : & 48.344406_dbl_kind, 55.145460_dbl_kind, 69.000539_dbl_kind, &
296 : & 11.071350_dbl_kind, 74.291298_dbl_kind, 11.047742_dbl_kind, &
297 : & 0.636717_dbl_kind, 12.844549_dbl_kind/)
298 :
299 : real (dbl_kind), parameter :: obphas(poblen) = & ! phases for obliquity cosine series
300 : & (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, &
301 : & 292.7252_dbl_kind, 15.3747_dbl_kind, 263.7951_dbl_kind, &
302 : & 308.4258_dbl_kind, 240.0099_dbl_kind, 222.9725_dbl_kind, &
303 : & 268.7809_dbl_kind, 316.7998_dbl_kind, 319.6024_dbl_kind, &
304 : & 143.8050_dbl_kind, 172.7351_dbl_kind, 28.9300_dbl_kind, &
305 : & 123.5968_dbl_kind, 20.2082_dbl_kind, 40.8226_dbl_kind, &
306 : & 123.4722_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, &
307 : & 267.2772_dbl_kind, 55.0196_dbl_kind, 152.5268_dbl_kind, &
308 : & 49.1382_dbl_kind, 204.6609_dbl_kind, 56.5233_dbl_kind, &
309 : & 200.3284_dbl_kind, 201.6651_dbl_kind, 213.5577_dbl_kind, &
310 : & 17.0374_dbl_kind, 164.4194_dbl_kind, 94.5422_dbl_kind, &
311 : & 131.9124_dbl_kind, 61.0309_dbl_kind, 296.2073_dbl_kind, &
312 : & 135.4894_dbl_kind, 114.8750_dbl_kind, 247.0691_dbl_kind, &
313 : & 256.6114_dbl_kind, 32.1008_dbl_kind, 143.6804_dbl_kind, &
314 : & 16.8784_dbl_kind, 160.6835_dbl_kind, 27.5932_dbl_kind, &
315 : & 348.1074_dbl_kind, 82.6496_dbl_kind/)
316 :
317 : ! Cosine/sine series data for computation of eccentricity and fixed vernal
318 : ! equinox longitude of perihelion (fvelp): amplitude,
319 : ! rate (arc seconds/year), phase (degrees).
320 :
321 : real (dbl_kind), parameter :: ecamp (pecclen) = & ! ampl for eccen/fvelp cos/sin series
322 : & (/ 0.01860798_dbl_kind, 0.01627522_dbl_kind, -0.01300660_dbl_kind, &
323 : & 0.00988829_dbl_kind, -0.00336700_dbl_kind, 0.00333077_dbl_kind, &
324 : & -0.00235400_dbl_kind, 0.00140015_dbl_kind, 0.00100700_dbl_kind, &
325 : & 0.00085700_dbl_kind, 0.00064990_dbl_kind, 0.00059900_dbl_kind, &
326 : & 0.00037800_dbl_kind, -0.00033700_dbl_kind, 0.00027600_dbl_kind, &
327 : & 0.00018200_dbl_kind, -0.00017400_dbl_kind, -0.00012400_dbl_kind, &
328 : & 0.00001250_dbl_kind/)
329 :
330 : real (dbl_kind), parameter :: ecrate(pecclen) = & ! rates for eccen/fvelp cos/sin series
331 : & (/ 4.2072050_dbl_kind, 7.3460910_dbl_kind, 17.8572630_dbl_kind, &
332 : & 17.2205460_dbl_kind, 16.8467330_dbl_kind, 5.1990790_dbl_kind, &
333 : & 18.2310760_dbl_kind, 26.2167580_dbl_kind, 6.3591690_dbl_kind, &
334 : & 16.2100160_dbl_kind, 3.0651810_dbl_kind, 16.5838290_dbl_kind, &
335 : & 18.4939800_dbl_kind, 6.1909530_dbl_kind, 18.8677930_dbl_kind, &
336 : & 17.4255670_dbl_kind, 6.1860010_dbl_kind, 18.4174410_dbl_kind, &
337 : & 0.6678630_dbl_kind/)
338 :
339 : real (dbl_kind), parameter :: ecphas(pecclen) = & ! phases for eccen/fvelp cos/sin series
340 : & (/ 28.620089_dbl_kind, 193.788772_dbl_kind, 308.307024_dbl_kind, &
341 : & 320.199637_dbl_kind, 279.376984_dbl_kind, 87.195000_dbl_kind, &
342 : & 349.129677_dbl_kind, 128.443387_dbl_kind, 154.143880_dbl_kind, &
343 : & 291.269597_dbl_kind, 114.860583_dbl_kind, 332.092251_dbl_kind, &
344 : & 296.414411_dbl_kind, 145.769910_dbl_kind, 337.237063_dbl_kind, &
345 : & 152.092288_dbl_kind, 126.839891_dbl_kind, 210.667199_dbl_kind, &
346 : & 72.108838_dbl_kind/)
347 :
348 : ! Sine series data for computation of moving vernal equinox longitude of
349 : ! perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees).
350 :
351 : real (dbl_kind), parameter :: mvamp (pmvelen) = & ! amplitudes for mvelp sine series
352 : & (/ 7391.0225890_dbl_kind, 2555.1526947_dbl_kind, 2022.7629188_dbl_kind, &
353 : & -1973.6517951_dbl_kind, 1240.2321818_dbl_kind, 953.8679112_dbl_kind, &
354 : & -931.7537108_dbl_kind, 872.3795383_dbl_kind, 606.3544732_dbl_kind, &
355 : & -496.0274038_dbl_kind, 456.9608039_dbl_kind, 346.9462320_dbl_kind, &
356 : & -305.8412902_dbl_kind, 249.6173246_dbl_kind, -199.1027200_dbl_kind, &
357 : & 191.0560889_dbl_kind, -175.2936572_dbl_kind, 165.9068833_dbl_kind, &
358 : & 161.1285917_dbl_kind, 139.7878093_dbl_kind, -133.5228399_dbl_kind, &
359 : & 117.0673811_dbl_kind, 104.6907281_dbl_kind, 95.3227476_dbl_kind, &
360 : & 86.7824524_dbl_kind, 86.0857729_dbl_kind, 70.5893698_dbl_kind, &
361 : & -69.9719343_dbl_kind, -62.5817473_dbl_kind, 61.5450059_dbl_kind, &
362 : & -57.9364011_dbl_kind, 57.1899832_dbl_kind, -57.0236109_dbl_kind, &
363 : & -54.2119253_dbl_kind, 53.2834147_dbl_kind, 52.1223575_dbl_kind, &
364 : & -49.0059908_dbl_kind, -48.3118757_dbl_kind, -45.4191685_dbl_kind, &
365 : & -42.2357920_dbl_kind, -34.7971099_dbl_kind, 34.4623613_dbl_kind, &
366 : & -33.8356643_dbl_kind, 33.6689362_dbl_kind, -31.2521586_dbl_kind, &
367 : & -30.8798701_dbl_kind, 28.4640769_dbl_kind, -27.1960802_dbl_kind, &
368 : & 27.0860736_dbl_kind, -26.3437456_dbl_kind, 24.7253740_dbl_kind, &
369 : & 24.6732126_dbl_kind, 24.4272733_dbl_kind, 24.0127327_dbl_kind, &
370 : & 21.7150294_dbl_kind, -21.5375347_dbl_kind, 18.1148363_dbl_kind, &
371 : & -16.9603104_dbl_kind, -16.1765215_dbl_kind, 15.5567653_dbl_kind, &
372 : & 15.4846529_dbl_kind, 15.2150632_dbl_kind, 14.5047426_dbl_kind, &
373 : & -14.3873316_dbl_kind, 13.1351419_dbl_kind, 12.8776311_dbl_kind, &
374 : & 11.9867234_dbl_kind, 11.9385578_dbl_kind, 11.7030822_dbl_kind, &
375 : & 11.6018181_dbl_kind, -11.2617293_dbl_kind, -10.4664199_dbl_kind, &
376 : & 10.4333970_dbl_kind, -10.2377466_dbl_kind, 10.1934446_dbl_kind, &
377 : & -10.1280191_dbl_kind, 10.0289441_dbl_kind, -10.0034259_dbl_kind/)
378 :
379 : real (dbl_kind), parameter :: mvrate(pmvelen) = & ! rates for mvelp sine series
380 : & (/ 31.609974_dbl_kind, 32.620504_dbl_kind, 24.172203_dbl_kind, &
381 : & 0.636717_dbl_kind, 31.983787_dbl_kind, 3.138886_dbl_kind, &
382 : & 30.973257_dbl_kind, 44.828336_dbl_kind, 0.991874_dbl_kind, &
383 : & 0.373813_dbl_kind, 43.668246_dbl_kind, 32.246691_dbl_kind, &
384 : & 30.599444_dbl_kind, 2.147012_dbl_kind, 10.511172_dbl_kind, &
385 : & 42.681324_dbl_kind, 13.650058_dbl_kind, 0.986922_dbl_kind, &
386 : & 9.874455_dbl_kind, 13.013341_dbl_kind, 0.262904_dbl_kind, &
387 : & 0.004952_dbl_kind, 1.142024_dbl_kind, 63.219948_dbl_kind, &
388 : & 0.205021_dbl_kind, 2.151964_dbl_kind, 64.230478_dbl_kind, &
389 : & 43.836462_dbl_kind, 47.439436_dbl_kind, 1.384343_dbl_kind, &
390 : & 7.437771_dbl_kind, 18.829299_dbl_kind, 9.500642_dbl_kind, &
391 : & 0.431696_dbl_kind, 1.160090_dbl_kind, 55.782177_dbl_kind, &
392 : & 12.639528_dbl_kind, 1.155138_dbl_kind, 0.168216_dbl_kind, &
393 : & 1.647247_dbl_kind, 10.884985_dbl_kind, 5.610937_dbl_kind, &
394 : & 12.658184_dbl_kind, 1.010530_dbl_kind, 1.983748_dbl_kind, &
395 : & 14.023871_dbl_kind, 0.560178_dbl_kind, 1.273434_dbl_kind, &
396 : & 12.021467_dbl_kind, 62.583231_dbl_kind, 63.593761_dbl_kind, &
397 : & 76.438310_dbl_kind, 4.280910_dbl_kind, 13.218362_dbl_kind, &
398 : & 17.818769_dbl_kind, 8.359495_dbl_kind, 56.792707_dbl_kind, &
399 : & 8.448301_dbl_kind, 1.978796_dbl_kind, 8.863925_dbl_kind, &
400 : & 0.186365_dbl_kind, 8.996212_dbl_kind, 6.771027_dbl_kind, &
401 : & 45.815258_dbl_kind, 12.002811_dbl_kind, 75.278220_dbl_kind, &
402 : & 65.241008_dbl_kind, 18.870667_dbl_kind, 22.009553_dbl_kind, &
403 : & 64.604291_dbl_kind, 11.498094_dbl_kind, 0.578834_dbl_kind, &
404 : & 9.237738_dbl_kind, 49.747842_dbl_kind, 2.147012_dbl_kind, &
405 : & 1.196895_dbl_kind, 2.133898_dbl_kind, 0.173168_dbl_kind/)
406 :
407 : real (dbl_kind), parameter :: mvphas(pmvelen) = & ! phases for mvelp sine series
408 : & (/ 251.9025_dbl_kind, 280.8325_dbl_kind, 128.3057_dbl_kind, &
409 : & 348.1074_dbl_kind, 292.7252_dbl_kind, 165.1686_dbl_kind, &
410 : & 263.7951_dbl_kind, 15.3747_dbl_kind, 58.5749_dbl_kind, &
411 : & 40.8226_dbl_kind, 308.4258_dbl_kind, 240.0099_dbl_kind, &
412 : & 222.9725_dbl_kind, 106.5937_dbl_kind, 114.5182_dbl_kind, &
413 : & 268.7809_dbl_kind, 279.6869_dbl_kind, 39.6448_dbl_kind, &
414 : & 126.4108_dbl_kind, 291.5795_dbl_kind, 307.2848_dbl_kind, &
415 : & 18.9300_dbl_kind, 273.7596_dbl_kind, 143.8050_dbl_kind, &
416 : & 191.8927_dbl_kind, 125.5237_dbl_kind, 172.7351_dbl_kind, &
417 : & 316.7998_dbl_kind, 319.6024_dbl_kind, 69.7526_dbl_kind, &
418 : & 123.5968_dbl_kind, 217.6432_dbl_kind, 85.5882_dbl_kind, &
419 : & 156.2147_dbl_kind, 66.9489_dbl_kind, 20.2082_dbl_kind, &
420 : & 250.7568_dbl_kind, 48.0188_dbl_kind, 8.3739_dbl_kind, &
421 : & 17.0374_dbl_kind, 155.3409_dbl_kind, 94.1709_dbl_kind, &
422 : & 221.1120_dbl_kind, 28.9300_dbl_kind, 117.1498_dbl_kind, &
423 : & 320.5095_dbl_kind, 262.3602_dbl_kind, 336.2148_dbl_kind, &
424 : & 233.0046_dbl_kind, 155.6977_dbl_kind, 184.6277_dbl_kind, &
425 : & 267.2772_dbl_kind, 78.9281_dbl_kind, 123.4722_dbl_kind, &
426 : & 188.7132_dbl_kind, 180.1364_dbl_kind, 49.1382_dbl_kind, &
427 : & 152.5268_dbl_kind, 98.2198_dbl_kind, 97.4808_dbl_kind, &
428 : & 221.5376_dbl_kind, 168.2438_dbl_kind, 161.1199_dbl_kind, &
429 : & 55.0196_dbl_kind, 262.6495_dbl_kind, 200.3284_dbl_kind, &
430 : & 201.6651_dbl_kind, 294.6547_dbl_kind, 99.8233_dbl_kind, &
431 : & 213.5577_dbl_kind, 154.1631_dbl_kind, 232.7153_dbl_kind, &
432 : & 138.3034_dbl_kind, 204.6609_dbl_kind, 106.5938_dbl_kind, &
433 : & 250.4676_dbl_kind, 332.3345_dbl_kind, 27.3039_dbl_kind/)
434 :
435 : !---------------------------Local variables----------------------------------
436 : integer(int_kind) :: i ! Index for series summations
437 60 : real (dbl_kind) :: obsum ! Obliquity series summation
438 60 : real (dbl_kind) :: cossum ! Cos series summation for eccentricity/fvelp
439 60 : real (dbl_kind) :: sinsum ! Sin series summation for eccentricity/fvelp
440 60 : real (dbl_kind) :: fvelp ! Fixed vernal equinox long of perihelion
441 60 : real (dbl_kind) :: mvsum ! mvelp series summation
442 60 : real (dbl_kind) :: beta ! Intermediate argument for lambm0
443 60 : real (dbl_kind) :: years ! Years to time of interest ( pos <=> future)
444 60 : real (dbl_kind) :: eccen2 ! eccentricity squared
445 60 : real (dbl_kind) :: eccen3 ! eccentricity cubed
446 60 : real (dbl_kind) :: degrad ! degrees to rad conversion
447 : integer (int_kind), parameter :: s_loglev = 0
448 : character(len=*),parameter :: subname='(shr_orb_params)'
449 :
450 : !-------------------------- Formats -----------------------------------------
451 : character(len=*),parameter :: F00 = "('(shr_orb_params) ',4a)"
452 : character(len=*),parameter :: F01 = "('(shr_orb_params) ',a,i9)"
453 : character(len=*),parameter :: F02 = "('(shr_orb_params) ',a,f6.3)"
454 : character(len=*),parameter :: F03 = "('(shr_orb_params) ',a,es14.6)"
455 :
456 : !----------------------------------------------------------------------------
457 : ! radinp and algorithms below will need a degree to radian conversion factor
458 :
459 : !call icepack_warnings_add(subname//' ')
460 160 : degrad = pi/180._dbl_kind ! degree to radian conversion factor
461 :
462 : if ( log_print .and. s_loglev > 0 ) then
463 : write(warnstr,F00) subname//'Calculate characteristics of the orbit:'
464 : call icepack_warnings_add(warnstr)
465 : end if
466 :
467 : ! Check for flag to use input orbit parameters
468 :
469 160 : IF ( iyear_AD == SHR_ORB_UNDEF_INT ) THEN
470 :
471 : ! Check input obliq, eccen, and mvelp to ensure reasonable
472 :
473 0 : if( obliq == SHR_ORB_UNDEF_REAL )then
474 0 : write(warnstr,F00) subname//' Have to specify orbital parameters:'
475 0 : call icepack_warnings_add(warnstr)
476 0 : write(warnstr,F00) subname//'Either set: iyear_AD, OR [obliq, eccen, and mvelp]:'
477 0 : call icepack_warnings_add(warnstr)
478 0 : write(warnstr,F00) subname//'iyear_AD is the year to simulate orbit for (ie. 1950): '
479 0 : call icepack_warnings_add(warnstr)
480 0 : write(warnstr,F00) subname//'obliq, eccen, mvelp specify the orbit directly:'
481 0 : call icepack_warnings_add(warnstr)
482 0 : write(warnstr,F00) subname//'The AMIP II settings (for a 1995 orbit) are: '
483 0 : call icepack_warnings_add(warnstr)
484 0 : write(warnstr,F00) subname//' obliq = 23.4441'
485 0 : call icepack_warnings_add(warnstr)
486 0 : write(warnstr,F00) subname//' eccen = 0.016715'
487 0 : call icepack_warnings_add(warnstr)
488 0 : write(warnstr,F00) subname//' mvelp = 102.7'
489 0 : call icepack_warnings_add(warnstr)
490 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
491 0 : call icepack_warnings_add(subname//' unreasonable oblip')
492 0 : else if ( log_print ) then
493 0 : write(warnstr,F00) subname//'Use input orbital parameters: '
494 0 : call icepack_warnings_add(warnstr)
495 : end if
496 0 : if( (obliq < SHR_ORB_OBLIQ_MIN).or.(obliq > SHR_ORB_OBLIQ_MAX) ) then
497 0 : write(warnstr,F03) subname//'Input obliquity unreasonable: ', obliq
498 0 : call icepack_warnings_add(warnstr)
499 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
500 0 : call icepack_warnings_add(subname//' unreasonable obliq')
501 : end if
502 0 : if( (eccen < SHR_ORB_ECCEN_MIN).or.(eccen > SHR_ORB_ECCEN_MAX) ) then
503 0 : write(warnstr,F03) subname//'Input eccentricity unreasonable: ', eccen
504 0 : call icepack_warnings_add(warnstr)
505 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
506 0 : call icepack_warnings_add(subname//' unreasonable eccen')
507 : end if
508 0 : if( (mvelp < SHR_ORB_MVELP_MIN).or.(mvelp > SHR_ORB_MVELP_MAX) ) then
509 0 : write(warnstr,F03) subname//'Input mvelp unreasonable: ' , mvelp
510 0 : call icepack_warnings_add(warnstr)
511 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
512 0 : call icepack_warnings_add(subname//' unreasonable mvelp')
513 : end if
514 0 : eccen2 = eccen*eccen
515 0 : eccen3 = eccen2*eccen
516 :
517 : ELSE ! Otherwise calculate based on years before present
518 :
519 : if ( log_print .and. s_loglev > 0) then
520 : write(warnstr,F01) subname//'Calculate orbit for year: ' , iyear_AD
521 : call icepack_warnings_add(warnstr)
522 : end if
523 160 : yb4_1950AD = 1950.0_dbl_kind - real(iyear_AD,dbl_kind)
524 160 : if ( abs(yb4_1950AD) .gt. 1000000.0_dbl_kind )then
525 0 : write(warnstr,F00) subname//'orbit only valid for years+-1000000'
526 0 : call icepack_warnings_add(warnstr)
527 0 : write(warnstr,F00) subname//'Relative to 1950 AD'
528 0 : call icepack_warnings_add(warnstr)
529 0 : write(warnstr,F03) subname//'# of years before 1950: ',yb4_1950AD
530 0 : call icepack_warnings_add(warnstr)
531 0 : write(warnstr,F01) subname//'Year to simulate was : ',iyear_AD
532 0 : call icepack_warnings_add(warnstr)
533 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
534 0 : call icepack_warnings_add(subname//' unreasonable year')
535 : end if
536 :
537 : ! The following calculates the earths obliquity, orbital eccentricity
538 : ! (and various powers of it) and vernal equinox mean longitude of
539 : ! perihelion for years in the past (future = negative of years past),
540 : ! using constants (see parameter section) given in the program of:
541 : !
542 : ! Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations
543 : ! of Daily Insolation. Contribution 18, Institute of Astronomy and
544 : ! Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
545 : !
546 : ! and formulas given in the paper (where less precise constants are also
547 : ! given):
548 : !
549 : ! Berger, Andre. 1978. Long-Term Variations of Daily Insolation and
550 : ! Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367
551 : !
552 : ! The algorithm is valid only to 1,000,000 years past or hence.
553 : ! For a solution valid to 5-10 million years past see the above author.
554 : ! Algorithm below is better for years closer to present than is the
555 : ! 5-10 million year solution.
556 : !
557 : ! Years to time of interest must be negative of years before present
558 : ! (1950) in formulas that follow.
559 :
560 160 : years = - yb4_1950AD
561 :
562 : ! In the summations below, cosine or sine arguments, which end up in
563 : ! degrees, must be converted to radians via multiplication by degrad.
564 : !
565 : ! Summation of cosine series for obliquity (epsilon in Berger 1978) in
566 : ! degrees. Convert the amplitudes and rates, which are in arc secs, into
567 : ! degrees via multiplication by psecdeg (arc seconds to degrees conversion
568 : ! factor). For obliq, first term is Berger 1978 epsilon star; second
569 : ! term is series summation in degrees.
570 :
571 160 : obsum = 0.0_dbl_kind
572 7680 : do i = 1, poblen
573 2820 : obsum = obsum + obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + &
574 7680 : & obphas(i))*degrad)
575 : end do
576 160 : obliq = 23.320556_dbl_kind + obsum
577 :
578 : ! Summation of cosine and sine series for computation of eccentricity
579 : ! (eccen; e in Berger 1978) and fixed vernal equinox longitude of
580 : ! perihelion (fvelp; pi in Berger 1978), which is used for computation
581 : ! of moving vernal equinox longitude of perihelion. Convert the rates,
582 : ! which are in arc seconds, into degrees via multiplication by psecdeg.
583 :
584 160 : cossum = 0.0_dbl_kind
585 3200 : do i = 1, pecclen
586 3200 : cossum = cossum+ecamp(i)*cos((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
587 : end do
588 :
589 160 : sinsum = 0.0_dbl_kind
590 3200 : do i = 1, pecclen
591 3200 : sinsum = sinsum+ecamp(i)*sin((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
592 : end do
593 :
594 : ! Use summations to calculate eccentricity
595 :
596 160 : eccen2 = cossum*cossum + sinsum*sinsum
597 160 : eccen = sqrt(eccen2)
598 160 : eccen3 = eccen2*eccen
599 :
600 : ! A series of cases for fvelp, which is in radians.
601 :
602 160 : if (abs(cossum) .le. 1.0E-8_dbl_kind) then
603 0 : if (sinsum .eq. 0.0_dbl_kind) then
604 0 : fvelp = 0.0_dbl_kind
605 0 : else if (sinsum .lt. 0.0_dbl_kind) then
606 0 : fvelp = 1.5_dbl_kind*pi
607 0 : else if (sinsum .gt. 0.0_dbl_kind) then
608 0 : fvelp = .5_dbl_kind*pi
609 : endif
610 160 : else if (cossum .lt. 0.0_dbl_kind) then
611 160 : fvelp = atan(sinsum/cossum) + pi
612 0 : else if (cossum .gt. 0.0_dbl_kind) then
613 0 : if (sinsum .lt. 0.0_dbl_kind) then
614 0 : fvelp = atan(sinsum/cossum) + 2.0_dbl_kind*pi
615 : else
616 0 : fvelp = atan(sinsum/cossum)
617 : endif
618 : endif
619 :
620 : ! Summation of sin series for computation of moving vernal equinox long
621 : ! of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp,
622 : ! first term is fvelp in degrees; second term is Berger 1978 psi bar
623 : ! times years and in degrees; third term is Berger 1978 zeta; fourth
624 : ! term is series summation in degrees. Convert the amplitudes and rates,
625 : ! which are in arc seconds, into degrees via multiplication by psecdeg.
626 : ! Series summation plus second and third terms constitute Berger 1978
627 : ! psi, which is the general precession.
628 :
629 160 : mvsum = 0.0_dbl_kind
630 12640 : do i = 1, pmvelen
631 4680 : mvsum = mvsum + mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + &
632 12640 : & mvphas(i))*degrad)
633 : end do
634 160 : mvelp = fvelp/degrad + 50.439273_dbl_kind*psecdeg*years + 3.392506_dbl_kind + mvsum
635 :
636 : ! Cases to make sure mvelp is between 0 and 360.
637 :
638 160 : do while (mvelp .lt. 0.0_dbl_kind)
639 160 : mvelp = mvelp + 360.0_dbl_kind
640 : end do
641 160 : do while (mvelp .ge. 360.0_dbl_kind)
642 160 : mvelp = mvelp - 360.0_dbl_kind
643 : end do
644 :
645 : END IF ! end of test on whether to calculate or use input orbital params
646 :
647 : ! Orbit needs the obliquity in radians
648 :
649 160 : obliqr = obliq*degrad
650 :
651 : ! 180 degrees must be added to mvelp since observations are made from the
652 : ! earth and the sun is considered (wrongly for the algorithm) to go around
653 : ! the earth. For a more graphic explanation see Appendix B in:
654 : !
655 : ! A. Berger, M. Loutre and C. Tricot. 1993. Insolation and Earth Orbital
656 : ! Periods. J. of Geophysical Research 98:10,341-10,362.
657 : !
658 : ! Additionally, orbit will need this value in radians. So mvelp becomes
659 : ! mvelpp (mvelp plus pi)
660 :
661 160 : mvelpp = (mvelp + 180._dbl_kind)*degrad
662 :
663 : ! Set up an argument used several times in lambm0 calculation ahead.
664 :
665 160 : beta = sqrt(1._dbl_kind - eccen2)
666 :
667 : ! The mean longitude at the vernal equinox (lambda m nought in Berger
668 : ! 1978; in radians) is calculated from the following formula given in
669 : ! Berger 1978. At the vernal equinox the true longitude (lambda in Berger
670 : ! 1978) is 0.
671 :
672 : lambm0 = 2._dbl_kind*((.5_dbl_kind*eccen + .125_dbl_kind*eccen3)*(1._dbl_kind + beta)*sin(mvelpp) &
673 : & - .250_dbl_kind*eccen2*(.5_dbl_kind + beta)*sin(2._dbl_kind*mvelpp) &
674 160 : & + .125_dbl_kind*eccen3*(1._dbl_kind/3._dbl_kind + beta)*sin(3._dbl_kind*mvelpp))
675 :
676 160 : if ( log_print ) then
677 0 : write(warnstr,F03) subname//'------ Computed Orbital Parameters ------'
678 0 : call icepack_warnings_add(warnstr)
679 0 : write(warnstr,F03) subname//'Eccentricity = ',eccen
680 0 : call icepack_warnings_add(warnstr)
681 0 : write(warnstr,F03) subname//'Obliquity (deg) = ',obliq
682 0 : call icepack_warnings_add(warnstr)
683 0 : write(warnstr,F03) subname//'Obliquity (rad) = ',obliqr
684 0 : call icepack_warnings_add(warnstr)
685 0 : write(warnstr,F03) subname//'Long of perh(deg) = ',mvelp
686 0 : call icepack_warnings_add(warnstr)
687 0 : write(warnstr,F03) subname//'Long of perh(rad) = ',mvelpp
688 0 : call icepack_warnings_add(warnstr)
689 0 : write(warnstr,F03) subname//'Long at v.e.(rad) = ',lambm0
690 0 : call icepack_warnings_add(warnstr)
691 0 : write(warnstr,F03) subname//'-----------------------------------------'
692 0 : call icepack_warnings_add(warnstr)
693 : end if
694 :
695 329581 : END SUBROUTINE shr_orb_params
696 :
697 : !===============================================================================
698 :
699 909228 : SUBROUTINE shr_orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr ,delta ,eccf)
700 :
701 : !-------------------------------------------------------------------------------
702 : !
703 : ! Compute earth/orbit parameters using formula suggested by
704 : ! Duane Thresher.
705 : !
706 : !---------------------------Code history----------------------------------------
707 : !
708 : ! Original version: Erik Kluzek
709 : ! Date: Oct/1997
710 : !
711 : !-------------------------------------------------------------------------------
712 :
713 : !------------------------------Arguments--------------------------------
714 : real (dbl_kind),intent(in) :: calday ! Calendar day, including fraction
715 : real (dbl_kind),intent(in) :: eccen ! Eccentricity
716 : real (dbl_kind),intent(in) :: obliqr ! Earths obliquity in radians
717 : real (dbl_kind),intent(in) :: lambm0 ! Mean long of perihelion at the
718 : ! vernal equinox (radians)
719 : real (dbl_kind),intent(in) :: mvelpp ! moving vernal equinox longitude
720 : ! of perihelion plus pi (radians)
721 : real (dbl_kind),intent(out) :: delta ! Solar declination angle in rad
722 : real (dbl_kind),intent(out) :: eccf ! Earth-sun distance factor (ie. (1/r)**2)
723 :
724 : !---------------------------Local variables-----------------------------
725 : real (dbl_kind),parameter :: dayspy = 365.0_dbl_kind ! days per year
726 : real (dbl_kind),parameter :: ve = 80.5_dbl_kind ! Calday of vernal equinox
727 : ! assumes Jan 1 = calday 1
728 :
729 329481 : real (dbl_kind) :: lambm ! Lambda m, mean long of perihelion (rad)
730 329481 : real (dbl_kind) :: lmm ! Intermediate argument involving lambm
731 329481 : real (dbl_kind) :: lamb ! Lambda, the earths long of perihelion
732 329481 : real (dbl_kind) :: invrho ! Inverse normalized sun/earth distance
733 329481 : real (dbl_kind) :: sinl ! Sine of lmm
734 :
735 : character(len=*),parameter :: subname='(shr_orb_decl)'
736 :
737 : ! Compute eccentricity factor and solar declination using
738 : ! day value where a round day (such as 213.0) refers to 0z at
739 : ! Greenwich longitude.
740 : !
741 : ! Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
742 : ! Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
743 : ! 35:2362-2367.
744 : !
745 : ! To get the earths true longitude (position in orbit; lambda in Berger
746 : ! 1978) which is necessary to find the eccentricity factor and declination,
747 : ! must first calculate the mean longitude (lambda m in Berger 1978) at
748 : ! the present day. This is done by adding to lambm0 (the mean longitude
749 : ! at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
750 : ! an increment (delta lambda m in Berger 1978) that is the number of
751 : ! days past or before (a negative increment) the vernal equinox divided by
752 : ! the days in a model year times the 2*pi radians in a complete orbit.
753 :
754 909228 : lambm = lambm0 + (calday - ve)*2._dbl_kind*pi/dayspy
755 909228 : lmm = lambm - mvelpp
756 :
757 : ! The earths true longitude, in radians, is then found from
758 : ! the formula in Berger 1978:
759 :
760 909228 : sinl = sin(lmm)
761 : lamb = lambm + eccen*(2._dbl_kind*sinl + eccen*(1.25_dbl_kind*sin(2._dbl_kind*lmm) &
762 909228 : & + eccen*((13.0_dbl_kind/12.0_dbl_kind)*sin(3._dbl_kind*lmm) - 0.25_dbl_kind*sinl)))
763 :
764 : ! Using the obliquity, eccentricity, moving vernal equinox longitude of
765 : ! perihelion (plus), and earths true longitude, the declination (delta)
766 : ! and the normalized earth/sun distance (rho in Berger 1978; actually inverse
767 : ! rho will be used), and thus the eccentricity factor (eccf), can be
768 : ! calculated from formulas given in Berger 1978.
769 :
770 909228 : invrho = (1._dbl_kind + eccen*cos(lamb - mvelpp)) / (1._dbl_kind - eccen*eccen)
771 :
772 : ! Set solar declination and eccentricity factor
773 :
774 909228 : delta = asin(sin(obliqr)*sin(lamb))
775 909228 : eccf = invrho*invrho
776 :
777 909228 : return
778 :
779 : END SUBROUTINE shr_orb_decl
780 : #endif
781 :
782 : !=======================================================================
783 :
784 : end module icepack_orbital
785 :
786 : !=======================================================================
|
