dslt_dep.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 dslt_dep (PSVT, PFSVT, PUSTAR, PRESA, PTA, PRHODREF,             &
                     PEMISSIG, PEMISRADIUS, KPMODE, PDENSITY, PMOLARWEIGHT, &
                     PCONVERTFACM0, PCONVERTFACM6, PCONVERTFACM3,           &
                     OVARSIG, ORGFIX, HVERMOD                               )  
!###########################################################                    
!!
!!    PURPOSE
!!    -------
!!      
!!    Compute dry deposition velocity for dust species 
!!
!!    AUTHOR
!!    ------
!!      P.Tulet      * CNRM *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original      20/02/05
!!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_csts, ONLY : xpi, xavogadro, xg
!
USE mode_dslt_surf
USE modi_dslt_velgrav1d
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
USE modd_surf_par , ONLY : xsurf_tiny
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
REAL, DIMENSION(:,:),   INTENT(IN)    :: PSVT       ! friction velocity
REAL, DIMENSION(:,:),   INTENT(INOUT) :: PFSVT      ! flux
REAL, DIMENSION(:),     INTENT(IN)    :: PUSTAR       ! friction velocity
REAL, DIMENSION(:),     INTENT(IN)    :: PRESA        ! aerodynamical resistance
REAL, DIMENSION(:),     INTENT(IN)    :: PTA          ! ait temperature
REAL, DIMENSION(:),     INTENT(IN)    :: PRHODREF     ! air density
REAL, DIMENSION(:),     INTENT(IN)    :: PEMISSIG  
REAL, DIMENSION(:),     INTENT(IN)    :: PEMISRADIUS
INTEGER,                INTENT(IN)    :: KPMODE
REAL,                   INTENT(IN)    :: PDENSITY
REAL,                   INTENT(IN)    :: PMOLARWEIGHT
REAL,                   INTENT(OUT)   :: PCONVERTFACM0
REAL,                   INTENT(OUT)   :: PCONVERTFACM6
REAL,                   INTENT(OUT)   :: PCONVERTFACM3
LOGICAL,                INTENT(IN)    :: OVARSIG
LOGICAL,                INTENT(IN)    :: ORGFIX
 CHARACTER(LEN=6),       INTENT(IN)    :: HVERMOD
!
!*       0.2   Declarations of local variables :
!
REAL , DIMENSION(SIZE(PSVT,1), KPMODE) :: ZSIG, ZRG, ZVG, ZDG
REAL , DIMENSION(SIZE(PSVT,1), KPMODE) :: ZSTN ! Stockes number
REAL , DIMENSION(SIZE(PSVT,1), KPMODE) :: ZSC  ! Schmidt number
REAL , DIMENSION(SIZE(PSVT,1), KPMODE) :: ZRD  ! surface  resistance
REAL , DIMENSION(SIZE(PSVT,1), KPMODE*3) :: ZVGK, ZDPK
REAL , DIMENSION(SIZE(PSVT,1), KPMODE*3) :: ZVD  ! [m/s] dry deposition velocity 
REAL , DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2)) :: ZSVT
REAL , DIMENSION(SIZE(PSVT,1)) :: ZUSTAR, ZRESA
REAL , DIMENSION(SIZE(PSVT,1)) :: ZNU
REAL , DIMENSION(SIZE(PSVT,1)) :: ZMU
INTEGER,DIMENSION(KPMODE) :: NM0                 ! [idx] index for Mode 0 in passed variables
INTEGER,DIMENSION(KPMODE) :: NM3                 ! [idx] indexes for Mode 3 in passed variables
INTEGER,DIMENSION(KPMODE) :: NM6                 ! [idx] indexes for Mode 6 in passed variables
INTEGER :: JN, J0
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!============================================================================
!
!            Primilary
!            ---------
!Default values
!--------------
! Cf Ackermann (all to black carbon except water)
IF (lhook) CALL dr_hook('DSLT_DEP',0,zhook_handle)
!
zustar(:) = max(pustar(:), 1.e-20)
zresa(:) = min(max(presa(:), 1.e-20), 9999.)
! Save scalars in local array
zsvt(:,:) = max(psvt(:,:), xsurf_tiny)
!
zmu(:)    = 0.
zvgk(:,:) = 0.
zvg(:,:) = 0.
zdpk(:,:) = 0.
!
IF (ovarsig) THEN
  DO jn=1,kpmode
    nm0(jn) = 1+(jn-1)*3
    nm3(jn) = 2+(jn-1)*3
    nm6(jn) = 3+(jn-1)*3
  END DO
ELSE IF (orgfix) THEN
  DO jn=1,kpmode
    nm3(jn) = jn
  END DO
ELSE
  DO jn=1,kpmode
    nm0(jn) = 1+(jn-1)*2
    nm3(jn) = 2+(jn-1)*2
  END DO
END IF
!
pconvertfacm0 = pmolarweight / xavogadro
pconvertfacm6 = pconvertfacm0 * 1.d6
pconvertfacm3 = 4./3. * xpi * pdensity / 1.d18
!
 CALL dsltmoment2size(zsvt, prhodref, pemisradius, pemissig, nm0, nm3, nm6, &
                     pconvertfacm0, pconvertfacm6, pconvertfacm3,          &
                     ovarsig, orgfix, psig1d=zsig, prg1d=zrg               )
!
 CALL dslt_velgrav1d(zsig, zrg, pta, prhodref, pdensity, zmu, zvgk, zdpk, zvg, zdg)
!
znu(:)   = zmu(:)/prhodref(:)
!
zvgk(:,:) = max(zvgk(:,:),1.e-20)
zdpk(:,:) = max(zdpk(:,:),1.e-40)
!
zvg(:,:) = max(zvg(:,:),1.e-20)
zdg(:,:) = max(zdg(:,:),1.e-40)
!
zstn(:,:) =0.
!
DO  jn = 1,kpmode        
  !
  !         deposition velocity for each cover type
  !         ----------------------------------------
  !Stoke's number, Seinfeld & Pandis, pp 965
  zstn(:,jn) = zvg(:,jn)*zustar(:)**2/(xg*znu(:))
  zstn(:,jn) = max(zstn(:,jn), 0.05)
  !
  !     compute Schmidt number
  !     ----------------------
  zsc(:,jn) = znu(:)/zdg(:,jn)
  !
  !Get nominator of equation 19.18 Seinfeld & Pandis==> 1/rd
  zrd(:,jn) = zustar(:) * (zsc(:,jn)**(-2./3.)+ 10**(-3./zstn(:,jn)))
  !
  !Limit to reasonable values
  !ZRD(:,JN) = MAX(ZRD(:,JN),1.E-10)
  !Get rd
  zrd(:,jn) = 1. / zrd(:,jn)
  !
ENDDO
!
DO  jn = 1,kpmode*3
  !
  j0 = (jn-1)/3+1
  !
  !         deposition velocity for each cover type
  !         ----------------------------------------
  zvd(:,jn) = 0.
  !
  !Get ra + rd + ra*rd*vg which is equal to 
  !getting nominator of equation 19.7 Seinfeld & Pandis 
  zvd(:,jn)= zresa(:) + zrd(:,j0) + zresa(:)*zrd(:,j0)*zvgk(:,jn)
  !
  !Limit to reasonable values
  zvd(:,jn)= max(zvd(:,jn), 1.e-10)
  !
  !Get the total dry dep velocity (Seinfeld & Pandis, eqn 19.7)
  IF (hvermod=='CMDVER') THEN 
    zvd(:,jn)= zvgk(:,jn)         &  ! Gravitation term 
               + 1./zvd(:,jn)        ! turbulence and surface resistance term
  ELSE    
!     ZVD(:,JN)=zvs(:,JN)  &  !Gravitation term 
!         + 1./ZVD(:,JN)     !turbulence and surface resistance term
      ! The gravitation term as been computed by MesoNH (see sedim_dust.f90)
    zvd(:,jn) =  1./zvd(:,jn) ! turbulence and surface resistance term
  END IF
  !
END DO

! Only M3 flux (mass) has been used ; flux for over moment has been made after
DO jn = 1,kpmode
  pfsvt(:,nm3(jn)) =  pfsvt(:,nm3(jn)) - psvt(:,nm3(jn))  * zvd(:,2+(jn-1)*3)
ENDDO
!
IF (lhook) CALL dr_hook('DSLT_DEP',1,zhook_handle)
!---------------------------------------------------------------------
!
END SUBROUTINE dslt_dep