#include "EXF_OPTIONS.h" C-- File exf_interp.F: Routines to interpolate input field on to model grid C-- Contents C-- o S/R EXF_INTERPOLATE C-- o FCT LAGRAN C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: LAGRAN C !INTERFACE: _RL FUNCTION LAGRAN(i,x,a,sp) C !DESCRIPTION: \bv C *==========================================================* C | FUNCTION LAGRAN C | o Provide linear (sp=2) and cubic (sp=4) interpolation C | weight as lagrange polynomial coefficient. C *==========================================================* C \ev C !USES: IMPLICIT NONE C !INPUT/OUTPUT PARAMETERS: INTEGER i _RS x _RL a(4) INTEGER sp C !LOCAL VARIABLES: INTEGER k _RL numer,denom numer = 1. _d 0 denom = 1. _d 0 #ifdef TARGET_NEC_SX !CDIR UNROLL=8 #endif /* TARGET_NEC_SX */ DO k=1,sp IF ( k .NE. i) THEN denom = denom*(a(i) - a(k)) numer = numer*(x - a(k)) ENDIF ENDDO LAGRAN = numer/denom RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: EXF_INTERPOLATE C !INTERFACE: SUBROUTINE EXF_INTERPOLATE( I inFile, irecord, method, I nxIn, nyIn, x_in, y_in, I arrayin, O arrayout, I xG, yG, I w_ind, s_ind, I bi, bj, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE EXF_INTERPOLATE C | o Interpolate a regular lat-lon input field C | on to the model grid location C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" C !INPUT/OUTPUT PARAMETERS: C inFile :: name of the binary input file (direct access) C irecord :: record number in input file C method :: 1,11,21 for bilinear; 2,12,22 for bicubic C :: 1,2 for tracer; 11,12 for U; 21,22 for V C nxIn,nyIn :: size in x & y direction of input field C x_in :: longitude vector defining input field grid C y_in :: latitude vector defining input field grid C arrayin :: input field array (loaded from file) C arrayout :: output array C xG, yG :: coordinates for output grid to interpolate to 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 myThid :: My Thread Id number CHARACTER*(*) inFile INTEGER irecord, method INTEGER nxIn, nyIn _RL x_in(-1:nxIn+2), y_in(-1:nyIn+2) _RL arrayin( -1:nxIn+2, -1:nyIn+2 ) _RL arrayout(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS xG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS yG (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) INTEGER w_ind(sNx,sNy), s_ind(sNx,sNy) INTEGER bi, bj, myThid C !FUNCTIONS: EXTERNAL LAGRAN _RL LAGRAN INTEGER ILNBLNK EXTERNAL ILNBLNK C !LOCAL VARIABLES: C px_ind :: local copy of longitude position around current model grid pt C py_ind :: local copy of latitude position around current model grid pt C ew_val :: intermediate field value after interpolation in East-West dir. C sp :: number of input-field values used in 1.d interpolation C i, j, k, l :: loop indices C msgBuf :: Informational/error message buffer _RL px_ind(4), py_ind(4), ew_val(4) INTEGER sp INTEGER i, j, k, l CHARACTER*(MAX_LEN_MBUF) msgBuf #ifdef TARGET_NEC_SX _RL ew_val1, ew_val2, ew_val3, ew_val4 #endif CEOP IF (method.EQ.1 .OR. method.EQ.11 .OR. method.EQ.21) THEN C-- Bilinear interpolation sp = 2 DO j=1,sNy DO i=1,sNx arrayout(i,j,bi,bj) = 0. DO l=0,1 px_ind(l+1) = x_in(w_ind(i,j)+l) py_ind(l+1) = y_in(s_ind(i,j)+l) ENDDO #ifndef TARGET_NEC_SX DO k=1,2 ew_val(k) = arrayin(w_ind(i,j) ,s_ind(i,j)+k-1) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-1) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj) & + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp) ENDDO #else ew_val1 = arrayin(w_ind(i,j) ,s_ind(i,j) ) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j) ) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) ew_val2 = arrayin(w_ind(i,j) ,s_ind(i,j)+1) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)+1) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) arrayout(i,j,bi,bj)= & +ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp) & +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp) #endif /* TARGET_NEC_SX defined */ ENDDO ENDDO ELSEIF (method .EQ. 2 .OR. method.EQ.12 .OR. method.EQ.22) THEN C-- Bicubic interpolation sp = 4 DO j=1,sNy DO i=1,sNx arrayout(i,j,bi,bj) = 0. DO l=-1,2 px_ind(l+2) = x_in(w_ind(i,j)+l) py_ind(l+2) = y_in(s_ind(i,j)+l) ENDDO #ifndef TARGET_NEC_SX DO k=1,4 ew_val(k) = arrayin(w_ind(i,j)-1,s_ind(i,j)+k-2) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j) ,s_ind(i,j)+k-2) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)+k-2) & *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+2,s_ind(i,j)+k-2) & *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp) arrayout(i,j,bi,bj) = arrayout(i,j,bi,bj) & + ew_val(k)*LAGRAN(k,yG(i,j,bi,bj),py_ind,sp) ENDDO #else ew_val1 = arrayin(w_ind(i,j)-1,s_ind(i,j)-1) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j) ,s_ind(i,j)-1) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)-1) & *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+2,s_ind(i,j)-1) & *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp) ew_val2 = arrayin(w_ind(i,j)-1,s_ind(i,j) ) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j) ,s_ind(i,j) ) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j) ) & *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+2,s_ind(i,j) ) & *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp) ew_val3 = arrayin(w_ind(i,j)-1,s_ind(i,j)+1) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j) ,s_ind(i,j)+1) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)+1) & *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+2,s_ind(i,j)+1) & *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp) ew_val4 = arrayin(w_ind(i,j)-1,s_ind(i,j)+2) & *LAGRAN(1,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j) ,s_ind(i,j)+2) & *LAGRAN(2,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+1,s_ind(i,j)+2) & *LAGRAN(3,xG(i,j,bi,bj),px_ind,sp) & + arrayin(w_ind(i,j)+2,s_ind(i,j)+2) & *LAGRAN(4,xG(i,j,bi,bj),px_ind,sp) arrayout(i,j,bi,bj) = & ew_val1*LAGRAN(1,yG(i,j,bi,bj),py_ind,sp) & +ew_val2*LAGRAN(2,yG(i,j,bi,bj),py_ind,sp) & +ew_val3*LAGRAN(3,yG(i,j,bi,bj),py_ind,sp) & +ew_val4*LAGRAN(4,yG(i,j,bi,bj),py_ind,sp) #endif /* TARGET_NEC_SX defined */ ENDDO ENDDO ELSE l = ILNBLNK(inFile) WRITE(msgBuf,'(3A,I6)') & 'EXF_INTERPOLATE: file="', inFile(1:l), '", rec=', irecord CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A,I8,A)') & 'EXF_INTERPOLATE: method=', method,' not supported' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R EXF_INTERPOLATE: invalid method' ENDIF RETURN END