SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
narp_module Module Reference

Functions/Subroutines

subroutine radmethod (netradiationmethod, snowuse, netradiationmethod_use, albedochoice, ldown_option)
 
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 narp_cal_sunposition (year, idectime, utc, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
 
subroutine julian_calculation (year, month, day, hour, min, sec, utc, juliancentury, julianday, julianephemeris_century, julianephemeris_day, julianephemeris_millenium)
 
subroutine earth_heliocentric_position_calculation (julianephemeris_millenium, earth_heliocentric_positionlatitude, earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
 
subroutine sun_geocentric_position_calculation (earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude, sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
 
subroutine nutation_calculation (julianephemeris_century, nutationlongitude, nutationobliquity)
 
subroutine corr_obliquity_calculation (julianephemeris_millenium, nutationobliquity, corr_obliquity)
 
subroutine abberation_correction_calculation (earth_heliocentric_positionradius, aberration_correction)
 
subroutine apparent_sun_longitude_calculation (sun_geocentric_positionlongitude, nutationlongitude, aberration_correction, apparent_sun_longitude)
 
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 sun_geocentric_declination_calculation (apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, sun_geocentric_declination)
 
subroutine observer_local_hour_calculation (apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)
 
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 topocentric_local_hour_calculate (observer_local_hour, topocentric_sun_positionrigth_ascension_parallax, topocentric_local_hour)
 
subroutine sun_topocentric_zenith_angle_calculate (locationlatitude, topocentric_sun_positiondeclination, topocentric_local_hour, sunazimuth, sunzenith)
 
real(kind(1d0)) function set_to_range (var)
 
real(kind(1d0)) function dewpoint_narp (temp_c, rh)
 
real(kind(1d0)) function prata_emis (temp_k, ea_hpa)
 
real(kind(1d0)) function emis_cloud (emis_a, fcld)
 
real(kind(1d0)) function emis_cloud_sq (emis_a, fcld)
 
real(kind(1d0)) function cloud_fraction (kdown, kclear)
 
real(kind(1d0)) function wc_fraction (rh, temp)
 
real(kind(1d0)) function isurface (doy, zenith)
 
real(kind(1d0)) function solar_esdist (doy)
 
real(kind(1d0)) function, dimension(365) smithlambda (lat)
 
real(kind(1d0)) function transmissivity (press_hpa, temp_c_dew, g, zenith)
 

Function/Subroutine Documentation

◆ abberation_correction_calculation()

subroutine narp_module::abberation_correction_calculation ( real(kind(1d0)), intent(in) earth_heliocentric_positionradius,
real(kind(1d0)), intent(out) aberration_correction )

Definition at line 1048 of file suews_phys_narp.f95.

1049 IMPLICIT NONE
1050
1051 REAL(KIND(1D0)), INTENT(out) :: aberration_correction
1052 REAL(KIND(1D0)), INTENT(in) :: earth_heliocentric_positionradius
1053
1054 ! This function compute the aberration_correction, as a function of the
1055 ! earth-sun distance.
1056
1057 aberration_correction = -20.4898/(3600*earth_heliocentric_positionradius)
1058

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ apparent_stime_at_greenwich_calculation()

subroutine narp_module::apparent_stime_at_greenwich_calculation ( real(kind(1d0)), intent(in) julianday,
real(kind(1d0)), intent(in) juliancentury,
real(kind(1d0)), intent(in) nutationlongitude,
real(kind(1d0)), intent(in) corr_obliquity,
real(kind(1d0)), intent(out) apparent_stime_at_greenwich )

Definition at line 1076 of file suews_phys_narp.f95.

1078 IMPLICIT NONE
1079
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
1089
1090 ! This function compute the apparent sideral time at Greenwich.
1091
1092 jd = julianday
1093 jc = juliancentury
1094
1095 ! Mean sideral time, in degrees
1096 mean_stime = 280.46061837d+0 + (360.98564736629d+0*(jd - 2451545.0d+0)) + (0.000387933d+0*jc**2) - (jc**3/38710000.0d+0)
1097
1098 ! Limit the range to [0-360];
1099 mean_stime = set_to_range(mean_stime)
1100
1101 apparent_stime_at_greenwich = mean_stime + (nutationlongitude*cos(corr_obliquity*pi/180))

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ apparent_sun_longitude_calculation()

subroutine narp_module::apparent_sun_longitude_calculation ( real(kind(1d0)), intent(in) sun_geocentric_positionlongitude,
real(kind(1d0)), intent(in) nutationlongitude,
real(kind(1d0)), intent(in) aberration_correction,
real(kind(1d0)), intent(out) apparent_sun_longitude )

Definition at line 1061 of file suews_phys_narp.f95.

1063 IMPLICIT NONE
1064
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
1069
1070 ! This function compute the sun apparent longitude
1071
1072 apparent_sun_longitude = sun_geocentric_positionlongitude + nutationlongitude + aberration_correction
1073

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ cloud_fraction()

real(kind(1d0)) function narp_module::cloud_fraction ( real(kind(1d0)) kdown,
real(kind(1d0)) kclear )

Definition at line 1338 of file suews_phys_narp.f95.

1339 REAL(KIND(1D0)) :: KDOWN, KCLEAR, FCLD
1340
1341 fcld = 1.-kdown/kclear
1342 IF (fcld > 1.) fcld = 1.
1343 IF (fcld < 0.) fcld = 0.

Referenced by narp().

Here is the caller graph for this function:

◆ corr_obliquity_calculation()

subroutine narp_module::corr_obliquity_calculation ( real(kind(1d0)), intent(in) julianephemeris_millenium,
real(kind(1d0)), intent(in) nutationobliquity,
real(kind(1d0)), intent(out) corr_obliquity )

Definition at line 1025 of file suews_phys_narp.f95.

1026 IMPLICIT NONE
1027
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
1034
1035 ! This function compute the true obliquity of the ecliptic.
1036
1037 p = (/2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 84381.448/)
1038 ! mean_obliquity = polyval(p, julian.ephemeris_millenium/10);
1039
1040 u = julianephemeris_millenium/10
1041 mean_obliquity =&
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)
1044
1045 corr_obliquity = (mean_obliquity/3600) + nutationobliquity

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ dewpoint_narp()

real(kind(1d0)) function narp_module::dewpoint_narp ( real(kind(1d0)) temp_c,
real(kind(1d0)) rh )

Definition at line 1301 of file suews_phys_narp.f95.

1302 ! ea = vapor pressure (hPa)
1303 ! td = dewpoint (oC)
1304 ! calculates dewpoint in degC from
1305 ! http://www.atd.ucar.edu/weather_fl/dewpoint.html
1306 ! dewpoint = (237.3 * ln(e_vp/6.1078)) / (17.27 - (ln(e_vp/6.1078)))
1307
1308 REAL(KIND(1D0)) :: rh, td, Temp_C, g
1309 !http://en.wikipedia.org/wiki/Dew_point
1310 g = ((17.27*temp_c)/(237.7 + temp_c)) + log(rh/100)
1311 td = (237.7*g)/(17.27 - g)
1312 !td = (237.3 * LOG(ea_hPa/6.1078)) / (17.27 - (LOG(ea_hPa/6.1078)))

Referenced by narp().

Here is the caller graph for this function:

◆ earth_heliocentric_position_calculation()

subroutine narp_module::earth_heliocentric_position_calculation ( real(kind(1d0)) julianephemeris_millenium,
real(kind(1d0)) earth_heliocentric_positionlatitude,
real(kind(1d0)) earth_heliocentric_positionlongitude,
real(kind(1d0)) earth_heliocentric_positionradius )

Definition at line 701 of file suews_phys_narp.f95.

703 IMPLICIT NONE
704
705 REAL(KIND(1D0)) :: julianephemeris_millenium
706
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
751 !REAL(KIND(1D0)) :: L0_terms !>
752 REAL(KIND(1D0)) :: L1
753 !REAL(KIND(1D0)) :: L1_terms !>
754 REAL(KIND(1D0)) :: L2
755 !REAL(KIND(1D0)) :: L2_terms !>
756 REAL(KIND(1D0)) :: L3
757 !REAL(KIND(1D0)) :: L3_terms !>
758 REAL(KIND(1D0)) :: L4
759 !REAL(KIND(1D0)) :: L4_terms !>
760 REAL(KIND(1D0)) :: L5
761 !REAL(KIND(1D0)), dimension(3) :: L5_terms !>
762 !REAL(KIND(1D0)) :: R0_terms !>
763 !REAL(KIND(1D0)) :: R1_terms !>
764 !REAL(KIND(1D0)) :: R2_terms !>
765 !REAL(KIND(1D0)) :: R3_terms !>
766 !REAL(KIND(1D0)), dimension(3) :: R4_terms !>
767 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793d+0
768
769 ! This function compute the earth position relative to the sun, using
770 ! tabulated values.
771
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,&
793 &5.3, 2.65, 4.67/)
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/)
806 a4 = (/114, 8, 1/)
807 b4 = (/3.1420, 4.1300, 3.8400/)
808 c4 = (/0., 6283.08, 12566.15/)
809 a5 = 1.
810 b5 = 3.1400
811 c5 = 0.
812
813 jme = julianephemeris_millenium
814
815 ! Compute the Earth Heliochentric longitude from the tabulated values.
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))
822
823 earth_heliocentric_positionlongitude = &
824 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4) + (l5*jme**5))/1e8
825 ! Convert the longitude to degrees.
826 earth_heliocentric_positionlongitude = earth_heliocentric_positionlongitude*180./pi
827 ! Limit the range to [0,360[;
828 earth_heliocentric_positionlongitude = set_to_range(earth_heliocentric_positionlongitude)
829
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/)
833 a1i = (/9, 6/)
834 b1i = (/3.90000000000000, 1.73000000000000/)
835 c1i = (/5507.55000000000, 5223.69000000000/)
836
837 l0 = sum(a0i*cos(b0i + (c0i*jme)))
838 l1 = sum(a1i*cos(b1i + (c1i*jme)))
839
840 earth_heliocentric_positionlatitude = (l0 + (l1*jme))/1e8
841 ! Convert the latitude to degrees.
842 earth_heliocentric_positionlatitude = earth_heliocentric_positionlatitude*180/pi
843 ! Limit the range to [0,360];
844 earth_heliocentric_positionlatitude = set_to_range(earth_heliocentric_positionlatitude)
845
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./)
863 a3j = (/145, 7/)
864 b3j = (/4.273, 3.92/)
865 c3j = (/6283.07600, 12566.1500/)
866 a4j = 4
867 b4j = 2.56
868 c4j = 6283.08000
869
870 ! Compute the Earth heliocentric radius vector
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))
876
877 ! Units are in AU
878 earth_heliocentric_positionradius = &
879 &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4))/1e8
880

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ emis_cloud()

real(kind(1d0)) function narp_module::emis_cloud ( real(kind(1d0)) emis_a,
real(kind(1d0)) fcld )

Definition at line 1324 of file suews_phys_narp.f95.

1325 !calculates adjusted emissivity due to clouds
1326 REAL(KIND(1D0)) :: EMIS_A, FCLD, em_adj
1327 !T. Loridan, removing the square for FCLD in the emissivity correction
1328 !em_adj=EMIS_A+(1.-EMIS_A)*FCLD*FCLD
1329 em_adj = emis_a + (1.-emis_a)*fcld

Referenced by narp().

Here is the caller graph for this function:

◆ emis_cloud_sq()

real(kind(1d0)) function narp_module::emis_cloud_sq ( real(kind(1d0)) emis_a,
real(kind(1d0)) fcld )

Definition at line 1332 of file suews_phys_narp.f95.

1333 !calculates adjusted emissivity due to clouds
1334 REAL(KIND(1D0)) :: EMIS_A, FCLD, em_adj
1335 em_adj = emis_a + (1.-emis_a)*fcld*fcld

Referenced by narp().

Here is the caller graph for this function:

◆ isurface()

real(kind(1d0)) function narp_module::isurface ( integer doy,
real(kind(1d0)) zenith )

Definition at line 1389 of file suews_phys_narp.f95.

1390 ! Calculates ground level solar irradiance clear sky
1391 ! assuming transmissivity = 1
1392 ! let it report zero if zenith >= 90
1393 REAL(KIND(1D0)) :: zenith, Isurf
1394 INTEGER :: doy
1395 REAL(KIND(1D0)) :: Rmean, Rse, cosZ, Itoa
1396 REAL(KIND(1D0)), PARAMETER :: DEG2RAD = 0.017453292
1397
1398 rmean = 149.6 !Stull 1998
1399 rse = solar_esdist(doy)
1400 IF (zenith < 90.*deg2rad) THEN
1401 cosz = cos(zenith)
1402 itoa = 1370.*(rmean/rse)**2 !top of the atmosphere
1403 isurf = itoa*cosz !ground level solar irradiance in W/m2
1404 ELSE
1405 isurf = 0.
1406 END IF
1407

References solar_esdist().

Referenced by narp().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ julian_calculation()

subroutine narp_module::julian_calculation ( real(kind(1d0)) year,
integer month,
integer day,
integer hour,
integer min,
real(kind(1d0)) sec,
real(kind(1d0)) utc,
real(kind(1d0)) juliancentury,
real(kind(1d0)) julianday,
real(kind(1d0)) julianephemeris_century,
real(kind(1d0)) julianephemeris_day,
real(kind(1d0)) julianephemeris_millenium )

Definition at line 631 of file suews_phys_narp.f95.

633 IMPLICIT NONE
634
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
643 !REAL(KIND(1D0)) :: time !>
644 REAL(KIND(1D0)) :: ut_time, Y !tt,
645 !
646 ! This function compute the julian day and julian century from the local
647 ! time and timezone information. Ephemeris are calculated with a delta_t=0
648 ! seconds.
649
650 IF (month == 1 .OR. month == 2) THEN
651 y = year - 1.
652 m = month + 12
653 ELSE
654 y = year
655 m = month
656 END IF
657 ut_time = ((float(hour) - utc)/24.) + (float(min)/(60.*24.)) + (sec/(60.*60.*24.)) ! time of day in UT time.
658 d = day + ut_time ! Day of month in decimal time, ex. 2sd day of month at 12:30:30UT, D=2.521180556
659
660 ! In 1582, the gregorian calendar was adopted
661 IF (year == 1582.) THEN
662 IF (month == 10) THEN
663 IF (day <= 4) THEN ! The Julian calendar ended on October 4, 1582
664 b = 0
665 ELSE IF (day >= 15) THEN ! The Gregorian calendar started on October 15, 1582
666 a = floor(y/100)
667 b = 2 - a + floor(a/4)
668 ELSE
669 !disp('This date never existed!. Date automatically set to October 4, 1582')
670 month = 10
671 day = 4
672 b = 0
673 END IF
674 ELSE IF (month < 10) THEN ! Julian calendar
675 b = 0
676 ELSE ! Gregorian calendar
677 a = floor(y/100)
678 b = 2 - a + floor(a/4)
679 END IF
680
681 ELSE IF (year < 1582.) THEN ! Julian calendar
682 b = 0
683 ELSE
684 a = floor(y/100) ! Gregorian calendar
685 b = 2 - a + floor(a/4)
686 END IF
687
688 julianday = floor(365.25*(y + 4716.)) + floor(30.6001*(m + 1)) + d + b - 1524.5
689
690 delta_t = 0. ! 33.184;
691 julianephemeris_day = julianday + (delta_t/86400)
692
693 juliancentury = (julianday - 2451545.)/36525.
694
695 julianephemeris_century = (julianephemeris_day - 2451545.)/36525.
696
697 julianephemeris_millenium = julianephemeris_century/10.
698

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ narp()

subroutine narp_module::narp ( integer, intent(in) storageheatmethod,
integer, intent(in) nsurf,
real(kind(1d0)), dimension(nsurf), intent(in) sfr_surf,
real(kind(1d0)), dimension(nsurf), intent(in) tsfc_surf,
real(kind(1d0)), dimension(nsurf), intent(in) snowfrac,
real(kind(1d0)), dimension(nsurf), intent(in) alb,
real(kind(1d0)), dimension(nsurf), intent(in) emis,
real(kind(1d0)), dimension(nsurf), intent(in) icefrac,
real(kind(1d0)), intent(in) narp_trans_site,
real(kind(1d0)), intent(in) narp_emis_snow,
real(kind(1d0)), intent(in) dtime,
real(kind(1d0)), intent(in) zenith_deg,
real(kind(1d0)), intent(in) tsurf_0,
real(kind(1d0)), intent(in) kdown,
real(kind(1d0)), intent(in) temp_c,
real(kind(1d0)), intent(in) rh,
real(kind(1d0)), intent(in) press_hpa,
real(kind(1d0)), intent(in) qn1_obs,
real(kind(1d0)), intent(in) ldown_obs,
real(kind(1d0)), intent(in) snowalb,
integer, intent(in) albedochoice,
integer, intent(in) ldown_option,
integer, intent(in) netradiationmethod_use,
integer, intent(in) diagqn,
real(kind(1d0)), dimension(nsurf), intent(out) qn_surf,
real(kind(1d0)), intent(out) qstarall,
real(kind(1d0)), intent(out) qstar_sf,
real(kind(1d0)), intent(out) qstar_s,
real(kind(1d0)), intent(out) kclear,
real(kind(1d0)), intent(out) kupall,
real(kind(1d0)), intent(out) ldown,
real(kind(1d0)), intent(out) lupall,
real(kind(1d0)), intent(out) fcld,
real(kind(1d0)), intent(out) tsurfall,
real(kind(1d0)), dimension(nsurf), intent(out) qn1_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(out) kup_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(out) tsurf_ind_snow,
real(kind(1d0)), dimension(nsurf), intent(out) tsurf_surf,
real(kind(1d0)), intent(out) albedo_snowfree,
real(kind(1d0)), intent(out) albedo_snow )

Definition at line 121 of file suews_phys_narp.f95.

133 !KCLEAR,FCLD,DTIME,KDOWN,QSTARall,KUPall,LDOWN,LUPall,TSURFall,&
134 !AlbedoChoice,ldown_option,Temp_C,Press_hPa,Ea_hPa,qn1_obs,RH,&
135 !,zenith_degnetRadiationChoice,
136
137 !SUBROUTINE NARP(QSTAR,DTIME,KDOWN,LDOWN,T,RH,PRESS,FCLD,SNOW)
138 !returns estimate of Q* given the meteorological fields and prior
139 !configuration call.
140
141 !OUTPUT FIELDS
142 !QSTARall (W/m2) = net all wave radiation
143 !KCLEAR = clear sky incoming solar radiation
144 !KUPall = reflect solar radiation
145 !LDOWN = incoming longwave radiation (observed or modelled depending on ldown_option)
146 !LUPall = outgoing longwaver radiation
147 !FCLD = estimated cloud fraction (used only for emissivity estimate)
148 ! FCLD USED AS INPUT ALSO
149 !TSURFall (DEG C)= estimated surface temperature
150 !QSTAR_SF = net all wave radiation for snow free surface
151 !QSTAR_S = net all wave radiation for SnowPack
152
153 !INPUT FIELDS
154 !DTIME (days) = local time, not daylight savings
155 !KDOWN (W/m2) = incoming solar radiation
156 !T (DEG C) = air temperature near model height
157 !RH (%) = relative humidity near model height
158 !PRESS (mb) = station pressure, use estimate if unavailable
159 !qn1_obs = Observed Q*
160 !ldown_obs = Observed Ldown
161
162 !INTERNAL FIELDS
163 ! TemP_K = air temperature in K
164 ! ES_hPs = saturation vapor pressure (hPa)
165 ! EA_Pa = vapor pressure (hPa)
166 ! TD = dewpoint (C)
167 ! ZENITH = solar zenith angle
168 ! albedo_snowfree = surface albedo
169 ! EMIS0 = surface emissivity
170 ! EMIS_A = atmospheric emissivity
171 ! TRANS = atmospheric transmissivity
172 ! LUPCORR = correction for surface heating by kdown (W/m2)
173 ! SIGMATK4 = energy flux density
174 ! KDOWN_HR = hourly average insolation
175 ! DOY = day of year
176
177 !Modified by LJ to calcuate snow free and SnowPack components (May 2013)
178 !Modified to define variables in data_in module
179 !-------------------------------------------------------------------------------
180 ! USE allocateArray
181 ! use gis_data
182 ! use data_in ! Included 20140701, FL
183 ! use moist ! Included 20140701, FL
184 ! use time ! Included 20140701, FL
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
191 ! REAL(KIND(1D0)),DIMENSION(365),INTENT(in) ::NARP_G
192
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
205
206 INTEGER, INTENT(in) :: nsurf
207 INTEGER, INTENT(in) :: NetRadiationMethod_use ! the one processed by RadMethod
208 INTEGER, INTENT(in) :: storageheatmethod ! needed for separate surface temperatures
209 INTEGER, INTENT(in) :: AlbedoChoice ! flag if correction to albedo of snow cover should be applied
210 INTEGER, INTENT(in) :: ldown_option ! flag for different ldown modelling options; 1 for obs; see code below for other parameterisations
211 INTEGER, INTENT(in) :: DiagQN
212
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
222
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
228
229 REAL(KIND(1D0)), INTENT(out) :: albedo_snowfree
230 REAL(KIND(1D0)), INTENT(out) :: albedo_snow
231
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
235
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
241
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 !,RH,DTIME,KDOWN
246 REAL(KIND(1D0)) :: LUPCORR, SIGMATK4, KDOWN_HR
247 INTEGER :: DOY, is
248
249 REAL(KIND(1D0)) :: qn1_cum, kup_cum, lup_cum, tsurf_cum, & !Cumulative radiation components
250 qn_is, kup_is, lup_is, tsurf_is, & !Sub-surface radiation components
251 sf_all
252
253 REAL(KIND(1D0)), PARAMETER :: DEG2RAD = 0.017453292
254 ! REAL(KIND(1D0)),PARAMETER ::RAD2DEG=57.29577951
255 REAL(KIND(1D0)), PARAMETER :: SIGMA_SB = 5.67e-8
256
257 ! NB: NARP_G is not assigned with a value in SUEWS_translate.
258 ! 3.0 is used here as annual average for mid-latitude areas. TS 24 Oct 2017
259 REAL(KIND(1D0)), DIMENSION(365), PARAMETER :: NARP_G = 3.0
260 !Initialize variables
261 ! RH=avrh
262 ! DTIME=dectime
263 ! KDOWN=avkdn
264 kdown_hr = 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
269 td = dewpoint_narp(temp_c, rh)
270 ! Sun postition is now calculated in the main loop, FL
271 !ZENITH=SOLAR_ZENITH(NARP_LAT,NARP_LONG,NARP_TZ,DTIME)
272 !call NARP_cal_SunPosition(NARP_YEAR,DTIME,NARP_TZ,NARP_LAT,NARP_LONG,Alt,AZIMUTH,ZENITH)
273 zenith = zenith_deg*deg2rad
274 doy = int(dtime)
275 IF (doy == 366) doy = 365
276
277 !===================================================================================
278 !Calculate radiation for each sub-surface
279
280 ! initialisation
281 qn1_cum = 0
282 kup_cum = 0
283 lup_cum = 0
284 tsurf_cum = 0
285
286 qstarall = 0
287 qstar_sf = 0
288 qstar_s = 0
289 kclear = 0
290 kupall = 0
291 IF (ldown_option .NE. 1) THEN
292 ldown = 0
293 END IF
294 lupall = 0
295 fcld = 0
296 tsurfall = 0
297
298 qn1_ind_snow = 0
299 kup_ind_snow = 0
300 lup_ind_snow = 0
301 tsurf_ind_snow = 0
302 tsurf_surf = 0
303 albedo_snowfree = 0.2 ! arbitrary non-zero value for initialisatoin
304 albedo_snow = snowalb
305
306 IF (ldown_option == 1) THEN !Observed ldown provided as forcing
307 ldown = ldown_obs
308 ELSE
309 ldown = 0 !to be filled in NARP
310 END IF
311
312 !Total snowfree surface fraction
313 sf_all = 0
314 DO is = 1, nsurf
315 IF (sfr_surf(is) /= 0) sf_all = sf_all + sfr_surf(is)*(1 - snowfrac(is))
316 END DO
317
318 !looping over the surfaces
319 DO is = 1, nsurf
320 IF (diagqn == 1) WRITE (*, *) 'is ', is
321
322 emis_a = prata_emis(temp_k, press_hpa)
323
324 !--------------------------------------------------
325 !-------SNOW FREE SURFACE--------------------------
326
327 !NARP
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 !AIDA 1982
330 ELSE
331 albedo_snowfree = alb(is)
332 END IF
333 emis0 = emis(is)
334
335 !Downward longwave radiation
336 IF ((ldown_option == 4) .OR. (ldown_option == 5)) THEN !Estimate FCLD from Kdown (Offerle et al. 2003)
337 IF (zenith < 1.5) THEN !DAYTIME CALCULATIONS
338 trans = transmissivity(press_hpa, td, narp_g(doy), zenith)
339 kclear = isurface(doy, zenith)*trans*narp_trans_site
340 IF (kclear > 50.) THEN
341 fcld = cloud_fraction(kdown, kclear)
342 emis_a = emis_cloud_sq(emis_a, fcld)
343 ELSE
344 IF (ldown_option == 5) THEN ! Use RH when Kdown can not be used
345 fcld = wc_fraction(rh, temp_c)
346 emis_a = emis_cloud(emis_a, fcld)
347 ELSE
348 !FCLD is left to the latest calculable value
349 emis_a = emis_cloud_sq(emis_a, fcld)
350 END IF
351 END IF
352 ELSE !NIGHT TIME CALCULATIONS
353 IF (ldown_option == 4) THEN
354 !FCLD is left to the latest calculable value
355 emis_a = emis_cloud_sq(emis_a, fcld)
356 ELSEIF ((ldown_option == 5)) THEN ! Use RH
357 fcld = wc_fraction(rh, temp_c)
358 emis_a = emis_cloud(emis_a, fcld)
359 END IF
360 END IF
361 ELSEIF (ldown_option == 3) THEN !Always use RH
362 fcld = wc_fraction(rh, temp_c)
363 emis_a = emis_cloud(emis_a, fcld)
364 ELSEIF (ldown_option == 2) THEN ! FCLD obs, came from input
365 emis_a = emis_cloud(emis_a, fcld)
366 END IF
367 IF (diagqn == 1) WRITE (*, *) 'ldown_option: ', ldown_option, 'FCLD:', fcld
368
369 IF (ldown_option > 1) THEN ! Forcing available if ldown_option=1, model otherwise
370 ldown = emis_a*sigmatk4
371 IF (diagqn == 1) WRITE (*, *) 'EMIS_A: ', emis_a, 'SIGMATK4:', sigmatk4, 'LDOWN: ', ldown
372 END IF
373
374 !Upward longwave radiation
375 !NARP
376 !Note that this is not averaged over the hour for cases where time step < 1hr
377 kdown_hr = kdown
378 !calculate Lup correction using Eq. 15 of Offerle et al.
379 IF (kdown_hr > 0) THEN
380 lupcorr = (1 - albedo_snowfree)*(0.08*kdown_hr)
381 ELSE
382 lupcorr = 0.
383 END IF
384 !calculate Lup
385 IF (netradiationmethod_use < 10) THEN
386 ! NARP method
387 tsurf_k = ((emis0*sigmatk4 + lupcorr)/(emis0*sigma_sb))**0.25 !Eqs. (14) and (15),
388 lup = emis0*sigmatk4 + lupcorr + (1 - emis0)*ldown !Eq (16) in Offerle et al. (2002)
389 ELSE
390 ! use iteration-based approach to calculate LUP and also TSURF; TS 20 Sep 2019
391 IF (storageheatmethod == 5) THEN
392 tsurf_k = tsfc_surf_k(is)
393 ELSE
394 tsurf_k = tsurf_0_k
395 END IF
396 ! TSURF = tsfc_surf_K(is)
397 lup = emis0*sigma_sb*tsurf_k**4 + (1 - emis0)*ldown
398 END IF
399
400 !Upward shortwave radiation of the snowfree surface fractions
401 kup = albedo_snowfree*kdown
402
403 !Net all wave
404 qstar = kdown - kup + ldown - lup
405 tsurf_c = tsurf_k - 273.16
406
407 !Define sub-surface radiation components
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
412
413 !!!!!!!!!!!!!!!! snow section !!!!!!!!!!!!!!!!!!!!!
414 !======================================================================
415 !Snow related parameters if snow pack existing
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 !AIDA 1982
419 ELSE
420 albedo_snow = snowalb
421 END IF
422
423 kup_snow = (albedo_snow*(snowfrac(is) - snowfrac(is)*icefrac(is)) + albedo_snowfree*snowfrac(is)*icefrac(is))*kdown
424
425 IF (netradiationmethod_use < 10) THEN
426 ! NARP method
427 tsurf_snow_k = ((narp_emis_snow*sigmatk4)/(narp_emis_snow*sigma_sb))**0.25 !Snow surface temperature
428 !IF (TSURF_SNOW>273.16) TSURF_SNOW=min(273.16,Temp_K)!Set this to 2 degrees (melted water on top)
429 !open(34,file='TestingSnowFrac.txt',position='append')
430 !write(34,*) dectime,is,albedo_snow,albedo_snowfree,SnowFrac(is),IceFrac(is),KDOWN,KUP_snow
431 !close(34)
432 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow_k**4 + (1 - narp_emis_snow)*ldown
433 ELSE
434 ! use iteration-based approach to calculate LUP and also TSURF; TS 20 Sep 2019
435 tsurf_snow_k = tsurf_0_k
436 lup_snow = narp_emis_snow*sigma_sb*tsurf_snow_k**4 + (1 - narp_emis_snow)*ldown
437 END IF
438
439 qstar_snow = kdown - kup_snow + ldown - lup_snow
440 tsurf_snow_c = tsurf_snow_k - 273.16
441
442 ELSE
443 kup_snow = 0
444 lup_snow = 0
445 tsurf_snow_c = 0
446 qstar_snow = 0
447 tsurf_snow_k = 0
448 !QSTAR_ICE = 0
449 !KUP_ICE = 0
450 END IF
451
452 !Define snow sub-surface radiation components
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
457
458 IF (sf_all /= 0) THEN
459 qstar_sf = qstar_sf + qstar*sfr_surf(is)*(1 - snowfrac(is))/sf_all
460 ELSE
461 qstar_sf = qstar_sf + qstar*sfr_surf(is)*(1 - snowfrac(is))
462 END IF
463
464 IF ((1 - sf_all) /= 0) THEN
465 qstar_s = qstar_s + qstar_snow*sfr_surf(is)*snowfrac(is)/(1 - sf_all)
466 ELSE
467 qstar_s = qstar_s + qstar_snow*sfr_surf(is)*snowfrac(is)
468 END IF
469 !!!!!!!!!!!!!!!! end snow section !!!!!!!!!!!!!!!!!!!!!
470
471 !---------------------------------------------------------------------
472 !Calculate snow weighted variables for each subsurface
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)
477
478 IF (diagqn == 1) WRITE (*, *) 'QSTAR', qstar, 'QSTAR_SNOW', qstar_snow, 'SnowFrac', snowfrac(is)
479
480 !Add up the radiation from each surface
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))
485
486 !Define sub-surface radiation components
487 qn_surf(is) = qn_is
488 kup_surf(is) = kup_is
489 lup_surf(is) = lup_is
490 tsurf_surf(is) = tsurf_is
491
492 IF (diagqn == 1) WRITE (*, *) 'qn1_is: ', qn_is
493
494 END DO !End of the surface types
495
496 !Set overall radiation components
497 IF (netradiationmethod_use /= 3000) THEN !Observed Q* used and snow is modeled
498 qstarall = qn1_cum
499 ELSE
500 qstarall = qn1_obs
501 END IF
502 kupall = kup_cum
503 lupall = lup_cum
504 tsurfall = tsurf_cum
505
506 ! qn1=QSTARall
507 ! kup=KUPall
508 !LDOWN has same name
509 ! LUP=LUPall
510 tsurf_k = tsurfall
511 !if (kup>500) then
512 ! write(*,*) Kdown, kup, kup_ind(1),kup_ind(2),kup_ind(3),kup_ind(4),kup_ind(5),kup_ind(6),SnowAlb
513 ! pause
514 !endif
515
516 IF (diagqn == 1) WRITE (*, *) 'kdown: ', kdown, 'kup:', kup, 'LDOWN: ', ldown, 'LUP: ', lup
517 IF (diagqn == 1) WRITE (*, *) 'Qn: ', qstarall
518

References cloud_fraction(), dewpoint_narp(), emis_cloud(), emis_cloud_sq(), isurface(), prata_emis(), transmissivity(), and wc_fraction().

Referenced by suews_driver::suews_cal_qn(), and suews_driver::suews_cal_qn_dts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ narp_cal_sunposition()

subroutine narp_module::narp_cal_sunposition ( real(kind(1d0)), intent(in) year,
real(kind(1d0)), intent(in) idectime,
real(kind(1d0)), intent(in) utc,
real(kind(1d0)), intent(in) locationlatitude,
real(kind(1d0)), intent(in) locationlongitude,
real(kind(1d0)), intent(in) locationaltitude,
real(kind(1d0)), intent(out) sunazimuth,
real(kind(1d0)), intent(out) sunzenith )

Definition at line 522 of file suews_phys_narp.f95.

525 IMPLICIT NONE
526
527 REAL(KIND(1D0)), INTENT(in) :: year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude
528 REAL(KIND(1D0)), INTENT(out) :: sunazimuth, sunzenith
529
530 REAL(KIND(1D0)) :: sec
531 INTEGER :: month, day, hour, min, seas, dayofyear, year_int
532
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
549 ! REAL(KIND(1D0)) :: sunazimuth,sunzenith
550
551 ! This function compute the sun position (zenith and azimuth angle (in degrees) at the observer
552 ! location) as a function of the observer local time and position.
553 !
554 ! Input lat and lng should be in degrees, alt in meters.
555 !
556 ! It is an implementation of the algorithm presented by Reda et Andreas in:
557 ! Reda, I., Andreas, A. (2003) Solar position algorithm for solar
558 ! radiation application. National Renewable Energy Laboratory (NREL)
559 ! Technical report NREL/TP-560-34302.
560 ! This document is available at www.osti.gov/bridge
561 ! Code is translated from matlab code by Fredrik Lindberg (fredrikl@gvc.gu.se)
562 ! Last modified: LJ 27 Jan 2016 - Tabs removed
563
564 ! Convert to timevectors from dectime and year
565 CALL dectime_to_timevec(idectime, hour, min, sec)
566 dayofyear = floor(idectime)
567 year_int = int(year)
568 CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)
569
570 ! 1. Calculate the Julian Day, and Century. Julian Ephemeris day, century
571 ! and millenium are calculated using a mean delta_t of 33.184 seconds.
572 CALL julian_calculation(year, month, day, hour, min, sec, utc, juliancentury, julianday, julianephemeris_century, &
573 julianephemeris_day, julianephemeris_millenium)
574
575 ! 2. Calculate the Earth heliocentric longitude, latitude, and radius
576 ! vector (L, B, and R)
577 CALL earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude,&
578 &earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
579
580 ! 3. Calculate the geocentric longitude and latitude
581 CALL sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude,&
582 & sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
583
584 ! 4. Calculate the nutation in longitude and obliquity (in degrees).
585 CALL nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)
586
587 ! 5. Calculate the true obliquity of the ecliptic (in degrees).
588 CALL corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)
589
590 ! 6. Calculate the aberration correction (in degrees)
591 CALL abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)
592
593 ! 7. Calculate the apparent sun longitude in degrees)
594 CALL apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude,&
595 & aberration_correction, apparent_sun_longitude)
596
597 ! 8. Calculate the apparent sideral time at Greenwich (in degrees)
598 CALL apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, &
599 &corr_obliquity, apparent_stime_at_greenwich)
600
601 ! 9. Calculate the sun rigth ascension (in degrees)
602 CALL sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
603 &sun_rigth_ascension)
604
605 ! 10. Calculate the geocentric sun declination (in degrees). Positive or
606 ! negative if the sun is north or south of the celestial equator.
607 CALL sun_geocentric_declination_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
608 &sun_geocentric_declination)
609
610 ! 11. Calculate the observer local hour angle (in degrees, westward from south).
611 CALL observer_local_hour_calculation(apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)
612
613 ! 12. Calculate the topocentric sun position (rigth ascension, declination and
614 ! rigth ascension parallax in degrees)
615 CALL topocentric_sun_position_calculate(topocentric_sun_positionrigth_ascension,&
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)
619
620 ! 13. Calculate the topocentric local hour angle (in degrees)
621 CALL topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax,&
622 & topocentric_local_hour)
623
624 ! 14. Calculate the topocentric zenith and azimuth angle (in degrees)
625 CALL sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination,&
626 & topocentric_local_hour, sunazimuth, sunzenith)
627
dectime_to_timevec
subroutine dectime_to_timevec(dectime, hours, mins, secs)
Definition suews_util_time.f95:203
day2month
subroutine day2month(b, mb, md, seas, year, latitude)
Definition suews_util_time.f95:24

References abberation_correction_calculation(), apparent_stime_at_greenwich_calculation(), apparent_sun_longitude_calculation(), corr_obliquity_calculation(), day2month(), dectime_to_timevec(), earth_heliocentric_position_calculation(), julian_calculation(), nutation_calculation(), observer_local_hour_calculation(), sun_geocentric_declination_calculation(), sun_geocentric_position_calculation(), sun_rigth_ascension_calculation(), sun_topocentric_zenith_angle_calculate(), topocentric_local_hour_calculate(), and topocentric_sun_position_calculate().

Referenced by beers_module::cal_ci_latenight(), metdisagg::disaggregatemet(), suews_driver::suews_cal_main(), suews_driver::suews_cal_main_dts(), and beers_module::tsurfbeers().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ nutation_calculation()

subroutine narp_module::nutation_calculation ( real(kind(1d0)), intent(in) julianephemeris_century,
real(kind(1d0)) nutationlongitude,
real(kind(1d0)) nutationobliquity )

Definition at line 904 of file suews_phys_narp.f95.

905 IMPLICIT NONE
906
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
925 INTEGER :: i
926 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793d+0
927
928 ! This function compute the nutation in longtitude and in obliquity, in
929 ! degrees.
930
931 ! All Xi are in degrees.
932 jce = julianephemeris_century
933
934 ! 1. Mean elongation of the moon from the sun
935 p0 = (/(1/189474.), -0.0019142, 445267.11148, 297.85036/)
936 ! X0 = polyval(p, JCE);
937 x0 = p0(1)*jce**3 + p0(2)*jce**2 + p0(3)*jce + p0(4) ! This is faster than polyval...
938
939 ! 2. Mean anomaly of the sun (earth)
940 p1 = (/-(1/300000.), -0.0001603, 35999.05034, 357.52772/)
941 ! X1 = polyval(p, JCE);
942 x1 = p1(1)*jce**3 + p1(2)*jce**2 + p1(3)*jce + p1(4)
943
944 ! 3. Mean anomaly of the moon
945 p2 = (/(1/56250.), 0.0086972, 477198.867398, 134.96298/)
946 ! X2 = polyval(p, JCE);
947 x2 = p2(1)*jce**3 + p2(2)*jce**2 + p2(3)*jce + p2(4)
948
949 ! 4. Moon argument of latitude
950 p3 = (/(1/327270.), -0.0036825, 483202.017538, 93.27191/)
951 ! X3 = polyval(p, JCE);
952 x3 = p3(1)*jce**3 + p3(2)*jce**2 + p3(3)*jce + p3(4)
953
954 ! 5. Longitude of the ascending node of the moon's mean orbit on the
955 ! ecliptic, measured from the mean equinox of the date
956 p4 = (/(1/450000.), 0.0020708, -1934.136261, 125.04452/)
957 ! X4 = polyval(p, JCE);
958 x4 = p4(1)*jce**3 + p4(2)*jce**2 + p4(3)*jce + p4(4)
959
960 ! Y tabulated terms from the original code
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, &
964 0, 0, 0, 0, 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, &
967 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., &
984 0.2, -895., 0.5, &
985 1426., -3.4, 54., -0.1, 712., 0.1, -7., 0., -517., 1.2, 224., -0.6, -386., -0.4, 200., 0., -301., &
986 0., 129., -0.1, &
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/))
1004 ! Using the tabulated values, compute the delta_longitude and
1005 ! delta_obliquity.
1006 xi = (/x0, x1, x2, x3, x4/)
1007
1008 DO i = 1, 63
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
1012 END DO
1013
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)
1016
1017 ! Nutation in longitude
1018 nutationlongitude = sum(delta_longitude)/36000000.0
1019
1020 ! Nutation in obliquity
1021 nutationobliquity = sum(delta_obliquity)/36000000.0
1022

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ observer_local_hour_calculation()

subroutine narp_module::observer_local_hour_calculation ( real(kind(1d0)), intent(in) apparent_stime_at_greenwich,
real(kind(1d0)), intent(in) locationlongitude,
real(kind(1d0)), intent(in) sun_rigth_ascension,
real(kind(1d0)), intent(out) observer_local_hour )

Definition at line 1144 of file suews_phys_narp.f95.

1146 IMPLICIT NONE
1147
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
1152
1153 observer_local_hour = apparent_stime_at_greenwich + locationlongitude - sun_rigth_ascension
1154 ! Set the range to [0-360]
1155 observer_local_hour = set_to_range(observer_local_hour)

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ prata_emis()

real(kind(1d0)) function narp_module::prata_emis ( real(kind(1d0)) temp_k,
real(kind(1d0)) ea_hpa )

Definition at line 1315 of file suews_phys_narp.f95.

1316 ! clear sky emissivity function Prata 1996
1317 REAL(KIND(1D0)) :: Temp_K, ea_hPa, EMIS_A
1318 REAL(KIND(1D0)) :: W
1319
1320 w = 46.5*(ea_hpa/temp_k)
1321 emis_a = 1.-(1.+w)*exp(-sqrt(1.2 + 3.*w))

Referenced by narp().

Here is the caller graph for this function:

◆ radmethod()

subroutine narp_module::radmethod ( integer, intent(in) netradiationmethod,
integer, intent(in) snowuse,
integer, intent(out) netradiationmethod_use,
integer, intent(out) albedochoice,
integer, intent(out) ldown_option )

Definition at line 44 of file suews_phys_narp.f95.

48 IMPLICIT NONE
49 INTEGER, INTENT(in) :: NetRadiationMethod ! the one from RunControl setting
50 INTEGER, INTENT(in) :: SnowUse
51 INTEGER, INTENT(out) :: NetRadiationMethod_use ! processed NetRadiationMethod to be used for other radiation calculations
52 INTEGER, INTENT(out) :: AlbedoChoice, ldown_option
53 !Determine what should be done with respect to radiation
54 albedochoice = 0
55 ldown_option = 0
56 IF (netradiationmethod == 0) THEN !Observed Q* from the met input file will be used
57 netradiationmethod_use = 0
58 ! ldown_option is not required if NetRadiationMethodX=0 as LDOWN calculations are skipped
59
60 IF (snowuse == 1) THEN !If snow is modelled, NARP is needed for surface temperature
61 ! NetRadiationMethod=3000
62 netradiationmethod_use = 3000
63 ldown_option = 3 !LDOWN will be modelled
64 !NetRadiationMethod=NetRadiationMethod/1000
65 END IF
66
67 ELSEIF (netradiationmethod > 0) THEN !Modelled Q* is used (NARP)
68 albedochoice = -9
69 IF (netradiationmethod < 100) THEN
70 albedochoice = 0
71 ! after the introduction of iteration-based tsurf, TS 20 Sep 2019
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
76 ! recover values before modulus calculation
77 netradiationmethod_use = netradiationmethod
78
79 ! prior to introduction of iteration-based tsurf, TS 20 Sep 2019
80 ! IF (NetRadiationMethod == 1) ldown_option = 1
81 ! IF (NetRadiationMethod == 2) ldown_option = 2
82 ! IF (NetRadiationMethod == 3) ldown_option = 3
83 ! NetRadiationMethod_use = NetRadiationMethod
84
85 ELSEIF (netradiationmethod >= 100 .AND. netradiationmethod < 1000) THEN
86 albedochoice = 1
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
91
92 !choose Ldown method for Spartacus
93 ELSEIF (netradiationmethod > 1000) THEN
94 albedochoice = 0
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
99 ! recover values before modulus calculation
100 netradiationmethod_use = netradiationmethod
101 !If bad NetRadiationMethod value
102 IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9) THEN
103 WRITE (*, *) 'NetRadiationMethod=', netradiationmethod_use
104 WRITE (*, *) 'Value not usable'
105 stop
106 END IF
107
108 END IF
109
110 !If bad NetRadiationMethod value
111 IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9) THEN
112 WRITE (*, *) 'NetRadiationMethod=', netradiationmethod_use
113 WRITE (*, *) 'Value not usable'
114 stop
115 END IF
116 END IF
117

Referenced by suews_driver::suews_cal_qn(), and suews_driver::suews_cal_qn_dts().

Here is the caller graph for this function:

◆ set_to_range()

real(kind(1d0)) function narp_module::set_to_range ( real(kind(1d0)) var)

Definition at line 1281 of file suews_phys_narp.f95.

1282 ! This function make sure the variable is in the specified range.
1283
1284 REAL(KIND(1D0)) :: max_interval
1285 REAL(KIND(1D0)) :: min_interval
1286 REAL(KIND(1D0)) :: var
1287 REAL(KIND(1D0)) :: vari
1288 !
1289 max_interval = 360.0
1290 min_interval = 0.0
1291
1292 vari = var - max_interval*floor(var/max_interval)
1293
1294 IF (vari < min_interval) THEN
1295 vari = vari + max_interval
1296 END IF
1297

Referenced by apparent_stime_at_greenwich_calculation(), earth_heliocentric_position_calculation(), observer_local_hour_calculation(), sun_geocentric_position_calculation(), sun_rigth_ascension_calculation(), and sun_topocentric_zenith_angle_calculate().

Here is the caller graph for this function:

◆ smithlambda()

real(kind(1d0)) function, dimension(365) narp_module::smithlambda ( integer lat)

Definition at line 1427 of file suews_phys_narp.f95.

1428 USE filename
1429 USE defaultnotused
1430 !read kriged data based on Smith 1966 (JAM)
1431 ! Smith, William L.
1432 ! "Note on the relationship between total precipitable water and surface dew point."
1433 ! Journal of Applied Meteorology 5.5 (1966): 726-727.
1434 INTEGER :: lat, ios, ilat
1435 REAL(KIND(1D0)), DIMENSION(365) :: G
1436
1437 !open(99,file="Smith1966.grd",access="direct",action="read",recl=365*4,iostat=ios)
1438 !read(99,rec=lat+1,iostat=ios) G
1439 OPEN (99, file=smithfile, iostat=ios)
1440 DO ilat = 1, lat
1441 READ (99, *)
1442 END DO
1443 READ (99, *, iostat=ios) ilat, g
1444 IF (ios /= 0) THEN
1445 CALL errorhint(11, 'reading Smith1966.grd (ios).', notused, notused, ios)
1446 END IF
1447 CLOSE (99)
defaultnotused::notused
real(kind(1d0)) notused
Definition suews_ctrl_const.f95:1394
filename::smithfile
character(len=90) smithfile
Definition suews_ctrl_const.f95:1694
errorhint
subroutine errorhint(errh, problemfile, value, value2, valuei)
Definition suews_ctrl_error.f95:2

References errorhint(), defaultnotused::notused, and filename::smithfile.

Here is the call graph for this function:

◆ solar_esdist()

real(kind(1d0)) function narp_module::solar_esdist ( integer doy)

Definition at line 1411 of file suews_phys_narp.f95.

1412 !from Stull, 1998 Keep! called from SOLWEIG_clearnessindex_2013b
1413 INTEGER :: doy
1414 REAL(KIND(1D0)) :: Rse
1415 REAL(KIND(1D0)) :: MA, nu, e, a
1416
1417 e = 0.0167
1418 a = 146.457
1419
1420 ma = 2.*3.141592654*(doy - 3)/365.25463 !Mean anomaly
1421 nu = ma + 0.0333988*sin(ma) + .0003486*sin(2.*ma) + 5e-6*sin(3.*ma) !true anomaly
1422 rse = a*(1 - e*e)/(1 + e*cos(nu))
1423

Referenced by isurface().

Here is the caller graph for this function:

◆ sun_geocentric_declination_calculation()

subroutine narp_module::sun_geocentric_declination_calculation ( real(kind(1d0)), intent(in) apparent_sun_longitude,
real(kind(1d0)), intent(in) corr_obliquity,
real(kind(1d0)), intent(in) sun_geocentric_positionlatitude,
real(kind(1d0)), intent(out) sun_geocentric_declination )

Definition at line 1127 of file suews_phys_narp.f95.

1129 IMPLICIT NONE
1130
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
1137
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))
1140
1141 sun_geocentric_declination = asin(argument)*180.0/pi

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ sun_geocentric_position_calculation()

subroutine narp_module::sun_geocentric_position_calculation ( real(kind(1d0)), intent(in) earth_heliocentric_positionlongitude,
real(kind(1d0)), intent(in) earth_heliocentric_positionlatitude,
real(kind(1d0)) sun_geocentric_positionlatitude,
real(kind(1d0)) sun_geocentric_positionlongitude )

Definition at line 883 of file suews_phys_narp.f95.

886 IMPLICIT NONE
887
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
892
893 ! This function compute the sun position relative to the earth.
894
895 sun_geocentric_positionlongitude = earth_heliocentric_positionlongitude + 180.0
896 ! Limit the range to [0,360];
897 sun_geocentric_positionlongitude = set_to_range(sun_geocentric_positionlongitude)
898
899 sun_geocentric_positionlatitude = -earth_heliocentric_positionlatitude
900 ! Limit the range to [0,360]
901 sun_geocentric_positionlatitude = set_to_range(sun_geocentric_positionlatitude)

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sun_rigth_ascension_calculation()

subroutine narp_module::sun_rigth_ascension_calculation ( real(kind(1d0)), intent(in) apparent_sun_longitude,
real(kind(1d0)), intent(in) corr_obliquity,
real(kind(1d0)), intent(in) sun_geocentric_positionlatitude,
real(kind(1d0)), intent(out) sun_rigth_ascension )

Definition at line 1104 of file suews_phys_narp.f95.

1106 IMPLICIT NONE
1107
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
1115
1116 ! This function compute the sun rigth ascension.
1117
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)
1121
1122 sun_rigth_ascension = atan2(argument_numerator, argument_denominator)*180.0/pi
1123 ! Limit the range to [0,360];
1124 sun_rigth_ascension = set_to_range(sun_rigth_ascension)

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sun_topocentric_zenith_angle_calculate()

subroutine narp_module::sun_topocentric_zenith_angle_calculate ( real(kind(1d0)), intent(in) locationlatitude,
real(kind(1d0)), intent(in) topocentric_sun_positiondeclination,
real(kind(1d0)), intent(in) topocentric_local_hour,
real(kind(1d0)) sunazimuth,
real(kind(1d0)) sunzenith )

Definition at line 1230 of file suews_phys_narp.f95.

1232 IMPLICIT NONE
1233
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
1246 ! This function compute the sun zenith angle, taking into account the
1247 ! atmospheric refraction. A default temperature of 283K and a
1248 ! default pressure of 1010 mbar are used.
1249
1250 ! Topocentric elevation, without atmospheric refraction
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
1255
1256 ! Atmospheric refraction correction (in degrees)
1257 argument = corr_elevation + (10.3/(corr_elevation + 5.11))
1258 refraction_corr = 1.02/(60*tan(argument*pi/180.0))
1259
1260 ! For exact pressure and temperature correction, use this,
1261 ! with P the pressure in mbar amd T the temperature in Kelvins:
1262 ! refraction_corr = (P/1010) * (283/T) * 1.02 / (60 * tan(argument * pi/180));
1263
1264 ! Apparent elevation
1265 apparent_elevation = corr_elevation + refraction_corr
1266
1267 sunzenith = 90.0 - apparent_elevation
1268
1269 ! Topocentric azimuth angle. The +180 conversion is to pass from astronomer
1270 ! notation (westward from south) to navigation notation (eastward from
1271 ! north);
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
1276 ! Set the range to [0-360]
1277 sunazimuth = set_to_range(sunazimuth)
1278

References set_to_range().

Referenced by narp_cal_sunposition().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ topocentric_local_hour_calculate()

subroutine narp_module::topocentric_local_hour_calculate ( real(kind(1d0)), intent(in) observer_local_hour,
real(kind(1d0)), intent(in) topocentric_sun_positionrigth_ascension_parallax,
real(kind(1d0)), intent(out) topocentric_local_hour )

Definition at line 1217 of file suews_phys_narp.f95.

1219 IMPLICIT NONE
1220
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
1224
1225 ! This function compute the topocentric local jour angle in degrees
1226
1227 topocentric_local_hour = observer_local_hour - topocentric_sun_positionrigth_ascension_parallax

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ topocentric_sun_position_calculate()

subroutine narp_module::topocentric_sun_position_calculate ( real(kind(1d0)) topocentric_sun_positionrigth_ascension,
real(kind(1d0)) topocentric_sun_positionrigth_ascension_parallax,
real(kind(1d0)) topocentric_sun_positiondeclination,
real(kind(1d0)), intent(in) locationaltitude,
real(kind(1d0)), intent(in) locationlatitude,
real(kind(1d0)), intent(in) observer_local_hour,
real(kind(1d0)), intent(in) sun_rigth_ascension,
real(kind(1d0)), intent(in) sun_geocentric_declination,
real(kind(1d0)), intent(in) earth_heliocentric_positionradius )

Definition at line 1158 of file suews_phys_narp.f95.

1162 IMPLICIT NONE
1163
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
1181
1182 ! topocentric_sun_positionrigth_ascension_parallax
1183 ! This function compute the sun position (rigth ascension and declination)
1184 ! with respect to the observer local position at the Earth surface.
1185
1186 ! Equatorial horizontal parallax of the sun in degrees
1187 eq_horizontal_parallax = 8.794/(3600*earth_heliocentric_positionradius)
1188
1189 ! Term u, used in the following calculations (in radians)
1190 u = atan(0.99664719*tan(locationlatitude*pi/180))
1191
1192 ! Term x, used in the following calculations
1193 x = cos(u) + ((locationaltitude/6378140)*cos(locationaltitude*pi/180))
1194
1195 ! Term y, used in the following calculations
1196 y = (0.99664719d+0*sin(u)) + ((locationaltitude/6378140)*sin(locationlatitude*pi/180))
1197
1198 ! Parallax in the sun rigth ascension (in radians)
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)
1203 ! Conversion to degrees.
1204 topocentric_sun_positionrigth_ascension_parallax = sun_rigth_ascension_parallax*180.0/pi
1205
1206 ! Topocentric sun rigth ascension (in degrees)
1207 topocentric_sun_positionrigth_ascension = sun_rigth_ascension + (sun_rigth_ascension_parallax*180.0/pi)
1208
1209 ! Topocentric sun declination (in degrees)
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

Referenced by narp_cal_sunposition().

Here is the caller graph for this function:

◆ transmissivity()

real(kind(1d0)) function narp_module::transmissivity ( real(kind(1d0)) press_hpa,
real(kind(1d0)) temp_c_dew,
real(kind(1d0)) g,
real(kind(1d0)) zenith )

Definition at line 1451 of file suews_phys_narp.f95.

1452 ! bulk atmospheric transmissivity (Crawford and Duchon, 1999)
1453 ! P = pressure (hPa)
1454 ! Td = dewpoint (C)
1455 ! G parameter is empirical value from Smith 1966 (JAM)
1456 ! zenith in radians
1457 ! if zenith > 80 use the value for 80.
1458
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
1463
1464 IF (zenith > 80.*deg2rad) THEN
1465 cosz = cos(80.*deg2rad)
1466 ELSE
1467 cosz = cos(zenith)
1468 END IF
1469
1470 tdf = temp_c_dew*1.8 + 32. !celsius to fahrenheit
1471 ! Transmission coefficients
1472 m = 35*cosz/sqrt(1224.*cosz*cosz + 1) !optical air mass at p=1013 mb
1473 !Rayleigh & permanent gases
1474 trtpg = 1.021 - 0.084*sqrt(m*(0.000949*press_hpa + 0.051)) !first two trans coeff
1475 u = exp(0.113 - log(g + 1) + 0.0393*tdf) !precipitable water
1476 tw = 1 - 0.077*(u*m)**0.3 !vapor transmission coe3ff.
1477 ta = 0.935**m !4th trans coeff
1478 trans = trtpg*tw*ta !bulk atmospheric transmissivity

Referenced by narp().

Here is the caller graph for this function:

◆ wc_fraction()

real(kind(1d0)) function narp_module::wc_fraction ( real(kind(1d0)), intent(in) rh,
real(kind(1d0)), intent(in) temp )

Definition at line 1346 of file suews_phys_narp.f95.

1347 ! Thomas Loridan, King's College London: June 2009
1348 ! Parameterisation of fraction water content using the relative humidity
1349 REAL(KIND(1D0)), INTENT(in) :: RH !Relative Humidity in %
1350 REAL(KIND(1D0)), INTENT(in) :: Temp !Temperature in degre C
1351
1352 REAL(KIND(1D0)) :: FWC !Fraction water content between 0 and 1
1353 REAL(KIND(1D0)) :: A, B !Parameters in the expo
1354
1355 !Parameters
1356 !A=0.078
1357 !B=0.026
1358
1359 a = 0.185
1360 b = 0.00019*temp + 0.015
1361
1362 !FWC parameterization
1363 fwc = a*(exp(b*rh) - 1)
1364 IF (fwc > 1.) fwc = 1.
1365 IF (fwc < 0.) fwc = 0.

Referenced by narp().

Here is the caller graph for this function:
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0