pkg/exf/exf_interp_read.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.10 2007/04/09 23:57:50 dimitri Exp $
C $Name:  $

#include "EXF_OPTIONS.h"

CBOP
C !ROUTINE: EXF_INTERP_READ
C !INTERFACE:
       SUBROUTINE EXF_INTERP_READ(
     I                infile, filePrec,
     O                arrayin,
     I                irecord, nx_in, ny_in, myThid )

C !DESCRIPTION:

C !USES:
      IMPLICIT NONE

C Global variables / common blocks
#include "SIZE.h"
#include "EEPARAMS.h"
#include "EXF_PARAM.h"
#ifdef ALLOW_USE_MPI
# include "EESUPPORT.h"
#endif /* ALLOW_USE_MPI */
#include "PARAMS.h"


C !INPUT/OUTPUT PARAMETERS:
C  infile      (string)  :: name of the binary input file (direct access)
C  filePrec    (integer) :: number of bits per word in file (32 or 64)
C  arrayin     ( _RL )   :: array to read file into
C  irecord     (integer) :: record number to read
C  nx_in,ny_in (integer) :: size in x & y direction of input file to read
C  myThid      (integer) :: My Thread Id number

      CHARACTER*(*) infile
      INTEGER       filePrec, irecord, nx_in, ny_in
       _RL          arrayin( -1:nx_in+2 , -1:ny_in+2 )
      INTEGER       myThid
CEOP

C !FUNCTIONS
      INTEGER MDS_RECLEN
      LOGICAL MASTER_CPU_IO
      EXTERNAL MDS_RECLEN
      EXTERNAL MASTER_CPU_IO

C !LOCAL VARIABLES
      INTEGER  i, j
      INTEGER  ioUnit, length_of_rec
#ifdef EXF_INTERP_USE_DYNALLOC
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
C     When using threads the address of the local automatic array
C     "buffer" is not visible to the other threads. So we create
C     a pointer to share that address here. This is presently
C     in an ifdef because it won't go through g77 and I'm not
C     currently sure what TAF would do with this.
      COMMON /EXF_IOPTR8/ glPtr8
      REAL*8, POINTER :: glPtr8(:,:)
      COMMON /EXF_IOPTR4/ glPtr4
      REAL*4, POINTER :: glPtr4(:,:)

      Real*8, target ::  buffer_r8(nx_in,ny_in)
      Real*4, target ::  buffer_r4(nx_in,ny_in)
#else  /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
      Real*8   buffer_r8(nx_in,ny_in)
      Real*4   buffer_r4(nx_in,ny_in)
#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
#else  /* ndef EXF_INTERP_USE_DYNALLOC */
      Real*8   buffer_r8(exf_interp_bufferSize)
      Real*4   buffer_r4(exf_interp_bufferSize)
      COMMON /EXF_INTERP_BUFFER/ buffer_r8, buffer_r4
      INTEGER ijs
#endif /* ndef EXF_INTERP_USE_DYNALLOC */
#ifdef ALLOW_USE_MPI
      INTEGER  ierr
#endif

C--   Check for consistency:
#ifdef EXF_INTERP_USE_DYNALLOC
#ifndef EXF_IREAD_USE_GLOBAL_POINTER
C     The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the
C     case of nThreads > 1. Stop IF it isnt.
      IF ( nThreads .GT. 1 ) THEN
      STOP
     &'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'
      ENDIF
#endif
#else  /* ndef EXF_INTERP_USE_DYNALLOC */
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
      STOP
     &'EXF_INTERP_READ: USE_GLOBAL_POINTER needs INTERP_USE_DYNALLOC'
#endif
      IF ( nx_in*ny_in .GT. exf_interp_bufferSize ) THEN
        STOP 'EXF_INTERP_READ: exf_interp_bufferSize too small'
      ENDIF
#endif /* ndef EXF_INTERP_USE_DYNALLOC */

C--   before starting to read, wait for everyone to finish
      _BARRIER

C---  read in input data

      IF ( MASTER_CPU_IO(myThid) ) THEN
C--   master thread of process 0, only, opens a global file

        CALL MDSFINDUNIT( ioUnit, myThid )
        length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, myThid )
        OPEN( ioUnit, file=infile, status='old', access='direct',
     &       recl=length_of_rec )
        IF ( filePrec .EQ. 32 ) THEN
#ifdef EXF_INTERP_USE_DYNALLOC
          READ(ioUnit,rec=irecord)  buffer_r4
#else
          READ(ioUnit,rec=irecord) (buffer_r4(i),i=1,nx_in*ny_in)
#endif
#ifdef _BYTESWAPIO
          CALL MDS_BYTESWAPR4(nx_in*ny_in,buffer_r4)
#endif /* _BYTESWAPIO */
        ELSE
#ifdef EXF_INTERP_USE_DYNALLOC
          READ(ioUnit,rec=irecord)  buffer_r8
#else
          READ(ioUnit,rec=irecord) (buffer_r8(i),i=1,nx_in*ny_in)
#endif
#ifdef _BYTESWAPIO
          CALL MDS_BYTESWAPR8(nx_in*ny_in,buffer_r8)
#endif /* _BYTESWAPIO */
        ENDIF
        CLOSE( ioUnit )
C--   end if MASTER_CPU_IO
      ENDIF

      _BEGIN_MASTER( myThid )
#ifdef ALLOW_USE_MPI
C--   broadcast to all processes
       IF ( useSingleCPUIO ) THEN
         IF ( filePrec .EQ. 32 ) THEN
           CALL MPI_BCAST(buffer_r4,nx_in*ny_in,MPI_REAL,
     &          0,MPI_COMM_MODEL,ierr)
         ELSE
           CALL MPI_BCAST(buffer_r8,nx_in*ny_in,MPI_DOUBLE_PRECISION,
     &          0,MPI_COMM_MODEL,ierr)
         ENDIF
       ENDIF
#endif /* ALLOW_USE_MPI */

#ifdef EXF_IREAD_USE_GLOBAL_POINTER
       IF ( filePrec .EQ. 32 ) THEN
         glPtr4 => buffer_r4
       ELSE
         glPtr8 => buffer_r8
       ENDIF
#endif
      _END_MASTER( myThid )
      _BARRIER

C---  Transfer buffer to "arrayin" array:
#ifdef EXF_INTERP_USE_DYNALLOC
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
      IF ( filePrec .EQ. 32 ) THEN
        DO j=1,ny_in
          DO i=1,nx_in
            arrayin(i,j)=glPtr4(i,j)
          ENDDO
        ENDDO
      ELSE
        DO j=1,ny_in
          DO i=1,nx_in
            arrayin(i,j)=glPtr8(i,j)
          ENDDO
        ENDDO
      ENDIF
#else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
      IF ( filePrec .EQ. 32 ) THEN
        DO j=1,ny_in
          DO i=1,nx_in
            arrayin(i,j)=buffer_r4(i,j)
          ENDDO
        ENDDO
      ELSE
        DO j=1,ny_in
          DO i=1,nx_in
            arrayin(i,j)=buffer_r8(i,j)
          ENDDO
        ENDDO
      ENDIF
#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
#else  /* ndef EXF_INTERP_USE_DYNALLOC */
      IF ( filePrec .EQ. 32 ) THEN
        DO j=1,ny_in
          ijs = (j-1)*nx_in
          DO i=1,nx_in
            arrayin(i,j)=buffer_r4(i+ijs)
          ENDDO
        ENDDO
      ELSE
        DO j=1,ny_in
          ijs = (j-1)*nx_in
          DO i=1,nx_in
            arrayin(i,j)=buffer_r8(i+ijs)
          ENDDO
        ENDDO
      ENDIF
#endif /* ndef EXF_INTERP_USE_DYNALLOC */

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.10 2007/04/09 23:57:50 dimitri Exp $
2	C $Name: $
3
4	#include "EXF_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: EXF_INTERP_READ
8	C !INTERFACE:
9	SUBROUTINE EXF_INTERP_READ(
10	I infile, filePrec,
11	O arrayin,
12	I irecord, nx_in, ny_in, myThid )
13
14	C !DESCRIPTION:
15
16	C !USES:
17	IMPLICIT NONE
18
19	C Global variables / common blocks
20	#include "SIZE.h"
21	#include "EEPARAMS.h"
22	#include "EXF_PARAM.h"
23	#ifdef ALLOW_USE_MPI
24	# include "EESUPPORT.h"
25	#endif /* ALLOW_USE_MPI */
26	#include "PARAMS.h"
27
28
29	C !INPUT/OUTPUT PARAMETERS:
30	C infile (string) :: name of the binary input file (direct access)
31	C filePrec (integer) :: number of bits per word in file (32 or 64)
32	C arrayin ( _RL ) :: array to read file into
33	C irecord (integer) :: record number to read
34	C nx_in,ny_in (integer) :: size in x & y direction of input file to read
35	C myThid (integer) :: My Thread Id number
36
37	CHARACTER() infile
38	INTEGER filePrec, irecord, nx_in, ny_in
39	_RL arrayin( -1:nx_in+2 , -1:ny_in+2 )
40	INTEGER myThid
41	CEOP
42
43	C !FUNCTIONS
44	INTEGER MDS_RECLEN
45	LOGICAL MASTER_CPU_IO
46	EXTERNAL MDS_RECLEN
47	EXTERNAL MASTER_CPU_IO
48
49	C !LOCAL VARIABLES
50	INTEGER i, j
51	INTEGER ioUnit, length_of_rec
52	#ifdef EXF_INTERP_USE_DYNALLOC
53	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
54	C When using threads the address of the local automatic array
55	C "buffer" is not visible to the other threads. So we create
56	C a pointer to share that address here. This is presently
57	C in an ifdef because it won't go through g77 and I'm not
58	C currently sure what TAF would do with this.
59	COMMON /EXF_IOPTR8/ glPtr8
60	REAL*8, POINTER :: glPtr8(:,:)
61	COMMON /EXF_IOPTR4/ glPtr4
62	REAL*4, POINTER :: glPtr4(:,:)
63
64	Real*8, target :: buffer_r8(nx_in,ny_in)
65	Real*4, target :: buffer_r4(nx_in,ny_in)
66	#else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
67	Real*8 buffer_r8(nx_in,ny_in)
68	Real*4 buffer_r4(nx_in,ny_in)
69	#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
70	#else /* ndef EXF_INTERP_USE_DYNALLOC */
71	Real*8 buffer_r8(exf_interp_bufferSize)
72	Real*4 buffer_r4(exf_interp_bufferSize)
73	COMMON /EXF_INTERP_BUFFER/ buffer_r8, buffer_r4
74	INTEGER ijs
75	#endif /* ndef EXF_INTERP_USE_DYNALLOC */
76	#ifdef ALLOW_USE_MPI
77	INTEGER ierr
78	#endif
79
80	C-- Check for consistency:
81	#ifdef EXF_INTERP_USE_DYNALLOC
82	#ifndef EXF_IREAD_USE_GLOBAL_POINTER
83	C The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the
84	C case of nThreads > 1. Stop IF it isnt.
85	IF ( nThreads .GT. 1 ) THEN
86	STOP
87	&'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'
88	ENDIF
89	#endif
90	#else /* ndef EXF_INTERP_USE_DYNALLOC */
91	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
92	STOP
93	&'EXF_INTERP_READ: USE_GLOBAL_POINTER needs INTERP_USE_DYNALLOC'
94	#endif
95	IF ( nx_in*ny_in .GT. exf_interp_bufferSize ) THEN
96	STOP 'EXF_INTERP_READ: exf_interp_bufferSize too small'
97	ENDIF
98	#endif /* ndef EXF_INTERP_USE_DYNALLOC */
99
100	C-- before starting to read, wait for everyone to finish
101	_BARRIER
102
103	C--- read in input data
104
105	IF ( MASTER_CPU_IO(myThid) ) THEN
106	C-- master thread of process 0, only, opens a global file
107
108	CALL MDSFINDUNIT( ioUnit, myThid )
109	length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, myThid )
110	OPEN( ioUnit, file=infile, status='old', access='direct',
111	& recl=length_of_rec )
112	IF ( filePrec .EQ. 32 ) THEN
113	#ifdef EXF_INTERP_USE_DYNALLOC
114	READ(ioUnit,rec=irecord) buffer_r4
115	#else
116	READ(ioUnit,rec=irecord) (buffer_r4(i),i=1,nx_in*ny_in)
117	#endif
118	#ifdef _BYTESWAPIO
119	CALL MDS_BYTESWAPR4(nx_in*ny_in,buffer_r4)
120	#endif /* _BYTESWAPIO */
121	ELSE
122	#ifdef EXF_INTERP_USE_DYNALLOC
123	READ(ioUnit,rec=irecord) buffer_r8
124	#else
125	READ(ioUnit,rec=irecord) (buffer_r8(i),i=1,nx_in*ny_in)
126	#endif
127	#ifdef _BYTESWAPIO
128	CALL MDS_BYTESWAPR8(nx_in*ny_in,buffer_r8)
129	#endif /* _BYTESWAPIO */
130	ENDIF
131	CLOSE( ioUnit )
132	C-- end if MASTER_CPU_IO
133	ENDIF
134
135	_BEGIN_MASTER( myThid )
136	#ifdef ALLOW_USE_MPI
137	C-- broadcast to all processes
138	IF ( useSingleCPUIO ) THEN
139	IF ( filePrec .EQ. 32 ) THEN
140	CALL MPI_BCAST(buffer_r4,nx_in*ny_in,MPI_REAL,
141	& 0,MPI_COMM_MODEL,ierr)
142	ELSE
143	CALL MPI_BCAST(buffer_r8,nx_in*ny_in,MPI_DOUBLE_PRECISION,
144	& 0,MPI_COMM_MODEL,ierr)
145	ENDIF
146	ENDIF
147	#endif /* ALLOW_USE_MPI */
148
149	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
150	IF ( filePrec .EQ. 32 ) THEN
151	glPtr4 => buffer_r4
152	ELSE
153	glPtr8 => buffer_r8
154	ENDIF
155	#endif
156	_END_MASTER( myThid )
157	_BARRIER
158
159	C--- Transfer buffer to "arrayin" array:
160	#ifdef EXF_INTERP_USE_DYNALLOC
161	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
162	IF ( filePrec .EQ. 32 ) THEN
163	DO j=1,ny_in
164	DO i=1,nx_in
165	arrayin(i,j)=glPtr4(i,j)
166	ENDDO
167	ENDDO
168	ELSE
169	DO j=1,ny_in
170	DO i=1,nx_in
171	arrayin(i,j)=glPtr8(i,j)
172	ENDDO
173	ENDDO
174	ENDIF
175	#else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
176	IF ( filePrec .EQ. 32 ) THEN
177	DO j=1,ny_in
178	DO i=1,nx_in
179	arrayin(i,j)=buffer_r4(i,j)
180	ENDDO
181	ENDDO
182	ELSE
183	DO j=1,ny_in
184	DO i=1,nx_in
185	arrayin(i,j)=buffer_r8(i,j)
186	ENDDO
187	ENDDO
188	ENDIF
189	#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
190	#else /* ndef EXF_INTERP_USE_DYNALLOC */
191	IF ( filePrec .EQ. 32 ) THEN
192	DO j=1,ny_in
193	ijs = (j-1)*nx_in
194	DO i=1,nx_in
195	arrayin(i,j)=buffer_r4(i+ijs)
196	ENDDO
197	ENDDO
198	ELSE
199	DO j=1,ny_in
200	ijs = (j-1)*nx_in
201	DO i=1,nx_in
202	arrayin(i,j)=buffer_r8(i+ijs)
203	ENDDO
204	ENDDO
205	ENDIF
206	#endif /* ndef EXF_INTERP_USE_DYNALLOC */
207
208	RETURN
209	END