source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_bilan_dyn.f90

MODULE lmdz_bilan_dyn
  IMPLICIT NONE; PRIVATE
  PUBLIC bilan_dyn

CONTAINS


  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
    USE lmdz_comgeom2

    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    USE lmdz_paramet
    IMPLICIT NONE




    !====================================================================

    !   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

END MODULE lmdz_bilan_dyn