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