! $Id: bilan_dyn.F90 5134 2024-07-26 15:56:37Z abarral $

SUBROUTINE bilan_dyn(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 comconst_mod, ONLY: pi, cpp
  USE comvert_mod, ONLY: presnivs
  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
  USE lmdz_iniprint, ONLY: lunout, prt_level


  IMPLICIT NONE

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  INCLUDE "comgeom2.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, jjp1)
  REAL :: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
  REAL :: flux_u(iip1, jjp1, llm)
  REAL :: flux_v(iip1, jjm, llm)
  REAL :: teta(iip1, jjp1, llm)
  REAL :: phi(iip1, jjp1, llm)
  REAL :: ucov(iip1, jjp1, llm)
  REAL :: vcov(iip1, jjm, llm)
  REAL :: trac(iip1, jjp1, llm, ntrac)

  !   Local :
  !   =======

  INTEGER :: icum, ncum
  LOGICAL :: first
  REAL :: zz, zqy, zfactv(jjm, llm)

  INTEGER :: nQ
  parameter (nQ = 7)


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

  CHARACTER(LEN = 10) :: file
  INTEGER :: ifile
  parameter (ifile = 4)

  INTEGER :: itemp, igeop, iecin, iang, iu, iovap, iun
  INTEGER :: i_sortie

  save first, icum, ncum
  save itemp, igeop, iecin, iang, iu, iovap, iun
  save i_sortie

  REAL :: time
  INTEGER :: itau
  save time, itau
  data time, itau/0., 0/

  data first/.TRUE./
  data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/
  data i_sortie/1/

  REAL :: ww

  !   variables dynamiques intermédiaires
  REAL :: vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
  REAL :: ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
  REAL :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
  REAL :: vorpot(iip1, jjm, llm)
  REAL :: w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
  REAL :: bern(iip1, jjp1, llm)

  !   champ contenant les scalaires advectés.
  REAL :: Q(iip1, jjp1, llm, nQ)

  !   champs cumulés
  REAL :: ps_cum(iip1, jjp1)
  REAL :: masse_cum(iip1, jjp1, llm)
  REAL :: flux_u_cum(iip1, jjp1, llm)
  REAL :: flux_v_cum(iip1, jjm, llm)
  REAL :: Q_cum(iip1, jjp1, llm, nQ)
  REAL :: flux_uQ_cum(iip1, jjp1, llm, nQ)
  REAL :: flux_vQ_cum(iip1, jjm, llm, nQ)
  REAL :: flux_wQ_cum(iip1, jjp1, llm, nQ)
  REAL :: dQ(iip1, jjp1, llm, nQ)

  save ps_cum, masse_cum, flux_u_cum, flux_v_cum
  save Q_cum, flux_uQ_cum, flux_vQ_cum

  !   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 :: iave, itot, immc, itrs, istn
  data iave, itot, immc, itrs, istn/1, 2, 3, 4, 5/
  CHARACTER(LEN = 3) :: ctrs(ntr)
  data ctrs/'  ', 'TOT', 'MMC', 'TRS', 'STN'/

  REAL :: zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
  REAL :: zavQ(jjm, ntr, nQ), psiQ(jjm, llm + 1, nQ)
  REAL :: zmasse(jjm, llm), zamasse(jjm)

  REAL :: zv(jjm, llm), psi(jjm, llm + 1)

  INTEGER :: i, j, l, iQ


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

  CHARACTER(LEN = 10) :: infile

  INTEGER :: fileid
  INTEGER :: thoriid, zvertiid
  save fileid

  INTEGER :: ndex3d(jjm * llm)

  !   Variables locales
  !
  INTEGER :: tau0
  REAL :: zjulian
  CHARACTER(LEN = 3) :: str
  CHARACTER(LEN = 10) :: ctrac
  INTEGER :: ii, jj
  INTEGER :: zan, dayref
  !
  REAL :: rlong(jjm), rlatg(jjm)



  !=====================================================================
  !   Initialisation
  !=====================================================================

  time = time + dt_app
  itau = itau + 1
  !IM
  ndex3d = 0

  IF (first) THEN
    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('bilan_dyn', 'stopped', 1)
    endif

    IF (i_sortie==1) THEN
      file = 'dynzon'
      CALL inigrads(ifile, 1 &
              , 0., 180. / pi, 0., 0., jjm, rlatv, -90., 90., 180. / pi &
              , llm, presnivs, 1. &
              , dt_cum, file, 'dyn_zon ')
    endif

    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

    CALL histbeg(infile, 1, rlong, jjm, rlatg, &
            1, 1, 1, jjm, &
            tau0, zjulian, dt_cum, thoriid, fileid)

    !
    !  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, jjm, 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, jjm, 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, jjm, thoriid, llm, 1, llm, zvertiid, &
            32, 'ave(X)', dt_cum, dt_cum)
    CALL histdef(fileid, 'v', 'v', &
            'm/s', 1, jjm, 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, jjm, 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, jjm, thoriid, 1, 1, 1, -99, &
                32, 'ave(X)', dt_cum, dt_cum)
      enddo
    enddo


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

  ENDIF


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

  !   énergie cinétique
  ucont(:, :, :) = 0
  CALL covcont(llm, ucov, vcov, ucont, vcont)
  CALL enercin(vcov, ucov, vcont, ucont, ecin)

  !   moment cinétique
  do l = 1, llm
    ang(:, :, l) = ucov(:, :, l) + constang(:, :)
    unat(:, :, l) = ucont(:, :, l) * cu(:, :)
  enddo

  Q(:, :, :, itemp) = teta(:, :, :) * pk(:, :, :) / cpp
  Q(:, :, :, igeop) = phi(:, :, :)
  Q(:, :, :, iecin) = ecin(:, :, :)
  Q(:, :, :, iang) = ang(:, :, :)
  Q(:, :, :, iu) = unat(:, :, :)
  Q(:, :, :, iovap) = trac(:, :, :, 1)
  Q(:, :, :, iun) = 1.


  !=====================================================================
  !   Cumul
  !=====================================================================
  !
  IF(icum==0) THEN
    ps_cum = 0.
    masse_cum = 0.
    flux_u_cum = 0.
    flux_v_cum = 0.
    Q_cum = 0.
    flux_vQ_cum = 0.
    flux_uQ_cum = 0.
  ENDIF

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

  !   accumulation des flux de masse horizontaux
  ps_cum = ps_cum + ps
  masse_cum = masse_cum + masse
  flux_u_cum = flux_u_cum + flux_u
  flux_v_cum = flux_v_cum + flux_v
  do iQ = 1, nQ
    Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ) * masse(:, :, :)
  enddo

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

  !   Flux longitudinal
  !   -----------------
  do iQ = 1, nQ
    do l = 1, llm
      do j = 1, jjp1
        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
  enddo

  !    flux méridien
  !    -------------
  do iQ = 1, nQ
    do l = 1, llm
      do j = 1, jjm
        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
  enddo


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

  !   convergence horizontale
  CALL  convflu(flux_uQ_cum, flux_vQ_cum, llm * nQ, dQ)

  !   calcul de la vitesse verticale
  CALL convmas(flux_u_cum, flux_v_cum, convm)
  CALL vitvert(convm, w)

  do iQ = 1, nQ
    do l = 1, llm - 1
      do j = 1, jjp1
        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
  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'

    !   Normalisation
    do iQ = 1, nQ
      Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum(:, :, :)
    enddo
    zz = 1. / REAL(ncum)
    ps_cum = ps_cum * zz
    masse_cum = masse_cum * zz
    flux_u_cum = flux_u_cum * zz
    flux_v_cum = flux_v_cum * zz
    flux_uQ_cum = flux_uQ_cum * zz
    flux_vQ_cum = flux_vQ_cum * zz
    dQ = dQ * zz


    !   A retravailler eventuellement
    !   division de dQ par la masse pour revenir aux bonnes grandeurs
    do iQ = 1, nQ
      dQ(:, :, :, iQ) = dQ(:, :, :, iQ) / masse_cum(:, :, :)
    enddo

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

    !   cumul zonal des masses des mailles
    !   ----------------------------------
    zv = 0.
    zmasse = 0.
    CALL massbar(masse_cum, massebx, masseby)
    do l = 1, llm
      do j = 1, jjm
        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

    ! 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
    !   ----------------------------------------

    zvQ = 0.
    psiQ = 0.
    do iQ = 1, nQ
      zvQtmp = 0.
      do l = 1, llm
        do j = 1, jjm
          ! 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
      !   fonction de courant meridienne pour la quantite Q
      do l = llm, 1, -1
        do j = 1, jjm
          psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
        enddo
      enddo
    enddo

    !   fonction de courant pour la circulation meridienne moyenne
    psi = 0.
    do l = llm, 1, -1
      do j = 1, jjm
        psi(j, l) = psi(j, l + 1) + zv(j, l)
        zv(j, l) = zv(j, l) * zfactv(j, l)
      enddo
    enddo

    ! PRINT*,'4OK'
    !   sorties proprement dites
    IF (i_sortie==1) THEN
      do iQ = 1, nQ
        do itr = 1, ntr
          CALL histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ) &
                  , jjm * llm, ndex3d)
        enddo
        CALL histwrite(fileid, 'psi' // nom(iQ), itau, psiQ(:, 1:llm, iQ) &
                , jjm * llm, ndex3d)
      enddo

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

    endif


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

    zamasse = 0.
    do l = 1, llm
      zamasse(:) = zamasse(:) + zmasse(:, l)
    enddo
    zavQ = 0.
    do iQ = 1, nQ
      do itr = 2, ntr
        do l = 1, llm
          zavQ(:, itr, iQ) = zavQ(:, itr, iQ) + zvQ(:, l, itr, iQ) * zmasse(:, l)
        enddo
        zavQ(:, itr, iQ) = zavQ(:, itr, iQ) / zamasse(:)
        CALL histwrite(fileid, 'a' // znom(itr, iQ), itau, zavQ(:, itr, iQ) &
                , jjm * llm, ndex3d)
      enddo
    enddo

    ! on doit pouvoir tracer systematiquement la fonction de courant.

    !=====================================================================
    !/////////////////////////////////////////////////////////////////////
    icum = 0                  !///////////////////////////////////////
  ENDIF ! icum.EQ.ncum    !///////////////////////////////////////
  !/////////////////////////////////////////////////////////////////////
  !=====================================================================


END SUBROUTINE  bilan_dyn