SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
      IMPLICIT NONE
c=======================================================================
c
c   Subject: computation of the surface drag coefficient using the
c   -------  approch developed by Loui for ECMWF.
c
c   Author: Frederic Hourdin  15 /10 /93
c   -------
c
c   Arguments:
c   ----------
c
c   inputs:
c   ------
c     ngrid            size of the horizontal grid
c     pg               gravity (m s -2)
c     pz(ngrid)        height of the first atmospheric layer
c     pu(ngrid)        u component of the wind in that layer
c     pv(ngrid)        v component of the wind in that layer
c     pts(ngrid)       surfacte temperature
c     ph(ngrid)        potential temperature T*(p/ps)^kappa
c
c   outputs:
c   --------
c     pcdv(ngrid)      Cd for the wind
c     pcdh(ngrid)      Cd for potential temperature
c
c=======================================================================
c
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

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

      INTEGER ngrid,nlay
      REAL pz0(ngrid)
      REAL pg,pz(ngrid)
      REAL pu(ngrid),pv(ngrid)
      REAL pts(ngrid),ph(ngrid)
      REAL pcdv(ngrid),pcdh(ngrid)

c   Local:
c   ------

      INTEGER ig

      REAL zu2,z1,zri,zcd0,zz

      REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2
      LOGICAL firstcal
      DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
      DATA firstcal/.true./
      SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2

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

c     DO ig=1,ngrid
c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
c        pcdh(ig)=pcdv(ig)
c     ENDDO
c     RETURN

      IF (firstcal) THEN
         c2b=2.*b
         c3bc=3.*b*c
         c3b=3.*b
         firstcal=.false.
      ENDIF

c!!!! WARNING, verifier la formule originale de Louis!
      DO ig=1,ngrid
         zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
         zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
         z1=1.+pz(ig)/pz0(ig)
         zcd0=karman/log(z1)
         zcd0=zcd0*zcd0*sqrt(zu2)
         IF(zri.LT.0.) THEN
            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
            pcdv(ig)=zcd0*(1.-2.*z1)
            pcdh(ig)=zcd0*(1.-3.*z1)
         ELSE
            zz=sqrt(1.+d*zri)
            pcdv(ig)=zcd0/(1.+c2b*zri/zz)
            pcdh(ig)=zcd0/(1.+c3b*zri*zz)
         ENDIF
      ENDDO

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

      RETURN
      END