dustflux_get.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(            &
         pustar,                    &
         prhoa,                     & 
         pwg,                       & 
         pz0,                       &
         pwsat,                     &
         pclay,                     &
         psand,                     &
         pwind10m,                  &
         psfdst,                    &
         ksize                      &
         )  
!
!PURPOSE: 
!Take input from ISBA land surface model and
!calculate a dust flux which is consistent with the input.

!THEORY:
!Based on Marticorena/Bergametti, 1995 and Zender et al 2003 (JGR)

!CODE HISTORY
!Code is a modified version of dstmbl.F90 in the DEAD model
!Original version was downloaded from the DEAD homepage
!http://dust.ess.uci.edu/dead/ on January 10th 2005

!AUTHOR (or rather "code modifyer")
!Alf Grini <alf.grini@cnrm.meteo.fr>
!
USE modd_dst_surf, ONLY :  xflx_mss_fdg_fct
USE modd_dstmbl,   ONLY : xrgh_mmn_smt, xcst_slt, xdmt_slt_opt
USE mode_dstmblutl               !Dust mobilization subroutines
!
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)   :: 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
!REAL,PARAMETER::flx_mss_fdg_fctm=4.0e-4 ! [frc] Global mass flux tuning factor (a posteriori) (proposez by Pierre)
!
!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) :: zgwc_sfc                  ! [kg/kg] Gravimetric water content
REAL,    DIMENSION(KSIZE) :: zgwc_thr                  ! [kg kg-1] Threshold 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) :: ZWND_RFR_THR_SLT          ! [m s-1] Threshold wind speed at reference level
REAL,    DIMENSION(KSIZE) :: zcoef
REAL,    DIMENSION(KSIZE) :: zflx_mss_hrz_slt_ttl_wbn  ! [kg m-1 s-1] Vertically integrated horizontal saltation soil flux for a wind bin 
REAL,    DIMENSION(KSIZE) :: 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  :: zclay(ksize)
INTEGER         :: i                    !Counter for number of points (used in loops)
REAL(KIND=JPRB) :: zhook_handle
    
!Allocate the local variables
IF (lhook) CALL dr_hook('DUSTFLUX_GET',0,zhook_handle)
!
zclay(:)=0.2
!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
!
! Old Alf Grini code: bug on DEAD detected ? 
! Modification proposed by M. Mokhtari.. Accepted for all cases
!  if (CVERMOD=='CMDVER') then
zgwc_thr(:) = min(0.14,max(0.02,3. * pclay(:) * (0.17 + 0.14 * pclay(:))))
!ZGWC_THR(:) = PCLAY(:) * (0.17d0 + 0.14d0 * PCLAY(:))
!  else 
! gwc_thr=mss_frc_cly*(0.17d0+0.14d0*mss_frc_cly) ! [m3 m-3] FMB99 p. 155 (14)
! fxm: 19991105 remove factor of mss_frc_cly from gwc_thr to improve large scale behavior
! Begin Old Alf code 
!   gwc_thr(lon_idx)=0.17d0+0.14d0*mss_frc_cly(lon_idx) ! [m3 m-3] 
!  endif
!
! 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
DO i = 1,SIZE(pz0)
  CALL frc_thr_ncr_drg_get(pz0(i), xrgh_mmn_smt, 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, zgwc_thr, zgwc_sfc, zfrc_thr_ncr_wtr)
!
! fxm: Use surface density not midlayer density
 CALL wnd_frc_thr_slt_get(prhoa, xdmt_slt_opt, zwnd_frc_thr_slt)
!
DO i = 1, ksize
  zwnd_frc_thr_slt(i) =   & ! [m s-1] Threshold friction velocity for saltation
    zwnd_frc_thr_slt(i) * & ! [m s-1] Threshold for dry, flat ground
    zfrc_thr_ncr_wtr(i)     ! [frc] Adjustment for moisture
ENDDO
!
! Threshold saltation wind speed
!Needed for the saltation feedback roughening effect
!do i=1,KSIZE
!   if (flg_mbl(i)) then
!      wnd_rfr_thr_slt(i)= & ! [m s-1] Threshold 10 m wind speed for saltation
!           wnd_rfr(i)*wnd_frc_thr_slt(i)/PUSTAR(i) !++alfgr 
!   endif                  ! endif flg_mbl
!end do                    ! end loop over lon

!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( &
!    flg_mbl, & ! I [flg] Mobilization candidate flag
!     PUSTAR, & ! I [m s-1] Surface friction velocity
!     wnd_frc_slt, & ! O [m s-1] Saltating friction velocity
!     wnd_rfr, & ! I [m s-1] Wind speed at reference height
!     wnd_rfr_thr_slt) ! I [m s-1] Threshold 10 m wind speed for saltation
!
!
!Skip the roughening of surface effect for now, and 
!just use the wind friction speed as it is modified
!by drag partitioning
zwnd_frc_slt(:) = pustar(:) / zfrc_thr_ncr_drg(:) 
! 
! Horizontal streamwise mass flux for old "bulk" formulation
zcoef(:) = xcst_slt
 CALL flx_mss_hrz_slt_ttl_whi79_get(zcoef, gflg_mbl, prhoa, zwnd_frc_slt, &
        zwnd_frc_thr_slt, zflx_mss_hrz_slt_ttl_wbn)
!
! Apply land surface and vegetation limitations and global tuning factor
DO i = 1, ksize
  zflx_mss_hrz_slt_ttl_wbn(i) = zflx_mss_hrz_slt_ttl_wbn(i) & ! [kg m-2 s-1]
            !*lnd_frc_mbl(i) & ! [frc] Bare ground fraction
                                * zmbl_bsn_fct(i)           & ! [frc] Erodibility factor
                                * xflx_mss_fdg_fct            ! [frc] Global mass flux tuning factor (empirical) 
ENDDO
!    
! Vertical dust mass flux
 CALL flx_mss_vrt_dst_ttl_mab95_get(gflg_mbl, zclay, zflx_mss_hrz_slt_ttl_wbn, & 
        zdst_slt_flx_rat_ttl, zflx_mss_vrt_dst_ttl_wbn)
!
!Assign the output vertical dust flux to the value calculated
!PSFDST(:) = flx_mss_vrt_dst_ttl_wbn(:)
psfdst(:) = zdst_slt_flx_rat_ttl(:) * zflx_mss_hrz_slt_ttl_wbn(:)
!
IF (lhook) CALL dr_hook('DUSTFLUX_GET',1,zhook_handle)
!
END SUBROUTINE dustflux_get