      SUBROUTINE SSALTAERO5(EXTA,  EXOMA, EXOMGA, FA, ABSA, 
     1                      EXTA055, OMA055, GA055,
     2                      EXTB, ODFB, SSAB, ABSB, EXTB055, ODF055, 
     3                      SSA055, ABS055, RH, ALOAD1,  ALOAD2, 
     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, TWO MODES ARE 0.732 AND
C     *                             6.13 UM AS THE WET SIZE, PLEASE NOTE
C     *                             THE BAND CALCULTION ARE BASED ON DRY
C     *                             RECALCUTE THE 4 BAND SHORTWAVE PART
C     *                             TO MAKE IT IN THE SAME SIZE AS 0.55
C     *                             KEEP THE LONGWAVE PART AS BEFORE (GLEN).
C     * FEB 13,2009 - M.LAZARE. PREVIOUS VERSION SSALTAERO4 FOR GCM15H/I:
C     *                         - RHMAX LOWERED FROM 0.98 TO 0.95. 
C     *                         - USE BOUNDED RH (SRH) INSTEAD OF RH TO
C     *                           CALCULATE GROWTH FACTOR X.
C     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION SSALTAERO3 FOR GCM15G:
C     *                         - SRH,FRH,SLOAD1,SLOAD2 NOW INTERNAL 
C     *                           WORK ARRAYS.
C     * JUN 19/06 -  M.LAZARE. PREVIOUS VERSION SSALTAERO2 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     * APR 25,2003 - J.LI.    PREVIOUS VERSION SSALTAERO FOR UP TO
C     *                        AND INCLUDING GCM15E.  
C----------------------------------------------------------------------C
C     CALCULATION OF OPTICAL PROPERTIES FOR SEA SALT WITH TWO MODES    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     SLOAD1: SEA SALT LOADING FOR FINE MODE                           C
C     SLOAD2: SEA SALT LOADING FOR COARSE MODE                         C
C----------------------------------------------------------------------C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9)
C
C     * OUTPUT ARRAYS.
C
      REAL, DIMENSION(ILG,LAY,NBS):: EXTA, EXOMA, EXOMGA, FA, 
     1                               ABSB, ODFB, EXTB, SSAB
      REAL ABSA(ILG,LAY,NBL)
C
      REAL, DIMENSION(ILG,LAY):: EXTA055, OMA055, GA055, EXTB055, 
     1                           ODF055, SSA055, ABS055
C
C     * INPUT ARRAYS.
C
      REAL, DIMENSION(ILG,LAY):: RH, ALOAD1, ALOAD2
C
C     * INTERNAL WORK ARRAYS.
C
      REAL, DIMENSION(ILG,LAY):: SRH, FRH, SLOAD1, SLOAD2
      REAL, DIMENSION(4,NBS,2):: SE, SW, SG
      REAL SA(4,NBL,2)
C
      REAL   SE055(4,2),CW055(2,2),G055(4,2)
      REAL   SLOAD551(ILG,LAY), SLOAD552(ILG,LAY)
      REAL   FACSS
C
CCCC Jiangnan corrected parameters for band calculation
C
        DATA (((SE(I,J,K), I = 1, 4), J = 1, NBS), K= 1, 2) /
     1  1.391E+00,-1.545E-01, 4.648E-02, 1.281E+00,
     1  1.309E+00, 5.348E-02, 5.152E-02, 1.226E+00,
     1  8.957E-01, 3.558E-01, 4.128E-02, 9.062E-01,
     1  4.099E-01, 3.625E-01, 1.496E-02, 4.504E-01,
     2  1.317E-01,-7.488E-03, 4.379E-03, 1.218E-01,
     2  1.362E-01,-9.357E-03, 4.498E-03, 1.262E-01,
     2  1.434E-01,-1.308E-02, 4.664E-03, 1.330E-01,
     2  1.554E-01,-1.798E-02, 5.028E-03, 1.330E-01           /
C
        DATA (((SW(I,J,K), I = 1, 4), J = 1, NBS), K= 1, 2) /
     1  3.938E-06,-2.809E-06, 5.888E-09, 8.927E-06,
     1  3.236E-04,-2.417E-04,-6.789E-07, 7.616E-04,
     1  3.045E-03,-1.918E-03,-1.590E-04, 4.597E-03,
     1  2.221E-01,-6.472E-03, 1.382E-03, 3.493E-02,
     2  3.277E-05,-2.335E-05, 9.002E-08, 7.038E-05,
     2  2.869E-03,-2.044E-03, 4.838E-06, 6.325E-03,
     2  1.950E-02,-8.772E-03,-1.112E-03, 2.904E-02,
     2  2.218E-01, 3.505E-02,-4.841E-03, 9.729E-02          /
C
        DATA (((SG(I,J,K), I = 1, 4), J = 1, NBS), K= 1, 2) /
     1  7.667E-01, 1.750E-02,-2.244E-03, 6.964E-01,
     1  7.764E-01, 2.025E-02,-3.875E-04, 7.049E-01,
     1  7.682E-01, 4.617E-02, 2.603E-04, 7.028E-01,
     1  7.245E-01, 7.342E-02,-3.145E-03, 6.550E-01,
     2  8.102E-01, 7.596E-02, 1.021E-02, 7.908E-01,
     2  8.030E-01, 7.067E-02, 3.140E-03, 7.784E-01,
     2  7.973E-01, 4.625E-02,-8.634E-04, 7.613E-01,
     2  8.404E-01, 3.119E-02,-2.783E-03, 7.470E-01          /
C
        DATA (((SA(I,J,K), I = 1, 4), J = 1, NBL), K= 1, 2) /
     1  1.300E-02, 2.315E-02,-3.508E-04, 2.517E-03,
     1  1.322E-02, 2.226E-02,-4.119E-04, 3.647E-03,
     1  4.858E-02, 7.799E-02, 7.454E-04, 1.185E-02,
     1  2.469E-02, 2.915E-02,-4.593E-04, 1.497E-02,
     1  2.434E-02, 3.196E-02,-3.737E-04, 1.115E-02,
     1  4.888E-02, 7.948E-02, 1.013E-03, 6.968E-03,
     1  9.355E-02, 1.607E-01, 3.343E-03, 1.769E-02,
     1  6.312E-02, 9.555E-02, 8.945E-05, 3.203E-02,
     1  2.530E-02, 2.432E-02,-4.512E-04, 1.908E-02,
     2  1.621E-02, 1.378E-02, 6.919E-04, 3.440E-03,
     2  1.691E-02, 1.455E-02, 7.138E-04, 5.132E-03,
     2  4.295E-02, 1.331E-02, 1.832E-03, 1.481E-02,
     2  2.903E-02, 1.505E-02, 1.161E-03, 2.070E-02,
     2  2.939E-02, 1.525E-02, 1.167E-03, 1.781E-02,
     2  4.336E-02, 1.393E-02, 1.806E-03, 1.142E-02,
     2  6.891E-02, 7.288E-03, 2.792E-03, 2.550E-02,
     2  6.488E-02, 1.362E-02, 2.775E-03, 4.389E-02,
     2  2.800E-02, 9.300E-03, 1.212E-03, 2.329E-02           /

C
        DATA ((SE055(I,K), I = 1, 4), K= 1, 2) /
     1  1.415E+00,-1.692E-01, 4.704E-02, 1.30086, 
     1  1.320E-01,-7.585E-03, 4.385E-03, 0.12199             /
C
        DATA ((CW055(I,K), I = 1, 2), K= 1, 2) /
     1  0.99999,  0.99999, 
     1  0.99999,  0.99999                                    /
C
        DATA ((G055(I,K), I = 1, 4), K= 1, 2)  /
     1  7.646E-01, 4.447E-02,-1.197E-03, 0.69225,
     1  7.731E-01,-6.938E-02, 1.272E-02, 0.79106             /
C
      DATA RHMAX /0.95/
C
C----------------------------------------------------------------------C
C     SEA SALT WITH TWO MODES: ACCUMULATION MODE WITH MODEL RADIUS     C
C     0.209 UM, SIGMA 2.03; COARSE MODE WITH MODEL RADIUS 1.75 UM,     C
C     SIGMA 2.03. DENSITY 2.24 KG/M^3                                  C
C                                                                      C
C     SOLAR, SE: SPECIFIC EXTINCTION, OMA: SINGLE SCATTERING ALBEDO,   C
C     GA: ASYMMETRY FACTOR, EXTA: EXTINCTION COEFFICIENT               C
C     FACTOR 10000 BECAUSE THE UNIT OF SPECIFIC EXTINCTION FOR AEROSOL C
C     IS M^2/GRAM, WHILE FOR GAS IS CM^2/GRAM, IN RADDRIV WE USE SAME  C
C     DP (AIR DENSITY * LAYER THICKNESS) IS USED FOR BOTH GA AND       C
C     AEROSOL. ALOAD S SALT AEROSOL LOADING IN UNIT G (AEROSOL) / G    C
C     (AIR). TAU = ALOAD * SPECIFIC EXTINCTION * AIR DENSITY *         C
C     LAYER THICKNESS                                                  C
C     DATA BY G. LENSINS, CODE BY J. LI (2002, 2)                      C
C----------------------------------------------------------------------C
C
      FACSS              =  0.90
      DO 100 K = 1, LAY
      DO 100 I = IL1, IL2
        SRH(I,K)         =  MIN (RH(I,K), RHMAX)
        FRH(I,K)         =  1.0 / (SRH(I,K) - 1.05)
        X                =  10000.0 * EXP(2.9115E-03 + 1.5673 * 
     1                      SRH(I,K) + 1.1818E-02 / ((SRH(I,K) - 1.06) *
     2                      (SRH(I,K) - 1.06)))
        SLOAD1(I,K)      =  X * ALOAD1(I,K)
        SLOAD2(I,K)      =  X * ALOAD2(I,K)
        SLOAD551(I,K)    =  X * ALOAD1(I,K) * FACSS
        SLOAD552(I,K)    =  X * (ALOAD2(I,K) + ALOAD1(I,K)*(1.0-FACSS)) 
  100 CONTINUE
C
      DO 200 J = 1, NBS
      DO 200 K = 1, LAY
      DO 200 I = IL1, IL2
        IF (RH(I,K) .GE. 0.45)                                      THEN
          EXT1           =  SLOAD551(I,K) * (SE(1,J,1) + SE(2,J,1) * 
     1                      SRH(I,K) + SE(3,J,1) * FRH(I,K))
          EXT2           =  SLOAD552(I,K) * (SE(1,J,2) + SE(2,J,2) * 
     1                      SRH(I,K) + SE(3,J,2) * FRH(I,K))
C
          OM1            =  1.0 - (SW(1,J,1) + SW(2,J,1) * SRH(I,K) + 
     1                      SW(3,J,1) * FRH(I,K))
          OM2            =  1.0 - (SW(1,J,2) + SW(2,J,2) * SRH(I,K) +
     1                      SW(3,J,2) * FRH(I,K))
          EXOM1          =  EXT1 * OM1
          EXOM2          =  EXT2 * OM2
C
          G1             =  SG(1,J,1) + SG(2,J,1) * SRH(I,K) + 
     1                      SG(3,J,1) * FRH(I,K)
          G2             =  SG(1,J,2) + SG(2,J,2) * SRH(I,K) + 
     1                      SG(3,J,2) * FRH(I,K)
          EXOMG1         =  EXOM1 * G1
          EXOMG2         =  EXOM2 * G2 
C 
          ABSB(I,K,J)    =  EXT1 * (1.0 - OM1) + EXT2 * (1.0 - OM2)   
C
        ELSE
          EXT1           =  SLOAD551(I,K) * SE(4,J,1)
          EXT2           =  SLOAD552(I,K) * SE(4,J,2)
          EXOM1          =  EXT1 * (1.0 - SW(4,J,1))
          EXOM2          =  EXT2 * (1.0 - SW(4,J,2))
          G1             =  SG(4,J,1)
          G2             =  SG(4,J,2)
          EXOMG1         =  EXOM1 * G1
          EXOMG2         =  EXOM2 * G2

          ABSB(I,K,J)    =  EXT1 * SW(4,J,1) + EXT2 * SW(4,J,2)
        ENDIF
C
        EXTA(I,K,J)      =  EXTA(I,K,J) + EXT1 + EXT2
        EXOMA(I,K,J)     =  EXOMA(I,K,J) + EXOM1 + EXOM2
        EXOMGA(I,K,J)    =  EXOMGA(I,K,J) + EXOMG1 + EXOMG2
        FA(I,K,J)        =  FA(I,K,J) + EXOMG1 * G1 + EXOMG2 * G2
C
        TMASS            =  10000.0 * (ALOAD1(I,K) + ALOAD2(I,K))
        IF (TMASS.GT.1.E-15) THEN
         EXTB(I,K,J)     =  (EXT1 + EXT2) / TMASS
        ELSE
         EXTB(I,K,J)     =  0.0
        ENDIF
        IF (ALOAD1(I,K).GT.1.E-19) THEN
         ODFB(I,K,J)     =  EXT1 / (10000.0 * ALOAD1(I,K))
        ELSE
         ODFB(I,K,J)     =  0.0
        ENDIF
C
        IF ((EXT1 + EXT2) .GT. 1.E-20) THEN
         SSAB(I,K,J)     =  (EXOM1 + EXOM2) / (EXT1 + EXT2)
        ELSE
         SSAB(I,K,J)     =  0.0
        ENDIF
C
  200  CONTINUE
C
      DO 250 K = 1, LAY
      DO 250 I = IL1, IL2
        IF (RH(I,K) .GE. 0.45) THEN
          EXT0551        =  SLOAD551(I,K) * (SE055(1,1) + SE055(2,1) *
     1                      SRH(I,K) + SE055(3,1) * FRH(I,K))
          EXT0552        =  SLOAD552(I,K) * (SE055(1,2) + SE055(2,2) *
     1                      SRH(I,K) + SE055(3,2) * FRH(I,K))
C
          OM0551         =  CW055(1,1)
          OM0552         =  CW055(1,2)
C
          GA0551         =  G055(1,1) + G055(2,1) * SRH(I,K) +
     1                      G055(3,1) * FRH(I,K)
          GA0552         =  G055(1,2) + G055(2,2) * SRH(I,K) +
     1                      G055(3,2) * FRH(I,K)
        ELSE
          EXT0551        =  SLOAD551(I,K) * SE055(4,1)
          EXT0552        =  SLOAD552(I,K) * SE055(4,2)
          OM0551         =  CW055(2,1)
          OM0552         =  CW055(2,2)
          GA0551         =  G055(4,1)
          GA0552         =  G055(4,2)
        ENDIF
C
        EXTA055(I,K)     =  EXTA055(I,K) + EXT0551 + EXT0552
        OMASUM           =  EXT0551 * OM0551 + EXT0552 * OM0552
        OMA055(I,K)      =  OMA055(I,K) + OMASUM
        GA055(I,K)       =  GA055(I,K) + (EXT0551 * OM0551 * GA0551 * 
     1                      GA0551 + EXT0552 * OM0552 * GA0552 * GA0552)
C
        ABS055(I,K)      =  EXT0551*(1.0 - OM0551) + 
     1                      EXT0552*(1.0 - OM0552)
        TMASS            =  10000.0 * (ALOAD1(I,K) + ALOAD2(I,K))
        IF (TMASS.GT.1.E-15) THEN
         EXTB055(I,K)    =  (EXT0551 + EXT0552) / TMASS
        ELSE
         EXTB055(I,K)    =  0.0
        ENDIF
        IF (ALOAD1(I,K).GT.1.E-19) THEN
         ODF055(I,K)     =  EXT0551 / (10000.0 * ALOAD1(I,K))
        ELSE
         ODF055(I,K)     =  0.0
        ENDIF
C
        IF ((EXT0551+EXT0552).GT.1.E-20) THEN
         SSA055(I,K)     =  OMASUM / (EXT0551 + EXT0552)
        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
C
        IF (RH(I,K) .GE. 0.45)                                      THEN
          ABS1           =  SLOAD551(I,K) * (SA(1,J,1) + SA(2,J,1) *
     1                      SRH(I,K) + SA(3,J,1) * FRH(I,K))
          ABS2           =  SLOAD552(I,K) * (SA(1,J,2) + SA(2,J,2) *
     1                      SRH(I,K) + SA(3,J,2) * FRH(I,K))
        ELSE
          ABS1           =  SLOAD551(I,K) * SA(4,J,1)
          ABS2           =  SLOAD552(I,K) * SA(4,J,2)
        ENDIF
C
        ABSA(I,K,J)      =  ABSA(I,K,J) + ABS1 + ABS2
 400  CONTINUE
C
      RETURN
      END
