SUBROUTINE Dust_Emission
c ****************************************************************************
c * This subroutine updates the surface mixing ratio of dust aerosols for
c * NMX size bins. The uplifting of dust depends in space on the source 
c * function, and in time and space on the soil moisture and surface
c * wind spedd (10 meters). Dust is uplifted if the wind speed is greater
c * than a threshold velocity which is calculated with the formula of
c * Marticorena et al.  (JGR, v.102, p 23277-23287, 1997).
c * To run this subroutine you need the source function which can be
c * obtained by contacting Paul Ginoux at ginoux@rondo.gsfc.nasa.gov
c * If you are not using GEOS DAS met fields, you will most likely need
c * to adapt the adjusting parameter.
c *
c * Contact: Paul Ginoux (ginoux@rondo.gsfc.nasa.gov) 
c *
c *
c *  Input:
c *         SRCE      Source function                               (-)
c *         DUSTDEN   Dust density                                  (kg/m3)
c *         DUSTREFF  Effective radius                              (um)
c *         DXY       Surface of each grid cell                     (m2)
c *         AIRMAS    Air mass for each grid box                    (kg)
c *         NTDT      Time step                                     (s)
c *         W10m      Velocity at the anemometer level (10meters)   (m/s)
c *         GWET      Surface wetness                               (-)
c *      
c *  Update:
c *         TC        Mixing ratio at first model level over Earth surface
c *
c *
c *  Parameters used in GOCART
c *
c *  Longitude: IMX = 144
c *  Latitude : JMX =  91
c *  Levels   : LMX =  20 (GEOS-1), 26 (GEOS-strat), 30 (GEOS-terra)
c *  Size bins: NMX =   4
c *
c *  Dust properties used in GOCART
c *
c *  Size classes: 01-1, 1-1.8, 1.8-3, 3-6 (um)
c *  Radius: 0.7, 1.5, 2.5, 4  (um)
c *  Density: 2500, 2650, 2650, 2650 (kg/m3)
c *
c *  References:
c *   1) Ginoux, P., M. Chin, I. Tegen, J. Prospero, B. Hoben, O. Dubovik,
c *      and S.-J. Lin, "Sources and distributions of dust aerosols simulated
c *      with the GOCART model", J. Geophys. Res., 2001
c *
c ****************************************************************************
c     
      parameter (imx=144,jmx=91,lmx=26,nmx=4)
c Input
      dimension GWET(imx,jmx),W10M(imx,jmx),SRCE(imx,jmx)
      dimension AIRMAS(imx,jmx,LMX),DXY(imx,jmx)
      dimension AIRDEN(imx,jmx,lmx)
c Mixing Ratio
      dimension TC(imx,jmx,lmx,nmx)
c
c Soil size distribution 
      dimension FRAC_S(nmx)
c
      dimension DUSTREFF(nmx),DUSTDEN(nmx)
c Particle radii (m)
      data DUSTREFF/0.75e-6,1.5e-6,2.5e-6,4e-6/
c Soil density (kg/m3)
      data DUSTDEN/2500.,2650.,2650.,2650./
c
c Fraction of each size classes (deduced from Tegen and Fung, 1994)
c
      data frac_s/0.1,0.25,0.25,0.25/
c
      Ch_dust    = 1.e-9   ! Transfert coeff for type natural source 
c                            ! (kg*s2/m5)
c Loop over each size bins
c
      do n = 1, nmx
c Threshold velocity as a function of the dust density and the diameter
c from Bagnold (1941), valid for particles larger than 10 um
c
c       u_ts0 = 6.5*sqrt(dustden(n)*g0*2.*dustreff(n))
c Threshold velocity from Greely and Iversen (1985)
c  Convert units to fit dimensional parameters
        den=dustden(n)*1.e-3
        diam=2.*dustreff(n)*1.e2
        g=9.81*1.e2

        TSRC = 0.
c
c Loop over each longitude
c
        do i = 1, IMX
c
c Loop over each latitude
c
          do j = 1, JMX
            rhoa=airden(i,j,1)*1.e-3

            u_ts0=0.13*1.e-2*sqrt(den*g*diam/rhoa)*
     >        sqrt(1.+0.006/den/g/diam**2.5)/
     >        sqrt(1.928*(1331*(diam)**1.56+0.38)**0.092-1)
c Case of surface dry enough to erode
            if (gwet(i,j).lt.0.5) then
              u_ts=amax1(0.,u_ts0*(1.2+0.2*alog10(amax1(1.e-3,
     >           gwet(i,j)))))
            else
c Case of wet surface, no erosion
              u_ts=100.
            endif
            SRCE_p=frac_s(n)*SRCE(i,j)*DXY(i,j)  ! (m2)
            DSRC=Ch_dust*SRCE_p*W10m(i,j)**2
     >                * (W10m(i,j)-u_ts)*NTDT  ! (kg)
            if (DSRC.lt.0.) DSRC=0.
c
c Update dust mixing ratio at first model level.
c
            TC(i,j,1,n) = TC(i,j,1,n) + DSRC / AIRMAS(i,j,1)
          enddo
        enddo
      enddo
c    
      RETURN
      END