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 (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 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 (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 988 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

988  IMPLICIT NONE
989 
990  REAL(KIND(1D0)), INTENT(out) :: aberration_correction
991  REAL(KIND(1D0)), INTENT(in) :: earth_heliocentric_positionradius
992 
993  ! This function compute the aberration_correction, as a function of the
994  ! earth-sun distance.
995 
996  aberration_correction = -20.4898/(3600*earth_heliocentric_positionradius)
997 
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 1017 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

1017  IMPLICIT NONE
1018 
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
1028 
1029  ! This function compute the apparent sideral time at Greenwich.
1030 
1031  jd = julianday
1032  jc = juliancentury
1033 
1034  ! Mean sideral time, in degrees
1035  mean_stime = 280.46061837d+0 + (360.98564736629d+0*(jd - 2451545.0d+0)) + (0.000387933d+0*jc**2) - (jc**3/38710000.0d+0)
1036 
1037  ! Limit the range to [0-360];
1038  mean_stime = set_to_range(mean_stime)
1039 
1040  apparent_stime_at_greenwich = mean_stime + (nutationlongitude*cos(corr_obliquity*pi/180))
real(kind(1d0)), parameter pi
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 1002 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

1002  IMPLICIT NONE
1003 
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
1008 
1009  ! This function compute the sun apparent longitude
1010 
1011  apparent_sun_longitude = sun_geocentric_positionlongitude + nutationlongitude + aberration_correction
1012 
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 1278 of file suews_phys_narp.f95.

Referenced by narp().

1278  REAL(KIND(1d0))::kdown, kclear, fcld
1279 
1280  fcld = 1.-kdown/kclear
1281  IF (fcld > 1.) fcld = 1.
1282  IF (fcld < 0.) fcld = 0.
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 965 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

965  IMPLICIT NONE
966 
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
972  REAL(KIND(1D0)) :: u
973 
974  ! This function compute the true obliquity of the ecliptic.
975 
976  p = (/2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 84381.448/)
977  ! mean_obliquity = polyval(p, julian.ephemeris_millenium/10);
978 
979  u = julianephemeris_millenium/10
980  mean_obliquity =&
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)
983 
984  corr_obliquity = (mean_obliquity/3600) + nutationobliquity
Here is the caller graph for this function:

◆ dewpoint()

real(kind(1d0)) function narp_module::dewpoint ( real(kind(1d0))  Temp_C,
real(kind(1d0))  rh 
)

Definition at line 1241 of file suews_phys_narp.f95.

Referenced by narp().

1241  ! ea = vapor pressure (hPa)
1242  ! td = dewpoint (oC)
1243  ! calculates dewpoint in degC from
1244  ! http://www.atd.ucar.edu/weather_fl/dewpoint.html
1245  ! dewpoint = (237.3 * ln(e_vp/6.1078)) / (17.27 - (ln(e_vp/6.1078)))
1246 
1247  REAL(KIND(1d0))::rh, td, temp_c, g
1248  !http://en.wikipedia.org/wiki/Dew_point
1249  g = ((17.27*temp_c)/(237.7 + temp_c)) + log(rh/100)
1250  td = (237.7*g)/(17.27 - g)
1251  !td = (237.3 * LOG(ea_hPa/6.1078)) / (17.27 - (LOG(ea_hPa/6.1078)))
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 642 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

642  IMPLICIT NONE
643 
644  REAL(KIND(1D0)) :: julianephemeris_millenium
645 
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
690  !REAL(KIND(1D0)) :: L0_terms !>
691  REAL(KIND(1D0)) :: l1
692  !REAL(KIND(1D0)) :: L1_terms !>
693  REAL(KIND(1D0)) :: l2
694  !REAL(KIND(1D0)) :: L2_terms !>
695  REAL(KIND(1D0)) :: l3
696  !REAL(KIND(1D0)) :: L3_terms !>
697  REAL(KIND(1D0)) :: l4
698  !REAL(KIND(1D0)) :: L4_terms !>
699  REAL(KIND(1D0)) :: l5
700  !REAL(KIND(1D0)), dimension(3) :: L5_terms !>
701  !REAL(KIND(1D0)) :: R0_terms !>
702  !REAL(KIND(1D0)) :: R1_terms !>
703  !REAL(KIND(1D0)) :: R2_terms !>
704  !REAL(KIND(1D0)) :: R3_terms !>
705  !REAL(KIND(1D0)), dimension(3) :: R4_terms !>
706  REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793d+0
707 
708  ! This function compute the earth position relative to the sun, using
709  ! tabulated values.
710 
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,&
732  &5.3, 2.65, 4.67/)
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/)
745  a4 = (/114, 8, 1/)
746  b4 = (/3.1420, 4.1300, 3.8400/)
747  c4 = (/0., 6283.08, 12566.15/)
748  a5 = 1.
749  b5 = 3.1400
750  c5 = 0.
751 
752  jme = julianephemeris_millenium
753 
754  ! Compute the Earth Heliochentric longitude from the tabulated values.
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))
761 
762  earth_heliocentric_positionlongitude = &
763  &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4) + (l5*jme**5))/1e8
764  ! Convert the longitude to degrees.
765  earth_heliocentric_positionlongitude = earth_heliocentric_positionlongitude*180./pi
766  ! Limit the range to [0,360[;
767  earth_heliocentric_positionlongitude = set_to_range(earth_heliocentric_positionlongitude)
768 
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/)
772  a1i = (/9, 6/)
773  b1i = (/3.90000000000000, 1.73000000000000/)
774  c1i = (/5507.55000000000, 5223.69000000000/)
775 
776  l0 = sum(a0i*cos(b0i + (c0i*jme)))
777  l1 = sum(a1i*cos(b1i + (c1i*jme)))
778 
779  earth_heliocentric_positionlatitude = (l0 + (l1*jme))/1e8
780  ! Convert the latitude to degrees.
781  earth_heliocentric_positionlatitude = earth_heliocentric_positionlatitude*180/pi
782  ! Limit the range to [0,360];
783  earth_heliocentric_positionlatitude = set_to_range(earth_heliocentric_positionlatitude)
784 
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./)
802  a3j = (/145, 7/)
803  b3j = (/4.273, 3.92/)
804  c3j = (/6283.07600, 12566.1500/)
805  a4j = 4
806  b4j = 2.56
807  c4j = 6283.08000
808 
809  ! Compute the Earth heliocentric radius vector
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))
815 
816  ! Units are in AU
817  earth_heliocentric_positionradius = &
818  &(l0 + (l1*jme) + (l2*jme**2) + (l3*jme**3) + (l4*jme**4))/1e8
819 
real(kind(1d0)), parameter pi
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 1264 of file suews_phys_narp.f95.

Referenced by narp().

1264  !calculates adjusted emissivity due to clouds
1265  REAL(KIND(1d0))::emis_a, fcld, em_adj
1266  !T. Loridan, removing the square for FCLD in the emissivity correction
1267  !em_adj=EMIS_A+(1.-EMIS_A)*FCLD*FCLD
1268  em_adj = emis_a + (1.-emis_a)*fcld
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 1272 of file suews_phys_narp.f95.

Referenced by narp().

1272  !calculates adjusted emissivity due to clouds
1273  REAL(KIND(1d0))::emis_a, fcld, em_adj
1274  em_adj = emis_a + (1.-emis_a)*fcld*fcld
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 1329 of file suews_phys_narp.f95.

References solar_esdist().

Referenced by narp().

1329  ! Calculates ground level solar irradiance clear sky
1330  ! assuming transmissivity = 1
1331  ! let it report zero if zenith >= 90
1332  REAL(KIND(1d0))::zenith, isurf
1333  INTEGER::doy
1334  REAL(KIND(1d0))::rmean, rse, cosz, itoa
1335  REAL(KIND(1D0)), PARAMETER :: deg2rad = 0.017453292
1336 
1337  rmean = 149.6 !Stull 1998
1338  rse = solar_esdist(doy)
1339  IF (zenith < 90.*deg2rad) THEN
1340  cosz = cos(zenith)
1341  itoa = 1370.*(rmean/rse)**2 !top of the atmosphere
1342  isurf = itoa*cosz !ground level solar irradiance in W/m2
1343  ELSE
1344  isurf = 0.
1345  ENDIF
1346 
real(kind(1d0)) function solar_esdist(doy)
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 572 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

572  IMPLICIT NONE
573 
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
582  !REAL(KIND(1D0)) :: time !>
583  REAL(KIND(1D0)) :: ut_time, y !tt,
584  !
585  ! This function compute the julian day and julian century from the local
586  ! time and timezone information. Ephemeris are calculated with a delta_t=0
587  ! seconds.
588 
589  IF (month == 1 .OR. month == 2) THEN
590  y = year - 1.
591  m = month + 12
592  ELSE
593  y = year
594  m = month
595  END IF
596  ut_time = ((float(hour) - utc)/24.) + (float(min)/(60.*24.)) + (sec/(60.*60.*24.)) ! time of day in UT time.
597  d = day + ut_time ! Day of month in decimal time, ex. 2sd day of month at 12:30:30UT, D=2.521180556
598 
599  ! In 1582, the gregorian calendar was adopted
600  IF (year == 1582.) THEN
601  IF (month == 10) THEN
602  IF (day <= 4) THEN ! The Julian calendar ended on October 4, 1582
603  b = 0
604  ELSE IF (day >= 15) THEN ! The Gregorian calendar started on October 15, 1582
605  a = floor(y/100)
606  b = 2 - a + floor(a/4)
607  ELSE
608  !disp('This date never existed!. Date automatically set to October 4, 1582')
609  month = 10
610  day = 4
611  b = 0
612  END IF
613  ELSE IF (month < 10) THEN ! Julian calendar
614  b = 0
615  ELSE ! Gregorian calendar
616  a = floor(y/100)
617  b = 2 - a + floor(a/4)
618  END IF
619 
620  ELSE IF (year < 1582.) THEN ! Julian calendar
621  b = 0
622  ELSE
623  a = floor(y/100) ! Gregorian calendar
624  b = 2 - a + floor(a/4)
625  END IF
626 
627  julianday = floor(365.25*(y + 4716.)) + floor(30.6001*(m + 1)) + d + b - 1524.5
628 
629  delta_t = 0. ! 33.184;
630  julianephemeris_day = julianday + (delta_t/86400)
631 
632  juliancentury = (julianday - 2451545.)/36525.
633 
634  julianephemeris_century = (julianephemeris_day - 2451545.)/36525.
635 
636  julianephemeris_millenium = julianephemeris_century/10.
637 
real(kind(1d0)), dimension(6) y
Here is the caller graph for this function:

◆ narp()

subroutine narp_module::narp ( integer, intent(in)  nsurf,
real(kind(1d0)), dimension(nsurf), intent(in)  sfr,
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)  SnowAlb,
integer, intent(in)  AlbedoChoice,
integer, intent(in)  ldown_option,
integer, intent(in)  NetRadiationMethod_use,
integer, intent(in)  DiagQN,
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_ind,
real(kind(1d0)), intent(out)  alb0,
real(kind(1d0)), intent(out)  alb1 
)

Definition at line 116 of file suews_phys_narp.f95.

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

Referenced by suews_driver::suews_cal_qn().

116  !KCLEAR,FCLD,DTIME,KDOWN,QSTARall,KUPall,LDOWN,LUPall,TSURFall,&
117  !AlbedoChoice,ldown_option,Temp_C,Press_hPa,Ea_hPa,qn1_obs,RH,&
118  !,zenith_degnetRadiationChoice,
119 
120  !SUBROUTINE NARP(QSTAR,DTIME,KDOWN,LDOWN,T,RH,PRESS,FCLD,SNOW)
121  !returns estimate of Q* given the meteorological fields and prior
122  !configuration call.
123 
124  !OUTPUT FIELDS
125  !QSTARall (W/m2) = net all wave radiation
126  !KCLEAR = clear sky incoming solar radiation
127  !KUPall = reflect solar radiation
128  !LDOWN = incoming longwave radiation (observed or modelled depending on ldown_option)
129  !LUPall = outgoing longwaver radiation
130  !FCLD = estimated cloud fraction (used only for emissivity estimate)
131  ! FCLD USED AS INPUT ALSO
132  !TSURFall (DEG C)= estimated surface temperature
133  !QSTAR_SF = net all wave radiation for snow free surface
134  !QSTAR_S = net all wave radiation for SnowPack
135 
136  !INPUT FIELDS
137  !DTIME (days) = local time, not daylight savings
138  !KDOWN (W/m2) = incoming solar radiation
139  !T (DEG C) = air temperature near model height
140  !RH (%) = relative humidity near model height
141  !PRESS (mb) = station pressure, use estimate if unavailable
142  !qn1_obs = Observed Q*
143 
144  !INTERNAL FIELDS
145  ! TemP_K = air temperature in K
146  ! ES_hPs = saturation vapor pressure (hPa)
147  ! EA_Pa = vapor pressure (hPa)
148  ! TD = dewpoint (C)
149  ! ZENITH = solar zenith angle
150  ! ALB0 = surface albedo
151  ! EMIS0 = surface emissivity
152  ! EMIS_A = atmospheric emissivity
153  ! TRANS = atmospheric transmissivity
154  ! LUPCORR = correction for surface heating by kdown (W/m2)
155  ! SIGMATK4 = energy flux density
156  ! KDOWN_HR = hourly average insolation
157  ! DOY = day of year
158 
159  !Modified by LJ to calcuate snow free and SnowPack components (May 2013)
160  !Modified to define variables in data_in module
161  !-------------------------------------------------------------------------------
162  ! USE allocateArray
163  ! use gis_data
164  ! use data_in ! Included 20140701, FL
165  ! use moist ! Included 20140701, FL
166  ! use time ! Included 20140701, FL
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
172  ! REAL(KIND(1D0)),DIMENSION(365),INTENT(in) ::NARP_G
173 
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
185 
186  INTEGER, INTENT(in) ::nsurf
187  INTEGER, INTENT(in) ::netradiationmethod_use ! the one processed by RadMethod
188  INTEGER, INTENT(in) ::albedochoice
189  INTEGER, INTENT(in) ::ldown_option
190  INTEGER, INTENT(in) ::diagqn
191 
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
201 
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
206 
207  REAL(KIND(1d0)), INTENT(out) ::alb0
208  REAL(KIND(1d0)), INTENT(out) ::alb1
209 
210  REAL(KIND(1d0)), DIMENSION(nsurf) ::qn1_ind
211  REAL(KIND(1d0)), DIMENSION(nsurf) ::kup_ind
212  REAL(KIND(1d0)), DIMENSION(nsurf) ::lup_ind
213 
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
219 
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!,RH,DTIME,KDOWN
223  REAL(KIND(1D0)) ::lupcorr, sigmatk4, kdown_hr = 0.
224  INTEGER ::doy, is
225 
226  REAL(KIND(1D0))::qn1_cum, kup_cum, lup_cum, tsurf_cum, & !Cumulative radiation components
227  qn1_is, kup_is, lup_is, tsurf_is, & !Sub-surface radiation components
228  sf_all
229 
230  REAL(KIND(1D0)), PARAMETER :: deg2rad = 0.017453292
231  ! REAL(KIND(1D0)),PARAMETER ::RAD2DEG=57.29577951
232  REAL(KIND(1D0)), PARAMETER :: sigma_sb = 5.67e-8
233 
234  ! NB: NARP_G is not assigned with a value in SUEWS_translate.
235  ! 3.0 is used here as annual average for mid-latitude areas. TS 24 Oct 2017
236  REAL(KIND(1D0)), DIMENSION(365), PARAMETER :: narp_g = 3.0
237  !Initialize variables
238  ! RH=avrh
239  ! DTIME=dectime
240  ! KDOWN=avkdn
241  tsurf_0_k = tsurf_0 + 273.16
242  temp_k = temp_c + 273.16
243  sigmatk4 = sigma_sb*temp_k**4
244  td = dewpoint(temp_c, rh)
245  ! Sun postition is now calculated in the main loop, FL
246  !ZENITH=SOLAR_ZENITH(NARP_LAT,NARP_LONG,NARP_TZ,DTIME)
247  !call NARP_cal_SunPosition(NARP_YEAR,DTIME,NARP_TZ,NARP_LAT,NARP_LONG,Alt,AZIMUTH,ZENITH)
248  zenith = zenith_deg*deg2rad
249  doy = int(dtime)
250  IF (doy == 366) doy = 365
251 
252  !===================================================================================
253  !Calculate radiation for each sub-surface
254  qn1_cum = 0
255  kup_cum = 0
256  lup_cum = 0
257  tsurf_cum = 0
258 
259  qstar_sf = 0
260  qstar_s = 0
261 
262  qn1_ind_snow = 0
263  kup_ind_snow = 0
264  lup_ind_snow = 0
265  tsurf_ind_snow = 0
266 
267  !Total snowfree surface fraction
268  sf_all = 0
269  DO is = 1, nsurf
270  IF (sfr(is) /= 0) sf_all = sf_all + sfr(is)*(1 - snowfrac(is))
271  ENDDO
272 
273  DO is = 1, nsurf
274  IF (diagqn == 1) WRITE (*, *) 'is ', is
275 
276  emis_a = prata_emis(temp_k, press_hpa)
277 
278  !--------------------------------------------------
279  !-------SNOW FREE SURFACE--------------------------
280 
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 !AIDA 1982
283  ELSE
284  alb0 = alb(is)
285  ENDIF
286  emis0 = emis(is)
287 
288  !Downward longwave radiation
289  IF ((ldown_option == 4) .OR. (ldown_option == 5)) THEN !Estimate FCLD from Kdown (Offerle et al. 2003)
290  IF (zenith < 1.5) THEN !DAYTIME CALCULATIONS
291  trans = transmissivity(press_hpa, td, narp_g(doy), zenith)
292  kclear = isurface(doy, zenith)*trans*narp_trans_site
293  IF (kclear > 50.) THEN
294  fcld = cloud_fraction(kdown, kclear)
295  emis_a = emis_cloud_sq(emis_a, fcld)
296  ELSE
297  IF (ldown_option == 5) THEN ! Use RH when Kdown can not be used
298  fcld = wc_fraction(rh, temp_c)
299  emis_a = emis_cloud(emis_a, fcld)
300  ELSE
301  !FCLD is left to the latest calculable value
302  emis_a = emis_cloud_sq(emis_a, fcld)
303  ENDIF
304  ENDIF
305  ELSE !NIGHT TIME CALCULATIONS
306  IF (ldown_option == 4) THEN
307  !FCLD is left to the latest calculable value
308  emis_a = emis_cloud_sq(emis_a, fcld)
309  ELSEIF ((ldown_option == 5)) THEN ! Use RH
310  fcld = wc_fraction(rh, temp_c)
311  emis_a = emis_cloud(emis_a, fcld)
312  ENDIF
313  ENDIF
314  ELSEIF (ldown_option == 3) THEN !Always use RH
315  fcld = wc_fraction(rh, temp_c)
316  emis_a = emis_cloud(emis_a, fcld)
317  ELSEIF (ldown_option == 2) THEN ! FCLD obs, came from input
318  emis_a = emis_cloud(emis_a, fcld)
319  ENDIF
320  IF (diagqn == 1) WRITE (*, *) 'ldown_option: ', ldown_option, 'FCLD:', fcld
321 
322  IF (ldown_option > 1) THEN ! Forcing available if ldown_option=1, model otherwise
323  ldown = emis_a*sigmatk4
324  IF (diagqn == 1) WRITE (*, *) 'EMIS_A: ', emis_a, 'SIGMATK4:', sigmatk4, 'LDOWN: ', ldown
325  ENDIF
326 
327  !----------------------------------------------------------------------------
328  !Note that this is not averaged over the hour for cases where time step < 1hr
329  kdown_hr = kdown
330  IF (kdown_hr > 0) THEN
331  lupcorr = (1 - alb0)*(0.08*kdown_hr)
332  ELSE
333  lupcorr = 0.
334  ENDIF
335 
336  kup = alb0*kdown
337 
338  if (netradiationmethod_use < 10) then
339  ! NARP method
340  tsurf = ((emis0*sigmatk4 + lupcorr)/(emis0*sigma_sb))**0.25 !Eqs. (14) and (15),
341  lup = emis0*sigmatk4 + lupcorr + (1 - emis0)*ldown !Eq (16) in Offerle et al. (2002)
342  else
343  ! use iteration-based approach to calculate LUP and also TSURF; TS 20 Sep 2019
344  tsurf = tsurf_0_k
345  lup = emis0*sigma_sb*tsurf**4 + (1 - emis0)*ldown
346 
347  end if
348 
349  qstar = kdown - kup + ldown - lup
350  tsurf = tsurf - 273.16
351 
352  !======================================================================
353  !Snow related parameters if snow pack existing
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 !AIDA 1982
357  ELSE
358  alb1 = snowalb
359  ENDIF
360 
361  kup_snow = (alb1*(snowfrac(is) - snowfrac(is)*icefrac(is)) + alb0*snowfrac(is)*icefrac(is))*kdown
362 
363  if (netradiationmethod_use < 10) then
364  ! NARP method
365  tsurf_snow = ((narp_emis_snow*sigmatk4)/(narp_emis_snow*sigma_sb))**0.25 !Snow surface temperature
366  !IF (TSURF_SNOW>273.16) TSURF_SNOW=min(273.16,Temp_K)!Set this to 2 degrees (melted water on top)
367  !open(34,file='TestingSnowFrac.txt',position='append')
368  !write(34,*) dectime,is,alb1,alb0,SnowFrac(is),IceFrac(is),KDOWN,KUP_snow
369  !close(34)
370  lup_snow = narp_emis_snow*sigma_sb*tsurf_snow**4 + (1 - narp_emis_snow)*ldown
371  else
372  ! use iteration-based approach to calculate LUP and also TSURF; TS 20 Sep 2019
373  tsurf_snow = tsurf_0_k
374  lup_snow = narp_emis_snow*sigma_sb*tsurf_snow**4 + (1 - narp_emis_snow)*ldown
375  end if
376 
377  qstar_snow = kdown - kup_snow + ldown - lup_snow
378  tsurf_snow = tsurf_snow - 273.16
379 
380  ELSE
381  kup_snow = 0
382  lup_snow = 0
383  tsurf_snow = 0
384  qstar_snow = 0
385  !QSTAR_ICE = 0
386  !KUP_ICE = 0
387  ENDIF
388 
389  qn1_ind_nosnow(is) = qstar !Define sub-surface radiation components
390  kup_ind_nosnow(is) = kup
391  lup_ind_nosnow(is) = lup
392  tsurf_ind_nosnow(is) = tsurf
393 
394  qn1_ind_snow(is) = qstar_snow !Define snow sub-surface radiation components
395  kup_ind_snow(is) = kup_snow
396  lup_ind_snow(is) = lup_snow
397  tsurf_ind_snow(is) = tsurf_snow
398 
399  IF (sf_all /= 0) THEN
400  qstar_sf = qstar_sf + qstar*sfr(is)*(1 - snowfrac(is))/sf_all
401  ELSE
402  qstar_sf = qstar_sf + qstar*sfr(is)*(1 - snowfrac(is))
403  ENDIF
404 
405  IF ((1 - sf_all) /= 0) THEN
406  qstar_s = qstar_s + qstar_snow*sfr(is)*snowfrac(is)/(1 - sf_all)
407  ELSE
408  qstar_s = qstar_s + qstar_snow*sfr(is)*snowfrac(is)
409  ENDIF
410 
411  !---------------------------------------------------------------------
412  !Calculate weighted variables for each subsurface
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)
417 
418  IF (diagqn == 1) WRITE (*, *) 'QSTAR', qstar, 'QSTAR_SNOW', qstar_snow, 'SnowFrac', snowfrac(is)
419 
420  qn1_cum = qn1_cum + (qn1_is*sfr(is)) !Calculate cumulative radiation components
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))
424 
425  qn1_ind(is) = qn1_is !Define sub-surface radiation components
426  kup_ind(is) = kup_is
427  lup_ind(is) = lup_is
428  tsurf_ind(is) = tsurf_is
429 
430  IF (diagqn == 1) WRITE (*, *) 'qn1_is: ', qn1_is
431 
432  ENDDO !End of the surface types
433 
434  !Set overall radiation components
435  IF (netradiationmethod_use /= 3000) THEN !Observed Q* used and snow is modeled
436  qstarall = qn1_cum
437  ELSE
438  qstarall = qn1_obs
439  ENDIF
440 
441  kupall = kup_cum
442  lupall = lup_cum
443  tsurfall = tsurf_cum
444 
445  ! qn1=QSTARall
446  ! kup=KUPall
447  !LDOWN has same name
448  ! LUP=LUPall
449  tsurf = tsurfall
450  !if (kup>500) then
451  ! write(*,*) Kdown, kup, kup_ind(1),kup_ind(2),kup_ind(3),kup_ind(4),kup_ind(5),kup_ind(6),SnowAlb
452  ! pause
453  !endif
454 
455  IF (diagqn == 1) WRITE (*, *) 'kdown: ', kdown, 'kup:', kup, 'LDOWN: ', ldown, 'LUP: ', lup
456  IF (diagqn == 1) WRITE (*, *) 'Qn: ', qstarall
457 
real(kind(1d0)) snowalb
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 464 of file suews_phys_narp.f95.

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 metdisagg::disaggregatemet(), suews_driver::suews_cal_main(), and suews_driver::suews_cal_sunposition().

464  IMPLICIT NONE
465 
466  REAL(KIND(1D0)), INTENT(in) :: year, idectime, utc, locationlatitude, locationlongitude, locationaltitude
467  REAL(KIND(1D0)), INTENT(out) ::sunazimuth, sunzenith
468 
469  REAL(KIND(1D0)):: sec
470  INTEGER :: month, day, hour, min, seas, dayofyear, year_int
471 
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
488  ! REAL(KIND(1D0)) :: sunazimuth,sunzenith
489 
490  ! This function compute the sun position (zenith and azimuth angle (in degrees) at the observer
491  ! location) as a function of the observer local time and position.
492  !
493  ! Input lat and lng should be in degrees, alt in meters.
494  !
495  ! It is an implementation of the algorithm presented by Reda et Andreas in:
496  ! Reda, I., Andreas, A. (2003) Solar position algorithm for solar
497  ! radiation application. National Renewable Energy Laboratory (NREL)
498  ! Technical report NREL/TP-560-34302.
499  ! This document is available at www.osti.gov/bridge
500  ! Code is translated from matlab code by Fredrik Lindberg (fredrikl@gvc.gu.se)
501  ! Last modified: LJ 27 Jan 2016 - Tabs removed
502 
503  ! Convert to timevectors from dectime and year
504  CALL dectime_to_timevec(idectime, hour, min, sec)
505  dayofyear = floor(idectime)
506  year_int = int(year)
507  CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)
508 
509  ! 1. Calculate the Julian Day, and Century. Julian Ephemeris day, century
510  ! and millenium are calculated using a mean delta_t of 33.184 seconds.
511  CALL julian_calculation(year, month, day, hour, min, sec, utc, juliancentury, julianday, julianephemeris_century, &
512  julianephemeris_day, julianephemeris_millenium)
513 
514  ! 2. Calculate the Earth heliocentric longitude, latitude, and radius
515  ! vector (L, B, and R)
516  CALL earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude,&
517  &earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)
518 
519  ! 3. Calculate the geocentric longitude and latitude
520  CALL sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude,&
521  & sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)
522 
523  ! 4. Calculate the nutation in longitude and obliquity (in degrees).
524  CALL nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)
525 
526  ! 5. Calculate the true obliquity of the ecliptic (in degrees).
527  CALL corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)
528 
529  ! 6. Calculate the aberration correction (in degrees)
530  CALL abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)
531 
532  ! 7. Calculate the apparent sun longitude in degrees)
533  CALL apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude,&
534  & aberration_correction, apparent_sun_longitude)
535 
536  ! 8. Calculate the apparent sideral time at Greenwich (in degrees)
537  CALL apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, &
538  &corr_obliquity, apparent_stime_at_greenwich)
539 
540  ! 9. Calculate the sun rigth ascension (in degrees)
541  CALL sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
542  &sun_rigth_ascension)
543 
544  ! 10. Calculate the geocentric sun declination (in degrees). Positive or
545  ! negative if the sun is north or south of the celestial equator.
546  CALL sun_geocentric_declination_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
547  &sun_geocentric_declination)
548 
549  ! 11. Calculate the observer local hour angle (in degrees, westward from south).
550  CALL observer_local_hour_calculation(apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)
551 
552  ! 12. Calculate the topocentric sun position (rigth ascension, declination and
553  ! rigth ascension parallax in degrees)
554  CALL topocentric_sun_position_calculate(topocentric_sun_positionrigth_ascension,&
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)
558 
559  ! 13. Calculate the topocentric local hour angle (in degrees)
560  CALL topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax,&
561  & topocentric_local_hour)
562 
563  ! 14. Calculate the topocentric zenith and azimuth angle (in degrees)
564  CALL sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination,&
565  & topocentric_local_hour, sunazimuth, sunzenith)
566 
subroutine day2month(b, mb, md, seas, year, latitude)
subroutine dectime_to_timevec(dectime, HOURS, MINS, SECS)
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 844 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

844  IMPLICIT NONE
845 
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
864  INTEGER :: i
865  REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793d+0
866 
867  ! This function compute the nutation in longtitude and in obliquity, in
868  ! degrees.
869 
870  ! All Xi are in degrees.
871  jce = julianephemeris_century
872 
873  ! 1. Mean elongation of the moon from the sun
874  p0 = (/(1/189474.), -0.0019142, 445267.11148, 297.85036/)
875  ! X0 = polyval(p, JCE);
876  x0 = p0(1)*jce**3 + p0(2)*jce**2 + p0(3)*jce + p0(4) ! This is faster than polyval...
877 
878  ! 2. Mean anomaly of the sun (earth)
879  p1 = (/-(1/300000.), -0.0001603, 35999.05034, 357.52772/)
880  ! X1 = polyval(p, JCE);
881  x1 = p1(1)*jce**3 + p1(2)*jce**2 + p1(3)*jce + p1(4)
882 
883  ! 3. Mean anomaly of the moon
884  p2 = (/(1/56250.), 0.0086972, 477198.867398, 134.96298/)
885  ! X2 = polyval(p, JCE);
886  x2 = p2(1)*jce**3 + p2(2)*jce**2 + p2(3)*jce + p2(4)
887 
888  ! 4. Moon argument of latitude
889  p3 = (/(1/327270.), -0.0036825, 483202.017538, 93.27191/)
890  ! X3 = polyval(p, JCE);
891  x3 = p3(1)*jce**3 + p3(2)*jce**2 + p3(3)*jce + p3(4)
892 
893  ! 5. Longitude of the ascending node of the moon's mean orbit on the
894  ! ecliptic, measured from the mean equinox of the date
895  p4 = (/(1/450000.), 0.0020708, -1934.136261, 125.04452/)
896  ! X4 = polyval(p, JCE);
897  x4 = p4(1)*jce**3 + p4(2)*jce**2 + p4(3)*jce + p4(4)
898 
899  ! Y tabulated terms from the original code
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, &
903  0, 0, 0, 0, 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, &
906  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., &
923  0.2, -895., 0.5, &
924  1426., -3.4, 54., -0.1, 712., 0.1, -7., 0., -517., 1.2, 224., -0.6, -386., -0.4, 200., 0., -301., &
925  0., 129., -0.1, &
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/))
943  ! Using the tabulated values, compute the delta_longitude and
944  ! delta_obliquity.
945  xi = (/x0, x1, x2, x3, x4/)
946 
947  DO i = 1, 63
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
951  END DO
952 
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)
955 
956  ! Nutation in longitude
957  nutationlongitude = sum(delta_longitude)/36000000.0
958 
959  ! Nutation in obliquity
960  nutationobliquity = sum(delta_obliquity)/36000000.0
961 
real(kind(1d0)), parameter pi
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 1085 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

1085  IMPLICIT NONE
1086 
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
1091 
1092  observer_local_hour = apparent_stime_at_greenwich + locationlongitude - sun_rigth_ascension
1093  ! Set the range to [0-360]
1094  observer_local_hour = set_to_range(observer_local_hour)
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 1255 of file suews_phys_narp.f95.

Referenced by narp().

1255  ! clear sky emissivity function Prata 1996
1256  REAL(KIND(1d0))::temp_k, ea_hpa, emis_a
1257  REAL(KIND(1d0))::w
1258 
1259  w = 46.5*(ea_hpa/temp_k)
1260  emis_a = 1.-(1.+w)*exp(-sqrt(1.2 + 3.*w))
real(kind(1d0)) ea_hpa
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 48 of file suews_phys_narp.f95.

Referenced by suews_driver::suews_cal_qn().

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  ! TODO: this can be wrapped into a subroutine, TS 20 Oct 2017
55  albedochoice = 0
56  ldown_option = 0
57  IF (netradiationmethod == 0) THEN !Observed Q* from the met input file will be used
58  netradiationmethod_use = 0
59  ! ldown_option is not required if NetRadiationMethodX=0 as LDOWN calculations are skipped
60 
61  IF (snowuse == 1) THEN !If snow is modelled, NARP is needed for surface temperature
62  ! NetRadiationMethod=3000
63  netradiationmethod_use = 3000
64  ldown_option = 3 !LDOWN will be modelled
65  !NetRadiationMethod=NetRadiationMethod/1000
66  ENDIF
67 
68  ELSEIF (netradiationmethod > 0) THEN !Modelled Q* is used (NARP)
69  albedochoice = -9
70  IF (netradiationmethod < 100) THEN
71  albedochoice = 0
72  ! after the introduction of iteration-based tsurf, TS 20 Sep 2019
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
77  ! recover values before modulus calculation
78  netradiationmethod_use = netradiationmethod
79 
80  ! prior to introduction of iteration-based tsurf, TS 20 Sep 2019
81  ! IF (NetRadiationMethod == 1) ldown_option = 1
82  ! IF (NetRadiationMethod == 2) ldown_option = 2
83  ! IF (NetRadiationMethod == 3) ldown_option = 3
84  ! NetRadiationMethod_use = NetRadiationMethod
85 
86  ELSEIF (netradiationmethod >= 100 .AND. netradiationmethod < 1000) THEN
87  albedochoice = 1
88  IF (netradiationmethod == 100) ldown_option = 1
89  IF (netradiationmethod == 200) ldown_option = 2
90  IF (netradiationmethod == 300) ldown_option = 3
91  ! NetRadiationMethod=NetRadiationMethod/100
92  netradiationmethod_use = netradiationmethod/100
93  ENDIF
94 
95  !If bad NetRadiationMethod value
96  IF (mod(netradiationmethod, 10) > 3 .OR. albedochoice == -9) THEN
97  WRITE (*, *) 'NetRadiationMethod=', netradiationmethod_use
98  WRITE (*, *) 'Value not usable'
99  stop
100  ENDIF
101  ENDIF
102 
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 1221 of file suews_phys_narp.f95.

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().

1221  ! This function make sure the variable is in the specified range.
1222 
1223  REAL(KIND(1D0)) :: max_interval
1224  REAL(KIND(1D0)) :: min_interval
1225  REAL(KIND(1D0)) :: var
1226  REAL(KIND(1D0)) :: vari
1227  !
1228  max_interval = 360.0
1229  min_interval = 0.0
1230 
1231  vari = var - max_interval*floor(var/max_interval)
1232 
1233  IF (vari < min_interval) THEN
1234  vari = vari + max_interval
1235  END IF
1236 
Here is the caller graph for this function:

◆ smithlambda()

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

Definition at line 1367 of file suews_phys_narp.f95.

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

1367  USE filename
1368  USE defaultnotused
1369  !read kriged data based on Smith 1966 (JAM)
1370  ! Smith, William L.
1371  ! "Note on the relationship between total precipitable water and surface dew point."
1372  ! Journal of Applied Meteorology 5.5 (1966): 726-727.
1373  INTEGER :: lat, ios, ilat
1374  REAL(KIND(1d0)), DIMENSION(365):: g
1375 
1376  !open(99,file="Smith1966.grd",access="direct",action="read",recl=365*4,iostat=ios)
1377  !read(99,rec=lat+1,iostat=ios) G
1378  OPEN (99, file=smithfile, iostat=ios)
1379  DO ilat = 1, lat
1380  READ (99, *)
1381  ENDDO
1382  READ (99, *, iostat=ios) ilat, g
1383  IF (ios /= 0) THEN
1384  CALL errorhint(11, 'reading Smith1966.grd (ios).', notused, notused, ios)
1385  ENDIF
1386  CLOSE (99)
real(kind(1d0)) notused
character(len=90) smithfile
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
Here is the call graph for this function:

◆ solar_esdist()

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

Definition at line 1351 of file suews_phys_narp.f95.

Referenced by isurface().

1351  !from Stull, 1998 Keep! called from SOLWEIG_clearnessindex_2013b
1352  INTEGER ::doy
1353  REAL(KIND(1d0)) ::rse
1354  REAL(KIND(1d0)) ::ma, nu, e, a
1355 
1356  e = 0.0167
1357  a = 146.457
1358 
1359  ma = 2.*3.141592654*(doy - 3)/365.25463 !Mean anomaly
1360  nu = ma + 0.0333988*sin(ma) + .0003486*sin(2.*ma) + 5e-6*sin(3.*ma) !true anomaly
1361  rse = a*(1 - e*e)/(1 + e*cos(nu))
1362 
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 1068 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

1068  IMPLICIT NONE
1069 
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
1076 
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))
1079 
1080  sun_geocentric_declination = asin(argument)*180.0/pi
real(kind(1d0)), parameter pi
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 825 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

825  IMPLICIT NONE
826 
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
831 
832  ! This function compute the sun position relative to the earth.
833 
834  sun_geocentric_positionlongitude = earth_heliocentric_positionlongitude + 180.0
835  ! Limit the range to [0,360];
836  sun_geocentric_positionlongitude = set_to_range(sun_geocentric_positionlongitude)
837 
838  sun_geocentric_positionlatitude = -earth_heliocentric_positionlatitude
839  ! Limit the range to [0,360]
840  sun_geocentric_positionlatitude = set_to_range(sun_geocentric_positionlatitude)
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 1045 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

1045  IMPLICIT NONE
1046 
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
1054 
1055  ! This function compute the sun rigth ascension.
1056 
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)
1060 
1061  sun_rigth_ascension = atan2(argument_numerator, argument_denominator)*180.0/pi
1062  ! Limit the range to [0,360];
1063  sun_rigth_ascension = set_to_range(sun_rigth_ascension)
real(kind(1d0)), parameter pi
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 1171 of file suews_phys_narp.f95.

References set_to_range().

Referenced by narp_cal_sunposition().

1171  IMPLICIT NONE
1172 
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
1185  ! This function compute the sun zenith angle, taking into account the
1186  ! atmospheric refraction. A default temperature of 283K and a
1187  ! default pressure of 1010 mbar are used.
1188 
1189  ! Topocentric elevation, without atmospheric refraction
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
1194 
1195  ! Atmospheric refraction correction (in degrees)
1196  argument = corr_elevation + (10.3/(corr_elevation + 5.11))
1197  refraction_corr = 1.02/(60*tan(argument*pi/180.0))
1198 
1199  ! For exact pressure and temperature correction, use this,
1200  ! with P the pressure in mbar amd T the temperature in Kelvins:
1201  ! refraction_corr = (P/1010) * (283/T) * 1.02 / (60 * tan(argument * pi/180));
1202 
1203  ! Apparent elevation
1204  apparent_elevation = corr_elevation + refraction_corr
1205 
1206  sunzenith = 90.0 - apparent_elevation
1207 
1208  ! Topocentric azimuth angle. The +180 conversion is to pass from astronomer
1209  ! notation (westward from south) to navigation notation (eastward from
1210  ! north);
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
1215  ! Set the range to [0-360]
1216  sunazimuth = set_to_range(sunazimuth)
1217 
real(kind(1d0)), parameter pi
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 1158 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

1158  IMPLICIT NONE
1159 
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
1163 
1164  ! This function compute the topocentric local jour angle in degrees
1165 
1166  topocentric_local_hour = observer_local_hour - topocentric_sun_positionrigth_ascension_parallax
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 1101 of file suews_phys_narp.f95.

Referenced by narp_cal_sunposition().

1101  IMPLICIT NONE
1102 
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
1120 
1121  ! topocentric_sun_positionrigth_ascension_parallax
1122  ! This function compute the sun position (rigth ascension and declination)
1123  ! with respect to the observer local position at the Earth surface.
1124 
1125  ! Equatorial horizontal parallax of the sun in degrees
1126  eq_horizontal_parallax = 8.794/(3600*earth_heliocentric_positionradius)
1127 
1128  ! Term u, used in the following calculations (in radians)
1129  u = atan(0.99664719*tan(locationlatitude*pi/180))
1130 
1131  ! Term x, used in the following calculations
1132  x = cos(u) + ((locationaltitude/6378140)*cos(locationaltitude*pi/180))
1133 
1134  ! Term y, used in the following calculations
1135  y = (0.99664719d+0*sin(u)) + ((locationaltitude/6378140)*sin(locationlatitude*pi/180))
1136 
1137  ! Parallax in the sun rigth ascension (in radians)
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)
1142  ! Conversion to degrees.
1143  topocentric_sun_positionrigth_ascension_parallax = sun_rigth_ascension_parallax*180.0/pi
1144 
1145  ! Topocentric sun rigth ascension (in degrees)
1146  topocentric_sun_positionrigth_ascension = sun_rigth_ascension + (sun_rigth_ascension_parallax*180.0/pi)
1147 
1148  ! Topocentric sun declination (in degrees)
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
real(kind(1d0)), dimension(6) y
real(kind(1d0)), parameter pi
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 1391 of file suews_phys_narp.f95.

Referenced by narp().

1391  ! bulk atmospheric transmissivity (Crawford and Duchon, 1999)
1392  ! P = pressure (hPa)
1393  ! Td = dewpoint (C)
1394  ! G parameter is empirical value from Smith 1966 (JAM)
1395  ! zenith in radians
1396  ! if zenith > 80 use the value for 80.
1397 
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
1402 
1403  IF (zenith > 80.*deg2rad) THEN
1404  cosz = cos(80.*deg2rad)
1405  ELSE
1406  cosz = cos(zenith)
1407  ENDIF
1408 
1409  tdf = temp_c_dew*1.8 + 32. !celsius to fahrenheit
1410  ! Transmission coefficients
1411  m = 35*cosz/sqrt(1224.*cosz*cosz + 1) !optical air mass at p=1013 mb
1412  !Rayleigh & permanent gases
1413  trtpg = 1.021 - 0.084*sqrt(m*(0.000949*press_hpa + 0.051)) !first two trans coeff
1414  u = exp(0.113 - log(g + 1) + 0.0393*tdf) !precipitable water
1415  tw = 1 - 0.077*(u*m)**0.3 !vapor transmission coe3ff.
1416  ta = 0.935**m !4th trans coeff
1417  trans = trtpg*tw*ta !bulk atmospheric transmissivity
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 1286 of file suews_phys_narp.f95.

Referenced by narp().

1286  ! Thomas Loridan, King's College London: June 2009
1287  ! Parameterisation of fraction water content using the relative humidity
1288  REAL(KIND(1d0)), INTENT(in) :: rh !Relative Humidity in %
1289  REAL(KIND(1d0)), INTENT(in) :: temp !Temperature in degre C
1290 
1291  REAL(KIND(1d0)) :: fwc !Fraction water content between 0 and 1
1292  REAL(KIND(1d0)) :: a, b !Parameters in the expo
1293 
1294  !Parameters
1295  !A=0.078
1296  !B=0.026
1297 
1298  a = 0.185
1299  b = 0.00019*temp + 0.015
1300 
1301  !FWC parameterization
1302  fwc = a*(exp(b*rh) - 1)
1303  IF (fwc > 1.) fwc = 1.
1304  IF (fwc < 0.) fwc = 0.
Here is the caller graph for this function: