pkg/exf/exf_interpolate.F

bc51b4ef50 Jean*0001 #include "EXF_OPTIONS.h"
                0002
                0003 C--  File exf_interp.F: Routines to interpolate input field on to model grid
                0004 C--   Contents
                0005 C--   o S/R EXF_INTERPOLATE
                0006 C--   o FCT LAGRAN
                0007
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0009
                0010 CBOP
                0011 C     !ROUTINE: LAGRAN
                0012 C     !INTERFACE:
                0013        _RL FUNCTION LAGRAN(i,x,a,sp)
                0014
                0015 C !DESCRIPTION: \bv
                0016 C  *==========================================================*
                0017 C  | FUNCTION LAGRAN
                0018 C  | o Provide linear (sp=2) and cubic (sp=4) interpolation
                0019 C  |   weight as lagrange polynomial coefficient.
                0020 C  *==========================================================*
                0021 C \ev
                0022
                0023 C !USES:
                0024        IMPLICIT NONE
                0025
                0026 C !INPUT/OUTPUT PARAMETERS:
                0027         INTEGER i
                0028         _RS x
                0029         _RL a(4)
                0030         INTEGER sp
                0031
                0032 C !LOCAL VARIABLES:
                0033         INTEGER k
                0034         _RL numer,denom
                0035
                0036         numer = 1. _d 0
                0037         denom = 1. _d 0
                0038
                0039 #ifdef TARGET_NEC_SX
                0040 !CDIR UNROLL=8
                0041 #endif /* TARGET_NEC_SX */
                0042         DO k=1,sp
                0043          IF ( k .NE. i) THEN
                0044           denom = denom*(a(i) - a(k))
                0045           numer = numer*(x    - a(k))
                0046          ENDIF
                0047         ENDDO
                0048
                0049         LAGRAN = numer/denom
                0050
                0051        RETURN
                0052        END
                0053
                0054 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0055
                0056 CBOP
                0057 C     !ROUTINE: EXF_INTERPOLATE
                0058 C     !INTERFACE:
                0059        SUBROUTINE EXF_INTERPOLATE(
                0060      I                inFile, irecord, method,
                0061      I                nxIn, nyIn, x_in, y_in,
                0062      I                arrayin,
                0063      O                arrayout,
                0064      I                xG, yG,
                0065      I                w_ind, s_ind,
                0066      I                bi, bj, myThid )
                0067
                0068 C !DESCRIPTION: \bv
                0069 C  *==========================================================*
                0070 C  | SUBROUTINE EXF_INTERPOLATE
                0071 C  | o Interpolate a regular lat-lon input field
                0072 C  |   on to the model grid location
                0073 C  *==========================================================*
                0074 C \ev
                0075
                0076 C !USES:
                0077       IMPLICIT NONE
                0078 C     === Global variables ===
                0079 #include "SIZE.h"
                0080 #include "EEPARAMS.h"
                0081 #include "PARAMS.h"
                0082
                0083 C !INPUT/OUTPUT PARAMETERS:
                0084 C   inFile     :: name of the binary input file (direct access)
                0085 C   irecord    :: record number in input file
                0086 C   method     :: 1,11,21 for bilinear; 2,12,22 for bicubic
                0087 C              :: 1,2 for tracer; 11,12 for U; 21,22 for V
                0088 C   nxIn,nyIn  :: size in x & y direction of input field
                0089 C    x_in      :: longitude vector defining input field grid
                0090 C    y_in      :: latitude  vector defining input field grid
                0091 C   arrayin    :: input field array (loaded from file)
                0092 C   arrayout   :: output array
                0093 C   xG, yG     :: coordinates for output grid to interpolate to
                0094 C    w_ind     :: input field longitudinal index, on western side of model grid pt
                0095 C    s_ind     :: input field latitudinal index, on southern side of model grid pt
                0096 C   bi, bj     :: tile indices
                0097 C   myThid     :: My Thread Id number
                0098
                0099       CHARACTER*(*) inFile
                0100       INTEGER       irecord, method
                0101       INTEGER       nxIn, nyIn
                0102       _RL           x_in(-1:nxIn+2), y_in(-1:nyIn+2)
                0103       _RL           arrayin( -1:nxIn+2, -1:nyIn+2 )
                0104       _RL           arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0105       _RS           xG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0106       _RS           yG      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0107       INTEGER       w_ind(sNx,sNy), s_ind(sNx,sNy)
                0108       INTEGER       bi, bj, myThid
                0109
                0110 C !FUNCTIONS:
                0111       EXTERNAL LAGRAN
                0112       _RL      LAGRAN
                0113       INTEGER  ILNBLNK
                0114       EXTERNAL ILNBLNK
                0115
                0116 C !LOCAL VARIABLES:
                0117 C   px_ind     :: local copy of longitude position around current model grid pt
                0118 C   py_ind     :: local copy of latitude  position around current model grid pt
                0119 C   ew_val     :: intermediate field value after interpolation in East-West dir.
                0120 C   sp         :: number of input-field values used in 1.d interpolation
                0121 C   i, j, k, l :: loop indices
                0122 C   msgBuf     :: Informational/error message buffer
                0123       _RL      px_ind(4), py_ind(4), ew_val(4)
                0124       INTEGER  sp
                0125       INTEGER  i, j, k, l
                0126       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0127 #ifdef TARGET_NEC_SX
                0128       _RL      ew_val1, ew_val2, ew_val3, ew_val4
                0129 #endif
                0130 CEOP
                0131
                0132       IF (method.EQ.1 .OR. method.EQ.11 .OR. method.EQ.21) THEN
                0133
                0134 C--   Bilinear interpolation
                0135          sp = 2
                0136          DO j=1,sNy
                0137           DO i=1,sNx
                0138            arrayout(i,j,bi,bj) = 0.
                0139            DO l=0,1
                0140             px_ind(l+1) = x_in(w_ind(i,j)+l)
                0141             py_ind(l+1) = y_in(s_ind(i,j)+l)
                0142            ENDDO
                0143 #ifndef TARGET_NEC_SX
                0144            DO k=1,2
                0145             ew_val(k) = arrayin(w_ind(i,j)  ,s_ind(i,j)+k-1)
                0146      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0147      &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-1)
                0148      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0149             arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
                0150      &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)
                0151            ENDDO
                0152 #else
                0153            ew_val1 = arrayin(w_ind(i,j)  ,s_ind(i,j)  )
                0154      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0155      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
                0156      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0157            ew_val2 = arrayin(w_ind(i,j)  ,s_ind(i,j)+1)
                0158      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0159      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
                0160      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0161            arrayout(i,j,bi,bj)=
                0162      &            +ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
                0163      &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
                0164 #endif /* TARGET_NEC_SX defined */
                0165           ENDDO
                0166          ENDDO
                0167       ELSEIF (method .EQ. 2 .OR. method.EQ.12 .OR. method.EQ.22) THEN
                0168
                0169 C--   Bicubic interpolation
                0170          sp = 4
                0171          DO j=1,sNy
                0172           DO i=1,sNx
                0173            arrayout(i,j,bi,bj) = 0.
                0174            DO l=-1,2
                0175             px_ind(l+2) = x_in(w_ind(i,j)+l)
                0176             py_ind(l+2) = y_in(s_ind(i,j)+l)
                0177            ENDDO
                0178 #ifndef TARGET_NEC_SX
                0179            DO k=1,4
                0180             ew_val(k) = arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2)
                0181      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0182      &                + arrayin(w_ind(i,j)  ,s_ind(i,j)+k-2)
                0183      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0184      &                + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-2)
                0185      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
                0186      &                + arrayin(w_ind(i,j)+2,s_ind(i,j)+k-2)
                0187      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
                0188             arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj)
                0189      &         + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp)
                0190            ENDDO
                0191 #else
                0192            ew_val1 = arrayin(w_ind(i,j)-1,s_ind(i,j)-1)
                0193      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0194      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)-1)
                0195      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0196      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)-1)
                0197      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
                0198      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)-1)
                0199      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
                0200            ew_val2 = arrayin(w_ind(i,j)-1,s_ind(i,j)  )
                0201      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0202      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)  )
                0203      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0204      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)  )
                0205      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
                0206      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)  )
                0207      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
                0208            ew_val3 = arrayin(w_ind(i,j)-1,s_ind(i,j)+1)
                0209      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0210      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+1)
                0211      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0212      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+1)
                0213      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
                0214      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+1)
                0215      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
                0216            ew_val4 = arrayin(w_ind(i,j)-1,s_ind(i,j)+2)
                0217      &                    *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp)
                0218      &             + arrayin(w_ind(i,j)  ,s_ind(i,j)+2)
                0219      &                    *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp)
                0220      &             + arrayin(w_ind(i,j)+1,s_ind(i,j)+2)
                0221      &                    *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp)
                0222      &             + arrayin(w_ind(i,j)+2,s_ind(i,j)+2)
                0223      &                    *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp)
                0224            arrayout(i,j,bi,bj) =
                0225      &             ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp)
                0226      &            +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp)
                0227      &            +ew_val3*LAGRAN(3,yG(i,j,bi,bj),py_ind,sp)
                0228      &            +ew_val4*LAGRAN(4,yG(i,j,bi,bj),py_ind,sp)
                0229 #endif /* TARGET_NEC_SX defined */
                0230           ENDDO
                0231          ENDDO
                0232       ELSE
                0233          l = ILNBLNK(inFile)
                0234          WRITE(msgBuf,'(3A,I6)')
                0235      &    'EXF_INTERPOLATE: file="', inFile(1:l), '", rec=', irecord
                0236          CALL PRINT_ERROR( msgBuf, myThid )
                0237          WRITE(msgBuf,'(A,I8,A)')
                0238      &    'EXF_INTERPOLATE: method=', method,' not supported'
                0239          CALL PRINT_ERROR( msgBuf, myThid )
                0240          STOP 'ABNORMAL END: S/R EXF_INTERPOLATE: invalid method'
                0241       ENDIF
                0242
                0243       RETURN
                0244       END