MODULE vdif_mod
  IMPLICIT NONE
  SAVE

CONTAINS

  SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
    USE planet, ONLY : karman
    !=======================================================================
    !
    !   Subject: computation of the surface drag coefficient using the
    !   -------  approch developed by Loui for ECMWF.
    !
    !   Author: Frederic Hourdin  15 /10 /93
    !   -------
    !
    !   Arguments:
    !   ----------
    !
    !   inputs:
    !   ------
    !     ngrid            size of the horizontal grid
    !     pg               gravity (m s -2)
    !     pz(ngrid)        height of the first atmospheric layer
    !     pu(ngrid)        u component of the wind in that layer
    !     pv(ngrid)        v component of the wind in that layer
    !     pts(ngrid)       surfacte temperature
    !     ph(ngrid)        potential temperature T*(p/ps)^kappa
    !
    !   outputs:
    !   --------
    !     pcdv(ngrid)      Cd for the wind
    !     pcdh(ngrid)      Cd for potential temperature
    !
    !=======================================================================
    !
    !-----------------------------------------------------------------------
    !   Declarations:
    !   -------------
    
    !   Arguments:
    !   ----------
    
    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)
    
    !   Local:
    !   ------
    
    INTEGER ig
    
    REAL zu2,z1,zri,zcd0,zz
    
    REAL b,c,d,c2b,c3bc,c3b,z0,umin2
    LOGICAL firstcal
    DATA b,c,d,umin2/5.,5.,5.,1.e-12/
    DATA firstcal/.true./
    SAVE b,c,d,c2b,c3bc,c3b,firstcal,umin2
    
    !-----------------------------------------------------------------------
    !   couche de surface:
    !   ------------------
    
    !     DO ig=1,ngrid
    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    !        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
    !        pcdh(ig)=pcdv(ig)
    !     ENDDO
    !     RETURN
    
    IF (firstcal) THEN
       c2b=2.*b
       c3bc=3.*b*c
       c3b=3.*b
       firstcal=.false.
    ENDIF
    
!!!! 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
    
    !-----------------------------------------------------------------------
    
    RETURN
  END SUBROUTINE vdif_cd
  
  
  SUBROUTINE vdif_k(ngrid,nlay,   &
       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh)
    ! FIXME : pkh :: pkv
    USE planet
    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
    
    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
    
    DO il=2,nlay
       DO ig=1,ngrid
          z1=pzlev(ig,il)+pz0(ig)
          lmix=karman*z1/(1.+karman*z1/lmixmin)
          !           lmix=lmixmin
          ! 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)) &
                  / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
             pkv(ig,il)= &
                  lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
          ENDIF
          pkh(ig,il)=pkv(ig,il)
          !           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
          !    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
          !    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
          !    s      ph(ig,il),ph(ig,il-1)
       ENDDO
    ENDDO
    
    RETURN
  END SUBROUTINE vdif_k
  
END MODULE vdif_mod