      SUBROUTINE SFCALB(ALBSUR, CSALG,                             ! OUTPUT 
     1                  SALB, CSAL,                                ! INPUT
     2                  WRKA, WRKB, ZSPD, CLDT, RMUG, GCG,         ! INPUT
     3                  GT, ISUN,                                  ! INPUT
     4                  KSALB, LENGATH, IL1, IL2, ILG, LAY, IB)    ! CONTROL
C
C     * JAN 19/2008 - M.LAZARE.  COMPUTE CLEAR AND ALL-SKY SURFACE ALBEDO.
C
      IMPLICIT NONE
C
      INTEGER NBS, NBL
      PARAMETER (NBS = 4, NBL = 9)
C
C     * OUTPUT DATA.
C
      REAL, INTENT(OUT) ::
     1 ALBSUR (ILG),
     2 CSALG  (ILG)  
C
C     * INPUT DATA.
C
      REAL, INTENT(IN) ::
     1 SALB   (ILG,NBS),
     2 CSAL   (ILG,NBS)

      REAL ::
     1 WRKA   (ILG),  ! OUTPUT FROM SW_CLD_PROPS
     2 WRKB   (ILG)  ! OUTPUT FROM SW_CLD_PROPS

      REAL, INTENT(IN) ::
     1 ZSPD   (ILG),
     2 CLDT   (ILG),
     3 RMUG   (ILG),
     4 GCG    (ILG),
     6 GT     (ILG)

      INTEGER, INTENT(IN) ::
     1 ISUN   (ILG)

      INTEGER, INTENT(IN) :: KSALB, LENGATH, IL1, IL2, ILG, LAY, IB
C
C     * LOCAL DATA.
C
      REAL, DIMENSION(ILG) :: CHL, ZSPDG, WRKC, WRKD
C
      REAL  ZSPDGUST, ZSPDMAX, TEMP, ALBALL, FWCAP, ALWCAP, X
      REAL  ACOEF, BCOEF, ARG 

      INTEGER I, J, K
C
C     * STATEMENT FUNCTIONS FOR WHITECAPS PARAMETERIZATIONS.
C
      ACOEF(ARG)=6.779E-03-1.83E-03*ARG+1.917E-04*ARG**2
     1          -3.778E-06*ARG**3
      BCOEF(ARG)=0.7566+6.096E-02*ARG-6.547E-03*ARG**2+1.276E-04*ARG**3
C======================================================================
      IF (KSALB.EQ.1) THEN
C
C       * WHEN ALBEDO IS TO BE DETERMINED FROM A TABLE SET THE MAX
C       * WIND SPEED TO THE LARGEST WIND SPEED NODAL VALUE IN THAT
C       * TABLE, TO AVOID EXTRAPOLATION.
C
        ZSPDMAX=18.
      ELSE
        ZSPDMAX=30.
      ENDIF
C
      DO I = 1, LENGATH
        J = ISUN(I)
        ZSPDGUST = ZSPD(J)
        ZSPDG(I) =  MAX(MIN(ZSPDGUST, ZSPDMAX), 0.)
      ENDDO
C
      IF (KSALB.EQ.1) THEN
          !--- Determine open water albedoes from NASA lookup table
          CHL(1:LENGATH)=0.3 ! no variability due to chlorophyl
          !--- Limit optical depths to 25 which is the largest nodal
          !--- value in the lookup table used by albval
          where (wrka(1:lengath)>25.0) wrka=25.0
          where (wrkb(1:lengath)>25.0) wrkb=25.0
          IF (ANY(WRKA(1:LENGATH).LT.0.0))  CALL XIT('SFCALB',-1)
          !--- on return from albval wrkc will contain clear sky albedos
          call albval(wrkc,ib,wrka,rmug,zspdg,chl,lengath)
          IF (ANY(WRKB(1:LENGATH).LT.0.0))  CALL XIT('SFCALB',-2)
          !--- on return from albval wrkd will contain all sky albedos
          call albval(wrkd,ib,wrkb,rmug,zspdg,chl,lengath)
          ALWCAP = 0.3
          DO I = 1, LENGATH
            J = ISUN(I)
            IF (ABS(GCG(I)) .LT. 0.5) THEN
              !--- Open water albedo adjusted for white caps
              TEMP=GT(J)-273.16
              FWCAP=MIN(ACOEF(TEMP)*ZSPDG(I)**BCOEF(TEMP),1.) 
              CSALG(I)  = (1. - FWCAP) * WRKC(I) +  FWCAP * ALWCAP
              CSALG(I)  = MIN( MAX(CSALG(I),0.), 1.)
              ALBALL    = CLDT(J) * (WRKD(I) - WRKC(I)) + WRKC(I)
              ALBSUR(I) = (1. - FWCAP) * ALBALL + FWCAP * ALWCAP
              ALBSUR(I) = MIN( MAX(ALBSUR(I),0. ), 1.)
            ELSE
              !--- Land or sea ice values as determined in class_
              CSALG(I)  = CSAL(J,IB)
              ALBSUR(I) = SALB(J,IB)
            ENDIF
          ENDDO
      ELSE IF (KSALB.EQ.2) THEN
          !--- Determine open water albedoes using Hadley center formula
          ALWCAP = 0.3
          DO I = 1, LENGATH
            J = ISUN(I)
            IF (ABS(GCG(I)) .LT. 0.5) THEN
              !--- Open water albedo adjusted for white caps
              TEMP=GT(J)-273.16
              FWCAP=MIN(ACOEF(TEMP)*ZSPDG(I)**BCOEF(TEMP),1.) 
              X         = 0.037/(1.1*(RMUG(I)**1.4) + 0.15)
              ALBSUR(I) = (1. - FWCAP) * X + FWCAP * ALWCAP
              ALBSUR(I) = MIN( MAX(ALBSUR(I),0. ), 1.)
              CSALG(I)  = ALBSUR(I)
            ELSE
              !--- Land or sea ice values as determined in class_
              CSALG(I)  = CSAL(J,IB)
              ALBSUR(I) = SALB(J,IB)
            ENDIF
          ENDDO
      ELSE
          !--- All albedo values are calculated in class_
          DO I = 1, LENGATH
            J = ISUN(I)
             CSALG(I) = CSAL(J,IB)
            ALBSUR(I) = SALB(J,IB)
          ENDDO
      ENDIF
      RETURN
      END
