! $Id$

SUBROUTINE addfi_loc(pdt, leapf, forward, &
        pucov, pvcov, pteta, pq, pps, &
        pdufi, pdvfi, pdhfi, pdqfi, pdpfi)
  USE parallel_lmdz
  USE infotrac, ONLY: nqtot
  USE control_mod, ONLY: planet_type
  USE lmdz_ssum_scopy, ONLY: ssum
  USE lmdz_comgeom

  IMPLICIT NONE
  !
  !=======================================================================
  !
  !    Addition of the physical tendencies
  !
  !    Interface :
  !    -----------
  !
  !  Input :
  !  -------
  !  pdt                    time step of integration
  !  leapf                  logical
  !  forward                logical
  !  pucov(ip1jmp1,llm)     first component of the covariant velocity
  !  pvcov(ip1ip1jm,llm)    second component of the covariant velocity
  !  pteta(ip1jmp1,llm)     potential temperature
  !  pts(ip1jmp1,llm)       surface temperature
  !  pdufi(ip1jmp1,llm)     |
  !  pdvfi(ip1jm,llm)       |   respective
  !  pdhfi(ip1jmp1)         |      tendencies
  !  pdtsfi(ip1jmp1)        |
  !
  !  Output :
  !  --------
  !  pucov
  !  pvcov
  !  ph
  !  pts
  !
  !
  !=======================================================================
  !
  !-----------------------------------------------------------------------
  !
  !    0.  Declarations :
  !    ------------------
  !
  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  !
  !    Arguments :
  !    -----------
  !
  REAL, INTENT(IN) :: pdt ! time step for the integration (s)
  !
  REAL, INTENT(INOUT) :: pvcov(ijb_v:ije_v, llm) ! covariant meridional wind
  REAL, INTENT(INOUT) :: pucov(ijb_u:ije_u, llm) ! covariant zonal wind
  REAL, INTENT(INOUT) :: pteta(ijb_u:ije_u, llm) ! potential temperature
  REAL, INTENT(INOUT) :: pq(ijb_u:ije_u, llm, nqtot) ! tracers
  REAL, INTENT(INOUT) :: pps(ijb_u:ije_u) ! surface pressure (Pa)
  ! respective tendencies (.../s) to add
  REAL, INTENT(IN) :: pdvfi(ijb_v:ije_v, llm)
  REAL, INTENT(IN) :: pdufi(ijb_u:ije_u, llm)
  REAL, INTENT(IN) :: pdqfi(ijb_u:ije_u, llm, nqtot)
  REAL, INTENT(IN) :: pdhfi(ijb_u:ije_u, llm)
  REAL, INTENT(IN) :: pdpfi(ijb_u:ije_u)
  !
  LOGICAL, INTENT(IN) :: leapf, forward ! not used
  !
  !
  !    Local variables :
  !    -----------------
  !
  REAL :: xpn(iim), xps(iim), tpn, tps
  INTEGER :: j, k, iq, ij
  REAL, PARAMETER :: qtestw = 1.0e-15
  REAL, PARAMETER :: qtestt = 1.0e-40

  INTEGER :: ijb, ije
  !
  !-----------------------------------------------------------------------

  ijb = ij_begin
  ije = ij_end

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO k = 1, llm
    DO j = ijb, ije
      pteta(j, k) = pteta(j, k) + pdhfi(j, k) * pdt
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  IF (pole_nord) THEN
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO  k = 1, llm
      DO  ij = 1, iim
        xpn(ij) = aire(ij) * pteta(ij, k)
      ENDDO
      tpn = SSUM(iim, xpn, 1) / apoln

      DO ij = 1, iip1
        pteta(ij, k) = tpn
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT
  ENDIF

  IF (pole_sud) THEN
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO  k = 1, llm
      DO  ij = 1, iim
        xps(ij) = aire(ij + ip1jm) * pteta(ij + ip1jm, k)
      ENDDO
      tps = SSUM(iim, xps, 1) / apols

      DO ij = 1, iip1
        pteta(ij + ip1jm, k) = tps
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT
  ENDIF
  !

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

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO k = 1, llm
    DO j = ijb, ije
      pucov(j, k) = pucov(j, k) + pdufi(j, k) * pdt
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  IF (pole_nord) ijb = ij_begin

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO k = 1, llm
    DO j = ijb, ije
      pvcov(j, k) = pvcov(j, k) + pdvfi(j, k) * pdt
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  !
  IF (pole_sud)  ije = ij_end
  !$OMP MASTER
  DO j = ijb, ije
    pps(j) = pps(j) + pdpfi(j) * pdt
  ENDDO
  !$OMP END MASTER

  IF (planet_type=="earth") THEN
    ! earth case, special treatment for first 2 tracers (water)
    DO iq = 1, 2
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO j = ijb, ije
          pq(j, k, iq) = pq(j, k, iq) + pdqfi(j, k, iq) * pdt
          pq(j, k, iq) = AMAX1(pq(j, k, iq), qtestw)
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT
    ENDDO

    DO iq = 3, nqtot
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO j = ijb, ije
          pq(j, k, iq) = pq(j, k, iq) + pdqfi(j, k, iq) * pdt
          pq(j, k, iq) = AMAX1(pq(j, k, iq), qtestt)
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT
    ENDDO
  else
    ! general case, treat all tracers equally)
    DO iq = 1, nqtot
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO j = ijb, ije
          pq(j, k, iq) = pq(j, k, iq) + pdqfi(j, k, iq) * pdt
          pq(j, k, iq) = AMAX1(pq(j, k, iq), qtestt)
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT
    ENDDO
  ENDIF ! of if (planet_type=="earth")

  !$OMP MASTER
  IF (pole_nord) THEN
    DO  ij = 1, iim
      xpn(ij) = aire(ij) * pps(ij)
    ENDDO

    tpn = SSUM(iim, xpn, 1) / apoln

    DO ij = 1, iip1
      pps (ij) = tpn
    ENDDO

  ENDIF

  IF (pole_sud) THEN
    DO  ij = 1, iim
      xps(ij) = aire(ij + ip1jm) * pps(ij + ip1jm)
    ENDDO

    tps = SSUM(iim, xps, 1) / apols

    DO ij = 1, iip1
      pps (ij + ip1jm) = tps
    ENDDO

  ENDIF
  !$OMP END MASTER

  IF (pole_nord) THEN
    DO iq = 1, nqtot
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO  k = 1, llm
        DO  ij = 1, iim
          xpn(ij) = aire(ij) * pq(ij, k, iq)
        ENDDO
        tpn = SSUM(iim, xpn, 1) / apoln

        DO ij = 1, iip1
          pq (ij, k, iq) = tpn
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT
    ENDDO
  ENDIF

  IF (pole_sud) THEN
    DO iq = 1, nqtot
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO  k = 1, llm
        DO  ij = 1, iim
          xps(ij) = aire(ij + ip1jm) * pq(ij + ip1jm, k, iq)
        ENDDO
        tps = SSUM(iim, xps, 1) / apols

        DO ij = 1, iip1
          pq (ij + ip1jm, k, iq) = tps
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT
    ENDDO
  ENDIF

END SUBROUTINE addfi_loc