mode_dstmblutl.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.
! ######spl
MODULE mode_dstmblutl ! [mdl] Mobilization utilities
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
 CONTAINS
!
!----------------------------------------------------------------------------------------
SUBROUTINE wnd_frc_thr_slt_get(PDNS_MDP, PDP, PWND_FRC_THR_SLT)  
! Purpose: Compute dry threshold friction velocity for saltation
!++grini use dstgrd ! [mdl] Dust grid sizes
!++grini use pmgrid ! [mdl] Spatial resolution parameters
!++grini use dstcst ! [mdl] Physical constants for dust routines
USE modd_dstmbl,    ONLY : xdns_slt
USE modd_csts,      ONLY : xg
USE modi_abor1_sfx
!
IMPLICIT NONE
!
! Input
!++grini real,intent(in)::dns_aer(dst_nbr) ! [kg m-3] Particle density
!++grini real,intent(in)::dmt_aer(dst_nbr) ! [m] Particle diameter
REAL, DIMENSION(:), INTENT(IN)   :: pdns_mdp        ! [kg m-3] Midlayer density
REAL, INTENT(IN)  :: pdp   
! Output
REAL, DIMENSION(:), INTENT(OUT) :: pwnd_frc_thr_slt ! [m s-1] Threshold friction velocity for saltation
!
! Local
REAL :: zryn, zryn1 ! [frc] Threshold friction Reynolds number approximation for optimal size
REAL :: zdns_fct    ! Density ratio factor for saltation calculation
REAL :: zicf_fct    ! Interparticle cohesive forces factor for saltation calculation
REAL :: ztmp        ! Factor in saltation computation
INTEGER :: i        ! [idx] Counting index
!
REAL(KIND=JPRB) :: zhook_handle   
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:WND_FRC_THR_SLT_GET',0,zhook_handle)
!
! Initialize some variables
! MaB95 pzn. for Re*t(D_opt) circumvents iterative solution
zryn = 0.38d0 + 1331.0d0*(100.0d0*pdp)**1.56d0                        ! [frc] "B" MaB95 p. 16417 (5)
! Given Re*t(D_opt), compute time independent factors contributing to u*t
zdns_fct                 = xdns_slt * xg * pdp                        ! IvW82 p. 115 (6) MaB95 p. 16417 (4)
zicf_fct                 = 1.0d0 + 6.0d-07/(xdns_slt*xg*(pdp**2.5d0)) ! [frc] IvW82 p. 115 (6) MaB95 p. 16417 (4) Interparticle cohesive forces
!
IF (zryn < 0.03d0) THEN
  CALL abor1_sfx('MODE_DSTMBLUTL:WND_FRC_THR_SLT_GET: RYN < 0.03')
ELSEIF (zryn < 10.0d0) THEN
  zryn1 = -1.0d0 + 1.928d0 * (zryn**0.0922d0) ! [frc] IvW82 p. 114 (3), MaB95 p. 16417 (6)
  zryn1 = 0.1291d0 * 0.1291d0 / zryn1          ! [frc] 
ELSE
  zryn1 = 1.0d0 - 0.0858d0 * exp(-0.0617d0*(zryn-10.0d0)) ! [frc] IvW82 p. 114 (3), MaB95 p. 16417 (7)
  zryn1 = 0.120d0 **2 * zryn1**2               ! [frc]
  !ryn1=0.129*0.129*ryn1*ryn1 dans le cas mm, à vérifier
ENDIF 
!
! This method minimizes the number of square root computations performed
ztmp = sqrt(zicf_fct*zdns_fct*zryn1)
DO i = 1, SIZE(pdns_mdp)
  pwnd_frc_thr_slt(i) = ztmp / sqrt(pdns_mdp(i))
ENDDO
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:WND_FRC_THR_SLT_GET',1,zhook_handle)
!
END SUBROUTINE wnd_frc_thr_slt_get
!
!----------------------------------------------------------------------------------------
SUBROUTINE vwc2gwc (OFLG_MBL, PVWC_SAT, PVWC_SFC, PGWC_SFC)
! Purpose: Convert volumetric water content to gravimetric water content
!++alfgr use pmgrid ! [mdl] Spatial resolution parameters
use modd_csts, only : xrholi  ! [kg/m3] density of liquid water
USE modd_dstmbl, ONLY : xdns_slt
!
!++alfgr use dstblm,only:dns_H2O_lqd_std ! [mdl] Boundary layer meteorology for non-vegetated land surfaces
IMPLICIT NONE
! Input
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl ! [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pvwc_sat ! [m3 m-3] Saturated volumetric water content (sand-dependent)
REAL,    DIMENSION(:), INTENT(IN) :: pvwc_sfc ! [m3 m-3] Volumetric water content
REAL,    DIMENSION(:), INTENT(OUT):: pgwc_sfc ! [kg kg-1] Gravimetric water content
! Local
REAL, DIMENSION(SIZE(PVWC_SAT)) :: zdns_blk_dry ! [kg m-3] Bulk density of dry surface soil
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle 
! Main Code
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:VWC2GWC',0,zhook_handle)
!
! Initialize output
DO i=1,SIZE(pvwc_sat)
  IF (oflg_mbl(i)) THEN
    ! Assume volume of air pores when dry equals saturated VWC
    ! This implies air pores are completely filled by water in saturated soil
    zdns_blk_dry(i) = xdns_slt * (1.0 - pvwc_sat(i)) ! [kg m-3] Bulk density of dry surface soil
    pgwc_sfc(i) = pvwc_sfc(i) * xrholi / zdns_blk_dry(i)
  ENDIF
ENDDO
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:VWC2GWC',1,zhook_handle)
END SUBROUTINE vwc2gwc
!
!----------------------------------------------------------------------------------------  
SUBROUTINE frc_thr_ncr_wtr_get(OFLG_MBL, PGWC_THR, PGWC_SFC, PFRC_THR_NCR_WTR)   
! Purpose: Compute factor by which soil moisture increases threshold friction velocity 
! This parameterization is based on FMB99
!++alfgr use pmgrid ! [mdl] Spatial resolution parameters
IMPLICIT NONE
! dans le calcul de l'humidit? du sol(a posteriori)
! Input
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl         ! I [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pgwc_thr         ! [kg kg-1] Threshold gravimetric water content
REAL,    DIMENSION(:), INTENT(IN) :: pgwc_sfc         ! [kg kg-1] Gravimetric water content
REAL,    DIMENSION(:), INTENT(OUT):: pfrc_thr_ncr_wtr ! [frc] Factor by which moisture increases threshold friction velocity
! Local
INTEGER :: i ! [idx] Counting index
REAL(KIND=JPRB) :: zhook_handle 

IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FRC_THR_NCR_WTR_GET',0,zhook_handle)
! Main Code
! Initialize output
pfrc_thr_ncr_wtr(:) = 1.0d0 ! [frc] Factor by which moisture increases threshold friction velocity
!
DO i = 1, SIZE(pgwc_sfc)
  IF (oflg_mbl(i)) THEN
    ! Adjust threshold velocity for inhibition by moisture
    ! frc_thr_ncr_wtr(lon_idx)=exp(22.7d0*vwc_sfc(lon_idx)) ! [frc] SRL96
    ! Compute threshold soil moisture based on clay content
    IF (pgwc_sfc(i) > pgwc_thr(i)) &
      pfrc_thr_ncr_wtr(i) = sqrt(1.0d0 + 1.21d0 * (100.0d0*(pgwc_sfc(i)-pgwc_thr(i)))**0.68d0) ! [frc] FMB99 p. 155 (15)
  ENDIF
ENDDO
! Uncomment following line to remove all dependence on gwc_sfc
! frc_thr_ncr_wtr(lon_idx)=1.0d0 ! [frc]
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FRC_THR_NCR_WTR_GET',1,zhook_handle)
END SUBROUTINE frc_thr_ncr_wtr_get
!  
!----------------------------------------------------------------------------------------  
!++alfgr fxm: Fix this so that we can actually use rgh_mmn_mbl different for grid cells
SUBROUTINE frc_thr_ncr_drg_get(PRGH_MMN_MBL, PRGH_MMN_SMT, PFRC_THR_NCR_DRG)
! Purpose: Compute factor by which surface roughness increases threshold friction velocity
! This parameterization is based on MaB95 and GMB98
!++alfgr use pmgrid ! [mdl] Spatial resolution parameters
USE modi_abor1_sfx
!
IMPLICIT NONE
! Input
REAL, INTENT(IN) :: prgh_mmn_mbl     ! [m] Roughness length momentum for erodible surfaces
REAL, INTENT(IN) :: prgh_mmn_smt                   ! [m] Smooth roughness length
! Output
REAL, INTENT(OUT):: pfrc_thr_ncr_drg ! [frc] Factor by which roughness increases threshold friction velocity
! Local
REAL :: zwnd_frc_fsh_frc     ! [frc] Efficient fraction of wind friction
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FRC_THR_NCR_DRG_GET',0,zhook_handle)
! Main Code
!
! Adjust threshold velocity for inhibition by roughness elements
zwnd_frc_fsh_frc = & ! [frc] MaB95 p. 16420, GMB98 p. 6207
  1.0d0 - log(prgh_mmn_mbl/prgh_mmn_smt) / log(0.35d0 * ((0.1d0/prgh_mmn_smt)**0.8d0))
zwnd_frc_fsh_frc = max(1.e-6, min(1., zwnd_frc_fsh_frc))
IF (zwnd_frc_fsh_frc <= 0.0d0 .OR. zwnd_frc_fsh_frc > 1.0d0) &
  CALL abor1_sfx("MODE_DSTMBLUTL:FRC_THR_NCR_DRG_GET0: WND_FRC_FSH_FRC OUT OF RANGE")
pfrc_thr_ncr_drg = 1.0d0 / zwnd_frc_fsh_frc! [frc]
! fxm: 19991012 Set frc_thr_ncr_drg=1.0, equivalent to assuming mobilization takes place at smooth roughness length
!++alfgr removed this frc_thr_ncr_drg=1.0d0
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FRC_THR_NCR_DRG_GET',1,zhook_handle)
END SUBROUTINE frc_thr_ncr_drg_get
!
!----------------------------------------------------------------------------------------
SUBROUTINE wnd_frc_slt_get(OFLG_MBL, PWND_FRC, PWND_RFR, PWND_RFR_THR_SLT, PWND_FRC_SLT)
! Purpose: Compute the saltating friction velocity
! Saltation increases friction speed by roughening surface, AKA "Owen's effect"
! This acts as a positive feedback to the friction speed
! GMB98 parameterized this feedback in terms of 10 m windspeeds
!++alfgr use pmgrid ! [mdl] Spatial resolution parameters
IMPLICIT NONE
! Input
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl         ! I [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_frc         ! I [m s-1] Surface friction velocity
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_rfr         ! I [m s-1] Wind speed at reference height
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_rfr_thr_slt ! I [m s-1] Threshold 10 m wind speed for saltation
REAL,    DIMENSION(:), INTENT(OUT):: pwnd_frc_slt     ! O [m s-1] Saltating friction velocity
! Local
REAL :: zwnd_rfr_dlt     ! [m s-1] Reference windspeed excess over threshold
REAL :: zwnd_frc_slt_dlt ! [m s-1] Friction velocity increase from saltation
INTEGER :: i             ! [idx] Counting index
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:WND_FRC_SLT_GET',0,zhook_handle)   
! Main Code
! Compute saltating friction velocity, accounting for "Owen's effect"
pwnd_frc_slt(:) = pwnd_frc(:) ! [m s-1] Saltating friction velocity
!
DO i = 1, SIZE(pwnd_frc)
  IF (oflg_mbl(i) .AND. pwnd_rfr(i) >= pwnd_rfr_thr_slt(i)) THEN
    ! Saltation roughens the boundary layer, AKA "Owen's effect"
    ! GMB98 p. 6206 Fig. 1 shows observed/computed u* dependence on observed U(1 m)
    ! GMB98 p. 6209 (12) has u* in cm s-1 and U, Ut in m s-1, personal communication, D. Gillette, 19990529
    ! With everything in MKS, the 0.3 coefficient in GMB98 (12) becomes 0.003 
    ! Increase in friction velocity due to saltation varies as square of 
    ! difference between reference wind speed and reference threshold speed 
    zwnd_rfr_dlt     = pwnd_rfr(i) - pwnd_rfr_thr_slt(i)
    zwnd_frc_slt_dlt = 0.003d0 * zwnd_rfr_dlt**2          ! [m s-1] Friction velocity increase from saltation GMB98 p. 6209
    pwnd_frc_slt(i)  = pwnd_frc_slt(i) + zwnd_frc_slt_dlt ! [m s-1] Saltating friction velocity
  ENDIF
ENDDO
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:WND_FRC_SLT_GET',1,zhook_handle)
END SUBROUTINE wnd_frc_slt_get
!
!----------------------------------------------------------------------------------------
SUBROUTINE flx_mss_hrz_slt_ttl_whi79_get(PCOEFF, OFLG_MBL, PDNS_MDP, PWND_FRC, &
                        pwnd_frc_thr_slt, pflx_mss_hrz_slt_ttl)
! Purpose: Compute vertically integrated streamwise mass flux of particles
! Theory: Uses method proposed by White (1979)
! fxm: use surface air density not midlayer density
!++alfgr use pmgrid ! [mdl] Spatial resolution parameters
!++alfgr use dstcst ! [mdl] Physical constants for dust routines
use modd_csts,      ONLY : xg  ! Gravitation constant
!
IMPLICIT NONE
! Input
REAL,    DIMENSION(:), INTENT(IN) :: pcoeff
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl             ! I [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pdns_mdp             ! I [kg m-3] Midlayer density
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_frc             ! I [m s-1] Surface friction velocity
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_frc_thr_slt     ! I [m s-1] Threshold friction speed for saltation
REAL,    DIMENSION(:), INTENT(OUT):: pflx_mss_hrz_slt_ttl ! O [kg m-1 s-1] Vertically integrated streamwise mass flux
! Local
REAL :: zwnd_frc_rat ! [frc] Ratio of wind friction threshold to wind friction
INTEGER :: i    ! [idx] Counting index for lon
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_HRZ_SLT_TTL_WHI79_GET',0,zhook_handle)   
! Main Code
! Initialize output
pflx_mss_hrz_slt_ttl(:) = 0.0d0 ! [kg m-1 s-1]
!
DO i = 1, SIZE(pdns_mdp)
  IF (oflg_mbl(i) .AND. pwnd_frc(i) > pwnd_frc_thr_slt(i)) THEN
    zwnd_frc_rat = pwnd_frc_thr_slt(i) / pwnd_frc(i)            ! [frc]
    pflx_mss_hrz_slt_ttl(i) =                                 & ! [kg m-1 s-1]
      pcoeff(i) * pdns_mdp(i) * (pwnd_frc(i)**3.0d0) *      &
      (1.0d0 - zwnd_frc_rat) * (1.0d0 + zwnd_frc_rat)**2 / xg   ! Whi79 p. 4648 (19), MaB97 p. 16422 (28)
  ENDIF
ENDDO

IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_HRZ_SLT_TTL_WHI79_GET',1,zhook_handle) 
END SUBROUTINE flx_mss_hrz_slt_ttl_whi79_get
! 
!----------------------------------------------------------------------------------------
SUBROUTINE flx_mss_vrt_dst_ttl_mab95_get(OFLG_MBL, PMSS_FRC_CLY, PFLX_MSS_HRZ_SLT_TTL, & 
                        pdst_slt_flx_rat_ttl, pflx_mss_vrt_dst_ttl)
! Purpose: Diagnose total vertical mass flux of dust from vertically integrated streamwise mass flux
! Theory: Uses clay-based method proposed by Marticorena & Bergametti (1995)
! Their parameterization is based only on data for mss_frc_cly < 0.20
! For clayier soils, dst_slt_flx_rat_ttl may behave dramatically differently
! Whether this behavior changes when mss_frc_cly > 0.20 is unknown
! Anecdotal evidence suggests vertical flux decreases for mss_frc_cly > 0.20
! Thus we use min[mss_frc_cly,0.20] in MaB95 parameterization
IMPLICIT NONE
! Input
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl             ! I [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pmss_frc_cly         ! I [frc] Mass fraction clay
REAL,    DIMENSION(:), INTENT(IN) :: pflx_mss_hrz_slt_ttl ! I [kg m-1 s-1] Vertically integrated streamwise mass flux
REAL,    DIMENSION(:), INTENT(OUT):: pdst_slt_flx_rat_ttl ! O [m-1] Ratio of vertical dust flux to streamwise mass flux
REAL,    DIMENSION(:), INTENT(OUT):: pflx_mss_vrt_dst_ttl ! O [kg m-2 s-1] Total vertical mass flux of dust
! Local
REAL :: zmss_frc_cly_vld ! [frc] Mass fraction clay limited to 0.20
INTEGER :: i             ! [idx] Counting index for lon
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_VRT_DST_TTL_MAB95_GET',0,zhook_handle)  
!
DO i = 1, SIZE(pmss_frc_cly)
  IF (oflg_mbl(i)) THEN
    ! 19990603: fxm: Dust production is EXTREMELY sensitive to this parameter, which changes flux by 3 orders of magnitude in 0.0 < mss_frc_cly < 0.20
    zmss_frc_cly_vld = min(pmss_frc_cly(i), 0.2) ! [frc]
    pdst_slt_flx_rat_ttl(i) = & ! [m-1]
      100.0d0 * exp(log(10.0d0)*(13.4d0*zmss_frc_cly_vld - 6.0d0))              ! MaB95 p. 16423 (47)
    pflx_mss_vrt_dst_ttl(i) = pflx_mss_hrz_slt_ttl(i) * pdst_slt_flx_rat_ttl(i) ! [kg m-1 s-1]
  ENDIF
ENDDO
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_VRT_DST_TTL_MAB95_GET',1,zhook_handle) 
END SUBROUTINE flx_mss_vrt_dst_ttl_mab95_get
!
!----------------------------------------------------------------------------------------
SUBROUTINE flx_mss_vrt_dst_aust_get(OFLG_MBL, PDNS_MDP, PWND_FRC_THR_SLT, &
        pflx_mss_hrz_slt_ttl, pdst_slt_flx_rat_ttl, pflx_mss_vrt_dst_ttl)
! Purpose: Diagnose total vertical mass flux of dust from vertically integrated streamwise mass flux
use modd_csts,      ONLY : xg  ! Gravitation constant
USE modd_dstmbl, ONLY : xdmt_slt_opt, xdmt_ero_opt, xdns_slt, xgama
!
IMPLICIT NONE
! Input
LOGICAL, DIMENSION(:), INTENT(IN) :: oflg_mbl             ! I [flg] Mobilization candidate flag
REAL,    DIMENSION(:), INTENT(IN) :: pdns_mdp             ! I [kg m-3] Midlayer density
REAL,    DIMENSION(:), INTENT(IN) :: pwnd_frc_thr_slt     ! I [m s-1] Threshold friction speed for saltation
REAL,    DIMENSION(:), INTENT(IN) :: pflx_mss_hrz_slt_ttl ! I [kg m-1 s-1] Vertically integrated streamwise mass flux
REAL,    DIMENSION(:), INTENT(OUT):: pdst_slt_flx_rat_ttl ! O [m-1] Ratio of vertical dust flux to streamwise mass flux
REAL,    DIMENSION(:), INTENT(OUT):: pflx_mss_vrt_dst_ttl ! O [kg m-2 s-1] Total vertical mass flux of dust
! Local
REAL, DIMENSION(SIZE(PDNS_MDP)) :: zrop_roa  ! [frc] ratio of particle densite to Midlayer density
REAL :: zexpdd, zlnds, zbeta, zbgxg
INTEGER :: i
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_VRT_DST_AUST_GET',0,zhook_handle)  
!
! Initialize some variables
zexpdd = exp(-140.7d0 * xdmt_ero_opt + 0.37d0)
zlnds  = (0.328*1.0d-4) + (0.125d0*1.0d-4*log(xdmt_slt_opt*1.0e+3))
zbeta  = zlnds * zexpdd ! [mm]
zbgxg  = zbeta * xgama * xg
!
DO i = 1, SIZE(pdns_mdp)
  IF (oflg_mbl(i)) THEN
    zrop_roa(i) = xdns_slt / pdns_mdp(i) 
    pdst_slt_flx_rat_ttl(i) =          & ! [mm/m]
      (2.0d0/3.0d0) * zbgxg * zrop_roa(i) / (pwnd_frc_thr_slt(i)**2.0d0)
    pflx_mss_vrt_dst_ttl(i) = 1.0d-3 * & ! [kg m-1 s-1]
      pflx_mss_hrz_slt_ttl(i) * pdst_slt_flx_rat_ttl(i)
  ENDIF
ENDDO
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:FLX_MSS_VRT_DST_AUST_GET',1,zhook_handle)
END SUBROUTINE flx_mss_vrt_dst_aust_get
!
!----------------------------------------------------------------------------------------
SUBROUTINE distribution(KMODE, PDLNDP, KTEXT, PDP, PDSRLV, PZS0)
!
USE modd_csts,      ONLY : xpi
USE modd_dstmbl, ONLY : xdns_slt
!
IMPLICIT NONE
!
!INPUT, set their dimensions to their passed lengths or to KSIZE ?
INTEGER, INTENT(IN) :: kmode  ![nbr] loop over mode
REAL,    INTENT(IN) :: pdlndp ! delta Dp [?m]
INTEGER, DIMENSION(:), INTENT(IN) :: ktext  ! texture
REAL, DIMENSION(:),   INTENT(OUT) :: pdp
REAL, DIMENSION(:,:), INTENT(OUT) :: pdsrlv ! [--] Output surface relative
REAL, DIMENSION(:),   INTENT(OUT) :: pzs0   ! [m] Output rugosit? de la surface lisse (Dp/30)
!Define local variables:
REAL, DIMENSION(KMODE,12) :: zsigma ! déviation géométrique standard du mode
REAL, DIMENSION(KMODE,12) :: zmfrac ! fraction massique du mode
REAL, DIMENSION(KMODE,12) :: zdmed  ! diametre median du mode
REAL, DIMENSION(SIZE(PDSRLV,1),SIZE(PDSRLV,2)) :: zds
REAL, DIMENSION(SIZE(PDSRLV,2)) :: zdpln, zsp, zdmln
REAL :: zdm1, zdm2
INTEGER :: imod  ! Counter for number of mode
INTEGER :: idp   ! Counter for number of particle
INTEGER :: itex  ! Counter for number of texture
REAL(KIND=JPRB) :: zhook_handle 
!
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:DISTRIBUTION',0,zhook_handle)  
!
!1 Sand
!2 Loamy Sand
!3 Sandy Loam
!4 Silt Loam
!5 Loam
!6 Sandy Clay Loam
!7 Silty Clay Loam
!8 Clay Loam
!9 Sandy Clay
!10 Silty Clay
!11 Clay
!12 Silt
zmfrac(1,:) = (/0.0, 0.1, 0.1, 0.15, 0.15, 0.2, 0.2, 0.3, 0.35, 0.4, 0.5, 0.15/)
zmfrac(2,:) = (/0.1, 0.3, 0.3, 0.35, 0.50, 0.5, 0.5, 0.5, 0.00, 0.0, 0.0, 0.40/)
zmfrac(3,:) = (/0.9, 0.6, 0.6, 0.50, 0.35, 0.3, 0.3, 0.2, 0.65, 0.6, 0.5, 0.45/)
!
zsigma(1,:) = (/1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8/)
zsigma(2,:) = (/1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.7/)
zsigma(3,:) = (/1.6, 1.6, 1.6, 1.6, 1.6, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.6/)
!
zdmed(1,:) = (/  10., 10.,  5.,  5., 2.5, 2.5, 2.5,  1.,  1., 0.5, 0.5, 2.5/)
zdmed(2,:) = (/ 100.,100.,100.,100., 75., 75., 50., 50., 10., 10., 10., 75./)
zdmed(3,:) = (/1000.,690.,520.,520.,520.,210.,210.,125.,100.,100.,100.,520./)
!
DO itex = 1, SIZE(pzs0) 
  pzs0(itex) =1e-6*zdmed(3,itex)/30.0
ENDDO
!
! Initialisation pour IDP = 1
pdp(1) = 0.1          ! ?m
zdpln(1) = log(pdp(1))  ! ?m
! 
DO idp = 2, SIZE(pdsrlv,2)
  zdpln(idp) = zdpln(idp-1) + pdlndp
  pdp(idp) = exp(zdpln(idp))
END DO
!
! calcul pour chaque particule 
! Boucle NDP 
! Initialisation
! surface basale totale
zsp(:) = 0.d0
! Calcul de la distribution massique des particules (dm(Dp)/dln(Dp)) pour chaque particule de diametre Dp
DO itex = 1, SIZE(pdsrlv,1)
  DO idp = 1, SIZE(pdsrlv,2)
    zdmln(idp) = 0.0d0
    DO imod = 1, kmode
      zdm1 = zmfrac(imod,itex) / (sqrt(2.d0*xpi) * log(zsigma(imod,itex)))
      zdm2 = exp((log(pdp(idp)) - log(zdmed(imod,itex)))**2. / (-2.d0 * (log(zsigma(imod,itex)))**2.d0))
      zdmln(idp) = zdmln(idp) + zdm1*zdm2
    END DO 
    ! Calcul de la surface basale occup?e par chaque particule de diam?tre Dp  
    zds(itex,idp) = 3.d0 * zdmln(idp) * pdlndp / (2.* xdns_slt * pdp(idp))
    ! Calcul de la surface basale totale Sp= (somme (ds(Dp)):
    zsp(itex) = zsp(itex) + zds(itex,idp)
  ENDDO
ENDDO 
!
! Calcul de la distribution de surface relative des particules
DO itex = 1, SIZE(pzs0)
  DO idp = 1, SIZE(pdsrlv,2)
    pdsrlv(itex,idp) = zds(itex,idp) / zsp(itex)
    ENDDO 
END DO 
! 
IF (lhook) CALL dr_hook('MODE_DSTMBLUTL:DISTRIBUTION',1,zhook_handle)  
END SUBROUTINE distribution
!
END MODULE mode_dstmblutl