pkg/exf/exf_interp_read.F

C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.8 2006/12/14 08:17:15 mlosch Exp $
C $Name:  $

#include "EXF_OPTIONS.h"

       SUBROUTINE exf_interp_read(
     I   infile, filePrec,
     O   arrayin,
     I   irecord, nx_in, ny_in, mythid)

      implicit none

C     infile       = name of the input file (direct access binary)
C     filePrec     = file precicision (currently not used, assumes real*4)
C     arrout       = output arrays (different for each processor)
C     irecord      = record number in global file
C     nx_in, ny_in = input x-grid and y-grid size
C     mythid       = thread id

#include "SIZE.h"
#include "EEPARAMS.h"
#ifdef ALLOW_USE_MPI
# include "EESUPPORT.h"
#endif /* ALLOW_USE_MPI */
#include "PARAMS.h"

#ifdef EXF_IREAD_USE_GLOBAL_POINTER
C     When using threads the address of the local automatic array
C     "global" 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_IOPTR/ glPtr
      REAL*4, POINTER :: glPtr(:,:)
      COMMON /EXF_IOPTR8/ glPtr8
      REAL*8, POINTER :: glPtr8(:,:)
#endif

C subroutine variables
      character*(*) infile
      integer       filePrec, irecord, nx_in, ny_in
      real*4        arrayin( -1:nx_in+2 , -1:ny_in+2 )
      integer       mythid

C Functions
      integer MDS_RECLEN

C local variables
      integer  ierr, length_of_rec
      real*8   ne_fac,nw_fac,se_fac,sw_fac
      integer  e_ind(snx,sny),w_ind(snx,sny)
      integer  n_ind(snx,sny),s_ind(snx,sny)
      real*8   px_ind(4), py_ind(4), ew_val(4)
      external lagran
      real*8   lagran
      integer  i, j, k, l, js, bi, bj, sp, interp_unit
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
      real*4, target ::   global(nx_in,ny_in)
      real*8, target ::   global8(nx_in,ny_in)
#else
      real*4   global(nx_in,ny_in)
      real*8   global8(nx_in,ny_in)
#endif

      _BARRIER
      _BEGIN_MASTER( myThid )

#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

C read in input data
#ifdef ALLOW_USE_MPI
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
       if (.FALSE.) then
#else
      if (useSingleCPUIO) then
#endif

C master thread of process 0, only, opens a global file
        IF( mpiMyId .EQ. 0 ) THEN
         call mdsfindunit( interp_unit, mythid)
         length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
         open(interp_unit,file=infile,status='old',access='direct',
     &        recl=length_of_rec)
         IF ( filePrec .EQ. 32 ) THEN
            read(interp_unit,rec=irecord)
     &           ((global(i,j),i=1,nx_in),j=1,ny_in)
#ifdef _BYTESWAPIO
            call MDS_BYTESWAPR4(nx_in*ny_in,global)
#endif /* _BYTESWAPIO */
         ELSE
            read(interp_unit,rec=irecord)
     &           ((global8(i,j),i=1,nx_in),j=1,ny_in)
#ifdef _BYTESWAPIO
            call MDS_BYTESWAPR8(nx_in*ny_in,global8)
#endif /* _BYTESWAPIO */
         ENDIF
         close(interp_unit)
        ENDIF

C broadcast to all processes
        IF ( filePrec .EQ. 32 ) THEN
           call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,
     &          0,MPI_COMM_MODEL,ierr)
        ELSE
           call MPI_BCAST(global8,nx_in*ny_in,MPI_DOUBLE_PRECISION,
     &          0,MPI_COMM_MODEL,ierr)
        ENDIF
       else
#endif /* ALLOW_USE_MPI */

        call mdsfindunit( interp_unit, mythid)
        length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
        open(interp_unit,file=infile,status='old',access='direct',
     &       recl=length_of_rec)
        IF ( filePrec .EQ. 32 ) THEN
           read(interp_unit,rec=irecord) global
#ifdef _BYTESWAPIO
           call MDS_BYTESWAPR4(nx_in*ny_in,global)
#endif /* _BYTESWAPIO */
        ELSE
           read(interp_unit,rec=irecord) global8
#ifdef _BYTESWAPIO
           call MDS_BYTESWAPR8(nx_in*ny_in,global8)
#endif /* _BYTESWAPIO */
        ENDIF
        close(interp_unit)

#ifdef ALLOW_USE_MPI
       endif
#endif /* ALLOW_USE_MPI */
#ifdef EXF_IREAD_USE_GLOBAL_POINTER
       IF ( filePrec .EQ. 32 ) THEN
          glPtr => global
       ELSE
          glPtr8 => global8
       ENDIF
#endif
      _END_MASTER( myThid )
      _BARRIER

#ifdef EXF_IREAD_USE_GLOBAL_POINTER
      IF ( filePrec .EQ. 32 ) THEN
         do j=1,ny_in
            do i=1,nx_in
               arrayin(i,j)=glPtr(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
      IF ( filePrec .EQ. 32 ) THEN
         do j=1,ny_in
            do i=1,nx_in
               arrayin(i,j)=global(i,j)
            enddo
         enddo
      ELSE
         do j=1,ny_in
            do i=1,nx_in
               arrayin(i,j)=global8(i,j)
            enddo
         enddo
      ENDIF
#endif

      RETURN
      END
1	dimitri	1.9	C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_interp_read.F,v 1.8 2006/12/14 08:17:15 mlosch Exp $
2	jmc	1.5	C $Name: $
3
4	heimbach	1.1	#include "EXF_OPTIONS.h"
5
6			SUBROUTINE exf_interp_read(
7	dimitri	1.4	I infile, filePrec,
8	heimbach	1.1	O arrayin,
9	dimitri	1.4	I irecord, nx_in, ny_in, mythid)
10	heimbach	1.1
11			implicit none
12
13			C infile = name of the input file (direct access binary)
14			C filePrec = file precicision (currently not used, assumes real*4)
15			C arrout = output arrays (different for each processor)
16			C irecord = record number in global file
17			C nx_in, ny_in = input x-grid and y-grid size
18			C mythid = thread id
19
20			#include "SIZE.h"
21			#include "EEPARAMS.h"
22			#ifdef ALLOW_USE_MPI
23			# include "EESUPPORT.h"
24			#endif /* ALLOW_USE_MPI */
25			#include "PARAMS.h"
26
27	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
28			C When using threads the address of the local automatic array
29			C "global" is not visible to the other threads. So we create
30			C a pointer to share that address here. This is presently
31			C in an ifdef because it won't go through g77 and I'm not
32			C currently sure what TAF would do with this.
33			COMMON /EXF_IOPTR/ glPtr
34			REAL*4, POINTER :: glPtr(:,:)
35	dimitri	1.9	COMMON /EXF_IOPTR8/ glPtr8
36			REAL*8, POINTER :: glPtr8(:,:)
37	cnh	1.3	#endif
38
39	heimbach	1.1	C subroutine variables
40			character() infile
41			integer filePrec, irecord, nx_in, ny_in
42	dimitri	1.9	real*4 arrayin( -1:nx_in+2 , -1:ny_in+2 )
43	dimitri	1.4	integer mythid
44	heimbach	1.1
45			C Functions
46			integer MDS_RECLEN
47
48			C local variables
49			integer ierr, length_of_rec
50			real*8 ne_fac,nw_fac,se_fac,sw_fac
51			integer e_ind(snx,sny),w_ind(snx,sny)
52			integer n_ind(snx,sny),s_ind(snx,sny)
53			real*8 px_ind(4), py_ind(4), ew_val(4)
54			external lagran
55			real*8 lagran
56			integer i, j, k, l, js, bi, bj, sp, interp_unit
57	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
58			real*4, target :: global(nx_in,ny_in)
59	dimitri	1.9	real*8, target :: global8(nx_in,ny_in)
60	cnh	1.3	#else
61	heimbach	1.1	real*4 global(nx_in,ny_in)
62	dimitri	1.9	real*8 global8(nx_in,ny_in)
63	cnh	1.3	#endif
64	heimbach	1.1
65	jmc	1.5	_BARRIER
66	cnh	1.3	_BEGIN_MASTER( myThid )
67	heimbach	1.1
68	cnh	1.3	#ifndef EXF_IREAD_USE_GLOBAL_POINTER
69			C The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the
70			C case of nThreads > 1. Stop if it isnt.
71			IF ( nThreads .GT. 1 ) THEN
72			STOP
73			&'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'
74			ENDIF
75			#endif
76	heimbach	1.1
77			C read in input data
78			#ifdef ALLOW_USE_MPI
79	dimitri	1.6	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
80	cnh	1.3	if (.FALSE.) then
81	dimitri	1.6	#else
82			if (useSingleCPUIO) then
83			#endif
84	heimbach	1.1
85			C master thread of process 0, only, opens a global file
86			IF( mpiMyId .EQ. 0 ) THEN
87			call mdsfindunit( interp_unit, mythid)
88			length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
89			open(interp_unit,file=infile,status='old',access='direct',
90			& recl=length_of_rec)
91	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
92			read(interp_unit,rec=irecord)
93			& ((global(i,j),i=1,nx_in),j=1,ny_in)
94			#ifdef _BYTESWAPIO
95			call MDS_BYTESWAPR4(nx_in*ny_in,global)
96			#endif /* _BYTESWAPIO */
97			ELSE
98			read(interp_unit,rec=irecord)
99			& ((global8(i,j),i=1,nx_in),j=1,ny_in)
100			#ifdef _BYTESWAPIO
101			call MDS_BYTESWAPR8(nx_in*ny_in,global8)
102			#endif /* _BYTESWAPIO */
103			ENDIF
104	heimbach	1.1	close(interp_unit)
105			ENDIF
106
107			C broadcast to all processes
108	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
109			call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,
110			& 0,MPI_COMM_MODEL,ierr)
111			ELSE
112			call MPI_BCAST(global8,nx_in*ny_in,MPI_DOUBLE_PRECISION,
113			& 0,MPI_COMM_MODEL,ierr)
114			ENDIF
115	heimbach	1.1	else
116			#endif /* ALLOW_USE_MPI */
117
118			call mdsfindunit( interp_unit, mythid)
119			length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )
120			open(interp_unit,file=infile,status='old',access='direct',
121			& recl=length_of_rec)
122	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
123			read(interp_unit,rec=irecord) global
124			#ifdef _BYTESWAPIO
125			call MDS_BYTESWAPR4(nx_in*ny_in,global)
126			#endif /* _BYTESWAPIO */
127			ELSE
128			read(interp_unit,rec=irecord) global8
129			#ifdef _BYTESWAPIO
130			call MDS_BYTESWAPR8(nx_in*ny_in,global8)
131			#endif /* _BYTESWAPIO */
132			ENDIF
133	heimbach	1.1	close(interp_unit)
134
135			#ifdef ALLOW_USE_MPI
136			endif
137			#endif /* ALLOW_USE_MPI */
138	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
139	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
140			glPtr => global
141			ELSE
142			glPtr8 => global8
143			ENDIF
144	cnh	1.3	#endif
145			_END_MASTER( myThid )
146			_BARRIER
147	dimitri	1.9
148	cnh	1.3	#ifdef EXF_IREAD_USE_GLOBAL_POINTER
149	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
150			do j=1,ny_in
151			do i=1,nx_in
152			arrayin(i,j)=glPtr(i,j)
153			enddo
154			enddo
155			ELSE
156			do j=1,ny_in
157			do i=1,nx_in
158			arrayin(i,j)=glPtr8(i,j)
159			enddo
160			enddo
161			ENDIF
162	cnh	1.3	#else
163	dimitri	1.9	IF ( filePrec .EQ. 32 ) THEN
164			do j=1,ny_in
165			do i=1,nx_in
166			arrayin(i,j)=global(i,j)
167			enddo
168			enddo
169			ELSE
170			do j=1,ny_in
171			do i=1,nx_in
172			arrayin(i,j)=global8(i,j)
173			enddo
174			enddo
175			ENDIF
176	cnh	1.3	#endif
177	heimbach	1.1
178	jmc	1.5	RETURN
179	heimbach	1.1	END