#include "EXF_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: EXF_INTERP_UV
C     !INTERFACE:
       SUBROUTINE EXF_INTERP_UV(
     I                inFileU, inFileV, filePrec, irecord,
     I                nxIn, nyIn,
     I                lon_0, lon_inc, lat_0, lat_inc,
#ifdef EXF_INTERP_USE_DYNALLOC
     O                arrUout, arrVout,
#else
     O                arrUout, arrVout, arrUin, arrVin,
#endif
     I                xG_in, yG,
     I                methodU, methodV, myIter, myThid )

C !DESCRIPTION: \bv
C  *==========================================================*
C  | SUBROUTINE EXF_INTERP_UV
C  | o Load from file the 2 vector components of regular
C  |   lat-lon input field and interpolate on to the model
C  |   grid location
C  *==========================================================*
C \ev

C !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EXF_INTERP_SIZE.h"
#ifdef ALLOW_DEBUG
# include "EXF_PARAM.h"
#endif

C !INPUT/OUTPUT PARAMETERS:
C  inFileU     (string)  :: input file name for the 1rst component (U)
C  inFileV     (string)  :: input file name for the 2nd  component (V)
C  filePrec    (integer) :: number of bits per word in file (32 or 64)
C  irecord     (integer) :: record number to read
C  nxIn, nyIn  (integer) :: size in x & y direction of input file to read
C   lon_0, lat_0         :: lon and lat of sw corner of global input grid
C   lon_inc              :: scalar x-grid increment
C   lat_inc              :: vector y-grid increments
C  arrUout     ( _RL )   :: 1rst component (U) output array
C  arrVout     ( _RL )   :: 2nd  component (V) output array
#ifndef EXF_INTERP_USE_DYNALLOC
C  arrUin      ( _RL )   :: 1rst component input field array (loaded from file)
C  arrVin      ( _RL )   :: 2nd  component input field array (loaded from file)
#endif
C   xG_in,  yG           :: coordinates for output grid to interpolate to
C  methodU, methodV      :: 1,11,21 for bilinear; 2,12,22 for bicubic
C                        :: 1,2 for tracer; 11,12 for U; 21,22 for V
C  myIter      (integer) :: current iteration number
C  myThid      (integer) :: My Thread Id number
      CHARACTER*(*) inFileU
      CHARACTER*(*) inFileV
      INTEGER       filePrec, irecord, nxIn, nyIn
      _RL           arrUout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL           arrVout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifndef EXF_INTERP_USE_DYNALLOC
      _RL           arrUin ( -1:nxIn+2, -1:nyIn+2 )
      _RL           arrVin ( -1:nxIn+2, -1:nyIn+2 )
#endif
      _RS           xG_in  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS           yG     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL           lon_0, lon_inc
c     _RL           lat_0, lat_inc(nyIn-1)
      _RL           lat_0, lat_inc(*)
      INTEGER       methodU, methodV, myIter, myThid

C !FUNCTIONS:
#ifdef ALLOW_DEBUG
      INTEGER  ILNBLNK
      EXTERNAL ILNBLNK
#endif

C !LOCAL VARIABLES:
C     x_in        :: longitude vector defining input field grid
C     y_in        :: latitude  vector defining input field grid
C     w_ind       :: input field longitudinal index, on western side of model grid pt
C     s_ind       :: input field latitudinal index, on southern side of model grid pt
C     bi, bj      :: tile indices
C     i, j, k, l  :: loop indices
C     msgBuf      :: Informational/error message buffer
#ifdef EXF_INTERP_USE_DYNALLOC
C     arrUin      :: 1rst component input field array (loaded from file)
C     arrVin      :: 2nd  component input field array (loaded from file)
      _RL      arrUin( -1:nxIn+2, -1:nyIn+2 )
      _RL      arrVin( -1:nxIn+2, -1:nyIn+2 )
      _RL      csLon(-1:nxIn+2), snLon(-1:nxIn+2)
      _RL      x_in (-1:nxIn+2), y_in (-1:nyIn+2)
#else /* EXF_INTERP_USE_DYNALLOC */
      _RL      csLon(-1:exf_max_nLon+2), snLon(-1:exf_max_nLon+2)
      _RL      x_in (-1:exf_max_nLon+2), y_in (-1:exf_max_nLat+2)
#endif /* EXF_INTERP_USE_DYNALLOC */
      INTEGER  w_ind(sNx,sNy), s_ind(sNx,sNy)
      INTEGER  bi, bj
      INTEGER  i, j, k, l
      INTEGER  nLoop
      _RL      tmpVar
      _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS      threeSixtyRS
      _RL      yPole, symSign, poleU, poleV
      _RL      csdLon, sndLon, pSign
      _RL      edgeFac, poleFac
      LOGICAL  calcLonCS
      PARAMETER ( threeSixtyRS = 360. )
      PARAMETER ( yPole = 90. )
      INTEGER  nxd2
      LOGICAL  xIsPeriodic, poleSymmetry
#ifdef ALLOW_DEBUG
      LOGICAL  debugFlag
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RS      prtPole(-1:4)
#endif
CEOP

#ifndef EXF_INTERP_USE_DYNALLOC
C--   Check size declaration:
      IF ( nxIn.GT.exf_max_nLon ) THEN
       STOP 'EXF_INTERP_UV: exf_max_nLon too small'
      ENDIF
      IF ( nyIn.GT.exf_max_nLat ) THEN
       STOP 'EXF_INTERP_UV: exf_max_nLat too small'
      ENDIF
      IF ( (nxIn+4)*(nyIn+4).GT.exf_interp_bufferSize ) THEN
       STOP 'EXF_INTERP_UV: exf_interp_bufferSize too small'
      ENDIF
#endif /* ndef EXF_INTERP_USE_DYNALLOC */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C---  Load inut field

      CALL EXF_INTERP_READ(
     I         inFileU, filePrec,
     O         arrUin,
     I         irecord, nxIn, nyIn, myThid )
      CALL EXF_INTERP_READ(
     I         inFileV, filePrec,
     O         arrVin,
     I         irecord, nxIn, nyIn, myThid )

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C---  Prepare input grid and input field

C--   setup input longitude grid
      DO i=-1,nxIn+2
       x_in(i) = lon_0 + (i-1)*lon_inc
      ENDDO
      xIsPeriodic = nxIn.EQ.NINT( threeSixtyRS / lon_inc )
      nxd2 = NINT( nxIn*0.5 )
      poleSymmetry = xIsPeriodic .AND. ( nxIn.EQ.2*nxd2 )

C--   setup input latitude grid
      y_in(1) = lat_0
      DO j=1,nyIn+1
       i = MIN(j,nyIn-1)
       y_in(j+1) = y_in(j) + lat_inc(i)
      ENDDO
      y_in(0) = y_in(1) - lat_inc(1)
      y_in(-1)= y_in(0) - lat_inc(1)
#ifdef ALLOW_DEBUG
      DO l=-1,4
        prtPole(l) = 0.
      ENDDO
#endif
C--   Add 1 row @ the pole (if last row is not) and will fill it
C     with the polarmost-row zonally averaged value.

C-    Add 2 row @ southern end; if one is beyond S.pole, put one @ S.pole
      j = 0
      IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
         y_in(j) = -yPole
         y_in(j-1) = -2.*yPole - y_in(j+1)
#ifdef ALLOW_DEBUG
         prtPole(j)   = 1.
         prtPole(j-1) = 2.
#endif
      ENDIF
      j = -1
      IF ( ABS(y_in(j+1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
         y_in(j) = -yPole
#ifdef ALLOW_DEBUG
         prtPole(j)   = 1.
#endif
      ENDIF
C-    Add 2 row @ northern end; if one is beyond N.pole, put one @ N.pole
      j = nyIn+1
      IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
         y_in(j) = yPole
         y_in(j+1) = 2.*yPole - y_in(j-1)
#ifdef ALLOW_DEBUG
         prtPole(3)   = 1.
         prtPole(3+1) = 2.
#endif
      ENDIF
      j = nyIn+2
      IF ( ABS(y_in(j-1)).LT.yPole .AND. ABS(y_in(j)).GT.yPole ) THEN
         y_in(j) = yPole
#ifdef ALLOW_DEBUG
         prtPole(4)   = 1.
#endif
      ENDIF

C--   Enlarge boundary
      IF ( xIsPeriodic ) THEN
C-    fill-in added column assuming periodicity
        DO j=1,nyIn
         arrUin( 0,j)     = arrUin(nxIn  ,j)
         arrUin(-1,j)     = arrUin(nxIn-1,j)
         arrUin(nxIn+1,j) = arrUin(1,j)
         arrUin(nxIn+2,j) = arrUin(2,j)
         arrVin( 0,j)     = arrVin(nxIn  ,j)
         arrVin(-1,j)     = arrVin(nxIn-1,j)
         arrVin(nxIn+1,j) = arrVin(1,j)
         arrVin(nxIn+2,j) = arrVin(2,j)
        ENDDO
      ELSE
C-    fill-in added column from nearest column
        DO j=1,nyIn
         arrUin( 0,j)     = arrUin(1,j)
         arrUin(-1,j)     = arrUin(1,j)
         arrUin(nxIn+1,j) = arrUin(nxIn,j)
         arrUin(nxIn+2,j) = arrUin(nxIn,j)
         arrVin( 0,j)     = arrVin(1,j)
         arrVin(-1,j)     = arrVin(1,j)
         arrVin(nxIn+1,j) = arrVin(nxIn,j)
         arrVin(nxIn+2,j) = arrVin(nxIn,j)
        ENDDO
      ENDIF
      symSign = -1. _d 0
      DO l=-1,2
       j = l
       IF ( l.GE.1 ) j = nyIn+l
       k = MAX(1,MIN(j,nyIn))
       IF ( poleSymmetry .AND. ABS(y_in(j)).GT.yPole ) THEN
        IF ( nyIn.GE.3 .AND. ABS(y_in(k)).EQ.yPole )
     &    k = MAX(2,MIN(j,nyIn-1))
C-    fill-in added row assuming pole-symmetry
        DO i=-1,nxd2
         arrUin(i,j) = symSign*arrUin(i+nxd2,k)
         arrVin(i,j) = symSign*arrVin(i+nxd2,k)
        ENDDO
        DO i=1,nxd2+2
         arrUin(i+nxd2,j) = symSign*arrUin(i,k)
         arrVin(i+nxd2,j) = symSign*arrVin(i,k)
        ENDDO
#ifdef ALLOW_DEBUG
        i = l + 2*( (l+1)/2 )
        prtPole(i) = prtPole(i) + 0.2
#endif
       ELSE
C-    fill-in added row from nearest column values
        DO i=-1,nxIn+2
         arrUin(i,j) = arrUin(i,k)
         arrVin(i,j) = arrVin(i,k)
        ENDDO
       ENDIF
      ENDDO

C     For vector, replace value at the pole with northernmost/southermost
C     zonal-mean value (poleU,poleV corresponds to Lon=0.E orientation).
      IF ( methodU.GE.10 .AND. methodV.GE.10 ) THEN
       calcLonCS = .TRUE.
       DO l=-1,4
        j = l
        IF ( l.GE.2 ) j = nyIn+l-2
        IF ( ABS(y_in(j)).EQ.yPole ) THEN
          pSign = SIGN( oneRL, y_in(j) )
          IF ( calcLonCS ) THEN
            csLon(-1) = cos(x_in(-1)*deg2rad)
            snLon(-1) = sin(x_in(-1)*deg2rad)
            csdLon = cos(lon_inc*deg2rad)
            sndLon = sin(lon_inc*deg2rad)
            DO i=-1,nxIn+1
              csLon(i+1) = csLon(i)*csdLon - snLon(i)*sndLon
              snLon(i+1) = csLon(i)*sndLon + snLon(i)*csdLon
            ENDDO
            calcLonCS = .FALSE.
c           write(0,'(a,1P3E13.5)') 'cs 1,nxIn+1,diff=',
c    &           csLon(1),csLon(nxIn+1),csLon(1)-csLon(nxIn+1)
c           write(0,'(a,1P3E13.5)') 'sn 1,nxIn+1,diff=',
c    &           snLon(1),snLon(nxIn+1),snLon(1)-snLon(nxIn+1)
          ENDIF
C     account for local orientation when averaging
          poleU = 0.
          poleV = 0.
          DO i=1,nxIn
            poleU = poleU
     &       + ( csLon(i)*arrUin(i,j) - pSign*snLon(i)*arrVin(i,j) )
            poleV = poleV
     &       + ( pSign*snLon(i)*arrUin(i,j) + csLon(i)*arrVin(i,j) )
          ENDDO
          poleU = poleU / nxIn
          poleV = poleV / nxIn
C     put back zonal-mean value but locally orientated
          DO i=-1,nxIn+2
           arrUin(i,j) =  csLon(i)*poleU + pSign*snLon(i)*poleV
           arrVin(i,j) = -pSign*snLon(i)*poleU + csLon(i)*poleV
          ENDDO
#ifdef ALLOW_DEBUG
          prtPole(l) = prtPole(l) + 0.1
#endif
        ENDIF
       ENDDO
C-    change first additional row from simple copy to linear interpolation
C     between nearest column values and pole value
       DO l=0,1
        k = l*(nyIn+3) -1
        IF ( ABS(y_in(k)).EQ.yPole ) THEN
         j = l*(nyIn+1)
         i = l*(nyIn-1) +1
         edgeFac = (y_in(j) - y_in(k)) / (y_in(i) - y_in(k))
         poleFac = (y_in(i) - y_in(j)) / (y_in(i) - y_in(k))
         DO i=-1,nxIn+2
           arrUin(i,j) = arrUin(i,j) * edgeFac
     &                 + arrUin(i,k) * poleFac
           arrVin(i,j) = arrVin(i,j) * edgeFac
     &                 + arrVin(i,k) * poleFac
         ENDDO
#ifdef ALLOW_DEBUG
         prtPole(3*l) = prtPole(3*l) + 0.3
#endif
        ENDIF
       ENDDO
      ENDIF

#ifdef ALLOW_DEBUG
      debugFlag = ( exf_debugLev.GE.debLevC )
     &       .OR. ( exf_debugLev.GE.debLevB .AND. myIter.LE.nIter0 )
C     prtPole(l)=0 : extended, =1 : changed to pole, =2 : changed to symetric
      IF ( debugFlag ) THEN
        l = ILNBLNK(inFileU)
        WRITE(msgBuf,'(3A,I6,A,2L5)')
     &   ' EXF_INTERP_UV: fileU="',inFileU(1:l),'", rec=', irecord,
     &   ' , x-Per,P.Sym=', xIsPeriodic, poleSymmetry
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' S.edge (j=-1,0,1) :',
     &   ' proc=', (prtPole(j),j=-1,1), ', yIn=', (y_in(j),j=-1,1)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(2A,3F4.1,A,3F12.6)') ' N.edge (j=+0,+1,+2)',
     &   ' proc=', (prtPole(j),j=2,4), ', yIn=',(y_in(j),j=nyIn,nyIn+2)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
c      IF ( exf_output_interp ) THEN
c       j = (nxIn+4)*(nyIn+4)
c       CALL WRITE_GLVEC_RL( inFileU,'_in', arrUin, j, myIter, myThid )
c       CALL WRITE_GLVEC_RL( inFileV,'_in', arrVin, j, myIter, myThid )
c      ENDIF
      ENDIF
#endif /* ALLOW_DEBUG */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C---  Prepare output grid and interpolate for each tile

C--   put xG in interval [ lon_0 , lon_0+360 [
      DO bj=myByLo(myThid),myByHi(myThid)
       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
#ifdef ALLOW_DEBUG
C--   Check validity of input/output coordinates
        IF ( debugFlag ) THEN
         DO j=1,sNy
          DO i=1,sNx
           IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
     &          xG(i,j,bi,bj) .GE. x_in(nxIn+1)   .OR.
     &          yG(i,j,bi,bj) .LT. y_in(0)        .OR.
     &          yG(i,j,bi,bj) .GE. y_in(nyIn+1) ) THEN
            l = ILNBLNK(inFileU)
            WRITE(msgBuf,'(3A,I6)')
     &        'EXF_INTERP_UV: fileU="',inFileU(1:l),'", rec=',irecord
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A)')
     &        'EXF_INTERP_UV: input grid must encompass output grid.'
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,2I8,2I6,A,1P2E14.6)') 'i,j,bi,bj=',
     &        i,j,bi,bj, ' , xG,yG=', xG(i,j,bi,bj), yG(i,j,bi,bj)
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nxIn=', nxIn,
     &        ' , x_in(0,nxIn+1)=', x_in(0) ,x_in(nxIn+1)
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,I9,A,1P2E14.6)') 'nyIn=', nyIn,
     &        ' , y_in(0,nyIn+1)=', y_in(0) ,y_in(nyIn+1)
            CALL PRINT_ERROR( msgBuf, myThid )
            STOP 'ABNORMAL END: S/R EXF_INTERP_UV'
           ENDIF
          ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_DEBUG */
       ENDDO
      ENDDO

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C--   Compute interpolation lon & lat index mapping
C--     latitude index
        DO j=1,sNy
         DO i=1,sNx
          s_ind(i,j) = 0
          w_ind(i,j) = nyIn+1
         ENDDO
        ENDDO
C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
        tmpVar = nyIn + 1. _d -3
        nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
        DO l=1,nLoop
         DO j=1,sNy
          DO i=1,sNx
           IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
            k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
            IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
              w_ind(i,j) = k
            ELSE
              s_ind(i,j) = k
            ENDIF
           ENDIF
          ENDDO
         ENDDO
        ENDDO
#ifdef ALLOW_DEBUG
        IF ( debugFlag ) THEN
C-      Check that we found the right lat. index
         DO j=1,sNy
          DO i=1,sNx
           IF ( w_ind(i,j).NE.s_ind(i,j)+1 ) THEN
            l = ILNBLNK(inFileU)
            WRITE(msgBuf,'(3A,I4,A,I4)')
     &      'EXF_INTERP_UV: file="', inFileU(1:l), '", rec=', irecord,
     &      ', nLoop=', nLoop
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A)')
     &      'EXF_INTERP_UV: did not find Latitude index for grid-pt:'
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,2I8,2I6,A,1PE16.8)')
     &      'EXF_INTERP_UV: i,j,bi,bj=',i,j,bi,bj,' , yG=',yG(i,j,bi,bj)
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,I8,A,1PE16.8)')
     &      'EXF_INTERP_UV: s_ind=',s_ind(i,j),', lat=',y_in(s_ind(i,j))
            CALL PRINT_ERROR( msgBuf, myThid )
            WRITE(msgBuf,'(A,I8,A,1PE16.8)')
     &      'EXF_INTERP_UV: n_ind=',w_ind(i,j),', lat=',y_in(w_ind(i,j))
            CALL PRINT_ERROR( msgBuf, myThid )
            STOP 'ABNORMAL END: S/R EXF_INTERP_UV'
           ENDIF
          ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_DEBUG */
C--     longitude index
        DO j=1,sNy
         DO i=1,sNx
           w_ind(i,j) = INT((xG(i,j,bi,bj)-x_in(-1))/lon_inc) - 1
         ENDDO
        ENDDO

C--   Do interpolation using lon & lat index mapping
        CALL EXF_INTERPOLATE(
     I                inFileU, irecord, methodU,
     I                nxIn, nyIn, x_in, y_in,
     I                arrUin,
     O                arrUout,
     I                xG, yG,
     I                w_ind, s_ind,
     I                bi, bj, myThid )
        CALL EXF_INTERPOLATE(
     I                inFileV, irecord, methodV,
     I                nxIn, nyIn, x_in, y_in,
     I                arrVin,
     O                arrVout,
     I                xG, yG,
     I                w_ind, s_ind,
     I                bi, bj, myThid )

       ENDDO
      ENDDO

      RETURN
      END
