/[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.33 by mlosch, Thu Apr 16 16:13:30 2015 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,
      I                filePrec,  
13       O                arrayout,       O                arrayout,
14       I                irecord, xG_in, yG,       I                irecord, xG_in, yG,
15       I                lon_0, lon_inc,       I                lon_0, lon_inc, lat_0, lat_inc,
16       I                lat_0, lat_inc,       I                nxIn, nyIn, method, myIter, myThid )
      I                nxIn, nyIn, method, myThid )  
17    
18  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
19  C  *==========================================================*  C  *==========================================================*
# Line 75  C     === Global variables === Line 29  C     === Global variables ===
29  #include "SIZE.h"  #include "SIZE.h"
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32    #ifdef ALLOW_DEBUG
33    # include "EXF_PARAM.h"
34    #endif
35    
36  C !INPUT/OUTPUT PARAMETERS:  C !INPUT/OUTPUT PARAMETERS:
37  C  inFile      (string)  :: name of the binary input file (direct access)  C  inFile      (string)  :: name of the binary input file (direct access)
# Line 89  C     lat_inc            :: vector y-gri Line 45  C     lat_inc            :: vector y-gri
45  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
46  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic  C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
47  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
48    C  myIter      (integer) :: current iteration number
49  C  myThid      (integer) :: My Thread Id number  C  myThid      (integer) :: My Thread Id number
50    
51        CHARACTER*(*) infile        CHARACTER*(*) inFile
52        INTEGER       filePrec, irecord, nxIn, nyIn        INTEGER       filePrec, irecord, nxIn, nyIn
53        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54        _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)
# Line 99  C  myThid      (integer) :: My Thread Id Line 56  C  myThid      (integer) :: My Thread Id
56        _RL           lon_0, lon_inc        _RL           lon_0, lon_inc
57  c     _RL           lat_0, lat_inc(nyIn-1)  c     _RL           lat_0, lat_inc(nyIn-1)
58        _RL           lat_0, lat_inc(*)        _RL           lat_0, lat_inc(*)
59        INTEGER       method, myThid        INTEGER       method, myIter, myThid
60    
61  C !FUNCTIONS:  C !FUNCTIONS:
62        EXTERNAL LAGRAN  #ifdef ALLOW_DEBUG
       _RL      LAGRAN  
63        INTEGER  ILNBLNK        INTEGER  ILNBLNK
64        EXTERNAL ILNBLNK        EXTERNAL ILNBLNK
65    #endif
66    
67  C !LOCAL VARIABLES:  C !LOCAL VARIABLES:
68  C     msgBuf      :: Informational/error message buffer  C     arrayin     :: input field array (loaded from file)
69    C     x_in        :: longitude vector defining input field grid
70    C     y_in        :: latitude  vector defining input field grid
71    C     w_ind       :: input field longitudinal index, on western side of model grid pt
72    C     s_ind       :: input field latitudinal index, on southern side of model grid pt
73  C     bi, bj      :: tile indices  C     bi, bj      :: tile indices
74        CHARACTER*(MAX_LEN_MBUF) msgBuf  C     i, j, k, l  :: loop indices
75        INTEGER  bi, bj  C     msgBuf      :: Informational/error message buffer
       INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)  
       _RL      px_ind(4), py_ind(4), ew_val(4)  
76        _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )        _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )
77        _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)        _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)
78        _RL      NorthValue        INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)
79        INTEGER  i, j, k, l, sp        INTEGER  bi, bj
80  #ifdef OLD_EXF_INTERP_LAT_INDEX        INTEGER  i, j, k, l
       INTEGER  js  
 #else  
81        INTEGER  nLoop        INTEGER  nLoop
82  #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  
83        _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. )  
84        _RS      threeSixtyRS        _RS      threeSixtyRS
85          _RL      yPole, symSign, poleValue
86        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
87        LOGICAL  xIsPeriodic        PARAMETER ( yPole = 90. )
88          INTEGER  nxd2
89          LOGICAL  xIsPeriodic, poleSymmetry
90    #ifdef ALLOW_DEBUG
91          LOGICAL  debugFlag
92          CHARACTER*(MAX_LEN_MBUF) msgBuf
93          _RS      prtPole(-1:4)
94    #endif
95  CEOP  CEOP
96    
97  C--   put xG in interval [ lon_0 , lon_0+360 [  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
98        DO bj=myByLo(myThid),myByHi(myThid)  C---  Load inut field
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0  
      &                  + threeSixtyRS*2.  
           xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
99    
 C--   Load inut field  
100        CALL EXF_INTERP_READ(        CALL EXF_INTERP_READ(
101       I         inFile, filePrec,       I         inFile, filePrec,
102       O         arrayin,       O         arrayin,
103       I         irecord, nxIn, nyIn, myThid )       I         irecord, nxIn, nyIn, myThid )
104    
105    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
106    C---  Prepare input grid and input field
107    
108  C--   setup input longitude grid  C--   setup input longitude grid
109        DO i=-1,nxIn+2        DO i=-1,nxIn+2
110         x_in(i) = lon_0 + (i-1)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
111        ENDDO        ENDDO
112        xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )        xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
113          nxd2 = NINT( nxIn*0.5 )
114          poleSymmetry = xIsPeriodic .AND. ( nxIn.EQ.2*nxd2 )
115    #ifdef EXF_USE_OLD_INTERP_POLE
116          poleSymmetry = .FALSE.
117    #endif
118    
119  C--   setup input latitude grid  C--   setup input latitude grid
120        y_in(1) = lat_0        y_in(1) = lat_0
# Line 167  C--   setup input latitude grid Line 122  C--   setup input latitude grid
122         i = MIN(j,nyIn-1)         i = MIN(j,nyIn-1)
123         y_in(j+1) = y_in(j) + lat_inc(i)         y_in(j+1) = y_in(j) + lat_inc(i)
124        ENDDO        ENDDO
 C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole  
125        y_in(0) = y_in(1) - lat_inc(1)        y_in(0) = y_in(1) - lat_inc(1)
126        y_in(-1)= y_in(0) - lat_inc(1)        y_in(-1)= y_in(0) - lat_inc(1)
127  c     IF ( y_in(1).GT.-ninety .AND. y_in(0).LT.-ninety ) THEN  #ifdef ALLOW_DEBUG
128  c       y_in(0) = -ninety        DO l=-1,4
129  c       y_in(-1) = -2.*ninety - y_in(1)          prtPole(l) = 0.
130  c     ENDIF        ENDDO
131  c     IF ( y_in(0).GT.-ninety .AND. y_in(-1).LT.-ninety ) THEN  #endif
132  c       y_in(-1) = -ninety  C--   For tracer (method=1,2) add 1 row @ the pole (if last row is not)
133  c     ENDIF  C     and will fill it with the polarmost-row zonally averaged value.
134    C     For vector component, cannot do much without the other component;
135    C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
136    #ifdef EXF_USE_OLD_INTERP_POLE
137          IF ( .TRUE. ) THEN
138    #else
139          IF ( method.LT.10 ) THEN
140    C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
141           j = 0
142           IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
143             y_in(j) = -yPole
144             y_in(j-1) = -2.*yPole - y_in(1)
145    #ifdef ALLOW_DEBUG
146             prtPole(j)   = 1.
147             prtPole(j-1) = 2.
148    #endif
149           ENDIF
150           j = -1
151           IF ( ABS(y_in(j+1)).GT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
152             y_in(j) = -yPole
153    #ifdef ALLOW_DEBUG
154             prtPole(j)   = 1.
155    #endif
156           ENDIF
157    #endif /* EXF_USE_OLD_INTERP_POLE */
158  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
159        j = nyIn+1         j = nyIn+1
160        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
161          y_in(j) = ninety           y_in(j) = yPole
162          y_in(j+1) = 2.*ninety - y_in(j-1)           y_in(j+1) = 2.*yPole - y_in(j-1)
163        ENDIF  #ifdef ALLOW_DEBUG
164        j = nyIn+2           prtPole(3)   = 1.
165        IF ( y_in(j-1).LT.ninety .AND. y_in(j).GT.ninety ) THEN           prtPole(3+1) = 2.
166          y_in(j) = ninety  #endif
167           ENDIF
168           j = nyIn+2
169           IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
170             y_in(j) = yPole
171    #ifdef ALLOW_DEBUG
172             prtPole(4)   = 1.
173    #endif
174           ENDIF
175        ENDIF        ENDIF
176    
177  C--   enlarge boundary  C--   Enlarge boundary
178        IF ( xIsPeriodic ) THEN        IF ( xIsPeriodic ) THEN
179    C-    fill-in added column assuming periodicity
180          DO j=1,nyIn          DO j=1,nyIn
181           arrayin( 0,j)     = arrayin(nxIn  ,j)           arrayin( 0,j)     = arrayin(nxIn  ,j)
182           arrayin(-1,j)     = arrayin(nxIn-1,j)           arrayin(-1,j)     = arrayin(nxIn-1,j)
# Line 197  C--   enlarge boundary Line 184  C--   enlarge boundary
184           arrayin(nxIn+2,j) = arrayin(2,j)           arrayin(nxIn+2,j) = arrayin(2,j)
185          ENDDO          ENDDO
186        ELSE        ELSE
187    C-    fill-in added column from nearest column
188          DO j=1,nyIn          DO j=1,nyIn
189           arrayin( 0,j)     = arrayin(1,j)           arrayin( 0,j)     = arrayin(1,j)
190           arrayin(-1,j)     = arrayin(1,j)           arrayin(-1,j)     = arrayin(1,j)
# Line 204  C--   enlarge boundary Line 192  C--   enlarge boundary
192           arrayin(nxIn+2,j) = arrayin(nxIn,j)           arrayin(nxIn+2,j) = arrayin(nxIn,j)
193          ENDDO          ENDDO
194        ENDIF        ENDIF
195        DO i=-1,nxIn+2        symSign = 1. _d 0
196           arrayin(i, 0)     = arrayin(i,1)        IF ( method.GE.10 ) symSign = -1. _d 0
197           arrayin(i,-1)     = arrayin(i,1)        DO l=-1,2
198           arrayin(i,nyIn+1) = arrayin(i,nyIn)         j = l
199           arrayin(i,nyIn+2) = arrayin(i,nyIn)         IF ( l.GE.1 ) j = nyIn+l
200           k = MAX(1,MIN(j,nyIn))
201           IF ( poleSymmetry .AND. ABS(y_in(j)).GT.yPole ) THEN
202    C-    fill-in added row assuming pole-symmetry
203            DO i=-1,nxd2
204             arrayin(i,j) = symSign*arrayin(i+nxd2,k)
205            ENDDO
206            DO i=1,nxd2+2
207             arrayin(i+nxd2,j) = symSign*arrayin(i,k)
208            ENDDO
209    #ifdef ALLOW_DEBUG
210            i = l + 2*( (l+1)/2 )
211            prtPole(i) = prtPole(i) + 0.2
212    #endif
213           ELSE
214    C-    fill-in added row from nearest column values
215            DO i=-1,nxIn+2
216             arrayin(i,j) = arrayin(i,k)
217            ENDDO
218           ENDIF
219        ENDDO        ENDDO
220    
221  C-    For tracer (method=1,2) set to northernmost zonal-mean value  C--   For tracer (method=1,2) set to northernmost zonal-mean value
222  C     at 90N to avoid sharp zonal gradients near the Pole.  C     at 90N to avoid sharp zonal gradients near the Pole.
223  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;
224  C     gradient at North Pole  C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
225  C     For V (method=11,12) set to northernmost zonal value at 90N,  #ifdef EXF_USE_OLD_INTERP_POLE
226  C     as is already done above in order to allow cross-PoleArctic flow        IF ( .TRUE. ) THEN
227        DO j=nyIn,nyIn+2         DO l= 3,4
228         IF (y_in(j).EQ.ninety) THEN  #else
229          IF (method.EQ.1 .OR. method.EQ.2) THEN        IF ( method.LT.10 ) THEN
230           NorthValue = 0.         DO l=-1,4
231           DO i=1,nxIn  #endif
232            NorthValue = NorthValue + arrayin(i,j)          j = l
233           ENDDO          IF ( l.GE.2 ) j = nyIn+l-2
234           NorthValue = NorthValue / nxIn          IF ( ABS(y_in(j)).EQ.yPole ) THEN
235           DO i=-1,nxIn+2           IF (method.EQ.1 .OR. method.EQ.2) THEN
236            arrayin(i,j) = NorthValue            poleValue = 0.
237           ENDDO            DO i=1,nxIn
238          ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN             poleValue = poleValue + arrayin(i,j)
239           DO i=-1,nxIn+2            ENDDO
240            arrayin(i,j) = 0.            poleValue = poleValue / nxIn
241           ENDDO            DO i=-1,nxIn+2
242               arrayin(i,j) = poleValue
243              ENDDO
244    #ifdef ALLOW_DEBUG
245              prtPole(l) = prtPole(l) + 0.1
246    #endif
247    #ifdef EXF_USE_OLD_INTERP_POLE
248             ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
249              DO i=-1,nxIn+2
250               arrayin(i,j) = 0.
251              ENDDO
252    #ifdef ALLOW_DEBUG
253              prtPole(l) = prtPole(l) + 0.9
254    #endif
255    #endif /* EXF_USE_OLD_INTERP_POLE */
256             ENDIF
257          ENDIF          ENDIF
258         ENDIF         ENDDO
259        ENDDO        ENDIF
260    
261        DO bj = myByLo(myThid), myByHi(myThid)  #ifdef ALLOW_DEBUG
262         DO bi = myBxLo(myThid), myBxHi(myThid)        debugFlag = ( exf_debugLev.GE.debLevC )
263         &       .OR. ( exf_debugLev.GE.debLevB .AND. myIter.LE.nIter0 )
264    C     prtPole(l)=0 : extended, =1 : changed to pole, =2 : changed to symetric
265          IF ( debugFlag ) THEN
266            l = ILNBLNK(inFile)
267            _BEGIN_MASTER(myThid)
268            WRITE(msgBuf,'(3A,I6,A,2L5)')
269         &   'EXF_INTERP: file="',inFile(1:l),'", rec=', irecord,
270         &   ' , x-Per,P.Sym=', xIsPeriodic, poleSymmetry
271            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
272         &                      SQUEEZE_RIGHT, myThid )
273            WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') 'S.edge (j=-1,0,1) :',
274         &   ' proc=', (prtPole(j),j=-1,1), ', yIn=', (y_in(j),j=-1,1)
275            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
276         &                      SQUEEZE_RIGHT, myThid )
277            WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') 'N.edge (j=+0,+1,+2)',
278         &   ' proc=', (prtPole(j),j=2,4), ', yIn=',(y_in(j),j=nyIn,nyIn+2)
279            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
280         &                      SQUEEZE_RIGHT, myThid )
281            _END_MASTER(myThid)
282          ENDIF
283    #endif /* ALLOW_DEBUG */
284    
285  C--   Check validity of input/output coordinates  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
286    C---  Prepare output grid and interpolate for each tile
287    
288    C--   put xG in interval [ lon_0 , lon_0+360 [
289          DO bj=myByLo(myThid),myByHi(myThid)
290           DO bi=myBxLo(myThid),myBxHi(myThid)
291            DO j=1-OLy,sNy+OLy
292             DO i=1-OLx,sNx+OLx
293              xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
294         &                  + threeSixtyRS*2.
295              xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)
296             ENDDO
297            ENDDO
298  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
299          IF ( debugLevel.GE.debLevC ) THEN  C--   Check validity of input/output coordinates
300            IF ( debugFlag ) THEN
301           DO j=1,sNy           DO j=1,sNy
302            DO i=1,sNx            DO i=1,sNx
303             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 326  C--   Check validity of input/output coo
326           ENDDO           ENDDO
327          ENDIF          ENDIF
328  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
329           ENDDO
330          ENDDO
331    
332  C--   Compute interpolation indices        DO bj = myByLo(myThid), myByHi(myThid)
333  #ifdef OLD_EXF_INTERP_LAT_INDEX         DO bi = myBxLo(myThid), myBxHi(myThid)
334          DO j=1,sNy  
335           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 */  
336  C--     latitude index  C--     latitude index
337          DO j=1,sNy          DO j=1,sNy
338           DO i=1,sNx           DO i=1,sNx
# Line 359  C--     latitude index Line 342  C--     latitude index
342          ENDDO          ENDDO
343  C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as  C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
344  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
345          nLoop = 1 + INT( LOG(DFLOAT(nyIn)+1. _d -3)/LOG(2. _d 0) )          tmpVar = nyIn + 1. _d -3
346            nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
347          DO l=1,nLoop          DO l=1,nLoop
348           DO j=1,sNy           DO j=1,sNy
349            DO i=1,sNx            DO i=1,sNx
# Line 375  C       1 + truncated log2(# interval -1 Line 359  C       1 + truncated log2(# interval -1
359           ENDDO           ENDDO
360          ENDDO          ENDDO
361  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
362          IF ( debugLevel.GE.debLevC ) THEN          IF ( debugFlag ) THEN
363  C-      Check that we found the right lat. index  C-      Check that we found the right lat. index
364           DO j=1,sNy           DO j=1,sNy
365            DO i=1,sNx            DO i=1,sNx
# Line 386  C-      Check that we found the right la Line 370  C-      Check that we found the right la
370       &        ', nLoop=', nLoop       &        ', nLoop=', nLoop
371              CALL PRINT_ERROR( msgBuf, myThid )              CALL PRINT_ERROR( msgBuf, myThid )
372              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
373       &        'EXF_INTERP: did not found Latitude index for grid-pt:'       &        'EXF_INTERP: did not find Latitude index for grid-pt:'
374              CALL PRINT_ERROR( msgBuf, myThid )              CALL PRINT_ERROR( msgBuf, myThid )
375              WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')              WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
376       &        '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 393  C--     longitude index
393             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
394           ENDDO           ENDDO
395          ENDDO          ENDDO
 #endif /* ndef OLD_EXF_INTERP_LAT_INDEX */  
   
         IF (method.EQ.1 .OR. method.EQ.11 .OR. method.EQ.21) THEN  
396    
397  C--   Bilinear interpolation  C--   Do interpolation using lon & lat index mapping
398           sp = 2          CALL EXF_INTERPOLATE(
399           DO j=1,sNy       I                inFile, irecord, method,
400            DO i=1,sNx       I                nxIn, nyIn, x_in, y_in,
401             arrayout(i,j,bi,bj) = 0.       O                arrayin,
402             DO l=0,1       O                arrayout,
403              px_ind(l+1) = x_in(w_ind(i,j)+l)       I                xG, yG,
404              py_ind(l+1) = y_in(s_ind(i,j)+l)       I                w_ind, s_ind,
405             ENDDO       I                bi, bj, myThid )
 #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  
406    
 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  
407         ENDDO         ENDDO
408        ENDDO        ENDDO
409    
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22