!
! $Id: addfi.f90 5246 2024-10-21 12:58:45Z 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

  RETURN
END SUBROUTINE addfi