! $Id: addfi.F90 5117 2024-07-24 14:23:34Z abarral $

SUBROUTINE addfi(pdt, leapf, forward, &
        pucov, pvcov, pteta, pq, pps, &
        pdufi, pdvfi, pdhfi, pdqfi, pdpfi)

  USE infotrac, ONLY: nqtot
  USE control_mod, ONLY: planet_type
  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"
  include "comgeom.h"
  !
  !    Arguments :
  !    -----------
  !
  REAL, INTENT(IN) :: pdt ! time step for the integration (s)
  !
  REAL, INTENT(INOUT) :: pvcov(ip1jm, llm) ! covariant meridional wind
  REAL, INTENT(INOUT) :: pucov(ip1jmp1, llm) ! covariant zonal wind
  REAL, INTENT(INOUT) :: pteta(ip1jmp1, llm) ! potential temperature
  REAL, INTENT(INOUT) :: pq(ip1jmp1, llm, nqtot) ! tracers
  REAL, INTENT(INOUT) :: pps(ip1jmp1) ! surface pressure (Pa)
  ! respective tendencies (.../s) to add
  REAL, INTENT(IN) :: pdvfi(ip1jm, llm)
  REAL, INTENT(IN) :: pdufi(ip1jmp1, llm)
  REAL, INTENT(IN) :: pdqfi(ip1jmp1, llm, nqtot)
  REAL, INTENT(IN) :: pdhfi(ip1jmp1, llm)
  REAL, INTENT(IN) :: pdpfi(ip1jmp1)
  !
  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

  REAL :: SSUM
  !
  !-----------------------------------------------------------------------

  DO k = 1, llm
    DO j = 1, ip1jmp1
      pteta(j, k) = pteta(j, k) + pdhfi(j, k) * pdt
    ENDDO
  ENDDO

  DO  k = 1, llm
    DO  ij = 1, iim
      xpn(ij) = aire(ij) * pteta(ij, k)
      xps(ij) = aire(ij + ip1jm) * pteta(ij + ip1jm, k)
    ENDDO
    tpn = SSUM(iim, xpn, 1) / apoln
    tps = SSUM(iim, xps, 1) / apols

    DO ij = 1, iip1
      pteta(ij, k) = tpn
      pteta(ij + ip1jm, k) = tps
    ENDDO
  ENDDO
  !

  DO k = 1, llm
    DO j = iip2, ip1jm
      pucov(j, k) = pucov(j, k) + pdufi(j, k) * pdt
    ENDDO
  ENDDO

  DO k = 1, llm
    DO j = 1, ip1jm
      pvcov(j, k) = pvcov(j, k) + pdvfi(j, k) * pdt
    ENDDO
  ENDDO

  !
  DO j = 1, ip1jmp1
    pps(j) = pps(j) + pdpfi(j) * pdt
  ENDDO

  IF (planet_type=="earth") THEN
    ! earth case, special treatment for first 2 tracers (water)
    DO iq = 1, 2
      DO k = 1, llm
        DO j = 1, ip1jmp1
          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
    ENDDO

    DO iq = 3, nqtot
      DO k = 1, llm
        DO j = 1, ip1jmp1
          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
    ENDDO
  else
    ! general case, treat all tracers equally)
    DO iq = 1, nqtot
      DO k = 1, llm
        DO j = 1, ip1jmp1
          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
    ENDDO
  ENDIF ! of if (planet_type=="earth")

  DO  ij = 1, iim
    xpn(ij) = aire(ij) * pps(ij)
    xps(ij) = aire(ij + ip1jm) * pps(ij + ip1jm)
  ENDDO
  tpn = SSUM(iim, xpn, 1) / apoln
  tps = SSUM(iim, xps, 1) / apols

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

  DO iq = 1, nqtot
    DO  k = 1, llm
      DO  ij = 1, iim
        xpn(ij) = aire(ij) * pq(ij, k, iq)
        xps(ij) = aire(ij + ip1jm) * pq(ij + ip1jm, k, iq)
      ENDDO
      tpn = SSUM(iim, xpn, 1) / apoln
      tps = SSUM(iim, xps, 1) / apols

      DO ij = 1, iip1
        pq (ij, k, iq) = tpn
        pq (ij + ip1jm, k, iq) = tps
      ENDDO
    ENDDO
  ENDDO


END SUBROUTINE addfi