      SUBROUTINE DUFCC (FLND,GT,SMFR,FN,ZSPD,USTARD, &
                        SUZ0,SLAI,SPOT,ST02,ST03,ST04,ST06, &
                        ST13,ST14,ST15,ST16,ST17, &
                        UTH,SREL,SRELV,SU_SRELV, &
                        TDFLX,TDSIZE,TDVAR,TDMASK,ILGA,ISEC, &
                        ISVDUST,DEFA,DEFC, &
                        FALL,FA10,FA2,FA1, &
                        DUWD,DUST,DUTH,USMK, &
                        DIAGMAS,DIAGNUM,ISDUST)
!
!     * THIS ROUTINE CALCULATES THE SURFACE DUST EMISSION FLUX
!     * IN TWO SECTIONS FOR PLA
!
!     * FEB 05/2015 - M.LAZARE/     NEW VERSION FOR GCM18:
!     *               K.VONSALZEN.  - FLND PASSED IN INSTEAD OF {GC,MASK},
!     *                               FOR FRACTIONAL LAND MASK SUPPORT.
!     *                               THIS ELIMINATES THE BUG WHERE HAD "MASK(IL).EQ.-0.5".
!     *                             - {ZSPD,USTARD} PASSED IN INSTEAD OF {SFCU,SFCV,GUST}
!     *                               BECAUSE THIS SCREEN-LEVEL QUANTITY
!     *                               ALREADY CALCULATED IN EARLIER ROUTINES.
!     *                               THIS CHANGES AND SIMPLIFIES THE CALCULATION OF
!     *                               USTAR.
!     *                             - CORRECT BUG WHERE FLUXMAS WAS SET TO ZERO FOR NON-DRY
!     *                               SOILS **BEFORE** ACTUAL CALCULATION (THUS NEVER "SEEN"),
!     *                               BY MOVING IT TO AFTER THE CALCULATION.
!     * NOV 14/2013 - M. LAZARE.    CUSCPLA NOT PASSED IN CALL FROM DUSTP
!     *                             AND REAL DECLARATION REMOVED.
!     * MAY 28/2013 - M. LAZARE     INCLUDE "CUSCPLA" IN CALL
!     *                             RATHER THAN THROUGH MODULE "DUPARM1"
!     * FEB 26/2010 - K. VONSALZEN  CORRECT USE OF DIAGMAS,DIAGNUM
!     * OCT 2009    - Y. PENG       REVISED FOR PLA SCHEME
!
!     * AUG 20/2007 - Y. PENG       NEW SCHEME INCL. SOIL TEXTURE:
!     *                             - USE ZOBLER SOIL PROPERTIES
!     *                             - USE MARTICORENA ET AL 1997 TO
!     *                               CALCULATE THE THRESHOLD WIND SPEED
!     *                             - WORK ARRAYS NOW LOCAL.
!     *
!
       USE DUPARM1
!
       IMPLICIT REAL (A-H,O-Z),INTEGER (I-N)
!
!      * I/O FIELDS.
!
       REAL, INTENT(OUT), DIMENSION(ILGA,ISDUST*2) :: DIAGMAS,DIAGNUM
       REAL, INTENT(IN),  DIMENSION(ILGA) :: FLND,GT,SMFR,FN, &
                                             ZSPD,USTARD,SUZ0,SLAI, &
                                             SPOT,ST02,ST03,ST04,ST06, &
                                             ST13,ST14,ST15,ST16,ST17
!
       REAL, DIMENSION(ILGA,ISEC) :: TDFLX,TDSIZE,TDVAR
       REAL, DIMENSION(ILGA) :: TDMASK
!
       REAL, INTENT(OUT), DIMENSION(ILGA) :: DEFA,DEFC,FALL,FA10,FA2, &
                                             FA1,DUWD,DUST,DUTH,USMK
!
!      * INTERNAL WORK FIELDS.
!
        REAL, DIMENSION(NCLASS) :: UTH
        REAL, DIMENSION(NATS,NCLASS) :: SREL,SRELV,SU_SRELV
!
        REAL, DIMENSION(NTRACE) :: DPK,DBMIN,DBMAX
        REAL, DIMENSION(ILGA) :: C_EFF,DUST_MASK,USTAR_D, &
                                 DU_WIND,DU_TH, &
                                 USTAR_ACRIT,DUZ01,DUZ02
        REAL, DIMENSION(ILGA) :: ST01
!
        REAL, DIMENSION(ILGA,NTRACE) :: FLUX_6H,FLUXBIN
        REAL, DIMENSION(ILGA,NCLASS) :: FLUXTYP, &
                                        FLUXDIAM1,FLUXDIAM2, &
                                        FLUXDIAM3,FLUXDIAM4, &
                                        FLUXDIAM6,FLUXDIAM_PF, &
                                        FLUXDIAM13,FLUXDIAM14, &
                                        FLUXDIAM15,FLUXDIAM16, &
                                        FLUXDIAM17
        REAL, DIMENSION(ILGA) :: FLUXALL,FLUXA10,FLUXA2,FLUXA1, &
                                 FLUX_AI,FLUX_CI
!
!      * AUXILIARY VARIABLES
!
        REAL DPB(NCLASS+1),DPC(NCLASS)
        REAL FLUXMAS(ILGA,NCLASS)
        REAL FLUXNUM(ILGA,NCLASS)
        REAL DRB(ISEC),DRE(ISEC),DRM(ISEC)
        REAL FLUX_NUM(ILGA,ISEC)
        INTEGER INB(ISEC),INE(ISEC)
!
        REAL AAA, BB, CCC, DDD, EE, FF
        REAL D1, FEFF
        REAL RDP, DLAST
        REAL USTAR,UTHP(ILGA,NATS)
        REAL ALPHA
        REAL WAA(NATS),WBB(ILGA)
        REAL FDP1, FDP2
!
        INTEGER IL,J, KK, NN, I_SOIL
        INTEGER KKMIN,KKK
        INTEGER DN(ILGA),DK(ILGA,NCLASS)
!
        INTEGER ZSCHM,WSCHM,SOLSCHM
!        z0 scheme, 1: new scheme Prigent et al 2002, 0: old scheme z01=z02=0.001
        DATA ZSCHM /1/
!        soil water scheme, 1: new scheme Fecan, 0: old scheme 0.99
        DATA WSCHM /1/
!        soil property scheme, 1: new scheme Cheng, 0: old scheme
        DATA SOLSCHM /0/
!    --------------------------------------------------------------
!     * INITIALIZATION OF DIAGNOSTIC ARRAYS.
!
      IF (ISVDUST.GT.0) THEN
        DEFA=0.
        DEFC=0.
        FALL=0.
        FA10=0.
        FA2=0.
        FA1=0.
        DUWD=0.
        DUST=0.
        DUTH=0.
        USMK=0.
      ENDIF

!     * INITIAL CALCULATION (C_EFF)
!
!      ! EFFECTIVE FRACTION CALCULATION (ROUGHNESS Z0 RELATED)
!      ! CONSTANT Z0 TEMPORARILY
!      ! MARTICORENA ET AL 1997 EQUATIONS (15)-(17)
!
!      ! initial values !D1 is set to zero temporarily
!
       D1 = 0.
       FEFF = 0.
       AAA = 0.
       BB = 0.
       CCC = 0.
       DDD = 0.
       EE = 0.
       FF = 0.
       DO IL=1,ILGA
        C_EFF(IL)       = 0.0
        DUZ01(IL)       = 0.0
        DUZ02(IL)       = 0.0
       END DO
!
       DO IL=1,ILGA
!
       IF (ZSCHM .EQ. 1) THEN
!       new Z0 scheme, monthly average
        DUZ01(IL)=SUZ0(IL)
        DUZ02(IL)=SUZ0(IL)
       ELSE
!       old z0 scheme
        DUZ01(IL)=Z01
        DUZ02(IL)=Z02
       ENDIF
!
       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_EFF(IL) = FEFF
!
       END DO  !IL
!
!     * INITIAL CALCULATION (C_EFF) DONE
!
!   --------------------------------------------------------------
!
!     * INITIAL CALCULATION (WINDS and MASKS)
!
!      ! initial values
!
       USTAR=0.0
       RDP=0.0
       DLAST=0.0
       DO KK=1,NCLASS
        DPB(KK)  = 0.0
        DPC(KK)  = 0.0
       END DO
       DO IL=1,ILGA
        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
       END DO
!
!      ! CALCULATE BIN MINIMUM/MAXIMUM RADIUS IN UM
!      ! RDP IS DIAMETER IN CM
!
       RDP=DMIN
       DPB(1)=RDP
       DLAST=DMIN
       DO NN=1,NTRACE
        DPK(NN)=0.
        DBMIN(NN)=0.
        DBMAX(NN)=0.
       END DO
!
       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)
        DPB(KK+1)=RDP
       END DO
!
!      ! CALCULATE THE FRICTION WIND SPEED USTAR
!
       DO IL=1,ILGA
!      ! set wind speed, use surface (10m) wind instead of ZSPD
         DU_WIND(IL)=ZSPD(IL)
         IF (DUZ02(IL) .EQ. 0.) THEN
          USTAR=0.
         ELSE
          USTAR=100.*USTARD(IL)
         ENDIF
         USTAR_D(IL)=USTAR
!
!         IF (C_EFF(IL) .EQ. 0.) THEN
          DU_TH(IL)=0.
!         ELSE
!          DU_TH(IL)=UMIN*CUSCPLA/C_EFF(IL)
!         ENDIF
!
!      ! check critical wind speed (wind stress threshold scaled)
          IF (USTAR.GT.0..AND.USTAR.GE.DU_TH(IL)) THEN
           IF (SLAI(IL).GT.0.0 .AND. DUZ02(IL).GT.0.0) THEN
!      ! set grid point as a potential dust source point
            DUST_MASK(IL)=1.
           ENDIF
!      ! store the time fraction with ustar > dust threshold for output
           USTAR_ACRIT(IL)=1.
          ELSE
           USTAR_ACRIT(IL)=0.
          ENDIF   ! USTAR
!
       END DO    ! IL
!
!
!     * INITIAL CALCULATION (WINDS and MASKS) DONE
!
!   --------------------------------------------------------------
!
!     * EMISSION FLUX CALCULATION
!
!      ! initialization
!
        ALPHA=0.0
        FDP1=0.0
        FDP2=0.0
       DO IS=1,ISEC
        DRB(IS)=0.0
        DRE(IS)=0.0
        DRM(IS)=0.0
        INB(IS)=0.0
        INE(IS)=0.0
       ENDDO
       DO J=1,NATS
        WAA(J)=0.0
       ENDDO
!
       DO IL=1,ILGA
        WBB(IL)=0.0
        DN(IL) =0
        ST01(IL)=1.
!
        DO IS=1,ISEC
         FLUX_NUM(IL,IS)=0.
        ENDDO
!
        DO NN=1,NTRACE
         FLUXBIN(IL,NN)=0.
         FLUX_6H(IL,NN)=0.
        END DO
!
        DO KK=1,NCLASS
         DK(IL,KK)     =0
         FLUXTYP(IL,KK)=0.
         FLUXMAS(IL,KK)=0.
         FLUXNUM(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
!
       END DO ! IL
!
!      ! calculation
!
       CD=1.*ROA/GRAVI
!
       DO KK=1,NCLASS
        DO IL=1,ILGA
         IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
!
!   ! soil moisture modification
!
          IF (WSCHM .EQ. 1) THEN
!          new ws scheme
            IF (FLND(IL).GT.0. .AND. GT(IL).GT.273.15) THEN
             WBB(IL)=MIN(SMFR(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* &
               (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
!          old ws scheme
            DO J=1,NATS
             UTHP(IL,J)=UTH(KK)
            ENDDO !J
          ENDIF !soil moisture
!
!    ! fluxdiam calculation (FDP1,FDP2 follows Marticorena)
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! add Asian soil types
           IF (SOLSCHM.EQ.1) THEN
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
!          ! 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) &
             *FDP1*FDP2*CD*USTAR_D(IL)**3*ALPHA)
            ENDIF
!
           ENDIF !Asian soil
!
!    ! initialize fluxtyp of all soil types at first bin (kk=1)
!
            IF (KK.EQ.1) THEN
             IF (SOLSCHM.EQ.1) THEN
             FLUXTYP(IL,KK)=FLUXTYP(IL,KK) &
                       +FLUXDIAM1(IL,KK)*(1.-SPOT(IL))* &
          (ST01(IL)-ST02(IL)-ST03(IL)-ST04(IL)-ST06(IL) &
          -ST13(IL)-ST14(IL)-ST15(IL)-ST16(IL)-ST17(IL)) &
                       +FLUXDIAM2(IL,KK)*(1.-SPOT(IL))*ST02(IL) &
                       +FLUXDIAM3(IL,KK)*(1.-SPOT(IL))*ST03(IL) &
                       +FLUXDIAM4(IL,KK)*(1.-SPOT(IL))*ST04(IL) &
                       +FLUXDIAM6(IL,KK)*(1.-SPOT(IL))*ST06(IL) &
                       +FLUXDIAM13(IL,KK)*(1.-SPOT(IL))*ST13(IL) &
                       +FLUXDIAM14(IL,KK)*(1.-SPOT(IL))*ST14(IL) &
                       +FLUXDIAM15(IL,KK)*(1.-SPOT(IL))*ST15(IL) &
                       +FLUXDIAM16(IL,KK)*(1.-SPOT(IL))*ST16(IL) &
                       +FLUXDIAM17(IL,KK)*(1.-SPOT(IL))*ST17(IL) &
                       +FLUXDIAM_PF(IL,KK)*SPOT(IL)
             ELSE
             FLUXTYP(IL,KK)=FLUXTYP(IL,KK) &
                       +FLUXDIAM1(IL,KK)*(1.-SPOT(IL))* &
          (ST01(IL)-ST02(IL)-ST03(IL)-ST04(IL)-ST06(IL)) &
                       +FLUXDIAM2(IL,KK)*(1.-SPOT(IL))*ST02(IL) &
                       +FLUXDIAM3(IL,KK)*(1.-SPOT(IL))*ST03(IL) &
                       +FLUXDIAM4(IL,KK)*(1.-SPOT(IL))*ST04(IL) &
                       +FLUXDIAM6(IL,KK)*(1.-SPOT(IL))*ST06(IL) &
                       +FLUXDIAM_PF(IL,KK)*SPOT(IL)
             ENDIF !Asian soil
            ELSE
             DN(IL)=DN(IL)+1           ! increase number of dislocated particle classes
             DK(IL,DN(IL))=KK          ! store dislocated particle class
            ENDIF  ! KK=1
!
         ENDIF     ! DUST_MASK=1 and C_EFF>1
        ENDDO      ! IL
       ENDDO       ! KK
!
!      ! calculate fluxtyp of all soil types for all 192 bins (kkk=2~192)
!
       KKMIN=1
       DO IL=1,ILGA
         IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
          DO N=1,DN(IL)
           DO KKK=1,DK(IL,N)           ! scaling with relative contribution of dust size fraction
             I_SOIL=1
             IF (SOLSCHM.EQ.1) THEN
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
         *(ST01(IL)-ST02(IL)-ST03(IL)-ST04(IL)-ST06(IL) &
          -ST13(IL)-ST14(IL)-ST15(IL)-ST16(IL)-ST17(IL)) &
         *FLUXDIAM1(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
         ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
             ELSE
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
         *(ST01(IL)-ST02(IL)-ST03(IL)-ST04(IL)-ST06(IL)) &
         *FLUXDIAM1(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
         ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
             ENDIF  ! Asian soil
!
             I_SOIL=2
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST02(IL) &
          *FLUXDIAM2(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=3
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST03(IL) &
          *FLUXDIAM3(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=4
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST04(IL) &
          *FLUXDIAM4(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=6
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST06(IL) &
          *FLUXDIAM6(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
            IF (DU_WIND(IL).GT. 8.)  I_SOIL=11
            IF (DU_WIND(IL).LE. 8.)  I_SOIL=10
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+SPOT(IL) &
          *FLUXDIAM_PF(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
          IF (SOLSCHM.EQ.1) THEN
             I_SOIL=13
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST13(IL) &
          *FLUXDIAM13(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=14
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST14(IL) &
          *FLUXDIAM14(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=15
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST15(IL) &
          *FLUXDIAM15(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=16
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST16(IL) &
          *FLUXDIAM16(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
             I_SOIL=17
             FLUXTYP(IL,KKK) = FLUXTYP(IL,KKK)+(1.-SPOT(IL)) &
          *ST17(IL) &
          *FLUXDIAM17(IL,DK(IL,N))*SRELV(I_SOIL,KKK)/ &
          ((SU_SRELV(I_SOIL,DK(IL,N))-SU_SRELV(I_SOIL,KKMIN)))
!
          ENDIF !Asian soil
!
           ENDDO  ! KKK
          ENDDO   ! N
         ENDIF    ! DUST_MASK
       ENDDO      ! IL
!
!      ! CALCULATE FLUXES IN EACH INTERNAL TRACER CLASSES (8)
!
       DO NN=1,NTRACE
        DO IL=1,ILGA
         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 (SMFR(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 (FN(IL).GT.0.1) FLUXBIN(IL,NN)=0.
           FLUX_6H(IL,NN)=FLUXBIN(IL,NN)*10000.*SLAI(IL)
         ENDIF ! DUST_MASK
        ENDDO  ! IL
       ENDDO   ! NN
!
!     * EMISSION FLUX CALCULATION DONE
!
!   --------------------------------------------------------------
!
!     * DIAGNOSTIC CALCULATION FOR MODEL OUTPUT
!
!      ! initialization for output arrays
!
       DO IL=1,ILGA
        FLUXALL(IL)=0.
        FLUXA10(IL)=0.
        FLUXA2(IL)=0.
        FLUXA1(IL)=0.
        FLUX_AI(IL)=0.
        FLUX_CI(IL)=0.
        TDFLX(IL,:)=0.
        TDSIZE(IL,:)=0.
        TDVAR(IL,:)=0.
        TDMASK(IL)=0.
!
        DO ISF=1,ISDUST*2
         DIAGMAS(IL,ISF)= 0.0
         DIAGNUM(IL,ISF)= 0.0
        ENDDO
!
       ENDDO
!
       DO IL=1,ILGA
       IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
!
!       ! unit of emission flux: FLUX_6H [g/m2/sec]
!
        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
!
!       ! SUM OVER INTERNAL TRACER CLASSES FOR AI AND CI MODES
!
!       ! accumulation mode
!
        DO NN=MINAI,MAXAI
         FLUX_AI(IL)=FLUX_AI(IL)+FLUX_6H(IL,NN)
        ENDDO
!
!       ! coarse mode
!
        DO NN=MINCI,MAXCI
         FLUX_CI(IL)=FLUX_CI(IL)+FLUX_6H(IL,NN)
        ENDDO
!
!       ! emission mass flux: g/m2/sec -> kg/m2/sec
!
        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.)
!
!       ! SAVE DUST EMISSION FLUXES
!
        IF (ISVDUST.GT.0) THEN
         DEFA(IL)=FLUX_AI(IL)
         DEFC(IL)=FLUX_CI(IL)
!
         FALL(IL)=FLUXALL(IL)
         FA10(IL)=FLUXA10(IL)
         FA2(IL) =FLUXA2(IL)
         FA1(IL) =FLUXA1(IL)
        ENDIF
!
       ENDIF ! (DUST_MASK and C_EFF)
!
!       ! SAVE DUST WIND COMPONENTS
!
        IF (ISVDUST.GT.0) THEN
         DUWD(IL)=DU_WIND(IL)
         DUST(IL)=USTAR_D(IL)
         DUTH(IL)=DU_TH(IL)
         USMK(IL)=USTAR_ACRIT(IL)
        ENDIF
!
       ENDDO ! IL
!
!   --------------------------------------------------------------
!
!     * DUST EMISSION FOR PLA SECTIONS AND DIAGNOSTIC FOR CHECKING
!
!       ! unit of radius
!
       DO KK=1,NCLASS+1
        DPB(KK)=DPB(KK)*0.5*1.e4          ! convert unit to um
       ENDDO
!
!       ! dust emission size distri. (for diagnostics)
!
       DO IL=1,ILGA
        IF (DUST_MASK(IL).EQ.1 .AND. C_EFF(IL).GT.0.) THEN
!
        DO KK=1,NCLASS
!       ! unit of dust emission mass flux
!       ! fluxtyp: g/cm2/sec; fluxmas: g/m2/sec
         FLUXMAS(IL,KK)=FLUXTYP(IL,KK)*10000.*SLAI(IL)
        ! mask out dust fluxes, where soil moisture threshold is reached
         IF (SMFR(IL).GT.W0) FLUXMAS(IL,KK) = 0.
        ! mask out dust fluxes where snow present
         IF (FN(IL).GT.0.1)  FLUXMAS(IL,KK)=0.
        ENDDO
!
!       ! convert mass flux to number conc.
!
        DO KK=1,NCLASS
         DPC(KK)=DPB(KK)+(DPB(KK+1)-DPB(KK))/2.
         FLUXNUM(IL,KK)=1.e12*FLUXMAS(IL,KK)/(ROP* &
                          4./3.*PI*DPC(KK)**3)
!
!       ! fluxmas: g/m2/sec; fluxnum: 1/m2/sec
!
         IF (ISVDUST.GT.0) THEN
          DIAGMAS(IL,KK)=FLUXMAS(IL,KK)/LOG(DPB(KK+1)/DPB(KK))
          DIAGNUM(IL,KK)=FLUXNUM(IL,KK)/LOG(DPB(KK+1)/DPB(KK))
         ENDIF
        ENDDO
!
!        ! diagmas: kg/m2/sec; diagnum: 1/m2/sec
!
        DO ISF=1,ISDUST*2
         DIAGMAS(IL,ISF)=MAX(DIAGMAS(IL,ISF)*1.e-3,0.)
         DIAGNUM(IL,ISF)=MAX(DIAGNUM(IL,ISF),0.)
        ENDDO
!
        ENDIF  ! (DUST_MASK and C_EFF)
       ENDDO   ! IL
!
!      ! write to dust emission size distri. (for diagnostics)
!
!      ! setup boundary size of each pla section
!
       IF ( ISEC==1 ) THEN
         INB(1)=MINAI
         INE(1)=MAXCI
       ELSE IF ( ISEC==2 ) THEN
         INB(1)=MINAI
         INE(1)=MAXAI
         INB(2)=MINCI
         INE(2)=MAXCI
       ELSE
         CALL XIT('DUFCC',-1)
       ENDIF
!
       DO IL=1,ILGA
!
        DO IS=1,ISEC
          DO KK=1,NCLASS
           IF (MOD(KK,NBIN) .EQ. 0 .AND. KK/NBIN .EQ. INE(IS)) THEN
            DRE(IS)=DPB(KK+1)
           ENDIF
!
           IF (KK .GT. (INB(IS)-1)*NBIN .AND. KK .LE. INE(IS)*NBIN) THEN
            FLUX_NUM(IL,IS) =FLUX_NUM(IL,IS) +FLUXNUM(IL,KK)
            DRM(IS)=DRM(IS)+FLUXNUM(IL,KK)*DPC(KK)
           ENDIF
          ENDDO  ! KK
        ENDDO    ! IS
!
         DRB(1) = DPB(1)
        DO IS=2,ISEC
         DRB(IS)=DRE(IS-1)
        ENDDO ! IS
!
!      ! dust emission size distri. parameters in pla sections
!
        INO=0
        DO IS=1,ISEC
          IF ( FLUX_NUM(IL,IS) <= 0. ) INO=1
        ENDDO
        IF (INO==0) THEN
         TDMASK(IL)=1.
         DO IS=1,ISEC
          DRM(IS)=DRM(IS)/FLUX_NUM(IL,IS)
          TDSIZE(IL,IS)=DRM(IS)*1.e-6
          TDFLX(IL,IS) =FLUX_NUM(IL,IS)/LOG(DRE(IS)/DRB(IS))
         ENDDO    ! IS
         IF ( ISEC==1 ) THEN
           TDVAR(IL,1)=2.3
         ELSE IF ( ISEC==2 ) THEN
           TDVAR(IL,1)=2.5
           TDVAR(IL,2)=1.5
         ENDIF
        ENDIF
!
       ENDDO     ! IL
!
!   --------------------------------------------------------------
!
      RETURN
!
      END SUBROUTINE DUFCC
