      SUBROUTINE OCNALB(SALB, CSAL,                                ! OUTPUT 
     1                  WRKA, WRKB,                                ! INPUT
     2                  ZSPD, CLDT, CSZ,                           ! INPUT
     3                  GC,   WCAP, GT,                            ! INPUT
     4                  FSD,  FSF,  CSD,  CSF,                     ! INPUT
     5                  ILG, IL1, IL2, NBS, ISNOALB, KSALB)        ! CONTROL
C
C     * APR 18/2018 - J. COLE    Adjusted call to lookup table to only 
C     *                          use the visible optical thickness (wrka,wrkb)
C     *                          instead of values for each band.
C     * DEC 02/2014 - M.LAZARE.  New version for gcm18+:
C     *                          - Band-specific whitecap albedoes
C     *                            introduced.
C     *                          - Combination of sfcalb and
C     *                            gc=0 section from oifprp9.
C     *                          - Whitecap effect limited
C     *                            to where SICN<SICN_CRT
C     *                            (WCAP passed in and used).
C
C     * Initialize surface albedo parameters
C     * KSALB = 0 ...surface albedo is calculated in from old CLASS_ formulation
C     * KSALB = 1 ...surface albedo is calculated from a lookup table
C     * KSALB = 2 ...surface albedo is calculated from a formula
C
      IMPLICIT NONE
C
C     * OUTPUT DATA.
C
      REAL, INTENT(OUT) ::
     1 SALB   (ILG,NBS),
     2 CSAL   (ILG,NBS)  
C
C     * INPUT DATA.
C
      REAL ::
     1 WRKA   (ILG,NBS), ! OUTPUT FROM SW_CLD_PROPS3
     2 WRKB   (ILG,NBS)  ! OUTPUT FROM SW_CLD_PROPS3

      REAL, INTENT(IN) ::
     1 ZSPD   (ILG),
     2 CLDT   (ILG),
     3 CSZ    (ILG),
     4 GC     (ILG),
     5 WCAP   (ILG),
     6 GT     (ILG)

      INTEGER, INTENT(IN) :: KSALB, IL1, IL2, ILG, NBS, ISNOALB

      INTEGER, PARAMETER :: NSALB = 4
C
C     * LOCAL DATA.
C
      REAL, DIMENSION(ILG,NBS) :: FSD, FSF, CSD, CSF
      REAL, DIMENSION(ILG,NBS) :: WDIRCTA, WDIRCTC, WDIFFUA, WDIFFUC
C
      REAL, DIMENSION(ILG) :: FSD_BB, FSF_BB, CSD_BB, CSF_BB
      REAL, DIMENSION(ILG) :: WDIRCTA_BB, WDIRCTC_BB,
     1                        WDIFFUA_BB, WDIFFUC_BB
      REAL, DIMENSION(ILG) :: FWCAP, CHL, WRKC, WRKD, VMOD
C
      REAL, PARAMETER :: ! BAND-SPECIFIC WHITECAP ALBEDOES
     1 ALWCAPB(NSALB) = (/0.216, 0.134, 0.044, 0.005/)
C
      REAL  ZSPDMAX, TEMP, ALBALL, ALWCAP, X
      REAL  ACOEF, BCOEF, ARG 
C
      REAL WS, FS, F1, F12, F13, F14, F15, F16, FI, FI2, FI3, WO, WTST
      REAL ENH, WO3, WTTT, FCT, WSS
      REAL ALBDIF, ALBDIR, ALBA, ALBC, ALA, ALC

      INTEGER I, IB
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
C     * CALCULATE WHITECAP FRACTION.
C     * OBSERVATIONAL STUDIES SUGGEST THAT A LARGE FETCH OF
C     * OPEN WATER IS REQUIRED FOR WHITECAPS TO DEVELOP.
C     * THIS IS LIKELY NOT POSSIBLE IF THERE IS A SUFFICIENT
C     * AMOUNT OF SEAICE IN THE GRID SQUARE.
C     * THEREFORE WE LIMIT THIS EFFECT TO OCCUR ONLY WHERE
C     * SICN<SICN_CRT (DEFINED BY WCAP=1. DEFINED IN SFCPROC).
C
      DO I = IL1,IL2
        FWCAP(I)=0.
        VMOD(I) =  MAX(MIN(ZSPD(I), ZSPDMAX), 0.)
        IF(GC(I).EQ.0. .AND. WCAP(I).EQ.1.) THEN
          TEMP=GT(I)-273.16
          FWCAP(I)=MIN(ACOEF(TEMP)*VMOD(I)**BCOEF(TEMP),1.) 
        ENDIF
      ENDDO
C
      IF (KSALB.EQ.1) THEN
        !--- Determine open water albedoes from NASA lookup table
        CHL(IL1:IL2)=0.3 ! no variability due to chlorophyl
        DO IB=1,NBS
          !--- Limit optical depths to 25 which is the largest nodal
          !--- value in the lookup table used by albval
          where (wrka(il1:il2,ib)>25.0) wrka(il1:il2,ib)=25.0
          where (wrkb(il1:il2,ib)>25.0) wrkb(il1:il2,ib)=25.0
          IF (ANY(WRKA(IL1:IL2,IB).LT.0.0))  CALL XIT('OCNALB',-1)
          !--- on return from albval wrkc will contain clear sky albedos
          !--- note for now only passing in the band 1 (visible) optical thickness
          call albval2(wrkc,ib,wrka(1,1),csz,vmod,chl,il1,il2,ilg)
          IF (ANY(WRKB(IL1:IL2,IB).LT.0.0))  CALL XIT('OCNALB',-2)
          !--- on return from albval wrkd will contain all sky albedos
          !--- note for now only passing in the visible (550 nm) optical thickness
          call albval2(wrkd,ib,wrkb(1,1),csz,vmod,chl,il1,il2,ilg)
          !--- whitecap albedo
          ALWCAP=ALWCAPB(IB)
          DO I = IL1,IL2
            IF(GC(I).EQ.0.) THEN
              CSAL(I,IB) = (1. - FWCAP(I)) * WRKC(I) + FWCAP(I) * ALWCAP
              CSAL(I,IB) = MIN( MAX(CSAL(I,IB),0.), 1.)
              ALBALL  = CLDT(I) * (WRKD(I) - WRKC(I)) + WRKC(I)
              SALB(I,IB) = (1. - FWCAP(I)) * ALBALL  + FWCAP(I) * ALWCAP
              SALB(I,IB) = MIN( MAX(SALB(I,IB),0.), 1.)
            ENDIF
          ENDDO
        ENDDO
      ELSE IF (KSALB.EQ.2) THEN
        !--- Determine open water albedoes using Hadley center formula
        DO IB=1,NBS
          !--- whitecap albedo
          ALWCAP=ALWCAPB(IB)
          DO I = IL1,IL2
            IF(GC(I).EQ.0.) THEN
              X         = 0.037/(1.1*(CSZ(I)**1.4) + 0.15)
              SALB(I,IB) = (1. - FWCAP(I)) * X + FWCAP(I) * ALWCAP
              SALB(I,IB) = MIN( MAX(SALB(I,IB),0.), 1.)
              CSAL(I,IB) = SALB(I,IB)
            ENDIF
          ENDDO
        ENDDO
      ELSE
        !--- Albedo values are calculated as in old class_  formula
C
C       * COMPUTE FRACTIONS OF DIRECT AND DIFFUSE BEAMS FOR BOTH
C       * EACH BAND AND OVER ALL BANDS (BOTH ALL-SKY AND CLEAR-SKY).
C
        FSD_BB(:) = 0.0
        FSF_BB(:) = 0.0
        CSD_BB(:) = 0.0
        CSF_BB(:) = 0.0

        DO IB = 1, NBS
          DO I = IL1,IL2
            WDIRCTA(I,IB) = FSD(I,IB)/(FSD(I,IB)+FSF(I,IB)+1.E-10)
            WDIRCTC(I,IB) = CSD(I,IB)/(CSD(I,IB)+CSF(I,IB)+1.E-10)
            WDIFFUA(I,IB) = 1.0 - WDIRCTA(I,IB)
            WDIFFUC(I,IB) = 1.0 - WDIRCTC(I,IB)
            FSD_BB(I)     = FSD_BB(I) + FSD(I,IB)
            FSF_BB(I)     = FSF_BB(I) + FSF(I,IB)
            CSD_BB(I)     = CSD_BB(I) + CSD(I,IB)
            CSF_BB(I)     = CSF_BB(I) + CSF(I,IB)

            IF (WDIRCTA(I,IB) .GT. 1.0 .OR. WDIRCTA(I,IB) .LT. 0.0) THEN
               WRITE(6,*) I,IB,WDIRCTA(I,IB),FSD(I,IB),FSF(I,IB)
               CALL XIT('OCNALB',-3)
            ENDIF

            IF (WDIRCTC(I,IB) .GT. 1.0 .OR. WDIRCTC(I,IB) .LT. 0.0) THEN
               WRITE(6,*) I,IB,WDIRCTC(I,IB),CSD(I,IB),CSF(I,IB)
               CALL XIT('OCNALB',-4)
            ENDIF
          ENDDO
        ENDDO
C
        IF (ISNOALB .EQ. 1) THEN ! COMPUTE FROM BAND MEANS
          DO I = IL1,IL2
            WDIRCTA_BB(I) = FSD_BB(I)/(FSD_BB(I)+FSF_BB(I)+1.E-10)
            WDIFFUA_BB(I) = 1.0 - WDIRCTA_BB(I)
            WDIRCTC_BB(I) = CSD_BB(I)/(CSD_BB(I)+CSF_BB(I)+1.E-10)
            WDIFFUC_BB(I) = 1.0 - WDIRCTC_BB(I)
          ENDDO 
        ELSE IF (ISNOALB .EQ. 0) THEN ! HAVE ONLY TOTAL RADIATION.  SHOULD HAVE SAME VALUE IN 
                                      ! EACH BAND SO JUST USE FIRST ONE.
          DO I = IL1,IL2
            WDIRCTA_BB(I) = FSD(I,1)/(FSD(I,1)+FSF(I,1)+1.E-10)
            WDIFFUA_BB(I) = 1.0 - WDIRCTA_BB(I)
            WDIRCTC_BB(I) = CSD(I,1)/(CSD(I,1)+CSF(I,1)+1.E-10)
            WDIFFUC_BB(I) = 1.0 - WDIRCTC_BB(I)
          ENDDO 
        ENDIF
C
        !--- whitecap albedo
        ALWCAP=0.3     ! old constant value.
        DO I = IL1,IL2
          IF(GC(I).EQ.0.) THEN
            IF(CSZ(I).GE.0.001)     THEN
C                                                                            
C             * OVER OPEN WATER, THE ALBEDOES ARE A COMPLICATED FUNCTION OF
C             * COSINE OF SOLAR ZENITH ANGLE AND WIND SPEED (FOR THE DIRECT
C             * BEAM) AND WIND SPEED FOR THE DIFFUSE BEAM, PLUS THE RESULT
C             * IS WEIGHTED BY THE FRACTION AND ALBEDO OF WHITE CAPS WHICH
C             * MAY BE PRESENT. WHEN THE SUN IS ON OR BELOW THE HORIZON,
C             * AN EQUAL WEIGHT OF DIRECT AND DIFFUSE BEAM IS ASSUMED, SINCE
C             * THERE IS NO RADIATION INFORMATION TO DETERMINE THE PROPER
C             * WEIGHTS AND THERE IS LITTLE SOLAR INPUT.
C
C             * NOTE THAT FOR NOW, THERE IS NO WAVELENGTH DEPENDENCE ON THE
C             * DIRECT OR DIFFUSE BEAMS. 
C
              F1   =  1. - CSZ(I)
              F12  =  F1 * F1
              F13  =  F12 * F1
              F14  =  F12 * F12
              F15  =  F14 * F1
              WS=ZSPD(I)
              F16 = F15*F1
              FI=CSZ(I)
              FI2= FI*FI
              FI3= FI2*FI
              FS=1.-FI2
              WO=180.*FI3*FS

              WTST=WS-WO
              IF(WTST.LT.0.) THEN
                ENH=5.4*FI2*FS*FS
                WO3=WO*WO*WO
                WTTT=WS-1.1*WO
                FCT=1.+ENH*WS*WTTT*WTTT/WO3
              ELSE
                FCT=1.
              ENDIF
              WSS=MAX(0.,WTST)
              ALBDIR=FCT*(0.021 + 0.0421*F12 + 0.128*F13 - 0.04*F16 + 
     1               (4./(5.68+WSS)+(0.074*F1)/(1.+3.*WSS))*F16)
              ALBDIF=.067*(0.55*EXP(-(WS/7.)**2)+1.45*EXP(-(WS/40.)**2)
     1              +1.)/3.
              ALBA=ALBDIR*WDIRCTA_BB(I) + ALBDIF*WDIFFUA_BB(I)
              ALBC=ALBDIR*WDIRCTC_BB(I) + ALBDIF*WDIFFUC_BB(I)
C
              ALA=(1.-FWCAP(I))*ALBA+FWCAP(I)*ALWCAP
              ALC=(1.-FWCAP(I))*ALBC+FWCAP(I)*ALWCAP
C
              SALB(I,1)=ALA
              SALB(I,2)=ALA
              SALB(I,3)=ALA
              SALB(I,4)=ALA
              CSAL(I,1)=ALC
              CSAL(I,2)=ALC
              CSAL(I,3)=ALC
              CSAL(I,4)=ALC
            ELSE
              SALB(I,1)=0.066
              SALB(I,2)=0.066
              SALB(I,3)=0.066
              SALB(I,4)=0.066
              CSAL(I,1)=0.066
              CSAL(I,2)=0.066
              CSAL(I,3)=0.066
              CSAL(I,4)=0.066
            ENDIF
          ENDIF
        ENDDO
      ENDIF !KSALB
C
      RETURN
      END
