source: LMDZ6/trunk/libf/dyn3d/vlsplt.F90 @ 5423

!
! $Id: vlsplt.F90 5285 2024-10-28 13:33:29Z evignon $
!

SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
  USE infotrac, ONLY: nqtot,tracers
  !
  ! 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 ....
  !
  !   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
  !
  !   --------------------------------------------------------------------
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !



  !
  !   Arguments:
  !   ----------
  REAL :: masse(ip1jmp1,llm),pente_max
   ! REAL masse(iip1,jjp1,llm),pente_max
  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
  REAL :: q(ip1jmp1,llm,nqtot)
   ! REAL q(iip1,jjp1,llm)
  REAL :: w(ip1jmp1,llm),pdt
  INTEGER :: iq ! CRisi
  !
  !  Local
  !   ---------
  !
  INTEGER :: ij,l
  !
  REAL :: zm(ip1jmp1,llm,nqtot)
  REAL :: mu(ip1jmp1,llm)
  REAL :: mv(ip1jm,llm)
  REAL :: mw(ip1jmp1,llm+1)
  REAL :: zq(ip1jmp1,llm,nqtot)
  REAL :: zzpbar, zzw
  INTEGER :: ifils,iq2 ! CRisi

  REAL :: qmin,qmax
  DATA qmin,qmax/0.,1.e33/

    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

  !print*,'Entree vlx1'
  ! call minmaxq(zq,qmin,qmax,'avant vlx     ')
  call vlx(zq,pente_max,zm,mu,iq)
  !print*,'Sortie vlx1'
  ! call minmaxq(zq,qmin,qmax,'apres vlx1    ')

  ! print*,'Entree vly1'

  call vly(zq,pente_max,zm,mv,iq)
  ! call minmaxq(zq,qmin,qmax,'apres vly1     ')
  !print*,'Sortie vly1'
  call vlz(zq,pente_max,zm,mw,iq)
  ! call minmaxq(zq,qmin,qmax,'apres vlz     ')


  call vly(zq,pente_max,zm,mv,iq)
  ! call minmaxq(zq,qmin,qmax,'apres vly     ')


  call vlx(zq,pente_max,zm,mu,iq)
  ! call minmaxq(zq,qmin,qmax,'apres vlx2    ')


  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

  RETURN
END SUBROUTINE vlsplt
RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi
        min_qParent,min_qMass,min_ratio ! MVals et CRisi

  ! 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 dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
  USE iniprint_mod_h
IMPLICIT NONE
  !


  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
  REAL :: u_m( ip1jmp1,llm )
  REAL :: q(ip1jmp1,llm,nqtot)
  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 sigu(ip1jmp1)
  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
  DATA first/.true./

  !   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
  !print*,'Ok calcul des pentes'

  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

  ! print*,'Bouclage en iip1'

  !   calcul des flux a gauche et a droite

#ifdef CRAY

  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.
  !print*,'Cumule ....'

  DO l=1,llm
   DO ij=iip2,ip1jm-1
  ! print*,'masse(',ij,')=',masse(ij,l,iq)
      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)*(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
#endif

  ! go to 9999
  !   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
  !print*,'Ok test 1'
  DO l=1,llm
   DO ij=iip1+iip1,ip1jm,iip1
      iadvplus(ij,l)=iadvplus(ij-iim,l)
   ENDDO
  ENDDO
  ! 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
  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
  if (prt_level > 2) PRINT *, &
        'Nombre de points pour lesquels on advect plus que le' &
        ,'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
  !9999    continue


  !   bouclage en latitude
  !print*,'cvant 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(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen

  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)
        ! !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
        if (q(ij,l,iq).gt.min_qParent) then
          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
        else
          Ratio(ij,l,iq2)=min_ratio
        endif
      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
        ! !MVals: veiller a ce qu'on 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=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 ! DO ij=ijb+iip1-1,ije,iip1
    enddo !DO l=1,llm
  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 vlx
RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
  USE comgeom_mod_h
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi
        min_qParent,min_qMass,min_ratio ! MVals et 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 ....
  ! dq         sont des arguments 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)
  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 appn apps
  ! REAL newq,oldmasse
  LOGICAL :: first
  SAVE 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
  SAVE airej2,airejjm

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

  !
  !
  REAL :: SSUM

  DATA first/.true./

  ! !write(*,*) 'vly 578: entree, iq=',iq

  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'

  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

  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
  DO ij=1,iip1
     dyq(ij,l)=0.
     dyq(ip1jm+ij,l)=0.
  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

  ! !write(*,*) 'vly 756'
  DO l=1,llm
   DO ij=1,ip1jm
      IF(masse_adv_v(ij,l).gt.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

  ! 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,nqDesc(iq)=',iq,tracers(iq)%nqDescen

  do ifils=1,tracers(iq)%nqDescen
    iq2=tracers(iq)%iqDescen(ifils)
    DO l=1,llm
      DO ij=1,ip1jmp1
        ! ! attention, chaque fils doit avoir son masseq, sinon, le 1er
        ! ! fils ecrase le masseq de ses freres.
        ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
  !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
        if (q(ij,l,iq).gt.min_qParent) then
          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
        else
          Ratio(ij,l,iq2)=min_ratio
        endif
      enddo
    enddo
  enddo

  do ifils=1,tracers(iq)%nqDescen
    iq2=tracers(iq)%iqDescen(ifils)
    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

     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

     ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
     ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols

     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

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

  ! 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 853: sortie'

  RETURN
END SUBROUTINE vly
RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi
        min_qParent,min_qMass,min_ratio ! MVals et CRisi
  !
  ! 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 dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !


  !
  !
  !   Arguments:
  !   ----------
  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
  REAL :: q(ip1jmp1,llm,nqtot)
  REAL :: w(ip1jmp1,llm+1)
  INTEGER :: iq
  !
  !  Local
  !   ---------
  !
  INTEGER :: ij,l
  !
  REAL :: wq(ip1jmp1,llm+1),newmasse

  REAL :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
  REAL :: sigw

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

  LOGICAL :: testcpu
  SAVE testcpu

#ifdef BIDON
  REAL :: temps0,temps1,second
  SAVE temps0,temps1

  DATA testcpu/.false./
  DATA temps0,temps1/0.,0./
#endif

  !    On oriente tout dans le sens de la pression c'est a dire dans le
  !    sens de W

  ! !write(*,*) 'vlz 923: entree'

#ifdef BIDON
  IF(testcpu) THEN
     temps0=second(0.)
  ENDIF
#endif
  DO l=2,llm
     DO ij=1,ip1jmp1
        dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
        adzqw(ij,l)=abs(dzqw(ij,l))
     ENDDO
  ENDDO

  DO l=2,llm-1
     DO ij=1,ip1jmp1
#ifdef CRAY
        dzq(ij,l)=0.5* &
              cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
#else
        IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
            dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
        ELSE
            dzq(ij,l)=0.
        ENDIF
#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

  ! !write(*,*) 'vlz 954'
  DO ij=1,ip1jmp1
     dzq(ij,1)=0.
     dzq(ij,llm)=0.
  ENDDO

#ifdef BIDON
  IF(testcpu) THEN
     temps1=temps1+second(0.)-temps0
  ENDIF
#endif
  ! ---------------------------------------------------------------
  !   .... calcul des termes d'advection verticale  .......
  ! ---------------------------------------------------------------

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

   ! !write(*,*) 'vlz 969'
   DO l = 1,llm-1
     do  ij = 1,ip1jmp1
      IF(w(ij,l+1).gt.0.) THEN
         sigw=w(ij,l+1)/masse(ij,l+1,iq)
         wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq) &
               +0.5*(1.-sigw)*dzq(ij,l+1))
      ELSE
         sigw=w(ij,l+1)/masse(ij,l,iq)
         wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
      ENDIF
     ENDDO
   ENDDO

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

  ! 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,nqChildren(iq)
  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)
        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
        if (q(ij,l,iq).gt.min_qParent) then
          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
        else
          Ratio(ij,l,iq2)=min_ratio
        endif
      enddo
    enddo
  enddo

  do ifils=1,tracers(iq)%nqChildren
    iq2=tracers(iq)%iqDescen(ifils)
    call vlz(Ratio,pente_max,masseq,wq,iq2)
  enddo
  ! end CRisi

  DO l=1,llm
     DO ij=1,ip1jmp1
        newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
        q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l)) &
              /newmasse
        masse(ij,l,iq)=newmasse
     ENDDO
  ENDDO

  ! 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(*,*) 'vlsplt 1032'

  RETURN
END SUBROUTINE vlz
 ! SUBROUTINE minmaxq(zq,qmin,qmax,comment)
!
!INCLUDE "dimensions_mod.f90"
!INCLUDE "paramet.h"

!  CHARACTER*(*) comment
!  real qmin,qmax
!  real zq(ip1jmp1,llm)

!  INTEGER jadrs(ip1jmp1), jbad, k, i


!  DO k = 1, llm
!     jbad = 0
!     DO i = 1, ip1jmp1
!     IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
!        jbad = jbad + 1
!        jadrs(jbad) = i
!     ENDIF
!     ENDDO
!     IF (jbad.GT.0) THEN
!     PRINT*, comment
!     DO i = 1, jbad
!c            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
!     ENDDO
!     ENDIF
!  ENDDO

!  return
!  end
subroutine minmaxq(zq,qmin,qmax,comment)
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
  USE paramet_mod_h

  character(len=20) :: comment
  real :: qmin,qmax
  real :: zq(ip1jmp1,llm)
  real :: zzq(iip1,jjp1,llm)

#ifdef isminmax
  integer :: imin,jmin,lmin,ijlmin
  integer :: imax,jmax,lmax,ijlmax

  integer :: ismin,ismax

  call scopy (ip1jmp1*llm,zq,1,zzq,1)

  ijlmin=ismin(ijp1llm,zq,1)
  lmin=(ijlmin-1)/ip1jmp1+1
  ijlmin=ijlmin-(lmin-1.)*ip1jmp1
  jmin=(ijlmin-1)/iip1+1
  imin=ijlmin-(jmin-1.)*iip1
  zqmin=zq(ijlmin,lmin)

  ijlmax=ismax(ijp1llm,zq,1)
  lmax=(ijlmax-1)/ip1jmp1+1
  ijlmax=ijlmax-(lmax-1.)*ip1jmp1
  jmax=(ijlmax-1)/iip1+1
  imax=ijlmax-(jmax-1.)*iip1
  zqmax=zq(ijlmax,lmax)

   if(zqmin.lt.qmin) &
  ! s     write(*,9999) comment,
         write(*,*) comment, &
         imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
   if(zqmax.gt.qmax) &
  ! s     write(*,9999) comment,
         write(*,*) comment, &
         imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)

#endif
  return
  !9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
end subroutine minmaxq