      SUBROUTINE SULFAERO5 (EXTA, EXOMA, EXOMGA, FA, ABSA,
     1                      EXTA055, OMA055, GA055, 
     2                      EXTB, ODFB, SSAB, ABSB, EXTB055,ODF055, 
     3                      SSA055, ABS055, RHIN, ALOAD, 
     4                      IL1, IL2, ILG, LAY) 
C
C     * APR 23,2012 - J.LI & Y.PENG. NEW VERSION FOR GCM16:
C     *                         - ADD THE OPTICAL DEPTH CALCULATION
C     *                           FOR 0.55 UM
C     * FEB 13,2009 - M.LAZARE. PREVIOUS VERSION SULFAERO4 FOR GCM15H/I:
C     *                         RHMAX LOWERED FROM 0.98 TO 0.95. 
C     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION SULFAERO3 FOR GCM15G:
C     *                         - HF,RS2,FR1,FR2,FR3,RE,RH NOW
C     *                           INTERNAL WORK ARRAYS.
C     * JUN 19/06 -  M.LAZARE. PREVIOUS VERSION SULFAERO2 FOR GCM15F:
C     *                        - COSMETIC: USE VARIABLE INSTEAD OF 
C     *                          CONSTANT IN INTRINSICS SUCH AS "MAX",
C     *                          SO THAT CAN COMPILE IN 32-BIT MODE
C     *                          WITH REAL*8.  
C     * MAY 14,2003 - J.LI.    PREVIOUS VERSION SULFAERO FOR UP TO
C     *                        AND INCLUDING GCM15E.  
C
C----------------------------------------------------------------------C
C     CALCULATION OF OPTICAL PROPERTIES FOR SULFATE AEROSOL, H2SO4,    C
C     (NH4)2SO4 AND NH4HSO4 (LI ET AL 2001 JAS P193; LI AND MIN 2002   C
C     JAS P3130)                                                       C
C     S1: H2SO4, S2: (NH4)2SO4, S3: NH4HSO4,                           C
C     9: 1-3 FOR RE = 0.166 UM; 4-6 FOR RE = 0.5 UM; 7-9 FOR RE = 1.0  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     RHIN:   INPUT RELATIVE HUMIDITY                                  C
C     ALOAD:  SULFATE AEROSOL LOADING FOR EACH LAYER                   C
C     RE:     EFFECTIVE RADIUS                                         C
C     RH:     RELATIVE HUMIDITY                                        C
C     FR1-3:  INTERPOLATION FACTORS                                    C
C----------------------------------------------------------------------C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9, NBLM1 = 8)
C
C     * OUTPUT ARRAYS.
C
      REAL   EXTA(ILG,LAY,NBS), EXOMA(ILG,LAY,NBS), EXOMGA(ILG,LAY,NBS),
     1       FA  (ILG,LAY,NBS)
      REAL   ABSB(ILG,LAY,NBS),  ODFB(ILG,LAY,NBS),
     1       EXTB(ILG,LAY,NBS),  SSAB(ILG,LAY,NBS)
      REAL   ABSA (ILG,LAY,NBL)
C
      REAL   EXTA055(ILG,LAY),OMA055(ILG,LAY),GA055(ILG,LAY), 
     1       EXTB055(ILG,LAY),ODF055(ILG,LAY),SSA055(ILG,LAY),
     2       ABS055(ILG,LAY)
C
C     * INPUT ARRAYS.
C
      REAL   RHIN(ILG,LAY), ALOAD(ILG,LAY)
C
C     * INTERNAL WORK ARRAYS.
C
      REAL   HF(ILG,LAY), RS2(ILG,LAY), FR1(ILG,LAY), FR2(ILG,LAY), 
     1       FR3(ILG,LAY), RE(ILG,LAY), RH(ILG,LAY)
C
      REAL   SE2(9,NBS), CW2(9,NBS), G2(9,NBS), SA2(9,NBLM1)
      REAL   SE055(3), CW055, G055(2)
C
       DATA ((SE2(I,J), I = 1, 9), J = 1, NBS)                         /
     1  3.94   , 9.46E-1, 7.56E-2, 2.44   , 1.34E-1, 8.80E-2,
     2  1.16   , 1.85E-2, 4.07E-2, 1.64   , 5.55E-1,-3.05E-2,
     3  2.04   , 3.24E-1, 6.64E-2, 1.24   , 5.11E-2, 4.49E-2,
     4  5.27E-1, 1.38E-1,-4.33E-2, 1.22   , 3.20E-1, 2.00E-2,
     5  1.07   , 1.40E-1, 3.58E-2, 1.86E-1, 8.73E-2,-1.75E-2,
     6  4.51E-1, 1.86E-1,-7.23E-3, 5.97E-1, 1.59E-1, 1.13E-2           /
C
       DATA ((CW2(I,J), I = 1, 9), J = 1, NBS)                         /
     1  6.71E-7,-4.52E-7, 7.34E-9, 1.36E-6,-8.83E-7, 9.43E-9,
     2  2.93E-6,-1.84E-6, 2.51E-8, 2.89E-6, 2.73E-6,-7.21E-8,
     3  3.65E-6, 1.25E-6,-5.96E-7, 6.88E-6, 1.09E-6,-1.48E-6,
     4  1.49E-3, 6.05E-4, 2.99E-5, 9.64E-4, 2.89E-4,-5.47E-5,
     5  1.38E-3, 7.31E-5,-1.54E-4, 6.64E-1, 3.33E-2, 1.73E-2,
     6  3.05E-1, 4.78E-2, 7.74E-3, 2.19E-1, 3.67E-2, 3.41E-3           /
C
       DATA ((G2(I,J), I = 1, 9), J = 1, NBS)                          /
     1  .639, .116,-2.73E-3, .671, .100,-1.56E-3, .687, .084,-2.99E-3,
     2  .591, .114,-5.56E-3, .675, .107,-1.13E-3, .683, .091,-1.62E-3,
     3  .512, .108,-9.00E-3, .658, .111,-2.76E-3, .689, .102,-1.11E-3,
     4  .349, .067,-1.41E-2, .597, .100,-6.89E-3, .697, .097,-3.16E-3  /
C
       DATA ((SA2(I,J), I = 1, 9), J = 1, NBLM1)                       /
     1  1.20E-2,  8.16E-3, -5.83E-4, 1.68E-2,  1.04E-2, -6.24E-4,
     2  2.22E-2,  1.09E-2, -3.54E-4, 3.04E-1, -1.62E-1,  6.23E-3,
     3  2.74E-1, -1.23E-1,  6.16E-3, 2.05E-1, -6.27E-2,  5.43E-3,
     4  3.15E-1, -1.35E-1,  5.29E-3, 2.94E-1, -1.05E-1,  5.20E-3,
     5  2.32E-1, -5.06E-2,  5.07E-3, 3.54E-1, -2.11E-1,  5.25E-3,
     6  3.22E-1, -1.78E-1,  4.83E-3, 2.56E-1, -1.17E-1,  4.45E-3,
     7  2.34E-1,  4.69E-2,  1.34E-2, 3.03E-1, -1.35E-2,  1.37E-2,
     8  3.25E-1, -6.14E-2,  1.21E-2, 2.35E-3,  3.21E-2, -1.34E-3,
     9  1.36E-3,  3.60E-2, -9.26E-4, 3.30E-3,  3.59E-2, -2.26E-4,
     A  3.03E-2,  1.37E-1, -6.75E-3, 2.63E-2,  1.62E-1, -4.78E-3,
     B  3.06E-2,  1.75E-1, -7.06E-4, 5.75E-3,  1.01E-1, -3.85E-3,
     C  7.97E-3,  1.08E-1, -4.22E-3, 8.59E-3,  1.24E-1, -2.95E-3      /
C
       DATA (SE055(I), I =1, 3)                                       /
     1  3.44744,  0.88011, 0.38768E-2                                 /
C
       DATA CW055 /0.99999/
C
       DATA (G055(I), I=1, 2)                                         /
     1  0.61204,  1.6612E-01                                          /
C
      DATA RHMAX /0.95/
C
C----------------------------------------------------------------------C
C     LI ET AL (2001) CAN HANDLE S1: H2SO4, S2: (NH4)2SO4, S3: NH4HSO4,C
C     WE ONLY CONSIDER 2 HERE                                          C
C     RS: RATIO OF WET TO DRY SULFATE AEROSOL CONTENT, RE: EFFECTIVE   C
C     RADIUS OF SULFATE AEROSOL, FR1,2,3 ARE COEFFICIENTS FOR          C
C     LAGRANGIAN INTERPOLATION BASED ON RESULTS OF RE = 0.166, 0.5 AND C
C     1 UM.                                                            C
C----------------------------------------------------------------------C
C
C----------------------------------------------------------------------C
C     SOLAR, SE: SPECIFIC EXTINCTION, OMAS: SINGLE SCATTERING ALBEDO,  C
C     GAS: ASYMMETRY FACTOR, EXTAS: EXTINCTION COEFFICIENT             C
C     10000 IN RS2, BECAUSE THE UNIT OF SPECIFIC EXTINCTION FOR AEROSOLC
C     IS M^2/GRAM, WHILE FOR GAS IS CM^2/GRAM, IN RADDRIV WE SAME      C
C     DP (AIR DENSITY * LAYER THICKNESS) IS USED FOR BOTH GAS AND      C
C     AEROSOL. ALOAD IS PURE SO4 LOADING IN UNIT G (SULFUR) / G        C
C     (AIR). TAU = ALOAD * SPECIFIC EXTINCTION * AIR DENSITY *         C
C     LAYER THICKNESS                                                  C
C     4.12: CONVERT THE MASS OF SULFUR TO MASS OF (NH4)2SO4            C
C           (NOT FOR PLA MASS!!)                                       C
C----------------------------------------------------------------------C
C
      DO 100 K = 1, LAY
      DO 100 I = IL1, IL2
        RE(I,K)         =  0.166
        RH(I,K)         =  MIN (RHIN(I,K), RHMAX)
        HF(I,K)         =  1.0 / (RH(I,K) - 1.05)
        RS2(I,K)        =  10000.0 * EXP(-0.2373 + 1.135 * RH(I,K) +
     1                     0.01048 / ((RH(I,K) - 1.06) *
     2                    (RH(I,K) - 1.06)))
        RE0166          =  RE(I,K) - 0.166
        RE05            =  RE(I,K) - 0.5
        RE10            =  RE(I,K) - 1.0
        FR1(I,K)        =  3.5899424 * RE05 * RE10
        FR2(I,K)        =  - 5.988024 * RE0166 * RE10
        FR3(I,K)        =  2.3980815 * RE0166 * RE05
  100 CONTINUE
C
      DO 200 J = 1, NBS
      DO 200 K = 1, LAY
      DO 200 I = IL1, IL2
        SE20166         =  SE2(1,J) + SE2(2,J) * RH(I,K) +
     1                     SE2(3,J) * HF(I,K)
        SE205           =  SE2(4,J) + SE2(5,J) * RH(I,K) +
     1                     SE2(6,J) * HF(I,K)
        SE210           =  SE2(7,J) + SE2(8,J) * RH(I,K) +
     1                     SE2(9,J) * HF(I,K)
        EXT             =  4.12 * ALOAD(I,K) * RS2(I,K) *
     1                    (FR1(I,K) * SE20166 + FR2(I,K) * SE205 +
     2                     FR3(I,K) * SE210)
C
        CW20166         =  CW2(1,J) + CW2(2,J) * RH(I,K) +
     1                     CW2(3,J) * HF(I,K)
        CW205           =  CW2(4,J) + CW2(5,J) * RH(I,K) +
     1                     CW2(6,J) * HF(I,K)
        CW210           =  CW2(7,J) + CW2(8,J) * RH(I,K) +
     1                     CW2(9,J) * HF(I,K)
        OM              =  1.0 - FR1(I,K) * CW20166 - FR2(I,K) * CW205 -
     1                     FR3(I,K) * CW210
        EXOM            =  EXT * OM
C
        G20166          =  G2(1,J) + G2(2,J) * RH(I,K) +
     1                     G2(3,J) * HF(I,K)
        G205            =  G2(4,J) + G2(5,J) * RH(I,K) +
     1                     G2(6,J) * HF(I,K)
        G210            =  G2(7,J) + G2(8,J) * RH(I,K) +
     1                     G2(9,J) * HF(I,K)
        G               =  FR1(I,K) * G20166 + FR2(I,K) * G205 +
     1                     FR3(I,K) * G210
        EXOMG           =  EXOM * G
        EXTA(I,K,J)     =  EXTA(I,K,J) + EXT
        EXOMA(I,K,J)    =  EXOMA(I,K,J) + EXOM
        EXOMGA(I,K,J)   =  EXOMGA(I,K,J) + EXOMG
        FA(I,K,J)       =  FA(I,K,J) + EXOMG * G
C
        ABSB(I,K,J)     =  EXT * (1.0 - OM)
        TMASS           =  4.12*10000.0*ALOAD(I,K)
        IF (TMASS .GT. 1.E-15) THEN
         EXTB(I,K,J)    =  EXT / TMASS
        ELSE
         EXTB(I,K,J)    =  0.0
        ENDIF
        ODFB(I,K,J)     =  EXTB(I,K,J)
        SSAB(I,K,J)     =  OM
C
 200  CONTINUE
C
      DO 250 K = 1, LAY
      DO 250 I = IL1, IL2
C
C * 550 NM, WET APPROACH, SAME AS BAND CAL., VM
C
        EXT055          =  4.12*ALOAD(I,K) * RS2(I,K) *
     1                    (SE055(1) + SE055(2) * RH(I,K) + 
     2                     SE055(3) * HF(I,K))
        EXTA055(I,K)    =  EXTA055(I,K) + EXT055
        OMA055(I,K)     =  OMA055(I,K) + EXT055 * CW055
        GG055           =  G055(1) + G055(2) * RH(I,K)
        GA055(I,K)      =  GA055(I,K) + EXT055 * CW055 * GG055 **2
C
        ABS055(I,K)     =  EXT055 * (1.0 - CW055)
        TMASS           =  4.12*10000.0*ALOAD(I,K)
        IF (TMASS .GT. 1.E-15) THEN
         EXTB055(I,K)   =  EXT055/TMASS
        ELSE
         EXTB055(I,K)   =  0.0
        ENDIF
        ODF055(I,K)     =  EXTB055(I,K)
        SSA055(I,K)     =  CW055
C
 250  CONTINUE
C
C----------------------------------------------------------------------C
C     INFRARED, SA: SPECIFIC ABSORPTANCE, ABSAS: ABSORPTION            C
C     COEFFICIENT THE CONTRIBUTION OF BAND 9 IS VERY SMALL, NEGLECTED  C
C----------------------------------------------------------------------C
C
      DO 400 J = 1, NBLM1
      DO 400 K = 1, LAY
      DO 400 I = IL1, IL2
        SA20166         =  SA2(1,J) + SA2(2,J) * RH(I,K) +
     1                     SA2(3,J) * HF(I,K)
        SA205           =  SA2(4,J) + SA2(5,J) * RH(I,K) +
     1                     SA2(6,J) * HF(I,K)
        SA210           =  SA2(7,J) + SA2(8,J) * RH(I,K) +
     1                     SA2(9,J) * HF(I,K)
        ABSA(I,K,J)     =  ABSA(I,K,J) + 4.12*ALOAD(I,K) * RS2(I,K) *
     1                    (FR1(I,K) * SA20166 + FR2(I,K) * SA205 +
     2                     FR3(I,K) * SA210)
 400  CONTINUE
C
      DO 600 K = 1, LAY
      DO 600 I = IL1, IL2
        ABSA(I,K,9)     =  0.0
 600  CONTINUE
C
      RETURN
      END
