C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_traj.F,v 1.14 2012/03/30 18:25:03 jmc Exp $
C $Name: $
#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
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