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