source: dynamico_lmdz/simple_physics/phyparam/param/vdif_k.F @ 4178

      SUBROUTINE vdif_k(ngrid,nlay,
     s   ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh)

      IMPLICIT NONE
#include "planete.h"

      INTEGER ngrid,nlay

      REAL ptimestep
      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
      REAL pz0(ngrid)
      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
      REAL pg,pcdv(ngrid)
      REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)

      INTEGER ig,il
      REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix

      REAL karman
      SAVE karman
      DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/


      print*,'LMIXMIN',lmixmin
      DO ig=1,ngrid
         pkv(ig,1)=0.
         pkh(ig,1)=0.
         pkv(ig,nlay+1)=0.
         pkh(ig,nlay+1)=0.
      ENDDO
c    s     ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)'
      DO il=2,nlay
         DO ig=1,ngrid
            z1=pzlev(ig,il)+pz0(ig)
            lmix=karman*z1/(1.+karman*z1/lmixmin)
c           lmix=lmixmin
c WARNING test lmix=lmixmin
            zdu=pu(ig,il)-pu(ig,il-1)
            zdv=pv(ig,il)-pv(ig,il-1)
            zdz=pzlay(ig,il)-pzlay(ig,il-1)
            zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
            IF(zdvodz2.LT.1.e-5) THEN
                pkv(ig,il)=lmix*sqrt(emin_turb)
            ELSE
               zri=2.*pg*(ph(ig,il)-ph(ig,il-1))
     s         / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
               pkv(ig,il)=
     s         lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
            ENDIF
            pkh(ig,il)=pkv(ig,il)
c           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
c    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
c    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
c    s      ph(ig,il),ph(ig,il-1)
         ENDDO
      ENDDO

      RETURN
      END