!!****f* ariane/posini()
!! NAME
!!   posini() (posini.f90)
!!
!! USAGE
!!   
!!
!! FUNCTION
!!   
!! 
!! AUTHOR
!!   * Origin  : Bruno Blanke  (1992)
!!   * F77toF90: Nicolas Grima (April-May 2005)
!! 
!! CREATION DATE
!!   * April-May 2005
!!
!! HISTORY
!!   Date (dd/mm/yyyy/) - Modification(s)
!!
!! RESULT
!!   
!!
!! EXAMPLES
!!   
!!
!! NOTES
!!   ROBODoc header style.
!!
!! TODO
!!   
!!
!! PORTABILITY
!!         Machine-OS    - Fortran90/95 compiler
!!   * i686-pc-linux-gnu -         ifort
!!   * i686-pc-linux-gnu -          g95
!!
!! SEE ALSO
!!
!! TODO
!!   * Should be transform to a module.
!!   * Should be migrate to F90 !!!
!!   * 3 similar loops, try to find a common base to reduce writing.
!!   * more modularity.
!!   * more comments.
!!
!! USES
!!   * USE mod_precision
!!   * USE mod_namelist
!!   * USE mod_input_grid
!!   * USE mod_input_data
!!   * USE mod_posin
!!   * USE mod_flags
!!   * USE mod_stats
!!   * USE mod_stati
!!   * USE mod_txt
!!   * USE mod_fx
!!   * USE mod_fy
!!   * USE mod_fz
!!   * USE mod_criter0
!!   * USE mod_sigma
!!   * USE mod_zinter
!!
!! USED BY
!!   * trajec
!!
!! SOURCE
!!=========================================================================

SUBROUTINE posini(nstat, trmin, maxsect,t0,nsegm,ntraj,nsect,trtot,nfnt,ibin) 2,172
  !-----------------------------------------------------------------------
  ! automatic positioning: particles are gathered where the transport is
  ! the largest, giving an equivalent individual tranport to all
  ! particles (not to exceed t0)
  !------------------------------------------------------------------------

  !------------------!
  ! USE ASSOCIAITION !
  !------------------!
  USE mod_precision
  USE mod_namelist
  USE mod_input_grid
  USE mod_input_data
  USE mod_posin
  USE mod_flags
  USE mod_stats
  USE mod_stati
  USE mod_txt
  USE mod_fx
  USE mod_fy
  USE mod_fz
  USE mod_criter0
  USE mod_sigma
  USE mod_zinter
  USE mod_reducmem
  USE mod_lun

  !-------------!
  ! Declaration !
  !-------------!
  IMPLICIT NONE

  !- arguments -!
  INTEGER(kind=iprec), INTENT(in)  :: ibin, nstat
  INTEGER(kind=iprec), INTENT(out) :: maxsect,nsegm, ntraj, nsect, nfnt
  REAL(kind=rprec)   , INTENT(in)  :: t0, trmin
  REAL(kind=rprec)   , INTENT(out) :: trtot

  !- local variables -!
  INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE :: &
       it1, & !
       it2, & !
       jt1, & ! 
       jt2, & !
       kt1, & !
       kt2, & !
       num, & !
       ndir   !

  INTEGER(kind=iprec) :: n     !
  INTEGER(kind=iprec) :: maxsegm, idum     !
  INTEGER(kind=iprec) :: i,j,k
  INTEGER(kind=iprec) :: is,js,ks
  INTEGER(kind=iprec) :: it0s,jt0s,kt0s
  INTEGER(kind=iprec) :: iu0s,jv0s,kw0s
  INTEGER(kind=iprec) :: l     !
  INTEGER(kind=iprec) :: iu0   !
  INTEGER(kind=iprec) :: it0   !
  INTEGER(kind=iprec) :: it00  !
  INTEGER(kind=iprec) :: n1    !
  INTEGER(kind=iprec) :: n2    !
  INTEGER(kind=iprec) :: i1    !
  INTEGER(kind=iprec) :: j1    !
  INTEGER(kind=iprec) :: k1    !
  INTEGER(kind=iprec) :: l1    !
  INTEGER(kind=iprec) :: jt0   !
  INTEGER(kind=iprec) :: jt00  !
  INTEGER(kind=iprec) :: ilim1 !
  INTEGER(kind=iprec) :: ilim2 !
  INTEGER(kind=iprec) :: jv0   !
  INTEGER(kind=iprec) :: kt0   !
  INTEGER(kind=iprec) :: kt00  !
  INTEGER(kind=iprec) :: kw0   !

  REAL(kind=rprec) :: t  ! temperature
  REAL(kind=rprec) :: s  ! salinity
  REAL(kind=rprec) :: r  ! density
  REAL(kind=rprec) :: u  ! zonal current
  REAL(kind=rprec) :: v  ! meridional current
  REAL(kind=rprec) :: w  ! vertical current
  REAL(kind=rprec) :: tr ! trajectory
  REAL(kind=rprec) :: prof ! W depth

  !-------------!
  ! Code begins !
  !-------------!
  !- Open the file param0 and compute the number of segments -!
  !- and the number of sections.                             -!
  !- (remember that a section could be composed by more than -!
  !-  one segment)                                           -!
  maxsegm = 1
  maxsect = 1
  OPEN(unit=lun_dummy, file='param0', action="read")
  DO 
    READ(lun_dummy,*, END=7777, ERR=7777) idum
    maxsegm = maxsegm + 1
    IF (idum > maxsect) THEN
      maxsect = maxsect + 1
    ENDIF
  ENDDO

7777 continue
  CLOSE(lun_dummy)

  !- Dynamic Allocations -!

  CALL sub_txt_alloc(maxsect, nstat)
  CALL sub_txt_init()

  CALL sub_flags_alloc(imt, jmt, kmt)
  CALL sub_flags_init(0)

  CALL sub_stati_alloc(maxsect)
  CALL sub_stati_init(0)

  CALL sub_stats_alloc(nstat, maxsect) 
  CALL sub_stats_init()

  ALLOCATE(it1(maxsegm))
  ALLOCATE(it2(maxsegm))
  ALLOCATE(jt1(maxsegm))
  ALLOCATE(jt2(maxsegm))
  ALLOCATE(kt1(maxsegm))
  ALLOCATE(kt2(maxsegm))
  ALLOCATE(num(maxsegm))
  ALLOCATE(ndir(maxsegm))

  !- Initializations -!
  ntraj = 0
  nfnt  = 0
  nsect = 0
  trtot = 0.

  OPEN(unit=lun_dummy, file='param0', action="read")
  !- DoLoop on segment number -!
  DO n=1,maxsegm+1

    READ(lun_dummy,*,END=7054,ERR=7054) &
         num(n),it1(n),it2(n),jt1(n),jt2(n),kt1(n),kt2(n),sname
    ndir(n)=0
    WRITE(58,*)num(n),' ',it1(n),' ',it2(n),' ',jt1(n),' ',jt2(n), &
         ' ',kt1(n),' ',kt2(n),' ',sname

    IF (num(n) < 0) THEN
      WRITE(58,*)'transparent section associated with segment #',n
      nfnt=max0(nfnt,-num(n))
    ENDIF

    IF (ABS(num(n)) > maxsect) THEN
      WRITE(0,*)'problem with segment #',n
      WRITE(0,*)'      maximum number of sections is: ',maxsect
      STOP
    ENDIF

    IF ((it1(n).NE.it2(n)).AND. &
         (jt1(n).NE.jt2(n)).AND. &
         (kt1(n).NE.kt2(n))) THEN
      WRITE(0,*)'problem with segment #',n
      WRITE(0,*)'      one pair of indices must be identical'
      STOP
    ENDIF

    IF (  ( (it1(n) == it2(n)).AND.(jt1(n) == jt2(n)) ).OR. &
         ( (it1(n) == it2(n)).AND.(kt1(n) == kt2(n)) ).OR. &
         ( (jt1(n) == jt2(n)).AND.(kt1(n) == kt2(n)) )  ) THEN
      WRITE(0,*)'problem with segment #',n
      WRITE(0,*)'      only one pair of indices can be identical'
      STOP
    ENDIF

    IF (num(n) > 0) THEN
      secname(num(n))=sname
      IF (num(n) > nsect) nsect=num(n)
    ELSE
      fntname(-num(n))=sname
    ENDIF

    IF (it1(n) == it2(n)) THEN
      ndir(n)=1
      i=it1(n)
      IF (i < 0) THEN
        i=-i
        ndir(n)=-1
      ENDIF
      DO k=kt1(n),kt2(n)
        DO j=jt1(n),jt2(n)
          mtfin(i,j,k)=num(n)
        END DO
      END DO
    ENDIF
    IF (it1(n) < 0) it1(n)=-it1(n)
    IF (it2(n) < 0) it2(n)=-it2(n)

    IF (jt1(n) == jt2(n)) THEN
      ndir(n)=2
      j=jt1(n)
      IF (j < 0) THEN
        j=-j
        ndir(n)=-2
      ENDIF
      DO k=kt1(n),kt2(n)
        DO i=it1(n),it2(n)
          mtfin(i,j,k)=num(n)
        END DO
      END DO
    ENDIF
    IF (jt1(n) < 0) jt1(n)=-jt1(n)
    IF (jt2(n) < 0) jt2(n)=-jt2(n)

    IF (kt1(n) == kt2(n)) THEN
      ndir(n)=3
      k=kt1(n)
      IF (k < 0) THEN
        k=-k
        ndir(n)=-3
      ENDIF
      DO j=jt1(n),jt2(n)
        DO i=it1(n),it2(n)
          mtfin(i,j,k)=num(n)
        END DO
      END DO
    ENDIF
    IF (kt1(n) < 0) kt1(n)=-kt1(n)
    IF (kt2(n) < 0) kt2(n)=-kt2(n)

  END DO

  WRITE(0,*)'# of segments to be read exceeds maxsegm'
  STOP

7054 nsegm=n-1
  WRITE(0,*)'# of segments that could be read: ',nsegm

  secname(0)='meanders'

  WRITE(58,*)'nsegm= ',nsegm

  IF (ibin < 1) THEN

    WRITE(0,*)'initial positions are automatically computed'

    DO n = 1, nsegm
      IF (num(n) == 1) THEN

        IF ((ndir(n) == 1).OR.(ndir(n) == -1)) THEN
          it0=it1(n)
          iu0=it0                !NG  Performance !
          it00=it0+1             !
          IF (ndir(n) < 0) THEN  ! iu0  = it0 +  INT( ( 1.- REAL(ndir(n), kind=rprec) ) * 0.5)
            iu0=it0-1            ! it00 = it0 + SIGN(1, ndir(n))
            it00=it0-1           !
          ENDIF                  !NG
          DO l=1,lmt
            DO k=kt1(n),kt2(n)
              DO j=jt1(n),jt2(n)

                CALL sub_shift_or_not_ind(it0, j, k, it0s, js, ks)

                IF ( (tmask(it0s,js,ks,1) > 0.5).AND.(mtfin(it00,j,k) <= 0) ) THEN
                  IF (key_nointerpolstats) THEN
                    IF (key_alltracers) THEN
                      t=tt(it0s,js,ks,l)
                      s=ss(it0s,js,ks,l)
                      r=rr(it0s,js,ks,l)
                    ENDIF
                  ENDIF

                  CALL sub_shift_or_not_ind(iu0, j, k, iu0s, js, ks)

                  u=uu(iu0s,js,ks,l)
                  !----------------------------------
                  ! test on the sign of the velocity
                  !----------------------------------
                  IF ((ABS(u) > trmin).AND.(float(ndir(n))*u > 0)) THEN

                    IF (lmt == 1) THEN
                      n1=1+INT( SQRT(ABS(u)/t0) )
                      n2=1
                    ELSE
                      n1=1+INT( (ABS(u)/t0)**(1./3.) )
                      n2=n1
                    ENDIF

                    tr=ABS(u)/n1/n1/n2

                    DO j1=1,n1
                      DO k1=1,n1
                        DO l1=1,n2
                          ntraj=ntraj+1
                          IF (ntraj > nmax) THEN
                            WRITE(0,*)'too many particles, maximum is: ',nmax
                            STOP
                          ENDIF
                          tfi(ntraj)=iu0
                          tfj(ntraj)=j-1.-.5/n1+float(j1)/n1
                          tfk(ntraj)=k-.5/n1+float(k1)/n1
                          tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
                          ttr(ntraj)=tr

                          IF (criter0(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
                               tfl(ntraj),it0,j,k,l)) THEN

                            trtot=trtot+tr
                            isn(-1)=isn(-1)+1
                            sn(-1)=sn(-1)+tr
                            IF (key_alltracers) THEN
                              IF (.NOT.key_nointerpolstats) THEN
                                t=zinter &
                                     (tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                s=zinter &
                                     (ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                IF (key_approximatesigma) THEN
                                  r=zinter &
                                       (rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                ELSE
                                  r=sigma(zsigma,s,t)
                                ENDIF
                              ENDIF
                              IF (t <= s0min(5,-1)) s0min(5,-1)=t
                              IF (s <= s0min(6,-1)) s0min(6,-1)=s
                              IF (r <= s0min(7,-1)) s0min(7,-1)=r
                              IF (t >= s0max(5,-1)) s0max(5,-1)=t
                              IF (s >= s0max(6,-1)) s0max(6,-1)=s
                              IF (r >= s0max(7,-1)) s0max(7,-1)=r
                              s0x(5,-1)=s0x(5,-1)+t*tr
                              s0x(6,-1)=s0x(6,-1)+s*tr
                              s0x(7,-1)=s0x(7,-1)+r*tr
                              s0x2(5,-1)=s0x2(5,-1)+t*t*tr
                              s0x2(6,-1)=s0x2(6,-1)+s*s*tr
                              s0x2(7,-1)=s0x2(7,-1)+r*r*tr

!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG                                     t,s,r
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
                                     t,s,r
!!NG                              ENDIF
                            ELSE
!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG                              ENDIF
                            ENDIF

                            IF (key_roms) THEN
                              prof=fz(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)
                            ELSE
                              prof=fz(gk=-tfk(ntraj))
                            ENDIF

                            IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
                                 s0min(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
                                 s0min(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof <= s0min(3,-1)) &
                                 s0min(3,-1)=-prof
                            IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
                                 s0max(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
                                 s0max(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof >= s0max(3,-1)) &
                                 s0max(3,-1)=-prof
                            s0x(1,-1)=s0x(1,-1)+fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x(2,-1)=s0x(2,-1)+fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x(3,-1)=s0x(3,-1)- prof *tr
                            s0x2(1,-1)=s0x2(1,-1)+fx(tfi(ntraj),tfj(ntraj))* &
                                 fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(2,-1)=s0x2(2,-1)+fy(tfi(ntraj),tfj(ntraj))* &
                                 fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(3,-1)=s0x2(3,-1)+ prof * &
                                 prof *tr

                          ELSE

                            ntraj=ntraj-1

                          ENDIF

                        END DO
                      END DO
                    END DO

7050                FORMAT(4(1x,f7.3),1x,f12.3,3(1x,f6.3))
                    !
                    ! fin de test sur signe de la vitesse
                  ENDIF
                  !
                  ! fin de test sur masque/critere
                ENDIF
                !
                ! fins de boucle sur j,k,l
              END DO
            END DO
          END DO
          !
          ! fin test it1=it2
        ENDIF
        !
        IF ((ndir(n) == 2).OR.(ndir(n) == -2)) THEN
          jt0=jt1(n)
          jv0=jt0
          jt00=jt0+1
          IF (ndir(n) < 0) THEN
            jv0=jt0-1
            jt00=jt0-1
          ENDIF
          DO l=1,lmt
            DO k=kt1(n),kt2(n)
              IF (key_periodic) THEN
                ilim1=MAX(2,it1(n))
                ilim2=MIN(imt-1,it2(n))
              ELSE
                ilim1=it1(n)
                ilim2=it2(n)
              ENDIF
              DO i=ilim1, ilim2

                CALL sub_shift_or_not_ind(i, jt0, k, is, jt0s, ks)

                IF ( (tmask(is,jt0s,ks,1) > 0.5).AND.(mtfin(i,jt00,k) <= 0) ) THEN
                  IF (key_nointerpolstats) THEN
                    IF (key_alltracers) THEN
                      t=tt(is,jt0s,ks,l)
                      s=ss(is,jt0s,ks,l)
                      r=rr(is,jt0s,ks,l)
                    ENDIF
                  ENDIF

                  CALL sub_shift_or_not_ind(i, jv0, k, is, jv0s, ks)

                  v=vv(is,jv0s,ks,l)
                  !----------------------------------
                  ! test on the sign of the velocity
                  !----------------------------------
                  IF ((ABS(v) > trmin).AND.(float(ndir(n))*v > 0)) THEN

                    IF (lmt == 1) THEN
                      n1=1+INT( SQRT(ABS(v)/t0) )
                      n2=1
                    ELSE
                      n1=1+INT( (ABS(v)/t0)**(1./3.) )
                      n2=n1
                    ENDIF

                    tr=ABS(v)/n1/n1/n2

                    DO i1=1,n1
                      DO k1=1,n1
                        DO l1=1,n2
                          ntraj=ntraj+1
                          IF (ntraj > nmax) THEN
                            WRITE(0,*)'too many particles, maximum is: ',nmax
                            STOP
                          ENDIF
                          tfi(ntraj)=i-1.-.5/n1+float(i1)/n1
                          tfj(ntraj)=jv0
                          tfk(ntraj)=k-.5/n1+float(k1)/n1
                          tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
                          ttr(ntraj)=tr

                          IF (criter0(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
                               tfl(ntraj),i,jt0,k,l)) THEN

                            trtot=trtot+tr
                            isn(-1)=isn(-1)+1
                            sn(-1)=sn(-1)+tr
                            IF (key_alltracers) THEN
                              IF (.NOT.key_nointerpolstats) THEN
                                t=zinter &
                                     (tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                s=zinter &
                                     (ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                IF (key_approximatesigma) THEN
                                  r=zinter &
                                       (rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                ELSE
                                  r=sigma(zsigma,s,t)
                                ENDIF
                              ENDIF
                              IF (t <= s0min(5,-1)) s0min(5,-1)=t
                              IF (s <= s0min(6,-1)) s0min(6,-1)=s
                              IF (r <= s0min(7,-1)) s0min(7,-1)=r
                              IF (t >= s0max(5,-1)) s0max(5,-1)=t
                              IF (s >= s0max(6,-1)) s0max(6,-1)=s
                              IF (r >= s0max(7,-1)) s0max(7,-1)=r
                              s0x(5,-1)=s0x(5,-1)+t*tr
                              s0x(6,-1)=s0x(6,-1)+s*tr
                              s0x(7,-1)=s0x(7,-1)+r*tr
                              s0x2(5,-1)=s0x2(5,-1)+t*t*tr
                              s0x2(6,-1)=s0x2(6,-1)+s*s*tr
                              s0x2(7,-1)=s0x2(7,-1)+r*r*tr

!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG                                     t,s,r
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
                                     t,s,r
!!NG                              ENDIF
                            ELSE
!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG                              ENDIF
                            ENDIF

                            IF (key_roms) THEN
                              prof=fz(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)
                            ELSE
                              prof=fz(gk=-tfk(ntraj))
                            ENDIF

                            IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
                                 s0min(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
                                 s0min(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof <= s0min(3,-1)) &
                                 s0min(3,-1)=-prof
                            IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
                                 s0max(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
                                 s0max(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof >= s0max(3,-1)) &
                                 s0max(3,-1)=-prof
                            s0x(1,-1)=s0x(1,-1)+fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x(2,-1)=s0x(2,-1)+fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x(3,-1)=s0x(3,-1)- prof *tr
                            s0x2(1,-1)=s0x2(1,-1)+fx(tfi(ntraj),tfj(ntraj))* &
                                 fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(2,-1)=s0x2(2,-1)+fy(tfi(ntraj),tfj(ntraj))* &
                                 fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(3,-1)=s0x2(3,-1)+ prof * &
                                 prof *tr

                          ELSE

                            ntraj=ntraj-1

                          ENDIF

                        END DO
                      END DO
                    END DO
                    !
                    ! fin de test sur signe de la vitesse
                  ENDIF
                  !
                  ! fin de test sur masque/critere
                ENDIF
                !
                ! fins de boucle sur i,k,l
              END DO
            END DO
          END DO
          !
          ! fin test jt1=jt2
        ENDIF
        !
        IF ((ndir(n) == 3).OR.(ndir(n) == -3)) THEN
          kt0=kt1(n)
          kw0=kt0+1
          kt00=kt0+1
          IF (ndir(n) < 0) THEN
            kw0=kt0
            kt00=kt0-1
          ENDIF
          IF (kt0 == 0) THEN
            kt0=1
          ENDIF
          DO l=1,lmt
            DO j=jt1(n),jt2(n)
              DO i=it1(n),it2(n)

                CALL sub_shift_or_not_ind(i, j, kt0, is, js, kt0s)

                IF ( (tmask(is,js,kt0s,1) > 0.5).AND.(mtfin(i,j,kt00) <= 0) ) THEN
                  IF (key_nointerpolstats) THEN
                    IF (key_alltracers) THEN

                      t=tt(is,js,kt0s,l)
                      s=ss(is,js,kt0s,l)
                      r=rr(is,js,kt0s,l)

                    ENDIF
                  ENDIF

                  CALL sub_shift_or_not_ind(i, j, kw0, is, js, kw0s)

                  w=ww(is,js,kw0s,l)
                  !----------------------------------
                  ! test on the sign of the velocity
                  !----------------------------------
                  IF ((ABS(w) > trmin).AND.(float(ndir(n))*w < 0)) THEN

                    IF (lmt == 1) THEN
                      n1=1+INT( SQRT(ABS(w)/t0) )
                      n2=1
                    ELSE
                      n1=1+INT( (ABS(w)/t0)**(1./3.) )
                      n2=n1
                    ENDIF

                    tr=ABS(w)/n1/n1/n2

                    DO j1=1,n1
                      DO i1=1,n1
                        DO l1=1,n2
                          ntraj=ntraj+1
                          IF (ntraj > nmax) THEN
                            WRITE(0,*)'too many particles, maximum is: ',nmax
                            STOP
                          ENDIF
                          tfi(ntraj)=i-1.-.5/n1+float(i1)/n1
                          tfj(ntraj)=j-1.-.5/n1+float(j1)/n1
                          tfk(ntraj)=kw0
                          tfl(ntraj)=l-.5-.5/n2+float(l1)/n2
                          ttr(ntraj)=tr

                          IF (criter0(tfi(ntraj),tfj(ntraj),-tfk(ntraj), &
                               tfl(ntraj),i,j,kt0,l)) THEN

                            trtot=trtot+tr
                            isn(-1)=isn(-1)+1
                            sn(-1)=sn(-1)+tr
                            IF (key_alltracers) THEN
                              IF (.NOT.key_nointerpolstats) THEN
                                t=zinter &
                                     (tt,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                s=zinter &
                                     (ss,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                IF (key_approximatesigma) THEN
                                  r=zinter &
                                       (rr,tfi(ntraj),tfj(ntraj),-tfk(ntraj),tfl(ntraj))
                                ELSE
                                  r=sigma(zsigma,s,t)
                                ENDIF
                              ENDIF
                              IF (t <= s0min(5,-1)) s0min(5,-1)=t
                              IF (s <= s0min(6,-1)) s0min(6,-1)=s
                              IF (r <= s0min(7,-1)) s0min(7,-1)=r
                              IF (t >= s0max(5,-1)) s0max(5,-1)=t
                              IF (s >= s0max(6,-1)) s0max(6,-1)=s
                              IF (r >= s0max(7,-1)) s0max(7,-1)=r
                              s0x(5,-1)=s0x(5,-1)+t*tr
                              s0x(6,-1)=s0x(6,-1)+s*tr
                              s0x(7,-1)=s0x(7,-1)+r*tr
                              s0x2(5,-1)=s0x2(5,-1)+t*t*tr
                              s0x2(6,-1)=s0x2(6,-1)+s*s*tr
                              s0x2(7,-1)=s0x2(7,-1)+r*r*tr
!!NG
!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr, &
!!NG                                     t,s,r
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr, &
                                     t,s,r
!!NG                              ENDIF
                            ELSE
!!NG                              IF (key_outgridt) THEN
!!NG                                WRITE(50,7050) &
!!NG                                     tfi(ntraj)+.5,tfj(ntraj)+.5,tfk(ntraj)-.5,tfl(ntraj),tr
!!NG                              ELSE
                                WRITE(50,7050) &
                                     tfi(ntraj),tfj(ntraj),tfk(ntraj),tfl(ntraj),tr
!!NG                              ENDIF
                            ENDIF

                            IF (key_roms) THEN

                              prof=fz(gi=tfi(ntraj), gj=tfj(ntraj), gk=-tfk(ntraj), il=l)

                            ELSE
                              prof=fz(gk=-tfk(ntraj))
                            ENDIF

                            IF (fx(tfi(ntraj),tfj(ntraj)) <= s0min(1,-1)) &
                                 s0min(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) <= s0min(2,-1)) &
                                 s0min(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof <= s0min(3,-1)) &
                                 s0min(3,-1)=-prof
                            IF (fx(tfi(ntraj),tfj(ntraj)) >= s0max(1,-1)) &
                                 s0max(1,-1)=fx(tfi(ntraj),tfj(ntraj))
                            IF (fy(tfi(ntraj),tfj(ntraj)) >= s0max(2,-1)) &
                                 s0max(2,-1)=fy(tfi(ntraj),tfj(ntraj))
                            IF (-prof>= s0max(3,-1)) &
                                 s0max(3,-1)=-prof
                            s0x(1,-1)=s0x(1,-1)+fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x(2,-1)=s0x(2,-1)+fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x(3,-1)=s0x(3,-1)- prof *tr
                            s0x2(1,-1)=s0x2(1,-1)+fx(tfi(ntraj),tfj(ntraj))* &
                                 fx(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(2,-1)=s0x2(2,-1)+fy(tfi(ntraj),tfj(ntraj))* &
                                 fy(tfi(ntraj),tfj(ntraj))*tr
                            s0x2(3,-1)=s0x2(3,-1)+ prof * &
                                 prof *tr

                          ELSE

                            ntraj=ntraj-1

                          ENDIF

                        END DO
                      END DO
                    END DO
                    !
                    ! fin de test sur signe de la vitesse
                  ENDIF
                  !
                  ! fin de test sur masque/critere
                ENDIF
                !
                ! fins de boucle sur i,j,l
              END DO
            END DO
          END DO
          !
          ! fin test kt1=kt2
        ENDIF
        !
        ! fin test num(n)
      ENDIF
      !
      ! fin du positionnement
    END DO

  ELSE
    WRITE(0,*) ' posini: ibin =', ibin
  ENDIF

  !- Deallocate dynamic memory array
  DEALLOCATE(it1)
  DEALLOCATE(it2)
  DEALLOCATE(jt1)
  DEALLOCATE(jt2)
  DEALLOCATE(kt1)
  DEALLOCATE(kt2)
  DEALLOCATE(num)
  DEALLOCATE(ndir)

  RETURN
END SUBROUTINE posini
!!***