! $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
  USE lmdz_iniprint, ONLY: lunout, prt_level
  USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4

  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"

  ! 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