122 storageheatmethod, & !input
123 nsurf, sfr_surf, tsfc_surf, SnowFrac, alb, emis, IceFrac, & ! input:
124 NARP_TRANS_SITE, NARP_EMIS_SNOW, &
125 DTIME, ZENITH_deg, tsurf_0, kdown, Temp_C, RH, Press_hPa, qn1_obs, ldown_obs, &
127 AlbedoChoice, ldown_option, &
128 NetRadiationMethod_use, DiagQN, &
130 QSTARall, QSTAR_SF, QSTAR_S, kclear, KUPall, LDOWN, LUPall, fcld, TSURFall, & ! output:
131 qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_surf, &
132 albedo_snowfree, albedo_snow)
185 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: sfr_surf
186 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: tsfc_surf
187 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: SnowFrac
188 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: alb
189 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: emis
190 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) :: IceFrac
193 REAL(KIND(1D0)),
INTENT(in) :: DTIME
194 REAL(KIND(1D0)),
INTENT(in) :: ZENITH_deg
195 REAL(KIND(1D0)),
INTENT(in) :: tsurf_0
196 REAL(KIND(1D0)),
INTENT(in) :: kdown
197 REAL(KIND(1D0)),
INTENT(in) :: Temp_C
198 REAL(KIND(1D0)),
INTENT(in) :: RH
199 REAL(KIND(1D0)),
INTENT(in) :: Press_hPa
200 REAL(KIND(1D0)),
INTENT(in) :: qn1_obs
201 REAL(KIND(1D0)),
INTENT(in) :: ldown_obs
202 REAL(KIND(1D0)),
INTENT(in) :: SnowAlb
203 REAL(KIND(1D0)),
INTENT(in) :: NARP_TRANS_SITE
204 REAL(KIND(1D0)),
INTENT(in) :: NARP_EMIS_SNOW
206 INTEGER,
INTENT(in) :: nsurf
207 INTEGER,
INTENT(in) :: NetRadiationMethod_use
208 INTEGER,
INTENT(in) :: storageheatmethod
209 INTEGER,
INTENT(in) :: AlbedoChoice
210 INTEGER,
INTENT(in) :: ldown_option
211 INTEGER,
INTENT(in) :: DiagQN
213 REAL(KIND(1D0)),
INTENT(out) :: QSTARall
214 REAL(KIND(1D0)),
INTENT(out) :: QSTAR_SF
215 REAL(KIND(1D0)),
INTENT(out) :: QSTAR_S
216 REAL(KIND(1D0)),
INTENT(out) :: kclear
217 REAL(KIND(1D0)),
INTENT(out) :: KUPall
218 REAL(KIND(1D0)),
INTENT(out) :: LDOWN
219 REAL(KIND(1D0)),
INTENT(out) :: LUPall
220 REAL(KIND(1D0)),
INTENT(out) :: fcld
221 REAL(KIND(1D0)),
INTENT(out) :: TSURFall
223 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(out) :: qn1_ind_snow
224 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(out) :: kup_ind_snow
225 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(out) :: Tsurf_ind_snow
226 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(out) :: Tsurf_surf
227 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(out) :: qn_surf
229 REAL(KIND(1D0)),
INTENT(out) :: albedo_snowfree
230 REAL(KIND(1D0)),
INTENT(out) :: albedo_snow
232 REAL(KIND(1D0)),
DIMENSION(nsurf) :: tsfc_surf_K
233 REAL(KIND(1D0)),
DIMENSION(nsurf) :: kup_surf
234 REAL(KIND(1D0)),
DIMENSION(nsurf) :: lup_surf
236 REAL(KIND(1D0)),
DIMENSION(nsurf) :: qn1_ind_nosnow
237 REAL(KIND(1D0)),
DIMENSION(nsurf) :: kup_ind_nosnow
238 REAL(KIND(1D0)),
DIMENSION(nsurf) :: lup_ind_nosnow
239 REAL(KIND(1D0)),
DIMENSION(nsurf) :: Tsurf_ind_nosnow
240 REAL(KIND(1D0)),
DIMENSION(nsurf) :: lup_ind_snow
242 REAL(KIND(1D0)) :: tsurf_0_K
243 REAL(KIND(1D0)) :: Temp_K, TD, ZENITH, QSTAR, QSTAR_SNOW, KUP_SNOW, LUP_SNOW
244 REAL(KIND(1D0)) :: TSURF_SNOW_K, TSURF_SNOW_C, KUP, LUP, TSURF_K, TSURF_C
245 REAL(KIND(1D0)) :: EMIS0, EMIS_A, TRANS
246 REAL(KIND(1D0)) :: LUPCORR, SIGMATK4, KDOWN_HR
249 REAL(KIND(1D0)) :: qn1_cum, kup_cum, lup_cum, tsurf_cum, &
250 qn_is, kup_is, lup_is, tsurf_is, &
253 REAL(KIND(1D0)),
PARAMETER :: DEG2RAD = 0.017453292
255 REAL(KIND(1D0)),
PARAMETER :: SIGMA_SB = 5.67e-8
259 REAL(KIND(1D0)),
DIMENSION(365),
PARAMETER :: NARP_G = 3.0
265 tsfc_surf_k = tsfc_surf + 273.16
266 tsurf_0_k = tsurf_0 + 273.16
267 temp_k = temp_c + 273.16
268 sigmatk4 = sigma_sb*temp_k**4
273 zenith = zenith_deg*deg2rad
275 IF (doy == 366) doy = 365
291 IF (ldown_option .NE. 1)
THEN
303 albedo_snowfree = 0.2
304 albedo_snow = snowalb
306 IF (ldown_option == 1)
THEN
315 IF (sfr_surf(is) /= 0) sf_all = sf_all + sfr_surf(is)*(1 - snowfrac(is))
320 IF (diagqn == 1)
WRITE (*, *)
'is ', is
328 IF (albedochoice == 1 .AND. 180*zenith/acos(0.0) < 90)
THEN
329 albedo_snowfree = alb(is) + 0.5e-16*(180*zenith/acos(0.0))**8
331 albedo_snowfree = alb(is)
336 IF ((ldown_option == 4) .OR. (ldown_option == 5))
THEN
337 IF (zenith < 1.5)
THEN
339 kclear =
isurface(doy, zenith)*trans*narp_trans_site
340 IF (kclear > 50.)
THEN
344 IF (ldown_option == 5)
THEN
353 IF (ldown_option == 4)
THEN
356 ELSEIF ((ldown_option == 5))
THEN
361 ELSEIF (ldown_option == 3)
THEN
364 ELSEIF (ldown_option == 2)
THEN
367 IF (diagqn == 1)
WRITE (*, *)
'ldown_option: ', ldown_option,
'FCLD:', fcld
369 IF (ldown_option > 1)
THEN
370 ldown = emis_a*sigmatk4
371 IF (diagqn == 1)
WRITE (*, *)
'EMIS_A: ', emis_a,
'SIGMATK4:', sigmatk4,
'LDOWN: ', ldown
379 IF (kdown_hr > 0)
THEN
380 lupcorr = (1 - albedo_snowfree)*(0.08*kdown_hr)
385 IF (netradiationmethod_use < 10)
THEN
387 tsurf_k = ((emis0*sigmatk4 + lupcorr)/(emis0*sigma_sb))**0.25
388 lup = emis0*sigmatk4 + lupcorr + (1 - emis0)*ldown
391 IF (storageheatmethod == 5)
THEN
392 tsurf_k = tsfc_surf_k(is)
397 lup = emis0*sigma_sb*tsurf_k**4 + (1 - emis0)*ldown
401 kup = albedo_snowfree*kdown
404 qstar = kdown - kup + ldown - lup
405 tsurf_c = tsurf_k - 273.16
408 qn1_ind_nosnow(is) = qstar
409 kup_ind_nosnow(is) = kup
410 lup_ind_nosnow(is) = lup
411 tsurf_ind_nosnow(is) = tsurf_c
416 IF (snowfrac(is) > 0)
THEN
417 IF (albedochoice == 1 .AND. 180*zenith/acos(0.0) < 90)
THEN
418 albedo_snow = snowalb + 0.5e-16*(180*zenith/acos(0.0))**8
420 albedo_snow = snowalb
423 kup_snow = (albedo_snow*(snowfrac(is) - snowfrac(is)*icefrac(is)) + albedo_snowfree*snowfrac(is)*icefrac(is))*kdown
425 IF (netradiationmethod_use < 10)
THEN
427 tsurf_snow_k = ((narp_emis_snow*sigmatk4)/(narp_emis_snow*sigma_sb))**0.25
432 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow_k**4 + (1 - narp_emis_snow)*ldown
435 tsurf_snow_k = tsurf_0_k
436 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow_k**4 + (1 - narp_emis_snow)*ldown
439 qstar_snow = kdown - kup_snow + ldown - lup_snow
440 tsurf_snow_c = tsurf_snow_k - 273.16
453 qn1_ind_snow(is) = qstar_snow
454 kup_ind_snow(is) = kup_snow
455 lup_ind_snow(is) = lup_snow
456 tsurf_ind_snow(is) = tsurf_snow_k
458 IF (sf_all /= 0)
THEN
459 qstar_sf = qstar_sf + qstar*sfr_surf(is)*(1 - snowfrac(is))/sf_all
461 qstar_sf = qstar_sf + qstar*sfr_surf(is)*(1 - snowfrac(is))
464 IF ((1 - sf_all) /= 0)
THEN
465 qstar_s = qstar_s + qstar_snow*sfr_surf(is)*snowfrac(is)/(1 - sf_all)
467 qstar_s = qstar_s + qstar_snow*sfr_surf(is)*snowfrac(is)
473 qn_is = qstar*(1 - snowfrac(is)) + qstar_snow*snowfrac(is)
474 kup_is = kup*(1 - snowfrac(is)) + kup_snow*snowfrac(is)
475 lup_is = lup*(1 - snowfrac(is)) + lup_snow*snowfrac(is)
476 tsurf_is = tsurf_c*(1 - snowfrac(is)) + tsurf_snow_c*snowfrac(is)
478 IF (diagqn == 1)
WRITE (*, *)
'QSTAR', qstar,
'QSTAR_SNOW', qstar_snow,
'SnowFrac', snowfrac(is)
481 qn1_cum = qn1_cum + (qn_is*sfr_surf(is))
482 kup_cum = kup_cum + (kup_is*sfr_surf(is))
483 lup_cum = lup_cum + (lup_is*sfr_surf(is))
484 tsurf_cum = tsurf_cum + (tsurf_is*sfr_surf(is))
488 kup_surf(is) = kup_is
489 lup_surf(is) = lup_is
490 tsurf_surf(is) = tsurf_is
492 IF (diagqn == 1)
WRITE (*, *)
'qn1_is: ', qn_is
497 IF (netradiationmethod_use /= 3000)
THEN
516 IF (diagqn == 1)
WRITE (*, *)
'kdown: ', kdown,
'kup:', kup,
'LDOWN: ', ldown,
'LUP: ', lup
517 IF (diagqn == 1)
WRITE (*, *)
'Qn: ', qstarall
523 locationlatitude, locationlongitude, locationaltitude, &
524 sunazimuth, sunzenith)
527 REAL(KIND(1D0)),
INTENT(in) :: year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude
528 REAL(KIND(1D0)),
INTENT(out) :: sunazimuth, sunzenith
530 REAL(KIND(1D0)) :: sec
531 INTEGER :: month, day, hour, min, seas, dayofyear, year_int
533 REAL(KIND(1D0)) :: juliancentury, julianday, julianephemeris_century, julianephemeris_day, &
534 julianephemeris_millenium
535 REAL(KIND(1D0)) :: earth_heliocentric_positionlatitude, earth_heliocentric_positionlongitude, &
536 earth_heliocentric_positionradius
537 REAL(KIND(1D0)) :: sun_geocentric_positionlatitude, sun_geocentric_positionlongitude
538 REAL(KIND(1D0)) :: nutationlongitude, nutationobliquity
539 REAL(KIND(1D0)) :: corr_obliquity
540 REAL(KIND(1D0)) :: aberration_correction
541 REAL(KIND(1D0)) :: apparent_sun_longitude
542 REAL(KIND(1D0)) :: apparent_stime_at_greenwich
543 REAL(KIND(1D0)) :: sun_rigth_ascension
544 REAL(KIND(1D0)) :: sun_geocentric_declination
545 REAL(KIND(1D0)) :: observer_local_hour
546 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension, topocentric_sun_positionrigth_ascension_parallax
547 REAL(KIND(1D0)) :: topocentric_sun_positiondeclination
548 REAL(KIND(1D0)) :: topocentric_local_hour
566 dayofyear = floor(idectime)
568 CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)
572 CALL julian_calculation(year, month, day, hour, min, sec, utc, juliancentury, julianday, julianephemeris_century, &
573 julianephemeris_day, julianephemeris_millenium)
578 &earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
582 & sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
595 & aberration_correction, apparent_sun_longitude)
599 &corr_obliquity, apparent_stime_at_greenwich)
603 &sun_rigth_ascension)
608 &sun_geocentric_declination)
616 &topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination, locationaltitude,&
617 &locationlatitude, observer_local_hour, sun_rigth_ascension, sun_geocentric_declination,&
618 &earth_heliocentric_positionradius)
622 & topocentric_local_hour)
626 & topocentric_local_hour, sunazimuth, sunzenith)
702 &, earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
705 REAL(KIND(1D0)) :: julianephemeris_millenium
707 REAL(KIND(1D0)),
DIMENSION(64) :: A0
708 REAL(KIND(1D0)),
DIMENSION(34) :: A1
709 REAL(KIND(1D0)),
DIMENSION(20) :: A2
710 REAL(KIND(1D0)),
DIMENSION(7) :: A3
711 REAL(KIND(1D0)),
DIMENSION(3) :: A4
712 REAL(KIND(1D0)) :: A5
713 REAL(KIND(1D0)),
DIMENSION(64) :: B0
714 REAL(KIND(1D0)),
DIMENSION(34) :: B1
715 REAL(KIND(1D0)),
DIMENSION(20) :: B2
716 REAL(KIND(1D0)),
DIMENSION(7) :: B3
717 REAL(KIND(1D0)),
DIMENSION(3) :: B4
718 REAL(KIND(1D0)) :: B5
719 REAL(KIND(1D0)),
DIMENSION(64) :: C0
720 REAL(KIND(1D0)),
DIMENSION(34) :: C1
721 REAL(KIND(1D0)),
DIMENSION(20) :: C2
722 REAL(KIND(1D0)),
DIMENSION(7) :: C3
723 REAL(KIND(1D0)),
DIMENSION(3) :: C4
724 REAL(KIND(1D0)) :: C5
725 REAL(KIND(1D0)),
DIMENSION(40) :: A0j
726 REAL(KIND(1D0)),
DIMENSION(10) :: A1j
727 REAL(KIND(1D0)),
DIMENSION(6) :: A2j
728 REAL(KIND(1D0)),
DIMENSION(2) :: A3j
729 REAL(KIND(1D0)) :: A4j
730 REAL(KIND(1D0)),
DIMENSION(40) :: B0j
731 REAL(KIND(1D0)),
DIMENSION(10) :: B1j
732 REAL(KIND(1D0)),
DIMENSION(6) :: B2j
733 REAL(KIND(1D0)),
DIMENSION(2) :: B3j
734 REAL(KIND(1D0)) :: B4j
735 REAL(KIND(1D0)),
DIMENSION(40) :: C0j
736 REAL(KIND(1D0)),
DIMENSION(10) :: C1j
737 REAL(KIND(1D0)),
DIMENSION(6) :: C2j
738 REAL(KIND(1D0)),
DIMENSION(2) :: C3j
739 REAL(KIND(1D0)) :: C4j
740 REAL(KIND(1D0)),
DIMENSION(5) :: A0i
741 REAL(KIND(1D0)),
DIMENSION(5) :: B0i
742 REAL(KIND(1D0)),
DIMENSION(5) :: C0i
743 REAL(KIND(1D0)),
DIMENSION(2) :: A1i
744 REAL(KIND(1D0)),
DIMENSION(2) :: B1i
745 REAL(KIND(1D0)),
DIMENSION(2) :: C1i
746 REAL(KIND(1D0)) :: earth_heliocentric_positionlatitude
747 REAL(KIND(1D0)) :: earth_heliocentric_positionlongitude
748 REAL(KIND(1D0)) :: earth_heliocentric_positionradius
749 REAL(KIND(1D0)) :: JME
750 REAL(KIND(1D0)) :: L0
752 REAL(KIND(1D0)) :: L1
754 REAL(KIND(1D0)) :: L2
756 REAL(KIND(1D0)) :: L3
758 REAL(KIND(1D0)) :: L4
760 REAL(KIND(1D0)) :: L5
767 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
772 a0 = (/175347046, 3341656, 34894, 3497, 3418, 3136, 2676, 2343, 1324, 1273, 1199, 990, 902, 857, 780, 753, 505,&
773 &492, 357, 317, 284, 271, 243, 206, 205, 202, 156, 132, 126, 115, 103, 102, 102, 99, 98, 86, 85, 85, 80, 79, 71,&
774 &74, 74, 70, 62, 61, 57, 56, 56, 52, 52, 51, 49, 41, 41, 39, 37, 37, 36, 36, 33, 30, 30, 25/)
775 b0 = (/0., 4.669256800, 4.626100, 2.744100, 2.828900, 3.627700, 4.418100, 6.135200, 0.7425000, 2.037100,&
776 &1.109600, 5.233000, 2.045000, 3.508000, 1.179000, 2.533000, 4.583000, 4.205000, 2.92, 5.849000,&
777 &1.899000, 0.315, 0.345, 4.806000, 1.869000, 2.445800, 0.833, 3.411000, 1.083000, 0.645, 0.636, 0.976,&
778 &4.267000, 6.21, 0.68, 5.98, 1.3, 3.67, 1.81, 3.04, 1.76, 3.5, 4.68, 0.83, 3.98, 1.82, 2.78, 4.39, 3.47, 0.19,&
779 &1.33, 0.28, 0.49, 5.37, 2.4, 6.17, 6.04, 2.57, 1.71, 1.78, 0.59, 0.44, 2.74, 3.16/)
780 c0 = (/0., 6283.075850, 12566.15170, 5753.384900, 3.523100, 77713.77150, 7860.419400, 3930.209700,&
781 &11506.76980, 529.6910, 1577.343500, 5884.927, 26.29800, 398.1490, 5223.694, 5507.553,&
782 &18849.22800, 775.5230, 0.067, 11790.62900, 796.2980, 10977.07900, 5486.778, 2544.314,&
783 &5573.143, 6069.777, 213.2990, 2942.463, 20.77500, 0.98, 4694.003, 15720.83900, 7.114000,&
784 &2146.170, 155.4200, 161000.6900, 6275.960, 71430.70, 17260.15, 12036.46, 5088.630, 3154.690,&
785 &801.8200, 9437.760, 8827.390, 7084.900, 6286.600, 14143.50, 6279.550, 12139.55, 1748.020,&
786 &5856.480, 1194.450, 8429.240, 19651.05, 10447.39, 10213.29, 1059.380, 2352.870, 6812.770,&
787 &17789.85, 83996.85, 1349.870, 4690.480/)
788 a1 = (/628331966747.000, 206059., 4303., 425., 119., 109., &
789 93., 72., 68., 67., 59., 56., 45., 36., 29., 21., 19., 19., 17., 16., &
790 16., 15., 12., 12., 12., 12., 11., 10., 10., 9., 9., 8., 6., 6./)
791 b1 = (/0., 2.678235, 2.635100, 1.59, 5.796000, 2.966000, 2.59, 1.14, 1.87, 4.41, 2.89, 2.17, 0.40, 0.47,&
792 &2.65, 5.34, 1.85, 4.97, 2.99, 0.030, 1.43, 1.21, 2.83, 3.26, 5.27, 2.08, 0.77, 1.3, 4.24, 2.7, 5.64,&
794 c1 = (/0., 6283.075850, 12566.15170, 3.523000, 26.29800, 1577.344, 18849.23, 529.6900, 398.1500,&
795 &5507.550, 5223.690, 155.4200, 796.3000, 775.5200, 7.11, 0.98, 5486.780, 213.3000, 6275.960,&
796 &2544.310, 2146.170, 10977.08, 1748.020, 5088.630, 1194.450, 4694., 553.5700, 3286.600,&
797 &1349.870, 242.7300, 951.7200, 2352.870, 9437.760, 4690.480/)
798 a2 = (/52919, 8720, 309, 27, 16, 16, 10, 9, 7, 5, 4, 4, 3, 3, 3, 3, 3, 3, 2, 2/)
799 b2 = (/0., 1.072100, 0.867, 0.050, 5.19, 3.68, 0.76, 2.06, 0.83, 4.66, 1.03, 3.44, 5.14, 6.05, 1.19,&
800 &6.12, 0.31, 2.28, 4.38, 3.75/)
801 c2 = (/0., 6283.075800, 12566.15200, 3.52, 26.3, 155.4200, 18849.23, 77713.77, 775.5200, 1577.340,&
802 &7.11, 5573.140, 796.3000, 5507.550, 242.7300, 529.6900, 398.1500, 553.5700, 5223.690, 0.98/)
803 a3 = (/289, 35, 17, 3, 1, 1, 1/)
804 b3 = (/5.8440, 0., 5.4900, 5.2000, 4.7200, 5.3000, 5.9700/)
805 c3 = (/6283.076, 0., 12566.15, 155.4200, 3.52, 18849.23, 242.7300/)
807 b4 = (/3.1420, 4.1300, 3.8400/)
808 c4 = (/0., 6283.08, 12566.15/)
813 jme = julianephemeris_millenium
816 l0 = sum(a0*cos(b0 + (c0*jme)))
817 l1 = sum(a1*cos(b1 + (c1*jme)))
818 l2 = sum(a2*cos(b2 + (c2*jme)))
819 l3 = sum(a3*cos(b3 + (c3*jme)))
820 l4 = sum(a4*cos(b4 + (c4*jme)))
821 l5 = a5*cos(b5 + (c5*jme))
823 earth_heliocentric_positionlongitude = &
824 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4) + (l5*jme**5))/1e8
826 earth_heliocentric_positionlongitude = earth_heliocentric_positionlongitude*180./pi
828 earth_heliocentric_positionlongitude =
set_to_range(earth_heliocentric_positionlongitude)
830 a0i = (/280, 102, 80, 44, 32/)
831 b0i = (/3.19900000000000, 5.42200000000000, 3.88000000000000, 3.70000000000000, 4./)
832 c0i = (/84334.6620000000, 5507.55300000000, 5223.69000000000, 2352.87000000000, 1577.34000000000/)
834 b1i = (/3.90000000000000, 1.73000000000000/)
835 c1i = (/5507.55000000000, 5223.69000000000/)
837 l0 = sum(a0i*cos(b0i + (c0i*jme)))
838 l1 = sum(a1i*cos(b1i + (c1i*jme)))
840 earth_heliocentric_positionlatitude = (l0 + (l1*jme))/1e8
842 earth_heliocentric_positionlatitude = earth_heliocentric_positionlatitude*180/pi
844 earth_heliocentric_positionlatitude =
set_to_range(earth_heliocentric_positionlatitude)
846 a0j = (/100013989, 1670700, 13956, 3084, 1628, 1576, 925, 542, 472, 346, 329, 307, 243, 212, 186, 175, 110,&
847 &98, 86, 86, 85, 63, 57, 56, 49, 47, 45, 43, 39, 38, 37, 37, 36, 35, 33, 32, 32, 28, 28, 26/)
848 b0j = (/0., 3.09846350, 3.05525000, 5.1985, 1.1739, 2.8469, 5.453, 4.564, 3.661, 0.964, 5.90, 0.299,&
849 &4.273, 5.847, 5.022, 3.012, 5.055, 0.890, 5.69, 1.27, 0.270, 0.920, 2.01, 5.24, 3.25, 2.58, 5.54,&
850 &6.01, 5.36, 2.39, 0.830, 4.90, 1.67, 1.84, 0.240, 0.180, 1.78, 1.21, 1.90, 4.59/)
851 c0j = (/0., 6283.07585, 12566.1517, 77713.7715, 5753.38490, 7860.41940, 11506.7700, 3930.21000,&
852 &5884.92700, 5507.55300, 5223.69400, 5573.14300, 11790.6290, 1577.34400, 10977.0790, 18849.2280,&
853 &5486.77800, 6069.78000, 15720.8400, 161000.690, 17260.1500, 529.69, 83996.8500, 71430.7000,&
854 &2544.31000, 775.52, 9437.76000, 6275.96000, 4694., 8827.39000, 19651.0500, 12139.5500,&
855 &12036.4600, 2942.46000, 7084.9, 5088.63000, 398.15, 6286.6, 6279.55000, 10447.3900/)
856 a1j = (/103019, 1721, 702, 32, 31, 25, 18, 10, 9, 9/)
857 b1j = (/1.10749000, 1.0644, 3.142, 1.02, 2.84, 1.32, 1.42, 5.91, 1.42, 0.270/)
858 c1j = (/6283.07585, 12566.1517, 0., 18849.2300, 5507.55000, 5223.69000, 1577.34000, 10977.0800,&
859 &6275.96000, 5486.78000/)
860 a2j = (/4359, 124, 12, 9, 6, 3/)
861 b2j = (/5.7846, 5.579, 3.14, 3.63, 1.87, 5.47/)
862 c2j = (/6283.07580, 12566.1520, 0., 77713.7700, 5573.14000, 18849./)
864 b3j = (/4.273, 3.92/)
865 c3j = (/6283.07600, 12566.1500/)
871 l0 = sum(a0j*cos(b0j + (c0j*jme)))
872 l1 = sum(a1j*cos(b1j + (c1j*jme)))
873 l2 = sum(a2j*cos(b2j + (c2j*jme)))
874 l3 = sum(a3j*cos(b3j + (c3j*jme)))
875 l4 = a4j*cos(b4j + (c4j*jme))
878 earth_heliocentric_positionradius = &
879 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4))/1e8
907 REAL(KIND(1D0)),
INTENT(in) :: julianephemeris_century
908 REAL(KIND(1D0)),
DIMENSION(63) :: delta_longitude
909 REAL(KIND(1D0)),
DIMENSION(63) :: delta_obliquity
910 REAL(KIND(1D0)) :: JCE
911 REAL(KIND(1D0)) :: nutationlongitude
912 REAL(KIND(1D0)) :: nutationobliquity
913 REAL(KIND(1D0)),
DIMENSION(4) :: p0, p1, p2, p3, p4
914 REAL(KIND(1D0)),
DIMENSION(63) :: tabulated_argument
915 REAL(KIND(1D0)) :: X0
916 REAL(KIND(1D0)) :: X1
917 REAL(KIND(1D0)) :: X2
918 REAL(KIND(1D0)) :: X3
919 REAL(KIND(1D0)) :: X4
920 REAL(KIND(1D0)),
DIMENSION(5) :: Xi
921 INTEGER,
DIMENSION(315) :: Y_terms1
922 INTEGER,
DIMENSION(5, 63) :: Y_terms
923 REAL(KIND(1D0)),
DIMENSION(252) :: nutation_terms1
924 REAL(KIND(1D0)),
DIMENSION(4, 63) :: nutation_terms
926 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
932 jce = julianephemeris_century
935 p0 = (/(1/189474.), -0.0019142, 445267.11148, 297.85036/)
937 x0 = p0(1)*jce**3 + p0(2)*jce**2 + p0(3)*jce + p0(4)
940 p1 = (/-(1/300000.), -0.0001603, 35999.05034, 357.52772/)
942 x1 = p1(1)*jce**3 + p1(2)*jce**2 + p1(3)*jce + p1(4)
945 p2 = (/(1/56250.), 0.0086972, 477198.867398, 134.96298/)
947 x2 = p2(1)*jce**3 + p2(2)*jce**2 + p2(3)*jce + p2(4)
950 p3 = (/(1/327270.), -0.0036825, 483202.017538, 93.27191/)
952 x3 = p3(1)*jce**3 + p3(2)*jce**2 + p3(3)*jce + p3(4)
956 p4 = (/(1/450000.), 0.0020708, -1934.136261, 125.04452/)
958 x4 = p4(1)*jce**3 + p4(2)*jce**2 + p4(3)*jce + p4(4)
961 y_terms1 = (/0, 0, 0, 0, 1, -2, 0, 0, 2, 2, 0, 0, 0, 2, 2, 0, 0, &
962 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 2, 2, 0, 0, 0, 2, 1, &
963 0, 0, 1, 2, 2, -2, -1, 0, 2, 2, -2, 0, 1, 0, 0, -2, 0, 0, 2, 1, 0, 0, -1, 2, 2, 2, 0, &
965 1, 0, 1, 2, 0, -1, 2, 2, &
966 0, 0, -1, 0, 1, 0, 0, 1, 2, 1, -2, 0, 2, 0, 0, 0, 0, -2, 2, 1, 2, 0, 0, 2, 2, 0, 0, 2, &
968 0, 0, -2, 0, 1, 2, 2, 0, &
969 0, 0, 2, 0, -2, 0, 0, 2, 0, 0, 0, -1, 2, 1, 0, 2, 0, 0, 0, &
970 2, 0, -1, 0, 1, -2, 2, 0, 2, 2, 0, 1, 0, 0, 1, &
971 -2, 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 0, &
972 2, -2, 0, 2, 0, -1, 2, 1, 2, 0, 1, 2, 2, 0, 1, 0, 2, 2, -2, &
973 1, 1, 0, 0, 0, -1, 0, 2, 2, 2, 0, 0, 2, 1, 2, &
974 0, 1, 0, 0, -2, 0, 2, 2, 2, -2, 0, 1, &
975 2, 1, 2, 0, -2, 0, 1, 2, 0, 0, 0, 1, 0, -1, 1, 0, 0, -2, -1, &
976 0, 2, 1, -2, 0, 0, 0, 1, 0, 0, 2, 2, 1, -2, &
977 0, 2, 0, 1, -2, 1, 0, 2, 1, 0, 0, 1, -2, &
978 0, -1, 0, 1, 0, 0, -2, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 2, &
979 0, 0, 0, -2, 2, 2, -1, -1, 1, 0, 0, 0, 1, &
980 1, 0, 0, 0, -1, 1, 2, 2, 2, -1, -1, 2, 2, &
981 0, 0, 3, 2, 2, 2, -1, 0, 2, 2/)
982 y_terms = reshape(y_terms1, (/5, 63/))
983 nutation_terms1 = (/-171996., -174.2, 92025., 8.9, -13187., -1.6, 5736., -3.1, -2274., -0.2, 977., -0.5, 2062., &
985 1426., -3.4, 54., -0.1, 712., 0.1, -7., 0., -517., 1.2, 224., -0.6, -386., -0.4, 200., 0., -301., &
987 217., -0.5, -95., 0.3, -158., 0., 0., 0., 129., 0.1, -70., 0., 123., 0., -53., 0., 63., 0., 0., 0., &
988 63., 0.1, -33., 0., -59., 0., 26., 0., &
989 -58., -0.1, 32., 0., -51., 0., 27., 0., 48., 0., 0., 0., 46., 0., -24., 0., -38., &
990 0., 16., 0., -31., 0., 13., 0., 29., 0., &
991 0., 0., 29., 0., -12., 0., 26., 0., 0., 0., -22., 0., 0., 0., 21., 0., -10., 0., 17., -0.1, 0., 0., 16., &
992 0., -8., 0., -16., 0.1, 7., 0., &
993 -15., 0., 9., 0., -13., 0., 7., 0., -12., 0., 6., 0., 11., 0., 0., 0., &
994 -10., 0., 5., 0., -8., 0., 3., 0., 7., &
995 0., -3., 0., -7., 0., 0., 0., &
996 -7., 0., 3., 0., -7., 0., 3., 0., 6., 0., 0., 0., 6., 0., -3., 0., 6., 0., &
997 -3., 0., -6., 0., 3., 0., -6., 0., 3., &
998 0., 5., 0., 0., 0., -5., 0., &
999 3., 0., -5., 0., 3., 0., -5., 0., 3., 0., 4., 0., 0., 0., 4., 0., 0., 0., 4., 0., 0., 0., -4., 0., 0., &
1000 0., -4., 0., 0., 0., -4., 0., 0., 0., &
1001 3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., &
1002 0., -3., 0., 0., 0., -3., 0., 0., 0./)
1003 nutation_terms = reshape(nutation_terms1, (/4, 63/))
1006 xi = (/x0, x1, x2, x3, x4/)
1009 tabulated_argument(i) = ((y_terms(1, i)*xi(1)) &
1010 + (y_terms(2, i)*xi(2)) + (y_terms(3, i)*xi(3)) &
1011 + (y_terms(4, i)*xi(4)) + (y_terms(5, i)*xi(5)))*pi/180
1014 delta_longitude = ((nutation_terms(1, :) + (nutation_terms(2, :)*jce)))*sin(tabulated_argument)
1015 delta_obliquity = ((nutation_terms(3, :) + (nutation_terms(4, :)*jce)))*cos(tabulated_argument)
1018 nutationlongitude = sum(delta_longitude)/36000000.0
1021 nutationobliquity = sum(delta_obliquity)/36000000.0