MODULE qf03
!
! Y. Shao, 29 Jan 2004 
!         
! JY Kang, 01 Dec 2008
! Modify the code for WRF_chem
!
! M. Klose, 2010-2015 Modifications
!-----------------------------------------------------------------------------------
! Calculate sediment flux for multi-particle size soils as a weighted average of Q(d) 
!           dust emission F(d) for covered and moisture soil
!
! Options for dust calculation:
! 1 Shao (2001)
! 2 Shao (2004)
! 3 Shao (2011): simplification of 2; added on 26 Sep 2009
!
!--------------------------------------------------------------------------------------
!
!      input: 
!      n:          number of particle size ranges.
!      dm:         median diameter of each particle size.                   [m]
!      m_fract:    Weight fraction of each particle range. Sum m_fract = 1
!      ustar:      Friciton velocity.                                       [m/s]
!      cf:         fraction area covered by roughness elements/vegetation cover
!      w:          surface soil moisture content [m^3/m^3]
!      wr:         air-dry soil moisture [m^3/m^3] from SOILPARM.TBL
!      c:          Owen's coefficient
!
!      output:  
!      ustart:     Mean threshold velocity of each particle range.
!      q:          Sand flux from each size range.
!      ffq:        Weighted sand flux from each size range.
!      qtotal:     Total sand flux (weighted average).
!      f:          Weighted dust flux ejected by single sand size range for given d_d. 
!      fff:        Weighted dust flux ejected by all sand size range from each dust size range.
!      ftotal:     Total dust flux (weighted average).
!--------------------------------------------------------------------------------------
!

  CONTAINS

  subroutine qf03_driver ( nmx, idst, g, rhop, rho, dt, &
                    ustar, w, cf, ust_min, imod, dz_lowest, &
                    soilc, tot_soilc, domsoilc,    &
                    tc, bems, rough_cor_in, smois_cor_in, wr, &
                    d0, dd, psd_m, dpsd_m, ppsd_m, psd_f, dpsd_f, ppsd_f, imax, stype)
       
       INTEGER, INTENT(IN   )                   :: nmx, imod, idst
       REAL   , INTENT(IN   )                   :: g, dt, rho, dz_lowest, ustar
       REAL(8), INTENT(IN   )                   :: w
       REAL(8), INTENT(INOUT)                   :: cf
       REAL(8), INTENT(  OUT)                   :: ust_min, rough_cor_in, smois_cor_in
       REAL   , INTENT(INOUT), DIMENSION(nmx)   :: tc
       REAL   , INTENT(  OUT), DIMENSION(nmx)   :: bems
       REAL   , INTENT(IN   ), DIMENSION(12)    :: wr

!kang [2009/01/07] soil class type
       REAL, INTENT(IN   )                      :: tot_soilc       
       REAL, DIMENSION(16), INTENT(IN)          :: soilc       
       INTEGER,  INTENT(IN  )                   :: domsoilc
       
       integer :: imax, stype
       real(8), dimension(0:imax), intent(in) :: d0
       real(8), dimension(imax), intent(in)   :: dd 
       real(8), dimension(imax,stype), intent(in)   :: psd_m, dpsd_m, ppsd_m
       real(8), dimension(imax,stype), intent(in)   :: psd_f, dpsd_f, ppsd_f
 
!local variables
!      particle-size distributions 
       real(8), dimension(imax)      :: psdm, dpsdm, ppsdm    ! minimally dispersed
       real(8), dimension(imax)      :: psdf, dpsdf, ppsdf    ! fully dispersed
       real(8), dimension(imax)      :: psds, dpsds, ppsds    ! sediment       

       real(8) :: smois_correc
       integer :: i, j, n, kk, ij, ffile
       real(8) :: total, qtotal, ftotal
       real(8) :: ftotalb, ftotalp
!
       real(8), dimension(imax)   :: beta_d, beta_s  ! beta1 and beta2 used by Shao et al. (1996) dust emission
       real(8), dimension(imax)   :: ustart, q, ffq, fff
       real(8), dimension(imax,imax) :: f
!
       real(8), parameter :: c_lambda=0.35
       real(8) :: h, lambda
       real(8) :: ghl, fc
       real(8) :: phl
       real(8) :: cys, u0, al0, rhos, smass, omega, rys
       real(8) :: ddm, a1, a2, a3  
       real(8) :: zeta, sigma_m                      ! u*sqrt(rhos/p), bombardment coefficient
!
       real(8) :: ustart0_out, qwhite_out, f_mb_out, f_hlys_out, pmass_out, vhlys_out
! 
       integer, parameter :: nbins = 5 
       real(8) :: sigma
       real, dimension(nbins)    :: dbin, fbin, cell_fbin
       integer, dimension(nbins) :: ibin
       data dbin/2.,3.6,6.,12.,20./ !size cut diameter (um) consistent with GOCART model
       real(8), intent(in) :: rhop
       real    :: cell_ftotal
       integer :: isl, cc
       character :: msg, soil
!*******************************************************************************************
!
! initialization
       cell_ftotal = 0.
       do n = 1, nbins
       cell_fbin(n) = 0.
       enddo
! 
       DO cc = 1, 12  ! soil category
           if (soilc(cc).eq.0.) then
              go to 103
           endif
!        
           psdm(:)=psd_m(:,cc)
           psdf(:)=psd_f(:,cc)
           dpsdm(:)=dpsd_m(:,cc)
           dpsdf(:)=dpsd_f(:,cc)
           ppsdm(:)=ppsd_m(:,cc)
           ppsdf(:)=ppsd_f(:,cc)

! default settings: 
       rhos = 1000.d0      ! bulk density of soil [kg/m3] 
       phl = 30000.        ! plastic pressure [N/m2]
       cys = 0.00001       ! cys : parameter

       sigma = rhop/rho    ! particle-air density ratio

!
!-------------------
! frontal area index
!-------------------
         if (cf .gt. 0.95) then 
            write(6,*) 'cover fraction too large, reset'
            cf = 0.95
         endif! as larger values lead to large lambda, which is problematic in Raupach scheme
         
         lambda = - c_lambda*dlog( 1.d0 - cf )
         call r_c(lambda, rough_cor_in)

! Matching WRF_soil class with Shao's class for moisture correction of Fecan (cc:WRF_soil class, isl:Shao's class)
! cc
! 1:sand, 2:loamy sand, 3:sandy loam, 4:silt loam, 5:silt, 6:loam, 7:sandy clay loam,
! 8:silty clay loam, 9:clay loam, 10:sandy clay, 11:silty clay, 12:clay
! isl
! 1:sand, 2:loamy sand, 3:sandy loam, 4:loam, 5:silt loam, 6: silt, 7:sandy clay loam, 
! 8:clay loam, 9:silty clay loam, 10:sandy clay, 11:silty clay, 12:clay
        if (cc.eq.1.or.cc.eq.2.or.cc.eq.3.or.cc.eq.7.or.cc.eq.10.or.  & 
        &   cc.eq.11.or.cc.eq.12) then
         isl = cc
        elseif (cc.eq.4) then
         isl = 5
        elseif (cc.eq.5) then
         isl = 6
        elseif (cc.eq.6) then
         isl = 4
        elseif (cc.eq.8) then
         isl = 9
        elseif (cc.eq.9) then
         isl = 8
        endif

         call h_c(w, wr(cc), isl, smois_correc)
!----------------------------------------------
! for each particle size group, estimate ustart
!----------------------------------------------
!
         ust_min = 999.0
         do i = 1, imax
            call ustart0(dd(i), sigma, g, rho, ustart0_out)
            ustart(i) = ustart0_out
            ustart(i) = rough_cor_in*smois_correc*ustart(i)
            ust_min = dmin1(ust_min, ustart(i))
            call qwhite(ustart(i), ustar, rho, g, qwhite_out)
            q(i) = qwhite_out
            q(i) = (1.d0-cf)*q(i)
         enddo
         if (cc.eq.domsoilc) then
          smois_cor_in = smois_correc
         endif
! 
! 
         IF ( ustar .le. ust_min ) THEN    ! no erosion goto 102
            q = 0.d0
            ffq = 0.d0
            qtotal = 0.d0
            fff = 0.d0
            ftotal = 0.d0
            fbin = 0.d0
            goto 102
         ELSE
            ghl = dexp( -(ustar - ust_min)**3.d0 )
            dpsds = ghl*dpsdm + (1.-ghl)*dpsdf
            psds  = ghl*psdm + (1.-ghl)*psdf
            ppsds = ghl*ppsdm + (1.-ghl)*ppsdf
!
            ffq = q*dpsds
            qtotal = sum(ffq)

!--------------
! dust emission
!--------------
! 
! size bin
            do n=1,nbins
             ibin(n)=0
             do i=imax,1,-1
              if(d0(i).ge.dbin(n)) ibin(n)=i
             enddo
             if(ibin(n).eq.0) stop 'wrong dust classes'
            enddo
! 
! 
! 
!--------------------------------
! Shao (2001) dust emission model
!--------------------------------
        IF (imod .eq. 1) THEN 
           do i = idst+1, imax
              ddm = dd(i)*1.d-6
              call pmass(rhop, ddm, pmass_out)                       ! mass of saltating particles
              smass = pmass_out
              u0 = 10*ustar
              al0 = 13.d0*3.14159d0/180.d0
              call vhlys(phl, 2, smass, al0, u0, ddm, vhlys_out)
              omega = vhlys_out     ![m3]
!               
              do j = 1, idst
                 rys = psdm(j)/psdf(j)
                 if (rys.gt.1.) then 
                     rys = 1.
                 endif
                 a1 = cys*( (1.-ghl) + ghl*rys )
                 a2 = ffq(i)*g/ustar**2/smass
                 
                 if ( dpsdf(j) .lt. dpsdm(j) ) then 
                    a3 = dpsdf(j)*rhos*omega
                 else 
                    a3 = dpsdf(j)*rhos*omega + (dpsdf(j)-dpsdm(j))*smass
                 endif
                 f(i,j) = a1*a2*a3    ![kg/m2/s]
              enddo
           enddo
!          
           ftotal = 0.0
           do j = 1, idst
              fff(j) = 0.
              do i = idst+1, imax
                 fff(j) = fff(j) + f(i,j)
              enddo
              fff(j) = (1.d0-cf)*fff(j)
              ftotal = ftotal + fff(j)
           enddo
! 
           do n=1,nbins
            j0=1
            if(n.gt.1) j0=ibin(n-1)+1
            fbin(n)=0
            do j=j0,ibin(n)
             fbin(n)=fbin(n)+fff(j)
            enddo
           enddo
! 
!--------------------------------
! Shao (2004) dust emission model
!--------------------------------
        ELSEIF (imod .eq. 2) THEN
          zeta = ustar*dsqrt( rhos/phl )
          sigma_m = 12.d0*zeta*zeta*(1.d0+14.d0*zeta)
! 
          do i = idst+1, imax
             do j = 1, idst
                rys = psdm(j)/psdf(j)
                if (rys .gt.1.) then
                    rys = 1.
                endif
                a1 = cys*dpsdf(j)*( (1.-ghl) + ghl*rys )
                a2 = (1.+sigma_m)
                a3 = ffq(i)*g/ustar**2
                f(i,j) = a1*a2*a3
             enddo
          enddo
!           
          ftotal = 0.0
          do j = 1, idst
             fff(j) = 0.
             do i = idst+1, imax
                fff(j) = fff(j) + f(i,j)
             enddo
             fff(j) = (1.d0-cf)*fff(j)             
             ftotal = ftotal + fff(j)
          enddo
!           
          do n=1,nbins
           j0=1
           if(n.gt.1) j0=ibin(n-1)+1
           fbin(n)=0
           do j=j0,ibin(n)
            fbin(n)=fbin(n)+fff(j) 
           enddo
          enddo
! 
! 
!--------------------------------------------------------------------------
! Shao (2011) minimal version, ghl = 1, Q independent of sand particle size
! 
! See Eq. (34) in 
! Shao, Y., M. Ishizuka, M. Mikami, J. Leys (2011): Parameterization of size-
! resolved dust emission and validation with measurements, JGR, 116, D08203, 
! doi: 10.1029/2010JD014527
!--------------------------------------------------------------------------
        ELSEIF (imod .eq. 3) THEN 
          zeta = ustar*dsqrt( rhos/phl )
          sigma_m = 12.d0*zeta*zeta*(1.d0+14.d0*zeta)
! 
          ftotal = 0.0          
! 
          do j = 1, idst
             a1 = cys*dpsdm(j)
             a2 = (1.+sigma_m)
             a3 = qtotal*g/ustar**2
             fff(j) = a1*a2*a3
             fff(j) = (1.d0-cf)*fff(j)
             ftotal = ftotal + fff(j)
          enddo
!           
          do n=1,nbins
           j0=1
           if(n.gt.1) j0=ibin(n-1)+1
           fbin(n)=0
           do j=j0,ibin(n)
            fbin(n)=fbin(n)+fff(j) 
           enddo
          enddo
! 
        ENDIF  ! dust scheme
      ENDIF    ! ust < ust_t
! 
! 
! 
102   continue

      do n = 1, nbins
      cell_fbin(n) = cell_fbin(n) + (soilc(cc)/tot_soilc)*fbin(n)
      enddo

103   continue  

      ENDDO  ! cc, soil category
      
      do n = 1, nbins
! fbin : [kg/m2/s],  dz_lowest : [m],  rho : [kg/m3],  dt : [s] -> tc : [kg/kg-dryair]
         tc(n) = tc(n) + cell_fbin(n)/dz_lowest/rho*dt     ![kg/kg-dryair]
         bems(n) = cell_fbin(n)     ![kg/m2/s]
      enddo


      END subroutine qf03_driver


!*****************************************************************************
      subroutine ustart0(dum, sigma, g, rho, ustart0_out)
!
! Y. Shao, 13 June 2000
!
! Calculate ustar0(d) using Shao and Lu (2000) for uncovered
! dry surface
!
! dum:    particle diameter			[um]
! ustar0: threshold friction velocity   	[m/s]
!
      real, intent(in)     :: g, rho
      real(8), intent(in)  :: dum, sigma
      real(8), intent(out) :: ustart0_out
      real(8) :: dm
      real(8), parameter :: gamma = 1.65d-4      ! a constant
      real(8), parameter :: f = 0.0123

      dm = dum*1d-6

      ustart0_out = f*(sigma*g*dm + gamma/(rho*dm) )
      ustart0_out = dsqrt( ustart0_out )
!      end function
      end subroutine ustart0
!*****************************************************************************
      subroutine qwhite(ust, ustar, rho, g, qwhite_out)
!
! Yaping Shao 17-07-99!
!
! White (1979) Sand Flux Equation
! Q = c*rho*u_*^3 over g (1 - u_*t over u_*)(1 + u_*t^2/u_*^2)
! qwhite: Streamwise Sand Flux;       [kg m-1 s-1]
! c     : 2.6
! ust   : threhold friction velocity  [m/s]
! ustar : friction velocity           [m/s]
!
      real(8) :: c
      real, intent(in) :: ustar, rho, g
      real(8), intent(in) :: ust
      real(8), intent(out) :: qwhite_out
      real(8) :: a, b
! default setting:       
      c = 2.3d0
      a = rho/g

      IF (ustar.lt.ust) THEN 
        qwhite_out = 0.d0
      ELSE
        b = ust/ustar 
        qwhite_out = c*a*ustar**3.*(1.-b)*(1.+b)**2.
      ENDIF

      END subroutine qwhite
!*****************************************************************************
      subroutine vhlys(p, k, xm, alpha, u, d, vhlys_out)
!
! Volume removal according to Lu and Shao (1999), Equation (8)
! alpha: impact angle             [^o]
! u    : impact velocity          [m/s]
! p    : plastic pressure         [N/m^2]
! xm   : particle mass            [kg]
! d    : particle diameter        [m]
!
!
      REAL(8),intent(in) :: alpha, xm, u, d, p
      REAL(8) :: beta
      REAL(8), PARAMETER :: pi=3.1415927d0
      REAL(8) :: t1, t2, t3
      INTEGER,intent(in) :: k
      real(8),intent(out) :: vhlys_out

      beta = dsqrt( p*k*d/xm )
      t1 = u*u/(beta*beta)*( dsin(2.d0*alpha) - 4.d0*dsin(alpha)*dsin(alpha) )
      t2 = u*dsin(alpha)/beta
      t3 = 7.5d0*pi*t2**3.d0/d
      vhlys_out = d*( t1 + t3 )
      
      END subroutine vhlys      
!*****************************************************************************
! A routine for correction of ust for soil moisture content
!
! w : volumetric soil moisture
! isl: soil texture type, ranging from 1 to 12
!
! Author: Yaping Shao, 5/05/2001
! Reference: Fecan et al. (1999), Ann. Geophysicae,17,149-157
!
! Data based on Shao and Jung, 2000, unpublished manuscript
! Data invented for sand, loamy sand, sandy loam, loam, clay loam, and clay
! Soil classes:
! 1 sand, 2 loamy sand, 3 sandy loam, 4 loam, 5 silty loam, 6 silt, 7 sandy clay loam, 8 clay loam, 
! 9 silty clay loam, 10 sandy clay, 11 silty clay, 12 clay
!----------------------------------------------------------------------
      subroutine h_c (w, wr, isl, h)

      real(8) :: a(12), b(12)!, thr(12) 
      real(8), intent(in)  :: w
      real(8), intent(out) :: h
      real, intent(in)     :: wr
      integer, intent(in)  :: isl
      character*100 :: msg  ! error message string

!       data thr/0.001, 0.003, 0.037, 0.049, 0.061, 0.072, 0.084, 0.095, 0.110, 0.126, 0.141, 0.156/
      data a /21.19, 33.03, 44.87, 17.79, 20.81, 23.83, 26.84, 29.86, 27.51, 25.17, 22.82, 20.47/
      data b / 0.68,  0.71,  0.85,  0.61,  0.66,  0.71,  0.75,  0.80,  0.75,  0.70,  0.64,  0.59/     
      
      if ( w.lt.0. ) then
        write(msg, *) 'soil moisture correction (h_c): w = ', w, ' <  0'
        call wrf_error_fatal(msg)
!         stop
      endif
      
      if ( w.le.wr ) then
         h = 1.0
      else
         h = sqrt( 1 + a(isl)*( w-wr )**b(isl) )         
      endif

      END subroutine h_c

!*****************************************************************************
      subroutine pmass(rhop, d, pmass_out)
!
! 	Particle Mass	
!     	rhop: particle density		[kg m^-3]
!     	d   : particle size		[m]
!	
      REAL(8), PARAMETER :: pi=3.1415927d0
      REAL(8),intent(in) :: rhop, d
      real(8),intent(out) :: pmass_out

      pmass_out = (pi*rhop*d**3.d0)/6.d0

      END subroutine pmass

!*****************************************************************************
      subroutine r_c (x,r)
!
!   Y. Shao 17-07-92
!   CORRECTION FUNCTION FOR UST(D) BASED ON Raupach et al. (1992)
!   x = frontal area index
!
!   R_C = (1 - sig m x)^{1/2} (1 + m beta x)^{1/2}
!         Note I deife R_C = u_{*tR}/u_{*tS}
!         While Raupach et al. defined 
!                      R_C = u_{*tS}/u_{*tR} and their R function is
!   R_C = (1 - sig m x)^{-1/2} (1 + m beta x)^{-1/2}
!                      
!   sig   : basal to frontal area; about 1
!   m     : parameter less than 1; about 0.5
!   beta  : a ratio of drag coef.; about 90. ; changed to 200 (recommendation by Y. Shao based on data by Gillies et al.)
!
      real(8) :: xc
      real(8), intent(in) :: x
      real(8), intent(out) :: r
      real(8), parameter :: sig=1., m=0.5, beta=200.
!
      xc = 1./(sig*m)
      IF (x.ge.xc) THEN
        r   = 999.           ! Full covered surface
      ELSE
        r   = dsqrt(1.-sig*m*x)*dsqrt(1.+m*beta*x)
      ENDIF
!
      END subroutine r_c



END MODULE qf03