MODULE uoc_dust
!----------------------------------------------------------------------------
! Dust emission module developed at the University of Cologne, Germany. 
! Dust emission schemes and framework developed by Y Shao (yshao@uni-koeln.de)
! Implementation into WRF and modifications by JY Kang (jy.kang@kiaps.org), 
! CL Wu (wuchenglai@mail.iap.ac.cn), and M Klose (mklose@uni-koeln.de).
! 
! For references and available schemes, see module_qf03.F
! Martina Klose, 29 May 2013
! 
!----------------------------------------------------------------------------
  USE module_data_gocart_dust
  USE qf03
  USE module_soilpsd
  USE module_sf_noahlsm, ONLY:DRYSMC
  USE NOAHMP_TABLES, ONLY: DRYSMC_nmp => SMCDRY_TABLE
  USE module_sf_ruclsm,  ONLY:DRYSMC_ruc => DRYSMC
  
  CONTAINS
  subroutine uoc_dust_driver(ktau,dt,config_flags,                         &
         chem,rho_phy,dz8w,smois,ust,                                      &
         isltyp,vegfra,g,emis_dust,                                        &
         ust_t_min, imod, rough_cor, smois_cor,                            &
         soil_top_cat, erod,                                               &
         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_model_constants, ONLY: mwdry
  IMPLICIT NONE
   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags

   INTEGER,      INTENT(IN   ) :: ktau, imod,                              &
                                  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   ) ::                                 isltyp 

   REAL, DIMENSION(ims:ime,1:config_flags%num_soil_cat,jms:jme) ,          &
          INTENT(IN   ) ::                           soil_top_cat   
          
   REAL,  DIMENSION( ims:ime , jms:jme, 3 ) ,                              &
          INTENT(IN   ) ::                           erod 

! ust_t_min is calculated value from qf03
   REAL,  DIMENSION( ims:ime , jms:jme ) ,                                 &
          INTENT(INOUT) ::                           ust_t_min,            &
                                                     rough_cor,            &   
                                                     smois_cor
   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(INOUT) ::                               smois
   REAL,  DIMENSION( ims:ime , jms:jme )  ,                                &
          INTENT(IN   ) ::                            ust,                 &
                                                      vegfra
   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
          INTENT(IN   ) ::                           dz8w,                 &
                                                     rho_phy         
 
  REAL, INTENT(IN   ) :: dt,g
!
! local variables
!
  integer, parameter :: imax=100              ! No. of particle size intervals for psd
  integer, parameter :: jmax=4                ! No. of log-normal distributions for constructing psd
  integer, parameter :: stype=12               ! No. of soil texture classes
  real(8), dimension(0:imax) :: d0
  real(8), dimension(imax)   :: dd 
  real(8), dimension(imax)   :: psdm, dpsdm, ppsdm
  real(8), dimension(imax)   :: psdf, dpsdf, ppsdf
  real(8), dimension(imax,stype)   :: psd_m, dpsd_m, ppsd_m
  real(8), dimension(imax,stype)   :: psd_f, dpsd_f, ppsd_f
  real(8), parameter :: dcut=20.d0            ! dust cutoff particle size
  real(8), parameter :: rhop = 2650.d0        ! particle density [kg/m3]
  
  
  integer             :: i,j,k,p,idst
  integer, parameter  :: nmx = 5  ! No. of dust bins
  real                :: ust_grid, airden, dz_lowest, smc
  real, dimension(nmx):: tc, bems
  real(8)             :: gwet, cf
  real(8)             :: conver,converi
  real(8)             :: ust_min, rough_cor_in, smois_cor_in
  real, dimension(16) :: soilc
  real                :: tot_soilc
  integer             :: domsoilc
  integer             :: cc    
! 
  real(8)             :: ustart0_out, lambda, sigma
  real(8), dimension(imax) :: ustart
  real, dimension(12) :: thr
  data thr/0.001, 0.003, 0.037, 0.061, 0.072, 0.049, 0.084, 0.110, 0.095, 0.126, 0.141, 0.156/
  
!   character*4, dimension(12):: s_type
!   data s_type /'sand', 'losa', 'salo', 'silo', 'silt', 'loam', 'sclo', 'sicl', 'cloa', &
!   &            'sacl', 'silc', 'clay'/
! saturation soil moisture from SOILPARM.TBL
!***************************************************************************
! initialization
   
  conver=1.e-9
  converi=1.e9

  k=kts   !in the bottom layer
  
! calculate soil-psd once for all 12 types
! initialize
   psd_m(:,:) = 0.
   psd_f(:,:) = 0.
   dpsd_m(:,:) = 0.
   dpsd_f(:,:) = 0.
   ppsd_m(:,:) = 0.
   ppsd_f(:,:) = 0.

   do p = 1, stype
!     write(6,*) p
      if (p.eq.1) then            ! sand
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csandm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csandf, jmax)
      elseif (p.eq.2) then        ! loamy sand
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, closam, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, closaf, jmax)
      elseif (p.eq.3) then        ! sandy loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csalom, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csalof, jmax)
      elseif (p.eq.4) then        ! silt loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csilom, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csilof, jmax)
      elseif (p.eq.5) then        ! silt
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csiltm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csiltf, jmax)
      elseif (p.eq.6) then        ! loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, cloamm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, cloamf, jmax)
      elseif (p.eq.7) then        ! sandy clay loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csclom, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csclof, jmax)
      elseif (p.eq.8) then        ! silty clay loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csiclm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csiclf, jmax)
      elseif (p.eq.9) then        ! clay loam
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, ccloam, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, ccloaf, jmax)
      elseif (p.eq.10) then        ! sandy clay
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csaclm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csaclf, jmax)
      elseif (p.eq.11) then        ! silty clay
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, csilcm, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, csilcf, jmax)         
      elseif (p.eq.12) then        ! clay
         call psd_create(d0, dd, psdm, dpsdm, ppsdm, imax, cclaym, jmax)
         call psd_create(d0, dd, psdf, dpsdf, ppsdf, imax, cclayf, jmax)       
      endif
         psd_m(:,p) = psdm(:)
         psd_f(:,p) = psdf(:)
         dpsd_m(:,p) = dpsdm(:)
         dpsd_f(:,p) = dpsdf(:)
         ppsd_m(:,p) = ppsdm(:)
         ppsd_f(:,p) = ppsdf(:)      
   enddo  
!    stop
!    
! Before calculating dust emission, some parameters should be set.
   j = 0
1  j = j+1
   if ( dd(j) .le. dcut ) then
      idst = j
      goto 1
   endif
! 
  do j=jts,jte
  do i=its,ite

! do dust over dust source area only
     if (sum(erod(i,j,:)).gt.0.) then    !use fraction of erodible surface area as dust source indicator    
        tc(1)=chem(i,kts,j,p_dust_1)*conver   ![kg/kg-dryair]
        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

        ust_grid=ust(i,j)   ! u* for one grid

! delta z for the lowest layer is needed for unit conversion
        dz_lowest = dz8w(i,1,j)
 
! for soil moisture using volumetric soil moisture (smois)
        if (smois(i,1,j).gt.1.) then 
           smois(i,1,j) = 1.
           CALL wrf_message('UoC CTDE WARNING: vol. soil moisture > 1, reset')
        endif 

        gwet=smois(i,1,j)
        airden=rho_phy(i,kts,j)                       ![kg/m3]
        cf=vegfra(i,j)                                ! in [%]
        cf = cf/100.d0                                ! fraction 

! initialization
         tot_soilc=0.
         soilc = 0.

         do cc = 1, 12 
            soilc(cc) = soil_top_cat(i,cc,j)
            tot_soilc = tot_soilc + soilc(cc)
         enddo
!
         domsoilc = isltyp(i,j)
         if (  config_flags%sf_surface_physics .eq. 3   ) then
             DRYSMC = DRYSMC_ruc
         elseif (  config_flags%sf_surface_physics .eq. 4   ) then
             DRYSMC = DRYSMC_nmp                       
         elseif ( config_flags%sf_surface_physics .eq. 1 .or. & 
         &        config_flags%sf_surface_physics .eq. 5 .or. & 
         &        config_flags%sf_surface_physics .eq. 7 .or. & 
         &        config_flags%sf_surface_physics .eq. 8 .or. & 
         &        config_flags%sf_surface_physics .eq. 0) then 
             DRYSMC(1:12) = thr             
             CALL wrf_message('UoC dust: DRYSMC reset for dust emission')
         endif           
!------------------------------------------------------------------------
         call qf03_driver( nmx, idst, g, rhop, airden, dt,        &
                  ust_grid, gwet, cf, ust_min, imod, dz_lowest, &
                  soilc, tot_soilc, domsoilc,                   &
                  tc, bems, rough_cor_in, smois_cor_in, DRYSMC(1:12),    & 
                  d0, dd, psd_m, dpsd_m, ppsd_m, psd_f, dpsd_f, ppsd_f, imax, stype)
                  
!                   
         chem(i,kts,j,p_dust_1)=tc(1)*converi   ![ug/kg-dryair]
         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
! for output diagnostics
         emis_dust(i,1,j,p_edust1)=bems(1)*converi 
         emis_dust(i,1,j,p_edust2)=bems(2)*converi
         emis_dust(i,1,j,p_edust3)=bems(3)*converi 
         emis_dust(i,1,j,p_edust4)=bems(4)*converi
         emis_dust(i,1,j,p_edust5)=bems(5)*converi      ![kg/m2/s] -> [ug/m2/s] 

     else   ! no dust source
         emis_dust(i,1,j,p_edust1)=0.
         emis_dust(i,1,j,p_edust2)=0.
         emis_dust(i,1,j,p_edust3)=0.
         emis_dust(i,1,j,p_edust4)=0.
         emis_dust(i,1,j,p_edust5)=0.
!          if (imod .ne. 4) then
            ust_min = -9999.d0
            rough_cor_in = 1.d0
            smois_cor_in = 1.d0
!          endif
     endif !dsource/erod

        ust_t_min(i,j) = ust_min
        rough_cor(i,j) = rough_cor_in
        smois_cor(i,j) = smois_cor_in
  enddo ! i loop
  enddo ! j loop
!

end subroutine UoC_dust_driver


!*****************************************************************************
      subroutine psd_create(d, dm, psd, dpsd, ppsd, imax, cmtrix, jmax)
!
!----------------------------------------------------------------------------
! Yaping Shao, 13 June 2000
!
! - Generate particle size distribution density function 
!   (both minimally-dispersed and fully-dispersed as the 
!   sum of four log-normal distributions.
!
! d(0,imax):	output, particle size at 0, 1, 2, ..., imax points	[um]
! dm(imax):	output, particle size at middle of 0-1, 1-2, etc      [um]
! psd(imax):    output, particle size distribution density at dm	[um^-1]
! dpsd(imax):   output, Delta P for sections 0-1, 1-2, etc.  		[ ]
! ppsd(imax):   output, P for sections 0-1, 1-2, etc.  		[ ]
! imax:         input, length dm, psd, dpsdm, ppsd, etc.
! cmtrix: 	jmaxx coefficient matrix 
!               e.g.
!		w1   = cmtrix(1, 1): weight for first log-normal distribution
!               dln1 = cmtrix(2, 1): mean log-particle size of first log-normal distribution 	
!               sig1 = cmtrix(3, 1): sigma of log-particle size for first log-normal distribution
!		etc. 
!		careful with the dimension of dln and sig
!----------------------------------------------------------------------------
!
      integer :: i, j, imax, jmax	
      real(8), dimension(3, jmax) :: cmtrix
      real(8) :: d(0:imax), dm(imax)                 
      real(8) :: psd(imax), dpsd(imax), ppsd(imax)   ! for p(d), Delta P(d) and P(d)
      real(8) :: p, pp, w, dln, sig
      real(8) :: cn
      real(8), parameter :: eps=1.d-7
      real(8), parameter :: dref=1000.d0
      real(8) :: fu, fd, phi
!
      cn = 1.d0/dsqrt(2.d0*3.14159d0)
!
! initialise psd, dpsd, ppsd
!
      psd = 0.d0
      dpsd = 0.d0
      ppsd = 0.d0
!
! Estimate d using phi scale. phi varies between from 9 to -1
! with increment 0.1. Reference particle size d0 = 1000 um
!
      fu = 10.d0
      fd = -1.d0
      do i = 0, imax
	phi = fu - i*(fu-fd)/imax
	d(i) = dref/2.d0**phi
      enddo 	
!
      do i = 1, imax
        dm(i) = (d(i)-d(i-1))/dlog(d(i)/d(i-1))           ! d = delta(d) / delta(log d), because del (log d) = (del d) / d   ! mklose [17122014]

        pp = 0.d0
	do j = 1, jmax
	  w = cmtrix(1, j)
          dln = cmtrix(2, j)
          sig = cmtrix(3, j)
	  if ( (w.gt.eps) .and. (sig.ne.0.) ) then
	    p = w*cn/sig*dexp( -(dlog(dm(i))-dln)**2/(2*sig**2) )
          else
            p=0.d0	
          endif
	  pp = pp + p
	enddo
!
	dpsd(i) = pp*( dlog(d(i)) - dlog(d(i-1)) )        ! Delta P over i
	if (i.eq.1) then 
          ppsd(i) = 0.d0 + dpsd(i)                        ! P(d), with P(0) = 0
	else
          ppsd(i) = ppsd(i-1) + dpsd(i)      
	endif
        psd(i) = pp/dm(i)                                 ! p(d), particle size distribution density

      enddo
!
! Normalize in case ppsd(imax) is not 1
!
      dpsd = dpsd/ppsd(imax)
      psd  = psd/ppsd(imax)
      ppsd = ppsd/ppsd(imax)

!
     end subroutine
!*****************************************************************************

END MODULE uoc_dust
