/[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.19 by jmc, Mon Apr 16 23:27:21 2007 UTC revision 1.22 by mlosch, Wed Jan 23 16:41:01 2008 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  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 12  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(
# Line 46  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Line 54  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
54    
55        implicit none        implicit none
56    
57  C     infile       = name of the input file (direct access binary)  C  infile      (string)  :: name of the binary input file (direct access)
58  C     filePrec     = file precicision (currently not used, assumes real*4)  C  filePrec    (integer) :: number of bits per word in file (32 or 64)
59  C     arrout       = output arrays (different for each processor)  C  arrout      ( _RL )   :: output array
60  C     irecord      = record number in global file  C  irecord     (integer) :: record number to read
61  C     xG,yG        = coordinates for output grid  C     xG,yG              :: coordinates for output grid to interpolate to
62  C     lon_0, lat_0 = lon and lat of sw corner of global input grid  C     lon_0, lat_0       :: lon and lat of sw corner of global input grid
63  C     lon_inc      = scalar x-grid increment  C     lon_inc            :: scalar x-grid increment
64  C     lat_inc      = vector y-grid increments  C     lat_inc            :: vector y-grid increments
65  C     nx_in, ny_in = input x-grid and y-grid size  C  nx_in,ny_in (integer) :: size in x & y direction of input file to read
66  C     method       = 1,11,21 for bilinear; 2,12,22 for bicubic  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  C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
68  C     mythid       = thread id  C  myThid      (integer) :: My Thread Id number
69  C  C
70    
71  #include "SIZE.h"  #include "SIZE.h"
# Line 68  C subroutine variables Line 76  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_in      (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 functions
86          external lagran
87          _RL      lagran
88    
89  C local variables  C local variables
90        integer  e_ind(snx,sny),w_ind(snx,sny)        integer  e_ind(snx,sny),w_ind(snx,sny)
91        integer  n_ind(snx,sny),s_ind(snx,sny)        integer  n_ind(snx,sny),s_ind(snx,sny)
92        real*8   px_ind(4), py_ind(4), ew_val(4)        _RL      px_ind(4), py_ind(4), ew_val(4)
93        external lagran        _RL      arrayin(-1:nx_in+2 ,      -1:ny_in+2)
94        real*8   lagran        _RL      NorthValue
95        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        _RL      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)  
       real*8   ninety  
       PARAMETER ( ninety = 90. )  
96        integer  i, j, k, l, js, bi, bj, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
97    #ifdef TARGET_NEC_SX
98          integer  ic, ii, icnt
99          integer  inx(snx*sny,2)
100    #endif
101        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
102        _RS      threeSixtyRS, NorthValue        _RL      ninety
103          PARAMETER ( ninety = 90. )
104          _RS      threeSixtyRS
105        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
106    
107  C     put xG in interval [ lon_0 , lon_0+360 [  C     put xG in interval [ lon_0 , lon_0+360 [
# Line 106  C     put xG in interval [ lon_0 , lon_0 Line 121  C     put xG in interval [ lon_0 , lon_0
121       I   infile, filePrec,       I   infile, filePrec,
122       O   arrayin,       O   arrayin,
123       I   irecord, nx_in, ny_in, mythid)       I   irecord, nx_in, ny_in, mythid)
       _BARRIER  
   
 C     _BEGIN_MASTER( myThid )  
124    
125  C setup input longitude grid  C setup input longitude grid
126        do i=-1,nx_in+2        do i=-1,nx_in+2
# Line 141  C enlarge boundary Line 153  C enlarge boundary
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        enddo
158    
159  C     For tracer (method=1,2) set to northernmost zonal-mean value  C     For tracer (method=1,2) set to northernmost zonal-mean value
# Line 153  C     as is already done above in order Line 165  C     as is already done above in order
165        do j=ny_in,ny_in+2        do j=ny_in,ny_in+2
166         if (y_in(j).eq.ninety) then         if (y_in(j).eq.ninety) then
167          if (method.eq.1 .or. method.eq.2) then          if (method.eq.1 .or. method.eq.2) then
168           NorthValue = 0           NorthValue = 0.
169           do i=1,nx_in           do i=1,nx_in
170            NorthValue = NorthValue + arrayin(i,j)            NorthValue = NorthValue + arrayin(i,j)
171           enddo           enddo
# Line 163  C     as is already done above in order Line 175  C     as is already done above in order
175           enddo           enddo
176          elseif (method.eq.11 .or. method.eq.12) then          elseif (method.eq.11 .or. method.eq.12) then
177           do i=-1,nx_in+2           do i=-1,nx_in+2
178            arrayin(i,j) = 0            arrayin(i,j) = 0.
179           enddo           enddo
180          endif          endif
181         endif         endif
182        enddo        enddo
183    
 C     _END_MASTER( myThid )  
         
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    
# Line 197  C check validity of input/output coordin Line 207  C check validity of input/output coordin
207          endif          endif
208  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
209    
210  C compute interpolation indices  C compute interpolation indices
211          do i=1,snx          do i=1,snx
212           do j=1,sny           do j=1,sny
213            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then
# Line 206  C compute interpolation indices Line 216  C compute interpolation indices
216             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
217            endif            endif
218            e_ind(i,j) = w_ind(i,j) + 1            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            js = ny_in*.5
227            do while (yG(i,j,bi,bj) .lt. y_in(js))            do while (yG(i,j,bi,bj) .lt. y_in(js))
228             js = (js - 1)*.5             js = (js - 1)*.5
# Line 214  C compute interpolation indices Line 231  C compute interpolation indices
231             js = js + 1             js = js + 1
232            enddo            enddo
233            s_ind(i,j) = js            s_ind(i,j) = js
234            n_ind(i,j) = js + 1           enddo
235            enddo
236    #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
242             do i=1,snx
243              s_ind(i,j) = ny_in*.5
244              icnt = icnt+1
245              inx(icnt,1) = i
246              inx(icnt,2) = j
247             enddo
248            enddo
249            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
292          enddo          enddo
293    
# Line 257  C bicubic interpolation Line 330  C bicubic interpolation
330       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
331       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)
332       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)
333       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)
334       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)
335       &             *lagran(4,xG(i,j,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(i,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
# Line 271  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.19  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22