!+ Program to calculate the solar flux at the top of the atmosphere.
!
! Method:
!	A solar constant and a day of the year are supplied. The
!	solar flux at the top of the atmosphere is calculated 
!	accordingly using the same algorithm as in the Unified Model.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date		Comment
!	1.0		12-04-95	First Version under RCS
!					(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77 with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      PROGRAM SOLAR_INC
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/astron_constants.ccf'
      INCLUDE 'incl/math_cnst.ccf'
      INCLUDE 'incl/def_std_io.icf'
!
!     DECLARATION OF VARIABLES.
      INTEGER
     &     IOS
!             I/O ERROR FLAG
      REAL*8
     &     DAY
!             DAY OF YEAR
     &   , SOLAR_CONSTANT
!             USER'S MEAN VALUE OF SOLAR CONSTANT
     &   , SOLAR_CONSTANT_LOCAL
!             LOCAL SOLAR CONSTANT
      REAL*8
     &     M
!             MEAN ANOMALY
     &   , V
!             TRUE ANOMALY
!
!
!     OBTAIN DAY OF YEAR AND SOLAR CONSTANT.
      WRITE(IU_USER, '(/A)')
     &   'ENTER DAY OF YEAR AND MEAN SOLAR CONSTANT.'
1     READ(IU_STDIN, *, IOSTAT=IOS) DAY, SOLAR_CONSTANT
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(A)') '+++ ILLEGAL INPUT: PLEASE RE-ENTER.'
         GOTO 1
      ENDIF
!
      M=2.0D+00*PI*(DAY-DAY_PERIHELION)/LENGTH_YEAR
      V=M+(2.0D+00*ECCENTRICITY-0.25D+00*ECCENTRICITY**3)*SIN(M)
     &   +0.8D+00*ECCENTRICITY**2*SIN(2.0D+00*M)
     &   +(13.0D+00*ECCENTRICITY**3/12.0D+00)*SIN(3.0D+00*M)
      SOLAR_CONSTANT_LOCAL=SOLAR_CONSTANT*(
     &    (1.0D+00+0.5D+00*ECCENTRICITY**2)
     &   *(1.0D+00+ECCENTRICITY*COS(V))/(1.0D+00-ECCENTRICITY**2))**2
!
      WRITE(IU_STDOUT, '(/A, 1PE12.5)') 'LOCAL SOLAR CONSTANT = '
     &   , SOLAR_CONSTANT_LOCAL
!
!
      STOP
      END
