      SUBROUTINE DUSTAERO3 (EXTA, EXOMA, EXOMGA, FA, ABSA, 
     1                      EXTA055, OMA055, GA055,
     2                      EXTB, ODFB, SSAB, ABSB, EXTB055,ODF055,
     3                      SSA055, ABS055,
     4                      ALOAD1, ALOAD2, 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     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION DUSTAERO2 FOR GCM15G/H/I:
C     *                         - SLOAD1,SLOAD2 NOW INTERNAL WORK ARRAYS. 
C     * APR 25,2003 - J.LI. PREVIOUS VERSION DUSTAERO FOR GCM15E/F.
C----------------------------------------------------------------------C
C     CALCULATION OF OPTICAL PROPERTIES FOR DUST 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     EXTA055:EXTINCTION COEFFICIENT AT 0.55 UM                        C
C     OMA055: SINGLE SCATTERING ALBEDO AT 0.55 UM                      C
C     GA055:  ASYMMETRY FACTOR AT 0.55 UM                              C
C     ALOAD1: DUST LOADING FOR FINE MODE                               C
C     ALOAD2: DUST LOADING FOR COARSE MODE                             C
C----------------------------------------------------------------------C
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9)
C
C     * OUTPUT ARRAYS.
C
      REAL   EXTA(ILG,LAY,NBS), EXOMA(ILG,LAY,NBS),
     1       EXOMGA(ILG,LAY,NBS), FA(ILG,LAY,NBS) 
      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),
     1       GA055(ILG,LAY)
      REAL   EXTB055(ILG,LAY), ODF055(ILG,LAY),
     1       SSA055(ILG,LAY), ABS055(ILG,LAY)
      REAL   ABSA(ILG,LAY,NBL)
C
C     * INPUT ARRAYS.
C
      REAL   ALOAD1(ILG,LAY), ALOAD2(ILG,LAY),
     1       ALOAD(ILG,LAY)
C
C     * INTERNAL WORK ARRAYS.
C
      REAL   SLOAD1(ILG,LAY), SLOAD2(ILG,LAY)
      REAL   SE(NBS,2), SW(NBS,2), SG(NBS,2), SA(NBL,2)
C
      REAL   SE055(2),SW055(2),G055(2)
      REAL   SLOAD551(ILG,LAY), SLOAD552(ILG,LAY)
      REAL   FACDU
C
c        DATA ((SE(J,K), J = 1, NBS), K= 1, 2)          /
c     1  0.5472E+00, 0.6015E+00, 0.6026E+00, 0.4516E+00,
c     2  0.7444E-01, 0.7657E-01, 0.7999E-01, 0.8605E-01 /
cC
c        DATA ((SW(J,K), J = 1, NBS), K= 1, 2)          /
c     1  0.8225, 0.9368, 0.9443, 0.8790,
c     2  0.6333, 0.7635, 0.7812, 0.7016                 / ! ????
cC
c        DATA ((SG(J,K), J = 1, NBS), K= 1, 2)          /
c     1  0.7552, 0.6976, 0.6887, 0.6941,
c     2  0.8999, 0.8396, 0.8077, 0.8048                 /
cC
c        DATA ((SA(J,K), J = 1, NBL), K= 1, 2)                    /
c     1  1.396E-02,  2.372E-02,  5.294E-02,  9.146E-02, 1.842E-01,
c     1  1.145E-01,  8.252E-02,  7.634E-02,  4.076E-02,
c     2  1.178E-02,  1.816E-02,  3.068E-02,  3.588E-02, 4.513E-02,
c     2  4.572E-02,  4.220E-02,  4.700E-02,  4.146E-02            /
C
CCCC Jiangnan corrected parameters for band calculation
C
c        DATA ((SE(J,K), J = 1, NBS), K= 1, 2)          /
c     1  1.9827E+00, 1.5467E+00, 5.8710E-01, 6.0982E-02,
c     2  3.5795E-01, 3.8706E-01, 3.7115E-01, 1.6996E-01 /
C
c        DATA ((SW(J,K), J = 1, NBS), K= 1, 2)          /
c     1  9.2771E-01, 9.4909E-01, 9.2426E-01, 6.6900E-01,
c     2  7.8074E-01, 8.5087E-01, 9.0524E-01, 8.7040E-01 /
C
c        DATA ((SG(J,K), J = 1, NBS), K= 1, 2)          /
c     1  6.9309E-01, 6.7709E-01, 6.8139E-01, 6.0267E-01,
c     2  7.9170E-01, 7.4094E-01, 7.8688E-01, 8.7738E-01 /
C
c        DATA ((SA(J,K), J = 1, NBL), K= 1, 2)                     /
c     1  1.5186E-02, 1.5452E-02, 2.7977E-02, 4.2728E-02, 5.7816E-02,
c     1  3.8332E-02, 2.5348E-02, 3.4011E-02, 3.9038E-02,           
c     2  1.8447E-02, 1.8589E-02, 2.8312E-02, 4.9108E-02, 7.9832E-02,
c     2  5.3812E-02, 3.4652E-02, 4.5685E-02, 4.6972E-02            /
C
        DATA ((SE(J,K), J = 1, NBS), K= 1, 2)          /
     1  0.1983E+01, 0.1547E+01, 0.5870E+00, 0.6095E-01,
     2  0.3579E+00, 0.3871E+00, 0.3711E+00, 0.1700E+00 / 
C
        DATA ((SW(J,K), J = 1, NBS), K= 1, 2)          /
     1  0.9277E+00, 0.9491E+00, 0.9243E+00, 0.6690E+00,
     2  0.7807E+00, 0.8509E+00, 0.9052E+00, 0.8704E+00 /
C
        DATA ((SG(J,K), J = 1, NBS), K= 1, 2)          /
     1  0.6931E+00, 0.6771E+00, 0.6814E+00, 0.6027E+00,
     2  0.7917E+00, 0.7409E+00, 0.7869E+00, 0.8774E+00 /
C
        DATA ((SA(J,K), J = 1, NBL), K= 1, 2)                    /
     1 0.1519E-01, 0.1545E-01, 0.2797E-01, 0.4272E-01, 0.5782E-01,
     1 0.3833E-01, 0.2535E-01, 0.3402E-01, 0.3904E-01,
     2 0.1847E-01, 0.1861E-01, 0.2836E-01, 0.4921E-01, 0.8002E-01,
     2 0.5393E-01, 0.3471E-01, 0.4577E-01, 0.4706E-01      /
C
C
        DATA (SE055(K),K = 1, 2)                       /
     1  1.95913, 0.36309                               /
C
        DATA (SW055(K),K = 1, 2)                       /
     1  0.93730, 0.79935                               /
C
        DATA (G055(K),K = 1, 2)                        /
     1  0.68922, 0.77746                               /
C
C----------------------------------------------------------------------C
C     DUST WITH TWO MODES: ACCUMULATION MODE WITH MODEL RADIUS         C
C     0.39 UM, SIGMA 2.0; COARSE MODE WITH MODEL RADIUS 1.9 UM,        C
C     SIGMA 2.15. DRY DENSITY 2.6 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 DUST 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
      FACDU             =  1.0
      DO 100 K = 1, LAY
      DO 100 I = IL1, IL2
        SLOAD1(I,K)     =  10000.0 * ALOAD1(I,K)
        SLOAD2(I,K)     =  10000.0 * ALOAD2(I,K)
        ALOAD(I,K)      =  ALOAD1(I,K) + ALOAD2(I,K)
        SLOAD551(I,K)   =  10000.0 * ALOAD1(I,K) * FACDU
        SLOAD552(I,K)   =  10000.0 * (ALOAD2(I,K) + 
     1                     ALOAD1(I,K)*(1.0-FACDU))
  100 CONTINUE
      DO 200 J = 1, NBS
      DO 200 K = 1, LAY
      DO 200 I = IL1, IL2
        EXT1            =  SLOAD551(I,K) * SE(J,1)
        EXT2            =  SLOAD552(I,K) * SE(J,2)
        EXOM1           =  EXT1 * SW(J,1)
        EXOM2           =  EXT2 * SW(J,2)
        EXOMG1          =  EXOM1 * SG(J,1)
        EXOMG2          =  EXOM2 * SG(J,2)
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 * SG(J,1) + 
     1                     EXOMG2 * SG(J,2)
C
        ABSB(I,K,J)     =  EXT1 * (1.0 - SW(J,1)) +
     1                     EXT2 * (1.0 - SW(J,2))
C
        TMASS           =  10000.0 * ALOAD(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
      DO 250 K = 1, LAY
      DO 250 I = IL1, IL2
        EXT0551         =  SLOAD551(I,K) * SE055(1)
        EXT0552         =  SLOAD552(I,K) * SE055(2) 
        EXTA055(I,K)    =  EXTA055(I,K) + EXT0551 + EXT0552
        OM0551          =  SW055(1) * EXT0551
        OM0552          =  SW055(2) * EXT0552
        OMA055(I,K)     =  OMA055(I,K) + OM0551 + OM0552
        GA0551          =  G055(1) * OM0551
        GA0552          =  G055(2) * OM0552
        GA055(I,K)      =  GA055(I,K) + GA0551*G055(1) + GA0552*G055(2)
C
        ABS055(I,K)     =  EXT0551*(1.0-SW055(1))+
     1                     EXT0552*(1.0-SW055(2))
        TMASS           =  10000.0 * ALOAD(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)    = (OM0551+OM0552) / (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
        ABSA(I,K,J)     =  ABSA(I,K,J) + SLOAD551(I,K) * SA(J,1) +
     1                     SLOAD552(I,K) * SA(J,2)
 400  CONTINUE
C
      RETURN
      END
