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
57 IF (netradiationmethod == 0)
THEN 58 netradiationmethod_use = 0
61 IF (snowuse == 1)
THEN 63 netradiationmethod_use = 3000
68 ELSEIF (netradiationmethod > 0)
THEN 70 IF (netradiationmethod < 100)
THEN 73 netradiationmethod_use = mod(netradiationmethod, 10)
74 IF (netradiationmethod_use == 1) ldown_option = 1
75 IF (netradiationmethod_use == 2) ldown_option = 2
76 IF (netradiationmethod_use == 3) ldown_option = 3
78 netradiationmethod_use = netradiationmethod
86 ELSEIF (netradiationmethod >= 100 .AND. netradiationmethod < 1000)
THEN 88 IF (netradiationmethod == 100) ldown_option = 1
89 IF (netradiationmethod == 200) ldown_option = 2
90 IF (netradiationmethod == 300) ldown_option = 3
92 netradiationmethod_use = netradiationmethod/100
96 IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9)
THEN 97 WRITE (*, *)
'NetRadiationMethod=', netradiationmethod_use
98 WRITE (*, *)
'Value not usable' 107 nsurf, sfr, SnowFrac, alb, emis, IceFrac, &! input:
108 NARP_TRANS_SITE, NARP_EMIS_SNOW, &
109 DTIME, ZENITH_deg, tsurf_0, kdown, Temp_C, RH, Press_hPa, qn1_obs, &
111 AlbedoChoice, ldown_option, &
112 NetRadiationMethod_use, DiagQN, &
113 QSTARall, QSTAR_SF, QSTAR_S, kclear, KUPall, LDOWN, LUPall, fcld, TSURFall, &! output:
114 qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_ind, &
167 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) ::sfr
168 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) ::SnowFrac
169 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) ::alb
170 REAL(KIND(1D0)),
DIMENSION(nsurf),
INTENT(in) ::emis
171 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(in) ::IceFrac
174 REAL(KIND(1D0)),
INTENT(in) ::DTIME
175 REAL(KIND(1D0)),
INTENT(in) ::ZENITH_deg
176 REAL(KIND(1D0)),
INTENT(in) ::tsurf_0
177 REAL(KIND(1D0)),
INTENT(in) ::kdown
178 REAL(KIND(1D0)),
INTENT(in) ::Temp_C
179 REAL(KIND(1D0)),
INTENT(in) ::RH
180 REAL(KIND(1D0)),
INTENT(in) ::Press_hPa
181 REAL(KIND(1D0)),
INTENT(in) ::qn1_obs
182 REAL(KIND(1D0)),
INTENT(in) ::SnowAlb
183 REAL(KIND(1D0)),
INTENT(in) ::NARP_TRANS_SITE
184 REAL(KIND(1D0)),
INTENT(in) ::NARP_EMIS_SNOW
186 INTEGER,
INTENT(in) ::nsurf
187 INTEGER,
INTENT(in) ::NetRadiationMethod_use
188 INTEGER,
INTENT(in) ::AlbedoChoice
189 INTEGER,
INTENT(in) ::ldown_option
190 INTEGER,
INTENT(in) ::DiagQN
192 REAL(KIND(1D0)),
INTENT(out) ::QSTARall
193 REAL(KIND(1D0)),
INTENT(out) ::QSTAR_SF
194 REAL(KIND(1D0)),
INTENT(out) ::QSTAR_S
195 REAL(KIND(1D0)),
INTENT(out) ::kclear
196 REAL(KIND(1D0)),
INTENT(out) ::KUPall
197 REAL(KIND(1D0)),
INTENT(out) ::LDOWN
198 REAL(KIND(1D0)),
INTENT(out) ::LUPall
199 REAL(KIND(1D0)),
INTENT(out) ::fcld
200 REAL(KIND(1D0)),
INTENT(out) ::TSURFall
202 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(out) ::qn1_ind_snow
203 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(out) ::kup_ind_snow
204 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(out) ::Tsurf_ind_snow
205 REAL(KIND(1d0)),
DIMENSION(nsurf),
INTENT(out) ::Tsurf_ind
207 REAL(KIND(1d0)),
INTENT(out) ::alb0
208 REAL(KIND(1d0)),
INTENT(out) ::alb1
210 REAL(KIND(1d0)),
DIMENSION(nsurf) ::qn1_ind
211 REAL(KIND(1d0)),
DIMENSION(nsurf) ::kup_ind
212 REAL(KIND(1d0)),
DIMENSION(nsurf) ::lup_ind
214 REAL(KIND(1d0)),
DIMENSION(nsurf) ::qn1_ind_nosnow
215 REAL(KIND(1d0)),
DIMENSION(nsurf) ::kup_ind_nosnow
216 REAL(KIND(1d0)),
DIMENSION(nsurf) ::lup_ind_nosnow
217 REAL(KIND(1d0)),
DIMENSION(nsurf) ::Tsurf_ind_nosnow
218 REAL(KIND(1d0)),
DIMENSION(nsurf) ::lup_ind_snow
220 REAL(KIND(1D0)) ::tsurf_0_K
221 REAL(KIND(1D0)) ::Temp_K, TD, ZENITH, QSTAR, QSTAR_SNOW, KUP_SNOW, LUP_SNOW, TSURF_SNOW, KUP, LUP, TSURF
222 REAL(KIND(1D0)) ::EMIS0, EMIS_A, TRANS
223 REAL(KIND(1D0)) ::LUPCORR, SIGMATK4, KDOWN_HR = 0.
226 REAL(KIND(1D0))::qn1_cum, kup_cum, lup_cum, tsurf_cum, &
227 qn1_is, kup_is, lup_is, tsurf_is, &
230 REAL(KIND(1D0)),
PARAMETER :: DEG2RAD = 0.017453292
232 REAL(KIND(1D0)),
PARAMETER :: SIGMA_SB = 5.67e-8
236 REAL(KIND(1D0)),
DIMENSION(365),
PARAMETER :: NARP_G = 3.0
241 tsurf_0_k = tsurf_0 + 273.16
242 temp_k = temp_c + 273.16
243 sigmatk4 = sigma_sb*temp_k**4
248 zenith = zenith_deg*deg2rad
250 IF (doy == 366) doy = 365
270 IF (sfr(is) /= 0) sf_all = sf_all + sfr(is)*(1 - snowfrac(is))
274 IF (diagqn == 1)
WRITE (*, *)
'is ', is
281 IF (albedochoice == 1 .AND. 180*zenith/acos(0.0) < 90)
THEN 282 alb0 = alb(is) + 0.5e-16*(180*zenith/acos(0.0))**8
289 IF ((ldown_option == 4) .OR. (ldown_option == 5))
THEN 290 IF (zenith < 1.5)
THEN 292 kclear =
isurface(doy, zenith)*trans*narp_trans_site
293 IF (kclear > 50.)
THEN 297 IF (ldown_option == 5)
THEN 306 IF (ldown_option == 4)
THEN 309 ELSEIF ((ldown_option == 5))
THEN 314 ELSEIF (ldown_option == 3)
THEN 317 ELSEIF (ldown_option == 2)
THEN 320 IF (diagqn == 1)
WRITE (*, *)
'ldown_option: ', ldown_option,
'FCLD:', fcld
322 IF (ldown_option > 1)
THEN 323 ldown = emis_a*sigmatk4
324 IF (diagqn == 1)
WRITE (*, *)
'EMIS_A: ', emis_a,
'SIGMATK4:', sigmatk4,
'LDOWN: ', ldown
330 IF (kdown_hr > 0)
THEN 331 lupcorr = (1 - alb0)*(0.08*kdown_hr)
338 if (netradiationmethod_use < 10)
then 340 tsurf = ((emis0*sigmatk4 + lupcorr)/(emis0*sigma_sb))**0.25
341 lup = emis0*sigmatk4 + lupcorr + (1 - emis0)*ldown
345 lup = emis0*sigma_sb*tsurf**4 + (1 - emis0)*ldown
349 qstar = kdown - kup + ldown - lup
350 tsurf = tsurf - 273.16
354 IF (snowfrac(is) > 0)
THEN 355 IF (albedochoice == 1 .AND. 180*zenith/acos(0.0) < 90)
THEN 356 alb1 = snowalb + 0.5e-16*(180*zenith/acos(0.0))**8
361 kup_snow = (alb1*(snowfrac(is) - snowfrac(is)*icefrac(is)) + alb0*snowfrac(is)*icefrac(is))*kdown
363 if (netradiationmethod_use < 10)
then 365 tsurf_snow = ((narp_emis_snow*sigmatk4)/(narp_emis_snow*sigma_sb))**0.25
370 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow**4 + (1 - narp_emis_snow)*ldown
373 tsurf_snow = tsurf_0_k
374 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow**4 + (1 - narp_emis_snow)*ldown
377 qstar_snow = kdown - kup_snow + ldown - lup_snow
378 tsurf_snow = tsurf_snow - 273.16
389 qn1_ind_nosnow(is) = qstar
390 kup_ind_nosnow(is) = kup
391 lup_ind_nosnow(is) = lup
392 tsurf_ind_nosnow(is) = tsurf
394 qn1_ind_snow(is) = qstar_snow
395 kup_ind_snow(is) = kup_snow
396 lup_ind_snow(is) = lup_snow
397 tsurf_ind_snow(is) = tsurf_snow
399 IF (sf_all /= 0)
THEN 400 qstar_sf = qstar_sf + qstar*sfr(is)*(1 - snowfrac(is))/sf_all
402 qstar_sf = qstar_sf + qstar*sfr(is)*(1 - snowfrac(is))
405 IF ((1 - sf_all) /= 0)
THEN 406 qstar_s = qstar_s + qstar_snow*sfr(is)*snowfrac(is)/(1 - sf_all)
408 qstar_s = qstar_s + qstar_snow*sfr(is)*snowfrac(is)
413 qn1_is = qstar*(1 - snowfrac(is)) + qstar_snow*snowfrac(is)
414 kup_is = kup*(1 - snowfrac(is)) + kup_snow*snowfrac(is)
415 lup_is = lup*(1 - snowfrac(is)) + lup_snow*snowfrac(is)
416 tsurf_is = tsurf*(1 - snowfrac(is)) + tsurf_snow*snowfrac(is)
418 IF (diagqn == 1)
WRITE (*, *)
'QSTAR', qstar,
'QSTAR_SNOW', qstar_snow,
'SnowFrac', snowfrac(is)
420 qn1_cum = qn1_cum + (qn1_is*sfr(is))
421 kup_cum = kup_cum + (kup_is*sfr(is))
422 lup_cum = lup_cum + (lup_is*sfr(is))
423 tsurf_cum = tsurf_cum + (tsurf_is*sfr(is))
428 tsurf_ind(is) = tsurf_is
430 IF (diagqn == 1)
WRITE (*, *)
'qn1_is: ', qn1_is
435 IF (netradiationmethod_use /= 3000)
THEN 455 IF (diagqn == 1)
WRITE (*, *)
'kdown: ', kdown,
'kup:', kup,
'LDOWN: ', ldown,
'LUP: ', lup
456 IF (diagqn == 1)
WRITE (*, *)
'Qn: ', qstarall
462 locationlatitude, locationlongitude, locationaltitude, &
463 sunazimuth, sunzenith)
466 REAL(KIND(1D0)),
INTENT(in) :: year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude
467 REAL(KIND(1D0)),
INTENT(out) ::sunazimuth, sunzenith
469 REAL(KIND(1D0)):: sec
470 INTEGER :: month, day, hour, min, seas, dayofyear, year_int
472 REAL(KIND(1D0)) :: juliancentury, julianday, julianephemeris_century, julianephemeris_day, &
473 julianephemeris_millenium
474 REAL(KIND(1D0)) :: earth_heliocentric_positionlatitude, earth_heliocentric_positionlongitude, &
475 earth_heliocentric_positionradius
476 REAL(KIND(1D0)) :: sun_geocentric_positionlatitude, sun_geocentric_positionlongitude
477 REAL(KIND(1D0)) :: nutationlongitude, nutationobliquity
478 REAL(KIND(1D0)) :: corr_obliquity
479 REAL(KIND(1D0)) :: aberration_correction
480 REAL(KIND(1D0)) :: apparent_sun_longitude
481 REAL(KIND(1D0)) :: apparent_stime_at_greenwich
482 REAL(KIND(1D0)) :: sun_rigth_ascension
483 REAL(KIND(1D0)) :: sun_geocentric_declination
484 REAL(KIND(1D0)) :: observer_local_hour
485 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension, topocentric_sun_positionrigth_ascension_parallax
486 REAL(KIND(1D0)) :: topocentric_sun_positiondeclination
487 REAL(KIND(1D0)) :: topocentric_local_hour
505 dayofyear = floor(idectime)
507 CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)
511 CALL julian_calculation(year, month, day, hour, min, sec, utc, juliancentury, julianday, julianephemeris_century, &
512 julianephemeris_day, julianephemeris_millenium)
517 &earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
521 & sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
534 & aberration_correction, apparent_sun_longitude)
538 &corr_obliquity, apparent_stime_at_greenwich)
542 &sun_rigth_ascension)
547 &sun_geocentric_declination)
555 &topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination, locationaltitude,&
556 &locationlatitude, observer_local_hour, sun_rigth_ascension, sun_geocentric_declination,&
557 &earth_heliocentric_positionradius)
561 & topocentric_local_hour)
565 & topocentric_local_hour, sunazimuth, sunzenith)
570 SUBROUTINE julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century&
571 &, julianephemeris_day, julianephemeris_millenium)
574 REAL(KIND(1D0)) :: A, B, D, delta_t
575 REAL(KIND(1D0)) :: juliancentury
576 REAL(KIND(1D0)) :: julianday
577 REAL(KIND(1D0)) :: julianephemeris_century
578 REAL(KIND(1D0)) :: julianephemeris_day
579 REAL(KIND(1D0)) :: julianephemeris_millenium
580 REAL(KIND(1D0)) :: M, sec, year, UTC
581 INTEGER :: day, hour, min, month
583 REAL(KIND(1D0)) :: ut_time, Y
589 IF (month == 1 .OR. month == 2)
THEN 596 ut_time = ((float(hour) - utc)/24.) + (float(min)/(60.*24.)) + (sec/(60.*60.*24.))
600 IF (year == 1582.)
THEN 601 IF (month == 10)
THEN 604 ELSE IF (day >= 15)
THEN 606 b = 2 - a + floor(a/4)
613 ELSE IF (month < 10)
THEN 617 b = 2 - a + floor(a/4)
620 ELSE IF (year < 1582.)
THEN 624 b = 2 - a + floor(a/4)
627 julianday = floor(365.25*(y + 4716.)) + floor(30.6001*(m + 1)) + d + b - 1524.5
630 julianephemeris_day = julianday + (delta_t/86400)
632 juliancentury = (julianday - 2451545.)/36525.
634 julianephemeris_century = (julianephemeris_day - 2451545.)/36525.
636 julianephemeris_millenium = julianephemeris_century/10.
641 &, earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
644 REAL(KIND(1D0)) :: julianephemeris_millenium
646 REAL(KIND(1D0)),
DIMENSION(64) :: A0
647 REAL(KIND(1D0)),
DIMENSION(34) :: A1
648 REAL(KIND(1D0)),
DIMENSION(20) :: A2
649 REAL(KIND(1D0)),
DIMENSION(7) :: A3
650 REAL(KIND(1D0)),
DIMENSION(3) :: A4
651 REAL(KIND(1D0)) :: A5
652 REAL(KIND(1D0)),
DIMENSION(64) :: B0
653 REAL(KIND(1D0)),
DIMENSION(34) :: B1
654 REAL(KIND(1D0)),
DIMENSION(20) :: B2
655 REAL(KIND(1D0)),
DIMENSION(7) :: B3
656 REAL(KIND(1D0)),
DIMENSION(3) :: B4
657 REAL(KIND(1D0)) :: B5
658 REAL(KIND(1D0)),
DIMENSION(64) :: C0
659 REAL(KIND(1D0)),
DIMENSION(34) :: C1
660 REAL(KIND(1D0)),
DIMENSION(20) :: C2
661 REAL(KIND(1D0)),
DIMENSION(7) :: C3
662 REAL(KIND(1D0)),
DIMENSION(3) :: C4
663 REAL(KIND(1D0)) :: C5
664 REAL(KIND(1D0)),
DIMENSION(40) :: A0j
665 REAL(KIND(1D0)),
DIMENSION(10) :: A1j
666 REAL(KIND(1D0)),
DIMENSION(6) :: A2j
667 REAL(KIND(1D0)),
DIMENSION(2) :: A3j
668 REAL(KIND(1D0)) :: A4j
669 REAL(KIND(1D0)),
DIMENSION(40) :: B0j
670 REAL(KIND(1D0)),
DIMENSION(10) :: B1j
671 REAL(KIND(1D0)),
DIMENSION(6) :: B2j
672 REAL(KIND(1D0)),
DIMENSION(2) :: B3j
673 REAL(KIND(1D0)) :: B4j
674 REAL(KIND(1D0)),
DIMENSION(40) :: C0j
675 REAL(KIND(1D0)),
DIMENSION(10) :: C1j
676 REAL(KIND(1D0)),
DIMENSION(6) :: C2j
677 REAL(KIND(1D0)),
DIMENSION(2) :: C3j
678 REAL(KIND(1D0)) :: C4j
679 REAL(KIND(1D0)),
DIMENSION(5) :: A0i
680 REAL(KIND(1D0)),
DIMENSION(5) :: B0i
681 REAL(KIND(1D0)),
DIMENSION(5) :: C0i
682 REAL(KIND(1D0)),
DIMENSION(2) :: A1i
683 REAL(KIND(1D0)),
DIMENSION(2) :: B1i
684 REAL(KIND(1D0)),
DIMENSION(2) :: C1i
685 REAL(KIND(1D0)) :: earth_heliocentric_positionlatitude
686 REAL(KIND(1D0)) :: earth_heliocentric_positionlongitude
687 REAL(KIND(1D0)) :: earth_heliocentric_positionradius
688 REAL(KIND(1D0)) :: JME
689 REAL(KIND(1D0)) :: L0
691 REAL(KIND(1D0)) :: L1
693 REAL(KIND(1D0)) :: L2
695 REAL(KIND(1D0)) :: L3
697 REAL(KIND(1D0)) :: L4
699 REAL(KIND(1D0)) :: L5
706 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
711 a0 = (/175347046, 3341656, 34894, 3497, 3418, 3136, 2676, 2343, 1324, 1273, 1199, 990, 902, 857, 780, 753, 505,&
712 &492, 357, 317, 284, 271, 243, 206, 205, 202, 156, 132, 126, 115, 103, 102, 102, 99, 98, 86, 85, 85, 80, 79, 71,&
713 &74, 74, 70, 62, 61, 57, 56, 56, 52, 52, 51, 49, 41, 41, 39, 37, 37, 36, 36, 33, 30, 30, 25/)
714 b0 = (/0., 4.669256800, 4.626100, 2.744100, 2.828900, 3.627700, 4.418100, 6.135200, 0.7425000, 2.037100,&
715 &1.109600, 5.233000, 2.045000, 3.508000, 1.179000, 2.533000, 4.583000, 4.205000, 2.92, 5.849000,&
716 &1.899000, 0.315, 0.345, 4.806000, 1.869000, 2.445800, 0.833, 3.411000, 1.083000, 0.645, 0.636, 0.976,&
717 &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,&
718 &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/)
719 c0 = (/0., 6283.075850, 12566.15170, 5753.384900, 3.523100, 77713.77150, 7860.419400, 3930.209700,&
720 &11506.76980, 529.6910, 1577.343500, 5884.927, 26.29800, 398.1490, 5223.694, 5507.553,&
721 &18849.22800, 775.5230, 0.067, 11790.62900, 796.2980, 10977.07900, 5486.778, 2544.314,&
722 &5573.143, 6069.777, 213.2990, 2942.463, 20.77500, 0.98, 4694.003, 15720.83900, 7.114000,&
723 &2146.170, 155.4200, 161000.6900, 6275.960, 71430.70, 17260.15, 12036.46, 5088.630, 3154.690,&
724 &801.8200, 9437.760, 8827.390, 7084.900, 6286.600, 14143.50, 6279.550, 12139.55, 1748.020,&
725 &5856.480, 1194.450, 8429.240, 19651.05, 10447.39, 10213.29, 1059.380, 2352.870, 6812.770,&
726 &17789.85, 83996.85, 1349.870, 4690.480/)
727 a1 = (/628331966747.000, 206059., 4303., 425., 119., 109., &
728 93., 72., 68., 67., 59., 56., 45., 36., 29., 21., 19., 19., 17., 16., &
729 16., 15., 12., 12., 12., 12., 11., 10., 10., 9., 9., 8., 6., 6./)
730 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,&
731 &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,&
733 c1 = (/0., 6283.075850, 12566.15170, 3.523000, 26.29800, 1577.344, 18849.23, 529.6900, 398.1500,&
734 &5507.550, 5223.690, 155.4200, 796.3000, 775.5200, 7.11, 0.98, 5486.780, 213.3000, 6275.960,&
735 &2544.310, 2146.170, 10977.08, 1748.020, 5088.630, 1194.450, 4694., 553.5700, 3286.600,&
736 &1349.870, 242.7300, 951.7200, 2352.870, 9437.760, 4690.480/)
737 a2 = (/52919, 8720, 309, 27, 16, 16, 10, 9, 7, 5, 4, 4, 3, 3, 3, 3, 3, 3, 2, 2/)
738 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,&
739 &6.12, 0.31, 2.28, 4.38, 3.75/)
740 c2 = (/0., 6283.075800, 12566.15200, 3.52, 26.3, 155.4200, 18849.23, 77713.77, 775.5200, 1577.340,&
741 &7.11, 5573.140, 796.3000, 5507.550, 242.7300, 529.6900, 398.1500, 553.5700, 5223.690, 0.98/)
742 a3 = (/289, 35, 17, 3, 1, 1, 1/)
743 b3 = (/5.8440, 0., 5.4900, 5.2000, 4.7200, 5.3000, 5.9700/)
744 c3 = (/6283.076, 0., 12566.15, 155.4200, 3.52, 18849.23, 242.7300/)
746 b4 = (/3.1420, 4.1300, 3.8400/)
747 c4 = (/0., 6283.08, 12566.15/)
752 jme = julianephemeris_millenium
755 l0 = sum(a0*cos(b0 + (c0*jme)))
756 l1 = sum(a1*cos(b1 + (c1*jme)))
757 l2 = sum(a2*cos(b2 + (c2*jme)))
758 l3 = sum(a3*cos(b3 + (c3*jme)))
759 l4 = sum(a4*cos(b4 + (c4*jme)))
760 l5 = a5*cos(b5 + (c5*jme))
762 earth_heliocentric_positionlongitude = &
763 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4) + (l5*jme**5))/1e8
765 earth_heliocentric_positionlongitude = earth_heliocentric_positionlongitude*180./pi
767 earth_heliocentric_positionlongitude =
set_to_range(earth_heliocentric_positionlongitude)
769 a0i = (/280, 102, 80, 44, 32/)
770 b0i = (/3.19900000000000, 5.42200000000000, 3.88000000000000, 3.70000000000000, 4./)
771 c0i = (/84334.6620000000, 5507.55300000000, 5223.69000000000, 2352.87000000000, 1577.34000000000/)
773 b1i = (/3.90000000000000, 1.73000000000000/)
774 c1i = (/5507.55000000000, 5223.69000000000/)
776 l0 = sum(a0i*cos(b0i + (c0i*jme)))
777 l1 = sum(a1i*cos(b1i + (c1i*jme)))
779 earth_heliocentric_positionlatitude = (l0 + (l1*jme))/1e8
781 earth_heliocentric_positionlatitude = earth_heliocentric_positionlatitude*180/pi
783 earth_heliocentric_positionlatitude =
set_to_range(earth_heliocentric_positionlatitude)
785 a0j = (/100013989, 1670700, 13956, 3084, 1628, 1576, 925, 542, 472, 346, 329, 307, 243, 212, 186, 175, 110,&
786 &98, 86, 86, 85, 63, 57, 56, 49, 47, 45, 43, 39, 38, 37, 37, 36, 35, 33, 32, 32, 28, 28, 26/)
787 b0j = (/0., 3.09846350, 3.05525000, 5.1985, 1.1739, 2.8469, 5.453, 4.564, 3.661, 0.964, 5.90, 0.299,&
788 &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,&
789 &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/)
790 c0j = (/0., 6283.07585, 12566.1517, 77713.7715, 5753.38490, 7860.41940, 11506.7700, 3930.21000,&
791 &5884.92700, 5507.55300, 5223.69400, 5573.14300, 11790.6290, 1577.34400, 10977.0790, 18849.2280,&
792 &5486.77800, 6069.78000, 15720.8400, 161000.690, 17260.1500, 529.69, 83996.8500, 71430.7000,&
793 &2544.31000, 775.52, 9437.76000, 6275.96000, 4694., 8827.39000, 19651.0500, 12139.5500,&
794 &12036.4600, 2942.46000, 7084.9, 5088.63000, 398.15, 6286.6, 6279.55000, 10447.3900/)
795 a1j = (/103019, 1721, 702, 32, 31, 25, 18, 10, 9, 9/)
796 b1j = (/1.10749000, 1.0644, 3.142, 1.02, 2.84, 1.32, 1.42, 5.91, 1.42, 0.270/)
797 c1j = (/6283.07585, 12566.1517, 0., 18849.2300, 5507.55000, 5223.69000, 1577.34000, 10977.0800,&
798 &6275.96000, 5486.78000/)
799 a2j = (/4359, 124, 12, 9, 6, 3/)
800 b2j = (/5.7846, 5.579, 3.14, 3.63, 1.87, 5.47/)
801 c2j = (/6283.07580, 12566.1520, 0., 77713.7700, 5573.14000, 18849./)
803 b3j = (/4.273, 3.92/)
804 c3j = (/6283.07600, 12566.1500/)
810 l0 = sum(a0j*cos(b0j + (c0j*jme)))
811 l1 = sum(a1j*cos(b1j + (c1j*jme)))
812 l2 = sum(a2j*cos(b2j + (c2j*jme)))
813 l3 = sum(a3j*cos(b3j + (c3j*jme)))
814 l4 = a4j*cos(b4j + (c4j*jme))
817 earth_heliocentric_positionradius = &
818 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4))/1e8
823 &earth_heliocentric_positionlatitude, sun_geocentric_positionlatitude, &
824 &sun_geocentric_positionlongitude)
827 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionlongitude
828 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionlatitude
829 REAL(KIND(1D0)) :: sun_geocentric_positionlatitude
830 REAL(KIND(1D0)) :: sun_geocentric_positionlongitude
834 sun_geocentric_positionlongitude = earth_heliocentric_positionlongitude + 180.0
836 sun_geocentric_positionlongitude =
set_to_range(sun_geocentric_positionlongitude)
838 sun_geocentric_positionlatitude = -earth_heliocentric_positionlatitude
840 sun_geocentric_positionlatitude =
set_to_range(sun_geocentric_positionlatitude)
846 REAL(KIND(1D0)),
INTENT(in) :: julianephemeris_century
847 REAL(KIND(1D0)),
DIMENSION(63) :: delta_longitude
848 REAL(KIND(1D0)),
DIMENSION(63) :: delta_obliquity
849 REAL(KIND(1D0)) :: JCE
850 REAL(KIND(1D0)) :: nutationlongitude
851 REAL(KIND(1D0)) :: nutationobliquity
852 REAL(KIND(1D0)),
DIMENSION(4) :: p0, p1, p2, p3, p4
853 REAL(KIND(1D0)),
DIMENSION(63) ::tabulated_argument
854 REAL(KIND(1D0)) :: X0
855 REAL(KIND(1D0)) :: X1
856 REAL(KIND(1D0)) :: X2
857 REAL(KIND(1D0)) :: X3
858 REAL(KIND(1D0)) :: X4
859 REAL(KIND(1D0)),
DIMENSION(5) :: Xi
860 INTEGER,
DIMENSION(315) :: Y_terms1
861 INTEGER,
DIMENSION(5, 63) ::Y_terms
862 REAL(KIND(1D0)),
DIMENSION(252) :: nutation_terms1
863 REAL(KIND(1D0)),
DIMENSION(4, 63) ::nutation_terms
865 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
871 jce = julianephemeris_century
874 p0 = (/(1/189474.), -0.0019142, 445267.11148, 297.85036/)
876 x0 = p0(1)*jce**3 + p0(2)*jce**2 + p0(3)*jce + p0(4)
879 p1 = (/-(1/300000.), -0.0001603, 35999.05034, 357.52772/)
881 x1 = p1(1)*jce**3 + p1(2)*jce**2 + p1(3)*jce + p1(4)
884 p2 = (/(1/56250.), 0.0086972, 477198.867398, 134.96298/)
886 x2 = p2(1)*jce**3 + p2(2)*jce**2 + p2(3)*jce + p2(4)
889 p3 = (/(1/327270.), -0.0036825, 483202.017538, 93.27191/)
891 x3 = p3(1)*jce**3 + p3(2)*jce**2 + p3(3)*jce + p3(4)
895 p4 = (/(1/450000.), 0.0020708, -1934.136261, 125.04452/)
897 x4 = p4(1)*jce**3 + p4(2)*jce**2 + p4(3)*jce + p4(4)
900 y_terms1 = (/0, 0, 0, 0, 1, -2, 0, 0, 2, 2, 0, 0, 0, 2, 2, 0, 0, &
901 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 2, 2, 0, 0, 0, 2, 1, &
902 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, &
904 1, 0, 1, 2, 0, -1, 2, 2, &
905 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, &
907 0, 0, -2, 0, 1, 2, 2, 0, &
908 0, 0, 2, 0, -2, 0, 0, 2, 0, 0, 0, -1, 2, 1, 0, 2, 0, 0, 0, &
909 2, 0, -1, 0, 1, -2, 2, 0, 2, 2, 0, 1, 0, 0, 1, &
910 -2, 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 0, &
911 2, -2, 0, 2, 0, -1, 2, 1, 2, 0, 1, 2, 2, 0, 1, 0, 2, 2, -2, &
912 1, 1, 0, 0, 0, -1, 0, 2, 2, 2, 0, 0, 2, 1, 2, &
913 0, 1, 0, 0, -2, 0, 2, 2, 2, -2, 0, 1, &
914 2, 1, 2, 0, -2, 0, 1, 2, 0, 0, 0, 1, 0, -1, 1, 0, 0, -2, -1, &
915 0, 2, 1, -2, 0, 0, 0, 1, 0, 0, 2, 2, 1, -2, &
916 0, 2, 0, 1, -2, 1, 0, 2, 1, 0, 0, 1, -2, &
917 0, -1, 0, 1, 0, 0, -2, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 2, &
918 0, 0, 0, -2, 2, 2, -1, -1, 1, 0, 0, 0, 1, &
919 1, 0, 0, 0, -1, 1, 2, 2, 2, -1, -1, 2, 2, &
920 0, 0, 3, 2, 2, 2, -1, 0, 2, 2/)
921 y_terms = reshape(y_terms1, (/5, 63/))
922 nutation_terms1 = (/-171996., -174.2, 92025., 8.9, -13187., -1.6, 5736., -3.1, -2274., -0.2, 977., -0.5, 2062., &
924 1426., -3.4, 54., -0.1, 712., 0.1, -7., 0., -517., 1.2, 224., -0.6, -386., -0.4, 200., 0., -301., &
926 217., -0.5, -95., 0.3, -158., 0., 0., 0., 129., 0.1, -70., 0., 123., 0., -53., 0., 63., 0., 0., 0., &
927 63., 0.1, -33., 0., -59., 0., 26., 0., &
928 -58., -0.1, 32., 0., -51., 0., 27., 0., 48., 0., 0., 0., 46., 0., -24., 0., -38., &
929 0., 16., 0., -31., 0., 13., 0., 29., 0., &
930 0., 0., 29., 0., -12., 0., 26., 0., 0., 0., -22., 0., 0., 0., 21., 0., -10., 0., 17., -0.1, 0., 0., 16., &
931 0., -8., 0., -16., 0.1, 7., 0., &
932 -15., 0., 9., 0., -13., 0., 7., 0., -12., 0., 6., 0., 11., 0., 0., 0., &
933 -10., 0., 5., 0., -8., 0., 3., 0., 7., &
934 0., -3., 0., -7., 0., 0., 0., &
935 -7., 0., 3., 0., -7., 0., 3., 0., 6., 0., 0., 0., 6., 0., -3., 0., 6., 0., &
936 -3., 0., -6., 0., 3., 0., -6., 0., 3., &
937 0., 5., 0., 0., 0., -5., 0., &
938 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., &
939 0., -4., 0., 0., 0., -4., 0., 0., 0., &
940 3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., 0., -3., 0., 0., &
941 0., -3., 0., 0., 0., -3., 0., 0., 0./)
942 nutation_terms = reshape(nutation_terms1, (/4, 63/))
945 xi = (/x0, x1, x2, x3, x4/)
948 tabulated_argument(i) = ((y_terms(1, i)*xi(1)) &
949 + (y_terms(2, i)*xi(2)) + (y_terms(3, i)*xi(3)) &
950 + (y_terms(4, i)*xi(4)) + (y_terms(5, i)*xi(5)))*pi/180
953 delta_longitude = ((nutation_terms(1, :) + (nutation_terms(2, :)*jce)))*sin(tabulated_argument)
954 delta_obliquity = ((nutation_terms(3, :) + (nutation_terms(4, :)*jce)))*cos(tabulated_argument)
957 nutationlongitude = sum(delta_longitude)/36000000.0
960 nutationobliquity = sum(delta_obliquity)/36000000.0
967 REAL(KIND(1D0)),
INTENT(out) :: corr_obliquity
968 REAL(KIND(1D0)),
INTENT(in) :: julianephemeris_millenium
969 REAL(KIND(1D0)),
INTENT(in) :: nutationobliquity
970 REAL(KIND(1D0)) :: mean_obliquity
971 REAL(KIND(1D0)),
DIMENSION(11) :: p
976 p = (/2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 84381.448/)
979 u = julianephemeris_millenium/10
981 &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 + &
982 &p(8)*u**3 + p(9)*u**2 + p(10)*u + p(11)
984 corr_obliquity = (mean_obliquity/3600) + nutationobliquity
990 REAL(KIND(1D0)),
INTENT(out) :: aberration_correction
991 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionradius
996 aberration_correction = -20.4898/(3600*earth_heliocentric_positionradius)
1001 & aberration_correction, apparent_sun_longitude)
1004 REAL(KIND(1D0)),
INTENT(in) :: aberration_correction
1005 REAL(KIND(1D0)),
INTENT(out) :: apparent_sun_longitude
1006 REAL(KIND(1D0)),
INTENT(in) :: nutationlongitude
1007 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlongitude
1011 apparent_sun_longitude = sun_geocentric_positionlongitude + nutationlongitude + aberration_correction
1016 & corr_obliquity, apparent_stime_at_greenwich)
1019 REAL(KIND(1D0)),
INTENT(out) :: apparent_stime_at_greenwich
1020 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1021 REAL(KIND(1D0)),
INTENT(in) :: julianday
1022 REAL(KIND(1D0)),
INTENT(in) :: juliancentury
1023 REAL(KIND(1D0)),
INTENT(in) :: nutationlongitude
1024 REAL(KIND(1D0)) :: JC
1025 REAL(KIND(1D0)) :: JD
1026 REAL(KIND(1D0)) :: mean_stime
1027 REAL(KIND(1D0)),
PARAMETER :: pi = 3.14159265358979d+0
1035 mean_stime = 280.46061837d+0 + (360.98564736629d+0*(jd - 2451545.0d+0)) + (0.000387933d+0*jc**2) - (jc**3/38710000.0d+0)
1040 apparent_stime_at_greenwich = mean_stime + (nutationlongitude*cos(corr_obliquity*pi/180))
1044 &sun_geocentric_positionlatitude, sun_rigth_ascension)
1047 REAL(KIND(1D0)),
INTENT(in) :: apparent_sun_longitude
1048 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1049 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlatitude
1050 REAL(KIND(1D0)),
INTENT(out) :: sun_rigth_ascension
1051 REAL(KIND(1D0)) :: argument_denominator
1052 REAL(KIND(1D0)) :: argument_numerator
1053 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1057 argument_numerator = (sin(apparent_sun_longitude*pi/180.0)*cos(corr_obliquity*pi/180.0)) - &
1058 (tan(sun_geocentric_positionlatitude*pi/180.0)*sin(corr_obliquity*pi/180.0))
1059 argument_denominator = cos(apparent_sun_longitude*pi/180.0)
1061 sun_rigth_ascension = atan2(argument_numerator, argument_denominator)*180.0/pi
1063 sun_rigth_ascension =
set_to_range(sun_rigth_ascension)
1067 &sun_geocentric_positionlatitude, sun_geocentric_declination)
1070 REAL(KIND(1D0)),
INTENT(in) :: apparent_sun_longitude
1071 REAL(KIND(1D0)),
INTENT(in) :: corr_obliquity
1072 REAL(KIND(1D0)),
INTENT(out) :: sun_geocentric_declination
1073 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_positionlatitude
1074 REAL(KIND(1D0)) :: argument
1075 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1077 argument = (sin(sun_geocentric_positionlatitude*pi/180.0)*cos(corr_obliquity*pi/180.0)) + &
1078 (cos(sun_geocentric_positionlatitude*pi/180.0)*sin(corr_obliquity*pi/180)*sin(apparent_sun_longitude*pi/180.0))
1080 sun_geocentric_declination = asin(argument)*180.0/pi
1084 &sun_rigth_ascension, observer_local_hour)
1087 REAL(KIND(1D0)),
INTENT(in) :: apparent_stime_at_greenwich
1088 REAL(KIND(1D0)),
INTENT(in) :: locationlongitude
1089 REAL(KIND(1D0)),
INTENT(out) :: observer_local_hour
1090 REAL(KIND(1D0)),
INTENT(in) :: sun_rigth_ascension
1092 observer_local_hour = apparent_stime_at_greenwich + locationlongitude - sun_rigth_ascension
1094 observer_local_hour =
set_to_range(observer_local_hour)
1098 &, topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination,&
1099 &locationaltitude, locationlatitude, observer_local_hour, sun_rigth_ascension,&
1100 &sun_geocentric_declination, earth_heliocentric_positionradius)
1103 REAL(KIND(1D0)),
INTENT(in) :: earth_heliocentric_positionradius
1104 REAL(KIND(1D0)),
INTENT(in) :: locationlatitude
1105 REAL(KIND(1D0)),
INTENT(in) :: locationaltitude
1106 REAL(KIND(1D0)),
INTENT(in) :: observer_local_hour
1107 REAL(KIND(1D0)),
INTENT(in) :: sun_geocentric_declination
1108 REAL(KIND(1D0)),
INTENT(in) :: sun_rigth_ascension
1109 REAL(KIND(1D0)) :: denominator
1110 REAL(KIND(1D0)) :: eq_horizontal_parallax
1111 REAL(KIND(1D0)) :: nominator
1112 REAL(KIND(1D0)) :: sun_rigth_ascension_parallax
1113 REAL(KIND(1D0)) :: topocentric_sun_positiondeclination
1114 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension
1115 REAL(KIND(1D0)) :: topocentric_sun_positionrigth_ascension_parallax
1116 REAL(KIND(1D0)) :: u
1117 REAL(KIND(1D0)) :: x
1118 REAL(KIND(1D0)) :: y
1119 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1126 eq_horizontal_parallax = 8.794/(3600*earth_heliocentric_positionradius)
1129 u = atan(0.99664719*tan(locationlatitude*pi/180))
1132 x = cos(u) + ((locationaltitude/6378140)*cos(locationaltitude*pi/180))
1135 y = (0.99664719d+0*sin(u)) + ((locationaltitude/6378140)*sin(locationlatitude*pi/180))
1138 nominator = -x*sin(eq_horizontal_parallax*pi/180.0)*sin(observer_local_hour*pi/180.0)
1139 denominator = cos(sun_geocentric_declination*pi/180.0) - &
1140 (x*sin(eq_horizontal_parallax*pi/180.0)*cos(observer_local_hour*pi/180.0))
1141 sun_rigth_ascension_parallax = atan2(nominator, denominator)
1143 topocentric_sun_positionrigth_ascension_parallax = sun_rigth_ascension_parallax*180.0/pi
1146 topocentric_sun_positionrigth_ascension = sun_rigth_ascension + (sun_rigth_ascension_parallax*180.0/pi)
1149 nominator = (sin(sun_geocentric_declination*pi/180.0) - (y*sin(eq_horizontal_parallax*pi/180.0)))&
1150 & *cos(sun_rigth_ascension_parallax)
1151 denominator = cos(sun_geocentric_declination*pi/180.0) - (y*sin(eq_horizontal_parallax*pi/180.0))&
1152 & *cos(observer_local_hour*pi/180.0)
1153 topocentric_sun_positiondeclination = atan2(nominator, denominator)*180.0/pi
1157 & topocentric_local_hour)
1160 REAL(KIND(1D0)),
INTENT(in) :: observer_local_hour
1161 REAL(KIND(1D0)),
INTENT(out) :: topocentric_local_hour
1162 REAL(KIND(1D0)),
INTENT(in) :: topocentric_sun_positionrigth_ascension_parallax
1166 topocentric_local_hour = observer_local_hour - topocentric_sun_positionrigth_ascension_parallax
1170 &topocentric_local_hour, sunazimuth, sunzenith)
1173 REAL(KIND(1D0)),
INTENT(in) :: locationlatitude
1174 REAL(KIND(1D0)),
INTENT(in) :: topocentric_local_hour
1175 REAL(KIND(1D0)),
INTENT(in) :: topocentric_sun_positiondeclination
1176 REAL(KIND(1D0)) :: corr_elevation
1177 REAL(KIND(1D0)) :: apparent_elevation
1178 REAL(KIND(1D0)) :: argument
1179 REAL(KIND(1D0)) :: denominator
1180 REAL(KIND(1D0)) :: nominator
1181 REAL(KIND(1D0)) :: refraction_corr
1182 REAL(KIND(1D0)) :: sunazimuth
1183 REAL(KIND(1D0)) :: sunzenith
1184 REAL(KIND(1D0)),
PARAMETER :: pi = 3.141592653589793d+0
1190 argument = (sin(locationlatitude*pi/180.0)*sin(topocentric_sun_positiondeclination*pi/180.0)) + &
1191 (cos(locationlatitude*pi/180.0)*cos(topocentric_sun_positiondeclination*pi/180.0)* &
1192 &cos(topocentric_local_hour*pi/180.0))
1193 corr_elevation = asin(argument)*180.0/pi
1196 argument = corr_elevation + (10.3/(corr_elevation + 5.11))
1197 refraction_corr = 1.02/(60*tan(argument*pi/180.0))
1204 apparent_elevation = corr_elevation + refraction_corr
1206 sunzenith = 90.0 - apparent_elevation
1211 nominator = sin(topocentric_local_hour*pi/180.0)
1212 denominator = (cos(topocentric_local_hour*pi/180.0)*sin(locationlatitude*pi/180.0)) - &
1213 (tan(topocentric_sun_positiondeclination*pi/180.0)*cos(locationlatitude*pi/180.0))
1214 sunazimuth = (atan2(nominator, denominator)*180.0/pi) + 180.0
1223 REAL(KIND(1D0)) :: max_interval
1224 REAL(KIND(1D0)) :: min_interval
1225 REAL(KIND(1D0)) :: var
1226 REAL(KIND(1D0)) :: vari
1228 max_interval = 360.0
1231 vari = var - max_interval*floor(var/max_interval)
1233 IF (vari < min_interval)
THEN 1234 vari = vari + max_interval
1240 FUNCTION dewpoint(Temp_C, rh)
RESULT(td)
1247 REAL(KIND(1d0))::rh, td, Temp_C, g
1249 g = ((17.27*temp_c)/(237.7 + temp_c)) + log(rh/100)
1250 td = (237.7*g)/(17.27 - g)
1254 FUNCTION prata_emis(Temp_K, EA_hPa)
RESULT(EMIS_A)
1256 REAL(KIND(1d0))::Temp_K, ea_hPa, EMIS_A
1259 w = 46.5*(ea_hpa/temp_k)
1260 emis_a = 1.-(1.+w)*exp(-sqrt(1.2 + 3.*w))
1263 FUNCTION emis_cloud(EMIS_A, FCLD)
RESULT(em_adj)
1265 REAL(KIND(1d0))::EMIS_A, FCLD, em_adj
1268 em_adj = emis_a + (1.-emis_a)*fcld
1273 REAL(KIND(1d0))::EMIS_A, FCLD, em_adj
1274 em_adj = emis_a + (1.-emis_a)*fcld*fcld
1278 REAL(KIND(1d0))::KDOWN, KCLEAR, FCLD
1280 fcld = 1.-kdown/kclear
1281 IF (fcld > 1.) fcld = 1.
1282 IF (fcld < 0.) fcld = 0.
1288 REAL(KIND(1d0)),
INTENT(in) :: RH
1289 REAL(KIND(1d0)),
INTENT(in) :: Temp
1291 REAL(KIND(1d0)) :: FWC
1292 REAL(KIND(1d0)) :: A, B
1299 b = 0.00019*temp + 0.015
1302 fwc = a*(exp(b*rh) - 1)
1303 IF (fwc > 1.) fwc = 1.
1304 IF (fwc < 0.) fwc = 0.
1328 FUNCTION isurface(doy, zenith)
RESULT(Isurf)
1332 REAL(KIND(1d0))::zenith, Isurf
1334 REAL(KIND(1d0))::Rmean, Rse, cosZ, Itoa
1335 REAL(KIND(1D0)),
PARAMETER :: DEG2RAD = 0.017453292
1339 IF (zenith < 90.*deg2rad)
THEN 1341 itoa = 1370.*(rmean/rse)**2
1353 REAL(KIND(1d0)) ::Rse
1354 REAL(KIND(1d0)) ::MA, nu, e, a
1359 ma = 2.*3.141592654*(doy - 3)/365.25463
1360 nu = ma + 0.0333988*sin(ma) + .0003486*sin(2.*ma) + 5e-6*sin(3.*ma)
1361 rse = a*(1 - e*e)/(1 + e*cos(nu))
1373 INTEGER :: lat, ios, ilat
1374 REAL(KIND(1d0)),
DIMENSION(365):: G
1382 READ (99, *, iostat=ios) ilat, g
1390 FUNCTION transmissivity(Press_hPa, Temp_C_dew, G, zenith)
RESULT(trans)
1398 REAL(KIND(1d0)) ::Press_hPa, TemP_C_dew, zenith, G, trans
1399 REAL(KIND(1d0))::m, TrTpg, u, Tw, Ta, cosZ
1400 REAL(KIND(1d0))::Tdf
1401 REAL(KIND(1D0)),
PARAMETER :: DEG2RAD = 0.017453292
1403 IF (zenith > 80.*deg2rad)
THEN 1404 cosz = cos(80.*deg2rad)
1409 tdf = temp_c_dew*1.8 + 32.
1411 m = 35*cosz/sqrt(1224.*cosz*cosz + 1)
1413 trtpg = 1.021 - 0.084*sqrt(m*(0.000949*press_hpa + 0.051))
1414 u = exp(0.113 - log(g + 1) + 0.0393*tdf)
1415 tw = 1 - 0.077*(u*m)**0.3
real(kind(1d0)) function transmissivity(Press_hPa, Temp_C_dew, G, zenith)
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)
subroutine nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)
subroutine narp_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
subroutine earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude, earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
real(kind(1d0)) function prata_emis(Temp_K, EA_hPa)
subroutine day2month(b, mb, md, seas, year, latitude)
subroutine radmethod(NetRadiationMethod, snowUse, NetRadiationMethod_use, AlbedoChoice, ldown_option)
real(kind(1d0)) function solar_esdist(doy)
real(kind(1d0)) function emis_cloud(EMIS_A, FCLD)
subroutine corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)
subroutine abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)
subroutine apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, corr_obliquity, apparent_stime_at_greenwich)
subroutine sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, sun_rigth_ascension)
subroutine observer_local_hour_calculation(apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)
subroutine apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude, aberration_correction, apparent_sun_longitude)
subroutine julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century, julianephemeris_day, julianephemeris_millenium)
subroutine sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude, sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
real(kind(1d0)) function, dimension(365) smithlambda(lat)
real(kind(1d0)) function emis_cloud_sq(EMIS_A, FCLD)
real(kind(1d0)) function wc_fraction(RH, Temp)
subroutine narp(nsurf, sfr, SnowFrac, alb, emis, IceFrac, NARP_TRANS_SITE, NARP_EMIS_SNOW, DTIME, ZENITH_deg, tsurf_0, kdown, Temp_C, RH, Press_hPa, qn1_obs, SnowAlb, AlbedoChoice, ldown_option, NetRadiationMethod_use, DiagQN, QSTARall, QSTAR_SF, QSTAR_S, kclear, KUPall, LDOWN, LUPall, fcld, TSURFall, qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_ind, alb0, alb1)
subroutine sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination, topocentric_local_hour, sunazimuth, sunzenith)
character(len=90) smithfile
real(kind(1d0)) function isurface(doy, zenith)
subroutine dectime_to_timevec(dectime, HOURS, MINS, SECS)
real(kind(1d0)) function cloud_fraction(KDOWN, KCLEAR)
real(kind(1d0)) function dewpoint(Temp_C, rh)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax, topocentric_local_hour)
real(kind(1d0)) function set_to_range(var)