/[MITgcm]/MITgcm/pkg/exf/exf_interp_read.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_interp_read.F

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

revision 1.10 by dimitri, Mon Apr 9 23:57:50 2007 UTC revision 1.11 by jmc, Thu May 10 22:21:55 2007 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6         SUBROUTINE exf_interp_read(  CBOP
7       I   infile, filePrec,  C !ROUTINE: EXF_INTERP_READ
8       O   arrayin,  C !INTERFACE:
9       I   irecord, nx_in, ny_in, mythid)         SUBROUTINE EXF_INTERP_READ(
10         I                infile, filePrec,
11        implicit none       O                arrayin,
12         I                irecord, nx_in, ny_in, myThid )
 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  
13    
14    C !DESCRIPTION:
15    
16    C !USES:
17          IMPLICIT NONE
18    
19    C Global variables / common blocks
20  #include "SIZE.h"  #include "SIZE.h"
21  #include "EEPARAMS.h"  #include "EEPARAMS.h"
22    #include "EXF_PARAM.h"
23  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
24  # include "EESUPPORT.h"  # include "EESUPPORT.h"
25  #endif /* ALLOW_USE_MPI */  #endif /* ALLOW_USE_MPI */
26  #include "PARAMS.h"  #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  #ifdef EXF_IREAD_USE_GLOBAL_POINTER
54  C     When using threads the address of the local automatic array  C     When using threads the address of the local automatic array
55  C     "global" is not visible to the other threads. So we create  C     "buffer" is not visible to the other threads. So we create
56  C     a pointer to share that address here. This is presently  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  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.  C     currently sure what TAF would do with this.
       COMMON /EXF_IOPTR/ glPtr  
       REAL*4, POINTER :: glPtr(:,:)  
59        COMMON /EXF_IOPTR8/ glPtr8        COMMON /EXF_IOPTR8/ glPtr8
60        REAL*8, POINTER :: glPtr8(:,:)        REAL*8, POINTER :: glPtr8(:,:)
61  #endif        COMMON /EXF_IOPTR4/ glPtr4
62          REAL*4, POINTER :: glPtr4(:,:)
63    
64  C subroutine variables        Real*8, target ::  buffer_r8(nx_in,ny_in)
65        character*(*) infile        Real*4, target ::  buffer_r4(nx_in,ny_in)
66        integer       filePrec, irecord, nx_in, ny_in  #else  /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
67        real*4        arrayin( -1:nx_in+2 , -1:ny_in+2 )        Real*8   buffer_r8(nx_in,ny_in)
68        integer       mythid        Real*4   buffer_r4(nx_in,ny_in)
69    #endif /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
70  C Functions  #else  /* ndef EXF_INTERP_USE_DYNALLOC */
71        integer MDS_RECLEN        Real*8   buffer_r8(exf_interp_bufferSize)
72          Real*4   buffer_r4(exf_interp_bufferSize)
73  C local variables        COMMON /EXF_INTERP_BUFFER/ buffer_r8, buffer_r4
74        integer  ierr, length_of_rec        INTEGER ijs
75        real*8   ne_fac,nw_fac,se_fac,sw_fac  #endif /* ndef EXF_INTERP_USE_DYNALLOC */
76        integer  e_ind(snx,sny),w_ind(snx,sny)  #ifdef ALLOW_USE_MPI
77        integer  n_ind(snx,sny),s_ind(snx,sny)        INTEGER  ierr
       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)  
78  #endif  #endif
79    
80        _BARRIER  C--   Check for consistency:
81        _BEGIN_MASTER( myThid )  #ifdef EXF_INTERP_USE_DYNALLOC
   
82  #ifndef EXF_IREAD_USE_GLOBAL_POINTER  #ifndef EXF_IREAD_USE_GLOBAL_POINTER
83  C     The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the  C     The CPP symbol EXF_IREAD_USE_GLOBAL_POINTER must be defined for the
84  C     case of nThreads > 1. Stop if it isnt.  C     case of nThreads > 1. Stop IF it isnt.
85        IF ( nThreads .GT. 1 ) THEN        IF ( nThreads .GT. 1 ) THEN
86        STOP        STOP
87       &'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'       &'EXF_INTERP_READ: nThreads > 1 needs EXF_IREAD_USE_GLOBAL_POINTER'
88        ENDIF        ENDIF
89  #endif  #endif
90    #else  /* ndef EXF_INTERP_USE_DYNALLOC */
 C read in input data  
 #ifdef ALLOW_USE_MPI  
91  #ifdef EXF_IREAD_USE_GLOBAL_POINTER  #ifdef EXF_IREAD_USE_GLOBAL_POINTER
92         if (.FALSE.) then        STOP
93  #else       &'EXF_INTERP_READ: USE_GLOBAL_POINTER needs INTERP_USE_DYNALLOC'
       if (useSingleCPUIO) then  
94  #endif  #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 master thread of process 0, only, opens a global file  C--   before starting to read, wait for everyone to finish
101          IF( mpiMyId .EQ. 0 ) THEN        _BARRIER
          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)  
         ENDIF  
102    
103  C broadcast to all processes  C---  read in input data
104          IF ( filePrec .EQ. 32 ) THEN  
105             call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,        IF ( MASTER_CPU_IO(myThid) ) THEN
106       &          0,MPI_COMM_MODEL,ierr)  C--   master thread of process 0, only, opens a global file
         ELSE  
            call MPI_BCAST(global8,nx_in*ny_in,MPI_DOUBLE_PRECISION,  
      &          0,MPI_COMM_MODEL,ierr)  
         ENDIF  
        else  
 #endif /* ALLOW_USE_MPI */  
107    
108          call mdsfindunit( interp_unit, mythid)          CALL MDSFINDUNIT( ioUnit, myThid )
109          length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, mythid )          length_of_rec=MDS_RECLEN( filePrec, nx_in*ny_in, myThid )
110          open(interp_unit,file=infile,status='old',access='direct',          OPEN( ioUnit, file=infile, status='old', access='direct',
111       &       recl=length_of_rec)       &       recl=length_of_rec )
112          IF ( filePrec .EQ. 32 ) THEN          IF ( filePrec .EQ. 32 ) THEN
113             read(interp_unit,rec=irecord) global  #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  #ifdef _BYTESWAPIO
119             call MDS_BYTESWAPR4(nx_in*ny_in,global)            CALL MDS_BYTESWAPR4(nx_in*ny_in,buffer_r4)
120  #endif /* _BYTESWAPIO */  #endif /* _BYTESWAPIO */
121          ELSE          ELSE
122             read(interp_unit,rec=irecord) global8  #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  #ifdef _BYTESWAPIO
128             call MDS_BYTESWAPR8(nx_in*ny_in,global8)            CALL MDS_BYTESWAPR8(nx_in*ny_in,buffer_r8)
129  #endif /* _BYTESWAPIO */  #endif /* _BYTESWAPIO */
130          ENDIF          ENDIF
131          close(interp_unit)          CLOSE( ioUnit )
132    C--   end if MASTER_CPU_IO
133          ENDIF
134    
135          _BEGIN_MASTER( myThid )
136  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
137         endif  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 */  #endif /* ALLOW_USE_MPI */
148    
149  #ifdef EXF_IREAD_USE_GLOBAL_POINTER  #ifdef EXF_IREAD_USE_GLOBAL_POINTER
150         IF ( filePrec .EQ. 32 ) THEN         IF ( filePrec .EQ. 32 ) THEN
151            glPtr => global           glPtr4 => buffer_r4
152         ELSE         ELSE
153            glPtr8 => global8           glPtr8 => buffer_r8
154         ENDIF         ENDIF
155  #endif  #endif
156        _END_MASTER( myThid )        _END_MASTER( myThid )
157        _BARRIER        _BARRIER
158    
159    C---  Transfer buffer to "arrayin" array:
160    #ifdef EXF_INTERP_USE_DYNALLOC
161  #ifdef EXF_IREAD_USE_GLOBAL_POINTER  #ifdef EXF_IREAD_USE_GLOBAL_POINTER
162        IF ( filePrec .EQ. 32 ) THEN        IF ( filePrec .EQ. 32 ) THEN
163           do j=1,ny_in          DO j=1,ny_in
164              do i=1,nx_in            DO i=1,nx_in
165                 arrayin(i,j)=glPtr(i,j)              arrayin(i,j)=glPtr4(i,j)
166              enddo            ENDDO
167           enddo          ENDDO
168        ELSE        ELSE
169           do j=1,ny_in          DO j=1,ny_in
170              do i=1,nx_in            DO i=1,nx_in
171                 arrayin(i,j)=glPtr8(i,j)              arrayin(i,j)=glPtr8(i,j)
172              enddo            ENDDO
173           enddo          ENDDO
174        ENDIF        ENDIF
175  #else  #else /* ndef EXF_IREAD_USE_GLOBAL_POINTER */
176        IF ( filePrec .EQ. 32 ) THEN        IF ( filePrec .EQ. 32 ) THEN
177           do j=1,ny_in          DO j=1,ny_in
178              do i=1,nx_in            DO i=1,nx_in
179                 arrayin(i,j)=global(i,j)              arrayin(i,j)=buffer_r4(i,j)
180              enddo            ENDDO
181           enddo          ENDDO
182        ELSE        ELSE
183           do j=1,ny_in          DO j=1,ny_in
184              do i=1,nx_in            DO i=1,nx_in
185                 arrayin(i,j)=global8(i,j)              arrayin(i,j)=buffer_r8(i,j)
186              enddo            ENDDO
187           enddo          ENDDO
188        ENDIF        ENDIF
189  #endif  #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        RETURN
209        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22