! $Header$

SUBROUTINE yamada4(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
    cd, q2, km, kn, kq, ustar, iflag_pbl)
  USE dimphy
  IMPLICIT NONE
  include "iniprint.h"
  ! .......................................................................
  ! ym#include "dimensions.h"
  ! ym#include "dimphy.h"
  ! .......................................................................

  ! dt : pas de temps
  ! g  : g
  ! zlev : altitude a chaque niveau (interface inferieure de la couche
  ! de meme indice)
  ! zlay : altitude au centre de chaque couche
  ! u,v : vitesse au centre de chaque couche
  ! (en entree : la valeur au debut du pas de temps)
  ! teta : temperature potentielle au centre de chaque couche
  ! (en entree : la valeur au debut du pas de temps)
  ! cd : cdrag
  ! (en entree : la valeur au debut du pas de temps)
  ! q2 : $q^2$ au bas de chaque couche
  ! (en entree : la valeur au debut du pas de temps)
  ! (en sortie : la valeur a la fin du pas de temps)
  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
  ! couche)
  ! (en sortie : la valeur a la fin du pas de temps)
  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
  ! (en sortie : la valeur a la fin du pas de temps)

  ! iflag_pbl doit valoir entre 6 et 9
  ! l=6, on prend  systematiquement une longueur d'equilibre
  ! iflag_pbl=6 : MY 2.0
  ! iflag_pbl=7 : MY 2.0.Fournier
  ! iflag_pbl=8/9 : MY 2.5
  ! iflag_pbl=8 with special obsolete treatments for convergence
  ! with Cmpi5 NPv3.1 simulations
  ! iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
  ! iflag_pbl=12 = 11 with vertical diffusion off q2

  ! 2013/04/01 (FH hourdin@lmd.jussieu.fr)
  ! Correction for very stable PBLs (iflag_pbl=10 and 11)
  ! iflag_pbl=8 converges numerically with NPv3.1
  ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l
  ! -> the model can run with longer time-steps.
  ! .......................................................................

  REAL dt, g, rconst
  REAL plev(klon, klev+1), temp(klon, klev)
  REAL ustar(klon)
  REAL kmin, qmin, pblhmin(klon), coriol(klon)
  REAL zlev(klon, klev+1)
  REAL zlay(klon, klev)
  REAL u(klon, klev)
  REAL v(klon, klev)
  REAL teta(klon, klev)
  REAL cd(klon)
  REAL q2(klon, klev+1), qpre
  REAL unsdz(klon, klev)
  REAL unsdzdec(klon, klev+1)

  REAL km(klon, klev+1)
  REAL kmpre(klon, klev+1), tmp2
  REAL mpre(klon, klev+1)
  REAL kn(klon, klev+1)
  REAL kq(klon, klev+1)
  REAL ff(klon, klev+1), delta(klon, klev+1)
  REAL aa(klon, klev+1), aa0, aa1
  INTEGER iflag_pbl, ngrid
  INTEGER nlay, nlev

  LOGICAL first
  INTEGER ipas
  SAVE first, ipas
  ! FH/IM     data first,ipas/.true.,0/
  DATA first, ipas/.FALSE., 0/
  !$OMP THREADPRIVATE( first,ipas)

  INTEGER ig, k


  REAL ri, zrif, zalpha, zsm, zsn
  REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev)

  REAL m2(klon, klev+1), dz(klon, klev+1), zq, n2(klon, klev+1)
  REAL dtetadz(klon, klev+1)
  REAL m2cstat, mcstat, kmcstat
  REAL l(klon, klev+1)
  REAL, ALLOCATABLE, SAVE :: l0(:)
  !$OMP THREADPRIVATE(l0)
  REAL sq(klon), sqz(klon), zz(klon, klev+1)
  INTEGER iter

  REAL ric, rifc, b1, kap
  SAVE ric, rifc, b1, kap
  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
  !$OMP THREADPRIVATE(ric,rifc,b1,kap)
  REAL frif, falpha, fsm
  REAL fl, zzz, zl0, zq2, zn2

  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)
  LOGICAL, SAVE :: firstcall = .TRUE.
  !$OMP THREADPRIVATE(firstcall)
  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)
  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
  fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &
    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &
    max(n2(ig,k),1.E-10))), 1.)


  nlay = klev
  nlev = klev + 1

  IF (firstcall) THEN
    ALLOCATE (l0(klon))
    firstcall = .FALSE.
  END IF


  IF (.NOT. (iflag_pbl>=6 .AND. iflag_pbl<=12)) THEN
    STOP 'probleme de coherence dans appel a MY'
  END IF

  ipas = ipas + 1


  ! .......................................................................
  ! les increments verticaux
  ! .......................................................................

  ! !!!!! allerte !!!!!c
  ! !!!!! zlev n'est pas declare a nlev !!!!!c
  ! !!!!! ---->
  DO ig = 1, ngrid
    zlev(ig, nlev) = zlay(ig, nlay) + (zlay(ig,nlay)-zlev(ig,nlev-1))
  END DO
  ! !!!!! <----
  ! !!!!! allerte !!!!!c

  DO k = 1, nlay
    DO ig = 1, ngrid
      unsdz(ig, k) = 1.E+0/(zlev(ig,k+1)-zlev(ig,k))
    END DO
  END DO
  DO ig = 1, ngrid
    unsdzdec(ig, 1) = 1.E+0/(zlay(ig,1)-zlev(ig,1))
  END DO
  DO k = 2, nlay
    DO ig = 1, ngrid
      unsdzdec(ig, k) = 1.E+0/(zlay(ig,k)-zlay(ig,k-1))
    END DO
  END DO
  DO ig = 1, ngrid
    unsdzdec(ig, nlay+1) = 1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
  END DO

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Computing M^2, N^2, Richardson numbers, stability functions
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! initialize arrays:
  m2(:, :) = 0.0
  sm(:, :) = 0.0
  rif(:, :) = 0.0

  DO k = 2, klev
    DO ig = 1, ngrid
      dz(ig, k) = zlay(ig, k) - zlay(ig, k-1)
      m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &
        k-1))**2)/(dz(ig,k)*dz(ig,k))
      dtetadz(ig, k) = (teta(ig,k)-teta(ig,k-1))/dz(ig, k)
      n2(ig, k) = g*2.*dtetadz(ig, k)/(teta(ig,k-1)+teta(ig,k))
      ! n2(ig,k)=0.
      ri = n2(ig, k)/max(m2(ig,k), 1.E-10)
      IF (ri<ric) THEN
        rif(ig, k) = frif(ri)
      ELSE
        rif(ig, k) = rifc
      END IF
      IF (rif(ig,k)<0.16) THEN
        alpha(ig, k) = falpha(rif(ig,k))
        sm(ig, k) = fsm(rif(ig,k))
      ELSE
        alpha(ig, k) = 1.12
        sm(ig, k) = 0.085
      END IF
      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
    END DO
  END DO


  ! ====================================================================
  ! Computing the mixing length
  ! ====================================================================

  ! Mise a jour de l0
  IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Iterative computation of l0
    ! This version is kept for iflag_pbl only for convergence
    ! with NPv3.1 Cmip5 simulations
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    DO ig = 1, ngrid
      sq(ig) = 1.E-10
      sqz(ig) = 1.E-10
    END DO
    DO k = 2, klev - 1
      DO ig = 1, ngrid
        zq = sqrt(q2(ig,k))
        sqz(ig) = sqz(ig) + zq*zlev(ig, k)*(zlay(ig,k)-zlay(ig,k-1))
        sq(ig) = sq(ig) + zq*(zlay(ig,k)-zlay(ig,k-1))
      END DO
    END DO
    DO ig = 1, ngrid
      l0(ig) = 0.2*sqz(ig)/sq(ig)
    END DO
    DO k = 2, klev
      DO ig = 1, ngrid
        l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
      END DO
    END DO
    ! print*,'L0 cas 8 ou 10 ',l0

  ELSE

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! In all other case, the assymptotic mixing length l0 is imposed (100m)
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    l0(:) = 150.
    DO k = 2, klev
      DO ig = 1, ngrid
        l(ig, k) = fl(zlev(ig,k), l0(ig), q2(ig,k), n2(ig,k))
      END DO
    END DO
    ! print*,'L0 cas autres ',l0

  END IF


  ! ====================================================================
  ! Yamada 2.0
  ! ====================================================================
  IF (iflag_pbl==6) THEN

    DO k = 2, klev
      q2(:, k) = l(:, k)**2*zz(:, k)
    END DO


  ELSE IF (iflag_pbl==7) THEN
    ! ====================================================================
    ! Yamada 2.Fournier
    ! ====================================================================

    ! Calcul de l,  km, au pas precedent
    DO k = 2, klev
      DO ig = 1, ngrid
        ! print*,'SMML=',sm(ig,k),l(ig,k)
        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
        kmpre(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
        mpre(ig, k) = sqrt(m2(ig,k))
        ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
      END DO
    END DO

    DO k = 2, klev - 1
      DO ig = 1, ngrid
        m2cstat = max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1, 1.E-12)
        mcstat = sqrt(m2cstat)

        ! print*,'M2 L=',k,mpre(ig,k),mcstat

        ! -----{puis on ecrit la valeur de q qui annule l'equation de m
        ! supposee en q3}

        IF (k==2) THEN
          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
            unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec &
            (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &
            -1))
        ELSE
          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &
            unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ &
            (unsdz(ig,k)+unsdz(ig,k-1))
        END IF
        ! print*,'T2 L=',k,tmp2
        tmp2 = kmcstat/(sm(ig,k)/q2(ig,k))/l(ig, k)
        q2(ig, k) = max(tmp2, 1.E-12)**(2./3.)
        ! print*,'Q2 L=',k,q2(ig,k)

      END DO
    END DO

  ELSE IF (iflag_pbl==8 .OR. iflag_pbl==9) THEN
    ! ====================================================================
    ! Yamada 2.5 a la Didi
    ! ====================================================================


    ! Calcul de l,  km, au pas precedent
    DO k = 2, klev
      DO ig = 1, ngrid
        ! print*,'SMML=',sm(ig,k),l(ig,k)
        delta(ig, k) = q2(ig, k)/(l(ig,k)**2*sm(ig,k))
        IF (delta(ig,k)<1.E-20) THEN
          ! print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
          delta(ig, k) = 1.E-20
        END IF
        km(ig, k) = l(ig, k)*sqrt(q2(ig,k))*sm(ig, k)
        aa0 = (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
        aa1 = (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
        ! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
        aa(ig, k) = aa1*dt/(delta(ig,k)*l(ig,k))
        ! print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
        qpre = sqrt(q2(ig,k))
        ! if (iflag_pbl.eq.8 ) then
        IF (aa(ig,k)>0.) THEN
          q2(ig, k) = (qpre+aa(ig,k)*qpre*qpre)**2
        ELSE
          q2(ig, k) = (qpre/(1.-aa(ig,k)*qpre))**2
        END IF
        ! else ! iflag_pbl=9
        ! if (aa(ig,k)*qpre.gt.0.9) then
        ! q2(ig,k)=(qpre*10.)**2
        ! else
        ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
        ! endif
        ! endif
        q2(ig, k) = min(max(q2(ig,k),1.E-10), 1.E4)
        ! print*,'Q2 L=',k,q2(ig,k),qpre*qpre
      END DO
    END DO

  ELSE IF (iflag_pbl>=10) THEN

    ! print*,'Schema mixte D'
    ! print*,'Longueur ',l(:,:)
    DO k = 2, klev - 1
      l(:, k) = max(l(:,k), 1.)
      km(:, k) = l(:, k)*sqrt(q2(:,k))*sm(:, k)
      q2(:, k) = q2(:, k) + dt*km(:, k)*m2(:, k)*(1.-rif(:,k))
      q2(:, k) = min(max(q2(:,k),1.E-10), 1.E4)
      q2(:, k) = 1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))
      q2(:, k) = q2(:, k)*q2(:, k)
    END DO


  ELSE
    STOP 'Cas nom prevu dans yamada4'

  END IF ! Fin du cas 8

  ! print*,'OK8'

  ! ====================================================================
  ! Calcul des coefficients de m�ange
  ! ====================================================================
  DO k = 2, klev
    ! print*,'k=',k
    DO ig = 1, ngrid
      ! abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
      zq = sqrt(q2(ig,k))
      km(ig, k) = l(ig, k)*zq*sm(ig, k)
      kn(ig, k) = km(ig, k)*alpha(ig, k)
      kq(ig, k) = l(ig, k)*zq*0.2
      ! print*,'KML=',km(ig,k),kn(ig,k)
    END DO
  END DO
    ! initialize near-surface and top-layer mixing coefficients
  kq(1:ngrid, 1) = kq(1:ngrid, 2) ! constant (ie no gradient) near the surface
  kq(1:ngrid, klev+1) = 0 ! zero at the top

  ! Transport diffusif vertical de la TKE.
  IF (iflag_pbl>=12) THEN
    ! print*,'YAMADA VDIF'
    q2(:, 1) = q2(:, 2)
    CALL vdif_q2(dt, g, rconst, ngrid, plev, temp, kq, q2)
  END IF

  ! Traitement des cas noctrunes avec l'introduction d'une longueur
  ! minilale.

  ! ====================================================================
  ! Traitement particulier pour les cas tres stables.
  ! D'apres Holtslag Boville.

  IF (prt_level>1) THEN
    PRINT *, 'YAMADA4 0'
  END IF !(prt_level>1) THEN
  DO ig = 1, ngrid
    coriol(ig) = 1.E-4
    pblhmin(ig) = 0.07*ustar(ig)/max(abs(coriol(ig)), 2.546E-5)
  END DO

  ! print*,'pblhmin ',pblhmin
  ! Test a remettre 21 11 02
  ! test abd 13 05 02      if(0.eq.1) then
  IF (1==1) THEN
    IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN

      DO k = 2, klev
        DO ig = 1, ngrid
          IF (teta(ig,2)>teta(ig,1)) THEN
            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
            kmin = kap*zlev(ig, k)*qmin
          ELSE
            kmin = -1. ! kmin n'est utilise que pour les SL stables.
          END IF
          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
            kn(ig, k) = kmin
            km(ig, k) = kmin
            kq(ig, k) = kmin
            ! la longueur de melange est suposee etre l= kap z
            ! K=l q Sm d'ou q2=(K/l Sm)**2
            q2(ig, k) = (qmin/sm(ig,k))**2
          END IF
        END DO
      END DO

    ELSE

      DO k = 2, klev
        DO ig = 1, ngrid
          IF (teta(ig,2)>teta(ig,1)) THEN
            qmin = ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
            kmin = kap*zlev(ig, k)*qmin
          ELSE
            kmin = -1. ! kmin n'est utilise que pour les SL stables.
          END IF
          IF (kn(ig,k)<kmin .OR. km(ig,k)<kmin) THEN
            ! print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
            ! s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
            kn(ig, k) = kmin
            km(ig, k) = kmin
            kq(ig, k) = kmin
            ! la longueur de melange est suposee etre l= kap z
            ! K=l q Sm d'ou q2=(K/l Sm)**2
            sm(ig, k) = 1.
            alpha(ig, k) = 1.
            q2(ig, k) = min((qmin/sm(ig,k))**2, 10.)
            zq = sqrt(q2(ig,k))
            km(ig, k) = l(ig, k)*zq*sm(ig, k)
            kn(ig, k) = km(ig, k)*alpha(ig, k)
            kq(ig, k) = l(ig, k)*zq*0.2
          END IF
        END DO
      END DO
    END IF

  END IF

  IF (prt_level>1) THEN
    PRINT *, 'YAMADA4 1'
  END IF !(prt_level>1) THEN
  ! Diagnostique pour stokage

  IF (1==0) THEN
    rino = rif
    smyam(1:ngrid, 1) = 0.
    styam(1:ngrid, 1) = 0.
    lyam(1:ngrid, 1) = 0.
    knyam(1:ngrid, 1) = 0.
    w2yam(1:ngrid, 1) = 0.
    t2yam(1:ngrid, 1) = 0.

    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)
    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)
    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)
    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)

    ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane

    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &
      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
      sqrt(q2(1:ngrid,2:klev))

    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
      lyam(1:ngrid, 2:klev)
  END IF

  ! print*,'OKFIN'
  first = .FALSE.
  RETURN
END SUBROUTINE yamada4
SUBROUTINE vdif_q2(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
  USE dimphy
  IMPLICIT NONE
  ! .......................................................................
  include "dimensions.h"
  ! ccc#include "dimphy.h"
  ! .......................................................................

  ! dt : pas de temps

  REAL plev(klon, klev+1)
  REAL temp(klon, klev)
  REAL timestep
  REAL gravity, rconst
  REAL kstar(klon, klev+1), zz
  REAL kmy(klon, klev+1)
  REAL q2(klon, klev+1)
  REAL deltap(klon, klev+1)
  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
  INTEGER ngrid

  INTEGER i, k

  ! print*,'RD=',rconst
  DO k = 1, klev
    DO i = 1, ngrid
      ! test
      ! print*,'i,k',i,k
      ! print*,'temp(i,k)=',temp(i,k)
      ! print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
        (plev(i,k)-plev(i,k+1))*timestep
    END DO
  END DO

  DO k = 2, klev
    DO i = 1, ngrid
      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
    END DO
  END DO
  DO i = 1, ngrid
    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
    denom(i, klev+1) = deltap(i, klev+1) + kstar(i, klev)
    alpha(i, klev+1) = deltap(i, klev+1)*q2(i, klev+1)/denom(i, klev+1)
    beta(i, klev+1) = kstar(i, klev)/denom(i, klev+1)
  END DO

  DO k = klev, 2, -1
    DO i = 1, ngrid
      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &
        kstar(i, k-1)
      ! correction d'un bug 10 01 2001
      alpha(i, k) = (q2(i,k)*deltap(i,k)+kstar(i,k)*alpha(i,k+1))/denom(i, k)
      beta(i, k) = kstar(i, k-1)/denom(i, k)
    END DO
  END DO

  ! Si on recalcule q2(1)
  IF (1==0) THEN
    DO i = 1, ngrid
      denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)
      q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)
    END DO
  END IF
  ! sinon, on peut sauter cette boucle...

  DO k = 2, klev + 1
    DO i = 1, ngrid
      q2(i, k) = alpha(i, k) + beta(i, k)*q2(i, k-1)
    END DO
  END DO

  RETURN
END SUBROUTINE vdif_q2
SUBROUTINE vdif_q2e(timestep, gravity, rconst, ngrid, plev, temp, kmy, q2)
  USE dimphy
  IMPLICIT NONE
  ! .......................................................................
  include "dimensions.h"
  ! ccc#include "dimphy.h"
  ! .......................................................................

  ! dt : pas de temps

  REAL plev(klon, klev+1)
  REAL temp(klon, klev)
  REAL timestep
  REAL gravity, rconst
  REAL kstar(klon, klev+1), zz
  REAL kmy(klon, klev+1)
  REAL q2(klon, klev+1)
  REAL deltap(klon, klev+1)
  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)
  INTEGER ngrid

  INTEGER i, k

  DO k = 1, klev
    DO i = 1, ngrid
      zz = (plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &
        (plev(i,k)-plev(i,k+1))*timestep
    END DO
  END DO

  DO k = 2, klev
    DO i = 1, ngrid
      deltap(i, k) = 0.5*(plev(i,k-1)-plev(i,k+1))
    END DO
  END DO
  DO i = 1, ngrid
    deltap(i, 1) = 0.5*(plev(i,1)-plev(i,2))
    deltap(i, klev+1) = 0.5*(plev(i,klev)-plev(i,klev+1))
  END DO

  DO k = klev, 2, -1
    DO i = 1, ngrid
      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &
        k))-kstar(i,k-1)*(q2(i,k)-q2(i,k-1)))/deltap(i, k)
    END DO
  END DO

  DO i = 1, ngrid
    q2(i, 1) = q2(i, 1) + (kstar(i,1)*(q2(i,2)-q2(i,1)))/deltap(i, 1)
    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &
      klev)))/deltap(i, klev+1)
  END DO

  RETURN
END SUBROUTINE vdif_q2e