 MODULE trackpoints
   !!======================================================================
   !!                   ***  MODULE  trackpoints ***
   !! Outputs : save data every time step at point locations
   !!
   !!   History :
   !!        !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !!
   !!======================================================================
   !! * Modules used

#include "toolcpp.h"

   USE parameters, ONLY : riosh,rsh,rlg,imin,imax,jmin,jmax,kmax,limaxp1,ljmaxp1, &
                          valmanq,rg_valmanq_io,lchain,limin,limax,ljmin,ljmax,   &
                          liminm1,liminm2,limaxp2,ljminm1,ljminm2,ljmaxp2
   USE ionc4, ONLY : ionc4_createfile_sta,ionc4_torigine,ionc4_createvar_sta,      &
                     ionc4_write_sta,ionc4_write_time,ionc4_sync,ionc4_write_stat, &
                     ionc4_write_zstat,ionc4_write_z,ionc4_write,ionc4_vatt_chain
#ifdef key_sedim_MUSTANG
   USE parameters, ONLY :  ksdmin,ksdmax
#endif
                          
#ifdef key_MPI_2D
   USE toolmpi, ONLY: ADD_TAG,MASTER
#endif

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC track_points                ! routine called by main.F90 step.F90

   !! * Shared module variables
   LOGICAL, PUBLIC                    :: l_points = .false.
   CHARACTER(LEN=lchain), PUBLIC      :: file_point


   !! * Private variables
   INTEGER                            :: nbtrack,unit_track
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: itrack,jtrack,ktrack
   REAL(KIND=rlg),ALLOCATABLE,DIMENSION(:) :: xtrack,ytrack
   CHARACTER(LEN=lchain)              :: filetrack
   CHARACTER(LEN=3)                   :: vardims
   !REAL(KIND=rlg)                     :: tdebdiag,tfindiag,t0
   REAL(KIND=riosh),ALLOCATABLE,DIMENSION(:,:) :: wtrackz
   REAL(KIND=riosh),ALLOCATABLE,DIMENSION(:)   :: wtrack
#ifdef key_sedim_MUSTANG
   REAL(KIND=riosh),ALLOCATABLE,DIMENSION(:,:) :: wtracksedz
   CHARACTER(LEN=lchain)              :: filetrack_sed
#endif
   LOGICAL                            :: l_out_temp_track,l_out_sal_track, &
                                         l_out_ubt_track,l_out_u3d_track,  &
                                         l_out_ssh_track

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

  SUBROUTINE track_points(u,v,uz,vz,xe,sal,temp        &
#ifdef key_substance
                                               ,cv_wat &
#endif
                                               ,ect,bz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn_nxt  ***
   !&E
   !&E ** Purpose : save data at several locations
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : main, step
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2006-10  (V. Garnier) Original code
   !&E       !  2008-02  (M. Jouan) modifications of output formats
   !&E       !  2009-03  (V. Garnier) MPI2D
   !&E       !  2010-10  (B. Thouvenin) Corrections to wtrackz and evolution for the sedimentology
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules
   USE comvars2d,    ONLY : nbouc,iwarnlog,iscreenlog,t,cdate,suffix,l_modele2d,      &
                            h0,h0_g,root,date_ref, name_out_temp,name_out_sal,name_out_xe, &
                            name_out_u,name_out_v,name_out_uz,name_out_vz,name_out_h0,&
                            lon2d_g,lat2d_g
   USE comvars3d,    ONLY : sig
   USE comdiag,      ONLY : tdebdiag,tfindiag,t0,pasordiag
#if defined key_siggen || defined key_gencoord
   USE comsiggen,   ONLY : hc_sig,Csu_sig,b_sig,theta_sig
#endif
#ifdef key_substance
   USE comsubstance, ONLY : irk_fil,name_var,l_out_log,unit_var,nv_tot,nv_adv
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : ksmi,ksma,dzs,cv_sed,poro
#endif
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: xe
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: u,v
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: sal,temp
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: uz,vz
   REAL(KIND=rsh), DIMENSION(kmax,limin:limax,ljmin:ljmax),         INTENT(in) :: ect,bz
#ifdef key_substance
   REAL(KIND=rsh),DIMENSION(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: cv_wat
#endif

   !! * Local declarations
   LOGICAL                                        :: l_pos_ij,ex
   CHARACTER(LEN=5)                               :: comment
   CHARACTER(LEN=1000)                            :: chain
   CHARACTER(LEN=19)                              :: date_startdiag,date_enddiag
   CHARACTER(LEN=lchain),ALLOCATABLE,DIMENSION(:) :: nametrack, filetrac
   INTEGER                                        :: i,j,k,l,eof,ir
   REAL(KIND=rsh)                                 :: tool_latlon2j,tool_latlon2i,x,y
   REAL(KIND=rlg)                                 :: g,phi
   REAL(KIND=rlg)                                 :: tool_datosec
   REAL(KIND=riosh),DIMENSION(1:kmax)             :: sigecrit

#ifdef key_MPI_2D
   INTEGER                            :: kk
   INTEGER,ALLOCATABLE,DIMENSION(:)   :: itrack_mpi,jtrack_mpi,ktrack_mpi
   REAL(KIND=rlg),ALLOCATABLE,DIMENSION(:) :: xtrack_mpi,ytrack_mpi
   CHARACTER(LEN=lchain),ALLOCATABLE,DIMENSION(:) :: nametrack_mpi
#endif

#if defined key_substance & defined key_sedim_MUSTANG
   REAL(KIND=riosh),DIMENSION(ksdmin:ksdmax)  :: niv
#endif

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

   PRINT_DBG*, 'enter TRACKPOINTS routine'

   ! INITIALIZATION

   IF (nbouc == 0) THEN

   ! save into simu.log
   !-------------------
     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) '************* POINTSUIVI.F90 *********************'
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'TRACKING FILE :',TRIM(file_point)
     ENDIF_MPI

     t0=t           ! reference date (pas besoin en netcdf)
     nbtrack=0
     eof=0
!cval     nrec_track=0

   ! read points.dat file
     INQUIRE(file=file_point,exist=ex)
     IF (ex) THEN
       OPEN(55,file=file_point,form='formatted')
       comment='debut'
       DO WHILE (comment /= '*****')
         READ(55,'(a)',iostat=eof) comment
       END DO
       DO WHILE(eof==0)
         READ(55,'(a)',iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         READ(55,*,iostat=eof)
         IF(eof==0) nbtrack=nbtrack+1
         READ(55,*,iostat=eof)
       END DO
       CLOSE(55)

   ! save into simu.log
       IF_MPI (MASTER) THEN
       WRITE(iscreenlog,*) ' '
       WRITE(iscreenlog,*) 'NB OF TRACKING POINTS :',nbtrack
       ENDIF_MPI

   ! allocate arrays
       ALLOCATE(itrack(nbtrack),jtrack(nbtrack),ktrack(nbtrack))
       ALLOCATE(xtrack(nbtrack),ytrack(nbtrack),wtrack(nbtrack))
       ALLOCATE(wtrackz(kmax,nbtrack))
#ifdef key_sedim_MUSTANG
       ALLOCATE(wtracksedz(ksdmax,nbtrack))
#endif
       ALLOCATE(nametrack(nbtrack),filetrac(nbtrack))

   ! READ DATA

       OPEN(55,file=file_point,form='formatted')
       comment='debut'
       DO WHILE (comment /= '*****')
         READ(55,'(a)',iostat=eof) comment
       END DO
       DO k=1,nbtrack
  ! file name
         READ(55,'(a)',iostat=eof) nametrack(k)
  ! dates
         READ(55,'(a)',iostat=eof) date_startdiag
         READ(55,'(a)',iostat=eof) date_enddiag
         !tdebdiag = tool_datosec(date_startdiag)
         tdebdiag = MAX(tdebdiag,t)
         !tfindiag = tool_datosec(date_enddiag)
  ! time step
         READ(55,*) pasordiag
         pasordiag = pasordiag*3600.0_rlg
  ! location of tracking point
         READ(55,*,iostat=eof) l_pos_ij
         IF (.NOT.l_pos_ij) THEN
           READ(55,*) g,phi,ktrack(k)
           xtrack(k)=g
           ytrack(k)=phi
           x=tool_latlon2i(g,phi)
           y=tool_latlon2j(g,phi)
           itrack(k)=NINT(x)
           jtrack(k)=NINT(y)
         ELSE
           READ(55,*) itrack(k),jtrack(k),ktrack(k)
           xtrack(k)=lon2d_g(itrack(k),jtrack(k))
           ytrack(k)=lat2d_g(itrack(k),jtrack(k))
         END IF
         IF(l_modele2d) ktrack(k)=1
         IF (itrack(k)<imin .OR. itrack(k)>imax) THEN
           PRINT*,'tracking ',TRIM(nametrack(k))
           PRINT*,'index itrack ',itrack(k),' is out of range'
           PRINT*,'simulation stopped'
           STOP
         END IF
         IF (jtrack(k)<jmin .OR. jtrack(k)>jmax) THEN
           PRINT*,'tracking ',TRIM(nametrack(k))
           PRINT*,'index jtrack ',jtrack,' is out of range'
           PRINT*,'simulation stopped'
           STOP
         END IF
   ! look for the nearest point into the water
         ir=0
         IF (h0_g(itrack(k),jtrack(k))==-valmanq) THEN
           ir=ir+1
           IF(ir==1) PRINT*,'tracking point ',TRIM(nametrack(k)), ' is moved'
           DO j=jtrack(k)-ir,jtrack(k)+ir
             DO i=itrack(k)-ir,itrack(k)+ir
               IF(h0_g(i,j) /=-valmanq) THEN
                 itrack(k)=i
                 jtrack(k)=j
               END IF
             END DO
           END DO
         END IF
         IF(ir/=0) THEN
           xtrack(k)=lon2d_g(itrack(k),jtrack(k))
           ytrack(k)=lat2d_g(itrack(k),jtrack(k))
           PRINT*,'to the cell (i,j) : ',itrack(k),jtrack(k)
           PRINT*,'at the location :', xtrack(k),ytrack(k)
           IF (h0_g(itrack(k),jtrack(k))==-valmanq) PRINT*,'but point is still at land !'
         END IF
  ! data to be saved
         READ(55,*,iostat=eof) l_out_ssh_track
         READ(55,*,iostat=eof) l_out_ubt_track
         READ(55,*,iostat=eof) l_out_u3d_track
         READ(55,*,iostat=eof) l_out_sal_track
         READ(55,*,iostat=eof) l_out_temp_track
  ! end tracking point
         READ(55,*,iostat=eof) comment
       END DO   ! end loop on nbtrack
       CLOSE(55)

   ! save into simu.log
       DO k=1,nbtrack
         IF_MPI (MASTER) THEN
         WRITE(iscreenlog,*) ' '
         WRITE(iscreenlog,*) TRIM(nametrack(k))
         WRITE(iscreenlog,*) 'LOCATION ',itrack(k),jtrack(k),ktrack(k)
         WRITE(iscreenlog,*) 'TRACK sea surface height : ', l_out_ssh_track
         WRITE(iscreenlog,*) 'TRACK barotropic velocity : ', l_out_ubt_track
         WRITE(iscreenlog,*) 'TRACK baroclinic velocity : ', l_out_u3d_track
         WRITE(iscreenlog,*) 'TRACK salinity : ', l_out_sal_track
         WRITE(iscreenlog,*) 'TRACK temperature : ', l_out_temp_track
         ENDIF_MPI
       END DO

     END IF   ! fin du test sur l existence du fichier

#ifdef key_MPI_2D
     ALLOCATE(itrack_mpi(nbtrack),jtrack_mpi(nbtrack),ktrack_mpi(nbtrack))
     ALLOCATE(xtrack_mpi(nbtrack),ytrack_mpi(nbtrack))
     ALLOCATE(nametrack_mpi(nbtrack))
     kk=0
     DO k=1,nbtrack
       IF((jtrack(k).GE.ljmin) .AND. (jtrack(k).LE.ljmax) .AND. (itrack(k).GE.limin) .AND. (itrack(k).LE.limax)) THEN
             kk=kk+1
             jtrack_mpi(kk)=jtrack(k)
             itrack_mpi(kk)=itrack(k)
             ktrack_mpi(kk)=ktrack(k)
             xtrack_mpi(kk)=xtrack(k)
             ytrack_mpi(kk)=ytrack(k)
             nametrack_mpi(kk)=nametrack(k)
       ENDIF
     ENDDO
     nbtrack=kk
     ! nbtrack differs with each CPU depending on the location of the tracking-point
     ! consequently, all the variables must be reallocated for each CPU
     DEALLOCATE(itrack,jtrack,ktrack,xtrack,ytrack,nametrack,wtrackz,wtrack,filetrac)
     ALLOCATE(itrack(nbtrack),jtrack(nbtrack),ktrack(nbtrack))
     ALLOCATE(filetrac(nbtrack),xtrack(nbtrack),ytrack(nbtrack))
     ALLOCATE(nametrack(nbtrack),wtrackz(kmax,nbtrack),wtrack(nbtrack))
#ifdef key_sedim_MUSTANG
     DEALLOCATE(wtracksedz)
     ALLOCATE(wtracksedz(ksdmax,nbtrack))
#endif
     jtrack(1:nbtrack)=jtrack_mpi(1:nbtrack)
     itrack(1:nbtrack)=itrack_mpi(1:nbtrack)
     ktrack(1:nbtrack)=ktrack_mpi(1:nbtrack)
     xtrack(1:nbtrack)=xtrack_mpi(1:nbtrack)
     ytrack(1:nbtrack)=ytrack_mpi(1:nbtrack)
     nametrack(1:nbtrack)=nametrack_mpi(1:nbtrack)
     DEALLOCATE(itrack_mpi,jtrack_mpi,ktrack_mpi,xtrack_mpi,ytrack_mpi,nametrack_mpi)
#endif

    IF (nbtrack/=0) THEN 

   ! CREATE OUTPUT FILE
   !-------------------

   ! ascii file
#if ! defined  key_MPI_2D
     DO k=1,nbtrack
       unit_track=255+k
       filetrac(k)=TRIM(root)//TRIM(nametrack(k))//'_tracking_'//TRIM(suffix)//'.txt'
       OPEN(unit_track,file=TRIM(filetrac(k)),form='formatted')
       WRITE(unit_track,156)'Nb_jours','Date_complete','Hauteur_eau','Elevation','Vzonale', &
#ifdef key_substance
                          'Vmeridienne','Salinite','Temperature', &
                         (TRIM(name_var(i)),i=1,nv_tot)
#else
                         'Vmeridienne','Salinite','Temperature'
#endif
     END DO
#endif

   ! netcdf file

     filetrack=TRIM(root)//'tracking_'//TRIM(suffix)//'.nc'
#ifdef key_sedim_MUSTANG
     filetrack_sed=TRIM(root)//'tracking_sed_'//TRIM(suffix)//'.nc'
#endif
#ifdef key_MPI_2D
     CALL ADD_TAG(filetrack,LEN(filetrack))
#ifdef key_sedim_MUSTANG
     CALL ADD_TAG(filetrack_sed,LEN(filetrack_sed))
#endif
#endif
     CALL ionc4_createfile_sta(TRIM(filetrack),1,kmax,1,kmax,nbtrack,0,xtrack,ytrack)
     CALL ionc4_torigine(TRIM(filetrack),date_ref)

     CALL ionc4_createvar_sta(TRIM(filetrack),'staname','station name',dims='s')
     CALL ionc4_write_sta(TRIM(filetrack),'staname',nametrack,1,nbtrack)

     sigecrit(1:kmax)=sig(1:kmax)
     CALL ionc4_createvar_sta(TRIM(filetrack),'SIG','',    &
                             'sigma variable',fill_value=rg_valmanq_io,dims='z')
     CALL ionc4_write_z(TRIM(filetrack),'SIG',sigecrit,kmax,fill_value=rg_valmanq_io)

#if defined key_siggen || defined key_gencoord
       sigecrit(1:kmax)=Csu_sig(1:kmax)
       CALL ionc4_createvar_sta(TRIM(filetrack),'Csu_sig',' ','C(s) law at mid layer',&
                                     standard_name='ocean_s_coordinate_function_at_midlayer',dims='z',&
                                     valid_min=-2.0_riosh,valid_max=2.0_riosh,fill_value=rg_valmanq_io)
       CALL ionc4_write_z(TRIM(filetrack),'Csu_sig',sigecrit,kmax,fill_value=rg_valmanq_io)
       CALL ionc4_vatt_chain(TRIM(filetrack),'Csu_sig','formula_definition', &
          'C(k) = (1-b)*sinh(a*s(k))/sinh(a) + b*[tanh(a*(s(k)+0.5))-tanh(0.5*a)]/(2*tanh(0.5*a))')
       
       CALL ionc4_createvar_sta(TRIM(filetrack),'hc',units='m',dims='st',                              &
                                     long_name='depth over which the resolution is increased',      &
                                     standard_name='ocean_s_coordinate_thick_of_surface_resolution',&
                                     valid_min=-20.0_riosh,valid_max=500.0_riosh,fill_value=-rg_valmanq_io)
       
       CALL ionc4_createvar_sta(TRIM(filetrack),'b',' ','bottom control parameter',&
                                     standard_name='ocean_s_coordinate_bottom_control',dims='s',&
                                     valid_min=0.0_riosh,valid_max=1.0_riosh,fill_value=rg_valmanq_io)
       CALL ionc4_write(TRIM(filetrack),'b',b_sig)
      
       CALL ionc4_createvar_sta(TRIM(filetrack),'theta',' ','surface control parameter',&
                                     standard_name='ocean_s_coordinate_surface_control',dims='s',&
                                     valid_min=0.0_riosh,valid_max=20.0_riosh,fill_value=rg_valmanq_io)
       CALL ionc4_write(TRIM(filetrack),'theta',theta_sig)
#endif

     CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_h0),'m',  &
                           'bathymetrie par rapport au niv moyen',fill_value=rg_valmanq_io,dims='s')
     DO k=1,nbtrack
       wtrack(k)=h0_g(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_sta(TRIM(filetrack),TRIM(name_out_h0),wtrack,1,nbtrack)
     IF (l_out_ssh_track) THEN
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_xe),'m',  &
                           'elevation de la surface libre',fill_value=rg_valmanq_io,dims='st')
     ENDIF

     IF (l_out_ubt_track) THEN
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_u),'m s-1',    &
                             'composante zonale du courant',fill_value=rg_valmanq_io,dims='st')
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_v),'m s-1',    &
                             'composante meridienne du courant',fill_value=rg_valmanq_io,dims='st')
     ENDIF

     IF (l_out_u3d_track.AND..NOT.l_modele2d) THEN
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_uz),'m s-1',    &
                             'composante zonale du courant',fill_value=rg_valmanq_io,dims='szt')
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_vz),'m s-1',    &
                             'composante meridienne du courant',fill_value=rg_valmanq_io,dims='szt')
     ENDIF

     IF (l_modele2d) THEN
       vardims = 'st'
     ELSE
       vardims = 'szt'
     ENDIF

     IF (l_out_temp_track) THEN
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_temp),           &
                            'Deg Celsius','temperature',fill_value=rg_valmanq_io,dims=TRIM(vardims))
     ENDIF

     IF (l_out_sal_track) THEN
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_out_sal),           &
                            'psu','salinite',fill_value=rg_valmanq_io,dims=TRIM(vardims))
     ENDIF

#ifdef key_substance
     DO l = 1,nv_tot
       IF (l_out_log) THEN
         unit_var(irk_fil(l))='log('//TRIM(unit_var(irk_fil(l)))//')'
       ELSE
         unit_var(irk_fil(l))='('//TRIM(unit_var(irk_fil(l)))//')'
       END IF
       CALL ionc4_createvar_sta(TRIM(filetrack),TRIM(name_var(irk_fil(l))),  &
                                TRIM(unit_var(irk_fil(l))),                  &
                                TRIM(name_var(irk_fil(l))),fill_value=rg_valmanq_io,dims=TRIM(vardims))
     END DO

#ifdef key_sedim_MUSTANG
     CALL ionc4_createfile_sta(TRIM(filetrack_sed),1,ksdmax,1,ksdmax,nbtrack,0,xtrack,ytrack)
     CALL ionc4_torigine(TRIM(filetrack_sed),date_ref)

     niv = (/ (i, i=ksdmin,ksdmax) /)
     CALL ionc4_createvar_sta(TRIM(filetrack_sed),'layer','layer',  &
                                  'Number of layer in the sediment',fill_value=rg_valmanq_io,dims='z')
     CALL ionc4_write_z(TRIM(filetrack_sed),'layer',niv,ksdmax,fill_value=rg_valmanq_io)

     CALL ionc4_createvar_sta(TRIM(filetrack_sed),'NBNIV','quantite',  &
                                  'Nbre de couche de sediment',fill_value=rg_valmanq_io,dims='st')

     CALL ionc4_createvar_sta(TRIM(filetrack_sed),'EPTOT','m',  &
                                  'Epaiss total du sediment',fill_value=rg_valmanq_io,dims='st')

     CALL ionc4_createvar_sta(TRIM(filetrack_sed),'POROSITE','fraction',  &
                                  'porosite',fill_value=rg_valmanq_io,dims='szt')
     CALL ionc4_createvar_sta(TRIM(filetrack_sed),'DZS','m',  &
                                  'Epaiss de la couche',fill_value=rg_valmanq_io,dims='szt')
     DO l = 1,nv_adv
         IF (l_out_log) THEN
             unit_var(irk_fil(l))='log('//TRIM(unit_var(irk_fil(l)))//')'
         ELSE
             unit_var(irk_fil(l))='('//TRIM(unit_var(irk_fil(l)))//')'
         END IF
         CALL ionc4_createvar_sta(TRIM(filetrack_sed),TRIM(name_var(irk_fil(l))),  &
                                  unit_var(irk_fil(l)), &
                                  name_var(irk_fil(l)),fill_value=rg_valmanq_io,dims='szt')
     END DO
#endif
#endif

     DEALLOCATE(nametrack,filetrac)

    END IF  ! end test on nbtrack

    tdebdiag = tdebdiag - pasordiag    ! to save the initial condition

   END IF   ! end test on nbouc

   ! DEROULEMENT NORMAL
   !-------------------

   IF ( nbtrack/=0 .AND. t >= tdebdiag .AND. t <= tfindiag ) THEN

     tdebdiag = tdebdiag + pasordiag

!cval     nrec_track=nrec_track+1

   ! netcdf file

     CALL ionc4_write_time(TRIM(filetrack),0,t)
     !tina i need to re-write EVERYTHIHNG here by do loop of k.....

#if defined key_siggen || defined key_gencoord
     DO k=1,nbtrack
         wtrack(k)=hc_sig(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_stat(TRIM(filetrack),'hc',wtrack,1,nbtrack,0)
#endif

     IF (l_out_ssh_track) THEN
     DO k=1,nbtrack
       wtrack(k)=xe(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_out_xe),wtrack,1,nbtrack,0)
     END IF

     IF (l_out_ubt_track) THEN
     DO k=1,nbtrack
       wtrack(k)=u(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_out_u),wtrack,1,nbtrack,0)
     DO k=1,nbtrack
       wtrack(k)=v(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_out_v),wtrack,1,nbtrack,0)
     END IF

     IF (l_out_u3d_track.AND..NOT.l_modele2d) THEN
     DO k=1,nbtrack
       wtrackz(:,k)=uz(:,itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_zstat(TRIM(filetrack),TRIM(name_out_uz),wtrackz,1,nbtrack,kmax,0)
     DO k=1,nbtrack
       wtrackz(:,k)=vz(:,itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_zstat(TRIM(filetrack),TRIM(name_out_vz),wtrackz,1,nbtrack,kmax,0)
     END IF

     IF (l_out_sal_track) THEN
       IF (l_modele2d) THEN
         DO k=1,nbtrack
           wtrack(k)=sal(1,itrack(k),jtrack(k))
         END DO
         CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_out_sal),wtrack,1,nbtrack,0)
       ELSE
         DO k=1,nbtrack
           wtrackz(:,k)=sal(:,itrack(k),jtrack(k))
         END DO
         CALL ionc4_write_zstat(TRIM(filetrack),TRIM(name_out_sal),wtrackz,1,nbtrack,kmax,0)
       END IF
     END IF

     IF (l_out_temp_track) THEN
       IF (l_modele2d) THEN
         DO k=1,nbtrack
           wtrack(k)=temp(1,itrack(k),jtrack(k))
         END DO
         CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_out_temp),wtrack,1,nbtrack,0)
       ELSE
         DO k=1,nbtrack
           wtrackz(:,k)=temp(:,itrack(k),jtrack(k))
         END DO
         CALL ionc4_write_zstat(TRIM(filetrack),TRIM(name_out_temp),wtrackz,1,nbtrack,kmax,0)
       END IF
     END IF

#ifdef key_substance
       DO l = 1,nv_adv
         IF (l_modele2d) THEN
           IF (l_out_log) THEN
             DO k=1,nbtrack
               wtrack(k)=LOG(MAX(cv_wat(l,1,itrack(k),jtrack(k)),1.0e-30))/LOG(10.0_rsh)
             END DO
           ELSE
             DO k=1,nbtrack
               !wtrack(k)=MAX(cv_wat(l,1,itrack(k),jtrack(k)),1.0e-30)
               wtrack(k)=cv_wat(l,1,itrack(k),jtrack(k))
             END DO
           END IF
           CALL ionc4_write_stat(TRIM(filetrack),TRIM(name_var(irk_fil(l))),wtrack,1,nbtrack,0)
         ELSE
           IF (l_out_log) THEN
             DO k=1,nbtrack
               wtrackz(:,k)=LOG(MAX(cv_wat(l,:,itrack(k),jtrack(k)),1.0e-30))/LOG(10.0_rsh)
             END DO
           ELSE
             DO k=1,nbtrack
               !wtrackz(:,k)=MAX(cv_wat(l,:,itrack(k),jtrack(k)),1.0e-30)
               wtrackz(:,k)=cv_wat(l,:,itrack(k),jtrack(k))
             END DO
           END IF
           CALL ionc4_write_zstat(TRIM(filetrack),TRIM(name_var(irk_fil(l))),wtrackz,1,nbtrack,kmax,0)
         END IF
       END DO
#endif

     ! To write the data on the disk and not loose data in case of run crash
     CALL ionc4_sync(TRIM(filetrack))


#if defined key_substance & defined key_sedim_MUSTANG
     CALL ionc4_write_time(TRIM(filetrack_sed),0,t)
     DO k=1,nbtrack
       wtrack(k)=ksma(itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_stat(TRIM(filetrack_sed),'NBNIV',wtrack,1,nbtrack,0)
     wtracksedz(:,:)=0.0_riosh
     DO k=1,nbtrack
       wtrack(k)=0.0_rsh
       DO ir = ksmi(itrack(k),jtrack(k)),ksma(itrack(k),jtrack(k))
             wtrack(k)=wtrack(k)+dzs(ir,itrack(k),jtrack(k))
             wtracksedz(ir,k)=poro(ir,itrack(k),jtrack(k))
       END DO
     END DO
     CALL ionc4_write_stat(TRIM(filetrack_sed),'EPTOT',wtrack,1,nbtrack,0)
     CALL ionc4_write_zstat(TRIM(filetrack_sed),'POROSITE',wtracksedz,1,nbtrack,ksdmax,0)
     wtracksedz(:,:)=0.0_riosh
     DO k=1,nbtrack
       wtracksedz(:,k)=dzs(:,itrack(k),jtrack(k))
     END DO
     CALL ionc4_write_zstat(TRIM(filetrack_sed),'DZS',wtracksedz,1,nbtrack,ksdmax,0)
     wtracksedz(:,:)=0.0_riosh
       DO l = 1,nv_adv
           IF (l_out_log) THEN
             DO k=1,nbtrack
              DO ir = ksmi(itrack(k),jtrack(k)),ksma(itrack(k),jtrack(k))
                  wtracksedz(ir,k)=LOG(MAX(cv_sed(l,ir,itrack(k),jtrack(k)),1.0e-30))/LOG(10.0_rsh)
               END DO
             END DO
           ELSE
             DO k=1,nbtrack
               DO ir = ksmi(itrack(k),jtrack(k)),ksma(itrack(k),jtrack(k))
                  !wtracksedz(ir,k)=MAX(cv_sed(l,ir,itrack(k),jtrack(k)),1.0e-30)
                  wtracksedz(ir,k)=cv_sed(l,ir,itrack(k),jtrack(k))
               END DO
             END DO
           END IF
           CALL ionc4_write_zstat(TRIM(filetrack_sed),TRIM(name_var(irk_fil(l))),wtracksedz,1,nbtrack,ksdmax,0)
       END DO

     ! To write the data on the disk and not loose data in case of run crash
     CALL ionc4_sync(TRIM(filetrack_sed))

#endif

  ! ascii file

#ifdef key_MPI_2D
     WRITE(iwarnlog,*) ''
     WRITE(iwarnlog,*) 'WARNING : TRACKPOINTS routine'
     WRITE(iwarnlog,*) 'With MPI, it is not possible to save an ascii file yet'
#else

     DO k=1,nbtrack                        ! suivi des points
       unit_track=k+255
       WRITE(chain,155) (t-t0)/86400.0_rlg,cdate,h0(itrack(k),jtrack(k))+ &
                                                 xe(itrack(k),jtrack(k)), &
                                                 xe(itrack(k),jtrack(k)), &
                                                 u(itrack(k),jtrack(k)),  &
                                                 v(itrack(k),jtrack(k)),  &
                                                 sal(ktrack(k),itrack(k),jtrack(k)), &
#ifdef key_substance
                                                 temp(ktrack(k),itrack(k),jtrack(k)), &
                                       (cv_wat(i,ktrack(k),itrack(k),jtrack(k)),i=1,nv_tot)
#else
                                                 temp(ktrack(k),itrack(k),jtrack(k))
#endif

       DO l=1,LEN_TRIM(chain)
         IF (chain(l:l)=='.') chain(l:l)=','
       END DO

       WRITE(unit_track,'(a)') TRIM(chain)
     END DO
#endif

   END IF

  155 FORMAT(e15.7,';',a19,100(';',e13.5))
  156 FORMAT(a15,';',a19,6(';',a13),100(';',a13))

   PRINT_DBG*, 'exit TRACKPOINTS routine'

  END SUBROUTINE track_points

 END MODULE
