pkg/exf/exf_interp.F

6d54cf9ca1 Ed H*0001 #include "EXF_OPTIONS.h"
981f32514d Dimi*0002
b3be92d084 Jean*0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
981f32514d Dimi*0004
ef1828fb24 Jean*0005 CBOP
                0006 C     !ROUTINE: EXF_INTERP
                0007 C     !INTERFACE:
                0008        SUBROUTINE EXF_INTERP(
1bdb364c26 Jean*0009      I                inFile, filePrec,
a9085e980c Jean*0010 #ifdef EXF_INTERP_USE_DYNALLOC
ef1828fb24 Jean*0011      O                arrayout,
a9085e980c Jean*0012 #else
                0013      O                arrayout, arrayin,
                0014 #endif
ef1828fb24 Jean*0015      I                irecord, xG_in, yG,
1bdb364c26 Jean*0016      I                lon_0, lon_inc, lat_0, lat_inc,
                0017      I                nxIn, nyIn, method, myIter, myThid )
ef1828fb24 Jean*0018
                0019 C !DESCRIPTION: \bv
                0020 C  *==========================================================*
                0021 C  | SUBROUTINE EXF_INTERP
                0022 C  | o Load from file a regular lat-lon input field
                0023 C  |   and interpolate on to the model grid location
                0024 C  *==========================================================*
                0025 C \ev
                0026
                0027 C !USES:
                0028       IMPLICIT NONE
                0029 C     === Global variables ===
                0030 #include "SIZE.h"
                0031 #include "EEPARAMS.h"
                0032 #include "PARAMS.h"
a9085e980c Jean*0033 #include "EXF_INTERP_SIZE.h"
9f46642c85 Jean*0034 #ifdef ALLOW_DEBUG
                0035 # include "EXF_PARAM.h"
                0036 #endif
981f32514d Dimi*0037
ef1828fb24 Jean*0038 C !INPUT/OUTPUT PARAMETERS:
                0039 C  inFile      (string)  :: name of the binary input file (direct access)
3e6ce1cd50 Jean*0040 C  filePrec    (integer) :: number of bits per word in file (32 or 64)
ef1828fb24 Jean*0041 C  arrayout    ( _RL )   :: output array
a9085e980c Jean*0042 #ifndef EXF_INTERP_USE_DYNALLOC
                0043 C  arrayin     ( _RL )   :: input field array (loaded from file)
                0044 #endif
3e6ce1cd50 Jean*0045 C  irecord     (integer) :: record number to read
ef1828fb24 Jean*0046 C     xG_in,yG           :: coordinates for output grid to interpolate to
3e6ce1cd50 Jean*0047 C     lon_0, lat_0       :: lon and lat of sw corner of global input grid
                0048 C     lon_inc            :: scalar x-grid increment
                0049 C     lat_inc            :: vector y-grid increments
ef1828fb24 Jean*0050 C  nxIn,nyIn   (integer) :: size in x & y direction of input file to read
3e6ce1cd50 Jean*0051 C     method             :: 1,11,21 for bilinear; 2,12,22 for bicubic
                0052 C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
1bdb364c26 Jean*0053 C  myIter      (integer) :: current iteration number
3e6ce1cd50 Jean*0054 C  myThid      (integer) :: My Thread Id number
1bdb364c26 Jean*0055       CHARACTER*(*) inFile
ef1828fb24 Jean*0056       INTEGER       filePrec, irecord, nxIn, nyIn
57bf62b584 Dimi*0057       _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
a9085e980c Jean*0058 #ifndef EXF_INTERP_USE_DYNALLOC
                0059       _RL           arrayin ( -1:nxIn+2, -1:nyIn+2 )
                0060 #endif
3e6ce1cd50 Jean*0061       _RS           xG_in   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57bf62b584 Dimi*0062       _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0063       _RL           lon_0, lon_inc
ef1828fb24 Jean*0064 c     _RL           lat_0, lat_inc(nyIn-1)
b3be92d084 Jean*0065       _RL           lat_0, lat_inc(*)
1bdb364c26 Jean*0066       INTEGER       method, myIter, myThid
981f32514d Dimi*0067
ef1828fb24 Jean*0068 C !FUNCTIONS:
1bdb364c26 Jean*0069 #ifdef ALLOW_DEBUG
ef1828fb24 Jean*0070       INTEGER  ILNBLNK
                0071       EXTERNAL ILNBLNK
1bdb364c26 Jean*0072 #endif
3e6ce1cd50 Jean*0073
ef1828fb24 Jean*0074 C !LOCAL VARIABLES:
1bdb364c26 Jean*0075 C     x_in        :: longitude vector defining input field grid
                0076 C     y_in        :: latitude  vector defining input field grid
                0077 C     w_ind       :: input field longitudinal index, on western side of model grid pt
                0078 C     s_ind       :: input field latitudinal index, on southern side of model grid pt
ef1828fb24 Jean*0079 C     bi, bj      :: tile indices
1bdb364c26 Jean*0080 C     i, j, k, l  :: loop indices
                0081 C     msgBuf      :: Informational/error message buffer
a9085e980c Jean*0082 #ifdef EXF_INTERP_USE_DYNALLOC
                0083 C     arrayin     :: input field array (loaded from file)
ef1828fb24 Jean*0084       _RL      arrayin( -1:nxIn+2, -1:nyIn+2 )
                0085       _RL      x_in(-1:nxIn+2), y_in(-1:nyIn+2)
a9085e980c Jean*0086 #else /* EXF_INTERP_USE_DYNALLOC */
                0087       _RL      x_in(-1:exf_max_nLon+2), y_in(-1:exf_max_nLat+2)
                0088 #endif /* EXF_INTERP_USE_DYNALLOC */
1bdb364c26 Jean*0089       INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)
                0090       INTEGER  bi, bj
                0091       INTEGER  i, j, k, l
ef1828fb24 Jean*0092       INTEGER  nLoop
1806bc85ed Jean*0093       _RL      tmpVar
8f3c818576 Patr*0094       _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
3e6ce1cd50 Jean*0095       _RS      threeSixtyRS
1bdb364c26 Jean*0096       _RL      yPole, symSign, poleValue
e2ffdbe062 Jean*0097       _RL      edgeFac, poleFac
8f3c818576 Patr*0098       PARAMETER ( threeSixtyRS = 360. )
1bdb364c26 Jean*0099       PARAMETER ( yPole = 90. )
                0100       INTEGER  nxd2
b374dc700b Mart*0101       LOGICAL  xIsPeriodic, poleSymmetry
1bdb364c26 Jean*0102 #ifdef ALLOW_DEBUG
b374dc700b Mart*0103       LOGICAL  debugFlag
1bdb364c26 Jean*0104       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0105       _RS      prtPole(-1:4)
                0106 #endif
ef1828fb24 Jean*0107 CEOP
9e9a5e4e0c Patr*0108
a9085e980c Jean*0109 #ifndef EXF_INTERP_USE_DYNALLOC
                0110 C--   Check size declaration:
                0111       IF ( nxIn.GT.exf_max_nLon ) THEN
                0112        STOP 'EXF_INTERP: exf_max_nLon too small'
                0113       ENDIF
                0114       IF ( nyIn.GT.exf_max_nLat ) THEN
                0115        STOP 'EXF_INTERP: exf_max_nLat too small'
                0116       ENDIF
                0117       IF ( (nxIn+4)*(nyIn+4).GT.exf_interp_bufferSize ) THEN
                0118        STOP 'EXF_INTERP: exf_interp_bufferSize too small'
                0119       ENDIF
                0120 #endif /* ndef EXF_INTERP_USE_DYNALLOC */
                0121
1bdb364c26 Jean*0122 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0123 C---  Load inut field
ef1828fb24 Jean*0124
                0125       CALL EXF_INTERP_READ(
                0126      I         inFile, filePrec,
                0127      O         arrayin,
                0128      I         irecord, nxIn, nyIn, myThid )
                0129
1bdb364c26 Jean*0130 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0131 C---  Prepare input grid and input field
                0132
ef1828fb24 Jean*0133 C--   setup input longitude grid
                0134       DO i=-1,nxIn+2
9da68b3e79 Dimi*0135        x_in(i) = lon_0 + (i-1)*lon_inc
ef1828fb24 Jean*0136       ENDDO
                0137       xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
1bdb364c26 Jean*0138       nxd2 = NINT( nxIn*0.5 )
                0139       poleSymmetry = xIsPeriodic .AND. ( nxIn.EQ.2*nxd2 )
                0140 #ifdef EXF_USE_OLD_INTERP_POLE
                0141       poleSymmetry = .FALSE.
                0142 #endif
9e9a5e4e0c Patr*0143
ef1828fb24 Jean*0144 C--   setup input latitude grid
9da68b3e79 Dimi*0145       y_in(1) = lat_0
ef1828fb24 Jean*0146       DO j=1,nyIn+1
                0147        i = MIN(j,nyIn-1)
                0148        y_in(j+1) = y_in(j) + lat_inc(i)
                0149       ENDDO
                0150       y_in(0) = y_in(1) - lat_inc(1)
                0151       y_in(-1)= y_in(0) - lat_inc(1)
1bdb364c26 Jean*0152 #ifdef ALLOW_DEBUG
                0153       DO l=-1,4
                0154         prtPole(l) = 0.
                0155       ENDDO
                0156 #endif
                0157 C--   For tracer (method=1,2) add 1 row @ the pole (if last row is not)
                0158 C     and will fill it with the polarmost-row zonally averaged value.
                0159 C     For vector component, cannot do much without the other component;
                0160 C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
                0161 #ifdef EXF_USE_OLD_INTERP_POLE
                0162       IF ( .TRUE. ) THEN
                0163 #else
                0164       IF ( method.LT.10 ) THEN
                0165 C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
                0166        j = 0
                0167        IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
                0168          y_in(j) = -yPole
c61b8a5e5a Jean*0169          y_in(j-1) = -2.*yPole - y_in(j+1)
1bdb364c26 Jean*0170 #ifdef ALLOW_DEBUG
                0171          prtPole(j)   = 1.
                0172          prtPole(j-1) = 2.
                0173 #endif
                0174        ENDIF
                0175        j = -1
c61b8a5e5a Jean*0176        IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
1bdb364c26 Jean*0177          y_in(j) = -yPole
                0178 #ifdef ALLOW_DEBUG
                0179          prtPole(j)   = 1.
                0180 #endif
                0181        ENDIF
                0182 #endif /* EXF_USE_OLD_INTERP_POLE */
ef1828fb24 Jean*0183 C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
1bdb364c26 Jean*0184        j = nyIn+1
                0185        IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
                0186          y_in(j) = yPole
                0187          y_in(j+1) = 2.*yPole - y_in(j-1)
                0188 #ifdef ALLOW_DEBUG
                0189          prtPole(3)   = 1.
                0190          prtPole(3+1) = 2.
                0191 #endif
                0192        ENDIF
                0193        j = nyIn+2
                0194        IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
                0195          y_in(j) = yPole
                0196 #ifdef ALLOW_DEBUG
                0197          prtPole(4)   = 1.
                0198 #endif
                0199        ENDIF
ef1828fb24 Jean*0200       ENDIF
                0201
1bdb364c26 Jean*0202 C--   Enlarge boundary
ef1828fb24 Jean*0203       IF ( xIsPeriodic ) THEN
1bdb364c26 Jean*0204 C-    fill-in added column assuming periodicity
ef1828fb24 Jean*0205         DO j=1,nyIn
                0206          arrayin( 0,j)     = arrayin(nxIn  ,j)
                0207          arrayin(-1,j)     = arrayin(nxIn-1,j)
                0208          arrayin(nxIn+1,j) = arrayin(1,j)
                0209          arrayin(nxIn+2,j) = arrayin(2,j)
                0210         ENDDO
                0211       ELSE
1bdb364c26 Jean*0212 C-    fill-in added column from nearest column
ef1828fb24 Jean*0213         DO j=1,nyIn
                0214          arrayin( 0,j)     = arrayin(1,j)
                0215          arrayin(-1,j)     = arrayin(1,j)
                0216          arrayin(nxIn+1,j) = arrayin(nxIn,j)
                0217          arrayin(nxIn+2,j) = arrayin(nxIn,j)
                0218         ENDDO
                0219       ENDIF
1bdb364c26 Jean*0220       symSign = 1. _d 0
                0221       IF ( method.GE.10 ) symSign = -1. _d 0
                0222       DO l=-1,2
                0223        j = l
                0224        IF ( l.GE.1 ) j = nyIn+l
                0225        k = MAX(1,MIN(j,nyIn))
                0226        IF ( poleSymmetry .AND. ABS(y_in(j)).GT.yPole ) THEN
e2ffdbe062 Jean*0227         IF ( nyIn.GE.3 .AND. ABS(y_in(k)).EQ.yPole )
                0228      &    k = MAX(2,MIN(j,nyIn-1))
1bdb364c26 Jean*0229 C-    fill-in added row assuming pole-symmetry
                0230         DO i=-1,nxd2
                0231          arrayin(i,j) = symSign*arrayin(i+nxd2,k)
                0232         ENDDO
                0233         DO i=1,nxd2+2
                0234          arrayin(i+nxd2,j) = symSign*arrayin(i,k)
                0235         ENDDO
                0236 #ifdef ALLOW_DEBUG
                0237         i = l + 2*( (l+1)/2 )
                0238         prtPole(i) = prtPole(i) + 0.2
                0239 #endif
                0240        ELSE
                0241 C-    fill-in added row from nearest column values
                0242         DO i=-1,nxIn+2
                0243          arrayin(i,j) = arrayin(i,k)
                0244         ENDDO
                0245        ENDIF
ef1828fb24 Jean*0246       ENDDO
                0247
1bdb364c26 Jean*0248 C--   For tracer (method=1,2) set to northernmost zonal-mean value
6519624399 Dimi*0249 C     at 90N to avoid sharp zonal gradients near the Pole.
1bdb364c26 Jean*0250 C     For vector component, cannot do much without the other component;
                0251 C     averaging properly done if uvInterp=T in S/R EXF_INTERP_UV
                0252 #ifdef EXF_USE_OLD_INTERP_POLE
                0253       IF ( .TRUE. ) THEN
                0254        DO l= 3,4
                0255 #else
                0256       IF ( method.LT.10 ) THEN
                0257        DO l=-1,4
                0258 #endif
                0259         j = l
                0260         IF ( l.GE.2 ) j = nyIn+l-2
                0261         IF ( ABS(y_in(j)).EQ.yPole ) THEN
                0262          IF (method.EQ.1 .OR. method.EQ.2) THEN
                0263           poleValue = 0.
                0264           DO i=1,nxIn
                0265            poleValue = poleValue + arrayin(i,j)
                0266           ENDDO
                0267           poleValue = poleValue / nxIn
                0268           DO i=-1,nxIn+2
                0269            arrayin(i,j) = poleValue
                0270           ENDDO
                0271 #ifdef ALLOW_DEBUG
                0272           prtPole(l) = prtPole(l) + 0.1
                0273 #endif
                0274 #ifdef EXF_USE_OLD_INTERP_POLE
                0275          ELSEIF (method.EQ.11 .OR. method.EQ.12) THEN
                0276           DO i=-1,nxIn+2
                0277            arrayin(i,j) = 0.
                0278           ENDDO
                0279 #ifdef ALLOW_DEBUG
                0280           prtPole(l) = prtPole(l) + 0.9
                0281 #endif
                0282 #endif /* EXF_USE_OLD_INTERP_POLE */
                0283          ENDIF
ef1828fb24 Jean*0284         ENDIF
1bdb364c26 Jean*0285        ENDDO
                0286       ENDIF
e2ffdbe062 Jean*0287 #ifndef EXF_USE_OLD_INTERP_POLE
                0288       IF (method.EQ.1 .OR. method.EQ.2) THEN
                0289 C-    change first additional row from simple copy to linear interpolation
                0290 C     between nearest column values and pole value
                0291        DO l=0,1
                0292         k = l*(nyIn+3) -1
                0293         IF ( ABS(y_in(k)).EQ.yPole ) THEN
                0294          j = l*(nyIn+1)
                0295          i = l*(nyIn-1) +1
                0296          edgeFac = (y_in(j) - y_in(k)) / (y_in(i) - y_in(k))
                0297          poleFac = (y_in(i) - y_in(j)) / (y_in(i) - y_in(k))
                0298          DO i=-1,nxIn+2
                0299            arrayin(i,j) = arrayin(i,j) * edgeFac
                0300      &                  + arrayin(i,k) * poleFac
                0301          ENDDO
                0302 #ifdef ALLOW_DEBUG
                0303          prtPole(3*l) = prtPole(3*l) + 0.3
                0304 #endif
                0305         ENDIF
                0306        ENDDO
                0307       ENDIF
                0308 #endif /* EXF_USE_OLD_INTERP_POLE */
981f32514d Dimi*0309
1bdb364c26 Jean*0310 #ifdef ALLOW_DEBUG
9f46642c85 Jean*0311       debugFlag = ( exf_debugLev.GE.debLevC )
                0312      &       .OR. ( exf_debugLev.GE.debLevB .AND. myIter.LE.nIter0 )
1bdb364c26 Jean*0313 C     prtPole(l)=0 : extended, =1 : changed to pole, =2 : changed to symetric
                0314       IF ( debugFlag ) THEN
                0315         l = ILNBLNK(inFile)
2c18da0d62 Jean*0316         _BEGIN_MASTER(myThid)
1bdb364c26 Jean*0317         WRITE(msgBuf,'(3A,I6,A,2L5)')
a72a954434 Jean*0318      &   ' EXF_INTERP: file="',inFile(1:l),'", rec=', irecord,
1bdb364c26 Jean*0319      &   ' , x-Per,P.Sym=', xIsPeriodic, poleSymmetry
                0320         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0321      &                      SQUEEZE_RIGHT, myThid )
a72a954434 Jean*0322         WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' S.edge (j=-1,0,1) :',
1bdb364c26 Jean*0323      &   ' proc=', (prtPole(j),j=-1,1), ', yIn=', (y_in(j),j=-1,1)
                0324         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0325      &                      SQUEEZE_RIGHT, myThid )
a72a954434 Jean*0326         WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' N.edge (j=+0,+1,+2)',
1bdb364c26 Jean*0327      &   ' proc=', (prtPole(j),j=2,4), ', yIn=',(y_in(j),j=nyIn,nyIn+2)
                0328         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0329      &                      SQUEEZE_RIGHT, myThid )
2c18da0d62 Jean*0330         _END_MASTER(myThid)
1bdb364c26 Jean*0331       ENDIF
                0332 #endif /* ALLOW_DEBUG */
57bf62b584 Dimi*0333
1bdb364c26 Jean*0334 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0335 C---  Prepare output grid and interpolate for each tile
                0336
                0337 C--   put xG in interval [ lon_0 , lon_0+360 [
                0338       DO bj=myByLo(myThid),myByHi(myThid)
                0339        DO bi=myBxLo(myThid),myBxHi(myThid)
                0340         DO j=1-OLy,sNy+OLy
                0341          DO i=1-OLx,sNx+OLx
                0342           xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
                0343      &                  + threeSixtyRS*2.
                0344           xG(i,j,bi,bj) = lon_0+MOD(xG(i,j,bi,bj),threeSixtyRS)
                0345          ENDDO
                0346         ENDDO
77d64cb6a9 Dimi*0347 #ifdef ALLOW_DEBUG
1bdb364c26 Jean*0348 C--   Check validity of input/output coordinates
                0349         IF ( debugFlag ) THEN
ef1828fb24 Jean*0350          DO j=1,sNy
                0351           DO i=1,sNx
b374dc700b Mart*0352            IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
                0353      &          xG(i,j,bi,bj) .GE. x_in(nxIn+1)   .OR.
                0354      &          yG(i,j,bi,bj) .LT. y_in(0)        .OR.
                0355      &          yG(i,j,bi,bj) .GE. y_in(nyIn+1) ) THEN
ef1828fb24 Jean*0356             l = ILNBLNK(inFile)
                0357             WRITE(msgBuf,'(3A,I6)')
                0358      &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
                0359             CALL PRINT_ERROR( msgBuf, myThid )
                0360             WRITE(msgBuf,'(A)')
                0361      &        'EXF_INTERP: input grid must encompass output grid.'
                0362             CALL PRINT_ERROR( msgBuf, myThid )
                0363             WRITE(msgBuf,'(A,2I8,2I6,A,1P2E14.6)') 'i,j,bi,bj=',
                0364      &        i,j,bi,bj, ' , xG,yG=', xG(i,j,bi,bj), yG(i,j,bi,bj)
                0365             CALL PRINT_ERROR( msgBuf, myThid )
                0366             WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nxIn=', nxIn,
                0367      &        ' , x_in(0,nxIn+1)=', x_in(0) ,x_in(nxIn+1)
                0368             CALL PRINT_ERROR( msgBuf, myThid )
                0369             WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nyIn=', nyIn,
                0370      &        ' , y_in(0,nyIn+1)=', y_in(0) ,y_in(nyIn+1)
                0371             CALL PRINT_ERROR( msgBuf, myThid )
                0372             STOP 'ABNORMAL END: S/R EXF_INTERP'
                0373            ENDIF
                0374           ENDDO
                0375          ENDDO
                0376         ENDIF
77d64cb6a9 Dimi*0377 #endif /* ALLOW_DEBUG */
1bdb364c26 Jean*0378        ENDDO
                0379       ENDDO
bda7c25f4a Dimi*0380
1bdb364c26 Jean*0381       DO bj = myByLo(myThid), myByHi(myThid)
                0382        DO bi = myBxLo(myThid), myBxHi(myThid)
                0383
b374dc700b Mart*0384 C--   Compute interpolation lon & lat index mapping
                0385 C--     latitude index
                0386         DO j=1,sNy
                0387          DO i=1,sNx
                0388           s_ind(i,j) = 0
                0389           w_ind(i,j) = nyIn+1
ef1828fb24 Jean*0390          ENDDO
b374dc700b Mart*0391         ENDDO
                0392 C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
                0393 C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
                0394         tmpVar = nyIn + 1. _d -3
                0395         nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
                0396         DO l=1,nLoop
f1d225df6a Mart*0397          DO j=1,sNy
                0398           DO i=1,sNx
b374dc700b Mart*0399            IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
                0400             k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
                0401             IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
f1d225df6a Mart*0402               w_ind(i,j) = k
b374dc700b Mart*0403             ELSE
                0404               s_ind(i,j) = k
f1d225df6a Mart*0405             ENDIF
b374dc700b Mart*0406            ENDIF
f1d225df6a Mart*0407           ENDDO
                0408          ENDDO
b374dc700b Mart*0409         ENDDO
ef1828fb24 Jean*0410 #ifdef ALLOW_DEBUG
1bdb364c26 Jean*0411         IF ( debugFlag ) THEN
ef1828fb24 Jean*0412 C-      Check that we found the right lat. index
                0413          DO j=1,sNy
                0414           DO i=1,sNx
b374dc700b Mart*0415            IF ( w_ind(i,j).NE.s_ind(i,j)+1 ) THEN
ef1828fb24 Jean*0416             l = ILNBLNK(inFile)
                0417             WRITE(msgBuf,'(3A,I4,A,I4)')
                0418      &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,
                0419      &        ', nLoop=', nLoop
                0420             CALL PRINT_ERROR( msgBuf, myThid )
                0421             WRITE(msgBuf,'(A)')
f1d225df6a Mart*0422      &        'EXF_INTERP: did not find Latitude index for grid-pt:'
ef1828fb24 Jean*0423             CALL PRINT_ERROR( msgBuf, myThid )
                0424             WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
                0425      &        'EXF_INTERP: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
                0426             CALL PRINT_ERROR( msgBuf, myThid )
                0427             WRITE(msgBuf,'(A,I8,A,1PE16.8)')
                0428      &        'EXF_INTERP: s_ind=',s_ind(i,j),', lat=',y_in(s_ind(i,j))
                0429             CALL PRINT_ERROR( msgBuf, myThid )
                0430             WRITE(msgBuf,'(A,I8,A,1PE16.8)')
                0431      &        'EXF_INTERP: n_ind=',w_ind(i,j),', lat=',y_in(w_ind(i,j))
                0432             CALL PRINT_ERROR( msgBuf, myThid )
                0433             STOP 'ABNORMAL END: S/R EXF_INTERP'
                0434            ENDIF
                0435           ENDDO
                0436          ENDDO
                0437         ENDIF
                0438 #endif /* ALLOW_DEBUG */
                0439 C--     longitude index
                0440         DO j=1,sNy
                0441          DO i=1,sNx
                0442            w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
                0443          ENDDO
                0444         ENDDO
981f32514d Dimi*0445
1bdb364c26 Jean*0446 C--   Do interpolation using lon & lat index mapping
                0447         CALL EXF_INTERPOLATE(
                0448      I                inFile, irecord, method,
                0449      I                nxIn, nyIn, x_in, y_in,
c2d594df43 Jean*0450      I                arrayin,
1bdb364c26 Jean*0451      O                arrayout,
                0452      I                xG, yG,
                0453      I                w_ind, s_ind,
                0454      I                bi, bj, myThid )
981f32514d Dimi*0455
ef1828fb24 Jean*0456        ENDDO
                0457       ENDDO
981f32514d Dimi*0458
3e6ce1cd50 Jean*0459       RETURN
981f32514d Dimi*0460       END