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