      SUBROUTINE BCAERO4 (EXTA, EXOMA, EXOMGA, FA, ABSA,
     1                    EXTA055, OMA055, GA055, 
     2                    EXTB, ODFB, SSAB, ABSB, EXTB055, ODF055,
     3                    SSA055, ABS055, BCOLOAD, BCYLOAD, CLD, WCDW,
     4                    IL1, IL2, ILG, LAY)
C
C     * JUN 14,2013 - J. LI.     NEW VERSION FOR GCM17+:
C     *                          - ADD SEMI-DIRECT EFFECT INSIDE CLOUD
C     *                            (REQUIRES PASSING AND USING CLD AND WCDW).
C     * APR 23,2012 - J. LI & Y. PENG. PREVIOUS VERSION BCAERO3 FOR GCM16:
C     *                                - ADD OPTICAL PROPERTY AT 0.55 UM 
C     *                                  AND RECALCULTE SOLAR BAND 
C     *                                  RESULT, IN THE PREVIOUS VERSION
C     *                                  THE ACC AND COA ARE SET OPPOSITE. 
C     * FEB 09,2009 - M.LAZARE. PREVIOUS VERSION BCAERO2 FOR GCM15H/I:
C     *                         - REMOVE UNUSED CODE RELATED TO RH EFFECT.
C     * SEP/13 2007 - X. MA/    PREVIOUS VERSION BCAERO FOR GCM15F:
C     *               J. LI.    - SEPERATE THE CODE FOR BC AND OC
C     *                           INDIVIDUALLY.
C     *                         - MODIFICATION TO MATCH MORE WITH
C     *                           OTHER PART OF AEROSOL AND SIMPLIFY 
C     *                           THE CODE.
C     * JUNE 08,2004 - DOMINIQUE BAEUMER: ORIGINAL CODING BASED ON: 
C                      D. BAEUMER, ET AL. J. GEOPHYS, RES. 2007 D10207
C
C----------------------------------------------------------------------C
C     CALCULATION OF OPTICAL PROPERTIES FOR BC (50% ACC, 50% CM)C
C     THIS FRACTION CAN BE CHANGED IF MORE INFO IS PROVIDED!           C
C                                                                      C
C     EXTA:      EXTINCTION COEFFICIENT                                C
C     EXOMA:     EXTINCTION COEFFICIENT TIMES SINGLE SCATTERING ALBEDO C
C     EXOMGA:    EXOMA TIMES ASYMMETRY FACTOR                          C
C     FA:        SQUARE OF ASYMMETRY FACTOR                            C
C     ABSA:      ABSORPTION COEFFICIENT                                C
C     BCLOAD:    BC (TOTAL) AEROSOL LOADING FOR EACH LAYER             C
C     BCOLOAD:   BCO (HYDROPHOBIC) AEROSOL LOADING FOR EACH LAYER      C
C     BCYLOAD:   BCY (HYDROPHYLIC) AEROSOL LOADING FOR EACH LAYER      C
C                                                                      C
C     BCTLDACC:  SCALED BCLOAD FOR ACCUM MODE                          C
C     BCTLDCOA:  SCALED BCLOAD FOR COARSE MODE                         C
C                                                                      C
C     SPHOBC:    SOLAR OPTICAL DATA FOR HYDROPHOBIC BC,                C
C                (NBS,6), INDEX NBS FOR SOLAR BAND, 6: 1, 2, 3 FOR EXT,C
C                SINGLE SCATTERING ALBEDO (OMEGA), ASYMMETRY FACTOR (G)C
C                OF ACCUM MODE, 4, 5, 6 FOR COARSE MODE                C
C     IPHOBC:    INFRARED OPTICAL DATA FOR HYDROPHOBIC BC,             C
C                (NBLM1,2), INDEX NBL FOR INFRARED BAND, 2: 1 FOR ACCUMC
C                2 FOR CORASE                                          C
C                                                                      C
C     FBCACC: FRACTION OF BC IN ACCUM MODE                             C
C----------------------------------------------------------------------C
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9, NBLM1 = NBL - 1, 
     1           FBCACC = 1.0, FBCCOA = 1.0 - FBCACC)
C
C     * OUTPUT ARRAYS.
C
      REAL EXTA(ILG,LAY,NBS), EXOMA(ILG,LAY,NBS), EXOMGA(ILG,LAY,NBS),
     1     FA  (ILG,LAY,NBS), ABSA (ILG,LAY,NBL)
      REAL ABSB(ILG,LAY,NBS), ODFB(ILG,LAY,NBS),
     1     EXTB(ILG,LAY,NBS), SSAB(ILG,LAY,NBS)
      REAL EXTA055(ILG,LAY), OMA055(ILG,LAY), GA055(ILG,LAY) 
      REAL EXTB055(ILG,LAY), ODF055(ILG,LAY), SSA055(ILG,LAY),
     1     ABS055( ILG,LAY), CLD(ILG,LAY), WCDW(ILG,LAY)
C
C     * INPUT ARRAYS.
C
      REAL BCLOAD(ILG,LAY),BCOLOAD(ILG,LAY),BCYLOAD(ILG,LAY)
C
C     * INTERNAL WORK ARRAYS.
C
      REAL BCTLDACC(ILG,LAY), BCTLDCOA(ILG,LAY)

      REAL SPHOBC(NBS,6), IPHOBC(NBLM1,2)
      REAL SE055(2), SW055(2), G055(2)

      COMMON /PARAM4/ ACTFRC
C      
C  *  HYDROPHOBIC DATA  *
C 
C  *  CORRECTED DATA BY JIANGNAN
C
      DATA ((SPHOBC(I,J), I = 1, 4), J = 1, 6)  /
     1 11.309423, 5.300023, 2.691659, 1.354601, 
     2  0.227982, 0.121912, 0.048024, 0.009951,
     3  0.357685, 0.252115, 0.160571, 0.076695,
     1  8.700115, 6.658090, 4.179159, 2.025891,
     2  0.444720, 0.401069, 0.327568, 0.208560,
     3  0.676934, 0.574040, 0.459321, 0.315548  /
C
      DATA ((IPHOBC(I,J), I = 1, NBLM1), J = 1, 2)                     /
     1 1.03056, .87149, .68762, .48330, .40322, .33417, .23639, .14938,
     2  .87573, .74959, .60427, .43853, .37152, .31304, .22652, .14601 /
C
      DATA (SE055(J), J = 1, 2)            /
     1 9.4789, 8.4650                      /
C
      DATA (SW055(J), J = 1, 2)            /
     1 0.2093, 0.4390                      / 
C
      DATA (G055(J), J = 1, 2)             /
     1 0.3348, 0.6577                      / 
C
C----------------------------------------------------------------------C
C     10000. IS THE UNIT SCALED DISCUSSED IN OTHER AEROSOL ROUTINES    C
C     FOR 0.55 UM THE ACCUMULATION MODE USING SIGMA = 2, GEOMETRY      C
C     RADIUS R0 = 0.0118 UM (REFF = 0.03922), COARSE MODE SIGMA = 2,   C
C     R0 = 0.06017 (REFF = 0.2 UM), ACCUMULATION MODE SIZE FOLLOWS     C
C     HESS ETC BUL OF AMER METEOR SOC 1998                             C
C----------------------------------------------------------------------C
C
      DO 100 K = 1, LAY
      DO 100 I = IL1, IL2
        BCLOAD(I,K)    =  BCOLOAD(I,K) + BCYLOAD(I,K)
        X              =  10000.0 * BCLOAD(I,K)
C
C----------------------------------------------------------------------C
C     ACTFRC% BC IS CONSIDERED IN THE INTERNAL MIXING OF CLOUD,        C
C     1-ACTFRC IS TREATED AS EXTERNAL MIXING INSIDE CLOUD,             C
C     1.0 - ACTFRC * CLD = (1.0 - ACTFRC) * CLD + 1.0 - CLD            C
C----------------------------------------------------------------------C
C
        IF (CLD(I,K) .GT. 0.01 .AND. WCDW(I,K) .GT. 0.001)          THEN
          Y            =  10000.0 * BCYLOAD(I,K)
          Z            =  10000.0 * BCOLOAD(I,K)
          EFCLD        =  1.0 - ACTFRC * CLD(I,K)
          BCTLDACC(I,K)=  EFCLD * FBCACC * Y + FBCACC * Z
          BCTLDCOA(I,K)=  EFCLD * FBCCOA * Y + FBCCOA * Z
        ELSE
          BCTLDACC(I,K)=  FBCACC * X
          BCTLDCOA(I,K)=  FBCCOA * X
        ENDIF
  100 CONTINUE
C
C----------------------------------------------------------------------C
C     SOLAR                                                            C
C     FOR BC THE HYDROPHYLIC EFFECT ON OPTICAL PROPERTY IS TOO SMALL   C
C     NOT ACCOUNTED FOR BOTH SOLAR AND INFRARED                        C
C----------------------------------------------------------------------C
C
      DO 200 J = 1, NBS
      DO 200 K = 1, LAY
      DO 200 I = IL1, IL2
C
C --- BC ACCUMULATION MODE ---
C
        EXT            =  BCTLDACC(I,K) * SPHOBC(J,1)
        EXTA(I,K,J)    =  EXTA(I,K,J) + EXT
        EXOM           =  EXT * SPHOBC(J,2)
        EXOMA(I,K,J)   =  EXOMA(I,K,J) + EXOM
        EXOMG          =  EXOM * SPHOBC(J,3)
        EXOMGA(I,K,J)  =  EXOMGA(I,K,J) + EXOMG
        EXOMGG         =  EXOMG * SPHOBC(J,3)
        FA(I,K,J)      =  FA(I,K,J) + EXOMGG
C
        ABS1           =  EXT * (1.0 - SPHOBC(J,2))
        EXT1           =  EXT
        SSA1           =  EXOM
C
C --- BC COARSE MODE ---
C
        EXT            =  BCTLDCOA(I,K) * SPHOBC(J,4)
        EXTA(I,K,J)    =  EXTA(I,K,J) + EXT
        EXOM           =  EXT * SPHOBC(J,5)
        EXOMA(I,K,J)   =  EXOMA(I,K,J) + EXOM
        EXOMG          =  EXOM * SPHOBC(J,6)
        EXOMGA(I,K,J)  =  EXOMGA(I,K,J) + EXOMG
        EXOMGG         =  EXOMG * SPHOBC(J,6)
        FA(I,K,J)      =  FA(I,K,J) + EXOMGG
C
        ABS2           =  EXT * (1.0 - SPHOBC(J,5))
        ABSB(I,K,J)    =  ABS1 + ABS2
        EXT2           =  EXT
        SSA2           =  EXOM
C
        TMASS          =  10000.0 * BCLOAD(I,K) 
        IF (TMASS .GT. 1.E-15) THEN
         EXTB(I,K,J)   = (EXT1+EXT2)/ TMASS
         ODFB(I,K,J)   =  EXT1 / BCTLDACC(I,K)
        ELSE
         EXTB(I,K,J)   =  0.0
         ODFB(I,K,J)   =  0.0
        ENDIF
C
        IF ((EXT1+EXT2) .GT. 1.E-20) THEN
         SSAB(I,K,J)   = (SSA1+SSA2)/(EXT1+EXT2)
        ELSE
         SSAB(I,K,J)   =  0.0
        ENDIF
C
  200 CONTINUE
C
      DO 250 K = 1, LAY
      DO 250 I = IL1, IL2
        EXT0551        =  BCTLDACC(I,K) * SE055(1) 
        EXT0552        =  BCTLDCOA(I,K) * SE055(2)
        EXTA055(I,K)   =  EXTA055(I,K) + EXT0551 + EXT0552
        OM0551         =  EXT0551 * SW055(1)
        OM0552         =  EXT0552 * SW055(2)
        OMA055(I,K)    =  OMA055(I,K) + OM0551 + OM0552
        GA0551         =  OM0551 * G055(1) * G055(1)
        GA0552         =  OM0552 * G055(2) * G055(2)
        GA055(I,K)     =  GA055(I,K) + GA0551 + GA0552
C
        ABS055(I,K)    =  EXT0551*(1.0-SW055(1))+
     1                    EXT0552*(1.0-SW055(2))
        TMASS          =  10000.0 * BCLOAD(I,K) 
        IF (TMASS .GT. 1.E-15) THEN
         EXTB055(I,K)  = (EXT0551+EXT0552) / TMASS
         ODF055(I,K)   =  EXT0551 / BCTLDACC(I,K)
        ELSE
         EXTB055(I,K)  =  0.0
         ODF055(I,K)   =  0.0
        ENDIF
C
        IF ((EXT0551+EXT0552) .GT. 1.E-20) THEN
         SSA055(I,K)   = (OM0551 + OM0552) / (EXT0551+EXT0552)
        ELSE
         SSA055(I,K)   =  0.0
        ENDIF
C
  250 CONTINUE
C
C----------------------------------------------------------------------C
C     INFRARED                                                         C
C----------------------------------------------------------------------C
C
      DO 300 J  = 1, NBLM1
      DO 300 K  = 1, LAY
      DO 300 I  = IL1, IL2
          ABSA(I,K,J)  =  ABSA(I,K,J) + 
     1                    BCTLDACC(I,K) * IPHOBC(J,1) +
     2                    BCTLDCOA(I,K) * IPHOBC(J,2)
 300  CONTINUE      
C
      RETURN
      END
