source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlsplt_loc.f90 @ 5133

! $Id: vlsplt_loc.f90 5128 2024-07-25 15:47:25Z abarral $

SUBROUTINE vlx_loc(q, pente_max, masse, u_m, ijb_x, ije_x, iq)

  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
  !
  !    ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  !    ********************************************************************
  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
  !
  !
  !   --------------------------------------------------------------------
  USE parallel_lmdz
  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
          min_qParent, min_qMass, min_ratio ! MVals et CRisi
  USE lmdz_iniprint, ONLY: lunout, prt_level
  IMPLICIT NONE
  !
  include "dimensions.h"
  include "paramet.h"
  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
  REAL :: u_m(ijb_u:ije_u, llm), pbarv(iip1, jjb_v:jje_v, llm)
  REAL :: q(ijb_u:ije_u, llm, nqtot) ! CRisi: ajout dimension nqtot
  REAL :: w(ijb_u:ije_u, llm)
  INTEGER :: iq ! CRisi
  !
  !  Local
  !   ---------
  !
  INTEGER :: ij, l, j, i, iju, ijq, indu(ijnb_u), niju
  INTEGER :: n0, iadvplus(ijb_u:ije_u, llm), nl(llm)
  !
  REAL :: new_m, zu_m, zdum(ijb_u:ije_u, llm)
  REAL :: sigu(ijb_u:ije_u), dxq(ijb_u:ije_u, llm), dxqu(ijb_u:ije_u)
  REAL :: zz(ijb_u:ije_u)
  REAL :: adxqu(ijb_u:ije_u), dxqmax(ijb_u:ije_u, llm)
  REAL :: u_mq(ijb_u:ije_u, llm)

  REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi
  INTEGER :: ifils, iq2 ! CRisi

  Logical :: extremum

  REAL :: z1, z2, z3

  INTEGER :: ijb, ije, ijb_x, ije_x

  !WRITE(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
  ! &   iq,ijb_x
  !   calcul de la pente a droite et a gauche de la maille

  ijb = ijb_x
  ije = ije_x

  IF (pole_nord.AND.ijb==1) ijb = ijb + iip1
  IF (pole_sud.AND.ije==ip1jmp1)  ije = ije - iip1

  IF (pente_max>-1.e-5) THEN
    ! IF (pente_max.gt.10) THEN

    !   calcul des pentes avec limitation, Van Leer scheme I:
    !   -----------------------------------------------------
    ! on a besoin de q entre ijb et ije
    !   calcul de la pente aux points u
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm

      DO ij = ijb, ije - 1
        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
        ! IF(u_m(ij,l).lt.0.) stop 'limx n admet pas les U<0'
        ! sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
      ENDDO
      DO ij = ijb + iip1 - 1, ije, iip1
        dxqu(ij) = dxqu(ij - iim)
        ! sigu(ij)=sigu(ij-iim)
      ENDDO

      DO ij = ijb, ije
        adxqu(ij) = abs(dxqu(ij))
      ENDDO

      !   calcul de la pente maximum dans la maille en valeur absolue

      DO ij = ijb + 1, ije
        dxqmax(ij, l) = pente_max * &
                min(adxqu(ij - 1), adxqu(ij))
        ! limitation subtile
        !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))

      ENDDO

      DO ij = ijb + iip1 - 1, ije, iip1
        dxqmax(ij - iim, l) = dxqmax(ij, l)
      ENDDO

      DO ij = ijb + 1, ije
        IF(dxqu(ij - 1) * dxqu(ij)>0) THEN
          dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
        ELSE
          !   extremum local
          dxq(ij, l) = 0.
        ENDIF
        dxq(ij, l) = 0.5 * dxq(ij, l)
        dxq(ij, l) = &
                sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
      ENDDO

    ENDDO ! l=1,llm
    !$OMP END DO NOWAIT
    ! PRINT*,'Ok calcul des pentes'

  ELSE ! (pente_max.lt.-1.e-5)

    !   Pentes produits:
    !   ----------------
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije - 1
        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
      ENDDO
      DO ij = ijb + iip1 - 1, ije, iip1
        dxqu(ij) = dxqu(ij - iim)
      ENDDO

      DO ij = ijb + 1, ije
        zz(ij) = dxqu(ij - 1) * dxqu(ij)
        zz(ij) = zz(ij) + zz(ij)
        IF(zz(ij)>0) THEN
          dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
        ELSE
          !   extremum local
          dxq(ij, l) = 0.
        ENDIF
      ENDDO

    ENDDO
    !$OMP END DO NOWAIT
  ENDIF ! (pente_max.lt.-1.e-5)

  !WRITE(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x

  !   bouclage de la pente en iip1:
  !   -----------------------------
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb + iip1 - 1, ije, iip1
      dxq(ij - iim, l) = dxq(ij, l)
    ENDDO
    DO ij = ijb, ije
      iadvplus(ij, l) = 0
    ENDDO

  ENDDO
  !$OMP END DO NOWAIT
  ! PRINT*,'Bouclage en iip1'

  !   calcul des flux a gauche et a droite


  !   on cumule le flux correspondant a toutes les mailles dont la masse
  !   au travers de la paroi pENDant le pas de temps.
  ! PRINT*,'Cumule ....'
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  ! on a besoin de masse entre ijb et ije
  DO l = 1, llm
    DO ij = ijb, ije - 1
      ! PRINT*,'masse(',ij,')=',masse(ij,l,iq)
      IF (u_m(ij, l)>0.) THEN
        zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq)
        u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) &
                + 0.5 * zdum(ij, l) * dxq(ij, l))
      ELSE
        zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq)
        u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) &
                - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
      ENDIF
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  !   detection des points ou on advecte plus que la masse de la
  !   maille
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije - 1
      IF(zdum(ij, l)<0) THEN
        iadvplus(ij, l) = 1
        u_mq(ij, l) = 0.
      ENDIF
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  ! PRINT*,'Ok test 1'

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb + iip1 - 1, ije, iip1
      iadvplus(ij, l) = iadvplus(ij - iim, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  ! PRINT*,'Ok test 2'


  !   traitement special pour le cas ou on advecte en longitude plus que le
  !   contenu de la maille.
  !   cette partie est mal vectorisee.

  !  calcul du nombre de maille sur lequel on advecte plus que la maille.

  n0 = 0
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    nl(l) = 0
    DO ij = ijb, ije
      nl(l) = nl(l) + iadvplus(ij, l)
    ENDDO
    n0 = n0 + nl(l)
  ENDDO
  !$OMP END DO NOWAIT
  !ym      IF(n0.gt.1) THEN
  !ym      IF(n0.gt.0) THEN

  ! PRINT*,'Nombre de points pour lesquels on advect plus que le'
  ! &       ,'contenu de la maille : ',n0
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)

  DO l = 1, llm
    IF(nl(l)>0) THEN
      iju = 0
      !   indicage des mailles concernees par le traitement special
      DO ij = ijb, ije
        IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN
          iju = iju + 1
          indu(iju) = ij
        ENDIF
      ENDDO
      niju = iju
      !PRINT*,'vlx 278, niju,nl',niju,nl(l)

      !  traitement des mailles
      DO iju = 1, niju
        ij = indu(iju)
        j = (ij - 1) / iip1 + 1
        zu_m = u_m(ij, l)
        u_mq(ij, l) = 0.
        IF(zu_m>0.) THEN
          ijq = ij
          i = ijq - (j - 1) * iip1
          !   accumulation pour les mailles completements advectees
          do while(zu_m>masse(ijq, l, iq))
            u_mq(ij, l) = u_mq(ij, l) &
                    + q(ijq, l, iq) * masse(ijq, l, iq)
            zu_m = zu_m - masse(ijq, l, iq)
            i = mod(i - 2 + iim, iim) + 1
            ijq = (j - 1) * iip1 + i
          ENDDO
          !   ajout de la maille non completement advectee
          u_mq(ij, l) = u_mq(ij, l) + zu_m * &
                  (q(ijq, l, iq) + 0.5 * &
                          (1. - zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
        ELSE
          ijq = ij + 1
          i = ijq - (j - 1) * iip1
          !   accumulation pour les mailles completements advectees
          do while(-zu_m>masse(ijq, l, iq))
            u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) &
                    * masse(ijq, l, iq)
            zu_m = zu_m + masse(ijq, l, iq)
            i = mod(i, iim) + 1
            ijq = (j - 1) * iip1 + i
          ENDDO
          !   ajout de la maille non completement advectee
          u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - &
                  0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
        ENDIF
      ENDDO
    ENDIF
  ENDDO
  !$OMP END DO NOWAIT
  !ym      ENDIF  ! n0.gt.0

  !   bouclage en latitude
  ! PRINT*,'Avant bouclage en latitude'
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb + iip1 - 1, ije, iip1
      u_mq(ij, l) = u_mq(ij - iim, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ! CRisi: appel récursif de l'advection sur les fils.
  ! Il faut faire ça avant d'avoir mis à jour q et masse

  do ifils = 1, tracers(iq)%nqDescen
    ! attention: comme Ratio est utilisé comme q dans l'appel
    ! recursif, il doit contenir à lui seul tous les indices de tous
    ! les descendants!
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        ! On a besoin de q et masse seulement entre ijb et ije. On ne
        ! les calcule donc que de ijb à ije
        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
        IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020
          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
        else
          Ratio(ij, l, iq2) = min_ratio
        endif
      enddo
    enddo
    !$OMP END DO NOWAIT
  enddo !do ifils=1,tracers(iq)%nqDescen
  do ifils = 1, tracers(iq)%nqChildren
    iq2 = tracers(iq)%iqDescen(ifils)
    CALL vlx_loc(Ratio, pente_max, masse, u_mq, ijb_x, ije_x, iq2)
  enddo
  ! end CRisi


  !   calcul des tENDances
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb + 1, ije
      !MVals: veiller a ce qu'on n'ait pas de denominateur nul
      new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass)
      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + &
              u_mq(ij - 1, l) - u_mq(ij, l)) &
              / new_m
      masse(ij, l, iq) = new_m
    ENDDO
    !   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    DO ij = ijb + iip1 - 1, ije, iip1
      q(ij - iim, l, iq) = q(ij, l, iq)
      masse(ij - iim, l, iq) = masse(ij, l, iq)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ! retablir les fils en rapport de melange par rapport a l'air:
  ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio
  ! puis on boucle en longitude
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb + 1, ije
        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
      enddo
      DO ij = ijb + iip1 - 1, ije, iip1
        q(ij - iim, l, iq2) = q(ij, l, iq2)
      enddo
    enddo
    !$OMP END DO NOWAIT
  enddo

END SUBROUTINE vlx_loc


SUBROUTINE vly_loc(q, pente_max, masse, masse_adv_v, iq)
  !
  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
  !
  !    ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  !    ********************************************************************
  ! q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
  ! dq         sont des arguments de sortie pour le s-pg ....
  !
  !
  !   --------------------------------------------------------------------
  USE parallel_lmdz
  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
          min_qParent, min_qMass, min_ratio ! MVals et CRisi
  USE comconst_mod, ONLY: pi
  USE lmdz_ssum_scopy, ONLY: ssum

  IMPLICIT NONE
  !
  include "dimensions.h"
  include "paramet.h"
  include "comgeom.h"
  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
  REAL :: masse_adv_v(ijb_v:ije_v, llm)
  REAL :: q(ijb_u:ije_u, llm, nqtot), dq(ijb_u:ije_u, llm)
  INTEGER :: iq ! CRisi
  !
  !  Local
  !   ---------
  !
  INTEGER :: i, ij, l
  !
  REAL :: airej2, airejjm, airescb(iim), airesch(iim)
  REAL :: dyq(ijb_u:ije_u, llm), dyqv(ijb_v:ije_v), zdvm(ijb_u:ije_u, llm)
  REAL :: adyqv(ijb_v:ije_v), dyqmax(ijb_u:ije_u)
  REAL :: qbyv(ijb_v:ije_v, llm)

  REAL :: qpns, qpsn, appn, apps, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
  ! REAL newq,oldmasse
  Logical :: extremum, first
  REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second
  SAVE temps0, temps1, temps2, temps3, temps4, temps5
  !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
  SAVE first
  !$OMP THREADPRIVATE(first)

  REAL :: convpn, convps, convmpn, convmps
  REAL :: massepn, masseps, qpn, qps
  REAL :: sinlon(iip1), sinlondlon(iip1)
  REAL :: coslon(iip1), coslondlon(iip1)
  SAVE sinlon, coslon, sinlondlon, coslondlon
  !$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
  SAVE airej2, airejjm
  !$OMP THREADPRIVATE(airej2,airejjm)

  REAL :: Ratio(ijb_u:ije_u, llm, nqtot) ! CRisi
  INTEGER :: ifils, iq2 ! CRisi

  DATA first/.TRUE./
  DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./
  INTEGER :: ijb, ije
  INTEGER :: ijbm, ijem

  ijb = ij_begin - 2 * iip1
  ije = ij_end + 2 * iip1
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud)  ije = ij_end

  IF(first) THEN
    PRINT*, 'Shema  Amont nouveau  appele dans  Vanleer   '
    first = .FALSE.
    do i = 2, iip1
      coslon(i) = cos(rlonv(i))
      sinlon(i) = sin(rlonv(i))
      coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi
      sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi
    ENDDO
    coslon(1) = coslon(iip1)
    coslondlon(1) = coslondlon(iip1)
    sinlon(1) = sinlon(iip1)
    sinlondlon(1) = sinlondlon(iip1)
    airej2 = SSUM(iim, aire(iip2), 1)
    airejjm = SSUM(iim, aire(ip1jm - iim), 1)
  ENDIF

  !
  ! PRINT*,'CALCUL EN LATITUDE'

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    !
    !   --------------------------------
    !  CALCUL EN LATITUDE
    !   --------------------------------

    !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
    !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
    !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.

    IF (pole_nord) THEN
      DO i = 1, iim
        airescb(i) = aire(i + iip1) * q(i + iip1, l, iq)
      ENDDO
      qpns = SSUM(iim, airescb, 1) / airej2
    endif

    IF (pole_sud) THEN
      DO i = 1, iim
        airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq)
      ENDDO
      qpsn = SSUM(iim, airesch, 1) / airejjm
    endif

    !   calcul des pentes aux points v

    ijb = ij_begin - 2 * iip1
    ije = ij_end + iip1
    IF (pole_nord) ijb = ij_begin
    IF (pole_sud)  ije = ij_end - iip1

    ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1
    ! Si pole sud, entre ij_begin-2*iip1 et ij_end
    ! Si pole Nord, entre ij_begin et ij_end+2*iip1
    DO ij = ijb, ije
      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
      adyqv(ij) = abs(dyqv(ij))
    ENDDO


    !   calcul des pentes aux points scalaires
    ijb = ij_begin - iip1
    ije = ij_end + iip1
    IF (pole_nord) ijb = ij_begin + iip1
    IF (pole_sud)  ije = ij_end - iip1

    DO ij = ijb, ije
      dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij))
      dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij))
      dyqmax(ij) = pente_max * dyqmax(ij)
    ENDDO

    !   calcul des pentes aux poles
    IF (pole_nord) THEN
      DO ij = 1, iip1
        dyq(ij, l) = qpns - q(ij + iip1, l, iq)
      ENDDO

      dyn1 = 0.
      dyn2 = 0.
      DO ij = 1, iim
        dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
        dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
      ENDDO
      DO ij = 1, iip1
        dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
      ENDDO

      DO ij = 1, iip1
        dyq(ij, l) = 0.
      ENDDO
      ! ym tout cela ne sert pas a grand chose
    ENDIF

    IF (pole_sud) THEN

      DO ij = 1, iip1
        dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn
      ENDDO

      dys1 = 0.
      dys2 = 0.

      DO ij = 1, iim
        dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
        dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
      ENDDO

      DO ij = 1, iip1
        dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
      ENDDO

      DO ij = 1, iip1
        dyq(ip1jm + ij, l) = 0.
      ENDDO
      ! ym tout cela ne sert pas a grand chose
    ENDIF

    !   filtrage de la derivee

    !   calcul des pentes limites aux poles
    ! ym partie inutile
    ! goto 8888
    ! fn=1.
    ! fs=1.
    ! DO ij=1,iim
    !    IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
    !       fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
    !    ENDIF
    ! IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
    !    fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
    !    ENDIF
    ! ENDDO
    ! DO ij=1,iip1
    !    dyq(ij,l)=fn*dyq(ij,l)
    !    dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
    ! ENDDO
    ! 8888    continue


    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    !  En memoire de dIFferents tests sur la
    !  limitation des pentes aux poles.
    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    ! PRINT*,dyq(1)
    ! PRINT*,dyqv(iip1+1)
    ! appn=abs(dyq(1)/dyqv(iip1+1))
    ! PRINT*,dyq(ip1jm+1)
    ! PRINT*,dyqv(ip1jm-iip1+1)
    ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    ! DO ij=2,iim
    !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
    !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    ! ENDDO
    ! appn=min(pente_max/appn,1.)
    ! apps=min(pente_max/apps,1.)
    !
    !
    !   cas ou on a un extremum au pole
    !
    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    !    &   appn=0.
    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    !    &   apps=0.
    !
    !   limitation des pentes aux poles
    ! DO ij=1,iip1
    !    dyq(ij)=appn*dyq(ij)
    !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    ! ENDDO
    !
    !   test
    !  DO ij=1,iip1
    !     dyq(iip1+ij)=0.
    !     dyq(ip1jm+ij-iip1)=0.
    !  ENDDO
    !  DO ij=1,ip1jmp1
    !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
    !  ENDDO
    !
    ! changement 10 07 96
    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    !    &   THEN
    !    DO ij=1,iip1
    !       dyqmax(ij)=0.
    !    ENDDO
    ! ELSE
    !    DO ij=1,iip1
    !       dyqmax(ij)=pente_max*abs(dyqv(ij))
    !    ENDDO
    ! ENDIF
    !
    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    !    &THEN
    !    DO ij=ip1jm+1,ip1jmp1
    !       dyqmax(ij)=0.
    !    ENDDO
    ! ELSE
    !    DO ij=ip1jm+1,ip1jmp1
    !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
    !    ENDDO
    ! ENDIF
    !   fin changement 10 07 96
    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

    !   calcul des pentes limitees
    ijb = ij_begin - iip1
    ije = ij_end + iip1
    IF (pole_nord) ijb = ij_begin + iip1
    IF (pole_sud)  ije = ij_end - iip1

    DO ij = ijb, ije
      IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN
        dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l))
      ELSE
        dyq(ij, l) = 0.
      ENDIF
    ENDDO

  ENDDO
  !$OMP END DO NOWAIT

  ijb = ij_begin - iip1
  ije = ij_end
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud)  ije = ij_end - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      IF(masse_adv_v(ij, l)>0) THEN
        qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * &
                0.5 * (1. - masse_adv_v(ij, l) &
                / masse(ij + iip1, l, iq))
      ELSE
        qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * &
                0.5 * (1. + masse_adv_v(ij, l) / masse(ij, l, iq))
      ENDIF
      qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l)
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ! CRisi: appel récursif de l'advection sur les fils.
  ! Il faut faire ça avant d'avoir mis à jour q et masse
  ! WRITE(*,*)'vly 689: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren

  ijb = ij_begin - 2 * iip1
  ije = ij_end + 2 * iip1
  ijbm = ij_begin - iip1
  ijem = ij_end + iip1
  IF (pole_nord) ijb = ij_begin
  IF (pole_sud)  ije = ij_end
  IF (pole_nord) ijbm = ij_begin
  IF (pole_sud)  ijem = ij_end

  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      ! modif des bornes: CRisi 16 nov 2020
      ! d'abord masse avec bornes corrigées
      DO ij = ijbm, ijem
        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
      enddo

      ! ensuite Ratio avec anciennes bornes
      DO ij = ijb, ije
        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        IF (q(ij, l, iq)>min_qParent) then ! modif 13 nov 2020
          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
        else
          Ratio(ij, l, iq2) = min_ratio
        endif
      enddo !DO ij=ijbm,ijem
    enddo !DO l=1,llm
    !$OMP END DO NOWAIT
  enddo

  do ifils = 1, tracers(iq)%nqChildren
    iq2 = tracers(iq)%iqDescen(ifils)
    CALL vly_loc(Ratio, pente_max, masse, qbyv, iq2)
  enddo
  ! end CRisi

  ijb = ij_begin
  ije = ij_end
  IF (pole_nord) ijb = ij_begin + iip1
  IF (pole_sud)  ije = ij_end - iip1

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      newmasse = masse(ij, l, iq) &
              + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l)

      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) &
              - qbyv(ij - iip1, l)) / newmasse

      masse(ij, l, iq) = newmasse

    ENDDO

    IF (pole_nord) THEN
      convpn = SSUM(iim, qbyv(1, l), 1)
      convmpn = ssum(iim, masse_adv_v(1, l), 1)
      massepn = ssum(iim, masse(1, l, iq), 1)
      qpn = 0.
      do ij = 1, iim
        qpn = qpn + masse(ij, l, iq) * q(ij, l, iq)
      enddo
      qpn = (qpn + convpn) / (massepn + convmpn)
      do ij = 1, iip1
        q(ij, l, iq) = qpn
      enddo
    endif

    IF (pole_sud) THEN
      convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1)
      convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1)
      masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1)
      qps = 0.
      do ij = ip1jm + 1, ip1jmp1 - 1
        qps = qps + masse(ij, l, iq) * q(ij, l, iq)
      enddo
      qps = (qps + convps) / (masseps + convmps)
      do ij = ip1jm + 1, ip1jmp1
        q(ij, l, iq) = qps
      enddo
    endif
  ENDDO
  !$OMP END DO NOWAIT

  ! retablir les fils en rapport de melange par rapport a l'air:
  ijb = ij_begin
  ije = ij_end
  ! if (pole_nord) ijb=ij_begin
  ! if (pole_sud)  ije=ij_end

  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
      enddo
    enddo
    !$OMP END DO NOWAIT
  enddo

END SUBROUTINE vly_loc


SUBROUTINE vlz_loc(q, pente_max, masse, w, ijb_x, ije_x, iq)
  !
  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
  !
  !    ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  !    ********************************************************************
  !    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
  ! dq         sont des arguments de sortie pour le s-pg ....
  !
  !
  !   --------------------------------------------------------------------
  USE parallel_lmdz
  USE vlz_mod
  USE infotrac, ONLY: nqtot, tracers, & ! CRisi                 &
          min_qParent, min_qMass, min_ratio ! MVals et CRisi
  USE lmdz_iniprint, ONLY: lunout, prt_level


  IMPLICIT NONE
  !
  include "dimensions.h"
  include "paramet.h"
  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ijb_u:ije_u, llm, nqtot), pente_max
  REAL :: q(ijb_u:ije_u, llm, nqtot)
  REAL :: w(ijb_u:ije_u, llm + 1, nqtot)
  INTEGER :: iq
  !
  !  Local
  !   ---------
  !
  INTEGER :: i, ij, l, j, ii

  REAL, DIMENSION(ijb_u:ije_u, llm + 1) :: wresi, morig, qorig, dzqorig
  INTEGER, DIMENSION(ijb_u:ije_u, llm + 1) :: lorig
  INTEGER, SAVE :: countcfl
  !$OMP THREADPRIVATE(countcfl)
  !
  REAL :: newmasse

  REAL :: dzqmax
  REAL :: sigw

  REAL :: temps0, temps1, temps2, temps3, temps4, temps5, second
  SAVE temps0, temps1, temps2, temps3, temps4, temps5
  !$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)

  DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./
  INTEGER :: ijb, ije, ijb_x, ije_x
  LOGICAL, SAVE :: first = .TRUE.
  !$OMP THREADPRIVATE(first)

  !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
  ! Ces varibles doivent être déclarées en pointer et en save dans
  ! vlz_loc si on veut qu'elles soient vues par tous les threads.
  INTEGER :: ifils, iq2 ! CRisi

  IF (first) THEN
    first = .FALSE.
  ENDIF
  !    On oriente tout dans le sens de la pression c'est a dire dans le
  !    sens de W

  !WRITE(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq

  ijb = ijb_x
  ije = ije_x

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 2, llm
    DO ij = ijb, ije
      dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq)
      adzqw(ij, l) = abs(dzqw(ij, l))
    ENDDO
  ENDDO
  !$OMP END DO

  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 2, llm - 1
    DO ij = ijb, ije
      IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN
        dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1))
      ELSE
        dzq(ij, l) = 0.
      ENDIF
      dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1))
      dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l))
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  !$OMP MASTER
  DO ij = ijb, ije
    dzq(ij, 1) = 0.
    dzq(ij, llm) = 0.
  ENDDO
  !$OMP END MASTER
  !$OMP BARRIER

  !--------------------------------------------------------
  ! On repere les points qui violent le CFL (|w| > masse)
  !--------------------------------------------------------

  countcfl = 0
  ! PRINT*,'vlz nouveau'
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 2, llm
    DO ij = ijb, ije
      IF((w(ij, l, iq)>0.AND.w(ij, l, iq)>masse(ij, l, iq)) &
              .OR. (w(ij, l, iq)<=0.AND.ABS(w(ij, l, iq))>masse(ij, l - 1, iq))) &
              countcfl = countcfl + 1
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  ! ---------------------------------------------------------------
  !  Identification des mailles ou on viole le CFL : w > masse
  ! ---------------------------------------------------------------

  IF (countcfl==0) THEN

    ! ---------------------------------------------------------------
    !   .... calcul des termes d'advection verticale  .......
    ! Dans le cas où le  |w| < masse partout.
    ! Version d'origine
    ! Pourrait etre enleve si on voit que le code plus general
    ! est aussi rapide
    ! ---------------------------------------------------------------

    ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq

    !  !WRITE(*,*) 'vlz 982,ijb,ije=',ijb,ije
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm - 1
      do  ij = ijb, ije
        IF(w(ij, l + 1, iq)>0.) THEN
          sigw = w(ij, l + 1, iq) / masse(ij, l + 1, iq)
          wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l + 1, iq) &
                  + 0.5 * (1. - sigw) * dzq(ij, l + 1))
        ELSE
          sigw = w(ij, l + 1, iq) / masse(ij, l, iq)
          wq(ij, l + 1, iq) = w(ij, l + 1, iq) * (q(ij, l, iq) &
                  - 0.5 * (1. + sigw) * dzq(ij, l))
        ENDIF
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT
    !WRITE(*,*) 'vlz 1001'

  ELSE ! countcfl>=1

    IF (prt_level>9) THEN
      WRITE(lunout, *)'vlz passage dans le non local'
    ENDIF
    ! ---------------------------------------------------------------
    !  Debut du traitement du cas ou on viole le CFL : w > masse
    ! ---------------------------------------------------------------

    ! Initialisation

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 2, llm
      DO ij = ijb, ije
        wresi(ij, l) = w(ij, l, iq)
        wq(ij, l, iq) = 0.
        IF(w(ij, l, iq)>0.) THEN
          lorig(ij, l) = l
          morig(ij, l) = masse(ij, l, iq)
          qorig(ij, l) = q(ij, l, iq)
          dzqorig(ij, l) = dzq(ij, l)
        ELSE
          lorig(ij, l) = l - 1
          morig(ij, l) = masse(ij, l - 1, iq)
          qorig(ij, l) = q(ij, l - 1, iq)
          dzqorig(ij, l) = dzq(ij, l - 1)
        ENDIF
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT

    ! Reindicage vertical en accumulant les flux sur
    !  les mailles qui viollent le CFL
    !  on itère jusqu'à ce que tous les poins satisfassent
    !  le critère
    DO WHILE (countcfl>=1)
      IF (prt_level>9) THEN
        WRITE(lunout, *)'On viole le CFL Vertical sur ', countcfl, ' pts'
      ENDIF
      countcfl = 0

      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 2, llm
        DO ij = ijb, ije
          IF (ABS(wresi(ij, l))>morig(ij, l)) THEN
            countcfl = countcfl + 1
            ! rm : les 8 lignes ci dessous pourraient sans doute s'ecrire
            ! avec la fonction sign
            IF(w(ij, l, iq)>0.) THEN
              wresi(ij, l) = wresi(ij, l) - morig(ij, l)
              wq(ij, l, iq) = wq(ij, l, iq) + morig(ij, l) * qorig(ij, l)
              lorig(ij, l) = lorig(ij, l) + 1
            ELSE
              wresi(ij, l) = wresi(ij, l) + morig(ij, l)
              wq(ij, l, iq) = wq(ij, l, iq) - morig(ij, l) * qorig(ij, l)
              lorig(ij, l) = lorig(ij, l) - 1
            ENDIF
            ! CRisi 24nov2020: ajout d'un message d'erreur clair au lieu d'un plantage
            ! pour seg fault
            IF (lorig(ij, l)==0) THEN
              CALL abort_gcm("vlz in vlsplt_loc", &
                      "unfixable violation of CFL", 1)
            endif
            morig(ij, l) = masse(ij, lorig(ij, l), iq)
            qorig(ij, l) = q(ij, lorig(ij, l), iq)
            dzqorig(ij, l) = dzq(ij, lorig(ij, l))
          ENDIF
        ENDDO
      ENDDO
      !$OMP END DO NOWAIT

    ENDDO ! WHILE (countcfl>=1)

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 2, llm
      do  ij = ijb, ije
        sigw = wresi(ij, l) / morig(ij, l)
        IF(w(ij, l, iq)>0.) THEN
          wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) &
                  + 0.5 * (1. - sigw) * dzqorig(ij, l))
        ELSE
          wq(ij, l, iq) = wq(ij, l, iq) + wresi(ij, l) * (qorig(ij, l) &
                  - 0.5 * (1. + sigw) * dzqorig(ij, l))
        ENDIF
      ENDDO
    ENDDO
    !$OMP END DO NOWAIT

  ENDIF ! councfl=0



  !$OMP MASTER
  DO ij = ijb, ije
    wq(ij, llm + 1, iq) = 0.
    wq(ij, 1, iq) = 0.
  ENDDO
  !$OMP END MASTER
  !$OMP BARRIER

  ! CRisi: appel récursif de l'advection sur les fils.
  ! Il faut faire ça avant d'avoir mis à jour q et masse
  ! WRITE(*,*)'vlsplt 942: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        masse(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
        IF (q(ij, l, iq)>min_qParent) THEN
          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
        else
          Ratio(ij, l, iq2) = min_ratio
        endif
        !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
        w(ij, l, iq2) = wq(ij, l, iq)
      enddo
    enddo
    !$OMP END DO NOWAIT
  enddo
  !$OMP BARRIER

  do ifils = 1, tracers(iq)%nqChildren
    iq2 = tracers(iq)%iqDescen(ifils)
    CALL vlz_loc(Ratio, pente_max, masse, w, ijb_x, ije_x, iq2)
  enddo
  ! end CRisi

  ! CRisi: On rajoute ici une barrière car on veut être sur que tous les
  ! wq soient synchronisés

  !$OMP BARRIER
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    DO ij = ijb, ije
      newmasse = masse(ij, l, iq) + w(ij, l + 1, iq) - w(ij, l, iq)
      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) &
              + wq(ij, l + 1, iq) - wq(ij, l, iq)) &
              / newmasse
      masse(ij, l, iq) = newmasse
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT


  ! retablir les fils en rapport de melange par rapport a l'air:
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      DO ij = ijb, ije
        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
      enddo
    enddo
    !$OMP END DO NOWAIT
  enddo

END SUBROUTINE vlz_loc