ch_aer_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 ch_aer_dep (PSVT, PFSVT,  PUSTAR, &
                      presa, pta, prhodref)  
!###########################################################
  !
  !!                   
  !!                       
  !!
  !!    PURPOSE
  !!    -------
  !!      
  !!    Compute dry deposition velocity for aerosol species 
  !!
  !!    AUTHOR
  !!    ------
  !!      P.Tulet      * CNRM *
  !!
  !!    MODIFICATIONS
  !!    -------------
  !!      Original      20/02/05
  !!      J.Escobar     06/2013  for REAL4/8 add EPSILON management
  !!
  !-------------------------------------------------------------------------------
  !
  !*       0.    DECLARATIONS
  !              ------------
  !
  USE mode_aer_surf
  USE modi_ch_aer_velgrav1d
  USE modd_chs_aerosol
  !
!
  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
  !
  !
  !*       0.2   Declarations of local variables :
  !
  REAL , DIMENSION(SIZE(PSVT,1), JPIN) :: zrd ! surface  resistance
  REAL , DIMENSION(SIZE(PSVT,1), JPIN) :: zvd  
  REAL , DIMENSION(SIZE(PSVT,1), JPIN) :: stn ! Stockes number
  REAL , DIMENSION(SIZE(PSVT,1), JPIN) :: sc
  REAL , DIMENSION(SIZE(PSVT,1))      :: wcn
  REAL , DIMENSION(SIZE(PSVT,1))      :: zustar, zresa
  REAL , DIMENSION(SIZE(PSVT,1), JPIN):: zwork
  REAL , DIMENSION(SIZE(PSVT,1),NSP+NCARB+NSOA,JPMODE):: zctota, zcctot
  REAL, DIMENSION(SIZE(PSVT,1),JPMODE):: zrhop
  REAL, DIMENSION(SIZE(PSVT,1)) :: znu
  REAL, DIMENSION(SIZE(PSVT,1),JPIN)     :: dg,zvs,zvsg, zdsg
  REAL, DIMENSION(SIZE(PSVT,1))        :: zmu
  REAL, DIMENSION(SIZE(PSVT,1),JPIN)   :: zvgk, zdpk
  REAL, DIMENSION(SIZE(PSVT,1),JPMODE) :: zsig, zrg, zn
  REAL, DIMENSION(SIZE(PSVT,1),JPMODE) :: zvg, zdg
  REAL, DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2)) :: zsvt
  REAL, DIMENSION(SIZE(PSVT,1))        :: zsum
  REAL, DIMENSION(NSP+NCARB+NSOA)      :: zrhoi
  INTEGER :: jt, jj, jsv, jn
  INTEGER :: m6i, m6j
  REAL :: zden2mol, zg, ztmp1, ztmp2, ztmp3, ztmp4
  REAL(KIND=JPRB) :: zhook_handle


  !
  !============================================================================
  !
  !            Primilary
  !            ---------
  !Default values
  !--------------
! Cf Ackermann (all to black carbon except water)
IF (lhook) CALL dr_hook('CH_AER_DEP',0,zhook_handle)
!
zrhoi(:) = 1.8e3
zrhoi(jp_aer_h2o) = 1.0e3   ! water

zden2mol = 1e-6 * xavogadro / xmd
zg       = 9.80665
zmu(:) = 0.
zvgk(:,:) = 0.
zvg(:,:) = 0.
zdpk(:,:) = 0.
zustar(:) = max(pustar(:), 1.e-20)
zresa(:) = min(max(presa(:),   1.e-20), 9999.)
! molec./m3 to part/part
DO jsv=1,SIZE(psvt,2)
zsvt(:,jsv) = psvt(:,jsv) * xmd / (xavogadro * prhodref(:))
ENDDO
zsvt(:,:) = max(zsvt(:,:),xsurf_tiny)
 CALL ppp2aero_surf(zsvt, prhodref, psig1d=zsig, prg1d=zrg, pn1d=zn, pctota=zctota)
zrhop(:,:) = 0.
DO jn=1,jpmode
  zsum(:)=0.
  DO jj=1,nsp+ncarb+nsoa
   zsum(:)=zsum(:)+zctota(:,jj,jn)/zrhoi(jj)
  END DO

  DO jj=1,nsp+ncarb+nsoa
   zcctot(:,jj,jn) = zctota(:,jj,jn)/zrhoi(jj)/zsum(:)
   zrhop(:,jn)=zrhop(:,jn)+zcctot(:,jj,jn)*zrhoi(jj)
  ENDDO
ENDDO
 CALL ch_aer_velgrav1d(zsig, zrg, pta, prhodref, zrhop, zmu, zvgk,zdpk, zvg, zdg)
dg(:,:)  = max(zdpk(:,:),1.e-40)
zvs(:,:) = max(zvgk(:,:),1.e-20)
znu(:)   = zmu(:)/prhodref(:)
DO jn=1,jpmode
DO jj= 0,2
zvsg(:,3*jn+jj-2) =  max(zvg(:,jn),1.e-20)
zdsg(:,3*jn+jj-2) =  max(zdg(:,jn),1.e-40)
END DO
END DO
!     compute Schmidt number
!     ----------------------
  DO  jn=1,jpin
      !Sc(:,JN)= ZNU(:)/Dg(:,JN)
      sc(:,jn)= znu(:)/zdsg(:,jn)
  END DO
!Scale for convective velocity 
! WCn(:) = MAX((PTKE(:) - 4.65* ZUSTAR(:)**2)/0.3, 1.E-20)


        ! verifier l'echelle convective sinon laisser la formulation de seinfeld
!    WHERE (WCn(:,:) /= XUNDEF)
!     WCn(:,:) = SQRT(WCn(:,:))
!    ELSEWHERE
     wcn(:) = 0.
!    END WHERE
!
!
stn(:,:) =0.
zvd(:,:) = 0.
zwork(:,:) = 0.
DO jt=1,SIZE(psvt,1)
  IF (zustar(jt).GE.1.e-10) THEN
    DO jn=1,jpin
      ztmp1=0.
      ztmp2=0.
      ztmp3=0.
      ztmp4=0.  
      stn(jt,jn)= zvsg(jt,jn)*zustar(jt)**2/(zg*znu(jt))
      ztmp1=sc(jt,jn)**(-2./3.)
      ztmp2=(-3./stn(jt,jn))
      IF (ztmp2.gt.-10) then
        ztmp3=10.**(ztmp2)
      ELSE
        ztmp3=0.
      ENDIF
      ztmp4=ztmp1+ztmp3
     
     !ZRD(:,JN) = (Sc(:,JN)**(-2./3.)+ 10**(-3./Stn(:,JN)))&
     !      * (1 + 0.24 * WCn(:)**2 /ZUSTAR(:)**2) &
     !      * ZUSTAR(:)  
     !ZRD(:,JN) = ZUSTAR(:) * (Sc(:,JN)**(-2./3.)+ &
     !                           10**(-3./Stn(:,JN)))   
     zrd(jt,jn) = zustar(jt) * ztmp4                                
     zrd(jt,jn) = max(zrd(jt,jn),1.e-10) 
     zrd(jt,jn) = 1. / zrd(jt,jn) 
     zwork(jt,jn)= zresa(jt) + zrd(jt,jn) + zresa(jt)*zrd(jt,jn)*zvs(jt,jn)  
     zwork(jt,jn)= max(zwork(jt,jn), 1.e-10)
     zwork(jt,jn)= zvs(jt,jn) + 1./ zwork(jt,jn)
        !         deposition velocity for each cover type
        !         ----------------------------------------
     zvd(jt,jn) = zvd(jt,jn) + zwork(jt,jn)
  END DO
 ELSE
   zvd(jt,:) = 0.
 END IF
ENDDO  

m6i=0
m6j=0
IF (lvarsigi) m6i=1
IF (lvarsigj) m6j=1
DO jsv=1,SIZE(psvt,2)-1-(jpmode+m6i+m6j),2 ! mass deposition for I mode
  pfsvt(:,jsv) =  pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,2)
ENDDO
DO jsv=2,SIZE(psvt,2)-(jpmode+m6i+m6j),2   ! mass deposition for J mode
  pfsvt(:,jsv) =  pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,5)
ENDDO
! number particles deposition I
jsv = SIZE(psvt,2)-(1+m6i+m6j)
pfsvt(:,jsv) = pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,1)
! number particles deposition J
jsv = SIZE(psvt,2)-(m6i+m6j)
pfsvt(:,jsv) = pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,4)
! m6 deposition I
jsv = SIZE(psvt,2)-m6j
IF (lvarsigi) pfsvt(:,jsv) = pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,3)
! m6 deposition J
jsv = SIZE(psvt,2)
IF (lvarsigj) pfsvt(:,jsv) = pfsvt(:,jsv) - psvt(:,jsv)  * zvd(:,6)
IF (lhook) CALL dr_hook('CH_AER_DEP',1,zhook_handle)
!
!---------------------------------------------------------------------
!
END SUBROUTINE ch_aer_dep