!
! $Id: vlspltqs.F90 5103 2024-07-23 13:29:36Z abarral $
!
SUBROUTINE vlspltqs (q, pente_max, masse, w, pbaru, pbarv, pdt, &
        p, pk, teta, iq)
  USE infotrac, ONLY: nqtot, tracers
  !
  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
  !
  !    ********************************************************************
  !      Shema  d'advection " pseudo amont " .
  !  + test sur humidite specifique: Q advecte< Qsat aval
  !               (F. Codron, 10/99)
  !    ********************************************************************
  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
  !
  ! pente_max facteur de limitation des pentes: 2 en general
  !                                            0 pour un schema amont
  ! pbaru,pbarv,w flux de masse en u ,v ,w
  ! pdt pas de temps
  !
  ! teta temperature potentielle, p pression aux interfaces,
  ! pk exner au milieu des couches necessaire pour calculer Qsat
  !   --------------------------------------------------------------------

  USE comconst_mod, ONLY: cpp
  USE logic_mod, ONLY: adv_qsat_liq
  IMPLICIT NONE
  !
  include "dimensions.h"
  include "paramet.h"

  !
  !   Arguments:
  !   ----------
  REAL :: masse(ip1jmp1, llm), pente_max
  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
  REAL :: q(ip1jmp1, llm, nqtot)
  REAL :: w(ip1jmp1, llm), pdt
  REAL :: p(ip1jmp1, llmp1), teta(ip1jmp1, llm), pk(ip1jmp1, llm)
  INTEGER :: iq ! CRisi
  !
  !  Local
  !   ---------
  !
  INTEGER :: i, ij, l, j, ii
  INTEGER :: ifils, iq2 ! CRisi
  !
  REAL :: qsat(ip1jmp1, llm)
  REAL :: zm(ip1jmp1, llm, nqtot)
  REAL :: mu(ip1jmp1, llm)
  REAL :: mv(ip1jm, llm)
  REAL :: mw(ip1jmp1, llm + 1)
  REAL :: zq(ip1jmp1, llm, nqtot)
  REAL :: temps1, temps2, temps3
  REAL :: zzpbar, zzw
  SAVE temps1, temps2, temps3

  REAL :: qmin, qmax
  DATA qmin, qmax/0., 1.e33/
  DATA temps1, temps2, temps3/0., 0., 0./

  !--pour rapport de melange saturant--

  REAL :: rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
  REAL :: ptarg, pdelarg, foeew, zdelta
  REAL :: tempe(ip1jmp1)

  !    fonction psat(T)

  FOEEW (PTARG, PDELARG) = EXP (&
          (R3LES * (1. - PDELARG) + R3IES * PDELARG) * (PTARG - RTT) &
                  / (PTARG - (R4LES * (1. - PDELARG) + R4IES * PDELARG)))

  r2es = 380.11733
  r3les = 17.269
  r3ies = 21.875
  r4les = 35.86
  r4ies = 7.66
  retv = 0.6077667
  rtt = 273.16

  !-- Calcul de Qsat en chaque point
  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
  !   pour eviter une exponentielle.
  DO l = 1, llm
    DO ij = 1, ip1jmp1
      tempe(ij) = teta(ij, l) * pk(ij, l) / cpp
    ENDDO
    DO ij = 1, ip1jmp1
      IF (adv_qsat_liq) THEN
        zdelta = 0.
      ELSE
        zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
      ENDIF
      play = 0.5 * (p(ij, l) + p(ij, l + 1))
      qsat(ij, l) = MIN(0.5, r2es * FOEEW(tempe(ij), zdelta) / play)
      qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
    ENDDO
  ENDDO

  ! PRINT*,'Debut vlsplt version debug sans vlyqs'

  zzpbar = 0.5 * pdt
  zzw = pdt
  DO l = 1, llm
    DO ij = iip2, ip1jm
      mu(ij, l) = pbaru(ij, l) * zzpbar
    ENDDO
    DO ij = 1, ip1jm
      mv(ij, l) = pbarv(ij, l) * zzpbar
    ENDDO
    DO ij = 1, ip1jmp1
      mw(ij, l) = w(ij, l) * zzw
    ENDDO
  ENDDO

  DO ij = 1, ip1jmp1
    mw(ij, llm + 1) = 0.
  ENDDO

  CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1)
  CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1)
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1)
  enddo

  ! CALL minmaxq(zq,qmin,qmax,'avant vlxqs     ')
  CALL vlxqs(zq, pente_max, zm, mu, qsat, iq)

  ! CALL minmaxq(zq,qmin,qmax,'avant vlyqs     ')

  CALL vlyqs(zq, pente_max, zm, mv, qsat, iq)

  ! CALL minmaxq(zq,qmin,qmax,'avant vlz     ')

  CALL vlz(zq, pente_max, zm, mw, iq)

  ! CALL minmaxq(zq,qmin,qmax,'avant vlyqs     ')
  ! CALL minmaxq(zm,qmin,qmax,'M avant vlyqs     ')

  CALL vlyqs(zq, pente_max, zm, mv, qsat, iq)

  ! CALL minmaxq(zq,qmin,qmax,'avant vlxqs     ')
  ! CALL minmaxq(zm,qmin,qmax,'M avant vlxqs     ')

  CALL vlxqs(zq, pente_max, zm, mu, qsat, iq)

  ! CALL minmaxq(zq,qmin,qmax,'apres vlxqs     ')
  ! CALL minmaxq(zm,qmin,qmax,'M apres vlxqs     ')

  DO l = 1, llm
    DO ij = 1, ip1jmp1
      q(ij, l, iq) = zq(ij, l, iq)
    ENDDO
    DO ij = 1, ip1jm + 1, iip1
      q(ij + iim, l, iq) = q(ij, l, iq)
    ENDDO
  ENDDO
  ! ! CRisi: aussi pour les fils
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    DO l = 1, llm
      DO ij = 1, ip1jmp1
        q(ij, l, iq2) = zq(ij, l, iq2)
      ENDDO
      DO ij = 1, ip1jm + 1, iip1
        q(ij + iim, l, iq2) = q(ij, l, iq2)
      ENDDO
    ENDDO
  enddo
  ! !write(*,*) 'vlspltqs 183: fin de la routine'

  RETURN
END SUBROUTINE vlspltqs
SUBROUTINE vlxqs(q, pente_max, masse, u_m, qsat, iq)
  USE infotrac, ONLY: nqtot, tracers ! CRisi

  !
  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
  !
  !    ********************************************************************
  ! Shema  d'advection " pseudo amont " .
  !    ********************************************************************
  !
  !   --------------------------------------------------------------------
  IMPLICIT NONE
  !
  include "dimensions.h"
  include "paramet.h"
  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
  REAL :: u_m(ip1jmp1, llm)
  REAL :: q(ip1jmp1, llm, nqtot)
  REAL :: qsat(ip1jmp1, llm)
  INTEGER :: iq ! CRisi
  !
  !  Local
  !   ---------
  !
  INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
  INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm)
  !
  REAL :: new_m, zu_m, zdum(ip1jmp1, llm)
  REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1)
  REAL :: zz(ip1jmp1)
  REAL :: adxqu(ip1jmp1), dxqmax(ip1jmp1, llm)
  REAL :: u_mq(ip1jmp1, llm)

  ! ! CRisi
  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot)
  INTEGER :: ifils, iq2 ! CRisi

  Logical :: first
  SAVE first

  REAL :: SSUM
  REAL :: temps0, temps1, temps2, temps3, temps4, temps5
  SAVE temps0, temps1, temps2, temps3, temps4, temps5

  DATA first/.TRUE./

  IF(first) THEN
    temps1 = 0.
    temps2 = 0.
    temps3 = 0.
    temps4 = 0.
    temps5 = 0.
    first = .FALSE.
  ENDIF

  !   calcul de la pente a droite et a gauche de la maille

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

    !   calcul des pentes avec limitation, Van Leer scheme I:
    !   -----------------------------------------------------

    !   calcul de la pente aux points u
    DO l = 1, llm
      DO ij = iip2, ip1jm - 1
        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
      ENDDO
      DO ij = iip1 + iip1, ip1jm, iip1
        dxqu(ij) = dxqu(ij - iim)
        ! sigu(ij)=sigu(ij-iim)
      ENDDO

      DO ij = iip2, ip1jm
        adxqu(ij) = abs(dxqu(ij))
      ENDDO

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

      DO ij = iip2 + 1, ip1jm
        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 = iip1 + iip1, ip1jm, iip1
        dxqmax(ij - iim, l) = dxqmax(ij, l)
      ENDDO

      DO ij = iip2 + 1, ip1jm
        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

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

    !   Pentes produits:
    !   ----------------

    DO l = 1, llm
      DO ij = iip2, ip1jm - 1
        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
      ENDDO
      DO ij = iip1 + iip1, ip1jm, iip1
        dxqu(ij) = dxqu(ij - iim)
      ENDDO

      DO ij = iip2 + 1, ip1jm
        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

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

  !   bouclage de la pente en iip1:
  !   -----------------------------

  DO l = 1, llm
    DO ij = iip1 + iip1, ip1jm, iip1
      dxq(ij - iim, l) = dxq(ij, l)
    ENDDO

    DO ij = 1, ip1jmp1
      iadvplus(ij, l) = 0
    ENDDO

  ENDDO


  !   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.
  !   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
  DO l = 1, llm
    DO ij = iip2, ip1jm - 1
      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) * &
                min(q(ij, l, iq) + 0.5 * zdum(ij, l) * dxq(ij, l), qsat(ij + 1, l))
      ELSE
        zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq)
        u_mq(ij, l) = u_m(ij, l) * &
                min(q(ij + 1, l, iq) - 0.5 * zdum(ij, l) * dxq(ij + 1, l), qsat(ij, l))
      ENDIF
    ENDDO
  ENDDO


  !   detection des points ou on advecte plus que la masse de la
  !   maille
  DO l = 1, llm
    DO ij = iip2, ip1jm - 1
      IF(zdum(ij, l)<0) THEN
        iadvplus(ij, l) = 1
        u_mq(ij, l) = 0.
      ENDIF
    ENDDO
  ENDDO
  DO l = 1, llm
    DO ij = iip1 + iip1, ip1jm, iip1
      iadvplus(ij, l) = iadvplus(ij - iim, l)
    ENDDO
  ENDDO



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

  !   pas d'influence de la pression saturante (pour l'instant)

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

  n0 = 0
  DO l = 1, llm
    nl(l) = 0
    DO ij = iip2, ip1jm
      nl(l) = nl(l) + iadvplus(ij, l)
    ENDDO
    n0 = n0 + nl(l)
  ENDDO

  IF(n0>0) THEN
    !cc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
    !cc     &       ,'contenu de la maille : ',n0

    DO l = 1, llm
      IF(nl(l)>0) THEN
        iju = 0
        !   indicage des mailles concernees par le traitement special
        DO ij = iip2, ip1jm
          IF(iadvplus(ij, l)==1.and.mod(ij, iip1)/=0) THEN
            iju = iju + 1
            indu(iju) = ij
          ENDIF
        ENDDO
        niju = iju
        ! PRINT*,'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
  ENDIF  ! n0.gt.0



  !   bouclage en latitude

  DO l = 1, llm
    DO ij = iip1 + iip1, ip1jm, iip1
      u_mq(ij, l) = u_mq(ij - iim, l)
    ENDDO
  ENDDO

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

  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    DO l = 1, llm
      DO ij = iip2, ip1jm
        ! ! On a besoin de q et masse seulement entre iip2 et ip1jm
        masseq(ij, l, iq2) = masse(ij, l, iq) * q(ij, l, iq)
        Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
      enddo
    enddo
  enddo
  do ifils = 1, tracers(iq)%nqChildren
    iq2 = tracers(iq)%iqDescen(ifils)
    CALL vlx(Ratio, pente_max, masseq, u_mq, iq2)
  enddo
  ! end CRisi

  !   calcul des tendances

  DO l = 1, llm
    DO ij = iip2 + 1, ip1jm
      new_m = masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l)
      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 = iip1 + iip1, ip1jm, iip1
      q(ij - iim, l, iq) = q(ij, l, iq)
      masse(ij - iim, l, iq) = masse(ij, l, iq)
    ENDDO
  ENDDO

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

  ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
  ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)

  RETURN
END SUBROUTINE vlxqs
SUBROUTINE vlyqs(q, pente_max, masse, masse_adv_v, qsat, iq)
  USE infotrac, ONLY: nqtot, tracers ! CRisi
  !
  ! 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 ....
  !     qsat 	       est   un argument de sortie pour le s-pg ....
  !
  !
  !   --------------------------------------------------------------------

  USE comconst_mod, ONLY: pi

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

  REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
  ! REAL newq,oldmasse
  Logical :: first
  REAL :: temps0, temps1, temps2, temps3, temps4, temps5
  SAVE temps0, temps1, temps2, temps3, temps4, temps5
  SAVE first

  REAL :: convpn, convps, convmpn, convmps
  REAL :: sinlon(iip1), sinlondlon(iip1)
  REAL :: coslon(iip1), coslondlon(iip1)
  SAVE sinlon, coslon, sinlondlon, coslondlon
  SAVE airej2, airejjm

  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
  INTEGER :: ifils, iq2 ! CRisi
  !
  !
  REAL :: SSUM

  DATA first/.TRUE./
  DATA temps0, temps1, temps2, temps3, temps4, temps5/0., 0., 0., 0., 0., 0./

  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

  !

  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.

    DO i = 1, iim
      airescb(i) = aire(i + iip1) * q(i + iip1, l, iq)
      airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq)
    ENDDO
    qpns = SSUM(iim, airescb, 1) / airej2
    qpsn = SSUM(iim, airesch, 1) / airejjm

    !   calcul des pentes aux points v

    DO ij = 1, ip1jm
      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
      adyqv(ij) = abs(dyqv(ij))
    ENDDO

    !   calcul des pentes aux points scalaires

    DO ij = iip2, ip1jm
      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

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

    !   filtrage de la derivee
    dyn1 = 0.
    dys1 = 0.
    dyn2 = 0.
    dys2 = 0.
    DO ij = 1, iim
      dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
      dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
      dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
      dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
    ENDDO
    DO ij = 1, iip1
      dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
      dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
    ENDDO

    !   calcul des pentes limites aux poles

    fn = 1.
    fs = 1.
    DO ij = 1, iim
      IF(pente_max * adyqv(ij)<abs(dyq(ij, l))) THEN
        fn = min(pente_max * adyqv(ij) / abs(dyq(ij, l)), fn)
      ENDIF
      IF(pente_max * adyqv(ij + ip1jm - iip1)<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

    !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

    DO ij = iip2, ip1jm
      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

  DO l = 1, llm
    DO ij = 1, ip1jm
      IF(masse_adv_v(ij, l)>0.) THEN
        qbyv(ij, l) = MIN(qsat(ij + iip1, 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) = MIN(qsat(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


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

  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    DO l = 1, llm
      DO ij = 1, ip1jmp1
        masseq(ij, l, iq2) = masse(ij, l, iq) * q(ij, l, iq)
        Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
      enddo
    enddo
  enddo
  do ifils = 1, tracers(iq)%nqChildren
    iq2 = tracers(iq)%iqDescen(ifils)
    ! !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
    CALL vly(Ratio, pente_max, masseq, qbyv, iq2)
  enddo

  DO l = 1, llm
    DO ij = iip2, ip1jm
      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
    !.-. ancienne version
    convpn = SSUM(iim, qbyv(1, l), 1) / apoln
    convmpn = ssum(iim, masse_adv_v(1, l), 1) / apoln
    DO ij = 1, iip1
      newmasse = masse(ij, l, iq) + convmpn * aire(ij)
      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + convpn * aire(ij)) / &
              newmasse
      masse(ij, l, iq) = newmasse
    ENDDO
    convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1) / apols
    convmps = -SSUM(iim, masse_adv_v(ip1jm - iim, l), 1) / apols
    DO ij = ip1jm + 1, ip1jmp1
      newmasse = masse(ij, l, iq) + convmps * aire(ij)
      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + convps * aire(ij)) / &
              newmasse
      masse(ij, l, iq) = newmasse
    ENDDO
    !.-. fin ancienne version

    !._. nouvelle version
    ! convpn=SSUM(iim,qbyv(1,l),1)
    ! convmpn=ssum(iim,masse_adv_v(1,l),1)
    ! oldmasse=ssum(iim,masse(1,l),1)
    ! newmasse=oldmasse+convmpn
    ! newq=(q(1,l)*oldmasse+convpn)/newmasse
    ! newmasse=newmasse/apoln
    ! DO ij = 1,iip1
    !    q(ij,l)=newq
    !    masse(ij,l,iq)=newmasse*aire(ij)
    ! ENDDO
    ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
    ! newmasse=oldmasse+convmps
    ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
    ! newmasse=newmasse/apols
    ! DO ij = ip1jm+1,ip1jmp1
    !    q(ij,l)=newq
    !    masse(ij,l,iq)=newmasse*aire(ij)
    ! ENDDO
    !._. fin nouvelle version
  ENDDO

  ! !write(*,*) 'vly 866'

  ! retablir les fils en rapport de melange par rapport a l'air:
  do ifils = 1, tracers(iq)%nqDescen
    iq2 = tracers(iq)%iqDescen(ifils)
    DO l = 1, llm
      DO ij = 1, ip1jmp1
        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
      enddo
    enddo
  enddo
  ! !write(*,*) 'vly 879'

  RETURN
END SUBROUTINE vlyqs