! $Id: traceurpole.f90 5134 2024-07-26 15:56:37Z abarral $

SUBROUTINE traceurpole(q, masse)
  USE lmdz_description, ONLY: descript
  USE lmdz_comdissip, ONLY: coefdis, tetavel, tetatemp, gamdissip, niterdis

  IMPLICIT NONE

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  INCLUDE "comgeom2.h"


  !   Arguments
  INTEGER :: iq
  REAL :: masse(iip1, jjp1, llm)
  REAL :: q(iip1, jjp1, llm)


  !   Locals
  INTEGER :: i, j, l
  REAL :: sommemassen(llm)
  REAL :: sommemqn(llm)
  REAL :: sommemasses(llm)
  REAL :: sommemqs(llm)
  REAL :: qpolen(llm), qpoles(llm)


  ! On impose une seule valeur au pôle Sud j=jjm+1=jjp1
  sommemasses = 0
  sommemqs = 0
  do l = 1, llm
    do i = 1, iip1
      sommemasses(l) = sommemasses(l) + masse(i, jjp1, l)
      sommemqs(l) = sommemqs(l) + masse(i, jjp1, l) * q(i, jjp1, l)
    enddo
    qpoles(l) = sommemqs(l) / sommemasses(l)
  enddo

  ! On impose une seule valeur du traceur au pôle Nord j=1
  sommemassen = 0
  sommemqn = 0
  do l = 1, llm
    do i = 1, iip1
      sommemassen(l) = sommemassen(l) + masse(i, 1, l)
      sommemqn(l) = sommemqn(l) + masse(i, 1, l) * q(i, 1, l)
    enddo
    qpolen(l) = sommemqn(l) / sommemassen(l)
  enddo

  ! On force le traceur à prendre cette valeur aux pôles
  do l = 1, llm
    do i = 1, iip1
      q(i, 1, l) = qpolen(l)
      q(i, jjp1, l) = qpoles(l)
    enddo
  enddo

  RETURN
END SUBROUTINE traceurpole