! $Id: interpre.f90 5116 2024-07-24 12:54:37Z abarral $

SUBROUTINE interpre(q, qppm, w, fluxwppm, masse, &
        apppm, bpppm, massebx, masseby, pbaru, pbarv, &
        unatppm, vnatppm, psppm)

  USE comconst_mod, ONLY: g
  USE comvert_mod, ONLY: ap, bp
  USE lmdz_description, ONLY: descript

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "comdissip.h"
  include "comgeom2.h"

  !---------------------------------------------------
  ! Arguments
  REAL :: apppm(llm + 1), bpppm(llm + 1)
  REAL :: q(iip1, jjp1, llm), qppm(iim, jjp1, llm)
  !---------------------------------------------------
  REAL :: masse(iip1, jjp1, llm)
  REAL :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
  REAL :: w(iip1, jjp1, llm)
  REAL :: fluxwppm(iim, jjp1, llm)
  REAL :: pbaru(iip1, jjp1, llm)
  REAL :: pbarv(iip1, jjm, llm)
  REAL :: unatppm(iim, jjp1, llm)
  REAL :: vnatppm(iim, jjp1, llm)
  REAL :: psppm(iim, jjp1)
  !---------------------------------------------------
  ! Local
  REAL :: vnat(iip1, jjp1, llm)
  REAL :: unat(iip1, jjp1, llm)
  REAL :: fluxw(iip1, jjp1, llm)
  REAL :: smass(iip1, jjp1)
  !----------------------------------------------------
  INTEGER :: l, ij, i, j

  ! CALCUL DE LA PRESSION DE SURFACE
  ! Les coefficients ap et bp sont passés en common
  ! Calcul de la pression au sol en mb optimisée pour
  ! la vectorialisation

  do j = 1, jjp1
    do i = 1, iip1
      smass(i, j) = 0.
    enddo
  enddo

  do l = 1, llm
    do j = 1, jjp1
      do i = 1, iip1
        smass(i, j) = smass(i, j) + masse(i, j, l)
      enddo
    enddo
  enddo

  do j = 1, jjp1
    do i = 1, iim
      psppm(i, j) = smass(i, j) / aire(i, j) * g * 0.01
    END DO
  END DO

  ! RECONSTRUCTION DES CHAMPS CONTRAVARIANTS
  ! Le programme ppm3d travaille avec les composantes
  ! de vitesse et pas les flux, on doit donc passer de l'un à l'autre
  ! Dans le même temps, on fait le changement d'orientation du vent en v
  do l = 1, llm
    do j = 1, jjm
      do i = 1, iip1
        vnat(i, j, l) = -pbarv(i, j, l) / masseby(i, j, l) * cv(i, j)
      enddo
    enddo
    do  i = 1, iim
      vnat(i, jjp1, l) = 0.
    enddo
    do j = 1, jjp1
      do i = 1, iip1
        unat(i, j, l) = pbaru(i, j, l) / massebx(i, j, l) * cu(i, j)
      enddo
    enddo
  enddo

  ! CALCUL DU FLUX MASSIQUE VERTICAL
  ! Flux en l=1 (sol) nul
  fluxw = 0.
  do l = 1, llm
    do j = 1, jjp1
      do i = 1, iip1
        fluxw(i, j, l) = w(i, j, l) * g * 0.01 / aire(i, j)
        ! PRINT*,i,j,l,'fluxw(i,j,l)=',fluxw(i,j,l),
        ! c                      'w(i,j,l)=',w(i,j,l)
      enddo
    enddo
  enddo

  ! INVERSION DES NIVEAUX
  ! le programme ppm3d travaille avec une 3ème coordonnée inversée par rapport
  ! de celle du LMDZ: z=1<=>niveau max, z=llm+1<=>surface
  ! On passe donc des niveaux du LMDZ à ceux de Lin

  do l = 1, llm + 1
    apppm(l) = ap(llm + 2 - l)
    bpppm(l) = bp(llm + 2 - l)
  enddo

  do l = 1, llm
    do j = 1, jjp1
      do i = 1, iim
        unatppm(i, j, l) = unat(i, j, llm - l + 1)
        vnatppm(i, j, l) = vnat(i, j, llm - l + 1)
        fluxwppm(i, j, l) = fluxw(i, j, llm - l + 1)
        qppm(i, j, l) = q(i, j, llm - l + 1)
      enddo
    enddo
  enddo

  RETURN
END SUBROUTINE interpre