/[MITgcm]/MITgcm/pkg/mdsio/mdsio_readfield.F
ViewVC logotype

Contents of /MITgcm/pkg/mdsio/mdsio_readfield.F

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


Revision 1.13 - (show annotations) (download)
Tue Apr 6 00:25:56 2004 UTC (20 years, 2 months ago) by dimitri
Branch: MAIN
Changes since 1.12: +57 -15 lines
o extending useSingleCpuIO option to work with new exch2 I/O format
  - old-style, missing-tile I/O is still accessible by defining CPP
    option MISSING_TILE_IO in pkg/mdsio/MDSIO_OPTIONS.h

1 C $Header: /u/gcmpack/MITgcm/pkg/mdsio/mdsio_readfield.F,v 1.12 2004/03/16 06:38:00 dimitri Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "MDSIO_OPTIONS.h"
6
7 SUBROUTINE MDSREADFIELD(
8 I fName,
9 I filePrec,
10 I arrType,
11 I nNz,
12 O arr,
13 I irecord,
14 I myThid )
15 C
16 C Arguments:
17 C
18 C fName string base name for file to read
19 C filePrec integer number of bits per word in file (32 or 64)
20 C arrType char(2) declaration of "arr": either "RS" or "RL"
21 C nNz integer size of third dimension: normally either 1 or Nr
22 C arr RS/RL array to read into, arr(:,:,nNz,:,:)
23 C irecord integer record number to read
24 C myThid integer thread identifier
25 C
26 C MDSREADFIELD first checks to see if the file "fName" exists, then
27 C if the file "fName.data" exists and finally the tiled files of the
28 C form "fName.xxx.yyy.data" exist. Currently, the meta-files are not
29 C read because it is difficult to parse files in fortran.
30 C The precision of the file is decsribed by filePrec, set either
31 C to floatPrec32 or floatPrec64. The precision or declaration of
32 C the array argument must be consistently described by the char*(2)
33 C string arrType, either "RS" or "RL". nNz allows for both 2-D and
34 C 3-D arrays to be handled. nNz=1 implies a 2-D model field and
35 C nNz=Nr implies a 3-D model field. irecord is the record number
36 C to be read and must be >= 1. The file data is stored in
37 C arr *but* the overlaps are *not* updated. ie. An exchange must
38 C be called. This is because the routine is sometimes called from
39 C within a MASTER_THID region.
40 C
41 C Created: 03/16/99 adcroft@mit.edu
42
43 implicit none
44 C Global variables / common blocks
45 #include "SIZE.h"
46 #include "EEPARAMS.h"
47 #include "PARAMS.h"
48 #include "EESUPPORT.h"
49 #ifdef ALLOW_EXCH2
50 #include "W2_EXCH2_TOPOLOGY.h"
51 #include "W2_EXCH2_PARAMS.h"
52 #endif /* ALLOW_EXCH2 */
53
54 C Routine arguments
55 character*(*) fName
56 integer filePrec
57 character*(2) arrType
58 integer nNz
59 Real arr(*)
60 integer irecord
61 integer myThid
62 C Functions
63 integer ILNBLNK
64 integer MDS_RECLEN
65 C Local variables
66 character*(80) dataFName,pfName
67 integer iG,jG,irec,bi,bj,i,j,k,dUnit,IL,pIL
68 logical exst
69 Real*4 r4seg(sNx)
70 Real*8 r8seg(sNx)
71 logical globalFile,fileIsOpen
72 integer x_size,y_size,iG_IO,jG_IO,length_of_rec
73 #if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
74 PARAMETER ( x_size = exch2_domain_nxt * sNx )
75 PARAMETER ( y_size = exch2_domain_nyt * sNy )
76 #else
77 PARAMETER ( x_size = Nx )
78 PARAMETER ( y_size = Ny )
79 #endif
80 Real*4 xy_buffer_r4(x_size,y_size)
81 Real*8 xy_buffer_r8(x_size,y_size)
82 character*(max_len_mbuf) msgbuf
83 Real*8 global (Nx,Ny)
84 _RL local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85 #ifdef ALLOW_EXCH2
86 integer domainHeight,domainLength,tby,tgx,tny,tnx,tn
87 #endif /* ALLOW_EXCH2 */
88 C ------------------------------------------------------------------
89
90 C Only do I/O if I am the master thread
91 _BEGIN_MASTER( myThid )
92
93 C Record number must be >= 1
94 if (irecord .LT. 1) then
95 write(msgbuf,'(a,i9.8)')
96 & ' MDSREADFIELD: argument irecord = ',irecord
97 call print_message( msgbuf, standardmessageunit,
98 & SQUEEZE_RIGHT , mythid)
99 write(msgbuf,'(a)')
100 & ' MDSREADFIELD: Invalid value for irecord'
101 call print_error( msgbuf, mythid )
102 stop 'ABNORMAL END: S/R MDSREADFIELD'
103 endif
104
105 C Assume nothing
106 globalFile = .FALSE.
107 fileIsOpen = .FALSE.
108 IL = ILNBLNK( fName )
109 pIL = ILNBLNK( mdsioLocalDir )
110
111 C Assign special directory
112 if ( mdsioLocalDir .NE. ' ' ) then
113 write(pFname(1:80),'(2a)') mdsioLocalDir(1:pIL), fName(1:IL)
114 else
115 pFname= fName
116 endif
117 pIL=ILNBLNK( pfName )
118
119 C Assign a free unit number as the I/O channel for this routine
120 call MDSFINDUNIT( dUnit, mythid )
121
122 C Check first for global file with simple name (ie. fName)
123 dataFName = fName
124 inquire( file=dataFname, exist=exst )
125 if (exst) then
126 if ( debugLevel .GE. debLevA ) then
127 write(msgbuf,'(a,a)')
128 & ' MDSREADFIELD: opening global file: ',dataFName
129 call print_message( msgbuf, standardmessageunit,
130 & SQUEEZE_RIGHT , mythid)
131 endif
132 globalFile = .TRUE.
133 endif
134
135 C If negative check for global file with MDS name (ie. fName.data)
136 if (.NOT. globalFile) then
137 write(dataFname(1:80),'(2a)') fName(1:IL),'.data'
138 inquire( file=dataFname, exist=exst )
139 if (exst) then
140 if ( debugLevel .GE. debLevA ) then
141 write(msgbuf,'(a,a)')
142 & ' MDSREADFIELD: opening global file: ',dataFName
143 call print_message( msgbuf, standardmessageunit,
144 & SQUEEZE_RIGHT , mythid)
145 endif
146 globalFile = .TRUE.
147 endif
148 endif
149
150 if ( .not. ( globalFile .and. useSingleCPUIO ) ) then
151
152 C If we are reading from a global file then we open it here
153 if (globalFile) then
154 length_of_rec=MDS_RECLEN( filePrec, sNx, mythid )
155 open( dUnit, file=dataFName, status='old',
156 & access='direct', recl=length_of_rec )
157 fileIsOpen=.TRUE.
158 endif
159
160 #ifdef ALLOW_EXCH2
161
162 if (globalFile) then
163 domainLength = exch2_domain_nxt
164 domainHeight = exch2_domain_nyt
165
166 cafe write(*,fmt='(1X,A,I3,A,I3)') 'L=', domainlength,
167 cafe & 'H=', domainheight
168 C Loop over all tiles
169 do bj=1,nSy
170 do bi=1,nSx
171 C If we are reading from a tiled MDS file then we open each one here
172 cafe if (.NOT. globalFile) then
173 cafe write(msgbuf,'(a)')
174 cafe & ' MDSREADFIELD: non-global input files not
175 cafe & implemented with exch2'
176 cafe call print_error( msgbuf, mythid )
177 cafe stop 'ABNORMAL END: S/R MDSREADFIELD'
178 cafe endif
179 tn = W2_myTileList(bi)
180 tby = exch2_tbasey(tn)
181 tgx = exch2_txglobalo(tn)
182 tny = exch2_tny(tn)
183 tnx = exch2_tnx(tn)
184 if (fileIsOpen) then
185 do k=1,nNz
186 do j=1,tNy
187 cafe write(*,fmt='(1X,A,I3,A,I3,A,I3,A,I3,A,I3,A,I3)') 'tby=', tby,
188 cafe & ', tgx=', tgx,
189 cafe & ', tnx=',tnx, ', tny=', tny, ', j=',j,', tn=',tn
190
191 irec = domainLength*tby + (tgx-1)/tnx + 1 +
192 & domainLength*(j-1) +
193 & domainLength*domainHeight*tny*(k-1) +
194 & domainLength*domainHeight*tny*nNz*(irecord-1)
195 cafe write(*,fmt='(1X,A,I6,A,I3)') 'record = ',irec,',thingy=',
196 cafe & (tgx-1)/tnx
197 cafe write(*,fmt='(1X,A,I6)') 'record = ',irec
198
199 if (filePrec .eq. precFloat32) then
200 read(dUnit,rec=irec) r4seg
201 #ifdef _BYTESWAPIO
202 call MDS_BYTESWAPR4( sNx, r4seg )
203 #endif
204 if (arrType .eq. 'RS') then
205 call MDS_SEG4toRS( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
206 elseif (arrType .eq. 'RL') then
207 call MDS_SEG4toRL( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
208 else
209 write(msgbuf,'(a)')
210 & ' MDSREADFIELD: illegal value for arrType'
211 call print_error( msgbuf, mythid )
212 stop 'ABNORMAL END: S/R MDSREADFIELD'
213 endif
214 elseif (filePrec .eq. precFloat64) then
215 read(dUnit,rec=irec) r8seg
216 #ifdef _BYTESWAPIO
217 call MDS_BYTESWAPR8( sNx, r8seg )
218 #endif
219 if (arrType .eq. 'RS') then
220 call MDS_SEG8toRS( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
221 elseif (arrType .eq. 'RL') then
222 call MDS_SEG8toRL( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
223 else
224 write(msgbuf,'(a)')
225 & ' MDSREADFIELD: illegal value for arrType'
226 call print_error( msgbuf, mythid )
227 stop 'ABNORMAL END: S/R MDSREADFIELD'
228 endif
229 else
230 write(msgbuf,'(a)')
231 & ' MDSREADFIELD: illegal value for filePrec'
232 call print_error( msgbuf, mythid )
233 stop 'ABNORMAL END: S/R MDSREADFIELD'
234 endif
235 C End of j loop
236 enddo
237 C End of k loop
238 enddo
239 if (.NOT. globalFile) then
240 close( dUnit )
241 fileIsOpen = .FALSE.
242 endif
243 endif
244 C End of bi,bj loops
245 enddo
246 enddo
247 else ! if globalFile
248 c#else /* .not. ALLOW_EXCH2 */
249 #endif /* ALLOW_EXCH2 */
250
251
252 C Loop over all tiles
253 do bj=1,nSy
254 do bi=1,nSx
255 C If we are reading from a tiled MDS file then we open each one here
256 if (.NOT. globalFile) then
257 iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
258 jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
259 write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)')
260 & pfName(1:pIL),'.',iG,'.',jG,'.data'
261 inquire( file=dataFname, exist=exst )
262 C Of course, we only open the file if the tile is "active"
263 C (This is a place-holder for the active/passive mechanism
264 if (exst) then
265 if ( debugLevel .GE. debLevA ) then
266 write(msgbuf,'(a,a)')
267 & ' MDSREADFIELD: opening file: ',dataFName
268 call print_message( msgbuf, standardmessageunit,
269 & SQUEEZE_RIGHT , mythid)
270 endif
271 length_of_rec=MDS_RECLEN( filePrec, sNx, mythid )
272 open( dUnit, file=dataFName, status='old',
273 & access='direct', recl=length_of_rec )
274 fileIsOpen=.TRUE.
275 else
276 fileIsOpen=.FALSE.
277 write(msgbuf,'(3a)')
278 & ' MDSREADFIELD: filename: ',dataFName, pfName
279 call print_message( msgbuf, standardmessageunit,
280 & SQUEEZE_RIGHT , mythid)
281 write(msgbuf,'(a)')
282 & ' MDSREADFIELD: File does not exist'
283 call print_error( msgbuf, mythid )
284 stop 'ABNORMAL END: S/R MDSREADFIELD'
285 endif
286 endif
287
288 if (fileIsOpen) then
289 do k=1,nNz
290 do j=1,sNy
291 if (globalFile) then
292 iG = myXGlobalLo-1 + (bi-1)*sNx
293 jG = myYGlobalLo-1 + (bj-1)*sNy
294 irec=1 + INT(iG/sNx) + nSx*nPx*(jG+j-1) + nSx*nPx*Ny*(k-1)
295 & + nSx*nPx*Ny*nNz*(irecord-1)
296 else
297 iG = 0
298 jG = 0
299 irec=j + sNy*(k-1) + sNy*nNz*(irecord-1)
300 endif
301 if (filePrec .eq. precFloat32) then
302 read(dUnit,rec=irec) r4seg
303 #ifdef _BYTESWAPIO
304 call MDS_BYTESWAPR4( sNx, r4seg )
305 #endif
306 if (arrType .eq. 'RS') then
307 call MDS_SEG4toRS( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
308 elseif (arrType .eq. 'RL') then
309 call MDS_SEG4toRL( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
310 else
311 write(msgbuf,'(a)')
312 & ' MDSREADFIELD: illegal value for arrType'
313 call print_error( msgbuf, mythid )
314 stop 'ABNORMAL END: S/R MDSREADFIELD'
315 endif
316 elseif (filePrec .eq. precFloat64) then
317 read(dUnit,rec=irec) r8seg
318 #ifdef _BYTESWAPIO
319 call MDS_BYTESWAPR8( sNx, r8seg )
320 #endif
321 if (arrType .eq. 'RS') then
322 call MDS_SEG8toRS( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
323 elseif (arrType .eq. 'RL') then
324 call MDS_SEG8toRL( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
325 else
326 write(msgbuf,'(a)')
327 & ' MDSREADFIELD: illegal value for arrType'
328 call print_error( msgbuf, mythid )
329 stop 'ABNORMAL END: S/R MDSREADFIELD'
330 endif
331 else
332 write(msgbuf,'(a)')
333 & ' MDSREADFIELD: illegal value for filePrec'
334 call print_error( msgbuf, mythid )
335 stop 'ABNORMAL END: S/R MDSREADFIELD'
336 endif
337 C End of j loop
338 enddo
339 C End of k loop
340 enddo
341 if (.NOT. globalFile) then
342 close( dUnit )
343 fileIsOpen = .FALSE.
344 endif
345 endif
346 C End of bi,bj loops
347 enddo
348 enddo
349
350 #ifdef ALLOW_EXCH2
351 endif ! globalFile
352
353 #endif /* ALLOW_EXCH2 */
354
355 C If global file was opened then close it
356 if (fileIsOpen .AND. globalFile) then
357 close( dUnit )
358 fileIsOpen = .FALSE.
359 endif
360
361 endif
362 c endif ( .not. ( globalFile .and. useSingleCPUIO ) )
363
364 _END_MASTER( myThid )
365
366 if ( globalFile .and. useSingleCPUIO ) then
367
368 C master thread of process 0, only, opens a global file
369 _BEGIN_MASTER( myThid )
370 #ifdef ALLOW_USE_MPI
371 IF( mpiMyId .EQ. 0 ) THEN
372 #else
373 IF ( .TRUE. ) THEN
374 #endif /* ALLOW_USE_MPI */
375 length_of_rec=MDS_RECLEN( filePrec, xy_size, mythid )
376 open( dUnit, file=dataFName, status='old',
377 & access='direct', recl=length_of_rec )
378 ENDIF
379 _END_MASTER( myThid )
380
381 DO k=1,nNz
382
383 _BEGIN_MASTER( myThid )
384 #ifdef ALLOW_USE_MPI
385 IF( mpiMyId .EQ. 0 ) THEN
386 #else
387 IF ( .TRUE. ) THEN
388 #endif /* ALLOW_USE_MPI */
389 irec = k+nNz*(irecord-1)
390 if (filePrec .eq. precFloat32) then
391 read(dUnit,rec=irec) xy_buffer_r4
392 #ifdef _BYTESWAPIO
393 call MDS_BYTESWAPR4( xy_size, xy_buffer_r4 )
394 #endif
395 #if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
396 bj=1
397 DO npe=1,nPx*nPy
398 DO bi=1,nSx
399 DO J=1,sNy
400 DO I=1,sNx
401 iG=mpi_myXGlobalLo(npe)-1+(bi-1)*sNx+i
402 jG=mpi_myYGlobalLo(npe)-1+(bj-1)*sNy+j
403 iG_IO=exch2_txglobalo(W2_mpi_myTileList(npe,bi))+i-1
404 jG_IO=exch2_tyglobalo(W2_mpi_myTileList(npe,bi))+j-1
405 global(iG,jG) = xy_buffer_r4(iG_IO,jG_IO)
406 ENDDO
407 ENDDO
408 ENDDO
409 ENDDO
410 #else /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
411 DO J=1,Ny
412 DO I=1,Nx
413 global(I,J) = xy_buffer_r4(I+(J-1)*Ny)
414 ENDDO
415 ENDDO
416 #endif /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
417 elseif (filePrec .eq. precFloat64) then
418 read(dUnit,rec=irec) xy_buffer_r8
419 #ifdef _BYTESWAPIO
420 call MDS_BYTESWAPR8( xy_size, xy_buffer_r8 )
421 #endif
422 #if defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO)
423 bj=1
424 DO npe=1,nPx*nPy
425 DO bi=1,nSx
426 DO J=1,sNy
427 DO I=1,sNx
428 iG=mpi_myXGlobalLo(npe)-1+(bi-1)*sNx+i
429 jG=mpi_myYGlobalLo(npe)-1+(bj-1)*sNy+j
430 iG_IO=exch2_txglobalo(W2_mpi_myTileList(npe,bi))+i-1
431 jG_IO=exch2_tyglobalo(W2_mpi_myTileList(npe,bi))+j-1
432 global(iG,jG) = xy_buffer_r8(iG_IO,jG_IO)
433 ENDDO
434 ENDDO
435 ENDDO
436 ENDDO
437 #else /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
438 DO J=1,Ny
439 DO I=1,Nx
440 global(I,J) = xy_buffer_r8(I+(J-1)*Ny)
441 ENDDO
442 ENDDO
443 #endif /* defined(ALLOW_EXCH2) && !defined(MISSING_TILE_IO) */
444 else
445 write(msgbuf,'(a)')
446 & ' MDSREADFIELD: illegal value for filePrec'
447 call print_error( msgbuf, mythid )
448 stop 'ABNORMAL END: S/R MDSREADFIELD'
449 endif
450 ENDIF
451 _END_MASTER( myThid )
452 CALL SCATTER_2D(global,local,mythid)
453 if (arrType .eq. 'RS') then
454 call PASStoRS( local,arr,k,nNz,mythid )
455 elseif (arrType .eq. 'RL') then
456 call PASStoRL( local,arr,k,nNz,mythid )
457 else
458 write(msgbuf,'(a)')
459 & ' MDSREADFIELD: illegal value for arrType'
460 call print_error( msgbuf, mythid )
461 stop 'ABNORMAL END: S/R MDSREADFIELD'
462 endif
463
464 ENDDO
465 c ENDDO k=1,nNz
466
467 _BEGIN_MASTER( myThid )
468 close( dUnit )
469 _END_MASTER( myThid )
470
471 endif
472 c endif ( globalFile .and. useSingleCPUIO )
473
474 return
475 end
476
477
478 C ==================================================================
479
480 subroutine passToRS(local,arr,k,nNz,mythid)
481 implicit none
482 #include "EEPARAMS.h"
483 #include "SIZE.h"
484 _RL local(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
485 integer i,j,k,bi,bj,nNz
486 _RS arr(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nNz,nSx,nSy)
487 integer mythid
488 DO bj = myByLo(myThid), myByHi(myThid)
489 DO bi = myBxLo(myThid), myBxHi(myThid)
490 DO J=1-Oly,sNy+Oly
491 DO I=1-Olx,sNx+Olx
492 arr(I,J,k,bi,bj) = local(I,J,bi,bj)
493 ENDDO
494 ENDDO
495 ENDDO
496 ENDDO
497 return
498 end
499
500 subroutine passToRL(local,arr,k,nNz,mythid)
501 implicit none
502 #include "EEPARAMS.h"
503 #include "SIZE.h"
504 _RL local(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
505 integer i,j,k,bi,bj,nNz
506 _RL arr(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nNz,nSx,nSy)
507 integer mythid
508 DO bj = myByLo(myThid), myByHi(myThid)
509 DO bi = myBxLo(myThid), myBxHi(myThid)
510 DO J=1-Oly,sNy+Oly
511 DO I=1-Olx,sNx+Olx
512 arr(I,J,k,bi,bj) = local(I,J,bi,bj)
513 ENDDO
514 ENDDO
515 ENDDO
516 ENDDO
517 return
518 end

  ViewVC Help
Powered by ViewVC 1.1.22