Returns the azimuth and zenith angles
Routine modified from NORTRIP
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(in) | :: | lat | |||
real, | intent(in) | :: | lon | |||
integer, | intent(in) | :: | date_a(6) | |||
real(kind=dp), | intent(in) | :: | date_num | |||
real, | intent(in) | :: | difutc_h | |||
real, | intent(out) | :: | azimuth_ang | |||
real, | intent(out) | :: | zenith_ang |
subroutine get_sun_angles(lat, lon, date_a, date_num, difutc_h, azimuth_ang, zenith_ang) !! Returns the azimuth and zenith angles !! !! Routine modified from NORTRIP real, intent(in) :: lat real, intent(in) :: lon integer, intent(in) :: date_a(6) real(dp), intent(in) :: date_num real, intent(in) :: difutc_h real, intent(out) :: azimuth_ang real, intent(out) :: zenith_ang ! Local variables real :: julian_day, time_s, dayang, dec, eqtime, solartime, hourang, azt, az real, parameter :: s0 = 1367.0 real, parameter :: pi = 3.14159/180.0 integer :: ref_year = 2000 if (date_a(1) .eq. 0) then julian_day = real(date_num) else julian_day = date_to_julian(date_a, ref_year) end if time_s = (julian_day - 1)*24.0*3600.0 dayang = 360.0/365.0*(julian_day - 1.0) dec = 0.396 - 22.91*cos(pi*dayang) + 4.025*sin(pi*dayang) eqtime = (1.03 + 25.7*cos(pi*dayang) - 440.0*sin(pi*dayang) - 201.0*cos(2.0*pi*dayang) - 562.0*sin(2.0*pi*dayang))/secphour solartime = mod(time_s + secpday + secphour*(lon/15.0 + difutc_h + eqtime), secpday) hourang = 15.0*(12.0 - solartime/secphour) ! Set zenith angle for atmospheric corrections azt = sin(pi*dec)*sin(pi*lat) + cos(pi*dec)*cos(pi*lat)*cos(pi*hourang) if (abs(azt) .lt. 1.0) then az = acos(azt)/pi else az = 0.0 end if azimuth_ang = 180.0 - hourang zenith_ang = az end subroutine get_sun_angles