45 NetRadiationMethod, & !input
47 NetRadiationMethod_use, AlbedoChoice, ldown_option)
49 INTEGER,
INTENT(in) :: NetRadiationMethod
50 INTEGER,
INTENT(in) :: SnowUse
51 INTEGER,
INTENT(out) :: NetRadiationMethod_use
52 INTEGER,
INTENT(out) :: AlbedoChoice, ldown_option
56 IF (netradiationmethod == 0)
THEN
57 netradiationmethod_use = 0
60 IF (snowuse == 1)
THEN
62 netradiationmethod_use = 3000
67 ELSEIF (netradiationmethod > 0)
THEN
69 IF (netradiationmethod < 100)
THEN
72 netradiationmethod_use = mod(netradiationmethod, 10)
73 IF (netradiationmethod_use == 1) ldown_option = 1
74 IF (netradiationmethod_use == 2) ldown_option = 2
75 IF (netradiationmethod_use == 3) ldown_option = 3
77 netradiationmethod_use = netradiationmethod
85 ELSEIF (netradiationmethod >= 100 .AND. netradiationmethod < 1000)
THEN
87 IF (netradiationmethod == 100) ldown_option = 1
88 IF (netradiationmethod == 200) ldown_option = 2
89 IF (netradiationmethod == 300) ldown_option = 3
90 netradiationmethod_use = netradiationmethod/100
93 ELSEIF (netradiationmethod > 1000)
THEN
95 netradiationmethod_use = mod(netradiationmethod, 10)
96 IF (netradiationmethod_use == 1) ldown_option = 1
97 IF (netradiationmethod_use == 2) ldown_option = 2
98 IF (netradiationmethod_use == 3) ldown_option = 3
100 netradiationmethod_use = netradiationmethod
102 IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9)
THEN
103 WRITE (*, *)
'NetRadiationMethod=', netradiationmethod_use
104 WRITE (*, *)
'Value not usable'
111 IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9)
THEN
112 WRITE (*, *)
'NetRadiationMethod=', netradiationmethod_use
113 WRITE (*, *)
'Value not usable'
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)
631 SUBROUTINE julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century&
632 &, julianephemeris_day, julianephemeris_millenium)
635 REAL(KIND(1D0)) :: A, B, D, delta_t
636 REAL(KIND(1D0)) :: juliancentury
637 REAL(KIND(1D0)) :: julianday
638 REAL(KIND(1D0)) :: julianephemeris_century
639 REAL(KIND(1D0)) :: julianephemeris_day
640 REAL(KIND(1D0)) :: julianephemeris_millenium
641 REAL(KIND(1D0)) :: M, sec, year, UTC
642 INTEGER :: day, hour, min, month
644 REAL(KIND(1D0)) :: ut_time, Y
650 IF (month == 1 .OR. month == 2)
THEN
657 ut_time = ((float(hour) - utc)/24.) + (float(min)/(60.*24.)) + (sec/(60.*60.*24.))
661 IF (year == 1582.)
THEN
662 IF (month == 10)
THEN
665 ELSE IF (day >= 15)
THEN
667 b = 2 - a + floor(a/4)
674 ELSE IF (month < 10)
THEN
678 b = 2 - a + floor(a/4)
681 ELSE IF (year < 1582.)
THEN
685 b = 2 - a + floor(a/4)
688 julianday = floor(365.25*(y + 4716.)) + floor(30.6001*(m + 1)) + d + b - 1524.5
691 julianephemeris_day = julianday + (delta_t/86400)
693 juliancentury = (julianday - 2451545.)/36525.
695 julianephemeris_century = (julianephemeris_day - 2451545.)/36525.
697 julianephemeris_millenium = julianephemeris_century/10.
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
884 &earth_heliocentric_positionlatitude, sun_geocentric_positionlatitude, &
885 &sun_geocentric_positionlongitude)
888 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionlongitude
889 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionlatitude
890 REAL(KIND(1D0)) :: sun_geocentric_positionlatitude
891 REAL(KIND(1D0)) :: sun_geocentric_positionlongitude
895 sun_geocentric_positionlongitude = earth_heliocentric_positionlongitude + 180.0
897 sun_geocentric_positionlongitude =
set_to_range(sun_geocentric_positionlongitude)
899 sun_geocentric_positionlatitude = -earth_heliocentric_positionlatitude
901 sun_geocentric_positionlatitude =
set_to_range(sun_geocentric_positionlatitude)
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
1028 REAL(KIND(1D0)),
INTENT(out) :: corr_obliquity
1029 REAL(KIND(1D0)),
INTENT(in) :: julianephemeris_millenium
1030 REAL(KIND(1D0)),
INTENT(in) :: nutationobliquity
1031 REAL(KIND(1D0)) :: mean_obliquity
1032 REAL(KIND(1D0)),
DIMENSION(11) :: p
1033 REAL(KIND(1D0)) :: U
1037 p = (/2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 84381.448/)
1040 u = julianephemeris_millenium/10
1042 &p(1)*u**10 + p(2)*u**9 + p(3)*u**8 + p(4)*u**7 + p(5)*u**6 + p(6)*u**5 + p(7)*u**4 + &
1043 &p(8)*u**3 + p(9)*u**2 + p(10)*u + p(11)
1045 corr_obliquity = (mean_obliquity/3600) + nutationobliquity
1051 REAL(KIND(1D0)),
INTENT(out) :: aberration_correction
1052 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionradius
1057 aberration_correction = -20.4898/(3600*earth_heliocentric_positionradius)
1062 & aberration_correction, apparent_sun_longitude)
1065 REAL(KIND(1D0)),
INTENT(in) :: aberration_correction
1066 REAL(KIND(1D0)),
INTENT(out) :: apparent_sun_longitude
1067 REAL(KIND(1D0)),
INTENT(in) :: nutationlongitude
1068 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlongitude
1072 apparent_sun_longitude = sun_geocentric_positionlongitude + nutationlongitude + aberration_correction
1077 & corr_obliquity, apparent_stime_at_greenwich)
1080 REAL(KIND(1D0)),
INTENT(out) :: apparent_stime_at_greenwich
1081 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1082 REAL(KIND(1D0)),
INTENT(in) :: julianday
1083 REAL(KIND(1D0)),
INTENT(in) :: juliancentury
1084 REAL(KIND(1D0)),
INTENT(in) :: nutationlongitude
1085 REAL(KIND(1D0)) :: JC
1086 REAL(KIND(1D0)) :: JD
1087 REAL(KIND(1D0)) :: mean_stime
1088 REAL(KIND(1D0)),
PARAMETER :: pi = 3.14159265358979d+0
1096 mean_stime = 280.46061837d+0 + (360.98564736629d+0*(jd - 2451545.0d+0)) + (0.000387933d+0*jc**2) - (jc**3/38710000.0d+0)
1101 apparent_stime_at_greenwich = mean_stime + (nutationlongitude*cos(corr_obliquity*pi/180))
1105 &sun_geocentric_positionlatitude, sun_rigth_ascension)
1108 REAL(KIND(1D0)),
INTENT(in) :: apparent_sun_longitude
1109 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1110 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlatitude
1111 REAL(KIND(1D0)),
INTENT(out) :: sun_rigth_ascension
1112 REAL(KIND(1D0)) :: argument_denominator
1113 REAL(KIND(1D0)) :: argument_numerator
1114 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1118 argument_numerator = (sin(apparent_sun_longitude*pi/180.0)*cos(corr_obliquity*pi/180.0)) - &
1119 (tan(sun_geocentric_positionlatitude*pi/180.0)*sin(corr_obliquity*pi/180.0))
1120 argument_denominator = cos(apparent_sun_longitude*pi/180.0)
1122 sun_rigth_ascension = atan2(argument_numerator, argument_denominator)*180.0/pi
1124 sun_rigth_ascension =
set_to_range(sun_rigth_ascension)
1128 &sun_geocentric_positionlatitude, sun_geocentric_declination)
1131 REAL(KIND(1D0)),
INTENT(in) :: apparent_sun_longitude
1132 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1133 REAL(KIND(1D0)),
INTENT(out) :: sun_geocentric_declination
1134 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlatitude
1135 REAL(KIND(1D0)) :: argument
1136 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1138 argument = (sin(sun_geocentric_positionlatitude*pi/180.0)*cos(corr_obliquity*pi/180.0)) + &
1139 (cos(sun_geocentric_positionlatitude*pi/180.0)*sin(corr_obliquity*pi/180)*sin(apparent_sun_longitude*pi/180.0))
1141 sun_geocentric_declination = asin(argument)*180.0/pi
1145 &sun_rigth_ascension, observer_local_hour)
1148 REAL(KIND(1D0)),
INTENT(in) :: apparent_stime_at_greenwich
1149 REAL(KIND(1D0)),
INTENT(in) :: locationlongitude
1150 REAL(KIND(1D0)),
INTENT(out) :: observer_local_hour
1151 REAL(KIND(1D0)),
INTENT(in) :: sun_rigth_ascension
1153 observer_local_hour = apparent_stime_at_greenwich + locationlongitude - sun_rigth_ascension
1155 observer_local_hour =
set_to_range(observer_local_hour)
1159 &, topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination,&
1160 &locationaltitude, locationlatitude, observer_local_hour, sun_rigth_ascension,&
1161 &sun_geocentric_declination, earth_heliocentric_positionradius)
1164 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionradius
1165 REAL(KIND(1D0)),
INTENT(in) :: locationlatitude
1166 REAL(KIND(1D0)),
INTENT(in) :: locationaltitude
1167 REAL(KIND(1D0)),
INTENT(in) :: observer_local_hour
1168 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_declination
1169 REAL(KIND(1D0)),
INTENT(in) :: sun_rigth_ascension
1170 REAL(KIND(1D0)) :: denominator
1171 REAL(KIND(1D0)) :: eq_horizontal_parallax
1172 REAL(KIND(1D0)) :: nominator
1173 REAL(KIND(1D0)) :: sun_rigth_ascension_parallax
1174 REAL(KIND(1D0)) :: topocentric_sun_positiondeclination
1175 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension
1176 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension_parallax
1177 REAL(KIND(1D0)) :: u
1178 REAL(KIND(1D0)) :: x
1179 REAL(KIND(1D0)) :: y
1180 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1187 eq_horizontal_parallax = 8.794/(3600*earth_heliocentric_positionradius)
1190 u = atan(0.99664719*tan(locationlatitude*pi/180))
1193 x = cos(u) + ((locationaltitude/6378140)*cos(locationaltitude*pi/180))
1196 y = (0.99664719d+0*sin(u)) + ((locationaltitude/6378140)*sin(locationlatitude*pi/180))
1199 nominator = -x*sin(eq_horizontal_parallax*pi/180.0)*sin(observer_local_hour*pi/180.0)
1200 denominator = cos(sun_geocentric_declination*pi/180.0) - &
1201 (x*sin(eq_horizontal_parallax*pi/180.0)*cos(observer_local_hour*pi/180.0))
1202 sun_rigth_ascension_parallax = atan2(nominator, denominator)
1204 topocentric_sun_positionrigth_ascension_parallax = sun_rigth_ascension_parallax*180.0/pi
1207 topocentric_sun_positionrigth_ascension = sun_rigth_ascension + (sun_rigth_ascension_parallax*180.0/pi)
1210 nominator = (sin(sun_geocentric_declination*pi/180.0) - (y*sin(eq_horizontal_parallax*pi/180.0)))&
1211 & *cos(sun_rigth_ascension_parallax)
1212 denominator = cos(sun_geocentric_declination*pi/180.0) - (y*sin(eq_horizontal_parallax*pi/180.0))&
1213 & *cos(observer_local_hour*pi/180.0)
1214 topocentric_sun_positiondeclination = atan2(nominator, denominator)*180.0/pi
1218 & topocentric_local_hour)
1221 REAL(KIND(1D0)),
INTENT(in) :: observer_local_hour
1222 REAL(KIND(1D0)),
INTENT(out) :: topocentric_local_hour
1223 REAL(KIND(1D0)),
INTENT(in) :: topocentric_sun_positionrigth_ascension_parallax
1227 topocentric_local_hour = observer_local_hour - topocentric_sun_positionrigth_ascension_parallax
1231 &topocentric_local_hour, sunazimuth, sunzenith)
1234 REAL(KIND(1D0)),
INTENT(in) :: locationlatitude
1235 REAL(KIND(1D0)),
INTENT(in) :: topocentric_local_hour
1236 REAL(KIND(1D0)),
INTENT(in) :: topocentric_sun_positiondeclination
1237 REAL(KIND(1D0)) :: corr_elevation
1238 REAL(KIND(1D0)) :: apparent_elevation
1239 REAL(KIND(1D0)) :: argument
1240 REAL(KIND(1D0)) :: denominator
1241 REAL(KIND(1D0)) :: nominator
1242 REAL(KIND(1D0)) :: refraction_corr
1243 REAL(KIND(1D0)) :: sunazimuth
1244 REAL(KIND(1D0)) :: sunzenith
1245 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1251 argument = (sin(locationlatitude*pi/180.0)*sin(topocentric_sun_positiondeclination*pi/180.0)) + &
1252 (cos(locationlatitude*pi/180.0)*cos(topocentric_sun_positiondeclination*pi/180.0)* &
1253 &cos(topocentric_local_hour*pi/180.0))
1254 corr_elevation = asin(argument)*180.0/pi
1257 argument = corr_elevation + (10.3/(corr_elevation + 5.11))
1258 refraction_corr = 1.02/(60*tan(argument*pi/180.0))
1265 apparent_elevation = corr_elevation + refraction_corr
1267 sunzenith = 90.0 - apparent_elevation
1272 nominator = sin(topocentric_local_hour*pi/180.0)
1273 denominator = (cos(topocentric_local_hour*pi/180.0)*sin(locationlatitude*pi/180.0)) - &
1274 (tan(topocentric_sun_positiondeclination*pi/180.0)*cos(locationlatitude*pi/180.0))
1275 sunazimuth = (atan2(nominator, denominator)*180.0/pi) + 180.0
1284 REAL(kind(1d0)) :: max_interval
1285 REAL(kind(1d0)) :: min_interval
1286 REAL(kind(1d0)) :: var
1287 REAL(kind(1d0)) :: vari
1289 max_interval = 360.0
1292 vari = var - max_interval*floor(var/max_interval)
1294 IF (vari < min_interval)
THEN
1295 vari = vari + max_interval
1308 REAL(kind(1d0)) :: rh, td, temp_c, g
1310 g = ((17.27*temp_c)/(237.7 + temp_c)) + log(rh/100)
1311 td = (237.7*g)/(17.27 - g)
1317 REAL(kind(1d0)) :: temp_k, ea_hpa, emis_a
1318 REAL(kind(1d0)) :: w
1320 w = 46.5*(ea_hpa/temp_k)
1321 emis_a = 1.-(1.+w)*exp(-sqrt(1.2 + 3.*w))
1326 REAL(kind(1d0)) :: emis_a, fcld, em_adj
1329 em_adj = emis_a + (1.-emis_a)*fcld
1334 REAL(kind(1d0)) :: emis_a, fcld, em_adj
1335 em_adj = emis_a + (1.-emis_a)*fcld*fcld
1339 REAL(kind(1d0)) :: kdown, kclear, fcld
1341 fcld = 1.-kdown/kclear
1342 IF (fcld > 1.) fcld = 1.
1343 IF (fcld < 0.) fcld = 0.
1349 REAL(kind(1d0)),
INTENT(in) :: rh
1350 REAL(kind(1d0)),
INTENT(in) :: temp
1352 REAL(kind(1d0)) :: fwc
1353 REAL(kind(1d0)) :: a, b
1360 b = 0.00019*temp + 0.015
1363 fwc = a*(exp(b*rh) - 1)
1364 IF (fwc > 1.) fwc = 1.
1365 IF (fwc < 0.) fwc = 0.
1393 REAL(kind(1d0)) :: zenith, isurf
1395 REAL(kind(1d0)) :: rmean, rse, cosz, itoa
1396 REAL(kind(1d0)),
PARAMETER :: deg2rad = 0.017453292
1400 IF (zenith < 90.*deg2rad)
THEN
1402 itoa = 1370.*(rmean/rse)**2
1414 REAL(kind(1d0)) :: rse
1415 REAL(kind(1d0)) :: ma, nu, e, a
1420 ma = 2.*3.141592654*(doy - 3)/365.25463
1421 nu = ma + 0.0333988*sin(ma) + .0003486*sin(2.*ma) + 5e-6*sin(3.*ma)
1422 rse = a*(1 - e*e)/(1 + e*cos(nu))
1434 INTEGER :: lat, ios, ilat
1435 REAL(kind(1d0)),
DIMENSION(365) :: g
1443 READ (99, *, iostat=ios) ilat, g
1459 REAL(kind(1d0)) :: press_hpa, temp_c_dew, zenith, g, trans
1460 REAL(kind(1d0)) :: m, trtpg, u, tw, ta, cosz
1461 REAL(kind(1d0)) :: tdf
1462 REAL(kind(1d0)),
PARAMETER :: deg2rad = 0.017453292
1464 IF (zenith > 80.*deg2rad)
THEN
1465 cosz = cos(80.*deg2rad)
1470 tdf = temp_c_dew*1.8 + 32.
1472 m = 35*cosz/sqrt(1224.*cosz*cosz + 1)
1474 trtpg = 1.021 - 0.084*sqrt(m*(0.000949*press_hpa + 0.051))
1475 u = exp(0.113 - log(g + 1) + 0.0393*tdf)
1476 tw = 1 - 0.077*(u*m)**0.3
character(len=90) smithfile
subroutine abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)
real(kind(1d0)) function wc_fraction(RH, Temp)
subroutine radmethod(NetRadiationMethod, SnowUse, NetRadiationMethod_use, AlbedoChoice, ldown_option)
subroutine topocentric_sun_position_calculate(topocentric_sun_positionrigth_ascension, topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination, locationaltitude, locationlatitude, observer_local_hour, sun_rigth_ascension, sun_geocentric_declination, earth_heliocentric_positionradius)
subroutine sun_geocentric_declination_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, sun_geocentric_declination)
real(kind(1d0)) function, dimension(365) smithlambda(lat)
real(kind(1d0)) function dewpoint(Temp_C, rh)
subroutine sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination, topocentric_local_hour, sunazimuth, sunzenith)
real(kind(1d0)) function emis_cloud_sq(EMIS_A, FCLD)
subroutine apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude, aberration_correction, apparent_sun_longitude)
subroutine corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)
subroutine earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude, earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
real(kind(1d0)) function isurface(doy, zenith)
subroutine apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, corr_obliquity, apparent_stime_at_greenwich)
real(kind(1d0)) function cloud_fraction(KDOWN, KCLEAR)
subroutine narp(storageheatmethod, nsurf, sfr_surf, tsfc_surf, SnowFrac, alb, emis, IceFrac, NARP_TRANS_SITE, NARP_EMIS_SNOW, DTIME, ZENITH_deg, tsurf_0, kdown, Temp_C, RH, Press_hPa, qn1_obs, ldown_obs, SnowAlb, AlbedoChoice, ldown_option, NetRadiationMethod_use, DiagQN, qn_surf, QSTARall, QSTAR_SF, QSTAR_S, kclear, KUPall, LDOWN, LUPall, fcld, TSURFall, qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_surf, albedo_snowfree, albedo_snow)
subroutine sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude, sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
real(kind(1d0)) function emis_cloud(EMIS_A, FCLD)
subroutine topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax, topocentric_local_hour)
subroutine sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, sun_rigth_ascension)
subroutine narp_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
subroutine observer_local_hour_calculation(apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)
subroutine julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century, julianephemeris_day, julianephemeris_millenium)
real(kind(1d0)) function prata_emis(Temp_K, EA_hPa)
real(kind(1d0)) function solar_esdist(doy)
subroutine nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)
real(kind(1d0)) function transmissivity(Press_hPa, Temp_C_dew, G, zenith)
real(kind(1d0)) function set_to_range(var)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine day2month(b, mb, md, seas, year, latitude)
subroutine dectime_to_timevec(dectime, HOURS, MINS, SECS)