! $Id: friction_p.F 1299 2010-01-20 14:27:21Z fairhead $

!=======================================================================
SUBROUTINE friction_loc(ucov,vcov,pdt)
  USE parallel_lmdz
  USE control_mod
  USE IOIPSL
  USE comconst_mod, ONLY: pi
  IMPLICIT NONE

  !=======================================================================

  !   Friction for the Newtonian case:
  !   --------------------------------
  !    2 possibilities (depending on flag 'friction_type'
  ! friction_type=0 : A friction that is only applied to the lowermost
  !                   atmospheric layer
  ! friction_type=1 : Friction applied on all atmospheric layer (but
  !   (default)       with stronger magnitude near the surface; see
  !                   iniacademic.F)
  !=======================================================================

  include "dimensions.h"
  include "paramet.h"
  include "comgeom2.h"
  include "iniprint.h"
  include "academic.h"

  ! arguments:
  REAL,INTENT(inout) :: ucov( iip1,jjb_u:jje_u,llm )
  REAL,INTENT(inout) :: vcov( iip1,jjb_v:jje_v,llm )
  REAL,INTENT(in) :: pdt ! time step

  ! local variables:

  REAL :: modv(iip1,jjb_u:jje_u),zco,zsi
  REAL :: vpn,vps,upoln,upols,vpols,vpoln
  REAL :: u2(iip1,jjb_u:jje_u),v2(iip1,jjb_v:jje_v)
  INTEGER :: i,j,l
  REAL,PARAMETER :: cfric=1.e-5
  LOGICAL,SAVE :: firstcall=.TRUE.
  INTEGER,SAVE :: friction_type=1
  CHARACTER(len=20) :: modname="friction_p"
  CHARACTER(len=80) :: abort_message
!$OMP THREADPRIVATE(firstcall,friction_type)
  integer :: jjb,jje

!$OMP SINGLE
  IF (firstcall) THEN
    ! set friction type
    CALL getin("friction_type",friction_type)
    if ((friction_type<0).or.(friction_type>1)) then
      abort_message="wrong friction type"
      write(lunout,*)'Friction: wrong friction type',friction_type
      CALL abort_gcm(modname,abort_message,42)
    endif
    firstcall=.FALSE.
  ENDIF
!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)

  if (friction_type==0) then ! friction on first layer only
!$OMP SINGLE
  !   calcul des composantes au carre du vent naturel
  jjb=jj_begin
  jje=jj_end+1
  if (pole_sud) jje=jj_end

  do j=jjb,jje
     do i=1,iip1
        u2(i,j)=ucov(i,j,1)*ucov(i,j,1)*unscu2(i,j)
     enddo
  enddo

  jjb=jj_begin-1
  jje=jj_end+1
  if (pole_nord) jjb=jj_begin
  if (pole_sud) jje=jj_end-1

  do j=jjb,jje
     do i=1,iip1
        v2(i,j)=vcov(i,j,1)*vcov(i,j,1)*unscv2(i,j)
     enddo
  enddo

  !   calcul du module de V en dehors des poles
  jjb=jj_begin
  jje=jj_end+1
  if (pole_nord) jjb=jj_begin+1
  if (pole_sud) jje=jj_end-1

  do j=jjb,jje
     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

  !   les deux composantes du vent au pole sont obtenues comme
  !   premiers modes de fourier de v pres du pole
  if (pole_nord) then

    upoln=0.
    vpoln=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)
       upoln=upoln+zco*vpn
       vpoln=vpoln+zsi*vpn
    enddo
    vpn=sqrt(upoln*upoln+vpoln*vpoln)/pi
    do i=1,iip1
       ! modv(i,1)=vpn
       modv(i,1)=modv(i,2)
    enddo

  endif

  if (pole_sud) then

    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))
       vps=vcov(i,jjm,1)/cv(i,jjm)
       upols=upols+zco*vps
       vpols=vpols+zsi*vps
    enddo
    vps=sqrt(upols*upols+vpols*vpols)/pi
    do i=1,iip1
     ! modv(i,jjp1)=vps
     modv(i,jjp1)=modv(i,jjm)
    enddo

  endif

  !   calcul du frottement au sol.

  jjb=jj_begin
  jje=jj_end
  if (pole_nord) jjb=jj_begin+1
  if (pole_sud) jje=jj_end-1

  do j=jjb,jje
     do i=1,iim
        ucov(i,j,1)=ucov(i,j,1) &
              -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

  jjb=jj_begin
  jje=jj_end
  if (pole_sud) jje=jj_end-1

  do j=jjb,jje
     do i=1,iip1
        vcov(i,j,1)=vcov(i,j,1) &
              -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
!$OMP END SINGLE
  endif ! of if (friction_type.eq.0)

  if (friction_type==1) then
   ! for ucov()
    jjb=jj_begin
    jje=jj_end
    if (pole_nord) jjb=jj_begin+1
    if (pole_sud) jje=jj_end-1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l=1,llm
      ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)* &
            (1.-pdt*kfrict(l))
    enddo
!$OMP END DO NOWAIT

   ! for vcoc()
    jjb=jj_begin
    jje=jj_end
    if (pole_sud) jje=jj_end-1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l=1,llm
      vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)* &
            (1.-pdt*kfrict(l))
    enddo
!$OMP END DO
  endif ! of if (friction_type.eq.1)


END SUBROUTINE friction_loc