1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: CALC_DES_DRAG_GS C 4 ! Purpose: This subroutine is only called from the DISCRETE side. C 5 ! It performs the following functions: C 6 ! - If interpolated, then it calculates the fluid-particle C 7 ! drag coefficient (F_GP) based on the particle velocity and C 8 ! interpolated fluid velocity. This F_GP is used to calculate C 9 ! the fluid-solids drag on the particle. C 10 ! - The total contact force on the particle is then updated to C 11 ! include the gas-solids drag force and gas pressure force C 12 ! C 13 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 14 SUBROUTINE DRAG_GS_DES0 15 16 !----------------------------------------------- 17 ! Modules 18 !----------------------------------------------- 19 USE param 20 USE param1 21 USE constant 22 USE physprop 23 USE fldvar 24 USE run 25 USE geometry 26 USE indices 27 USE bc 28 USE compar 29 USE sendrecv 30 USE discretelement 31 USE drag 32 USE interpolation 33 use desmpi 34 USE cutcell 35 USE mfix_pic 36 USE functions 37 IMPLICIT NONE 38 !----------------------------------------------- 39 ! Local variables 40 !----------------------------------------------- 41 ! general i, j, k indices 42 INTEGER :: I, J, K, IJK, cur_ijk 43 INTEGER :: II, JJ, KK 44 INTEGER :: IPJK, IJPK, IJKP,& 45 IPJPK, IPJKP, IJPKP, IPJPKP 46 ! indices used for interpolation stencil (unclear why IE, JN, KTP are 47 ! needed) 48 INTEGER :: IW, IE, JS, JN, KB, KTP 49 ! i,j,k indices of the fluid cell the particle resides in minus 1 50 ! (e.g., shifted 1 in west, south, bottom direction) 51 INTEGER, DIMENSION(3) :: PCELL 52 ! order of interpolation set in the call to set_interpolation_scheme 53 ! unless it is re/set later through the call to set_interpolation_stencil 54 INTEGER :: ONEW 55 ! index of solid phase that particle NP belongs to 56 INTEGER :: M 57 ! particle number index, used for looping 58 INTEGER :: NP, nindx 59 ! constant whose value depends on dimension of system 60 ! avg_factor=0.125 (in 3D) or =0.25 (in 2D) 61 ! avg_factor=0.250 (in 3D) or =0.50 (in 2D) 62 DOUBLE PRECISION :: AVG_FACTOR 63 ! for error messages 64 INTEGER :: IER 65 66 !Handan Liu added temporary variables on April 20 2012 67 DOUBLE PRECISION, DIMENSION(2,2,2,3) :: gst_tmp,vst_tmp 68 DOUBLE PRECISION, DIMENSION(2,2,2) :: weight_ft 69 DOUBLE PRECISION :: velfp(3), desposnew(3) 70 DOUBLE PRECISION :: D_FORCE(3) 71 DOUBLE PRECISION, DIMENSION(3) :: VEL_NEW 72 73 74 ! INTERPOLATED fluid-solids drag (the rest of this routine): 75 ! Calculate the gas solids drag coefficient using the particle 76 ! velocity and the fluid velocity interpolated to particle 77 ! position. 78 !----------------------------------------------------------------->>> 79 ! avg_factor=0.25 (in 3D) or =0.5 (in 2D) 80 AVG_FACTOR = merge(0.50d0, 0.25d0, NO_K) 81 82 ! sets several quantities including interp_scheme, scheme, and 83 ! order and allocates arrays necessary for interpolation 84 call set_interpolation_scheme(2) 85 86 ! There is some issue associated to gstencil, vstencil which are 87 ! allocatable variables 88 89 !$omp parallel do default(none) & 90 !$omp shared(ijkstart3,ijkend3,pinc,i_of,j_of,k_of,no_k,interp_scheme, & 91 !$omp funijk_map_c,xe,yn,dz,zt,avg_factor,do_k,pic,pea,des_pos_new, & 92 !$omp des_vel_new, mppic, mppic_pdrag_implicit,p_force, & 93 !$omp u_g,v_g,w_g,model_b,pvol,fc,f_gp,ep_g) & 94 !$omp private(ijk, i, j, k, pcell, iw, ie, js, jn, kb, ktp, & 95 !$omp onew, ii, jj, kk,cur_ijk, ipjk, ijpk, ipjpk, & 96 !$omp gst_tmp, vst_tmp, velfp, desposnew, ijpkp, ipjkp, & 97 !$omp ipjpkp,ijkp,nindx,np,weight_ft,d_force, vel_new) 98 DO ijk = ijkstart3,ijkend3 99 if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle 100 i = i_of(ijk) 101 j = j_of(ijk) 102 k = k_of(ijk) 103 104 ! generally a particle may not exist in a ghost cell. however, if the 105 ! particle is adjacent to the west, south or bottom boundary, then pcell 106 ! may be assigned indices of a ghost cell which will be passed to 107 ! set_interpolation_stencil 108 pcell(1) = i-1 109 pcell(2) = j-1 110 pcell(3) = merge(1, k-1, NO_K) 111 112 ! setup the stencil based on the order of interpolation and factoring in 113 ! whether the system has any periodic boundaries. sets onew to order. 114 call set_interpolation_stencil(pcell,iw,ie,js,jn,kb,& 115 ktp,interp_scheme,dimn,ordernew = onew) 116 117 ! Compute velocity at grid nodes and set the geometric stencil 118 DO k = 1,merge(1, ONEW, NO_K) 119 DO j = 1,onew 120 DO i = 1,onew 121 ii = iw + i-1 122 jj = js + j-1 123 kk = kb + k-1 124 cur_ijk = funijk_map_c(ii,jj,kk) 125 ipjk = funijk_map_c(ii+1,jj,kk) 126 ijpk = funijk_map_c(ii,jj+1,kk) 127 ipjpk = funijk_map_c(ii+1,jj+1,kk) 128 129 gst_tmp(i,j,k,1) = xe(ii) 130 gst_tmp(i,j,k,2) = yn(jj) 131 gst_tmp(i,j,k,3) = merge(DZ(1), zt(kk), NO_K) 132 vst_tmp(i,j,k,1) = avg_factor*(u_g(cur_ijk)+u_g(ijpk)) 133 vst_tmp(i,j,k,2) = avg_factor*(v_g(cur_ijk)+v_g(ipjk)) 134 135 if(DO_K) then 136 ijpkp = funijk_map_c(ii,jj+1,kk+1) 137 ipjkp = funijk_map_c(ii+1,jj,kk+1) 138 ipjpkp = funijk_map_c(ii+1,jj+1,kk+1) 139 ijkp = funijk_map_c(ii,jj,kk+1) 140 vst_tmp(i,j,k,1) = vst_tmp(i,j,k,1) + avg_factor*(u_g(ijkp) + u_g(ijpkp)) 141 vst_tmp(i,j,k,2) = vst_tmp(i,j,k,2) + avg_factor*(v_g(ijkp) + v_g(ipjkp)) 142 vst_tmp(i,j,k,3) = avg_factor*(w_g(cur_ijk)+& 143 w_g(ijpk)+w_g(ipjk)+w_g(ipjpk)) 144 else 145 vst_tmp(i,j,k,3) = 0.d0 146 ENDIF 147 ENDDO 148 ENDDO 149 ENDDO 150 ! loop through particles in the cell 151 ! interpolate the fluid velocity (VELFP) to the particle's position. 152 DO nindx = 1,PINC(IJK) 153 NP = PIC(ijk)%p(nindx) 154 ! skipping indices that do not represent particles and ghost particles 155 if(.not.pea(np,1)) cycle 156 if(pea(np,4)) cycle 157 158 desposnew(:) = des_pos_new(:,np) 159 call DRAG_INTERPOLATION(gst_tmp,vst_tmp,desposnew,velfp,weight_ft) 160 161 ! Calculate the particle centered drag coefficient (F_GP) using the 162 ! particle velocity and the interpolated gas velocity. Note F_GP 163 ! obtained from des_drag_gp subroutine is given as: 164 ! f_gp=beta*vol_p/eps where vol_p is the particle volume. 165 ! The drag force on each particle is equal to: 166 ! beta(u_g-u_s)*vol_p/eps. 167 ! Therefore, the drag force = f_gp*(u_g - u_s) 168 VEL_NEW(:) = DES_VEL_NEW(:,NP) 169 CALL DES_DRAG_GP(NP, VEL_NEW, VELFP, EP_G(IJK)) 170 171 ! Calculate the gas-solids drag force on the particle 172 IF(MPPIC .AND. MPPIC_PDRAG_IMPLICIT) THEN 173 ! implicit treatment of the drag term for mppic 174 D_FORCE(1:3) = F_GP(NP)*(VELFP) 175 ELSE 176 ! default case 177 D_FORCE(1:3) = F_GP(NP)*(VELFP-VEL_NEW) 178 ENDIF 179 180 ! Update the contact forces (FC) on the particle to include gas 181 ! pressure and gas-solids drag 182 FC(:3,NP) = FC(:3,NP) + D_FORCE(:3) 183 184 IF(.NOT.MODEL_B) THEN 185 ! P_force is evaluated as -dp/dx 186 FC(:3,NP) = FC(:3,NP) + P_FORCE(:,IJK)*PVOL(NP) 187 ENDIF 188 ENDDO ! end do (nindx = 1,pinc(ijk)) 189 190 ENDDO ! end do (ijk=ijkstart3,ijkend3) 191 !$omp end parallel do 192 193 194 RETURN 195 END SUBROUTINE DRAG_GS_DES0 196 197 198 199 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 200 ! C 201 ! Subroutine: DES_DRAG_GS C 202 ! Purpose: This subroutine is only called from the CONTINUUM side. C 203 ! It performs the following functions: C 204 ! - If interpolated then, it calculates the fluid-particle C 205 ! drag coefficient (F_GP) based on the particle velocity and C 206 ! interpolated fluid velocity. It then determines the C 207 ! the contributions of fluid-particle drag to the center C 208 ! coefficient of the A matrix and the b (source) vector in the C 209 ! matrix equation (A*VELFP=b) equation for the fluid phase C 210 ! x, y and z momentum balances using F_GP. C 211 ! C 212 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 213 SUBROUTINE DRAG_GS_GAS0 214 215 !----------------------------------------------- 216 ! Modules 217 !----------------------------------------------- 218 USE param 219 USE param1 220 USE constant 221 USE physprop 222 USE fldvar 223 USE run 224 USE geometry 225 USE indices 226 USE bc 227 USE compar 228 USE sendrecv 229 USE discretelement 230 USE drag 231 USE interpolation 232 use desmpi 233 USE cutcell 234 USE mfix_pic 235 USE functions 236 use mpi_node_des, only: des_addnodevalues 237 238 IMPLICIT NONE 239 !----------------------------------------------- 240 ! Local variables 241 !----------------------------------------------- 242 ! local variable used for debugging 243 LOGICAL :: FOCUS 244 ! general i, j, k indices 245 INTEGER :: I, J, K, IJK, cur_ijk 246 INTEGER :: II, JJ, KK 247 INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,& 248 IPJPK, IPJKP, IJPKP, IPJPKP, & 249 IMJMK, IMJKM, IJMKM, IMJMKM 250 ! indices used for interpolation stencil (unclear why IE, JN, KTP are 251 ! needed) 252 INTEGER :: IW, IE, JS, JN, KB, KTP 253 ! i,j,k indices of the fluid cell the particle resides in minus 1 254 ! (e.g., shifted 1 in west, south, bottom direction) 255 INTEGER, DIMENSION(3) :: PCELL 256 ! order of interpolation set in the call to set_interpolation_scheme 257 ! unless it is re/set later through the call to set_interpolation_stencil 258 INTEGER :: ONEW 259 ! index of solid phase that particle NP belongs to 260 INTEGER :: M 261 ! particle number index, used for looping 262 INTEGER :: NP, nindx 263 ! one over the volume of fluid cell 264 DOUBLE PRECISION :: OVOL 265 ! volume of fluid cell particle resides in 266 DOUBLE PRECISION :: VCELL 267 ! constant whose value depends on dimension of system 268 ! avg_factor=0.125 (in 3D) or =0.25 (in 2D) 269 ! avg_factor=0.250 (in 3D) or =0.50 (in 2D) 270 DOUBLE PRECISION :: AVG_FACTOR 271 ! Statistical weight of the particle. Equal to one for DEM 272 DOUBLE PRECISION :: WTP 273 !Handan Liu added temporary variables on April 20 2012 274 DOUBLE PRECISION, DIMENSION(2,2,2,3) :: gst_tmp,vst_tmp 275 DOUBLE PRECISION, DIMENSION(2,2,2) :: weight_ft 276 DOUBLE PRECISION :: velfp(3), desposnew(3) 277 DOUBLE PRECISION, DIMENSION(3) :: VEL_NEW 278 279 !----------------------------------------------- 280 281 282 283 ! INTERPOLATED fluid-solids drag (the rest of this routine): 284 ! Calculate the fluid solids drag coefficient using the particle 285 ! velocity and the fluid velocity interpolated to particle 286 ! position. 287 !----------------------------------------------------------------->>> 288 ! initializations 289 drag_am = ZERO 290 drag_bm = ZERO 291 292 ! avg_factor=0.25 (in 3D) or =0.5 (in 2D) 293 AVG_FACTOR = merge(0.50d0, 0.25d0, NO_K) 294 295 ! sets several quantities including interp_scheme, scheme, and 296 ! order and allocates arrays necessary for interpolation 297 call set_interpolation_scheme(2) 298 ! There is some issue associated to gstencil, vstencil which are 299 ! allocatable variables 300 301 302 !!!$omp parallel default(shared) & 303 !!!$omp private(ijk,i,j,k,pcell,iw,ie,js,jn,kb,ktp,onew, & 304 !!!$omp ii,jj,kk,cur_ijk,ipjk,ijpk,ipjpk, & 305 !!!$omp gst_tmp,vst_tmp,velfp,desposnew,ijpkp,ipjkp, & 306 !!!$omp ipjpkp,ijkp,nindx,focus,np,wtp,m,weight_ft, & 307 !!!$omp vcell,ovol) 308 !!!$omp do reduction(+:drag_am) reduction(+:drag_bm) 309 DO IJK = IJKSTART3,IJKEND3 310 IF(.NOT.FLUID_AT(IJK) .OR. PINC(IJK)==0) cycle 311 i = i_of(ijk) 312 j = j_of(ijk) 313 k = k_of(ijk) 314 315 ! generally a particle may not exist in a ghost cell. however, if the 316 ! particle is adjacent to the west, south or bottom boundary, then pcell 317 ! may be assigned indices of a ghost cell which will be passed to 318 ! set_interpolation_stencil 319 pcell(1) = i-1 320 pcell(2) = j-1 321 pcell(3) = merge(1, k-1, NO_K) 322 323 ! setup the stencil based on the order of interpolation and factoring in 324 ! whether the system has any periodic boundaries. sets onew to order. 325 call set_interpolation_stencil(pcell,iw,ie,js,jn,kb,& 326 ktp,interp_scheme,dimn,ordernew = onew) 327 328 ! Compute velocity at grid nodes and set the geometric stencil 329 DO k = 1, merge(1, ONEW, NO_K) 330 DO j = 1,onew 331 DO i = 1,onew 332 ii = iw + i-1 333 jj = js + j-1 334 kk = kb + k-1 335 cur_ijk = funijk_map_c(ii,jj,kk) 336 ipjk = funijk_map_c(ii+1,jj,kk) 337 ijpk = funijk_map_c(ii,jj+1,kk) 338 ipjpk = funijk_map_c(ii+1,jj+1,kk) 339 GST_TMP(I,J,K,1) = XE(II) 340 GST_TMP(I,J,K,2) = YN(JJ) 341 GST_TMP(I,J,K,3) = merge(DZ(1), ZT(KK), NO_K) 342 VST_TMP(I,J,K,1) = AVG_FACTOR*(U_G(CUR_IJK)+U_G(IJPK)) 343 VST_TMP(I,J,K,2) = AVG_FACTOR*(V_G(CUR_IJK)+V_G(IPJK)) 344 345 IF(DO_K) THEN 346 IJPKP = FUNIJK_MAP_C(II,JJ+1,KK+1) 347 IPJKP = FUNIJK_MAP_C(II+1,JJ,KK+1) 348 IPJPKP = FUNIJK_MAP_C(II+1,JJ+1,KK+1) 349 IJKP = FUNIJK_MAP_C(II,JJ,KK+1) 350 VST_TMP(I,J,K,1) = VST_TMP(I,J,K,1) + & 351 AVG_FACTOR*(U_G(IJKP) + U_G(IJPKP)) 352 353 VST_TMP(I,J,K,2) = VST_TMP(I,J,K,2) + & 354 AVG_FACTOR*(V_G(IJKP) + V_G(IPJKP)) 355 356 VST_TMP(I,J,K,3) = AVG_FACTOR*(W_G(CUR_IJK)+& 357 W_G(IJPK)+W_G(IPJK)+W_G(IPJPK)) 358 ELSE 359 VST_TMP(I,J,K,3) = 0.D0 360 ENDIF 361 ENDDO 362 ENDDO 363 ENDDO 364 ! loop through particles in the cell 365 ! interpolate the fluid velocity (VELFP) to the particle's position. 366 DO nindx = 1,PINC(IJK) 367 NP = PIC(ijk)%p(nindx) 368 ! skipping indices that do not represent particles and ghost particles 369 if(.not.pea(np,1)) cycle 370 if(pea(np,4)) cycle 371 desposnew(:) = des_pos_new(:,np) 372 call DRAG_INTERPOLATION(gst_tmp,vst_tmp,desposnew,velfp,weight_ft) 373 ! 374 ! Calculate the particle centered drag coefficient (F_GP) using the 375 ! particle velocity and the interpolated gas velocity. Note F_GP 376 ! obtained from des_drag_gp subroutine is given as: 377 ! f_gp=beta*vol_p/eps where vol_p is the particle volume. 378 ! The drag force on each particle is equal to: 379 ! beta(u_g-u_s)*vol_p/eps. 380 ! Therefore, the drag force = f_gp*(u_g - u_s) 381 VEL_NEW(:) = DES_VEL_NEW(:,NP) 382 CALL DES_DRAG_GP(NP, VEL_NEW, VELFP, EP_G(IJK)) 383 !-----------------------------------------------------------------<<< 384 ! Calculate the corresponding gas solids drag force that is used in 385 ! the gas phase momentum balances. 386 !----------------------------------------------------------------->>> 387 focus = .false. 388 M = pijk(np,5) 389 WTP = ONE 390 IF(MPPIC) WTP = DES_STAT_WT(NP) 391 392 DO k = 1, merge(1, ONEW, NO_K) 393 DO j = 1, onew 394 DO i = 1, onew 395 ! shift loop index to new variables for manipulation 396 ii = iw + i-1 397 jj = js + j-1 398 kk = kb + k-1 399 ! The interpolation is done using node. so one should use consistent 400 ! numbering system. in the current version imap_c is used instead of 401 ! ip_of or im_of 402 cur_ijk = funijk_map_c(ii, jj, kk) 403 404 ! Replacing the volume of cell to volume at the node 405 vcell = des_vol_node(cur_ijk) 406 ovol = one/vcell 407 408 !!!$omp critical 409 drag_am(cur_ijk) = drag_am(cur_ijk) + & 410 f_gp(np)*weight_ft(i,j,k)*ovol*wtp 411 412 drag_bm(cur_ijk,1:3) = & 413 drag_bm(cur_ijk,1:3) + & 414 f_gp(np) * vel_new(1:3) * & 415 weight_ft(i,j,k)*ovol*wtp 416 !!!$omp end critical 417 ENDDO 418 ENDDO 419 ENDDO 420 ENDDO ! end do (nindx = 1,pinc(ijk)) 421 422 ENDDO ! end do (ijk=ijkstart3,ijkend3) 423 !!!$omp end parallel 424 425 426 ! At the interface drag_am and drag_bm have to be added 427 ! send recv will be called and the node values will be added 428 ! at the junction. drag_am are drag_bm are altered by the 429 ! routine when periodic boundaries are invoked. so both 430 ! quantities are needed at the time of this call. 431 call des_addnodevalues 432 !-----------------------------------------------------------------<<< 433 ! Calculate/update the cell centered drag coefficient F_GDS for use 434 ! in the pressure correction equation 435 !----------------------------------------------------------------->>> 436 ! avg_factor=0.125 (in 3D) or =0.25 (in 2D) 437 AVG_FACTOR = merge(0.25d0, 0.125D0, NO_K) 438 439 !!$omp parallel do default(shared) & 440 !!$omp private(ijk,i,j,k,imjk,ijmk,imjmk,ijkm,imjkm,ijmkm, & 441 !!$omp imjmkm) & 442 !!$omp schedule (guided,20) 443 DO ijk = ijkstart3, ijkend3 444 IF(FLUID_AT(IJK)) THEN 445 446 i = i_of(ijk) 447 j = j_of(ijk) 448 k = k_of(ijk) 449 if (i.lt.istart2 .or. i.gt.iend2) cycle 450 if (j.lt.jstart2 .or. j.gt.jend2) cycle 451 if (k.lt.kstart2 .or. k.gt.kend2) cycle 452 imjk = funijk_map_c(i-1,j,k) 453 ijmk = funijk_map_c(i,j-1,k) 454 imjmk = funijk_map_c(i-1,j-1,k) 455 456 f_gds(ijk) = avg_factor*& 457 (drag_am(ijk) + drag_am(ijmk) +& 458 drag_am(imjmk) + drag_am(imjk)) 459 460 IF(DO_K) THEN 461 ijkm = funijk_map_c(i,j,k-1) 462 imjkm = funijk_map_c(i-1,j,k-1) 463 ijmkm = funijk_map_c(i,j-1,k-1) 464 imjmkm = funijk_map_c(i-1,j-1,k-1) 465 f_gds(ijk) = f_gds(ijk) + avg_factor*& 466 (drag_am(ijkm) + drag_am(ijmkm) +& 467 drag_am(imjmkm)+drag_am(imjkm) ) 468 ENDIF ! end if 469 ELSE ! else branch of if (fluid_at(ijk)) 470 F_GDS(IJK) = ZERO 471 ENDIF ! end if/else (fluid_at(ijk)) 472 473 ENDDO ! end do loop (ijk=ijkstart3,ijkend3) 474 !!$omp end parallel do 475 476 RETURN 477 END SUBROUTINE DRAG_GS_GAS0 478 479 480 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 481 ! Subroutine: DRAG_INTERPOLATION C 482 ! Purpose: DES - Calculate the fluid velocity interpolated at the C 483 ! particle's location and weights. Replace 'interpolator' C 484 ! interface for OpenMP implementation. C 485 ! C 486 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 487 488 SUBROUTINE DRAG_INTERPOLATION(GSTEN,VSTEN,DESPOS,VELFP,WEIGHTFACTOR) 489 490 use geometry, only: NO_K 491 492 IMPLICIT NONE 493 494 !----------------------------------------------- 495 ! Local Variables 496 !----------------------------------------------- 497 DOUBLE PRECISION, DIMENSION(2,2,2,3), INTENT(IN):: GSTEN 498 DOUBLE PRECISION, DIMENSION(2,2,2,3), INTENT(IN):: VSTEN 499 DOUBLE PRECISION, DIMENSION(3), INTENT(IN):: DESPOS 500 DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: VELFP 501 DOUBLE PRECISION, DIMENSION(2,2,2), INTENT(OUT) :: WEIGHTFACTOR 502 INTEGER :: II, JJ, KK 503 504 DOUBLE PRECISION, DIMENSION(2) :: XXVAL, YYVAL, ZZVAL 505 DOUBLE PRECISION :: DXX, DYY, DZZ 506 DOUBLE PRECISION, DIMENSION(3) :: ZETAA 507 508 DXX = GSTEN(2,1,1,1) - GSTEN(1,1,1,1) 509 DYY = GSTEN(1,2,1,2) - GSTEN(1,1,1,2) 510 511 ZETAA(1:2) = DESPOS(1:2) - GSTEN(1,1,1,1:2) 512 513 ZETAA(1) = ZETAA(1)/DXX 514 ZETAA(2) = ZETAA(2)/DYY 515 516 XXVAL(1)=1-ZETAA(1) 517 YYVAL(1)=1-ZETAA(2) 518 XXVAL(2)=ZETAA(1) 519 YYVAL(2)=ZETAA(2) 520 521 VELFP(:) = 0.D0 522 523 IF(NO_K) THEN 524 DO JJ=1,2 525 DO II=1,2 526 WEIGHTFACTOR(II,JJ,1) = XXVAL(II)*YYVAL(JJ) 527 VELFP(1:2) = VELFP(1:2) + VSTEN(II,JJ,1,1:2)*WEIGHTFACTOR(II,JJ,1) 528 ENDDO 529 ENDDO 530 ELSE 531 DZZ = GSTEN(1,1,2,3) - GSTEN(1,1,1,3) 532 ZETAA(3) = DESPOS(3) - GSTEN(1,1,1,3) 533 ZETAA(3) = ZETAA(3)/DZZ 534 ZZVAL(1)=1-ZETAA(3) 535 ZZVAL(2)=ZETAA(3) 536 DO KK=1,2 537 DO JJ=1,2 538 DO II=1,2 539 WEIGHTFACTOR(II,JJ,KK) = XXVAL(II)*YYVAL(JJ)*ZZVAL(KK) 540 VELFP(1:3) = VELFP(1:3) + VSTEN(II,JJ,KK,1:3)*WEIGHTFACTOR(II,JJ,KK) 541 ENDDO 542 ENDDO 543 ENDDO 544 ENDIF 545 546 END SUBROUTINE DRAG_INTERPOLATION 547