Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:39:55 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!