! $Header$

SUBROUTINE advect_new_loc(ucov, vcov, teta, w, massebx, masseby, &
        du, dv, dteta)
  USE parallel_lmdz
  USE write_field_loc
  USE advect_new_mod
  USE comconst_mod, ONLY: daysec
  USE logic_mod, ONLY: conser
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO

  IMPLICIT NONE
  !=======================================================================
  !
  !   Auteurs:  P. Le Van , Fr. Hourdin  .
  !   -------
  !
  !   Objet:
  !   ------
  !
  !   *************************************************************
  !   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
  !   *************************************************************
  !    ces termes sont ajoutes a du,dv,dteta et dq .
  !  Modif F.Forget 03/94 : on retire q de advect
  !
  !=======================================================================
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

  include "dimensions.h"
  include "paramet.h"
  include "comgeom.h"

  !   Arguments:
  !   ----------

  REAL :: vcov(ijb_v:ije_v, llm), ucov(ijb_u:ije_u, llm)
  REAL :: teta(ijb_u:ije_u, llm)
  REAL :: massebx(ijb_u:ije_u, llm), masseby(ijb_v:ije_v, llm)
  REAL :: w(ijb_u:ije_u, llm)
  REAL :: dv(ijb_v:ije_v, llm), du(ijb_u:ije_u, llm)
  REAL :: dteta(ijb_u:ije_u, llm)
  !   Local:
  !   ------

  REAL :: wsur2(ijb_u:ije_u)
  REAL :: unsaire2(ijb_u:ije_u), ge(ijb_u:ije_u)
  REAL :: deuxjour, ww, gt, uu, vv

  INTEGER :: ij, l, ijb, ije

  !-----------------------------------------------------------------------
  !   2. Calculs preliminaires:
  !   -------------------------

  IF (conser.AND.1==0)  THEN
    deuxjour = 2. * daysec

    DO   ij = 1, ip1jmp1
      unsaire2(ij) = unsaire(ij) * unsaire(ij)
    END DO
  END IF


  !------------------  -yy ----------------------------------------------
  !   .  Calcul de     u

  !$OMP MASTER
  ijb = ij_begin
  ije = ij_end
  IF (pole_nord) ijb = ijb + iip1
  IF (pole_sud)  ije = ije - iip1

  DO ij = ijb, ije
    du2(ij, 1) = 0.
    du1(ij, llm) = 0.
  ENDDO

  ijb = ij_begin
  ije = ij_end
  IF (pole_sud)  ije = ij_end - iip1

  DO ij = ijb, ije
    dv2(ij, 1) = 0.
    dv1(ij, llm) = 0.
  ENDDO

  ijb = ij_begin
  ije = ij_end

  DO ij = ijb, ije
    dteta2(ij, 1) = 0.
    dteta1(ij, llm) = 0.
  ENDDO
  !$OMP END MASTER
  !$OMP BARRIER

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO  l = 1, llm

    ijb = ij_begin
    ije = ij_end
    IF (pole_nord) ijb = ijb + iip1
    IF (pole_sud)  ije = ije - iip1

    ! DO    ij     = iip2, ip1jmp1
    !    uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
    ! ENDDO

    ! DO    ij     = iip2, ip1jm
    !    uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
    ! ENDDO

    DO    ij = ijb, ije

      uav(ij, l) = 0.25 * (ucov(ij, l) + ucov(ij - iip1, l)) &
              + 0.25 * (ucov(ij + iip1, l) + ucov(ij, l))
    ENDDO

    IF (pole_nord) THEN
      DO      ij = 1, iip1
        uav(ij, l) = 0.
      ENDDO
    endif

    IF (pole_sud) THEN
      DO      ij = 1, iip1
        uav(ip1jm + ij, l) = 0.
      ENDDO
    endif

  ENDDO
  !$OMP END DO
  ! CALL write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))

  !------------------  -xx ----------------------------------------------
  !   .  Calcul de     v

  ijb = ij_begin
  ije = ij_end
  IF (pole_sud)  ije = ij_end - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO  l = 1, llm

    DO    ij = ijb + 1, ije
      vav(ij, l) = 0.25 * (vcov(ij, l) + vcov(ij - 1, l))
    ENDDO

    DO    ij = ijb, ije, iip1
      vav(ij, l) = vav(ij + iim, l)
    ENDDO

    DO    ij = ijb, ije - 1
      vav(ij, l) = vav(ij, l) + vav(ij + 1, l)
    ENDDO

    DO    ij = ijb, ije, iip1
      vav(ij + iim, l) = vav(ij, l)
    ENDDO

  ENDDO
  !$OMP END DO
  ! CALL write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))

  !-----------------------------------------------------------------------
  !$OMP BARRIER

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llmm1


    ! ......   calcul de  - w/2.    au niveau  l+1   .......
    ijb = ij_begin
    ije = ij_end + iip1
    IF (pole_sud)  ije = ij_end

    DO ij = ijb, ije
      wsur2(ij) = - 0.5 * w(ij, l + 1)
    END DO


    ! .....................     calcul pour  du     ..................

    ijb = ij_begin
    ije = ij_end
    IF (pole_nord) ijb = ijb + iip1
    IF (pole_sud)  ije = ije - iip1

    DO ij = ijb, ije - 1
      ww = wsur2 (ij) + wsur2(ij + 1)
      uu = 0.5 * (ucov(ij, l) + ucov(ij, l + 1))
      du1(ij, l) = ww * (uu - uav(ij, l)) / massebx(ij, l)
      du2(ij, l + 1) = ww * (uu - uav(ij, l + 1)) / massebx(ij, l + 1)
    END DO

    ! .................    calcul pour   dv      .....................
    ijb = ij_begin
    ije = ij_end
    IF (pole_sud)  ije = ij_end - iip1

    DO ij = ijb, ije
      ww = wsur2(ij + iip1) + wsur2(ij)
      vv = 0.5 * (vcov(ij, l) + vcov(ij, l + 1))
      dv1(ij, l) = ww * (vv - vav(ij, l)) / masseby(ij, l)
      dv2(ij, l + 1) = ww * (vv - vav(ij, l + 1)) / masseby(ij, l + 1)
    END DO

    !

    ! ............................................................
    ! ...............    calcul pour   dh      ...................
    ! ............................................................

    !                   ---z
    !   calcul de  - d( teta  * w )      qu'on ajoute a   dh
    !               ...............
    ijb = ij_begin
    ije = ij_end

    DO ij = ijb, ije
      ww = wsur2(ij) * (teta(ij, l) + teta(ij, l + 1))
      dteta1(ij, l) = ww
      dteta2(ij, l + 1) = ww
    END DO

    ! ym ---> conser a voir plus tard

    ! IF( conser)  THEN
    !
    !    DO 17 ij = 1,ip1jmp1
    !    ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
    !  17    CONTINUE
    !    gt       = SSUM( ip1jmp1,ge,1 )
    !    gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
    !  END IF

  END DO
  !$OMP END DO

  ijb = ij_begin
  ije = ij_end
  IF (pole_nord) ijb = ijb + iip1
  IF (pole_sud)  ije = ije - iip1
  IF (CPPKEY_DEBUGIO) THEN
    CALL WriteField_u('du_bis', du)
  END IF
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije - 1
      du(ij, l) = du(ij, l) + du2(ij, l) - du1(ij, l)
    ENDDO

    DO   ij = ijb + iip1 - 1, ije, iip1
      du(ij, l) = du(ij - iim, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  IF (CPPKEY_DEBUGIO) THEN
    CALL WriteField_u('du1', du1)
    CALL WriteField_u('du2', du2)
    CALL WriteField_u('du_bis', du)
  END IF
  ijb = ij_begin
  ije = ij_end
  IF (pole_sud)  ije = ij_end - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      dv(ij, l) = dv(ij, l) + dv2(ij, l) - dv1(ij, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  ijb = ij_begin
  ije = ij_end

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      dteta(ij, l) = dteta(ij, l) + dteta2(ij, l) - dteta1(ij, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT


END SUBROUTINE advect_new_loc