File: RELATIVE:/../../../mfix.git/model/des/drag_gs_des0.f

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