SUBROUTINE zenith(times,dayi,hour)
C
C  ! $Author: CPTEC $
C  ! $Date: 2000/02/24 13:15:43 $
C  ! $Revision: 1.0 $
C  !
C  ! subprogram:    zenith      compute the solar zenith angle
C  !   prgrmmr: black           org: w/nmc2     date: 93-10-28
C  ! abstract:
C  !     zenith calculates the cosine of the solar zenith angles        
C  !     at each point for use in swrad
C  ! usage: call zenith from subroutine radtn
C  !   input argument list:
C  !       times  - the forecast time in seconds
C  !   output argument list:
C  !       dayi   - the day of the year
C  !       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"
       INCLUDE "parm.tbl"
      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(30) :: CSIDX,CLONX

       INCLUDE "CTLBLK.comm"
       INCLUDE "PHYS.comm"

       INTEGER,DIMENSION(12) :: month 

       DATA month/31,28,31,30,31,30,31,31,30,31,30,31/
       SAVE month
C  !*** read annual values of the greenwich mean sidereal time
C  !*** for 0000 utc december 31 and of the geometric mean
C  !*** 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,6
        ll1=ll2+1
        ll2=ll1+4
        READ(nsolar,25)(csidx(l),l=ll1,ll2)
  25    FORMAT(5f9.6)
      ENDDO
      ll2=0
      DO n=1,6
       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

C       !calculate exact number of days from beginning of year to
C       !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

C       !find celestial longitude of the sun then the solar 
C       !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

C       !find the greenwich sidereal time then the local solar
C       !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)
C             !the zenith angle is the complement of the altitude
C             !thus the cosine of the zenith angle equals  
C             !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 CONTINUE


C       !if the forecast is in a different year than the start time,
C       !reset dayi to the proper day of the new year (it must not be
C       !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