vlspltqs.F90 Source File

!
! $Id: vlspltqs.F90 5285 2024-10-28 13:33:29Z 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
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !



  !
  !   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
  LOGICAL :: testcpu
  SAVE testcpu
  SAVE temps1,temps2,temps3

  REAL :: qmin,qmax
  DATA qmin,qmax/0.,1.e33/
  DATA testcpu/.false./
  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 " .
  !    ********************************************************************
  !
  !   --------------------------------------------------------------------
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !


  !
  !
  !   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,testcpu
  SAVE first,testcpu

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


  DATA first,testcpu/.true.,.false./

  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.gt.-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
#ifdef CRAY
           dxq(ij,l)= &
                 cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
#else
           IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
              dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
           ELSE
  !   extremum local
              dxq(ij,l)=0.
           ENDIF
#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).gt.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

#ifdef CRAY
  !--pas encore modification sur Qsat
  DO l=1,llm
   DO ij=iip2,ip1jm-1
      zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), &
            1.+u_m(ij,l)/masse(ij+1,l,iq), &
            u_m(ij,l))
      zdum(ij,l)=0.5*zdum(ij,l)
      u_mq(ij,l)=cvmgp( &
            q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), &
            q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), &
            u_m(ij,l))
      u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
   ENDDO
  ENDDO
#else
  !   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).gt.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
#endif


  !   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).lt.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.gt.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).gt.0) THEN
           iju=0
  !   indicage des mailles concernees par le traitement special
           DO ij=iip2,ip1jm
              IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.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.gt.0.) THEN
                 ijq=ij
                 i=ijq-(j-1)*iip1
  !   accumulation pour les mailles completements advectees
                 do while(zu_m.gt.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.gt.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 comgeom_mod_h
  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

  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !


  !
  !
  !   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,testcpu
  REAL :: temps0,temps1,temps2,temps3,temps4,temps5
  SAVE temps0,temps1,temps2,temps3,temps4,temps5
  SAVE first,testcpu

  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,testcpu/.true.,.false./
  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).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

  !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).gt.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).GT.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