source: LMDZ6/branches/contrails/libf/phylmd/Dust/inscav_spl.f90 @ 5445

SUBROUTINE inscav_spl(pdtime,it,masse,henry,kk,qliq, &
        flxr,flxs,zrho,zdz,t,x, &
        his_dh)
  USE dimphy
  USE yomcst_mod_h
  USE yoecumf_mod_h
  USE chem_mod_h
IMPLICIT NONE
  !=====================================================================
  ! Objet : depot humide de traceurs
  ! Date : mars 1998
  ! Auteur: O. Boucher (LOA)
  !=====================================================================

  !
  INTEGER :: it
  REAL :: pdtime              ! pas de temps (s)
  REAL :: masse               ! molar mass (except for BC/OM/IF/DUST=Nav)
  REAL :: henry               ! constante de Henry en mol/l/atm
  REAL :: kk                  ! coefficient de dependence en T (K)
  REAL :: qliq                ! contenu en eau liquide dans le nuage (kg/kg)
   ! REAL flxr(klon,klev+1)   ! flux precipitant de pluie
   ! REAL flxs(klon,klev+1)   ! flux precipitant de neige
  REAL :: flxr(klon,klev)   ! flux precipitant de pluie   ! Titane
  REAL :: flxs(klon,klev)   ! flux precipitant de neige   ! Titane
  REAL :: flxr_aux(klon,klev+1)
  REAL :: flxs_aux(klon,klev+1)
  REAL :: zrho(klon,klev)
  REAL :: zdz(klon,klev)
  REAL :: t(klon,klev)
  REAL :: x(klon,klev)        ! q de traceur
  REAL :: his_dh(klon)        ! tendance de traceur integre verticalement
  !
  !--variables locales
  INTEGER :: i, k
  !
  REAL :: dx      ! tendance de traceur
  REAL :: f_a     !--rapport de la phase aqueuse a la phase gazeuse
  REAL :: beta    !--taux de conversion de l'eau en pluie
  REAL :: henry_t !--constante de Henry a T t  (mol/l/atm)
  REAL :: scav(klon,klev)    !--fraction aqueuse du constituant
  REAL :: K1, K2, ph, frac
  REAL :: frac_gas, frac_aer !-cste pour la reevaporation
  PARAMETER (ph=5., frac_gas=1.0, frac_aer=0.5)
  !---cste de dissolution pour le depot humide
  REAL :: frac_fine_scav,frac_coar_scav
  !---added by nhl
  REAL :: aux_cte

  PARAMETER (frac_fine_scav=0.7)
  PARAMETER (frac_coar_scav=0.7)

  !--101.325  m3/l x Pa/atm
  !--R        Pa.m3/mol/K
  !
  !------------------------------------------
  !
  !nhl      IF (it.EQ.2.OR.it.EQ.3) THEN !--aerosol  ! AS IT WAS FIRST
  IF (it.EQ.2.OR.it.EQ.3.OR.it.EQ.4) THEN !--aerosol
    frac=frac_aer
  ELSE                                                !--gas
    frac=frac_gas
  ENDIF
  !
  IF (it.EQ.1) THEN
  DO k=1, klev
  DO i=1, klon
    henry_t=henry*exp(-kk*(1./298.-1./t(i,k)))    !--mol/l/atm
    K1=1.2e-2*exp(-2010*(1/298.-1/t(i,k)))
    K2=6.6e-8*exp(-1510*(1/298.-1/t(i,k)))
    henry_t=henry_t*(1 + K1/10.**(-ph) + K1*K2/(10.**(-ph))**2)
    f_a=henry_t/101.325*R*t(i,k)*qliq*zrho(i,k)/rho_water
    scav(i,k)=f_a/(1.+f_a)
  ENDDO
  ENDDO
  ELSEIF (it.EQ.2) THEN
  DO k=1, klev
  DO i=1, klon
    scav(i,k)=frac_fine_scav
  ENDDO
  ENDDO
  ELSEIF (it.EQ.3) THEN
  DO k=1, klev
  DO i=1, klon
    scav(i,k)=frac_coar_scav
  ENDDO
  ENDDO
  ELSEIF (it.EQ.4) THEN
  DO k=1, klev
  DO i=1, klon
    scav(i,k)=frac_coar_scav
  ENDDO
  ENDDO
  ELSE
    PRINT *,'it non pris en compte'
    STOP
  ENDIF
  !
  ! NHL
  ! Auxiliary variables defined to deal with the fact that precipitation
  ! fluxes are defined on klev levels only.
  ! NHL
  !
  flxr_aux(:,klev+1)=0.0
  flxs_aux(:,klev+1)=0.0
  flxr_aux(:,1:klev)=flxr(:,:)
  flxs_aux(:,1:klev)=flxs(:,:)
  DO k=klev, 1, -1
  DO i=1, klon
  !--scavenging
    beta=flxr_aux(i,k)-flxr_aux(i,k+1)+flxs_aux(i,k)-flxs_aux(i,k+1)
    beta=beta/zdz(i,k)/qliq/zrho(i,k)
    beta=MAX(0.0,beta)
    dx=x(i,k)*(exp(-scav(i,k)*beta*pdtime)-1.)
    x(i,k)=x(i,k)+dx
    his_dh(i)=his_dh(i)-dx/RNAVO* &
          masse*1.e3*1.e6*zdz(i,k)/pdtime !--mgS/m2/s
  !--reevaporation
    beta=flxr_aux(i,k)-flxr_aux(i,k+1)+flxs_aux(i,k)-flxs_aux(i,k+1)
    IF (beta.LT.0.) beta=beta/(flxr_aux(i,k+1)+flxs_aux(i,k+1))
    IF (flxr_aux(i,k)+flxs_aux(i,k).EQ.0) THEN  !--reevaporation totale
      beta=MIN(MAX(0.0,-beta),1.0)
    ELSE                          !--reevaporation non totale pour aerosols
      ! !print *,'FRAC USED IN INSCAV_SPL'
      beta=MIN(MAX(0.0,-beta)*frac,1.0)
    ENDIF
    dx=beta*his_dh(i)*RNAVO/masse/1.e3/1.e6/zdz(i,k)*pdtime !ORIG LINE
  ! funny line for TL/AD
  ! AD test works without (x) and for xd = dxd*1.e5 : 2.79051851638 times the 0.
  ! AD test does not work with the line : 754592404.083 times the 0.
  ! problem seems to be linked to the largest dx wrt x
    ! x(i, k) = x(i, k) + dx
    !  x(i, k) = x(i, k) + dx         ! THIS LINE WAS COMMENTED OUT ORIGINALY !nhl
    his_dh(i)=(1.-beta)*his_dh(i)
  ENDDO
  ENDDO
  !
  RETURN
END SUBROUTINE inscav_spl