SUBROUTINE cltrac_spl(dtime, coef, yu1, yv1, t, tr, flux, paprs, pplay, d_tr)

  USE dimphy
  USE lmdz_yomcst
  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm

  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 second)
  ! coef-----input-R- le coefficient d'echange (m**2/s) l>1
  ! yu1------input-R- le vent dans le 1iere couche
  ! yv1------input-R- le vent dans le 1iere couche
  ! 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
  ! flux_tr--output-R- flux de tr
  !======================================================================
  REAL :: dtime
  REAL :: coef(klon, klev)
  REAL :: yu1(klon), yv1(klon)
  REAL :: t(klon, klev), tr(klon, klev)
  REAL :: paprs(klon, klev + 1), pplay(klon, klev), delp(klon, klev)
  REAL :: d_tr(klon, klev)
  REAL :: flux(klon), cdrag(klon), tr0(klon)
  ! REAL flux_tr(klon,klev)

  INTEGER :: i, k
  REAL :: zx_ctr(klon, 2:klev)
  REAL :: zx_dtr(klon, 2:klev)
  REAL :: zx_buf(klon)
  REAL :: zx_coef(klon, klev)
  REAL :: local_tr(klon, klev)
  REAL :: zx_alf1(klon), zx_alf2(klon), zx_flux(klon)
  !======================================================================
  ! CHECKING VALUES
  ! PRINT *,'CHECKING VALUES IN CLTRAC (INI)'
  ! PRINT *,'d_tr = ',sum(d_tr),MINVAL(d_tr),MAXVAL(d_tr)
  ! PRINT *,'flux = ',sum(flux),MINVAL(flux),MAXVAL(flux)
  ! PRINT *,'tr = ',sum(tr),MINVAL(tr),MAXVAL(tr)
  !======================================================================
  DO k = 1, klev
    DO i = 1, klon
      local_tr(i, k) = tr(i, k)
      delp(i, k) = paprs(i, k) - paprs(i, k + 1)
    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)
    zx_flux(i) = -flux(i) * dtime * RG
    !--pour le moment le flux est prescrit
    cdrag(i) = 0.0
    ! cdrag(i) =  coef(i,1) * (1.0+SQRT(yu1(i)**2+yv1(i)**2))
    ! .                * pplay(i,1)/(RD*t(i,1))
    tr0(i) = 0.0
    zx_coef(i, 1) = cdrag(i) * dtime * RG
  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)
  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
  !======================================================================
  ! PRINT *,'CHECKING VALUES IN CLTRAC (FIN)'
  ! PRINT *,'local_tr = ',sum(local_tr),MINVAL(local_tr),
  ! .                                    MAXVAL(local_tr)
  !  PRINT *,'zx_ctr = ',sum(zx_ctr),MINVAL(zx_ctr),MAXVAL(zx_ctr)
  !  PRINT *,'zx_dtr = ',sum(zx_dtr),MINVAL(zx_dtr),MAXVAL(zx_dtr)
  !  PRINT *,'tr = ',sum(tr),MINVAL(tr),MAXVAL(tr)
  !======================================================================
  !== 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
  ! PRINT *,'CHECKING VALUES IN CLTRAC (END)'
  ! PRINT *,'d_tr = ',sum(d_tr),MINVAL(d_tr),MAXVAL(d_tr)

END SUBROUTINE cltrac_spl