pkg/flt/flt_traj.F

C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_traj.F,v 1.4 2009/02/01 21:10:51 jmc Exp $
C $Name:  $

#include "FLT_OPTIONS.h"


      SUBROUTINE FLT_TRAJ (
     I                      myTime, myIter, myThid )

C     ==================================================================
C     SUBROUTINE FLT_TRAJ
C     ==================================================================
C     o This routine samples the model state at float position every
C       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 "GRID.h"
#include "DYNVARS.h"
#include "FLT.h"

C     == routine arguments ==
      _RL myTime
      INTEGER myIter, myThid

C     === Functions ==
      INTEGER ILNBLNK
      _RL global2local_i
      _RL global2local_j

C     == local variables ==
      INTEGER bi, bj, imax
      PARAMETER (imax=10)
      INTEGER ip, kp, ii
      _RL xx, yy
      _RL uu,vv,tt,ss, pp

      INTEGER ioUnit, irecord
      _RL tmp(imax)
      _RL npart_read,npart_times
      CHARACTER*(MAX_LEN_FNAM) fn
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      CHARACTER*(80) dataFName
      INTEGER iG,jG,IL
      LOGICAL exst
      LOGICAL globalFile

C     == end of interface ==

      fn = 'float_trajectories'

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

C (1) read actual number floats from file (if exists)
         IL=ILNBLNK( fn )
         iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
         jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
         WRITE(dataFname(1:80),'(2a,i3.3,a,i3.3,a)')
     &              fn(1:IL),'.',iG,'.',jG,'.data'
         INQUIRE( file=dataFname, exist=exst )
         IF (exst) THEN
            CALL FLT_MDSREADVECTOR(fn,globalFile,precFloat64,'RL',
     &                             imax,tmp,bi,bj,1,myThid)
            npart_read  = tmp(1)
            npart_times = tmp(5)
         ELSE
            npart_read  = 0.
            npart_times = 0.
            tmp(2)      = myTime
         ENDIF

C the standard routine mds_writevec_loc can be used here
C (2) WRITE new actual number floats and time axis into file
C
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
         DO ip=7,imax
            tmp(ip) = 0.
         ENDDO

         ioUnit = -1
         CALL MDS_WRITEVEC_LOC( fn, precFloat64, ioUnit,
     &                          'RL', imax, tmp,
     &                           bi,bj,-1, myIter, myThid )

         DO ip=1,npart_tile(bi,bj)

            xx=global2local_i(xpart(ip,bi,bj),bi,bj,myThid)
            yy=global2local_j(ypart(ip,bi,bj),bi,bj,myThid)
            kp = NINT(kpart(ip,bi,bj))
            tmp(1)  = npart(ip,bi,bj)
            tmp(2)  = myTime
            tmp(3)  = xpart(ip,bi,bj)
            tmp(4)  = ypart(ip,bi,bj)
            tmp(5)  = kpart(ip,bi,bj)

            IF ( ( 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: ', (tmp(ii),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  (xx,yy,uu,uVel,  kp,1,bi,bj,myThid)
              CALL FLT_BILINEAR  (xx,yy,vv,vVel,  kp,2,bi,bj,myThid)
              CALL FLT_BILINEAR2D(xx,yy,pp,etaN,     0,bi,bj,myThid)
              CALL FLT_BILINEAR  (xx,yy,tt,theta, kp,0,bi,bj,myThid)
              CALL FLT_BILINEAR  (xx,yy,ss,salt,  kp,0,bi,bj,myThid)

              tmp(6)   = uu
              tmp(7)   = vv
              tmp(8)   = tt
              tmp(9)   = ss
              tmp(10)  = pp
            ELSE
              tmp(6)   = flt_nan
              tmp(7)   = flt_nan
              tmp(8)   = flt_nan
              tmp(9)   = flt_nan
              tmp(10)  = flt_nan
            ENDIF

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, precFloat64, ioUnit,
     &                            'RL', imax, tmp,
     &                             bi,bj,irecord, myIter, myThid )

         ENDDO
         CLOSE( ioUnit )

       ENDDO
      ENDDO

      RETURN
      END
1	jmc	1.5	C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_traj.F,v 1.4 2009/02/01 21:10:51 jmc Exp $
2	jmc	1.2	C $Name: $
3	adcroft	1.1
4	jmc	1.2	#include "FLT_OPTIONS.h"
5	adcroft	1.1
6
7	jmc	1.2	SUBROUTINE FLT_TRAJ (
8	jmc	1.3	I myTime, myIter, myThid )
9	adcroft	1.1
10	jmc	1.2	C ==================================================================
11			C SUBROUTINE FLT_TRAJ
12			C ==================================================================
13			C o This routine samples the model state at float position every
14			C flt_int_traj time steps and writes output.
15			C ==================================================================
16	jmc	1.4
17			C !USES:
18			IMPLICIT NONE
19	adcroft	1.1
20	jmc	1.2	C == global variables ==
21			#include "SIZE.h"
22	adcroft	1.1	#include "EEPARAMS.h"
23	jmc	1.2	#include "PARAMS.h"
24			#include "GRID.h"
25	adcroft	1.1	#include "DYNVARS.h"
26			#include "FLT.h"
27
28	jmc	1.2	C == routine arguments ==
29	jmc	1.3	_RL myTime
30	jmc	1.2	INTEGER myIter, myThid
31	jmc	1.3
32	jmc	1.4	C === Functions ==
33			INTEGER ILNBLNK
34			_RL global2local_i
35			_RL global2local_j
36
37	jmc	1.3	C == local variables ==
38	adcroft	1.1	INTEGER bi, bj, imax
39	jmc	1.2	PARAMETER (imax=10)
40			INTEGER ip, kp, ii
41	jmc	1.3	_RL xx, yy
42			_RL uu,vv,tt,ss, pp
43	jmc	1.2
44	jmc	1.5	INTEGER ioUnit, irecord
45	adcroft	1.1	_RL tmp(imax)
46			_RL npart_read,npart_times
47			CHARACTER*(MAX_LEN_FNAM) fn
48	jmc	1.2	CHARACTER*(MAX_LEN_MBUF) msgBuf
49			CHARACTER*(80) dataFName
50			INTEGER iG,jG,IL
51			LOGICAL exst
52			LOGICAL globalFile
53	adcroft	1.1
54	jmc	1.2	C == end of interface ==
55	adcroft	1.1
56			fn = 'float_trajectories'
57
58			DO bj=myByLo(myThid),myByHi(myThid)
59	jmc	1.4	DO bi=myBxLo(myThid),myBxHi(myThid)
60	adcroft	1.1
61	jmc	1.2	C (1) read actual number floats from file (if exists)
62	adcroft	1.1	IL=ILNBLNK( fn )
63			iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
64			jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
65	jmc	1.2	WRITE(dataFname(1:80),'(2a,i3.3,a,i3.3,a)')
66	adcroft	1.1	& fn(1:IL),'.',iG,'.',jG,'.data'
67	jmc	1.2	INQUIRE( file=dataFname, exist=exst )
68			IF (exst) THEN
69	jmc	1.5	CALL FLT_MDSREADVECTOR(fn,globalFile,precFloat64,'RL',
70	jmc	1.2	& imax,tmp,bi,bj,1,myThid)
71	adcroft	1.1	npart_read = tmp(1)
72			npart_times = tmp(5)
73	jmc	1.2	ELSE
74	adcroft	1.1	npart_read = 0.
75			npart_times = 0.
76	jmc	1.2	tmp(2) = myTime
77			ENDIF
78	adcroft	1.1
79	jmc	1.5	C the standard routine mds_writevec_loc can be used here
80	jmc	1.2	C (2) WRITE new actual number floats and time axis into file
81			C
82			C total number of records in this file
83	adcroft	1.1	tmp(1) = DBLE(npart_tile(bi,bj))+npart_read
84	jmc	1.2	C first time of writing floats (do not change when written)
85	adcroft	1.1	c tmp(2) = tmp(2)
86	jmc	1.2	C current time
87			tmp(3) = myTime
88			C timestep
89	adcroft	1.1	tmp(4) = flt_int_traj
90	jmc	1.2	C total number of timesteps
91	adcroft	1.1	tmp(5) = npart_times + 1.
92	jmc	1.2	C total number of floats
93	adcroft	1.1	tmp(6) = max_npart
94	jmc	1.2	DO ip=7,imax
95	adcroft	1.1	tmp(ip) = 0.
96	jmc	1.2	ENDDO
97	adcroft	1.1
98	jmc	1.5	ioUnit = -1
99			CALL MDS_WRITEVEC_LOC( fn, precFloat64, ioUnit,
100			& 'RL', imax, tmp,
101			& bi,bj,-1, myIter, myThid )
102	adcroft	1.1
103	jmc	1.2	DO ip=1,npart_tile(bi,bj)
104	adcroft	1.1
105	jmc	1.2	xx=global2local_i(xpart(ip,bi,bj),bi,bj,myThid)
106			yy=global2local_j(ypart(ip,bi,bj),bi,bj,myThid)
107	jmc	1.4	kp = NINT(kpart(ip,bi,bj))
108	adcroft	1.1	tmp(1) = npart(ip,bi,bj)
109	jmc	1.2	tmp(2) = myTime
110			tmp(3) = xpart(ip,bi,bj)
111			tmp(4) = ypart(ip,bi,bj)
112	adcroft	1.1	tmp(5) = kpart(ip,bi,bj)
113
114	jmc	1.4	IF ( ( myTime.GE.tstart(ip,bi,bj)) .AND.
115	jmc	1.2	& ( tend(ip,bi,bj).EQ.-1. .OR. myTime.LE.tend(ip,bi,bj))
116			& ) THEN
117
118			IF ( kp.LT.1 .OR. kp.GT.Nr ) THEN
119			WRITE(msgBuf,'(2A,I8)') ' WARNING FLT_TRAJ: ',
120			& ' illegal value for kp=',kp
121			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
122			& SQUEEZE_RIGHT, myThid )
123			WRITE(msgBuf,'(A,1P5E20.13)')
124			& ' FLT_TRAJ: ', (tmp(ii),ii=1,5)
125			CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
126			& SQUEEZE_RIGHT, myThid )
127			c CALL PRINT_ERROR( msgBuf, myThid )
128			c STOP 'ABNORMAL END: S/R FLT_TRAJ'
129	jmc	1.4	C-- jmc: not sure if this is right but added to avoid Pb in FLT_BILINEAR:
130			kp = MIN( MAX(kp,1), Nr)
131	jmc	1.2	ENDIF
132
133	jmc	1.4	CALL FLT_BILINEAR (xx,yy,uu,uVel, kp,1,bi,bj,myThid)
134			CALL FLT_BILINEAR (xx,yy,vv,vVel, kp,2,bi,bj,myThid)
135			CALL FLT_BILINEAR2D(xx,yy,pp,etaN, 0,bi,bj,myThid)
136			CALL FLT_BILINEAR (xx,yy,tt,theta, kp,0,bi,bj,myThid)
137			CALL FLT_BILINEAR (xx,yy,ss,salt, kp,0,bi,bj,myThid)
138	jmc	1.2
139			tmp(6) = uu
140			tmp(7) = vv
141			tmp(8) = tt
142			tmp(9) = ss
143			tmp(10) = pp
144			ELSE
145			tmp(6) = flt_nan
146			tmp(7) = flt_nan
147			tmp(8) = flt_nan
148			tmp(9) = flt_nan
149			tmp(10) = flt_nan
150			ENDIF
151
152			C (3) WRITE float positions into file
153	jmc	1.5	irecord = npart_read+ip+1
154			IF ( ip.NE.npart_tile(bi,bj) ) irecord = -irecord
155			CALL MDS_WRITEVEC_LOC( fn, precFloat64, ioUnit,
156			& 'RL', imax, tmp,
157			& bi,bj,irecord, myIter, myThid )
158	adcroft	1.1
159	jmc	1.2	ENDDO
160	jmc	1.5	CLOSE( ioUnit )
161	adcroft	1.1
162			ENDDO
163	jmc	1.4	ENDDO
164	adcroft	1.1
165	jmc	1.2	RETURN
166			END