/[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.27 by jmc, Sun Jan 1 01:24:54 2012 UTC revision 1.37 by jmc, Fri Mar 10 00:16:11 2017 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
 #undef OLD_EXF_INTERP_LAT_INDEX  
   
 C==========================================*  
 C Flux Coupler using                       |  
 C Bilinear interpolation of forcing fields |  
 C                                          |  
 C B. Cheng (12/2002)                       |  
 C                                          |  
 C added Bicubic (bnc 1/2003)               |  
 C                                          |  
 C==========================================*  
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
        _RL FUNCTION LAGRAN(i,x,a,sp)  
   
        IMPLICIT NONE  
   
         INTEGER i  
         _RS x  
         _RL a(4)  
         INTEGER sp  
   
 C-      local variables:  
         INTEGER k  
         _RL numer,denom  
   
         numer = 1. _d 0  
         denom = 1. _d 0  
   
 #ifdef TARGET_NEC_SX  
 !CDIR UNROLL=8  
 #endif /* TARGET_NEC_SX */  
         DO k=1,sp  
          IF ( k .NE. i) THEN  
           denom = denom*(a(i) - a(k))  
           numer = numer*(x    - a(k))  
          ENDIF  
         ENDDO  
   
         LAGRAN = numer/denom  
   
        RETURN  
        END  
5    
6  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
# Line 53  CBOP Line 9  CBOP
9  C     !ROUTINE: EXF_INTERP  C     !ROUTINE: EXF_INTERP
10  C     !INTERFACE:  C     !INTERFACE:
11         SUBROUTINE EXF_INTERP(         SUBROUTINE EXF_INTERP(
12       I                inFile,       I                inFile, filePrec,
13       I                filePrec,  #ifdef EXF_INTERP_USE_DYNALLOC
14       O                arrayout,       O                arrayout,
15    #else
16         O                arrayout, arrayin,
17    #endif
18       I                irecord, xG_in, yG,       I                irecord, xG_in, yG,
19       I                lon_0, lon_inc,       I                lon_0, lon_inc, lat_0, lat_inc,
20       I                lat_0, lat_inc,       I                nxIn, nyIn, method, myIter, myThid )
      I                nxIn, nyIn, method, myThid )  
21    
22  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
23  C  *==========================================================*  C  *==========================================================*
# Line 75  C     === Global variables === Line 33  C     === Global variables ===
33  #include "SIZE.h"  #include "SIZE.h"
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36    #include "EXF_INTERP_SIZE.h"
37    #ifdef ALLOW_DEBUG
38    # include "EXF_PARAM.h"
39    #endif
40    
41  C !INPUT/OUTPUT PARAMETERS:  C !INPUT/OUTPUT PARAMETERS:
42  C  inFile      (string)  :: name of the binary input file (direct access)  C  inFile      (string)  :: name of the binary input file (direct access)
43  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)
44  C  arrayout    ( _RL )   :: output array  C  arrayout    ( _RL )   :: output array
45    #ifndef EXF_INTERP_USE_DYNALLOC
46    C  arrayin     ( _RL )   :: input field array (loaded from file)
47    #endif
48  C  irecord     (integer) :: record number to read  C  irecord     (integer) :: record number to read
49  C     xG_in,yG           :: coordinates for output grid to interpolate to  C     xG_in,yG           :: coordinates for output grid to interpolate to
50  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
# Line 89  C     lat_inc            :: vector y-gri Line 53  C     lat_inc            :: vector y-gri
53  C  nxIn,nyIn   (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
54  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
55  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
56    C  myIter      (integer) :: current iteration number
57  C  myThid      (integer) :: My Thread Id number  C  myThid      (integer) :: My Thread Id number
58          CHARACTER*(*) inFile
       CHARACTER*(*) infile  
59        INTEGER       filePrec, irecord, nxIn, nyIn        INTEGER       filePrec, irecord, nxIn, nyIn
60        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61    #ifndef EXF_INTERP_USE_DYNALLOC
62          _RL           arrayin ( -1:nxIn+2, -1:nyIn+2 )
63    #endif
64        _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)
65        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
67  c     _RL           lat_0, lat_inc(nyIn-1)  c     _RL           lat_0, lat_inc(nyIn-1)
68        _RL           lat_0, lat_inc(*)        _RL           lat_0, lat_inc(*)
69        INTEGER       method, myThid        INTEGER       method, myIter, myThid
70    
71  C !FUNCTIONS:  C !FUNCTIONS:
72        EXTERNAL LAGRAN  #ifdef ALLOW_DEBUG
       _RL      LAGRAN  
73        INTEGER  ILNBLNK        INTEGER  ILNBLNK
74        EXTERNAL ILNBLNK        EXTERNAL ILNBLNK
75    #endif
76    
77  C !LOCAL VARIABLES:  C !LOCAL VARIABLES:
78  C     msgBuf      :: Informational/error message buffer  C     x_in        :: longitude vector defining input field grid
79    C     y_in        :: latitude  vector defining input field grid
80    C     w_ind       :: input field longitudinal index, on western side of model grid pt
81    C     s_ind       :: input field latitudinal index, on southern side of model grid pt
82  C     bi, bj      :: tile indices  C     bi, bj      :: tile indices
83        CHARACTER*(MAX_LEN_MBUF) msgBuf  C     i, j, k, l  :: loop indices
84        INTEGER  bi, bj  C     msgBuf      :: Informational/error message buffer
85        INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)  #ifdef EXF_INTERP_USE_DYNALLOC
86        _RL      px_ind(4), py_ind(4), ew_val(4)  C     arrayin     :: input field array (loaded from file)
87        _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )        _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )
88        _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)        _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)
89        _RL      NorthValue  #else /* EXF_INTERP_USE_DYNALLOC */
90        INTEGER  i, j, k, l, sp        _RL      x_in(-1:exf_max_nLon+2), y_in(-1:exf_max_nLat+2)
91  #ifdef OLD_EXF_INTERP_LAT_INDEX  #endif /* EXF_INTERP_USE_DYNALLOC */
92        INTEGER  js        INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)
93  #else        INTEGER  bi, bj
94          INTEGER  i, j, k, l
95        INTEGER  nLoop        INTEGER  nLoop
96  #endif        _RL      tmpVar
 #ifdef TARGET_NEC_SX  
       INTEGER  ic, ii, icnt  
       INTEGER  inx(sNx*sNy,2)  
       _RL      ew_val1, ew_val2, ew_val3, ew_val4  
 #endif  
97        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL      ninety  
       PARAMETER ( ninety = 90. )  
98        _RS      threeSixtyRS        _RS      threeSixtyRS
99          _RL      yPole, symSign, poleValue
100        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
101        LOGICAL  xIsPeriodic        PARAMETER ( yPole = 90. )
102          INTEGER  nxd2
103          LOGICAL  xIsPeriodic, poleSymmetry
104    #ifdef ALLOW_DEBUG
105          LOGICAL  debugFlag
106          CHARACTER*(MAX_LEN_MBUF) msgBuf
107          _RS      prtPole(-1:4)
108    #endif
109  CEOP  CEOP
110    
111  C--   put xG in interval [ lon_0 , lon_0+360 [  #ifndef EXF_INTERP_USE_DYNALLOC
112        DO bj=myByLo(myThid),myByHi(myThid)  C--   Check size declaration:
113         DO bi=myBxLo(myThid),myBxHi(myThid)        IF ( nxIn.GT.exf_max_nLon ) THEN
114          DO j=1-OLy,sNy+OLy         STOP 'EXF_INTERP: exf_max_nLon too small'
115           DO i=1-OLx,sNx+OLx        ENDIF
116            xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0        IF ( nyIn.GT.exf_max_nLat ) THEN
117       &                  + threeSixtyRS*2.         STOP 'EXF_INTERP: exf_max_nLat too small'
118            xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)        ENDIF
119           ENDDO        IF ( (nxIn+4)*(nyIn+4).GT.exf_interp_bufferSize ) THEN
120          ENDDO         STOP 'EXF_INTERP: exf_interp_bufferSize too small'
121         ENDDO        ENDIF
122        ENDDO  #endif /* ndef EXF_INTERP_USE_DYNALLOC */
123    
124    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
125    C---  Load inut field
126    
 C--   Load inut field  
127        CALL EXF_INTERP_READ(        CALL EXF_INTERP_READ(
128       I         inFile, filePrec,       I         inFile, filePrec,
129       O         arrayin,       O         arrayin,
130       I         irecord, nxIn, nyIn, myThid )       I         irecord, nxIn, nyIn, myThid )
131    
132    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
133    C---  Prepare input grid and input field
134    
135  C--   setup input longitude grid  C--   setup input longitude grid
136        DO i=-1,nxIn+2        DO i=-1,nxIn+2
137         x_in(i) = lon_0 + (i-1)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
138        ENDDO        ENDDO
139        xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )        xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
140          nxd2 = NINT( nxIn*0.5 )
141          poleSymmetry = xIsPeriodic .AND. ( nxIn.EQ.2*nxd2 )
142    #ifdef EXF_USE_OLD_INTERP_POLE
143          poleSymmetry = .FALSE.
144    #endif
145    
146  C--   setup input latitude grid  C--   setup input latitude grid
147        y_in(1) = lat_0        y_in(1) = lat_0
# Line 167  C--   setup input latitude grid Line 149  C--   setup input latitude grid
149         i = MIN(j,nyIn-1)         i = MIN(j,nyIn-1)
150         y_in(j+1) = y_in(j) + lat_inc(i)         y_in(j+1) = y_in(j) + lat_inc(i)
151        ENDDO        ENDDO
 C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole  
152        y_in(0) = y_in(1) - lat_inc(1)        y_in(0) = y_in(1) - lat_inc(1)
153        y_in(-1)= y_in(0) - lat_inc(1)        y_in(-1)= y_in(0) - lat_inc(1)
154  c     IF ( y_in(1).GT.-ninety .AND. y_in(0).LT.-ninety ) THEN  #ifdef ALLOW_DEBUG
155  c       y_in(0) = -ninety        DO l=-1,4
156  c       y_in(-1) = -2.*ninety - y_in(1)          prtPole(l) = 0.
157  c     ENDIF        ENDDO
158  c     IF ( y_in(0).GT.-ninety .AND. y_in(-1).LT.-ninety ) THEN  #endif
159  c       y_in(-1) = -ninety  C--   For tracer (method=1,2) add 1 row @ the pole (if last row is not)
160  c     ENDIF  C     and will fill it with the polarmost-row zonally averaged value.
161    C     For vector component, cannot do much without the other component;
162    C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
163    #ifdef EXF_USE_OLD_INTERP_POLE
164          IF ( .TRUE. ) THEN
165    #else
166          IF ( method.LT.10 ) THEN
167    C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
168           j = 0
169           IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
170             y_in(j) = -yPole
171             y_in(j-1) = -2.*yPole - y_in(j+1)
172    #ifdef ALLOW_DEBUG
173             prtPole(j)   = 1.
174             prtPole(j-1) = 2.
175    #endif
176           ENDIF
177           j = -1
178           IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
179             y_in(j) = -yPole
180    #ifdef ALLOW_DEBUG
181             prtPole(j)   = 1.
182    #endif
183           ENDIF
184    #endif /* EXF_USE_OLD_INTERP_POLE */
185  C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole  C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
186        j = nyIn+1         j = nyIn+1
187        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN         IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
188          y_in(j) = ninety           y_in(j) = yPole
189          y_in(j+1) = 2.*ninety - y_in(j-1)           y_in(j+1) = 2.*yPole - y_in(j-1)
190        ENDIF  #ifdef ALLOW_DEBUG
191        j = nyIn+2           prtPole(3)   = 1.
192        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN           prtPole(3+1) = 2.
193          y_in(j) = ninety  #endif
194           ENDIF
195           j = nyIn+2
196           IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
197             y_in(j) = yPole
198    #ifdef ALLOW_DEBUG
199             prtPole(4)   = 1.
200    #endif
201           ENDIF
202        ENDIF        ENDIF
203    
204  C--   enlarge boundary  C--   Enlarge boundary
205        IF ( xIsPeriodic ) THEN        IF ( xIsPeriodic ) THEN
206    C-    fill-in added column assuming periodicity
207          DO j=1,nyIn          DO j=1,nyIn
208           arrayin( 0,j)     = arrayin(nxIn  ,j)           arrayin( 0,j)     = arrayin(nxIn  ,j)
209           arrayin(-1,j)     = arrayin(nxIn-1,j)           arrayin(-1,j)     = arrayin(nxIn-1,j)
# Line 197  C--   enlarge boundary Line 211  C--   enlarge boundary
211           arrayin(nxIn+2,j) = arrayin(2,j)           arrayin(nxIn+2,j) = arrayin(2,j)
212          ENDDO          ENDDO
213        ELSE        ELSE
214    C-    fill-in added column from nearest column
215          DO j=1,nyIn          DO j=1,nyIn
216           arrayin( 0,j)     = arrayin(1,j)           arrayin( 0,j)     = arrayin(1,j)
217           arrayin(-1,j)     = arrayin(1,j)           arrayin(-1,j)     = arrayin(1,j)
# Line 204  C--   enlarge boundary Line 219  C--   enlarge boundary
219           arrayin(nxIn+2,j) = arrayin(nxIn,j)           arrayin(nxIn+2,j) = arrayin(nxIn,j)
220          ENDDO          ENDDO
221        ENDIF        ENDIF
222        DO i=-1,nxIn+2        symSign = 1. _d 0
223           arrayin(i, 0)     = arrayin(i,1)        IF ( method.GE.10 ) symSign = -1. _d 0
224           arrayin(i,-1)     = arrayin(i,1)        DO l=-1,2
225           arrayin(i,nyIn+1) = arrayin(i,nyIn)         j = l
226           arrayin(i,nyIn+2) = arrayin(i,nyIn)         IF ( l.GE.1 ) j = nyIn+l
227           k = MAX(1,MIN(j,nyIn))
228           IF ( poleSymmetry .AND. ABS(y_in(j)).GT.yPole ) THEN
229    C-    fill-in added row assuming pole-symmetry
230            DO i=-1,nxd2
231             arrayin(i,j) = symSign*arrayin(i+nxd2,k)
232            ENDDO
233            DO i=1,nxd2+2
234             arrayin(i+nxd2,j) = symSign*arrayin(i,k)
235            ENDDO
236    #ifdef ALLOW_DEBUG
237            i = l + 2*( (l+1)/2 )
238            prtPole(i) = prtPole(i) + 0.2
239    #endif
240           ELSE
241    C-    fill-in added row from nearest column values
242            DO i=-1,nxIn+2
243             arrayin(i,j) = arrayin(i,k)
244            ENDDO
245           ENDIF
246        ENDDO        ENDDO
247    
248  C-    For tracer (method=1,2) set to northernmost zonal-mean value  C--   For tracer (method=1,2) set to northernmost zonal-mean value
249  C     at 90N to avoid sharp zonal gradients near the Pole.  C     at 90N to avoid sharp zonal gradients near the Pole.
250  C     For U (method=11,12) set to zero at 90N to minimize velocity  C     For vector component, cannot do much without the other component;
251  C     gradient at North Pole  C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
252  C     For V (method=11,12) set to northernmost zonal value at 90N,  #ifdef EXF_USE_OLD_INTERP_POLE
253  C     as is already done above in order to allow cross-PoleArctic flow        IF ( .TRUE. ) THEN
254        DO j=nyIn,nyIn+2         DO l= 3,4
255         IF (y_in(j).EQ.ninety) THEN  #else
256          IF (method.EQ.1 .OR. method.EQ.2) THEN        IF ( method.LT.10 ) THEN
257           NorthValue = 0.         DO l=-1,4
258           DO i=1,nxIn  #endif
259            NorthValue = NorthValue + arrayin(i,j)          j = l
260           ENDDO          IF ( l.GE.2 ) j = nyIn+l-2
261           NorthValue = NorthValue / nxIn          IF ( ABS(y_in(j)).EQ.yPole ) THEN
262           DO i=-1,nxIn+2           IF (method.EQ.1 .OR. method.EQ.2) THEN
263            arrayin(i,j) = NorthValue            poleValue = 0.
264           ENDDO            DO i=1,nxIn
265          ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN             poleValue = poleValue + arrayin(i,j)
266           DO i=-1,nxIn+2            ENDDO
267            arrayin(i,j) = 0.            poleValue = poleValue / nxIn
268           ENDDO            DO i=-1,nxIn+2
269               arrayin(i,j) = poleValue
270              ENDDO
271    #ifdef ALLOW_DEBUG
272              prtPole(l) = prtPole(l) + 0.1
273    #endif
274    #ifdef EXF_USE_OLD_INTERP_POLE
275             ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
276              DO i=-1,nxIn+2
277               arrayin(i,j) = 0.
278              ENDDO
279    #ifdef ALLOW_DEBUG
280              prtPole(l) = prtPole(l) + 0.9
281    #endif
282    #endif /* EXF_USE_OLD_INTERP_POLE */
283             ENDIF
284          ENDIF          ENDIF
285         ENDIF         ENDDO
286        ENDDO        ENDIF
287    
288        DO bj = myByLo(myThid), myByHi(myThid)  #ifdef ALLOW_DEBUG
289         DO bi = myBxLo(myThid), myBxHi(myThid)        debugFlag = ( exf_debugLev.GE.debLevC )
290         &       .OR. ( exf_debugLev.GE.debLevB .AND. myIter.LE.nIter0 )
291    C     prtPole(l)=0 : extended, =1 : changed to pole, =2 : changed to symetric
292          IF ( debugFlag ) THEN
293            l = ILNBLNK(inFile)
294            _BEGIN_MASTER(myThid)
295            WRITE(msgBuf,'(3A,I6,A,2L5)')
296         &   ' EXF_INTERP: file="',inFile(1:l),'", rec=', irecord,
297         &   ' , x-Per,P.Sym=', xIsPeriodic, poleSymmetry
298            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
299         &                      SQUEEZE_RIGHT, myThid )
300            WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' S.edge (j=-1,0,1) :',
301         &   ' proc=', (prtPole(j),j=-1,1), ', yIn=', (y_in(j),j=-1,1)
302            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
303         &                      SQUEEZE_RIGHT, myThid )
304            WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' N.edge (j=+0,+1,+2)',
305         &   ' proc=', (prtPole(j),j=2,4), ', yIn=',(y_in(j),j=nyIn,nyIn+2)
306            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
307         &                      SQUEEZE_RIGHT, myThid )
308            _END_MASTER(myThid)
309          ENDIF
310    #endif /* ALLOW_DEBUG */
311    
312  C--   Check validity of input/output coordinates  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313    C---  Prepare output grid and interpolate for each tile
314    
315    C--   put xG in interval [ lon_0 , lon_0+360 [
316          DO bj=myByLo(myThid),myByHi(myThid)
317           DO bi=myBxLo(myThid),myBxHi(myThid)
318            DO j=1-OLy,sNy+OLy
319             DO i=1-OLx,sNx+OLx
320              xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
321         &                  + threeSixtyRS*2.
322              xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)
323             ENDDO
324            ENDDO
325  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
326          IF ( debugLevel.GE.debLevC ) THEN  C--   Check validity of input/output coordinates
327            IF ( debugFlag ) THEN
328           DO j=1,sNy           DO j=1,sNy
329            DO i=1,sNx            DO i=1,sNx
330             IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.             IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
# Line 270  C--   Check validity of input/output coo Line 353  C--   Check validity of input/output coo
353           ENDDO           ENDDO
354          ENDIF          ENDIF
355  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
356           ENDDO
357          ENDDO
358    
359  C--   Compute interpolation indices        DO bj = myByLo(myThid), myByHi(myThid)
360  #ifdef OLD_EXF_INTERP_LAT_INDEX         DO bi = myBxLo(myThid), myBxHi(myThid)
361          DO j=1,sNy  
362           DO i=1,sNx  C--   Compute interpolation lon & lat index mapping
           IF (xG(i,j,bi,bj)-x_in(1) .GE. 0.) THEN  
            w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc) + 1  
           ELSE  
            w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(1))/lon_inc)  
           ENDIF  
          ENDDO  
         ENDDO  
 #ifndef TARGET_NEC_SX  
 C-    use the original and more readable variant of the algorithm that  
 C     has unvectorizable while-loops for each (i,j)  
         DO j=1,sNy  
          DO i=1,sNx  
           js = nyIn*.5  
           DO WHILE (yG(i,j,bi,bj) .LT. y_in(js))  
            js = (js - 1)*.5  
           ENDDO  
           DO WHILE (yG(i,j,bi,bj) .GE. y_in(js+1))  
            js = js + 1  
           ENDDO  
           s_ind(i,j) = js  
          ENDDO  
         ENDDO  
 #else /* TARGET_NEC_SX defined */  
 C-    this variant vectorizes more efficiently than the original one because  
 C     it moves the while loops out of the i,j loops (loop pushing) but  
 C     it is ugly and incomprehensible  
         icnt = 0  
         DO j=1,sNy  
          DO i=1,sNx  
           s_ind(i,j) = nyIn*.5  
           icnt = icnt+1  
           inx(icnt,1) = i  
           inx(icnt,2) = j  
          ENDDO  
         ENDDO  
         DO WHILE (icnt .GT. 0)  
          ii = 0  
 !CDIR NODEP  
          DO ic=1,icnt  
           i = inx(ic,1)  
           j = inx(ic,2)  
           IF (yG(i,j,bi,bj) .LT. y_in(s_ind(i,j))) THEN  
            s_ind(i,j) = (s_ind(i,j) - 1)*.5  
            ii = ii+1  
            inx(ii,1) = i  
            inx(ii,2) = j  
           ENDIF  
          ENDDO  
          icnt = ii  
         ENDDO  
         icnt = 0  
         DO j=1,sNy  
          DO i=1,sNx  
           icnt = icnt+1  
           inx(icnt,1) = i  
           inx(icnt,2) = j  
          ENDDO  
         ENDDO  
         DO WHILE (icnt .GT. 0)  
          ii = 0  
 !CDIR NODEP  
          DO ic=1,icnt  
           i = inx(ic,1)  
           j = inx(ic,2)  
           IF (yG(i,j,bi,bj) .GE. y_in(s_ind(i,j)+1)) THEN  
            s_ind(i,j) = s_ind(i,j) + 1  
            ii = ii+1  
            inx(ii,1) = i  
            inx(ii,2) = j  
           ENDIF  
          ENDDO  
          icnt = ii  
         ENDDO  
 #endif /* TARGET_NEC_SX defined */  
 #else /* OLD_EXF_INTERP_LAT_INDEX */  
363  C--     latitude index  C--     latitude index
364          DO j=1,sNy          DO j=1,sNy
365           DO i=1,sNx           DO i=1,sNx
# Line 359  C--     latitude index Line 369  C--     latitude index
369          ENDDO          ENDDO
370  C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as  C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
371  C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey  C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
372          nLoop = 1 + INT( LOG(DFLOAT(nyIn)+1. _d -3)/LOG(2. _d 0) )          tmpVar = nyIn + 1. _d -3
373            nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
374          DO l=1,nLoop          DO l=1,nLoop
375           DO j=1,sNy           DO j=1,sNy
376            DO i=1,sNx            DO i=1,sNx
# Line 375  C       1 + truncated log2(# interval -1 Line 386  C       1 + truncated log2(# interval -1
386           ENDDO           ENDDO
387          ENDDO          ENDDO
388  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
389          IF ( debugLevel.GE.debLevC ) THEN          IF ( debugFlag ) THEN
390  C-      Check that we found the right lat. index  C-      Check that we found the right lat. index
391           DO j=1,sNy           DO j=1,sNy
392            DO i=1,sNx            DO i=1,sNx
# Line 386  C-      Check that we found the right la Line 397  C-      Check that we found the right la
397       &        ', nLoop=', nLoop       &        ', nLoop=', nLoop
398              CALL PRINT_ERROR( msgBuf, myThid )              CALL PRINT_ERROR( msgBuf, myThid )
399              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
400       &        'EXF_INTERP: did not found Latitude index for grid-pt:'       &        'EXF_INTERP: did not find Latitude index for grid-pt:'
401              CALL PRINT_ERROR( msgBuf, myThid )              CALL PRINT_ERROR( msgBuf, myThid )
402              WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')              WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
403       &        'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)       &        'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
# Line 409  C--     longitude index Line 420  C--     longitude index
420             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1             w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
421           ENDDO           ENDDO
422          ENDDO          ENDDO
 #endif /* ndef OLD_EXF_INTERP_LAT_INDEX */  
423    
424          IF (method.EQ.1 .OR. method.EQ.11 .OR. method.EQ.21) THEN  C--   Do interpolation using lon & lat index mapping
425            CALL EXF_INTERPOLATE(
426  C--   Bilinear interpolation       I                inFile, irecord, method,
427           sp = 2       I                nxIn, nyIn, x_in, y_in,
428           DO j=1,sNy       I                arrayin,
429            DO i=1,sNx       O                arrayout,
430             arrayout(i,j,bi,bj) = 0.       I                xG, yG,
431             DO l=0,1       I                w_ind, s_ind,
432              px_ind(l+1) = x_in(w_ind(i,j)+l)       I                bi, bj, myThid )
             py_ind(l+1) = y_in(s_ind(i,j)+l)  
            ENDDO  
 #ifndef TARGET_NEC_SX  
            DO k=1,2  
             ew_val(k) = arrayin(w_ind(i,j)  ,s_ind(i,j)+k-1)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-1)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
             arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)  
      &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)  
            ENDDO  
 #else  
            ew_val1 = arrayin(w_ind(i,j)  ,s_ind(i,j)  )  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
            ew_val2 = arrayin(w_ind(i,j)  ,s_ind(i,j)+1)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
            arrayout(i,j,bi,bj)=  
      &            +ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)  
      &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)  
 #endif /* TARGET_NEC_SX defined */  
           ENDDO  
          ENDDO  
         ELSEIF (method .EQ. 2 .OR. method.EQ.12 .OR. method.EQ.22) THEN  
433    
 C--   Bicubic interpolation  
          sp = 4  
          DO j=1,sNy  
           DO i=1,sNx  
            arrayout(i,j,bi,bj) = 0.  
            DO l=-1,2  
             px_ind(l+2) = x_in(w_ind(i,j)+l)  
             py_ind(l+2) = y_in(s_ind(i,j)+l)  
            ENDDO  
 #ifndef TARGET_NEC_SX  
            DO k=1,4  
             ew_val(k) = arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &                + arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
      &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-2)  
      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)  
      &                + arrayin(w_ind(i,j)+2,s_ind(i,j)+k-2)  
      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)  
             arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)  
      &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)  
            ENDDO  
 #else  
            ew_val1 = arrayin(w_ind(i,j)-1,s_ind(i,j)-1)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)-1)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)-1)  
      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)-1)  
      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)  
            ew_val2 = arrayin(w_ind(i,j)-1,s_ind(i,j)  )  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)  )  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )  
      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)  )  
      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)  
            ew_val3 = arrayin(w_ind(i,j)-1,s_ind(i,j)+1)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+1)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)  
      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+1)  
      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)  
            ew_val4 = arrayin(w_ind(i,j)-1,s_ind(i,j)+2)  
      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+2)  
      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+2)  
      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)  
      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+2)  
      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)  
            arrayout(i,j,bi,bj) =  
      &             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)  
 #endif /* TARGET_NEC_SX defined */  
           ENDDO  
          ENDDO  
         ELSE  
          l = ILNBLNK(inFile)  
          WRITE(msgBuf,'(3A,I6)')  
      &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord  
          CALL PRINT_ERROR( msgBuf, myThid )  
          WRITE(msgBuf,'(A,I8,A)')  
      &        'EXF_INTERP: method=', method,' not supported'  
          CALL PRINT_ERROR( msgBuf, myThid )  
          STOP 'ABNORMAL END: S/R EXF_INTERP: invalid method'  
         ENDIF  
434         ENDDO         ENDDO
435        ENDDO        ENDDO
436    
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22