solar_calculations.F90

Go to the documentation of this file.

module solar_calculations
 
  use iso_c_binding
  use constants_and_conversions
  use meteorological_calculations
  use exceptions
  implicit none
 
  real (c_float) :: earth_sun_dist_dr
  real (c_float) :: solar_declination_delta
 
contains
 
  !> Calculate the number of daylight hours at a location.
  !!
  !! @param[in]  dOmega_s   Sunset hour angle in Radians.
  !! @retval           rN   Number of daylight hours.
  !!
  !! @note Implementation follows equation 34, Allen and others (1998).
  !!
  !! @note Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!       "Crop Evapotranspiration (Guidelines for computing crop water
  !!       requirements)", Food and Agriculture Organization, Rome, Italy.
  function daylight_hours( dOmega_s )    result(dN)   bind(c)
 
    real (c_double), intent(in)   :: domega_s
    real (c_double)               :: dn
 
    dn = 24_c_double / pi * domega_s
 
  end function daylight_hours
 
 
! ***** TO ACCESS THE ABOVE FUNCTION FROM PYTHON VIA CTYPES (once the proper shared library has been created):
! from ctypes import *
! libf = CDLL("libsntemp_lib.dylib")
! libf.daylight_hours.restype = c_double
! libf.daylight_hours(byref(c_double(0.8)))
 
! IF one does not wish to add the 'byref' function call as shown above, the Fortran code itself
! may be modified to force values to be passed by value rather than by reference:
!
!  subroutine daylight_hours_c(dOmega_s, dN)  bind(c, name="daylight_hours_c")
!
!    real (c_double), intent(in), value        :: dOmega_s
!    real (c_double), intent(out)              :: dN
!
!    dN = daylight_hours(dOmega_s)
!
!  end subroutine daylight_hours_c
 
  !------------------------------------------------------------------------------------------------
 
  !>  Calculate extraterrestrial radiation given latitude and time of year.
  !!
  !! @param[in]  dLatitude   Latitude of grid cell in RADIANS.
  !! @param[in]     dDelta   Solar declination in RADIANS.
  !! @param[in]   dOmega_s   Sunset hour angle in RADIANS.
  !! @param[in]     dDsubR   Inverse relative distance Earth-Sun.
  !!
  !! @retval           rRa   Extraterrestrial radiation in MJ / m**2 / day.
  !!
  !! @note  1 MJ = 1e6 Joules; 1 Joule = 1 Watt / sec.
  !! @note   Therefore, multiply by 1e6 and divide by 86400 to get W/m*2-day.
  !!
  !! @note Source
  !!        Equation 21, Allen and others (1998).
  !!
  !! @note Reference
  !!        Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!        "Crop Evapotranspiration (Guidelines for computing crop water
  !!        requirements)", Food and Agriculture Organization, Rome, Italy.
  !!
  !! @sa http://www.fao.org/docrep/x0490e/x0490e07.htm#solar%20radiation
 
  elemental function extraterrestrial_radiation__ra(dLatitude, dDelta, dOmega_s, dDsubR)    result(dRa)
 
    ! [ ARGUMENTS ]
    real (c_double), intent(in) :: dlatitude
    real (c_double), intent(in) :: ddelta
    real (c_double), intent(in) :: domega_s
    real (c_double), intent(in) :: ddsubr
 
    ! [ LOCALS ]
    real (c_double) :: dra
    real (c_double) :: dparta, dpartb
    real (c_double), parameter :: dgsc = 0.0820_c_double  ! MJ / m**2 / min
 
    dparta = domega_s * sin( dlatitude ) * sin( ddelta )
    dpartb = cos( dlatitude ) * cos( ddelta ) * sin( domega_s )
 
    dra = 24_c_double * 60_c_double * dgsc * ddsubr * ( dparta + dpartb ) / pi
 
  end function extraterrestrial_radiation__ra
 
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate net shortwave radiation
  !!
  !! @param[in]      dRs   Incoming shortwave solar radiation, in MJ / m**2 / day
  !! @param[in]  dAlbedo   Albedo or canopy reflection coefficient; 0.23 for grass reference crop
  !! @retval        dRns   Net shortwave radiation, in MJ / m**2 / day
  !!
  !! @note Implementation follows equation 38, Allen and others (1998).
  !!
  !! @note   Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  function net_shortwave_radiation__rns(dRs, dAlbedo)  result(dRns)
 
    real(c_double), intent(in) :: drs
    real(c_double), intent(in) :: dalbedo
 
    ! [ LOCALS ]
    real(c_double) :: drns
 
    drns = (1_c_double - dalbedo) * drs
 
  end function net_shortwave_radiation__rns
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate the solar declination for a given day of the year.
  !!
  !! @param[in]      iDayOfYear   Integer day of the year (January 1 = 1)
  !! @param[in]  iNumDaysInYear   Number of days in the current year
  !! @retval             dDelta   Solar declination, in RADIANS
  !!
  !! @note Implementation follows equation XXX? in:
  !!
  !! @note Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  elemental function solar_declination_simple__delta(iDayOfYear, iNumDaysInYear) result(dDelta)
 
    integer (c_int), intent(in) :: idayofyear
    integer (c_int), intent(in) :: inumdaysinyear
    real (c_double) :: ddelta
 
    ddelta = 0.409_c_double &
             * sin( (twopi * real(idayofyear, c_double) / real(inumdaysinyear, c_double) ) &
                  - 1.39_c_double)
 
  end function solar_declination_simple__delta
 
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate the solar declination for a given day of the year.
  !!
  !! @param[in]      iDayOfYear   Integer day of the year (January 1 = 1)
  !! @param[in]  iNumDaysInYear   Number of days in the current year
  !! @retval             dDelta   Solar declination, in RADIANS
  !!
  !! @note Implementation follows equation 1.3.1 in Iqbal (1983).
  !!
  !! @note Iqbal (1983) reports maximum error of 0.0006 radians; if the last two terms are omitted,
  !!       the reported accuracy drops to 0.0035 radians.
  !!
  !! @note Reference:
  !!       Iqbal, Muhammad (1983-09-28). An Introduction To Solar Radiation (p. 10). Elsevier Science. Kindle Edition.
 
  elemental function solar_declination__delta(iDayOfYear, iNumDaysInYear) result(dDelta)
 
    integer (c_int), intent(in) :: idayofyear
    integer (c_int), intent(in) :: inumdaysinyear
    real (c_double) :: ddelta
 
    ! [ LOCALS ]
    real (c_double) :: dgamma
 
    dgamma = day_angle__gamma( idayofyear, inumdaysinyear )
 
    ddelta =   0.006918_c_double                                         &
             - 0.399912_c_double * cos( dgamma )                         &
             + 0.070257_c_double * sin( dgamma )                         &
             - 0.006758_c_double * cos( 2.0_c_double * dgamma )          &
             + 0.000907_c_double * sin( 2.0_c_double * dgamma )          &
             - 0.002697_c_double * cos( 3.0_c_double * dgamma )          &
             + 0.00148_c_double  * sin( 3.0_c_double * dgamma )
 
 
  end function solar_declination__delta
 
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate the inverse relative Earth-Sun distance for a given day of the year.
  !!
  !! @param[in]      iDayOfYear   Integer day of the year (January 1 = 1)
  !! @param[in]  iNumDaysInYear   Number of days in the current year
  !! @retval             dDsubR   Relative Earth-Sun distance
  !!
  !! @note Implementation follows equation 23, Allen and others (1998): <BR>
  !! @f$ d_r = 1 + 0.033 \cos \left( \frac{ 2 \pi }{365} J \right) @f$  <BR><BR>
  !! where: <BR>
  !!     @f$ d_r @f$ is the inverse relative distance between Earth and the Sun <BR>
  !!     @f$ J @f$ is the current day of the year <BR>
  !!
  !! @note References:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
  !!
  !! @note Duffie, J.A., Beckman, W.A. Solar Engineering of Thermal Processes. New York: Wiley, 1980.
  !!
  !! @note Equation 1.2.3 in Iqbal, Muhammad (1983-09-28). An Introduction To Solar Radiation (p. 28).
  !!       Elsevier Science. Kindle Edition.
 
  elemental function relative_earth_sun_distance__d_r( iDayOfYear, iNumDaysInYear )   result( dDsubR )
 
    ! [ ARGUMENTS ]
    integer (c_int), intent(in) :: idayofyear
    integer (c_int), intent(in) :: inumdaysinyear
    real (c_double) :: ddsubr
 
    ddsubr = 1.0_c_double + 0.033_c_double                                   &
             * cos( twopi * real( idayofyear, c_double )                     &
                                      / real( inumdaysinyear, c_double ) )
 
  end function relative_earth_sun_distance__d_r
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate sunrise/sunset angle, in RADIANS.
  !!
  !! @param[in]  dLatitude   Latitude, in RADIANS
  !! @param[in]     dDelta   Solar declination, in RADIANS
  !! @retval      dOmega_s   Sunset angle, in RADIANS
  !!
  !! @note Implementation follows equation 25, Allen and others (1998).
  !!
  !! @note Hour angle is zero at solar noon. Definition in Iqbal (1983) yields positive values before solar noon,
  !!       and negative values following solar noon.
  !!
  !! @note   Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
   elemental function sunrise_sunset_angle__omega_s( dLatitude, dDelta ) result( dOmega_s )
 
    real (c_double), intent(in) :: dlatitude
    real (c_double), intent(in) :: ddelta
    real (c_double) :: domega_s
 
    domega_s = acos( - tan(dlatitude) * tan(ddelta) )
 
  end function sunrise_sunset_angle__omega_s
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate shortwave solar radiation.
  !!
  !! Calculates the solar radiation using Hargreave's radiation formula.
  !! For use when percent possible daily sunshine value is not available.
  !!
  !! @param[in]    dRa   Extraterrestrial radiation, in MJ / m**2 / day
  !! @param[in]  rTMin   Minimum daily air temperature, in &deg;C
  !! @param[in]  rTMax   Maximum daily air temperature, in &deg;C
  !! @retval       dRa   Solar radiation, in MJ / m**2 / day
  !!
  !! @note Implementation follows equation 50, Allen and others (1998).
  !!
  !! @note  Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  function solar_radiation_hargreaves__rs( dRa, fTMin, fTMax )   result( dRs )   bind(c)
 
    real (c_double), intent(in) :: dra
    real (c_float), intent(in)  :: ftmin
    real (c_float), intent(in)  :: ftmax
    real (c_double)             :: drs
 
    ! [ LOCALS ]
    real (c_double), parameter :: dkrs = 0.175
    real (c_double)            :: tmink
    real (c_double)            :: tmaxk
 
    tmink = c_to_k(real(ftmin, c_double))
    tmaxk = c_to_k(real(ftmax, c_double))
 
    drs = dkrs * sqrt( tmaxk - tmink ) * dra
 
  end function solar_radiation_hargreaves__rs
 
  !------------------------------------------------------------------------------------------------
 
  !> Estimate percent of possible sunshine.
  !!
  !! This function follows from equation 5 in "The Rational Use of the FAO Blaney-Criddle
  !! substituting the rearranged Hargreaves solar radiation formula into
  !! equation 5 results in the formulation below
  !!
  !! @param[in]  fTMax   Maximum daily air temperature, in &deg;K
  !! @param[in]  fTMin   Minimum daily air temperature, in &deg;K
  !! @retval     fPsun   Percentage of possible sunshine, dimensionless percentage
  !!
  !! @todo [Need to add reference here...]
 
  elemental function estimate_percent_of_possible_sunshine__psun(fTMax, fTMin)  result(fPsun)
 
    real (c_float), intent(in) :: ftmax   ! Kelvin degrees
    real (c_float), intent(in) :: ftmin   ! Kelvin degrees
    real (c_float) :: fpsun
 
    ! [ LOCALS ]
    real (c_float), parameter :: fkrs = 0.175
 
    fpsun = ( 2_c_float * fkrs * sqrt( ftmax - ftmin ) ) - 0.5_c_float
 
    if ( fpsun < 0_c_float ) then
      fpsun = 0_c_float
    elseif ( fpsun > 1.0_c_float ) then
      fpsun = 100_c_float
    else
      fpsun = fpsun * 100_c_float
    endif
 
  end function estimate_percent_of_possible_sunshine__psun
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate clear sky solar radiation.
  !!
  !! Calculate the clear sky solar radiation (i.e. when rPctSun = 100,
  !!   n/N=1.  Required for computing net longwave radiation.
  !!
  !! @param[in]  dRa   Extraterrestrial radiation, in MJ / m**2 / day
  !! @param[in]  dAs   Solar radiation regression constant, expressing the fraction
  !!                     of extraterrestrial radiation that reaches earth on OVERCAST days.
  !! @param[in]  sBs   Solar radiation regression constant. As + Bs express the fraction
  !!                     of extraterrestrial radiation that reaches earth on CLEAR days.
  !! @retval    dRso   Clear sky solar radiation, in MJ / m**2 / day
  !!
  !! @note Implementation follows equation 36, Allen and others (1998).
  !!
  !! @note Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  function clear_sky_solar_radiation__rso( dRa, dAs, dBs )   result( dRso )   bind(c)
 
    real (c_double), intent(in)           :: dra
    real (c_double), intent(in), optional :: das
    real (c_double), intent(in),optional  :: dbs
    real (c_double)                       :: drso
 
    ! [ LOCALS ]
    real (c_double) :: das_l
    real (c_double) :: dbs_l
 
    ! assign default value to As if none is provided
    if ( present( das ) ) then
      das_l = das
    else
      das_l = 0.25_c_double
    end if
 
    ! assign default value to Bs if none is provided
    if ( present( dbs ) ) then
      dbs_l = dbs
    else
      dbs_l = 0.5_c_double
    end if
 
    drso = (das_l + dbs_l * dra)
 
  end function clear_sky_solar_radiation__rso
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate the clear sky solar radiation.
  !!
  !! Calculate the clear sky solar radiation (i.e. when rPctSun = 100,
  !!   n/N=1.  Required for computing net longwave radiation.
  !!   For use when no regression coefficients (A, B) are known.
  !!
  !! @param[in]         dRa   Extraterrestrial radiation, in MJ / m**2 / day
  !! @param[in]  fElevation   Elevation, in METERS above sea level
  !! @retval           dRso   Clear sky solar radiation, in MJ / m**2 / day
  !!
  !! @note Implementation follows equation 37, Allen and others (1998).
  !!
  !! @note Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  function clear_sky_solar_radiation_noab__rso(dRa, fElevation) result(dRso)
 
    real (c_double), intent(in) :: dra
    real (c_float), intent(in)  :: felevation
    real (c_double)             :: drso
 
    drso = ( 0.75_c_double + 1.0e-5_c_double * felevation ) * dra
 
  end function clear_sky_solar_radiation_noab__rso
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate solar radiation by means of the Angstrom formula.
  !!
  !! @param[in]      dRa   Extraterrestrial radiation in MJ / m**2 / day
  !! @param[in]      dAs   Solar radiation regression constant, expressing the fraction
  !!                         of extraterrestrial radiation that reaches earth on OVERCAST days.
  !! @param[in]      dBs   Solar radiation regression constant. As + Bs express the fraction
  !!                         of extraterrestrial radiation that reaches earth on CLEAR days.
  !! @param[in]  fPctSun   Percent of TOTAL number of sunshine hours during which the
  !!                         sun actually shown.
  !! @retval         fRs   Solar radiation in MJ / m**2 / day
  !!
  !! @notes
  !!  Implementation follows equation 35, Allen and others (1998).
  !!
  !!   Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  elemental function solar_radiation__rs(dRa, dAs, dBs, fPctSun) result(dRs)
 
    real (c_double), intent(in) :: dra
    real (c_double), intent(in) :: das
    real (c_double), intent(in) :: dbs
    real (c_double), intent(in) :: fpctsun
    real (c_double)             :: drs
 
    drs = ( das + (dbs * fpctsun / 100_c_double ) ) * dra
 
  end function solar_radiation__rs
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate net longwave radiation flux.
  !!
  !! @param[in]  fTMin   Minimum daily air temperature, in &deg;C
  !! @param[in]  fTMax   Maximum daily air temperature, in &deg;C
  !! @param[in]    dRs   Measured or calculated shortwave solar radiation, in MJ / m**2 / day
  !! @param[in]   dRso   Calculated clear-sky radiation, in MJ / m**2 / day
  !! @retval      dRnl   Net longwave solar radiation flux (incoming minus outgoing), in MJ / m**2 / day
  !!
  !! @note
  !! Implementation follows equation 39, Allen and others (1998).
  !!
  !! @note
  !!   Reference:
  !!   Allen, R.G., and others, 1998, FAO Irrigation and Drainage Paper No. 56,
  !!    "Crop Evapotranspiration (Guidelines for computing crop water
  !!    requirements)", Food and Agriculture Organization, Rome, Italy.
 
  elemental function net_longwave_radiation__rnl(fTMin, fTMax, dRs, dRso)  result(dRnl)
 
    real(c_float), intent(in)  :: ftmin
    real(c_float), intent(in)  :: ftmax
    real(c_double), intent(in) :: drs
    real(c_double), intent(in) :: drso
    real(c_double)             :: drnl
 
    ! [ LOCALS ]
    real(c_double)             :: dtavg_k
    real(c_double)             :: dtavg_4
    real (c_double)            :: d_ea
    real (c_double)            :: dcloudfrac
    real (c_double), parameter :: dsigma = 4.903e-9_c_double
 
    dtavg_k = c_to_k((real(ftmin, c_double) + real(ftmax, c_double) ) / 2.0_c_double )
 
    dtavg_4 = dtavg_k * dtavg_k * dtavg_k * dtavg_k * dsigma
    d_ea = dewpoint_vapor_pressure__e_a( ftmin )
 
    dcloudfrac = min( drs / drso, 1.0_c_double )
 
    drnl = dtavg_4 * ( 0.34_c_double - 0.14_c_double * sqrt( d_ea ) ) &
            * ( 1.35_c_double * dcloudfrac - 0.35_c_double )
 
  end function net_longwave_radiation__rnl
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate the day angle in RADIANS.
  !!
  !! This function expresses the integer day number as an angle (in
  !! radians). Output values range from 0 (on January 1st) to
  !! just less than @f$ 2\pi @f$ on December 31st.
  !!
  !! @param[in]        iDayOfYear   Current day of the year.
  !! @param[in]    iNumDaysInYear   Number of days in the current year.
  !! @retval               dGamma   Day angle in RADIANS.
  !!
  !! @note Implementation follows equation 1.2.2 in Iqbal (1983)
  !!
  !! @note Reference:
  !!       Iqbal, Muhammad (1983-09-28). An Introduction To Solar Radiation (p. 3). Elsevier Science. Kindle Edition.
 
  elemental function day_angle__gamma(iDayOfYear, iNumDaysInYear)     result(dGamma)
 
    integer (c_int), intent(in)   :: idayofyear
    integer (c_int), intent(in)   :: inumdaysinyear
    real (c_double)               :: dgamma
 
    dgamma = twopi * ( real(idayofyear, c_double) - 1.0_c_double ) / real(inumdaysinyear, c_double)
 
  end function day_angle__gamma
 
  !------------------------------------------------------------------------------------------------
 
 
  !> Calculate the solar altitude given the zenith angle.
  !!
  !! @param[in]      rTheta_z   Solar zenith angle for given location and time, in RADIANS
  !!
  !! @retval           rAlpha   Solar altitude angle for given location and time, in RADIANS
 
  function solar_altitude__alpha( dTheta_z )    result( dAlpha )   bind(c)
 
    real (c_double), intent(in)   :: dtheta_z
    real (c_double)               :: dalpha
 
    call assert( dtheta_z >= 0.0_c_double .and. dtheta_z <= halfpi, &
      "Internal programming error: solar zenith angle must be in radians and in the range 0 to pi/2", &
      __file__, __line__)
 
    dalpha = halfpi - dalpha
 
  end function solar_altitude__alpha
 
 
  !------------------------------------------------------------------------------------------------
 
    !> Calculate solar zenith angle given latitude, declination, and (optionally) hour angle.
    !!
  !! Calculate solar zenith. Solar zenith angle is the angle between a point directly overhead and
    !! the center of the sun's disk.
    !!
  !! @param[in]     rLatitude   Latitude of location for which estimate is being made, in RADIANS
    !! @param[in]        rDelta   Solar declination angle, in RADIANS
  !! @param[in]        rOmega   [OPTIONAL] Hour angle, measured at the celestial pole between the observer's meridian
  !!                            and the solar meridian, in RADIANS.
  !!
    !! @retval         rTheta_z   Solar zenith angle for given location and time, in RADIANS
    !!
  !! @note Implementation follows equation 9.67, Jacobson, 2005
  !!
  !! @note
  !!   Reference:
  !!   Jacobson, M.Z., 2005, Fundamentals of atmospheric modeling, Second Edition:
  !!   Cambridge University Press.
 
    function zenith_angle__theta_z( dLatitude, dDelta, dOmega ) result( dTheta_z )     bind(c)
 
      real (c_double), intent(in)            :: dlatitude
      real (c_double), intent(in)            :: ddelta
    real (c_double), intent(in), optional  :: domega
 
      ! [ LOCALS ]
      real (c_double) :: dtheta_z
    real (c_double) :: domega_l
 
    if ( present( domega ) ) then
      domega_l = domega_l
    else
      domega_l = 0.0_c_double
    endif
 
      call assert( dlatitude <= halfpi .and. dlatitude >= -halfpi , &
        "Internal programming error: Latitude must be expressed in RADIANS and range from -pi/2 to pi/2", &
        __file__, __line__ )
 
 
      dtheta_z = acos( sin(dlatitude) * sin(ddelta) + cos(dlatitude) * cos(ddelta) * cos(domega_l ) )
 
    end function zenith_angle__theta_z
 
  !------------------------------------------------------------------------------------------------
 
  !> Calculate solar azimuth angle.
  !!
  !! @param[in]     rAlpha   Solar altitude angle, in RADIANS.
  !! @param[in]  rLatitude   Latitude of location for which estimate is being made, in RADIANS.
  !! @param[in]     rDelta   Solar declination angle, in RADIANS.
  !!
  !! @retval          rPsi   Solar azimuth angle, in RADIANS.
  !!
  !! @note   Implementation follows equation 1.5.2a in Iqbal (1983).
  !!
  !! @note Reference:
  !!       Iqbal, Muhammad (1983-09-28). An Introduction To Solar Radiation (p. 15). Elsevier Science. Kindle Edition.
 
  function azimuth_angle__psi(rAlpha, rLatitude, rDelta)   result(rPsi)    bind(c)
 
    real (c_double), intent(in)      :: ralpha
    real (c_double), intent(in)      :: rlatitude
    real (c_double), intent(in)      :: rdelta
    real (c_double)                  :: rpsi
 
    ! [ LOCALS ]
    real (c_double) :: rtempval
 
    rtempval = ( sin( ralpha ) * sin( rlatitude ) - sin( rdelta ) )          &
              /   ( cos( ralpha ) * cos( rlatitude ) )
 
    ! avoid a NaN; cap value at 1.0
    if ( abs( rtempval ) > 1.0_c_double ) rtempval = sign( 1.0_c_double, rtempval )
 
    rpsi = acos( rtempval )
 
  end function azimuth_angle__psi
 
end module solar_calculations