! $Id $

SUBROUTINE cltrac(dtime, coef, t, tr, flux, paprs, pplay, delp, &
        d_tr, d_tr_dry, flux_tr_dry)                    !jyg

  USE dimphy
  USE lmdz_yomcst

  IMPLICIT NONE
  !======================================================================
  ! Auteur(s): O. Boucher (LOA/LMD) date: 19961127
  !            inspire de clvent
  ! Objet: diffusion verticale de traceurs avec flux fixe a la surface
  !        ou/et flux du type c-drag

  ! Arguments:
  !-----------
  ! dtime.......input-R- intervalle du temps (en secondes)
  ! coef........input-R- le coefficient d'echange (m**2/s) l>1
  ! t...........input-R- temperature (K)
  ! tr..........input-R- la q. de traceurs
  ! flux........input-R- le flux de traceurs a la surface
  ! paprs.......input-R- pression a inter-couche (Pa)
  ! pplay.......input-R- pression au milieu de couche (Pa)
  ! delp........input-R- epaisseur de couche (Pa)
  ! cdrag.......input-R- cdrag pour le flux de surface (non active)
  ! tr0.........input-R- traceurs a la surface ou dans l'ocean (non active)
  ! d_tr........output-R- le changement de tr
  ! d_tr_dry....output-R- le changement de tr du au depot sec (1st layer)
  ! flux_tr_dry.output-R- depot sec
  !!! flux_tr..output-R- flux de tr
  !======================================================================

  ! Entree

  REAL, INTENT(IN) :: dtime
  REAL, DIMENSION(klon, klev), INTENT(IN) :: coef
  REAL, DIMENSION(klon, klev), INTENT(IN) :: t, tr
  REAL, DIMENSION(klon), INTENT(IN) :: flux !(at/s/m2)
  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
  REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay, delp

  ! Sorties

  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_tr
  REAL, DIMENSION(klon), INTENT(OUT) :: d_tr_dry          !jyg
  REAL, DIMENSION(klon), INTENT(OUT) :: flux_tr_dry       !jyg
  !  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: flux_tr

  ! Local

  INTEGER :: i, k
  REAL, DIMENSION(klon) :: cdrag, tr0
  REAL, DIMENSION(klon, klev) :: zx_ctr
  REAL, DIMENSION(klon, klev) :: zx_dtr
  REAL, DIMENSION(klon) :: zx_buf
  REAL, DIMENSION(klon, klev) :: zx_coef
  REAL, DIMENSION(klon, klev) :: local_tr
  REAL, DIMENSION(klon) :: zx_alf1, zx_alf2, zx_flux

  !======================================================================

  DO k = 1, klev
    DO i = 1, klon
      local_tr(i, k) = tr(i, k)
    ENDDO
  ENDDO

  !======================================================================

  DO i = 1, klon
    zx_alf1(i) = (paprs(i, 1) - pplay(i, 2)) / (pplay(i, 1) - pplay(i, 2))
    zx_alf2(i) = 1.0 - zx_alf1(i)
    flux_tr_dry(i) = -flux(i) * dtime                              !jyg
    zx_flux(i) = flux_tr_dry(i) * RG                              !jyg
    !!     zx_flux(i) =  -flux(i)*dtime*RG                            !jyg
    ! Pour le moment le flux est prescrit cdrag et zx_coef(1) vaut 0
    cdrag(i) = 0.0
    tr0(i) = 0.0
    zx_coef(i, 1) = cdrag(i) * dtime * RG
    zx_ctr(i, 1) = 0.
    zx_dtr(i, 1) = 0.
  ENDDO

  !======================================================================

  DO k = 2, klev
    DO i = 1, klon
      zx_coef(i, k) = coef(i, k) * RG / (pplay(i, k - 1) - pplay(i, k))   &
              * (paprs(i, k) * 2 / (t(i, k) + t(i, k - 1)) / RD)**2
      zx_coef(i, k) = zx_coef(i, k) * dtime * RG
    ENDDO
  ENDDO

  !======================================================================

  DO i = 1, klon
    zx_buf(i) = delp(i, 1) + zx_coef(i, 1) * zx_alf1(i) + zx_coef(i, 2)

    zx_ctr(i, 2) = (local_tr(i, 1) * delp(i, 1) + &
            zx_coef(i, 1) * tr0(i) - zx_flux(i)) / zx_buf(i)

    zx_dtr(i, 2) = (zx_coef(i, 2) - zx_alf2(i) * zx_coef(i, 1)) / &
            zx_buf(i)
    d_tr_dry(i) = -zx_flux(i) / zx_buf(i)                          !jyg
  ENDDO

  DO k = 3, klev
    DO i = 1, klon
      zx_buf(i) = delp(i, k - 1) + zx_coef(i, k)      &
              + zx_coef(i, k - 1) * (1. - zx_dtr(i, k - 1))
      zx_ctr(i, k) = (local_tr(i, k - 1) * delp(i, k - 1)  &
              + zx_coef(i, k - 1) * zx_ctr(i, k - 1)) / zx_buf(i)
      zx_dtr(i, k) = zx_coef(i, k) / zx_buf(i)
    ENDDO
  ENDDO

  DO i = 1, klon
    local_tr(i, klev) = (local_tr(i, klev) * delp(i, klev) &
            + zx_coef(i, klev) * zx_ctr(i, klev))             &
            / (delp(i, klev) + zx_coef(i, klev)            &
                    - zx_coef(i, klev) * zx_dtr(i, klev))
  ENDDO

  DO k = klev - 1, 1, -1
    DO i = 1, klon
      local_tr(i, k) = zx_ctr(i, k + 1) + zx_dtr(i, k + 1) * local_tr(i, k + 1)
    ENDDO
  ENDDO

  !======================================================================
  !== flux_tr est le flux de traceur (positif vers bas)
  !      DO i = 1, klon
  !         flux_tr(i,1) = zx_coef(i,1)/(RG*dtime)
  !      ENDDO
  !      DO k = 2, klev
  !      DO i = 1, klon
  !         flux_tr(i,k) = zx_coef(i,k)/(RG*dtime)
  !     .               * (local_tr(i,k)-local_tr(i,k-1))
  !      ENDDO
  !      ENDDO
  !======================================================================
  DO k = 1, klev
    DO i = 1, klon
      d_tr(i, k) = local_tr(i, k) - tr(i, k)
    ENDDO
  ENDDO

END SUBROUTINE cltrac