/[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.3 by edhill, Thu Oct 9 04:19:19 2003 UTC revision 1.10 by heimbach, Mon May 2 15:31:46 2005 UTC
# Line 41  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 41  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
41       I   lat_0, lat_inc,       I   lat_0, lat_inc,
42       I   nx_in, ny_in, method, mythid)       I   nx_in, ny_in, method, mythid)
43    
44  C        implicit none
45    
46  C     infile       = name of the input file (direct access binary)  C     infile       = name of the input file (direct access binary)
47  C     filePrec     = file precicision (currently not used, assumes real*4)  C     filePrec     = file precicision (currently not used, assumes real*4)
48  C     arrout       = output arrays (different for each processor)  C     arrout       = output arrays (different for each processor)
# Line 57  C Line 58  C
58    
59  #include "SIZE.h"  #include "SIZE.h"
60  #include "EEPARAMS.h"  #include "EEPARAMS.h"
61    #include "PARAMS.h"
62    
63  C subroutine variables  C subroutine variables
64        character*(*) infile        character*(*) infile
# Line 69  C subroutine variables Line 71  C subroutine variables
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
# Line 80  C local variables Line 81  C local variables
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, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
83    
84  C check input arguments         call exf_interp_read(
85        if ( .NOT. (filePrec .EQ. 32) )       I   infile,
86       & stop 'stop in exf_interp.F: value of filePrec not allowed'       I   filePrec,
87         O   arrayin,
88  C read in input data       I   irecord, xG, yG,
89        call mdsfindunit( interp_unit, mythid)       I   lon_0, lon_inc,
90        open(interp_unit,file=infile,status='old',access='direct',       I   lat_0, lat_inc,
91       & recl=nx_in*ny_in*4)       I   nx_in, ny_in, method, mythid)
92        read(interp_unit,rec=irecord) ((arrayin(i,j),i=1,nx_in),j=1,ny_in)  
93  #ifdef _BYTESWAPIO        _BEGIN_MASTER( myThid )
       call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )  
 #endif  
       close(interp_unit)  
94    
95  C setup input grid  C setup input grid
96        do i=-1,nx_in+2         do i=-1,nx_in+2
97         x_in(i) = lon_0 + (i-1.)*lon_inc          x_in(i) = lon_0 + (i-1.)*lon_inc
98        enddo         enddo
99        y_in(0) = lat_0 - lat_inc(1)         y_in(0) = lat_0 - lat_inc(1)
100        y_in(-1)= lat_0 - 2.*lat_inc(1)         y_in(-1)= lat_0 - 2.*lat_inc(1)
101        y_in(1) = lat_0         y_in(1) = lat_0
102        do j=2,ny_in         do j=2,ny_in
103        y_in(j) = y_in(j-1) + lat_inc(j-1)          y_in(j) = y_in(j-1) + lat_inc(j-1)
104        enddo         enddo
105        y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)         y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)
106        y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)         y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)
107    
108  C enlarge boundary  C enlarge boundary
109        do j=1,ny_in         do j=1,ny_in
110           arrayin(0,j)       = arrayin(nx_in,j)          arrayin(0,j)       = arrayin(nx_in,j)
111           arrayin(-1,j)      = arrayin(nx_in-1,j)          arrayin(-1,j)      = arrayin(nx_in-1,j)
112           arrayin(nx_in+1,j) = arrayin(1,j)          arrayin(nx_in+1,j) = arrayin(1,j)
113           arrayin(nx_in+2,j) = arrayin(2,j)          arrayin(nx_in+2,j) = arrayin(2,j)
114        enddo         enddo
115        do i=-1,nx_in+2         do i=-1,nx_in+2
116           arrayin(i,0)       = arrayin(i,1)          arrayin(i,0)       = arrayin(i,1)
117           arrayin(i,-1)      = arrayin(i,1)          arrayin(i,-1)      = arrayin(i,1)
118           arrayin(i,ny_in+1) = arrayin(i,ny_in)          arrayin(i,ny_in+1) = arrayin(i,ny_in)
119           arrayin(i,ny_in+2) = arrayin(i,ny_in)          arrayin(i,ny_in+2) = arrayin(i,ny_in)
120        enddo         enddo
121    
122          _END_MASTER( myThid )
123    
124        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
125         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
126    
127  C check validity of input/output coordinates  C check validity of input/output coordinates
128          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.  #ifdef ALLOW_DEBUG
129       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.          if ( debugLevel .ge. debLevB ) then
130       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.           do i=1,snx
131       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then            do j=1,sny
132             print*,'ERROR in S/R EXF_INTERP:'             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
133             print*,'   input grid must encompass output grid.'       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
134             STOP   ' ABNORMAL END: S/R EXF_INTERP'       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
135         &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
136                  print*,'ERROR in S/R EXF_INTERP:'
137                  print*,'   input grid must encompass output grid.'
138                  print*,'i,j,bi,bj'      ,i,j,bi,bj
139                  print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)
140                  print*,'nx_in,ny_in'    ,nx_in        ,ny_in
141                  print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
142                  print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
143                  STOP   ' ABNORMAL END: S/R EXF_INTERP'
144                 endif
145               enddo
146             enddo
147          endif          endif
148    #endif /* ALLOW_DEBUG */
149    
150  C compute interpolation indices  C compute interpolation indices
151          do i=1,snx          do i=1,snx
152           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then           do j=1,sny
153              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
154           else             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
155              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)            else
156           endif             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
157           e_ind(i) = w_ind(i) + 1            endif
158          enddo            e_ind(i,j) = w_ind(i,j) + 1
159          js = ny_in/2            js = ny_in*.5
160          do j=1,sny            do while (yG(i,j,bi,bj) .lt. y_in(js))
161           do while (yG(1,j,bi,bj) .lt. y_in(js))             js = (js - 1)*.5
162              js = (js + 1)/2            enddo
163           enddo            do while (yG(i,j,bi,bj) .ge. y_in(js+1))
164           do while (yG(1,j,bi,bj) .ge. y_in(js+1))             js = js + 1
165              js = js + 1            enddo
166              s_ind(i,j) = js
167              n_ind(i,j) = js + 1
168           enddo           enddo
          s_ind(j) = js  
          n_ind(j) = js + 1  
169          enddo          enddo
170    
171          if (method .eq. 1) then          if (method .eq. 1) then
# Line 163  C bilinear interpolation Line 176  C bilinear interpolation
176            do i=1,snx            do i=1,snx
177             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
178             do l=0,1             do l=0,1
179              px_ind(l+1) = x_in(w_ind(i)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
180              py_ind(l+1) = y_in(s_ind(j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
181             enddo             enddo
182             do k=1,2             do k=1,2
183              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)
184       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
185       &             +arrayin(e_ind(i),s_ind(j)+k-1)       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
186       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
187              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
188       &             +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)
189             enddo             enddo
190            enddo            enddo
191           enddo           enddo
# Line 184  C bicubic interpolation Line 197  C bicubic interpolation
197            do i=1,snx            do i=1,snx
198             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
199             do l=-1,2             do l=-1,2
200              px_ind(l+2) = x_in(w_ind(i)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
201              py_ind(l+2) = y_in(s_ind(j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
202             enddo             enddo
203             do k=1,4             do k=1,4
204              ew_val(k) =              ew_val(k) =
205       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
206       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
207       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
208       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
209       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
210       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
211       &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
212       &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
213              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
214       &             +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)
215             enddo             enddo
216            enddo            enddo
217           enddo           enddo
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22