SUBROUTINE coefdifv( np,nlev,
     $               pu,pv,ph,pphi,
     $               pcdzv,pcdzh)
      USE constlim
      USE phys_const
      IMPLICIT NONE

c=======================================================================
c
c   calcul des coeficients de la diffusion verticale
c
c=======================================================================
c
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

#include "dimensions.h"

c   Arguments:
c   ----------

      INTEGER np,nlev
      REAL pu(np,nlev),pv(np,nlev)
      REAL ph(np,nlev),pphi(np,nlev)
      REAL pcdzv(np,nlev),pcdzh(np,nlev)

c   Local:
c   ------

      INTEGER ilev,ip
      REAL z1(np),z2(np)
      REAL z1odz(np),zdzu2(np),zhbar(np)

c-----------------------------------------------------------------------
c   initialisations:
c   ----------------

c-----------------------------------------------------------------------
c   couche de surface:
c   ------------------

      DO 210 ip=1,np
         z1(ip)=pu(ip,1)*pu(ip,1)+pv(ip,1)*pv(ip,1)
         pcdzv(ip,1)=dgcdrag(ip)*(SQRT(z1(ip))+1.)  / ph(ip,1)
         pcdzh(ip,1)=pcdzv(ip,1)
210   CONTINUE
c
c
c-----------------------------------------------------------------------
c   autres couches:
c   ---------------

      PRINT*,'Coeff k'
      DO 310 ilev=1,nlev-1
         DO 320 ip=1,np
            z1(ip)=pu(ip,ilev+1)-pu(ip,ilev)
            z2(ip)=pv(ip,ilev+1)-pv(ip,ilev)
            z1(ip)=z1(ip)*z1(ip)+z2(ip)*z2(ip)
            z1odz(ip)=g/(pphi(ip,ilev+1)-pphi(ip,ilev))
            zdzu2(ip)=z1(ip)*z1odz(ip)*z1odz(ip)
            zhbar(ip)=0.5*(ph(ip,ilev)+ph(ip,ilev+1))
            z1(ip)=( (ph(ip,ilev+1)-ph(ip,ilev))*z1odz(ip)-
     $             cpgam ) * ccdzh/zhbar(ip)
            pcdzv(ip,ilev+1)=cdzconst(ilev+1)*
     $            SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin)) /
c    $            SQRT(cdzmin) /
     $            (zhbar(ip)*zhbar(ip))
            IF(ip.eq.np/2+1) PRINT*,lmixmin*lmixmin*
     s     SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin))
            pcdzh(ip,ilev+1)=pcdzv(ip,ilev+1)
320      CONTINUE
310   CONTINUE

c-----------------------------------------------------------------------

      RETURN
      END