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