MODULE GOCART_DUST_AFWA
!
! AFWA dust routine
! Created by Sandra Jones (AER and AFWA) and Glenn Creighton (AFWA).
!

  USE module_data_gocart_dust

  IMPLICIT NONE

  INTRINSIC max, min

CONTAINS
  SUBROUTINE gocart_dust_afwa_driver(ktau,dt,config_flags,julday,alt,t_phy,moist,u_phy,  &
         v_phy,chem,rho_phy,dz8w,smois,u10,v10,p8w,erod,erod_dri,dustin,snowh,zs,   &
         ivgtyp,isltyp,vegfra,lai_vegmask,xland,xlat,xlong,gsw,dx,g,emis_dust,      &
         ust,znt,clay_wrf,sand_wrf,clay_nga,sand_nga,afwa_dustloft,                 &
         tot_dust,tot_edust,vis_dust,alpha,gamma,smtune,ustune,            &
         ids,ide, jds,jde, kds,kde,                                        &
         ims,ime, jms,jme, kms,kme,                                        &
         its,ite, jts,jte, kts,kte                                         )
  USE module_configure
  USE module_state_description
  USE module_data_sorgam, ONLY: factnuma,factnumc,soilfac

   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags

   INTEGER,      INTENT(IN   ) :: julday, ktau,                            &
                                  ids,ide, jds,jde, kds,kde,               &
                                  ims,ime, jms,jme, kms,kme,               &
                                  its,ite, jts,jte, kts,kte
   INTEGER,DIMENSION( ims:ime , jms:jme )                  ,               &
          INTENT(IN   ) ::                                                 &
                                                     ivgtyp,               &
                                                     isltyp
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),                &
         INTENT(IN ) ::                              moist
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                 &
         INTENT(INOUT ) ::                           chem
   REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL,           &
         INTENT(INOUT ) ::                                                 &
         emis_dust
   REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) ,     &
         INTENT(IN   ) ::                            smois
   REAL, DIMENSION( config_flags%num_soil_layers ) ,                       &
         INTENT(IN   ) ::                            zs
   REAL, DIMENSION( ims:ime , jms:jme, ndcls )             ,               &
         INTENT(IN   ) ::                            erod,erod_dri
   REAL, DIMENSION( ims:ime , jms:jme, 5 )                 ,               &
         INTENT(INOUT) ::                            dustin
   REAL, DIMENSION( ims:ime , jms:jme )                    ,               &
         INTENT(IN   ) ::                                                  &
                                                     u10,                  &
                                                     v10,                  &
                                                     gsw,                  &
                                                     vegfra,               &
                                                     lai_vegmask,          &
                                                     xland,                &
                                                     xlat,                 &
                                                     xlong,                &
                                                     ust,                  &
                                                     znt,                  &
                                                     clay_wrf,             &
                                                     sand_wrf,             &
                                                     clay_nga,             &
                                                     sand_nga,             &
                                                     snowh
   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                         &
         INTENT(IN   ) ::                                                  &
                                                     alt,                  &
                                                     t_phy,                &
                                                     dz8w,p8w,             &
                                                     u_phy,v_phy,rho_phy
  REAL, DIMENSION( ims:ime , jms:jme )                    ,                &
         INTENT(  OUT) ::                            afwa_dustloft,        &
                                                     tot_edust
  REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),                          &
         INTENT(  OUT) ::                            tot_dust,             &
                                                     vis_dust
  REAL, INTENT(IN   ) :: dt,dx,g

  ! Local variables

  INTEGER :: nmx,smx,i,j,k,imx,jmx,lmx,lhave
  INTEGER,DIMENSION (1,1) :: ilwi
  REAL*8, DIMENSION (1,1) :: erodtot
  REAL*8, DIMENSION (1,1) :: vegmask
  REAL*8, DIMENSION (1,1) :: gravsm
  REAL, DIMENSION( ims:ime , jms:jme ) :: clay,sand
  REAL*8, DIMENSION (1,1) :: drylimit
  REAL*8, DIMENSION (5)   :: tc,bems
  REAL*8, DIMENSION (1,1) :: airden,airmas,ustar
  REAL*8, DIMENSION (1) :: dxy
  REAL*8, DIMENSION (3) :: massfrac
  REAL*8                :: volsm
  REAL :: conver,converi
  REAL :: psi,ustart,w10
  REAL*8 :: zwant
  REAL, INTENT(IN   ) :: alpha, gamma, smtune, ustune
  INTEGER :: smois_opt

  conver=1.e-9
  converi=1.e9

  ! Number of dust bins

  imx=1
  jmx=1
  lmx=1
  nmx=ndust
  smx=ngsalt

  k=kts
  DO j=jts,jte
  DO i=its,ite

    ! Masked value for afwa_dustloft

    afwa_dustloft(i,j)=-99.

    ! Don't do dust over water!!!

    IF (xland(i,j) .lt. 1.5) THEN
      ilwi(1,1)=1

      ! Total concentration at lowest model level. This is still hardcoded
      ! for 5 bins.

      IF(config_flags%chem_opt == CB05_SORG_AQ_KPP .or. &
         config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP ) then
        tc(1)=0.0
        tc(2)=0.0
        tc(3)=0.0
        tc(4)=0.0
        tc(5)=0.0
      ELSE
        tc(1)=chem(i,kts,j,p_dust_1)*conver
        tc(2)=chem(i,kts,j,p_dust_2)*conver
        tc(3)=chem(i,kts,j,p_dust_3)*conver
        tc(4)=chem(i,kts,j,p_dust_4)*conver
        tc(5)=chem(i,kts,j,p_dust_5)*conver
      END IF

      ! Air mass and density at lowest model level.

      airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*dx*dx/g
      airden(1,1)=rho_phy(i,kts,j)
      ustar(1,1)=ust(i,j)
      dxy(1)=dx*dx

      ! Friction velocity tuning constant (Note: recommend 0.7 for PXLSM,
      ! else use 1.0.  This was created due to make the scheme compatible
      ! with the much stronger friction velocities coming out of PXLSM).

      ustar(1,1)=ustar(1,1) * ustune

      ! Total erodibility.  Determine what DSR we're using.
      ! Note, the DRI erodibility dataset is an optional 1km resolution dataset
      ! but currently is only available over southwest Asia. If running
      ! a domain outside of SWA, this will fill in the missing data from
      ! DRI with the Ginoux dataset, but there will be inconsistencies. GAC
      ! Erodibility is broken up into 3 bins, sum them up here.

      IF (config_flags%dust_dsr .eq. 1) then  ! DRI DSR
          IF (erod_dri(i,j,1).ge.0) THEN  ! Where DRI is defined
             erodtot(1,1) = SUM(erod_dri(i,j,:))
          ELSE ! Outside where DRI not defined, use Ginoux
             erodtot(1,1)=SUM(erod(i,j,:))
          ENDIF
      ELSE ! Ginoux DSR
          erodtot(1,1)=SUM(erod(i,j,:))
      ENDIF

      ! Set the vegmask variable to the desired vegation mask at this gridpoint

      IF (config_flags%dust_veg .eq. 1) then

         ! 12-month vegetation fraction
         ! If user chose this 12-month greenfrac vegetation mask option,
         ! cut off everything above 5%

         IF (vegfra(i,j) .ge. 5.) then
             vegmask(1,1) = 0.0
         ELSE
             vegmask(1,1) = 1.0
         ENDIF
      ELSE IF (config_flags%dust_veg .eq. 2) then

         ! 8-day MODIS LAI vegmask
         ! 1 = no veg, produce dust; 0 = vegation, do not produce dust

         vegmask(1,1) = lai_vegmask(i,j)
      ELSE

         ! Default choice = static ginoux vegmask

         IF (erod(i,j,1) .eq. 0) THEN
            vegmask(1,1) = 0.0
         ELSE
            vegmask(1,1) = 1.0
         ENDIF
      ENDIF

      ! Remove vegetated areas (vegmask=0) from total erodibility.

      erodtot(1,1) = erodtot(1,1) * vegmask(1,1)

      ! Option to use an optional high resolution soil type database from NGA.
      ! Option 0 = WRF (default); Option 1 = NGA
      ! Note NGA dataset is currently only available over southwest Asia
      ! Until a global dataset becomes available, option 0 is recommended
      ! for consistency.  If option 1 is chosen for a domain outside of
      ! SWA, this logic will fill in the areas missing from NGA with the
      ! defaults from WRF. It will work, but it will be inconsistent. GAC

      IF (config_flags%dust_soils .eq. 1) then
         IF (clay_nga(i,j) .ge.0) then
             clay(i,j) = clay_nga(i,j)
             sand(i,j) = sand_nga(i,j)
         ELSE
             clay(i,j) =clay_wrf(i,j)
             sand(i,j) =sand_wrf(i,j)
         ENDIF
      ELSE
         clay(i,j) =clay_wrf(i,j)
         sand(i,j) =sand_wrf(i,j)
      ENDIF

      ! Mass fractions of clay, silt, and sand.

      massfrac(1)=clay(i,j)
      massfrac(2)=1-(clay(i,j)+sand(i,j))
      massfrac(3)=sand(i,j)

      ! Don't allow roughness lengths greater than 20 cm to be lofted.
      ! This kludge accounts for land use types like urban areas and
      ! forests which would otherwise show up as high dust emitters.
      ! This is a placeholder for a more widely accepted kludge
      ! factor in the literature, which reduces lofting for rough areas.
      ! Forthcoming...

      IF (znt(i,j) .gt. 0.2) then
        ilwi(1,1)=0
      ENDIF

      ! Do not allow areas with bedrock, lava, or land-ice to loft.

      IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. &
          isltyp(i,j) .eq. 18) then
        ilwi(1,1)=0
      ENDIF

      ! Another hack to ensure dust does not loft from areas with snow
      ! cover...because, well, that doesn't make sense.

      IF (snowh(i,j) .gt. 0.01) then
        ilwi(1,1)=0
      ENDIF

      ! Volumetric soil moisture can be tuned here with a dust_smtune
      ! set in the namelist.

      volsm=max(smois(i,1,j)*smtune,0.)

      ! Calculate gravimetric soil moisture.

      gravsm(1,1)=100*volsm/((1.-porosity(isltyp(i,j)))*(2.65*(1-clay(i,j))+2.50*clay(i,j)))

      ! Choose an LSM option and drylimit option.
      ! Drylimit calculations based on look-up table in
      ! Clapp and Hornberger (1978) for RUC and PXLSM and
      ! Cosby et al. (1984) for Noah LSM.

      smois_opt = 0
      IF (config_flags%dust_smois == 1) then
         sfc_select: SELECT CASE(config_flags%sf_surface_physics)
            CASE ( RUCLSMSCHEME, PXLSMSCHEME )
               drylimit(1,1) =0.035*(13.52*clay(i,j)+3.53)**2.68
               smois_opt = 1
            CASE ( LSMSCHEME )
               drylimit(1,1) =0.0756*(15.127*clay(i,j)+3.09)**2.3211
               smois_opt = 1
            CASE DEFAULT

               ! Don't currently support volumetric soil moisture
               ! for this scheme, use drylimit based on gravimetric

               drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
         END SELECT sfc_select
      ELSE

         ! use drylimit based on gravimetric soil moisture

         drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
      END IF

      ! Call dust emission routine.

      call source_dust(imx, jmx, lmx, nmx, smx, dt, tc, ustar, massfrac, &
                       erodtot, ilwi, dxy, gravsm, volsm, airden, airmas, &
                       bems, ustart, g, drylimit, alpha, gamma, smois_opt)

      IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
        dustin(i,j,1:5)=tc(1:5)*converi
      ELSE IF(config_flags%chem_opt == CB05_SORG_AQ_KPP .or.  &
              config_flags%chem_opt == CB05_SORG_VBS_AQ_KPP) then
!KW            chem(i,kts,j,p_p25i)=chem(i,kts,j,p_p25i) &
!                         +.25*(tc(1)+.286*tc(2))*converi
!            chem(i,kts,j,p_p25i)=max(chem(i,kts,j,p_p25i),1.e-16)
!            chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) &
!                         +.75*(tc(1)+.286*tc(2))*converi
!            chem(i,kts,j,p_p25j)=max(chem(i,kts,j,p_p25j),1.e-16)
!            chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) &
!                         +(.714*tc(2)+tc(3))*converi
!            chem(i,kts,j,p_soila)=max(chem(i,kts,j,p_soila),1.e-16)
         chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) + 0.03*sum(tc(1:5))*converi
         chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) + 0.97*1.02*sum(tc(1:5))*converi

         chem(i,kts,j,p_ac0) =  chem(i,kts,j,p_ac0) +  0.03*sum(tc(1:5))*converi*factnuma*soilfac
         chem(i,kts,j,p_corn) =  chem(i,kts,j,p_corn) + 0.97*1.02*sum(tc(1:5))*converi*factnumc*soilfac
      ELSE
        chem(i,kts,j,p_dust_1)=tc(1)*converi
        chem(i,kts,j,p_dust_2)=tc(2)*converi
        chem(i,kts,j,p_dust_3)=tc(3)*converi
        chem(i,kts,j,p_dust_4)=tc(4)*converi
        chem(i,kts,j,p_dust_5)=tc(5)*converi
      ENDIF

      ! Diagnostic dust lofting potential diagnostic calculation

      psi=0.
      w10=(u10(i,j)**2.+v10(i,j)**2.)**0.5
      IF (ustar(1,1) .ne. 0. .and. znt(i,j) .ne. 0.) THEN
        psi=0.4*w10/ustar(1,1)-LOG(10.0/znt(i,j))
      ENDIF
      IF (erodtot(1,1) .gt. 0.) then
        afwa_dustloft(i,j)=ustune*w10-ustart*(LOG(10.0/znt(i,j))+psi)/0.4
      ENDIF

      ! For output diagnostics (g m^-2 s^-1)

      emis_dust(i,1,j,p_edust1)=bems(1)
      emis_dust(i,1,j,p_edust2)=bems(2)
      emis_dust(i,1,j,p_edust3)=bems(3)
      emis_dust(i,1,j,p_edust4)=bems(4)
      emis_dust(i,1,j,p_edust5)=bems(5)

      ! Diagnostic total emitted dust (g m^-2 s^-1)

      tot_edust(i,j)=(bems(1)+bems(2)+bems(3)+bems(4)+bems(5))

    ENDIF

    ! Cumulative dust concentration (ug m^-3) and visibility (m) diagnostics
    ! Note visibility is capped at 20 km.  Simple visibility algorithm based
    ! on mean particle diameter for each dust bin - perfect spheres - perfect
    ! blackbodies.

    DO k=kts,kte
      tot_dust(i,k,j)=(chem(i,k,j,p_dust_1)+chem(i,k,j,p_dust_2)+       &
                       chem(i,k,j,p_dust_3)+chem(i,k,j,p_dust_4)+       &
                       chem(i,k,j,p_dust_5))*rho_phy(i,k,j)
      IF ( tot_dust(i,k,j) .gt. 0. ) THEN
        vis_dust(i,k,j)=MIN(3.912/                                      &
                          ((1.470E-6*chem(i,k,j,p_dust_1)+              &
                            7.877E-7*chem(i,k,j,p_dust_2)+              &
                            4.623E-7*chem(i,k,j,p_dust_3)+              &
                            2.429E-7*chem(i,k,j,p_dust_4)+              &
                            1.387E-7*chem(i,k,j,p_dust_5))*             &
                           rho_phy(i,k,j)),999999.)
      ELSE
        vis_dust(i,k,j)=999999.
      ENDIF
    ENDDO

  ENDDO
  ENDDO

  END SUBROUTINE gocart_dust_afwa_driver

  SUBROUTINE source_dust(imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac,&
                         erod, ilwi, dxy, gravsm, volsm, airden, airmas, &
                         bems, ustart, g0, drylimit, alpha, gamma, smois_opt)

  ! ****************************************************************************
  ! *  Evaluate the source of each dust particles size bin by soil emission
  ! *
  ! *  Input:
  ! *         EROD      Fraction of erodible grid cell                (-)
  ! *         ILWI      Land/water flag                               (-)
  ! *         GRAVSM    Gravimetric soil moisture                     (g/g)
  ! *         VOLSM     Volumetric soil moisture                      (g/g)
  ! *         SOILM_OPT Soil moisture option (1:Use GRAVSM 2:VOLSM)   (-)
  ! *         DRYLIMIT  Upper GRAVSM (VOLSM) limit for air-dry soil   (g/g)
  ! *         ALPHA     Constant to fudge the total emission of dust  (1/m)
  ! *         GAMMA     Exponential tuning constant for erodibility   (-)
  ! *         DXY       Surface of each grid cell                     (m2)
  ! *         AIRMAS    Mass of air for each grid box                 (kg)
  ! *         AIRDEN    Density of air for each grid box              (kg/m3)
  ! *         USTAR     Friction velocity                             (m/s)
  ! *         MASSFRAC  Fraction of mass in each of 3 soil classes    (-)
  ! *         DT1       Time step                                     (s)
  ! *         NMX       Number of dust bins                           (-)
  ! *         SMX       Number of saltation bins                      (-)
  ! *         IMX       Number of I points                            (-)
  ! *         JMX       Number of J points                            (-)
  ! *         LMX       Number of L points                            (-)
  ! *
  ! *  Data (see module_data_gocart_dust):
  ! *         SPOINT    Pointer to 3 soil classes                     (-)
  ! *         DEN_DUST  Dust density                                  (kg/m3)
  ! *         DEN_SALT  Saltation particle density                    (kg/m3)
  ! *         REFF_SALT Reference saltation particle diameter         (m)
  ! *         REFF_DUST Reference dust particle diameter              (m)
  ! *         LO_DUST   Lower diameter limits for dust bins           (m)
  ! *         UP_DUST   Upper diameter limits for dust bins           (m)
  ! *         FRAC_SALT Soil class mass fraction for saltation bins   (-)
  ! *
  ! *  Parameters:
  ! *         BETAMAX   Maximum sandblasting mass efficiency          (-)
  ! *         CMB       Constant of proportionality                   (-)
  ! *         MMD_DUST  Mass median diameter of dust                  (m)
  ! *         GSD_DUST  Geometric standard deviation of dust          (-)
  ! *         LAMBDA    Side crack propogation length                 (m)
  ! *         CV        Normalization constant                        (-)
  ! *         G0        Gravitational acceleration                    (m/s2)
  ! *         G         Gravitational acceleration in cgs             (cm/s2)
  ! *
  ! *  Working:
  ! *         BETA      Sandblasting mass efficiency                  (-)
  ! *         U_TS0     "Dry" threshold friction velocity             (m/s)
  ! *         U_TS      Moisture-adjusted threshold friction velocity (m/s)
  ! *         RHOA      Density of air in cgs                         (g/cm3)
  ! *         DEN       Dust density in cgs                           (g/cm3)
  ! *         DIAM      Dust diameter in cgs                          (cm)
  ! *         DMASS     Saltation mass distribution                   (-)
  ! *         DSURFACE  Saltation surface area per unit mass          (m2/kg)
  ! *         DS_REL    Saltation surface area distribution           (-)
  ! *         SALT      Saltation flux                                (kg/m/s)
  ! *         DLNDP     Dust bin width                                (-)
  ! *         EMIT      Total vertical mass flux                      (kg/m2/s)
  ! *         EMIT_VOL  Total vertical volume flux                    (m/s)
  ! *         DSRC      Mass of emitted dust               (kg/timestep/cell)
  ! *
  ! *  Output:
  ! *         TC        Total concentration of dust        (kg/kg/timestep/cell)
  ! *         BEMS      Source of each dust type           (kg/timestep/cell)
  ! *         USTART    Threshold friction vel. (bin 7)               (m/s)
  ! *
  ! ****************************************************************************

  INTEGER, INTENT(IN)   :: nmx,imx,jmx,lmx,smx
  INTEGER, INTENT(IN)   :: ilwi(imx,jmx)
  REAL*8, INTENT(IN)    :: erod(imx,jmx)
  REAL*8, INTENT(IN)    :: ustar(imx,jmx)
  REAL*8, INTENT(IN)    :: gravsm(imx,jmx)
  REAL*8, INTENT(IN)    :: drylimit(imx,jmx)
  REAL*8, INTENT(IN)    :: dxy(jmx)
  REAL*8, INTENT(IN)    :: airden(imx,jmx,lmx), airmas(imx,jmx,lmx)
  REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
  REAL*8, INTENT(OUT)   :: bems(imx,jmx,nmx)
  REAL, INTENT(IN)    :: g0,dt1
  REAL, INTENT(OUT)   :: ustart
  INTEGER, INTENT(IN) :: smois_opt
  REAL*8, INTENT(IN)    :: volsm

  REAL*8    :: den(smx), diam(smx)
  REAL*8    :: dvol(nmx), dlndp(nmx)
! REAL*8    :: distr_dust(nmx)
  REAL*8    :: dsurface(smx), ds_rel(smx)
  REAL*8    :: massfrac(3)
  REAL*8    :: u_ts0, u_ts, dsrc, srce, dmass, dvol_tot
  REAL*8    :: emit, emit_vol
  REAL      :: rhoa, g
  REAL*8    :: salt, stotal
  INTEGER   :: i, j, m, n, s

! Global tuning constant, alpha.  Sandblasting mass efficiency, beta.
! Beta maxes out for clay fractions above 0.2 => betamax.

  REAL, INTENT(IN)  :: alpha
  REAL, PARAMETER :: betamax=5.25E-4
  REAL*8 :: beta

! Experimental optional exponential tuning constant for erodibility.
! 0 < gamma < 1 -> more relative impact by low erodibility regions.
! 1 < gamma -> accentuates high erodibility regions.  Recommend this
! be set to 1 (default) unless looking for a way to increase spread
! within an ensemble framework.

  REAL, INTENT(IN) :: gamma

! Constant of proportionality from Marticorena et al, 1997 (unitless)
! Arguably more ~consistent~ fudge than alpha, which has many walnuts
! sprinkled throughout the literature.

  REAL, PARAMETER :: cmb=1.0    ! Marticorena et al,1997
! REAL, PARAMETER :: cmb=2.61   ! White,1979

! Parameters used in Kok distribution function. Advise not to play with
! these without the expressed written consent of someone who knows what
! they're doing. (See Kok, 2010 PNAS for details on this scheme).

! REAL, PARAMETER :: mmd_dust=3.4D-6  ! median mass diameter (m)
! REAL, PARAMETER :: gsd_dust=3.0     ! geometric standard deviation
! REAL, PARAMETER :: lambda=12.0D-6   ! crack propogation length (m)
! REAL, PARAMETER :: cv=12.62D-6      ! normalization constant

! Calculate saltation surface area distribution from sand, silt, and clay
! mass fractions and saltation bin fraction. This will later become a
! modifier to the total saltation flux.  The reasoning here is that the
! size and availability of saltators affects saltation efficiency. Based
! on Eqn. (32) in Marticorena & Bergametti, 1995 (hereon, MB95).

  DO n=1,smx
    dmass=massfrac(spoint(n))*frac_salt(n)
    dsurface(n)=0.75*dmass/(den_salt(n)*reff_salt(n))
  ENDDO

! The following equation yields relative surface area fraction.  It will only
! work if you are representing the "full range" of all three soil classes.
! For this reason alone, we have incorporated particle sizes that encompass
! the clay class, to account for its relative area over the basal
! surface, even though these smaller bins would be unlikely to play any large
! role in the actual saltation process.

  stotal=SUM(dsurface(:))
  DO n=1,smx
    ds_rel(n)=dsurface(n)/stotal
  ENDDO

! Calculate total dust emission due to saltation of sand sized particles.
! Begin by calculating DRY threshold friction velocity (u_ts0).  Next adjust
! u_ts0 for moisture to get threshold friction velocity (u_ts). Then
! calculate saltation flux (salt) where ustar has exceeded u_ts.  Finally,
! calculate total dust emission (tot_emit), taking into account erodibility.

  g = g0*1.0E2                          ! (cm s^-2)
  emit=0.0

  DO n = 1, smx                         ! Loop over saltation bins
    den(n) = den_salt(n)*1.0D-3         ! (g cm^-3)
    diam(n) = 2.0*reff_salt(n)*1.0D2    ! (cm)
    DO i = 1,imx
      DO j = 1,jmx
        rhoa = airden(i,j,1)*1.0D-3     ! (g cm^-3)

        ! Threshold friction velocity as a function of the dust density and
        ! diameter from Bagnold (1941) (m s^-1).

        u_ts0 = 0.13*1.0D-2*SQRT(den(n)*g*diam(n)/rhoa)* &
                SQRT(1.0+0.006/den(n)/g/(diam(n))**2.5)/ &
                SQRT(1.928*(1331.0*(diam(n))**1.56+0.38)**0.092-1.0)

        ! Friction velocity threshold correction function based on physical
        ! properties related to moisture tension. Soil moisture greater than
        ! dry limit serves to increase threshold friction velocity (making
        ! it more difficult to loft dust). When soil moisture has not reached
        ! dry limit, treat as dry (no correction to threshold friction
        ! velocity). GAC

        ! Calculate threshold friction velocity. If volumetric (gravimetric)
        ! water content is less than the drylimit, then the threshold friction
        ! velocity (u_ts) will be equal to the dry threshold friction velocity
        ! (u_ts0). EDH

        IF (smois_opt .EQ. 1) THEN
          IF (100.*volsm > drylimit(i,j)) THEN
            u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*((100.*volsm)-drylimit(i,j))**0.68)))
          ELSE
            u_ts = u_ts0
          ENDIF
        ELSE
          IF (gravsm(i,j) > drylimit(i,j)) THEN
            u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*(gravsm(i,j)-drylimit(i,j))**0.68)))
          ELSE
            u_ts = u_ts0
          END IF
        END IF

        ! Bin 7 threshold friction velocity for diagnostic dust lofting
        ! potential product

        IF (n .eq. 7) THEN  ! Saltation bin 7 is small sand
          ustart = u_ts
        ENDIF

        ! Saltation flux (kg m^-1 s^-1) from MB95
        ! ds_rel is the relative surface area distribution

        IF (ustar(i,j) .gt. u_ts .and. erod(i,j) .gt. 0.0 .and. ilwi(i,j) == 1) THEN
          salt = cmb*ds_rel(n)*(airden(i,j,1)/g0)*(ustar(i,j)**3)* &
                 (1. + u_ts/ustar(i,j))*(1. - (u_ts**2)/(ustar(i,j)**2))
        ELSE
          salt = 0.D0
        ENDIF

        ! Calculate total vertical mass flux (note beta has units of m^-1)
        ! Beta acts to tone down dust in areas with so few dust-sized particles
        ! that the lofting efficiency decreases.  Otherwise, super sandy zones
        ! would be huge dust producers, which is generally not the case.
        ! Equation derived from wind-tunnel experiments (see MB95).

        beta=10**(13.6*massfrac(1)-6.0)  ! (unitless)
        IF (beta .gt. betamax) THEN
          beta=betamax
        ENDIF
        emit=emit+salt*(erod(i,j)**gamma)*alpha*beta    ! (kg m^-2 s^-1)
      END DO
    END DO
  END DO                                 ! End do over saltation bins

  ! Now that we have the total dust emission, distribute into dust bins using
  ! lognormal distribution (Dr. Jasper Kok, 2010), and
  ! calculate total mass emitted over the grid box over the timestep.
  !
  ! In calculating the Kok distribution, we assume upper and lower limits to
  ! each bin. For reff_dust=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/) (default),
  ! lower limits were ASSUMED at lo_dust=(/0.1D-6,1.0D-6,1.8D-6,3.0D-6,6.0D-6/)
  ! upper limits were ASSUMED at up_dust=(/1.0D-6,1.8D-6,3.0D-6,6.0D-6,10.0D-6/)
  ! These may be changed within module_data_gocart_dust.F, but make sure it is
  ! consistent with reff_dust values.  These values were taken from the original
  ! GOCART bin configuration. We use them here to calculate dust bin width,
  ! dlndp. dvol is the volume distribution. GAC
  !
  ! UPDATE: We bypass the calculation below and instead hardcode distr_dust for
  ! the five dust bins we are using here since this distribution is static and
  ! unnecessary to calculate at every time step. Keeping everything here to
  ! document the steps in obtaining distr_dust.  GAC 20140320

!  dvol_tot=0.
!  DO n=1,nmx  ! Loop over all dust bins
!    dlndp(n)=LOG(up_dust(n)/lo_dust(n))
!    dvol(n)=(2.0*reff_dust(n)/cv)*(1.+ERF(LOG(2.0*reff_dust(n)/mmd_dust)/(SQRT(2.)*LOG(gsd_dust))))*&
!          EXP(-(2.0*reff_dust(n)/lambda)**3.0)*dlndp(n)
!    dvol_tot=dvol_tot+dvol(n)
!
!   ! Convert mass flux to volume flux
!    emit_vol=emit/den_dust(n) ! (m s^-1)
!  END DO
!  DO n=1,nmx  ! Loop over all dust bins
!    distr_dust(n)=dvol(n)/dvol_tot
!  END DO

  ! Now distribute total vertical emission into dust bins and update
  ! concentration.

  DO n=1,nmx  ! Loop over all dust bins
    DO i=1,imx
      DO j=1,jmx

        ! Calculate total mass emitted

        dsrc = emit*distr_dust(n)*dxy(j)*dt1  ! (kg)
        IF (dsrc < 0.0) dsrc = 0.0

        ! Update dust mixing ratio at first model level.

        tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1) ! (kg/kg)
        bems(i,j,n) = 1000.*dsrc/(dxy(j)*dt1) ! diagnostic (g/m2/s)
      END DO
    END DO
  END DO

  END SUBROUTINE source_dust

END MODULE GOCART_DUST_AFWA