      SUBROUTINE DUINIT1(UTH,SREL,SRELV,SU_SRELV,CUSCAL)
!-----------------------------------------------------------------------
!     PURPOSE:
!     --------
!     INITIALIZE SOIL FIELDS FOR DUST EMISSION.
!
!     HISTORY:
!     --------
!     * MAY 28/2013 - M. LAZARE. NEW VERSION FOR GCM17:
!     *                          - "CUSCAL" PASSED IN (FROM DUSTEMI FOR
!     *                            BULK, AS "CUSCALE", OR FROM DUFCC FOR
!     *                            PLA, AS "CUSCPLA") RATHER THAN BEING
!     *                            OBTAINED FROM MODULE "DUPARM", SO THAT
!     *                            THIS ROUTINE WORKS WITH BULK OR PLA
!     *                            APPROACHES.
!     * JAN 29/2008 - Y. PENG    PREVIOUS VERSION DUINIT.
!
!-----------------------------------------------------------------------
!
      USE DUPARM1
!
!      * INTERNAL WORK FIELDS.
!
        REAL UTH(NCLASS)
        REAL SREL(NATS,NCLASS),SRELV(NATS,NCLASS)
        REAL SU_SRELV(NATS,NCLASS)
!
!      * AUXILIARY VARIABLES
!
        REAL AAA, BB, CCC, DDD, EE, FF
        REAL RDP, RSIZE(NCLASS), DLASTJ(NCLASS)
        REAL XKU, XK, XL, XM, XN, XNV
        REAL STOTAL,STOTALV,SU,SUV,SU_LOC,SU_LOCV
        REAL SU_CLASS(NCLASS),SU_CLASSV(NCLASS)
        REAL UTEST(NATS)
        REAL CUSCAL
!
        INTEGER I,J,KK
        INTEGER NS
        INTEGER ND,NSI,NP
!
!     * INITIALIZATION
!
!      ! THRESHOLD FRICTION VELOCITY UTH
!      ! MARTICORENA AND BERGAMETTI 1995 EQUATIONS (3) and (4)
!
       I=0
!
       RDP=DMIN
       DO WHILE(RDP .LE. DMAX+1.E-5)
        I=I+1
        RSIZE(I)=RDP
        BB = A_RNOLDS * (RDP ** X_RNOLDS) + B_RNOLDS
        XKU = SQRT(ROP * GRAVI * RDP / ROA)
        CCC = SQRT(1. + D_THRSLD /(RDP ** 2.5))
         IF (BB.LT.10.) THEN
          DDD=SQRT(1.928 * (BB ** 0.092) - 1.)
          UTH(I) = 0.129 * XKU * CCC / DDD
         ELSE
          EE = -0.0617 * (BB - 10.)
          FF = 1. -0.0858 * EXP(EE)
          UTH(I) = 0.12 * XKU * CCC * FF
         ENDIF
        RDP = RDP * EXP(DSTEP)
       END DO
!
!      ! SCALE THE WIND STRESS THRESHOLD
!
       DO I=1, NCLASS
        UTH(I)=UTH(I)*CUSCAL
       END DO
!
!      ! LOOP OVER SOIL TYPES NS
!      ! CALCULATE SOIL DISTRIBUTION AND RELATED SURFACES
!
       DO NS=1,NATS
        UTEST(NS)=0.
!      ! initial surface values
        RDP = DMIN
        KK = 0
        STOTAL = 0.
        STOTALV = 0.
        DO I=1,NCLASS
         SU_CLASS(I) = 0.
         SU_CLASSV(I) = 0.
        END DO  ! I NCLASS
!      ! surface calculation
        DO WHILE (RDP.LE.DMAX+1.E-5)
         KK = KK + 1
         SU = 0.
         SUV = 0.
!
         DO I = 1, NMODE
          ND  = ((I - 1) *3 ) + 1
          NSI = ND + 1
          NP  = ND + 2
          IF (SOLSPE(ND,NS).EQ.0.) THEN
           SU_LOC =0.
           SU_LOCV=0.
          ELSE
           XK = SOLSPE(NP,NS)/(SQRT(2.*PI)*LOG(SOLSPE(NSI,NS)))
           XL = ((LOG(RDP)-LOG(SOLSPE(ND,NS)))**2)/ &
                (2.*(LOG(SOLSPE(NSI,NS)))**2)
           XM = XK * EXP(-XL)
           XN = ROP*(2./3.)*(RDP/2.)
           XNV = 1.
           SU_LOC = (XM*DSTEP/XN)
           SU_LOCV = (XM*DSTEP/XNV)
          ENDIF
          SU  = SU  + SU_LOC
          SUV = SUV + SU_LOCV
         END DO
!
         SU_CLASS(KK) = SU
         SU_CLASSV(KK)= SUV
         STOTAL = STOTAL + SU
         STOTALV= STOTALV+SUV
         DLASTJ(KK)=RDP*5000.
         RDP=RDP*EXP(DSTEP)
        END DO     ! RDP<=DMAX
!      ! related surface
        DO J = 1,NCLASS
         IF (STOTAL.EQ.0.) THEN
          SREL(NS,J) = 0.
          SRELV(NS,J)= 0.
         ELSE
          SREL(NS,J) = SU_CLASS(J) /STOTAL
          SRELV(NS,J)= SU_CLASSV(J)/STOTALV
          UTEST(NS) = UTEST(NS)+SRELV(NS,J)
          SU_SRELV(NS,J) = UTEST(NS)
         ENDIF
        END DO    ! J NCLASS
!
       END DO     ! NS NATS
!
!      ! LOOP OVER SOIL TYPES NS ENDS
!
      END SUBROUTINE DUINIT1
