MODULE trajinitsave
   !!======================================================================
   !!                   ***  MODULE trajinitsave  ***
   !! Initialisation and saving of particles position
   !! Read circle or rectangle boundaries to calculate initial position, 
   !! or directly real position in netcdf file
   !! Save in a netcdf file
   !&E       !  01-2012 (M. Huret) separation from traject3d
   !&E
   !&E Description :
   !&E    traj_init3d
   !&E    traj_save3d
   !&E
   !&E ** History :
   !&E       !  2011-10 (T. Odaka, M. Huret, V. Garnier) Introduction of MPI
   !&E       !  2012-01 (M. Huret) Read netcdf
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,  ONLY : imin,imax,jmin,jmax,kmax,ljmin,ljmax,ljmaxp1,  &
                           limaxp1,liminm1,limaxp2,ljminm1,ljmaxp2,limax, &
                           rsh,rlg,riosh,lchain,limin,liminm2,ljminm2,    &
                           rg_valmanq_io,dg_valmanq_io
   USE ionc4,       ONLY : ionc4_createfile_traj, ionc4_createvar_traj,   &
                           ionc4_write_time,ionc4_sync,ionc4_close,       &
                           ionc4_write_trajt,ionc4_read_dimtraj,          &
                           ionc4_read_trajt,ionc4_openr
   USE_MPI toolmpi, ONLY : iproc, nprocs

   IMPLICIT NONE
   PRIVATE

#if defined key_ibm || defined key_trajec3d

   !! * Accessibility
   PUBLIC   :: traj_save3d          ! routine called by traject3d
   PUBLIC   :: traj_init3d          ! routine called by step.F90, ibm_init.F90
   PUBLIC   :: init_patch
   PUBLIC   :: indices_loc2glob

   !! * Shared module variables

   !! * Private variables

 CONTAINS
 
   !!======================================================================

   SUBROUTINE indices_loc2glob(in_start, nb_part, out_start, out_end)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE indices_loc2glob  ***
   !&E
   !&E ** Purpose : .
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_init3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2015-01  (M. Honnorat) IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE_MPI toolmpi, ONLY : exchange_vectcpu_int

   !! * Arguments
   INTEGER, INTENT(in)  :: in_start
   INTEGER, INTENT(in)  :: nb_part
   INTEGER, INTENT(out) :: out_start
   INTEGER, INTENT(out) :: out_end

   !! * Local declarations
#ifdef key_MPI_2D
   INTEGER, DIMENSION(0:nprocs-1)   :: nb_part_percpu
#endif

   !!----------------------------------------------------------------------
   !! * Executable part
 
   out_start = in_start
#ifdef key_MPI_2D
   nb_part_percpu(iproc) = nb_part
   CALL exchange_vectcpu_int(nb_part_percpu)

   IF ( iproc /= 0 ) THEN
      out_start = out_start + SUM(nb_part_percpu(0:iproc-1))
   END IF
#endif
   out_end = out_start-1 + nb_part

   END SUBROUTINE indices_loc2glob

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

   SUBROUTINE init_patch(patch, nb_part)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE init_patch  ***
   !&E
   !&E ** Purpose : .
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_init3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2015-01  (M. Honnorat) IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comtraj,     ONLY : type_patch
   USE_MPI toolmpi, ONLY : ADD_ALL_MPI_INT

   !! * Arguments
   TYPE(type_patch), INTENT(inout)  :: patch
   INTEGER,          INTENT(in)     :: nb_part

   !! * Local declarations
   INTEGER                          :: nb_part_total
   INTEGER                          :: nb_part_alloc

   !!----------------------------------------------------------------------
   !! * Executable part
 
   IF ( nb_part > 0 ) THEN
       nb_part_alloc = patch % nb_part_batch * ( INT(nb_part/patch%nb_part_batch)+1 )
   ELSE
       nb_part_alloc = 0
   ENDIF
   patch % nb_part_alloc = nb_part_alloc

   ! Allocation of particles of each patch
   ALLOCATE( patch % particles(patch % nb_part_alloc) )

   patch % particles(:) = patch % init_particle

   nb_part_total = nb_part
#ifdef key_MPI_2D
   CALL ADD_ALL_MPI_INT(nb_part_total)
#endif
   patch % nb_part_total = nb_part_total
   patch % nb_part_max   = nb_part_total    ! Will be overwritten in ibm_init

   END SUBROUTINE init_patch

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

   SUBROUTINE traj_init3d(xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_init3d  ***
   !&E
   !&E ** Purpose : Read file traject.dat to initialize trajectories variables.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  12-2006 (M. Chiffle, V. Garnier)  Original code
   !&E       !  05-2010 (M. Sourisseau)  Introduction of derivate type
   !&E       !  01-2011 (M. Huret)  Adaptation to siggen
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,  ONLY : pi
   USE comvars2d,   ONLY : iscreenlog,tdeb3d,tfin,h0_g,radearth,lat2d_g,lon2du_g,lat2dv_g,ierrorlog
   USE trajectools, ONLY : h0int,xeint,loc_h0
#if defined key_siggen
   USE trajectools, ONLY : ztosiggen,hc_sigint
#endif
   USE comtraj, ONLY : patch_list_append, patches, type_patch, file_trajec, itypetraj, ndtz
   USE_MPI comtraj, ONLY : init_mpi_type_particle
   USE_MPI toolmpi, ONLY : MASTER, ADD_ALL_MPI_INT

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

   !! * Local declarations
   LOGICAL                                     :: ex,l_posit
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=19)                           :: dateread, tool_sectodat,species
   CHARACTER(LEN=lchain)                       :: rec
   INTEGER                                     :: eof,i,j,k,nn,n,m,mm,kk,l
   INTEGER                                     :: nb_part
   INTEGER                                     :: ierr_mpi
   REAL(KIND=rsh)                              :: dxc,dyc  ! dx and dy at the center of the circle
   REAL(KIND=rsh)                              :: fi,radeg,ray_patch
   REAL(KIND=rsh)                              :: tool_latlon2i,tool_latlon2j,rx,ry,dl, &
                                                  d3,xtemp,ytemp,xe_lag,h0_lag
#if defined key_siggen
   REAL(KIND=rsh)                              :: kint, spos, hc_sig_lag
#endif
   REAL(KIND=rlg)          :: t_traj_beg, t_traj_end, dt_traj

   INTEGER  :: imin_patch, imax_patch,   &
               jmin_patch, jmax_patch,   &
               kmin_patch, kmax_patch
   INTEGER  :: istep_patch, jstep_patch, kstep_patch
   INTEGER  :: nb_patch, nb_part_nc, nb_part_intro, nbpart_patch
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)    :: depth_nc
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)    :: lon_nc,lat_nc
   REAL(KIND=rlg)                               :: tool_datosec,gmin,gmax,phimin,phimax
   TYPE(type_patch), POINTER                    :: new_patch, patch
!  parameters for interpolation at particle location
   REAL(KIND=rsh)                       :: px,py
   INTEGER                              :: igg,idd,jbb,jhh,hlb,hlt,hrb,hrt
   INTEGER  :: idx_s, idx_e

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

   PRINT_DBG*, 'enter TRAJ_INIT3D routine'

   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) '**************************************************'
      WRITE(iscreenlog,*) '***************** TRAJ_INIT3D.F90 ****************'
      WRITE(iscreenlog,*) '**************************************************'
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) 'fichier definissant les caracteristiques des trajectoires : ',trim(file_trajec)
   ENDIF_MPI
   ! Make sure the type of trajectory is correct
   INQUIRE(file=file_trajec, exist=ex)
   IF (.NOT. ex) THEN
     print*, "Trajectory file '" // trim(file_trajec) // "' does not exist."
     print*, "Check in 'paraspec.txt' or 'paraibm.txt' if you use key_ibm."
     print*, "Simulation stopped."
     CALL_MPI MPI_FINALIZE(ierr_mpi)
     stop
   END IF
   IF (itypetraj/=1 .AND. itypetraj/=2 .AND. itypetraj/=3) THEN
      print*, "Type of trajectory is not defined correctly."
      print*, "Must be 1 (circle patch), 2 (rectangle patch) or 3 (Netcdf)"
      print*, "Simulation stopped."
      CALL_MPI MPI_FINALIZE(ierr_mpi)
      stop
   END IF

   OPEN(49, file=file_trajec, form='formatted')
   comment='debut'
   DO WHILE ( comment /= '*****' )
      READ(49,'(a)',iostat=eof) comment
   END DO

   nb_patch = 0
   DO WHILE ( eof == 0 )
      READ(49,'(a)',iostat=eof) comment
      IF ( (comment == '*****') .AND. (eof == 0) ) nb_patch = nb_patch+1
   END DO
   IF_MPI (MASTER) THEN
      WRITE(iscreenlog,*) 'Number of patches for estimation of trajectories = ',nb_patch
   ENDIF_MPI

   ! Rewind
   REWIND(49)
   comment='debut'
   DO WHILE (comment /= '*****')
      READ(49,'(a)',iostat=eof) comment
   END DO

   DO n = 1,nb_patch

      ! Create new patch data structure
      new_patch => patch_list_append(patches)

      ! Read name of patch
      READ(49,'(a)',iostat=eof) ! name_patch

      ! Read starting date of trajectory
      READ(49,'(a)',iostat=eof) dateread
      t_traj_beg = tool_datosec(dateread)
      IF( t_traj_beg < tdeb3d .OR. t_traj_beg > tfin ) THEN
         IF_MPI (MASTER) THEN
            write(iscreenlog,*) ' '
            write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
            write(iscreenlog,*) 'WARNING : Starting date ',trim(dateread),' is uncorrect (before 3D simulation starting).'
            write(iscreenlog,*) 'WARNING : Consequently trajectory starts with departure of 3D run on ', tool_sectodat(tdeb3d)
         ENDIF_MPI
         t_traj_beg = tdeb3d
      ENDIF

      ! Read ending date of trajectory
      READ(49,'(a)',iostat=eof) dateread
      t_traj_end = tool_datosec(dateread)
      IF( t_traj_end < tdeb3d .OR. t_traj_end > tfin ) THEN
         IF_MPI (MASTER) THEN
            WRITE(iscreenlog,*) ' '
            WRITE(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
            WRITE(iscreenlog,*) 'WARNING : Ending date ',trim(dateread),' is uncorrect (after 3D simulation ending).'
            WRITE(iscreenlog,*) 'WARNING : Consequently trajectory ends with run on ', tool_sectodat(tfin)
         ENDIF_MPI
         t_traj_end = tfin
      ENDIF

      ! Read time step for outputs
      READ(49,*,iostat=eof) dt_traj

      new_patch % t_beg   = t_traj_beg
      new_patch % t_end   = t_traj_end
      new_patch % t_save  = t_traj_beg
      new_patch % dt_save = dt_traj

      IF_MPI (MASTER) THEN
         WRITE(iscreenlog,*) 'PATCH NUMBER : ',n, new_line(''),                                      &
                             '   trajectory from ' // trim(tool_sectodat(t_traj_beg)), new_line(''), &
                             '                to ' // trim(tool_sectodat(t_traj_end)), new_line(''), &
                             '   with a ', dt_traj, 'hours time step.'
      ENDIF_MPI

      IF ( itypetraj == 1 ) THEN

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! Initialize circle patches
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! Spatial extension of initial patch in indexes or long. lati. ?
         READ(49,*,iostat=eof) l_posit

         ! Read location of initial patch
         IF (l_posit) THEN
            READ(49,*,iostat=eof) imin_patch, jmin_patch
         ELSE
            READ(49,*,iostat=eof) gmin,phimin
            imin_patch = NINT(tool_latlon2i(gmin,phimin))
            jmin_patch = NINT(tool_latlon2j(gmin,phimin))
            IF ( imin_patch < imin .OR. imin_patch > imax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               imin_patch = MIN(MAX(imin_patch,imin),imax)
               print*,' its longitude is :', imin_patch
            END IF
            IF ( jmin_patch < jmin .OR. jmin_patch > jmax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               jmin_patch = MIN(MAX(jmin_patch,jmin),jmax)
               print*,' its latitude is :', jmin_patch
            END IF
         END IF

         ! Read ray of initial patch
         READ(49,*,iostat=eof) ray_patch

         ! Read nbpart_patch
         READ(49,*,iostat=eof) nbpart_patch

         ! Read depth of initial patch (read as immersion in meters)
         READ(49,*,iostat=eof) kmin_patch,kmax_patch

         ! Read resolution depth of initial patch
         READ(49,*,iostat=eof) kstep_patch

         ! Correct if depth not set as imersion
         kstep_patch = ABS(kstep_patch)
         kmin_patch  = ABS(kmin_patch)
         kmax_patch  = ABS(kmax_patch)

         ! Number of particles set at each initial position
         READ(49,*,iostat=eof) nb_part_intro

         ! Type of vertical behavior (integer):
             ! itypevert = 0 if constant depth
             ! itypevert < 0 if no random walk (vertical advection only)
             ! itypevert > 0 if random walk (advection + diffusion)
             ! abs(itypevert) = 1 if no vertical swimming 
             ! abs(itypevert) > 1 if vertical swimming (larval behavior):
             !                     = 2 for nycthemeral migration
             !                     = 3 for ontogenic migration (sakina), ...
         READ(49,*,iostat=eof) new_patch % init_particle % itypevert

         ! Read output file
         READ(49,'(a)',iostat=eof) rec
         kk = index(rec,',|')
         IF (kk > 0 ) THEN
            new_patch % file_out = rec(1:kk-1)
         ELSE
            new_patch % file_out = rec
         END IF
         READ(49,*,iostat=eof)
         ! == End of file reading

         ! Estimate number of particle inside the circle patch
         ! mean spatial size of the cell 
         radeg = pi/180.0_rsh
         fi = REAL(lat2d_g(imin_patch,jmin_patch),rsh)
         IF ( imin_patch == imin ) THEN
            dxc = radearth*COS(fi*radeg)*REAL(lon2du_g(imin_patch+1,jmin_patch)  &
                                             -lon2du_g(imin_patch,  jmin_patch),rsh)*radeg
         ELSE
            dxc = radearth*COS(fi*radeg)*REAL(lon2du_g(imin_patch,  jmin_patch)  &
                                             -lon2du_g(imin_patch-1,jmin_patch),rsh)*radeg
         END IF
         IF (jmin_patch==jmin) THEN
            dyc = radearth*REAL(lat2dv_g(imin_patch,jmin_patch+1)    &
                               -lat2dv_g(imin_patch,jmin_patch),rsh)*radeg
         ELSE
            dyc = radearth*REAL(lat2dv_g(imin_patch,jmin_patch)      &
                               -lat2dv_g(imin_patch,jmin_patch-1),rsh)*radeg
         END IF

         nn = INT(sqrt(4/pi*nbpart_patch))+1
         dl = 2*ray_patch/nn
         IF ( h0_g(imin_patch,jmin_patch) < kmin_patch ) THEN
            print*,'Patch on land or too deep'
            print*,'Check patch number', new_patch % id, 'in file "',trim(file_trajec),'"'
            print*,'Simulation stopped.'
            stop
         END IF

         nb_part = 0
         DO j=1,nn
         DO i=1,nn
            rx = -ray_patch+REAL(i-1)*dl
            ry = -ray_patch+REAL(j-1)*dl
            IF ( sqrt(rx*rx+ry*ry) <= ray_patch ) THEN
               xtemp = imin_patch+rx/dxc
               ytemp = jmin_patch+ry/dyc
               IF_MPI ( limin<=NINT(xtemp) .AND. NINT(xtemp)<=limax .AND. ljmin<=NINT(ytemp) .AND. NINT(ytemp)<=ljmax ) THEN
                  DO k=kmin_patch,kmax_patch,kstep_patch
                     IF ( h0_g(NINT(xtemp),NINT(ytemp)) > k ) THEN
                        CALL loc_h0(xtemp,ytemp,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        h0_lag = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        d3 = h0_lag+xe_lag
                        IF (d3 > k) THEN
                           nb_part = nb_part+nb_part_intro
                        END IF
                     END IF
                  END DO
               ENDIF_MPI
            END IF
         END DO
         END DO

         CALL init_patch(new_patch, nb_part)
         CALL indices_loc2glob(0, nb_part, idx_s, idx_e)
         m = 0
         DO j=1,nn
         DO i=1,nn
            rx = -ray_patch+REAL(i-1)*dl
            ry = -ray_patch+REAL(j-1)*dl
            IF ( sqrt(rx*rx+ry*ry) <= ray_patch ) THEN
               xtemp = imin_patch+rx/dxc
               ytemp = jmin_patch+ry/dyc
               IF_MPI ( limin<=NINT(xtemp) .AND. NINT(xtemp)<=limax .AND. ljmin<=NINT(ytemp) .AND. NINT(ytemp)<=ljmax ) THEN
                  DO k=kmin_patch,kmax_patch,kstep_patch
                     IF ( h0_g(NINT(xtemp),NINT(ytemp)) > k ) THEN
                        mm = m + 1
                        m  = m + nb_part_intro
                        new_patch % particles(mm:m) % xpos = imin_patch + rx/dxc ! i index of initial location
                        new_patch % particles(mm:m) % ypos = jmin_patch + ry/dyc ! j index of initial location
                        ! total depth at particle s location
                        CALL loc_h0( new_patch % particles(m) % xpos,  &
                                     new_patch % particles(m) % ypos,  &
                                     px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        h0_lag = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                        d3 = h0_lag+xe_lag
                        IF (d3 > k) THEN
#ifdef key_MPI_2D
                           new_patch % particles(mm:m) % active = .TRUE.
#endif 
                           new_patch % particles(mm:m) % zpos = k      ! immersion depth
#ifdef key_siggen        
                           kint = -float(k)+xe_lag ! switch from immersion to real z
                           hc_sig_lag = hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                           CALL ztosiggen(kint,spos,xe_lag,h0_lag,hc_sig_lag)
                           new_patch % particles(mm:m) % spos = spos
                           new_patch % particles(mm:m) % hc   = hc_sig_lag
#else
                           new_patch % particles(mm:m) % spos = -(k/d3)
#endif
                           new_patch % particles(mm:m) % d3 = d3
                           new_patch % particles(mm:m) % h0 = h0_lag
                           new_patch % particles(mm:m) % xe = xe_lag
                           DO l=0,nb_part_intro-1
                              new_patch % particles(mm+l) % num = idx_s+mm+l
                           ENDDO
                        ELSE
                           m = m-nb_part_intro
                        END IF
                     END IF
                  END DO
               ENDIF_MPI
            END IF
         END DO
         END DO

      ELSEIF ( itypetraj == 2 ) THEN

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! Initialize rectangular patches
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! Spatial extension of initial patch in indexes or long. lati. ?
         READ(49,*,iostat=eof) l_posit

         ! Read spatial extension of initial patch
         IF (l_posit) THEN
            READ(49,*,iostat=eof) imin_patch, imax_patch, jmin_patch, jmax_patch
         ELSE
            READ(49,*,iostat=eof) gmin,gmax,phimin,phimax
            imin_patch = NINT(tool_latlon2i(gmin,phimin))
            imax_patch = NINT(tool_latlon2i(gmax,phimax))
            jmin_patch = NINT(tool_latlon2j(gmin,phimin))
            jmax_patch = NINT(tool_latlon2j(gmax,phimax))
            IF ( imin_patch < imin .OR. imin_patch > imax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               imin_patch = MIN(MAX(imin_patch,imin),imax)
               print*,' its western longitude is :', imin_patch
            END IF
            IF ( jmin_patch < jmin .OR. jmin_patch > jmax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               jmin_patch = MIN(MAX(jmin_patch,jmin),jmax)
               print*,' its southern latitude is :', jmin_patch
            END IF
            IF ( imax_patch < imin .OR. imax_patch > imax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               imax_patch=MIN(MAX(imax_patch,imin),imax)
               print*,' its eastern longitude is :', imax_patch
            END IF
            IF ( jmax_patch < jmin .OR. jmax_patch > jmax ) THEN
               print*,'LOCATION OUT OF THE DOMAIN'
               print*,'Have a look at file ',trim(file_trajec),' patch number',new_patch % id
               jmax_patch=MIN(MAX(jmax_patch,jmin),jmax)
               print*,' its northern latitude is :', jmax_patch
            END IF
         END IF

         ! Read spatial dispersion of particles inside initial patch
         READ(49,*,iostat=eof) istep_patch, jstep_patch

         ! Read depth of initial patch (read as immersion in meters)
         READ(49,*,iostat=eof) kmin_patch,kmax_patch

         ! Read resolution depth of initial patch
         READ(49,*,iostat=eof) kstep_patch

         ! Number of particles set at each initial position
         READ(49,*,iostat=eof) nb_part_intro

         ! Type of vertical behavior (integer):
             ! itypevert = 0 if constant depth
             ! itypevert < 0 if no random walk (vertical advection only)
             ! itypevert > 0 if random walk (advection + diffusion)
             ! abs(itypevert) = 1 if no vertical swimming 
             ! abs(itypevert) > 1 if vertical swimming (larval behavior):
             !                     = 2 for nycthemeral migration
             !                     = 3 for ontogenic migration (sakina), ...
         READ(49,*,iostat=eof) new_patch % init_particle % itypevert

         ! Read output file
         READ(49,'(a)',iostat=eof) rec
         kk = index(rec,',|')
         IF (kk > 0 ) THEN
            new_patch % file_out = rec(1:kk-1)
         ELSE
            new_patch % file_out = rec
         END IF
         READ(49,*,iostat=eof)
         ! == End of file reading

         ! Estimate number of particle inside the rectangular patch
         nb_part = 0
         DO j=MAX0(ljmin,jmin_patch),MIN0(ljmax,jmax_patch),jstep_patch
         DO i=MAX0(limin,imin_patch),MIN0(limax,imax_patch),istep_patch
         DO k=kmin_patch,kmax_patch,kstep_patch

            IF (h0_g(i,j) > k) THEN
               xtemp=FLOAT(i)
               ytemp=FLOAT(j)
               CALL loc_h0(xtemp,ytemp,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               h0_lag = h0int(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               d3 = h0_lag+xe_lag
               IF (d3 > k) THEN
                  nb_part = nb_part + nb_part_intro
               END IF
            END IF
         END DO
         END DO
         END DO

         CALL init_patch(new_patch, nb_part)
         CALL indices_loc2glob(0, nb_part, idx_s, idx_e)
         m = 0
         DO j=MAX0(ljmin,jmin_patch),MIN0(ljmax,jmax_patch),jstep_patch
         DO i=MAX0(limin,imin_patch),MIN0(limax,imax_patch),istep_patch
         DO k=kmin_patch,kmax_patch,kstep_patch
            IF (h0_g(i,j) > k) THEN
               mm = m + 1
               m  = m + nb_part_intro
               new_patch % particles(mm:m) % xpos = FLOAT(i)   ! i index of initial location    
               new_patch % particles(mm:m) % ypos = FLOAT(j)   ! j index                                        
               ! total depth at particle s location
               CALL loc_h0( new_patch % particles(m) % xpos,  &
                            new_patch % particles(m) % ypos,  &
                            px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               h0_lag = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
               d3 = h0_lag+xe_lag
               IF (d3 > k) THEN
                  new_patch % particles(mm:m) % active = .TRUE.
                  new_patch % particles(mm:m) % zpos = k      ! immersion depth
#ifdef key_siggen        
                  kint = -float(k)+xe_lag ! switch from immersion to real z
                  hc_sig_lag = hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  CALL ztosiggen(kint,spos,xe_lag,h0_lag,hc_sig_lag)
                  new_patch % particles(mm:m) % spos = spos
                  new_patch % particles(mm:m) % hc   = hc_sig_lag
#else
                  new_patch % particles(mm:m) % spos = -(k/d3)
#endif
                  new_patch % particles(mm:m) % d3 = d3
                  new_patch % particles(mm:m) % h0 = h0_lag
                  new_patch % particles(mm:m) % xe = xe_lag
                  DO l=0,nb_part_intro-1
                     new_patch % particles(mm+l) % num = idx_s+mm+l
                  ENDDO
               ELSE
                  m = m-nb_part_intro
               END IF
            END IF
         END DO
         END DO
         END DO

      ELSEIF ( itypetraj == 3 ) THEN

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! Initialize netcdf patches
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! Read name of input netcdf file
         READ(49,'(a)',iostat=eof) rec
         kk = index(rec,',|')
         IF (kk > 0 ) THEN
            new_patch % file_inp = rec(1:kk-1)
         ELSE
            new_patch % file_inp = rec
         END IF

         ! Number of particles set at each initial position
         READ(49,*,iostat=eof) nb_part_intro

         ! Type of vertical behavior (integer):
             ! itypevert = 0 if constant depth
             ! itypevert < 0 if no random walk (vertical advection only)
             ! itypevert > 0 if random walk (advection + diffusion)
             ! abs(itypevert) = 1 if no vertical swimming 
             ! abs(itypevert) > 1 if vertical swimming (larval behavior):
             !                     = 2 for nycthemeral migration
             !                     = 3 for ontogenic migration (sakina), ...
         READ(49,*,iostat=eof) new_patch % init_particle % itypevert

         READ(49,'(a)',iostat=eof) species
#ifdef key_ibm_species
         new_patch % species = species
#endif
         ! Read output file
         READ(49,'(a)',iostat=eof) rec
         kk = index(rec,',|')
         IF (kk > 0 ) THEN
            new_patch % file_out = rec(1:kk-1)
         ELSE
            new_patch % file_out = rec
         END IF
         READ(49,*,iostat=eof)
         ! == End of file reading

         ! Estimate/correct number of particle inside the patch

         ! Open input file and read number of particles
         CALL ionc4_openr(trim(new_patch % file_inp),.false.)
         CALL ionc4_read_dimtraj(trim(new_patch % file_inp), nb_part_nc)

         ALLOCATE( lon_nc(nb_part_nc), lat_nc(nb_part_nc), depth_nc(nb_part_nc) )

         ! read lat,lon,depth of particles in file
         CALL ionc4_read_trajt(trim(new_patch%file_inp), "longitude",lon_nc,  1,nb_part_nc,1)
         CALL ionc4_read_trajt(trim(new_patch%file_inp), "latitude", lat_nc,  1,nb_part_nc,1)
         CALL ionc4_read_trajt(trim(new_patch%file_inp), "DEPTH",    depth_nc,1,nb_part_nc,1)

         nb_part = 0
         DO nn = 1,nb_part_nc

            IF ( depth_nc(nn) < 0.0_rsh ) THEN
               WRITE(ierrorlog,*) 'Function INIT_TRAJ : depth of particle has to be > 0'
               WRITE(ierrorlog,*) 'Check the netcdf traj file :', new_patch % file_inp
               WRITE(ierrorlog,*) 'The simulation is stopped'
               CALL_MPI MPI_FINALIZE(ierr_mpi)
               STOP
            END IF

            i = NINT(tool_latlon2i(lon_nc(nn),lat_nc(nn)))
            j = NINT(tool_latlon2j(lon_nc(nn),lat_nc(nn)))
            IF ( i >= limin .and. i<=limax .and. j>=ljmin .and. j<=ljmax ) THEN
               IF ( h0_g(i,j) > depth_nc(nn) ) THEN
                  xtemp = tool_latlon2i(lon_nc(nn),lat_nc(nn))
                  ytemp = tool_latlon2j(lon_nc(nn),lat_nc(nn))
                  CALL loc_h0(xtemp,ytemp,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  h0_lag = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  d3 = h0_lag+xe_lag
            !      IF (d3 > depth_nc(nn)) THEN   ! patch tempo pour restart (martin)
                     nb_part = nb_part + nb_part_intro
            !      END IF
               END IF
            END IF
         END DO

         CALL init_patch(new_patch, nb_part)
         CALL indices_loc2glob(0, nb_part, idx_s, idx_e)
         m = 0
         DO nn = 1,nb_part_nc
            i = NINT(tool_latlon2i(lon_nc(nn),lat_nc(nn)))
            j = NINT(tool_latlon2j(lon_nc(nn),lat_nc(nn)))
            IF ( i >= limin .and. i<=limax .and. j>=ljmin .and. j<=ljmax ) THEN
               IF ( h0_g(i,j) > depth_nc(nn) ) THEN
                  mm = m + 1
                  m  = m + nb_part_intro
                  new_patch % particles(mm:m) % xpos = tool_latlon2i(lon_nc(nn),lat_nc(nn))
                  new_patch % particles(mm:m) % ypos = tool_latlon2j(lon_nc(nn),lat_nc(nn))                                   

                  ! total depth at particle s location
                  CALL loc_h0( new_patch % particles(m) % xpos,  &
                               new_patch % particles(m) % ypos,  &
                               px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  xe_lag = xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  h0_lag = h0int(   px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                  d3 = h0_lag+xe_lag
                  IF ( d3 < depth_nc(nn) ) depth_nc(nn)=d3-0.1_rsh ! patch tempo pour restart (martin)
                  IF ( d3 > depth_nc(nn) ) THEN
                     new_patch % particles(mm:m) % active = .TRUE.
                     new_patch % particles(mm:m) % zpos = depth_nc(nn)   ! immersion depth
#ifdef key_siggen        
                     kint = -depth_nc(nn)+xe_lag ! switch from immersion to real z
                     hc_sig_lag = hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
                     CALL ztosiggen(kint,spos,xe_lag,h0_lag,hc_sig_lag)
                     new_patch % particles(mm:m) % spos = spos
                     new_patch % particles(mm:m) % hc   = hc_sig_lag
#else
                     new_patch % particles(mm:m) % spos = -(depth_nc(nn)/d3)
#endif		
                     new_patch % particles(mm:m) % d3 = d3
                     new_patch % particles(mm:m) % h0 = h0_lag
                     new_patch % particles(mm:m) % xe = xe_lag
                     DO l=0,nb_part_intro-1
                        new_patch % particles(mm+l) % num = idx_s+mm+l
                     ENDDO
                  ELSE
                     m = m-nb_part_intro
                  END IF
               END IF
            END IF
         END DO
         DEALLOCATE(lon_nc,lat_nc,depth_nc)

         ! close netcdf file
         CALL ionc4_close(new_patch % file_inp)

      END IF  ! end test on itypetraj

   END DO  ! loop on patches

   CLOSE(49)

   IF_MPI (MASTER) THEN
      INQUIRE(file=file_trajec,exist=ex)
      IF (ex) THEN
         patch => patches % first
         write(iscreenlog,*) 'VERTICAL COMPONENT OF THE DISPLACEMENT:  '
         DO n = 1,patches%nb
            write(iscreenlog,*) ' ndtz : ', ndtz
            IF     ( patch % init_particle % itypevert == 0 ) THEN
               WRITE(iscreenlog,*) ' Patch number',n,': Trajectories at constant depth ' 
            ELSEIF ( patch % init_particle % itypevert < 0  ) THEN
               WRITE(iscreenlog,*) ' Patch number',n,': Advection only (no random walk) '
            ELSE
               WRITE(iscreenlog,*) ' Patch number',n,': With random walk (Advection + Diffusion) '
            ENDIF
            patch => patch % next
         END DO
      ENDIF
   ENDIF_MPI

#ifdef key_MPI_2D
   CALL init_mpi_type_particle
#endif

#ifndef key_ibm 
   ! save initialization
   CALL traj_save3d
#endif

   PRINT_DBG*, 'exit TRAJ_INIT3D routine'

   END SUBROUTINE traj_init3d

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

   SUBROUTINE traj_save3d
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_save3d  ***
   !&E
   !&E ** Purpose : Save trajectories
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2006-12 (M. Chiffle, V. Garnier)  Original code
   !&E       !  2010-05 (M. Sourisseau, M. Huret) Coupling with environment variables
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,   ONLY : t
   USE comtraj,     ONLY : patches, type_patch, type_particle

   USE_MPI toolmpi, ONLY : ADD_ALL_MPI_INT

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=lchain)    :: file_out
   LOGICAL                                     :: l_out_nc4par
   INTEGER                                     :: n,m,num1,num2,nb_part,nb_part_nc,p
   REAL(kind=rsh)                              :: fillval
   REAL(KIND=rlg)                              :: tool_ind2lat,tool_ind2lon
   INTEGER,        ALLOCATABLE, DIMENSION(:)   :: num_out
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)   :: lat_out,lon_out
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: xpos_out,ypos_out,flag_out,spos_out,zpos_out,h0pos_out
   TYPE(type_patch),    POINTER     :: patch
   TYPE(type_particle), POINTER     :: particle
   INTEGER  :: idx_s, idx_e

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

   PRINT_DBG*, 'enter TRAJ_SAVE3D routine'

   IF(rsh==8) THEN 
      fillval = dg_valmanq_io
   ELSE
      fillval = rg_valmanq_io
   ENDIF

#ifdef key_MPI_2D
   l_out_nc4par = .true.
#else
   l_out_nc4par = .false.
#endif
   patch => patches % first

   DO n = 1,patches % nb

      IF ( (t < patch % t_save) .OR. (t > patch % t_end) ) THEN
         patch => patch % next
         CYCLE
      END IF

      file_out = trim(patch % file_out)

      IF ( .NOT. patch % file_out_init ) THEN
         ! Create output file
         patch % file_out_init = .TRUE.
         nb_part_nc = patch % nb_part_total
         CALL ionc4_createfile_traj(file_out, nb_part_nc, 0, 0, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "latitude", "degrees_north","latitude",                     &
                                             fill_value=dg_valmanq_io, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "longitude","degrees_east", "longitude",                    &
                                             fill_value=dg_valmanq_io, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DEPTH","m","depth",fill_value=-fillval,l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "H0","m","h0",fill_value=-fillval,l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "NUM","","number of the particle",fill_value=0,l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "flag","nbr","flag",fill_value=fillval,l_out_nc4par=l_out_nc4par)

      END IF    ! (.NOT. patch % file_out_init)

      nb_part = patch % nb_part_alloc
      ALLOCATE( lat_out(nb_part),  lon_out(nb_part) )
      ALLOCATE( xpos_out(nb_part), ypos_out(nb_part) )
      ALLOCATE( spos_out(nb_part), zpos_out(nb_part) )
      ALLOCATE( num_out(nb_part),  h0pos_out(nb_part) )
      ALLOCATE( flag_out(nb_part) )

      lat_out(:)  = dg_valmanq_io ; lon_out(:) = dg_valmanq_io
      xpos_out(:) = fillval ; ypos_out(:) = fillval ; spos_out(:) = fillval ; zpos_out(:) = -fillval
      h0pos_out(:)=-fillval ; flag_out(:) = fillval ; num_out(:)  = 0

      p = 0
      DO m = 1,nb_part
         particle => patch % particles(m)
         IF ( .NOT. particle % active ) CYCLE
         p = p+1
         xpos_out(p)  = particle % xpos
         ypos_out(p)  = particle % ypos
         lat_out(p)   = tool_ind2lat(particle % xpos, particle % ypos)
         lon_out(p)   = tool_ind2lon(particle % xpos, particle % ypos)
         spos_out(p)  = particle % spos
         zpos_out(p)  = particle % zpos ! for output as immersion
         h0pos_out(p) = particle % d3
         flag_out(p)  = particle % flag
         num_out(p)   = particle % num
      END DO

      CALL ionc4_write_time(file_out,0,t)

      nb_part = p
      CALL indices_loc2glob(1, nb_part, idx_s, idx_e)

#ifdef key_MPI_2D
      num1 = idx_s
      num2 = idx_e
#else
      num1 = 1
      num2 = nb_part
#endif
      CALL ionc4_write_trajt(file_out, 'latitude',  lat_out(1:nb_part),num1,num2,0,dg_valmanq_io)
      CALL ionc4_write_trajt(file_out, 'longitude', lon_out(1:nb_part),num1,num2,0,dg_valmanq_io)
      CALL ionc4_write_trajt(file_out, 'DEPTH',    zpos_out(1:nb_part),num1,num2,0,-fillval)
      CALL ionc4_write_trajt(file_out, 'H0',      h0pos_out(1:nb_part),num1,num2,0,-fillval)
      CALL ionc4_write_trajt(file_out, 'NUM',       num_out(1:nb_part),num1,num2,0,0)
      CALL ionc4_write_trajt(file_out, 'flag',     flag_out(1:nb_part),num1,num2,0,fillval)

      ! To write the data on the disk and not loose data in case of run crash
      CALL ionc4_sync(file_out)
      DEALLOCATE(xpos_out, ypos_out, zpos_out, spos_out)
      DEALLOCATE(lat_out, lon_out, h0pos_out, flag_out, num_out)

      patch % t_save = patch % t_save + patch % dt_save*3600.0_rlg
      patch => patch % next
      
   END DO

   PRINT_DBG*, 'exit TRAJ_SAVE3D routine'

   END SUBROUTINE traj_save3d

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

#endif
END MODULE
