/[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.11 - (show annotations) (download)
Mon Mar 15 01:37:22 2004 UTC (20 years, 3 months ago) by cnh
Branch: MAIN
Changes since 1.10: +2 -1 lines
o Fixes for ALLOW_EXCH2 as a packages_conf entry.

o Validation configurations for testing a variety of cube tilings including
  tiles that are left out over land. This second set of items needs
  mods to genmake2 to be automatically included in regression tests - until
  then it can be tested by hand. The file
  verification/global_ocean.cs32x15/code_alt/README
  shows how to do this manually with genmake2

1 C $Header: /u/gcmpack/MITgcm/pkg/mdsio/mdsio_readfield.F,v 1.10 2004/03/08 21:20:02 adcroft 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,j,k,dUnit,IL,pIL
68 logical exst
69 Real*4 r4seg(sNx)
70 Real*8 r8seg(sNx)
71 logical globalFile,fileIsOpen
72 integer length_of_rec
73 character*(max_len_mbuf) msgbuf
74 integer i
75 Real*4 global_r4(Nx,Ny)
76 Real*8 global (Nx,Ny)
77 _RL local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78 #ifdef ALLOW_EXCH2
79 integer domainHeight,domainLength,tby,tgx,tny,tnx,tn
80 #endif /* ALLOW_EXCH2 */
81 C ------------------------------------------------------------------
82
83 C Only do I/O if I am the master thread
84 _BEGIN_MASTER( myThid )
85
86 C Record number must be >= 1
87 if (irecord .LT. 1) then
88 write(msgbuf,'(a,i9.8)')
89 & ' MDSREADFIELD: argument irecord = ',irecord
90 call print_message( msgbuf, standardmessageunit,
91 & SQUEEZE_RIGHT , mythid)
92 write(msgbuf,'(a)')
93 & ' MDSREADFIELD: Invalid value for irecord'
94 call print_error( msgbuf, mythid )
95 stop 'ABNORMAL END: S/R MDSREADFIELD'
96 endif
97
98 C Assume nothing
99 globalFile = .FALSE.
100 fileIsOpen = .FALSE.
101 IL = ILNBLNK( fName )
102 pIL = ILNBLNK( mdsioLocalDir )
103
104 C Assign special directory
105 if ( mdsioLocalDir .NE. ' ' ) then
106 write(pFname(1:80),'(2a)') mdsioLocalDir(1:pIL), fName(1:IL)
107 else
108 pFname= fName
109 endif
110 pIL=ILNBLNK( pfName )
111
112 C Assign a free unit number as the I/O channel for this routine
113 call MDSFINDUNIT( dUnit, mythid )
114
115 C Check first for global file with simple name (ie. fName)
116 dataFName = fName
117 inquire( file=dataFname, exist=exst )
118 if (exst) then
119 if ( debugLevel .GE. debLevA ) then
120 write(msgbuf,'(a,a)')
121 & ' MDSREADFIELD: opening global file: ',dataFName
122 call print_message( msgbuf, standardmessageunit,
123 & SQUEEZE_RIGHT , mythid)
124 endif
125 globalFile = .TRUE.
126 endif
127
128 C If negative check for global file with MDS name (ie. fName.data)
129 if (.NOT. globalFile) then
130 write(dataFname(1:80),'(2a)') fName(1:IL),'.data'
131 inquire( file=dataFname, exist=exst )
132 if (exst) then
133 if ( debugLevel .GE. debLevA ) then
134 write(msgbuf,'(a,a)')
135 & ' MDSREADFIELD: opening global file: ',dataFName
136 call print_message( msgbuf, standardmessageunit,
137 & SQUEEZE_RIGHT , mythid)
138 endif
139 globalFile = .TRUE.
140 endif
141 endif
142
143 if ( .not. ( globalFile .and. useSingleCPUIO ) ) then
144
145 C If we are reading from a global file then we open it here
146 if (globalFile) then
147 length_of_rec=MDS_RECLEN( filePrec, sNx, mythid )
148 open( dUnit, file=dataFName, status='old',
149 & access='direct', recl=length_of_rec )
150 fileIsOpen=.TRUE.
151 endif
152
153 #ifdef ALLOW_EXCH2
154
155 if (globalFile) then
156 domainLength = exch2_domain_nxt
157 domainHeight = exch2_domain_nyt
158
159 cafe write(*,fmt='(1X,A,I3,A,I3)') 'L=', domainlength,
160 cafe & 'H=', domainheight
161 C Loop over all tiles
162 do bj=1,nSy
163 do bi=1,nSx
164 C If we are reading from a tiled MDS file then we open each one here
165 cafe if (.NOT. globalFile) then
166 cafe write(msgbuf,'(a)')
167 cafe & ' MDSREADFIELD: non-global input files not
168 cafe & implemented with exch2'
169 cafe call print_error( msgbuf, mythid )
170 cafe stop 'ABNORMAL END: S/R MDSREADFIELD'
171 cafe endif
172 tn = W2_myTileList(bi)
173 tby = exch2_tbasey(tn)
174 tgx = exch2_txglobalo(tn)
175 tny = exch2_tny(tn)
176 tnx = exch2_tnx(tn)
177 if (fileIsOpen) then
178 do k=1,nNz
179 do j=1,tNy
180 cafe write(*,fmt='(1X,A,I3,A,I3,A,I3,A,I3,A,I3,A,I3)') 'tby=', tby,
181 cafe & ', tgx=', tgx,
182 cafe & ', tnx=',tnx, ', tny=', tny, ', j=',j,', tn=',tn
183
184 irec = domainLength*tby + (tgx-1)/tnx + 1 +
185 & domainLength*(j-1) +
186 & domainLength*domainHeight*tny*(k-1) +
187 & domainLength*domainHeight*tny*nNz*(irecord-1)
188 cafe write(*,fmt='(1X,A,I6,A,I3)') 'record = ',irec,',thingy=',
189 cafe & (tgx-1)/tnx
190 cafe write(*,fmt='(1X,A,I6)') 'record = ',irec
191
192 if (filePrec .eq. precFloat32) then
193 read(dUnit,rec=irec) r4seg
194 #ifdef _BYTESWAPIO
195 call MDS_BYTESWAPR4( sNx, r4seg )
196 #endif
197 if (arrType .eq. 'RS') then
198 call MDS_SEG4toRS( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
199 elseif (arrType .eq. 'RL') then
200 call MDS_SEG4toRL( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
201 else
202 write(msgbuf,'(a)')
203 & ' MDSREADFIELD: illegal value for arrType'
204 call print_error( msgbuf, mythid )
205 stop 'ABNORMAL END: S/R MDSREADFIELD'
206 endif
207 elseif (filePrec .eq. precFloat64) then
208 read(dUnit,rec=irec) r8seg
209 #ifdef _BYTESWAPIO
210 call MDS_BYTESWAPR8( sNx, r8seg )
211 #endif
212 if (arrType .eq. 'RS') then
213 call MDS_SEG8toRS( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
214 elseif (arrType .eq. 'RL') then
215 call MDS_SEG8toRL( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
216 else
217 write(msgbuf,'(a)')
218 & ' MDSREADFIELD: illegal value for arrType'
219 call print_error( msgbuf, mythid )
220 stop 'ABNORMAL END: S/R MDSREADFIELD'
221 endif
222 else
223 write(msgbuf,'(a)')
224 & ' MDSREADFIELD: illegal value for filePrec'
225 call print_error( msgbuf, mythid )
226 stop 'ABNORMAL END: S/R MDSREADFIELD'
227 endif
228 C End of j loop
229 enddo
230 C End of k loop
231 enddo
232 if (.NOT. globalFile) then
233 close( dUnit )
234 fileIsOpen = .FALSE.
235 endif
236 endif
237 C End of bi,bj loops
238 enddo
239 enddo
240 else ! if globalFile
241 c#else /* .not. ALLOW_EXCH2 */
242 #endif /* ALLOW_EXCH2 */
243
244
245 C Loop over all tiles
246 do bj=1,nSy
247 do bi=1,nSx
248 C If we are reading from a tiled MDS file then we open each one here
249 if (.NOT. globalFile) then
250 iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
251 jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
252 write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)')
253 & pfName(1:pIL),'.',iG,'.',jG,'.data'
254 inquire( file=dataFname, exist=exst )
255 C Of course, we only open the file if the tile is "active"
256 C (This is a place-holder for the active/passive mechanism
257 if (exst) then
258 if ( debugLevel .GE. debLevA ) then
259 write(msgbuf,'(a,a)')
260 & ' MDSREADFIELD: opening file: ',dataFName
261 call print_message( msgbuf, standardmessageunit,
262 & SQUEEZE_RIGHT , mythid)
263 endif
264 length_of_rec=MDS_RECLEN( filePrec, sNx, mythid )
265 open( dUnit, file=dataFName, status='old',
266 & access='direct', recl=length_of_rec )
267 fileIsOpen=.TRUE.
268 else
269 fileIsOpen=.FALSE.
270 write(msgbuf,'(3a)')
271 & ' MDSREADFIELD: filename: ',dataFName, pfName
272 call print_message( msgbuf, standardmessageunit,
273 & SQUEEZE_RIGHT , mythid)
274 write(msgbuf,'(a)')
275 & ' MDSREADFIELD: File does not exist'
276 call print_error( msgbuf, mythid )
277 stop 'ABNORMAL END: S/R MDSREADFIELD'
278 endif
279 endif
280
281 if (fileIsOpen) then
282 do k=1,nNz
283 do j=1,sNy
284 if (globalFile) then
285 iG = myXGlobalLo-1 + (bi-1)*sNx
286 jG = myYGlobalLo-1 + (bj-1)*sNy
287 irec=1 + INT(iG/sNx) + nSx*nPx*(jG+j-1) + nSx*nPx*Ny*(k-1)
288 & + nSx*nPx*Ny*nNz*(irecord-1)
289 else
290 iG = 0
291 jG = 0
292 irec=j + sNy*(k-1) + sNy*nNz*(irecord-1)
293 endif
294 if (filePrec .eq. precFloat32) then
295 read(dUnit,rec=irec) r4seg
296 #ifdef _BYTESWAPIO
297 call MDS_BYTESWAPR4( sNx, r4seg )
298 #endif
299 if (arrType .eq. 'RS') then
300 call MDS_SEG4toRS( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
301 elseif (arrType .eq. 'RL') then
302 call MDS_SEG4toRL( j,bi,bj,k,nNz, r4seg, .TRUE., arr )
303 else
304 write(msgbuf,'(a)')
305 & ' MDSREADFIELD: illegal value for arrType'
306 call print_error( msgbuf, mythid )
307 stop 'ABNORMAL END: S/R MDSREADFIELD'
308 endif
309 elseif (filePrec .eq. precFloat64) then
310 read(dUnit,rec=irec) r8seg
311 #ifdef _BYTESWAPIO
312 call MDS_BYTESWAPR8( sNx, r8seg )
313 #endif
314 if (arrType .eq. 'RS') then
315 call MDS_SEG8toRS( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
316 elseif (arrType .eq. 'RL') then
317 call MDS_SEG8toRL( j,bi,bj,k,nNz, r8seg, .TRUE., arr )
318 else
319 write(msgbuf,'(a)')
320 & ' MDSREADFIELD: illegal value for arrType'
321 call print_error( msgbuf, mythid )
322 stop 'ABNORMAL END: S/R MDSREADFIELD'
323 endif
324 else
325 write(msgbuf,'(a)')
326 & ' MDSREADFIELD: illegal value for filePrec'
327 call print_error( msgbuf, mythid )
328 stop 'ABNORMAL END: S/R MDSREADFIELD'
329 endif
330 C End of j loop
331 enddo
332 C End of k loop
333 enddo
334 if (.NOT. globalFile) then
335 close( dUnit )
336 fileIsOpen = .FALSE.
337 endif
338 endif
339 C End of bi,bj loops
340 enddo
341 enddo
342
343 #ifdef ALLOW_EXCH2
344 endif ! globalFile
345
346 #endif /* ALLOW_EXCH2 */
347
348 C If global file was opened then close it
349 if (fileIsOpen .AND. globalFile) then
350 close( dUnit )
351 fileIsOpen = .FALSE.
352 endif
353
354 endif
355 c endif ( .not. ( globalFile .and. useSingleCPUIO ) )
356
357 _END_MASTER( myThid )
358
359 if ( globalFile .and. useSingleCPUIO ) then
360
361 #ifdef ALLOW_EXCH2
362 write(msgbuf,'(a)')
363 & ' MDSREADFIELD: useSingleCPUIO not
364 & implemented with exch2'
365 call print_error( msgbuf, mythid )
366 stop 'ABNORMAL END: S/R MDSREADFIELD'
367 #endif /* ALLOW_EXCH2 */
368
369 C master thread of process 0, only, opens a global file
370 _BEGIN_MASTER( myThid )
371 #ifdef ALLOW_USE_MPI
372 IF( mpiMyId .EQ. 0 ) THEN
373 #else
374 IF ( .TRUE. ) THEN
375 #endif /* ALLOW_USE_MPI */
376 length_of_rec=MDS_RECLEN( filePrec, Nx*Ny, mythid )
377 open( dUnit, file=dataFName, status='old',
378 & access='direct', recl=length_of_rec )
379 ENDIF
380 _END_MASTER( myThid )
381
382 DO k=1,nNz
383
384 _BEGIN_MASTER( myThid )
385 #ifdef ALLOW_USE_MPI
386 IF( mpiMyId .EQ. 0 ) THEN
387 #else
388 IF ( .TRUE. ) THEN
389 #endif /* ALLOW_USE_MPI */
390 irec = k+nNz*(irecord-1)
391 if (filePrec .eq. precFloat32) then
392 read(dUnit,rec=irec) global_r4
393 #ifdef _BYTESWAPIO
394 call MDS_BYTESWAPR4( Nx*Ny, global_r4 )
395 #endif
396 DO J=1,Ny
397 DO I=1,Nx
398 global(I,J) = global_r4(I,J)
399 ENDDO
400 ENDDO
401 elseif (filePrec .eq. precFloat64) then
402 read(dUnit,rec=irec) global
403 #ifdef _BYTESWAPIO
404 call MDS_BYTESWAPR8( Nx*Ny, global )
405 #endif
406 else
407 write(msgbuf,'(a)')
408 & ' MDSREADFIELD: illegal value for filePrec'
409 call print_error( msgbuf, mythid )
410 stop 'ABNORMAL END: S/R MDSREADFIELD'
411 endif
412 ENDIF
413 _END_MASTER( myThid )
414
415 CALL SCATTER_2D(global,local,mythid)
416 if (arrType .eq. 'RS') then
417 call PASStoRS( local,arr,k,nNz,mythid )
418 elseif (arrType .eq. 'RL') then
419 call PASStoRL( local,arr,k,nNz,mythid )
420 else
421 write(msgbuf,'(a)')
422 & ' MDSREADFIELD: illegal value for arrType'
423 call print_error( msgbuf, mythid )
424 stop 'ABNORMAL END: S/R MDSREADFIELD'
425 endif
426
427 ENDDO
428 c ENDDO k=1,nNz
429
430 _BEGIN_MASTER( myThid )
431 close( dUnit )
432 _END_MASTER( myThid )
433
434 endif
435 c endif ( globalFile .and. useSingleCPUIO )
436
437 return
438 end
439
440
441 C ==================================================================
442
443 subroutine passToRS(local,arr,k,nNz,mythid)
444 implicit none
445 #include "EEPARAMS.h"
446 #include "SIZE.h"
447 _RL local(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
448 integer i,j,k,bi,bj,nNz
449 _RS arr(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nNz,nSx,nSy)
450 integer mythid
451 DO bj = myByLo(myThid), myByHi(myThid)
452 DO bi = myBxLo(myThid), myBxHi(myThid)
453 DO J=1-Oly,sNy+Oly
454 DO I=1-Olx,sNx+Olx
455 arr(I,J,k,bi,bj) = local(I,J,bi,bj)
456 ENDDO
457 ENDDO
458 ENDDO
459 ENDDO
460 return
461 end
462
463 subroutine passToRL(local,arr,k,nNz,mythid)
464 implicit none
465 #include "EEPARAMS.h"
466 #include "SIZE.h"
467 _RL local(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
468 integer i,j,k,bi,bj,nNz
469 _RL arr(1-oLx:sNx+oLx,1-oLy:sNy+oLy,nNz,nSx,nSy)
470 integer mythid
471 DO bj = myByLo(myThid), myByHi(myThid)
472 DO bi = myBxLo(myThid), myBxHi(myThid)
473 DO J=1-Oly,sNy+Oly
474 DO I=1-Olx,sNx+Olx
475 arr(I,J,k,bi,bj) = local(I,J,bi,bj)
476 ENDDO
477 ENDDO
478 ENDDO
479 ENDDO
480 return
481 end

  ViewVC Help
Powered by ViewVC 1.1.22