dustflux_get_mb.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE dustflux_get_mb(     &
        pustar                  &  !I [m/s] Wind friction speed 
       ,prhoa                   &  !I [kg/m3] air density at 2m height 
       ,pwg                     &  !I [m3/m3] volumetric water content 
       ,pz0                     &  !I [m] roughness length of surface
       ,pwsat                   &  !I [m3 m-3] saturation liquid water content
       ,pclay                   &  !I [frc] mass fraction clay
       ,psand                   &  !I [frc] mass fraction sand
       ,pdst_erod               &  !I [frc] erodible surface
       ,pwind10m                &  !I [m/s] wind at 10m altitude
       ,psfdst                  &  !O [kg/m2/sec] Vertical dust flux
       ,ksize                   &  !I [nbr] number of points for calculation
       )
!
USE mode_dstmblutl                     !Dust mobilization subroutines
USE modd_dst_surf, ONLY :  cvermod
USE modd_dstmbl, ONLY : xflx_mss_fdg_fctm, ntex, nmode, ndp, nbin, xcst_slt, &
                        xdmt_slt_opt
USE modd_csts, only: xpi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!INPUT, set their dimensions to their passed lengths or to KSIZE ?
INTEGER, INTENT(IN)                  :: ksize    ![nbr] length of passed arrays
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pustar   ![m/s] wind friction speed
REAL, INTENT(IN), DIMENSION(KSIZE)   :: prhoa    ![kg/m3] air density at 2m height
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pclay    ![frc] mass fraction clay
REAL, INTENT(IN), DIMENSION(KSIZE)   :: psand    ![frc] mass fraction sand
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pdst_erod![frc]
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pwg      ![m3 m-3] volumetric water fraction
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pwsat    ![m3 m-3] saturation water content
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pz0      ![m] surface roughness length
REAL, INTENT(IN), DIMENSION(KSIZE)   :: pwind10m ![m/s] wind at 10m altitude
!OUTPUT the flux of dust
REAL, INTENT(OUT), DIMENSION(KSIZE)  :: psfdst   ! [kg m-2 s-1] Output flux of atmospheric dust
!
!!!!!!!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&!!!!!!

!#ifdef AlG01
!    real,parameter::XFLX_MSS_FDG_FCT=28. ! [frc] Global mass flux tuning factor (a posteriori)
!#else
!    real,parameter::XFLX_MSS_FDG_FCT=7.0e-4 ! [frc] Global mass flux tuning factor (a posteriori)
!    real,parameter::XFLX_MSS_FDG_FCT=21.0e-4 ! [frc] Global mass flux tuning factor (a posteriori)
!    real,parameter::XFLX_MSS_FDG_FCT=12.0e-4 ! [frc] values used in Masdev47
!
!Define local variables:
LOGICAL, DIMENSION(KSIZE) :: gflg_mbl                  ! [frc] Mobilization candidate flag
!REAL,    DIMENSION(KSIZE) :: ZMBL_BSN_FCT              ! [frc] enhancement factor for grid cells with higher erodibility
REAL,    DIMENSION(KSIZE) :: zwnd_rfr                  ! [m s-1] wind speed at reference level 
REAL,    DIMENSION(KSIZE) :: zwnd_frc_thr_slt          ! [m/s] Threshold wind friction speed when all effects taken into account
REAL,    DIMENSION(KSIZE) :: zwnd_rfr_thr_slt          ! [m s-1] Threshold wind speed at reference level
REAL,    DIMENSION(KSIZE) :: zgwc_sfc                  ! [kg/kg] Gravimetric water content
REAL,    DIMENSION(KSIZE) :: zfrc_thr_ncr_wtr          ! [frc] Fraction by which soil wetness increases threshold wind
REAL,    DIMENSION(KSIZE) :: zfrc_thr_ncr_drg          ! [frc] fraction by which drag partitioning increases threshold wind
REAL,    DIMENSION(KSIZE) :: zwnd_frc_slt              ! [m/s] wind friction speed after modified for saltation feedbacks
REAL,    DIMENSION(KSIZE,NBIN) :: zflx_mss_hrz_slt_ttl_wbn  ! [kg m-1 s-1] Vertically integrated horizontal saltation soil flux for a wind bin 
REAL,    DIMENSION(KSIZE,NBIN) :: zflx_mss_vrt_dst_ttl_wbn  ! [kg m-2 s-1]
REAL,    DIMENSION(KSIZE) :: zdst_slt_flx_rat_ttl      ! [m-1] ratio of vertical to horizontal flux (alpha in several papers)
REAL,    DIMENSION(KSIZE) :: zsilt                     ! [frc] dummy for fraction of silt 
REAL,    DIMENSION(KSIZE) :: zwprm                     ! threshold soil wetness
INTEGER, DIMENSION(KSIZE) :: itext                     ! soil texture
REAL, DIMENSION(NTEX,NDP) :: zdsrlv                    ! Surface relative des grains du sol
REAL, DIMENSION(KSIZE,NBIN) :: zdsbin                    !
REAL, DIMENSION(KSIZE,NBIN) :: zgamma                    !
INTEGER, DIMENSION(NBIN+1)   :: iseuil
REAL, DIMENSION(NBIN,2)   :: zpcen
REAL, DIMENSION(NTEX)     :: zzs0                      ! rugosité de la surface lisse
REAL, DIMENSION(NDP)      :: zdp                       ! [µm] diamètre du particule
REAL :: zdlndp, zrgh_z0
INTEGER :: i                  !Counter for number of points (used in loops)
INTEGER :: idp                !Counter for number of particle
INTEGER :: itex               !Counter for number of texture
INTEGER :: is
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('DUSTFLUX_GET_MB',0,zhook_handle)
!
!Initialize mobilization candidate flag
gflg_mbl(:) = .true.
!fxm: Get erodibility limitation factor, use something connected to amount of sand
!Discuss with Valery Masson
!ZMBL_BSN_FCT(:) = PSAND(:)
! utilisé dans le calcul de l'effet Owen 
zwnd_rfr(:) = pwind10m(:)
!
! Initialize vertical dust flux
zflx_mss_vrt_dst_ttl_wbn(:,:) = 0.d0
zflx_mss_hrz_slt_ttl_wbn(:,:) = 0.d0
!
psfdst(:) = 0.0d0
!
! CLAY : CLAY >= 0.40 SILT < 0.40 SAND < 0.45
! SANDY CLAY : CLAY >= 0.36 SAND >= 0.45 
! SILTY CLAY : CLAY >= 0.40 SILT >= 0.40 
! SILT : SILT >= 0.8 CLAY < 0.12
! SAND : SAND >= 0.3*CLAY + 0.87
! SANDY CLAY LOAM : CLAY >= 0.28  CLAY < 0.36 SAND >= 0.45 | CLAY >= 0.20 CLAY < 0.28 SILT < 0.28
! SILTY CLAY LOAM : CLAY >= 0.28 CLAY < 0.40 SAND < 0.20
! CLAY LOAM : CLAY >= 0.28 CLAY < 0.40 SAND >= 0.20 SAND < 0.45
! SILT LOAM : SILT >= 0.8 CLAY >= 0.12 | SILT >= 0.5 SILT < 0.8 CLAY < 0.28
! LOAMY SAND : SAND >= CLAY + 0.7 SAND < 0.3*CLAY + 0.87
! SANDY LOAM : SAND >= 0.52  CLAY < 0.20 | SAND >= (0.5 - CLAY)  CLAY < 0.07
! LOAM : CLAY >= 0.20 CLAY < 0.28 SILT >= 0.28 SILT < 0.5 | SAND >= (0.5 - CLAY)  CLAY < 0.20
DO i = 1, ksize
  !
  zsilt(i) = 1.0 - pclay(i) - psand(i)
  IF (zsilt(i) <= 0.) zsilt(i) = 0.0
  zwprm(i) = 3.0*pclay(i) * (0.17 + 0.0014 * pclay(i))
  zwprm(i) = min( 0.15, max(0.053, zwprm(i)) )
  ! 
  !tous les cas CLAY >= 0.28 sont traités  IF ( PCLAY(I) >= 0.28 ) THEN
    IF ( PSAND(I) >= 0.45 ) THEN      
      IF (PCLAY(I) >= 0.36 ) THEN     ! Sandy Clay 
        ITEXT(I) = 9
        !ZWPRM(I) = 10.0E-2
      ELSE                            ! Sandy Clay Loam 
        ITEXT(I) = 6 
        !ZWPRM(I) = 6.0E-2
      ENDIF
    ELSEIF ( PCLAY(I) >= 0.40 ) THEN
      IF ( ZSILT(I) >= 0.40 ) THEN    ! Silty Clay 
        ITEXT(I) = 10
        !ZWPRM(I) = 10.5E-2
      ELSE                            ! Clay 
        ITEXT(I) = 11
        !ZWPRM(I) = 11.5E-2
      ENDIF
    ELSEIF (PSAND(I) >= 0.20 ) THEN  ! Clay Loam 
      ITEXT(I) = 8
      !ZWPRM(I) = 6.8E-2
    ELSE                             ! Silty Clay Loam 
      ITEXT(I) = 7
      !ZWPRM(I) = 6.8E-2
    ENDIF
  ENDIF
  ! les cas ZSILT >= 0.5 .AND. PCLAY < 0.28 sont traités
  ! les cas ZSILT < 0.5 .AND. PCLAY >= 0.20 sont traités 
  ! ne sont pas traités les cas PCLAY < 0.20 .AND. ZSILT < 0.5 => SAND >= 0.3
  IF ( ZSILT(I) >= 0.8 .AND. PCLAY(I) < 0.12 ) THEN ! Silt 
    ITEXT(I) = 12
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( PCLAY(I) < 0.28 ) THEN    ! ( clay est forcément < 0.28 )
    IF ( ZSILT(I) >= 0.5 ) THEN      ! Silt Loam 
      ITEXT(I) = 4
      !ZWPRM(I) = 5.0E-2
    ELSEIF ( PCLAY(I) >= 0.20 ) THEN
      IF ( ZSILT(I) >= 0.28 ) THEN   ! Loam  
        ITEXT(I) = 5
        !ZWPRM(I) = 4.0E-2
      ELSE                           ! Sandy Clay Loam 
        ITEXT(I) = 6
        !ZWPRM(I) = 4.0E-2
      ENDIF
    ENDIF
  ENDIF
  ! les cas SAND >= 0.87 sont traités: silt < 0.13, clay < 0.1 => entrent dans les cas non encore traités
  ! les cas SAND >= 0.7 => CLAY < 0.15 , SILT < 0.3 sont traités ) => ""
  ! les cas SAND >=0.52 .AND. CLAY < 0.20 sont traités SILT < 0.48 => ""
  ! les cas restants considèrement pclay < 0.20, sand >= 0.5 - pclay => sand >=  0.3
  ! => clay + sand >= 0.5 => silt < 0.5
  IF ( PSAND(I) >= (0.3*PCLAY(I) + 0.87) ) THEN   ! Sand 
    ITEXT(I) = 1
    !ZWPRM(I) = 1.5E-2
  ELSEIF ( PSAND(I) >= (PCLAY(I) + 0.7) ) THEN    ! Loamy Sand
    ITEXT(I) = 2
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( PSAND(I) >= 0.52 .AND. PCLAY(I) < 0.20 ) THEN ! Sandy Loam            
    ITEXT(I) = 3
    !ZWPRM(I) = 3.0E-2   
  ELSEIF ( PSAND(I) >= (0.5 - PCLAY(I)) ) THEN
    IF ( PCLAY(I) < 0.07 ) THEN                   ! Sandy Loam
      ITEXT(I) = 3
      !ZWPRM(I) = 3.0E-2
    ELSEIF ( PCLAY(I) < 0.20 ) THEN               ! Loam
      ITEXT(I) = 5
      !ZWPRM(I) = 4.0E-2
    ENDIF
  ENDIF
  !
ENDDO
!
ZDLNDP = 0.1d0        ! [µm]  Dln(DP)
 CALL DISTRIBUTION (NMODE, ZDLNDP, ITEXT, ZDP, ZDSRLV, ZZS0)
!
ISEUIL = (/0, 30, 40, 65, NDP/)
ZPCEN(:,1) = (/0.005, 0.006, 0.1, 0.10/)
ZPCEN(:,2) = (/0.005, 0.006, 1.0, 0.12/)
!
ZDSBIN(:,:) = 0.0
!
DO IS = 1, NBIN
  DO IDP = ISEUIL(IS)+1, ISEUIL(IS+1)
    DO I = 1, KSIZE
      ITEX = ITEXT(I)
      ZDSBIN(I,IS) = ZDSBIN(I,IS) + ZDSRLV(ITEX,IDP) / FLOAT(ISEUIL(IS+1)-ISEUIL(IS))
      IF (ITEX==4) THEN
        ZGAMMA(I,IS) = ZPCEN(IS,1)
      ELSE
        ZGAMMA(I,IS) = ZPCEN(IS,2)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
! Adjust threshold velocity for inhibition by roughness elements
DO I = 1, SIZE(ZFRC_THR_NCR_DRG)
  ITEX    = ITEXT(I)
  ZRGH_Z0 = PZ0(I)
  IF (ZRGH_Z0 <= ZZS0(ITEX)) ZRGH_Z0 = ZZS0(ITEX)
  ! Factor by which surface roughness increases threshold friction velocity 
  !++grini: fxm: USE WHOLE ARRAY OF Z0 INSTEAD OF ONLY RGH_MMN_MBL AS IN OLD CODE 
  CALL FRC_THR_NCR_DRG_GET(ZRGH_Z0, ZZS0(ITEX), ZFRC_THR_NCR_DRG(I))
ENDDO
!
! Convert volumetric water content to gravimetric water content
 CALL VWC2GWC(GFLG_MBL, PWSAT, PWG, ZGWC_SFC)
! Factor by which soil moisture increases threshold friction velocity 
 CALL FRC_THR_NCR_WTR_GET (GFLG_MBL, ZWPRM, ZGWC_SFC, ZFRC_THR_NCR_WTR)
ZFRC_THR_NCR_WTR(:) = MAX(ZFRC_THR_NCR_WTR(:), 1.0)
!
!
 CALL WND_FRC_THR_SLT_GET(PRHOA, XDMT_SLT_OPT, ZWND_FRC_THR_SLT)
!
DO I = 1,KSIZE
  IF (GFLG_MBL(I)) THEN
    ZWND_FRC_THR_SLT(I) = ZWND_FRC_THR_SLT(I) * ZFRC_THR_NCR_WTR(I) &  ! [frc] Adjustment for moisture
                                              * ZFRC_THR_NCR_DRG(I)    ! I [frc] Adjustment for roughness
    ZWND_RFR_THR_SLT(I) = ZWND_RFR(I) * ZWND_FRC_THR_SLT(I) / PUSTAR(I)
  ENDIF
ENDDO
!
! CHECK IF THIS CAN BE USED EASILY
! NEEDS 10M WIND SPEED WHICH IS MAYBE KNOWN MAYBE NOT !
! Saltation increases friction speed by roughening surface
 CALL WND_FRC_SLT_GET(GFLG_MBL, PUSTAR, ZWND_RFR, ZWND_RFR_THR_SLT, ZWND_FRC_SLT)
!
DO IS = 1,NBIN
  CALL FLX_MSS_HRZ_SLT_TTL_WHI79_GET(ZGAMMA(:,IS)*ZDSBIN(:,IS), GFLG_MBL, PRHOA, &
         ZWND_FRC_SLT, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS))
ENDDO
!
! Vertical dust mass flux
DO IS = 1,NBIN
  CALL FLX_MSS_VRT_DST_AUST_GET(GFLG_MBL, PRHOA, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS), & 
                      ZDST_SLT_FLX_RAT_TTL, ZFLX_MSS_VRT_DST_TTL_WBN(:,IS))
  !
  DO I = 1, KSIZE
    ! Vertically integrated streamwise mass flux in wind bin
    PSFDST(I) =  PSFDST(I) + PDST_EROD(I) * XFLX_MSS_FDG_FCTM * ZFLX_MSS_VRT_DST_TTL_WBN(I,IS)
  ENDDO
ENDDO
!
END SUBROUTINE DUSTFLUX_GET_MB 

  IF ( pclay(i) >= 0.28 ) THEN
    IF ( psand(i) >= 0.45 ) THEN      
      IF (pclay(i) >= 0.36 ) THEN     ! Sandy Clay 
        itext(i) = 9
        !ZWPRM(I) = 10.0E-2
      ELSE                            ! Sandy Clay Loam 
        itext(i) = 6 
        !ZWPRM(I) = 6.0E-2
      ENDIF
    ELSEIF ( pclay(i) >= 0.40 ) THEN
      IF ( zsilt(i) >= 0.40 ) THEN    ! Silty Clay 
        itext(i) = 10
        !ZWPRM(I) = 10.5E-2
      ELSE                            ! Clay 
        itext(i) = 11
        !ZWPRM(I) = 11.5E-2
      ENDIF
    ELSEIF (psand(i) >= 0.20 ) THEN  ! Clay Loam 
      itext(i) = 8
      !ZWPRM(I) = 6.8E-2
    ELSE                             ! Silty Clay Loam 
      itext(i) = 7
      !ZWPRM(I) = 6.8E-2
    ENDIF
  ENDIF
  ! les cas ZSILT >= 0.5 .AND. PCLAY < 0.28 sont traités  ! les cas ZSILT < 0.5 .AND. PCLAY >= 0.20 sont traités 
  ! ne sont pas traités les cas PCLAY < 0.20 .AND. ZSILT < 0.5 => SAND >= 0.3
  IF ( ZSILT(I) >= 0.8 .AND. PCLAY(I) < 0.12 ) THEN ! Silt 
    ITEXT(I) = 12
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( PCLAY(I) < 0.28 ) THEN    ! ( clay est forcément < 0.28 )
    IF ( ZSILT(I) >= 0.5 ) THEN      ! Silt Loam 
      ITEXT(I) = 4
      !ZWPRM(I) = 5.0E-2
    ELSEIF ( PCLAY(I) >= 0.20 ) THEN
      IF ( ZSILT(I) >= 0.28 ) THEN   ! Loam  
        ITEXT(I) = 5
        !ZWPRM(I) = 4.0E-2
      ELSE                           ! Sandy Clay Loam 
        ITEXT(I) = 6
        !ZWPRM(I) = 4.0E-2
      ENDIF
    ENDIF
  ENDIF
  ! les cas SAND >= 0.87 sont traités: silt < 0.13, clay < 0.1 => entrent dans les cas non encore traités
  ! les cas SAND >= 0.7 => CLAY < 0.15 , SILT < 0.3 sont traités ) => ""
  ! les cas SAND >=0.52 .AND. CLAY < 0.20 sont traités SILT < 0.48 => ""
  ! les cas restants considèrement pclay < 0.20, sand >= 0.5 - pclay => sand >=  0.3
  ! => clay + sand >= 0.5 => silt < 0.5
  IF ( PSAND(I) >= (0.3*PCLAY(I) + 0.87) ) THEN   ! Sand 
    ITEXT(I) = 1
    !ZWPRM(I) = 1.5E-2
  ELSEIF ( PSAND(I) >= (PCLAY(I) + 0.7) ) THEN    ! Loamy Sand
    ITEXT(I) = 2
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( PSAND(I) >= 0.52 .AND. PCLAY(I) < 0.20 ) THEN ! Sandy Loam            
    ITEXT(I) = 3
    !ZWPRM(I) = 3.0E-2   
  ELSEIF ( PSAND(I) >= (0.5 - PCLAY(I)) ) THEN
    IF ( PCLAY(I) < 0.07 ) THEN                   ! Sandy Loam
      ITEXT(I) = 3
      !ZWPRM(I) = 3.0E-2
    ELSEIF ( PCLAY(I) < 0.20 ) THEN               ! Loam
      ITEXT(I) = 5
      !ZWPRM(I) = 4.0E-2
    ENDIF
  ENDIF
  !
ENDDO
!
ZDLNDP = 0.1d0        ! [µm]  Dln(DP)
 CALL DISTRIBUTION (NMODE, ZDLNDP, ITEXT, ZDP, ZDSRLV, ZZS0)
!
ISEUIL = (/0, 30, 40, 65, NDP/)
ZPCEN(:,1) = (/0.005, 0.006, 0.1, 0.10/)
ZPCEN(:,2) = (/0.005, 0.006, 1.0, 0.12/)
!
ZDSBIN(:,:) = 0.0
!
DO IS = 1, NBIN
  DO IDP = ISEUIL(IS)+1, ISEUIL(IS+1)
    DO I = 1, KSIZE
      ITEX = ITEXT(I)
      ZDSBIN(I,IS) = ZDSBIN(I,IS) + ZDSRLV(ITEX,IDP) / FLOAT(ISEUIL(IS+1)-ISEUIL(IS))
      IF (ITEX==4) THEN
        ZGAMMA(I,IS) = ZPCEN(IS,1)
      ELSE
        ZGAMMA(I,IS) = ZPCEN(IS,2)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
! Adjust threshold velocity for inhibition by roughness elements
DO I = 1, SIZE(ZFRC_THR_NCR_DRG)
  ITEX    = ITEXT(I)
  ZRGH_Z0 = PZ0(I)
  IF (ZRGH_Z0 <= ZZS0(ITEX)) ZRGH_Z0 = ZZS0(ITEX)
  ! Factor by which surface roughness increases threshold friction velocity 
  !++grini: fxm: USE WHOLE ARRAY OF Z0 INSTEAD OF ONLY RGH_MMN_MBL AS IN OLD CODE 
  CALL FRC_THR_NCR_DRG_GET(ZRGH_Z0, ZZS0(ITEX), ZFRC_THR_NCR_DRG(I))
ENDDO
!
! Convert volumetric water content to gravimetric water content
 CALL VWC2GWC(GFLG_MBL, PWSAT, PWG, ZGWC_SFC)
! Factor by which soil moisture increases threshold friction velocity 
 CALL FRC_THR_NCR_WTR_GET (GFLG_MBL, ZWPRM, ZGWC_SFC, ZFRC_THR_NCR_WTR)
ZFRC_THR_NCR_WTR(:) = MAX(ZFRC_THR_NCR_WTR(:), 1.0)
!
!
 CALL WND_FRC_THR_SLT_GET(PRHOA, XDMT_SLT_OPT, ZWND_FRC_THR_SLT)
!
DO I = 1,KSIZE
  IF (GFLG_MBL(I)) THEN
    ZWND_FRC_THR_SLT(I) = ZWND_FRC_THR_SLT(I) * ZFRC_THR_NCR_WTR(I) &  ! [frc] Adjustment for moisture
                                              * ZFRC_THR_NCR_DRG(I)    ! I [frc] Adjustment for roughness
    ZWND_RFR_THR_SLT(I) = ZWND_RFR(I) * ZWND_FRC_THR_SLT(I) / PUSTAR(I)
  ENDIF
ENDDO
!
! CHECK IF THIS CAN BE USED EASILY
! NEEDS 10M WIND SPEED WHICH IS MAYBE KNOWN MAYBE NOT !
! Saltation increases friction speed by roughening surface
 CALL WND_FRC_SLT_GET(GFLG_MBL, PUSTAR, ZWND_RFR, ZWND_RFR_THR_SLT, ZWND_FRC_SLT)
!
DO IS = 1,NBIN
  CALL FLX_MSS_HRZ_SLT_TTL_WHI79_GET(ZGAMMA(:,IS)*ZDSBIN(:,IS), GFLG_MBL, PRHOA, &
         ZWND_FRC_SLT, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS))
ENDDO
!
! Vertical dust mass flux
DO IS = 1,NBIN
  CALL FLX_MSS_VRT_DST_AUST_GET(GFLG_MBL, PRHOA, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS), & 
                      ZDST_SLT_FLX_RAT_TTL, ZFLX_MSS_VRT_DST_TTL_WBN(:,IS))
  !
  DO I = 1, KSIZE
    ! Vertically integrated streamwise mass flux in wind bin
    PSFDST(I) =  PSFDST(I) + PDST_EROD(I) * XFLX_MSS_FDG_FCTM * ZFLX_MSS_VRT_DST_TTL_WBN(I,IS)
  ENDDO
ENDDO
!
END SUBROUTINE DUSTFLUX_GET_MB 

  ! les cas ZSILT < 0.5 .AND. PCLAY >= 0.20 sont traités 
  ! ne sont pas traités les cas PCLAY < 0.20 .AND. ZSILT < 0.5 => SAND >= 0.3
  IF ( zsilt(i) >= 0.8 .AND. pclay(i) < 0.12 ) THEN ! Silt 
    itext(i) = 12
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( pclay(i) < 0.28 ) THEN    ! ( clay est forcément < 0.28 )
    IF ( zsilt(i) >= 0.5 ) THEN      ! Silt Loam 
      itext(i) = 4
      !ZWPRM(I) = 5.0E-2
    ELSEIF ( pclay(i) >= 0.20 ) THEN
      IF ( zsilt(i) >= 0.28 ) THEN   ! Loam  
        itext(i) = 5
        !ZWPRM(I) = 4.0E-2
      ELSE                           ! Sandy Clay Loam 
        itext(i) = 6
        !ZWPRM(I) = 4.0E-2
      ENDIF
    ENDIF
  ENDIF
  ! les cas SAND >= 0.87 sont traités: silt < 0.13, clay < 0.1 => entrent dans les cas non encore traités  ! les cas SAND >= 0.7 => CLAY < 0.15 , SILT < 0.3 sont traités ) => ""
  ! les cas SAND >=0.52 .AND. CLAY < 0.20 sont traités SILT < 0.48 => ""
  ! les cas restants considèrement pclay < 0.20, sand >= 0.5 - pclay => sand >=  0.3
  ! => clay + sand >= 0.5 => silt < 0.5
  IF ( PSAND(I) >= (0.3*PCLAY(I) + 0.87) ) THEN   ! Sand 
    ITEXT(I) = 1
    !ZWPRM(I) = 1.5E-2
  ELSEIF ( PSAND(I) >= (PCLAY(I) + 0.7) ) THEN    ! Loamy Sand
    ITEXT(I) = 2
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( PSAND(I) >= 0.52 .AND. PCLAY(I) < 0.20 ) THEN ! Sandy Loam            
    ITEXT(I) = 3
    !ZWPRM(I) = 3.0E-2   
  ELSEIF ( PSAND(I) >= (0.5 - PCLAY(I)) ) THEN
    IF ( PCLAY(I) < 0.07 ) THEN                   ! Sandy Loam
      ITEXT(I) = 3
      !ZWPRM(I) = 3.0E-2
    ELSEIF ( PCLAY(I) < 0.20 ) THEN               ! Loam
      ITEXT(I) = 5
      !ZWPRM(I) = 4.0E-2
    ENDIF
  ENDIF
  !
ENDDO
!
ZDLNDP = 0.1d0        ! [µm]  Dln(DP)
 CALL DISTRIBUTION (NMODE, ZDLNDP, ITEXT, ZDP, ZDSRLV, ZZS0)
!
ISEUIL = (/0, 30, 40, 65, NDP/)
ZPCEN(:,1) = (/0.005, 0.006, 0.1, 0.10/)
ZPCEN(:,2) = (/0.005, 0.006, 1.0, 0.12/)
!
ZDSBIN(:,:) = 0.0
!
DO IS = 1, NBIN
  DO IDP = ISEUIL(IS)+1, ISEUIL(IS+1)
    DO I = 1, KSIZE
      ITEX = ITEXT(I)
      ZDSBIN(I,IS) = ZDSBIN(I,IS) + ZDSRLV(ITEX,IDP) / FLOAT(ISEUIL(IS+1)-ISEUIL(IS))
      IF (ITEX==4) THEN
        ZGAMMA(I,IS) = ZPCEN(IS,1)
      ELSE
        ZGAMMA(I,IS) = ZPCEN(IS,2)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
! Adjust threshold velocity for inhibition by roughness elements
DO I = 1, SIZE(ZFRC_THR_NCR_DRG)
  ITEX    = ITEXT(I)
  ZRGH_Z0 = PZ0(I)
  IF (ZRGH_Z0 <= ZZS0(ITEX)) ZRGH_Z0 = ZZS0(ITEX)
  ! Factor by which surface roughness increases threshold friction velocity 
  !++grini: fxm: USE WHOLE ARRAY OF Z0 INSTEAD OF ONLY RGH_MMN_MBL AS IN OLD CODE 
  CALL FRC_THR_NCR_DRG_GET(ZRGH_Z0, ZZS0(ITEX), ZFRC_THR_NCR_DRG(I))
ENDDO
!
! Convert volumetric water content to gravimetric water content
 CALL VWC2GWC(GFLG_MBL, PWSAT, PWG, ZGWC_SFC)
! Factor by which soil moisture increases threshold friction velocity 
 CALL FRC_THR_NCR_WTR_GET (GFLG_MBL, ZWPRM, ZGWC_SFC, ZFRC_THR_NCR_WTR)
ZFRC_THR_NCR_WTR(:) = MAX(ZFRC_THR_NCR_WTR(:), 1.0)
!
!
 CALL WND_FRC_THR_SLT_GET(PRHOA, XDMT_SLT_OPT, ZWND_FRC_THR_SLT)
!
DO I = 1,KSIZE
  IF (GFLG_MBL(I)) THEN
    ZWND_FRC_THR_SLT(I) = ZWND_FRC_THR_SLT(I) * ZFRC_THR_NCR_WTR(I) &  ! [frc] Adjustment for moisture
                                              * ZFRC_THR_NCR_DRG(I)    ! I [frc] Adjustment for roughness
    ZWND_RFR_THR_SLT(I) = ZWND_RFR(I) * ZWND_FRC_THR_SLT(I) / PUSTAR(I)
  ENDIF
ENDDO
!
! CHECK IF THIS CAN BE USED EASILY
! NEEDS 10M WIND SPEED WHICH IS MAYBE KNOWN MAYBE NOT !
! Saltation increases friction speed by roughening surface
 CALL WND_FRC_SLT_GET(GFLG_MBL, PUSTAR, ZWND_RFR, ZWND_RFR_THR_SLT, ZWND_FRC_SLT)
!
DO IS = 1,NBIN
  CALL FLX_MSS_HRZ_SLT_TTL_WHI79_GET(ZGAMMA(:,IS)*ZDSBIN(:,IS), GFLG_MBL, PRHOA, &
         ZWND_FRC_SLT, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS))
ENDDO
!
! Vertical dust mass flux
DO IS = 1,NBIN
  CALL FLX_MSS_VRT_DST_AUST_GET(GFLG_MBL, PRHOA, ZWND_FRC_THR_SLT, ZFLX_MSS_HRZ_SLT_TTL_WBN(:,IS), & 
                      ZDST_SLT_FLX_RAT_TTL, ZFLX_MSS_VRT_DST_TTL_WBN(:,IS))
  !
  DO I = 1, KSIZE
    ! Vertically integrated streamwise mass flux in wind bin
    PSFDST(I) =  PSFDST(I) + PDST_EROD(I) * XFLX_MSS_FDG_FCTM * ZFLX_MSS_VRT_DST_TTL_WBN(I,IS)
  ENDDO
ENDDO
!
END SUBROUTINE DUSTFLUX_GET_MB 

  ! les cas SAND >= 0.7 => CLAY < 0.15 , SILT < 0.3 sont traités ) => ""
  ! les cas SAND >=0.52 .AND. CLAY < 0.20 sont traités SILT < 0.48 => ""
  ! les cas restants considèrement pclay < 0.20, sand >= 0.5 - pclay => sand >=  0.3
  ! => clay + sand >= 0.5 => silt < 0.5
  IF ( psand(i) >= (0.3*pclay(i) + 0.87) ) THEN   ! Sand 
    itext(i) = 1
    !ZWPRM(I) = 1.5E-2
  ELSEIF ( psand(i) >= (pclay(i) + 0.7) ) THEN    ! Loamy Sand
    itext(i) = 2
    !ZWPRM(I) = 2.5E-2
  ELSEIF ( psand(i) >= 0.52 .AND. pclay(i) < 0.20 ) THEN ! Sandy Loam            
    itext(i) = 3
    !ZWPRM(I) = 3.0E-2   
  ELSEIF ( psand(i) >= (0.5 - pclay(i)) ) THEN
    IF ( pclay(i) < 0.07 ) THEN                   ! Sandy Loam
      itext(i) = 3
      !ZWPRM(I) = 3.0E-2
    ELSEIF ( pclay(i) < 0.20 ) THEN               ! Loam
      itext(i) = 5
      !ZWPRM(I) = 4.0E-2
    ENDIF
  ENDIF
  !
ENDDO
!
zdlndp = 0.1d0        ! [µm]  Dln(DP)
 CALL distribution(nmode, zdlndp, itext, zdp, zdsrlv, zzs0)
!
iseuil = (/0, 30, 40, 65, ndp/)
zpcen(:,1) = (/0.005, 0.006, 0.1, 0.10/)
zpcen(:,2) = (/0.005, 0.006, 1.0, 0.12/)
!
zdsbin(:,:) = 0.0
!
DO is = 1, nbin
  DO idp = iseuil(is)+1, iseuil(is+1)
    DO i = 1, ksize
      itex = itext(i)
      zdsbin(i,is) = zdsbin(i,is) + zdsrlv(itex,idp) / float(iseuil(is+1)-iseuil(is))
      IF (itex==4) THEN
        zgamma(i,is) = zpcen(is,1)
      ELSE
        zgamma(i,is) = zpcen(is,2)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
! Adjust threshold velocity for inhibition by roughness elements
DO i = 1, SIZE(zfrc_thr_ncr_drg)
  itex    = itext(i)
  zrgh_z0 = pz0(i)
  IF (zrgh_z0 <= zzs0(itex)) zrgh_z0 = zzs0(itex)
  ! Factor by which surface roughness increases threshold friction velocity 
  !++grini: fxm: USE WHOLE ARRAY OF Z0 INSTEAD OF ONLY RGH_MMN_MBL AS IN OLD CODE 
  CALL frc_thr_ncr_drg_get(zrgh_z0, zzs0(itex), zfrc_thr_ncr_drg(i))
ENDDO
!
! Convert volumetric water content to gravimetric water content
 CALL vwc2gwc(gflg_mbl, pwsat, pwg, zgwc_sfc)
! Factor by which soil moisture increases threshold friction velocity 
 CALL frc_thr_ncr_wtr_get(gflg_mbl, zwprm, zgwc_sfc, zfrc_thr_ncr_wtr)
zfrc_thr_ncr_wtr(:) = max(zfrc_thr_ncr_wtr(:), 1.0)
!
!
 CALL wnd_frc_thr_slt_get(prhoa, xdmt_slt_opt, zwnd_frc_thr_slt)
!
DO i = 1,ksize
  IF (gflg_mbl(i)) THEN
    zwnd_frc_thr_slt(i) = zwnd_frc_thr_slt(i) * zfrc_thr_ncr_wtr(i) &  ! [frc] Adjustment for moisture
                                              * zfrc_thr_ncr_drg(i)    ! I [frc] Adjustment for roughness
    zwnd_rfr_thr_slt(i) = zwnd_rfr(i) * zwnd_frc_thr_slt(i) / pustar(i)
  ENDIF
ENDDO
!
! CHECK IF THIS CAN BE USED EASILY
! NEEDS 10M WIND SPEED WHICH IS MAYBE KNOWN MAYBE NOT !
! Saltation increases friction speed by roughening surface
 CALL wnd_frc_slt_get(gflg_mbl, pustar, zwnd_rfr, zwnd_rfr_thr_slt, zwnd_frc_slt)
!
DO is = 1,nbin
  CALL flx_mss_hrz_slt_ttl_whi79_get(zgamma(:,is)*zdsbin(:,is), gflg_mbl, prhoa, &
         zwnd_frc_slt, zwnd_frc_thr_slt, zflx_mss_hrz_slt_ttl_wbn(:,is))
ENDDO
!
! Vertical dust mass flux
DO is = 1,nbin
  CALL flx_mss_vrt_dst_aust_get(gflg_mbl, prhoa, zwnd_frc_thr_slt, zflx_mss_hrz_slt_ttl_wbn(:,is), & 
                      zdst_slt_flx_rat_ttl, zflx_mss_vrt_dst_ttl_wbn(:,is))
  !
  DO i = 1, ksize
    ! Vertically integrated streamwise mass flux in wind bin
    psfdst(i) =  psfdst(i) + pdst_erod(i) * xflx_mss_fdg_fctm * zflx_mss_vrt_dst_ttl_wbn(i,is)
  ENDDO
ENDDO
!
END SUBROUTINE dustflux_get_mb