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

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

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

revision 1.4 by dimitri, Sun Dec 7 07:41:45 2003 UTC revision 1.12 by heimbach, Fri Jan 13 16:48:03 2006 UTC
# Line 58  C Line 58  C
58    
59  #include "SIZE.h"  #include "SIZE.h"
60  #include "EEPARAMS.h"  #include "EEPARAMS.h"
61  #ifdef ALLOW_USE_MPI  #include "PARAMS.h"
 # include "EESUPPORT.h"  
 # include "PARAMS.h"  
 #endif /* ALLOW_USE_MPI */  
62    
63  C subroutine variables  C subroutine variables
64        character*(*) infile        character*(*) infile
65        integer       filePrec, irecord, nx_in, ny_in        integer       filePrec, irecord, nx_in, ny_in
66        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67        _RS           xG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           xG_in      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
69        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
70        _RL           lat_0, lat_inc(ny_in-1)        _RL           lat_0, lat_inc(ny_in-1)
71        integer       method, mythid        integer       method, mythid
72    
73  C local variables  C local variables
       integer  ierr  
74        real*8   ne_fac,nw_fac,se_fac,sw_fac        real*8   ne_fac,nw_fac,se_fac,sw_fac
75        integer  e_ind(snx),w_ind(snx)        integer  e_ind(snx,sny),w_ind(snx,sny)
76        integer  n_ind(sny),s_ind(sny)        integer  n_ind(snx,sny),s_ind(snx,sny)
77        real*8   px_ind(4), py_ind(4), ew_val(4)        real*8   px_ind(4), py_ind(4), ew_val(4)
78        external lagran        external lagran
79        real*8   lagran        real*8   lagran
80        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
81        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
82        integer  i, j, k, l, js, bi, bj, interp_unit, sp        integer  i, j, k, l, js, bi, bj, sp, interp_unit
83        real*4   global(nx_in,ny_in)        _RS      xg(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
84    
85        _BEGIN_MASTER( myThid )  !     nomalize xg
86          do bj=myByLo(myThid),myByHi(myThid)
87  C check input arguments         do bi=myBxLo(myThid),myBxHi(myThid)
88         if ( .NOT. (filePrec .EQ. 32) )          do j=1-OLy,sNy+OLy
89       &     stop 'stop in exf_interp.F: value of filePrec not allowed'           do i=1-OLx,sNx+OLx
90              xg(i,j,bi,bj) = xg_in(i,j,bi,bj)
91  C read in input data            xg(i,j,bi,bj) = mod(xg(i,j,bi,bj),360.)
92  #ifdef ALLOW_USE_MPI            if ( xg(i,j,bi,bj) .lt. 0. )
93         if (useSingleCPUIO) then       &     xg(i,j,bi,bj) = 360.+xg(i,j,bi,bj)
   
 C master thread of process 0, only, opens a global file  
         IF( mpiMyId .EQ. 0 ) THEN  
          call mdsfindunit( interp_unit, mythid)  
          open(interp_unit,file=infile,status='old',access='direct',  
      &        recl=nx_in*ny_in*4)  
          read(interp_unit,rec=irecord)  
      &        ((global(i,j),i=1,nx_in),j=1,ny_in)  
          close(interp_unit)  
         ENDIF  
   
 C broadcast to all processes  
         call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,  
      &       0,MPI_COMM_MODEL,ierr)  
         do j=1,ny_in  
          do i=1,nx_in  
           arrayin(i,j)=global(i,j)  
94           enddo           enddo
95          enddo          enddo
96           enddo
97          enddo
98    
99         else         call exf_interp_read(
100  #endif /* ALLOW_USE_MPI */       I   infile,
101         I   filePrec,
102         O   arrayin,
103         I   irecord, xG, yG,
104         I   lon_0, lon_inc,
105         I   lat_0, lat_inc,
106         I   nx_in, ny_in, method, mythid)
107          _BARRIER
108    
109          call mdsfindunit( interp_unit, mythid)  C     _BEGIN_MASTER( myThid )
         open(interp_unit,file=infile,status='old',access='direct',  
      &       recl=nx_in*ny_in*4)  
         read(interp_unit,rec=irecord)  
      &       ((arrayin(i,j),i=1,nx_in),j=1,ny_in)  
         close(interp_unit)  
   
 #ifdef ALLOW_USE_MPI  
        endif  
 #endif /* ALLOW_USE_MPI */  
   
 #ifdef _BYTESWAPIO  
        call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )  
 #endif /* _BYTESWAPIO */  
110    
111  C setup input grid  C setup input grid
112         do i=-1,nx_in+2         do i=-1,nx_in+2
113          x_in(i) = lon_0 + (i-1.)*lon_inc          x_in(i) = lon_0 + (i-1.)*lon_inc
114            x_in(i) = mod(x_in(i),360.)
115            if ( x_in(i) .lt. 0. )
116         &   x_in(i) = 360.+x_in(i)
117         enddo         enddo
118    
119         y_in(0) = lat_0 - lat_inc(1)         y_in(0) = lat_0 - lat_inc(1)
120         y_in(-1)= lat_0 - 2.*lat_inc(1)         y_in(-1)= lat_0 - 2.*lat_inc(1)
121         y_in(1) = lat_0         y_in(1) = lat_0
# Line 160  C enlarge boundary Line 139  C enlarge boundary
139          arrayin(i,ny_in+2) = arrayin(i,ny_in)          arrayin(i,ny_in+2) = arrayin(i,ny_in)
140         enddo         enddo
141    
142        _END_MASTER( myThid )  C     _END_MASTER( myThid )
143    
144        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
145         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
146    
147  C check validity of input/output coordinates  C check validity of input/output coordinates
148          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.  #ifdef ALLOW_DEBUG
149       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.          if ( debugLevel .ge. debLevB ) then
150       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.           do i=1,snx
151       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then            do j=1,sny
152             print*,'ERROR in S/R EXF_INTERP:'             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
153             print*,'   input grid must encompass output grid.'       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
154             STOP   ' ABNORMAL END: S/R EXF_INTERP'       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
155         &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
156                  print*,'ERROR in S/R EXF_INTERP:'
157                  print*,'   input grid must encompass output grid.'
158                  print*,'i,j,bi,bj'      ,i,j,bi,bj
159                  print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)
160                  print*,'nx_in,ny_in'    ,nx_in        ,ny_in
161                  print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
162                  print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
163                  STOP   ' ABNORMAL END: S/R EXF_INTERP'
164                 endif
165               enddo
166             enddo
167          endif          endif
168    #endif /* ALLOW_DEBUG */
169    
170  C compute interpolation indices  C compute interpolation indices
171          do i=1,snx          do i=1,snx
172           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then           do j=1,sny
173              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc) + 1            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
174           else             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
175              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)            else
176           endif             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
177           e_ind(i) = w_ind(i) + 1            endif
178          enddo            e_ind(i,j) = w_ind(i,j) + 1
179          js = ny_in/2            js = ny_in*.5
180          do j=1,sny            do while (yG(i,j,bi,bj) .lt. y_in(js))
181           do while (yG(1,j,bi,bj) .lt. y_in(js))             js = (js - 1)*.5
182              js = (js + 1)/2            enddo
183           enddo            do while (yG(i,j,bi,bj) .ge. y_in(js+1))
184           do while (yG(1,j,bi,bj) .ge. y_in(js+1))             js = js + 1
185              js = js + 1            enddo
186              s_ind(i,j) = js
187              n_ind(i,j) = js + 1
188           enddo           enddo
          s_ind(j) = js  
          n_ind(j) = js + 1  
189          enddo          enddo
190    
191          if (method .eq. 1) then          if (method .eq. 1) then
# Line 204  C bilinear interpolation Line 196  C bilinear interpolation
196            do i=1,snx            do i=1,snx
197             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
198             do l=0,1             do l=0,1
199              px_ind(l+1) = x_in(w_ind(i)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
200              py_ind(l+1) = y_in(s_ind(j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
201             enddo             enddo
202             do k=1,2             do k=1,2
203              ew_val(k) = arrayin(w_ind(i),s_ind(j)+k-1)              ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)
204       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
205       &             +arrayin(e_ind(i),s_ind(j)+k-1)       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
206       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
207              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
208       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
209             enddo             enddo
210            enddo            enddo
211           enddo           enddo
# Line 225  C bicubic interpolation Line 217  C bicubic interpolation
217            do i=1,snx            do i=1,snx
218             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
219             do l=-1,2             do l=-1,2
220              px_ind(l+2) = x_in(w_ind(i)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
221              py_ind(l+2) = y_in(s_ind(j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
222             enddo             enddo
223             do k=1,4             do k=1,4
224              ew_val(k) =              ew_val(k) =
225       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
226       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
227       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
228       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
229       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
230       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
231       &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
232       &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
233              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
234       &             +ew_val(k)*lagran(k,yG(1,j,bi,bj),py_ind,sp)       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)
235             enddo             enddo
236            enddo            enddo
237           enddo           enddo
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22