    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
!>              PHYS0
!>              TEMPCOM
!>              TOPO
!>
!> DRIVER     : INIT
!>              INITS
!>              RADTN
!>              RDTEMP
!>
!> CALLS      : GETDATE
!>--------------------------------------------------------------------------------------------------
    USE CTLBLK
    USE F77KINDS
    USE GLB_TABLE
    USE MAPPINGS
    USE MPPCOM
    USE PARMETA
    USE PARM_TBL
    USE PHYS0
    USE TEMPCOM
    USE TOPO
! 
    REAL   (KIND=R4KIND), PARAMETER :: GSTC1  =   24110.54841
    REAL   (KIND=R4KIND), PARAMETER :: GSTC2  = 8640184.812866
    REAL   (KIND=R4KIND), PARAMETER :: GSTC3  =       9.3104E-2
    REAL   (KIND=R4KIND), PARAMETER :: GSTC4  =      -6.2E-6 
    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
    REAL   (KIND=R4KIND), PARAMETER :: ZEROJD = 2451545.0
!
#include "sp.h"
!
    LOGICAL(KIND=L4KIND)                                                                        ::&
    & LEAP
!
    INTEGER(KIND=I4KIND), DIMENSION(12)                                                         ::&
    & MONTH
!
    DATA MONTH /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
!     
    INTEGER(KIND=I4KIND)                                                                        ::&
    & IYR     , IMO     , IDY     , CYR     , CMO     , CDY     , UTC     , CUTC    , I       ,   &
    & NLEAP   , NYEARS
!
    SAVE MONTH
!-------------------------------------------   
! MODIFICA O CALENDARIO DE 360 PARA 365 DIAS
!-------------------------------------------  
    IYR = IDAT(3)
    IMO = IDAT(1)
    IDY = IDAT(2)
!
    UTC = TIMES / 3600
!
    CALL GETDATE(IYR, IMO, IDY, UTC, CYR, CMO, CDY, CUTC)
!
    NLEAP = 0
!
    DO I=IYR, CYR 
        IF (MOD(I,4) == 0) NLEAP = NLEAP + 1
    END DO
!
    NYEARS = CYR - IYR
!
    DAY  = 0.
    LEAP = .FALSE.
!
    IF (MOD(CYR,4) == 0) THEN
        MONTH(2) = 29
        LEAP     = .TRUE.
    END IF
!
    IF (IDAT(1) > 1) THEN
        KMNTH = IDAT(1) - 1
        DO 10 KNT=1,KMNTH
            DAY = DAY + REAL(MONTH(KNT))
     10 CONTINUE
    END IF
!-----------------------------------------------------------------------------------
!C ALCULATE EXACT NUMBER OF DAYS FROM BEGINNING OF YEAR TO FORECAST TIME OF INTEREST 
!-----------------------------------------------------------------------------------
    DAY = DAY + REAL(IDAT(2) - 1) + (REAL(IHRST) + TIMES / 3600.) / 24.
!    
    DAYI  = REAL(INT(DAY) + 1)
    HOUR  = (DAY - DAYI + 1.) * 24.
    YFCTR = 2000. - CYR
!------------------------------------------------------------------------------------
! FIND CELESTIAL LONGITUDE OF THE SUN THEN THE SOLAR DECLINATION AND RIGHT ASCENSION.
!------------------------------------------------------------------------------------
    IDIFYR = IDAT(3) - 2000
!----------------------------------------------------------------------------- 
! FIND JULIAN DATE OF START OF THE RELEVANT YEAR ADDING IN LEAP DAYS AS NEEDED
!----------------------------------------------------------------------------- 
    IF (IDIFYR < 0) THEN
        ADDDAY = REAL(IDIFYR / 4)
    ELSE
        ADDDAY = REAL((IDIFYR + 3) / 4)
    END IF
!
    STARTYR = ZEROJD + IDIFYR * 365. + ADDDAY - 0.5
!---------------------------------------- 
! THE JULIAN DATE OF THE TIME IN QUESTION
!---------------------------------------- 
    DATJUL = STARTYR + DAY
!------------------------------------------------------------------------
! DIFFERENCE OF ACTUAL JULIAN DATE FROM JULIAN DATE AT 00H 1 JANUARY 2000
!------------------------------------------------------------------------
    DIFJD = DATJUL - ZEROJD
!------------------------------------ 
! MEAN GEOMETRIC LONGITUDE OF THE SUN
!------------------------------------
    SLONM = (280.460 + 0.9856474 * DIFJD) * DEG2RD + YFCTR * PI2
!-----------------
! THE MEAN ANOMOLY
!-----------------
    ANOM = (357.528 + 0.9856003 * DIFJD) * DEG2RD
!----------------------------------------
! APPARENT GEOMETRIC LONGITUDE OF THE SUN
!----------------------------------------
    SLON = SLONM + (1.915 * SIN(ANOM) + 0.020 * SIN(2. * ANOM)) * DEG2RD
!
    IF (SLON > PI2) SLON = SLON - PI2
!--------------------------------
! DECLINATION AND RIGHT ASCENSION
!--------------------------------
    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.
!------------------------------------------------------------------ 
    DATJ0  = STARTYR + DAYI - 1.
    TU     = (DATJ0 - 2451545.) / 36525.
    STIM0  = GSTC1  + GSTC2 * TU + GSTC3 * TU ** 2 + GSTC4 * TU ** 3
    SIDTIM = STIM0  / 3600. + YFCTR * 24. + 1.00273791 * HOUR
    SIDTIM = SIDTIM * 15.   * DEG2RD
!
    IF (SIDTIM < 0. ) 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 < 0.) SINALT = 0.
!
        CZEN(I,J) = SINALT
!
100 CONTINUE
!--------------------------------------------------------------------------------------------------
! 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 - ((365. * (NYEARS-1)) + (NLEAP))
        ELSE IF (LEAP.AND.DAYI > 366.) THEN
            DAYI = DAYI-((365. * (NYEARS)) + (NLEAP - 1))
        END IF
    END IF
!
    RETURN
!
    END SUBROUTINE ZENITH
