!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! - 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/trajec
!! NAME
!!   trajec (trajec.f90 - Main program - Fortran 90)
!!
!! USAGE
!!   
!!
!! FUNCTION
!! --- Particle trajectories from the analytical computation of 3D
!!     streamlines in a 3D velocity fiels.
!!
!! --- Qualitative or quantitative diagnostics
!!
!! --- Computations done over periods equal to the available sampling
!!     for the velocity field, and over which the velocity is assumed
!!     constant
!!
!! --- The program works with the OPA tensorial formalism
!!
!! --- input: simplified meshmask, 3D TRANSPORT field (velocity x
!!     section), parameters (on various files)
!!
!! --- output: individual trajectories or quantitative outputs
!!
!! --- no-vectorized version
!!
!! --- some comments (yes !) in English
!!
!! --- copyright: Bruno Blanke (Bruno.Blanke@univ-brest.fr)
!! first version: Summer 1992
!! this version: March 2007
!! http://www.ifremer.fr/lpo/blanke/ARIANE
!! 
!! 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
!!
!! (see mod_lun)
!! current input files:
!! -------------------

!! current output files:
!! ---------------------
!!
!! TODO
!!   
!!
!! PORTABILITY
!!         Machine-OS    - Fortran90/95 compiler
!!   * i686-pc-linux-gnu -         ifort
!!
!! SEE ALSO
!!
!! SOURCE
!!=========================================================================

SUBROUTINE trajec_seq() 1,342

  USE mod_precision
  USE mod_cst
  USE mod_namelist
  USE mod_input_grid
  USE mod_seq
  USE mod_input_data
  USE mod_output_data
  USE mod_trajec_subs
  USE mod_init_particules
  USE mod_posin
  USE mod_flags
  USE mod_stats
  USE mod_stati
  USE mod_txt
  USE mod_lun
  USE mod_fx
  USE mod_fy
  USE mod_fz
  USE mod_zinter
  USE mod_sigma
  USE mod_quant
  USE mod_reducmem
  USE mod_orca
  USE mod_save_netcdf

  !--------------!
  ! DECLARATIONS !
  !--------------!
  IMPLICIT NONE 

  INTEGER(kind=iprec) :: &
       nsect2      

  INTEGER(kind=iprec) :: &
       i, j, n         , & !
       is, js, ks      , & !
       nis, njs, nks   , & !
       i1s, j1s, k1s   , & !
       i2s, j2s, k2s   , & !
       ndone           , & !
       iter            , & !
       lref     = 1    , & !
       kk0             , & !
       ifl             , & !
       !!       cnt      = 1    , & ! counter
       iter_min        , & !
       iter_max        , & ! 
       imin_tmp        , & !
       imax_tmp            !

  INTEGER(kind=iprec), DIMENSION(:), ALLOCATABLE ::  &
       igo   , & !
       i0    , & !
       j0    , & !
       k0    , & !
       l0    , & !
       i1    , & !
       i2    , & !
       j1    , & !
       j2    , & !
       k1    , & !
       k2    , & !
       iflag , & !
       iperfl, & !
       jperfl, & !
       ind_traj, & !
       nfin

  REAL(kind=rprec) :: &
       subtime =   -1._rprec  , & !
       tfic    =    0._rprec  , & !
       tlap    =    0._rprec  , & !
       tfin    =    0._rprec  , & !
       trtot1  =    0._rprec  , & !
       prof    =    0._rprec  , & !
       x0   , & !
       y0   , & !
       z0   , & !
       ft   , & !
       xt   , & !
       yt   , & !
       zt   , & !
       st   , & !
       r_lmt, & !
       dumm

  REAL(kind=rprec), DIMENSION(:), ALLOCATABLE :: & !
       age  , & !
       agenew, & !
       truage, & !
       gi   , & !
       gj   , & !
       gk   , & !
       gl   , & !
       hi   , & !
       hj   , & !
       hk   , & !
       hl   , & !
       tnext, & !
       tswap, & !
       zbuoy, & !
       ti   , & !
       si   , & !
       ri   , & !
       u    , & !
       v    , & !
       w    , & !
       du   , & !
       dv   , & !
       dw   , & !
       tx   , & !
       ty   , & !
       tz   , & !
       tf   , & !
       sf   , & !
       rf   , & !
       tfac , & !
       init_depth , & ! !! NG: 15_09_2008
       final_depth !! NG: 15_09_2008

  REAL(kind = qprec), DIMENSION(:), ALLOCATABLE :: ttt

  REAL(kind=rprec), DIMENSION(:,:,:), ALLOCATABLE :: & !
       array_qual_traj

  LOGICAL :: open_door=.TRUE., first_pass=.TRUE., first_call=.TRUE.

  !-------------------------------------------!
  !- DYNAMIC ALLOCATIONS AND INITIALIZATIONS -!
  !-------------------------------------------!
  CALL sub_posin_alloc(nmax)

  ALLOCATE(igo(nmax))   ; igo(:)   = 0
  CALL sub_memory(size(igo),'i','igo','trajec_seq')
  ALLOCATE(i0(nmax))    ; i0(:)    = 0
  CALL sub_memory(size(i0),'i','i0','trajec_seq')
  ALLOCATE(j0(nmax))    ; j0(:)    = 0 
  CALL sub_memory(size(j0),'i','j0','trajec_seq')
  ALLOCATE(k0(nmax))    ; k0(:)    = 0
  CALL sub_memory(size(k0),'i','k0','trajec_seq')
  ALLOCATE(l0(nmax))    ; l0(:)    = 0
  CALL sub_memory(size(l0),'i','l0','trajec_seq')
  ALLOCATE(i1(nmax))    ; i1(:)    = 0
  CALL sub_memory(size(i1),'i','i1','trajec_seq')
  ALLOCATE(i2(nmax))    ; i2(:)    = 0
  CALL sub_memory(size(i2),'i','i2','trajec_seq')
  ALLOCATE(j1(nmax))    ; j1(:)    = 0
  CALL sub_memory(size(j1),'i','j1','trajec_seq')
  ALLOCATE(j2(nmax))    ; j2(:)    = 0
  CALL sub_memory(size(j2),'i','j2','trajec_seq')
  ALLOCATE(k1(nmax))    ; k1(:)    = 0
  CALL sub_memory(size(k1),'i','k1','trajec_seq')
  ALLOCATE(k2(nmax))    ; k2(:)    = 0
  CALL sub_memory(size(k2),'i','k2','trajec_seq')
  ALLOCATE(iflag(nmax)) ; iflag(:) = 0
  CALL sub_memory(size(iflag),'i','iflag','trajec_seq')
  ALLOCATE(iperfl(nmax)); iperfl(:)= 0
  CALL sub_memory(size(iperfl),'i','iperfl','trajec_seq')
  ALLOCATE(jperfl(nmax)); jperfl(:)= 0
  CALL sub_memory(size(jperfl),'i','jperfl','trajec_seq')
  ALLOCATE(nfin(nmax))  ; nfin(:)  = -1
  CALL sub_memory(size(nfin),'i','nfin','trajec_seq')

  ALLOCATE(age(nmax))   ; age(:)   = 0._rprec
  CALL sub_memory(size(age),'r','age','trajec_seq')
  ALLOCATE(agenew(nmax)); agenew(:)= 0._rprec
  CALL sub_memory(size(agenew),'r','agenew','trajec_seq')
  ALLOCATE(truage(nmax)); truage(:)= 0._rprec
  CALL sub_memory(size(truage),'r','truage','trajec_seq')

  tage(:) = 0._rprec ! define is mod_posin.f90

  ALLOCATE(gi(nmax))   ; gi(:)   = 0._rprec
  CALL sub_memory(size(gi),'r','gi','trajec_seq')
  ALLOCATE(gj(nmax))   ; gj(:)   = 0._rprec
  CALL sub_memory(size(gj),'r','gj','trajec_seq')
  ALLOCATE(gk(nmax))   ; gk(:)   = 0._rprec
  CALL sub_memory(size(gk),'r','gk','trajec_seq')
  ALLOCATE(gl(nmax))   ; gl(:)   = 0._rprec
  CALL sub_memory(size(gl),'r','gl','trajec_seq')
  ALLOCATE(hi(nmax))   ; hi(:)   = 0._rprec
  CALL sub_memory(size(hi),'r','hi','trajec_seq')
  ALLOCATE(hj(nmax))   ; hj(:)   = 0._rprec
  CALL sub_memory(size(hj),'r','hj','trajec_seq')
  ALLOCATE(hk(nmax))   ; hk(:)   = 0._rprec
  CALL sub_memory(size(hk),'r','hk','trajec_seq')
  ALLOCATE(hl(nmax))   ; hl(:)   = 0._rprec
  CALL sub_memory(size(hl),'r','hl','trajec_seq')
  ALLOCATE(tnext(nmax)); tnext(:)= 0._rprec
  CALL sub_memory(size(tnext),'r','tnext','trajec_seq')
  ALLOCATE(tswap(nmax)); tswap(:)= 0._rprec
  CALL sub_memory(size(tswap),'r','tswap','trajec_seq')
  ALLOCATE(zbuoy(nmax)); zbuoy(:)= 1.e9_rprec
  CALL sub_memory(size(zbuoy),'r','zbuoy','trajec_seq')
  ALLOCATE(u(nmax))    ; u(:)    = 0._rprec
  CALL sub_memory(size(u),'r','u','trajec_seq')
  ALLOCATE(v(nmax))    ; v(:)    = 0._rprec
  CALL sub_memory(size(v),'r','v','trajec_seq')
  ALLOCATE(w(nmax))    ; w(:)    = 0._rprec
  CALL sub_memory(size(w),'r','w','trajec_seq')
  ALLOCATE(du(nmax))   ; du(:)   = 0._rprec
  CALL sub_memory(size(du),'r','du','trajec_seq')
  ALLOCATE(dv(nmax))   ; dv(:)   = 0._rprec
  CALL sub_memory(size(dv),'r','dv','trajec_seq')
  ALLOCATE(dw(nmax))   ; dw(:)   = 0._rprec
  CALL sub_memory(size(dw),'r','dw','trajec_seq')
  ALLOCATE(tx(nmax))   ; tx(:)   = 0._rprec
  CALL sub_memory(size(tx),'r','tx','trajec_seq')
  ALLOCATE(ty(nmax))   ; ty(:)   = 0._rprec
  CALL sub_memory(size(ty),'r','ty','trajec_seq')
  ALLOCATE(tz(nmax))   ; tz(:)   = 0._rprec
  CALL sub_memory(size(tz),'r','tz','trajec_seq')
  ALLOCATE(tf(nmax))   ; tf(:)   = 0._rprec
  CALL sub_memory(size(tf),'r','tf','trajec_seq')
  ALLOCATE(sf(nmax))   ; sf(:)   = 0._rprec
  CALL sub_memory(size(sf),'r','sf','trajec_seq')
  ALLOCATE(rf(nmax))   ; rf(:)   = 0._rprec
  CALL sub_memory(size(rf),'r','rf','trajec_seq')
  ALLOCATE(tfac(nmax)) ; tfac(:) = 0._rprec
  CALL sub_memory(size(tfac),'r','tfac','trajec_seq')
  ALLOCATE(ttt(nmax))  ; ttt(:)  = 0._qprec
  CALL sub_memory(size(ttt),'r','ttt','trajec_seq')

  !! NG: 15_09_2008
  ALLOCATE(init_depth(nmax))  ; init_depth(:)  = 0._rprec
  CALL sub_memory(size(init_depth),'r','init_depth','trajec_seq')
  ALLOCATE(final_depth(nmax)) ; final_depth(:) = 0._rprec
  CALL sub_memory(size(final_depth),'r','final_depth','trajec_seq')

  IF (key_alltracers) THEN
     ALLOCATE(ti(nmax)); ti(:) = 0._rprec
     CALL sub_memory(size(ti),'r','ti','trajec_seq')
     ALLOCATE(si(nmax)); si(:) = 0._rprec
     CALL sub_memory(size(si),'r','si','trajec_seq')
     ALLOCATE(ri(nmax)); ri(:) = 0._rprec
     CALL sub_memory(size(ri),'r','ri','trajec_seq')
  ENDIF

  !======================!
  !- Read region limits -!
  !======================!----------------------------!
  ! This is done in qualitative and quantitative mode !
  !---------------------------------------------------!
  CALL sub_reducmem_read_reg_lim()

  IF (key_ascii_outputs) THEN
     !==========================!
     ! Open ASCII output files -!
     !==========================!
     CALL sub_output_data_open_files()
  ENDIF

  !-------------------!
  !- READ INPUT MESH -!
  !-------------------!---------------------------------------------------!
  !-- Allocate and read coordinates and scale factors 
  !-- xx_tt, xx_uu, xx_vv, xx_ff
  !-- yy_tt, yy_uu, yy_vv, yy_ff
  !-- zz_tt, zz_ww
  !-- e2u
  !-- e1v
  !-- e1t, e2t, e3t
  !-- tmask
  !-----------------------------------------------------------------------!
  CALL sub_input_grid()

  !==========================
  !- QUALITATIVE VARIABLES -!
  !=================================================================
  ! sampling time (in seconds) for the available transport field
  !=================================================================
  tfic = tunit   * REAL(ntfic, kind=rprec)

  !=================================================================
  ! sampling time (in seconds) between two successive output positions
  !=================================================================
  tlap = ABS(delta_t) * REAL(frequency, kind=rprec)

  !=================================================================
  ! maximum integration time (im seconds) for individual trajectories
  !=================================================================
  tfin = tlap    * REAL(ABS(nb_output), kind=rprec)

  !==============================================
  ! mask is written (for graphical trajectories)
  !==============================================
  IF ((TRIM(mode) == 'qualitative').AND.(mask)) THEN
     DO j = 1, jmt
        DO i = 1, imt

           IF (tmask(i,j,1,1) == 0._rprec) THEN

              IF (key_alltracers) THEN
                 WRITE(lun_traj,7055)0,xx_tt(i,j,1,1),yy_tt(i,j,1,1),0.,0.,0.,0.,0.
              ELSE
                 WRITE(lun_traj,7055)0,xx_tt(i,j,1,1),yy_tt(i,j,1,1),0.,0.
              ENDIF

           ENDIF

        END DO
     END DO
  ENDIF

  !!NG: BEGIN SEQUENTIAL DATA ACCESS
  !=============================================!
  !- Allocate and Initialize input data arrays -!
  !=============================================!
  CALL sub_input_data_alloc_init_seq()
  !!NG: END SEQUENTIAL DATA ACCESS

  !=================================================================
  ! we initialize the *new* particle position 
  !=================================================================
  ! BBL __________________________________________________________________________
  ! BBL BBL: Initialisation a zero des prochaines positions a calculer sur la
  ! BBL        grille spatiale (I, J, K) et temporelle (L)
  ! BBL
  hi(:) = 0._rprec
  hj(:) = 0._rprec
  hk(:) = 0._rprec
  hl(:) = 0._rprec

  !=======================================!
  !- Read or Compute particule positions -!
  !=======================================!
  CALL sub_init_particules_positions()

  !=======================!
  !- NetCDF data outputs -!
  !=======================!
  !---------------------------!
  !- Open Netcdf output file -!
  !---------------------------!
  CALL sub_save_netcdf_data_init()

  !!NG  !------------------------------------------------!
  !!NG  !- Save Initial Positions in netcdf output file -!
  !!NG  !------------------------------------------------!
  !!NG  CALL sub_save_netcdf_init_pos()
  !!NG
  !!NG  !- Deallocate init_temp, init_salt, init_dens (mod_posin) -!
  !!NG  CALL sub_posin_dealloc_init_tracers()

  IF (TRIM(mode) == 'quantitative') THEN
     CALL sub_quant_alloc(imt, jmt, kmt)
     CALL sub_quant_init() !  initialize a few arrays

     !-------------------------------------------------------------------------
     ! we do not yet know whether a final criterion (criter1) will be satisfied
     !-------------------------------------------------------------------------
     icrit1 = 0

     secname(nsect+1) = 'Criter1'

  ELSEIF (TRIM(mode) == 'qualitative') THEN

     ALLOCATE(ind_traj(ntraj))
     CALL sub_memory(size(ind_traj),'i','ind_traj','trajec_seq')
     ALLOCATE(array_qual_traj(ntraj,nb_output+1,nstat))
     CALL sub_memory(size(array_qual_traj),'r','array_qual_traj','trajec_seq')
     array_qual_traj(:,:,:)=mask_value

  ELSE
     STOP

  ENDIF

  !=============================================================!
  != First loop on the number of trajectories; initializations =!
  !=============================================================!
  DO n = 1, ntraj

     !====================================================================
     ! time index must be within [0.5, lmt+.5[; otherwise time translation
     !NG time index must be within ]0.5, lmt+.5]; otherwise time translation
     !====================================================================
     !NG !!! To work "lmt" must be greater than zero !!! (else infinit loop)
     !NG !!! "lmt" must be tested before to use it.  !!!
     !NG !!! If "lmt" is tested, GOTO could be removed. !!!
     ! BBL __________________________________________________________________________
     ! BBL BBL: Recadrage du parametre temporel FL dans la fenetre [0.5, lmt+0.5[
     ! BBL
100  IF (tfl(n) < 0.5_rprec) THEN
        tfl(n) = tfl(n) + REAL(lmt,kind=rprec)
        GOTO 100
     ENDIF ! (fl(n) < 0.5)

     IF (tfl(n) >= (REAL(lmt,kind=rprec) + 0.5_rprec)) THEN
        tfl(n) = tfl(n) - REAL(lmt,kind=rprec)
        GOTO 100
     ENDIF

     !====================================================================
     ! the program is able to deal with constant-depth particles
     !  (qualitative experiments only) [QUALI]
     !====================================================================
     ! BBL __________________________________________________________________________
     ! BBL BBL: La variable ZBUOY, initialisee par defaut a 1.e9, n'est active que
     ! BBL        dans le cas de calculs en mode QUALITATIVE, et seulement alors pour
     ! BBL        les particules dont on veut calculer des trajectoires a profondeur
     ! BBL        CONSTANTE
     ! BBL       Sa valeur est alors egale a l'immersion initiale de chaque particule
     ! BBL        concernee.
     ! BBL       NOTE: dans le cas d'un maillage variable dans le temps (evolution
     ! BBL        d'une surface libre par exemple), le calcul de l'immersion de
     ! BBL        reference n'est pas parfait car le maillage vertical utilise dans
     ! BBL        la fonction FZ est un maillage statique (ZETA = 0) [la dynamique du
     ! BBL        modele, incluant ZETA, n'a en effet pas encore ete lue a ce stade
     ! BBL        du code]
     ! BBL
     zbuoy(n)=1.e9_rprec

     IF (tfk(n) <= 0._rprec) THEN

        IF (TRIM(mode) == 'quantitative') THEN 
           WRITE(lun_error,*)''
           WRITE(lun_error,*)'---    ERROR    ---'
           WRITE(lun_error,*)'Isobaric particles allowed in QUALITATIVE experiments ONLY '
           WRITE(lun_error,*)'Particule number n      = ', n
           WRITE(lun_error,*)'                 tfk(n) = ', tfk(n)
           WRITE(lun_error,*)'--- STOP ARIANE ---'
           STOP
        ENDIF

        IF (key_roms) THEN
           If (first_pass) THEN
              CALL sub_input_approx_zz_ww_roms()
              first_pass=.FALSE.
           ENDIF
           zbuoy(n)= fz(gi=tfi(n), gj=tfj(n), gk=tfk(n), il=1)
        ELSEIF (key_symphonie) THEN
           WRITE(lun_error,*)'       Isobaric particles'
           WRITE(lun_error,*)' NOT AVAILABLE FOR SYMPHONIE - STOP'
           STOP
        ELSE
           zbuoy(n)= fz(gk=tfk(n))
        ENDIF

        tfk(n)   = -tfk(n)

        !!NG: to debug___ write(*,*) zbuoy(n)

     ENDIF

     !====================================================================
     ! time (tfl(n)) and position (VELOCITY grid: fi, tfj(n), tfk(n))
     ! indices used to  detect the relevant embedding temperature
     !  gridcell (i, j, k, l).
     !====================================================================
     ! BBL __________________________________________________________________________
     ! BBL BBL: Definition de la maille TEMPERATURE (I0, J0, K0, L0) contenant la
     ! BBL        position initiale de chaque particule
     ! BBL
     i0(n) =  INT(tfi(n)) + 1
     j0(n) =  INT(tfj(n)) + 1
     k0(n) =  INT(tfk(n))
     l0(n) = NINT(tfl(n))


     !===========================!
     !- Periodicity in OPA-ORCA -!
     !===========================!
     IF (.NOT.(key_roms.OR.key_symphonie)) THEN

        CALL sub_orca_north_pole(i0(n), j0(n))

        CALL sub_orca_east_west_periodic(i0(n))

     ENDIF

     !===================================================================
     ! Compute diagnostics for initial tracers values
     !===================================================================
     ! BBL __________________________________________________________________________
     ! BBL BBL: La maille TEMPERATURE contenant la position initiale ne saurait
     ! BBL        etre une maille TERRE
     ! BBL
     CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is,js,ks)

     IF (tmask(is,js,ks,1) <= .5_rprec) THEN

        WRITE(lun_error,7000)'   false start', &
             n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)

        STOP
     ENDIF ! (tmask(i0(n),j0(n),k0(n),1)

     !------------------------------------------------------------------!
     !- fi, fj, fk and fl refer to the initial position                -!
     !- gi, gj, gk and gl refer to the current position                -!
     !- BEWARE: we use a reverse vertical axis for easier computations -!
     !------------------------------------------------------------------!
     ! BBL __________________________________________________________________________
     ! BBL BBL: Positionnement initial de chaque particule a partir des valeurs
     ! BBL        FI, FJ, FK et FL lues sur fichier ou calculees par POSINI
     ! BBL       NOTE: les calculs faits dans le coeur du code utilisent des indices
     ! BBL        verticaux NEGATIFS, d'ou le changement de signe opere sur gk
     gi(n) =  tfi(n)
     gj(n) =  tfj(n)
     gk(n) = -tfk(n)
     gl(n) =  tfl(n)

     !----------------------!
     !- age initialization -!
     !----------------------!
     ! BBL __________________________________________________________________________
     ! BBL BBL: La variable AGE definit l'age (en secondes) de chaque particule
     ! BBL        depuis le debut de son integration temporelle
     ! BBL
     age(n) = tage(n)

     ! BBL __________________________________________________________________________
     ! BBL BBL: La variable IGO est utilisee pour referencer les etats de chaque
     ! BBL        particule, et en particulier les conditions suivantes:
     ! BBL         - pas encore partie
     ! BBL         - en cours d'integration
     ! BBL         - en attente d'une valeur actualisee du champ de transport
     !             - arrivee a destination
     igo(n) = 0

     !
     ! next output instant (for positions in qualitative experiments)
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: En mode QUALITATIVE la variable TNEXT definit le prochain instant
     ! BBL        (en secondes) ou la sortie d'une position sur le fichier TRAJ.QL
     ! BBL        sera necessaire
     ! BBL       NOTE: ce temps sera utilise par comparaison a l'AGE de chaque
     ! BBL        particule; dans la mesure ou on echantillonne les trajectoires de
     ! BBL        toutes les particules avec la meme frequence, la variable TNEXT est
     ! BBL        initialisee a la meme valeur pour toutes les particules
     ! BBL
     tnext(n) = tlap

     !------------------------------------------------------------------
     !- next time swap instant (between two samples of velocity field) -
     !------------------------------------------------------------------
     ! BBL __________________________________________________________________________
     ! BBL BBL: La variable TSWAP definit le prochain instant (en secondes) ou une
     ! BBL        actualisation du champ de transport sera necessaire
     ! BBL       NOTE: ce temps sera utilise par comparaison avec l'AGE de chaque
     ! BBL        particule; dans la mesure ou les particules ne sont pas toutes
     ! BBL        initialisees au meme moment sur l'axe temporel, la valeur de TSWAP
     ! BBL        differe selon la particule consideree; son calcul depend aussi de
     ! BBL        l'orientation temporelle chosie (BACKWARD ou FORWARD)
     ! BBL
     IF (TRIM(forback) == 'forward') THEN
        tswap(n)=tfic*ABS( tfl(n)-(l0(n)+.5_rprec) ) + tage(n)
     ELSE
        tswap(n)=tfic*ABS( tfl(n)-(l0(n)-.5_rprec) ) + tage(n)
     ENDIF

     !!BBL - march 2007
     !!NG: option to follow an experience interrupted. Age is stored
     !!NG: in the final position file.
     !!NG: Dont USE this option in backward mode !!!
     IF (key_read_age) THEN
        dumm = tswap(n) / tfic
        tswap(n) = NINT(tswap(n)-NINT(dumm)*tfic)+NINT(dumm)*tfic
     ENDIF
     !
     ! end initializations
     !
  ENDDO ! n=1,ntraj

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

  !-------------------------------------------!
  !- boucle sur les pas de temps disponibles -!
  !-------------------------------------------!
  ! BBL __________________________________________________________________________
  ! BBL BBL: Pour des calculs SEQUENTIELS en mode CLIMATOLOGIQUE, il est utile de
  ! BBL       definir le nombre maximal de cycles autorises sur l'axe temporel, de
  ! BBL       maniere a pouvoir filtrer les particules *trop lentes*
  ! BBL
  !! TODO: maxcyl should be a namelist parameter...
  !!NG : namelist => maxcycl0 = 300
  !!NG : namelist => maxcycl  = maxcycl0
  ! BBL __________________________________________________________________________
  ! BBL BBL: Dans le cas LMT=1, tout le champ de transport est charge en memoire,
  ! BBL       et la restriction precedente peut etre relachee
  ! BBL
  !!NG : IF (lmt == 1) maxcycl=12599999
  ! BBL __________________________________________________________________________
  ! BBL BBL: Boucle sur la duree definie par l'archivage utilise pour le champ de
  ! BBL      vitesse
  ! BBL


  !!========================================================================!!
  !!============================== DO LOOP: CYCLES =========================!!
  !!========================================================================!!
  !- Comments -!
  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'=================================='
  WRITE(lun_standard,*)'= ENTER >>>>> MAIN LOOP ON CYCLE ='
  WRITE(lun_standard,*)'=================================='

  !! TODO: infini do loop and exit.
  DO iter = 1, maxcycles

     if (iter==1) THEN
        first_call=.TRUE.
     ELSE
        first_call=.FALSE.
     ENDIF

     !- Initialize some variables need to the sequential mode -!
     iter_min=1
     iter_max=lmt
     IF (iter == 1) THEN
        IF ( (lmin /= 1).AND.(TRIM(forback) == 'forward')) THEN
           CALL sub_seq_init(lmin)
           iter_min=lmin
        ELSEIF (TRIM(forback) == 'backward') THEN
           CALL sub_seq_init(lmax)
           iter_max=lmax
        ELSE
           CALL sub_seq_init()
        ENDIF
     ELSE
        CALL sub_seq_init()
     ENDIF

     ! BBL __________________________________________________________________________
     ! BBL BBL: On utlise NDONE pour comptabiliser le nombre de particules ayant
     ! BBL       atteint leur etat final; le calcul n'est fait que pour informer
     ! BBL       l'utilisateur sur les progres effectifs realises par les calculs,
     ! BBL       cycle temporel apres cycle temporel
     ! BBL      Le calcul est fait a chaque iteration ITER si LMT > 1, et seulement
     ! BBL       au debut de la premiere iteration si LMT=1 (dans ce dernier cas en
     ! BBL       effet tout le champ de transport peut etre charge en memoire; les
     ! BBL       calculs sont donc beaucoup plus rapides, et il n'est pas necessaire
     ! BBL       d'informer l'utilisateur...)
     ! BBL
     IF ( ((iter == 1).AND.(lmt == 1)).OR. &
          ((iter /= 1) .AND.(lmt > 1)) ) THEN

        IF (TRIM(mode) == 'quantitative') THEN
           ndone = 0 ! initialize to 0 the number of trajectorie(s) done

           DO  n = 1, ntraj

              ! BBL ___________________________________________________________________
              ! BBL -- BBL: IGO=3 ou -3 signifie que l'integration de la particule est
              ! BBL          terminee
              ! BBL
              IF (ABS(igo(n)) >= 3) ndone = ndone + 1

           ENDDO

           WRITE(lun_standard,*)
           WRITE(lun_standard,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
           WRITE(lun_standard,*)'HOW MANY PARTICLES REACH A SECTION: ', &
                INT(100._rprec*REAL(ndone,kind=rprec)/REAL(ntraj,kind=rprec)), '%'
           WRITE(lun_standard,*)'    - number of particules         : ', ntraj
           WRITE(lun_standard,*)'    - done                         : ', ndone
           WRITE(lun_standard,*)'    - rest                         : ', ntraj - ndone
           WRITE(lun_standard,*)'    - cycle                        : ', iter-1
           WRITE(lun_standard,*)'<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'

        ELSE

           WRITE(lun_standard,*)
           WRITE(lun_standard,*)'>>>>>>>>>>>>>>>>>>>>>>>>'
           WRITE(lun_standard,*)'>>  Cycle : ', iter-1
           WRITE(lun_standard,*)'<<<<<<<<<<<<<<<<<<<<<<<<'

        ENDIF

        id_comments = .FALSE.

     ELSE

        WRITE(lun_standard,*)
        WRITE(lun_standard,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
        WRITE(lun_standard,*)' PARTICULE TRAJECTORIES will be computed:'
        WRITE(lun_standard,*)'    - number of trajectories         : ', ntraj
        WRITE(lun_standard,*)'<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'

        !- To add comments when input data are reading -!
        !- comments only at the first reading -!
        id_comments = .TRUE.

     ENDIF ! ( ((iter == 1).AND.(lmt == 1)).OR.((iter /= 1) .AND.(lmt > 1)) )

     ! BBL __________________________________________________________________________
     ! BBL BBL: LREF designe l'indice actif (entre 1 et LMT) sur l'axe temporel,
     ! BBL        c'est-a-dire l'indice temporel du champ de transport a charger en
     ! BBL        memoire en debut de cycle
     ! BBL       LREF vaut donc 1 en mode FORWARD et LMT en mode BACKWARD
     ! BBL

     lref = iter_min

     IF (TRIM(forback) == 'backward') lref = iter_max

     !!    cnt = iter_min                                    ! counter
     !!    IF (TRIM(forback) == 'backward') cnt = iter_max ! Backward mode

     ! BBL __________________________________________________________________________
     ! BBL BBL: Le chargement effectif en memoire s'effectue obligatoirement a la
     ! BBL        premiere iteration (il faut faire au moins une lecture!), et a
     ! BBL        chaque modification de l'indice temporel si jamais LMT > 1
     ! BBL

9998 CONTINUE

     IF ( (iter == 1).OR.(lmt > 1) ) THEN

        !!--------------------------------!!
        !!- READ INPUT DATA SEQUENTIALLY -!!
        !!--------------------------------!!

        !- Comments -!
        WRITE(lun_standard,*)'---------------------------------'
        WRITE(lun_standard,*)'- READ INPUT DATA :', lref, '-'

        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(:,:)            , &
                first_call                   )

           !!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 ! OPA-NEMO
           CALL sub_input_data_seq_opa( &
                ncids(:)              , &
                varids(:)             , &
                new_file(:)           , &
                ind_file(:)           , &
                ind_time(:)           , &
                ind_time_size(:)      , &
                sdimsorders(:,:)       )
        ENDIF

        IF (TRIM(forback) == 'backward') THEN
           !!        cnt = cnt - 1
           delta_t = - delta_t
           uu(:,:,:,:) = -uu(:,:,:,:)
           vv(:,:,:,:) = -vv(:,:,:,:)
           ww(:,:,:,:) = -ww(:,:,:,:)
           !!      ELSE
           !!        cnt = cnt + 1
        ENDIF

        IF (open_door) THEN
           !------------------------------------------------!
           !- Save Initial Positions in netcdf output file -!
           !------------------------------------------------!
           CALL sub_save_netcdf_init_pos()

           !- Deallocate init_temp, init_salt, init_dens (mod_posin) -!
           CALL sub_posin_dealloc_init_tracers()

           open_door=.FALSE.
        ENDIF


     ENDIF ! ((iter == 1).OR.(lmt > 1))

     !!----------------------------------------------------------------------------!!
     !!---1111111111111111111---- DO LOOP : ON TRAJECTORIES ----111111111111111----!!
     !!----------------------------------------------------------------------------!!
     ! marquage des particules
     ! code utilise:
     !  igo=-3: integration de la particule terminee
     !  igo=-2: attente de la disponibilite d'un nouveau champ de transport
     !  igo=0: pas encore partie
     !  igo=1: active
     !  igo=2: interception temporelle (swap) venant d'etre realisee
     !  igo=3: interception spatiale/criter1 venant d'etre realisee
     !

     DO n = 1, ntraj
        ! BBL __________________________________________________________________________
        ! BBL BBL: Si l'integration de la particule s'est terminee lors de
        ! BBL        l'utilisation du champ de transport precedent (IGO=3), on la place
        ! BBL        definitivement dans l'etat IGO=-3
        ! BBL
        IF (ABS(igo(n)) == 3) igo(n) = -3

        ! BBL __________________________________________________________________________
        ! BBL BBL: Si l'integration de la particule avait ete stoppee du fait d'une
        ! BBL        mise a jour necesaire du champ de transport (IGO=+/-2), on la remet
        ! BBL        en activite (IGO=1)
        ! BBL
        IF (ABS(igo(n)) == 2) igo(n) = 1

        IF ((igo(n) == 0).AND.(l0(n) == lref)) THEN

           ! BBL __________________________________________________________________________
           ! BBL BBL: Si la particule n'est pas encore partie (IGO=0) et si l'instant
           ! BBL        LREF qui vient d'etre charge en en memoire l'autorisee, elle est
           ! BBL        mise en activite (IGO=1)
           ! BBL
           igo(n)=1
           !
           ! diagnostic for initial tracers values
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On calcule la valeur des traceurs sur la position initiale

           CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is,js,ks)

           IF (key_alltracers) THEN

              IF (key_nointerpolstats) THEN
                 ti(n)=tt(is,js,ks,1)
                 si(n)=ss(is,js,ks,1)
                 ri(n)=rr(is,js,ks,1)
              ELSE
                 ti(n)=zinter(tt,gi(n),gj(n),gk(n),1._rprec)
                 si(n)=zinter(ss,gi(n),gj(n),gk(n),1._rprec)

                 IF (key_approximatesigma) THEN
                    ri(n)=zinter(rr,gi(n),gj(n),gk(n),1._rprec)

                 ELSE
                    ri(n)=sigma(zsigma,si(n),ti(n))

                 ENDIF ! key_approximatesigma

              ENDIF ! key_nointerpolstats

           ENDIF ! key_alltracers

           !-------------------------------!
           !-QUANTITATIVE INITIALIZATIONS -!
           !-------------------------------!
           IF (TRIM(mode) == 'quantitative') THEN

              ! BBL __________________________________________________________________________
              ! BBL -- BBL: En mode QUANTITATIVE, on met a jour si necessaire la valeur des
              ! BBL          champs de transport 2D au voisinage de la section initiale
              ! BBL

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

              !! NG: 15_09_2008
              !! NG: Here we stored the depth of the particle at the intial position.
              !! NG: Because with ROMS or symphonie the zz_ww in mod_fz is different at
              !! NG: time step.

              init_depth(n)=prof


              CALL sub_reducmem_shift_or_not_ind( NINT(tfi(n)), NINT(tfj(n)), NINT(tfk(n)), &
                   nis, njs, nks)

              IF ((ABS(ANINT(tfi(n))-tfi(n)) <= 1.e-5_rprec).AND.(uu(nis, js, ks,1) > 0._rprec)) THEN
                 IF (key_alltracers) THEN
                    CALL sub_quant_initside_u(NINT(tfi(n)), j0(n), k0(n), 1, &
                         prof, ttr(n), ti(n), si(n), ri(n))
                 ELSE
                    CALL sub_quant_initside_u(NINT(tfi(n)), j0(n), k0(n), 1, &
                         prof, ttr(n))
                 ENDIF
              ENDIF

              IF ((ABS(ANINT(tfj(n))-tfj(n)) <= 1.e-5_rprec).AND.(vv(is,njs,ks,1) > 0._rprec)) THEN
                 IF (key_alltracers) THEN
                    CALL sub_quant_initside_v(i0(n), NINT(tfj(n)), k0(n), 1, &
                         prof, ttr(n), ti(n), si(n), ri(n))
                 ELSE
                    CALL sub_quant_initside_v(i0(n), NINT(tfj(n)), k0(n), 1, &
                         prof, ttr(n))
                 ENDIF
              ENDIF

              IF ((ABS(ANINT(tfj(n))-tfj(n)) <= 1.e-5_rprec).AND.(ww(is,js,nks,1) < 0._rprec)) THEN
                 CALL sub_quant_initside_w(i0(n), j0(n), NINT(tfk(n)), 1, ttr(n))
              ENDIF

           ENDIF

           !---------------!
           !- QUALITATIVE -!
           !---------------!
           IF (TRIM(mode) == 'qualitative') THEN

              !------------------------------------------------!
              !- the initial position is written on "traj.ql" -!
              !------------------------------------------------!
              ! BBL ________________________________________________________________________
              ! BBL -- BBL: Calcul de la position geographique initiale
              ! BBL         La position verticale n'est pas calculee dans le cas de flotteurs
              ! BBL          isobares, mais elle est prise egale a la valeur de reference
              ! BBL          calculee a l'initialisation
              ! BBL
              x0 = fx(gi(n),gj(n))
              y0 = fy(gi(n),gj(n))

              IF (zbuoy(n) > 1.e8_rprec) THEN
                 IF (key_roms.OR.key_symphonie) THEN
                    z0= fz(gi=gi(n), gj=gj(n), gk=gk(n), il=1)
                 ELSE
                    z0=fz(gi(n),gj(n),gk(n))
                 ENDIF
              ELSE
                 z0=zbuoy(n)
              ENDIF

              ind_traj(n) = 1

              array_qual_traj(n,ind_traj(n),1) = x0
              array_qual_traj(n,ind_traj(n),2) = y0
              array_qual_traj(n,ind_traj(n),3) = z0
              array_qual_traj(n,ind_traj(n),4) = 0._rprec

              IF (key_alltracers) THEN

                 !!            WRITE(lun_traj,7055)n,x0,y0,z0,0.,ti(n),si(n),ri(n)

                 array_qual_traj(n,ind_traj(n),5) = ti(n)
                 array_qual_traj(n,ind_traj(n),6) = si(n)
                 array_qual_traj(n,ind_traj(n),7) = ri(n)

                 !!          ELSE
                 !!           WRITE(lun_traj,7055)n,x0,y0,z0,0.

              ENDIF ! key_alltracers

           ENDIF

        ENDIF

     ENDDO

2222 CONTINUE

     !!----------------------------------------------------------------------------!!
     !!---2222222222222222222---- DO LOOP : ON TRAJECTORIES ----222222222222222----!!
     !!----------------------------------------------------------------------------!!

     DO  n = 1, ntraj

        IF (igo(n) == 1) THEN
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
           ! BBL         (etat IGO=1)
           !
           ! linear interpolation for the local velocity
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Le champ de transport (UU, VV, WW) est interpole lineairement,
           ! BBL         composante par composante, sur la position courante (GI, GJ, GK)
           ! BBL
           CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is,js,ks)

           u(n) = uu(is-1,js,ks,1) + &
                (gi(n) - REAL(i0(n)-1,kind=rprec)) * (uu(is,js,ks,1)   - uu(is-1,js,ks,1))

           v(n) = vv(is,js-1,ks,1) + &
                (gj(n) - REAL(j0(n)-1,kind=rprec)) * (vv(is,js,ks,1)   - vv(is,js-1,ks,1))

           w(n) = ww(is,js,ks,1)   - &
                (gk(n) + REAL(k0(n),kind=rprec))   * (ww(is,js,ks+1,1) - ww(is,js,ks,1))
           !
           ! entry and exit indices (depending on the sign of the velocity)
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Pour chaque composante du transport, on calcule les indices des
           ! BBL         facettes *ENTREE* et *SORTIE*, a partir du signe de la composante
           ! BBL         du transport interpolee localement
           ! BBL
           i1(n) = i0(n) - 1
           i2(n) = i0(n)
           IF (u(n) < 0._rprec) THEN
              i2(n) = i0(n) - 1
              i1(n) = i0(n)
           ENDIF

           j1(n) = j0(n) - 1
           j2(n) = j0(n)
           IF (v(n) < 0._rprec) THEN
              j2(n) = j0(n) - 1
              j1(n) = j0(n)
           ENDIF

           k1(n) = k0(n) + 1
           k2(n) = k0(n)
           IF (w(n) < 0._rprec) THEN
              k2(n) = k0(n) + 1
              k1(n) = k0(n)
           ENDIF

           !!NG      ENDIF
           !!NG    ENDDO

           !!----------------------------------------------------------------------------!!
           !!---3333333333333333333---- DO LOOP : ON TRAJECTORIES ---3333333333333333----!!
           !!----------------------------------------------------------------------------!!

           !!NG    DO  n=1,ntraj
           !!NG      IF (igo(n) == 1) THEN

           ! BBL __________________________________________________________________________
           ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
           ! BBL         (etat IGO=1)
           !
           ! times to get to each edge of the (i, j, k) grid cell:
           ! - infinite, if Uexit * Ulocal < 0
           ! - linear relationship, if Uexit = Ulocal
           ! - logarithmic formulation, in all other cases
           !  BEWARE: use of log(u)-log(v) instead of log (u/v),
           !   better for numerical results
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: DU, DV et DW sont des tableaux intermediaires
           ! BBL
           CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is ,js ,ks )
           CALL sub_reducmem_shift_or_not_ind(i1(n),j1(n),k1(n),i1s,j1s,k1s)
           CALL sub_reducmem_shift_or_not_ind(i2(n),j2(n),k2(n),i2s,j2s,k2s)

           du(n) = ( uu(i2s,js ,ks ,1) - uu(i1s,js ,ks ,1) ) * REAL(i2(n)-i1(n),kind=rprec)
           dv(n) = ( vv(is ,j2s,ks ,1) - vv(is ,j1s,ks ,1) ) * REAL(j2(n)-j1(n),kind=rprec)
           dw(n) = ( ww(is ,js ,k2s,1) - ww(is ,js ,k1s,1) ) * REAL(k1(n)-k2(n),kind=rprec)

           ! BBL __________________________________________________________________________
           ! BBL - BBL: Les temps ARIANE (TX, TY, TZ, TTT) sont normalises, et le facteur
           ! BBL         de normalisation pour obtenir des temps *physiques* fait
           ! BBL         intervenir le volume de la maille TEMPERATURE contenant chaque
           ! BBL         particule
           ! BBL
           !!NG        tfac(n)=e3t(is, js, ks,1)*e1t(is, js,1,1)*e2t(is, js,1,1)

           IF (key_roms.OR.key_symphonie) THEN

              tfac(n) = e3t(is,js,ks,1) * e1t(is,js,1,1) * e2t(is,js,1,1)

           ELSE
              !- OPA case => E3T is not dependent of time -!
              IF (key_partialsteps) THEN
                 tfac(n) = e3t(is,js,ks,1) * e1t(is,js,1,1) * e2t(is,js,1,1)
              ELSE
                 tfac(n) = e3t(1,1,ks,1) * e1t(is,js,1,1) * e2t(is,js,1,1)
              ENDIF

           ENDIF

           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Si le signe de la composante zonale du transport est different du
           ! BBL         signe de la composante zonale du transport sur la facette de
           ! BBL         sortie, le temps d'acces vers cette facette est infini
           ! BBL
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Si la composante zonale du transport est la meme sur les facettes
           ! BBL         d'entree et de sortie, le temps d'acces est defini par une simple
           ! BBL         relation lineaire
           ! BBL
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Temps d'acces (vers la facette de sortie) dans le cas general
           ! BBL

           IF (( u(n) * uu(i2s,js,ks,1)) <= 0._rprec) THEN
              tx(n) = 1.e35_rprec
           ELSEIF (du(n) == 0._rprec) THEN
              tx(n) = ( REAL(i2(n),  kind=rprec) - gi(n) ) / u(n)
           ELSE
              tx(n) = ( LOG(ABS(uu(i2s,js,ks,1))) - LOG(ABS(u(n))) ) / du(n)
           ENDIF

           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Les calculs sont equivalents pour la composante meridienne du
           ! BBL         transport
           ! BBL

           IF ((v(n) * vv(is,j2s,ks,1)) <= 0._rprec) THEN
              ty(n) = 1.e35_rprec
           ELSEIF (dv(n) == 0._rprec) THEN
              ty(n) = ( REAL(j2(n), kind=rprec) - gj(n) ) / v(n)
           ELSE
              ty(n) = ( LOG(ABS(vv(is,j2s,ks,1))) - LOG(ABS(v(n))) ) / dv(n)
           ENDIF

           ! BBL __________________________________________________________________________
           ! BBL - BBL: Les calculs sont equivalents pour la composante verticale du
           ! BBL         transport, avec un traitement specifique des flotteurs isobares
           ! BBL         (pour lesquels les temps d'acces vers les facettes de sortie sont
           ! BBL         forces a une valeur *infinie*)
           ! BBL
           !!NG: pour BBL
           IF ((key_2dquant).OR.(zbuoy(n) <= 1.e8_rprec).OR.((w(n)*ww(is,js,k2s,1)) <= 0._rprec)) THEN
              tz(n) = 1.e35_rprec
           ELSEIF (dw(n) == 0._rprec) THEN
              tz(n) = -( REAL(k2(n), kind=rprec) + gk(n) ) / w(n)
           ELSE
              tz(n) = ( LOG(ABS(ww(is,js,k2s,1))) - LOG(ABS(w(n))) ) / dw(n)
           ENDIF

           !-----------------------------------------------------------------------------
           ! the true exit time is given by the shortest times (of tx, ty, tz)
           ! entry and exit indices are then updated
           !-----------------------------------------------------------------------------

           ! BBL __________________________________________________________________________
           ! BBL - BBL: Si une particule sort de la maille temperature qui la contient,
           ! BBL         elle le fait en un temps TTT, egal au minimum des 3 temps d'acces
           ! BBL         precedents (TX, TY, TZ)
           ! BBL
           ttt(n)=MIN(tx(n),ty(n),tz(n))

           !!NG: pour BBL 
           IF (.NOT.key_2dquant) THEN
              !
              ! BBL __________________________________________________________________________
              ! BBL - BBL: Dans certains cas, les 3 temps d'acces (TX, TY, TZ) sont infinis...
              ! BBL
              IF (ttt(n) >= 1.e34_qprec) THEN
                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: ... cela peut notamment se produire pour des particules isobares
                 ! BBL         (pour lesquelles on a force TZ a etre infini); dans ce cas, c'est
                 ! BBL         *normal* et on continue les calculs
                 ! BBL
                 !! NG: 12/05/2009 IF (zbuoy(n) <= 1.e8_rprec) GOTO 204
                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: ... si la particule n'est pas isobare, c'est qu'elle se trouve
                 ! BBL         piegee dans une maille ou les 6 composantes du transport sont
                 ! BBL         *entrantes*; cela peut se produire dans le cas d'un maillage
                 ! BBL         vertical qui *respire* (avec la montee et la descente de la
                 ! BBL         surface libre);
                 ! BBL        La valeur de TTT est ramenee a une valeur *finie*, choisie de telle
                 ! BBL         facon qu'elle excede la duree documentee par le champ de
                 ! BBL         transport present en memoire
                 ! BBL        NOTE: cette precaution est en fait inutile, et on pourrait tout
                 ! BBL         aussi bien continuer avec la valeur *infinie* de TTT, comme pour
                 ! BBL         les particules isobares; a la place de la ligne de code suivante,
                 ! BBL         il faudrait en fait introduire un test du type:
                 ! BBL          - si le maillage vertical peut evoluer dans le temps, une valeur
                 ! BBL           de TTT *infinie*  n'est pas incongrue, et on continue les
                 ! BBL           calculs sans modifier TTT
                 ! BBL          - si le maillage vertical est statique, TTT *infini* est une
                 ! BBL           aberration, et il faut stopper les calculs avec un message
                 ! BBL           d'erreur 'negative time'
                 ! BBL
                 !! NG: 12/05/2009 ttt(n)=2._rprec*(tswap(n)-age(n))/tfac(n)
                 IF (zbuoy(n) > 1.e8_rprec)  THEN
                    ttt(n)=2._qprec*(tswap(n)-age(n))/tfac(n)
                 ENDIF

              ENDIF

           ENDIF
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Dans de tres rares cas, TTT prend une valeur negative!
           ! BBL
204        IF (ttt(n) < 0._qprec) THEN
              ! BBL __________________________________________________________________________
              ! BBL - BBL: Si TTT est *significativement* negatif (> 1 seconde), il y a un
              ! BBL         probleme et on stoppe l'integration avec un message d'erreur
              ! BBL
              IF (ABS(tfac(n)*ttt(n)) > 1._rprec) THEN
                 WRITE(lun_error,7000)' negative time', &
                      n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
                 STOP
              ENDIF
              ! BBL __________________________________________________________________________
              ! BBL - BBL: Sinon, il s'agit d'un probleme d'arrondi machine, et on force TTT
              ! BBL         a etre egal a 0
              ! BBL
              ttt(n)=0._qprec
           ENDIF


           ! BBL __________________________________________________________________________
           ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
           ! BBL         (etat IGO=1)
           !
           ! do we need to update the velocity field WITHIN the gridcell ?
           !
           ! update needed
           !
           IF ((lmt > 1).AND.((age(n)+ttt(n)*tfac(n)) > tswap(n))) THEN
              ! BBL __________________________________________________________________________
              ! BBL - BBL: Si LMT > 1, le temps d'acces vers la facette de sortie NE SAURAIT
              ! BBL         exceder la duree documentee par le champ de transport present en
              ! BBL         memoire
              ! BBL        Si tel est le cas, on seuille TTT au temps que le champ present en
              ! BBL         memoire permet encore de documenter, et la particule passe dans
              ! BBL         l'etat IGO=2: une mise a jour du champ de transport s'avere
              ! BBL         necessaire avant de la faire *repartir*
              ! BBL        NOTE: pour comparer les temps ARIANE aux temps physiques (AGE,
              ! BBL         TSWAP, TFIN, ...), l'utilisation du facteur de normalisation est
              ! BBL         necessaire
              ! BBL
              ttt(n)=(tswap(n)-age(n))/tfac(n)
              igo(n)=2
           ENDIF
           !
           ! end of the integration
           !
           IF ((TRIM(mode) == 'qualitative').AND.((age(n)+ttt(n)*tfac(n)) > tfin)) THEN
              ! BBL __________________________________________________________________________
              ! BBL - BBL: Dans un cas QUALITATIVE, la valeur retenue pour TTT ne doit pas
              ! BBL         conduire a un age de la particule superieur au temps maximal
              ! BBL         d'integration demande
              ! BBL        Si tel est le cas, on seuille TTT au temps maximal admissible, et
              ! BBL         la particule passe dans l'etat IGO=3 (integration terminee)
              ! BBL
              ttt(n)=(tfin-age(n))/tfac(n)
              igo(n)=3
           ENDIF

           !====================================================================
           ! case 1: we do not reach the side of the cell (for a given direction)
           ! diagnostic of the final position
           ! - linear formulation, if Uexit = Ulocal
           ! - exponential formulation, in other cases
           !====================================================================
           IF (tx(n) > ttt(n)) THEN

              ! BBL __________________________________________________________________________
              ! BBL - BBL: On traite le cas ou le temps de sortie dans la direction zonale est
              ! BBL          plus grand que le temps TTT retenu pour les calculs
              ! BBL        On calcule donc une nouvelle position dans la maille TEMPERATURE,
              ! BBL         a partir de la valeur TTT (avec l'expression analytique de la
              ! BBL         trajectoire dans la direction zonale)
              ! BBL
              CALL sub_dont_reachside(hi(n), du(n), gi(n), u(n), ttt(n), i1(n), i2(n))

           ENDIF


           IF (ty(n) > ttt(n)) THEN

              ! BBL __________________________________________________________________________
              ! BBL - BBL: On applique un traitement equivalent dans la direction meridienne
              ! BBL
              CALL sub_dont_reachside(hj(n), dv(n), gj(n), v(n), ttt(n), j1(n), j2(n))

           ENDIF


           IF (TRIM(mode) == 'quantitative') THEN
              !!NG: BBL
              IF ((zbuoy(n) > 1.e8_rprec).AND.(.NOT.key_2dquant)) THEN

                 IF (tz(n) > ttt(n)) THEN
                    ! BBL __________________________________________________________________________
                    ! BBL -- BBL: On applique un traitement equivalent dans la direction verticale,
                    ! BBL          si la particule n'est pas isobare
                    ! BBL
                    CALL sub_dont_reachside(hk(n), dw(n), gk(n), w(n), ttt(n), -k1(n), -k2(n))

                 ENDIF
              ELSEIF(key_2dquant) THEN

                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: Si la particule est isobare, on ne fait aucun calcul et on force la
                 ! BBL         conservation de l'indice de maillage vertical
                 ! BBL
                 hk(n)=gk(n)

              ENDIF

           ELSE !! QUALITATIVE !!
              IF (zbuoy(n) > 1.e8_rprec) THEN

                 IF (tz(n) > ttt(n)) THEN
                    ! BBL __________________________________________________________________________
                    ! BBL -- BBL: On applique un traitement equivalent dans la direction verticale,
                    ! BBL          si la particule n'est pas isobare
                    ! BBL
                    CALL sub_dont_reachside(hk(n), dw(n), gk(n), w(n), ttt(n), -k1(n), -k2(n))

                 ENDIF
              ELSE

                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: Si la particule est isobare, on ne fait aucun calcul et on force la
                 ! BBL         conservation de l'indice de maillage vertical
                 ! BBL
                 hk(n)=gk(n)

              ENDIF
           ENDIF

           !====================================================================
           ! case 2: we DO reach the side of the cell (for a given direction)
           ! no need for time update within the gridcell
           ! - we record the exiting transport
           ! - we update the positions index (hi, hj, hk)
           !====================================================================

           IF (tx(n) <= ttt(n)) THEN

              ! BBL __________________________________________________________________________
              ! BBL -- BBL: On traite le cas ou le temps d'acces vers la section de sortie
              ! BBL         *zonale* (TX) est inferieur ou egal au temps TTT retenu pour
              ! BBL         l'integration de la particules
              ! BBL         NOTE: TX ne peut etre en fait STRICTEMENT inferieur a TTT (par
              ! BBL          construction), et le test capture en fait les cas TX=TTT
              ! BBL __________________________________________________________________________
              ! BBL -- BBL: Le nouvel indice de positionnement zonal est EXACTEMENT l'indice
              ! BBL          de la facette de sortie
              ! BBL
              hi(n)=REAL(i2(n),kind=rprec)

              IF (TRIM(mode) == 'quantitative') THEN

                 ! BBL __________________________________________________________________________
                 ! BBL BBL BBL: Dans le cas QUANTITATIVE, on memorise le passage de la particule
                 ! BBL           dans les champs de transport lagrangien 2D appropriés
                 ! BBL

                 CALL sub_quant_reachside_u( &
                      i2(n), j0(n), k0(n), 1, &
                      hi(n), hj(n), hk(n), 1._rprec, &
                      ttr(n))

              ENDIF

              !
              ! BBL __________________________________________________________________________
              ! BBL -- BBL: On effectue la mise a jour des indices des facettes *ENTREE* et
              ! BBL          *SORTIE* dans la direction zonale; cette mise a jour depend bien
              ! BBL           sur du sens du deplacement
              ! BBL
              IF (i2(n) > i1(n)) THEN
                 i1(n) = i2(n)
                 i2(n) = i2(n)+1
              ELSE
                 i1(n) = i2(n)
                 i2(n) = i2(n)-1
              ENDIF

           ENDIF

           IF (ty(n) <= ttt(n)) THEN

              ! BBL __________________________________________________________________________
              ! BBL -- BBL: On applique un traitement equivalent dans la direction meridienne
              ! BBL
              hj(n) = REAL(j2(n),kind=rprec)

              IF (TRIM(mode) == 'quantitative') THEN

                 CALL sub_quant_reachside_v(         &
                      i0(n), j2(n), k0(n), 1,        &
                      hi(n), hj(n), hk(n), 1._rprec, &
                      ttr(n))

              ENDIF

              IF (j2(n) > j1(n)) THEN
                 j1(n) = j2(n)
                 j2(n) = j2(n) + 1
              ELSE
                 j1(n) = j2(n)
                 j2(n) = j2(n) - 1
              ENDIF

           ENDIF


           IF (tz(n) <= ttt(n)) THEN

              ! BBL __________________________________________________________________________
              ! BBL -- BBL: On applique un traitement equivalent dans la direction verticale,
              ! BBL          sauf pour les particules isobares pour lesquelles aucun
              ! BBL          franchissement vertical de maille n'est necessaire
              ! BBL
              IF (TRIM(mode) == 'quantitative') THEN

                 CALL sub_quant_reachside_w(i0(n), j0(n), k2(n), 1, ttr(n))

              ENDIF

              !!NG: BBL
              IF (.NOT.key_2dquant) THEN

                 IF (zbuoy(n) > 1.e8_rprec) THEN

                    hk(n) = -REAL(k2(n),kind=rprec)

                    IF (k2(n) > k1(n)) THEN
                       k1(n) = k2(n)
                       k2(n) = k2(n) + 1
                    ELSE
                       k1(n) = k2(n)
                       k2(n) = k2(n) - 1
                    ENDIF

                 ENDIF

              ENDIF

           ENDIF


           ! BBL __________________________________________________________________________
           ! BBL - BBL: Le cas des particules isobares necessite un traitement particulier
           ! BBL        En effet, pour un systeme de coordonnees suivant la topographie, si
           ! BBL         l'on se deplace dans une couche donnee sans modifier l'indice
           ! BBL         vertical (HK), on ne conserve pas forcement l'immersion de
           ! BBL         reference
           ! BBL        Ainsi, il convient de verifier que la maille TEMPERATURE dans
           ! BBL         laquelle la particule evolue couvre toujours un domaine de
           ! BBL         profondeur incluant l'immersion de reference a laquelle on
           ! BBL         souhaite faire voyager la particule isobare
           ! BBL
           IF (zbuoy(n) < 1.e8_rprec) THEN
              iflag(n)=0
              ! BBL __________________________________________________________________________
              ! BBL -- BBL: Dans le contexte le plus simple, l'extension verticale de la
              ! BBL          maille ambiante inclut bien l'immersion de reference, et il n'est
              ! BBL          pas necessaire de modifier l'indice de positionnement vertical
              ! BBL 

              !NG: potential bug : change <0 by <=0._rprec (15/05/2009)

              IF (                                                                      &
                   (zbuoy(n) - fz( gi=hi(n), gj=hj(n), gk=-AINT(-hk(n)), il=1))       * &
                   (zbuoy(n) - fz( gi=hi(n), gj=hj(n), gk=-1._rprec-AINT(-hk(n)), il=1))&
                   <= 0._rprec) THEN
                 hk(n)=gk(n)
                 iflag(n)=1
              ELSE
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: Dans le cas contraire, il convient de scanner la verticale pour
                 ! BBL          trouver la nouvelle maille TEMPERATURE dont l'extension verticale
                 ! BBL          inclut l'immersion de reference de la particule isobare
                 ! BBL

                 DO kk0=1,kmt-1

                    !NG: potential bug : change <0 by <=0._rprec (15/05/2009)

                    IF (                                                                           &
                         (zbuoy(n) - fz(gi=hi(n),gj=hj(n),gk=-REAL(kk0, kind=rprec),il=1))       * &
                         (zbuoy(n) - fz(gi=hi(n),gj=hj(n),gk=-1._rprec-REAL(kk0,kind=rprec),il=1)) &
                         <= 0._rprec) THEN
                       hk(n)=-REAL(kk0, kind=rprec)-.5_rprec
                       k1(n)=kk0
                       k2(n)=kk0+1
                       iflag(n)=1
                    ENDIF
                 END DO

              ENDIF

              IF (iflag(n) == 0) THEN
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: La particule qui ne peut plus etre positionnee dans une maille
                 ! BBL          TEMPERATURE appropriee (incluant son immersion de reference) voit
                 ! BBL           son integration se finir, et son etat passe a IGO=3
                 ! BBL
                 WRITE(lun_error,7000)'out of zdomain', &
                      n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
                 igo(n)=3
              ENDIF
           ENDIF
           !
        ENDIF

     ENDDO

     !!----------------------------------------------------------------------------!!
     !!---5555555555555555555---- DO LOOP : ON TRAJECTORIES ---5555555555555555----!!
     !!----------------------------------------------------------------------------!!

     DO n=1,ntraj
        IF (igo(n) > 0) THEN
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
           ! BBL         (etat IGO=1,2,3)
           !
           ! update of the T-grid index of the new embedding gridcell
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Il convient de mettre a jour les indices (I0, J0, K0) de la
           ! BBL         nouvelle maille TEMPERATURE contenant la particule, a partir de
           ! BBL         la connaissance des nouveaux indices des facettes d'entree et de
           ! BBL         sortie
           ! BBL
           i0(n)    = max0(i1(n), i2(n))
           j0(n)    = max0(j1(n), j2(n))
           k0(n)    = min0(k1(n), k2(n))

           !===========================!
           !- Periodicity in OPA-ORCA -!
           !===========================!
           IF (.NOT.(key_roms.OR.key_symphonie)) THEN

              CALL sub_orca_north_pole(i1(n), j1(n), j0(n)+1)
              CALL sub_orca_north_pole(i2(n), j2(n), j0(n)+1)
              CALL sub_orca_north_pole(hi(n), hj(n), j0(n)+1)
              CALL sub_orca_north_pole(i0(n), j0(n), j0(n)+1, jperfl(n))

              !!NG: 26/06/2009: BUGGED CALL sub_orca_east_west_periodic(i1(n), min0(i1(n),i2(n))+1, i0(n))
              !!NG: 26/06/2009: BUGGED CALL sub_orca_east_west_periodic(i2(n), min0(i1(n),i2(n))+1, i0(n))
              !!NG: 26/06/2009: BUGGED CALL sub_orca_east_west_periodic(hi(n), min0(i1(n),i2(n))+1, i0(n))
              !!NG: 26/06/2009: BUGGED CALL sub_orca_east_west_periodic(i0(n), min0(i1(n),i2(n))+1, i0(n), iperfl(n))

              !!NG: 26/06/2009: Replace by:
              imin_tmp = min0(i1(n),i2(n))
              imax_tmp = max0(i1(n),i2(n))
              CALL sub_orca_east_west_periodic(i1(n), imin_tmp+1, imax_tmp)
              CALL sub_orca_east_west_periodic(i2(n), imin_tmp+1, imax_tmp)
              CALL sub_orca_east_west_periodic(hi(n), imin_tmp+1, imax_tmp)
              CALL sub_orca_east_west_periodic(i0(n), imin_tmp+1, imax_tmp, iperfl(n))

           ENDIF

           !====================================================================
           ! age update
           !====================================================================
           !
           ! BBL __________________________________________________________________________
           ! BBL - BBL: L'age des particules (en secondes) est mis a jour a partir de la
           ! BBL         valeur de TTT
           ! BBL
           agenew(n)=age(n)+ttt(n)*tfac(n)

           ! BBL __________________________________________________________________________
           ! BBL - BBL: On met a jour l'indice temporel des particules, selon le sens de
           ! BBL         l'integration (BACKWARD ou FORWARD)
           ! BBL
           IF (TRIM(forback) == 'forward') THEN
              hl(n) = gl(n) + ttt(n) * tfac(n) / tfic
           ELSE
              hl(n) = gl(n) - ttt(n) * tfac(n) / tfic
           ENDIF

           ! BBL __________________________________________________________________________
           ! BBL - BBL: Si une particule a ete stoppee dans l'attente d'une mise a jour du
           ! BBL         champ de transport charge en memoire, on force la valeur de
           ! BBL         l'indice temporel a l'instant exact de la transition
           ! BBL
           IF (igo(n) == 2)  THEN
              hl(n)     = NINT(hl(n)+.5_rprec) - .5_rprec
              agenew(n) = tswap(n)
           ENDIF

        ENDIF

     ENDDO


     !!================================================================================!!
     !!============================== QUALITATIVE MODE ================================!!
     !!================================================================================!!
     !
     !====================================================================
     ! we stop the integration if we reach a limit of the domain (it may
     !  happen in QUALITATIVE experiments) [QUALI]
     !====================================================================
     IF (TRIM(mode) == 'qualitative') THEN

        !!-------------------------------------!!
        !!------ DO LOOP: ON TRAJECTORIES -----!!
        !!-------------------------------------!!
        DO  n=1,ntraj

           IF (igo(n) >= 1) THEN
              ! BBL _________________________________________________________________
              ! BBL - BBL: On ne fait les calculs suivants que 
              ! BBL        pour les particules *actives* (etat IGO=1)
              !

              IF (k0(n) == 0) THEN
                 ! BBL ______________________________________________________________
                 ! BBL -- BBL: Traitement des sorties de domaine *par la surface*
                 ! BBL
                 WRITE(lun_error,7000)'out of zdomain', &
                      n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
                 igo(n)=3
              ENDIF

              IF ((igo(n) /= 3).AND.((j0(n) == dims_reg(2,1)).OR.(j0(n) == dims_reg(2,2)))) THEN
                 ! BBL _______________________________________________________________
                 ! BBL -- BBL: Traitement des sorties de domaine 
                 ! BBL         *par une frontiere meridienne*
                 ! BBL
                 WRITE(lun_error,7000)'out of ydomain', &
                      n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
                 igo(n)=3
              ENDIF

              !!NG          IF (.NOT.key_jfold) THEN
              !!NG            IF ((igo(n) /= 3).AND.(j0(n) == dims_reg(2,2))) THEN
              !!NG              ! BBL ______________________________________________________________
              !!NG              ! BBL -- BBL: Il y a un BUG ici, car le test suivant n'est jamais
              !!NG              ! BBL         verifie (le test precedent a deja fait le travail)
              !!NG              ! BBL         Il faudrait reecrire les tests de maniere plus propre 
              !!NG              ! BBL         dans le cas d'une condition de repliement meridien
              !!NG              ! BBL
              !!NG              WRITE(lun_error,7000)'out of Ydomain', &
              !!NG                   n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
              !!NG              igo(n)=3
              !!NG            ENDIF
              !!NG          ENDIF ! .NOT.key_jfold

              IF (.NOT.key_periodic) THEN
                 IF ((igo(n) /= 3).AND.((i0(n) == dims_reg(1,1)).OR.(i0(n) == dims_reg(1,2)))) THEN
                    ! BBL _________________________________________________________
                    ! BBL -- BBL: Traitement des sorties de domaine 
                    ! BBL         *par une frontiere zonale*
                    ! BBL
                    WRITE(lun_error,7000)'out of xdomain', &
                         n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
                    igo(n)=3
                 ENDIF
              ENDIF ! .NOT.key_periodic

           ENDIF
        ENDDO
        !
     ENDIF ! (TRIM(mode) == 'qualitative') => QUALITATIVE MODE

     !!================================================================================!!
     !!============================== QUANTITATIVE MODE ===============================!!
     !!================================================================================!!
     !
     !====================================================================
     ! quantitative analysis
     !
     ! in QUANTITATIVE experiments, the 2D projections of the transport
     !  field associated with 1 particle is summed depending on the final
     !  section [QUANT]
     !====================================================================
     ! BBL __________________________________________________________________________
     ! BBL BBL: En mode QUANTITATIVE, on teste l'atteinte des sections d'interception
     ! BBL
     IF (TRIM(mode) == 'quantitative') THEN

        !!-------------------------------------!!
        !!------ DO LOOP: ON TRAJECTORIES -----!!
        !!-------------------------------------!!
        DO n=1,ntraj

           IF (igo(n) > 0) THEN
              ! BBL __________________________________________________________________________
              ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
              ! BBL         (etat IGO=1,2,3)
              !------------------------------------------------------------------
              ! mtfin informs about the sections used in quantitative experiments
              !------------------------------------------------------------------
              nfin(n) = mtfin(i0(n),j0(n),k0(n))

              !!NG: REMARK - nothing here about criter2 like in trajec.f90 ?????

              !===============================================================
              ! a final criterion may define a "final section" (criter1)
              !===============================================================
              IF ((nfin(n) > 0).OR. &
                   (criter1(tfi(n),tfj(n),tfk(n),tfl(n), &
                   hi(n),hj(n),hk(n),hl(n),i0(n),j0(n),k0(n),l0(n)))) THEN

                 !----------------------------------------------------------
                 ! a final criterion uses "1 + nsect" as the section index
                 ! we flag the existence of such a criterion
                 !----------------------------------------------------------
                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: Quand une particule atteint une section d'interception ou verifie
                 ! BBL        critere d'interception (CRITER1), son integration s'arrete et son
                 ! BBL        etat passe a IGO=3
                 ! BBL
                 igo(n)=3
                 IF (nfin(n) <= 0) THEN
                    nfin(n) = 1 + nsect
                    icrit1 = 1
                 ENDIF
                 !-----------------------------------------------------------------
                 ! we distinguish meandering and recirculating particles (criter0)
                 !-----------------------------------------------------------------
                 IF ((TRIM(bin) == 'nobin').AND.(nfin(n) == 1) .AND. &
                      criter0(hi(n),hj(n),hk(n),hl(n),i0(n),j0(n),k0(n),l0(n))) THEN

                    nfin(n) = 0

                 ENDIF
                 !
                 ! hereafter,
                 ! nfin=0, if the initial section (with criter0) has been reached
                 ! nfin=1, if the initial section (no criter0) or a final section has
                 !  been reached
                 ! nfin=1+nsect, if a final criterion (criter1) is satisfied
                 !
                 sectrans(nfin(n))=sectrans(nfin(n))+ttr(n)

                 !???????????????????????????????????????????????????????????????????
                 ! there IS NOW storage of transport projections for "criter1" cases
                 !  (but never try to compute 2D streamfunctions with these fields !)
                 !???????????????????????????????????????????????????????????????????

                 !==================================================================
                 ! statistics...
                 !==================================================================

                 IF (key_alltracers) THEN
                    CALL sub_quant_statistics(               &
                         nfin(n), ttr(n), agenew(n), tcyc  , &
                         hi(n)  , hj(n) , hk(n) , 1._rprec , &
                         tfi(n) , tfj(n), tfk(n), 1._rprec , &
                         i0(n)  , j0(n) , k0(n) , 1        , &
                         tf(n)  , sf(n) , rf(n) , truage(n), &
                         ti(n), si(n), ri(n))
                 ELSE
                    CALL sub_quant_statistics(               &
                         nfin(n), ttr(n), agenew(n), tcyc  , &
                         hi(n)  , hj(n) , hk(n) , 1._rprec , &
                         tfi(n) , tfj(n), tfk(n), 1._rprec , &
                         i0(n)  , j0(n) , k0(n) , 1        , &
                         tf(n)  , sf(n) , rf(n) , truage(n)  )
                 endif

                 !! NG: 15_09_2008
                 !! NG: Here we store the final depth which depend of time for ROMS
                 !! NG: and symphonie.
                 IF (key_roms.OR.key_symphonie) THEN
                    prof=fz(gi=hi(n), gj=hj(n), gk=hk(n), il=1)
                 ELSE
                    prof=fz(gk=hk(n))
                 ENDIF
                 final_depth(n)=prof

              ENDIF

           ENDIF

        ENDDO

     ENDIF ! end of quantitative analysis

     !!----------------------------------------------------------------------------!!
     !!---6666666666666666666---- DO LOOP : ON TRAJECTORIES ---6666666666666666----!!
     !!----------------------------------------------------------------------------!!
     DO n=1,ntraj

        IF (igo(n) == 2) THEN
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On applique ici un traitement specifique aux particules dans
           ! BBL         l'attente d'une mise a jour en memoire du champ de transport
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On incremente (ou decremente) d'une unite l'indice temporel
           ! BBL
           IF (TRIM(forback) == 'forward') THEN
              l0(n) = l0(n) + 1
           ELSE
              l0(n) = l0(n) - 1
           ENDIF
           ! BBL __________________________________________________________________________
           ! BBL - BBL: Si necessaire, on boucle sur le cycle periodique
           ! BBL
           IF (l0(n) < 1)   l0(n) = l0(n) + lmt
           IF (l0(n) > lmt) l0(n) = l0(n) - lmt
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On met a jour pour chaque particule l'instant a utiliser pour le
           ! BBL         prochain swap temporel en memoire (prochaine mise a jour du champ
           ! BBL         de transport)
           ! BBL
           tswap(n)=tswap(n)+tfic
        ENDIF

     ENDDO

     !!----------------------------------------------------------------------------!!
     !!---7777777777777777777---- DO LOOP : ON TRAJECTORIES ---7777777777777777----!!
     !!----------------------------------------------------------------------------!!
     DO  n=1,ntraj

        IF ((igo(n) > 0).AND.(igo(n) /= 3)) THEN
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On applique ici un traitement specifique aux particules dont
           ! BBL         l'integration n'est pas encore terminee (mais deja commencee)
           ! BBL        On verifie que l'on n'a pas atteint un point TERRE

           CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is,js,ks)

           IF (tmask(is,js,ks,1) == 0) THEN
              WRITE(lun_error,7000)'   coast crash', &
                   n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
              WRITE(lun_output,*)' '
              WRITE(lun_output,*)'coast crash for particle # ',n
              WRITE(lun_output,*)'i j k: ',i0(n),' ',j0(n),' ',k0(n)
              WRITE(lun_output,*)'is js ks: ',is,' ',js,' ',ks
              WRITE(lun_output,*)'gi gj gk:',gi(n),' ',gj(n),' ',gk(n)
              WRITE(lun_output,*)'hi hj hk:',hi(n),' ',hj(n),' ',hk(n)
              !! NG: Bug          WRITE(lun_output,*)'u: ',uu(i0(n)-1,j0(n),k0(n),1),uu(i0(n),j0(n),k0(n),1)
              !! NG: Bug          WRITE(lun_output,*)'v: ',vv(i0(n),j0(n)-1,k0(n),1),vv(i0(n),j0(n),k0(n),1)
              !! NG: Bug          WRITE(lun_output,*)'w: ',ww(i0(n),j0(n),k0(n)+1,1),ww(i0(n),j0(n),k0(n),1)
              WRITE(lun_output,*)'u: ',uu(is-1,js,ks,1),uu(is,js,ks,1)
              WRITE(lun_output,*)'v: ',vv(is,js-1,ks,1),vv(is,js,ks,1)
              WRITE(lun_output,*)'w: ',ww(is,js,ks+1,1),ww(is,js,ks,1)
              WRITE(lun_output,*)'tx ty tz ttt: ',tx(n),ty(n),tz(n),ttt(n)
              STOP
           ENDIF
           !
        ENDIF

     ENDDO !  n=1,ntraj

     !!================================================================================!!
     !!============================== QUALITATIVE MODE ===============================!!
     !!================================================================================!!
     !
     ! do we need an output of the position (qualitative experiment) ?
     ! (linear interpolation)
     ! - age: previous age
     ! - agenew: current age
     ! - tnext: output time
     !
     IF (TRIM(mode) == 'qualitative') THEN
        ! BBL __________________________________________________________________________
        ! BBL BBL: En mode QUALITATIVE, on s'interesse ici au calcul des positions
        ! BBL        GEOGRAPHIQUES des particules
        ! BBL       ARIANE a en memoire une position ancienne (GI, GJ, GK, GL) et une
        ! BBL        position nouvelle (HI, HJ, HK, HL) pour chaque particule active, et
        ! BBL         on sait qu'il faut calculer une position geographique toutes les
        ! BBL        TNEXT secondes
        ! BBL       Les valeurs de GL et HL n'etant pas les memes pour toutes les
        ! BBL        particules, on travaille ici de maniere incrementale, et les
        ! BBL        calculs se font par incrementation temporelle, *tant que cela
        ! BBL        s'avere necessaire*
        ! BBL       Par defaut, toutes les particules sont *bonnes pour le service*
        ! BBL        (IFLAG=1)
        ! BBL

        iflag(:) = 1 ! array notation

        !----------
        ! BBL __________________________________________________________________________
        ! BBL BBL: La variable IFL sera mise a la valeur 1 si une particule au moins a
        ! BBL       besoin de voir echantillonner sa trajectoire GEOGRAPHIQUE
        ! BBL
1111    ifl=0

        !!-------------------------------------!!
        !!------ DO LOOP: ON TRAJECTORIES -----!!
        !!-------------------------------------!!
        DO n=1,ntraj

           IF (igo(n) > 0) THEN
              ! BBL __________________________________________________________________________
              ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
              ! BBL         (etat IGO=1,2,3)
              !
              IF ((iflag(n) == 1).AND.(agenew(n) >= tnext(n))) THEN
                 ! BBL __________________________________________________________________________
                 ! BBL - BBL: Si l'age *nouveau* AGENEW d'une particule excede l'instant TNEXT
                 ! BBL         correspondant a la demande d'une sortie GEOGRAPHIQUE, on calcule
                 ! BBL         une position par interpolation entre les positions (GI, GJ, GK) et
                 ! BBL         (HI, HJ, HK), avec un facteur de ponderation FT lie a
                 ! BBL         l'emplacement temporel de TNEXT entre AGE et AGENEW
                 ! BBL        NOTE: XT, YT et ZT ne sont pas des positions GEOGRAPHIQUES, mais
                 ! BBL         des indices sur la grille du modele, tout comme GI, GJ et GK, et
                 ! BBL         HI, HJ et HK
                 ! BBL
                 ft=(tnext(n)-age(n))/(agenew(n)-age(n))
                 xt=gi(n)+(hi(n)-gi(n))*ft
                 yt=gj(n)+(hj(n)-gj(n))*ft
                 zt=gk(n)+(hk(n)-gk(n))*ft

                 IF (TRIM(forback) == 'forward') THEN
                    st=gl(n)+(tnext(n)-age(n))/tfic
                 ELSE
                    st=gl(n)-(tnext(n)-age(n))/tfic
                 ENDIF

                 IF (key_periodic) THEN
                    ! BBL _____________________________________________________________________
                    ! BBL -- BBL: La variable IPERFL est utilisee pour corriger le cas echeant
                    ! BBL          l'interpolation lineaire precedente (donnant un resultat 
                    ! BBL          faux si les deux positions de reference, GI et HI, sont 
                    ! BBL          situes de part et d'autre de la periodicite
                    ! BBL
                    IF (iperfl(n) == 1) THEN
                       xt=gi(n)+(hi(n)-REAL(imt-2,kind=rprec)-gi(n))*ft
                    ENDIF
                    IF (iperfl(n) == 2) THEN
                       xt=gi(n)+(hi(n)+REAL(imt-2,kind=rprec)-gi(n))*ft
                    ENDIF

                 ENDIF


                 IF (key_jfold) THEN

                    ! extrapolation not yet implemented
                    IF (jperfl(n) == 1) THEN
                       xt=gi(n)
                       yt=gj(n)
                       zt=gk(n)
                    ENDIF

                 ENDIF


                 !---------------------------!
                 !- position output on file -!
                 !---------------------------!
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: On effectue ici la conversion entre indice de grille et
                 ! BBL          coordonnees geographiques
                 ! BBL
                 x0=fx(xt,yt)
                 y0=fy(xt,yt)

                 IF (zbuoy(n) > 1.e8_rprec) THEN
                    IF (key_roms.OR.key_symphonie) THEN
                       z0 = fz(gi=xt, gj=yt, gk=zt, il = 1)
                    ELSE
                       z0 = fz(gk=zt)
                    ENDIF
                 ELSE

                    ! BBL ________________________________________________________________________
                    ! BBL -- BBL: Pour les particules isobares, aucun calcul n'est necessaire pour
                    ! BBL          l'immersion...
                    ! BBL
                    z0=zbuoy(n)
                 ENDIF

                 IF (tnext(n) <= tfin) THEN
                    ! BBL _______________________________________________________________________
                    ! BBL -- BBL: L'ecriture sur fichier se fait si on n'a pas atteint le temps
                    ! BBL          limite TFIN

                    ind_traj(n) = ind_traj(n) + 1

                    array_qual_traj(n,ind_traj(n),1) = x0
                    array_qual_traj(n,ind_traj(n),2) = y0
                    array_qual_traj(n,ind_traj(n),3) = z0
                    array_qual_traj(n,ind_traj(n),4) = tnext(n)/tcyc

                    !!NG              To print the vertical velocity
                    !!NG              IF (n == 1) THEN
                    !!NG              CALL sub_reducmem_shift_or_not_ind(i0(n),j0(n),k0(n),is,js,ks)
                    !!NG              WRITE(0,*) z0, ( ww(is,js,ks,1)   - &
                    !!NG                   ( gk(n) + REAL(k0(n), kind = rprec))   * & 
                    !!NG                   (ww(is,js,ks+1,1) - ww(is,js,ks,1) )) / &
                    !!NG                   (e1t(is,js,1,1)*e2t(is,js,1,1))
                    !!NG              ENDIF

                    IF (key_alltracers) THEN

                       tf(n) = zinter(tt,xt,yt,zt, 1._rprec)
                       sf(n) = zinter(ss,xt,yt,zt, 1._rprec)

                       IF (key_approximatesigma) THEN

                          rf(n) = zinter(rr,xt,yt,zt, 1._rprec)

                       ELSE 

                          rf(n)=sigma(zsigma,sf(n),tf(n))

                       ENDIF

                       array_qual_traj(n,ind_traj(n),5) = tf(n)
                       array_qual_traj(n,ind_traj(n),6) = sf(n)
                       array_qual_traj(n,ind_traj(n),7) = rf(n)

                    ENDIF ! key_alltracers

                 ENDIF
                 !
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: On met a jour la date d'ecriture de la position suivante
                 ! BBL
                 tnext(n) = tnext(n) + tlap
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: Puisqu'une particule est *passee par ce morceau de code* il faut
                 ! BBL          continuer les calculs et on utilise IFL=1
                 ! BBL
                 ifl=1
              ELSE
                 ! BBL __________________________________________________________________________
                 ! BBL -- BBL: Cette particule n'a plus de raison de voir un calcul de positions
                 ! BBL          geographiques
                 ! BBL
                 iflag(n)=0
              ENDIF
              !
           ENDIF
        ENDDO
        !
        ! BBL __________________________________________________________________________
        ! BBL BBL: Si IFL=1, une particule au moins necessite toujours un calcul de
        ! BBL       positions geographiques
        ! BBL
        IF (ifl == 1) GOTO 1111
        !
     ENDIF

     !!----------------------------------------------------------------------------!!
     !!--88888888888888888888---- DO LOOP : ON TRAJECTORIES ---8888888888888888----!!
     !!----------------------------------------------------------------------------!!
     !
     ! swap for space and time index positions
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: Mise a jour des parametres de position spatiale et temporelle
     ! BBL
     DO n=1,ntraj

        IF (igo(n) > 0) THEN
           ! BBL __________________________________________________________________________
           ! BBL - BBL: On ne fait les calculs suivants que pour les particules *actives*
           ! BBL         (etat IGO=1,2,3)
           ! BBL
           gi(n)  = hi(n)
           gj(n)  = hj(n)
           gk(n)  = hk(n)
           gl(n)  = hl(n)
           age(n) = agenew(n)
        ENDIF

     ENDDO

     !!===================================!!
     !!= TESTS TO VERIFY IF GAME IS OVER =!!
     !!===================================!!

     !!----------------------------------------------------------------------------!!
     !!--88888888888888888888---- DO LOOP : ON TRAJECTORIES ---8888888888888888----!!
     !!----------------------------------------------------------------------------!!
     ! BBL __________________________________________________________________________
     ! BBL BBL: IFL controle ici la necessite (IFL=1) ou non (IFL=0) de continuer les
     ! BBL       calculs avec le MEME champ de transport charge en memoire
     ! BBL
     ifl=0

     DO  n=1,ntraj
        ! BBL __________________________________________________________________________
        ! BBL BBL: On *masque* d'abord les particules dans l'attente de la mise a jour
        ! BBL       du champ de transport (IGO=-2)
        ! BBL
        IF (igo(n) == 2) igo(n) = -2
        ! BBL __________________________________________________________________________
        ! BBL BBL: On *neutralise* definitivement les particules ayant atteint le terme
        ! BBL       de leur trajectoire (IGO=-3)
        ! BBL
        IF (igo(n) == 3) igo(n) = -3
        ! BBL __________________________________________________________________________
        ! BBL BBL: Si une particule est toujours active (IGO=1), il faudra continuer
        ! BBL       (IFL=1)
        ! BBL
        IF (igo(n) == 1) ifl = 1
     ENDDO
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: On continue donc avec le meme champ de transport (IFL=1)
     ! BBL
     IF (ifl == 1) THEN
        GOTO 2222 !! GAME IS NOT OVER !!
     ENDIF

     !!----------------------------------------------------------------------------!!
     !!--99999999999999999999---- DO LOOP : ON TRAJECTORIES ---9999999999999999----!!
     !!----------------------------------------------------------------------------!!
     ! BBL __________________________________________________________________________
     ! BBL BBL: IFL controle maintenant la necessite (IFL=1) ou non (IFL=0) de
     ! BBL       continuer les calculs, en chargeant en memoire un nouveau champ de
     ! BBL       transport
     ! BBL
     ifl=0

     DO  n=1,ntraj
        ! BBL __________________________________________________________________________
        ! BBL BBL: Si une particule n'a pas atteint le terme de sa trajectoire
        ! BBL        (IGO NE -3), il faudra continuer (IFL=1)
        ! BBL
        IF (igo(n) /= -3) ifl=1
     ENDDO
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: Pas la peine de continuer...
     ! BBL
     IF (ifl == 0) THEN
        GOTO 3333 !! GAME IS OVER !!
     ENDIF
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: Il faut continuer, et mettre a jour l'indice temporel du champ de
     ! BBL       transport a lire, en prenant garde au sens d'integration (FORWARD ou
     ! BBL        BACKWARD)
     ! BBL      Tant que l'on reste dans le meme cycle temporel, on relance vers le
     ! BBL       label 9998, sinon on enchaine sur le label 9999 qui correspond a
     ! BBL       l'incrementation d'un nouveau cycle temporel
     ! BBL
     IF (TRIM(forback) == 'forward') THEN
        lref = lref + 1
        IF (lref <= lmt) THEN
           GOTO 9998
        ENDIF

     ELSE ! backward
        lref = lref - 1
        IF (lref >= 1) THEN
           GOTO 9998
        ENDIF

     ENDIF

  ENDDO ! DO LOOP: CYCLES (iter=1,maxcycles)

  !!==================================================================================!!
  !!============================ MAIN LOOP IS FINISHED ===============================!!
  !!==================================================================================!!

3333 CONTINUE

  IF (TRIM(mode) == 'quantitative') THEN
     ndone=0

     DO  n = 1, ntraj

        ! BBL ___________________________________________________________________
        ! BBL -- BBL: IGO=3 ou -3 signifie que l'integration de la particule est
        ! BBL          terminee
        ! BBL
        IF (ABS(igo(n)) >= 3) ndone = ndone + 1

     ENDDO

     WRITE(lun_standard,*)
     WRITE(lun_standard,*)'>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
     WRITE(lun_standard,*)'HOW MANY PARTICLES REACH A SECTION: ', &
          INT(100._rprec*REAL(ndone,kind=rprec)/REAL(ntraj,kind=rprec)), '%'
     WRITE(lun_standard,*)'    - number of particules         : ', ntraj
     WRITE(lun_standard,*)'    - done                         : ', ndone
     WRITE(lun_standard,*)'    - rest                         : ', ntraj - ndone
     WRITE(lun_standard,*)'    - cycle                        : ', iter-1
     WRITE(lun_standard,*)'<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'
  ENDIF

  !- Comments -!
  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'=================================='
  WRITE(lun_standard,*)'= MAIN LOOP ON CYCLE >>>>>> EXIT ='
  WRITE(lun_standard,*)'=================================='

  !!----------------------------------------------------------------------------!!
  !!--10-10-10-10-10-10-10---- DO LOOP : ON TRAJECTORIES ---10-10-10-10-10-10---!!
  !!----------------------------------------------------------------------------!!

  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'-----------------------'
  WRITE(lun_standard,*)'= Writing Output Data ='
  WRITE(lun_standard,*)'-----------------------'

  !!===============================================================================!!
  !!============================== QUALITATIVE MODE ===============================!!
  !!===============================================================================!!
  IF (TRIM(mode) == 'qualitative' ) THEN

     !----------------!
     !- Trajectories -!
     !----------------!
     !- Netcdf -!
     !----------!

     IF (key_alltracers) THEN

        CALL sub_save_netcdf_trajectories(                   &
             array_qual_traj(:,:,1), array_qual_traj(:,:,2), &
             array_qual_traj(:,:,3), array_qual_traj(:,:,4), &
             array_qual_traj(:,:,5), array_qual_traj(:,:,6), &
             array_qual_traj(:,:,7)                          )
     ELSE
        CALL sub_save_netcdf_trajectories(                   &
             array_qual_traj(:,:,1), array_qual_traj(:,:,2), &
             array_qual_traj(:,:,3), array_qual_traj(:,:,4)  )
     ENDIF

     !---------!
     !- ASCII -!
     !---------!
     IF (key_ascii_outputs) THEN

        DO n = 1, ntraj
           DO i = 1, ind_traj(n)
              WRITE(lun_traj,7055)n,array_qual_traj(n,i,1:nstat)
           ENDDO
        ENDDO

     ENDIF

7055 FORMAT (i0,3(1x,f0.5),1x,f0.5,1x,3(1x,f0.5))

     !!================================================================================!!
     !!============================== QUANTITATIVE MODE ===============================!!
     !!================================================================================!!
     !
     ! BBL __________________________________________________________________________
     ! BBL BBL: On trouve ensuite la suite logique de l'archivage (mode QUANTITATIVE)
     ! BBL
  ELSE

     !-------------------!
     !- Final positions -!
     !-------------------! 
     !---------!
     !- ASCII -!
     !---------!
     IF (key_ascii_outputs) THEN
        DO  n=1,ntraj

           IF (nfin(n) >= 0) THEN
              WRITE(lun_output,7000)secname(nfin(n)), &
                   n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
           ELSE
              WRITE(lun_output,7000)'meanders (-1)', &
                   n,i0(n),j0(n),k0(n),l0(n),tfi(n),tfj(n),tfk(n),tfl(n)
           ENDIF

           IF (key_alltracers) THEN
              WRITE(lun_fin_pos,7059)hi(n),hj(n),-hk(n),hl(n), &
                   ttr(n),truage(n),nfin(n),tf(n),sf(n),rf(n)
           ELSE
              WRITE(lun_fin_pos,7059)hi(n),hj(n),-hk(n),hl(n), &
                   ttr(n),truage(n),nfin(n)
           ENDIF

        ENDDO

     ENDIF

7000 FORMAT(a14,', #',i0,' en:', 1x,i0,1x,i0,1x,i0,1x,i0,',  init.=',4(1x,f0.2),': ',f0.2)

7059 FORMAT(3(1x,f0.3),1x,f0.3,1x,f0.3,1x,f0.3,1x,i0,3(1x,f0.3))

     !----------------------------!
     !- statistics file in ASCII -!
     !----------------------------!
     WRITE(lun_stats,*)' '

     IF (key_unitm3) THEN
        WRITE(lun_stats,*)'total transport (in m3/s): ', trtot/REAL(lmt,kind=rprec), &
             ' ( x lmt =' ,trtot, ')'
        WRITE(lun_stats,*)'max_transport (in m3/s)  : ', max_transport
     ELSE
        WRITE(lun_stats,*)'total transport (in Sv)  : ', trtot/1.e6_rprec
        WRITE(lun_stats,*)'max_transport (in Sv)    : ', max_transport/1.e6_rprec
     ENDIF
     WRITE(lun_stats,*)'# particles              : '  , ntraj


     IF ((TRIM(bin) == 'nobin').AND.(isn(-1) > 0)) THEN

        !------------------------------------------
        ! we mask uninitialized min. and max. ages
        !------------------------------------------
        DO i=1,nstat
           IF (ABS(s1min(i,-1)) > 1.e10_rprec) s1min(i,-1)=0._rprec
           IF (ABS(s1max(i,-1)) > 1.e10_rprec) s1max(i,-1)=0._rprec
           IF (ABS(s0min(i,-1)) > 1.e10_rprec) s0min(i,-1)=0._rprec
           IF (ABS(s0max(i,-1)) > 1.e10_rprec) s0max(i,-1)=0._rprec
        END DO
        !
        WRITE(lun_stats,*)' '
        WRITE(lun_stats,7057)'initial state ',' ','#',isn(-1)
        WRITE(lun_stats,7157)'stats. for: ',(chpstat(i),i=1,nstat)
        WRITE(lun_stats,7257)' min: ',(s0min(i,-1),i=1,nstat)
        WRITE(lun_stats,7257)' max: ',(s0max(i,-1),i=1,nstat)
        WRITE(lun_stats,7257)'mean: ',(s0x(i,-1)/sn(-1),i=1,nstat)
        WRITE(lun_stats,7257)'std. dev.: ', &
             ( SQRT(ABS((s0x2(i,-1)-s0x(i,-1)*s0x(i,-1)/sn(-1))/(sn(-1)-1._rprec))),i=1,nstat)

     ENDIF

     IF (TRIM(bin) /= 'nobin') THEN
        WRITE(lun_stats,*)' '
        WRITE(lun_stats,*)'no stats for initial state available since'
        WRITE(lun_stats,*)'   automatic positioning was by-passed'
     ENDIF

     !================================
     ! 2D flow projection and z stats
     !================================
     subtime = REAL(lmt, kind=rprec) / REAL(1+lmax-lmin, kind =rprec)

     !NG Performances (sub_quant_store_transp => * r_lmt)
     r_lmt =  1._rprec / REAL(lmt, kind=rprec)

     !!NG    WRITE(lun_xy_zonal)((subtime*uxy(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_merid)((subtime*vxy(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_yz_merid)((subtime*vyz(j,k,1)*r_lmt,j=1,jmt),k=1,kmt)
     !!NG    WRITE(lun_yz_verti)((subtime*wyz(j,k,1)*r_lmt,j=1,jmt),k=1,kmt)
     !!NG    WRITE(lun_xz_zonal)((subtime*uxz(i,k,1)*r_lmt,i=1,imt),k=1,kmt)
     !!NG    WRITE(lun_xz_verti)((subtime*wxz(i,k,1)*r_lmt,i=1,imt),k=1,kmt)
     !!NG
     !!NG    IF (key_alltracers) THEN
     !!NG: seq-roms       WRITE(lun_yz_merid)((subtime*vyr(j,k)*r_lmt,j=1,jmt),k=1,nrclas)
     !!NG: seq-roms       WRITE(lun_yz_verti)((subtime*wyr(j,k)*r_lmt,j=1,jmt),k=1,nrclas)
     !!NG: seq-roms       WRITE(lun_xz_zonal)((subtime*uxr(i,k)*r_lmt,i=1,imt),k=1,nrclas)
     !!NG: seq-roms       WRITE(lun_xz_verti)((subtime*wxr(i,k)*r_lmt,i=1,imt),k=1,nrclas)
     !!NG: #endif
     !!NG    ENDIF
     !!NG
     !!NG    WRITE(lun_xy_stats)((uh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_stats)((vh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_stats)((zuh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_stats)((zvh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_stats)((z2uh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG    WRITE(lun_xy_stats)((z2vh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG
     !!NG    IF (key_alltracers) THEN
     !!NG      WRITE(lun_xy_stats)((tuh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((tvh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((t2uh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((t2vh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((suh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((svh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((s2uh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((s2vh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((ruh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((rvh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((r2uh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG      WRITE(lun_xy_stats)((r2vh(i,j,1)*r_lmt,i=1,imt),j=1,jmt)
     !!NG
     !!NG    ENDIF

     uxy(:,:,:)  = uxy(:,:,:) * subtime * r_lmt
     vxy(:,:,:)  = vxy(:,:,:) * subtime * r_lmt
     vyz(:,:,:)  = vyz(:,:,:) * subtime * r_lmt
     wyz(:,:,:)  = wyz(:,:,:) * subtime * r_lmt
     uxz(:,:,:)  = uxz(:,:,:) * subtime * r_lmt
     wxz(:,:,:)  = wxz(:,:,:) * subtime * r_lmt
     uh(:,:,:)   = uh(:,:,:)            * r_lmt
     vh(:,:,:)   = vh(:,:,:)            * r_lmt
     zuh(:,:,:)  = zuh(:,:,:)           * r_lmt
     zvh(:,:,:)  = zvh(:,:,:)           * r_lmt
     z2uh(:,:,:) = z2uh(:,:,:)          * r_lmt
     z2vh(:,:,:) = z2vh(:,:,:)          * r_lmt

     IF (key_alltracers) THEN
        tuh(:,:,:)  = tuh(:,:,:)         * r_lmt
        tvh(:,:,:)  = tvh(:,:,:)         * r_lmt
        t2uh(:,:,:) = t2uh(:,:,:)        * r_lmt
        t2vh(:,:,:) = t2vh(:,:,:)        * r_lmt
        suh(:,:,:)  = suh(:,:,:)         * r_lmt
        svh(:,:,:)  = svh(:,:,:)         * r_lmt
        s2uh(:,:,:) = s2uh(:,:,:)        * r_lmt
        s2vh(:,:,:) = s2vh(:,:,:)        * r_lmt
        ruh(:,:,:)  = ruh(:,:,:)         * r_lmt
        rvh(:,:,:)  = rvh(:,:,:)         * r_lmt
        r2uh(:,:,:) = r2uh(:,:,:)        * r_lmt
        r2vh(:,:,:) = r2vh(:,:,:)        * r_lmt
     ENDIF

     !!---------------------------!
     !!= Quantitative statistics =!
     !!---------------------------!
     trtot1=0._rprec

     WRITE(lun_stats,*)' '
     nsect2 = nsect
     IF (icrit1 == 1) THEN
        nsect2 = 1 + nsect
     ENDIF

     DO n = 0, nsect2
        IF (key_unitm3) THEN
           WRITE(lun_stats,7357)secname(n),subtime*sectrans(n)*r_lmt,n
        ELSE
           WRITE(lun_stats,7357)secname(n),subtime*sectrans(n)/1.e6_rprec*r_lmt,n
        ENDIF
        trtot1=trtot1+sectrans(n)
     END DO

     IF (key_unitm3) THEN
        WRITE(lun_stats,7357)'total',subtime*trtot*r_lmt
        WRITE(lun_stats,7357)'except mnds', &
             subtime*(trtot-sectrans(0))*r_lmt
        WRITE(lun_stats,7357)'lost',subtime*(trtot-trtot1)*r_lmt
     ELSE
        WRITE(lun_stats,7357)'total',subtime*trtot/1.e6_rprec*r_lmt
        WRITE(lun_stats,7357)'except mnds', &
             subtime*(trtot-sectrans(0))/1.e6_rprec*r_lmt
        WRITE(lun_stats,7357)'lost',subtime*(trtot-trtot1)/1.e6_rprec*r_lmt
     ENDIF

7357 FORMAT(a15,1x,f0.4,1x,i0)

     DO n=0,nsect2
        IF (isn(n) > 0) THEN
           WRITE(lun_stats,*)' '
           WRITE(lun_stats,7057)'final state ',secname(n),' #',isn(n)
           WRITE(lun_stats,7157)'stats. ini: ',(chpstat(i),i=1,nstat)
           WRITE(lun_stats,7257)' min: ',(s0min(i,n),i=1,nstat)
           WRITE(lun_stats,7257)' max: ',(s0max(i,n),i=1,nstat)
           WRITE(lun_stats,7257)'mean: ',(s0x(i,n)/sn(n),i=1,nstat)
           WRITE(lun_stats,7257)'std. dev.: ',( SQRT( &
                ABS((s0x2(i,n)-s0x(i,n)*s0x(i,n)/sn(n))/(sn(n)-1._rprec))),i=1,nstat)
           WRITE(lun_stats,7157)'stats. fin: ',(chpstat(i),i=1,nstat)
           WRITE(lun_stats,7257)' min: ',(s1min(i,n),i=1,nstat)
           WRITE(lun_stats,7257)' max: ',(s1max(i,n),i=1,nstat)
           WRITE(lun_stats,7257)'mean: ',(s1x(i,n)/sn(n),i=1,nstat)
           WRITE(lun_stats,7257)'std. dev.: ',( SQRT( &
                ABS((s1x2(i,n)-s1x(i,n)*s1x(i,n)/sn(n))/(sn(n)-1._rprec))),i=1,nstat)
        ENDIF
     END DO

     !----------------------------------------------------------------------!
     !- Define statistic variables and their attributes in quantative mode -!
     !----------------------------------------------------------------------!
     CALL sub_save_netcdf_data_init_stats()

     !------------------------------------!
     !- Save statistics in Netcdf Format -!
     !------------------------------------!
     Call sub_save_netcdf_stats()

     !----------------------------!
     !- Close Netcdf output file -!
     !----------------------------!
     WRITE(lun_standard,*)''
     WRITE(lun_standard,*)'--------------------------------'
     WRITE(lun_standard,*)'= Close Statistics NetCDF File ='
     WRITE(lun_standard,*)'--------------------------------'
     CALL sub_save_netcdf_data_close(ncid_data_stats)

  ENDIF

7057 FORMAT(a13,a15,a2,i7)
7157 FORMAT(a13,7a10)
7257 FORMAT(a13,7f10.3)


  !=====================================================!
  !- OLD binary storage of initial and final positions -!
  !=====================================================!
  IF (key_ascii_outputs) THEN
     DO n=1,ntraj
        ! BBL __________________________________________________________________________
        ! BBL BBL: Il est temps de sortir sur fichier les positions initiales et
        ! BBL        finales de toutes les particules
        ! BBL
        WRITE(lun_final)hi(n),hj(n),-hk(n),hl(n),ttr(n),agenew(n)
        WRITE(lun_init)tfi(n),tfj(n),tfk(n),tfl(n),ttr(n), tage(n)
     ENDDO
  ENDIF

  !! NG: 15_09_2008
  IF (TRIM(mode) == 'quantitative') THEN 
     !-----------------!
     !- Initial Depth -!
     !-----------------!
     !- Netcdf -!
     !----------!

     CALL sub_save_netcdf_init_depth(init_depth(:))

  ENDIF

  !-------------------!
  !- Final positions -!
  !-------------------! 
  !- Netcdf -!
  !----------!
  !! NG: 15_09_2008
  IF (key_alltracers) THEN
     CALL sub_save_netcdf_final_pos( &
          hi(:),hj(:),-hk(:),hl(:) , &
          agenew(:), nfin(:)       , &
          final_depth(:)           , &
          tf(:),sf(:),rf(:)          )
  ELSE
     CALL sub_save_netcdf_final_pos(                           &
          hi(1:ntraj),hj(1:ntraj),-hk(1:ntraj),hl(1:ntraj)   , &
          agenew(1:ntraj), nfin(1:ntraj),final_depth(1:ntraj)  )
  ENDIF


  !----------------------------!
  !- Close Netcdf output file -!
  !----------------------------!
  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'-------------------------------'
  WRITE(lun_standard,*)'= Close Positions NetCDF File ='
  WRITE(lun_standard,*)'-------------------------------'
  CALL sub_save_netcdf_data_close(ncid_data_pos)

  !=====================!
  !- DEALLOCATE memory -!
  !=====================!

  WRITE(lun_standard,*)''
  WRITE(lun_standard,*)'---------------------'
  WRITE(lun_standard,*)'= Deallocate Memory ='
  WRITE(lun_standard,*)'---------------------'

  IF (ALLOCATED(igo)) THEN
     CALL sub_memory(-size(igo),'i','igo','trajec_seq')
     DEALLOCATE(igo)
  END IF
  IF (ALLOCATED(i0)) THEN
     CALL sub_memory(-size(i0),'i','i0','trajec_seq')
     DEALLOCATE(i0)
  END IF
  IF (ALLOCATED(j0)) THEN
     CALL sub_memory(-size(j0),'i','j0','trajec_seq')
     DEALLOCATE(j0)
  END IF
  IF (ALLOCATED(k0)) THEN
     CALL sub_memory(-size(k0),'i','k0','trajec_seq')
     DEALLOCATE(k0)
  END IF
  IF (ALLOCATED(l0)) THEN
     CALL sub_memory(-size(l0),'i','l0','trajec_seq')
     DEALLOCATE(l0)
  END IF
  IF (ALLOCATED(i1)) THEN
     CALL sub_memory(-size(i1),'i','i1','trajec_seq')
     DEALLOCATE(i1)
  END IF
  IF (ALLOCATED(i2)) THEN
     CALL sub_memory(-size(i2),'i','i2','trajec_seq')
     DEALLOCATE(i2)
  END IF
  IF (ALLOCATED(j1)) THEN
     CALL sub_memory(-size(j1),'i','j1','trajec_seq')
     DEALLOCATE(j1)
  END IF
  IF (ALLOCATED(j2)) THEN
     CALL sub_memory(-size(j2),'i','j2','trajec_seq')
     DEALLOCATE(j2)
  END IF
  IF (ALLOCATED(k1)) THEN
     CALL sub_memory(-size(k1),'i','k1','trajec_seq')
     DEALLOCATE(k1)
  END IF
  IF (ALLOCATED(k2)) THEN
     CALL sub_memory(-size(k2),'i','k2','trajec_seq')
     DEALLOCATE(k2)
  END IF
  IF (ALLOCATED(iflag)) THEN
     CALL sub_memory(-size(iflag),'i','iflag','trajec_seq')
     DEALLOCATE(iflag)
  END IF
  IF (ALLOCATED(iperfl)) THEN
     CALL sub_memory(-size(iperfl),'i','iperfl','trajec_seq')
     DEALLOCATE(iperfl)
  END IF
  IF (ALLOCATED(jperfl)) THEN
     CALL sub_memory(-size(jperfl),'i','jperfl','trajec_seq')
     DEALLOCATE(jperfl)
  END IF
  IF (ALLOCATED(nfin)) THEN
     CALL sub_memory(-size(nfin),'i','nfin','trajec_seq')
     DEALLOCATE(nfin)
  END IF

  IF (ALLOCATED(age)) THEN
     CALL sub_memory(-size(age),'r','age','trajec_seq')
     DEALLOCATE(age)
  END IF
  IF (ALLOCATED(agenew)) THEN
     CALL sub_memory(-size(agenew),'r','agenew','trajec_seq')
     DEALLOCATE(agenew)
  END IF
  IF (ALLOCATED(truage)) THEN
     CALL sub_memory(-size(truage),'r','truage','trajec_seq')
     DEALLOCATE(truage)

  END IF
  IF (ALLOCATED(gi)) THEN
     CALL sub_memory(-size(gi),'r','gi','trajec_seq')
     DEALLOCATE(gi)
  END IF
  IF (ALLOCATED(gj)) THEN
     CALL sub_memory(-size(gj),'r','gj','trajec_seq')
     DEALLOCATE(gj)
  END IF
  IF (ALLOCATED(gk)) THEN
     CALL sub_memory(-size(gk),'r','gk','trajec_seq')
     DEALLOCATE(gk)
  END IF
  IF (ALLOCATED(gl)) THEN
     CALL sub_memory(-size(gl),'r','gl','trajec_seq')
     DEALLOCATE(gl)
  END IF
  IF (ALLOCATED(hi)) THEN
     CALL sub_memory(-size(hi),'r','hi','trajec_seq')
     DEALLOCATE(hi)
  END IF
  IF (ALLOCATED(hj)) THEN
     CALL sub_memory(-size(hj),'r','hj','trajec_seq')
     DEALLOCATE(hj)
  END IF
  IF (ALLOCATED(hk)) THEN
     CALL sub_memory(-size(hk),'r','hk','trajec_seq')
     DEALLOCATE(hk)
  END IF
  IF (ALLOCATED(hl)) THEN
     CALL sub_memory(-size(hl),'r','hl','trajec_seq')
     DEALLOCATE(hl)
  END IF
  IF (ALLOCATED(tnext)) THEN
     CALL sub_memory(-size(tnext),'r','tnext','trajec_seq')
     DEALLOCATE(tnext)
  END IF
  IF (ALLOCATED(tswap)) THEN
     CALL sub_memory(-size(tswap),'r','tswap','trajec_seq')
     DEALLOCATE(tswap)
  END IF
  IF (ALLOCATED(zbuoy)) THEN
     CALL sub_memory(-size(zbuoy),'r','zbuoy','trajec_seq')
     DEALLOCATE(zbuoy)
  END IF
  IF (ALLOCATED(u)) THEN
     CALL sub_memory(-size(u),'r','u','trajec_seq')
     DEALLOCATE(u)
  END IF
  IF (ALLOCATED(v)) THEN
     CALL sub_memory(-size(v),'r','v','trajec_seq')
     DEALLOCATE(v)
  END IF
  IF (ALLOCATED(w)) THEN
     CALL sub_memory(-size(w),'r','w','trajec_seq')
     DEALLOCATE(w)
  END IF
  IF (ALLOCATED(du)) THEN
     CALL sub_memory(-size(du),'r','du','trajec_seq')
     DEALLOCATE(du)
  END IF
  IF (ALLOCATED(dv)) THEN
     CALL sub_memory(-size(dv),'r','dv','trajec_seq')
     DEALLOCATE(dv)
  END IF
  IF (ALLOCATED(dw)) THEN
     CALL sub_memory(-size(dw),'r','dw','trajec_seq')
     DEALLOCATE(dw)
  END IF
  IF (ALLOCATED(tx)) THEN
     CALL sub_memory(-size(tx),'r','tx','trajec_seq')
     DEALLOCATE(tx)
  END IF
  IF (ALLOCATED(ty)) THEN
     CALL sub_memory(-size(ty),'r','ty','trajec_seq')
     DEALLOCATE(ty)
  END IF
  IF (ALLOCATED(tz)) THEN
     CALL sub_memory(-size(tz),'r','tz','trajec_seq')
     DEALLOCATE(tz)
  END IF
  IF (ALLOCATED(tf)) THEN
     CALL sub_memory(-size(tf),'r','tf','trajec_seq')
     DEALLOCATE(tf)
  END IF
  IF (ALLOCATED(sf)) THEN
     CALL sub_memory(-size(sf),'r','sf','trajec_seq')
     DEALLOCATE(sf)
  END IF
  IF (ALLOCATED(rf)) THEN
     CALL sub_memory(-size(rf),'r','rf','trajec_seq')
     DEALLOCATE(rf)
  END IF
  IF (ALLOCATED(tfac)) THEN
     CALL sub_memory(-size(tfac),'r','tfac','trajec_seq')
     DEALLOCATE(tfac)
  END IF
  IF (ALLOCATED(ttt)) THEN
     CALL sub_memory(-size(ttt),'r','ttt','trajec_seq')
     DEALLOCATE(ttt)
  END IF

  !! NG: 15_09_2008
  IF (ALLOCATED(init_depth)) THEN
     CALL sub_memory(-size(init_depth),'r','init_depth','trajec_seq')
     DEALLOCATE(init_depth)
  END IF
  IF (ALLOCATED(final_depth)) THEN
     CALL sub_memory(-size(final_depth),'r','final_depth','trajec_seq')
     DEALLOCATE(final_depth)
  END IF

  IF (ALLOCATED(ti)) THEN
     CALL sub_memory(-size(ti),'r','ti','trajec_seq')
     DEALLOCATE(ti)
  END IF
  IF (ALLOCATED(si)) THEN
     CALL sub_memory(-size(si),'r','si','trajec_seq')
     DEALLOCATE(si)
  END IF
  IF (ALLOCATED(ri)) THEN
     CALL sub_memory(-size(ri),'r','ri','trajec_seq')
     DEALLOCATE(ri)
  END IF

  IF (ALLOCATED(ind_traj)) THEN
     CALL sub_memory(-size(ind_traj),'i','ind_traj','trajec_seq')
     DEALLOCATE(ind_traj)
  END IF
  IF (ALLOCATED(array_qual_traj)) THEN
     CALL sub_memory(-size(array_qual_traj),'r','array_qual_traj','trajec_seq')
     DEALLOCATE(array_qual_traj)
  END IF


  CALL sub_netcdf_dealloc_mem()
  CALL sub_input_data_dealloc_mem()
  CALL sub_coord_dealloc()
  CALL sub_scalef_dealloc()
  IF (key_roms)  CALL sub_h_roms_dealloc()
  CALL sub_transp_dealloc()
  IF (key_alltracers) THEN
     CALL sub_tracer_dealloc()
  ENDIF
  CALL sub_posin_dealloc()
  CALL sub_tmask_dealloc()
  CALL sub_seq_dealloc()

  IF (TRIM(mode) == 'quantitative') THEN ! [QUANT]
     CALL sub_txt_dealloc()
     CALL sub_flags_dealloc()
     CALL sub_stats_dealloc()
     CALL sub_stati_dealloc()
     CALL sub_quant_dealloc()
  ENDIF

END SUBROUTINE trajec_seq