      SUBROUTINE OCAERO4 (EXTA, EXOMA, EXOMGA, FA, ABSA,
     1                    EXTA055, OMA055, GA055,
     2                    EXTB, ODFB, SSAB, ABSB, EXTB055, ODF055,
     3                    SSA055, ABS055, RH, OCOLOAD, OCYLOAD,
     4                    IL1, IL2, ILG, LAY)
C
C     * FEB 08, 2013 - J.LI.  NEW VERSION FOR GCM17 BASED ON OCAERO3.
C     *                       - UPDATE TO READING DATA, 
C     *                         SIMPLIFY THE CODE AND TOOK OUT
C     *                         THE COARSE MODE, SINCE NO INPUT 
C     *                         FROM GCM.
C     * APR 23, 2012 - J.LI & Y.PENG. PREVIOUS VERSION OCAERO3 FOR GCM16:
C                                     - REVISED CODE FOR OPTICAL OF OC.
C     * FEB 10, 2009 - J. LI/ PREVIOUS VERSION OCAERO2 FOR GCM15H/I:
C     *                X. MA. - INCLUDE RH EFFECT FOR HYDROPHYLIC OC.
C     * SEP 26,2007 - X. MA/  PREVIOUS VERSION OCAERO FOR GCM15F:
C     *               J. LI.  - CORRECT THE WRONG OC INPUT FOR HYDROPHOBIC
C     *                         AT INFRARED
C     *                       - SEPARATE THE CODE FOR BC AND OC.
C     *                       - MODIFICATION TO MATCH MORE WITH OTHER PART
C     *                         OF AEROSOL AND SIMPLIFY THE CODE
C----------------------------------------------------------------------C
C     CALCULATION OF OPTICAL PROPERTIES FOR OC WITH TWO MODES, FINE    C
C     MODE: SIGMA = 2.24, R0 = 0.0212, COARSE MODE: SIGMA = 2., R0 =   C
C     0.0602. NOTE IN GCM4 WE ONLY HAVE THE ACCUM MODE, WE KEEP THE    C
C     COARSE MODE DATA, BUT NOT USE IT NOW                             C
C     THE REFRACTIVE INDEX IS FROM READING                             C
C     THE GROWTH IS FOLLOWING KNUT'S SCHEME USED IN AEROMIX3           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     RH:     RELATIVE HUMIDITY                                        C
C     OCOLOAD:OCO AEROSOL LOADING FOR EACH LAYER                       C
C     OCYLOAD:OCY AEROSOL LOADING FOR EACH LAYER                       C
C     EXTA055, OMA055, GA055 ARE THE INFO FOR 0.55                     C
C     FOLLOWING INFO IS NO USEFUL EXCEPT 0.55  (I WILL TAKE OUT ALL    C
C     EXTB, SSAB, ABSB ARE BAND INFO, ODFB IS EXTB FOR FINE MODE       C
C     IN SE1 FIRST 4 IS FOR BANDS, 5 IS FOR MASS, 6 IS FOR 0.55, 7 IS  C
C     0.83. IN SW1 AND SG1, THE 5 IS FOR 0.55                          C
C----------------------------------------------------------------------C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9, NH = 7)
C
      REAL, DIMENSION(ILG,LAY,NBS):: EXTA, EXOMA, EXOMGA, FA, ABSB, 
     1                               ODFB, EXTB, SSAB
      REAL, DIMENSION(ILG,LAY,NBL):: ABSA
      REAL, DIMENSION(NBS):: SOMG, SG, SEXTA, SOMGA, SGA
      REAL, DIMENSION(7,NH):: SE1
      REAL, DIMENSION(5,NH):: SW1, SG1
      REAL, DIMENSION(NBL,NH):: SA1
      REAL, DIMENSION(ILG,LAY):: RH, OCYLOAD, OCOLOAD, OCOLDACC,
     1                           OCYLDACC, OCLDACC, DELTA, RMASS, 
     2                           EXTA055, OMA055, GA055, EXTB055, 
     3                           ODF055, SSA055, ABS055
      INTEGER NUMRH(ILG,LAY)
      REAL RHSTD(7)
C
        DATA ((SE1(I,J), I = 1, 7), J = 1, NH)                       /
     1  5.27833, 2.24282, 0.74013, 0.10034, 1.01570, 4.74178, 2.21780,
     2  5.17586, 2.21064, 0.73347, 0.12568, 1.14027, 4.64903, 2.18660,
     3  5.00877, 2.16995, 0.72961, 0.16993, 1.41475, 4.50327, 2.14763,
     4  4.76548, 2.15794, 0.75460, 0.24868, 2.20339, 4.31103, 2.13947,
     5  4.61216, 2.21680, 0.81762, 0.31564, 3.43786, 4.21564, 2.20286,
     6  4.50408, 2.29541, 0.89280, 0.36684, 4.97494, 4.16079, 2.28606,
     7  4.39900, 2.36639, 0.96723, 0.40893, 6.76721, 4.10382, 2.36177/

        DATA ((SW1(I,J), I = 1, 5), J = 1, NH)     /
     1  0.82337, 0.99565, 0.99390, 0.93582, 0.89469,
     2  0.83723, 0.99608, 0.99423, 0.75423, 0.90468,
     3  0.86112, 0.99679, 0.99483, 0.57845, 0.92101,
     4  0.90214, 0.99791, 0.99594, 0.45321, 0.94700,
     5  0.93271, 0.99868, 0.99682, 0.43287, 0.96498,
     6  0.95105, 0.99910, 0.99736, 0.44739, 0.97523,
     7  0.96239, 0.99935, 0.99770, 0.47112, 0.98137/

        DATA ((SG1(I,J), I = 1, 5), J = 1, NH)     /
     1  0.60207, 0.53575, 0.47566, 0.30760, 0.58512,
     2  0.62043, 0.55494, 0.49143, 0.31751, 0.60415,
     3  0.64999, 0.58655, 0.51907, 0.33807, 0.63505,
     4  0.69608, 0.63859, 0.56991, 0.38434, 0.68404,
     5  0.72850, 0.67871, 0.61442, 0.43241, 0.71933,
     6  0.74817, 0.70558, 0.64714, 0.47315, 0.74124,
     7  0.76078, 0.72447, 0.67172, 0.50713, 0.75557/

        DATA ((SA1(I,J), I = 1, NBL), J = 1, NH)                     /
     1  0.07194, 0.06362, 0.15064, 0.33950, 0.26811, 0.08958, 0.15957,
     1  0.13965, 0.13655,
     2  0.06767, 0.06088, 0.14525, 0.32087, 0.27992, 0.09750, 0.17296,
     2  0.14665, 0.14213,
     3  0.06049, 0.05596, 0.13554, 0.28726, 0.28377, 0.10968, 0.19419,
     3  0.15527, 0.14563,
     4  0.04907, 0.04758, 0.11899, 0.22296, 0.24680, 0.12643, 0.22548,
     4  0.16216, 0.13838,
     5  0.04147, 0.04177, 0.10748, 0.16983, 0.19475, 0.13612, 0.24529,
     5  0.16288, 0.12664,
     6  0.03733, 0.03861, 0.10115, 0.13646, 0.15764, 0.14105, 0.25618,
     6  0.16219, 0.11824,
     7  0.03502, 0.03687, 0.09761, 0.11546, 0.13323, 0.14384, 0.26269,
     7  0.16155, 0.11274                                              /

      DATA RHSTD /0.1, 0.5, 0.75, 0.9, 0.95, 0.97, 0.98/
C
C----------------------------------------------------------------------C
      DO 100 K = 1, LAY
      DO 100 I = IL1, IL2
        RHS       =  MAX( MIN (RH(I,K), 0.98), 0.1)
        IF (RHS .GE. 0.1 .AND. RHS .LT. 0.50)                       THEN
          L           =  1
        ELSEIF (RHS .GE. 0.50 .AND. RHS .LT. 0.75)                  THEN
          L           =  2
        ELSEIF (RHS .GE. 0.75 .AND. RHS .LT. 0.90)                  THEN
          L           =  3
        ELSEIF (RHS .GE. 0.90 .AND. RHS .LT. 0.95)                  THEN
          L           =  4
        ELSEIF (RHS .GE. 0.95 .AND. RHS .LT. 0.97)                  THEN
          L           =  5
        ELSEIF (RHS .GE. 0.97 .AND. RHS .LE. 0.98)                  THEN
          L           =  6
        ENDIF
        NUMRH(I,K)    =  L
        LP1           =  L + 1
        DELTA(I,K)    = (RH(I,K) - RHSTD(L)) / (RHSTD(LP1) - RHSTD(L))
C
C *   MASS RATIO
        RMASS(I,K)    =  SE1(5,L) + DELTA(I,K) * (SE1(5,LP1) - SE1(5,L))
C
        OCOLDACC(I,K) =  10000.0 * OCOLOAD(I,K)
        OCYLDACC(I,K) =  10000.0 * OCYLOAD(I,K)
        OCLDACC(I,K)  =  OCOLDACC(I,K) + OCYLDACC(I,K)
  100 CONTINUE
C
      DO 200 J = 1, NBS
      DO 200 K = 1, LAY
      DO 200 I = IL1, IL2
        IF (OCLDACC(I,K) .LT. 1.E-09)                           GOTO 200
        L             =  NUMRH(I,K)
        LP1           =  L + 1
        EXTY          =  OCYLDACC(I,K) * RMASS(I,K) * (SE1(J,L) + 
     1                   DELTA(I,K) * (SE1(J,LP1) - SE1(J,L)))
        EXTO          =  OCOLDACC(I,K) * SE1(J,1)
C
        OMY           =  SW1(J,L) + DELTA(I,K) * (SW1(J,LP1) - SW1(J,L))
        OMO           =  SW1(J,1)
C
        GY            =  SG1(J,L) + DELTA(I,K) * (SG1(J,LP1) - SG1(J,L))
        GO            =  SG1(J,1)
        EXOMY         =  EXTY * OMY
        EXOMO         =  EXTO * OMO
        EXOMGY        =  EXTY * OMY * GY
        EXOMGO        =  EXTO * OMO * GO
C
        EXTA(I,K,J)   =  EXTA(I,K,J) + EXTY + EXTO
        EXOMA(I,K,J)  =  EXOMA(I,K,J) + EXOMY + EXOMO
        EXOMGA(I,K,J) =  EXOMGA(I,K,J) + EXOMGY + EXOMGO
        FA(I,K,J)     =  FA(I,K,J) + EXOMGY * GY + EXOMGO * GO
C
        ABSB(I,K,J)   =  EXTY * (1.0 - OMY) + EXTO * (1.0 - OMO)
        SSAB(I,K,J)   =  EXOMY + EXOMO
        EXTB(I,K,J)   =  EXTY  + EXTO
        ODFB(I,K,J)   =  EXTB(I,K,J)
 200  CONTINUE
C
      DO 250 K = 1, LAY
      DO 250 I = IL1, IL2
        IF (OCLDACC(I,K) .LT. 1.E-09)                           GOTO 250
        L             =  NUMRH(I,K)
        LP1           =  L + 1
        EXTY          =  OCYLDACC(I,K) * RMASS(I,K) * (SE1(6,L) +
     1                   DELTA(I,K) * (SE1(6,LP1) - SE1(6,L)))
        EXTO          =  OCOLDACC(I,K) * SE1(6,1)
        EXTA055(I,K)  =  EXTA055(I,K) + EXTY + EXTO
        ODF055(I,K)   =  EXTA055(I,K)
C
        OMY           =  SW1(5,L) + DELTA(I,K) * (SW1(5,LP1) - SW1(5,L))
        OMO           =  SW1(5,1)
        OMA055(I,K)   =  OMA055(I,K) + OMY + OMO
C
        GY            =  SG1(5,L) + DELTA(I,K) * (SG1(5,LP1) - SG1(5,L))
        GO            =  SG1(5,1)
        GA055(I,K)    =  GA055(I,K) + GY + GO
C
        ODF055(I,K)   =  EXTA055(I,K)
        ABS055(I,K)   = (1.0 - OMY) * EXTY + (1.0 - OMO) * EXTO

        TMASS         =  OCYLDACC(I,K) + OCOLDACC(I,K)
        IF (TMASS .GT. 1.E-15) THEN
          EXTB055(I,K)= (EXTY + EXTO) / TMASS
          ODF055(I,K) =  EXTB055(I,K)
        ELSE
          EXTB055(I,K)=  0.0
          ODF055(I,K) =  0.0
        ENDIF
C
        IF (EXTY + EXTO .GT. 1.E-20)                                THEN
         SSA055(I,K)  =  (EXTY * OMY + EXTO * OMO) / (EXTY + EXTO)
        ELSE
         SSA055(I,K)  =  0.0
        ENDIF
C
 250  CONTINUE
C
C----------------------------------------------------------------------C
C     INFRARED, SA: SPECIFIC ABSORPTANCE, ABS: ABSORPTION              C
C----------------------------------------------------------------------C
C
      DO 400 J = 1, NBL 
      DO 400 K = 1, LAY
      DO 400 I = IL1, IL2
        IF (OCLDACC(I,K) .LT. 1.E-09)                           GOTO 400
        L             =  NUMRH(I,K)
        LP1           =  L + 1
        ABSY          =  OCYLDACC(I,K) * (SA1(J,L) +
     1                   DELTA(I,K) * (SA1(J,LP1) - SA1(J,L)))
        ABSO          =  OCOLDACC(I,K) * SA1(J,1)

        ABSA(I,K,J)   =  ABSA(I,K,J) + ABSY + ABSO
 400  CONTINUE
C
      RETURN
      END
