FUNCTION FLOAT Perez (h, slope, alpha, theta, Id, Ibn)

/* Calculate the solar anisotropic sky model of Perez  */

LANGUAGE F77

INPUT
      FLOAT  h      /* height above sea level                    */,
             slope  /* surface slope                             */,
             alpha  /* solar elevation (altitude)                */,
             theta  /* angle of incidence                        */,
             Id     /* diffuse irradiance on the horizotal plane */,
             Ibn    /* surface slope                             */;


CODE

      DOUBLE PRECISION FUNCTION PEREZ(h, slope, alpha, theta, Id, Ibn)

      DOUBLE PRECISION F11(8), F12(8), F13(8), F21(8), F22(8), F23(8)
      DATA F11 / -0.008,  0.130,  0.330,  0.568,
     &            0.873,  1.132,  1.060,  0.678 /
      DATA F12 /  0.588,  0.683,  0.487,  0.187,
     &           -0.392, -1.237, -1.600, -0.372 /
      DATA F13 / -0.062, -0.151, -0.221, -0.295,
     &           -0.362, -0.412, -0.359, -0.250 /
      DATA F21 / -0.060, -0.019,  0.055,  0.109,
     &            0.226,  0.288,  0.264,  0.156 /
      DATA F22 /  0.072,  0.066, -0.064, -0.152,
     &           -0.462, -0.823, -1.127, -1.377 /
      DATA F23 / -0.022, -0.029, -0.026, -0.014,
     &            0.001,  0.056,  0.131,  0.251 /

      DOUBLE PRECISION A, B, F1, F2, Epsilon, Delta, m, h, alpha, theta
      DOUBLE PRECISION PI, Id, Ibn, Idt, Ion, slope, zenith

      PI = DACOS(-1.0D0)

      zenith = 90 - alpha

      A = MAX(0D0,DCOS(Theta*PI/180.0))

      B = MAX(0.087D0,DCOS(zenith*PI/180.0))

      IF (Id .LT. 1D-06) THEN
        Epsilon = ((1D-06+Ibn)/1D-06 + 1.041*(zenith*PI/180.0)**3) / 
     &            (1 + 1.041*(zenith*PI/180.0)**3)
      ELSE
        Epsilon = ((Id+Ibn)/Id + 1.041*(Theta*PI/180.0)**3) / 
     &            (1 + 1.041*(zenith*PI/180.0)**3)
      END IF

      IF (alpha .LT. 1D-06) THEN
        m = 0
      ELSE
        m = EXP(h/1000*(-0.1174-0.0017*h/1000)) / 
     &      (SIN(alpha*PI/180.0) + 0.50572*(alpha+6.07995)**(-1.6364))
      END IF


C     Check value below from Duffie & Beckman, 1991
      Ion = 1370

      Delta = m * Id  / Ion

      IF (Epsilon .LT. 1.065) THEN
        F1 = F11(1) + F12(1) * delta + F13(1) * (Zenith*PI/180.0)
        F2 = F21(1) + F22(1) * delta + F23(1) * (Zenith*PI/180.0)
      ELSE IF (Epsilon .LT. 1.230) THEN
        F1 = F11(2) + F12(2) * delta + F13(2) * (Zenith*PI/180.0)
        F2 = F21(2) + F22(2) * delta + F23(2) * (Zenith*PI/180.0)
      ELSE IF (Epsilon .LT. 1.500) THEN
        F1 = F11(3) + F12(3) * delta + F13(3) * (Zenith*PI/180.0)
        F2 = F21(3) + F22(3) * delta + F23(3) * (Zenith*PI/180.0)
      ELSE IF (Epsilon .LT. 1.950) THEN
        F1 = F11(4) + F12(4) * delta + F13(4) * (Zenith*PI/180.0)
        F2 = F21(4) + F22(4) * delta + F23(4) * (Zenith*PI/180.0)
      ELSE IF (Epsilon .LT. 2.800) THEN
        F1 = F11(5) + F12(5) * delta + F13(5) * (Zenith*PI/180.0)
        F2 = F21(5) + F22(5) * delta + F23(5) * (Zenith*PI/180.0)
      ELSE IF (Epsilon .LT. 4.500) THEN
        F1 = F11(6) + F12(6) * delta + F13(6) * (Zenith*PI/180.0)
        F2 = F21(6) + F22(6) * delta + F23(6) * (Theta*PI/180.0)
      ELSE IF (Epsilon .LT. 6.200) THEN
        F1 = F11(7) + F12(7) * delta + F13(7) * (Zenith*PI/180.0)
        F2 = F21(7) + F22(7) * delta + F23(7) * (Zenith*PI/180.0)
      ELSE
        F1 = F11(8) + F12(8) * delta + F13(8) * (Zenith*PI/180.0)
        F2 = F21(8) + F22(8) * delta + F23(8) * (Zenith*PI/180.0)
      END IF


      Idt = Id *(F1*a/b + (1-F1)*(1+COS(Slope*PI/180))/2 +
     &           F2*SIN(Slope*PI/180))

      Perez = MAX(0D0,Idt)

      END


END_CODE


END_FUNCTION