/[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.22 by mlosch, Wed Jan 23 16:41:01 2008 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
7  C Flux Coupler using                       C  C Flux Coupler using                       C
8  C Bilinear interpolation of forcing fields C  C Bilinear interpolation of forcing fields C
# Line 9  C added Bicubic (bnc 1/2003) Line 13  C added Bicubic (bnc 1/2003)
13  C                                          C  C                                          C
14  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15    
16          real*8 function lagran(i,x,a,sp)         _RL FUNCTION LAGRAN(i,x,a,sp)
17    
18          INTEGER i,k,sp          INTEGER i
19          _RS x          _RS x
20          real*8 a(4)          _RL a(4)
21          real*8 numer,denom          INTEGER sp
   
         numer = 1.D0  
         denom = 1.D0  
22    
23    C-      local variables:
24            INTEGER k
25            _RL numer,denom
26    
27            numer = 1. _d 0
28            denom = 1. _d 0
29    
30    #ifdef TARGET_NEC_SX
31    !CDIR UNROLL=8
32    #endif /* TARGET_NEC_SX */
33          do k=1,sp          do k=1,sp
34          if ( k .ne. i) then           if ( k .ne. i) then
35            denom = denom*(a(i) - a(k))            denom = denom*(a(i) - a(k))
36            numer = numer*(x    - a(k))            numer = numer*(x    - a(k))
37          endif           endif
38          enddo          enddo
39    
40          lagran = numer/denom          lagran = numer/denom
41    
42          return         RETURN
43          end         END
44    
45    
46         SUBROUTINE exf_interp(         SUBROUTINE exf_interp(
47       I   infile,       I   infile,
48       I   filePrec,       I   filePrec,
49       O   arrayout,       O   arrayout,
50       I   irecord, xG, yG,       I   irecord, xG_in, yG,
51       I   lon_0, lon_inc,       I   lon_0, lon_inc,
52       I   lat_0, lat_inc,       I   lat_0, lat_inc,
53       I   nx_in, ny_in, method, mythid)       I   nx_in, ny_in, method, mythid)
54    
55  C        implicit none
56  C     infile       = name of the input file (direct access binary)  
57  C     filePrec     = file precicision (currently not used, assumes real*4)  C  infile      (string)  :: name of the binary input file (direct access)
58  C     arrout       = output arrays (different for each processor)  C  filePrec    (integer) :: number of bits per word in file (32 or 64)
59  C     irecord      = record number in global file  C  arrout      ( _RL )   :: output array
60  C     xG,yG        = coordinates for output grid  C  irecord     (integer) :: record number to read
61  C     lon_0, lat_0 = lon and lat of sw corner of global input grid  C     xG,yG              :: coordinates for output grid to interpolate to
62  C     lon_inc      = scalar x-grid increment  C     lon_0, lat_0       :: lon and lat of sw corner of global input grid
63  C     lat_inc      = vector y-grid increments  C     lon_inc            :: scalar x-grid increment
64  C     nx_in, ny_in = input x-grid and y-grid size  C     lat_inc            :: vector y-grid increments
65  C     method       = 1 for bilinear 2 for bicubic  C  nx_in,ny_in (integer) :: size in x & y direction of input file to read
66  C     mythid       = thread id  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
67    C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
68    C  myThid      (integer) :: My Thread Id number
69  C  C
70    
71  #include "SIZE.h"  #include "SIZE.h"
72  #include "EEPARAMS.h"  #include "EEPARAMS.h"
73    #include "PARAMS.h"
74    
75  C subroutine variables  C subroutine variables
76        character*(*) infile        character*(*) infile
77        integer       filePrec, irecord, nx_in, ny_in        integer       filePrec, irecord, nx_in, ny_in
78        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79        _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)
80        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
82        _RL           lat_0, lat_inc(ny_in-1)        _RL           lat_0, lat_inc(ny_in-1)
83        integer       method, mythid        integer       method, mythid
84    
85  C local variables  C functions
       integer  ierr  
       real*8   ne_fac,nw_fac,se_fac,sw_fac  
       integer  e_ind(snx),w_ind(snx)  
       integer  n_ind(sny),s_ind(sny)  
       real*8   px_ind(4), py_ind(4), ew_val(4)  
86        external lagran        external lagran
87        real*8   lagran        _RL      lagran
       real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)  
       real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)  
       integer  i, j, k, l, js, bi, bj, sp, interp_unit  
88    
89  C check input arguments  C local variables
90        if ( .NOT. (filePrec .EQ. 32) )        integer  e_ind(snx,sny),w_ind(snx,sny)
91       & stop 'stop in exf_interp.F: value of filePrec not allowed'        integer  n_ind(snx,sny),s_ind(snx,sny)
92          _RL      px_ind(4), py_ind(4), ew_val(4)
93  C read in input data        _RL      arrayin(-1:nx_in+2 ,      -1:ny_in+2)
94        call mdsfindunit( interp_unit, mythid)        _RL      NorthValue
95        open(interp_unit,file=infile,status='old',access='direct',        _RL      x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
96       & recl=nx_in*ny_in*4)        integer  i, j, k, l, js, bi, bj, sp, interp_unit
97        read(interp_unit,rec=irecord) ((arrayin(i,j),i=1,nx_in),j=1,ny_in)  #ifdef TARGET_NEC_SX
98  #ifdef _BYTESWAPIO        integer  ic, ii, icnt
99        call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )        integer  inx(snx*sny,2)
100  #endif  #endif
101        close(interp_unit)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
102          _RL      ninety
103          PARAMETER ( ninety = 90. )
104          _RS      threeSixtyRS
105          PARAMETER ( threeSixtyRS = 360. )
106    
107    C     put xG in interval [ lon_0 , lon_0+360 [
108          do bj=myByLo(myThid),myByHi(myThid)
109           do bi=myBxLo(myThid),myBxHi(myThid)
110            do j=1-OLy,sNy+OLy
111             do i=1-OLx,sNx+OLx
112              xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
113         &                  + threeSixtyRS*2.
114              xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
115             enddo
116            enddo
117           enddo
118          enddo
119    
120  C setup input grid         call exf_interp_read(
121         I   infile, filePrec,
122         O   arrayin,
123         I   irecord, nx_in, ny_in, mythid)
124    
125    C setup input longitude grid
126        do i=-1,nx_in+2        do i=-1,nx_in+2
127         x_in(i) = lon_0 + (i-1.)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
128        enddo        enddo
129    
130    C setup input latitude grid
131        y_in(0) = lat_0 - lat_inc(1)        y_in(0) = lat_0 - lat_inc(1)
132        y_in(-1)= lat_0 - 2.*lat_inc(1)        y_in(-1)= lat_0 - 2.*lat_inc(1)
133        y_in(1) = lat_0        y_in(1) = lat_0
134        do j=2,ny_in        do j=2,ny_in
135        y_in(j) = y_in(j-1) + lat_inc(j-1)         y_in(j) = y_in(j-1) + lat_inc(j-1)
136          enddo
137          do j=ny_in+1,ny_in+2
138           if (y_in(j-1).eq.ninety) then
139            y_in(j) = 2 * ninety - y_in(j-2)
140           else
141            y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
142           endif
143        enddo        enddo
       y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)  
       y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)  
144    
145  C enlarge boundary  C enlarge boundary
146        do j=1,ny_in        do j=1,ny_in
147           arrayin(0,j)       = arrayin(nx_in,j)         arrayin(0,j)       = arrayin(nx_in,j)
148           arrayin(-1,j)      = arrayin(nx_in-1,j)         arrayin(-1,j)      = arrayin(nx_in-1,j)
149           arrayin(nx_in+1,j) = arrayin(1,j)         arrayin(nx_in+1,j) = arrayin(1,j)
150           arrayin(nx_in+2,j) = arrayin(2,j)         arrayin(nx_in+2,j) = arrayin(2,j)
151        enddo        enddo
152        do i=-1,nx_in+2        do i=-1,nx_in+2
153           arrayin(i,0)       = arrayin(i,1)         arrayin(i,0)       = arrayin(i,1)
154           arrayin(i,-1)      = arrayin(i,1)         arrayin(i,-1)      = arrayin(i,1)
155           arrayin(i,ny_in+1) = arrayin(i,ny_in)         arrayin(i,ny_in+1) = arrayin(i,ny_in)
156           arrayin(i,ny_in+2) = arrayin(i,ny_in)         arrayin(i,ny_in+2) = arrayin(i,ny_in)
157          enddo
158    
159    C     For tracer (method=1,2) set to northernmost zonal-mean value
160    C     at 90N to avoid sharp zonal gradients near the Pole.
161    C     For U (method=11,12) set to zero at 90N to minimize velocity
162    C     gradient at North Pole
163    C     For V (method=11,12) set to northernmost zonal value at 90N,
164    C     as is already done above in order to allow cross-PoleArctic flow
165          do j=ny_in,ny_in+2
166           if (y_in(j).eq.ninety) then
167            if (method.eq.1 .or. method.eq.2) then
168             NorthValue = 0.
169             do i=1,nx_in
170              NorthValue = NorthValue + arrayin(i,j)
171             enddo
172             NorthValue = NorthValue / nx_in
173             do i=-1,nx_in+2
174              arrayin(i,j) = NorthValue
175             enddo
176            elseif (method.eq.11 .or. method.eq.12) then
177             do i=-1,nx_in+2
178              arrayin(i,j) = 0.
179             enddo
180            endif
181           endif
182        enddo        enddo
183    
184        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
185         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
186    
187  C check validity of input/output coordinates  C check validity of input/output coordinates
188          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.  #ifdef ALLOW_DEBUG
189       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.          if ( debugLevel .ge. debLevB ) then
190       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.           do i=1,snx
191       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then            do j=1,sny
192             print*,'ERROR in S/R EXF_INTERP:'             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
193             print*,'   input grid must encompass output grid.'       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
194             STOP   ' ABNORMAL END: S/R EXF_INTERP'       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
195         &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
196                  print*,'ERROR in S/R EXF_INTERP:'
197                  print*,'   input grid must encompass output grid.'
198                  print*,'i,j,bi,bj'      ,i,j,bi,bj
199                  print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)
200                  print*,'nx_in,ny_in'    ,nx_in        ,ny_in
201                  print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
202                  print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
203                  STOP   ' ABNORMAL END: S/R EXF_INTERP'
204                 endif
205               enddo
206             enddo
207          endif          endif
208    #endif /* ALLOW_DEBUG */
209    
210  C compute interpolation indices  C compute interpolation indices
211          do i=1,snx          do i=1,snx
212           if (xG(i,1,bi,bj)-x_in(1) .ge. 0.) then           do j=1,sny
213              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
214           else             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
215              w_ind(i) = int((xG(i,1,bi,bj)-x_in(1))/lon_inc)            else
216           endif             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
217           e_ind(i) = w_ind(i) + 1            endif
218              e_ind(i,j) = w_ind(i,j) + 1
219             enddo
220            enddo
221    #ifndef TARGET_NEC_SX
222    C     use the original and more readable variant of the algorithm that
223    C     has unvectorizable while-loops for each (i,j)
224            do i=1,snx
225             do j=1,sny
226              js = ny_in*.5
227              do while (yG(i,j,bi,bj) .lt. y_in(js))
228               js = (js - 1)*.5
229              enddo
230              do while (yG(i,j,bi,bj) .ge. y_in(js+1))
231               js = js + 1
232              enddo
233              s_ind(i,j) = js
234             enddo
235          enddo          enddo
236          js = ny_in/2  #else /* TARGET_NEC_SX defined */
237    C     this variant vectorizes more efficiently than the original one because
238    C     it moves the while loops out of the i,j loops (loop pushing) but
239    C     it is ugly and incomprehensible
240            icnt = 0
241          do j=1,sny          do j=1,sny
242           do while (yG(1,j,bi,bj) .lt. y_in(js))           do i=1,snx
243              js = (js + 1)/2            s_ind(i,j) = ny_in*.5
244              icnt = icnt+1
245              inx(icnt,1) = i
246              inx(icnt,2) = j
247           enddo           enddo
248           do while (yG(1,j,bi,bj) .ge. y_in(js+1))          enddo
249              js = js + 1          do while (icnt .gt. 0)
250             ii = 0
251    !CDIR NODEP
252             do ic=1,icnt
253              i = inx(ic,1)
254              j = inx(ic,2)
255              if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then
256               s_ind(i,j) = (s_ind(i,j) - 1)*.5
257               ii = ii+1
258               inx(ii,1) = i
259               inx(ii,2) = j
260              endif
261             enddo
262             icnt = ii
263            enddo
264            icnt = 0
265            do j=1,sny
266             do i=1,snx
267              icnt = icnt+1
268              inx(icnt,1) = i
269              inx(icnt,2) = j
270             enddo
271            enddo
272            do while (icnt .gt. 0)
273             ii = 0
274    !CDIR NODEP
275             do ic=1,icnt
276              i = inx(ic,1)
277              j = inx(ic,2)
278              if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then
279               s_ind(i,j) = s_ind(i,j) + 1
280               ii = ii+1
281               inx(ii,1) = i
282               inx(ii,2) = j
283              endif
284             enddo
285             icnt = ii
286            enddo
287    #endif /* TARGET_NEC_SX defined */
288            do i=1,snx
289             do j=1,sny
290              n_ind(i,j) = s_ind(i,j) + 1
291           enddo           enddo
          s_ind(j) = js  
          n_ind(j) = js + 1  
292          enddo          enddo
293    
294          if (method .eq. 1) then          if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then
295    
296  C bilinear interpolation  C bilinear interpolation
297           sp = 2           sp = 2
# Line 163  C bilinear interpolation Line 299  C bilinear interpolation
299            do i=1,snx            do i=1,snx
300             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
301             do l=0,1             do l=0,1
302              px_ind(l+1) = x_in(w_ind(i)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
303              py_ind(l+1) = y_in(s_ind(j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
304             enddo             enddo
305             do k=1,2             do k=1,2
306              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)
307       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
308       &             +arrayin(e_ind(i),s_ind(j)+k-1)       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)
309       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
310              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
311       &             +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)
312             enddo             enddo
313            enddo            enddo
314           enddo           enddo
315          elseif (method .eq. 2) then          elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then
316    
317  C bicubic interpolation  C bicubic interpolation
318           sp = 4           sp = 4
# Line 184  C bicubic interpolation Line 320  C bicubic interpolation
320            do i=1,snx            do i=1,snx
321             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
322             do l=-1,2             do l=-1,2
323              px_ind(l+2) = x_in(w_ind(i)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
324              py_ind(l+2) = y_in(s_ind(j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
325             enddo             enddo
326             do k=1,4             do k=1,4
327              ew_val(k) =              ew_val(k) =
328       &             arrayin(w_ind(i)-1,s_ind(j)+k-2)       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
329       &             *lagran(1,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)
330       &             +arrayin(w_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
331       &             *lagran(2,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
332       &             +arrayin(e_ind(i)  ,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
333       &             *lagran(3,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
334       &             +arrayin(e_ind(i)+1,s_ind(j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
335       &             *lagran(4,xG(i,1,bi,bj),px_ind,sp)       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)
336              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)
337       &             +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)
338             enddo             enddo
339            enddo            enddo
340           enddo           enddo
# Line 208  C bicubic interpolation Line 344  C bicubic interpolation
344         enddo         enddo
345        enddo        enddo
346    
347          RETURN
348        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22