! $Id: bilan_dyn_p.F 1299 2010-01-20 14:27:21Z fairhead $

SUBROUTINE bilan_dyn_loc (ntrac,dt_app,dt_cum, &
        ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)

  !   AFAIRE
  !   Prevoir en champ nq+1 le diagnostique de l'energie
  !   en faisant Qzon=Cv T + L * ...
  !             vQ..A=Cp T + L * ...

  USE IOIPSL
  USE parallel_lmdz
  USE mod_hallo
  use misc_mod
  USE write_field_loc
  USE comconst_mod, ONLY: cpp, pi
  USE comvert_mod, ONLY: presnivs
  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "comgeom2.h"
  include "iniprint.h"

  !====================================================================
  !
  !   Sous-programme consacre à des diagnostics dynamiques de base
  !
  !
  !   De facon generale, les moyennes des scalaires Q sont ponderees par
  !   la masse.
  !
  !   Les flux de masse sont eux simplement moyennes.
  !
  !====================================================================

  !   Arguments :
  !   ===========

  integer :: ntrac
  real :: dt_app,dt_cum
  real :: ps(iip1,jjb_u:jje_u)
  real :: masse(iip1,jjb_u:jje_u,llm),pk(iip1,jjb_u:jje_u,llm)
  real :: flux_u(iip1,jjb_u:jje_u,llm)
  real :: flux_v(iip1,jjb_v:jje_v,llm)
  real :: teta(iip1,jjb_u:jje_u,llm)
  real :: phi(iip1,jjb_u:jje_u,llm)
  real :: ucov(iip1,jjb_u:jje_u,llm)
  real :: vcov(iip1,jjb_v:jje_v,llm)
  real :: trac(iip1,jjb_u:jje_u,llm,ntrac)

  !   Local :
  !   =======

  integer,SAVE :: icum,ncum
!$OMP THREADPRIVATE(icum,ncum)
  LOGICAL,SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)

  real :: zz,zqy
  REAl,SAVE,ALLOCATABLE :: zfactv(:,:)

  INTEGER,PARAMETER :: nQ=7


  !ym      character*6 nom(nQ)
  !ym      character*6 unites(nQ)
  character(len=6),save :: nom(nQ)
  character(len=6),save :: unites(nQ)

  character(len=10) file
  integer :: ifile
  parameter (ifile=4)

  integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
  INTEGER,PARAMETER :: iovap=6,iun=7
  integer,PARAMETER :: i_sortie=1

  real,SAVE :: time=0.
  integer,SAVE :: itau=0.
!$OMP THREADPRIVATE(time,itau)

  real :: ww

  !   variables dynamiques intermédiaires
  REAL,SAVE,ALLOCATABLE :: vcont(:,:,:),ucont(:,:,:)
  REAL,SAVE,ALLOCATABLE :: ang(:,:,:),unat(:,:,:)
  REAL,SAVE,ALLOCATABLE :: massebx(:,:,:),masseby(:,:,:)
  REAL,SAVE,ALLOCATABLE :: vorpot(:,:,:)
  REAL,SAVE,ALLOCATABLE :: w(:,:,:),ecin(:,:,:),convm(:,:,:)
  REAL,SAVE,ALLOCATABLE :: bern(:,:,:)

  !   champ contenant les scalaires advectés.
  real,SAVE,ALLOCATABLE :: Q(:,:,:,:)

  !   champs cumulés
  real,SAVE,ALLOCATABLE ::  ps_cum(:,:)
  real,SAVE,ALLOCATABLE ::  masse_cum(:,:,:)
  real,SAVE,ALLOCATABLE ::  flux_u_cum(:,:,:)
  real,SAVE,ALLOCATABLE ::  flux_v_cum(:,:,:)
  real,SAVE,ALLOCATABLE ::  Q_cum(:,:,:,:)
  real,SAVE,ALLOCATABLE ::  flux_uQ_cum(:,:,:,:)
  real,SAVE,ALLOCATABLE ::  flux_vQ_cum(:,:,:,:)
  real,SAVE,ALLOCATABLE ::  flux_wQ_cum(:,:,:,:)
  real,SAVE,ALLOCATABLE ::  dQ(:,:,:,:)


  !   champs de tansport en moyenne zonale
  integer :: ntr,itr
  parameter (ntr=5)

  !ym      character*10 znom(ntr,nQ)
  !ym      character*20 znoml(ntr,nQ)
  !ym      character*10 zunites(ntr,nQ)
  character*10,save :: znom(ntr,nQ)
  character*20,save :: znoml(ntr,nQ)
  character*10,save :: zunites(ntr,nQ)

  INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5

  character(len=3) :: ctrs(ntr)
  data ctrs/'  ','TOT','MMC','TRS','STN'/

  real,SAVE,ALLOCATABLE ::  zvQ(:,:,:,:),zvQtmp(:,:)
  real,SAVE,ALLOCATABLE ::  zavQ(:,:,:),psiQ(:,:,:)
  real,SAVE,ALLOCATABLE ::  zmasse(:,:),zamasse(:)

  real,SAVE,ALLOCATABLE ::  zv(:,:),psi(:,:)

  integer :: i,j,l,iQ


  !   Initialisation du fichier contenant les moyennes zonales.
  !   ---------------------------------------------------------

  character(len=10) :: infile

  integer, save :: fileid
  integer :: thoriid, zvertiid

  INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)

  !   Variables locales
  !
  integer :: tau0
  real :: zjulian
  character(len=3) :: str
  character(len=10) :: ctrac
  integer :: ii,jj
  integer :: zan, dayref
  !
  real,SAVE,ALLOCATABLE :: rlong(:),rlatg(:)
  integer :: jjb,jje,jjn,ijb,ije
  type(Request),SAVE :: Req
!$OMP THREADPRIVATE(Req)

  ! definition du domaine d'ecriture pour le rebuild

  INTEGER,DIMENSION(1) :: ddid
  INTEGER,DIMENSION(1) :: dsg
  INTEGER,DIMENSION(1) :: dsl
  INTEGER,DIMENSION(1) :: dpf
  INTEGER,DIMENSION(1) :: dpl
  INTEGER,DIMENSION(1) :: dhs
  INTEGER,DIMENSION(1) :: dhe

  INTEGER :: bilan_dyn_domain_id

  !=====================================================================
  !   Initialisation
  !=====================================================================
  if (adjust) return

  time=time+dt_app
  itau=itau+1

  if (first) then
!$OMP BARRIER
!$OMP MASTER
  ALLOCATE(zfactv(jjb_v:jje_v,llm))
  ALLOCATE(vcont(iip1,jjb_v:jje_v,llm))
  ALLOCATE(ucont(iip1,jjb_u:jje_u,llm))
  ALLOCATE(ang(iip1,jjb_u:jje_u,llm))
  ALLOCATE(unat(iip1,jjb_u:jje_u,llm))
  ALLOCATE(massebx(iip1,jjb_u:jje_u,llm))
  ALLOCATE(masseby(iip1,jjb_v:jje_v,llm))
  ALLOCATE(vorpot(iip1,jjb_v:jje_v,llm))
  ALLOCATE(w(iip1,jjb_u:jje_u,llm))
  ALLOCATE(ecin(iip1,jjb_u:jje_u,llm))
  ALLOCATE(convm(iip1,jjb_u:jje_u,llm))
  ALLOCATE(bern(iip1,jjb_u:jje_u,llm))
  ALLOCATE(Q(iip1,jjb_u:jje_u,llm,nQ))
  ALLOCATE(ps_cum(iip1,jjb_u:jje_u))
  ALLOCATE(masse_cum(iip1,jjb_u:jje_u,llm))
  ALLOCATE(flux_u_cum(iip1,jjb_u:jje_u,llm))
  ALLOCATE(flux_v_cum(iip1,jjb_v:jje_v,llm))
  ALLOCATE(Q_cum(iip1,jjb_u:jje_u,llm,nQ))
  ALLOCATE(flux_uQ_cum(iip1,jjb_u:jje_u,llm,nQ))
  ALLOCATE(flux_vQ_cum(iip1,jjb_v:jje_v,llm,nQ))
  ALLOCATE(flux_wQ_cum(iip1,jjb_u:jje_u,llm,nQ))
  ALLOCATE(dQ(iip1,jjb_u:jje_u,llm,nQ))
  ALLOCATE(zvQ(jjb_v:jje_v,llm,ntr,nQ))
  ALLOCATE(zvQtmp(jjb_v:jje_v,llm))
  ALLOCATE(zavQ(jjb_v:jje_v,ntr,nQ))
  ALLOCATE(psiQ(jjb_v:jje_v,llm+1,nQ))
  ALLOCATE(zmasse(jjb_v:jje_v,llm))
  ALLOCATE(zamasse(jjb_v:jje_v))
  ALLOCATE(zv(jjb_v:jje_v,llm))
  ALLOCATE(psi(jjb_v:jje_v,llm+1))
  ALLOCATE(ndex3d(jjb_v:jje_v*llm))
  ndex3d=0
  ALLOCATE(rlong(1))
  ALLOCATE(rlatg(jjm))

!$OMP END MASTER
!$OMP BARRIER
    icum=0
    ! initialisation des fichiers
    first=.FALSE.
  !   ncum est la frequence de stokage en pas de temps
    ncum=dt_cum/dt_app
    if (abs(ncum*dt_app-dt_cum)>1.e-5*dt_app) then
       WRITE(lunout,*) &
             'Pb : le pas de cumule doit etre multiple du pas'
       WRITE(lunout,*)'dt_app=',dt_app
       WRITE(lunout,*)'dt_cum=',dt_cum
       CALL abort_gcm("conf_gcmbilan_dyn_loc","stopped",1)
    endif

!$OMP MASTER
    nom(itemp)='T'
    nom(igeop)='gz'
    nom(iecin)='K'
    nom(iang)='ang'
    nom(iu)='u'
    nom(iovap)='ovap'
    nom(iun)='un'

    unites(itemp)='K'
    unites(igeop)='m2/s2'
    unites(iecin)='m2/s2'
    unites(iang)='ang'
    unites(iu)='m/s'
    unites(iovap)='kg/kg'
    unites(iun)='un'


  !   Initialisation du fichier contenant les moyennes zonales.
  !   ---------------------------------------------------------

  infile='dynzon'

  zan = annee_ref
  dayref = day_ref
  CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
  tau0 = itau_dyn

  rlong=0.
  rlatg=rlatv*180./pi

  jjb=jj_begin
  jje=jj_end
  jjn=jj_nb
  IF (pole_sud) THEN
    jjn=jj_nb-1
    jje=jj_end-1
  ENDIF

  ddid=(/ 2 /)
  dsg=(/ jjm /)
  dsl=(/ jjn /)
  dpf=(/ jjb /)
  dpl=(/ jje /)
  dhs=(/ 0 /)
  dhe=(/ 0 /)

  CALL flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, &
        'box',bilan_dyn_domain_id)

  CALL histbeg(trim(infile), &
        1, rlong, jjn, rlatg(jjb:jje), &
        1, 1, 1, jjn, &
        tau0, zjulian, dt_cum, thoriid, fileid, &
        bilan_dyn_domain_id)

  !
  !  Appel a histvert pour la grille verticale
  !
  CALL histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
        llm, presnivs, zvertiid)
  !
  !  Appels a histdef pour la definition des variables a sauvegarder
  do iQ=1,nQ
     do itr=1,ntr
        if(itr==1) then
           znom(itr,iQ)=nom(iQ)
           znoml(itr,iQ)=nom(iQ)
           zunites(itr,iQ)=unites(iQ)
        else
           znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
           znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
           zunites(itr,iQ)='m/s * '//unites(iQ)
        endif
     enddo
  enddo

  !   Declarations des champs avec dimension verticale
   ! PRINT*,'1HISTDEF'
  do iQ=1,nQ
     do itr=1,ntr
  IF (prt_level > 5) &
        WRITE(lunout,*)'var ',itr,iQ &
        ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
        CALL histdef(fileid,znom(itr,iQ),znoml(itr,iQ), &
              zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid, &
              32,'ave(X)',dt_cum,dt_cum)
     enddo
  !   Declarations pour les fonctions de courant
   ! PRINT*,'2HISTDEF'
      CALL histdef(fileid,'psi'//nom(iQ) &
            ,'stream fn. '//znoml(itot,iQ), &
            zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid, &
            32,'ave(X)',dt_cum,dt_cum)
  enddo


  !   Declarations pour les champs de transport d'air
   ! PRINT*,'3HISTDEF'
  CALL histdef(fileid, 'masse', 'masse', &
        'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid, &
        32, 'ave(X)', dt_cum, dt_cum)
  CALL histdef(fileid, 'v', 'v', &
        'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid, &
        32, 'ave(X)', dt_cum, dt_cum)
  !   Declarations pour les fonctions de courant
   ! PRINT*,'4HISTDEF'
      CALL histdef(fileid,'psi','stream fn. MMC ','mega t/s', &
            1,jjn,thoriid,llm,1,llm,zvertiid, &
            32,'ave(X)',dt_cum,dt_cum)


  !   Declaration des champs 1D de transport en latitude
   ! PRINT*,'5HISTDEF'
  do iQ=1,nQ
     do itr=2,ntr
        CALL histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ), &
              zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99, &
              32,'ave(X)',dt_cum,dt_cum)
     enddo
  enddo


   ! PRINT*,'8HISTDEF'
           CALL histend(fileid)

!$OMP END MASTER
  endif


  !=====================================================================
  !   Calcul des champs dynamiques
  !   ----------------------------

  jjb=jj_begin
  jje=jj_end

  !   énergie cinétique
   ! ucont(:,jjb:jje,:)=0

  CALL Register_Hallo_u(ucov,llm,1,1,1,1,Req)
  CALL Register_Hallo_v(vcov,llm,1,1,1,1,Req)
  CALL SendRequest(Req)
!$OMP BARRIER
  CALL WaitRequest(Req)

  CALL covcont_loc(llm,ucov,vcov,ucont,vcont)
  CALL enercin_loc(vcov,ucov,vcont,ucont,ecin)

  !   moment cinétique
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  do l=1,llm
     ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
     unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
  enddo
!$OMP END DO NOWAIT

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    Q(:,jjb:jje,l,itemp)=teta(:,jjb:jje,l)*pk(:,jjb:jje,l)/cpp
    Q(:,jjb:jje,l,igeop)=phi(:,jjb:jje,l)
    Q(:,jjb:jje,l,iecin)=ecin(:,jjb:jje,l)
    Q(:,jjb:jje,l,iang)=ang(:,jjb:jje,l)
    Q(:,jjb:jje,l,iu)=unat(:,jjb:jje,l)
    Q(:,jjb:jje,l,iovap)=trac(:,jjb:jje,l,1)
    Q(:,jjb:jje,l,iun)=1.
  ENDDO
!$OMP END DO NOWAIT

  !=====================================================================
  !   Cumul
  !=====================================================================
  !
  if(icum==0) then
     jjb=jj_begin
     jje=jj_end

!$OMP MASTER
     ps_cum(:,jjb:jje)=0.
!$OMP END MASTER


!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l=1,llm
      masse_cum(:,jjb:jje,l)=0.
      flux_u_cum(:,jjb:jje,l)=0.
      Q_cum(:,jjb:jje,l,:)=0.
      flux_uQ_cum(:,jjb:jje,l,:)=0.
      if (pole_sud) jje=jj_end-1
      flux_v_cum(:,jjb:jje,l)=0.
      flux_vQ_cum(:,jjb:jje,l,:)=0.
    ENDDO
!$OMP END DO NOWAIT
  endif

  IF (prt_level > 5) &
        WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
  icum=icum+1

  !   accumulation des flux de masse horizontaux
  jjb=jj_begin
  jje=jj_end

!$OMP MASTER
  ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
!$OMP END MASTER


!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)+masse(:,jjb:jje,l)
    flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l) &
          +flux_u(:,jjb:jje,l)
  ENDDO
!$OMP END DO NOWAIT

  if (pole_sud) jje=jj_end-1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
   flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l) &
         +flux_v(:,jjb:jje,l)
  ENDDO
!$OMP END DO NOWAIT

  jjb=jj_begin
  jje=jj_end

  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l=1,llm
      Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ) &
            +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
    ENDDO
!$OMP END DO NOWAIT
  enddo

  !=====================================================================
  !  FLUX ET TENDANCES
  !=====================================================================

  !   Flux longitudinal
  !   -----------------
  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     do l=1,llm
        do j=jjb,jje
           do i=1,iim
              flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ) &
                    +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
           enddo
           flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
        enddo
     enddo
!$OMP END DO NOWAIT
  enddo

  !    flux méridien
  !    -------------
  do iQ=1,nQ
    CALL Register_Hallo_u(Q(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
  enddo
  CALL SendRequest(Req)
!$OMP BARRIER
  CALL WaitRequest(Req)

  jjb=jj_begin
  jje=jj_end
  if (pole_sud) jje=jj_end-1

  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     do l=1,llm
        do j=jjb,jje
           do i=1,iip1
              flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ) &
                    +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
           enddo
        enddo
     enddo
!$OMP ENDDO NOWAIT
!$OMP BARRIER
  enddo

  !    tendances
  !    ---------

  !   convergence horizontale
  CALL Register_Hallo_u(flux_uQ_cum,llm,2,2,2,2,Req)
  CALL Register_Hallo_v(flux_vQ_cum,llm,2,2,2,2,Req)
  CALL SendRequest(Req)
!$OMP BARRIER
  CALL WaitRequest(Req)

  CALL  convflu_loc(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)

  !   calcul de la vitesse verticale
  CALL Register_Hallo_u(flux_u_cum,llm,2,2,2,2,Req)
  CALL Register_Hallo_v(flux_v_cum,llm,2,2,2,2,Req)
  CALL SendRequest(Req)
!$OMP BARRIER
  CALL WaitRequest(Req)

  CALL convmas_loc(flux_u_cum,flux_v_cum,convm)
  CALL vitvert_loc(convm,w)
!$OMP BARRIER


  jjb=jj_begin
  jje=jj_end

   ! do iQ=1,nQ
   !    do l=1,llm-1
   !       do j=jjb,jje
   !          do i=1,iip1
   !             ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
   !             dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
   !             dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
   !          enddo
   !       enddo
   !     enddo
   !  enddo

  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     do l=1,llm
        IF (l<llm) THEN
          do j=jjb,jje
             do i=1,iip1
                ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
                dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
                dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
             enddo
          enddo
        ENDIF
        IF (l>2) THEN
          do j=jjb,jje
            do i=1,iip1
              ww=-0.5*w(i,j,l)*(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
              dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
            enddo
          enddo
        ENDIF
     enddo
!$OMP ENDDO NOWAIT
  enddo
  IF (prt_level > 5) &
        WRITE(lunout,*)'Apres les calculs fait a chaque pas'
  !=====================================================================
  !   PAS DE TEMPS D'ECRITURE
  !=====================================================================
  if (icum==ncum) then
  !=====================================================================

  IF (prt_level > 5) &
        WRITE(lunout,*)'Pas d ecriture'

  jjb=jj_begin
  jje=jj_end

  !   Normalisation
  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l=1,llm
      Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ) &
            /masse_cum(:,jjb:jje,l)
    enddo
!$OMP ENDDO NOWAIT
  enddo

  zz=1./REAL(ncum)

!$OMP MASTER
    ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
!$OMP END MASTER

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
    flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
    flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
    dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
  ENDDO
!$OMP ENDDO NOWAIT

  IF (pole_sud) jje=jj_end-1
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
    flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)*zz
    flux_vQ_cum(:,jjb:jje,l,:)=flux_vQ_cum(:,jjb:jje,l,:)*zz
  ENDDO
!$OMP ENDDO NOWAIT
!$OMP BARRIER

  jjb=jj_begin
  jje=jj_end


  !   A retravailler eventuellement
  !   division de dQ par la masse pour revenir aux bonnes grandeurs
  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l=1,llm
       dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
    ENDDO
!$OMP ENDDO NOWAIT
  enddo

  !=====================================================================
  !   Transport méridien
  !=====================================================================

  !   cumul zonal des masses des mailles
  !   ----------------------------------
  jjb=jj_begin
  jje=jj_end
  if (pole_sud) jje=jj_end-1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l=1,llm
      zv(jjb:jje,l)=0.
      zmasse(jjb:jje,l)=0.
    ENDDO
!$OMP ENDDO NOWAIT
!$OMP BARRIER

  CALL Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
  do iQ=1,nQ
    CALL Register_Hallo_u(Q_cum(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
  enddo

  CALL SendRequest(Req)
!$OMP BARRIER
  CALL WaitRequest(Req)

  CALL massbar_loc(masse_cum,massebx,masseby)

  jjb=jj_begin
  jje=jj_end
  if (pole_sud) jje=jj_end-1

!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  do l=1,llm
     do j=jjb,jje
        do i=1,iim
           zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
           zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
        enddo
        zfactv(j,l)=cv(1,j)/zmasse(j,l)
     enddo
  enddo
!$OMP ENDDO NOWAIT
!$OMP BARRIER

  ! PRINT*,'3OK'
  !   --------------------------------------------------------------
  !   calcul de la moyenne zonale du transport :
  !   ------------------------------------------
  !
  !                                 --
  ! TOT : la circulation totale       [ vq ]
  !
  !                                  -     -
  ! MMC : mean meridional circulation [ v ] [ q ]
  !
  !                                 ----      --       - -
  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
  !
  !                                 - * - *       - -       -     -
  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
  !
  !                                          - -
  !    on utilise aussi l'intermediaire TMP :  [ v q ]
  !
  !    la variable zfactv transforme un transport meridien cumule
  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
  !
  !   --------------------------------------------------------------


  !   ----------------------------------------
  !   Transport dans le plan latitude-altitude
  !   ----------------------------------------

  jjb=jj_begin
  jje=jj_end
  if (pole_sud) jje=jj_end-1

  zvQ=0.
  psiQ=0.
  do iQ=1,nQ
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     do l=1,llm
        zvQtmp(:,l)=0.
        do j=jjb,jje
           ! PRINT*,'j,l,iQ=',j,l,iQ
  !   Calcul des moyennes zonales du transort total et de zvQtmp
           do i=1,iim
              zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ) &
                    +flux_vQ_cum(i,j,l,iQ)
              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
                    Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
              zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy &
                    /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
              zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
           enddo
           ! PRINT*,'aOK'
  !   Decomposition
           zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
           zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
           zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
           zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
           zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
           zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
        enddo
     enddo
!$OMP ENDDO NOWAIT
  !   fonction de courant meridienne pour la quantite Q
!$OMP BARRIER
!$OMP MASTER
     do l=llm,1,-1
        do j=jjb,jje
           psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
        enddo
     enddo
!$OMP END MASTER
!$OMP BARRIER
  enddo

  !   fonction de courant pour la circulation meridienne moyenne
!$OMP BARRIER
!$OMP MASTER
  psi(jjb:jje,:)=0.
  do l=llm,1,-1
     do j=jjb,jje
        psi(j,l)=psi(j,l+1)+zv(j,l)
        zv(j,l)=zv(j,l)*zfactv(j,l)
     enddo
  enddo
!$OMP END MASTER
!$OMP BARRIER

  ! PRINT*,'4OK'
  !   sorties proprement dites
!$OMP MASTER
  if (i_sortie==1) then
  jjb=jj_begin
  jje=jj_end
  jjn=jj_nb
  if (pole_sud) jje=jj_end-1
  if (pole_sud) jjn=jj_nb-1
  do iQ=1,nQ
     do itr=1,ntr
        CALL histwrite(fileid,znom(itr,iQ),itau, &
              zvQ(jjb:jje,:,itr,iQ) &
              ,jjn*llm,ndex3d)
     enddo
     CALL histwrite(fileid,'psi'//nom(iQ), &
           itau,psiQ(jjb:jje,1:llm,iQ) &
           ,jjn*llm,ndex3d)
  enddo

  CALL histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm) &
        ,jjn*llm,ndex3d)
  CALL histwrite(fileid,'v',itau,zv(jjb:jje,1:llm) &
        ,jjn*llm,ndex3d)
  psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
  CALL histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm), &
        jjn*llm,ndex3d)

  endif


  !   -----------------
  !   Moyenne verticale
  !   -----------------

  zamasse(jjb:jje)=0.
  do l=1,llm
     zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
  enddo

  zavQ(jjb:jje,:,:)=0.
  do iQ=1,nQ
     do itr=2,ntr
        do l=1,llm
           zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ) &
                 +zvQ(jjb:jje,l,itr,iQ) &
                 *zmasse(jjb:jje,l)
        enddo
        zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
        CALL histwrite(fileid,'a'//znom(itr,iQ),itau, &
              zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
     enddo
  enddo
!$OMP END MASTER
  ! on doit pouvoir tracer systematiquement la fonction de courant.

  !=====================================================================
  !/////////////////////////////////////////////////////////////////////
  icum=0                  !///////////////////////////////////////
  endif ! icum.eq.ncum    !///////////////////////////////////////
  !/////////////////////////////////////////////////////////////////////
  !=====================================================================
!$OMP MASTER
  CALL histsync(fileid)
!$OMP END MASTER



end subroutine bilan_dyn_loc