!
! $Header$
!
      SUBROUTINE vdif_kcay(ngrid,dt,g,rconst,plev,temp
     s   ,zlev,zlay,u,v,teta,cd,q2,q2diag,km,kn,ustar
     s   ,l_mix)
      use dimphy
      IMPLICIT NONE
c.......................................................................
cym#include "dimensions.h"
cym#include "dimphy.h"
c.......................................................................
c
c dt : pas de temps
c g  : g
c zlev : altitude a chaque niveau (interface inferieure de la couche
c        de meme indice)
c zlay : altitude au centre de chaque couche
c u,v : vitesse au centre de chaque couche
c       (en entree : la valeur au debut du pas de temps)
c teta : temperature potentielle au centre de chaque couche
c        (en entree : la valeur au debut du pas de temps)
c cd : cdrag
c      (en entree : la valeur au debut du pas de temps)
c q2 : $q^2$ au bas de chaque couche
c      (en entree : la valeur au debut du pas de temps)
c      (en sortie : la valeur a la fin du pas de temps)
c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
c      couche)
c      (en sortie : la valeur a la fin du pas de temps)
c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
c      (en sortie : la valeur a la fin du pas de temps)
c 
c.......................................................................
      REAL dt,g,rconst
      real plev(klon,klev+1),temp(klon,klev)
      real ustar(klon),snstable
      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),q2s(klon,klev+1)
      REAL q2diag(klon,klev+1)
      REAL km(klon,klev+1)
      REAL kn(klon,klev+1)
      real sq(klon),sqz(klon),zz(klon,klev+1),zq,long0(klon)

      integer l_mix,iii
c.......................................................................
c
c nlay : nombre de couches        
c nlev : nombre de niveaux
c ngrid : nombre de points de grille       
c unsdz : 1 sur l'epaisseur de couche
c unsdzdec : 1 sur la distance entre le centre de la couche et le
c            centre de la couche inferieure
c q : echelle de vitesse au bas de chaque couche
c     (valeur a la fin du pas de temps)
c
c.......................................................................
      INTEGER nlay,nlev,ngrid
      REAL unsdz(klon,klev)
      REAL unsdzdec(klon,klev+1)
      REAL q(klon,klev+1)

c.......................................................................
c
c kmpre : km au debut du pas de temps
c qcstat : q : solution stationnaire du probleme couple
c          (valeur a la fin du pas de temps)
c q2cstat : q2 : solution stationnaire du probleme couple
c           (valeur a la fin du pas de temps)
c
c.......................................................................
      REAL kmpre(klon,klev+1)
      REAL qcstat
      REAL q2cstat
      real sss,sssq
c.......................................................................
c
c long : longueur de melange calculee selon Blackadar
c
c.......................................................................
      REAL long(klon,klev+1)
c.......................................................................
c
c kmq3 : terme en q^3 dans le developpement de km
c        (valeur au debut du pas de temps)
c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
c           (valeur a la fin du pas de temps)
c knq3 : terme en q^3 dans le developpement de kn
c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
c          (valeur a la fin du pas de temps)
c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
c           (valeur a la fin du pas de temps)
c m : valeur a la fin du pas de temps
c mpre : valeur au debut du pas de temps
c m2 : valeur a la fin du pas de temps
c n2 : valeur a la fin du pas de temps
c 
c.......................................................................
      REAL kmq3
      REAL kmcstat
      REAL knq3
      REAL mcstat
      REAL m2cstat
      REAL m(klon,klev+1)
      REAL mpre(klon,klev+1)
      REAL m2(klon,klev+1)
      REAL n2(klon,klev+1)
c.......................................................................
c
c gn : intermediaire pour les coefficients de stabilite
c gnmin : borne inferieure de gn (-0.23 ou -0.28)
c gnmax : borne superieure de gn (0.0233)
c gninf : vrai si gn est en dessous de sa borne inferieure
c gnsup : vrai si gn est en dessus de sa borne superieure
c gm : drole d'objet bien utile
c ri : nombre de Richardson
c sn : coefficient de stabilite pour n
c snq2 : premier terme du developement limite de sn en q2
c sm : coefficient de stabilite pour m
c smq2 : premier terme du developement limite de sm en q2
c
c.......................................................................
      REAL gn
      REAL gnmin
      REAL gnmax
      LOGICAL gninf
      LOGICAL gnsup
      REAL gm
c      REAL ri(klon,klev+1)
      REAL sn(klon,klev+1)
      REAL snq2(klon,klev+1)
      REAL sm(klon,klev+1)
      REAL smq2(klon,klev+1)
c.......................................................................
c
c kappa : consatnte de Von Karman (0.4)
c long00 : longueur de reference pour le calcul de long (160)
c a1,a2,b1,b2,c1 : constantes d'origine pour les  coefficients
c                  de stabilite (0.92/0.74/16.6/10.1/0.08)
c cn1,cn2 : constantes pour sn
c cm1,cm2,cm3,cm4 : constantes pour sm
c
c.......................................................................
      REAL kappa
      REAL long00
      REAL a1,a2,b1,b2,c1
      REAL cn1,cn2
      REAL cm1,cm2,cm3,cm4
c.......................................................................
c
c termq : termes en $q$ dans l'equation de q2
c termq3 : termes en $q^3$ dans l'equation de q2
c termqm2 : termes en $q*m^2$ dans l'equation de q2
c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
c
c.......................................................................
      REAL termq
      REAL termq3
      REAL termqm2
      REAL termq3m2
c.......................................................................
c
c q2min : borne inferieure de q2
c q2max : borne superieure de q2
c
c.......................................................................
      REAL q2min
      REAL q2max
c.......................................................................
c knmin : borne inferieure de kn
c kmmin : borne inferieure de km
c.......................................................................
      REAL knmin
      REAL kmmin
c.......................................................................
      INTEGER ilay,ilev,igrid
      REAL tmp1,tmp2
c.......................................................................
      PARAMETER (kappa=0.4E+0)
      PARAMETER (long00=160.E+0)
c     PARAMETER (gnmin=-10.E+0)
      PARAMETER (gnmin=-0.28)
      PARAMETER (gnmax=0.0233E+0)
      PARAMETER (a1=0.92E+0)
      PARAMETER (a2=0.74E+0)
      PARAMETER (b1=16.6E+0)
      PARAMETER (b2=10.1E+0)
      PARAMETER (c1=0.08E+0)
      PARAMETER (knmin=1.E-5)
      PARAMETER (kmmin=1.E-5)
      PARAMETER (q2min=1.e-5)
      PARAMETER (q2max=1.E+2)
cym      PARAMETER (nlay=klev)
cym      PARAMETER (nlev=klev+1)
c
      PARAMETER (
     &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
     &          )
      PARAMETER (
     &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
     &          )
      PARAMETER (
     &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
     &          )
      PARAMETER (
     &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
     &          -3.E+0 *c1*(b2+6.E+0 *a1)))
     &          )
      PARAMETER (
     &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
     &          )
      PARAMETER (
     &  cm4=-9.E+0 *a1*a2
     &          )

      logical first
      save first
      data first/.true./
c$OMP THREADPRIVATE(first)
c.......................................................................
c  traitment des valeur de q2 en entree
c.......................................................................
c
c   Initialisation de q2
      nlay=klev
      nlev=klev+1
       
      call yamada(ngrid,dt,g,rconst,plev,temp
     s   ,zlev,zlay,u,v,teta,cd,q2diag,km,kn,ustar
     s   ,l_mix)
      if (first.and.1.eq.1) then
      first=.false.
      q2=q2diag
      endif

      DO ilev=1,nlev
                                                      DO igrid=1,ngrid 
        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
        q(igrid,ilev)=sqrt(q2(igrid,ilev))
                                                      ENDDO
      ENDDO
c
                                                      DO igrid=1,ngrid 
      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
      q2(igrid,1)=amax1(q2(igrid,1),q2min)
      q(igrid,1)=sqrt(q2(igrid,1))
                                                      ENDDO
c
c.......................................................................
c  les increments verticaux
c.......................................................................
c
c!!!!! allerte !!!!!c
c!!!!! zlev n'est pas declare a nlev !!!!!c
c!!!!! ---->
                                                      DO igrid=1,ngrid 
            zlev(igrid,nlev)=zlay(igrid,nlay)
     &             +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
                                                      ENDDO            
c!!!!! <----
c!!!!! allerte !!!!!c
c
      DO ilay=1,nlay
                                                      DO igrid=1,ngrid 
        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
                                                      ENDDO
      ENDDO
                                                      DO igrid=1,ngrid 
      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
                                                      ENDDO
      DO ilay=2,nlay
                                                      DO igrid=1,ngrid 
        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
                                                      ENDDO
      ENDDO
                                                      DO igrid=1,ngrid 
      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
                                                      ENDDO
c
c.......................................................................
c  le cisaillement et le gradient de temperature
c.......................................................................
c
                                                      DO igrid=1,ngrid 
      m2(igrid,1)=(unsdzdec(igrid,1)
     &                   *u(igrid,1))**2
     &                 +(unsdzdec(igrid,1)
     &                   *v(igrid,1))**2
      m(igrid,1)=sqrt(m2(igrid,1))
      mpre(igrid,1)=m(igrid,1)
                                                      ENDDO
c
c-----------------------------------------------------------------------
      DO ilev=2,nlev-1
                                                      DO igrid=1,ngrid 
c-----------------------------------------------------------------------
c
        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
     &                   *(teta(igrid,ilev)-teta(igrid,ilev-1))
     &                   /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
c       n2(igrid,ilev)=0.
c
c --->
c       on ne sais traiter que les cas stratifies. et l'ajustement
c       convectif est cense faire en sorte que seul des configurations
c       stratifiees soient rencontrees en entree de cette routine.
c       mais, bon ... on sait jamais (meme on sait que n2 prends
c       quelques valeurs negatives ... parfois) alors : 
c<---
c
        IF (n2(igrid,ilev).lt.0.E+0) THEN
          n2(igrid,ilev)=0.E+0
        ENDIF
c
        m2(igrid,ilev)=(unsdzdec(igrid,ilev)
     &                     *(u(igrid,ilev)-u(igrid,ilev-1)))**2
     &                   +(unsdzdec(igrid,ilev)
     &                     *(v(igrid,ilev)-v(igrid,ilev-1)))**2
        m(igrid,ilev)=sqrt(m2(igrid,ilev))
        mpre(igrid,ilev)=m(igrid,ilev)
c
c-----------------------------------------------------------------------
                                                      ENDDO
      ENDDO
c-----------------------------------------------------------------------
c
                                                      DO igrid=1,ngrid 
      m2(igrid,nlev)=m2(igrid,nlev-1)
      m(igrid,nlev)=m(igrid,nlev-1)
      mpre(igrid,nlev)=m(igrid,nlev)
                                                      ENDDO
c
c.......................................................................
c  calcul des fonctions de stabilite
c.......................................................................
c
      if (l_mix.eq.4) then
                                                      DO igrid=1,ngrid 
         sqz(igrid)=1.e-10
         sq(igrid)=1.e-10
                                                      ENDDO
         do ilev=2,nlev-1
                                                      DO igrid=1,ngrid 
           zq=sqrt(q2(igrid,ilev))
           sqz(igrid)
     .     =sqz(igrid)+zq*zlev(igrid,ilev)
     .     *(zlay(igrid,ilev)-zlay(igrid,ilev-1))
           sq(igrid)=sq(igrid)+zq*(zlay(igrid,ilev)-zlay(igrid,ilev-1))
                                                      ENDDO
         enddo
                                                      DO igrid=1,ngrid 
         long0(igrid)=0.2*sqz(igrid)/sq(igrid)
                                                      ENDDO
      else if (l_mix.eq.3) then
         long0(igrid)=long00
      endif

c (abd 5 2)      print*,'LONG0=',long0

c-----------------------------------------------------------------------
      DO ilev=2,nlev-1
                                                      DO igrid=1,ngrid 
c-----------------------------------------------------------------------
c
        tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
        if (l_mix.ge.10) then
            long(igrid,ilev)=l_mix
        else
           long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0(igrid))
        endif
        long(igrid,ilev)=max(min(long(igrid,ilev)
     s    ,0.5*sqrt(q2(igrid,ilev))/sqrt(max(n2(igrid,ilev),1.e-10)))
     s    ,5.)

        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
     &                                           * n2(igrid,ilev)
        gm=long(igrid,ilev)**2 / q2(igrid,ilev)
     &                                           * m2(igrid,ilev)
c
        gninf=.false.
        gnsup=.false.
        long(igrid,ilev)=long(igrid,ilev)
        long(igrid,ilev)=long(igrid,ilev)
c
        IF (gn.lt.gnmin) THEN
          gninf=.true.
          gn=gnmin
        ENDIF
c
        IF (gn.gt.gnmax) THEN
          gnsup=.true.
          gn=gnmax
        ENDIF
c
        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
        sm(igrid,ilev)=
     &    (cm1+cm2*gn)
     &   /( (1.E+0 +cm3*gn)
     &     *(1.E+0 +cm4*gn) )
c
        IF ((gninf).or.(gnsup)) THEN
          snq2(igrid,ilev)=0.E+0
          smq2(igrid,ilev)=0.E+0
        ELSE
          snq2(igrid,ilev)=
     &     -gn
     &     *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
          smq2(igrid,ilev)=
     &     -gn
     &     *( cm2*(1.E+0 +cm3*gn)
     &           *(1.E+0 +cm4*gn)
     &       -( cm3*(1.E+0 +cm4*gn)
     &         +cm4*(1.E+0 +cm3*gn) )
     &       *(cm1+cm2*gn)            )
     &     /( (1.E+0 +cm3*gn)
     &       *(1.E+0 +cm4*gn) )**2
        ENDIF
c
c abd
c        if(ilev.le.57.and.ilev.ge.37) then
c            print*,'L=',ilev,'   GN=',gn,'  SM=',sm(igrid,ilev)
c        endif
c --->
c       la decomposition de Taylor en q2 n'a de sens que
c       dans les cas stratifies ou sn et sm sont quasi
c       proportionnels a q2. ailleurs on laisse le meme
c       algorithme car l'ajustement convectif fait le travail.
c       mais c'est delirant quand sn et snq2 n'ont pas le meme
c       signe : dans ces cas, on ne fait pas la decomposition.
c<---
c
        IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
     &      snq2(igrid,ilev)=0.E+0
        IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
     &      smq2(igrid,ilev)=0.E+0
c
C   Correction pour les couches stables.
C   Schema repris de JHoltzlag Boville, lui meme venant de...

        if (1.eq.1) then
        snstable=1.-zlev(igrid,ilev)
     s     /(700.*max(ustar(igrid),0.0001))
        snstable=1.-zlev(igrid,ilev)/400.
        snstable=max(snstable,0.)
        snstable=snstable*snstable

c abde       print*,'SN ',ilev,sn(1,ilev),snstable
        if (sn(igrid,ilev).lt.snstable) then
           sn(igrid,ilev)=snstable
           snq2(igrid,ilev)=0.
        endif

        if (sm(igrid,ilev).lt.snstable) then
           sm(igrid,ilev)=snstable
           smq2(igrid,ilev)=0.
        endif

        endif

c sn : coefficient de stabilite pour n
c snq2 : premier terme du developement limite de sn en q2
c-----------------------------------------------------------------------
                                                      ENDDO
      ENDDO
c-----------------------------------------------------------------------
c
c.......................................................................
c  calcul de km et kn au debut du pas de temps
c.......................................................................
c
                                                      DO igrid=1,ngrid 
      kn(igrid,1)=knmin
      km(igrid,1)=kmmin
      kmpre(igrid,1)=km(igrid,1)
                                                      ENDDO
c
c-----------------------------------------------------------------------
      DO ilev=2,nlev-1
                                                      DO igrid=1,ngrid 
c-----------------------------------------------------------------------
c
        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
     &                                         *sn(igrid,ilev)
        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
     &                                         *sm(igrid,ilev)
        kmpre(igrid,ilev)=km(igrid,ilev)
c
c-----------------------------------------------------------------------
                                                      ENDDO
      ENDDO
c-----------------------------------------------------------------------
c
                                                      DO igrid=1,ngrid 
      kn(igrid,nlev)=kn(igrid,nlev-1)
      km(igrid,nlev)=km(igrid,nlev-1)
      kmpre(igrid,nlev)=km(igrid,nlev)
                                                      ENDDO
c
c.......................................................................
c  boucle sur les niveaux 2 a nlev-1
c.......................................................................
c
c---->
      DO 10001 ilev=2,nlev-1
c---->
      DO 10002 igrid=1,ngrid 
c
c.......................................................................
c
c  calcul des termes sources et puits de l'equation de q2
c  ------------------------------------------------------
c
        knq3=kn(igrid,ilev)*snq2(igrid,ilev)
     &                                    /sn(igrid,ilev)
        kmq3=km(igrid,ilev)*smq2(igrid,ilev)
     &                                    /sm(igrid,ilev)
c
        termq=0.E+0
        termq3=0.E+0
        termqm2=0.E+0
        termq3m2=0.E+0
c
        tmp1=dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
        tmp2=dt*2.E+0 *kmq3*m2(igrid,ilev)
        termqm2=termqm2
     &    +dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
     &    -dt*2.E+0 *kmq3*m2(igrid,ilev)
        termq3m2=termq3m2
     &    +dt*2.E+0 *kmq3*m2(igrid,ilev)
c 
        termq=termq
     &    -dt*2.E+0 *kn(igrid,ilev)*n2(igrid,ilev)
     &    +dt*2.E+0 *knq3*n2(igrid,ilev)
        termq3=termq3
     &    -dt*2.E+0 *knq3*n2(igrid,ilev)
c
        termq3=termq3
     &    -dt*2.E+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
c
c.......................................................................
c
c  resolution stationnaire couplee avec le gradient de vitesse local
c  -----------------------------------------------------------------
c
c  -----{on cherche le cisaillement qui annule l'equation de q^2
c        supposee en q3}
c
        tmp1=termq+termq3
        tmp2=termqm2+termq3m2
        m2cstat=m2(igrid,ilev)
     &      -(tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev))
        mcstat=sqrt(m2cstat)

c  abde      print*,'M2 L=',ilev,mpre(igrid,ilev),mcstat
c
c  -----{puis on ecrit la valeur de q qui annule l'equation de m
c        supposee en q3}
c
        IF (ilev.eq.2) THEN
          kmcstat=1.E+0 / mcstat
     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
     &                        *mpre(igrid,ilev+1)
     &      +unsdz(igrid,ilev-1)
     &              *cd(igrid)
     &              *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
     &                -mcstat/unsdzdec(igrid,ilev)
     &                -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
        ELSE
          kmcstat=1.E+0 / mcstat
     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
     &                        *mpre(igrid,ilev+1)
     &      +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
     &                          *mpre(igrid,ilev-1) )
     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
        ENDIF
        tmp2=kmcstat
     &      /( sm(igrid,ilev)/q2(igrid,ilev) )
     &      /long(igrid,ilev)
        qcstat=tmp2**(1.E+0/3.E+0)
        q2cstat=qcstat**2
c
c.......................................................................
c
c  choix de la solution finale
c  ---------------------------
c
          q(igrid,ilev)=qcstat
          q2(igrid,ilev)=q2cstat
          m(igrid,ilev)=mcstat
c abd       if(ilev.le.57.and.ilev.ge.37) then
c           print*,'L=',ilev,'   M2=',m2(igrid,ilev),m2cstat,
c     s     'N2=',n2(igrid,ilev)
c abd       endif
          m2(igrid,ilev)=m2cstat
c
c --->
c       pour des raisons simples q2 est minore 
c<---
c
        IF (q2(igrid,ilev).lt.q2min) THEN
          q2(igrid,ilev)=q2min
          q(igrid,ilev)=sqrt(q2min)
        ENDIF
c
c.......................................................................
c
c  calcul final de kn et km
c  ------------------------
c
        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
     &                                           * n2(igrid,ilev)
        IF (gn.lt.gnmin) gn=gnmin
        IF (gn.gt.gnmax) gn=gnmax
        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
        sm(igrid,ilev)=
     &    (cm1+cm2*gn)
     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
     &                 *sn(igrid,ilev)
        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
     &                 *sm(igrid,ilev)
c abd
c        if(ilev.le.57.and.ilev.ge.37) then
c            print*,'L=',ilev,'   GN=',gn,'  SM=',sm(igrid,ilev)
c        endif
c
c.......................................................................
c
10002 CONTINUE
c
10001 CONTINUE
c
c.......................................................................
c
c
                                                      DO igrid=1,ngrid 
      kn(igrid,1)=knmin
      km(igrid,1)=kmmin
c     kn(igrid,1)=cd(igrid)
c     km(igrid,1)=cd(igrid)
      q2(igrid,nlev)=q2(igrid,nlev-1)
      q(igrid,nlev)=q(igrid,nlev-1)
      kn(igrid,nlev)=kn(igrid,nlev-1)
      km(igrid,nlev)=km(igrid,nlev-1)
                                                      ENDDO
c
c  CALCUL DE LA DIFFUSION VERTICALE DE Q2
      if (1.eq.1) then

        do ilev=2,klev-1
           sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
           sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
        enddo
c        print*,'Q2moy avant',sssq/sss
c       print*,'Q2q20 ',(q2(1,ilev),ilev=1,10)
c       print*,'Q2km0 ',(km(1,ilev),ilev=1,10)
c   ! C'est quoi ca qu'etait dans l'original???
c       do igrid=1,ngrid
c          q2(igrid,1)=10.
c       enddo
c        q2s=q2
c       do iii=1,10
c       call vdif_q2(dt,g,rconst,plev,temp,km,q2)
c       do ilev=1,klev+1
c          write(iii+49,*) q2(1,ilev),zlev(1,ilev)
c       enddo
c       enddo
c       stop
c       do ilev=1,klev
c          print*,zlev(1,ilev),q2s(1,ilev),q2(1,ilev)
c       enddo
c        q2s=q2-q2s
c       do ilev=1,klev
c          print*,q2s(1,ilev),zlev(1,ilev)
c       enddo
        do ilev=2,klev-1
           sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
           sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
        enddo
        print*,'Q2moy apres',sssq/sss
c
c
        do ilev=1,nlev
           do igrid=1,ngrid
              q2(igrid,ilev)=max(q2(igrid,ilev),q2min)
              q(igrid,ilev)=sqrt(q2(igrid,ilev))

c.......................................................................
c
c  calcul final de kn et km
c  ------------------------
c
        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
     &                                           * n2(igrid,ilev)
        IF (gn.lt.gnmin) gn=gnmin
        IF (gn.gt.gnmax) gn=gnmax
        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
        sm(igrid,ilev)=
     &    (cm1+cm2*gn)
     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
C   Correction pour les couches stables.
C   Schema repris de JHoltzlag Boville, lui meme venant de...

        if (1.eq.1) then
        snstable=1.-zlev(igrid,ilev)
     s     /(700.*max(ustar(igrid),0.0001))
        snstable=1.-zlev(igrid,ilev)/400.
        snstable=max(snstable,0.)
        snstable=snstable*snstable

c abde      print*,'SN ',ilev,sn(1,ilev),snstable
        if (sn(igrid,ilev).lt.snstable) then
           sn(igrid,ilev)=snstable
           snq2(igrid,ilev)=0.
        endif

        if (sm(igrid,ilev).lt.snstable) then
           sm(igrid,ilev)=snstable
           smq2(igrid,ilev)=0.
        endif

        endif

c sn : coefficient de stabilite pour n
        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
     &                 *sn(igrid,ilev)
        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
c
           enddo
        enddo
c       print*,'Q2km1 ',(km(1,ilev),ilev=1,10)

      endif

      RETURN
      END