MODULE traject3d
   !!======================================================================
   !!                   ***  MODULE traject3d  ***
   !! Trajectories dynamics:  time evolution of a patch in 3D
   !&E
   !&E Description :
   !&E
   !&E       ! At the location of ssh, the middle of the cell, the position of the particule
   !&E       ! (xpos ypos variables) is equal to the indexes i j.
   !&E       ! location of ssh : xpos=REAL(i)
   !&E       ! in a cell, REAL(i)-0.5 <= xpos <= REAL(i) + 0.5 - epsilon
   !&E       ! in a cell, REAL(j)-0.5 <= ypos <= REAL(j) + 0.5 - epsilon
   !&E
   !&E ** History :
   !&E       !  2010-05 (M. Sourisseau) Introduction of derivate type
   !&E       !  2011-01 (M. Huret) Review and introduction of generalized sigma
   !&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E       !  2011-01 (M. Huret) call from step (if key_trajec3d) or from
   !&E       !           module ibm (if key_ibm in case of particle tracking with biology)
   !&E       !  2011-09 (T. Odaka, V. Garnier) MPI
   !&E       !  2011-10 (T. Odaka, M. Huret, V. Garnier) Introduction of MPI
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,  ONLY : imin,imax,jmin,jmax,kmax,ljmin,ljmax,ljmaxp1,  &
                           limaxp1,liminm1,limaxp2,ljminm1,ljmaxp2,limax, &
                           rsh,rlg,riosh,lchain,limin,liminm2,ljminm2,valmanq

   IMPLICIT NONE
   PRIVATE

#if defined key_ibm || defined key_trajec3d

   !! * Accessibility
   PUBLIC :: traj_3d    ! routine called by step3dxy.F90 step3dyx.F90, ibm_3d

   !! * Shared module variables

   !! * Private variables
   REAL(KIND=rsh), parameter :: r = 1.0_rsh/3.0_rsh

 CONTAINS

!!======================================================================

  SUBROUTINE traj_3d(xe,uz,vz)
 
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_3d  ***
   !&E
   !&E ** Purpose : calcul de l advection dispersion d un nuage de
   !&E              particules.
   !&E
   !&E ** Description : Methode de random walk version non naive
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference : Visser A. W., Marine ecology progress series,
   !&E                vol 158, pp275-281, annee 1997 : Using random
   !&E                walk models to simulate the vertical distribution of
   !&E                particles in a turbulent water column
   !&E
   !&E ** History :
   !&E       !          (P. Lazure, F. Dumas)  Original code
   !&E       !  2008-08 (F. Dumas ?)
   !&E       !  2007-12 (M. Huret) ajout modif du profil kz (smooth + spline)
   !&E       !  see (Ross and Sharples, 2004, Limnol.Oceanogr.)and 
   !&E       !  North et al. 2006. to garantee the Well Mixed Condition test
   !&E       !  where sharp kz gradients occur
   !&E       !  2010-01  (F Batifoulier, M Sourisseau) 
   !&E       !  2010-05 (M. Sourisseau)  Introduction of derivate type
   !&E       !  2011-01 (M. Huret) Move of biology and behaviour in module ibm
   !&E       !  2011-01 (M. Huret) Adaptation to siggen, RW in z
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : t,dt,h0_g
   USE comvars3d,    ONLY : sigw
   USE trajinitsave, ONLY : traj_save3d
   USE trajectools,  ONLY : h0int,xeint,ksupkinf,loc_h0,wint,dksdzint,splint,kzprofile
   USE comtraj,      ONLY : patches, type_patch, type_particle, ndtz
   USE_MPI toolmpi,  ONLY : ex_traj
   USE_MPI comtraj,  ONLY : down_give, up_give, right_give, left_give
#if defined key_siggen
   USE trajectools,  ONLY : siggentoz, ztosiggen, hc_sigint
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),      INTENT( in )    :: xe
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )    :: uz,vz

   !! * Local declarations
   INTEGER              :: n,m,nt,l,nb_patch
   REAL(KIND=rsh)       :: x_old,y_old,s_old,z_old,z_int,s_int,slag,     &
                           x_mid,y_mid,z_mid,d3_mid,xe_mid,h0_mid,       &
                           tir,dksdz,kzz,ds_adv,ds_dif,d3,d3_final,      &
                           xe_final,h0_final,                            &
                           lb,lt,lbs,lts,ds_bio,wbio,d3avt,xeavt,hcavt,h0avt
#if defined key_siggen
   REAL(KIND=rsh)       :: hc_sig_mid, hc_sig_final, zposf, zlag
   INTEGER              :: kw, kwm
#endif
   REAL(KIND=rlg)       :: dtm,dtz
   REAL(KIND=rsh), DIMENSION(0:kmax) :: dkz2,kzsmth
   REAL(KIND=rsh), DIMENSION(0:kmax) :: zpos  
   TYPE(type_patch),    POINTER :: patch    => NULL()
   TYPE(type_particle), POINTER :: particle => NULL()

!  parameters for interpolation at particle location
   REAL(KIND=rsh)       :: px,py
   INTEGER              :: igg,idd,jbb,jhh,hlb,hlt,hrb,hrt,ipos,jpos

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter TRAJ_3D routine'

   ! Initialization
   dtm = dt*0.5_rlg
   nb_patch = patches % nb

   IF ( nb_patch /= 0 ) dtz = dtm/real(ndtz)

  ! We modify spos from zpos because it was first calculated at simulation intial time step with another xe...
   ! to be removed when reading of positions at debtraj time step and not beginning of run
   patch => patches % first
   DO n = 1,nb_patch
   
      IF ( (t < patch % t_beg) .OR. (t > patch % t_end) ) THEN
         patch => patch % next
         CYCLE
      END IF

      IF ( (t-patch%t_beg) <= dtm ) THEN ! first path for this patch
         ! We modify spos from zpos because it was first calculated at simulation intial time step with another xe...
         DO m = 1, patch % nb_part_alloc
            particle => patch % particles(m)
            IF ( .NOT. particle % active ) CYCLE

            ! total depth at particle s location
            CALL loc_h0(particle % xpos, particle % ypos,  &
                        px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            particle % xe = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            particle % d3 = particle % h0 + particle % xe
            IF ( particle % d3 > particle % zpos ) THEN                
#ifdef key_siggen        
               zlag = -(particle % zpos) + particle % xe ! switch from immersion to real z
               CALL ztosiggen(zlag, slag, particle % xe, particle % h0, particle % hc)
#else
               slag = -(particle % zpos) / particle % d3
#endif		
               particle  % spos = slag
            ELSE
               particle  % zpos = particle  % h0
               particle  % spos = -1.0_rsh
            END IF            
         END DO
      END IF

      patch => patch % next
            
   END DO

!**********************************
! Ecriture des sorties (before any advection to get the exact release location in the output file
!**********************************
#ifdef key_agrifs
   IF(Agrif_Root()) CALL traj_save3d
#else
# ifdef key_trajec3d 
   CALL traj_save3d
# endif
#endif

#ifdef key_MPI_2D
   down_give  = 0
   up_give    = 0
   right_give = 0
   left_give  = 0
#endif

   !***************************
   ! Start loop on patches 
   !***************************
   patch => patches % first
   DO n = 1,nb_patch

      IF ( (t < patch % t_beg) .OR. (t > patch % t_end) ) THEN
         patch => patch % next
         CYCLE
      END IF
   
      DO m = 1, patch % nb_part_alloc

         particle => patch % particles(m)

         ! Skip if particle is inactive.
         IF ( .NOT. particle % active ) CYCLE
#ifdef key_species
         ! Skip if species stage not appropriate
         IF ( particle % stage >= 5 .OR. particle % super <= 0.0_rsh )  CYCLE
#endif
         ! Skip if flag is missing...
         IF ( particle % flag == -valmanq )  CYCLE

         ! depth of cell in which particle is located
         d3 = h0_g(NINT(particle%xpos),NINT(particle%ypos))+xe(NINT(particle%xpos),NINT(particle%ypos))

         IF ( d3 <= 0.0_rsh ) CYCLE

      !***************************
      ! horizontal advection 
      !***************************

         ! Save former position
         x_old = particle % xpos
         y_old = particle % ypos
         s_old = particle % spos 
         d3avt = particle % d3
         xeavt = particle % xe
         h0avt = particle % h0
         hcavt = particle % hc

         ! along-sigma advection
         CALL avance(uz,vz,xe,dtm,particle%xpos,particle%ypos,particle%spos,particle%flag) !FB flag

         CALL loc_h0(particle%xpos,particle%ypos, &
                      px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         xe_final = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         h0_final = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         d3_final = h0_final+xe_final
#ifdef key_siggen 
         hc_sig_final = hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
#endif
         x_mid = 0.5_rsh*(particle%xpos+x_old) ! Mid-position between old(s_old) and new (particle%xpos)
         y_mid = 0.5_rsh*(particle%ypos+y_old)
         CALL loc_h0(x_mid,y_mid,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         xe_mid = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         h0_mid = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
         d3_mid = h0_mid+xe_mid
#ifdef key_siggen 
         hc_sig_mid  =hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
#endif
         IF ( (d3_final > 0.0_rsh).AND.(d3_mid > 0.0_rsh) ) THEN
            particle % d3 = d3_final
            particle % xe = xe_final
            particle % h0 = h0_final
#ifdef key_siggen 
            particle % hc = hc_sig_final
#endif
         ELSE
            write(*,*) 'WARNING: Particule a la cote'
            write(*,*) d3_final,h0_final,d3_mid,h0_mid,n,m
            particle % xpos = x_old ! on remet dans l eau
            particle % ypos = y_old
            stop
         ENDIF

#ifdef key_MPI_2D
       ! check if we need to exchange or not
       ! -----------------------------------
         ipos = NINT(particle%xpos)
         jpos = NINT(particle%ypos)

       ! Here, we verify if the particle stays in your zone
       ! although it gets out of its zone; it will only move max to i+-1 and j+-1 
       ! thus we can compute the particle s new value in this time step without pb
       !
       ! If the particle went out the zone, we check on which side it went out and
       ! count it in the proper variable 'down_give', 'up_give', 'right_give' or 'left_give'
       ! (used to call ex_traj)
       !
       ! The neighbor domain where the particle went is coded in variable 'limitbye' according
       ! to this scheme :
       !         7  |  6  |  5
       !       -----+-----+-----
       !         8  |  0  |  4
       !       -----+-----+-----
       !         1  |  2  |  3

         IF ( jpos < ljmin ) THEN
            particle % limitbye = 2
            down_give = down_give+1
            IF ( ipos < limin ) particle % limitbye = 1
            IF ( ipos > limax ) particle % limitbye = 3
         ELSE IF ( jpos > ljmax ) THEN
            particle % limitbye = 6
            up_give = up_give+1
            IF ( ipos < limin ) particle % limitbye = 7
            IF ( ipos > limax ) particle % limitbye = 5
         ELSE IF ( ipos < limin ) THEN
            particle % limitbye = 8
            left_give = left_give+1
         ELSE IF ( ipos > limax ) THEN
            particle % limitbye = 4
            right_give = right_give+1
         END IF
#endif

      !**********************************
      ! Vertical advection
      ! is estimated from the environment at the averaged location x_mid,y_mid
      !**********************************    

         IF ( (d3_final <= 0.0_rsh).OR.(d3_mid <= 0.0_rsh) )  CYCLE

         IF ( particle%itypevert == 0 ) THEN

          !**********************************
          ! No vertical advection (float trajectory at constant immersion)
          !**********************************    
#ifdef key_siggen 
            particle % zpos = MIN(particle%zpos,d3_final) ! when particle wants to go in shallower areas than immersion 
            zlag = -particle%zpos + xe_final  ! switch from immersion to real z
            CALL ztosiggen(zlag,slag,xe_final,h0_final,hc_sig_final)
#else
            slag = -particle%zpos / d3_final  ! new sigma 
#endif 
            particle%spos = MAX(slag,-1.0_rsh) ! when particle wants to go in shallower areas than immersion  

         ELSE   

          !**************************************
          ! Vertical advection + random walk
          !**************************************

          ! vertical displacement induced by advection : ds_adv

            ds_adv = dtm*wint(s_old,xe_mid,h0_mid,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)

            IF ( particle%itypevert > 0 ) THEN ! Random walk (turbulent diffusion)

             !*********************************************************
             ! + Random walk (turbulent diffusion)
             !*********************************************************
             ! Switch to z coordinate for RW
             ! allow gain in computing time when in siggen
             ! because then no switch of coordinate in the RW loop

               DO l=0,kmax
#ifdef key_siggen              
                  CALL siggentoz(zpos(l),sigw(l),xe_mid,h0_mid,hc_sig_mid)
#else
                  zpos(l) = sigw(l)*d3_mid+xe_mid 
#endif    
               END DO

             ! smooth and spline the kz profile at x_mid,y_mid location
               CALL kzprofile(kzsmth,dkz2,zpos,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)

             ! vertical displacement induced by diffusion : ds_dif
             ! displacement linked to kz intensity and structure
             ! ( kz*(dkz/dz) see Visser)
             ! loop on t to use a smaller time step
             ! sous-boucle de calcul du deplacement vertical
             ! lie a kz + structure de  kz (dkz/dz)
             ! deplacement net du centre de masse (voir Visser) : ds_dif
             ! s_int positions intermediaires successives
             ! s_mid position mediane entre deplacement du centre de masse et s_int

#ifdef key_siggen              
               CALL siggentoz(z_old,s_old,xe_mid,h0_mid,hc_sig_mid)
#else
               z_old = s_old*d3_mid + xe_mid
#endif  
               z_int = z_old

               IF ( d3_mid > 5.0_rsh ) THEN   ! attention : lb+lt for boundary layers should
                                              ! be lower than criteria here
                  DO nt=1,ndtz

                   ! get the first derivative of Kz at z_int location
                     dksdz  = dksdzint(z_int,zpos,kzsmth,dkz2)                                       
                     ds_dif = dksdz*dtz                                                                     
            
                     z_mid = MIN(MAX(z_int+ds_dif*0.5_rsh,zpos(0)),zpos(kmax))
                   ! get the kz value at s_mid location entre son centre de masse et sa position initiale
                   ! ancien calcul         kzz=kzint(x_mid,y_mid,s_mid)    
                     CALL splint(zpos,kzsmth,dkz2,kmax+1,z_mid,kzz)
                     IF ( kzz < 0.0 ) THEN                   ! kzz should be > 0 but...
                        kzz = 0.0
                     ENDIF
            
                   ! random walk component of vertical diffusion (between 0 and 1)
                     CALL random_number(tir)
            
                   ! s_int=s_int+ds_dif+(2.0_rsh*tir-1.0_rsh)*sqrt(2.0_rsh*dtz*kzz/(r*d3_mid*d3_mid))  en sigma
                     z_int=z_int+ds_dif+(2.0_rsh*tir-1.0_rsh)*sqrt(2.0_rsh*dtz*kzz/r) 
                   ! handle boundary conditions
                   ! s_int=MIN(MAX(s_int,-1.0_rsh),0.0_rsh)         ! particle set on boundary
                   ! s_int=MIN(MAX(s_int,-(2_rsh+s_int)),-s_int)  ! reflective boundary condition
                   ! random mixed layer to avoid accumulation (see Ross and Sharples, 2004)
                   ! lb=0.5_rsh        ! distance in meters for the random boundary layer
                   ! lt=0.5_rsh
                   ! lts=lt/d3_mid
                   ! lbs=lb/d3_mid
                   ! CALL random_number(tir)
                   ! IF (s_int > -lts) s_int=-tir*lts
                   ! IF (s_int < -1_rsh+lbs) s_int=tir*lbs-1_rsh
             
                   ! handle boundary conditions
                   ! z_int=MIN(MAX(z_int,-h0_mid),xe_mid)         ! particle set on boundary
                   ! z_int=MIN(MAX(z_int,-(2_rsh*h0_mid+z_int)),2_rsh*xe_mid-z_int)  ! reflective boundary condition
                   ! random mixed layer to avoid accumulation (see Ross and Sharples, 2004)
                     lb=2.0_rsh        ! distance in meters for the random boundary layer
                     lt=2.0_rsh        ! twice the boundary layer should be ok
                     CALL random_number(tir)
                     IF (z_int > zpos(kmax)-lt) z_int=zpos(kmax)-tir*lt
                     IF (z_int < zpos(0)+lb) z_int=zpos(0)+tir*lb

                  END DO   ! end loop on RW time
         
               ENDIF  ! if d3_mid > 5.0

             ! back to sigma position
#ifdef key_siggen
               kw  = kmax
               kwm = 0
               CALL ksupkinf(z_int,zpos,kmax+1,kw,kwm,2)
               particle%spos = ((z_int-zpos(kwm))*sigw(kw)+ &
                                 (zpos(kw)-z_int)*sigw(kwm))/(zpos(kw)-zpos(kwm))  !! approximation
#else
               particle%spos = (z_int-xe_mid)/d3_mid                     
#endif

            ENDIF ! end of RW (itypevert > 0)

          !**************************************
          ! Vertical behavior (larval behavior)
          !**************************************

          ! Behaviour: speed of vertical migration
            wbio = 0.0_rsh ! may be done here (if simple) or in module IBM

!#ifdef key_ibm_behav
!                  CALL w_behav_part(n,m,sal,temp)
!                  wbio = particle%w
!#endif
            ds_bio = wbio*dtm/d3_mid  ! approximation
          ! End of vertical migration
          !**************************************

          ! Estimate of the new vertical position
            particle%spos = MIN(MAX(particle%spos+ds_adv+ds_bio,-1.0_rsh),0.0_rsh) 
!                  particle%spos = MIN(MAX(particle%spos+ds_bio,-1.0_rsh),0.0_rsh) 

#ifdef key_siggen

            CALL siggentoz(zposf,particle%spos,xe_final,h0_final,hc_sig_final)
            particle%zpos = -zposf+xe_final           ! immersion
#else
            particle%zpos = -particle%spos*d3_final   ! z is positive = immersion
#endif
         END IF     ! ends test on vertical movement (itypevert)

        ! if (abs(particle%zpos-test).gt.1.0_rsh) then
     !  !  if (n.eq.1 .and. m.eq.1) then
        !    write(*,*) 'vit vert', ds_adv, particle%spos, s_old
        !    write(*,*) test, particle%zpos
        !    write(*,*) abs(particle%zpos-test)
        !    write(*,*) n,m, particle%xpos, particle%ypos
        !    write(*,*) x_old, y_old
        !    write(*,*) d3avt, xeavt, hcavt,h0avt
        !    write(*,*) d3_final, xe_final, hc_sig_final,h0_final
       !  endif

      END DO   ! ends loop on number of particles

      patch => patch % next

   END DO   ! Ends loop on number of patches

   CALL_MPI ex_traj(down_give, up_give, right_give, left_give)

   PRINT_DBG*, 'exit TRAJ_3D routine'

   END SUBROUTINE traj_3d


   !!======================================================================

   SUBROUTINE avance(u,v,xe,deltat,xx,yy,s0,statp)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE avance  ***
   !&E
   !&E ** Purpose : calcul du deplacement cinematique d une particule
   !&E              sous l action d un courant et d une composante
   !&E              aleatoire de deplacement (random walk), dans
   !&E              le plan horizontal.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  08-2002 (F. Dumas, P. Lazure)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : pi,valmanq
   USE comvars2d,  ONLY : dy,dx,h0_g
   USE trajectools, ONLY: uint,vint

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )             :: s0
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: u,v
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rlg), INTENT( in )             :: deltat
   REAL(KIND=rsh), INTENT( inout )          :: xx,yy
   REAL(KIND=rsh), INTENT( inout ),OPTIONAL :: statp

   !! * Local declarations
   INTEGER                                  :: j0,jst,i0,ist,i1,j1
   REAL(KIND=rsh)                           :: uyp,ux0,uy0,xst,yst,a,b,tir1
   REAL(KIND=rsh)                           :: x0,y0,x1,y1,xp,yp,uxp,tir2
   REAL(KIND=rsh)                           :: dkx,disp
   REAL(KIND=rsh)                           :: dmaill,tetha

   !!----------------------------------------------------------------------
   !! * Executable part

   dmaill=0.01_rsh ! plus eleve plus jump en z, car advecte le long des sigmas, martin
   disp=0.0_rsh

   x0 = xx
   y0 = yy
   ux0=uint(x0,y0,s0,u,xe)
   uy0=vint(x0,y0,s0,v,xe)
   xp=x0+ux0*deltat/dx(nint(x0),nint(y0))
   yp=y0+uy0*deltat/dy(nint(x0),nint(y0))
   uxp=uint(xp,yp,s0,u,xe)
   uyp=vint(xp,yp,s0,v,xe)

   x1=x0+0.5_rsh*deltat*(ux0/dx(nint(x0),nint(y0))+uxp/dx(nint(xp),nint(yp)))
   y1=y0+0.5_rsh*deltat*(uy0/dy(nint(x0),nint(y0))+uyp/dy(nint(xp),nint(yp)))

   if (disp /= 0.0_rsh) then
      call random_number(tir1)
      tetha=2.0_rsh*pi*tir1
      call random_number(tir2)
      dkx=sqrt(2.0_rsh*disp*deltat*3.0_rsh)*tir2
      IF (0.5_rsh*(ux0+uxp) /= 0.0_rsh .and. 0.5_rsh*(uy0+uyp) /= 0.0_rsh) THEN
         x1=x1+dkx*cos(tetha)/dx(nint(x1),nint(y1))
         y1=y1+dkx*sin(tetha)/dy(nint(x1),nint(y1))
      END IF
   endif

   !     detection d un franchissement de mur hy
   !--------------------------------------------------

   j0 = nint(y0)
   j1 = nint(y1)
   IF (j0 /= j1) THEN
     jst = min (j0,j1)
     yst = real(jst) + 0.5_rsh
     xst = x0
     IF (x1 /= x0) THEN
       a = (y1-y0)/(x1-x0)
       b = y0 -a*x0
       xst = (yst-b)/a
     END IF
     ist = nint(xst)
     IF (v(1,ist,jst) == 0.0_rsh) THEN

     !------------------------------------------------------------
     !    real(j1-j0) est juste la pour donner le signe
     !    et savoir si l on s arrete au-dessus ou au-dessous
     !    du mur.
     !    Pourquoi xst n est-il pas modifie ?
     !    Reflection a la normale de la position finale sur terre
     !------------------------------------------------------------

       y1 = yst - dmaill*real(j1-j0)
       statp = statp+1
     END IF
   END IF

   !     detection d un franchissement de mur hx
   !--------------------------------------------------

   i0 = nint(x0)
   i1 = nint(x1)
   IF (i0 /= i1) THEN
     ist = min (i0,i1)
     xst = real(ist) + 0.5_rsh
     yst = y0
     IF (y1 /= y0) THEN
       a = (x1-x0)/(y1-y0)
       b = x0 -a*y0
       yst = (xst-b)/a
     END IF
     jst = nint(yst)
     IF (u(1,ist,jst) == 0.0_rsh) THEN
     !------------------------------------------------------------
     !    real(i1-i0) est juste la pour donner le signe
     !    et savoir si l on s arrete au-dessus ou au-dessous
     !    du mur.
     !    Pourquoi yst n est-il pas modifie ?
     !------------------------------------------------------------

       x1 = xst - dmaill*real(i1-i0)
       statp = statp+1
     END IF
   END IF


! detection d un franchissement limite de domaine
   IF (y1 <= 1.0_rsh.or.y1 >= real(jmax)) THEN
    !  x1 = 0.0_rsh  ! plante en MPI car fait passer d un CPU a un autre non adjacent
    !  y1 = 0.0_rsh  ! limitby n est plus coherent dans traject.
      statp=-valmanq ! on prefere flagger la particule et la laisser ou elle est (ne bougera plus) 
   END IF
   IF (x1 <= 1.0_rsh.or.x1 >= real(imax)) THEN
     ! x1 = 0.0_rsh
     ! y1 = 0.0_rsh
      statp=-valmanq ! on prefere flagger la particule et la laisser ou elle est (ne bougera plus) 
   END IF

   xx = x1
   yy = y1

   if (h0_g(NINT(xx),NINT(yy)) == -valmanq) then ! pour cas special ou passage exact sur coin de cellule, avec un mur unique partant de la section, tests precedents ne marche pas forcement
      xx=x0
      yy=y0
      statp=statp+1
   endif

   END SUBROUTINE avance

#endif

END MODULE
