      SUBROUTINE SNOSICA(SALB,CSAL,
     1                   FSD,FSF,CSD,CSF,
     2                   SNO,TS,ALBS,GC,GTA,SIC,CSZ,
     3                   RHON,GT,REFSN,BCSN,
     4                   ILG,IS,IF,NBS,ISNOALB)
C
C     * Aug 06/19 - S.Kharin/ Modify definiton of SNOFAC to handle
C     *             M.Lazare. RHON=0. Cosmetic fix for CMIP6; will
C     *                       be looked at to avoid RHON=0 in next
C     *                       development cycle.
C     * APR 22/18 - M.Lazare. Replace fractional snow coverage
C     *                       by that from CLASSA (ie linear based
C     *                       on depth, instead of sqrt based on SWE).
C     * APR 30/18 - M.Lazare. Remove unused {ACOEF,BCOEF}. Calculation
C     *                       now in OCNALB.
C     * DEC 01/14 - M.Lazare. New version for gcm18+:
C     *                       - Based on old OIFPRP9 but without
C     *                         simpler calculations moved to new
C     *                         OIFPRP10.
C     *                       - Leads fraction effect removed
C     *                         since now applies to ice-covered
C     *                         portion only.
C     *                       - Unused {UL,VL,VMOD} removed.
C     *                       - Initialize C_FLAG to zero (bugfix).
C
C     * PERFORMS ALBEDO CALCULATIONS OVER SEA-ICE.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C 
C     * OUTPUT OR I/O FIELDS:
C
      REAL   SALB (ILG,NBS), CSAL (ILG,NBS)
C
C     * INPUT FIELDS:
C
      REAL   FSD(ILG,NBS), FSF(ILG,NBS), CSD(ILG,NBS), CSF(ILG,NBS)
C
      REAL   SNO    (ILG), TS     (ILG), ALBS   (ILG), GC     (ILG),
     1       GTA    (ILG), SIC    (ILG), CSZ    (ILG),
     2       RHON   (ILG), GT     (ILG), REFSN  (ILG), BCSN   (ILG)
C
C     * LOCAL FIELDS ("A"=ALL-SKY, "C"=CLEAR-SKY):
C
      REAL ALBBI(ILG,NBS), ALSNODIF(ILG,NBS), ALSNODIR(ILG,NBS)
C
      REAL ALSNOA (ILG,NBS), ALSNOC (ILG,NBS)
      REAL WDIRCTA(ILG,NBS), WDIRCTC(ILG,NBS)
      REAL WDIFFUA(ILG,NBS), WDIFFUC(ILG,NBS)
C
      REAL WDIRCTA_BB(ILG), WDIRCTC_BB(ILG)
      REAL WDIFFUA_BB(ILG), WDIFFUC_BB(ILG)
      REAL FSD_BB(ILG), FSF_BB(ILG), CSD_BB(ILG), CSF_BB(ILG)
C
      REAL SNOFAC(ILG), REFSNO(ILG), BCSNO(ILG)
      INTEGER C_FLAG(ILG), SUM_C_FLAG
C
      COMMON /PARAM1/ PI,     RVORD, TFREZ,  HS,   HV,   DAYLNT
      COMMON /PARAM3/ CSNO,   CPACK, GTFSW,  RKHI, SBC,  SNOMAX
      COMMON /PARAM5/ CONI,DENI,SICMIN,DFSIC,CONF
C
C     * MELT POND ALBEDOES FOR 4 SPECTRAL BANDS (BASED ON
C     * EBERT AND CURRY (1993), ASSUMING MELT POND DEPTH OF 0.2M
C
      DATA ALBP1,ALBP2,ALBP3,ALBP4 /0.27, 0.05, 0.05, 0.03/
      DATA PUDDLE_MAX  / 0.25   /
C-----------------------------------------------------------------------
C     * CHECK THAT THE ASSUMED BAND STRUCTURE, CURRENTLY 4 BANDS, 
C     * MATCHES THAT OF THE EXTERNAL RADIATIVE MODEL.
      
      IF (NBS .NE. 4) THEN
         CALL XIT("SNOSICA",-1)
      END IF
 
      SUM_C_FLAG = 0
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 = IS, IF
            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("SNOSICA",-2)
            END IF

            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("SNOSICA",-3)
            END IF
         END DO
      END DO
C
      IF (ISNOALB .EQ. 1) THEN ! COMPUTE FROM BAND MEANS
         DO I = IS, IF
            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)
         END DO 
      ELSE IF (ISNOALB .EQ. 0) THEN ! HAVE ONLY TOTAL RADIATION.  SHOULD HAVE SAME VALUE IN 
                                    ! EACH BAND SO JUST USE FIRST ONE.
         DO I = IS, IF
            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)
         END DO 
      END IF

      DO 400 I=IS,IF

C Compute the sea-ice albedos, snow water equivalent and run-off in this loop. Needed
C to perform this loop first so that appropriate fields can be passed to snow albedo
C look up tables.

C
C       * INITIALIZE VARIABLES
C
         ALSNOA(I,:) = 0.0
         ALSNOC(I,:) = 0.0
         C_FLAG(I)   = 0
C
C       * FRACTIONAL SNOW COVERAGE.
C
        ZSNOLIM=0.10
        IF(SNO(I).GT.0.0) THEN
          IF(RHON(I).EQ.0.)THEN
            SNOFAC(I)=1.0
          ELSE
            ZSNOW=SNO(I)/RHON(I)
            IF(ZSNOW.GE.(ZSNOLIM-0.00001)) THEN
              SNOFAC(I)=1.0
            ELSE
              SNOFAC(I)=ZSNOW/ZSNOLIM
            ENDIF
          ENDIF
        ELSE
           SNOFAC(I)=0.0
        ENDIF

        IF(GC(I).EQ.1.)                                           THEN
C                                                                              
C         * PACK (BARE ICE) ALBEDOES HAVE VALUES IN EACH SPECTRAL INTERVAL   
C         * WHICH FOLLOW FROM EBERT AND CURRY).
C         * THESE ARE IN EFFECT UNTIL "PUDDLING" OCCURS, WHICH                
C         * IS INCLUDED TO REDUCE ALBEDO IN EACH SPECTRAL INTERVAL,            
C         * ONLY IF GTA LIES BETWEEN GTFSW AND GTMSI. THE EFFECT               
C         * OF SNOW IS THEN SUPERIMPOSED, BASED ON THE ABOVE VALUE, THE        
C         * EFFECTIVE SNOW COVERAGE (SNOFAC) AND THE ABOVE-CALCULATED DEEP     
C         * SNOW ALBEDOES.                                                     
C                                                        
            HICE=SIC(I)/DENI
            IF (HICE .GE. 2.)                                      THEN
               ALBBI(I,1)    =  0.778
               ALBBI(I,2)    =  0.443
            ELSE IF (HICE .LT. 2. .AND. HICE .GE. 1.)              THEN
               ALBBI(I,1)    =  0.760 + 0.018 * (HICE - 1.0)
               ALBBI(I,2)    =  0.247 + 0.196 * (HICE - 1.0)
            ELSE
               X = LOG(MAX(HICE,0.01) + 1.E-10)
               ALBBI(I,1)    =  0.760 + 0.140 * X
               ALBBI(I,2)    =  0.247 + 0.029 * X
            ENDIF
            ALBBI(I,3)      =  0.055
            ALBBI(I,4)      =  0.036
            
            IF (ISNOALB .EQ. 1)                                    THEN
               IF (SNO(I) .GT. 0.0)                                THEN
                  C_FLAG(I) = 1
               ELSE
                  C_FLAG(I) = 0
               END IF
               SUM_C_FLAG = SUM_C_FLAG + C_FLAG(I)
            ELSE IF (ISNOALB .EQ. 0) THEN
C                                                                             
C         * ALBEDOES OVER SNOW ARE AGED AS IN "CLASS".
C         * FOR NOW, UNTIL A 4-BAND SNOW ALBEDO MODEL IS DEVELOPED,
C         * IT IS ASSUMED THAT THE ALBEDO FOR THE LAST 3 BANDS
C         * (NEAR-IR) ARE THE SAME.
C     
               IF(ALBS(I).GT.0.70) THEN                                          
                  TIMFAC=(ALBS(I)-0.70)/0.14                                       
                  ALSNOA(I,1)=0.135*TIMFAC+0.84                                         
                  ALSNOA(I,2)=0.272*TIMFAC+0.56                             
               ELSE                                                                
                  TIMFAC=(ALBS(I)-0.50)/0.34                                       
                  ALSNOA(I,1)=0.365*TIMFAC+0.61                                         
                  ALSNOA(I,2)=0.452*TIMFAC+0.38                          
               ENDIF 
               ALSNOA(I,3)=0.25
               ALSNOA(I,4)=0.025
C
               ALSNOC(I,1)=ALSNOA(I,1)         
               ALSNOC(I,2)=ALSNOA(I,2)         
               ALSNOC(I,3)=ALSNOA(I,3)       
               ALSNOC(I,4)=ALSNOA(I,4)         
                
            END IF
         END IF
 400   CONTINUE
       
       IF (ISNOALB .EQ. 1 .AND. SUM_C_FLAG .GT. 0) THEN
! Convert the units of the snow grain size and BC mixing ratio
! Snow grain size from meters to microns and BC from kg BC/kg SNOW to ng BC/kg SNOW
          DO I = IS,IF
             REFSNO(I) = REFSN(I)*1.0E6
             IF (RHON(I) .GT. 0.0) THEN
                BCSNO(I)  = (BCSN(I)/RHON(I))*1.0E12 
             ELSE
                BCSNO(I)  = BCSN(I)
             END IF
          END DO ! I
          
          CALL SNOW_ALBVAL(ALSNODIF, ! OUTPUT
     1                     ALSNODIR,
     2                     CSZ,      ! INPUT
     3                     ALBBI,
     4                     BCSNO,
     5                     REFSNO,
     6                     SNO,
     7                     C_FLAG,
     7                     IS,
     8                     IF,
     9                     ILG,
     1                     NBS )
          DO 420 IB = 1, NBS
             DO 450 I=IS,IF
                IF (C_FLAG(I) .EQ. 1) THEN
                   ALSNOA(I,IB) = ALSNODIF(I,IB)*WDIFFUA(I,IB)
     +                          + ALSNODIR(I,IB)*WDIRCTA(I,IB)
                   ALSNOC(I,IB) = ALSNODIF(I,IB)*WDIFFUC(I,IB)
     +                          + ALSNODIR(I,IB)*WDIRCTC(I,IB)
                ELSE
                   ALSNOA(I,IB) = 0.0
                   ALSNOC(I,IB) = 0.0
                END IF
 450         CONTINUE
 420      CONTINUE
       END IF
C
      DO 500 I=IS,IF
C
C       * ALBEDO CALCULATIONS (4 BANDS).
C
        IF(GC(I).EQ.1.)                                    THEN
C
C         * ESTIMATE MELT POND FRACTION ('PUDDLE') FROM SURFACE AIR TEMP.
C         * (NOTE MAXIMUM POND FRACTION IS 0.25, AS IN EBERT AND CURRY, 1993)
C
          PUDDLE = MAX(0., (MIN((GTA(I)-GTFSW)/4.,PUDDLE_MAX)))
C
C         * NOW TAKE TOTAL ICE ALBEDO TO BE MIX OF BARE ICE AND POND
C         * ALBEDO WEIGHTED BY MELT POND FRACTION.
C
          DPUDDLE = 1.-PUDDLE
          ALICE1  = DPUDDLE*ALBBI(I,1) + PUDDLE*ALBP1
          ALICE2  = DPUDDLE*ALBBI(I,2) + PUDDLE*ALBP2
          ALICE3  = DPUDDLE*ALBBI(I,3) + PUDDLE*ALBP3
          ALICE4  = DPUDDLE*ALBBI(I,4) + PUDDLE*ALBP4
C
C         * NOW OBTAIN AN EFFECTIVE ICE-SNOW ALBEDO BASED ON "SNOFAC"
C         * AS A FRACTIONAL COVERAGE (FOR BOTH ALL-SKY AND CLEAR-SKY).
C
          DSNOFAC   = 1.-SNOFAC(I)  
          SALB(I,1) = ALICE1*DSNOFAC + ALSNOA(I,1)*SNOFAC(I)
          SALB(I,2) = ALICE2*DSNOFAC + ALSNOA(I,2)*SNOFAC(I)
          SALB(I,3) = ALICE3*DSNOFAC + ALSNOA(I,3)*SNOFAC(I)
          SALB(I,4) = ALICE4*DSNOFAC + ALSNOA(I,4)*SNOFAC(I)
          CSAL(I,1) = ALICE1*DSNOFAC + ALSNOC(I,1)*SNOFAC(I)
          CSAL(I,2) = ALICE2*DSNOFAC + ALSNOC(I,2)*SNOFAC(I)
          CSAL(I,3) = ALICE3*DSNOFAC + ALSNOC(I,3)*SNOFAC(I)
          CSAL(I,4) = ALICE4*DSNOFAC + ALSNOC(I,4)*SNOFAC(I)
        ENDIF
  500 CONTINUE
C
      RETURN
      END 
