pkg/exf/exf_interp_read.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.15 2013/12/19 01:06:10 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 "PARAMS.h"
#include "EXF_INTERP_SIZE.h"
#include "EXF_PARAM.h"
#ifdef ALLOW_USE_MPI
# include "EESUPPORT.h"
#endif /* ALLOW_USE_MPI */

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  ILNBLNK
      INTEGER MDS_RECLEN
      LOGICAL MASTER_CPU_IO
      EXTERNAL ILNBLNK
      EXTERNAL MDS_RECLEN
      EXTERNAL MASTER_CPU_IO

C !LOCAL VARIABLES
      INTEGER  i, j
      INTEGER  ioUnit, length_of_rec, IL
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      LOGICAL  exst
#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

        IL  = ILNBLNK( infile )
        INQUIRE( file=infile, exist=exst )
        IF (exst) THEN
         IF ( debugLevel.GE.debLevB ) THEN
          WRITE(msgBuf,'(A,A)')
     &         ' EXF_INTERP_READ: opening file: ',infile(1:IL)
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &         SQUEEZE_RIGHT , myThid)
         ENDIF
        ELSE
         WRITE(msgBuf,'(2A)')
     &        ' EXF_INTERP_READ: filename: ', infile(1:IL)
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT , myThid)
         CALL PRINT_ERROR( msgBuf, myThid )
         WRITE(msgBuf,'(A)')
     &        ' EXF_INTERP_READ: File does not exist'
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT , myThid)
         CALL PRINT_ERROR( msgBuf, myThid )
         STOP 'ABNORMAL END: S/R EXF_INTERP_READ'
        ENDIF

        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 ( useSingleCpuInput ) 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	jmc	1.16	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.15 2013/12/19 01:06:10 dimitri Exp $
2	jmc	1.5	C $Name: $
3
4	heimbach	1.1	#include "EXF_OPTIONS.h"
5
6	jmc	1.11	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	heimbach	1.1
14	jmc	1.11	C !DESCRIPTION:
15
16			C !USES:
17			IMPLICIT NONE
18
19			C Global variables / common blocks
20	heimbach	1.1	#include "SIZE.h"
21			#include "EEPARAMS.h"
22	jmc	1.16	#include "PARAMS.h"
23			#include "EXF_INTERP_SIZE.h"
24	jmc	1.11	#include "EXF_PARAM.h"
25	heimbach	1.1	#ifdef ALLOW_USE_MPI
26			# include "EESUPPORT.h"
27			#endif /* ALLOW_USE_MPI */
28	jmc	1.11
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			CHARACTER() infile
37			INTEGER filePrec, irecord, nx_in, ny_in
38			_RL arrayin( -1:nx_in+2 , -1:ny_in+2 )
39			INTEGER myThid
40			CEOP
41
42			C !FUNCTIONS
43	mlosch	1.12	INTEGER ILNBLNK
44	jmc	1.11	INTEGER MDS_RECLEN
45			LOGICAL MASTER_CPU_IO
46	mlosch	1.12	EXTERNAL ILNBLNK
47	jmc	1.11	EXTERNAL MDS_RECLEN
48			EXTERNAL MASTER_CPU_IO
49
50			C !LOCAL VARIABLES
51			INTEGER i, j
52	mlosch	1.12	INTEGER ioUnit, length_of_rec, IL
53			CHARACTER*(MAX_LEN_MBUF) msgBuf
54			LOGICAL exst
55	jmc	1.11	#ifdef EXF_INTERP_USE_DYNALLOC
56	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
57			C When using threads the address of the local automatic array
58	jmc	1.11	C "buffer" is not visible to the other threads. So we create
59	cnh	1.3	C a pointer to share that address here. This is presently
60			C in an ifdef because it won't go through g77 and I'm not
61			C currently sure what TAF would do with this.
62	dimitri	1.9	COMMON /EXF_IOPTR8/ glPtr8
63			REAL*8, POINTER :: glPtr8(:,:)
64	jmc	1.11	COMMON /EXF_IOPTR4/ glPtr4
65			REAL*4, POINTER :: glPtr4(:,:)
66	cnh	1.3
67	jmc	1.11	Real*8, target :: buffer_r8(nx_in,ny_in)
68			Real*4, target :: buffer_r4(nx_in,ny_in)
69			#else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
70			Real*8 buffer_r8(nx_in,ny_in)
71			Real*4 buffer_r4(nx_in,ny_in)
72			#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
73			#else /* ndef EXF_INTERP_USE_DYNALLOC */
74			Real*8 buffer_r8(exf_interp_bufferSize)
75			Real*4 buffer_r4(exf_interp_bufferSize)
76			COMMON /EXF_INTERP_BUFFER/ buffer_r8, buffer_r4
77			INTEGER ijs
78			#endif /* ndef EXF_INTERP_USE_DYNALLOC */
79			#ifdef ALLOW_USE_MPI
80			INTEGER ierr
81	cnh	1.3	#endif
82	heimbach	1.1
83	jmc	1.11	C-- Check for consistency:
84			#ifdef EXF_INTERP_USE_DYNALLOC
85	cnh	1.3	#ifndef EXF_IREAD_USE_GLOBAL_POINTER
86	jmc	1.11	C The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the
87			C case of nThreads > 1. Stop IF it isnt.
88	cnh	1.3	IF ( nThreads .GT. 1 ) THEN
89	jmc	1.11	STOP
90	cnh	1.3	&'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'
91			ENDIF
92			#endif
93	jmc	1.11	#else /* ndef EXF_INTERP_USE_DYNALLOC */
94	dimitri	1.6	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
95	jmc	1.11	STOP
96			&'EXF_INTERP_READ: USE_GLOBAL_POINTER needs INTERP_USE_DYNALLOC'
97	dimitri	1.6	#endif
98	jmc	1.11	IF ( nx_in*ny_in .GT. exf_interp_bufferSize ) THEN
99			STOP 'EXF_INTERP_READ: exf_interp_bufferSize too small'
100			ENDIF
101			#endif /* ndef EXF_INTERP_USE_DYNALLOC */
102
103			C-- before starting to read, wait for everyone to finish
104			_BARRIER
105	heimbach	1.1
106	jmc	1.11	C--- read in input data
107	heimbach	1.1
108	jmc	1.11	IF ( MASTER_CPU_IO(myThid) ) THEN
109			C-- master thread of process 0, only, opens a global file
110	heimbach	1.1
111	mlosch	1.12	IL = ILNBLNK( infile )
112			INQUIRE( file=infile, exist=exst )
113			IF (exst) THEN
114	jmc	1.14	IF ( debugLevel.GE.debLevB ) THEN
115	mlosch	1.12	WRITE(msgBuf,'(A,A)')
116			& ' EXF_INTERP_READ: opening file: ',infile(1:IL)
117			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
118			& SQUEEZE_RIGHT , myThid)
119			ENDIF
120			ELSE
121			WRITE(msgBuf,'(2A)')
122			& ' EXF_INTERP_READ: filename: ', infile(1:IL)
123			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
124			& SQUEEZE_RIGHT , myThid)
125			CALL PRINT_ERROR( msgBuf, myThid )
126			WRITE(msgBuf,'(A)')
127			& ' EXF_INTERP_READ: File does not exist'
128			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
129			& SQUEEZE_RIGHT , myThid)
130			CALL PRINT_ERROR( msgBuf, myThid )
131			STOP 'ABNORMAL END: S/R EXF_INTERP_READ'
132			ENDIF
133
134	jmc	1.11	CALL MDSFINDUNIT( ioUnit, myThid )
135			length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, myThid )
136			OPEN( ioUnit, file=infile, status='old', access='direct',
137			& recl=length_of_rec )
138	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
139	jmc	1.11	#ifdef EXF_INTERP_USE_DYNALLOC
140			READ(ioUnit,rec=irecord) buffer_r4
141			#else
142			READ(ioUnit,rec=irecord) (buffer_r4(i),i=1,nx_in*ny_in)
143			#endif
144	dimitri	1.9	#ifdef _BYTESWAPIO
145	jmc	1.11	CALL MDS_BYTESWAPR4(nx_in*ny_in,buffer_r4)
146	dimitri	1.9	#endif /* _BYTESWAPIO */
147			ELSE
148	jmc	1.11	#ifdef EXF_INTERP_USE_DYNALLOC
149			READ(ioUnit,rec=irecord) buffer_r8
150			#else
151			READ(ioUnit,rec=irecord) (buffer_r8(i),i=1,nx_in*ny_in)
152			#endif
153	dimitri	1.9	#ifdef _BYTESWAPIO
154	jmc	1.11	CALL MDS_BYTESWAPR8(nx_in*ny_in,buffer_r8)
155	dimitri	1.9	#endif /* _BYTESWAPIO */
156			ENDIF
157	jmc	1.11	CLOSE( ioUnit )
158			C-- end if MASTER_CPU_IO
159			ENDIF
160	heimbach	1.1
161	jmc	1.11	_BEGIN_MASTER( myThid )
162	heimbach	1.1	#ifdef ALLOW_USE_MPI
163	jmc	1.11	C-- broadcast to all processes
164	dimitri	1.15	IF ( useSingleCpuInput ) THEN
165	jmc	1.11	IF ( filePrec .EQ. 32 ) THEN
166			CALL MPI_BCAST(buffer_r4,nx_in*ny_in,MPI_REAL,
167			& 0,MPI_COMM_MODEL,ierr)
168			ELSE
169			CALL MPI_BCAST(buffer_r8,nx_in*ny_in,MPI_DOUBLE_PRECISION,
170			& 0,MPI_COMM_MODEL,ierr)
171			ENDIF
172			ENDIF
173	heimbach	1.1	#endif /* ALLOW_USE_MPI */
174	jmc	1.11
175	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
176	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
177	jmc	1.11	glPtr4 => buffer_r4
178	dimitri	1.9	ELSE
179	jmc	1.11	glPtr8 => buffer_r8
180	dimitri	1.9	ENDIF
181	cnh	1.3	#endif
182			_END_MASTER( myThid )
183			_BARRIER
184	dimitri	1.9
185	jmc	1.11	C--- Transfer buffer to "arrayin" array:
186			#ifdef EXF_INTERP_USE_DYNALLOC
187	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
188	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
189	jmc	1.11	DO j=1,ny_in
190			DO i=1,nx_in
191			arrayin(i,j)=glPtr4(i,j)
192			ENDDO
193			ENDDO
194			ELSE
195			DO j=1,ny_in
196			DO i=1,nx_in
197			arrayin(i,j)=glPtr8(i,j)
198			ENDDO
199			ENDDO
200			ENDIF
201			#else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
202			IF ( filePrec .EQ. 32 ) THEN
203			DO j=1,ny_in
204			DO i=1,nx_in
205			arrayin(i,j)=buffer_r4(i,j)
206			ENDDO
207			ENDDO
208	dimitri	1.9	ELSE
209	jmc	1.11	DO j=1,ny_in
210			DO i=1,nx_in
211			arrayin(i,j)=buffer_r8(i,j)
212			ENDDO
213			ENDDO
214	dimitri	1.9	ENDIF
215	jmc	1.11	#endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
216			#else /* ndef EXF_INTERP_USE_DYNALLOC */
217	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
218	jmc	1.11	DO j=1,ny_in
219			ijs = (j-1)*nx_in
220			DO i=1,nx_in
221			arrayin(i,j)=buffer_r4(i+ijs)
222			ENDDO
223			ENDDO
224	dimitri	1.9	ELSE
225	jmc	1.11	DO j=1,ny_in
226			ijs = (j-1)*nx_in
227			DO i=1,nx_in
228			arrayin(i,j)=buffer_r8(i+ijs)
229			ENDDO
230			ENDDO
231	dimitri	1.9	ENDIF
232	jmc	1.11	#endif /* ndef EXF_INTERP_USE_DYNALLOC */
233	heimbach	1.1
234	jmc	1.5	RETURN
235	heimbach	1.1	END