source: LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F90 @ 5322

!
! $Id: vlsplt_loc.F90 5285 2024-10-28 13:33:29Z abarral $
!
RECURSIVE 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 iniprint_mod_h
  USE parallel_lmdz
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi                 &
        min_qParent,min_qMass,min_ratio ! MVals et CRisi
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !


  !
  !
  !   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 :: SSUM
  EXTERNAL  SSUM

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

#ifdef CRAY
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l=1,llm
   DO ij=ijb,ije-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,iq))
      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
!$OMP END DO NOWAIT
#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 ....'
!$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).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
!$OMP END DO NOWAIT
#endif

  ! go to 9999
  !   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).lt.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).gt.0) THEN
           iju=0
  !   indicage des mailles concernees par le traitement special
           DO ij=ijb,ije
              IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.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.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
!$OMP END DO NOWAIT
  !ym      ENDIF  ! n0.gt.0
9999   continue

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

  ! !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x
  ! 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_loc


RECURSIVE 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 comgeom_mod_h
  USE parallel_lmdz
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi                 &
        min_qParent,min_qMass,min_ratio ! MVals et CRisi
  USE comconst_mod, ONLY: pi
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE
  !


  !
  !
  !   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,testcpu
  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,testcpu
!$OMP THREADPRIVATE(first,testcpu)

  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
  !
  !
  REAL :: SSUM
  EXTERNAL  SSUM

  DATA first,testcpu/.true.,.false./
  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).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
!$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).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
!$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).gt.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


  !.-. ancienne version
     ! convpn=SSUM(iim,qbyv(1,l),1)/apoln
     ! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
     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

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

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


  RETURN
END SUBROUTINE vly_loc



RECURSIVE 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 iniprint_mod_h
  USE parallel_lmdz
  USE vlz_mod
  USE infotrac, ONLY : nqtot,tracers, & ! CRisi                 &
        min_qParent,min_qMass,min_ratio ! MVals et CRisi

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


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

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

  REAL :: SSUM
  EXTERNAL  SSUM

  DATA testcpu/.false./
  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
#ifdef BIDON
  IF(testcpu) THEN
     temps0=second(0.)
  ENDIF
#endif

  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
#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
!$OMP END DO NOWAIT

!$OMP MASTER
  DO ij=ijb,ije
     dzq(ij,1)=0.
     dzq(ij,llm)=0.
  ENDDO
!$OMP END MASTER
!$OMP BARRIER
#ifdef BIDON
  IF(testcpu) THEN
     temps1=temps1+second(0.)-temps0
  ENDIF
#endif

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

  RETURN
END SUBROUTINE vlz_loc
 ! 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