pkg/flt/flt_interp_linear.F

51ec3c32fe Jean*0001 #include "FLT_OPTIONS.h"
                0002
                0003 C--   Contents
                0004 C--   o FLT_BILINEAR
                0005 C--   o FLT_TRILINEAR
                0006 C--   o FLT_BILINEAR2D
                0007
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0009
                0010       SUBROUTINE FLT_BILINEAR(
                0011      I                         ix, jy,
                0012      O                         uu,
                0013      I                         var,
                0014      I                         kl, nu, bi, bj, myThid )
                0015
                0016 C     ==================================================================
                0017 C     SUBROUTINE FLT_BILINEAR
                0018 C     ==================================================================
                0019 C     o Bilinear scheme to interpolate variable to particle position
                0020 C       given by its fractional (real) index ix,jy location
                0021 C     ==================================================================
                0022
                0023 C     !USES:
                0024       IMPLICIT NONE
                0025
                0026 C     == global variables ==
                0027 #include "SIZE.h"
                0028
                0029 C     == routine arguments ==
                0030       _RL ix, jy
                0031       _RL uu
                0032       _RL var(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0033       INTEGER kl, nu, bi, bj, myThid
                0034
                0035 C     == local variables ==
                0036       INTEGER i1, j1, i2, j2, klp
                0037       _RL ddx, ddy
                0038       _RL u11, u12, u22, u21
                0039
                0040 C     == end of interface ==
                0041
                0042 C--   to choose the u box in which the particle is found
                0043 C nu=0 for T, S
                0044 C nu=1 for u
                0045 C nu=2 for v
                0046 C nu=3 for Vorticity
                0047 C nu=4 for w
                0048
                0049       IF ( kl.LT.1 .OR. kl.GT.Nr ) THEN
                0050 c           WRITE(msgbuf,'(A,I8)')
                0051 c    &        ' FLT_BILINEAR: illegal value for kl=',kl
                0052 c           CALL PRINT_ERROR( msgbuf, myThid )
                0053             STOP 'ABNORMAL END: S/R FLT_BILINEAR'
                0054       ENDIF
                0055
                0056 C--   find x-index according to grid-location of variable
                0057       IF ( MOD(nu,2).EQ.0 ) THEN
                0058         i1 = INT(ix)
                0059         ddx = ix - DFLOAT(i1)
                0060       ELSE
                0061         i1 = NINT(ix)
                0062         ddx = 0.5 _d 0 + ix - DFLOAT(i1)
                0063       ENDIF
                0064 C--   find y-index according to grid-location of variable
                0065       IF ( MOD(nu,4).LE.1 ) THEN
                0066         j1 = INT(jy)
                0067         ddy = jy - DFLOAT(j1)
                0068       ELSE
                0069         j1 = NINT(jy)
                0070         ddy = 0.5 _d 0 + jy - DFLOAT(j1)
                0071       ENDIF
                0072
                0073 C--   Set the higher index for interpolation
                0074       i2 = i1 + 1
                0075       j2 = j1 + 1
                0076
                0077 C--   No need to change start/end index : use array overlap if needed
                0078
                0079 C--   bilinear interpolation (from numerical recipes)
                0080       IF (nu.LE.3) THEN
                0081         uu =          ( (1.-ddx)*(1.-ddy)*var(i1,j1,kl,bi,bj)
                0082      &                +     ddx * ddy    *var(i2,j2,kl,bi,bj) )
                0083      &              + (     ddx *(1.-ddy)*var(i2,j1,kl,bi,bj)
                0084      &                + (1.-ddx)* ddy    *var(i1,j2,kl,bi,bj) )
                0085       ELSE
                0086         klp = MIN(kl+1,Nr)
                0087         u11 = ( var(i1,j1,kl,bi,bj)+var(i1,j1,klp,bi,bj) )*0.5 _d 0
                0088         u21 = ( var(i2,j1,kl,bi,bj)+var(i2,j1,klp,bi,bj) )*0.5 _d 0
                0089         u22 = ( var(i2,j2,kl,bi,bj)+var(i2,j2,klp,bi,bj) )*0.5 _d 0
                0090         u12 = ( var(i1,j2,kl,bi,bj)+var(i1,j2,klp,bi,bj) )*0.5 _d 0
                0091         uu =          ( (1.-ddx)*(1.-ddy)*u11
                0092      &                +     ddx * ddy    *u22 )
                0093      &              + (     ddx *(1.-ddy)*u21
                0094      &                + (1.-ddx)* ddy    *u12 )
                0095       ENDIF
                0096
                0097
                0098       RETURN
                0099       END
                0100
                0101 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0102
                0103       SUBROUTINE FLT_TRILINEAR(
                0104      I                         ix, jy, kz,
                0105      O                         uu,
                0106      I                         var,
                0107      I                         nu, bi, bj, myThid )
                0108
                0109 C     ==================================================================
                0110 C     SUBROUTINE FLT_TRILINEAR
                0111 C     ==================================================================
                0112 C     o Trilinear scheme to interpolate variable to particle position
                0113 C       given by its fractional (real) index ix,jy,kz location
                0114 C       This routine is a straight forward generalization of the
                0115 C       bilinear interpolation scheme.
                0116 C
                0117 C     started: 2004.05.28 Antti Westerlund (antti.westerlund@fimr.fi)
                0118 C              and Sergio Jaramillo (sju@eos.ubc.ca).
                0119 C              (adopted from SUBROUTINE bilinear by Arne Biastoch)
                0120 C     ==================================================================
                0121
                0122 C     !USES:
                0123       IMPLICIT NONE
                0124
                0125 C     == global variables ==
                0126 #include "SIZE.h"
                0127
                0128 C     == routine arguments ==
                0129       _RL ix, jy, kz
                0130       _RL uu
                0131       _RL var(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0132       INTEGER nu, bi, bj, myThid
                0133
                0134 C     == local variables ==
                0135       INTEGER i1, j1, k1, i2, j2, k2
                0136       _RL ddx, ddy, ddz
                0137
                0138 C     == end of interface ==
                0139
                0140 C to choose the u box in which the particle is found
                0141 C nu=0 for T, S
                0142 C nu=1 for u
                0143 C nu=2 for v
                0144 C nu=3 for Vorticity
                0145 C nu=4 for w
                0146
                0147 C--   find x-index according to grid-location of variable
                0148       IF ( MOD(nu,2).EQ.0 ) THEN
                0149         i1 = INT(ix)
                0150         ddx = ix - DFLOAT(i1)
                0151       ELSE
                0152         i1 = NINT(ix)
                0153         ddx = 0.5 _d 0 + ix - DFLOAT(i1)
                0154       ENDIF
                0155 C--   find y-index according to grid-location of variable
                0156       IF ( MOD(nu,4).LE.1 ) THEN
                0157         j1 = INT(jy)
                0158         ddy = jy - DFLOAT(j1)
                0159       ELSE
                0160         j1 = NINT(jy)
                0161         ddy = 0.5 _d 0 + jy - DFLOAT(j1)
                0162       ENDIF
                0163 C--   find z-index according to grid-location of variable
                0164       IF ( nu.LE.3 ) THEN
                0165         k1 = INT(kz)
                0166         ddz = kz - DFLOAT(k1)
                0167       ELSE
                0168         k1 = NINT(kz)
                0169         ddz = 0.5 _d 0 + kz - DFLOAT(k1)
                0170       ENDIF
                0171
                0172 C--   Set the higher index for interpolation
                0173       i2 = i1 + 1
                0174       j2 = j1 + 1
                0175       k2 = k1 + 1
                0176
                0177 C--   No need to change start/end horizontal index : use array overlap if needed
                0178 C--   Need to adjust start/end vertical index :
                0179       k1 = MIN( MAX( k1, 1 ), Nr )
                0180       k2 = MIN( MAX( k2, 1 ), Nr )
                0181
                0182 C     Trilinear interpolation, a straight forward generalization
                0183 C     of the bilinear interpolation scheme.
                0184       uu = (1.-ddz)*( ( (1.-ddx)*(1.-ddy)*var(i1,j1,k1,bi,bj)
                0185      &                +     ddx * ddy    *var(i2,j2,k1,bi,bj) )
                0186      &              + (     ddx *(1.-ddy)*var(i2,j1,k1,bi,bj)
                0187      &                + (1.-ddx)* ddy    *var(i1,j2,k1,bi,bj) ) )
                0188      &   +   ddz   *( ( (1.-ddx)*(1.-ddy)*var(i1,j1,k2,bi,bj)
                0189      &                +     ddx * ddy    *var(i2,j2,k2,bi,bj) )
                0190      &              + (     ddx*(1.-ddy) *var(i2,j1,k2,bi,bj)
                0191      &                + (1.-ddx)* ddy    *var(i1,j2,k2,bi,bj) ) )
                0192
                0193       RETURN
                0194       END
                0195
                0196 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0197
                0198       SUBROUTINE FLT_BILINEAR2D(
                0199      I                           ix, jy,
                0200      O                           uu,
                0201      I                           var,
                0202      I                           nu, bi, bj, myThid )
                0203
                0204 C     ==================================================================
                0205 C     SUBROUTINE FLT_BILINEAR2D
                0206 C     ==================================================================
                0207 C     o Bilinear scheme to interpolate 2-D variable to particle position
                0208 C       given by its fractional (real) index ix,jy location
                0209 C
                0210 C     started: Arne Biastoch abiastoch@ucsd.edu 13-Jan-2000
                0211 C              (adopted from SUBROUTINE bilinear)
                0212 C     ==================================================================
                0213
                0214 C     !USES:
                0215       IMPLICIT NONE
                0216
                0217 C     == global variables ==
                0218 #include "SIZE.h"
                0219
                0220 C     == routine arguments ==
                0221       _RL ix, jy
                0222       _RL uu
                0223       _RL var(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0224       INTEGER nu, bi, bj, myThid
                0225
                0226 C     == local variables ==
                0227       INTEGER i1, j1, i2, j2
                0228       _RL ddx, ddy
                0229
                0230 C     == end of interface ==
                0231
                0232 C to choose the u box in which the particle is found
                0233 C nu=0 for T, S
                0234 C nu=1 for u
                0235 C nu=2 for v
                0236 C nu=3 for Vorticity
                0237 C nu=4 for w
                0238
                0239 C--   find x-index according to grid-location of variable
                0240       IF ( MOD(nu,2).EQ.0 ) THEN
                0241         i1 = INT(ix)
                0242         ddx = ix - DFLOAT(i1)
                0243       ELSE
                0244         i1 = NINT(ix)
                0245         ddx = 0.5 _d 0 + ix - DFLOAT(i1)
                0246       ENDIF
                0247 C--   find y-index according to grid-location of variable
                0248       IF ( MOD(nu,4).LE.1 ) THEN
                0249         j1 = INT(jy)
                0250         ddy = jy - DFLOAT(j1)
                0251       ELSE
                0252         j1 = NINT(jy)
                0253         ddy = 0.5 _d 0 + jy - DFLOAT(j1)
                0254       ENDIF
                0255
                0256 C--   Set the higher index for interpolation
                0257       i2 = i1 + 1
                0258       j2 = j1 + 1
                0259
                0260 C--   No need to change start/end index : use array overlap if needed
                0261
                0262 C bilinear interpolation (from numerical recipes)
                0263       uu =            ( (1.-ddx)*(1.-ddy)*var(i1,j1,bi,bj)
                0264      &                +     ddx * ddy    *var(i2,j2,bi,bj) )
                0265      &              + (     ddx *(1.-ddy)*var(i2,j1,bi,bj)
                0266      &                + (1.-ddx)* ddy    *var(i1,j2,bi,bj) )
                0267
                0268       RETURN
                0269       END