 MODULE traject2d
   !!======================================================================
   !!                   ***  MODULE traj2d  ***
   !! Trajectories dynamics:  time evolution of a patch in 3D
   !!
   !!   History :
   !!       !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !!
   !!======================================================================
   !! * Modules used
   USE parameters, ONLY       : rsh,riosh,rlg,imin,imax,jmin,jmax,kmax,ljmin,ljmax,  &
                                lchain,valmanq,ngridmax,nbtachmax,nbpartmax,         &
                                liminm2,limaxp2,ljminm2,ljmaxp2
   USE ionc4, ONLY           : ionc4_createvar,ionc4_vatt_missing,ionc4_write_xy,    &
                               ionc4_createfile_traj,ionc4_torigine,ionc4_write_time,&
                               ionc4_write_trajt

   USE comtraj,    ONLY       : file_trajec

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility

   !! * Shared module variables
#ifdef key_agrif
   real(kind=rlg),dimension(0:ngridmax-1,7), PUBLIC  :: data_grid
   integer,dimension(nbpartmax+1,nbtachmax+1), PUBLIC :: id_part
#endif

#ifdef key_trajec2d
   !! * Accessibility
   PUBLIC traj_2d               ! routine called by stepu.F90 stepv.F90
   PUBLIC traj_init2d           ! routine called by step.F90

   !! * Private variables
                   ! name of trajectories output files
   CHARACTER(LEN=lchain), ALLOCATABLE, DIMENSION(:)  :: filetraj2dout
                   ! number of particles inside each patch
   INTEGER, ALLOCATABLE, DIMENSION(:)           :: nb_part
                   ! starting and endind date of trajectories
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)    :: tbegtraj,tendtraj,tsavtraj,dtsave_traj
                   ! i,j,k indexes of particle location
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:)  :: xpos,ypos
                   ! longitude and latitude of particle location
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:,:)  :: lon_pos,lat_pos
                   ! number of patches
   INTEGER                                      :: nb_patch
                   ! maximum number of particles
   INTEGER                                      :: nbmax_part
                   ! trajectories type
   INTEGER                                      :: itypetraj

                   ! estimate of residual  currents
                   ! M2 period (secondes)
   REAL(KIND=rlg), PARAMETER                    ::  m2per=44712.0_rlg
                   ! nb of hours for residual estimate
   INTEGER,        PARAMETER                    ::  nres=13
                   ! i,j,k indexes of particle location
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:)  :: xposres,yposres,xmoyres,ymoyres
                   ! residual currents
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:)    :: ures,vres,vitres
   INTEGER,        ALLOCATABLE, DIMENSION(:)      :: iswitj
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)      :: tbegres,tprectraj
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:)    :: nomb
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:)  :: uit,vit

 CONTAINS

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

   SUBROUTINE traj_init2d

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_init2d  ***
   !&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
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : pi
   USE comvars2d, ONLY  : iscreenlog,tdeb,tfin,h0,dy,dx

   !! * Arguments

   !! * Local declarations
   LOGICAL                                     :: ex,l_posit
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=lchain)                       :: rec
   CHARACTER(LEN=19)                           :: dateread,tool_sectodat,comdate
   CHARACTER(LEN=lchain), ALLOCATABLE, DIMENSION(:)  :: name_patch
   INTEGER                                     :: eof,i,j,nn,n,m,kk
   REAL(KIND=rsh)                              :: tool_latlon2i,tool_latlon2j,rx,ry,dl
   INTEGER, ALLOCATABLE, DIMENSION(:)          :: imin_patch,imax_patch,   &
                                                  jmin_patch,jmax_patch,   &
                                                  nbpart_patch
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: ray_patch
   INTEGER, ALLOCATABLE, DIMENSION(:)          :: istep_patch,jstep_patch
   REAL(KIND=rlg)                              :: tool_datosec,tool_ind2lat,tool_ind2lon, &
                                                  gmin,gmax,phimin,phimax

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

   ! save into simu.log
   !-------------------
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '***************** TRAJ_INIT2D.F90 ****************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'fichier definissant les caracteristiques des trajectoires : ',trim(file_trajec)

   ! Make sure the type of trajectory is correct
   eof = 0
   INQUIRE(file=file_trajec,exist=ex)
   IF (ex) THEN
     OPEN(49,file = file_trajec,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,*,iostat=eof) itypetraj
     IF (itypetraj/=1 .AND. itypetraj/=2 .AND. itypetraj/=3 .AND. itypetraj/=4) THEN
       print*,'Type of trajectory is not defined correctly'
       print*,'Must be 1 (circle patch) or 2 (rectangle patch)'
       print*,'or 3 () or 4 (discharge every hours)'
       print*,'Simulation stopped.'
       stop
     END IF
     CLOSE(49)
   END IF


   ! Read trajec.dat file

   eof = 0
   nb_patch = 0
   INQUIRE(file=file_trajec,exist=ex)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Initialize circle patches
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   IF (itypetraj==1 .AND. ex) THEN
     OPEN(49,file = file_trajec,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,*,iostat=eof) itypetraj
     READ(49,*,iostat=eof)
     DO WHILE(eof==0)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,'(a)',iostat=eof)
       IF (eof==0) nb_patch = nb_patch+1
       READ(49,'(a)',iostat=eof)
     END DO
     CLOSE(49)
     WRITE(iscreenlog,*) 'Number of circular patches for estimation of trajectories = ',nb_patch

   ! allocate arrays for trajectories variables
     ALLOCATE( name_patch(nb_patch),filetraj2dout(nb_patch) )
     ALLOCATE( nb_part(nb_patch), dtsave_traj(nb_patch) )
     ALLOCATE( tbegtraj(nb_patch),tendtraj(nb_patch),tsavtraj(nb_patch) )
     ALLOCATE( imin_patch(nb_patch),jmin_patch(nb_patch) )
     ALLOCATE( nbpart_patch(nb_patch),ray_patch(nb_patch) )

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

     DO n = 1,nb_patch

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

     ! Read starting date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tbegtraj(n) = tool_datosec(dateread)
       IF(tbegtraj(n)<tdeb .OR. tbegtraj(n)>tfin) THEN
         comdate=tool_sectodat(tdeb)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Starting date ,',trim(dateread),', is uncorrect,'
         write(iscreenlog,*) 'WARNING : consequently trajectory starts with departure of 3D run on', comdate
         tbegtraj(n)=tdeb
       ENDIF
       tsavtraj(n)=tbegtraj(n)

     ! Read ending date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tendtraj(n) = tool_datosec(dateread)
       IF(tendtraj(n)<tdeb .OR. tendtraj(n)>tfin) THEN
         comdate=tool_sectodat(tfin)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Ending date ,',trim(dateread),', is uncorrect, '
         write(iscreenlog,*) 'WARNING : consequently trajectory ends with run on', comdate
         tendtraj(n)=tfin
       ENDIF

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

     ! Location 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(n),jmin_patch(n)
       ELSE
         READ(49,*,iostat=eof) gmin,phimin
         imin_patch(n)=nint(tool_latlon2i(gmin,phimin))
         jmin_patch(n)=nint(tool_latlon2j(gmin,phimin))
         IF ( imin_patch(n)<imin .OR. imin_patch(n)>imax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           imin_patch(n)=MAX(imin_patch(n),imin)
           imin_patch(n)=MIN(imin_patch(n),imax) 
           print*,'The longitude of patch number ',n,'is :',imin_patch(n)
         END IF
         IF ( jmin_patch(n)<jmin .OR. jmin_patch(n)>jmax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           jmin_patch(n)=MAX(jmin_patch(n),jmin)
           jmin_patch(n)=MIN(jmin_patch(n),jmax) 
           print*,'The latitude of patch number ',n,'is :',jmin_patch(n)
         END IF
       END IF

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

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

     ! Read output file
       READ(49,'(a)',iostat=eof)  rec
       kk=index(rec,',|')
       IF (kk /= 0 ) THEN
         filetraj2dout(n)=rec(1:kk-1)
       ELSE
         filetraj2dout(n)=rec
       END IF
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)

     ! Estimate number of particle inside the circle patch
     nbmax_part=0
     DO n = 1,nb_patch
       nn=int(sqrt(4/pi*nbpart_patch(n)))+1
       dl=2*ray_patch(n)/nn
       m=0                               ! initialize nb_part for each patch
       IF (h0(imin_patch(n),jmin_patch(n)) == -valmanq) THEN
         nb_part(n)=0
         print*,'Patch on land'
         print*,'Check patch number',n
         print*,'in file "',trim(file_trajec),'"'
         print*,'Simulation stopped.'
         stop
       ELSE
         DO j=1,nn
         DO i=1,nn
           rx=-ray_patch(n)+real(i-1,rsh)*dl
           ry=-ray_patch(n)+real(j-1,rsh)*dl
           IF (sqrt(rx*rx+ry*ry) <= ray_patch(n)) m=m+1
         END DO
         END DO
       ENDIF
       nb_part(n)=m
       nbmax_part=max(nbmax_part,m)    ! maximum number of particles
       WRITE(iscreenlog,*) 'Number of particles inside circle number ',n,' = ',nb_part(n)
       IF ( m==0 ) THEN
         print*,'Problem with patch number',n
         print*,'in file "',trim(file_trajec),'"'
         print*,'Simulation stopped.'
         stop
       END IF
     END DO
     WRITE(iscreenlog,*) 'The maximum number of particles is ',nbmax_part

     ! Location of particles of each patch
     ALLOCATE( xpos(nb_patch,nbmax_part),ypos(nb_patch,nbmax_part) )
     ALLOCATE( lon_pos(nb_patch,nbmax_part),lat_pos(nb_patch,nbmax_part) )

     DO n = 1,nb_patch
       nn=int(sqrt(4/pi*nbpart_patch(n)))+1
       dl=2*ray_patch(n)/nn
       m=0
       DO j=1,nn
       DO i=1,nn
         rx=-ray_patch(n)+real(i-1,rsh)*dl
         ry=-ray_patch(n)+real(j-1,rsh)*dl

         IF (sqrt(rx*rx+ry*ry) <= ray_patch(n)) THEN
           m=m+1
                     ! i index of initial location
           xpos(n,m)=imin_patch(n)+rx/dx(imin_patch(n),jmin_patch(n))
                     ! j index of initial location
           ypos(n,m)=jmin_patch(n)+ry/dy(imin_patch(n),jmin_patch(n))
                     ! longitude of initial location
           lon_pos(n,m)=tool_ind2lon(xpos(n,m),ypos(n,m))
                     ! latitude of initial location
           lat_pos(n,m)=tool_ind2lat(xpos(n,m),ypos(n,m))
         END IF
       END DO
       END DO
     END DO

     DEALLOCATE(name_patch,imin_patch,jmin_patch,nbpart_patch,ray_patch)


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Initialize rectangular patches
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ELSEIF (itypetraj==2 .AND. ex) THEN
     OPEN(49,file = file_trajec,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,*,iostat=eof) itypetraj
     READ(49,*,iostat=eof)
     DO WHILE(eof==0)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,'(a)',iostat=eof)
       IF (eof==0) nb_patch = nb_patch+1
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)
     WRITE(iscreenlog,*) 'Number of rectangular patches for estimation of trajectories = ',nb_patch

   ! allocate arrays for trajectories variables
     ALLOCATE( name_patch(nb_patch),filetraj2dout(nb_patch) )
     ALLOCATE( nb_part(nb_patch), dtsave_traj(nb_patch) )
     ALLOCATE( tbegtraj(nb_patch),tendtraj(nb_patch),tsavtraj(nb_patch) )
     ALLOCATE( imin_patch(nb_patch),imax_patch(nb_patch) )
     ALLOCATE( jmin_patch(nb_patch),jmax_patch(nb_patch) )
     ALLOCATE( istep_patch(nb_patch),jstep_patch(nb_patch) )

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

     DO n = 1,nb_patch

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

     ! Read starting date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tbegtraj(n) = tool_datosec(dateread)
       IF(tbegtraj(n)<tdeb .OR. tbegtraj(n)>tfin) THEN
         comdate=tool_sectodat(tdeb)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Starting date ,',trim(dateread),', is uncorrect (before 3D simulations starting) . '
         write(iscreenlog,*) 'WARNING : consequently trajectory starts with departure of 3D run on', comdate
         tbegtraj(n)=tdeb
       ENDIF
       tsavtraj(n)=tbegtraj(n)

     ! Read ending date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tendtraj(n) = tool_datosec(dateread)
       IF(tendtraj(n)<tdeb .OR. tendtraj(n)>tfin) THEN
         comdate=tool_sectodat(tfin)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Ending date ,',trim(dateread),', is uncorrect, '
         write(iscreenlog,*) 'WARNING : consequently trajectory ends with run on', comdate
         tendtraj(n)=tfin
       ENDIF

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

     ! 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(n),imax_patch(n),jmin_patch(n),jmax_patch(n)
       ELSE
         READ(49,*,iostat=eof) gmin,gmax,phimin,phimax
         imin_patch(n)=nint(tool_latlon2i(gmin,phimin))
         imax_patch(n)=nint(tool_latlon2i(gmax,phimax))
         jmin_patch(n)=nint(tool_latlon2j(gmin,phimin))
         jmax_patch(n)=nint(tool_latlon2j(gmax,phimax))
         IF ( imin_patch(n)<imin .OR. imin_patch(n)>imax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           imin_patch(n)=MAX(imin_patch(n),imin)
           imin_patch(n)=MIN(imin_patch(n),imax)
           print*,'The western longitude of patch number ',n,'is :',imin_patch(n)
         END IF
         IF ( jmin_patch(n)<jmin .OR. jmin_patch(n)>jmax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           jmin_patch(n)=MAX(jmin_patch(n),jmin)
           jmin_patch(n)=MIN(jmin_patch(n),jmax)
           print*,'The southern latitude of patch number ',n,'is :',jmin_patch(n)
         END IF
         IF ( imax_patch(n)<imin .OR. imax_patch(n)>imax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           imax_patch(n)=MAX(imax_patch(n),imin)
           imax_patch(n)=MIN(imax_patch(n),imax)
           print*,'The eastern longitude of patch number ',n,'is :',imax_patch(n)
         END IF
         IF ( jmax_patch(n)<jmin .OR. jmax_patch(n)>jmax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           jmax_patch(n)=MAX(jmax_patch(n),jmin)
           jmax_patch(n)=MIN(jmax_patch(n),jmax)
           print*,'The northern latitude of patch number ',n,'is :',jmax_patch(n)
         END IF
       END IF

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

     ! Read output file
       READ(49,'(a)',iostat=eof)  rec
       kk=index(rec,',|')
       IF (kk /= 0 ) THEN
         filetraj2dout(n)=rec(1:kk-1)
       ELSE
         filetraj2dout(n)=rec
       END IF
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)

     ! Estimate number of particle inside the rectangular patch
     nbmax_part=0
     DO n = 1,nb_patch
       m=0                               ! initialize nb_part for each patch
       DO j=jmin_patch(n),jmax_patch(n),jstep_patch(n)
       DO i=imin_patch(n),imax_patch(n),istep_patch(n)
         IF (h0(i,j).gt.-valmanq) m=m+1
       END DO
       END DO
       nb_part(n)=m
       nbmax_part=max(nbmax_part,m)    ! maximum number of particles
       WRITE(iscreenlog,*) 'Number of particles inside rectangle number ',n
       IF (m==0 ) THEN
         print*,'Patch on land'
         print*,'Check patch number',n
         print*,'in file "',trim(file_trajec),'"'
         print*,'Simulation stopped.'
         stop
       END IF
     END DO
     WRITE(iscreenlog,*) 'The maximum number of particles is ',nbmax_part

     ! Location of particles of each patch
     ALLOCATE( xpos(nb_patch,nbmax_part),ypos(nb_patch,nbmax_part) )
     ALLOCATE( lon_pos(nb_patch,nbmax_part),lat_pos(nb_patch,nbmax_part) )
     DO n = 1,nb_patch
       m=0
       DO j=jmin_patch(n),jmax_patch(n),jstep_patch(n)
       DO i=imin_patch(n),imax_patch(n),istep_patch(n)
         IF (h0(i,j).gt.-valmanq) THEN
           m=m+1
           xpos(n,m)=FLOAT(i)                           ! i index of initial location
           ypos(n,m)=FLOAT(j)                           ! j index
           lon_pos(n,m)=tool_ind2lon(xpos(n,m),ypos(n,m))    ! longitude of initial location
           lat_pos(n,m)=tool_ind2lat(xpos(n,m),ypos(n,m))    ! latitude of initial location
         END IF
       END DO
       END DO
     END DO

   DEALLOCATE(name_patch,imin_patch,imax_patch,jmin_patch,jmax_patch)
   DEALLOCATE(istep_patch,jstep_patch)



   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Initialize a unique rectangular patch for estimating of residual currents
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ELSEIF (itypetraj==3 .AND. ex) THEN

   ! allocate arrays for residual estimates
     ALLOCATE( tbegtraj(1) )
     ALLOCATE( xposres(nres+1,imin:imax,jmin:jmax) )
     ALLOCATE( yposres(nres+1,imin:imax,jmin:jmax) )
     ALLOCATE( xmoyres(nres+1,imin:imax,jmin:jmax) )
     ALLOCATE( ymoyres(nres+1,imin:imax,jmin:jmax) )
     ALLOCATE( tbegres(nres+1),tprectraj(nres+1) )
     ALLOCATE( iswitj(0:nres+1) )
     ALLOCATE( ures(imin:imax,jmin:jmax),vres(imin:imax,jmin:jmax) )
     ALLOCATE( vitres(imin:imax,jmin:jmax),nomb(imin:imax,jmin:jmax) )
     ALLOCATE( uit(nres+1,imin:imax,jmin:jmax),vit(nres+1,imin:imax,jmin:jmax) )
     iswitj(:)=0

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

   ! Read starting date of trajectory
     READ(49,'(a)',iostat=eof) dateread
     tbegtraj(1) = tool_datosec(dateread)
     IF(tbegtraj(1)<tdeb .OR. tbegtraj(1)>tfin) THEN
       comdate=tool_sectodat(tdeb)
       write(iscreenlog,*) ' '
       write(iscreenlog,*) 'WARNING : Starting date ,',trim(dateread),', is uncorrect, '
       write(iscreenlog,*) 'WARNING : consequently residual currents estimates starts with departure of 2D run on', comdate
       tbegtraj(1)=tdeb
     ENDIF
  !   tsavtraj(1)=tbegtraj(1)+m2per

   ! Initialize variables
     DO j=jmin+2,jmax-1
       DO i=imin+2,imax-1
         xposres(:,i,j)=REAL(i,rsh)
         yposres(:,i,j)=REAL(j,rsh)
         xmoyres(:,i,j)=0.0_rsh
         ymoyres(:,i,j)=0.0_rsh
       END DO
     END DO


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! local discharge at every hours
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ELSEIF (itypetraj==4 .AND. ex) THEN
     OPEN(49,file = file_trajec,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,*,iostat=eof) itypetraj
     READ(49,*,iostat=eof)
     DO WHILE(eof==0)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,'(a)',iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,*,iostat=eof)
       READ(49,'(a)',iostat=eof)
       IF (eof==0) nb_patch = nb_patch+1
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)
     WRITE(iscreenlog,*) 'Number of local discharges for estimation of trajectories = ',nb_patch

   ! allocate arrays for trajectories variables
     ALLOCATE( name_patch(nb_patch),filetraj2dout(nb_patch) )
     ALLOCATE( nb_part(nb_patch), dtsave_traj(nb_patch) )
     ALLOCATE( tbegtraj(nb_patch),tendtraj(nb_patch),tsavtraj(nb_patch) )
     ALLOCATE( imin_patch(nb_patch) )
     ALLOCATE( jmin_patch(nb_patch) )
     nbmax_part=13
     ALLOCATE( xpos(nb_patch,nbmax_part),ypos(nb_patch,nbmax_part) )
     ALLOCATE( lon_pos(nb_patch,nbmax_part),lat_pos(nb_patch,nbmax_part) )
     nb_part(:)=nbmax_part

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

     DO n = 1,nb_patch

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

     ! Read starting date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tbegtraj(n) = tool_datosec(dateread)
       IF(tbegtraj(n)<tdeb .OR. tbegtraj(n)>tfin) THEN
         comdate=tool_sectodat(tdeb)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Starting date ,',trim(dateread),', is uncorrect (before 3D simulations starting) . '
         write(iscreenlog,*) 'WARNING : consequently trajectory starts with departure of 3D run on', comdate
         tbegtraj(n)=tdeb
       ENDIF
       tsavtraj(n)=tbegtraj(n)

     ! Read ending date of trajectory
       READ(49,'(a)',iostat=eof) dateread
       tendtraj(n) = tool_datosec(dateread)
       IF(tendtraj(n)<tdeb .OR. tendtraj(n)>tfin) THEN
         comdate=tool_sectodat(tfin)
         write(iscreenlog,*) ' '
         write(iscreenlog,*) 'WARNING : PATCH NUMBER : ',n
         write(iscreenlog,*) 'WARNING : Ending date ,',trim(dateread),', is uncorrect, '
         write(iscreenlog,*) 'WARNING : consequently trajectory ends with run on', comdate
         tendtraj(n)=tfin
       ENDIF

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

     ! 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(n),jmin_patch(n)
       ELSE
         READ(49,*,iostat=eof) gmin,phimin
         imin_patch(n)=nint(tool_latlon2i(gmin,phimin))
         jmin_patch(n)=nint(tool_latlon2j(gmin,phimin))
         IF ( imin_patch(n)<imin .OR. imin_patch(n)>imax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           imin_patch(n)=MAX(imin_patch(n),imin)
           imin_patch(n)=MIN(imin_patch(n),imax)
           print*,'The longitude of patch number ',n,'is :',imin_patch(n)
         END IF
         IF ( jmin_patch(n)<jmin .OR. jmin_patch(n)>jmax ) THEN
           print*,'LOCATION OUT OF THE DOMAIN'
           print*,'Have a look at file ',trim(file_trajec),' trajectory number ',n
           jmin_patch(n)=MAX(jmin_patch(n),jmin)
           jmin_patch(n)=MIN(jmin_patch(n),jmax)
           print*,'The latitude of patch number ',n,'is :',jmin_patch(n)
         END IF
       END IF

     ! Read output file
       READ(49,'(a)',iostat=eof)  rec
       kk=index(rec,',|')
       IF (kk /= 0 ) THEN
         filetraj2dout(n)=rec(1:kk-1)
       ELSE
         filetraj2dout(n)=rec
       END IF
       READ(49,*,iostat=eof)
     END DO
     CLOSE(49)

     DO n = 1,nb_patch
       IF (h0(imin_patch(n),jmin_patch(n)) == -valmanq) THEN
         nb_part(n)=0
         print*,'Patch on land'
         print*,'Check patch number',n
         print*,'in file "',trim(file_trajec),'"'
         print*,'Simulation stopped.'
         stop
       ELSE
         xpos(n,:)=imin_patch(n)                      ! i index of initial location
         ypos(n,:)=jmin_patch(n)                      ! j index
         lon_pos(n,:)=tool_ind2lon(xpos(n,:),ypos(n,:))    ! longitude of initial location
         lat_pos(n,:)=tool_ind2lat(xpos(n,:),ypos(n,:))    ! latitude of initial location
       END IF
     END DO

     DEALLOCATE(name_patch,imin_patch,jmin_patch)

   END IF  ! end test on itypetraj

   ! save initialization
   IF (itypetraj /= 3) CALL traj_save2d 

   END SUBROUTINE traj_init2d

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

   SUBROUTINE traj_2d(u,v)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_2d  ***
   !&E
   !&E ** Purpose : calcul de l advection dispersion d un nuage de
   !&E              particules (2D).
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : step2dxy.F90, step2dyx.F90,
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (P. Lazure, F. Dumas)  Original code
   !&E       !  01-2007 (V. Garnier dissociation partie initialisation / evolution)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY  : t,dt,dy,dx,h0
   USE output, ONLY     : fileresu
   USE trajectools

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

   !! * Local declarations
   INTEGER             :: n,m,i,j,k,ii,jj
   REAL(KIND=rsh)      :: s0
   REAL(KIND=rlg)      :: dtm,tool_ind2lon,tool_ind2lat
   REAL(KIND=rsh), DIMENSION(kmax,imin:imax,jmin:jmax)  :: ui,vi
   REAL(KIND=riosh), ALLOCATABLE, DIMENSION(:,:)        :: wres


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

   ! Initialization
   s0=0.0_rsh


   ! Residual estimates (itypetraj == 3) or usual trajectories (other itypetraj) ?

   IF (itypetraj == 3) THEN  ! residual current estimates

     ! Starts loop on hours
     DO k = 1,nres

       ! Test on starting date for residual estimates
       IF (t >= tbegtraj(1)+(dble(k-1)*3600.0_rlg)) THEN

         ! Initializes switch and starting times
         IF (iswitj(k) == 0) THEN
           iswitj(k)=1
           tbegres(k)=t
           tprectraj(k)=t
         END IF

         ! Residual estimates
         IF (iswitj(k) == 1) THEN

           ! Initializes time step
           IF (t < tbegres(k)+m2per) THEN
            dtm=t-tprectraj(k)
           ELSE
            dtm=tbegres(k)+m2per-tprectraj(k)
            iswitj(k)=2
           END IF

           ui(1,:,:)=u(:,:)
           vi(1,:,:)=v(:,:)
           DO j = jmin+2,jmax-1
           DO i = imin+2,imax-1
           ! horizontal advection
             IF (xposres(k,i,j) /= 0.0_rsh .and. yposres(k,i,j) /= 0.0_rsh) &
               call avance(ui,vi,dtm,xposres(k,i,j),yposres(k,i,j),s0)
           ! average incrementation
             IF (xposres(k,i,j) /= 0.0_rsh .and. yposres(k,i,j) /= 0.0_rsh) THEN
               xmoyres(k,i,j)=xmoyres(k,i,j)+xposres(k,i,j)*dtm/m2per
               ymoyres(k,i,j)=ymoyres(k,i,j)+yposres(k,i,j)*dtm/m2per
             ELSE
               xmoyres(k,i,j)=0.0_rsh
               ymoyres(k,i,j)=0.0_rsh
             END IF
           END DO
           END DO
           tprectraj(k)=t
         END IF

       END IF   ! Ends test on starting date
     END DO   ! Ends loop on hours (nres)

     ! Save residual estimates
     IF (iswitj(nres) == 2) THEN

       nomb(:,:)= 0.0_rsh
       ures(:,:)= 0.0_rsh
       vres(:,:)= 0.0_rsh
       vitres(:,:)= 0.0_rsh

       DO j=jmin+2,jmax-1
       DO i=imin+2,imax-1
         DO k=1,nres
           ii=nint(xmoyres(k,i,j))
           jj=nint(ymoyres(k,i,j))
           IF (ii>=imin+1.and.ii<=imax.and.jj>=jmin+1.and.jj<=jmax) THEN
             ures(ii,jj)=ures(ii,jj)+(xposres(k,i,j)-real(i,rsh))* &
                 dx(nint(0.5_rsh*(xposres(k,i,j)+real(i,rsh))),nint(0.5_rsh*(yposres(k,i,j)+real(j,rsh))))/m2per
             vres(ii,jj)=vres(ii,jj)+(yposres(k,i,j)-real(j,rsh))* &
                  dy(nint(0.5_rsh*(xposres(k,i,j)+real(i,rsh))),nint(0.5_rsh*(yposres(k,i,j)+real(j,rsh))))/m2per
             nomb(ii,jj)= nomb(ii,jj)+1.0_rsh
           END IF
         END DO
       END DO
       END DO

       ! bathymetry mask
       DO j=jmin+1,jmax
       DO i=imin+1,imax
         IF (nomb(i,j) /= 0.0_rsh .and. h0(i,j) > 0.0_rsh) THEN
           ures(i,j)=ures(i,j)/nomb(i,j)
           vres(i,j)=vres(i,j)/nomb(i,j)
           vitres(i,j)=sqrt(vres(i,j)**2+ures(i,j)**2)
         ELSE
           ures(i,j)=0.0_rsh
           vres(i,j)=0.0_rsh
           vitres(i,j)=0.0_rsh
         END IF
       END DO
       END DO
!
       call ionc4_createvar(fileresu(1),'ures','m s-1','u residuel lagrangien',fill_value=0.0_riosh,dims='xy')
       call ionc4_createvar(fileresu(1),'vres','m s-1','v residuel lagrangien',fill_value=0.0_riosh,dims='xy')
       call ionc4_createvar(fileresu(1),'vitres','m s-1','norme residuel lagrangien',fill_value=0.0_riosh,dims='xy')
       ALLOCATE(wres(imin:imax,jmin:jmax))
       wres(:,:)=ures(:,:)
       call ionc4_write_xy(fileresu(1),'ures',wres,imin,imax,jmin,jmax)
       wres(:,:)=vres(:,:)
       call ionc4_write_xy(fileresu(1),'vres',wres,imin,imax,jmin,jmax)
       wres(:,:)=vitres(:,:)
       call ionc4_write_xy(fileresu(1),'vitres',wres,imin,imax,jmin,jmax)
       DEALLOCATE(wres)
       iswitj(nres)=4

     END IF   ! Ends test on iswitj=2

   ELSE   ! trajectory estimates

     !***************************
     ! Start loop on patches
     !***************************
     DO n = 1,nb_patch
       IF ((t >= tbegtraj(n)).and.(t <= tendtraj(n))) THEN
         DO m = 1,nb_part(n)

     !*************************************
     ! horizontal advection at every time step
     !*************************************


           IF (itypetraj <= 2) THEN  ! general case

             ! time step
             dtm=0.5_rlg*dt
             !xpos,ypos: variables in/out of routine avance
             ui(1,:,:)=u(:,:)
             vi(1,:,:)=v(:,:)
             CALL avance(ui,vi,dtm,xpos(n,m),ypos(n,m),s0)

             ! new indexes after particle s advection
             lon_pos(n,m)=tool_ind2lon(xpos(n,m),ypos(n,m))
             lat_pos(n,m)=tool_ind2lat(xpos(n,m),ypos(n,m))

           ELSE IF (itypetraj==4 .AND. t >= tbegtraj(n)+(m-1)*3600.0_rlg) THEN

             ! time step
             dtm=0.5_rlg*dt
             !xpos,ypos: variables in/out of routine avance
             ui(1,:,:)=u(:,:)
             vi(1,:,:)=v(:,:)
             CALL avance(ui,vi,dtm,xpos(n,m),ypos(n,m),s0)

             ! new indexes after particle s advection
             lon_pos(n,m)=tool_ind2lon(xpos(n,m),ypos(n,m))
             lat_pos(n,m)=tool_ind2lat(xpos(n,m),ypos(n,m))

           END IF  ! Ends test on trajectory s type

         END DO   ! ends loop on number of particles
       END IF   ! Ends test on trajectory s duration
     END DO   ! Ends loop on number of patches

   END IF   ! Ends swith on residual estimates or trajectory estimates

   !*********
   ! 0utputs
   !*********
   CALL traj_save2d

   END SUBROUTINE traj_2d


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

   SUBROUTINE traj_save2d

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE traj_save2d  ***
   !&E
   !&E ** Purpose : Save trajectories
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_init2d, step.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  01-2007 (V. Garnier)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY  : nbouc,date_ref,t

   !! * Arguments

   !! * Local declarations
   INTEGER                                     :: n,m
   REAL(KIND=rlg)                              :: tool_ind2lat,tool_ind2lon
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)   :: xpos_out,ypos_out,   &
                                                  lat_out,lon_out
   CHARACTER(LEN=19)                           :: dateread,tool_sectodat

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

   ALLOCATE( xpos_out(nbmax_part),ypos_out(nbmax_part) )
   ALLOCATE( lat_out(nbmax_part),lon_out(nbmax_part) )

   ! Create output file
   IF (nbouc == 0) THEN
     DO n=1,nb_patch
       CALL ionc4_createfile_traj(filetraj2dout(n),nb_part(n),0,0)
       CALL ionc4_torigine(filetraj2dout(n),date_ref)
       CALL ionc4_createvar_traj(filetraj2dout(n), 'lat','degrees_north','latitude')
       CALL ionc4_createvar_traj(filetraj2dout(n), 'lon','degrees_east', 'longitude')
     END DO
   END IF

   DO n=1,nb_patch
     IF ((t >= tsavtraj(n)).and.(t <= tendtraj(n))) THEN
       DO m = 1,nb_part(n)
         xpos_out(m)=xpos(n,m)
         ypos_out(m)=ypos(n,m)
         lat_out(m)=tool_ind2lat(xpos(n,m),ypos(n,m))
         lon_out(m)=tool_ind2lon(xpos(n,m),ypos(n,m))
       END DO

       CALL ionc4_write_time(filetraj2dout(n),0,t)
       CALL ionc4_write_trajt(filetraj2dout(n),'lat',lat_out,1,nb_part(n),0)
       CALL ionc4_write_trajt(filetraj2dout(n),'lon',lon_out,1,nb_part(n),0)
       tsavtraj(n)=tsavtraj(n)+dtsave_traj(n)*3600.0_rlg
     END IF
   END DO

   END SUBROUTINE traj_save2d
#endif

 END MODULE
