/[MITgcm]/MITgcm/pkg/flt/flt_traj.F
ViewVC logotype

Diff of /MITgcm/pkg/flt/flt_traj.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.12 by jmc, Mon Dec 27 19:21:23 2010 UTC revision 1.13 by jmc, Wed Aug 31 21:41:55 2011 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "FLT_OPTIONS.h"  #include "FLT_OPTIONS.h"
5    
6    CBOP 0
7    C !ROUTINE: FLT_TRAJ
8    
9    C !INTERFACE:
10        SUBROUTINE FLT_TRAJ (        SUBROUTINE FLT_TRAJ (
11       I                      myTime, myIter, myThid )       I                      myTime, myIter, myThid )
12    
13  C     ==================================================================  C     !DESCRIPTION:
14  C     SUBROUTINE FLT_TRAJ  C     *==========================================================*
15  C     ==================================================================  C     | SUBROUTINE FLT_TRAJ
16  C     o This routine samples the model state at float position every  C     | o This routine samples the model state at float position
17  C       flt_int_traj time steps and writes output.  C     |   every flt_int_traj time steps and writes output.
18  C     ==================================================================  C     *==========================================================*
19    
20  C     !USES:  C     !USES:
21        IMPLICIT NONE        IMPLICIT NONE
   
22  C     == global variables ==  C     == global variables ==
23  #include "SIZE.h"  #include "SIZE.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 24  C     == global variables == Line 26  C     == global variables ==
26  #include "DYNVARS.h"  #include "DYNVARS.h"
27  #include "FLT_SIZE.h"  #include "FLT_SIZE.h"
28  #include "FLT.h"  #include "FLT.h"
29    #include "FLT_BUFF.h"
30  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
31  #include "W2_EXCH2_SIZE.h"  #include "W2_EXCH2_SIZE.h"
32  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
33  #endif  #endif
34    
35  C     == routine arguments ==  C     !INPUT PARAMETERS:
36    C     myTime :: current time in simulation
37    C     myIter :: current iteration number
38    C     myThid :: my Thread Id number
39        _RL myTime        _RL myTime
40        INTEGER myIter, myThid        INTEGER myIter, myThid
41    
42  C     === Functions ==  C     !FUNCTIONS:
43        _RL FLT_MAP_K2R        _RL FLT_MAP_K2R
44        EXTERNAL FLT_MAP_K2R        EXTERNAL FLT_MAP_K2R
45    
46  C     == local variables ==  C     !LOCAL VARIABLES:
47        INTEGER bi, bj, imax        INTEGER bi, bj, nFlds
       PARAMETER (imax=13)  
48        INTEGER ip, kp, ii        INTEGER ip, kp, ii
49        _RL ix, jy, i0x, j0y, xx, yy, zz        _RL ix, jy, i0x, j0y, xx, yy, zz
50        _RL uu, vv, tt, ss, pp        _RL uu, vv, tt, ss, pp
51    
52          INTEGER imax
53          PARAMETER (imax=13)
54        _RL tmp(imax)        _RL tmp(imax)
55        _RL npart_read,npart_times        _RL npart_read, npart_times
56        _RS dummyRS(1)        _RS dummyRS(1)
57        INTEGER fp, ioUnit, irecord        INTEGER fp, ioUnit, irecord
58        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
# Line 53  C     == local variables == Line 60  C     == local variables ==
60  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
61        INTEGER nT        INTEGER nT
62  #endif  #endif
63    CEOP
64    
65  C     == end of interface ==  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
66    
67        fn = 'float_trajectories'  C--   set number of fields to write
68        fp = writeBinaryPrec        nFlds = 0
69          IF ( flt_selectTrajOutp.GE.1 ) nFlds = nFlds + 8
70          IF ( flt_selectTrajOutp.GE.2 ) nFlds = nFlds + 5
71    
72    C--   check buffer size
73          IF ( nFlds.GT.fltBufDim ) THEN
74             WRITE(msgBuf,'(3(A,I4))') ' FLT_TRAJ: fltBufDim=', fltBufDim,
75         &                             ' too small (<', nFlds, ' )'
76             CALL PRINT_ERROR( msgBuf, myThid )
77             WRITE(msgBuf,'(2A)')     ' FLT_TRAJ: => increase fltBufDim',
78         &                            ' in "FLT_SIZE.h" & recompile'
79             CALL PRINT_ERROR( msgBuf, myThid )
80             CALL ALL_PROC_DIE( myThid )
81    c        IF (myThid.EQ.1) STOP 'ABNORMAL END: S/R FLT_TRAJ'
82    c        _BARRIER
83             STOP 'ABNORMAL END: S/R FLT_TRAJ'
84          ENDIF
85    
86        DO bj=myByLo(myThid),myByHi(myThid)        IF ( myIter.EQ.nIter0 .OR. flt_selectTrajOutp.LE.0 ) RETURN
        DO bi=myBxLo(myThid),myBxHi(myThid)  
87    
88  C (1) read actual number floats from file (if exists)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89           ioUnit = -2  C--   Calculate position + other fields at float position and fill up IO-buffer
          CALL MDS_READVEC_LOC(  fn, fp, ioUnit,  
      &                          'RL', imax, tmp, dummyRS,  
      &                           bi, bj, 1, myThid )  
          IF ( ioUnit.GT. 0 ) THEN  
             npart_read  = tmp(1)  
             npart_times = tmp(5)  
 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  
90    
91  C the standard routine mds_writevec_loc can be used here        DO bj=myByLo(myThid),myByHi(myThid)
92  C (2) WRITE new actual number floats and time axis into file         DO bi=myBxLo(myThid),myBxHi(myThid)
 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 ii=7,imax  
             tmp(ii) = 0.  
          ENDDO  
   
          ioUnit = -1  
          CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit,  
      &                          'RL', imax, tmp, dummyRS,  
      &                           bi,bj,-1, myIter, myThid )  
93    
94  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
95           nT = W2_myTileList(bi,bj)           nT = W2_myTileList(bi,bj)
# Line 118  C total number of floats Line 107  C total number of floats
107       I                               ix, jy, bi,bj, myThid )       I                               ix, jy, bi,bj, myThid )
108              zz = FLT_MAP_K2R( kpart(ip,bi,bj),bi,bj,myThid )              zz = FLT_MAP_K2R( kpart(ip,bi,bj),bi,bj,myThid )
109              kp = NINT(kpart(ip,bi,bj))              kp = NINT(kpart(ip,bi,bj))
110              tmp(1)  = npart(ip,bi,bj)              tmp(1) = npart(ip,bi,bj)
111              tmp(2)  = myTime              tmp(2) = myTime
112              tmp(3)  = xx              tmp(3) = xx
113              tmp(4)  = yy              tmp(4) = yy
114              tmp(5)  = zz              tmp(5) = zz
115              tmp(6)  = ix + i0x              tmp(6) = ix + i0x
116              tmp(7)  = jy + j0y              tmp(7) = jy + j0y
117              tmp(8)  = kpart(ip,bi,bj)              tmp(8) = kpart(ip,bi,bj)
118    
119              IF ( ( myTime.GE.tstart(ip,bi,bj)) .AND.              IF ( ( flt_selectTrajOutp.GE.2 )   .AND.
120         &           ( myTime.GE.tstart(ip,bi,bj)) .AND.
121       &           ( tend(ip,bi,bj).EQ.-1. .OR. myTime.LE.tend(ip,bi,bj))       &           ( tend(ip,bi,bj).EQ.-1. .OR. myTime.LE.tend(ip,bi,bj))
122       &         ) THEN       &         ) THEN
   
123                IF ( kp.LT.1 .OR. kp.GT.Nr ) THEN                IF ( kp.LT.1 .OR. kp.GT.Nr ) THEN
124                  WRITE(msgBuf,'(2A,I8)') '** WARNING ** FLT_TRAJ: ',                  WRITE(msgBuf,'(2A,I8)') '** WARNING ** FLT_TRAJ: ',
125       &            ' illegal value for kp=',kp       &            ' illegal value for kp=',kp
126                  CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,                  CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
127       &                              SQUEEZE_RIGHT, myThid )       &                              SQUEEZE_RIGHT, myThid )
128                  WRITE(msgBuf,'(A,1P5E20.13)')                  WRITE(msgBuf,'(A,1P5E20.13)')
129       &            ' FLT_TRAJ: ', (tmp(ii),ii=1,5)       &            ' FLT_TRAJ: ', (flt_io_buff(ii,ip,bi,bj),ii=1,5)
130                  CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,                  CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
131       &                              SQUEEZE_RIGHT, myThid )       &                              SQUEEZE_RIGHT, myThid )
132  c               CALL PRINT_ERROR( msgBuf, myThid )  c               CALL PRINT_ERROR( msgBuf, myThid )
# Line 145  c               STOP 'ABNORMAL END: S/R Line 134  c               STOP 'ABNORMAL END: S/R
134  C-- jmc: not sure if this is right but added to avoid Pb in FLT_BILINEAR:  C-- jmc: not sure if this is right but added to avoid Pb in FLT_BILINEAR:
135                  kp = MIN( MAX(kp,1), Nr)                  kp = MIN( MAX(kp,1), Nr)
136                ENDIF                ENDIF
   
137                CALL FLT_BILINEAR  (ix,jy,uu,uVel,  kp,1,bi,bj,myThid)                CALL FLT_BILINEAR  (ix,jy,uu,uVel,  kp,1,bi,bj,myThid)
138                CALL FLT_BILINEAR  (ix,jy,vv,vVel,  kp,2,bi,bj,myThid)                CALL FLT_BILINEAR  (ix,jy,vv,vVel,  kp,2,bi,bj,myThid)
139                CALL FLT_BILINEAR2D(ix,jy,pp,etaN,     0,bi,bj,myThid)                CALL FLT_BILINEAR2D(ix,jy,pp,etaN,     0,bi,bj,myThid)
140                CALL FLT_BILINEAR  (ix,jy,tt,theta, kp,0,bi,bj,myThid)                CALL FLT_BILINEAR  (ix,jy,tt,theta, kp,0,bi,bj,myThid)
141                CALL FLT_BILINEAR  (ix,jy,ss,salt,  kp,0,bi,bj,myThid)                CALL FLT_BILINEAR  (ix,jy,ss,salt,  kp,0,bi,bj,myThid)
   
142                tmp( 9) = pp                tmp( 9) = pp
143                tmp(10) = uu                tmp(10) = uu
144                tmp(11) = vv                tmp(11) = vv
145                tmp(12) = tt                tmp(12) = tt
146                tmp(13) = ss                tmp(13) = ss
147              ELSE              ELSEIF ( flt_selectTrajOutp.GE.2 ) THEN
148                tmp( 9) = flt_nan                tmp( 9) = flt_nan
149                tmp(10) = flt_nan                tmp(10) = flt_nan
150                tmp(11) = flt_nan                tmp(11) = flt_nan
# Line 165  C-- jmc: not sure if this is right but a Line 152  C-- jmc: not sure if this is right but a
152                tmp(13) = flt_nan                tmp(13) = flt_nan
153              ENDIF              ENDIF
154    
155  C (3) WRITE float positions into file              DO ii=1,nFlds
156                  flt_io_buff(ii,ip,bi,bj) = tmp(ii)
157                ENDDO
158    
159             ENDDO
160    
161           ENDDO
162          ENDDO
163    
164    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
165    
166    C--   Write shared buffer to file
167    
168          _BARRIER
169          _BEGIN_MASTER(myThid)
170    
171          fn = 'float_trajectories'
172          fp = writeBinaryPrec
173    
174          DO bj=1,nSy
175           DO bi=1,nSx
176    
177    C (1) read actual number floats from file (if exists)
178             ioUnit = -2
179             CALL MDS_READVEC_LOC(  fn, fp, ioUnit, 'RL', nFlds,
180         O                          tmp, dummyRS,
181         I                          bi, bj, 1, myThid )
182             IF ( ioUnit.GT. 0 ) THEN
183                npart_read  = tmp(1)
184                npart_times = tmp(5)
185                ii = NINT(tmp(7))
186    C-       for backward compatibility with old trajectory files:
187                IF ( ii.EQ.0 ) ii = 13
188                IF ( ii.NE.nFlds ) THEN
189                  WRITE(msgBuf,'(A,I4,A)')
190         &            'FLT_TRAJ: nFlds=', nFlds,' different from'
191                  CALL PRINT_ERROR( msgBuf, myThid )
192                  WRITE(msgBuf,'(3A,I4,A)')
193         &            'previous file (',fn(1:18),') value =',ii
194                  CALL PRINT_ERROR( msgBuf, myThid )
195                  CALL ALL_PROC_DIE( myThid )
196                  STOP 'ABNORMAL END: S/R FLT_TRAJ'
197                ENDIF
198    C-       close the read-unit (safer to use a different unit for writing)
199                CLOSE( ioUnit )
200             ELSE
201                npart_read  = 0.
202                npart_times = 0.
203                tmp(2)      = myTime
204             ENDIF
205    
206    C (2) write new actual number floats and time axis into file
207    C-    the standard routine mds_writevec_loc can be used here
208    
209    C     total number of records in this file
210             tmp(1) = DBLE(npart_tile(bi,bj))+npart_read
211    C     first time of writing floats (do not change when written)
212    c        tmp(2) = tmp(2)
213    C     current time
214             tmp(3) = myTime
215    C     timestep
216             tmp(4) = flt_int_traj
217    C     total number of timesteps
218             tmp(5) = npart_times + 1.
219    C     total number of floats
220             tmp(6) = max_npart
221    C     total number of fields
222             tmp(7) = nFlds
223             DO ii=8,nFlds
224               tmp(ii) = 0.
225             ENDDO
226             ioUnit = -1
227             CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit, 'RL', nFlds,
228         &                          tmp, dummyRS,
229         &                          bi, bj, -1, myIter, myThid )
230    
231             DO ip=1,npart_tile(bi,bj)
232    C (3) write float positions into file
233              irecord = npart_read+ip+1              irecord = npart_read+ip+1
234              IF ( ip.NE.npart_tile(bi,bj) ) irecord = -irecord              IF ( ip.NE.npart_tile(bi,bj) ) irecord = -irecord
235              CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit,              CALL MDS_WRITEVEC_LOC( fn, fp, ioUnit, 'RL', nFlds,
236       &                            'RL', imax, tmp, dummyRS,       I                             flt_io_buff(1,ip,bi,bj), dummyRS,
237       &                             bi,bj,irecord, myIter, myThid )       I                             bi, bj, irecord, myIter, myThid )
   
238           ENDDO           ENDDO
239           CLOSE( ioUnit )           CLOSE( ioUnit )
240    
241         ENDDO         ENDDO
242        ENDDO        ENDDO
243    
244          _END_MASTER(myThid)
245          _BARRIER
246    
247        RETURN        RETURN
248        END        END
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22