/[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.26 by jmc, Thu Dec 22 19:03:41 2011 UTC revision 1.27 by jmc, Sun Jan 1 01:24:54 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-|--+----|  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  c     _RL           lat_0, lat_inc(ny_in-1)  c     _RL           lat_0, lat_inc(nyIn-1)
101        _RL           lat_0, lat_inc(*)        _RL           lat_0, lat_inc(*)
102        integer       method, mythid        INTEGER       method, myThid
103    
104  C functions  C !FUNCTIONS:
105        external lagran        EXTERNAL LAGRAN
106        _RL      lagran        _RL      LAGRAN
107          INTEGER  ILNBLNK
108  C local variables        EXTERNAL ILNBLNK
109        integer  e_ind(snx,sny),w_ind(snx,sny)  
110        integer  n_ind(snx,sny),s_ind(snx,sny)  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  #ifdef OLD_EXF_INTERP_LAT_INDEX
122          INTEGER  js
123    #else
124          INTEGER  nLoop
125    #endif
126  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
127        integer  ic, ii, icnt        INTEGER  ic, ii, icnt
128        integer  inx(snx*sny,2)        INTEGER  inx(sNx*sNy,2)
129        _RL      ew_val1, ew_val2, ew_val3, ew_val4        _RL      ew_val1, ew_val2, ew_val3, ew_val4
130  #endif  #endif
131        _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 106  C local variables Line 133  C local variables
133        PARAMETER ( ninety = 90. )        PARAMETER ( ninety = 90. )
134        _RS      threeSixtyRS        _RS      threeSixtyRS
135        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
136          LOGICAL  xIsPeriodic
137    CEOP
138    
139  C     put xG in interval [ lon_0 , lon_0+360 [  C--   put xG in interval [ lon_0 , lon_0+360 [
140        do bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
141         do bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
142          do j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
143           do i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
144            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
145       &                  + threeSixtyRS*2.       &                  + threeSixtyRS*2.
146            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)
147           enddo           ENDDO
148          enddo          ENDDO
149         enddo         ENDDO
150        enddo        ENDDO
151    
152         call exf_interp_read(  C--   Load inut field
153       I   infile, filePrec,        CALL EXF_INTERP_READ(
154       O   arrayin,       I         inFile, filePrec,
155       I   irecord, nx_in, ny_in, mythid)       O         arrayin,
156         I         irecord, nxIn, nyIn, myThid )
157    
158  C setup input longitude grid  C--   setup input longitude grid
159        do i=-1,nx_in+2        DO i=-1,nxIn+2
160         x_in(i) = lon_0 + (i-1)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
161        enddo        ENDDO
162          xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
163    
164  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)  
165        y_in(1) = lat_0        y_in(1) = lat_0
166        do j=2,ny_in        DO j=1,nyIn+1
167         y_in(j) = y_in(j-1) + lat_inc(j-1)         i = MIN(j,nyIn-1)
168        enddo         y_in(j+1) = y_in(j) + lat_inc(i)
169        do j=ny_in+1,ny_in+2        ENDDO
170         if (y_in(j-1).eq.ninety) then  C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
171          y_in(j) = 2 * ninety - y_in(j-2)        y_in(0) = y_in(1) - lat_inc(1)
172         else        y_in(-1)= y_in(0) - lat_inc(1)
173          i = max(1,ny_in-1)  c     IF ( y_in(1).GT.-ninety .AND. y_in(0).LT.-ninety ) THEN
174          y_in(j) = min( y_in(j-1)+lat_inc(i), ninety )  c       y_in(0) = -ninety
175         endif  c       y_in(-1) = -2.*ninety - y_in(1)
176        enddo  c     ENDIF
177    c     IF ( y_in(0).GT.-ninety .AND. y_in(-1).LT.-ninety ) THEN
178  C enlarge boundary  c       y_in(-1) = -ninety
179        do j=1,ny_in  c     ENDIF
180         arrayin(0,j)       = arrayin(nx_in,j)  C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
181         arrayin(-1,j)      = arrayin(nx_in-1,j)        j = nyIn+1
182         arrayin(nx_in+1,j) = arrayin(1,j)        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN
183         arrayin(nx_in+2,j) = arrayin(2,j)          y_in(j) = ninety
184        enddo          y_in(j+1) = 2.*ninety - y_in(j-1)
185        do i=-1,nx_in+2        ENDIF
186         arrayin(i,0)       = arrayin(i,1)        j = nyIn+2
187         arrayin(i,-1)      = arrayin(i,1)        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN
188         arrayin(i,ny_in+1) = arrayin(i,ny_in)          y_in(j) = ninety
189         arrayin(i,ny_in+2) = arrayin(i,ny_in)        ENDIF
190        enddo  
191    C--   enlarge boundary
192          IF ( xIsPeriodic ) THEN
193            DO j=1,nyIn
194             arrayin( 0,j)     = arrayin(nxIn  ,j)
195             arrayin(-1,j)     = arrayin(nxIn-1,j)
196             arrayin(nxIn+1,j) = arrayin(1,j)
197             arrayin(nxIn+2,j) = arrayin(2,j)
198            ENDDO
199          ELSE
200            DO j=1,nyIn
201             arrayin( 0,j)     = arrayin(1,j)
202             arrayin(-1,j)     = arrayin(1,j)
203             arrayin(nxIn+1,j) = arrayin(nxIn,j)
204             arrayin(nxIn+2,j) = arrayin(nxIn,j)
205            ENDDO
206          ENDIF
207          DO i=-1,nxIn+2
208             arrayin(i, 0)     = arrayin(i,1)
209             arrayin(i,-1)     = arrayin(i,1)
210             arrayin(i,nyIn+1) = arrayin(i,nyIn)
211             arrayin(i,nyIn+2) = arrayin(i,nyIn)
212          ENDDO
213    
214  C     For tracer (method=1,2) set to northernmost zonal-mean value  C-    For tracer (method=1,2) set to northernmost zonal-mean value
215  C     at 90N to avoid sharp zonal gradients near the Pole.  C     at 90N to avoid sharp zonal gradients near the Pole.
216  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
217  C     gradient at North Pole  C     gradient at North Pole
218  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,
219  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
220        do j=ny_in,ny_in+2        DO j=nyIn,nyIn+2
221         if (y_in(j).eq.ninety) then         IF (y_in(j).EQ.ninety) THEN
222          if (method.eq.1 .or. method.eq.2) then          IF (method.EQ.1 .OR. method.EQ.2) THEN
223           NorthValue = 0.           NorthValue = 0.
224           do i=1,nx_in           DO i=1,nxIn
225            NorthValue = NorthValue + arrayin(i,j)            NorthValue = NorthValue + arrayin(i,j)
226           enddo           ENDDO
227           NorthValue = NorthValue / nx_in           NorthValue = NorthValue / nxIn
228           do i=-1,nx_in+2           DO i=-1,nxIn+2
229            arrayin(i,j) = NorthValue            arrayin(i,j) = NorthValue
230           enddo           ENDDO
231          elseif (method.eq.11 .or. method.eq.12) then          ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
232           do i=-1,nx_in+2           DO i=-1,nxIn+2
233            arrayin(i,j) = 0.            arrayin(i,j) = 0.
234           enddo           ENDDO
235          endif          ENDIF
236         endif         ENDIF
237        enddo        ENDDO
238    
239        do bj = mybylo(mythid), mybyhi(mythid)        DO bj = myByLo(myThid), myByHi(myThid)
240         do bi = mybxlo(mythid), mybxhi(mythid)         DO bi = myBxLo(myThid), myBxHi(myThid)
241    
242  C check validity of input/output coordinates  C--   Check validity of input/output coordinates
243  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
244          if ( debugLevel.GE.debLevC ) then          IF ( debugLevel.GE.debLevC ) THEN
245           do j=1,sny           DO j=1,sNy
246            do i=1,snx            DO i=1,sNx
247             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.             IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
248       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.       &          xG(i,j,bi,bj) .GE. x_in(nxIn+1)   .OR.
249       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.       &          yG(i,j,bi,bj) .LT. y_in(0)        .OR.
250       &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then       &          yG(i,j,bi,bj) .GE. y_in(nyIn+1) ) THEN
251                print*,'ERROR in S/R EXF_INTERP:'              l = ILNBLNK(inFile)
252                print*,'   input grid must encompass output grid.'              WRITE(msgBuf,'(3A,I6)')
253                print*,'i,j,bi,bj'      ,i,j,bi,bj       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
254                print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)              CALL PRINT_ERROR( msgBuf, myThid )
255                print*,'nx_in,ny_in'    ,nx_in        ,ny_in              WRITE(msgBuf,'(A)')
256                print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)       &        'EXF_INTERP: input grid must encompass output grid.'
257                print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)              CALL PRINT_ERROR( msgBuf, myThid )
258                STOP   ' ABNORMAL END: S/R EXF_INTERP'              WRITE(msgBuf,'(A,2I8,2I6,A,1P2E14.6)') 'i,j,bi,bj=',
259             endif       &        i,j,bi,bj, ' , xG,yG=', xG(i,j,bi,bj), yG(i,j,bi,bj)
260            enddo              CALL PRINT_ERROR( msgBuf, myThid )
261           enddo              WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nxIn=', nxIn,
262          endif       &        ' , x_in(0,nxIn+1)=', x_in(0) ,x_in(nxIn+1)
263                CALL PRINT_ERROR( msgBuf, myThid )
264                WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nyIn=', nyIn,
265         &        ' , y_in(0,nyIn+1)=', y_in(0) ,y_in(nyIn+1)
266                CALL PRINT_ERROR( msgBuf, myThid )
267                STOP 'ABNORMAL END: S/R EXF_INTERP'
268               ENDIF
269              ENDDO
270             ENDDO
271            ENDIF
272  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
273    
274  C compute interpolation indices  C--   Compute interpolation indices
275          do j=1,sny  #ifdef OLD_EXF_INTERP_LAT_INDEX
276           do i=1,snx          DO j=1,sNy
277            if (xG(i,j,bi,bj)-x_in(1) .ge. 0.) then           DO i=1,sNx
278             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
279            else             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1
280             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)            ELSE
281            endif             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc)
282            e_ind(i,j) = w_ind(i,j) + 1            ENDIF
283           enddo           ENDDO
284          enddo          ENDDO
285  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
286  C     use the original and more readable variant of the algorithm that  C-    use the original and more readable variant of the algorithm that
287  C     has unvectorizable while-loops for each (i,j)  C     has unvectorizable while-loops for each (i,j)
288          do j=1,sny          DO j=1,sNy
289           do i=1,snx           DO i=1,sNx
290            js = ny_in*.5            js = nyIn*.5
291            do while (yG(i,j,bi,bj) .lt. y_in(js))            DO WHILE (yG(i,j,bi,bj) .LT. y_in(js))
292             js = (js - 1)*.5             js = (js - 1)*.5
293            enddo            ENDDO
294            do while (yG(i,j,bi,bj) .ge. y_in(js+1))            DO WHILE (yG(i,j,bi,bj) .GE. y_in(js+1))
295             js = js + 1             js = js + 1
296            enddo            ENDDO
297            s_ind(i,j) = js            s_ind(i,j) = js
298           enddo           ENDDO
299          enddo          ENDDO
300  #else /* TARGET_NEC_SX defined */  #else /* TARGET_NEC_SX defined */
301  C     this variant vectorizes more efficiently than the original one because  C-    this variant vectorizes more efficiently than the original one because
302  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
303  C     it is ugly and incomprehensible  C     it is ugly and incomprehensible
304          icnt = 0          icnt = 0
305          do j=1,sny          DO j=1,sNy
306           do i=1,snx           DO i=1,sNx
307            s_ind(i,j) = ny_in*.5            s_ind(i,j) = nyIn*.5
308            icnt = icnt+1            icnt = icnt+1
309            inx(icnt,1) = i            inx(icnt,1) = i
310            inx(icnt,2) = j            inx(icnt,2) = j
311           enddo           ENDDO
312          enddo          ENDDO
313          do while (icnt .gt. 0)          DO WHILE (icnt .GT. 0)
314           ii = 0           ii = 0
315  !CDIR NODEP  !CDIR NODEP
316           do ic=1,icnt           DO ic=1,icnt
317            i = inx(ic,1)            i = inx(ic,1)
318            j = inx(ic,2)            j = inx(ic,2)
319            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
320             s_ind(i,j) = (s_ind(i,j) - 1)*.5             s_ind(i,j) = (s_ind(i,j) - 1)*.5
321             ii = ii+1             ii = ii+1
322             inx(ii,1) = i             inx(ii,1) = i
323             inx(ii,2) = j             inx(ii,2) = j
324            endif            ENDIF
325           enddo           ENDDO
326           icnt = ii           icnt = ii
327          enddo          ENDDO
328          icnt = 0          icnt = 0
329          do j=1,sny          DO j=1,sNy
330           do i=1,snx           DO i=1,sNx
331            icnt = icnt+1            icnt = icnt+1
332            inx(icnt,1) = i            inx(icnt,1) = i
333            inx(icnt,2) = j            inx(icnt,2) = j
334           enddo           ENDDO
335          enddo          ENDDO
336          do while (icnt .gt. 0)          DO WHILE (icnt .GT. 0)
337           ii = 0           ii = 0
338  !CDIR NODEP  !CDIR NODEP
339           do ic=1,icnt           DO ic=1,icnt
340            i = inx(ic,1)            i = inx(ic,1)
341            j = inx(ic,2)            j = inx(ic,2)
342            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
343             s_ind(i,j) = s_ind(i,j) + 1             s_ind(i,j) = s_ind(i,j) + 1
344             ii = ii+1             ii = ii+1
345             inx(ii,1) = i             inx(ii,1) = i
346             inx(ii,2) = j             inx(ii,2) = j
347            endif            ENDIF
348           enddo           ENDDO
349           icnt = ii           icnt = ii
350          enddo          ENDDO
351  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
352          do j=1,sny  #else /* OLD_EXF_INTERP_LAT_INDEX */
353           do i=1,snx  C--     latitude index
354            n_ind(i,j) = s_ind(i,j) + 1          DO j=1,sNy
355           enddo           DO i=1,sNx
356          enddo            s_ind(i,j) = 0
357              w_ind(i,j) = nyIn+1
358             ENDDO
359            ENDDO
360    C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
361    C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
362            nLoop = 1 + INT( LOG(DFLOAT(nyIn)+1. _d -3)/LOG(2. _d 0) )
363            DO l=1,nLoop
364             DO j=1,sNy
365              DO i=1,sNx
366               IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
367                k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
368                IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
369                  w_ind(i,j) = k
370                ELSE
371                  s_ind(i,j) = k
372                ENDIF
373               ENDIF
374              ENDDO
375             ENDDO
376            ENDDO
377    #ifdef ALLOW_DEBUG
378            IF ( debugLevel.GE.debLevC ) THEN
379    C-      Check that we found the right lat. index
380             DO j=1,sNy
381              DO i=1,sNx
382               IF ( w_ind(i,j).NE.s_ind(i,j)+1 ) THEN
383                l = ILNBLNK(inFile)
384                WRITE(msgBuf,'(3A,I4,A,I4)')
385         &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,
386         &        ', nLoop=', nLoop
387                CALL PRINT_ERROR( msgBuf, myThid )
388                WRITE(msgBuf,'(A)')
389         &        'EXF_INTERP: did not found Latitude index for grid-pt:'
390                CALL PRINT_ERROR( msgBuf, myThid )
391                WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
392         &        'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
393                CALL PRINT_ERROR( msgBuf, myThid )
394                WRITE(msgBuf,'(A,I8,A,1PE16.8)')
395         &        'EXF_INTERP: s_ind=',s_ind(i,j),', lat=',y_in(s_ind(i,j))
396                CALL PRINT_ERROR( msgBuf, myThid )
397                WRITE(msgBuf,'(A,I8,A,1PE16.8)')
398         &        'EXF_INTERP: n_ind=',w_ind(i,j),', lat=',y_in(w_ind(i,j))
399                CALL PRINT_ERROR( msgBuf, myThid )
400                STOP 'ABNORMAL END: S/R EXF_INTERP'
401               ENDIF
402              ENDDO
403             ENDDO
404            ENDIF
405    #endif /* ALLOW_DEBUG */
406    C--     longitude index
407            DO j=1,sNy
408             DO i=1,sNx
409               w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
410             ENDDO
411            ENDDO
412    #endif /* ndef OLD_EXF_INTERP_LAT_INDEX */
413    
414          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
415    
416  C bilinear interpolation  C--   Bilinear interpolation
417           sp = 2           sp = 2
418           do j=1,sny           DO j=1,sNy
419            do i=1,snx            DO i=1,sNx
420             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
421             do l=0,1             DO l=0,1
422              px_ind(l+1) = x_in(w_ind(i,j)+l)              px_ind(l+1) = x_in(w_ind(i,j)+l)
423              py_ind(l+1) = y_in(s_ind(i,j)+l)              py_ind(l+1) = y_in(s_ind(i,j)+l)
424             enddo             ENDDO
425  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
426             do k=1,2             DO k=1,2
427              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)
428       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
429       &             +arrayin(e_ind(i,j),s_ind(i,j)+k-1)       &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-1)
430       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
431              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)              arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
432       &             +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)
433             enddo             ENDDO
434  #else  #else
435             ew_val1 = arrayin(w_ind(i,j),s_ind(i,j)+1-1)             ew_val1 = arrayin(w_ind(i,j)  ,s_ind(i,j)  )
436       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
437       &             +arrayin(e_ind(i,j),s_ind(i,j)+1-1)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
438       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
439             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)
440       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
441       &             +arrayin(e_ind(i,j),s_ind(i,j)+2-1)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
442       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
443             arrayout(i,j,bi,bj)=             arrayout(i,j,bi,bj)=
444       &             +ew_val1*lagran(1,yG(i,j,bi,bj),py_ind,sp)       &            +ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
445       &             +ew_val2*lagran(2,yG(i,j,bi,bj),py_ind,sp)       &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
446  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
447            enddo            ENDDO
448           enddo           ENDDO
449          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
450    
451  C bicubic interpolation  C--   Bicubic interpolation
452           sp = 4           sp = 4
453           do j=1,sny           DO j=1,sNy
454            do i=1,snx            DO i=1,sNx
455             arrayout(i,j,bi,bj) = 0.             arrayout(i,j,bi,bj) = 0.
456             do l=-1,2             DO l=-1,2
457              px_ind(l+2) = x_in(w_ind(i,j)+l)              px_ind(l+2) = x_in(w_ind(i,j)+l)
458              py_ind(l+2) = y_in(s_ind(i,j)+l)              py_ind(l+2) = y_in(s_ind(i,j)+l)
459             enddo             ENDDO
460  #ifndef TARGET_NEC_SX  #ifndef TARGET_NEC_SX
461             do k=1,4             DO k=1,4
462              ew_val(k) =              ew_val(k) = arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
463       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
464       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
465       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
466       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-2)
467       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+k-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
468       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &                + arrayin(w_ind(i,j)+2,s_ind(i,j)+k-2)
469       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+k-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
470       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)              arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
471              arrayout(i,j,bi,bj)=arrayout(i,j,bi,bj)       &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)
472       &             +ew_val(k)*lagran(k,yG(i,j,bi,bj),py_ind,sp)             ENDDO
            enddo  
473  #else  #else
474             ew_val1 =             ew_val1 = arrayin(w_ind(i,j)-1,s_ind(i,j)-1)
475       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+1-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
476       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)-1)
477       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+1-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
478       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)-1)
479       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+1-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
480       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)-1)
481       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+1-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
482       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)             ew_val2 = arrayin(w_ind(i,j)-1,s_ind(i,j)  )
483              ew_val2 =       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
484       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)  )
485       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
486       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
487       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
488       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+2-2)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)  )
489       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
490       &             +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)
491       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
492              ew_val3 =       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+1)
493       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+3-2)       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
494       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
495       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+3-2)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
496       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+1)
497       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+3-2)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
498       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)             ew_val4 = arrayin(w_ind(i,j)-1,s_ind(i,j)+2)
499       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+3-2)       &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
500       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+2)
501              ew_val4 =       &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
502       &             arrayin(w_ind(i,j)-1,s_ind(i,j)+4-2)       &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+2)
503       &             *lagran(1,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
504       &             +arrayin(w_ind(i,j)  ,s_ind(i,j)+4-2)       &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+2)
505       &             *lagran(2,xG(i,j,bi,bj),px_ind,sp)       &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
506       &             +arrayin(e_ind(i,j)  ,s_ind(i,j)+4-2)             arrayout(i,j,bi,bj) =
507       &             *lagran(3,xG(i,j,bi,bj),px_ind,sp)       &             ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
508       &             +arrayin(e_ind(i,j)+1,s_ind(i,j)+4-2)       &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
509       &             *lagran(4,xG(i,j,bi,bj),px_ind,sp)       &            +ew_val3*LAGRAN(3,yG(i,j,bi,bj),py_ind,sp)
510              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)  
511  #endif /* TARGET_NEC_SX defined */  #endif /* TARGET_NEC_SX defined */
512            enddo            ENDDO
513           enddo           ENDDO
514          else          ELSE
515           stop 'stop in exf_interp.F: interpolation method not supported'           l = ILNBLNK(inFile)
516          endif           WRITE(msgBuf,'(3A,I6)')
517         enddo       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
518        enddo           CALL PRINT_ERROR( msgBuf, myThid )
519             WRITE(msgBuf,'(A,I8,A)')
520         &        'EXF_INTERP: method=', method,' not supported'
521             CALL PRINT_ERROR( msgBuf, myThid )
522             STOP 'ABNORMAL END: S/R EXF_INTERP: invalid method'
523            ENDIF
524           ENDDO
525          ENDDO
526    
527        RETURN        RETURN
528        END        END
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22