/[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.23 by mlosch, Thu Jan 24 08:29:51 2008 UTC revision 1.28 by jmc, Sun Jan 1 15:25:23 2012 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    #undef OLD_EXF_INTERP_LAT_INDEX
6    
7  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  C==========================================*
8  C Flux Coupler using                       C  C Flux Coupler using                       |
9  C Bilinear interpolation of forcing fields C  C Bilinear interpolation of forcing fields |
10  C                                          C  C                                          |
11  C B. Cheng (12/2002)                       C  C B. Cheng (12/2002)                       |
12  C                                          C  C                                          |
13  C added Bicubic (bnc 1/2003)               C  C added Bicubic (bnc 1/2003)               |
14  C                                          C  C                                          |
15  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  C==========================================*
16    
17    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
18    
19         _RL FUNCTION LAGRAN(i,x,a,sp)         _RL FUNCTION LAGRAN(i,x,a,sp)
20    
21           IMPLICIT NONE
22    
23          INTEGER i          INTEGER i
24          _RS x          _RS x
25          _RL a(4)          _RL a(4)
# Line 30  C-      local variables: Line 35  C-      local variables:
35  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
36  !CDIR UNROLL=8  !CDIR UNROLL=8
37  #endif /* TARGET_NEC_SX */  #endif /* TARGET_NEC_SX */
38          do k=1,sp          DO k=1,sp
39           if ( k .ne. i) then           IF ( k .NE. i) THEN
40            denom = denom*(a(i) - a(k))            denom = denom*(a(i) - a(k))
41            numer = numer*(x    - a(k))            numer = numer*(x    - a(k))
42           endif           ENDIF
43          enddo          ENDDO
44    
45          lagran = numer/denom          LAGRAN = numer/denom
46    
47         RETURN         RETURN
48         END         END
49    
50    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
51    
52         SUBROUTINE exf_interp(  CBOP
53       I   infile,  C     !ROUTINE: EXF_INTERP
54       I   filePrec,  C     !INTERFACE:
55       O   arrayout,         SUBROUTINE EXF_INTERP(
56       I   irecord, xG_in, yG,       I                inFile,
57       I   lon_0, lon_inc,       I                filePrec,
58       I   lat_0, lat_inc,       O                arrayout,
59       I   nx_in, ny_in, method, mythid)       I                irecord, xG_in, yG,
60         I                lon_0, lon_inc,
61         I                lat_0, lat_inc,
62         I                nxIn, nyIn, method, myThid )
63    
64    C !DESCRIPTION: \bv
65    C  *==========================================================*
66    C  | SUBROUTINE EXF_INTERP
67    C  | o Load from file a regular lat-lon input field
68    C  |   and interpolate on to the model grid location
69    C  *==========================================================*
70    C \ev
71    
72    C !USES:
73          IMPLICIT NONE
74    C     === Global variables ===
75    #include "SIZE.h"
76    #include "EEPARAMS.h"
77    #include "PARAMS.h"
78    
       implicit none  
79    
80  C  infile      (string)  :: name of the binary input file (direct access)  C !INPUT/OUTPUT PARAMETERS:
81    C  inFile      (string)  :: name of the binary input file (direct access)
82  C  filePrec    (integer) :: number of bits per word in file (32 or 64)  C  filePrec    (integer) :: number of bits per word in file (32 or 64)
83  C  arrout      ( _RL )   :: output array  C  arrayout    ( _RL )   :: output array
84  C  irecord     (integer) :: record number to read  C  irecord     (integer) :: record number to read
85  C     xG,yG              :: coordinates for output grid to interpolate to  C     xG_in,yG           :: coordinates for output grid to interpolate to
86  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
87  C     lon_inc            :: scalar x-grid increment  C     lon_inc            :: scalar x-grid increment
88  C     lat_inc            :: vector y-grid increments  C     lat_inc            :: vector y-grid increments
89  C  nx_in,ny_in (integer) :: size in x & y direction of input file to read  C  nxIn,nyIn   (integer) :: size in x & y direction of input file to read
90  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
91  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
92  C  myThid      (integer) :: My Thread Id number  C  myThid      (integer) :: My Thread Id number
 C  
   
 #include "SIZE.h"  
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
93    
94  C subroutine variables        CHARACTER*(*) infile
95        character*(*) infile        INTEGER       filePrec, irecord, nxIn, nyIn
       integer       filePrec, irecord, nx_in, ny_in  
96        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
97        _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)
98        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
99        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
100        _RL           lat_0, lat_inc(ny_in-1)  c     _RL           lat_0, lat_inc(nyIn-1)
101        integer       method, mythid        _RL           lat_0, lat_inc(*)
102          INTEGER       method, myThid
103  C functions  
104        external lagran  C !FUNCTIONS:
105        _RL      lagran        EXTERNAL LAGRAN
106          _RL      LAGRAN
107  C local variables        INTEGER  ILNBLNK
108        integer  e_ind(snx,sny),w_ind(snx,sny)        EXTERNAL ILNBLNK
109        integer  n_ind(snx,sny),s_ind(snx,sny)  
110    C !LOCAL VARIABLES:
111    C     msgBuf      :: Informational/error message buffer
112    C     bi, bj      :: tile indices
113          CHARACTER*(MAX_LEN_MBUF) msgBuf
114          INTEGER  bi, bj
115          INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)
116        _RL      px_ind(4), py_ind(4), ew_val(4)        _RL      px_ind(4), py_ind(4), ew_val(4)
117        _RL      arrayin(-1:nx_in+2 ,      -1:ny_in+2)        _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )
118          _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)
119        _RL      NorthValue        _RL      NorthValue
120        _RL      x_in   (-1:nx_in+2), y_in(-1:ny_in+2)        INTEGER  i, j, k, l, sp
121        integer  i, j, k, l, js, bi, bj, sp, interp_unit  #ifdef OLD_EXF_INTERP_LAT_INDEX
122          INTEGER  js
123    #else
124          INTEGER  nLoop
125          _RL      tmpVar
126    #endif
127  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
128        integer  ic, ii, icnt        INTEGER  ic, ii, icnt
129        integer  inx(snx*sny,2)        INTEGER  inx(sNx*sNy,2)
130        _RL      ew_val1, ew_val2, ew_val3, ew_val4        _RL      ew_val1, ew_val2, ew_val3, ew_val4
131  #endif  #endif
132        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
# Line 104  C local variables Line 134  C local variables
134        PARAMETER ( ninety = 90. )        PARAMETER ( ninety = 90. )
135        _RS      threeSixtyRS        _RS      threeSixtyRS
136        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
137          LOGICAL  xIsPeriodic
138    CEOP
139    
140  C     put xG in interval [ lon_0 , lon_0+360 [  C--   put xG in interval [ lon_0 , lon_0+360 [
141        do bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
142         do bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
143          do j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
144           do i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
145            xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0            xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
146       &                  + threeSixtyRS*2.       &                  + threeSixtyRS*2.
147            xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)            xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)
148           enddo           ENDDO
149          enddo          ENDDO
150         enddo         ENDDO
151        enddo        ENDDO
152    
153         call exf_interp_read(  C--   Load inut field
154       I   infile, filePrec,        CALL EXF_INTERP_READ(
155       O   arrayin,       I         inFile, filePrec,
156       I   irecord, nx_in, ny_in, mythid)       O         arrayin,
157         I         irecord, nxIn, nyIn, myThid )
158    
159  C setup input longitude grid  C--   setup input longitude grid
160        do i=-1,nx_in+2        DO i=-1,nxIn+2
161         x_in(i) = lon_0 + (i-1)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
162        enddo        ENDDO
163          xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
164    
165  C setup input latitude grid  C--   setup input latitude grid
       y_in(0) = lat_0 - lat_inc(1)  
       y_in(-1)= lat_0 - 2.*lat_inc(1)  
166        y_in(1) = lat_0        y_in(1) = lat_0
167        do j=2,ny_in        DO j=1,nyIn+1
168         y_in(j) = y_in(j-1) + lat_inc(j-1)         i = MIN(j,nyIn-1)
169        enddo         y_in(j+1) = y_in(j) + lat_inc(i)
170        do j=ny_in+1,ny_in+2        ENDDO
171         if (y_in(j-1).eq.ninety) then  C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
172          y_in(j) = 2 * ninety - y_in(j-2)        y_in(0) = y_in(1) - lat_inc(1)
173         else        y_in(-1)= y_in(0) - lat_inc(1)
174          y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )  c     IF ( y_in(1).GT.-ninety .AND. y_in(0).LT.-ninety ) THEN
175         endif  c       y_in(0) = -ninety
176        enddo  c       y_in(-1) = -2.*ninety - y_in(1)
177    c     ENDIF
178  C enlarge boundary  c     IF ( y_in(0).GT.-ninety .AND. y_in(-1).LT.-ninety ) THEN
179        do j=1,ny_in  c       y_in(-1) = -ninety
180         arrayin(0,j)       = arrayin(nx_in,j)  c     ENDIF
181         arrayin(-1,j)      = arrayin(nx_in-1,j)  C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
182         arrayin(nx_in+1,j) = arrayin(1,j)        j = nyIn+1
183         arrayin(nx_in+2,j) = arrayin(2,j)        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN
184        enddo          y_in(j) = ninety
185        do i=-1,nx_in+2          y_in(j+1) = 2.*ninety - y_in(j-1)
186         arrayin(i,0)       = arrayin(i,1)        ENDIF
187         arrayin(i,-1)      = arrayin(i,1)        j = nyIn+2
188         arrayin(i,ny_in+1) = arrayin(i,ny_in)        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN
189         arrayin(i,ny_in+2) = arrayin(i,ny_in)          y_in(j) = ninety
190        enddo        ENDIF
191    
192    C--   enlarge boundary
193          IF ( xIsPeriodic ) THEN
194            DO j=1,nyIn
195             arrayin( 0,j)     = arrayin(nxIn  ,j)
196             arrayin(-1,j)     = arrayin(nxIn-1,j)
197             arrayin(nxIn+1,j) = arrayin(1,j)
198             arrayin(nxIn+2,j) = arrayin(2,j)
199            ENDDO
200          ELSE
201            DO j=1,nyIn
202             arrayin( 0,j)     = arrayin(1,j)
203             arrayin(-1,j)     = arrayin(1,j)
204             arrayin(nxIn+1,j) = arrayin(nxIn,j)
205             arrayin(nxIn+2,j) = arrayin(nxIn,j)
206            ENDDO
207          ENDIF
208          DO i=-1,nxIn+2
209             arrayin(i, 0)     = arrayin(i,1)
210             arrayin(i,-1)     = arrayin(i,1)
211             arrayin(i,nyIn+1) = arrayin(i,nyIn)
212             arrayin(i,nyIn+2) = arrayin(i,nyIn)
213          ENDDO
214    
215  C     For tracer (method=1,2) set to northernmost zonal-mean value  C-    For tracer (method=1,2) set to northernmost zonal-mean value
216  C     at 90N to avoid sharp zonal gradients near the Pole.  C     at 90N to avoid sharp zonal gradients near the Pole.
217  C     For U (method=11,12) set to zero at 90N to minimize velocity  C     For U (method=11,12) set to zero at 90N to minimize velocity
218  C     gradient at North Pole  C     gradient at North Pole
219  C     For V (method=11,12) set to northernmost zonal value at 90N,  C     For V (method=11,12) set to northernmost zonal value at 90N,
220  C     as is already done above in order to allow cross-PoleArctic flow  C     as is already done above in order to allow cross-PoleArctic flow
221        do j=ny_in,ny_in+2        DO j=nyIn,nyIn+2
222         if (y_in(j).eq.ninety) then         IF (y_in(j).EQ.ninety) THEN
223          if (method.eq.1 .or. method.eq.2) then          IF (method.EQ.1 .OR. method.EQ.2) THEN
224           NorthValue = 0.           NorthValue = 0.
225           do i=1,nx_in           DO i=1,nxIn
226            NorthValue = NorthValue + arrayin(i,j)            NorthValue = NorthValue + arrayin(i,j)
227           enddo           ENDDO
228           NorthValue = NorthValue / nx_in           NorthValue = NorthValue / nxIn
229           do i=-1,nx_in+2           DO i=-1,nxIn+2
230            arrayin(i,j) = NorthValue            arrayin(i,j) = NorthValue
231           enddo           ENDDO
232          elseif (method.eq.11 .or. method.eq.12) then          ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
233           do i=-1,nx_in+2           DO i=-1,nxIn+2
234            arrayin(i,j) = 0.            arrayin(i,j) = 0.
235           enddo           ENDDO
236          endif          ENDIF
237         endif         ENDIF
238        enddo        ENDDO
239    
240        do bj = mybylo(mythid), mybyhi(mythid)        DO bj = myByLo(myThid), myByHi(myThid)
241         do bi = mybxlo(mythid), mybxhi(mythid)         DO bi = myBxLo(myThid), myBxHi(myThid)
242    
243  C check validity of input/output coordinates  C--   Check validity of input/output coordinates
244  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
245          if ( debugLevel .ge. debLevB ) then          IF ( debugLevel.GE.debLevC ) THEN
246           do i=1,snx           DO j=1,sNy
247            do j=1,sny            DO i=1,sNx
248             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.             IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
249       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.       &          xG(i,j,bi,bj) .GE. x_in(nxIn+1)   .OR.
250       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.       &          yG(i,j,bi,bj) .LT. y_in(0)        .OR.
251       &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then       &          yG(i,j,bi,bj) .GE. y_in(nyIn+1) ) THEN
252                print*,'ERROR in S/R EXF_INTERP:'              l = ILNBLNK(inFile)
253                print*,'   input grid must encompass output grid.'              WRITE(msgBuf,'(3A,I6)')
254                print*,'i,j,bi,bj'      ,i,j,bi,bj       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
255                print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)              CALL PRINT_ERROR( msgBuf, myThid )
256                print*,'nx_in,ny_in'    ,nx_in        ,ny_in              WRITE(msgBuf,'(A)')
257                print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)       &        'EXF_INTERP: input grid must encompass output grid.'
258                print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)              CALL PRINT_ERROR( msgBuf, myThid )
259                STOP   ' ABNORMAL END: S/R EXF_INTERP'              WRITE(msgBuf,'(A,2I8,2I6,A,1P2E14.6)') 'i,j,bi,bj=',
260               endif       &        i,j,bi,bj, ' , xG,yG=', xG(i,j,bi,bj), yG(i,j,bi,bj)
261             enddo              CALL PRINT_ERROR( msgBuf, myThid )
262           enddo              WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nxIn=', nxIn,
263          endif       &        ' , x_in(0,nxIn+1)=', x_in(0) ,x_in(nxIn+1)
264                CALL PRINT_ERROR( msgBuf, myThid )
265                WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nyIn=', nyIn,
266         &        ' , y_in(0,nyIn+1)=', y_in(0) ,y_in(nyIn+1)
267                CALL PRINT_ERROR( msgBuf, myThid )
268                STOP 'ABNORMAL END: S/R EXF_INTERP'
269               ENDIF
270              ENDDO
271             ENDDO
272            ENDIF
273  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
274    
275  C compute interpolation indices  C--   Compute interpolation indices
276          do i=1,snx  #ifdef OLD_EXF_INTERP_LAT_INDEX
277           do j=1,sny          DO j=1,sNy
278            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then           DO i=1,sNx
279             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1            IF (xG(i,j,bi,bj)-x_in(1) .GE. 0.) THEN
280            else             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
281             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)            ELSE
282            endif             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc)
283            e_ind(i,j) = w_ind(i,j) + 1            ENDIF
284           enddo           ENDDO
285          enddo          ENDDO
286  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
287  C     use the original and more readable variant of the algorithm that  C-    use the original and more readable variant of the algorithm that
288  C     has unvectorizable while-loops for each (i,j)  C     has unvectorizable while-loops for each (i,j)
289          do i=1,snx          DO j=1,sNy
290           do j=1,sny           DO i=1,sNx
291            js = ny_in*.5            js = nyIn*.5
292            do while (yG(i,j,bi,bj) .lt. y_in(js))            DO WHILE (yG(i,j,bi,bj) .LT. y_in(js))
293             js = (js - 1)*.5             js = (js - 1)*.5
294            enddo            ENDDO
295            do while (yG(i,j,bi,bj) .ge. y_in(js+1))            DO WHILE (yG(i,j,bi,bj) .GE. y_in(js+1))
296             js = js + 1             js = js + 1
297            enddo            ENDDO
298            s_ind(i,j) = js            s_ind(i,j) = js
299           enddo           ENDDO
300          enddo          ENDDO
301  #else /* TARGET_NEC_SX defined */  #else /* TARGET_NEC_SX defined */
302  C     this variant vectorizes more efficiently than the original one because  C-    this variant vectorizes more efficiently than the original one because
303  C     it moves the while loops out of the i,j loops (loop pushing) but  C     it moves the while loops out of the i,j loops (loop pushing) but
304  C     it is ugly and incomprehensible  C     it is ugly and incomprehensible
305          icnt = 0          icnt = 0
306          do j=1,sny          DO j=1,sNy
307           do i=1,snx           DO i=1,sNx
308            s_ind(i,j) = ny_in*.5            s_ind(i,j) = nyIn*.5
309            icnt = icnt+1            icnt = icnt+1
310            inx(icnt,1) = i            inx(icnt,1) = i
311            inx(icnt,2) = j            inx(icnt,2) = j
312           enddo           ENDDO
313          enddo          ENDDO
314          do while (icnt .gt. 0)          DO WHILE (icnt .GT. 0)
315           ii = 0           ii = 0
316  !CDIR NODEP  !CDIR NODEP
317           do ic=1,icnt           DO ic=1,icnt
318            i = inx(ic,1)            i = inx(ic,1)
319            j = inx(ic,2)            j = inx(ic,2)
320            if (yG(i,j,bi,bj) .lt. y_in(s_ind(i,j))) then            IF (yG(i,j,bi,bj) .LT. y_in(s_ind(i,j))) THEN
321             s_ind(i,j) = (s_ind(i,j) - 1)*.5             s_ind(i,j) = (s_ind(i,j) - 1)*.5
322             ii = ii+1             ii = ii+1
323             inx(ii,1) = i             inx(ii,1) = i
324             inx(ii,2) = j             inx(ii,2) = j
325            endif            ENDIF
326           enddo           ENDDO
327           icnt = ii           icnt = ii
328          enddo          ENDDO
329          icnt = 0          icnt = 0
330          do j=1,sny          DO j=1,sNy
331           do i=1,snx           DO i=1,sNx
332            icnt = icnt+1            icnt = icnt+1
333            inx(icnt,1) = i            inx(icnt,1) = i
334            inx(icnt,2) = j            inx(icnt,2) = j
335           enddo           ENDDO
336          enddo          ENDDO
337          do while (icnt .gt. 0)          DO WHILE (icnt .GT. 0)
338           ii = 0           ii = 0
339  !CDIR NODEP  !CDIR NODEP
340           do ic=1,icnt           DO ic=1,icnt
341            i = inx(ic,1)            i = inx(ic,1)
342            j = inx(ic,2)            j = inx(ic,2)
343            if (yG(i,j,bi,bj) .ge. y_in(s_ind(i,j)+1)) then            IF (yG(i,j,bi,bj) .GE. y_in(s_ind(i,j)+1)) THEN
344             s_ind(i,j) = s_ind(i,j) + 1             s_ind(i,j) = s_ind(i,j) + 1
345             ii = ii+1             ii = ii+1
346             inx(ii,1) = i             inx(ii,1) = i
347             inx(ii,2) = j             inx(ii,2) = j
348            endif            ENDIF
349           enddo           ENDDO
350           icnt = ii           icnt = ii
351          enddo          ENDDO
352  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
353          do i=1,snx  #else /* OLD_EXF_INTERP_LAT_INDEX */
354           do j=1,sny  C--     latitude index
355            n_ind(i,j) = s_ind(i,j) + 1          DO j=1,sNy
356           enddo           DO i=1,sNx
357          enddo            s_ind(i,j) = 0
358              w_ind(i,j) = nyIn+1
359             ENDDO
360            ENDDO
361    C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
362    C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
363            tmpVar = nyIn + 1. _d -3
364            nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
365            DO l=1,nLoop
366             DO j=1,sNy
367              DO i=1,sNx
368               IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
369                k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
370                IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
371                  w_ind(i,j) = k
372                ELSE
373                  s_ind(i,j) = k
374                ENDIF
375               ENDIF
376              ENDDO
377             ENDDO
378            ENDDO
379    #ifdef ALLOW_DEBUG
380            IF ( debugLevel.GE.debLevC ) THEN
381    C-      Check that we found the right lat. index
382             DO j=1,sNy
383              DO i=1,sNx
384               IF ( w_ind(i,j).NE.s_ind(i,j)+1 ) THEN
385                l = ILNBLNK(inFile)
386                WRITE(msgBuf,'(3A,I4,A,I4)')
387         &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,
388         &        ', nLoop=', nLoop
389                CALL PRINT_ERROR( msgBuf, myThid )
390                WRITE(msgBuf,'(A)')
391         &        'EXF_INTERP: did not found Latitude index for grid-pt:'
392                CALL PRINT_ERROR( msgBuf, myThid )
393                WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
394         &        'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
395                CALL PRINT_ERROR( msgBuf, myThid )
396                WRITE(msgBuf,'(A,I8,A,1PE16.8)')
397         &        'EXF_INTERP: s_ind=',s_ind(i,j),', lat=',y_in(s_ind(i,j))
398                CALL PRINT_ERROR( msgBuf, myThid )
399                WRITE(msgBuf,'(A,I8,A,1PE16.8)')
400         &        'EXF_INTERP: n_ind=',w_ind(i,j),', lat=',y_in(w_ind(i,j))
401                CALL PRINT_ERROR( msgBuf, myThid )
402                STOP 'ABNORMAL END: S/R EXF_INTERP'
403               ENDIF
404              ENDDO
405             ENDDO
406            ENDIF
407    #endif /* ALLOW_DEBUG */
408    C--     longitude index
409            DO j=1,sNy
410             DO i=1,sNx
411               w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
412             ENDDO
413            ENDDO
414    #endif /* ndef OLD_EXF_INTERP_LAT_INDEX */
415    
416          if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then          IF (method.EQ.1 .OR. method.EQ.11 .OR. method.EQ.21) THEN
417    
418  C bilinear interpolation  C--   Bilinear interpolation
419           sp = 2           sp = 2
420           do j=1,sny           DO j=1,sNy
421            do i=1,snx            DO i=1,sNx
422             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
423             do l=0,1             DO l=0,1
424              px_ind(l+1) = x_in(w_ind(i,j)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
425              py_ind(l+1) = y_in(s_ind(i,j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
426             enddo             ENDDO
427  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
428             do k=1,2             DO k=1,2
429              ew_val(k) = arrayin(w_ind(i,j),s_ind(i,j)+k-1)              ew_val(k) = arrayin(w_ind(i,j)  ,s_ind(i,j)+k-1)
430       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
431       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)       &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-1)
432       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
433              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
434       &             +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)
435             enddo             ENDDO
436  #else  #else
437             ew_val1 = arrayin(w_ind(i,j),s_ind(i,j)+1-1)             ew_val1 = arrayin(w_ind(i,j)  ,s_ind(i,j)  )
438       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
439       &             +arrayin(e_ind(i,j),s_ind(i,j)+1-1)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
440       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
441             ew_val2 = arrayin(w_ind(i,j),s_ind(i,j)+2-1)             ew_val2 = arrayin(w_ind(i,j)  ,s_ind(i,j)+1)
442       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
443       &             +arrayin(e_ind(i,j),s_ind(i,j)+2-1)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
444       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
445             arrayout(i,j,bi,bj)=             arrayout(i,j,bi,bj)=
446       &             +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)       &            +ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
447       &             +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)       &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
448  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
449            enddo            ENDDO
450           enddo           ENDDO
451          elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then          ELSEIF (method .EQ. 2 .OR. method.EQ.12 .OR. method.EQ.22) THEN
452    
453  C bicubic interpolation  C--   Bicubic interpolation
454           sp = 4           sp = 4
455           do j=1,sny           DO j=1,sNy
456            do i=1,snx            DO i=1,sNx
457             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
458             do l=-1,2             DO l=-1,2
459              px_ind(l+2) = x_in(w_ind(i,j)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
460              py_ind(l+2) = y_in(s_ind(i,j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
461             enddo             ENDDO
462  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
463             do k=1,4             DO k=1,4
464              ew_val(k) =              ew_val(k) = arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
465       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
466       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
467       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
468       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-2)
469       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
470       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)+2,s_ind(i,j)+k-2)
471       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
472       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)              arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
473              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)       &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)
474       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)             ENDDO
            enddo  
475  #else  #else
476             ew_val1 =             ew_val1 = arrayin(w_ind(i,j)-1,s_ind(i,j)-1)
477       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+1-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
478       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)-1)
479       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+1-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
480       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)-1)
481       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+1-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
482       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)-1)
483       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+1-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
484       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)             ew_val2 = arrayin(w_ind(i,j)-1,s_ind(i,j)  )
485              ew_val2 =       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
486       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)  )
487       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
488       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
489       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
490       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)  )
491       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
492       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+2-2)             ew_val3 = arrayin(w_ind(i,j)-1,s_ind(i,j)+1)
493       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
494              ew_val3 =       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+1)
495       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+3-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
496       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
497       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+3-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
498       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+1)
499       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+3-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
500       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)             ew_val4 = arrayin(w_ind(i,j)-1,s_ind(i,j)+2)
501       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+3-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
502       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+2)
503              ew_val4 =       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
504       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+4-2)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+2)
505       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
506       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+4-2)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+2)
507       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
508       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+4-2)             arrayout(i,j,bi,bj) =
509       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
510       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+4-2)       &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
511       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &            +ew_val3*LAGRAN(3,yG(i,j,bi,bj),py_ind,sp)
512              arrayout(i,j,bi,bj)=       &            +ew_val4*LAGRAN(4,yG(i,j,bi,bj),py_ind,sp)
      &             +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)  
      &             +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)  
      &             +ew_val3*lagran(3,yG(i,j,bi,bj),py_ind,sp)  
      &             +ew_val4*lagran(4,yG(i,j,bi,bj),py_ind,sp)  
513  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
514            enddo            ENDDO
515           enddo           ENDDO
516          else          ELSE
517           stop 'stop in exf_interp.F: interpolation method not supported'           l = ILNBLNK(inFile)
518          endif           WRITE(msgBuf,'(3A,I6)')
519         enddo       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
520        enddo           CALL PRINT_ERROR( msgBuf, myThid )
521             WRITE(msgBuf,'(A,I8,A)')
522         &        'EXF_INTERP: method=', method,' not supported'
523             CALL PRINT_ERROR( msgBuf, myThid )
524             STOP 'ABNORMAL END: S/R EXF_INTERP: invalid method'
525            ENDIF
526           ENDDO
527          ENDDO
528    
529        RETURN        RETURN
530        END        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22