[1793] | 1 | ! Copyright Université Reims Champagnne-Ardenne (2010-2015) |
---|
| 2 | ! contributor: Jérémie Burgalat |
---|
| 3 | ! |
---|
| 4 | ! jeremie.burgalat@univ-reims.fr |
---|
| 5 | ! |
---|
| 6 | ! This software is a computer program whose purpose is to compute multi-variate |
---|
| 7 | ! linear interpolation. |
---|
| 8 | ! |
---|
| 9 | ! This software is governed by the CeCILL-B license under French law and |
---|
| 10 | ! abiding by the rules of distribution of free software. You can use, |
---|
| 11 | ! modify and/ or redistribute the software under the terms of the CeCILL-B |
---|
| 12 | ! license as circulated by CEA, CNRS and INRIA at the following URL |
---|
| 13 | ! "http://www.cecill.info". |
---|
| 14 | ! |
---|
| 15 | ! As a counterpart to the access to the source code and rights to copy, |
---|
| 16 | ! modify and redistribute granted by the license, users are provided only |
---|
| 17 | ! with a limited warranty and the software's author, the holder of the |
---|
| 18 | ! economic rights, and the successive licensors have only limited |
---|
| 19 | ! liability. |
---|
| 20 | ! |
---|
| 21 | ! In this respect, the user's attention is drawn to the risks associated |
---|
| 22 | ! with loading, using, modifying and/or developing or reproducing the |
---|
| 23 | ! software by the user in light of its specific status of free software, |
---|
| 24 | ! that may mean that it is complicated to manipulate, and that also |
---|
| 25 | ! therefore means that it is reserved for developers and experienced |
---|
| 26 | ! professionals having in-depth computer knowledge. Users are therefore |
---|
| 27 | ! encouraged to load and test the software's suitability as regards their |
---|
| 28 | ! requirements in conditions enabling the security of their systems and/or |
---|
| 29 | ! data to be ensured and, more generally, to use and operate it in the |
---|
| 30 | ! same conditions as regards security. |
---|
| 31 | ! |
---|
| 32 | ! The fact that you are presently reading this means that you have had |
---|
| 33 | ! knowledge of the CeCILL-B license and that you accept its terms. |
---|
| 34 | |
---|
| 35 | !! file: lintgen.f90 |
---|
| 36 | !! summary: Linear interpolation generic interfaces |
---|
| 37 | !! author: burgalat |
---|
| 38 | !! date: 2010-2014 |
---|
| 39 | |
---|
| 40 | MODULE LINTGEN |
---|
| 41 | !! Generic linear interpolation module |
---|
| 42 | USE LINT_PREC |
---|
| 43 | USE LINTCORE |
---|
| 44 | IMPLICIT NONE |
---|
| 45 | |
---|
| 46 | PUBLIC :: lint_gen, hdcd_lint_gen |
---|
| 47 | PRIVATE |
---|
| 48 | |
---|
| 49 | !> Interface to multi-variate linear interpolation |
---|
| 50 | !! |
---|
| 51 | !! For __n__, the number of function's variables, user must provide: |
---|
| 52 | !! |
---|
| 53 | !! - __n__ scalar(s) with the coordinate(s) of the point to interpolate. |
---|
| 54 | !! - __n__ vector(s) with the tabulated values of each coordinate. |
---|
| 55 | !! - a single __n__-D array with the value of tabulated function at each |
---|
| 56 | !! tabulated coordinate. |
---|
| 57 | !! - A function, satisfying [[lintcore(module):locate(interface)]] interface, |
---|
| 58 | !! that locates the neighbourhood of the point to interpolate. |
---|
| 59 | !! - A scalar that stores the output value. |
---|
| 60 | INTERFACE lint_gen |
---|
| 61 | MODULE PROCEDURE l1d_gen_sc,l2d_gen_sc,l3d_gen_sc,l4d_gen_sc,l5d_gen_sc |
---|
| 62 | MODULE PROCEDURE l1d_gen_ve,l2d_gen_ve,l3d_gen_ve,l4d_gen_ve,l5d_gen_ve |
---|
| 63 | END INTERFACE |
---|
| 64 | |
---|
| 65 | !> Interface to multi-variate hard-coded linear interpolation |
---|
| 66 | !! |
---|
| 67 | !! Same remarks as in [[lintgen(module):lint_gen(interface)]] apply here. |
---|
| 68 | INTERFACE hdcd_lint_gen |
---|
| 69 | MODULE PROCEDURE hdcdl1d_gen_sc, hdcdl2d_gen_sc, hdcdl3d_gen_sc, & |
---|
| 70 | hdcdl4d_gen_sc, hdcdl5d_gen_sc |
---|
| 71 | MODULE PROCEDURE hdcdl1d_gen_ve, hdcdl2d_gen_ve, hdcdl3d_gen_ve, & |
---|
| 72 | hdcdl4d_gen_ve, hdcdl5d_gen_ve |
---|
| 73 | END INTERFACE |
---|
| 74 | |
---|
| 75 | CONTAINS |
---|
| 76 | |
---|
| 77 | ! GENERIC versions |
---|
| 78 | ! ---------------- |
---|
| 79 | |
---|
| 80 | FUNCTION l1d_gen_sc(x,cx,cv, locator,res) RESULT(ret) |
---|
| 81 | !! Wrapper interface for 1D linear interpolation (scalar) |
---|
| 82 | !! |
---|
| 83 | !! @warning |
---|
| 84 | !! On error, __res__ output value is undefined. |
---|
| 85 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 86 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 87 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values |
---|
| 88 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 89 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 90 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 91 | REAL(kind=wp), DIMENSION(2,2) :: g1d |
---|
| 92 | INTEGER :: i,ix |
---|
| 93 | ret = .false. |
---|
| 94 | ix = locator(x,cx) ; IF (ix == 0) RETURN |
---|
| 95 | DO i=0,1 |
---|
| 96 | g1d(i+1,1) = cx(ix+i) |
---|
| 97 | g1d(i+1,2) = cv(ix+i) |
---|
| 98 | ENDDO |
---|
| 99 | res = lintc_((/x/),g1d) |
---|
| 100 | ret = .true. |
---|
| 101 | END FUNCTION l1d_gen_sc |
---|
| 102 | |
---|
| 103 | FUNCTION l2d_gen_sc(x,y,cx,cy,cv,locator,res) RESULT(ret) |
---|
| 104 | !! Wrapper interface for 2D linear interpolation (scalar) |
---|
| 105 | !! |
---|
| 106 | !! @warning |
---|
| 107 | !! On error, __res__ output value is undefined. |
---|
| 108 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 109 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 110 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 111 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 112 | REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values |
---|
| 113 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 114 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 115 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 116 | REAL(kind=wp), DIMENSION(4,3) :: g2d |
---|
| 117 | INTEGER :: ix0,iy0,i,a,b |
---|
| 118 | ret = .false. |
---|
| 119 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 120 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 121 | DO i=1,4 |
---|
| 122 | a=ix0+MOD((i-1),2) ; g2d(i,1) = cx(a) |
---|
| 123 | b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = cy(b) |
---|
| 124 | g2d(i,3) = cv(a,b) |
---|
| 125 | ENDDO |
---|
| 126 | res = lintc_((/x,y/),g2d) |
---|
| 127 | ret = .true. |
---|
| 128 | END FUNCTION l2d_gen_sc |
---|
| 129 | |
---|
| 130 | FUNCTION l3d_gen_sc(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret) |
---|
| 131 | !! Wrapper interface for 3D linear interpolation (scalar) |
---|
| 132 | !! |
---|
| 133 | !! @warning |
---|
| 134 | !! On error, __res__ output value is undefined. |
---|
| 135 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 136 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 137 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 138 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 139 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 140 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 141 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values |
---|
| 142 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 143 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 144 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 145 | REAL(kind=wp), DIMENSION(8,4) :: g3d |
---|
| 146 | INTEGER :: ix0,iy0,iz0,i,a,b,c |
---|
| 147 | ret = .false. |
---|
| 148 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 149 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 150 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 151 | DO i=1,8 |
---|
| 152 | a=ix0+MOD((i-1),2) ; g3d(i,1) = cx(a) |
---|
| 153 | b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = cy(b) |
---|
| 154 | c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = cz(c) |
---|
| 155 | g3d(i,4) = cv(a,b,c) |
---|
| 156 | ENDDO |
---|
| 157 | res = lintc_((/x,y,z/),g3d) |
---|
| 158 | ret = .true. |
---|
| 159 | END FUNCTION l3d_gen_sc |
---|
| 160 | |
---|
| 161 | FUNCTION l4d_gen_sc(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret) |
---|
| 162 | !! Wrapper interface for 4D linear interpolation (scalar) |
---|
| 163 | !! |
---|
| 164 | !! @warning |
---|
| 165 | !! On error, __res__ output value is undefined. |
---|
| 166 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 167 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 168 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 169 | REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate |
---|
| 170 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 171 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 172 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 173 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 174 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values |
---|
| 175 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 176 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 177 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 178 | REAL(kind=wp), DIMENSION(16,5) :: g4d |
---|
| 179 | INTEGER :: ix0,iy0,iz0,it0 |
---|
| 180 | INTEGER :: i,a,b,c,d |
---|
| 181 | ret = .false. |
---|
| 182 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 183 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 184 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 185 | it0 = locator(t,ct) ; IF (it0 == 0) RETURN |
---|
| 186 | DO i=1,16 |
---|
| 187 | a=ix0+MOD((i-1),2) ; g4d(i,1) = cx(a) |
---|
| 188 | b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = cy(b) |
---|
| 189 | c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = cz(c) |
---|
| 190 | d=it0+MOD((i-1)/8,2) ; g4d(i,4) = ct(d) |
---|
| 191 | g4d(i,5) = cv(a,b,c,d) |
---|
| 192 | ENDDO |
---|
| 193 | res = lintc_((/x,y,z,t/),g4d) |
---|
| 194 | ret = .true. |
---|
| 195 | END FUNCTION l4d_gen_sc |
---|
| 196 | |
---|
| 197 | FUNCTION l5d_gen_sc(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret) |
---|
| 198 | !! Wrapper interface for 5D linear interpolation (scalar) |
---|
| 199 | !! |
---|
| 200 | !! @warning |
---|
| 201 | !! On error, __res__ output value is undefined. |
---|
| 202 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 203 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 204 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 205 | REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate |
---|
| 206 | REAL(kind=wp), INTENT(in) :: w !! W coordinate of the point to interpolate |
---|
| 207 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 208 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 209 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 210 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 211 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates |
---|
| 212 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values |
---|
| 213 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 214 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 215 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 216 | REAL(kind=wp), DIMENSION(32,6) :: g5d |
---|
| 217 | INTEGER :: ix0,iy0,iz0,it0,iw0 |
---|
| 218 | INTEGER :: i,a,b,c,d,e |
---|
| 219 | ret = .false. |
---|
| 220 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 221 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 222 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 223 | it0 = locator(t,ct) ; IF (it0 == 0) RETURN |
---|
| 224 | iw0 = locator(w,cw) ; IF (iw0 == 0) RETURN |
---|
| 225 | DO i=1,32 |
---|
| 226 | a=ix0+MOD((i-1),2) ; g5d(i,1) = cx(a) |
---|
| 227 | b=iy0+MOD((i-1)/2,2) ; g5d(i,2) = cy(b) |
---|
| 228 | c=iz0+MOD((i-1)/4,2) ; g5d(i,3) = cz(c) |
---|
| 229 | d=it0+MOD((i-1)/8,2) ; g5d(i,4) = ct(d) |
---|
| 230 | e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = cw(e) |
---|
| 231 | g5d(i,6) = cv(a,b,c,d,e) |
---|
| 232 | ENDDO |
---|
| 233 | res = lintc_((/x,y,z,t,w/),g5d) |
---|
| 234 | ret = .true. |
---|
| 235 | END FUNCTION l5d_gen_sc |
---|
| 236 | |
---|
| 237 | ! HARD-CODED versions |
---|
| 238 | ! ------------------- |
---|
| 239 | |
---|
| 240 | FUNCTION hdcdl1d_gen_sc(x,cx,cv,locator,res) RESULT(ret) |
---|
| 241 | !! Hard-coded 1D linear interpolation (scalar) |
---|
| 242 | !! |
---|
| 243 | !! @warning |
---|
| 244 | !! On error, __res__ output value is undefined. |
---|
| 245 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 246 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 247 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values |
---|
| 248 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 249 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 250 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 251 | REAL(kind=wp) :: xd |
---|
| 252 | INTEGER :: ix0,ix1 |
---|
| 253 | ret = .false. |
---|
| 254 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 255 | ix1 = ix0+1 |
---|
| 256 | xd = (x-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 257 | res = cv(ix0)*(1d0-xd)+cv(ix1)*xd |
---|
| 258 | ret = .true. |
---|
| 259 | END FUNCTION hdcdl1d_gen_sc |
---|
| 260 | |
---|
| 261 | FUNCTION hdcdl2d_gen_sc(x,y,cx,cy,cv,locator,res) RESULT(ret) |
---|
| 262 | !! Hard-coded 2D linear interpolation (scalar) |
---|
| 263 | !! |
---|
| 264 | !! @warning |
---|
| 265 | !! On error, __res__ output value is undefined. |
---|
| 266 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 267 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 268 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 269 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 270 | REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values |
---|
| 271 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 272 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 273 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 274 | REAL(kind=wp) :: xd,yd |
---|
| 275 | REAL(kind=wp) :: f0,f1 |
---|
| 276 | INTEGER :: ix0,iy0,ix1,iy1 |
---|
| 277 | ret = .false. |
---|
| 278 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 279 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 280 | ix1 = ix0 + 1 ; iy1 = iy0 + 1 |
---|
| 281 | xd = (x-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 282 | yd = (y-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 283 | f0 = cv(ix0,iy0)*(1d0-xd)+cv(ix1,iy0)*xd |
---|
| 284 | f1 = cv(ix0,iy1)*(1d0-xd)+cv(ix1,iy1)*xd |
---|
| 285 | res = f0*(1d0-yd)+f1*yd |
---|
| 286 | ret = .true. |
---|
| 287 | END FUNCTION hdcdl2d_gen_sc |
---|
| 288 | |
---|
| 289 | FUNCTION hdcdl3d_gen_sc(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret) |
---|
| 290 | !! Hard-coded 3D linear interpolation (scalar) |
---|
| 291 | !! |
---|
| 292 | !! @warning |
---|
| 293 | !! On error, __res__ output value is undefined. |
---|
| 294 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 295 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 296 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 297 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 298 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 299 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 300 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values |
---|
| 301 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 302 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 303 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 304 | REAL(kind=wp) :: xd,yd,zd |
---|
| 305 | REAL(kind=wp) :: f00,f10,f01,f11,f0,f1 |
---|
| 306 | INTEGER :: ix0,iy0,iz0,ix1,iy1,iz1 |
---|
| 307 | ret = .false. |
---|
| 308 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 309 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 310 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 311 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 |
---|
| 312 | xd = (x-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 313 | yd = (y-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 314 | zd = (z-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 315 | |
---|
| 316 | f00 = cv(ix0,iy0,iz0)*(1d0-xd)+cv(ix1,iy0,iz0)*xd |
---|
| 317 | f10 = cv(ix0,iy1,iz0)*(1d0-xd)+cv(ix1,iy1,iz0)*xd |
---|
| 318 | f01 = cv(ix0,iy0,iz1)*(1d0-xd)+cv(ix1,iy0,iz1)*xd |
---|
| 319 | f11 = cv(ix0,iy1,iz1)*(1d0-xd)+cv(ix1,iy1,iz1)*xd |
---|
| 320 | |
---|
| 321 | f0 = f00 *(1d0-yd)+f10*yd |
---|
| 322 | f1 = f00 *(1d0-yd)+f10*yd |
---|
| 323 | |
---|
| 324 | res = f0*(1d0-zd)+f1*zd |
---|
| 325 | ret = .true. |
---|
| 326 | END FUNCTION hdcdl3d_gen_sc |
---|
| 327 | |
---|
| 328 | FUNCTION hdcdl4d_gen_sc(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret) |
---|
| 329 | !! Hard-coded 4D linear interpolation (scalar) |
---|
| 330 | !! |
---|
| 331 | !! @warning |
---|
| 332 | !! On error, __res__ output value is undefined. |
---|
| 333 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 334 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 335 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 336 | REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate |
---|
| 337 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 338 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 339 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 340 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 341 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values |
---|
| 342 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 343 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 344 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 345 | REAL(kind=wp) :: xd,yd,zd,td |
---|
| 346 | REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, & |
---|
| 347 | f00,f10,f01,f11,f0,f1 |
---|
| 348 | INTEGER :: ix0,iy0,iz0,it0,ix1,iy1,iz1,it1 |
---|
| 349 | ret = .false. |
---|
| 350 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 351 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 352 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 353 | it0 = locator(t,ct) ; IF (it0 == 0) RETURN |
---|
| 354 | |
---|
| 355 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1 |
---|
| 356 | xd = (x-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 357 | yd = (y-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 358 | zd = (z-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 359 | td = (t-ct(it0))/(ct(it1)-ct(it0)) |
---|
| 360 | |
---|
| 361 | f000 = cv(ix0,iy0,iz0,it0)*(1d0-xd)+cv(ix1,iy0,iz0,it0)*xd |
---|
| 362 | f100 = cv(ix0,iy1,iz0,it0)*(1d0-xd)+cv(ix1,iy1,iz0,it0)*xd |
---|
| 363 | f010 = cv(ix0,iy0,iz1,it0)*(1d0-xd)+cv(ix1,iy0,iz1,it0)*xd |
---|
| 364 | f110 = cv(ix0,iy1,iz1,it0)*(1d0-xd)+cv(ix1,iy1,iz1,it0)*xd |
---|
| 365 | f001 = cv(ix0,iy0,iz0,it1)*(1d0-xd)+cv(ix1,iy0,iz0,it1)*xd |
---|
| 366 | f101 = cv(ix0,iy1,iz0,it1)*(1d0-xd)+cv(ix1,iy1,iz0,it1)*xd |
---|
| 367 | f011 = cv(ix0,iy0,iz1,it1)*(1d0-xd)+cv(ix1,iy0,iz1,it1)*xd |
---|
| 368 | f111 = cv(ix0,iy1,iz1,it1)*(1d0-xd)+cv(ix1,iy1,iz1,it1)*xd |
---|
| 369 | |
---|
| 370 | f00 = f000*(1d0-yd)+f100*yd |
---|
| 371 | f10 = f010*(1d0-yd)+f110*yd |
---|
| 372 | f01 = f001*(1d0-yd)+f101*yd |
---|
| 373 | f11 = f011*(1d0-yd)+f111*yd |
---|
| 374 | |
---|
| 375 | f0 = f00 *(1d0-zd)+f10*zd |
---|
| 376 | f1 = f01 *(1d0-zd)+f11*zd |
---|
| 377 | |
---|
| 378 | res = f0*(1d0-td)+f1*td |
---|
| 379 | ret = .true. |
---|
| 380 | END FUNCTION hdcdl4d_gen_sc |
---|
| 381 | |
---|
| 382 | FUNCTION hdcdl5d_gen_sc(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret) |
---|
| 383 | !! Hard-coded 5D linear interpolation (scalar) |
---|
| 384 | !! |
---|
| 385 | !! @warning |
---|
| 386 | !! On error, __res__ output value is undefined. |
---|
| 387 | REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate |
---|
| 388 | REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate |
---|
| 389 | REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate |
---|
| 390 | REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate |
---|
| 391 | REAL(kind=wp), INTENT(in) :: w !! W coordinate of the point to interpolate |
---|
| 392 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 393 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 394 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 395 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 396 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates |
---|
| 397 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values |
---|
| 398 | REAL(kind=wp), INTENT(out) :: res !! Interpolated value |
---|
| 399 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 400 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 401 | REAL(kind=wp) :: xd,yd,zd,td,wd |
---|
| 402 | REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, & |
---|
| 403 | f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, & |
---|
| 404 | f000,f100,f010,f110,f001,f101,f011,f111, & |
---|
| 405 | f00,f10,f01,f11,f0,f1 |
---|
| 406 | INTEGER :: ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1 |
---|
| 407 | ret = .false. |
---|
| 408 | ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN |
---|
| 409 | iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN |
---|
| 410 | iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN |
---|
| 411 | it0 = locator(t,ct) ; IF (it0 == 0) RETURN |
---|
| 412 | iw0 = locator(w,cw) ; IF (iw0 == 0) RETURN |
---|
| 413 | |
---|
| 414 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1 ; iw1=iw0+1 |
---|
| 415 | xd = (x-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 416 | yd = (y-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 417 | zd = (z-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 418 | td = (t-ct(it0))/(ct(it1)-ct(it0)) |
---|
| 419 | wd = (w-cw(iw0))/(cw(iw1)-cw(iw0)) |
---|
| 420 | |
---|
| 421 | f0000 = cv(ix0,iy0,iz0,it0,iw0)*(1d0-xd)+cv(ix1,iy0,iz0,it0,iw0)*xd |
---|
| 422 | f1000 = cv(ix0,iy1,iz0,it0,iw0)*(1d0-xd)+cv(ix0,iy1,iz0,it0,iw0)*xd |
---|
| 423 | f0100 = cv(ix0,iy0,iz1,it0,iw0)*(1d0-xd)+cv(ix1,iy0,iz1,it0,iw0)*xd |
---|
| 424 | f1100 = cv(ix0,iy1,iz1,it0,iw0)*(1d0-xd)+cv(ix1,iy1,iz1,it0,iw0)*xd |
---|
| 425 | f0010 = cv(ix0,iy0,iz0,it1,iw0)*(1d0-xd)+cv(ix1,iy0,iz0,it1,iw0)*xd |
---|
| 426 | f1010 = cv(ix0,iy1,iz0,it1,iw0)*(1d0-xd)+cv(ix1,iy1,iz0,it1,iw0)*xd |
---|
| 427 | f0110 = cv(ix0,iy0,iz1,it1,iw0)*(1d0-xd)+cv(ix1,iy0,iz1,it1,iw0)*xd |
---|
| 428 | f1110 = cv(ix0,iy1,iz1,it1,iw0)*(1d0-xd)+cv(ix1,iy1,iz1,it1,iw0)*xd |
---|
| 429 | f0001 = cv(ix0,iy0,iz0,it0,iw1)*(1d0-xd)+cv(ix1,iy0,iz0,it0,iw1)*xd |
---|
| 430 | f1001 = cv(ix0,iy1,iz0,it0,iw1)*(1d0-xd)+cv(ix0,iy1,iz0,it0,iw1)*xd |
---|
| 431 | f0101 = cv(ix0,iy0,iz1,it0,iw1)*(1d0-xd)+cv(ix1,iy0,iz1,it0,iw1)*xd |
---|
| 432 | f1101 = cv(ix0,iy1,iz1,it0,iw1)*(1d0-xd)+cv(ix1,iy1,iz1,it0,iw1)*xd |
---|
| 433 | f0011 = cv(ix0,iy0,iz0,it1,iw1)*(1d0-xd)+cv(ix1,iy0,iz0,it1,iw1)*xd |
---|
| 434 | f1011 = cv(ix0,iy1,iz0,it1,iw1)*(1d0-xd)+cv(ix1,iy1,iz0,it1,iw1)*xd |
---|
| 435 | f0111 = cv(ix0,iy0,iz1,it1,iw1)*(1d0-xd)+cv(ix1,iy0,iz1,it1,iw1)*xd |
---|
| 436 | f1111 = cv(ix0,iy1,iz1,it1,iw1)*(1d0-xd)+cv(ix1,iy1,iz1,it1,iw1)*xd |
---|
| 437 | |
---|
| 438 | f000 = f0000*(1d0-yd) + f1000*yd |
---|
| 439 | f100 = f0100*(1d0-yd) + f1100*yd |
---|
| 440 | f010 = f0010*(1d0-yd) + f1010*yd |
---|
| 441 | f110 = f0110*(1d0-yd) + f1110*yd |
---|
| 442 | f101 = f0101*(1d0-yd) + f1101*yd |
---|
| 443 | f011 = f0011*(1d0-yd) + f1011*yd |
---|
| 444 | f111 = f0111*(1d0-yd) + f1111*yd |
---|
| 445 | |
---|
| 446 | f00 = f000*(1d0-zd)+f100*zd |
---|
| 447 | f10 = f010*(1d0-zd)+f110*zd |
---|
| 448 | f01 = f001*(1d0-zd)+f101*zd |
---|
| 449 | f11 = f011*(1d0-zd)+f111*zd |
---|
| 450 | |
---|
| 451 | f0 = f00 *(1d0-td)+f10*td |
---|
| 452 | f1 = f01 *(1d0-td)+f11*td |
---|
| 453 | |
---|
| 454 | res = f0*(1d0-wd)+f1*wd |
---|
| 455 | ret = .true. |
---|
| 456 | END FUNCTION hdcdl5d_gen_sc |
---|
| 457 | |
---|
| 458 | !-------------------- |
---|
| 459 | ! VECTORIZED VERSIONS |
---|
| 460 | !-------------------- |
---|
| 461 | |
---|
| 462 | ! GENERIC versions |
---|
| 463 | ! ---------------- |
---|
| 464 | |
---|
| 465 | FUNCTION l1d_gen_ve(x,cx,cv,locator,res) RESULT(ret) |
---|
| 466 | !! Wrapper interface for 1D linear interpolation (vector) |
---|
| 467 | !! |
---|
| 468 | !! @warning |
---|
| 469 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 470 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 471 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 472 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values |
---|
| 473 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 474 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 475 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 476 | REAL(kind=wp), DIMENSION(2,2) :: g1d |
---|
| 477 | INTEGER :: n,j,i,ix |
---|
| 478 | ret = .false. |
---|
| 479 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 480 | DO j=1,n |
---|
| 481 | ix = locator(x(j),cx) |
---|
| 482 | IF (ix == 0) THEN |
---|
| 483 | DEALLOCATE(res) ; RETURN |
---|
| 484 | ENDIF |
---|
| 485 | DO i=0,1 |
---|
| 486 | g1d(i+1,1) = cx(ix+i) |
---|
| 487 | g1d(i+1,2) = cv(ix+i) |
---|
| 488 | ENDDO |
---|
| 489 | res(j) = lintc_((/x(j)/),g1d) |
---|
| 490 | ENDDO |
---|
| 491 | ret = .true. |
---|
| 492 | END FUNCTION l1d_gen_ve |
---|
| 493 | |
---|
| 494 | FUNCTION l2d_gen_ve(x,y,cx,cy,cv,locator,res) RESULT(ret) |
---|
| 495 | !! Wrapper interface for 2D linear interpolation (vector) |
---|
| 496 | !! |
---|
| 497 | !! @warning |
---|
| 498 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 499 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 500 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 501 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 502 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 503 | REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values |
---|
| 504 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 505 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 506 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 507 | REAL(kind=wp), DIMENSION(4,3) :: g2d |
---|
| 508 | INTEGER :: n,j,ix0,iy0,i,a,b |
---|
| 509 | ret = .false. |
---|
| 510 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 511 | DO j=1,n |
---|
| 512 | ix0 = locator(x(j),cx) |
---|
| 513 | iy0 = locator(y(j),cy) |
---|
| 514 | IF (ix0 == 0 .OR. iy0 == 0) THEN |
---|
| 515 | DEALLOCATE(res) ; RETURN |
---|
| 516 | ENDIF |
---|
| 517 | DO i=1,4 |
---|
| 518 | a=ix0+MOD((i-1),2) ; g2d(i,1) = cx(a) |
---|
| 519 | b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = cy(b) |
---|
| 520 | g2d(i,3) = cv(a,b) |
---|
| 521 | ENDDO |
---|
| 522 | res(j) = lintc_((/x(j),y(j)/),g2d) |
---|
| 523 | ENDDO |
---|
| 524 | ret = .true. |
---|
| 525 | END FUNCTION l2d_gen_ve |
---|
| 526 | |
---|
| 527 | FUNCTION l3d_gen_ve(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret) |
---|
| 528 | !! Wrapper interface for 3D linear interpolation (vector) |
---|
| 529 | !! |
---|
| 530 | !! @warning |
---|
| 531 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 532 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 533 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 534 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 535 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 536 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 537 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 538 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values |
---|
| 539 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 540 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 541 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 542 | REAL(kind=wp), DIMENSION(8,4) :: g3d |
---|
| 543 | INTEGER :: ix0,iy0,iz0 |
---|
| 544 | INTEGER :: n,j,i,a,b,c |
---|
| 545 | ret = .false. |
---|
| 546 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 547 | DO j=1,n |
---|
| 548 | ix0 = locator(x(j),cx) |
---|
| 549 | iy0 = locator(y(j),cy) |
---|
| 550 | iz0 = locator(z(j),cz) |
---|
| 551 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0) THEN |
---|
| 552 | DEALLOCATE(res) ; RETURN |
---|
| 553 | ENDIF |
---|
| 554 | DO i=1,8 |
---|
| 555 | a=ix0+MOD((i-1),2) ; g3d(i,1) = cx(a) |
---|
| 556 | b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = cy(b) |
---|
| 557 | c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = cz(c) |
---|
| 558 | g3d(i,4) = cv(a,b,c) |
---|
| 559 | ENDDO |
---|
| 560 | res(j) = lintc_((/x(j),y(j),z(j)/),g3d) |
---|
| 561 | ENDDO |
---|
| 562 | ret = .true. |
---|
| 563 | END FUNCTION l3d_gen_ve |
---|
| 564 | |
---|
| 565 | FUNCTION l4d_gen_ve(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret) |
---|
| 566 | !! Wrapper interface for 4D linear interpolation (vector) |
---|
| 567 | !! |
---|
| 568 | !! @warning |
---|
| 569 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 570 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 571 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 572 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 573 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate |
---|
| 574 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 575 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 576 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 577 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 578 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values |
---|
| 579 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 580 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 581 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 582 | REAL(kind=wp), DIMENSION(16,5) :: g4d |
---|
| 583 | INTEGER :: ix0,iy0,iz0,it0 |
---|
| 584 | INTEGER :: n,j,i,a,b,c,d |
---|
| 585 | ret = .false. |
---|
| 586 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 587 | DO j=1,n |
---|
| 588 | ix0 = locator(x(j),cx) |
---|
| 589 | iy0 = locator(y(j),cy) |
---|
| 590 | iz0 = locator(z(j),cz) |
---|
| 591 | it0 = locator(t(j),ct) |
---|
| 592 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0) THEN |
---|
| 593 | DEALLOCATE(res) ; RETURN |
---|
| 594 | ENDIF |
---|
| 595 | DO i=1,16 |
---|
| 596 | a=ix0+MOD((i-1),2) ; g4d(i,1) = cx(a) |
---|
| 597 | b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = cy(b) |
---|
| 598 | c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = cz(c) |
---|
| 599 | d=it0+MOD((i-1)/8,2) ; g4d(i,4) = ct(d) |
---|
| 600 | g4d(i,5) = cv(a,b,c,d) |
---|
| 601 | ENDDO |
---|
| 602 | res(j) = lintc_((/x(j),y(j),z(j),t(j)/),g4d) |
---|
| 603 | ENDDO |
---|
| 604 | ret = .true. |
---|
| 605 | END FUNCTION l4d_gen_ve |
---|
| 606 | |
---|
| 607 | FUNCTION l5d_gen_ve(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret) |
---|
| 608 | !! Wrapper interface for 5D linear interpolation (vector) |
---|
| 609 | !! |
---|
| 610 | !! @warning |
---|
| 611 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 612 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 613 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 614 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 615 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate |
---|
| 616 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: w !! W coordinate of the points to interpolate |
---|
| 617 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 618 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 619 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 620 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 621 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates |
---|
| 622 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values |
---|
| 623 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 624 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 625 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 626 | REAL(kind=wp), DIMENSION(32,6) :: g5d |
---|
| 627 | INTEGER :: ix0,iy0,iz0,it0,iw0 |
---|
| 628 | INTEGER :: n,j,i,a,b,c,d,e |
---|
| 629 | ret = .false. |
---|
| 630 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 631 | DO j=1,n |
---|
| 632 | ix0 = locator(x(j),cx) |
---|
| 633 | iy0 = locator(y(j),cy) |
---|
| 634 | iz0 = locator(z(j),cz) |
---|
| 635 | it0 = locator(t(j),ct) |
---|
| 636 | iw0 = locator(w(j),cw) |
---|
| 637 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0 .OR. iw0 == 0) THEN |
---|
| 638 | DEALLOCATE(res) ; RETURN |
---|
| 639 | ENDIF |
---|
| 640 | DO i=1,32 |
---|
| 641 | a=ix0+MOD((i-1),2) ; g5d(i,1) = cx(a) |
---|
| 642 | b=iy0+MOD((i-1)/2,2) ; g5d(i,2) = cy(b) |
---|
| 643 | c=iz0+MOD((i-1)/4,2) ; g5d(i,3) = cz(c) |
---|
| 644 | d=it0+MOD((i-1)/8,2) ; g5d(i,4) = ct(d) |
---|
| 645 | e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = cw(e) |
---|
| 646 | g5d(i,6) = cv(a,b,c,d,e) |
---|
| 647 | ENDDO |
---|
| 648 | res(j) = lintc_((/x,y,z,t,w/),g5d) |
---|
| 649 | ENDDO |
---|
| 650 | ret = .true. |
---|
| 651 | END FUNCTION l5d_gen_ve |
---|
| 652 | |
---|
| 653 | ! HARD-CODED versions |
---|
| 654 | ! ------------------- |
---|
| 655 | |
---|
| 656 | FUNCTION hdcdl1d_gen_ve(x,cx,cv,locator,res) RESULT(ret) |
---|
| 657 | !! Hard-coded 1D linear interpolation (vector) |
---|
| 658 | !! |
---|
| 659 | !! @warning |
---|
| 660 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 661 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 662 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 663 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values |
---|
| 664 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 665 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 666 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 667 | REAL(kind=wp) :: xd |
---|
| 668 | INTEGER :: ix0,ix1 |
---|
| 669 | INTEGER :: n,j |
---|
| 670 | ret = .false. |
---|
| 671 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 672 | DO j=1,n |
---|
| 673 | ix0 = locator(x(j),cx) |
---|
| 674 | IF (ix0 == 0) THEN |
---|
| 675 | DEALLOCATE(res) ; RETURN |
---|
| 676 | ENDIF |
---|
| 677 | ix1 = ix0+1 |
---|
| 678 | xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 679 | res(j) = cv(ix0)*(1d0-xd)+cv(ix1)*xd |
---|
| 680 | ENDDO |
---|
| 681 | ret = .true. |
---|
| 682 | END FUNCTION hdcdl1d_gen_ve |
---|
| 683 | |
---|
| 684 | FUNCTION hdcdl2d_gen_ve(x,y,cx,cy,cv,locator,res) RESULT(ret) |
---|
| 685 | !! Hard-coded 2D linear interpolation (vector) |
---|
| 686 | !! |
---|
| 687 | !! @warning |
---|
| 688 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 689 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 690 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 691 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 692 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 693 | REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values |
---|
| 694 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 695 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 696 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 697 | REAL(kind=wp) :: xd,yd |
---|
| 698 | REAL(kind=wp) :: f0,f1 |
---|
| 699 | INTEGER :: ix0,iy0,ix1,iy1 |
---|
| 700 | INTEGER :: n,j |
---|
| 701 | ret = .false. |
---|
| 702 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 703 | DO j=1,n |
---|
| 704 | ix0 = locator(x(j),cx) |
---|
| 705 | iy0 = locator(y(j),cy) |
---|
| 706 | IF (ix0 == 0 .OR. iy0 == 0) THEN |
---|
| 707 | DEALLOCATE(res) ; RETURN |
---|
| 708 | ENDIF |
---|
| 709 | ix1 = ix0 + 1 ; iy1 = iy0 + 1 |
---|
| 710 | xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 711 | yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 712 | f0 = cv(ix0,iy0)*(1d0-xd)+cv(ix1,iy0)*xd |
---|
| 713 | f1 = cv(ix0,iy1)*(1d0-xd)+cv(ix1,iy1)*xd |
---|
| 714 | res(j) = f0*(1d0-yd)+f1*yd |
---|
| 715 | ENDDO |
---|
| 716 | ret = .true. |
---|
| 717 | END FUNCTION hdcdl2d_gen_ve |
---|
| 718 | |
---|
| 719 | FUNCTION hdcdl3d_gen_ve(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret) |
---|
| 720 | !! Hard-coded 3D linear interpolation (vector) |
---|
| 721 | !! |
---|
| 722 | !! @warning |
---|
| 723 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 724 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 725 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 726 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 727 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 728 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 729 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 730 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values |
---|
| 731 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 732 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 733 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 734 | REAL(kind=wp) :: xd,yd,zd |
---|
| 735 | REAL(kind=wp) :: f00,f10,f01,f11,f0,f1 |
---|
| 736 | INTEGER :: ix0,iy0,iz0,ix1,iy1,iz1 |
---|
| 737 | INTEGER :: n,j |
---|
| 738 | ret = .false. |
---|
| 739 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 740 | DO j=1,n |
---|
| 741 | ix0 = locator(x(j),cx) |
---|
| 742 | iy0 = locator(y(j),cy) |
---|
| 743 | iz0 = locator(z(j),cz) |
---|
| 744 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0) THEN |
---|
| 745 | DEALLOCATE(res) ; RETURN |
---|
| 746 | ENDIF |
---|
| 747 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 |
---|
| 748 | xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 749 | yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 750 | zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 751 | |
---|
| 752 | f00 = cv(ix0,iy0,iz0)*(1d0-xd)+cv(ix1,iy0,iz0)*xd |
---|
| 753 | f10 = cv(ix0,iy1,iz0)*(1d0-xd)+cv(ix1,iy1,iz0)*xd |
---|
| 754 | f01 = cv(ix0,iy0,iz1)*(1d0-xd)+cv(ix1,iy0,iz1)*xd |
---|
| 755 | f11 = cv(ix0,iy1,iz1)*(1d0-xd)+cv(ix1,iy1,iz1)*xd |
---|
| 756 | |
---|
| 757 | f0 = f00 *(1d0-yd)+f10*yd |
---|
| 758 | f1 = f00 *(1d0-yd)+f10*yd |
---|
| 759 | |
---|
| 760 | res(j) = f0*(1d0-zd)+f1*zd |
---|
| 761 | ENDDO |
---|
| 762 | ret = .true. |
---|
| 763 | END FUNCTION hdcdl3d_gen_ve |
---|
| 764 | |
---|
| 765 | FUNCTION hdcdl4d_gen_ve(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret) |
---|
| 766 | !! Hard-coded 4D linear interpolation (vector) |
---|
| 767 | !! |
---|
| 768 | !! @warning |
---|
| 769 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 770 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 771 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 772 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 773 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate |
---|
| 774 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 775 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 776 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 777 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 778 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values |
---|
| 779 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 780 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 781 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 782 | REAL(kind=wp) :: xd,yd,zd,td |
---|
| 783 | REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, & |
---|
| 784 | f00,f10,f01,f11,f0,f1 |
---|
| 785 | INTEGER :: ix0,iy0,iz0,it0,ix1,iy1,iz1,it1 |
---|
| 786 | INTEGER :: n,j |
---|
| 787 | ret = .false. |
---|
| 788 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 789 | DO j=1,n |
---|
| 790 | ix0 = locator(x(j),cx) |
---|
| 791 | iy0 = locator(y(j),cy) |
---|
| 792 | iz0 = locator(z(j),cz) |
---|
| 793 | it0 = locator(t(j),ct) |
---|
| 794 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0) THEN |
---|
| 795 | DEALLOCATE(res) ; RETURN |
---|
| 796 | ENDIF |
---|
| 797 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1 |
---|
| 798 | xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 799 | yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 800 | zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 801 | td = (t(j)-ct(it0))/(ct(it1)-ct(it0)) |
---|
| 802 | |
---|
| 803 | f000 = cv(ix0,iy0,iz0,it0)*(1d0-xd)+cv(ix1,iy0,iz0,it0)*xd |
---|
| 804 | f100 = cv(ix0,iy1,iz0,it0)*(1d0-xd)+cv(ix1,iy1,iz0,it0)*xd |
---|
| 805 | f010 = cv(ix0,iy0,iz1,it0)*(1d0-xd)+cv(ix1,iy0,iz1,it0)*xd |
---|
| 806 | f110 = cv(ix0,iy1,iz1,it0)*(1d0-xd)+cv(ix1,iy1,iz1,it0)*xd |
---|
| 807 | f001 = cv(ix0,iy0,iz0,it1)*(1d0-xd)+cv(ix1,iy0,iz0,it1)*xd |
---|
| 808 | f101 = cv(ix0,iy1,iz0,it1)*(1d0-xd)+cv(ix1,iy1,iz0,it1)*xd |
---|
| 809 | f011 = cv(ix0,iy0,iz1,it1)*(1d0-xd)+cv(ix1,iy0,iz1,it1)*xd |
---|
| 810 | f111 = cv(ix0,iy1,iz1,it1)*(1d0-xd)+cv(ix1,iy1,iz1,it1)*xd |
---|
| 811 | |
---|
| 812 | f00 = f000*(1d0-yd)+f100*yd |
---|
| 813 | f10 = f010*(1d0-yd)+f110*yd |
---|
| 814 | f01 = f001*(1d0-yd)+f101*yd |
---|
| 815 | f11 = f011*(1d0-yd)+f111*yd |
---|
| 816 | |
---|
| 817 | f0 = f00 *(1d0-zd)+f10*zd |
---|
| 818 | f1 = f01 *(1d0-zd)+f11*zd |
---|
| 819 | |
---|
| 820 | res(j) = f0*(1d0-td)+f1*td |
---|
| 821 | ENDDO |
---|
| 822 | ret = .true. |
---|
| 823 | END FUNCTION hdcdl4d_gen_ve |
---|
| 824 | |
---|
| 825 | FUNCTION hdcdl5d_gen_ve(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret) |
---|
| 826 | !! Hard-coded 5D linear interpolation (vector) |
---|
| 827 | !! |
---|
| 828 | !! @warning |
---|
| 829 | !! On error, __res__ output vector is undefined (i.e. not allocated). |
---|
| 830 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate |
---|
| 831 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate |
---|
| 832 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate |
---|
| 833 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate |
---|
| 834 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: w !! W coordinate of the points to interpolate |
---|
| 835 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates |
---|
| 836 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates |
---|
| 837 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates |
---|
| 838 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates |
---|
| 839 | REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates |
---|
| 840 | REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values |
---|
| 841 | REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values |
---|
| 842 | PROCEDURE(locate) :: locator !! Locator function |
---|
| 843 | LOGICAL :: ret !! .true. on success, .false. otherwise |
---|
| 844 | REAL(kind=wp) :: xd,yd,zd,td,wd |
---|
| 845 | REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, & |
---|
| 846 | f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, & |
---|
| 847 | f000,f100,f010,f110,f001,f101,f011,f111, & |
---|
| 848 | f00,f10,f01,f11,f0,f1 |
---|
| 849 | INTEGER :: ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1 |
---|
| 850 | INTEGER :: n,j |
---|
| 851 | ret = .false. |
---|
| 852 | n = SIZE(x) ; ALLOCATE(res(n)) |
---|
| 853 | DO j=1,n |
---|
| 854 | ix0 = locator(x(j),cx) |
---|
| 855 | iy0 = locator(y(j),cy) |
---|
| 856 | iz0 = locator(z(j),cz) |
---|
| 857 | it0 = locator(t(j),ct) |
---|
| 858 | iw0 = locator(w(j),cw) |
---|
| 859 | IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0 .OR. iw0 == 0) THEN |
---|
| 860 | DEALLOCATE(res) ; RETURN |
---|
| 861 | ENDIF |
---|
| 862 | ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; |
---|
| 863 | xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0)) |
---|
| 864 | yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0)) |
---|
| 865 | zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0)) |
---|
| 866 | td = (t(j)-ct(it0))/(ct(it1)-ct(it0)) |
---|
| 867 | wd = (w(j)-cw(it0))/(cw(it1)-cw(it0)) |
---|
| 868 | |
---|
| 869 | f0000 = cv(ix0,iy0,iz0,it0,iw0)*(1d0-xd)+cv(ix1,iy0,iz0,it0,iw0)*xd |
---|
| 870 | f1000 = cv(ix0,iy1,iz0,it0,iw0)*(1d0-xd)+cv(ix0,iy1,iz0,it0,iw0)*xd |
---|
| 871 | f0100 = cv(ix0,iy0,iz1,it0,iw0)*(1d0-xd)+cv(ix1,iy0,iz1,it0,iw0)*xd |
---|
| 872 | f1100 = cv(ix0,iy1,iz1,it0,iw0)*(1d0-xd)+cv(ix1,iy1,iz1,it0,iw0)*xd |
---|
| 873 | f0010 = cv(ix0,iy0,iz0,it1,iw0)*(1d0-xd)+cv(ix1,iy0,iz0,it1,iw0)*xd |
---|
| 874 | f1010 = cv(ix0,iy1,iz0,it1,iw0)*(1d0-xd)+cv(ix1,iy1,iz0,it1,iw0)*xd |
---|
| 875 | f0110 = cv(ix0,iy0,iz1,it1,iw0)*(1d0-xd)+cv(ix1,iy0,iz1,it1,iw0)*xd |
---|
| 876 | f1110 = cv(ix0,iy1,iz1,it1,iw0)*(1d0-xd)+cv(ix1,iy1,iz1,it1,iw0)*xd |
---|
| 877 | f0001 = cv(ix0,iy0,iz0,it0,iw1)*(1d0-xd)+cv(ix1,iy0,iz0,it0,iw1)*xd |
---|
| 878 | f1001 = cv(ix0,iy1,iz0,it0,iw1)*(1d0-xd)+cv(ix0,iy1,iz0,it0,iw1)*xd |
---|
| 879 | f0101 = cv(ix0,iy0,iz1,it0,iw1)*(1d0-xd)+cv(ix1,iy0,iz1,it0,iw1)*xd |
---|
| 880 | f1101 = cv(ix0,iy1,iz1,it0,iw1)*(1d0-xd)+cv(ix1,iy1,iz1,it0,iw1)*xd |
---|
| 881 | f0011 = cv(ix0,iy0,iz0,it1,iw1)*(1d0-xd)+cv(ix1,iy0,iz0,it1,iw1)*xd |
---|
| 882 | f1011 = cv(ix0,iy1,iz0,it1,iw1)*(1d0-xd)+cv(ix1,iy1,iz0,it1,iw1)*xd |
---|
| 883 | f0111 = cv(ix0,iy0,iz1,it1,iw1)*(1d0-xd)+cv(ix1,iy0,iz1,it1,iw1)*xd |
---|
| 884 | f1111 = cv(ix0,iy1,iz1,it1,iw1)*(1d0-xd)+cv(ix1,iy1,iz1,it1,iw1)*xd |
---|
| 885 | |
---|
| 886 | f000 = f0000*(1d0-yd) + f1000*yd |
---|
| 887 | f100 = f0100*(1d0-yd) + f1100*yd |
---|
| 888 | f010 = f0010*(1d0-yd) + f1010*yd |
---|
| 889 | f110 = f0110*(1d0-yd) + f1110*yd |
---|
| 890 | f101 = f0101*(1d0-yd) + f1101*yd |
---|
| 891 | f011 = f0011*(1d0-yd) + f1011*yd |
---|
| 892 | f111 = f0111*(1d0-yd) + f1111*yd |
---|
| 893 | |
---|
| 894 | f00 = f000*(1d0-zd)+f100*zd |
---|
| 895 | f10 = f010*(1d0-zd)+f110*zd |
---|
| 896 | f01 = f001*(1d0-zd)+f101*zd |
---|
| 897 | f11 = f011*(1d0-zd)+f111*zd |
---|
| 898 | |
---|
| 899 | f0 = f00 *(1d0-td)+f10*td |
---|
| 900 | f1 = f01 *(1d0-td)+f11*td |
---|
| 901 | |
---|
| 902 | res(j) = f0*(1d0-wd)+f1*wd |
---|
| 903 | ENDDO |
---|
| 904 | ret = .true. |
---|
| 905 | END FUNCTION hdcdl5d_gen_ve |
---|
| 906 | |
---|
| 907 | END MODULE LINTGEN |
---|