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

   implicit none

  include "dimensions.h"
  include "paramet.h"
  include "comdissip.h"
  include "comgeom2.h"
  include "description.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