!
! $Id: friction.F 1299 2010-01-20 14:27:21Z acaubel $
!
c=======================================================================
      SUBROUTINE friction(ucov,vcov,pdt)

      USE control_mod
 
      IMPLICIT NONE

c=======================================================================
c
c
c   Objet:
c   ------
c
c  ***********
c    Friction
c  ***********
c
c=======================================================================

#include "dimensions.h"
#include "paramet.h"
#include "comgeom2.h"
#include "comconst.h"

      REAL pdt
      REAL modv(iip1,jjp1),zco,zsi
      REAL vpn,vps,upoln,upols,vpols,vpoln
      REAL u2(iip1,jjp1),v2(iip1,jjm)
      REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
      INTEGER  i,j
      REAL cfric
      parameter (cfric=1.e-5)


c   calcul des composantes au carre du vent naturel
      do j=1,jjp1
         do i=1,iip1
            u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
         enddo
      enddo
      do j=1,jjm
         do i=1,iip1
            v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
         enddo
      enddo

c   calcul du module de V en dehors des poles
      do j=2,jjm
         do i=2,iip1
            modv(i,j)=sqrt(0.5*(u2(i-1,j)+u2(i,j)+v2(i,j-1)+v2(i,j)))
         enddo
         modv(1,j)=modv(iip1,j)
      enddo

c   les deux composantes du vent au pole sont obtenues comme
c   premiers modes de fourier de v pres du pole
      upoln=0.
      vpoln=0.
      upols=0.
      vpols=0.
      do i=2,iip1
         zco=cos(rlonv(i))*(rlonu(i)-rlonu(i-1))
         zsi=sin(rlonv(i))*(rlonu(i)-rlonu(i-1))
         vpn=vcov(i,1,1)/cv(i,1)
         vps=vcov(i,jjm,1)/cv(i,jjm)
         upoln=upoln+zco*vpn
         vpoln=vpoln+zsi*vpn
         upols=upols+zco*vps
         vpols=vpols+zsi*vps
      enddo
      vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
      vps=sqrt(upols*upols+vpols*vpols)/pi
      do i=1,iip1
c        modv(i,1)=vpn
c        modv(i,jjp1)=vps
         modv(i,1)=modv(i,2)
         modv(i,jjp1)=modv(i,jjm)
      enddo

c   calcul du frottement au sol.
      do j=2,jjm
         do i=1,iim
            ucov(i,j,1)=ucov(i,j,1)
     s      -cfric*pdt*0.5*(modv(i+1,j)+modv(i,j))*ucov(i,j,1)
         enddo
         ucov(iip1,j,1)=ucov(1,j,1)
      enddo
      do j=1,jjm
         do i=1,iip1
            vcov(i,j,1)=vcov(i,j,1)
     s      -cfric*pdt*0.5*(modv(i,j+1)+modv(i,j))*vcov(i,j,1)
         enddo
         vcov(iip1,j,1)=vcov(1,j,1)
      enddo

      RETURN
      END