!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! - Ariane - (May - 2007)
!! 
!! bruno.blanke@univ-brest.fr and nicolas.grima@univ-brest.fr
!! 
!! This software is a computer program whose purpose is 
!! the computation of 3D streamlines in a given velocity field 
!! (as the output of an Ocean General Circulation Model) and 
!! subsequent water masses analyses.
!! 
!! This software is governed by the CeCILL license under French law and
!! abiding by the rules of distribution of free software.  You can  use, 
!! modify and/ or redistribute the software under the terms of the CeCILL
!! license as circulated by CEA, CNRS and INRIA at the following URL
!! "http://www.cecill.info". 
!! 
!! As a counterpart to the access to the source code and  rights to copy,
!! modify and redistribute granted by the license, users are provided only
!! with a limited warranty  and the software's author,  the holder of the
!! economic rights,  and the successive licensors  have only  limited
!! liability. 
!! 
!! In this respect, the user's attention is drawn to the risks associated
!! with loading,  using,  modifying and/or developing or reproducing the
!! software by the user in light of its specific status of free software,
!! that may mean  that it is complicated to manipulate,  and  that  also
!! therefore means  that it is reserved for developers  and  experienced
!! professionals having in-depth computer knowledge. Users are therefore
!! encouraged to load and test the software's suitability as regards their
!! requirements in conditions enabling the security of their systems and/or 
!! data to be ensured and,  more generally, to use and operate it in the 
!! same conditions as regards security. 
!! 
!! The fact that you are presently reading this means that you have had
!! knowledge of the CeCILL license and that you accept its terms.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!****f* ariane/posini_seq()
!! NAME
!!   posini_seq() (posini-seq.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_seq(trmin, nb_traj, nb_sect, trans_total, pos_nfnt) 1,107
  !-----------------------------------------------------------------------
  ! automatic positioning: particles are gathered where the transport is
  ! the largest, giving an equivalent individual tranport to all
  ! particles (not to exceed max_transport)
  !------------------------------------------------------------------------

  !------------------!
  ! USE ASSOCIAITION !
  !------------------!
  USE mod_precision
  USE mod_namelist
  USE mod_input_grid
  USE mod_seq
  USE mod_input_data
  USE mod_init_particules
  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(out) :: nb_traj, nb_sect, pos_nfnt
  REAL(kind=rprec)   , INTENT(in)  :: trmin
  REAL(kind=rprec)   , INTENT(out) :: trans_total

  INTEGER(kind=iprec) :: &
       born1, & !
       born2, &
       alloc_size, & !
       n         , & !
       i         , & !
       k         , & !
       j         , & !
       l         , & !
       it0       , & !
       iu0       , & !
       it00      , & !
       n1        , & !
       n2        , & !
       j1        , & !
       k1        , & !
       l1        , & !
       jt0       , & !
       jv0       , & !
       jt00      , & !
       i1        , & !
       kt0       , & !
       kw0       , & !
       kt00      , & !
       is,js,ks  , & !
       it0s,jt0s,kt0s , & !
       iu0s,jv0s,kw0s 
  !
  INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE :: &
       num , & !
       it1, & !
       it2, & !
       jt1, & !
       jt2, & !
       kt1, & !
       kt2, & !
       ndir

  REAL(kind=rprec) :: u  ! zonal current
  REAL(kind=rprec) :: v  ! meridional current
  REAL(kind=rprec) :: w  ! vertical current

  REAL(kind=rprec) :: init_trans ! transport

  REAL(kind=rprec) :: prof ! depth

  LOGICAL :: switch=.FALSE.

  !- Comments -!
  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'  ---------------------------'
  WRITE(lun_standard,*)'  = ENTER >>>> POSINI (SEQ) ='
  WRITE(lun_standard,*)'  ---------------------------'

  !-------------!
  ! Code begins !
  !-------------!

  !- Dynamic Allocations -!

  CALL sub_txt_alloc()
  CALL sub_txt_init()

  CALL sub_flags_alloc(imt, jmt, kmt)
  CALL sub_flags_init(-1)

  CALL sub_stati_alloc()
  CALL sub_stati_init(0)

  CALL sub_stats_alloc() 
  CALL sub_stats_init()

  ALLOCATE(num(maxsegm)) ; num(:)  = 0
  CALL sub_memory(size(num),'i','num','posini_seq')
  ALLOCATE(it1(maxsegm)) ; it1(:)  = 0
  CALL sub_memory(size(it1),'i','it1','posini_seq')
  ALLOCATE(it2(maxsegm)) ; it2(:)  = 0
  CALL sub_memory(size(it2),'i','it2','posini_seq')
  ALLOCATE(jt1(maxsegm)) ; jt1(:)  = 0
  CALL sub_memory(size(jt1),'i','jt1','posini_seq')
  ALLOCATE(jt2(maxsegm)) ; jt2(:)  = 0
  CALL sub_memory(size(jt2),'i','jt2','posini_seq')
  ALLOCATE(kt1(maxsegm)) ; kt1(:)  = 0
  CALL sub_memory(size(kt1),'i','kt1','posini_seq')
  ALLOCATE(kt2(maxsegm)) ; kt2(:)  = 0
  CALL sub_memory(size(kt2),'i','kt2','posini_seq')
  ALLOCATE(ndir(maxsegm)); ndir(:) = 0
  CALL sub_memory(size(ndir),'i','ndir','posini_seq')

  !- Initializations -!
  nb_traj     = 0 
  pos_nfnt    = 0
  nb_sect     = 0
  trans_total = 0._rprec

  !- Open sections.txt file -!
  OPEN(unit=lun_dummy, file='sections.txt', action="read", &
       POSITION='rewind')

  WRITE(lun_standard,*)'--- sections.txt is opened: '

  !- DoLoop on segment number -!
  DO n = 1, maxsegm

     READ(lun_dummy,*) &
          num(n),it1(n),it2(n),jt1(n),jt2(n),kt1(n),kt2(n),sname

     WRITE(lun_standard,*)'    --- Reading section: ', sname

     ndir(n) = 0
     
     WRITE(lun_output,*)         &
          num(n),' '           , &
          it1(n),' ',it2(n),' ', &
          jt1(n),' ',jt2(n),' ', &
          kt1(n),' ',kt2(n),' ',sname

     IF (num(n) < 0) THEN
        WRITE(lun_error,*)'    transparent sections are not allowed'
        STOP
     ENDIF

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

     IF ( (it1(n) /= it2(n)).AND. &
          (jt1(n) /= jt2(n)).AND. &
          (kt1(n) /= kt2(n))) THEN
        WRITE(lun_error,*)'    problem with segment #',n
        WRITE(lun_error,*)'      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(lun_error,*)'    problem with segment #',n
        WRITE(lun_error,*)'      only one pair of indices can be identical'
        STOP
     ENDIF

     secname(num(n)) = sname

     IF (num(n) > nb_sect) nb_sect = num(n)


     !-----------------------!
     !- Latitudinal section -!
     !-----------------------!
     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)

     !------------------------!
     !- Longitudinal section -!
     !------------------------!
     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)


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

  secname(0)='meanders'

  !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!
  !!===============================================================================!!
  !!_______________________________________________________________________________!!

  IF (TRIM(bin) == 'nobin') THEN

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

     !!NG: to have the results in  same order than the old version of Ariane
     IF (TRIM(forback) == 'backward') THEN
        forback='forward'
        switch=.TRUE.
     ELSE
        switch=.FALSE.
     ENDIF

     IF (key_alltracers) THEN
        IF ((key_roms).OR.(key_symphonie))THEN
           alloc_size = 7
        ELSE
           IF (TRIM(w_surf_option) == 'E-P-R') THEN
              alloc_size = 7
           ELSE
              alloc_size = 6
           ENDIF
        ENDIF
     ELSE
        IF ((key_roms).OR.(key_symphonie))THEN
           alloc_size = 3
        ELSE
           IF (TRIM(w_surf_option) == 'E-P-R') THEN
              alloc_size = 4
           ELSE
              alloc_size = 3
           ENDIF
        ENDIF
     ENDIF

     CALL sub_seq_alloc(alloc_size)

     IF ( (lmin /= 1).AND.(TRIM(forback) == 'forward')) THEN
        CALL sub_seq_init(lmin)
     ELSEIF (TRIM(forback) == 'backward') THEN  !! Never used !!
        CALL sub_seq_init(lmax)
     ELSE
        CALL sub_seq_init()
     ENDIF

     DO l = lmin, lmax

        !! READ INPUT DATA SEQUENTIALLY !!
        IF (key_roms) THEN
           CALL sub_input_data_seq_roms( &
                ncids(:)               , &
                varids(:)              , &
                new_file(:)            , &
                ind_file(:)            , &
                ind_time(:)            , &
                ind_time_size(:)       , &
                sdimsorders(:,:)        )
        ELSEIF(key_symphonie) THEN
           CALL sub_input_data_seq_symphonie( &
                ncids(:)                    , &
                varids(:)                   , &
                new_file(:)                 , &
                ind_file(:)                 , &
                ind_time(:)                 , &
                ind_time_size(:)            , &
                sdimsorders(:,:)              )

        ELSEIF (key_B2C_grid) THEN

          CALL sub_input_data_B2C_gridz_seq( &
                ncids(:)              , &
                varids(:)             , &
                new_file(:)           , &
                ind_file(:)           , &
                ind_time(:)           , &
                ind_time_size(:)      , &
                sdimsorders(:,:)       )

          !!NG: 16/04/2009 WRITE(lun_error,*)''
          !!NG: 16/04/2009 WRITE(lun_error,*)'STOP === key_B2C_grid in sequential mode is not available === STOP'
          !!NG: 16/04/2009 STOP

       ELSE
          CALL sub_input_data_seq_opa( &
                ncids(:)              , &
                varids(:)             , &
                new_file(:)           , &
                ind_file(:)           , &
                ind_time(:)           , &
                ind_time_size(:)      , &
                sdimsorders(:,:)       )
        ENDIF

        !!NG: Bug fixed 27/12/2006
        !!NG:    IF(TRIM(forback) == 'backward') THEN
        IF (switch) THEN
           uu(:,:,:,:) = -uu(:,:,:,:)
           vv(:,:,:,:) = -vv(:,:,:,:)
           ww(:,:,:,:) = -ww(:,:,:,:)
        ENDIF

        DO n = 1, maxsegm

           IF (num(n) == 1) THEN

              IF ((ndir(n) == 1).OR.(ndir(n) == -1)) THEN
                 it0  = it1(n)
                 iu0  = it0
                 it00 = it0 + 1

                 IF (ndir(n) < 0) THEN
                    iu0  = it0 - 1
                    it00 = it0 - 1
                 ENDIF

                 DO k = kt1(n), kt2(n)
                    DO j = jt1(n), jt2(n)

                       CALL sub_reducmem_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

                          !!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
                          !!NG: bug 14/04/2009   IF(key_alltracers) THEN
                          !!NG: bug 14/04/2009      init_temp(nb_traj) = tt(it0s,js,ks,1)
                          !!NG: bug 14/04/2009      init_salt(nb_traj) = ss(it0s,js,ks,1)
                          !!NG: bug 14/04/2009      init_dens(nb_traj) = rr(it0s,js,ks,1)
                          !!NG: bug 14/04/2009   ENDIF ! key_alltracers
                          !!NG: bug 14/04/2009 ENDIF ! key_nointerpolstats

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

                          u = uu(iu0s, js, ks, 1)

                          !----------------------------------
                          ! test on the sign of the velocity
                          !----------------------------------
                          IF ((ABS(u) > trmin).AND.(REAL(ndir(n),kind=rprec) * u > 0._rprec)) THEN

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

                             init_trans=ABS(u)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)

                             DO j1 = 1, n1
                                DO k1 = 1, n1
                                   DO l1 = 1, n2

                                      nb_traj = nb_traj + 1

                                      IF (nb_traj > nmax) THEN
                                         WRITE(lun_error,*)'too many particles, maximum is: ',nmax
                                         STOP
                                      ENDIF

                                      tfi(nb_traj) = REAL(iu0,kind=rprec)

                                      tfj(nb_traj) = REAL(j,kind=rprec) - 1._rprec - &
                                           .5_rprec/REAL(n1,kind=rprec) + &
                                           REAL(j1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfk(nb_traj) = REAL(k,kind=rprec)            - &
                                           .5_rprec/REAL(n1,kind=rprec) + &
                                           REAL(k1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfl(nb_traj) = REAL(l,kind=rprec) - .5_rprec - &
                                           .5_rprec/REAL(n2,kind=rprec) + &
                                           REAL(l1,kind=rprec)/REAL(n2,kind=rprec)

                                      ttr(nb_traj) = init_trans

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

                                         trans_total = trans_total + init_trans

                                         isn(-1) = isn(-1) + 1
                                         sn(-1)  =  sn(-1) + init_trans

                                         IF (key_alltracers) THEN

                                            IF (.NOT.key_nointerpolstats) THEN
                                               init_temp(nb_traj) = zinter(tt, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               init_salt(nb_traj) = zinter(ss, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               IF (key_approximatesigma) THEN
                                                  init_dens(nb_traj) = zinter(rr, & 
                                                       tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                       1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               ELSE ! key_approximatesigma
                                                  init_dens(nb_traj) = sigma(zsigma, init_salt(nb_traj), init_temp(nb_traj))

                                               ENDIF ! key_approximatesigma

                                            ELSE                                     !!NG: bug correction 14/04/2009
                          
                                               init_temp(nb_traj) = tt(it0s,js,ks,1) !!NG: bug correction 14/04/2009
                                               init_salt(nb_traj) = ss(it0s,js,ks,1) !!NG: bug correction 14/04/2009
                                               init_dens(nb_traj) = rr(it0s,js,ks,1) !!NG: bug correction 14/04/2009

                                            ENDIF ! key_nointerpolstats

                                            IF (init_temp(nb_traj) <= s0min(5,-1)) s0min(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) <= s0min(6,-1)) s0min(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) <= s0min(7,-1)) s0min(7,-1) = init_dens(nb_traj)
                                            IF (init_temp(nb_traj) >= s0max(5,-1)) s0max(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) >= s0max(6,-1)) s0max(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) >= s0max(7,-1)) s0max(7,-1) = init_dens(nb_traj)

                                            s0x(5,-1)  = s0x(5,-1)  + init_temp(nb_traj)        * init_trans
                                            s0x(6,-1)  = s0x(6,-1)  + init_salt(nb_traj)        * init_trans
                                            s0x(7,-1)  = s0x(7,-1)  + init_dens(nb_traj)        * init_trans
                                            s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
                                            s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
                                            s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans

                                         ENDIF !  key_alltracers

                                         IF (key_roms.OR.key_symphonie) THEN
                                            prof = &
                                                 fz(gi=tfi(nb_traj), gj=tfj(nb_traj), gk=-tfk(nb_traj), il=1)
                                         ELSE
                                            prof = fz(gk=-tfk(nb_traj))
                                         ENDIF

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
                                              s0min(1,-1) = fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
                                              s0min(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof <= s0min(3,-1)) s0min(3,-1) = -prof

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
                                              s0max(1,-1) = fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
                                              s0max(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof >= s0max(3,-1)) s0max(3,-1) = -prof

                                         s0x(1,-1) = s0x(1,-1) + fx(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x(2,-1) = s0x(2,-1) + fy(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x(3,-1) = s0x(3,-1) - prof * init_trans

                                         s0x2(1,-1) = s0x2(1,-1) + &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(2,-1) = s0x2(2,-1) + &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(3,-1) = s0x2(3,-1) + prof * prof * init_trans

                                      ELSE

                                         nb_traj = nb_traj - 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 j,k,l
                    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 k = kt1(n), kt2(n)

                    IF (key_periodic) THEN
                       born1 = MAX(2,it1(n))
                       born2 = MIN(imt-1,it2(n))

                    ELSE ! key_periodic
                       born1 = it1(n)
                       born2 = it2(n)

                    ENDIF ! key_periodic

                    DO i = born1, born2

                       CALL sub_reducmem_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

                          !!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
                          !!NG: bug 14/04/2009   IF (key_alltracers) THEN
                          !!NG: bug 14/04/2009      init_temp(nb_traj) = tt(is,jt0s,ks,1)
                          !!NG: bug 14/04/2009      init_salt(nb_traj) = ss(is,jt0s,ks,1)
                          !!NG: bug 14/04/2009      init_dens(nb_traj) = rr(is,jt0s,ks,1)
                          !!NG: bug 14/04/2009   ENDIF
                          !!NG: bug 14/04/2009 ENDIF

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

                          v = vv(is,jv0s,ks,1)

                          !------------------------------------!
                          !- test on the sign of the velocity -!
                          !------------------------------------!
                          IF ((ABS(v) > trmin).AND.(REAL(ndir(n),kind=rprec) * v > 0._rprec)) THEN
                             !
                             IF (lmt == 1) THEN
                                n1 = 1 + INT( SQRT(ABS(v)/max_transport) )
                                n2 = 1
                             ELSE
                                n1=  1 + INT( (ABS(v)/max_transport)**(1._rprec/3._rprec) )
                                n2= n1
                             ENDIF

                             init_trans = ABS(v)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)

                             DO i1 = 1, n1
                                DO k1 = 1, n1
                                   DO l1 = 1, n2

                                      nb_traj = nb_traj + 1

                                      IF (nb_traj > nmax) THEN
                                         WRITE(lun_error,*)'too many particles, maximum is: ',nmax
                                         STOP
                                      ENDIF

                                      tfi(nb_traj) = REAL(i,kind=rprec) -1._rprec - &
                                           .5_rprec/REAL(n1,kind=rprec)           + &
                                           REAL(i1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfj(nb_traj) = REAL(jv0,kind=rprec)

                                      tfk(nb_traj) = REAL(k,kind=rprec)           - &
                                           .5_rprec/REAL(n1,kind=rprec)           + &
                                           REAL(k1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfl(nb_traj) = REAL(l,kind=rprec) -.5_rprec - &
                                           .5_rprec/REAL(n2,kind=rprec)           + &
                                           REAL(l1,kind=rprec)/REAL(n2,kind=rprec)

                                      ttr(nb_traj) = init_trans

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

                                         trans_total = trans_total + init_trans

                                         isn(-1) = isn(-1) + 1
                                         sn(-1)  =  sn(-1) + init_trans

                                         IF (key_alltracers) THEN

                                            IF (.NOT.key_nointerpolstats) THEN
                                               init_temp(nb_traj) = zinter(tt, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               init_salt(nb_traj) = zinter(ss, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               IF (key_approximatesigma) THEN
                                                  init_dens(nb_traj) = zinter(rr, &
                                                       tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                       1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               ELSE
                                                  init_dens(nb_traj) = sigma(zsigma, init_salt(nb_traj), init_temp(nb_traj))

                                               ENDIF ! key_approximatesigma

                                            ELSE                                     !!NG: bug correction 14/04/2009

                                               init_temp(nb_traj) = tt(is,jt0s,ks,1) !!NG: bug correction 14/04/2009
                                               init_salt(nb_traj) = ss(is,jt0s,ks,1) !!NG: bug correction 14/04/2009
                                               init_dens(nb_traj) = rr(is,jt0s,ks,1) !!NG: bug correction 14/04/2009

                                            ENDIF !  key_nointerpolstats

                                            IF (init_temp(nb_traj) <= s0min(5,-1)) s0min(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) <= s0min(6,-1)) s0min(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) <= s0min(7,-1)) s0min(7,-1) = init_dens(nb_traj)
                                            IF (init_temp(nb_traj) >= s0max(5,-1)) s0max(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) >= s0max(6,-1)) s0max(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) >= s0max(7,-1)) s0max(7,-1) = init_dens(nb_traj)

                                            s0x(5,-1)  = s0x(5,-1)  + init_temp(nb_traj)        * init_trans
                                            s0x(6,-1)  = s0x(6,-1)  + init_salt(nb_traj)        * init_trans
                                            s0x(7,-1)  = s0x(7,-1)  + init_dens(nb_traj)        * init_trans
                                            s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
                                            s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
                                            s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans

                                         ENDIF ! key_alltracers

                                         IF (key_roms.OR.key_symphonie) THEN

                                            prof = &
                                                 fz(gi=tfi(nb_traj), gj=tfj(nb_traj), gk=-tfk(nb_traj), il=1)

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

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
                                              s0min(1,-1) = fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
                                              s0min(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof <= s0min(3,-1)) s0min(3,-1) = -prof

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
                                              s0max(1,-1)=fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
                                              s0max(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof >= s0max(3,-1)) s0max(3,-1) = -prof

                                         s0x(1,-1) = s0x(1,-1) + fx(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x(2,-1) = s0x(2,-1) + fy(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x(3,-1) = s0x(3,-1) - prof * init_trans

                                         s0x2(1,-1) = s0x2(1,-1) + &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(2,-1) = s0x2(2,-1) + &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(3,-1) = s0x2(3,-1) + prof * prof * init_trans

                                      ELSE

                                         nb_traj = nb_traj - 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
                 !
                 ! 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 j = jt1(n), jt2(n)
                    DO i = it1(n), it2(n)

                       CALL sub_reducmem_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

                          !!NG: bug 14/04/2009 IF (key_nointerpolstats) THEN
                          !!NG: bug 14/04/2009    IF (key_alltracers) THEN
                          !!NG: bug 14/04/2009       init_temp(nb_traj) = tt(is,js,kt0s,1)
                          !!NG: bug 14/04/2009       init_salt(nb_traj) = ss(is,js,kt0s,1)
                          !!NG: bug 14/04/2009       init_dens(nb_traj) = rr(is,js,kt0s,1)
                          !!NG: bug 14/04/2009    ENDIF
                          !!NG: bug 14/04/2009 ENDIF

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

                          w=ww(is,js,kw0s,1)

                          !------------------------------------!
                          !- test on the sign of the velocity -!
                          !------------------------------------!
                          IF ((ABS(w) > trmin).AND.(REAL(ndir(n),kind=rprec) * w < 0._rprec)) THEN

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

                             init_trans = ABS(w)/REAL(n1,kind=rprec)/REAL(n1,kind=rprec)/REAL(n2,kind=rprec)

                             DO j1 = 1, n1
                                DO i1 = 1, n1
                                   DO l1 = 1, n2

                                      nb_traj = nb_traj + 1

                                      IF (nb_traj > nmax) THEN
                                         WRITE(lun_error,*)'too many particles, maximum is: ', nmax
                                         STOP
                                      ENDIF

                                      tfi(nb_traj) = REAL(i,kind=rprec)   -1._rprec - &
                                           .5_rprec/REAL(n1,kind=rprec)             + &
                                           REAL(i1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfj(nb_traj) = REAL(j,kind=rprec)   -1._rprec - &
                                           .5_rprec/REAL(n1,kind=rprec)             + &
                                           REAL(j1,kind=rprec)/REAL(n1,kind=rprec)

                                      tfk(nb_traj) = REAL(kw0,kind=rprec)

                                      tfl(nb_traj) = REAL(l,kind=rprec)   -.5_rprec - &
                                           .5_rprec/REAL(n2,kind=rprec)             + &
                                           REAL(l1,kind=rprec)/REAL(n2,kind=rprec)

                                      ttr(nb_traj) = init_trans
                                      !
                                      IF (criter0(tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                           tfl(nb_traj),i,j,kt0,l)) THEN
                                         !
                                         trans_total = trans_total + init_trans

                                         isn(-1) = isn(-1) + 1
                                         sn(-1)  =  sn(-1) + init_trans

                                         IF (key_alltracers) THEN

                                            IF (.NOT.key_nointerpolstats) THEN
                                               init_temp(nb_traj) = zinter(tt, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )
                                               init_salt(nb_traj) = zinter(ss, &
                                                    tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                    1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               IF (key_approximatesigma) THEN
                                                  init_dens(nb_traj) = zinter(rr, &
                                                       tfi(nb_traj),tfj(nb_traj),-tfk(nb_traj), &
                                                       1._rprec + (tfl(nb_traj) - INT(tfl(nb_traj))) )

                                               ELSE
                                                  init_dens(nb_traj) = sigma(zsigma, init_salt(nb_traj), init_temp(nb_traj))

                                               ENDIF

                                            ELSE                                     !!NG: bug correction 14/04/2009

                                               init_temp(nb_traj) = tt(is,js,kt0s,1) !!NG: bug correction 14/04/2009
                                               init_salt(nb_traj) = ss(is,js,kt0s,1) !!NG: bug correction 14/04/2009
                                               init_dens(nb_traj) = rr(is,js,kt0s,1) !!NG: bug correction 14/04/2009

                                            ENDIF

                                            IF (init_temp(nb_traj) <= s0min(5,-1)) s0min(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) <= s0min(6,-1)) s0min(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) <= s0min(7,-1)) s0min(7,-1) = init_dens(nb_traj)
                                            IF (init_temp(nb_traj) >= s0max(5,-1)) s0max(5,-1) = init_temp(nb_traj)
                                            IF (init_salt(nb_traj) >= s0max(6,-1)) s0max(6,-1) = init_salt(nb_traj)
                                            IF (init_dens(nb_traj) >= s0max(7,-1)) s0max(7,-1) = init_dens(nb_traj)

                                            s0x(5,-1)  = s0x(5,-1)  + init_temp(nb_traj)        * init_trans
                                            s0x(6,-1)  = s0x(6,-1)  + init_salt(nb_traj)        * init_trans
                                            s0x(7,-1)  = s0x(7,-1)  + init_dens(nb_traj)        * init_trans
                                            s0x2(5,-1) = s0x2(5,-1) + init_temp(nb_traj) * init_temp(nb_traj) * init_trans
                                            s0x2(6,-1) = s0x2(6,-1) + init_salt(nb_traj) * init_salt(nb_traj) * init_trans
                                            s0x2(7,-1) = s0x2(7,-1) + init_dens(nb_traj) * init_dens(nb_traj) * init_trans


                                         ENDIF

                                         IF (key_roms.OR.key_symphonie) THEN
                                            prof = &
                                                 fz(gi=tfi(nb_traj), gj=tfj(nb_traj), gk=-tfk(nb_traj), il=1)
                                         ELSE
                                            prof = fz(gk=-tfk(nb_traj))
                                         ENDIF

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) <= s0min(1,-1)) &
                                              s0min(1,-1) = fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) <= s0min(2,-1)) &
                                              s0min(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof <= s0min(3,-1)) s0min(3,-1) = -prof

                                         IF (fx(tfi(nb_traj),tfj(nb_traj)) >= s0max(1,-1)) &
                                              s0max(1,-1) = fx(tfi(nb_traj),tfj(nb_traj))

                                         IF (fy(tfi(nb_traj),tfj(nb_traj)) >= s0max(2,-1)) &
                                              s0max(2,-1) = fy(tfi(nb_traj),tfj(nb_traj))

                                         IF (-prof >= s0max(3,-1)) s0max(3,-1) = -prof

                                         s0x(1,-1)  = s0x(1,-1) + fx(tfi(nb_traj),tfj(nb_traj))*init_trans

                                         s0x(2,-1)  = s0x(2,-1) + fy(tfi(nb_traj),tfj(nb_traj))*init_trans

                                         s0x(3,-1)  = s0x(3,-1) - prof * init_trans

                                         s0x2(1,-1) = s0x2(1,-1) + &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * &
                                              fx(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(2,-1) = s0x2(2,-1) + &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * &
                                              fy(tfi(nb_traj),tfj(nb_traj)) * init_trans

                                         s0x2(3,-1) = s0x2(3,-1) + prof * prof * init_trans

                                      ELSE

                                         nb_traj = nb_traj - 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
                 !
                 ! fin test kt1=kt2
              ENDIF
              !
              ! fin test num(n)
           ENDIF
           !
           ! fin du positionnement
        END DO
     END DO

     IF (switch) THEN
        forback='backward'
     ENDIF

  ELSE

     WRITE(lun_standard,*)'     - in Posini just to initialize some arrays...'
     WRITE(lun_standard,*)'       bin = ', TRIM(bin)

  ENDIF

  IF (ALLOCATED(num)) THEN
     CALL sub_memory(-size(num),'i','num','posini_seq')
     DEALLOCATE(num)
  END IF
  IF (ALLOCATED(it1)) THEN
     CALL sub_memory(-size(it1),'i','it1','posini_seq')
     DEALLOCATE(it1)
  END IF
  IF (ALLOCATED(it2)) THEN
     CALL sub_memory(-size(it2),'i','it2','posini_seq')
     DEALLOCATE(it2)
  END IF
  IF (ALLOCATED(jt1)) THEN
     CALL sub_memory(-size(jt1),'i','jt1','posini_seq')
     DEALLOCATE(jt1)
  END IF
  IF (ALLOCATED(jt2)) THEN
     CALL sub_memory(-size(jt2),'i','jt2','posini_seq')
     DEALLOCATE(jt2)
  END IF
  IF (ALLOCATED(kt1)) THEN
     CALL sub_memory(-size(kt1),'i','kt1','posini_seq')
     DEALLOCATE(kt1)
  END IF
  IF (ALLOCATED(kt2)) THEN
     CALL sub_memory(-size(kt2),'i','kt2','posini_seq')
     DEALLOCATE(kt2)
  END IF
  IF ( ALLOCATED(ndir)) THEN
     CALL sub_memory(-size(ndir),'i','ndir','posini_seq')
     DEALLOCATE(ndir)
  END IF

  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'  ----------------------------'
  WRITE(lun_standard,*)'  =  POSINI (SEQ)  >>>> EXIT ='
  WRITE(lun_standard,*)'  ----------------------------'

  RETURN

END SUBROUTINE posini_seq