C **********************************************************************
C *                                                                    *
C *                                SOFTWARE LICENSING                  *
C *                                                                    *
C *                This program is free software; you can redistribute *
C *                it and/or modify it under the terms of the GNU      *
C *                General Public License as published by the Free     *
C *                Software Foundation, either Version 2 of the        *
C *                license, or (at your option) any later version.     *
C *                                                                    *
C *                This program is distributed in the hope that it     *
C *                will be useful, but without any warranty; without   *
C *                even the implied warranty of merchantability or     *
C *                fitness for a particular purpose. See the GNU       *
C *                General Public License for more details.            *
C *                                                                    *
C *                A copy of the GNU General Public License is         *
C *                available at http://www.gnu.org/copyleft/gpl.html   *
C *                or by writing to the Free Software Foundation, Inc.,*
C *                59 Temple Place - Suite 330, Boston, MA 02111, USA. *
C *                                                                    *
C **********************************************************************

      Subroutine N_CLOUD(CUTOFF,ALAT,TIME,SWOBS,RAD,ATC,CLD)

c       A. Plueddemann 6-12-92
c       Modified from R. Weller's original routine CLOUD c. 6-18-87
c
c       This routine accepts:
c         ALAT = latitude of  observations, decimal degrees.
c         TIME = time in decimal yearday (example, day 1.5 is local
c              noon on the 2nd day).
c         SWOBS = observed shortwave radiation, it is assumed that 
c              SWOBS has been corrected for nighttime bias.
c         ATC = atmospheric transmission coeff., it is expected that
c              an appropriate ATC for the site has been empirically
c              derived.
c       and returns:
c         RAD = theoretical estimate of total shortwave radiation
c              reaching the earths surface
c         CLD = cloud cover fraction (0.0 to 1.0) for use in
c              correcting longwave radiation, nighttime values of
c              CLD are set to a flag value assuming that they will
c              later be interpolated.
c
c       routine to abandon attempts to compute cloud cover just prior
c       to sunset and just after dawn to avoid erroneous estmiates 
c       caused by persistent dawn/dusk differences in theoretical vs.
c       observed short wave.  The value of this cutoff must be "tuned"
c       for particular data sets by trial and error.

c       --------------------------------------------------------------

c       solar constant
      SC = 1395
c       conversion of degrees to radians
      DTOR = 6.283185 / 360.

c       compute altitude of the sun (angular elevation above the
c       horizon) from Smithsonian Met. Tables formulation.
c         SINA = sine of suns altitude
c         THTA = declination of sun 
c         HR = hour angle of the sun

      YRPH = (6.283185/365.) * (357.-TIME)
c       yrph uses 357 instead of 355.90826
      THTA = (-23.45*COS(YRPH)) * DTOR
c       thta uses 23.45 instead of 23.26658
      HR = (TIME+0.5) * 6.283185
      SINA = COS(THTA) * COS(ALAT*DTOR) * COS(HR) + SIN(THTA) * 
     *    SIN(ALAT*DTOR)

c       if low sun angle go to nighttime loop 
      If (SINA.GT.CUTOFF) Then
c       conditional cut off before R calculation (after R in old version)

c       estimate cloud cover using scheme described in Smithsonian
c       Meteorological Tables (R.J. List, 1984).  
c         R = radius vector of the earth, equal to the distance from
c              the center of the earth to the center of the sun
c         RSS = direct E solar radiation at the earths surface.
c         RDIFF = diffuse sky radiation.
c         RAD = total clear sky radiation at earths surface, can be
c              compared to SWOBS to determine ATC.
c         CLD = cloud cover fraction (0 to 1) used for correction to
c              longwave radiation (Fung et al, 1984).

        TRAD = 6.283185 * ((TIME-186.)/365.)
c       trad uses 186 instead of 186.72252
        R = 1.0002 + 0.01671 * COS(TRAD)
c       r uses 1.0002 instead of 1.0001531
        RSS = (SC/(R**2.)) * SINA * (ATC**(1./SINA))
        RDIFF = .5 * (SC/(R**2.)) * SINA * (0.91-(ATC**(1./SINA)))
        RAD = RSS + RDIFF
        CLD = (1.-(SWOBS/RAD)) * (10./7.)
c
c       dont allow cloud fraction greater than 1 or less than 0
        If (CLD.GT.1.) CLD = 1.
        If (CLD.LE.0.) CLD = 0.
      Else

c       nighttime loop, set theoretical radiation 
c       to zero.
        RAD = 0.
c redundant, but included to emphasize that the CLD passed INTO the 
c subroutine is passed OUT of the subroutine, i.e. cloud cover is
c PERSISTED at night.
        CLD = CLD
      End If

      Return
      End