/[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.32 by mlosch, Fri Mar 13 14:16:36 2015 UTC revision 1.33 by mlosch, Thu Apr 16 16:13:30 2015 UTC
# Line 86  C     msgBuf      :: Informational/error Line 86  C     msgBuf      :: Informational/error
86        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
87        PARAMETER ( yPole = 90. )        PARAMETER ( yPole = 90. )
88        INTEGER  nxd2        INTEGER  nxd2
89        LOGICAL  xIsPeriodic, poleSymmetry, latFromS2N        LOGICAL  xIsPeriodic, poleSymmetry
90  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
91        LOGICAL  debugFlag, gridNotRight        LOGICAL  debugFlag
92        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
93        _RS      prtPole(-1:4)        _RS      prtPole(-1:4)
94  #endif  #endif
# Line 173  C-    Add 2 row @ northern end; if one i Line 173  C-    Add 2 row @ northern end; if one i
173  #endif  #endif
174         ENDIF         ENDIF
175        ENDIF        ENDIF
 C     store orientation of latitude  
       latFromS2N = y_in(0) .LT. y_in(nyIn+1)  
176    
177  C--   Enlarge boundary  C--   Enlarge boundary
178        IF ( xIsPeriodic ) THEN        IF ( xIsPeriodic ) THEN
# Line 302  C--   Check validity of input/output coo Line 300  C--   Check validity of input/output coo
300          IF ( debugFlag ) THEN          IF ( debugFlag ) THEN
301           DO j=1,sNy           DO j=1,sNy
302            DO i=1,sNx            DO i=1,sNx
303             IF ( latFromS2N ) THEN             IF ( xG(i,j,bi,bj) .LT. x_in(0)        .OR.
304  C     y_in(0) is the southern edge       &          xG(i,j,bi,bj) .GE. x_in(nxIn+1)   .OR.
305              gridNotRight = (       &          yG(i,j,bi,bj) .LT. y_in(0)        .OR.
306       &           xG(i,j,bi,bj) .LT. x_in(0)        .OR.       &          yG(i,j,bi,bj) .GE. y_in(nyIn+1) ) THEN
      &           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) )  
            ELSE  
 C     y_in(0) is the northern edge  
             gridNotRight = (  
      &           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) .GE. y_in(0)        .OR.  
      &           yG(i,j,bi,bj) .LT. y_in(nyIn+1) )  
            ENDIF  
            IF ( gridNotRight ) THEN  
307              l = ILNBLNK(inFile)              l = ILNBLNK(inFile)
308              WRITE(msgBuf,'(3A,I6)')              WRITE(msgBuf,'(3A,I6)')
309       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord
# Line 343  C     y_in(0) is the northern edge Line 329  C     y_in(0) is the northern edge
329         ENDDO         ENDDO
330        ENDDO        ENDDO
331    
 C--   Compute interpolation lon & lat index mapping  
 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) )  
 C  
332        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
333         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
334    
335          IF ( latFromS2N ) THEN  C--   Compute interpolation lon & lat index mapping
336  C--   Case 1: y_in(0) is the southern edge  C--     latitude index
337  C--      latitude index          DO j=1,sNy
338           DO j=1,sNy           DO i=1,sNx
339            DO i=1,sNx            s_ind(i,j) = 0
340             s_ind(i,j) = 0            w_ind(i,j) = nyIn+1
            w_ind(i,j) = nyIn+1  
           ENDDO  
          ENDDO  
          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  
341           ENDDO           ENDDO
342          ELSE          ENDDO
343  C--   Case 2: y_in(0) is the northern edge  C       # of pts = nyIn+2 ; # of interval = nyIn+1 ; evaluate nLoop as
344  C--   In this case the roles of s_ind and w_ind are exchanged: now  C       1 + truncated log2(# interval -1); add epsil=1.e-3 for safey
345  C--   s_ind is the index on the northern side grid box          tmpVar = nyIn + 1. _d -3
346  C--      latitude index          nLoop = 1 + INT( LOG(tmpVar)/LOG(2. _d 0) )
347            DO l=1,nLoop
348           DO j=1,sNy           DO j=1,sNy
349            DO i=1,sNx            DO i=1,sNx
350             s_ind(i,j) = nyIn+1             IF ( w_ind(i,j).GT.s_ind(i,j)+1 ) THEN
351             w_ind(i,j) = 0              k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )
352            ENDDO              IF ( yG(i,j,bi,bj) .LT. y_in(k) ) THEN
          ENDDO  
          DO l=1,nLoop  
           DO j=1,sNy  
            DO i=1,sNx  
             IF ( s_ind(i,j).GT.w_ind(i,j)+1 ) THEN  
              k = NINT( (s_ind(i,j)+w_ind(i,j))*0.5 )  
              IF ( yG(i,j+1,bi,bj) .GT. y_in(k) ) THEN  
               s_ind(i,j) = k  
              ELSE  
353                w_ind(i,j) = k                w_ind(i,j) = k
354               ENDIF              ELSE
355                  s_ind(i,j) = k
356              ENDIF              ENDIF
357             ENDDO             ENDIF
358            ENDDO            ENDDO
359           ENDDO           ENDDO
360          ENDIF          ENDDO
361  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
362          IF ( debugFlag ) 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
366             IF ( latFromS2N ) THEN             IF ( w_ind(i,j).NE.s_ind(i,j)+1 ) THEN
             gridNotRight = w_ind(i,j).NE.s_ind(i,j)+1  
            ELSE  
             gridNotRight = s_ind(i,j).NE.w_ind(i,j)+1  
            ENDIF  
            IF ( gridNotRight ) THEN  
367              l = ILNBLNK(inFile)              l = ILNBLNK(inFile)
368              WRITE(msgBuf,'(3A,I4,A,I4)')              WRITE(msgBuf,'(3A,I4,A,I4)')
369       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,       &        'EXF_INTERP: file="', inFile(1:l), '", rec=', irecord,
Legend:
Removed from v.1.32  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22