    SUBROUTINE ZENITH(TIMES, DAYI, HOUR)
!>--------------------------------------------------------------------------------------------------
!> SUBROUTINE ZENITH
!>
!> SUBPROGRAM: ZENITH - COMPUTE THE SOLAR ZENITH ANGLE
!> PROGRAMMER: BLACK   
!> ORG: W/NMC2 
!> DATE: 93-10-28
!>
!> ABSTRACT:
!> ZENITH CALCULATES THE COSINE OF THE SOLAR ZENITH ANGLES AT EACH POINT FOR USE IN SWRAD
!>
!> PROGRAM HISTORY LOG:
!> 93-10-28  ?????      - ORIGINATOR
!> 18-01-15  LUCCI      - MODERNIZATION OF THE CODE, INCLUDING:
!>                        * F77 TO F90/F95
!>                        * INDENTATION & UNIFORMIZATION CODE
!>                        * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!>                        * DOCUMENTATION WITH DOXYGEN
!>                        * OPENMP FUNCTIONALITY
!>
!> INPUT ARGUMENT LIST:
!> TIMES  - THE FORECAST TIME IN SECONDS
!>
!> OUTPUT ARGUMENT LIST:
!> DAYI   - THE DAY OF THE YEAR
!> HOUR   - THE HOUR OF THE DAY
!>
!> INPUT/OUTPUT ARGUMENT LIST:
!> NONE
!>
!> USE MODULES: CTLBLK
!>              F77KINDS
!>              GLB_TABLE
!>              MAPPINGS
!>              MPPCOM
!>              PARMETA
!>              PARM_TBL
!>              PHYS
!>              TEMPCOM
!>              TOPO
!>
!> DRIVER     : INIT
!>              INITS
!>              RADTN
!>              RDTEMP
!>
!> CALLS      : -----
!>--------------------------------------------------------------------------------------------------
    USE CTLBLK
    USE F77KINDS
    USE GLB_TABLE
    USE MAPPINGS
    USE MPPCOM
    USE PARMETA
    USE PARM_TBL
    USE PHYS
    USE TEMPCOM
    USE TOPO
!
    IMPLICIT NONE
!
    REAL   (KIND=R4KIND), PARAMETER :: D00    =    0.0E0 
    REAL   (KIND=R4KIND), PARAMETER :: D5     =    0.5E0
    REAL   (KIND=R4KIND), PARAMETER :: H15    =   15.E0 
    REAL   (KIND=R4KIND), PARAMETER :: H24    =   24.E0 
    REAL   (KIND=R4KIND), PARAMETER :: H3600  = 3600.E0
    REAL   (KIND=R4KIND), PARAMETER :: CSID2  =    6.570982E-2
    REAL   (KIND=R4KIND), PARAMETER :: CSID3  =    1.002738E0
    REAL   (KIND=R4KIND), PARAMETER :: CLON2  =    9.856474E-1
    REAL   (KIND=R4KIND), PARAMETER :: PI     =    3.1415926 
    REAL   (KIND=R4KIND), PARAMETER :: PI2    =    2.    * PI
    REAL   (KIND=R4KIND), PARAMETER :: PIH    =    0.5   * PI
    REAL   (KIND=R4KIND), PARAMETER :: DEG2RD =    1.745329E-2
    REAL   (KIND=R4KIND), PARAMETER :: OBLIQ  =   23.440 * DEG2RD
!
#include "sp.h"
!
    INTEGER(KIND=I4KIND), PARAMETER :: LDA    = LM +  9 
    INTEGER(KIND=I4KIND), PARAMETER :: LA     =      13
    INTEGER(KIND=I4KIND), PARAMETER :: IMJM   = IM * JM - JM  / 2 
    INTEGER(KIND=I4KIND), PARAMETER :: JAM    = 6  +  2 * (JM - 10)
    INTEGER(KIND=I4KIND), PARAMETER :: LSCRCH = 4  * LM + 1 + LA + 1
    INTEGER(KIND=I4KIND), PARAMETER :: L1     = LA +  1
    INTEGER(KIND=I4KIND), PARAMETER :: L2     = LA + LM + 1 
    INTEGER(KIND=I4KIND), PARAMETER :: L3     = LA +  2 * LM + 1
    INTEGER(KIND=I4KIND), PARAMETER :: L4     = LA +  3 * LM + 1
!
    REAL   (KIND=R4KIND)                                                  , INTENT(IN)          ::&
    & TIMES
!
    REAL   (KIND=R4KIND)                                                  , INTENT(OUT)         ::&
    & DAYI
!
    REAL   (KIND=R4KIND)                                                  , INTENT(OUT)         ::&
    & HOUR
!
    LOGICAL(KIND=L4KIND)                                                                        ::&
    & LEAP
!
    INTEGER(KIND=I4KIND)                                                                        ::&
    & IYR     , INDX    , I       , J       , KMNTH   , KNT
!
    INTEGER(KIND=I4KIND)                                                                        ::&
    & LL1     , LL2
!
    REAL   (KIND=R4KIND)                                                                        ::&
    & CSID1   , CLON1   , DAY     , SLON    , DEC     , RA      , SIDTIM  , HRANG   , SINALT  ,   &
    & HRLCL
!
    REAL   (KIND=R4KIND), DIMENSION(50)                                                         ::&
    & CSIDX   , CLONX
!
    INTEGER(KIND=I4KIND), DIMENSION(12)                                                         ::&
    & MONTH
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
    INTEGER(KIND=I4KIND)                                                                        ::&
    & N       , K       , NSOLAR
!
    DATA MONTH /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30,31/
    SAVE MONTH
!-------------------------------------------------------------------------------------------------- 
! READ ANNUAL VALUES OF THE GREENWICH MEAN SIDEREAL TIME FOR 0000 UTC DECEMBER 31 AND OF THE 
! GEOMETRIC MEAN CELESTIAL LONGITUDE OF THE SUN FOR 0000 TDT DECEMBER 31
!-------------------------------------------------------------------------------------------------- 
    NSOLAR = 19
    OPEN (NSOLAR, FILE='./SOLAR', FORM='FORMATTED', STATUS='OLD')
!
    LL2 = 0
    DO N=1,10
        LL1 = LL2 + 1
        LL2 = LL1 + 4
        READ(NSOLAR,25) (CSIDX(K), K=LL1,LL2)
        25 FORMAT(5F9.6)
    END DO
!
    ll2 = 0
!
    DO N=1,10
        LL1 = LL2 + 1
        LL2 = LL1 + 4
        READ(NSOLAR,30) (CLONX(K), K=LL1,LL2)
        30 FORMAT(5F11.6)
    END DO
!
    CLOSE(NSOLAR)
!
    IYR   = IDAT(3) - 1900
    INDX  = IYR     -   85
    CSID1 = CSIDX(INDX)
    CLON1 = CLONX(INDX)
    DAY   = D00
    LEAP  = .FALSE. 
!
    IF (MOD(IDAT(3),4) == 0) THEN
        MONTH(2) = 29
        LEAP     = .TRUE. 
    END IF
!
    IF (IDAT(1) > 1) THEN
        KMNTH = IDAT(1) - 1
        DO  KNT=1,KMNTH
            DAY = DAY + REAL(MONTH(KNT))
        END DO
    END IF
!-----------------------------------------------------------------------------------
! CALCULATE EXACT NUMBER OF DAYS FROM BEGINNING OF YEAR TO FORECAST TIME OF INTEREST
!-----------------------------------------------------------------------------------
    DAY  = DAY + REAL(IDAT(2) - 1) + (REAL(IHRST)+ TIMES/H3600) / H24
    DAYI = REAL(INT(DAY) + 1)
    HOUR = (DAY - DAYI) * H24 + H24
!------------------------------------------------------------------------------------
! FIND CELESTIAL LONGITUDE OF THE SUN THEN THE SOLAR DECLINATION AND RIGHT ASCENSION.
!------------------------------------------------------------------------------------
    SLON = (CLON1 + CLON2 * DAY) * DEG2RD
!
    IF (SLON > PI2) SLON = SLON - PI2
!
    DEC = ASIN(SIN(SLON) * SIN(OBLIQ))
    RA  = ACOS(COS(SLON) / COS(DEC))
!
    IF (SLON > PI) RA = PI2 - RA
!------------------------------------------------------------------
! FIND THE GREENWICH SIDEREAL TIME THEN THE LOCAL SOLAR HOUR ANGLE.
!------------------------------------------------------------------
    SIDTIM = CSID1  + CSID2 * DAYI + CSID3 * HOUR
    SIDTIM = SIDTIM * H15   * DEG2RD
!
    IF (SIDTIM < D00) SIDTIM = SIDTIM + PI2
    IF (SIDTIM > PI2) SIDTIM = SIDTIM - PI2
!
    HRANG = SIDTIM - RA
!
    DO 100 J=MYJS,MYJE
        DO 100 I=MYIS,MYIE
            HRLCL = HRANG - GLON(I,J)
!--------------------------------------------------------------------------------------------------
! THE ZENITH ANGLE IS THE COMPLEMENT OF THE ALTITUDE THUS THE COSINE OF THE ZENITH ANGLE EQUALS THE
! SINE OF THE ALTITUDE.
!--------------------------------------------------------------------------------------------------
            SINALT = SIN(DEC) * SIN(GLAT(I,J)) + COS(DEC) * COS(HRLCL) * COS(GLAT(I,J))
!
            IF (SINALT < D00) SINALT = D00
!
            CZEN(I,J) = SINALT
!
100 END DO
!-------------------------------------------------------------------------------------------------- 
! IF THE FORECAST IS IN A DIFFERENT YEAR THAN THE START TIME, RESET DAYI TO THE PROPER DAY OF THE 
! NEW YEAR (IT MUST NOT BE RESET BEFORE THE SOLAR ZENITH ANGLE IS COMPUTED).
!-------------------------------------------------------------------------------------------------- 
    IF (DAYI > 365.) THEN
        IF (.NOT. LEAP) THEN
            DAYI = DAYI - 365.
        ELSE IF (LEAP .AND. DAYI > 366.) THEN
            DAYI = DAYI - 366.
        END IF
    END IF
!
    RETURN
!
    END SUBROUTINE ZENITH
