source: LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90 @ 5116

! $Id: integrd.F90 5116 2024-07-24 12:54:37Z abarral $

SUBROUTINE integrd &
        (nq, vcovm1, ucovm1, tetam1, psm1, massem1, &
        dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis & !,finvmaold
        )

  use control_mod, ONLY: planet_type
  use comconst_mod, ONLY: pi
  USE logic_mod, ONLY: leapf
  use comvert_mod, ONLY: ap, bp
  USE temps_mod, ONLY: dt

  IMPLICIT NONE


  !=======================================================================
  !
  !   Auteur:  P. Le Van
  !   -------
  !
  !   objet:
  !   ------
  !
  !   Incrementation des tendances dynamiques
  !
  !=======================================================================
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

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

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

  integer, intent(in) :: nq ! number of tracers to handle in this routine
  real, intent(inout) :: vcov(ip1jm, llm) ! covariant meridional wind
  real, intent(inout) :: ucov(ip1jmp1, llm) ! covariant zonal wind
  real, intent(inout) :: teta(ip1jmp1, llm) ! potential temperature
  real, intent(inout) :: q(ip1jmp1, llm, nq) ! advected tracers
  real, intent(inout) :: ps(ip1jmp1) ! surface pressure
  real, intent(inout) :: masse(ip1jmp1, llm) ! atmospheric mass
  real, intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
  ! values at previous time step
  real, intent(inout) :: vcovm1(ip1jm, llm)
  real, intent(inout) :: ucovm1(ip1jmp1, llm)
  real, intent(inout) :: tetam1(ip1jmp1, llm)
  real, intent(inout) :: psm1(ip1jmp1)
  real, intent(inout) :: massem1(ip1jmp1, llm)
  ! the tendencies to add
  real, intent(in) :: dv(ip1jm, llm)
  real, intent(in) :: du(ip1jmp1, llm)
  real, intent(in) :: dteta(ip1jmp1, llm)
  real, intent(in) :: dp(ip1jmp1)
  real, intent(in) :: dq(ip1jmp1, llm, nq) !!! unused
  ! real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused

  !   Local:
  !   ------

  REAL :: vscr(ip1jm), uscr(ip1jmp1), hscr(ip1jmp1), pscr(ip1jmp1)
  REAL :: massescr(ip1jmp1, llm)
  ! REAL finvmasse(ip1jmp1,llm)
  REAL :: p(ip1jmp1, llmp1)
  REAL :: tpn, tps, tppn(iim), tpps(iim)
  REAL :: qpn, qps, qppn(iim), qpps(iim)
  REAL :: deltap(ip1jmp1, llm)

  INTEGER :: l, ij, iq, i, j

  REAL :: SSUM

  !-----------------------------------------------------------------------

  DO  l = 1, llm
    DO  ij = 1, iip1
      ucov(ij, l) = 0.
      ucov(ij + ip1jm, l) = 0.
      uscr(ij) = 0.
      uscr(ij + ip1jm) = 0.
    ENDDO
  ENDDO


  !    ............    integration  de       ps         ..............

  CALL SCOPY(ip1jmp1 * llm, masse, 1, massescr, 1)

  DO ij = 1, ip1jmp1
    pscr (ij) = ps(ij)
    ps (ij) = psm1(ij) + dt * dp(ij)
  ENDDO
  !
  DO ij = 1, ip1jmp1
    IF(ps(ij)<0.) THEN
      WRITE(lunout, *) "integrd: negative surface pressure ", ps(ij)
      WRITE(lunout, *) " at node ij =", ij
      ! since ij=j+(i-1)*jjp1 , we have
      j = modulo(ij, jjp1)
      i = 1 + (ij - j) / jjp1
      WRITE(lunout, *) " lon = ", rlonv(i) * 180. / pi, " deg", &
              " lat = ", rlatu(j) * 180. / pi, " deg"
      CALL abort_gcm("integrd", "", 1)
    ENDIF
  ENDDO
  !
  DO  ij = 1, iim
    tppn(ij) = aire(ij) * ps(ij)
    tpps(ij) = aire(ij + ip1jm) * ps(ij + ip1jm)
  ENDDO
  tpn = SSUM(iim, tppn, 1) / apoln
  tps = SSUM(iim, tpps, 1) / apols
  DO ij = 1, iip1
    ps(ij) = tpn
    ps(ij + ip1jm) = tps
  ENDDO
  !
  !  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
  !
  CALL pression (ip1jmp1, ap, bp, ps, p)
  CALL massdair (p, masse)

  ! Ehouarn : we don't use/need finvmaold and finvmasse,
  ! so might as well not compute them
  ! CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
  ! CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
  !

  !    ............   integration  de  ucov, vcov,  h     ..............

  DO l = 1, llm

    DO ij = iip2, ip1jm
      uscr(ij) = ucov(ij, l)
      ucov(ij, l) = ucovm1(ij, l) + dt * du(ij, l)
    ENDDO

    DO ij = 1, ip1jm
      vscr(ij) = vcov(ij, l)
      vcov(ij, l) = vcovm1(ij, l) + dt * dv(ij, l)
    ENDDO

    DO ij = 1, ip1jmp1
      hscr(ij) = teta(ij, l)
      teta (ij, l) = tetam1(ij, l) * massem1(ij, l) / masse(ij, l) &
              + dt * dteta(ij, l) / masse(ij, l)
    ENDDO

    !   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
    !
    !
    DO  ij = 1, iim
      tppn(ij) = aire(ij) * teta(ij, l)
      tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
    ENDDO
    tpn = SSUM(iim, tppn, 1) / apoln
    tps = SSUM(iim, tpps, 1) / apols

    DO ij = 1, iip1
      teta(ij, l) = tpn
      teta(ij + ip1jm, l) = tps
    ENDDO
    !

    IF(leapf)  THEN
      CALL SCOPY (ip1jmp1, uscr(1), 1, ucovm1(1, l), 1)
      CALL SCOPY (ip1jm, vscr(1), 1, vcovm1(1, l), 1)
      CALL SCOPY (ip1jmp1, hscr(1), 1, tetam1(1, l), 1)
    END IF

  ENDDO ! of DO l = 1,llm


  !
  !   .......  integration de   q   ......
  !
  !$$$      IF( iadv(1).NE.3.AND.iadv(2).NE.3 )    THEN
  !$$$c
  !$$$       IF( forward .OR.  leapf )  THEN
  !$$$        DO iq = 1,2
  !$$$        DO  l = 1,llm
  !$$$        DO ij = 1,ip1jmp1
  !$$$        q(ij,l,iq) = ( q(ij,l,iq)*finvmaold(ij,l) + dtvr *dq(ij,l,iq) )/
  !$$$     $                            finvmasse(ij,l)
  !$$$        ENDDO
  !$$$        ENDDO
  !$$$        ENDDO
  !$$$       ELSE
  !$$$         DO iq = 1,2
  !$$$         DO  l = 1,llm
  !$$$         DO ij = 1,ip1jmp1
  !$$$         q( ij,l,iq ) = q( ij,l,iq ) * finvmaold(ij,l) / finvmasse(ij,l)
  !$$$         ENDDO
  !$$$         ENDDO
  !$$$         ENDDO
  !$$$
  !$$$       END IF
  !$$$c
  !$$$      ENDIF

  if (planet_type=="earth") THEN
    ! Earth-specific treatment of first 2 tracers (water)
    DO l = 1, llm
      DO ij = 1, ip1jmp1
        deltap(ij, l) = p(ij, l) - p(ij, l + 1)
      ENDDO
    ENDDO

    CALL qminimum(q, nq, deltap)

    !
    !    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
    !

    DO iq = 1, nq
      DO l = 1, llm

        DO ij = 1, iim
          qppn(ij) = aire(ij) * q(ij, l, iq)
          qpps(ij) = aire(ij + ip1jm) * q(ij + ip1jm, l, iq)
        ENDDO
        qpn = SSUM(iim, qppn, 1) / apoln
        qps = SSUM(iim, qpps, 1) / apols

        DO ij = 1, iip1
          q(ij, l, iq) = qpn
          q(ij + ip1jm, l, iq) = qps
        ENDDO

      ENDDO
    ENDDO

    ! Ehouarn: forget about finvmaold
    ! CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )

  endif ! of if (planet_type.eq."earth")
  !
  !
  ! .....   FIN  de l'integration  de   q    .......

  !    .................................................................

  IF(leapf)  THEN
    CALL SCOPY (ip1jmp1, pscr, 1, psm1, 1)
    CALL SCOPY (ip1jmp1 * llm, massescr, 1, massem1, 1)
  END IF


END SUBROUTINE integrd