      SUBROUTINE DUSTEMI3
     1    (FLNDROW,GTROW,PRESSG,DSHJ,SMFRAC,FNROL,BSFRAC,
     2    ZSPDSBS,USTARBS,SUZ0ROW,PDSFROW, 
     3    SPOTROW,ST01ROW,ST02ROW,ST03ROW,ST04ROW,ST06ROW,
     4    ST13ROW,ST14ROW,ST15ROW,ST16ROW,ST17ROW,
     5    ESDROW,XEMIS,IDUA,IDUC,SAVERAD,ISVDUST,
     6    UTH,SREL,SRELV,SU_SRELV,FALLROW,FA10ROW,FA2ROW,FA1ROW,
     7    DUWDROW,DUSTROW,DUTHROW,USMKROW,
     8    ILG,IL1,IL2,ILEV,LEV,NTRAC,ZTMST,CUSCALE)
C
C     * THIS ROUTINE CALCULATES THE SURFACE DUST EMISSION FLUX 
C     * IN ACCUMULATION AND COARSE MODES  
C 
C     * SEP 13/2014 - K. VONSALZEN/ NEW VERSION FOR GCM18+:
C     *               M. LAZARE:    - USTARBS PASSED IN AND USED
C     *                               TO DEFINE "USTAR", INSTEAD OF
C     *                               PREVIOUS DEFINITION.
C     *                             - PASSES IN AND USES FLNDROW
C     *                               INSTEAD OF {GCROW,MASKROW}.
C     *                             - BARE SOIL SURFACE WIND SPEED
C     *                               (WITH GUSTINESS EFFECTS) PASSED IN
C     *                               (ZSPDSBS) INSTEAD OF CALCULATING
C     *                               VIA {SFCUROW,SFCVROW,GUSTROL}
C     *                               WHICH ARE NOW REMOVED.
C     * JUN 18/2013 - K. VONSALZEN/ PREVIOUS VERSION DUSTEMI2 FOR GCM17+:
C     *               M. LAZARE:    - USES "DUPARM1" INSTED OF "DUPARM",
C     *                               WHICH HAS "CUSCALE" MOVED ELSEWHERE
C     *                               IN THE CODE ALLOWING "DUPARM1" TO
C     *                               BE USED BY BOTH BULK AND PLA AEROSOLS.
C     *                               "CUSCALE" IS NOW PASSED IN THE CALL.                 
C     * AUG 20/2007 - Y. PENG.      PREVIOUS VERSION DUSTEMI FOR GCM16:
C     *                             NEW SCHEME INCL. SOIL TEXTURE:
C     *                             - USE ZOBLER SOIL PROPERTIES 
C     *                             - USE MARTICORENA ET AL 1997 TO 
C     *                               CALCULATE THE THRESHOLD WIND SPEED
C     *                             - WORK ARRAYS NOW LOCAL.
C     *
C
       USE DUPARM1 
C
       IMPLICIT REAL (A-H,O-Z),INTEGER (I-N)
C 
C      * I/O FIELDS.
C 
       REAL FLNDROW(ILG),GTROW(ILG),SMFRAC(ILG)
       REAL PRESSG(ILG),DSHJ(ILG,ILEV)
       REAL FNROL(ILG),BSFRAC(ILG),ZSPDSBS(ILG),USTARBS(ILG)
       REAL SUZ0ROW(ILG),PDSFROW(ILG)
       REAL SPOTROW(ILG)
       REAL ST01ROW(ILG),ST02ROW(ILG)
       REAL ST03ROW(ILG),ST04ROW(ILG),ST06ROW(ILG)
       REAL ST13ROW(ILG),ST14ROW(ILG)
       REAL ST15ROW(ILG),ST16ROW(ILG),ST17ROW(ILG)
       REAL ESDROW(ILG),XEMIS(ILG,ILEV,NTRAC)   
C
       REAL FALLROW(ILG),FA10ROW(ILG)
       REAL FA2ROW(ILG),FA1ROW(ILG)
       REAL DUWDROW(ILG),USMKROW(ILG)
       REAL DUSTROW(ILG),DUTHROW(ILG)
C
C      * INTERNAL WORK FIELDS.
C
        REAL UTH(NCLASS),SREL(NATS,NCLASS)
        REAL SRELV(NATS,NCLASS)
        REAL SU_SRELV(NATS,NCLASS)
C
        REAL C_EFF(ILG)
        REAL DPK(NTRACE),DBMIN(NTRACE),DBMAX(NTRACE)
        REAL BAREFRAC(ILG),DUST_MASK(ILG),USTAR_D(ILG)
        REAL DU_WIND(ILG),DU_TH(ILG)
        REAL USTAR_ACRIT(ILG)
        REAL DUZ01(ILG),DUZ02(ILG)
C
        REAL FLUX_6H(ILG,NTRACE)
        REAL FLUXBIN(ILG,NTRACE)
        REAL FLUXALL(ILG),FLUXA10(ILG)
        REAL FLUXA2(ILG),FLUXA1(ILG)
        REAL FLUX_AI(ILG),FLUX_CI(ILG)
        REAL FLUXTYP(ILG,NCLASS)
        REAL FLUXDIAM1(ILG,NCLASS)
        REAL FLUXDIAM2(ILG,NCLASS)
        REAL FLUXDIAM3(ILG,NCLASS)
        REAL FLUXDIAM4(ILG,NCLASS)
        REAL FLUXDIAM6(ILG,NCLASS)
        REAL FLUXDIAM_PF(ILG,NCLASS)
        REAL FLUXDIAM13(ILG,NCLASS)
        REAL FLUXDIAM14(ILG,NCLASS)
        REAL FLUXDIAM15(ILG,NCLASS)
        REAL FLUXDIAM16(ILG,NCLASS)
        REAL FLUXDIAM17(ILG,NCLASS)
C
C      * AUXILIARY VARIABLES
C
        REAL AAA, BB, CCC, DDD, EE, FF
        REAL D1, FEFF
        REAL RDP, DLAST
        REAL USTAR,UTHP(ILG,NATS)
        REAL ALPHA
        REAL WAA(NATS),WBB(ILG)
        REAL FDP1, FDP2
C
        INTEGER IL,J, KK, NN, I_SOIL
        INTEGER KKMIN,KKK
        INTEGER DK(NCLASS-1)
C
       COMMON /PARAMS/ WW,TW,RAYON,ASQ,GRAV,RGAS,RGOCP,RGOASQ,CPRES        
C
        INTEGER ZSCHM,WSCHM,SOLSCHM
c        z0 scheme, 1: new scheme Prigent et al 2002, 0: old scheme z01=z02=0.001
        DATA ZSCHM /1/
c        soil water scheme, 1: new scheme Fecan, 0: old scheme 0.99
        DATA WSCHM /1/
c        soil property scheme, 1: new scheme Cheng, 0: old scheme 
        DATA SOLSCHM /1/
C
C    --------------------------------------------------------------
C
C     * INITIAL CALCULATION (C_EFF)
C
C      ! EFFECTIVE FRACTION CALCULATION (ROUGHNESS Z0 RELATED)
C      ! CONSTANT Z0 TEMPORARILY
C      ! MARTICORENA ET AL 1997 EQUATIONS (15)-(17)
C
C      ! initial values !D1 is set to zero temporarily
C
       D1 = 0.
       FEFF = 0.
       AAA = 0.
       BB = 0.
       CCC = 0.
       DDD = 0.
       EE = 0.
       FF = 0.
       DO IL=IL1,IL2
        C_EFF(IL)       = 0.0
        DUZ01(IL)       = 0.0
        DUZ02(IL)       = 0.0
       END DO
C
       DO IL=IL1,IL2
C
       IF (ZSCHM .EQ. 1) THEN
C       new Z0 scheme, monthly average
        DUZ01(IL)=SUZ0ROW(IL)
        DUZ02(IL)=SUZ0ROW(IL)  
       ELSE
C       old z0 scheme
        DUZ01(IL)=Z01
        DUZ02(IL)=Z02
       ENDIF
C  
       IF (DUZ01(IL).EQ.0..OR.DUZ01(IL).LT.Z0S.OR.DUZ02(IL).LT.Z0S) THEN
        FEFF = 0.
       ELSE
        AAA = LOG(DUZ01(IL)/Z0S)
        BB = LOG(AEFF*(XEFF/Z0S)**0.8)
        CCC = 1.- AAA/BB
         IF (D1 .EQ. 0.) THEN
          FF = 1.
         ELSE
          DDD = LOG(DUZ02(IL)/DUZ01(IL))
          EE = LOG(AEFF*(D1/DUZ01(IL))**0.8)
          FF = 1.-DDD/EE
         ENDIF  ! D1=0
        FEFF = FF*CCC
        IF (FEFF .LT. 0.) FEFF = 0.
        IF (FEFF .GT. 1.) FEFF = 1.
       ENDIF  ! Z01=0
C
       C_EFF(IL) = FEFF
C       
       END DO  !IL 
C
C     * INITIAL CALCULATION (C_EFF) DONE
C
C   --------------------------------------------------------------
C
C     * INITIAL CALCULATION (WINDS and MASKS)
C
C      ! initial values
C 
        USTAR=0.0
        RDP=0.0
        DLAST=0.0
       DO IL=IL1,IL2 
        UTHP(IL,:)      = 0.0
        DUST_MASK(IL)   = 0.0
        USTAR_D(IL)     = 0.0
        DU_WIND(IL)     = 0.0
        DU_TH(IL)       = 0.0
        USTAR_ACRIT(IL) = 0.0
C
C      BSFRAC IS THE BARE SOIL FRACTION PASSED FROM CCCMA GCM LAND SCHEME CLASS (DEFAULT). 
C      BSFRAC IS COMBINED WITH THE SOIL MOISTURE CRITERION (W0) TO OBTAIN THE POTENTIAL DUST 
C      SOURCE AREA.  
C      PDSFROW IS A PRESCRIBED POTENTIAL DUST SOURCE AREAL FRACTION OBTAINED FROM AN OFFLINE 
C      TERRESTRIAL MODEL, IN WHICH THE SOIL MOISTURE & SNOW COVER CRITERIA WERE SATISFIED.
C      PDSFROW IS ONLY APPLIED FOR THE PURPOSE OF DUST PARAMETERIZATION DEVELOPMENT.
C   
        IF (CCCMABF) THEN
         BAREFRAC(IL)=BSFRAC(IL)
        ELSE
         BAREFRAC(IL)=PDSFROW(IL)
        ENDIF
C
       END DO
C
C      ! CALCULATE BIN MINIMUM/MAXIMUM RADIUS IN UM
C      ! RDP IS DIAMETER IN CM 
C
       RDP=DMIN
       DLAST=DMIN
       DO NN=1,NTRACE
        DPK(NN)=0.
        DBMIN(NN)=0.
        DBMAX(NN)=0.
       END DO   
C
       NN=1
       DO KK=1,NCLASS
        IF (MOD(KK,NBIN).EQ.0) THEN
         DBMAX(NN)=RDP*5000.
         DBMIN(NN)=DLAST*5000.
         DPK(NN)=SQRT(DBMAX(NN)*DBMIN(NN))
         NN=NN+1
         DLAST=RDP
        ENDIF 
        RDP=RDP*EXP(DSTEP)
       END DO
C
C      ! CALCULATE THE FRICTION WIND SPEED USTAR
C
       DO IL=IL1,IL2
         DU_WIND(IL)=ZSPDSBS(IL)
         IF (DUZ02(IL) .EQ. 0.) THEN
          USTAR=0.
         ELSE
          USTAR=100.*USTARBS(IL)
         ENDIF
         USTAR_D(IL)=USTAR
C         
C         IF (C_EFF(IL) .EQ. 0.) THEN
          DU_TH(IL)=0.
C         ELSE
C          DU_TH(IL)=UMIN*CUSCALE/C_EFF(IL)
C         ENDIF
C
C      ! check critical wind speed (wind stress threshold scaled)
          IF (USTAR.GT.0..AND.USTAR.GE.DU_TH(IL)) THEN
C      ! check if the grid cell is a potential dust source (exclude vegetated area)
           IF (BAREFRAC(IL).GT.0.0.AND.DUZ02(IL).GT.0.0) THEN
C      ! set grid point as a potential dust source point
            DUST_MASK(IL)=1.
           ENDIF  ! bare fraction
C      ! store the time fraction with ustar > dust threshold for output 
           USTAR_ACRIT(IL)=1.
          ELSE
           USTAR_ACRIT(IL)=0.
          ENDIF   ! USTAR
C
       END DO    ! IL
C
C     * INITIAL CALCULATION (WINDS and MASKS) DONE
C
C   --------------------------------------------------------------
C
C     * EMISSION FLUX CALCULATION
C
C      ! initialization
C
        ALPHA=0.0
        FDP1=0.0
        FDP2=0.0
       DO J=1,NATS
        WAA(J)=0.0
       ENDDO
       DO IL=IL1,IL2
        WBB(IL)=0.0
C
        DO NN=1,NTRACE
         FLUXBIN(IL,NN)=0.
         FLUX_6H(IL,NN)=0.
        END DO
C
        DO KK=1,NCLASS
         FLUXTYP(IL,KK)=0.
         FLUXDIAM1(IL,KK)=0.
         FLUXDIAM2(IL,KK)=0.
         FLUXDIAM3(IL,KK)=0.
         FLUXDIAM4(IL,KK)=0.
         FLUXDIAM6(IL,KK)=0.
         FLUXDIAM_PF(IL,KK)=0.
         FLUXDIAM13(IL,KK)=0.
         FLUXDIAM14(IL,KK)=0.
         FLUXDIAM15(IL,KK)=0.
         FLUXDIAM16(IL,KK)=0.
         FLUXDIAM17(IL,KK)=0.
        END DO
C
       END DO ! IL
C
C      ! calculation
C
       CD=1.*ROA/GRAVI
C
       DO KK=1,NCLASS
        DO IL=IL1,IL2
         IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
C
C   ! soil moisture modification      
C
          IF (WSCHM .EQ. 1) THEN
C          new ws scheme
            IF (FLNDROW(IL).GT.0. .AND. GTROW(IL).GT.273.15) THEN
             WBB(IL)=MIN(SMFRAC(IL)/ROP,1.)*100.           
            ELSE
             WBB(IL)=0.
            ENDIF
            DO J=1,NATS
             WAA(J)=0.0014*(SOLSPE(NSPE-2,J)*100.)**2+0.17*
     1         (SOLSPE(NSPE-2,J)*100.)
             IF (WBB(IL).LE.WAA(J).OR.WAA(J).EQ.0.) THEN
              UTHP(IL,J)=UTH(KK)
             ELSE
              UTHP(IL,J)=UTH(KK)*SQRT(1.+1.21*(WBB(IL)-WAA(J))**0.68)
             ENDIF
            ENDDO !J 
          ELSE !soil moisture
C          old ws scheme
            DO J=1,NATS
             UTHP(IL,J)=UTH(KK) 
            ENDDO !J    
          ENDIF !soil moisture
C
C    ! fluxdiam calculation (FDP1,FDP2 follows Marticorena)
C
C          ! flux for soil type #1
            I_SOIL=1  
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM1(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #2
            I_SOIL=2
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM2(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #3
            I_SOIL=3
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM3(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #4
            I_SOIL=4 
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM4(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #6
            I_SOIL=6
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM6(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type in preferential source area
            I_SOIL=10
             FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
             FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM_PF(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! add Asian soil types
           IF (SOLSCHM.EQ.1) THEN
C          ! flux for soil type #13
            I_SOIL=13
            FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
            FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM13(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #14
            I_SOIL=14
            FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
            FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM14(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #15
            I_SOIL=15
            FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
            FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM15(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #16
            I_SOIL=16
            FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
            FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM16(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
C          ! flux for soil type #17
            I_SOIL=17
            FDP1=(1.-(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))
            FDP2=(1.+(UTHP(IL,I_SOIL)/(C_EFF(IL)*USTAR_D(IL))))**2
            IF (FDP1 .GT. 0.) THEN
             ALPHA=SOLSPE(NMODE*3+1,I_SOIL)
             FLUXDIAM17(IL,KK)=MAX(0.,SREL(I_SOIL,KK)
     1       *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
C
           ENDIF !Asian soil
C
C    ! initialize fluxtyp of all soil types at first bin (kk=1)
C
            IF (KK.EQ.1) THEN
             IF (SOLSCHM.EQ.1) THEN
             FLUXTYP(IL,KK)=FLUXTYP(IL,KK)
     1                 +FLUXDIAM1(IL,KK)*ST01ROW(IL)
     1                 +FLUXDIAM2(IL,KK)*ST02ROW(IL)
     1                 +FLUXDIAM3(IL,KK)*ST03ROW(IL)
     1                 +FLUXDIAM4(IL,KK)*ST04ROW(IL)
     1                 +FLUXDIAM6(IL,KK)*ST06ROW(IL)
     1                 +FLUXDIAM13(IL,KK)*ST13ROW(IL)
     1                 +FLUXDIAM14(IL,KK)*ST14ROW(IL)
     1                 +FLUXDIAM15(IL,KK)*ST15ROW(IL)
     1                 +FLUXDIAM16(IL,KK)*ST16ROW(IL)
     1                 +FLUXDIAM17(IL,KK)*ST17ROW(IL)
     1                 +FLUXDIAM_PF(IL,KK)*SPOTROW(IL)
             ELSE
             FLUXTYP(IL,KK)=FLUXTYP(IL,KK)
     1                 +FLUXDIAM1(IL,KK)*ST01ROW(IL)
     1                 +FLUXDIAM2(IL,KK)*ST02ROW(IL)
     1                 +FLUXDIAM3(IL,KK)*ST03ROW(IL)
     1                 +FLUXDIAM4(IL,KK)*ST04ROW(IL)
     1                 +FLUXDIAM6(IL,KK)*ST06ROW(IL)
     1                 +FLUXDIAM_PF(IL,KK)*SPOTROW(IL)
             ENDIF !Asian soil
            ELSE 
             DK(KK-1)=KK          ! store dislocated particle class
            ENDIF  ! KK=1
C
         ENDIF     ! DUST_MASK=1 and C_EFF>1
        ENDDO      ! IL
       ENDDO       ! KK
C
C      ! calculate fluxtyp of all soil types for all 192 bins (kkk=2~192)
C
       KKMIN=1
       DO IL=IL1,IL2
         IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
          DO N=1,NCLASS-1
           DO KKK=1,DK(N)           ! scaling with relative contribution of dust size fraction
             I_SOIL=1
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST01ROW(IL)
     1   *FLUXDIAM1(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1   ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN))) 
C 
             I_SOIL=2
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST02ROW(IL) 
     1    *FLUXDIAM2(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=3
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST03ROW(IL) 
     1    *FLUXDIAM3(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=4
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST04ROW(IL) 
     1    *FLUXDIAM4(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=6
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST06ROW(IL) 
     1    *FLUXDIAM6(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
            IF (DU_WIND(IL).GT. 10.)  I_SOIL=11
            IF (DU_WIND(IL).LE. 10.)  I_SOIL=10
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+SPOTROW(IL)
     1    *FLUXDIAM_PF(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
          IF (SOLSCHM.EQ.1) THEN
             I_SOIL=13
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST13ROW(IL) 
     1    *FLUXDIAM13(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=14
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST14ROW(IL) 
     1    *FLUXDIAM14(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=15
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST15ROW(IL) 
     1    *FLUXDIAM15(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=16
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST16ROW(IL) 
     1    *FLUXDIAM16(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
             I_SOIL=17
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+ST17ROW(IL) 
     1    *FLUXDIAM17(IL,DK(N))*SRELV(I_SOIL,KKK)/
     1    ((SU_SRELV(I_SOIL,DK(N))-SU_SRELV(I_SOIL,KKMIN)))
C
          ENDIF !Asian soil       
C
           ENDDO  ! KKK
          ENDDO   ! N 
         ENDIF    ! DUST_MASK
       ENDDO      ! IL
C
C      ! CALCULATE FLUXES IN EACH INTERNAL TRACER CLASSES (8)
C
       DO NN=1,NTRACE
        DO IL=IL1,IL2
         IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
           DO KK=(NN-1)*NBIN+1,NN*NBIN
            FLUXBIN(IL,NN)=FLUXBIN(IL,NN)+FLUXTYP(IL,KK)
           ENDDO
           ! mask out dust fluxes, where soil moisture threshold is reached
           IF (SMFRAC(IL).GT.W0) FLUXBIN(IL,NN)=0.
           ! fluxbin: g/cm2/sec; flux_6h: g/m2/sec
           ! exclude snow covered region and weighted by vegetation cover
           IF (FNROL(IL).GT.0.1) FLUXBIN(IL,NN)=0.
           FLUX_6H(IL,NN)=FLUXBIN(IL,NN)*10000.*BAREFRAC(IL)
         ENDIF ! DUST_MASK
        ENDDO  ! IL 
       ENDDO   ! NN
C
C     * EMISSION FLUX CALCULATION DONE
C  
C   --------------------------------------------------------------
C   
C     * DIAGNOSTIC CALCULATION FOR MODEL OUTPUT
C
C      ! initialization for output arrays
C
       DO IL=IL1,IL2
        FLUXALL(IL)=0.
        FLUXA10(IL)=0.
        FLUXA2(IL)=0.
        FLUXA1(IL)=0.
        FLUX_AI(IL)=0.
        FLUX_CI(IL)=0.
       ENDDO
C
       DO IL=IL1,IL2
       IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
C
C       ! unit of emission flux: FLUX_6H [g/m2/sec]
C
        DO NN=1,NTRACE
         FLUXALL(IL)=FLUXALL(IL)+FLUX_6H(IL,NN)
         IF (DBMAX(NN).LE.10.) FLUXA10(IL)=FLUXA10(IL)+FLUX_6H(IL,NN)
         IF (DBMAX(NN).LE.2.) FLUXA2(IL)=FLUXA2(IL)+FLUX_6H(IL,NN)
         IF (DBMAX(NN).LE.1.) FLUXA1(IL)=FLUXA1(IL)+FLUX_6H(IL,NN)
        ENDDO
C
C       ! SUM OVER INTERNAL TRACER CLASSES FOR AI AND CI MODES
C
C       ! accumulation mode
C
        DO NN=MINAI,MAXAI
         FLUX_AI(IL)=FLUX_AI(IL)+FLUX_6H(IL,NN)
        ENDDO
C
C       ! coarse mode
C
        DO NN=MINCI,MAXCI
         FLUX_CI(IL)=FLUX_CI(IL)+FLUX_6H(IL,NN)
        ENDDO
C
C       ! emission mass flux: g/m2/sec -> kg/m2/sec
C
        FLUX_AI(IL)=MAX(FLUX_AI(IL)*1.e-3,0.) 
        FLUX_CI(IL)=MAX(FLUX_CI(IL)*1.e-3,0.)
        FLUXALL(IL)=MAX(FLUXALL(IL)*1.e-3,0.)
        FLUXA10(IL)=MAX(FLUXA10(IL)*1.e-3,0.)
        FLUXA2(IL)=MAX(FLUXA2(IL)*1.e-3,0.)
        FLUXA1(IL)=MAX(FLUXA1(IL)*1.e-3,0.)
C
C       ! SAVE DUST EMISSION FLUXES
C
        FAC= ZTMST*GRAV*FLNDROW(IL)/(DSHJ(IL,ILEV)*PRESSG(IL))
        XEMIS(IL,ILEV,IDUA)=XEMIS(IL,ILEV,IDUA)+FLUX_AI(IL)*FAC
        XEMIS(IL,ILEV,IDUC)=XEMIS(IL,ILEV,IDUC)+FLUX_CI(IL)*FAC
C
        IF (ISVDUST.GT.0) THEN
         ESDROW (IL)=ESDROW (IL)+FLUX_AI(IL)*FLNDROW(IL)*SAVERAD
         ESDROW (IL)=ESDROW (IL)+FLUX_CI(IL)*FLNDROW(IL)*SAVERAD
C
         FALLROW(IL)=FALLROW(IL)+FLUXALL(IL)*FLNDROW(IL)*SAVERAD
         FA10ROW(IL)=FA10ROW(IL)+FLUXA10(IL)*FLNDROW(IL)*SAVERAD
         FA2ROW(IL) =FA2ROW (IL)+FLUXA2 (IL)*FLNDROW(IL)*SAVERAD
         FA1ROW(IL) =FA1ROW (IL)+FLUXA1 (IL)*FLNDROW(IL)*SAVERAD
C
        ENDIF
c        write(0,*) IL,SAVERAD,XEMIS(IL,ILEV,IDUA)
c        write(0,*) IL,XEMIS(IL,ILEV,IDUA),XEMIS(IL,ILEV,IDUC),
c     &         FLUX_AI(IL),FLUX_CI(IL),FAC,ESDROW(IL),SAVERAD
c       write(0,*) "DUSTEM: ",IL,FAC,ZTMST,GRAV
C
       ENDIF ! (DUST_MASK and C_EFF)
C
C       ! SAVE DUST WIND COMPONENTS     
C
        IF (ISVDUST.GT.0) THEN
         DUWDROW(IL)=DUWDROW(IL)+DU_WIND(IL)*SAVERAD
         DUSTROW(IL)=DUSTROW(IL)+USTAR_D(IL)*SAVERAD
         DUTHROW(IL)=DUTHROW(IL)+DU_TH(IL)*SAVERAD
         USMKROW(IL)=USMKROW(IL)+USTAR_ACRIT(IL)*SAVERAD
        ENDIF
C 
       ENDDO ! IL
C
C   --------------------------------------------------------------
C
      RETURN
C
      END
