#ifdef ISO
#ifdef ISOTRAC
! $Id: $

MODULE isotrac_routines_mod
! on créé ce module pour éviter dépendances circulaires.
! isotopes_routines_mod a besoin de isotrac et isotopes_verif
! isotopes_verif a besoin de isotopes et isotrac
! isotrac n'a besoin que de isotopes
    USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, index_trac=>itZonIso, ntraceurs_zone=>nzone
IMPLICIT NONE

CONTAINS


         SUBROUTINE uncompress_commun_zone(ncas, cas, &
         xtp_cas,xtp,xtwater_cas,xtwater,xtevap_cas,xtevap, &
                 ncum,izone)

    USE isotopes_mod, ONLY: ridicule,iso_eau

         implicit none 

         ! decompression des outputs communs à tous les cas dans
         ! appel_stewart
         ! cas des traceurs. Ici, aucun risque de revap franche.

         ! inputs
         integer ncas,ncum
         integer cas(ncum)
         real xtevap_cas(niso,ncum)
         real xtp_cas(niso,ncum)
         real xtwater_cas(niso,ncum)
         integer izone
         ! outputs
         real xtwater(ntraciso,ncum)
         real xtp(ntraciso,ncum)
         real xtevap(ntraciso,ncum)

         ! locals
         integer il,ixt,iiso,ixt_revap


         do il=1,ncas
             do iiso=1,niso 
               ixt=index_trac(izone,iiso) 
               xtevap(ixt,cas(il))=xtevap_cas(iiso,il)
               xtp(ixt,cas(il))=xtp_cas(iiso,il)
               xtwater(ixt,cas(il))=xtwater_cas(iiso,il)
           enddo !do iiso=1,niso
         enddo !do il=1,ncas

         END SUBROUTINE  uncompress_commun_zone


         SUBROUTINE uncompress_commun_zone_revap(ncas, cas, &
                xtp_cas,xtp,xtwater_cas,xtwater, &
                xtevap_cas,xtevap, &
                ncum,izone,Eqi_stewart,fac_ftmr_cas, &
#ifdef ISOVERIF
                 Exi_cas,Exi,    &
#endif
                xtp_avantevap_cas,liq,hdiag)

    USE isotopes_mod, ONLY: ridicule,iso_eau,iso_HDO,ridicule_evap
    USE isotrac_mod, only: option_revap,evap_franche,izone_revap, &
&       ridicule_trac
#ifdef ISOVERIF
    USE isotopes_verif_mod
#endif

         implicit none 

         ! decompression des outputs communs à tous les cas dans
         ! appel_stewart
         ! cas des traceurs: mais ici, risque de révap franche -> on fat
         ! plus attention

         ! inputs
         integer ncas,ncum
         integer cas(ncum)
         real xtevap_cas(niso,ncum)
         real xtp_cas(niso,ncum)
         real xtwater_cas(niso,ncum)
         integer izone
         real Eqi_stewart(ncum)
         real xtp_avantevap_cas(niso,ncum)
         real fac_ftmr_cas(ncum) 
         integer liq  
         real hdiag(ncas)      
         
         ! outputs
         real xtwater(ntraciso,ncum)
         real xtp(ntraciso,ncum)
         real xtevap(ntraciso,ncum)

         ! locals
         integer il,ixt,iiso,ixt_revap
         real xtaddp_tag(niso,ncum)
#ifdef ISOVERIF
         integer ieau,iHDO
         real Exi_cas(niso,ncum)
         real Exi(ntraciso,ncum)
         !USE isotopes_verif, ONLY: 
#endif         

!        write(*,*) 'compress_stewart 315 tmp: ',
!     :           'entrée dans uncompress_commun_zone_revap'

         do il=1,ncas       
          if ((option_revap.eq.1).and. &
              ((((liq.eq.1).and. &
              (Eqi_stewart(il)*fac_ftmr_cas(il).gt.evap_franche) &
              .and.(hdiag(il).lt.0.99)).or. &
              ((liq.eq.0).and. &
              (Eqi_stewart(il)*fac_ftmr_cas(il).ge.0.0))).or. &
              (izone.eq.izone_revap))) then
!          if ((option_revap.eq.1).and.
!     :           (((Eqi_stewart(il)*fac_ftmr_cas(il).gt.evap_franche)
!     :           .or.((liq.eq.0).and.
!     :           (Eqi_stewart(il)*fac_ftmr_cas(il).ge.0.0))).or.
!     :           (izone.eq.izone_revap))) then
            ! on met la revap dans izone_revap si option_revap=1 et si:
            ! * evap glace (non fractionnante)
            ! * ou evap liq suffisemment forte pour que pas de flux
            ! d'isotopes négatifs.
            ! si option_revap=1 et izone=izone_revap, on met aussi dans izone_revap
!#ifdef ISOVERIF            
!            if (il.eq.1) then
!            write(*,*) 'compress tmp 341: revap dans izone_revap'
!            write(*,*) 'Eqi_stewart(il),fac_ftmr_cas(il),evap_franche=',
!     :           Eqi_stewart(il),fac_ftmr_cas(il),evap_franche
!            write(*,*) 'il,xtp_cas(iso_eau,il)=',il,xtp_cas(iso_eau,il)
!            write(*,*) 'il,xtp_avantevap_cas(iso_eau,il)=',il,
!     :            xtp_avantevap_cas(iso_eau,il)
!            write(*,*) 'xtp(ixt_revap,cas(il))=',
!     :            xtp(index_trac(izone_revap,iso_eau),cas(il))
!            endif
!#endif           

            ! toute la révap franche va dans izone_revap
            do iiso=1,niso
               ixt=index_trac(izone,iiso)           
               ixt_revap=index_trac(izone_revap,iiso)  
               ! le terme d'évap qui était pour la zone izone devient
               ! nul, et à la place on le met dans izone_revap&               
               xtevap(ixt_revap,cas(il))=xtevap(ixt_revap,cas(il)) &
                         +xtevap_cas(iiso,il)

               ! ce qui a été ajouté à xtp par rapport à xtp_avantevap
               ! est ajouté à izone_revap, au lieu de izone
               xtaddp_tag(iiso,il)=xtp_cas(iiso,il) &
                         -xtp_avantevap_cas(iiso,il)
               xtp(ixt_revap,cas(il))= &
                        xtp(ixt_revap,cas(il)) &
                        +xtaddp_tag(iiso,il)
            enddo !do iiso=1,niso
#ifdef ISOVERIF
            if (iso_HDO.gt.0) then
                if (xtevap_cas(iso_eau,il).gt.ridicule_evap) then
                    CALL iso_verif_aberrant_choix( &
                       xtevap_cas(iso_HDO,il),xtevap_cas(iso_eau,il), &
                       ridicule_trac,deltalimtrac*2,'compress 344a')
                endif
                ieau=index_trac(izone_revap,iso_eau)
                iHDO=index_trac(izone_revap,iso_HDO) 
                CALL iso_verif_aberrant_choix(xtevap(iHDO,cas(il)), &
                    xtevap(ieau,cas(il)),ridicule_trac,deltalimtrac*2, &
                    'compress 344b')
            endif !if (iso_HDO.gt.0) then
#endif            

            ! l'évap des autres zones devient nulle
            if (izone.ne.izone_revap) then
                do iiso=1,niso
                  ixt=index_trac(izone,iiso)
                  xtevap(ixt,cas(il))=0.0
                  xtp(ixt,cas(il))=xtp_avantevap_cas(iiso,il)
                enddo
            endif

!#ifdef ISOVERIF            
!            if (il.eq.1) then
!            write(*,*) 'compress tmp 341: revap dans izone_revap'
!            write(*,*) 'xtp(ixt_revap,cas(il))=',
!     :           xtp(index_trac(izone_revap,iso_eau),cas(il))
!            write(*,*) 'xtp(ixt,cas(il))=',
!     :            xtp(index_trac(izone,iso_eau),cas(il))
!            endif
!#endif            

           else !if ((Eqi_stewart(il).gt.ridicule_evap*100)

             do iiso=1,niso 
               ixt=index_trac(izone,iiso) 
               xtevap(ixt,cas(il))=xtevap_cas(iiso,il)
               xtp(ixt,cas(il))=xtp_cas(iiso,il)
             enddo !do iiso=1,niso
           endif !if ((Eqi_stewart(il).gt.ridicule_evap*100)
        enddo !do il=1,ncas

        do il=1,ncas
           do iiso=1,niso
             ixt=index_trac(izone,iiso)
             xtwater(ixt,cas(il))=xtwater_cas(iiso,il)
#ifdef ISOVERIF
             Exi(ixt,cas(il))=Exi_cas(iiso,il)
#endif             
           enddo !do iiso=1,niso
         enddo !do il=1,ncas
!         write(*,*) 'compress_stewart 361 tmp: ',
!     :           'sortie de uncompress_commun_zone_revap'
         
         END SUBROUTINE  uncompress_commun_zone_revap



         SUBROUTINE compress_cond_facftmr_zone( &
          ncas,  cas, &
          Eqi_prime_cas,Eqi_prime, &
          Pqisup_cas,Pqisup,  &
          Pxtisup_cas,Pxtisup,   &
          qp_avantevap_cas,qp_avantevap, &
          xtp_avantevap_cas,xtp_avantevap,  &
          xtevapsup_cas,xtevap,&
          water_cas,water,&
#ifdef ISOVERIF        
          evap_cas,evap, &
#endif        
          nloc,ncum,nd,i,izone)

    USE isotopes_mod, ONLY: iso_eau
#ifdef ISOVERIF        
    USE isotopes_verif_mod  
#endif

         implicit none

         ! compression dans le cas condensation_facftmr
         ! inputs
         integer nd,ncum,nloc
         integer ncas
         integer cas(ncum)
         integer i
         integer izone
         real  &
          xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
          water_cas(ncum),water(ncum)
         real Eqi_prime_cas(ncum),Eqi_prime(ncum), &
          Pqisup_cas(ncum),Pqisup(ncum),  &
          Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum),  &
          qp_avantevap_cas(ncum),qp_avantevap(ncum), &
          xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
#ifdef ISOVERIF      
         real evap_cas(ncum),evap(ncum)
#endif         
         ! locals
         integer il,ixt,ieau,iiso

          ieau=index_trac(izone,iso_eau)
          do il=1,ncas
            if (qp_avantevap(cas(il)).gt.0.0) then
               Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
              *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
            else !if (qp_avantevap_cas(cas(il)).gt.0.0) then
#ifdef ISOVERIF
                CALL iso_verif_egalite_choix( &
                         (Eqi_prime(cas(il))),0.0, &
                         'compress_stewart 495',errmax,errmaxrel)
#endif
                Eqi_prime_cas(il)=0.0
            endif !if (qp_avantevap_cas(cas(il)).gt.0.0) thens
            
            if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
               water_cas(il)=water(cas(il)) &
                 *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
                 /(Pqisup(cas(il))-Eqi_prime(cas(il))))
            else !if (Pqisup(cas(il)).gt.0.0) then
#ifdef ISOVERIF
                CALL iso_verif_egalite_choix(water(cas(il)),0.0, &
                         'compress_stewart 507',errmax,errmaxrel)
#endif
                water_cas(il)=0.0
            endif !if (Pqisup(cas(il)).gt.0.0) then

            Pqisup_cas(il)=Pxtisup(ieau,cas(il))
            qp_avantevap_cas(il)=xtp_avantevap(ieau,cas(il))

!#ifdef ISOVERIF
!            CALL iso_verif_noNaN(water_cas(il),'compress_stewart 518')
!            evap_cas(il)=evap(cas(il))
!     :        *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
!!           qp_cas(il)=xtp(ieau,cas(il))     
!#endif
            do iiso=1,niso
              ixt=index_trac(izone,iiso)
              Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
              xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
              xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
            enddo
          enddo 

         END SUBROUTINE  compress_cond_facftmr_zone

         SUBROUTINE compress_cond_nofftmr_zone( &
          ncas,  cas, &
          Eqi_prime_cas,Eqi_prime, &
          Pqisup_cas,Pqisup,  &
          Pxtisup_cas,Pxtisup, &
          water_cas,water,   &
          qp_avantevap_cas,qp_avantevap, &
          xtp_avantevap_cas,xtp_avantevap, &
          xt_cas,xt,q_cas,q,  &
          xtevapsup_cas,xtevap, &
#ifdef ISOVERIF
          evap_cas,evap, &
#endif      
          nloc,ncum,nd,i,izone)

    USE isotopes_mod, ONLY: iso_eau
#ifdef ISOVERIF
    USE isotopes_verif_mod
#endif

         implicit none

         ! compression dans le cas condensation_facftmr
         ! inputs
         integer nloc,nd,ncum
         integer ncas
         integer cas(ncum)
         integer i
         integer izone
         real  &
          xt_cas(niso,ncum),q_cas(ncum),xt(ntraciso,ncum),q(ncum),  &
          xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
          water_cas(ncum),water(ncum)
         real Eqi_prime_cas(ncum),Eqi_prime(ncum), &
          Pqisup_cas(ncum),Pqisup(ncum),  &
          Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum), &
          qp_avantevap_cas(ncum),qp_avantevap(ncum), &
          xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
#ifdef ISOVERIF
         real evap_cas(ncum),evap(ncum)
#endif         
         ! locals
          integer il,ixt,ieau,iiso

          ieau=index_trac(izone,iso_eau)
          do il=1,ncas
            if (qp_avantevap(cas(il)).gt.0) then
              Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
                 *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
            else
                Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
                 *(xt(ieau,cas(il))/q(cas(il)))
            endif
            
            if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
              water_cas(il)=water(cas(il)) &
                 *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
                 /(Pqisup(cas(il))-Eqi_prime(cas(il))))
            else !if (Pqisup(cas(il)).gt.0.0) then
#ifdef ISOVERIF
                CALL iso_verif_egalite_choix(water(cas(il)),0.0, &
                         'compress_stewart 654',errmax,errmaxrel)
#endif
                water_cas(il)=0.0
            endif !if (Pqisup(cas(il)).gt.0.0) then  

            Pqisup_cas(il)=Pxtisup(ieau,cas(il))
            qp_avantevap_cas(il)=xtp_avantevap(ieau,cas(il))
            q_cas(il)=xt(ieau,cas(il))
!#ifdef ISOVERIF
!            if (qp_avantevap(cas(il)).gt.0.0) then
!              evap_cas(il)=evap(cas(il))
!     :           *(xtp_avantevap(ieau,cas(il))/qp_avantevap(cas(il)))
!            else
!                evap_cas(il)=evap(cas(il))
!     :           *(xt(ieau,cas(il))/q(cas(il)))
!            endif
!            qp_cas(il)=xtp(ieau,cas(il))
!#endif            
            do iiso=1,niso 
              ixt=index_trac(izone,iiso)            
              Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
              xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
              xt_cas(iiso,il)=xt(ixt,cas(il))
              xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
            enddo
          enddo 

         END SUBROUTINE  compress_cond_nofftmr_zone

         SUBROUTINE compress_noevap_zone( &
          ncas,  cas, &
          Pqisup_cas,Pqisup,  &
          Pxtisup_cas,Pxtisup,   &
          xtp_avantevap_cas,xtp_avantevap, &
          xtevapsup_cas,xtevap, &
          water_cas,water, &
#ifdef ISOVERIF         
          evap_cas,evap, &
#endif 
          nloc,ncum,nd,i,izone)

    USE isotopes_mod, ONLY: ridicule,iso_eau
#ifdef ISOVERIF
    USE isotopes_verif_mod
#endif
         implicit none

         ! compression dans le cas condensation_facftmr
         integer nloc,nd,ncum
         integer ncas
         integer cas(ncum)
         integer i
         integer izone
         real xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum), &
          water_cas(ncum),water(ncum)
         real Pqisup_cas(ncum),Pqisup(ncum),  &
          Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum),   &
          xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum)
#ifdef ISOVERIF        
        real evap_cas(ncum),evap(ncum)
#endif 
         integer il,ixt,ieau,iiso

          ieau=index_trac(izone,iso_eau)
          do il=1,ncas
            Pqisup_cas(il)=Pxtisup(ieau,cas(il))
            if (Pqisup(cas(il)).gt.0.0) then
               water_cas(il)=water(cas(il)) &
                 *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
            else
                water_cas(il)=0.0
#ifdef ISOVERIF                 
                CALL iso_verif_egalite_choix(water(cas(il)), &
                     0.0,'compress_stewart 709',errmax,errmaxrel)
#endif                
            endif
#ifdef ISOVERIF 
!            evap_cas(il)=evap(cas(il))
!     :           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
!            qp_cas(il)=xtp(ieau,cas(il))
#endif
            do iiso=1,niso
              ixt=index_trac(izone,iiso)            
              Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
              xtp_avantevap_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
              xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
            enddo
          enddo 

         END SUBROUTINE  compress_noevap_zone

         SUBROUTINE compress_evap_liq_zone(iflag_con,ncas,   &
          cas, &
          Pqisup_cas,Pqisup,  &
          Pxtisup_cas,Pxtisup,     &
          xtp_avantevap_cas,xtp_avantevap, &
          xtp_avantevaptrac_cas,qp_avantevaptrac_cas, &
          xtevapsup_cas,xtevap, &
          water_cas,water, &
          Eqi_stewart,Pqiinf_stewart,Eqi_prime_cas,&
          Pqiinf,Eqi_par,Pqiinf_par,Eqi_prime,ptrac, &
          Eqi,Eqi_cas,  &
!     &    qp_cas,
#ifdef ISOVERIF           
          evap_cas,evap,   &
#endif         
          nloc,ncum,nd,izone)

    USE isotopes_mod, ONLY: ridicule,iso_eau
#ifdef ISOVERIF
    USE isotopes_verif_mod
#endif
         implicit none

         ! compression dans le cas condensation_facftmr
        ! inputs et outputs 
         integer iflag_con
         integer izone    
         integer nloc,nd,ncum
         integer ncas
         integer cas(ncum)
!         integer i
         real xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum)
         real xtp_avantevap_cas(niso,ncum), &
                 xtp_avantevap(ntraciso,ncum)
         real water_cas(ncum),water(ncum)
         real xtp_avantevaptrac_cas(niso,ncum), &
                 qp_avantevaptrac_cas(ncum)
         real ptrac(ncum)
!         real qp_cas(ncum)
#ifdef ISOVERIF           
         real evap_cas(ncum),evap(ncum)
#endif         
         real &
          Eqi_stewart(ncum),Pqiinf_stewart(ncum),Eqi_prime_cas(ncum), &
          Pqiinf(ncum),Eqi_par(ncum),Pqiinf_par(ncum), &
          Eqi_prime(ncum),Pqisup(ncum),Pqisup_cas(ncum), &
          Pxtisup(ntraciso,ncum),Pxtisup_cas(niso,ncum), &
          Eqi(ncum),Eqi_cas(ncum)
         ! locals
          integer il,ixt,iiso,ieau

!          write(*,*) 'compress 910: xtp_avantevap(iso_eau,cas(1))=',
!     :           xtp_avantevap(iso_eau,cas(1))
!        write(*,*) 'compress_evap_liq_zone 510: ncas,ncum=',ncas,ncum
        ptrac(:)=0. ! CR 31 mars 2023: initialisation de ptrac

          ieau=index_trac(izone,iso_eau)
          do il=1,ncas            
            Pqisup_cas(il)=Pxtisup(ieau,cas(il))

            if (Pqisup(cas(il)).gt.0.0) then
              ptrac(il)=Pxtisup(ieau,cas(il))/Pqisup(cas(il))
              Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
                 *ptrac(il)
            else
#ifdef ISOVERIF 
              CALL iso_verif_egalite(( &
                         Eqi_prime(cas(il))),0.0, &
                         'compress_stewart 979')
#endif                
              Eqi_prime_cas(il)=0.0 
            endif

            if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
               water_cas(il)=water(cas(il)) &
                 *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
                 /(Pqisup(cas(il))-Eqi_prime(cas(il))))
            else !if (Pqisup(cas(il)).gt.0.0) then
#ifdef ISOVERIF
                CALL iso_verif_egalite_choix(water(cas(il)),0.0, &
                         'compress_stewart 507',errmax,errmaxrel)
#endif
                water_cas(il)=0.0
            endif !if (Pqisup(cas(il)).gt.0.0) then

            
!            qp_cas(il)=qp(cas(il))
#ifdef ISOVERIF              
!            evap_cas(il)=evap(cas(il))
!     &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
!            ! ce calcul est faux& normalement, 
!        evap_cas(il)=Eqi_prime_cas(il)//100.0/delP_cas(il)/SIGD*g*2.0
!     &         -(ieau,cas(il))
#endif            
            qp_avantevaptrac_cas(il)=xtp_avantevap(ieau,cas(il))
            do iiso=1,niso  
              ixt=index_trac(izone,iiso)            
              Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
              xtp_avantevap_cas(iiso,il)=xtp_avantevap(iiso,cas(il))
              xtp_avantevaptrac_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
              xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
            enddo
          enddo !do il=1,ncas
!         write(*,*) 'compress_stewart 931: ',
!     &           'xtp_avantevap_cas(iso_eau,1)=', 
!     &           xtp_avantevap_cas(iso_eau,1)
!         write(*,*) 'xtp_avantevap(iso_eau,cas(1))=',
!     &           xtp_avantevap(iso_eau,cas(1))

      ! calculs des flux de masses à mettre en argument de stewart:
      ! comme l'eau n'est pas bien concervée dans les ddfts, on est
      ! obligé de bidouillé.
      ! 1) soit on considère Pqisup, Eqi, et Pqiinf_par=Pqisup-Eqi
      !    et on suppose que dans la réalité les compositions de
      !    Pqiinf sont les même que Pqiinf_par
      ! 2) soit on considère Pqisup, Eqi_par=Pqisup-Pqiinf, et Pqiinf,
      !    et on suppose que dans la réalité les compositions de
      !    Eqi_prime sont les même que Eqi_par
          do il=1,ncas  
           if (Pqisup(cas(il)).gt.0.0) then            
            if ((water(cas(il)).gt.ridicule/100).and. &
                  (Pqiinf_par(cas(il)).le.0.0)) then
             ! on ne peut pas utiliser la méthode 1, car KE prédit de l'eau
             ! alors que le bilan de masse n'enprédit pas.
             ! Peut-on utiliser la méthode 2?
             Pqiinf_stewart(il)=Pqiinf(cas(il))*ptrac(il)
             Eqi_stewart(il)=Eqi_par(cas(il)) *ptrac(il)
           else !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
             ! il n'y a pas d'obstacles à l'utilisation de 1)
             Pqiinf_stewart(il)=Pqiinf_par(cas(il))*ptrac(il)
             if (iflag_con.eq.30) then
                Eqi_stewart(il)=Eqi_prime(cas(il))*ptrac(il)
             else
                Eqi_stewart(il)=Eqi(cas(il))*ptrac(il)
             endif
           endif !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
          else ! if (Pqisup(cas(il).gt.0.0) then
#ifdef ISOVERIF
               CALL iso_verif_egalite((Pqiinf(cas(il))), &
                 0.0,'compress_stewart 1047a')
               CALL iso_verif_egalite(( &
                 Eqi_prime(cas(il))),0.0,'compress_stewart 1047b')
               CALL iso_verif_egalite(( &
                 Pqiinf_par(cas(il))),0.0,'compress_stewart 1047c')
               CALL iso_verif_egalite((Eqi_par(cas(il))), &
                 0.0,'compress_stewart 1047d')
#endif              
               Pqiinf_stewart(il)=0.0
               Eqi_stewart(il)=0.0
          endif ! if (Pqisup(cas(il).gt.0.0) then 
         enddo !do il=1,ncas_evap_glace

         ! petite vérif
#ifdef ISOVERIF        
        do il=1,ncas
          if ((abs(water_cas(il)).ge.ridicule/10.) &
                 .and.(Pqiinf_stewart(il).le.0.0)) then
              write(*,*) 'compress_stewart 498: evap liq:'
              write(*,*) 'water(il,i)=', water_cas(il)
              write(*,*) 'Pqiinf=',Pqiinf(cas(il))
              write(*,*) 'Pqiinf_par=',Pqiinf_par(cas(il))
              write(*,*) 'Pqiinf_stewart=',Pqiinf_stewart(il)
              stop                   
          endif
        enddo !do il=1,ncas_evap_glace
#endif

         END SUBROUTINE  compress_evap_liq_zone

         SUBROUTINE compress_evap_glace_zone(iflag_con, &
          ncas, cas, &
          water_cas,water,     &
          Pqisup_cas,Pqisup,  &
          Pxtisup_cas,Pxtisup,  &
          xtp_avantevap_cas,xtp_avantevap,  &
          xtp_avantevaptrac_cas,qp_avantevaptrac_cas, &
          xtevapsup_cas,xtevap, &
          Eqi_stewart,Pqiinf_stewart,Eqi_prime_cas,Eqi_cas, &
          Pqiinf,Eqi_par,Pqiinf_par,Eqi_prime,Eqi, &
!     &    qp_cas,
#ifdef ISOVERIF           
          evap_cas,evap, &
#endif         
          nloc,ncum,nd,i,frac_sublim,izone)

    USE isotopes_mod, ONLY: ridicule,iso_eau
#ifdef ISOVERIF
    USE isotopes_verif_mod
#endif
         implicit none

         ! compression dans le cas condensation_facftmr
         integer iflag_con         
         integer nloc,nd,ncum
         integer ncas
         integer cas(ncum)
         integer i
         integer izone
         real  &
          water_cas(ncum),water(ncum), &
          xtevapsup_cas(niso,ncum),xtevap(ntraciso,ncum)
!         real qp_cas(ncum)    
#ifdef ISOVERIF           
         real evap_cas(ncum),evap(ncum)
#endif         
         real  &
          Pqisup_cas(ncum),Pqisup(ncum),  &
          Pxtisup_cas(niso,ncum),Pxtisup(ntraciso,ncum),   &
          xtp_avantevap_cas(niso,ncum),xtp_avantevap(ntraciso,ncum), &
          Eqi_stewart(ncum),Pqiinf_stewart(ncum),Eqi_prime_cas(ncum), &
          Pqiinf(ncum),Eqi_par(ncum),Pqiinf_par(ncum),Eqi_prime(ncum), &
          Eqi(ncum),Eqi_cas(ncum)

         real xtp_avantevaptrac_cas(niso,ncum), &
                 qp_avantevaptrac_cas(ncum)
          integer frac_sublim
          ! locals
          integer il,ixt,ieau,iiso

          ieau=index_trac(izone,iso_eau)
          do il=1,ncas
            Pqisup_cas(il)=Pxtisup(ieau,cas(il))  

            if (Pqisup(cas(il)).gt.0.0) then
              Eqi_prime_cas(il)=Eqi_prime(cas(il)) &
                 *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
              Eqi_cas(il)=Eqi(cas(il)) & ! corr bug Camille 15 juin 2024
                 *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
            else
#ifdef ISOVERIF 
                CALL iso_verif_egalite(( &
                         Eqi_prime(cas(il))),0.0, &
                         'compress_stewart 979b')
#endif                
              Eqi_prime_cas(il)=0.0 
              Eqi_cas(il)=0.0
            endif

            if (Pqisup(cas(il))-Eqi_prime(cas(il)).gt.0.0) then
               water_cas(il)=water(cas(il)) &
                 *((Pxtisup(ieau,cas(il))-Eqi_prime_cas(il)) &
                 /(Pqisup(cas(il))-Eqi_prime(cas(il))))
            else !if (Pqisup(cas(il)).gt.0.0) then
#ifdef ISOVERIF
                CALL iso_verif_egalite_choix(water(cas(il)),0.0, &
                         'compress_stewart 507',errmax,errmaxrel)
#endif
                water_cas(il)=0.0
            endif !if (Pqisup(cas(il)).gt.0.0) then


            qp_avantevaptrac_cas(il)=xtp_avantevap(ieau,cas(il))
!            qp_cas(il)=xtp(ieau,cas(il))
#ifdef ISOVERIF              
!            evap_cas(il)=evap(cas(il))
!     &           *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
            ! ce calcul est faux& faire plutot:
!            evap_cas(il)=Eqi_prime_cas(il)//100.0/delP_cas(il)/SIGD*g*2.0
!     &         -(ieau,cas(il))
#endif            
            do iiso=1,niso  
              ixt=index_trac(izone,iiso)            
              Pxtisup_cas(iiso,il)=Pxtisup(ixt,cas(il))
              xtp_avantevap_cas(iiso,il)=xtp_avantevap(iiso,cas(il))
              xtp_avantevaptrac_cas(iiso,il)=xtp_avantevap(ixt,cas(il))
              xtevapsup_cas(iiso,il)=xtevap(ixt,cas(il))
            enddo
          enddo  !do il=1,ncas     

          ! calculs des flux de masses à mettre en argument de stewart:
      ! comme l'eau n'est pas bien concervée dans les ddfts, on est
      ! obligé de bidouillé.
      ! 1) soit on considère Pqisup, Eqi, et Pqiinf_par=Pqisup-Eqi
      !    et on suppose que dans la réalité les compositions de
      !    Pqiinf sont les même que Pqiinf_par
      ! 2) soit on considère Pqisup, Eqi_par=Pqisup-Pqiinf, et Pqiinf,
      !    et on suppose que dans la réalité les compositions de
      !    Eqi_prime sont les même que Eqi_par
          do il=1,ncas  
           if (Pqisup(cas(il)).gt.0.0) then  

            if ((water(cas(il)).gt.ridicule/100).and. &
                  (Pqiinf_par(cas(il)).le.0.0)) then
             ! on ne peut pas utiliser la méthode 1, car KE prédit de l'eau
             ! alors que le bilan de masse n'enprédit pas.
             ! Peut-on utiliser la méthode 2?
             Pqiinf_stewart(il)=Pqiinf(cas(il)) &
                *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
             Eqi_stewart(il)=Eqi_par(cas(il)) &
                *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
           else !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then
             ! il n'y a pas d'obstacles à l'utilisation de 1)
             Pqiinf_stewart(il)=Pqiinf_par(cas(il)) &
                  *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
             if (iflag_con.eq.30) then
                Eqi_stewart(il)=Eqi_prime(cas(il)) &
                  *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
             else
                Eqi_stewart(il)=Eqi(cas(il)) &
                  *(Pxtisup(ieau,cas(il))/Pqisup(cas(il)))
             endif
           endif !if ((water(il,i).gt.ridicule/100).and.(Pqiinf_par.le.0.0)) then

          else !if (Pqisup(cas(il).gt.0.0) then  

#ifdef ISOVERIF
              CALL iso_verif_egalite((Pqiinf(cas(il))), &
                 0.0,'compress_stewart 1347a')
               CALL iso_verif_egalite(( &
                 Eqi_prime(cas(il))),0.0,'compress_stewart 1347b')
               CALL iso_verif_egalite(( &
                 Pqiinf_par(cas(il))),0.0,'compress_stewart 1347c')
               CALL iso_verif_egalite((Eqi_par(cas(il))), &
                 0.0,'compress_stewart 1347d')
#endif              
               Pqiinf_stewart(il)=0.0
               Eqi_stewart(il)=0.0              
          endif !if (Pqisup(cas(il).gt.0.0) then  
         enddo !do il=1,ncas_evap_glace

        ! petite vérif
#ifdef ISOVERIF        
        do il=1,ncas
          if ((abs(water_cas(il)).ge.ridicule/10.) &
                 .and.(Pqiinf_stewart(il).le.0.0)) then
              write(*,*) 'compress_stewart 498: evap glace:'
              write(*,*) 'water(il,i)=', water_cas(il)
              write(*,*) 'Pqiinf=',Pqiinf(cas(il))
              write(*,*) 'Pqiinf_par=',Pqiinf_par(cas(il))
              write(*,*) 'Pqiinf_stewart=',Pqiinf_stewart(il)
              stop                   
          endif
        enddo !do il=1,ncas_evap_glace
#endif             

         END SUBROUTINE  compress_evap_glace_zone

         SUBROUTINE uncompress_ilp_zone( &
             ncas,cas, &
             zxtrfln_cas,zxt_cas,zxtrfl,zxtrfln,zxt,klon, &
             izone,Eqi,Exi,fac_ftmr, &
             xtrevap_tag,liq,hdiag)

    USE isotopes_mod, ONLY: ridicule,iso_eau
    USE isotrac_mod, only: option_revap,evap_franche
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif

         implicit none

        ! inputs         
         integer ncas
         integer cas(ncas)
         integer klon
         integer izone                  
         real zxt_cas(niso,ncas),zxtrfln_cas(niso,ncas)
         real Exi(niso,ncas)
         real Eqi(ncas)
         real fac_ftmr(ncas)
         integer liq
         real hdiag(ncas)

         ! output
         real zxt(ntraciso,klon)
         real zxtrfl(ntraciso,klon),zxtrfln(ntraciso,klon)
         real xtrevap_tag(ntraciso,ncas)

         ! locals
         integer il,ixt,iiso

!         write(*,*) 'uncompress_stewart 1420 tmp: zxt=',
!     :            zxt(iso_eau:ntraciso:3,cas(1))
!         write(*,*) 'Exi,fac_ftmr=',
!     :            Exi(iso_eau,1),fac_ftmr(1)
         do il=1,ncas
          do iiso=1,niso
            ixt=index_trac(izone,iiso)
            zxtrfln(ixt,cas(il))=zxtrfln_cas(iiso,il)
            zxtrfl(ixt,cas(il))=zxtrfln_cas(iiso,il)            
          enddo !do iiso=1,niso
!#ifdef ISOVERIF
!          if (il.eq.9) then
!            write(*,*) 'uncompress 1521'
!            write(*,*) 'il,Eqi,fac_ftmr,evap_franche,Exi(2,il)=',
!     :         il,Eqi(il),fac_ftmr(il),evap_franche,Exi(2,il)
!          endif
!#endif

          if ((option_revap.eq.1).and. &
               (((liq.eq.1).and.(Eqi(il)*fac_ftmr(il).gt.evap_franche) &
               .and.(hdiag(il).lt.0.99)).or. &
               ((liq.eq.0).and. &
               (Eqi(il)*fac_ftmr(il).ge.0.0)))) then
          ! le flux d'évap va dans un tag particulier
          ! -> zxt est inchangé mais xtrevap_tag(ixt,il) est incrémenté
             do iiso=1,niso
              ixt=index_trac(izone,iiso)
              xtrevap_tag(ixt,il)=fac_ftmr(il)*Exi(iiso,il)
!              zxt(ixt,cas(il))=zxt_cas(iiso,il)
!     :                   -xtrevap_tag(ixt,il)
             enddo !do iiso=1,niso
                
          else !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
            ! reequilibration standard
            do iiso=1,niso
              ixt=index_trac(izone,iiso)
              zxt(ixt,cas(il))=zxt(ixt,cas(il)) &
                         +fac_ftmr(il)*Exi(iiso,il)
              zxt(ixt,cas(il))=max(0.0,zxt(ixt,cas(il)))
              xtrevap_tag(ixt,il)=0.0
#ifdef ISOVERIF
              CALL iso_verif_positif_choix(zxt(ixt,cas(il)), &
                    0.0,'compress 1508')
#endif              
            enddo ! do iiso=1,niso           
          endif !if ((Eqi(il)*fac_ftmr(il).gt.evap_franche).and.
         enddo !do il=1,ncas
!         write(*,*) 'compress_stewart 1453 tmp: zxt=',
!     :            zxt(iso_eau:ntraciso:3,cas(1))

         END SUBROUTINE  uncompress_ilp_zone

         SUBROUTINE compress_ilp_evap_liq_zone( &
             ncas,cas, &
             zxt_cas,zxt, &
             zxtrfl_cas,zxtrfl_ancien, &
             zrfln_cas,zrfln,   &
             zrfl_cas,zrfl_ancien,     &
             zqev_diag_cas,zqev_diag,  &
             klon,izone,ptrac)

    USE isotopes_mod, ONLY: ridicule,iso_eau
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
         implicit none

        ! inputs         
         integer ncas
         integer cas(ncas)
         integer klon
         real zxt(ntraciso,klon)
         real zxtrfl(ntraciso,klon)
         real zrfl_ancien(klon)
         real zqev_diag(klon)
         real zrfln(klon)
         real zxtrfl_ancien(ntraciso,klon)
         integer izone

         ! outputs
         real zxt_cas(niso,ncas)         
         real zxtrfl_cas(niso,ncas)         
         real zqev_diag_cas(ncas)
         real zrfln_cas(ncas)         
         real zrfl_cas(ncas)
         real ptrac(ncas)
         
         ! locals
         integer il,ixt,ieau,iiso

         ieau=index_trac(izone,iso_eau)
         do il=1,ncas
          do iiso=1,niso
            ixt=index_trac(izone,iiso)
            ! la compo de la vap à l'extérieure reste la vapeur totale
            zxt_cas(iiso,il)=zxt(iiso,cas(il))
            ! le flux de pluie est celui le flux de pluie lié à la zone
            zxtrfl_cas(iiso,il)=zxtrfl_ancien(ixt,cas(il))
          enddo
          zrfl_cas(il)=zxtrfl_ancien(ieau,cas(il))

          if (zrfl_ancien(cas(il)).gt.0.0) then
          ! proportion de izone dans l'évap = celle dans la goutte          
            ptrac(il)=zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il))
            zrfln_cas(il)=zrfln(cas(il))*ptrac(il)
            zqev_diag_cas(il)=zqev_diag(cas(il))*ptrac(il)
!#ifdef ISOVERIF              
!            if (il.eq.9) then
!              write(*,*) 'compress tmp: il, ptrac=',il,ptrac(il)
!              write(*,*) 'ieau,zxtrfl_ancien(ieau,cas(il))=',
!     :           ieau,zxtrfl_ancien(ieau,cas(il))
!              write(*,*) 'zrfl_ancien(cas(il))=',zrfl_ancien(cas(il))
!              write(*,*) 'zrfl_cas(il)=',zrfl_cas(il)
!            endif
!#endif

          else !if (zrfl_ancien(cas(il)).gt.0.0) then
#ifdef ISOVERIF  
!              write(*,*) 'il,cas(il),zrfln,zrfl_ancien,zqev_diag=',
!     :              il,cas(il),zrfln(cas(il)),zrfl_ancien(cas(il)),
!     :              zqev_diag(cas(il))   
              CALL iso_verif_egalite(zqev_diag(cas(il)),0.0, &
                 'compress_stewart 1591a')
              CALL iso_verif_egalite(zrfln(cas(il)),0.0, &
                 'compress_stewart 1591b')
#endif               
                zrfln_cas(il)=0.0
                zqev_diag_cas(il)=0.0
               
          endif !if (zrfl_ancien(cas(il)).gt.0.0) then
          ! les lignes suvantes ne sont pas à recalculer
!          zt_cas(il)=zt(cas(il))
!          delP(il)=paprs(cas(il),k)-paprs(cas(il),k+1)
         enddo !do il=1,ncas

         END SUBROUTINE  compress_ilp_evap_liq_zone


         SUBROUTINE compress_ilp_evap_glace_zone( &
             ncas,cas, &
             zxt_cas,zxt, &
             zxtrfl_cas,zxtrfl_ancien, &
             zrfln_cas,zrfln,   &
             zrfl_cas, zrfl_ancien, &
             zqev_diag_cas,zqev_diag,  &
             klon,izone)

    USE isotopes_mod, ONLY: ridicule,iso_eau
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif

         implicit none

        ! inputs         
         integer ncas
         integer cas(ncas)
         integer klon
         real zxt(ntraciso,klon)
         real zxtrfl_ancien(ntraciso,klon)
         real zqev_diag(klon)
         real zrfln(klon)
         integer izone
         real zrfl_ancien(klon)

         ! outputs
         real zxt_cas(niso,ncas)         
         real zxtrfl_cas(niso,ncas)         
         real zqev_diag_cas(ncas)
         real zrfln_cas(ncas)         
         real zrfl_cas(ncas)
         
         
         ! locals
         integer il,ixt,ieau,iiso     

         ieau=index_trac(izone,iso_eau)
         do il=1,ncas
          do iiso=1,niso
            ixt=index_trac(izone,iiso)
            zxt_cas(iiso,il)=zxt(iiso,cas(il))
            zxtrfl_cas(iiso,il)=zxtrfl_ancien(ixt,cas(il))
          enddo          
          zrfl_cas(il)=zxtrfl_ancien(ieau,cas(il))

          if (zrfl_ancien(cas(il)).gt.0.0) then
             zrfln_cas(il)=zrfln(cas(il)) &
                 *(zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il)))
            ! car la proportion de traceurs dans zqev_diag et la même
            ! que dans zrfl_ancien. Comme zrfln=zrfl-zqev_diag*fac_ftmr
            ! alors cette proportion de traceurs est aussi la même dans
            ! zrfln
             zqev_diag_cas(il)=zqev_diag(cas(il)) &
               *zxtrfl_ancien(ieau,cas(il))/zrfl_ancien(cas(il))
          else !if (zrfl_ancien(cas(il)).gt.0.0) then
#ifdef ISOVERIF  
              CALL iso_verif_egalite(zqev_diag(cas(il)),0.0, &
                 'compress_stewart 1791a')
              CALL iso_verif_egalite(zrfln(cas(il)),0.0, &
                 'compress_stewart 1791b')
#endif               
                zrfln_cas(il)=0.0
                zqev_diag_cas(il)=0.0
          endif !if (zrfl_ancien(cas(il)).gt.0.0) then
         enddo

         END SUBROUTINE  compress_ilp_evap_glace_zone


      SUBROUTINE ajoute_revap(ncas,cas, &
          klon,izone,zxt,xtrevap_tag)

#ifdef ISOVERIF
USE isotopes_verif_mod
#endif
USE isotrac_mod, only: izone_revap
      implicit none

      ! ajoute xtrevap_tag (evaps des différents traceurs d'isotopes)
      ! dans la vapeur qui est taggée par la révap des gouttes.

      ! input
      integer ncas
      integer cas(ncas)
      integer klon
      integer izone      
      real xtrevap_tag(ntraciso,ncas)

      ! inout
      real zxt(ntraciso,klon)

      ! local
      integer i,ixt,iiso,ixt_revap
!#ifdef ISOVERIF      
!      integer iso_verif_positif_choix_nostop
!#endif          

#ifdef ISOVERIF
      do i=1,ncas
        do ixt=1,ntraciso
            CALL iso_verif_positif_choix(zxt(ixt,cas(i)),0.0, &
                         'ajoute_revap 29')
        enddo
      enddo !do i=1,ncas
#endif    

      do i=1,ncas
         do iiso=1,niso
            ixt_revap=index_trac(izone_revap,iiso)
            do izone=1,ntraceurs_zone
               ixt=index_trac(izone,iiso)
               zxt(ixt_revap,cas(i))=zxt(ixt_revap,cas(i)) &
                   +xtrevap_tag(ixt,i)
#ifdef ISOVERIF
               if (iso_verif_positif_choix_nostop(zxt(ixt_revap,cas(i)), &
                         0.0,'ajoute_revap 46').eq.1) then
                  write(*,*) 'i,iiso,izone,ixt=',i,iiso,izone,ixt
                  write(*,*) 'xtrevap_tag(ixt,i)=',xtrevap_tag(ixt,i)
!                  stop
               endif
#endif               
               zxt(ixt_revap,cas(i))=max(0.0,zxt(ixt_revap,cas(i)))
            enddo !do izone=1,ntraceurs_zone
         enddo !do iiso=1,niso
      enddo !do i=1,ncas_evap_liq

#ifdef ISOVERIF
      do i=1,ncas
        do ixt=1,ntraciso
            CALL iso_verif_positif_choix(zxt(ixt,cas(i)),0.0, &
                         'ajoute_revap 40')
        enddo
      enddo !do i=1,ncas
#endif      

      return
      END SUBROUTINE  ajoute_revap


      function is_in_bassin(lat,lon,bassin)
      USE isotrac_mod, only: use_bassin_atlantic,use_bassin_medit, &
&       use_bassin_indian,use_bassin_austral,use_bassin_pacific, &
&       use_bassin_MerArabie,use_bassin_BengalGolf,use_bassin_SouthIndian, &
&       use_bassin_tropics,use_bassin_midlats,use_bassin_HighLats, &
&       bassin_atlantic,bassin_medit, &
&       bassin_indian,bassin_austral,bassin_pacific, &
&       bassin_MerArabie,bassin_BengalGolf,bassin_SouthIndian, &
&       bassin_tropics,bassin_midlats,bassin_HighLats
      implicit none
      ! répond true si lat,lon se trouve dans le bassin numéroté bassin

      ! input
      integer bassin
      real lat,lon

      ! outut
      logical is_in_bassin

      ! locals
      !logical is_in_rectangle
      !logical is_in_triangle
      
      is_in_bassin=.FALSE.
#ifdef ISOVERIF      
      write(*,*) 'is_in_basin 84: entree,bassin=',bassin
#endif
      if (use_bassin_atlantic .and. bassin==bassin_atlantic) then
#ifdef ISOVERIF           
          write(*,*) 'bassin Atlantique?'
#endif          
          if (is_in_rectangle(lon,lat,-67.0,28.0,20.0,-45.0)) then
            ! boite sud
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_rectangle(lon,lat,-100.0,40.0,-5.3,28.0)) then
            ! ouest gibraltar    
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_rectangle(lon,lat,-100.0,48.0,0.0,40.0)) then
            ! Ouest France
            is_in_bassin=.TRUE.
            return
          endif 
          if (is_in_rectangle(lon,lat,-90.0,80.0,10.0,46.0)) then
            ! Atlantic Nord
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_triangle(lon,lat, &
                 -62.0,0.0,-62.0,30.0,-112.0,30.0)) then
            ! golfe du Mexique
            is_in_bassin=.TRUE.
            return
          endif

      else if (use_bassin_medit .and. bassin==bassin_medit) then
#ifdef ISOVERIF           
          write(*,*) 'bassin Medit?'
#endif          
          if (is_in_rectangle(lon,lat,0.0,48.0,45.0,29.0)) then
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_rectangle(lon,lat,-5.3,42.0,45.0,29.0)) then
            is_in_bassin=.TRUE.
            return
          endif

      else if (use_bassin_indian .and. bassin==bassin_indian) then
#ifdef ISOVERIF           
          write(*,*) 'bassin indian?'
#endif           
          if (is_in_rectangle(lon,lat,20.0,30.0,110.0,-45.0)) then
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_triangle(lon,lat, &
                 90.0,30.0,90.0,-45.0,150.0,-45.0)) then
            ! Ouest Australie
            is_in_bassin=.TRUE.
            return
          endif   

      else if (use_bassin_SouthIndian .and. bassin==bassin_SouthIndian) then
#ifdef ISOVERIF           
          write(*,*) 'bassin indian hemisphere Sud?'
#endif           
          if (is_in_rectangle(lon,lat,20.0,0.0,120.0,-45.0)) then
            is_in_bassin=.TRUE.
            return
          endif
          
      else if (use_bassin_MerArabie .and. bassin==bassin_MerArabie) then
#ifdef ISOVERIF           
          write(*,*) 'bassin Mer d''Arabie?'
#endif           
          if (is_in_rectangle(lon,lat,20.0,30.0,76.0,0.0)) then
            is_in_bassin=.TRUE.
            return
          endif 

      else if (use_bassin_BengalGolf .and. bassin==bassin_BengalGolf) then
#ifdef ISOVERIF           
          write(*,*) 'bassin Golfe du Bengale?'
#endif           
          if (is_in_rectangle(lon,lat,76.0,30.0,110.0,0.0)) then
            is_in_bassin=.TRUE.
            return
          endif         

      else if (use_bassin_pacific .and. bassin==bassin_pacific) then
#ifdef ISOVERIF           
          write(*,*) 'bassin Pacific?'
#endif          
         if (is_in_rectangle(lon,lat,-180.0,80.0,-100.0,-45.0)) then
             ! pacifique Est
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_rectangle(lon,lat,110.0,80.0,180.0,28.0)) then
            ! Pacifique Nord Ouest
            is_in_bassin=.TRUE.
            return
          endif 
          if (is_in_rectangle(lon,lat,120.0,80.0,180.0,-45.0)) then
            ! Pacifique central Sud
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_triangle(lon,lat, &
                  90.0,28.0,150.0,-45.0,150.0,28.0)) then
            ! Pacifque Sud Ouest
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_triangle(lon,lat, &
                  -62.0,0.0,-112.0,30.0,-112.0,0.0)) then
            ! Ouest Amérique centrale
            is_in_bassin=.TRUE.
            return
          endif
          if (is_in_rectangle(lon,lat,-180.0,0.0,-67.0,-45.0)) then
            ! Ouest Chili
            is_in_bassin=.TRUE.
            return
          endif

      else if (use_bassin_austral .and. bassin==bassin_austral) then  
#ifdef ISOVERIF           
          write(*,*) 'bassin austral?'
#endif           
          if (lat.lt.-45.0+0.2) then
                is_in_bassin=.TRUE.
                return    
          endif  

      else if (use_bassin_HighLats .and. bassin==bassin_HighLats) then  
#ifdef ISOVERIF           
          write(*,*) 'bassin hautes lats?'
#endif           
          if (abs(lat).gt.35.0) then
                is_in_bassin=.TRUE.
                return    
          endif

      else if (use_bassin_tropics .and. bassin==bassin_tropics) then  
#ifdef ISOVERIF           
          write(*,*) 'bassin tropics?'
#endif           
          if (abs(lat).lt.15.0) then
                is_in_bassin=.TRUE.
                return    
          endif

       else if (use_bassin_midlats .and. bassin==bassin_midlats) then  
#ifdef ISOVERIF           
          write(*,*) 'bassin mid lats?'
#endif           
          if ((abs(lat).ge.15.0).and. &
                 (abs(lat).le.35.0)) then
                is_in_bassin=.TRUE.
                return    
          endif    

      else
          write(*,*) 'iso_traceurs_routines 59: bassin inconnu'
          write(*,*) 'bassin_atlantic=' ,bassin_atlantic  
          write(*,*) 'bassin_medit=' ,bassin_medit
          write(*,*) 'bassin_indian=' ,bassin_indian
          write(*,*) 'bassin_austral=' ,bassin_austral
          write(*,*) 'bassin_MerArabie=' ,bassin_MerArabie
          write(*,*) 'bassin_BengalGolf=' ,bassin_BengalGolf
          write(*,*) 'bassin_SouthIndian=' ,bassin_SouthIndian
          write(*,*) 'use_bassin_atlantic=' ,use_bassin_atlantic  
          write(*,*) 'use_bassin_medit=' ,use_bassin_medit
          write(*,*) 'use_bassin_indian=' ,use_bassin_indian
          write(*,*) 'use_bassin_austral=' ,use_bassin_austral
          write(*,*) 'use_bassin_MerArabie=' ,use_bassin_MerArabie
          write(*,*) 'use_bassin_BengalGolf=' ,use_bassin_BengalGolf
          write(*,*) 'use_bassin_SouthIndian=' ,use_bassin_SouthIndian
          stop
      endif
      
      return
      end function is_in_bassin

      SUBROUTINE find_bassin(lat,lon,bassin)
      use isotrac_mod, only: izone_poubelle,ntraceurs_zone=>ntiso,option_traceurs, &
&        bassin_map
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif
      implicit none

      ! inputs
      real lat,lon
      ! output
      integer bassin
      !locals
      logical continu
      !logical is_in_bassin

      continu=.TRUE.
      bassin=1
     
#ifdef ISOVERIF       
      write(*,*) ''
      write(*,*) 'find bassin lat,lon=',lat,lon
#endif       
      do while (continu)
!#ifdef ISOVERIF       
!!      write(*,*) 'find_bassin 169: lat,lon,bassin=',lat,lon,bassin
!#endif       
         if (is_in_bassin(lat,lon,bassin)) then
                continu=.FALSE.
#ifdef ISOVERIF                 
                write(*,*) 'find_bassin 173: trouve: bassin=',bassin
#endif                 
         else
!#ifdef ISOVERIF              
!             write(*,*) 'find_bassin 175: pas trouve: bassin=',bassin
!#endif              
             bassin=bassin+1
         endif
         if (bassin.eq.izone_poubelle) then
                continu=.FALSE.
                bassin=izone_poubelle
!#ifdef ISOVERIF                 
!                write(*,*) 'find_bassin 179: poubelle: bassin=',bassin
!#endif  
         endif
       enddo

       ! normalement, le bassin est soit un bassin oce, soit un résidu
       ! donc bassin<=ntraceurs_zone-1
#ifdef ISOVERIF
       CALL iso_verif_positif(float(ntraceurs_zone-1-bassin), &
                 'find_bassin 195')
#endif

        return
        END SUBROUTINE  find_bassin

        SUBROUTINE initialise_bassins_boites(presnivs)
        USE dimphy, only: klev
        USE geometry_mod, ONLY: longitude_deg, latitude_deg
        use isotrac_mod, only: bassin_map,option_traceurs,boite_map
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif
        implicit none
        real presnivs(klev)

       if (option_traceurs.eq.3) then
           ! initialisation de bassin_map
           CALL bassin_map_init(latitude_deg,longitude_deg,bassin_map)
       else if (option_traceurs.eq.20) then
           ! initialisation de bassin_map selon < ou > 35° lat
           write(*,*) 'physiq 1681: init de la map pour tag 20'
           CALL bassin_map_init_opt20(latitude_deg,bassin_map)
       else if (option_traceurs.eq.5) then
           CALL boite_AMMA_init(latitude_deg,longitude_deg,presnivs,boite_map)
       else if (option_traceurs.eq.21) then
           CALL boite_UT_extra_init(latitude_deg,longitude_deg,presnivs,boite_map)
       endif

        return
        END SUBROUTINE  initialise_bassins_boites

        SUBROUTINE bassin_map_init(lat,lon,bassin_map)
        USE dimphy, only: klon
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif

        implicit none

        ! inputs
        real lat(klon),lon(klon)
        ! output
        integer bassin_map(klon)
        ! locals
        integer i

        do i=1,klon
             CALL find_bassin(lat(i),lon(i),bassin_map(i))
#ifdef ISOVERIF              
             write(*,*) 'init 233: i,lat,lon,bassin=',i,lat(i),lon(i), &
                  bassin_map(i)
#endif             
        enddo

        return
        END SUBROUTINE  bassin_map_init

        function is_in_rectangle(x,y,x1,y1,x2,y2)
        implicit none
        ! inputs 
        real x,y
        ! point en haut à gauche         
        real x1,y1
        ! point en bas à droite
        real x2,y2
        ! output
        logical is_in_rectangle

!#ifdef ISOVERIF        
!        write(*,*) 'is_in_rectange 237: x,y=',x,y
!        write(*,*) 'x1,y1,x2,y2=',x1,y1,x2,y2
!#endif         
        if ((x-x2.lt.0.1).and.(x-x1.gt.-0.1).and. &
                  (y-y1.lt.0.1).and.(y-y2.gt.-0.1)) then
          is_in_rectangle=.TRUE.
        else
          is_in_rectangle=.FALSE.
        endif
!#ifdef ISOVERIF        
!        write(*,*) 'is_in_rectangle=',is_in_rectangle
!#endif        
        return
        end function is_in_rectangle

        function is_in_triangle(x,y,x1,y1,x2,y2,x3,y3)
        implicit none
        ! inputs
        real x,y
        ! points dans le sens trigo
        ! à gauche
        real x1,y1
        ! en bas
        real x2,y2
        ! à droite
        real x3,y3
        ! output
        logical is_in_triangle
        ! locals
        real det1
        real det2
        real det3
!#ifdef ISOVERIF
!        write(*,*) 'is_in_triange 271: x,y=',x,y
!        write(*,*) 'x1,y1,x2,y2,x3,y3=',x1,y1,x2,y2,x3,y3
!#endif        
        det1=(x1-x)*(y2-y)-(y1-y)*(x2-x)
        det2=(x2-x)*(y3-y)-(y2-y)*(x3-x)
        det3=(x3-x)*(y1-y)-(y3-y)*(x1-x)
        is_in_triangle=.FALSE.
        if ((det1*det2.gt.0.0).and.(det2*det3.gt.0.0)) then
          is_in_triangle=.TRUE.
        else
          is_in_triangle=.FALSE.
        endif
!#ifdef ISOVERIF        
!        write(*,*) 'det1,det2,det3,is_in_triangle',
!     :         det1,det2,det3,is_in_triangle
!#endif
        return
        end function is_in_triangle


        SUBROUTINE isotrac_recolorise_tmin(xt,t)
        USE dimphy, only: klon, klev
        USE isotrac_mod, only: zone_temp,nzone_temp
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        real t(klon,klev)
        ! locals
        integer izone_temp 
        integer ixt,ixt_recoit
        integer k,i,izone,iiso
        !integer find_index
      

        do k=1,klev
          do i=1,klon
!#ifdef ISOVERIF
!            if (i.eq.1) then
!                write(*,*) 'recolorise 396: i,k,t=',i,k,t(i,k)
!                write(*,*) 'zone_temp=',zone_temp     
!            endif
!#endif           
            ! trouver la zone de cette température
            izone_temp=find_index(t(i,k),nzone_temp,zone_temp)

!#ifdef ISOVERIF
!            if (i.eq.1) then
!                write(*,*) 'recolorise 414: izone_temp=',izone_temp    
!            endif
!#endif            
            do izone=1,nzone_temp-1
               ! tous les tags de zone < nzone_temp se trouvant à des
               ! températures plus basses sont convertis
!#ifdef ISOVERIF
!            if (i.eq.1) then
!                write(*,*) 'recolorise 405: izone,xt_eau=',
!     :               izone,xt(index_trac(izone,iso_eau),i,k)  
!            endif
!#endif                 
               if (izone.lt.izone_temp) then                   
                  do iiso=1,niso
                   ixt=index_trac(izone,iiso) ! emmetteur
                   ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
                   xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
                   xt(ixt,i,k)=0.0 
                  enddo !do iiso=1,niso
               endif  !if (izone.lt.izone_temp) then
!#ifdef ISOVERIF
!            if (i.eq.1) then
!                write(*,*) 'recolorise 419: xt_eau,xt_recoit=',
!     :               xt(index_trac(izone,iso_eau),i,k),
!     :               xt(index_trac(izone_temp,iso_eau),i,k)
!            endif
!#endif                
            enddo !do izone=1,zone_pot(k)-1

            ! conversion de l'évap en surf et de la revap des gouttes
            do izone=nzone_temp+1,ntraceurs_zone 
              do iiso=1,niso 
                   ixt=index_trac(izone,iiso) ! emmetteur
                   ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
                   xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
                   xt(ixt,i,k)=0.0           
              enddo !do iiso=1,niso 
            enddo !do izone=nzone_temp+1,ntraceurs_zone 
          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 403')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_tmin

        SUBROUTINE isotrac_recolorise_tmin_sfrev(xt,t)
        USE dimphy, only: klon,klev
        USE isotrac_mod, only: nzone_temp,zone_temp
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! recolorise selon la température minimum, mais les tags de
        ! revap sont laissés en revap

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        real t(klon,klev)
        ! locals
        integer izone_temp 
        integer ixt,ixt_recoit
        integer k,i,izone,iiso
        !integer find_index

        do k=1,klev
          do i=1,klon
            ! trouver la zone de cette température
            izone_temp=find_index(t(i,k),nzone_temp,zone_temp)

            do izone=1,nzone_temp-1
               ! tous les tags de zone < nzone_temp se trouvant à des
               ! températures plus basses sont convertis
               ! sauf la revap
               ! le tag de la revap est nzone_temp+1=ntraceurs_zone
               if (izone.lt.izone_temp) then                   
                  do iiso=1,niso
                   ixt=index_trac(izone,iiso) ! emmetteur
                   ixt_recoit=index_trac(izone_temp,iiso) ! recepteur
                   xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
                   xt(ixt,i,k)=0.0 
                  enddo !do iiso=1,niso
               endif  !if (izone.lt.izone_temp) then
            enddo !do izone=1,zone_pot(k)-1

          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 594')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_tmin_sfrev


        SUBROUTINE isotrac_recolorise_saturation(xt,rh,lat,pres)
        USE dimphy, only: klon,klev
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! recolorise selon la température minimum, mais les tags de
        ! revap sont laissés en revap

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        real rh(klon,klev)
        real lat(klon)
        real pres(klev)
        ! locals
        integer izone_recoit
        integer ixt,ixt_recoit
        integer k,i,izone,iiso
        logical continu
        real rh_seuil
        parameter (rh_seuil=0.90)
        !integer index_zone_latpres

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 612')
          enddo !do i=1,klon
        enddo !do k=1,klev        
#endif

        ! on ne sature pas les 2 premières couches: on les laisse se
        ! recharger en evap de surface
        do k=3,klev
          do i=1,klon
            if (rh(i,k).gt.rh_seuil) then
                izone_recoit=index_zone_latpres(lat(i),pres(k))
                do izone=1,ntraceurs_zone
                 if (izone.ne.izone_recoit) then
                  do iiso=1,niso
                   ixt=index_trac(izone,iiso) ! emmetteur
                   ixt_recoit=index_trac(izone_recoit,iiso) ! recepteur
                   xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k)+xt(ixt,i,k)
                   xt(ixt,i,k)=0.0 
                  enddo !do iiso=1,niso
                 endif
                enddo !do izone=1,ntraceurs_zone
            endif !if (rh(i,k).gt.rh_seuil) then
          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 637')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_saturation

        SUBROUTINE isotrac_recolorise_boite(xt,boite_map)
        USE dimphy, only: klon,klev
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! SUBROUTINE écrite à la base pour tagguer 3 boites AMMA.
        ! Mais ça peut être générique, selon comment est initialisée boite_map

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        integer boite_map(klon,klev)
        ! locals
        integer i,k
        integer izone_recoit,izone,iiso
        integer ixt,ixt_recoit
        
        do k=1,klev
          do i=1,klon
            izone_recoit=boite_map(i,k)
            if (izone_recoit.gt.0) then
                ! on est dans une boite connue
                ! toutes les molécules sont converties en cete couleur
               do izone=1,ntraceurs_zone
                  if (izone.ne.izone_recoit) then
                      ! on met les traceurs izone dans izone_recoit
                      do iiso=1,niso
                        ixt=index_trac(izone,iiso)
                        ixt_recoit=index_trac(izone_recoit,iiso)
                        xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
                               +xt(ixt,i,k)
                        xt(ixt,i,k)=0.0
                      enddo
                  endif !if (izone.ne.izone_recoit) then
               enddo !do izone=2,ntraceurs_zone
            endif !if (izone_recoit.gt.0) then
          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 514')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_boite

        SUBROUTINE isotrac_recolorise_extra(xt,rlat)
        USE dimphy, only: klon,klev
        usE isotrac_mod, only: lim_tag20,izone_trop,izone_extra
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! SUBROUTINE écrite pour l'option de taggage 20
        ! permet de retagguer la vapeur tropicale en vapeur
        ! extratropicale dès qu'elle atteint 35° de latitude

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        real rlat(klon)
        ! locals
        integer i,k
        integer iiso,ixt,ixt_recoit
       
!        write(*,*) 'iso_traceurs_routines 723: lim_tag20=',lim_tag20
        do k=1,klev
          do i=1,klon
            if (abs(rlat(i)).gt.lim_tag20) then
                ! on met les traceurs izone_trop dans izone_extra
                do iiso=1,niso
                    ixt=index_trac(izone_trop,iiso)
                    ixt_recoit=index_trac(izone_extra,iiso)
                        xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
                               +xt(ixt,i,k)
                        xt(ixt,i,k)=0.0
                 enddo                  
            endif ! if (abs(rlat(i)).lt.35.0) then
          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 741')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_extra

        SUBROUTINE isotrac_recolorise_conv(xt,rlat,presnivs,rain_con)
        USE dimphy, only: klon,klev
        use isotrac_mod, only: lim_precip_tag22, &
&       izone_conv_BT,izone_conv_UT
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! SUBROUTINE écrite pour l'option de taggage 20
        ! permet de retagguer la vapeur tropicale en vapeur
        ! extratropicale dès qu'elle atteint 35° de latitude

        ! inout
        real xt(ntraciso,klon,klev)
        ! input
        real rlat(klon)
        real presnivs(klev)
        real rain_con(klon)
        ! locals
        integer i,k
        integer iiso,ixt,ixt_recoit,izone
       
!        write(*,*) 'iso_traceurs_routines 723: lim_tag20=',lim_tag20
!        write(*,*) 'presnivs=',presnivs
!        stop
        do k=1,klev
          do i=1,klon
#ifdef ISOVERIF
          if ((abs(rlat(i)).lt.30.0).and.(k.eq.1)) then
          endif
#endif          
            if ((abs(rlat(i)).lt.30.0).and. &
                  (rain_con(i)*86400.gt.lim_precip_tag22)) then
            
                ! on met les traceurs izone_trop dans izone_conv fn z                
                do iiso=1,niso
                  if (presnivs(k).gt.650.0*100.0) then
                    ixt_recoit=index_trac(izone_conv_BT,iiso) 
                  else
                    ixt_recoit=index_trac(izone_conv_UT,iiso) 
                  endif
                  do izone=1,ntraceurs_zone
                    ixt=index_trac(izone,iiso)
                    if (ixt.ne.ixt_recoit) then                    
                        xt(ixt_recoit,i,k)=xt(ixt_recoit,i,k) &
                               +xt(ixt,i,k)
                        xt(ixt,i,k)=0.0
                    endif !if (ixt.ne.ixt_recoit) then
                  enddo !do izone=1,ntraceurs_zone
                enddo !do iiso=1,niso
!                write(*,*) 'k,presnivs,ixt,ixt_recoit=',k,presnivs(k),
!     :                   ixt,ixt_recoit
!                write(*,*) 'xt(:,i,k)=',xt(:,i,k)
            endif ! if (abs(rlat(i)).lt.35.0) then
          enddo !do i=1,klon
        enddo !do k=1,klev

#ifdef ISOVERIF
        do k=1,klev
          do i=1,klon
            CALL iso_verif_traceur(xt(1,i,k),'recolorise 741')
          enddo !do i=1,klon
        enddo !do k=1,klev
#endif

        return
        END SUBROUTINE  isotrac_recolorise_conv


        SUBROUTINE boite_AMMA_init(lat,lon,presnivs,boite_map)
        USE dimphy, only: klon,klev
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        USE isotrac_mod, only: izone_aej,izone_mousson,izone_harmattan
        implicit none

        real lat(klon),lon(klon)
        real presnivs(klev)
        ! output
        integer boite_map(klon,klev)
        ! locals
        integer i,k

!        write(*,*) 'izone_aej,izone_mousson,izone_harmattan=',
!     :       izone_aej,izone_mousson,izone_harmattan    
        do k=1,klev
          do i=1,klon
                boite_map(i,k)=0.0
!                write(*,*) 'i,k,lat,lon,pres=',
!     :                   i,k,lat(i),lon(i),presnivs(k)
                if ((presnivs(k).le.700.0*100.0).and. &
                    (presnivs(k).gt.400.0*100.0).and. &
                    (lat(i).gt.8.0).and. &
                    (lat(i).lt.20.0).and. &
                    (lon(i).gt.10.0).and. &
                    (lon(i).lt.30.0)) then
                   boite_map(i,k)=izone_aej
!                   write(*,*) '   -> zone AEJ'
                else if ((presnivs(k).ge.850.0*100.0).and. &
                    (lat(i).gt.-5.0).and. &
                    (lat(i).le.8.0).and. &
                    (lon(i).gt.-40.0).and. &
                    (lon(i).lt.15.0)) then
                   boite_map(i,k)=izone_mousson
!                   write(*,*) '   -> zone flux de mousson'
                else if ((presnivs(k).gt.700.0*100.0).and. &
                    (lat(i).ge.20.0).and. &
                    (lat(i).lt.30.0).and. &
                    (lon(i).gt.-10.0).and. &
                    (lon(i).lt.40.0)) then
                   boite_map(i,k)=izone_harmattan
!                   write(*,*) '   -> zone Harmattan'
                endif
!                write(*,*) '   ** boite_map=',boite_map(i,k)
          enddo
        enddo

        return
        END SUBROUTINE  boite_AMMA_init

        
        SUBROUTINE boite_UT_extra_init(lat,lon,presnivs,boite_map)
        USE dimphy, only: klon,klev
        use isotrac_mod, only: izone_extra,izone_trop
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        real lat(klon),lon(klon)
        real presnivs(klev)
        ! output
        integer boite_map(klon,klev)
        ! locals
        integer i,k

!        write(*,*) 'izone_trop,izone_extra=',
!     :       izone_trop,izone_extra    
        do k=1,klev
          do i=1,klon
                boite_map(i,k)=0.0
!                write(*,*) 'i,k,lat,lon,pres=',
!     :                   i,k,lat(i),lon(i),presnivs(k)
                if ((presnivs(k).le.500.0*100.0) &
                    .and.(abs(lat(i)).lt.15.0)) then
                   boite_map(i,k)=izone_trop
!                   write(*,*) '   -> zone trop'
                else if (abs(lat(i)).gt.35.0) then
                   boite_map(i,k)=izone_extra
!                   write(*,*) '   -> zone extratropiques'
                endif
!                write(*,*) '   ** boite_map=',boite_map(i,k)
          enddo
        enddo

        return
        END SUBROUTINE  boite_UT_extra_init


        function index_zone_lat(lat)
        use isotrac_mod, only: lattag_min,dlattag,nzone_lat
        implicit none

        ! inputs
        real lat,pres
        ! output 
        integer index_zone_lat

        if (lat.lt.lattag_min) then
                index_zone_lat=1
        else
                index_zone_lat=int((lat-lattag_min)/dlattag)+2
                index_zone_lat=min(index_zone_lat,nzone_lat)
        endif

        return
        end function index_zone_lat


        function index_zone_pres(pres)
        use isotrac_mod, only: nzone_pres,zone_pres
        implicit none

        ! inputs
        real lat,pres
        ! output 
        integer index_zone_pres
        !integer find_index

        
        index_zone_pres=find_index(pres,nzone_pres,zone_pres) 
        write(*,*) 'iso_traceurs_routines 802: pres,index_zone_pres=', &
                 pres,index_zone_pres
        write(*,*) 'zone_pres=',zone_pres(1:nzone_pres-1)

        return
        end function index_zone_pres

        function find_index(pres,nzone_pres,zone_pres)
        implicit none

        ! inputs
        real pres
        integer nzone_pres
        real zone_pres(nzone_pres)
        ! output 
        integer find_index
        logical continu

        if (nzone_pres.gt.1) then
          if (pres.ge.zone_pres(1)) then
                find_index=1
          else if (pres.lt.zone_pres(nzone_pres-1)) then
                find_index=nzone_pres
          else !if (t(i,k).ge.zone_temp1) then
                continu=.TRUE.
                find_index=2
                do while (continu)
                  if (pres.ge.zone_pres(find_index)) then
                     continu=.FALSE.
                     ! c'est izone_temp, zone trouvée
                  else   
                     find_index=find_index+1
                  endif     
                enddo !do while (continu)
           endif !if (t(i,k).ge.zone_temp1) then

        else !if (nzone_pres.gt.1) then
            find_index=1
        endif !if (nzone_pres.gt.1) then

        return
        end function find_index

        function index_zone_latpres(lat,pres)
        use isotrac_mod, only: nzone_lat
        implicit none

        ! inputs
        real lat,pres
        ! output 
        integer index_zone_latpres
        ! locals
        integer index_lat
        integer index_pres
        !integer index_zone_lat
        !integer index_zone_pres

        index_lat=index_zone_lat(lat)
        index_pres=index_zone_pres(pres)
        index_zone_latpres=index_lat+(index_pres-1)*nzone_lat

        return
        end function index_zone_latpres

        SUBROUTINE iso_recolorise_condensation(qt,cond, &
                 xt,zxtcond,tcond,ep,xtres, &
                 seuil_in)
        USE dimphy, only: klon,klev
        USE isotopes_mod, only: bidouille_anti_divergence,iso_eau
        use isotrac_mod, only: option_seuil_tag_tmin,izone_cond, &
&       nzone_temp,zone_temp
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! on recolorise la vapeur résiduelle selon la température de condensation
        ! on supose qu'une vapeur xt,q condense en cond,zxtcond, à une
        ! température tcond. A ce stade, la vapeur initiale n'est pas
        ! retranchée de son condensat. On calcule les tags dans la vepru
        ! résiduelle xtres qu'on aurait si on retranchait un fraction ep du
        ! condensat

        ! inputs
        real qt
        real cond
        real tcond
        real ep
        real xt(ntraciso)
        real zxtcond(ntraciso)
        real seuil_in
        ! outputs
        real xtres(ntraciso)
        ! locals
        integer izone_temp,izone 
        integer ixt,ixt_recoit
        integer iiso
        !integer find_index
        real fcond, qmicro
!        real f

        if ((cond.gt.0.0).and.(qt.gt.0.0)) then        

            izone_temp=find_index(tcond,nzone_temp,zone_temp)
!            WRITE(*,*) 'pgm 901 tmp: izone_temp=',izone_temp    
#ifdef ISOVERIF
           do ixt=1,ntraciso
              CALL iso_verif_positif(xt(ixt)-zxtcond(ixt), &
                 'iso_trac 898')
           enddo  !do ixt=1,ntraciso  
           CALL iso_verif_traceur_justmass(xt, &
                'iso_trac_routines 906')
           CALL iso_verif_traceur_justmass(zxtcond, &
                'iso_trac_routines 908')
#endif           
          ! bidouille
          if (bidouille_anti_divergence) then
              CALL iso_verif_traceur_jbidouille(xt)
              CALL iso_verif_traceur_jbidouille(zxtcond)
          endif 

          do ixt=1,niso
            xtres(ixt)=xt(ixt)-ep*zxtcond(ixt)
          enddo
          do ixt=1+niso,ntraciso
            xtres(ixt)=0.0
          enddo
!          write(*,*) 'iso_trac_routines tmp 916: xtres=',xtres

#ifdef ISOVERIF 
        do ixt=1,ntraciso
          CALL iso_verif_positif(xtres(ixt), &
                    'iso_trac_routines 921')
        enddo
#endif          

             ! cas de izone sfc et izone precip et izone cond et izone< izone_temp

!             write(*,*) 'iso_trac 940: cond/qt,seuil_in,izone_temp=',
!     :                   cond/qt,seuil_in,izone_temp

             if (option_seuil_tag_tmin.eq.2) then
                 qmicro=0.0
                 do izone=nzone_temp+1,ntraceurs_zone 
                   ixt= index_trac(izone,iso_eau)  
                   qmicro=qmicro+xt(ixt)
                 enddo !do izone=nzone_temp+1,ntraceurs_zone
                 if (qt-qmicro.gt.0.0) then 
                     fcond=(cond-qmicro)/(qt-qmicro)
                 else
                     fcond=0.0   
                 endif
             else
                 fcond=cond/qt
             endif
       
             if (fcond.gt.seuil_in) then

             ! on les transfert à izone_temp
             do izone=1,ntraceurs_zone
               if ((izone.gt.nzone_temp).or.(izone.lt.izone_temp)) then
!                 ieau=index_trac(izone,iso_eau)
                 do iiso=1,niso  
                   ixt= index_trac(izone,iiso)    
                   ixt_recoit=index_trac(izone_temp,iiso) ! recepteur                 
                   xtres(ixt_recoit)=xtres(ixt_recoit) &
                         +(xt(ixt)-zxtcond(ixt))
                   xtres(ixt)=0.0   
!                   write(*,*) 'iso_trac 920: izone,ixt,',
!     :                 'ixt_recoit=',
!     :                 izone,ixt,ixt_recoit     
!                   write(*,*) 'isotrac 924: xt=',xt
!                   write(*,*) 'isotrac 925: zxtcond=',zxtcond
                 enddo !do iiso=1,niso
!                 write(*,*) 'iso_trac tmp 944: izone,xtres=',izone,xtres
                endif !if (izone.ne.izone_cond) then  
              enddo !do izone=nzones_temp+1,ntraceurs_zone

             else !if (cond/qt.gt.seuil_in) then
                
                ! on les laisse sur place
               do izone=1,ntraceurs_zone
                if ((izone.gt.nzone_temp).or.(izone.lt.izone_temp)) then
                 do iiso=1,niso  
                   ixt= index_trac(izone,iiso)    
                   xtres(ixt)=(xt(ixt)-zxtcond(ixt))
                 enddo !do iiso=1,niso
                endif !if (izone.ne.izone_cond) then  
               enddo !do izone=nzones_temp+1,ntraceurs_zone

             endif !if (cond/qt.gt.seuil_in) then

              ! izone_temp est conservé, on lui enlève juste son
              ! condesat
              do iiso=1,niso  
                   ixt_recoit=index_trac(izone_temp,iiso) ! recepteur                 
                   xtres(ixt_recoit)=xtres(ixt_recoit) &
                         +(xt(ixt_recoit)-zxtcond(ixt_recoit))
              enddo !do iiso=1,niso  

#ifdef ISOVERIF 
        do ixt=1,ntraciso
          CALL iso_verif_positif(xtres(ixt), &
                    'iso_trac_routines 940')
        enddo
#endif

             ! cas des zones > izone temp
             ! on conserve le condensat résiduel
             do izone=izone_temp+1,nzone_temp   
               do iiso=1,niso
                   ixt= index_trac(izone,iiso)
                   xtres(ixt)=xt(ixt)-zxtcond(ixt)
!                   write(*,*) 'iso_trac 931: izone,ixt,ixt_recoit=',
!     :                 izone,ixt,ixt_recoit       
!                   write(*,*) 'isotrac 934: xt=',xt
!                   write(*,*) 'isotrac 935: zxtice=',zxtice
               enddo !do iiso=1,niso  
!               write(*,*) 'iso_trac tmp 965: izone,xtres=',izone,xtres
             enddo !do izone=izone_temp+1,nzones_temp

             ! on rajoute le condensat qui ne precipite pas
             if (ep.lt.1.0) then
               do iiso=1,niso    
                   ixt= index_trac(izone_cond,iiso) 
                   xtres(ixt)=xtres(ixt)+(1.0-ep)*zxtcond(iiso)
!                   write(*,*) 'iso_trac 940: izone,ixt,ixt_recoit=',
!     :                 izone,ixt,ixt_recoit     
!                   write(*,*) 'isotrac 1014: xt=',xt
!                   write(*,*) 'isotrac 945: zxtice=',zxtice
                enddo  !do iiso=1,niso      
            endif !if (ep.lt.0.0) then

        else
            ! si cond=0 ou qt=0, tot reste pareil
           do ixt=1,ntraciso
                xtres(ixt)=xt(ixt)
           enddo  !do ixt=1,ntraciso     

        endif ! if (qt.gt.0.0) then    

#ifdef ISOVERIF
        if (iso_verif_traceur_jm_nostop(xtres, &
                'iso_trac_routines 166').eq.1) then
          write(*,*) 'isotrac 1024: xt=',xt
          write(*,*) 'zxtcond=',zxtcond
          write(*,*) 'xtres=',xtres
          write(*,*) 'ep=',ep
          stop
        endif 
#endif          
#ifdef ISOVERIF              
        do ixt=1,ntraciso
          CALL iso_verif_positif(xtres(ixt),'iso_trac_routines 953')
        enddo
        if (nzone_temp.ge.5) then 
        if (iso_verif_tag17_q_deltaD_chns(xtres, &
                 'iso_trac_routines 1025').eq.1) then
            write(*,*) 'xt=',xt
            write(*,*) 'zxtcond=',zxtcond
            write(*,*) 'xtres=',xtres
            write(*,*) 'ep=',ep
            write(*,*) 'tcond=',tcond
            write(*,*) 'izone_temp=',izone_temp
            stop
        endif
        endif
!       write(*,*) 'isotrac 1048: sortie de iso_recolorise_condensation'
#endif

            return
            END SUBROUTINE  iso_recolorise_condensation

        SUBROUTINE bassin_map_init_opt20(lat,bassin_map)
        USE dimphy, only: klon
        use isotrac_mod, only: izone_cont,izone_trop,lim_tag20
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none

        ! inputs
        real lat(klon)
        ! output
        integer bassin_map(klon)
        ! locals
        integer i

        write(*,*) 'iso_traceurs_routines 1142: lim_tag20=',lim_tag20
        do i=1,klon
         if (abs(lat(i)).gt.lim_tag20) then   
             bassin_map(i)=izone_cont
         else  
             bassin_map(i)=izone_trop
         endif
        enddo !do i=1,klon
        
        return
        END SUBROUTINE  bassin_map_init_opt20

        SUBROUTINE isotrac_recolorise_general(xt_seri,t_seri,zx_rh,presnivs)
        USE geometry_mod, ONLY: latitude_deg
        USE dimphy, only: klon,klev
        use isotrac_mod, only: option_traceurs,boite_map
        implicit none

        ! inputs
        real, dimension(ntraciso,klon,klev), intent(in) :: xt_seri
        real, dimension(klon,klev), intent(in) :: t_seri
        real, dimension(klon,klev), intent(in) :: zx_rh
        real, dimension(klev), intent(in) :: presnivs

              if (option_traceurs.eq.4) then
         CALL isotrac_recolorise_tmin(xt_seri,t_seri)
      elseif ((option_traceurs.eq.5).or. &
                  (option_traceurs.eq.21)) then
        CALL isotrac_recolorise_boite(xt_seri,boite_map)
      elseif (option_traceurs.eq.13) then
        CALL isotrac_recolorise_tmin_sfrev(xt_seri,t_seri)
      elseif (option_traceurs.eq.14) then
        CALL isotrac_recolorise_saturation(xt_seri,zx_rh,latitude_deg,presnivs)
      elseif (option_traceurs.eq.20) then      
        CALL isotrac_recolorise_extra(xt_seri,latitude_deg)
      endif !if (option_traceurs.eq.4) then

        return
        END SUBROUTINE  isotrac_recolorise_general


        

        SUBROUTINE iso_verif_traceur_jbid_vect(x,n,m)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
        !use isotrac_mod, only: ntraceurs_zone=>nzone
        USE infotrac_phy, ONLY: ntraceurs_zone=>nzone
        implicit none
        
        ! version vectrisée de iso_verif_traceur_jbidouille
            ! inputs
       integer n,m
       real x(ntraciso,n,m)

       ! locals
       integer iiso,izone,ixt,i,j
       real xtractot(n,m)
                
        if (bidouille_anti_divergence) then        
        do iiso=1,niso

          do j=1,m
           do i=1,n          
            xtractot(i,j)=0.0
           enddo !do j=1,m
          enddo !do j=1,m

          do izone=1,ntraceurs_zone  
            ixt=index_trac(izone,iiso) 
            do j=1,m
             do i=1,n
              xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)
             enddo !do j=1,m
            enddo !do j=1,m
           enddo !do izone=1,ntraceurs_zone
         
              ! on réajuste pour que les traceurs fasses bien la somme
              ! des traceurs
           do izone=1,ntraceurs_zone
            ixt=index_trac(izone,iiso)
            do j=1,m
             do i=1,n
!              if (abs(xtractot(i,j)).gt.ridicule*10) then 
               if (abs(xtractot(i,j)).gt.ridicule) then 
                   ! modif le 19 fev 2011
                  x(ixt,i,j)=x(ixt,i,j)/xtractot(i,j)*x(iiso,i,j) 
              endif !if (abs(xtractot(i,j)).gt.ridicule*10) then
             enddo !do i=1,n
            enddo !do j=1,m   
           enddo !do izone=1,ntraceurs_zone 
 
!           ! ajout le 19 fev 2011
!           ! on rend plutot les vérifs plus strictes
!           ixt=index_trac(izone_poubelle,iiso)
!           do j=1,m
!             do i=1,n
!              if ((abs(xtractot(i,j)).lt.1e-18).and.
!     :           (x(iiso,i,j).gt.ridicule)) then
!                  x(ixt,i,j)=x(iiso,i,j)
!              endif !if (abs(xtractot(i,j)).gt.ridicule*10) then
!             enddo ! do i=1,n
!           enddo !do j=1,m

        enddo !do iiso=1,ntraceurs_iso  
        endif !if (bidouille_anti_divergence) then    

        END SUBROUTINE  iso_verif_traceur_jbid_vect

        SUBROUTINE iso_verif_traceur_jbidouille(x)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
        implicit none
        
        ! on réajuste aussi les valeurs des traceurs pour la
        ! conservation de la masse, dans le cas bidouille
       
       ! inputs
       real x(ntraciso)

       ! locals
       integer iiso,izone,ixt
       real xtractot
        
        if (bidouille_anti_divergence) then

        do iiso=1,niso

          xtractot=0.0
          do izone=1,ntraceurs_zone  
            ixt=index_trac(izone,iiso) 
            xtractot=xtractot+x(ixt)
          enddo !do izone=1,ntraceurs_zone
         
              ! on réajuste pour que les traceurs fasses bien la somme
              ! des traceurs
              if (abs(xtractot).gt.ridicule*10) then
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt)=x(ixt)/xtractot*x(iiso)
                enddo !do izone=1,ntraceurs_zone                
              endif         

        enddo !do iiso=1,ntraceurs_iso   
        endif !if (bidouille_anti_divergence) then     

        END SUBROUTINE  iso_verif_traceur_jbidouille


        SUBROUTINE iso_verif_traceur_jbid_pos(x)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
!#ifdef ISOVERIF
!        use isotopes_verif_mod, only: iso_verif_traceur_pbidouille
!#endif
        implicit none
        
        ! on réajuste les valeurs des traceurs pour qu'il n'y ai pas de
        ! valeurs négatives. Si valeurs négatives -> on pompe les autres
        ! traceurs
        ! attention: fait la même chose pour tous les isos -> peut
        ! induire des fractionnements.
        ! Pour ne pas induire des fractionnements, prendre 
        ! iso_verif_traceur_jbid_pos2
        ! avantage de cette subroutine: conserve la masse en isotopes
        ! légers, ce qui nest pas le cas de pos2
          
       ! inputs
       real x(ntraciso)

       ! locals
       integer iiso,izone,ixt
       real xtractot,xtractotprec
        
        if (bidouille_anti_divergence) then

!            write(*,*) 'pgm 532 tmp: x=',x
        do iiso=1,niso

          xtractot=0.0
          xtractotprec=0.0
          do izone=1,ntraceurs_zone  
            ixt=index_trac(izone,iiso) 
            xtractotprec=xtractotprec+x(ixt)
            x(ixt)=max(x(ixt),0.0) 
            xtractot=xtractot+x(ixt)
          enddo !do izone=1,ntraceurs_zone
!          write(*,*) 'iiso,xtractotprec,xtractot=',
!     :          iiso,xtractotprec,xtractot
         
          if (xtractot.gt.xtractotprec) then
              ! on réajuste pour que les traceurs fasses bien la somme
              ! des traceurs
              if (abs(xtractot).gt.ridicule) then
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt)=x(ixt)*xtractotprec/xtractot
                enddo !do izone=1,ntraceurs_zone  
                ! on modifie aussi l'isotope de base si lui aussi était
                ! négatif   
!                x(iiso)=xtractot
              else  !if (abs(xtractot).gt.ridicule) then
                ! normallement, valeurs restantes très faibles
                ! on ne fait rien. 
                ! on met juste un max
                x(iiso)=max(x(iiso),0.0)
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt)=max(x(ixt),0.0)
                enddo !do izone=1,ntraceurs_zone
              endif !if (abs(xtractot).gt.ridicule) then
          endif !if (xtractot.gt.xtractotprec) then      

        enddo !do iiso=1,ntraceurs_iso
#ifdef ISOVERIF        
        CALL iso_verif_traceur_pbidouille(x,'iso_verif_trac 558')
#else
        CALL iso_verif_traceur_jbidouille(x)
#endif
        endif !if (bidouille_anti_divergence) then     

        END SUBROUTINE  iso_verif_traceur_jbid_pos

        SUBROUTINE iso_verif_traceur_jbid_pos_vect(n,m,x)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
        implicit none
       
       ! inputs
       integer n,m
       real x(ntraciso,n,m)

       ! locals
       integer iiso,izone,ixt
       real xtractot(n,m),xtractotprec(n,m)
       integer i,j
        
        if (bidouille_anti_divergence) then

!            write(*,*) 'pgm 532 tmp: x=',x
               
        do iiso=1,niso          
          do j=1,m 
          do i=1,n
          xtractot(i,j)=0.0
          xtractotprec(i,j)=0.0
          enddo !do j=1,m 
          enddo !do i=1,n
          do izone=1,ntraceurs_zone 
            ixt=index_trac(izone,iiso)  

            do j=1,m 
            do i=1,n
            xtractotprec(i,j)=xtractotprec(i,j)+x(ixt,i,j)
            x(ixt,i,j)=max(x(ixt,i,j),0.0) 
            xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)
            enddo !do i=1,n
            enddo !do j=1,m 
            
          enddo !do izone=1,ntraceurs_zone
!          write(*,*) 'iiso,xtractotprec,xtractot=',
!     :          iiso,xtractotprec,xtractot
                  
         do j=1,m 
         do i=1,n
          if (xtractot(i,j).gt.xtractotprec(i,j)) then
              ! on réajuste pour que les traceurs fasses bien la somme
              ! des traceurs
              if (abs(xtractot(i,j)).gt.ridicule) then
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt,i,j)=x(ixt,i,j)*xtractotprec(i,j)/xtractot(i,j)
                enddo !do izone=1,ntraceurs_zone  
                ! on modifie aussi l'isotope de base si lui aussi était
                ! négatif   
!                x(iiso)=xtractot
              else  !if (abs(xtractot).gt.ridicule) then
                ! normallement, valeurs restantes très faibles
                ! on ne fait rien. 
                ! on met juste un max
                x(iiso,i,j)=max(x(iiso,i,j),0.0)
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt,i,j)=max(x(ixt,i,j),0.0)
                enddo !do izone=1,ntraceurs_zone
              endif !if (abs(xtractot).gt.ridicule) then
          endif !if (xtractot.gt.xtractotprec) then
         enddo !do i=1,n      
         enddo !do j=1,m 
         

        enddo !do iiso=1,ntraceurs_iso
#ifdef ISOVERIF        
        CALL iso_verif_traceur_pbid_vect(x,n,m,'iso_verif_trac 558')
#else
        CALL iso_verif_traceur_jbid_vect(x,n,m)
#endif
        endif !if (bidouille_anti_divergence) then     

        END SUBROUTINE  iso_verif_traceur_jbid_pos_vect

        SUBROUTINE iso_verif_traceur_jbid_pos2(x,q)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
#ifdef ISOVERIF
        use isotopes_verif_mod
#endif        
        implicit none
        
        ! même but que iso_verif_traceur_jbid_pos, mais n'induit
        ! pas de fractionnement.
        ! on regarde si xteau est positif. S'il ne l'est pas, on pompe
        ! dans les autres tags pour le mettre à 0. On conserve la compo
        ! iso.
        ! Pb: ne conserve pas la masse d'isotopes légers.
    
       ! inputs
       real x(ntraciso),q

       ! locals
       integer iiso,izone,ixt,ieau
       real dqtmp,factmp
        
        if (bidouille_anti_divergence) then
!        write(*,*) 'iso_verif_trac 578 tmp: q,xt=',
!     :                q,x(1:ntraciso)
        if (q.gt.0.0) then
              dqtmp=0.0
              do izone=1,ntraceurs_zone
                ieau=index_trac(izone,iso_eau)
                if (x(ieau).lt.0.0) then
!                    write(*,*) 'local_x<0 pour izone=',izone
                    dqtmp=dqtmp-x(ieau)
                    do iiso=1,niso    
                      ixt=index_trac(izone,iiso)
                      x(ixt) =0.0
                    enddo !do iiso=1,niso  
                endif  !if (local_xt(ieau,i,k).lt.0.0) then                
              enddo !do izone=1,ntraceurs_zone
!              write(*,*) 'dqtmp=',dqtmp
              if (dqtmp.gt.0.0) then  
!                  write(*,*) 'iso_verif_trac 593 warning: q,dqtmp,xt=',
!     :                q,dqtmp,x(1:ntraciso)
                    ! on redistribue la négativité des traceurs dans les
                    ! traceurs positifs
!                    factmp=(1.0-dqtmp/(local_q(i,k)+dqtmp))
                    ! correction janv 2010
                    factmp=(q/(q+dqtmp))
!                    write(*,*) 'factmp=',factmp
                     do izone=1,ntraceurs_zone
                       ieau=index_trac(izone,iso_eau)
                       if (x(ieau).gt.0.0) then
                          do iiso=1,niso    
                             ixt=index_trac(izone,iiso)
                             x(ixt)=x(ixt)*factmp
                          enddo !do iiso=1,niso
                       endif !if (local_xt(ieau,i,k).gt.0.0) then
                     enddo ! do izone=1,ntraceurs_zone
!                     write(*,*) 'apres bidouille: xt=',x(1:ntraciso)
              endif !if (dqtmp.gt.0.0) then   
#ifdef ISOVERIF
              CALL iso_verif_traceur(x,'iso_verif_traceurs 612')
#endif
          endif !if (local_q(i,k).lt.0.0) then

#ifdef ISOVERIF
        CALL iso_verif_traceur_pbidouille(x,'iso_verif_trac 625')
#endif
        endif ! if (bidouille_anti_divergence) then  

        END SUBROUTINE  iso_verif_traceur_jbid_pos2

        SUBROUTINE iso_verif_traceur_jbid_vect1D(x,n)
        USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
        implicit none
        
        ! version vectrisée de iso_verif_traceur_jbidouille
           
       ! inputs
       integer n
       real x(ntraciso,n)

       ! locals
       integer iiso,izone,ixt,i
       real xtractot
        
        if (bidouille_anti_divergence) then

        do i=1,n
        do iiso=1,niso

          xtractot=0.0
          do izone=1,ntraceurs_zone  
            ixt=index_trac(izone,iiso) 
            xtractot=xtractot+x(ixt,i)
          enddo !do izone=1,ntraceurs_zone
         
              ! on réajuste pour que les traceurs fasses bien la somme
              ! des traceurs
              if (abs(xtractot).gt.ridicule*10) then
                do izone=1,ntraceurs_zone
                  ixt=index_trac(izone,iiso) 
                  x(ixt,i)=x(ixt,i)/xtractot*x(iiso,i)
                enddo !do izone=1,ntraceurs_zone                
              endif         

        enddo !do iiso=1,ntraceurs_iso  
        enddo  !do i=1,n
        endif !if (bidouille_anti_divergence) then             

        END SUBROUTINE  iso_verif_traceur_jbid_vect1D
      
! on met ces routines ici pour éviter dépendances circulaires   
#ifdef ISOVERIF

        SUBROUTINE iso_verif_traceur_pbidouille(x,err_msg)
        use isotopes_verif_mod
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant
        ! on réajuste aussi les valeurs des traceurs pour la
        ! conservation de la masse, dans le cas bidouille

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! local
       !integer iso_verif_traceur_pbid_ns

        if (iso_verif_traceur_pbid_ns(x,err_msg).eq.1) then
            stop
        endif

        END SUBROUTINE  iso_verif_traceur_pbidouille

        function iso_verif_traceur_pbid_ns(x,err_msg)
        use isotopes_mod, ONLY: iso_HDO,bidouille_anti_divergence
        use isotrac_mod, only: ridicule_trac
        use isotopes_verif_mod
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant
        ! on réajuste aussi les valeurs des traceurs pour la
        ! conservation de la masse, dans le cas bidouille

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! output
       integer iso_verif_traceur_pbid_ns

       ! locals
       !integer iso_verif_traceur_noNaN_nostop
       !integer iso_verif_tracm_choix_nostop  
       !integer iso_verif_tracdD_choix_nostop  
       integer iiso,izone,ixt
       real xtractot
        
        ! verif noNaN
        iso_verif_traceur_pbid_ns=0

        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
!             stop
            iso_verif_traceur_pbid_ns=1
        endif

        ! verif masse
        if (iso_verif_tracm_choix_nostop(x,err_msg, &
                 errmax*10,errmaxrel*50).eq.1) then
             ! on est plus laxiste car ça vient en général après une
             ! bidouille pour iso_eau normal
!             stop
             iso_verif_traceur_pbid_ns=1   
        endif  

        if (bidouille_anti_divergence) then
            ! on réajuste pour que les traceurs fasses bien la somme
            ! des traceurs
            CALL iso_verif_traceur_jbidouille(x)
        endif !if (bidouille_anti_divergence) then    

       ! verif deltaD
       if (iso_HDO.gt.0) then
        if (iso_verif_tracdD_choix_nostop(x,err_msg, &
                 ridicule_trac,deltalimtrac).eq.1) then
!             stop
              iso_verif_traceur_pbid_ns=1
        endif  
       endif !if (iso_HDO.gt.0) then

        end function iso_verif_traceur_pbid_ns

        SUBROUTINE iso_verif_traceur_pbid_vect(x,n,m,err_msg)
        use isotopes_mod, ONLY: iso_HDO,bidouille_anti_divergence
        use isotopes_verif_mod
        implicit none
       
       ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher       

       ! locals
       integer iiso,izone,ixt
       real xtractot
        
        ! verif noNaN
        CALL iso_verif_traceur_noNaN_vect(x,n,m,err_msg)

        ! verif masse
        CALL iso_verif_trac_masse_vect(x,n,m,err_msg,errmax*10, &
                 errmaxrel*50)

        if (bidouille_anti_divergence) then
            ! on réajuste pour que les traceurs fasses bien la somme
            ! des traceurs
            CALL iso_verif_traceur_jbid_vect(x,n,m)
        endif !if (bidouille_anti_divergence) then    

       ! verif deltaD
       if (iso_HDO.gt.0) then
        CALL iso_verif_tracdd_vect(x,n,m,err_msg)
       endif

        END SUBROUTINE  iso_verif_traceur_pbid_vect
#endif

END MODULE isotrac_routines_mod
#endif
#endif
