SUBROUTINE zenith(times,dayi,hour)

!  ! $Author: CPTEC $
!  ! $Date: 2000/02/24 13:15:43 $
!  ! $Revision: 1.0 $
!  !
!  ! subprogram:    zenith      compute the solar zenith angle
!  !   prgrmmr: 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
!  ! usage: call zenith from subroutine radtn
!  !   input argument list:
!  !       times  - the forecast time in seconds
!  !   output argument list:
!  !       dayi   - the day of the year
!  !       hour   - the hour of the day


    REAL, PARAMETER :: d00=0.e0,d5=0.5e0
    REAL, PARAMETER :: h15=15.e0,h24=24.e0,h3600=3600.e0
    REAL, PARAMETER :: csid2=6.570982e-2,csid3=1.002738e0
    REAL, PARAMETER :: clon2=9.856474e-1
    REAL, PARAMETER :: pi=3.1415926,pi2=2.*pi,pih=0.5*pi
    REAL, PARAMETER :: deg2rd=1.745329e-2
    REAL, PARAMETER :: obliq=23.440*deg2rd

    INCLUDE "parmeta.f90"
    INCLUDE "parm.tbl.f90"
    INCLUDE "mpp.h"
#include "sp.h"


    INTEGER, PARAMETER :: lda=lm+9,la=13
    INTEGER, PARAMETER :: imjm=im*jm-jm/2,jam=6+2*(jm-10)
    INTEGER, PARAMETER :: lm1=lm-1,lp1=lm+1
    INTEGER, PARAMETER :: lscrch=4*lm+1+la+1,l1=la+1
    INTEGER, PARAMETER :: l2=la+lm+1,l3=la+2*lm+1
    INTEGER, PARAMETER :: l4=la+3*lm+1
    REAL, INTENT(in)                         :: times
    REAL, INTENT(out)                        :: dayi
    REAL, INTENT(out)                        :: hour
    LOGICAL :: leap
    INTEGER :: iyr,indx,i,j,kmnth,knt
    INTEGER :: l,ll1,ll2
    REAL :: csid1,clon1,day,slon,dec,ra,sidtim,hrang,sinalt,hrlcl
    REAL,DIMENSION(50) :: CSIDX,CLONX

!    USE COMM_CTLBLK
!    USE COMM_PHYS
    INCLUDE "COMM_CTLBLK.f90"
    INCLUDE "COMM_PHYS.f90"

    INTEGER,DIMENSION(12) :: month

    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(l),l=ll1,ll2)
        25 FORMAT(5f9.6)
    ENDDO
    ll2=0
    DO n=1,10
        ll1=ll2+1
        ll2=ll1+4
        READ(nsolar,30)(clonx(l),l=ll1,ll2)
        30 FORMAT(5f11.6)
    ENDDO
    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