#include "FLT_OPTIONS.h"

CBOP 0
C !ROUTINE: FLT_TRAJ

C !INTERFACE:
      SUBROUTINE FLT_TRAJ (
     I                      myTime, myIter, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE FLT_TRAJ
C     | o This routine samples the model state at float position
C     |   every flt_int_traj time steps and writes output.
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
C     == global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "FLT_SIZE.h"
#include "FLT.h"
#include "FLT_BUFF.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_SIZE.h"
#include "W2_EXCH2_TOPOLOGY.h"
#endif

C     !INPUT PARAMETERS:
C     myTime :: current time in simulation
C     myIter :: current iteration number
C     myThid :: my Thread Id number
      _RL myTime
      INTEGER myIter, myThid

C     !FUNCTIONS:
      _RL FLT_MAP_K2R
      EXTERNAL FLT_MAP_K2R

C     !LOCAL VARIABLES:
      INTEGER bi, bj, nFlds
      INTEGER ip, kp, ii
      _RL ix, jy, i0x, j0y, xx, yy, zz
      _RL uu, vv, tt, ss, pp

      INTEGER imax
      PARAMETER (imax=13)
      _RL tmp(imax)
      _RL npart_read, npart_times
      _RS dummyRS(1)
      INTEGER fp, ioUnit, irecord
      CHARACTER*(MAX_LEN_FNAM) fn
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_EXCH2
      INTEGER nT
#endif
CEOP

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   set number of fields to write
      nFlds = 0
      IF ( flt_selectTrajOutp.GE.1 ) nFlds = nFlds + 8
      IF ( flt_selectTrajOutp.GE.2 ) nFlds = nFlds + 5

C--   check buffer size
      IF ( nFlds.GT.fltBufDim ) THEN
         _BEGIN_MASTER(myThid)
         WRITE(msgBuf,'(3(A,I4))') ' FLT_TRAJ: fltBufDim=', fltBufDim,
     &                             ' too small (<', nFlds, ' )'
         CALL PRINT_ERROR( msgBuf, myThid )
         WRITE(msgBuf,'(2A)')     ' FLT_TRAJ: => increase fltBufDim',
     &                            ' in "FLT_SIZE.h" & recompile'
         CALL PRINT_ERROR( msgBuf, myThid )
         _END_MASTER(myThid)
         CALL ALL_PROC_DIE( myThid )
         STOP 'ABNORMAL END: S/R FLT_TRAJ'
      ENDIF

      IF ( myIter.EQ.nIter0 .OR. flt_selectTrajOutp.LE.0 ) RETURN

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Calculate position + other fields at float position and fill up IO-buffer

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_EXCH2
         nT = W2_myTileList(bi,bj)
         i0x = DFLOAT( exch2_txGlobalo(nT) - 1 )
         j0y = DFLOAT( exch2_tyGlobalo(nT) - 1 )
#else
         i0x = DFLOAT( myXGlobalLo-1 + (bi-1)*sNx )
         j0y = DFLOAT( myYGlobalLo-1 + (bj-1)*sNy )
#endif
         DO ip=1,npart_tile(bi,bj)

            ix = ipart(ip,bi,bj)
            jy = jpart(ip,bi,bj)
            CALL FLT_MAP_IJLOCAL2XY( xx, yy,
     I                               ix, jy, bi,bj, myThid )
            zz = FLT_MAP_K2R( kpart(ip,bi,bj),bi,bj,myThid )
            kp = NINT(kpart(ip,bi,bj))
            tmp(1) = npart(ip,bi,bj)
            tmp(2) = myTime
            tmp(3) = xx
            tmp(4) = yy
            tmp(5) = zz
            tmp(6) = ix + i0x
            tmp(7) = jy + j0y
            tmp(8) = kpart(ip,bi,bj)

            IF ( ( flt_selectTrajOutp.GE.2 )   .AND.
     &           ( myTime.GE.tstart(ip,bi,bj)) .AND.
     &           ( tend(ip,bi,bj).EQ.-1. .OR. myTime.LE.tend(ip,bi,bj))
     &         ) THEN
              IF ( kp.LT.1 .OR. kp.GT.Nr ) THEN
                WRITE(msgBuf,'(2A,I8)') '** WARNING ** FLT_TRAJ: ',
     &            ' illegal value for kp=',kp
                CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                              SQUEEZE_RIGHT, myThid )
                WRITE(msgBuf,'(A,1P5E20.13)')
     &            ' FLT_TRAJ: ', (flt_io_buff(ii,ip,bi,bj),ii=1,5)
                CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                              SQUEEZE_RIGHT, myThid )
c               CALL PRINT_ERROR( msgBuf, myThid )
c               STOP 'ABNORMAL END: S/R FLT_TRAJ'
C-- jmc: not sure if this is right but added to avoid Pb in FLT_BILINEAR:
                kp = MIN( MAX(kp,1), Nr)
              ENDIF
              CALL FLT_BILINEAR  (ix,jy,uu,uVel,  kp,1,bi,bj,myThid)
              CALL FLT_BILINEAR  (ix,jy,vv,vVel,  kp,2,bi,bj,myThid)
              CALL FLT_BILINEAR2D(ix,jy,pp,etaN,     0,bi,bj,myThid)
              CALL FLT_BILINEAR  (ix,jy,tt,theta, kp,0,bi,bj,myThid)
              CALL FLT_BILINEAR  (ix,jy,ss,salt,  kp,0,bi,bj,myThid)
              tmp( 9) = pp
              tmp(10) = uu
              tmp(11) = vv
              tmp(12) = tt
              tmp(13) = ss
            ELSEIF ( flt_selectTrajOutp.GE.2 ) THEN
              tmp( 9) = flt_nan
              tmp(10) = flt_nan
              tmp(11) = flt_nan
              tmp(12) = flt_nan
              tmp(13) = flt_nan
            ENDIF

            DO ii=1,nFlds
              flt_io_buff(ii,ip,bi,bj) = tmp(ii)
            ENDDO

         ENDDO

       ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Write shared buffer to file

      _BARRIER
      _BEGIN_MASTER(myThid)

      fn = 'float_trajectories'
      fp = writeBinaryPrec

      DO bj=1,nSy
       DO bi=1,nSx

C (1) read actual number floats from file (if exists)
         ioUnit = -2
         CALL MDS_READVEC_LOC(  fn, fp, ioUnit, 'RL', nFlds,
     O                          tmp, dummyRS,
     I                          bi, bj, 1, myThid )
         IF ( ioUnit.GT. 0 ) THEN
            npart_read  = tmp(1)
            npart_times = tmp(5)
            ii = NINT(tmp(7))
C-       for backward compatibility with old trajectory files:
            IF ( ii.EQ.0 ) ii = 13
            IF ( ii.NE.nFlds ) THEN
              WRITE(msgBuf,'(A,I4,A)')
     &            'FLT_TRAJ: nFlds=', nFlds,' different from'
              CALL PRINT_ERROR( msgBuf, myThid )
              WRITE(msgBuf,'(3A,I4,A)')
     &            'previous file (',fn(1:18),') value =',ii
              CALL PRINT_ERROR( msgBuf, myThid )
              CALL ALL_PROC_DIE( 0 )
              STOP 'ABNORMAL END: S/R FLT_TRAJ'
            ENDIF
C-       close the read-unit (safer to use a different unit for writing)
            CLOSE( ioUnit )
         ELSE
            npart_read  = 0.
            npart_times = 0.
            tmp(2)      = myTime
         ENDIF

C (2) write new actual number floats and time axis into file
C-    the standard routine mds_writevec_loc can be used here

C     total number of records in this file
         tmp(1) = DBLE(npart_tile(bi,bj))+npart_read
C     first time of writing floats (do not change when written)
c        tmp(2) = tmp(2)
C     current time
         tmp(3) = myTime
C     timestep
         tmp(4) = flt_int_traj
C     total number of timesteps
         tmp(5) = npart_times + 1.
C     total number of floats
         tmp(6) = max_npart
C     total number of fields
         tmp(7) = nFlds
         DO ii=8,nFlds
           tmp(ii) = 0.
         ENDDO
         ioUnit = -1
         CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit, 'RL', nFlds,
     &                          tmp, dummyRS,
     &                          bi, bj, -1, myIter, myThid )

         DO ip=1,npart_tile(bi,bj)
C (3) write float positions into file
            irecord = npart_read+ip+1
            IF ( ip.NE.npart_tile(bi,bj) ) irecord = -irecord
            CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit, 'RL', nFlds,
     I                             flt_io_buff(1,ip,bi,bj), dummyRS,
     I                             bi, bj, irecord, myIter, myThid )
         ENDDO
         CLOSE( ioUnit )

       ENDDO
      ENDDO

      _END_MASTER(myThid)
      _BARRIER

      RETURN
      END
