File: N:\mfix\model\des\interpolation_mod.f

1     ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: INTERPOLATION                                          C
4     !  Purpose: Contains number of generic interfaces that are used if
5     !           DES_INTERP_ON is TRUE
6     !                                                                      C
7     !                                                                      C
8     !  Author: Chidambaram Narayanan and Rahul Garg       Date: 01-Aug-07  C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10     
11     MODULE interpolation
12     
13       USE constant
14       USE discretelement
15       USE param1, only: half, one, zero
16       IMPLICIT NONE
17     
18       PRIVATE
19     
20       PUBLIC:: set_interpolation_stencil,  set_interpolation_scheme, set_interpolation_pstencil, array_dot_product
21       PUBLIC:: INTERPOLATE_CC
22     
23     
24       PUBLIC:: interpolator
25       INTERFACE interpolator
26          MODULE PROCEDURE interp_oned_scalar
27          MODULE PROCEDURE interp_oned_vector
28          MODULE PROCEDURE interp_threed_scalar
29          MODULE PROCEDURE interp_threed_vector
30          MODULE PROCEDURE interp_twod_scalar
31          MODULE PROCEDURE interp_twod_vector
32          MODULE PROCEDURE calc_weightderiv_threed
33          MODULE PROCEDURE calc_weightderiv_twod
34       END INTERFACE
35     
36       INTERFACE array_dot_product
37           Module procedure array_dot_product_2d
38           Module procedure array_dot_product_3d
39       END INTERFACE
40     
41       SAVE
42     
43     ! the interpolator subroutine essentially returns the value of interp_scl/interp_vec
44     ! which is the interpolated value of the quantity passed through 'scalar/vector'.
45     ! the value is interpolated based on grid (x,y,z position of grid)  passed
46     ! through coor at the x,y,z position of the particle passed through ppos
47     ! the interpolation order is specified by order and the interpolation scheme is
48     ! specified by isch/fun
49     
50     ! interp_oned_scalar  (coor,scalar,ppos,interp_scl,order, isch,weight_pointer)
51     !    REAL coor(:), scalar(:), PPOS, INTERP_SCL
52     !    INTEGER ORDER; CHARACTER ISCH; REAL WEIGHT_POINTER(:)
53     ! interp_oned_vector  (coor,vector,ppos,interp_vec,order, fun, weight_pointer)
54     !    REAL coor(:), vector(:,:), PPOS, INTERP_VEC(:)
55     !    INTEGER ORDER; CHARACTER FUN; REAL WEIGHT_POINTER(:)
56     ! interp_twod_scalar  (coor,scalar,ppos,interp_scl,order, isch,weight_pointer)
57     !    REAL coor(:,:,:), scalar(:,:), PPOS(2), INTERP_SCL
58     !    INTEGER ORDER; CHARACTER ISCH; REAL WEIGHT_POINTER(:,:,:)
59     ! interp_twod_vector  (coor,vector,ppos,interp_vec,order, fun, weight_pointer)
60     !    REAL coor(:,:,:), vector(:,:,:), PPOS(2), INTERP_VEC(:)
61     !    INTEGER ORDER; CHARACTER FUN; REAL WEIGHT_POINTER(:,:,:)
62     ! interp_threed_scalar(coor,scalar,ppos,interp_scl,order, isch,weight_pointer)
63     !    REAL coor(:,:,:,:), scalar(:,:,:), PPOS(3), INTERP_SCL
64     !    INTEGER ORDER; CHARACTER ISCH; REAL WEIGHT_POINTER(:,:,:)
65     ! interp_threed_vector(coor,vector,ppos,interp_vec,order, fun, weight_pointer)
66     !    REAL coor(:,:,:,:), vector(:,:,:,:), PPOS(3), INTERP_VEC(:)
67     !    INTEGER ORDER; CHARACTER FUN; REAL WEIGHT_POINTER(:,:,:)
68     
69       INTEGER, PARAMETER :: prcn=8
70       INTEGER, PARAMETER :: iprec=8
71       !double precision, Parameter:: zero = 0.0_iprec
72       !double precision, Parameter:: one  = 1.0_iprec
73       DOUBLE PRECISION, PARAMETER  :: two = 2.0_iprec
74       DOUBLE PRECISION, PARAMETER  :: three = 3.0_iprec
75       DOUBLE PRECISION, PARAMETER  :: four = 4.0_iprec
76       DOUBLE PRECISION, PARAMETER :: six = 6.0_iprec
77     
78       !double precision, Parameter:: half = 0.5_iprec
79       DOUBLE PRECISION, PARAMETER  :: fourth = 0.25_iprec
80     
81       INTEGER, PARAMETER :: maxorder=6
82       DOUBLE PRECISION, DIMENSION(maxorder), TARGET :: xval, yval, zval
83       DOUBLE PRECISION, DIMENSION(maxorder-1) :: dx, dy, dz
84       DOUBLE PRECISION, DIMENSION(maxorder,maxorder,maxorder), TARGET :: weights
85     
86     
87     
88      CONTAINS
89     
90             !----------
91       ! A . B where A and B are arrays of the same dimensions
92       !----------
93       !----------
94       FUNCTION array_dot_product_3d(A,B)  !3D Arrays
95         Real(prcn):: array_dot_product_3d
96     
97         Integer:: j, k, s1, s2, s3
98         Real(prcn), Dimension(:,:,:), Intent(in):: A, B
99     
100         s1 = SIZE(A,1)
101         s2 = SIZE(A,2)
102         s3 = SIZE(A,3)
103     
104         array_dot_product_3d = zero
105         DO k = 1,s3
106            DO j = 1,s2
107               array_dot_product_3d = array_dot_product_3d + dot_product(A(:,j,k), B(1:s1,j,k))
108            ENDDO
109         ENDDO
110       END FUNCTION array_dot_product_3d
111     
112       FUNCTION array_dot_product_2d(A,B)  !2D Arrays
113         Real(prcn):: array_dot_product_2d
114     
115         Integer:: j, s2
116         Real(prcn), Dimension(:,:), Intent(in):: A, B
117     
118         s2 = SIZE(A,2)
119     
120         array_dot_product_2d = zero
121         DO j = 1,s2
122            array_dot_product_2d = array_dot_product_2d + dot_product(A(:,j), B(:,j))
123         ENDDO
124       END FUNCTION array_dot_product_2d
125     
126     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
127       SUBROUTINE set_interpolation_stencil(PC, IW, IE, JS, JN, KB, KTP,&
128           isch,dimprob,ordernew)
129     
130     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
131            !USE discretelement, ONLY : order,ob2l,ob2r, des_periodic_walls_x,y,z
132     
133         USE geometry
134     
135         IMPLICIT NONE
136     !-----------------------------------------------
137     ! Local variables
138     !-----------------------------------------------
139         INTEGER, DIMENSION(3), INTENT(in):: pc   ! i,j,k indices of particle - 1
140         INTEGER, INTENT(out):: IW, IE, JS, JN, KB, KTP
141         CHARACTER(LEN=5), INTENT(in) :: isch   ! interpolation scheme
142         INTEGER, INTENT(in) :: dimprob   ! dimension of system = DIMN
143         INTEGER, OPTIONAL :: ordernew   ! interpolation order
144     
145         INTEGER :: im, jm, km   ! local variables assigned maximum i,j,k fluid cell indices
146         INTEGER :: ob2rtmp, ob2ltmp, ordertemp, ordernewtmp
147     !-----------------------------------------------
148     
149     ! reassigning ob2l and ob2r ?
150         ob2l = (order+1)/2
151         ob2r = order/2
152     ! local variables for maximum fluid cell indices
153         im = imax1
154         jm = jmax1
155         km = kmax1
156     
157         SELECT CASE(isch)
158         CASE('csi')
159     
160            ob2rtmp = ob2r
161            ob2ltmp = ob2l
162            ordertemp = order
163     
164     ! lowest IW will be assigned is 1 (ghost cell)
165            IW = MAX(1 ,pc(1) - (ob2l - 1))
166     ! highest IE will be assigned is maximum fluid cell index
167            IE = MIN(im,pc(1) + ob2r)
168     ! if IW is west ghost cell and/or IE is maximum fluid cell
169     ! reassign IW and/or IE accordingly
170            IF(.NOT.DES_PERIODIC_WALLS_X) THEN
171               IF (IW.EQ.1 ) IE = IW + order - 1
172               IF (IE.EQ.im) IW = IE - order + 1
173            ELSE
174               IF (IW.EQ.1 ) IW = IE - order + 1
175               IF (IE.EQ.im) IE = IW + order - 1
176            ENDIF
177     
178            JS = MAX(1 ,pc(2) - (ob2l - 1)) !non-periodic
179            JN = MIN(jm,pc(2) + ob2r)
180            IF(.NOT.DES_PERIODIC_WALLS_Y) THEN
181               IF (JS.EQ.1 ) JN = JS + order - 1
182               IF (JN.EQ.jm) JS = JN - order + 1
183            ELSE
184               IF (JS.EQ.1 ) JS = JN - order + 1
185               IF (JN.EQ.jm) JN = JS + order - 1
186            ENDIF
187     
188            KB = MAX(1 ,pc(3) - (ob2l - 1)) !non-periodic
189            KTP = MIN(km,pc(3) + ob2r)
190            IF(.NOT.DES_PERIODIC_WALLS_Z) THEN
191               IF (KB.EQ.1 ) KTP = KB + order - 1
192               IF (KTP.EQ.km) KB = KTP - order + 1
193            ELSE
194               IF (KB.EQ.1 ) KB = KTP - order + 1
195               IF (KTP.EQ.km) KTP = KB + order - 1
196            ENDIF
197     
198            ob2r =  ob2rtmp
199            ob2l = ob2ltmp
200            ordernewtmp = order
201            order = ordertemp !reset the order
202     
203         CASE('lpi')
204     
205            IW = MAX(1 ,pc(1) - (ob2l - 1)) !non-periodic
206            IE = MIN(im,pc(1) + ob2r)
207            IF(.NOT.DES_PERIODIC_WALLS_X) THEN
208               IF (IW.EQ.1 ) IE = IW + order - 1
209               IF (IE.EQ.im) IW = IE - order + 1
210            ELSE
211               IF (IW.EQ.1 ) IW = IE - order + 1
212               IF (IE.EQ.im) IE = IW + order - 1
213            ENDIF
214     
215            JS = MAX(1 ,pc(2) - (ob2l - 1)) !non-periodic
216            JN = MIN(jm,pc(2) + ob2r)
217            IF(.NOT.DES_PERIODIC_WALLS_Y) THEN
218               IF (JS.EQ.1 ) JN = JS + order - 1
219               IF (JN.EQ.jm) JS = JN - order + 1
220            ELSE
221               IF (JS.EQ.1 ) JS = JN - order + 1
222               IF (JN.EQ.jm) JN = JS + order - 1
223            ENDIF
224     
225            KB = MAX(1 ,pc(3) - (ob2l - 1)) !non-periodic
226            KTP = MIN(km,pc(3) + ob2r)
227            IF(.NOT.DES_PERIODIC_WALLS_Z) THEN
228               IF (KB.EQ.1 ) KTP = KB + order - 1
229               IF (KTP.EQ.km) KB = KTP - order + 1
230            ELSE
231               IF (KB.EQ.1 ) KB = KTP - order + 1
232               IF (KTP.EQ.km) KTP = KB + order - 1
233            ENDIF
234            ordernewtmp = order
235     
236         END SELECT
237     
238         IF(dimprob == 2) THEN
239            KB = pc(3)
240            KTP = pc(3)
241         ENDIF
242     
243     ! for debugging
244         !print*, 'order in set stencil = ', order,pc(3)
245         !Print*,'IW, IE = ', pc(1), IW, IE
246     
247         IF(PRESENT(ordernew)) ordernew = ordernewtmp
248     
249       END SUBROUTINE set_interpolation_stencil
250     
251     
252     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
253     !  Subroutine: SET_INTERPOLATION_STENCIL_CC                            !
254     !                                                                      !
255     !  Purpose:                                                            !
256     !                                                                      !
257     !                                                                      !
258     !  Author: J.Musser                                   Date:            !
259     !                                                                      !
260     !  Comments:                                                           !
261     !                                                                      !
262     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
263           SUBROUTINE SET_INTERPOLATION_STENCIL_CC(NP, PC, IW, JS, KB, FOCUS)
264     
265           USE discretelement
266           USE geometry
267           USE param1
268     
269           IMPLICIT NONE
270     
271     ! I, J, K indices of the cell containing the particle.
272           INTEGER, INTENT(IN) :: PC(3)
273     ! particle number
274           INTEGER, INTENT(IN) :: NP
275     
276     ! These values indicate the I,J,K starting indices used to select the
277     ! nodes for interpolation.
278           INTEGER, INTENT(OUT) :: IW, JS, KB
279     ! flag to tell whether local debugging was called in drag_fgs
280           LOGICAL, INTENT(IN) :: FOCUS
281     
282     !-----------------------------------------------
283     ! Local variables
284     !-----------------------------------------------
285     ! indices
286           INTEGER I_SHIFT, J_SHIFT, K_SHIFT, IE, JN, KTP
287     
288     ! Determine the shift index necessary to account for the quadrant of
289     ! the fluid cell where the particle is located.
290           IF( DES_POS_NEW(NP,1) <= XE(PC(1))-HALF*DX(PC(1)))THEN
291              I_SHIFT = 0
292           ELSE
293              I_SHIFT = 1
294           ENDIF
295           IF( DES_POS_NEW(NP,2) <= YN(PC(2))-HALF*DY(PC(2)))THEN
296              J_SHIFT = 0
297           ELSE
298              J_SHIFT = 1
299           ENDIF
300           IF(DIMN == 2) THEN
301              K_SHIFT = 0
302           ELSE
303              IF( DES_POS_NEW(NP,3) <= ZT(PC(3))-HALF*DZ(PC(3)))THEN
304                 K_SHIFT = 0
305              ELSE
306                 K_SHIFT = 1
307              ENDIF
308           ENDIF
309     
310     ! Calculate the west and east I indices to be used.
311     !-----------------------------------------------------------------------
312     ! Restrict the west I index to a minimum of 1
313           IW = MAX( 1, (PC(1) + I_SHIFT-1))
314     ! Restrict the east I index to a maximum of IMAX1
315           IE = MIN( IMAX2, PC(1) + I_SHIFT)
316           IF(.NOT.DES_PERIODIC_WALLS_X) THEN
317     ! If the YZ walls are not periodic, ensure that the indices used for
318     ! interpolation stay within the domain. This may require a shift in the
319     ! indices for cells near the walls and interpolation schemes of high
320     ! order.
321              IF (IE .EQ. IMAX2) IW = IMAX2 - 1
322           ELSE
323     ! If the YZ walls are periodic, the starting index can been less than 1.
324     ! This is accounted for when filling GSTENCIL. corrected later.
325              IF (IW .EQ. 1) IW = IE - 1
326           ENDIF
327     
328     ! Calculate the south and north J indices to be used.
329     !-----------------------------------------------------------------------
330     ! Restrict the south J index to a minimum of 1
331           JS = MAX( 1, PC(2)+ J_SHIFT - 1)
332     ! Restrict the north J index to a maximum of JMAX1
333           JN = MIN( JMAX2, PC(2) + J_SHIFT)
334           IF(.NOT.DES_PERIODIC_WALLS_Y) THEN
335     ! Shift indices for non periodic boundaries if needed.
336              IF (JN .EQ. JMAX2) JS = JMAX2 - 1
337           ELSE
338     ! Shift indices for periodic boundaries if needed.
339              IF (JS .EQ. 1 ) JS = JN - 1
340           ENDIF
341     
342     ! Calculate the bottom and top K indices to be used.
343     !-----------------------------------------------------------------------
344           IF(DIMN == 2) THEN
345     ! If 2 dimension simulation, set the bottom index to 1.
346              KB = 1
347           ELSE
348     ! Restric the bottom K index to a minimum  of 1
349              KB = MAX( 1, PC(3)+K_SHIFT-1)
350     ! Restric the top K index to a maximum  of KMAX1
351              KTP = MIN( KMAX2, PC(3) + K_SHIFT)
352              IF(.NOT.DES_PERIODIC_WALLS_Z) THEN
353     ! Shift indices for non periodic boundaries if needed.
354                 IF (KTP .EQ. KMAX1) KB = KTP - 1
355              ELSE
356     ! Shift indices for periodic boundaries if needed.
357                 IF (KB .EQ. 1 ) KB = KTP - 1
358              ENDIF
359           ENDIF
360     
361           RETURN
362           END SUBROUTINE SET_INTERPOLATION_STENCIL_CC
363     
364     
365     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
366     !  SUBROUTINE: INTERPOLATE_CC                                          !
367     !                                                                      !
368     !  Purpose:                                                            !
369     !                                                                      !
370     !  Author: J.Musser                                   Date: 13-Jan-11  !
371     !                                                                      !
372     !  Comments:                                                           !
373     !                                                                      !
374     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
375           SUBROUTINE INTERPOLATE_CC(NP, INTP_IJK, INTP_WEIGHTS, FOCUS)
376     
377           USE compar
378           USE discretelement
379           USE geometry
380           USE indices
381           USE param
382           USE param1
383           USE functions
384     
385           IMPLICIT NONE
386     
387     ! Particle index number
388           INTEGER, INTENT(IN)  :: NP
389     ! Local debugging flag
390           LOGICAL, INTENT(IN)  :: FOCUS
391     
392     ! IJK values of cells. These have already been adjusted for periodic
393     ! and non-periodic boundary conditions.
394           INTEGER, INTENT(INOUT)  :: INTP_IJK(2**DIMN)
395     ! Weights associated with the IJK cells in INTP_IJK
396           DOUBLE PRECISION, INTENT(INOUT)  :: INTP_WEIGHTS(2**DIMN)
397     
398     !-----------------------------------------------
399     ! Local variables
400     !-----------------------------------------------
401     
402     ! Starting grid indices used in interpolation
403           INTEGER IW, JS, KB
404     
405     ! Weights associated with the nodes used for interpolation. These values
406     ! only need to be returned with interpolating routine is invoked from
407     ! the continuum phase.
408           DOUBLE PRECISION, POINTER :: WEIGHTS_CC (:,:,:)
409     ! indices
410           INTEGER I, J, K, II, JJ, KK
411     ! Loop counter
412           INTEGER LC
413     ! geometery and field stencils for cell-center interpolations
414           DOUBLE PRECISION STNCL_GEMTY (2, 2, MAX(1, 2*(DIMN-2)), DIMN)
415           DOUBLE PRECISION STNCL_FEILD (2, 2, MAX(1, 2*(DIMN-2)))
416     ! Generic value. In the original interpolation routine, this would be
417     ! the value that was interpolated within the code. However, this routine
418     ! only needs to return the IJK indicies of the cells and the interpolation
419     ! weights. In order to satisfy the function interface of "interpolator"
420     ! is to send it is dummy value.
421           DOUBLE PRECISION INTERP_scalar
422     ! Interpolatin scheme. It is 'hard coded' so that only a second order
423     ! Lagrange polynomial is used. Higher order isn't necessary.
424           CHARACTER(LEN=5) :: INTP_SCHM = 'LPI'
425     !-----------------------------------------------
426     
427     ! Obtain the starting cell index values for interpolation around
428     ! particle NP
429           CALL SET_INTERPOLATION_STENCIL_CC(NP, PIJK(NP,1:3), IW, JS, KB, &
430              FOCUS)
431     
432           LC = 0
433           DO K = 1, (3-DIMN)*1+(DIMN-2)*2
434              DO J = 1, 2
435                 DO I = 1, 2
436     ! increment the loop counter
437                    LC = LC + 1
438     ! Shift loop indicies to the correct values
439                    II = IW + I-1
440                    JJ = JS + J-1
441                    KK = KB + K-1
442     ! Adjust for boundary conditions.
443                    IF (DES_PERIODIC_WALLS_X) THEN
444                       IF(II .LT. IMIN1)THEN
445     ! Shift from the east side of the domain to the west side.
446                          II = IMAX2 + (II-IMIN1)
447                          STNCL_GEMTY(I,J,K,1) = XE(II) - HALF*DX(II) - XLENGTH
448                       ELSEIF(II .GT. IMAX1) THEN
449     ! Shift for the west side of the domain to the east side.
450                          II = II - IMAX
451                          STNCL_GEMTY(I,J,K,1) = XLENGTH + XE(II) - HALF*DX(II)
452                       ELSE
453     ! No shift needed. Use actual node location.
454                          STNCL_GEMTY(I,J,K,1) = XE(II) - HALF*DX(II)
455                       ENDIF
456                    ELSE
457     ! If the YZ plane is not periodic, calculate the position of the node.
458                       STNCL_GEMTY(I,J,K,1) = XE(II) - HALF*DX(II)
459     ! Due to the restrictions on non periodic 1 <= II and II <= IMAX2.
460     ! Shift the II east/west to reflect the cell value at the boundary.
461                       IF(II .EQ. 1)THEN
462                          II = IMIN1
463                       ELSEIF(II .EQ. IMAX2)THEN
464                          II = IMAX1
465                       ENDIF
466                    ENDIF
467                    IF (DES_PERIODIC_WALLS_Y) THEN
468                       IF(JJ .LT. JMIN1) THEN
469     ! Shift from the north side of the domain to the south side.
470                          JJ = JMAX2 + (JJ-JMIN1)
471                          STNCL_GEMTY(I,J,K,2) = YN(JJ) - HALF*DY(JJ) - YLENGTH
472                       ELSEIF(JJ .GT. JMAX1) THEN
473     ! Shift from the south side of the domain to the north side.
474                          JJ = JJ - JMAX
475                          STNCL_GEMTY(I,J,K,2) = YLENGTH + YN(JJ) - HALF*DY(JJ)
476                       ELSE
477     ! No shift needed. Use the nodes actual position.
478                          STNCL_GEMTY(I,J,K,2) = YN(JJ) - HALF*DY(JJ)
479                       ENDIF
480                    ELSE
481     ! If the XZ plane is not periodic, calculate the position of the node.
482                       STNCL_GEMTY(I,J,K,2) = YN(JJ) - HALF*DY(JJ)
483     ! Due to the restrictions on non periodic 1 <= II and II <= IMAX2.
484     ! Shift the II east/west to reflect the cell value at the boundary.
485                       IF(JJ .EQ. 1)THEN
486                          JJ = JMIN1
487                       ELSEIF(JJ .EQ. JMAX2)THEN
488                          JJ = JMAX1
489                       ENDIF
490                    ENDIF
491                    IF (DIMN .EQ. 3) THEN
492                       IF (DES_PERIODIC_WALLS_Z) THEN
493                          IF(KK .LT. KMIN1) THEN
494     ! Shift from the top side of the domain to the bottom side.
495                             KK = KMAX1 + (KK-KMIN1)
496                             STNCL_GEMTY(I,J,K,3) = ZT(KK) - HALF*DZ(KK) - ZLENGTH
497                          ELSEIF(KK .GT. KMAX1) THEN
498     ! Shift from the bottom side of the domain to the top side.
499                             KK = KK - KMAX
500                             STNCL_GEMTY(I,J,K,3) = ZLENGTH + ZT(KK) - HALF*DZ(KK)
501                          ELSE
502     ! No shift needed. Use actual node location.
503                             STNCL_GEMTY(I,J,K,3) = ZT(KK) - HALF*DZ(KK)
504                          ENDIF
505                       ELSE
506     ! If the XY plane is not periodic, calculate the position of the node.
507                          STNCL_GEMTY(I,J,K,3) = ZT(KK) - HALF*DZ(KK)
508     ! Due to the restrictions on non periodic 1 <= II and II <= IMAX2.
509     ! Shift the II east/west to reflect the cell value at the boundary.
510                          IF(KK .EQ. 1)THEN
511                             KK = KMIN1
512                          ELSEIF(KK .EQ. KMAX2)THEN
513                             KK = KMAX1
514                          ENDIF
515                       ENDIF
516                    ELSE
517                       KK = 1
518                    ENDIF
519     ! Store the IJK in the interp array
520                    INTP_IJK(LC) = FUNIJK(II,JJ,KK)
521     ! Assign a dummy value for the feild variable
522                    STNCL_FEILD(I,J,K) = 1.0d0
523                 ENDDO
524              ENDDO
525           ENDDO
526     
527     ! Call the interpolation routine. If the call to this routine originates
528     ! from the continuum phase, the interpolation weights are of interest.
529     ! Otherwise, the interpolated scalar value (INTERP_scalar) is desired.
530           IF(DIMN.EQ.2) THEN
531              CALL interpolator( STNCL_GEMTY(1:2, 1:2, 1, 1:DIMN), &
532                 STNCL_FEILD(1:2, 1:2, 1), DES_POS_NEW(NP,1:2), INTERP_scalar,&
533                 2, INTP_SCHM, WEIGHTS_CC )
534           ELSE
535              CALL interpolator( STNCL_GEMTY(1:2, 1:2, 1:2, 1:DIMN), &
536                 STNCL_FEILD(1:2, 1:2, 1:2), DES_POS_NEW(NP,:), INTERP_scalar,&
537                 2, INTP_SCHM, WEIGHTS_CC )
538           ENDIF
539     ! Store interpolation weights in an array for calling routine.
540           LC = 0
541           DO K = 1, (3-DIMN)*1+(DIMN-2)*2
542              DO J = 1, 2
543                 DO I = 1, 2
544     ! increment the loop counter
545                    LC = LC + 1
546                    INTP_WEIGHTS(LC) = WEIGHTS_CC(I,J,K)
547                 ENDDO
548              ENDDO
549           ENDDO
550     
551     ! Local debugging
552           IF(FOCUS .AND. DEBUG_DES)THEN
553              WRITE(*,"(3X,A,1X,A3,4X,A,I2)") &
554                 'Interpolation Scheme:','LPI','Interpolation Order:',2
555              WRITE(*,"(3X,A,I5,3X,A,I6)")'Particle: ',NP,'IJK: ',PIJK(NP,4)
556              WRITE(*,"(/5X,'|',29('-'),'|')")
557              WRITE(*,"(5X,A)")'| Fluid Cell | Interp. Weight |'
558              WRITE(*,"(5X,'|',12('-'),'|',16('-'),'|')")
559              DO LC = 1, 2**DIMN
560                 WRITE(*,"(5X,'|',2X,I8,2X,'|',2X,F13.10,2X,'|')")&
561                    INTP_IJK(LC), INTP_WEIGHTS(LC)
562              ENDDO
563              WRITE(*,"(5X,'|',29('-'),'|'//)")
564           ENDIF ! Local Debugging
565     
566           END SUBROUTINE INTERPOLATE_CC
567     
568     
569     
570     
571       SUBROUTINE set_interpolation_pstencil(pc, ib,ie,jb,je,kb,ke, isch&
572            &,dimprob, ordernew)
573            !USE discretelement, ONLY : order,ob2l,ob2r,  intx_per, inty_per, intz_per
574         USE geometry
575         USE param1, only: half
576     
577         IMPLICIT NONE
578         INTEGER, DIMENSION(3), INTENT(in):: pc
579         INTEGER :: im,jm,km
580         INTEGER, INTENT(out):: ib, ie, jb, je, kb, ke
581         INTEGER, INTENT(in) :: dimprob
582         INTEGER, OPTIONAL :: ordernew
583         CHARACTER(LEN=5), INTENT(in) :: isch
584         INTEGER :: ob2rtmp, ob2ltmp, ordertemp, ordernewtmp
585         ob2l = (order+1)/2
586         ob2r = order/2
587         im = imax1
588         jm = jmax1
589         km = kmax1
590     
591         SELECT CASE(isch)
592         CASE('csi')
593     
594            ob2rtmp = ob2r
595            ob2ltmp = ob2l
596            ordertemp = order
597     !!$       IF(order.NE.3.AND.(pc(1).EQ.1.OR.pc(1).EQ.im&
598     !!$            &-1.OR.pc(2).EQ.1.OR.pc(2).EQ.jm&
599     !!$            &-1.OR.pc(3).EQ.1.OR.pc(3).EQ.km-1)) order = 3
600            !To make the order at boundary cells for csi equal to 3.
601            !print*,'ob2l = ',ob2l,ob2r
602            ib = MAX(1 ,pc(1) - (ob2l - 1)) !non-periodic
603            ie = MIN(im,pc(1) + ob2r)
604               IF (ib.EQ.1 ) ie = ib + order - 1
605               IF (ie.EQ.im) ib = ie - order + 1
606     
607     
608            jb = MAX(1 ,pc(2) - (ob2l - 1)) !non-periodic
609            je = MIN(jm,pc(2) + ob2r)
610     
611               IF (jb.EQ.1 ) je = jb + order - 1
612               IF (je.EQ.jm) jb = je - order + 1
613     
614     
615            kb = MAX(1 ,pc(3) - (ob2l - 1)) !non-periodic
616            ke = MIN(km,pc(3) + ob2r)
617     
618               IF (kb.EQ.1 ) ke = kb + order - 1
619               IF (ke.EQ.km) kb = ke - order + 1
620     
621     
622            ob2r =  ob2rtmp
623            ob2l = ob2ltmp
624            ordernewtmp = order
625     
626            !Print*,'ib,ie .... in processing =', ib,ie,jb,je,kb,ke,&
627            !     & ordernewtmp
628            !PRINT*, 'pc = ',pc(1),pc(2),pc(3)!
629            order = ordertemp !reset the order
630     
631         CASE('lpi')
632            !print*, 'order in set stencil = ', order,pc(3)
633            !Print*,'ib, ie = ', pc(1), ib, ie
634     
635     
636            ib = MAX(1 ,pc(1) - (ob2l - 1)) !non-periodic
637            ie = MIN(im,pc(1) + ob2r)
638     
639               IF (ib.EQ.1 ) ie = ib + order - 1
640               IF (ie.EQ.im) ib = ie - order + 1
641     
642     
643            jb = MAX(1 ,pc(2) - (ob2l - 1)) !non-periodic
644            je = MIN(jm,pc(2) + ob2r)
645     
646               IF (jb.EQ.1 ) je = jb + order - 1
647               IF (je.EQ.jm) jb = je - order + 1
648     
649     
650            kb = MAX(1 ,pc(3) - (ob2l - 1)) !non-periodic
651            ke = MIN(km,pc(3) + ob2r)
652     
653               IF (kb.EQ.1 ) ke = kb + order - 1
654               IF (ke.EQ.km) kb = ke - order + 1
655     
656            ordernewtmp = order
657         END SELECT
658         IF(dimprob == 2) THEN
659            kb = pc(3)
660            ke = pc(3)
661            !Print*,'kb,ke  =', kb,ke
662         ENDIF
663         IF(PRESENT(ordernew)) ordernew = ordernewtmp
664       END SUBROUTINE set_interpolation_pstencil
665     
666     
667     
668     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
669       SUBROUTINE interp_oned_scalar(coor,scalar,ppos,interp_scl,order, &
670           isch, weight_pointer)
671     
672     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
673     
674     !-----------------------------------------------
675     ! Local variables
676     !-----------------------------------------------
677         REAL(prcn), DIMENSION(:), INTENT(in):: coor, scalar
678         REAL(prcn), INTENT(in):: ppos
679         REAL(prcn), INTENT(out):: interp_scl
680         INTEGER, INTENT(in):: order
681         CHARACTER(LEN=5), INTENT(in) :: isch
682         REAL(prcn), DIMENSION(:), POINTER, OPTIONAL:: weight_pointer
683     
684         REAL(prcn), DIMENSION(:), ALLOCATABLE:: zetacsi
685         INTEGER ::  iorig, i
686         REAL(prcn):: zeta
687     !-----------------------------------------------
688     
689         DO i = 1,order-1
690            dx(i) = coor(i+1) - coor(i)
691         ENDDO
692     
693     
694         SELECT CASE(isch)
695         CASE('lpi')
696     
697            !order = SIZE(coor)
698            iorig = order/2
699     
700            zeta = (ppos - coor(iorig))/dx(iorig)
701     
702            !-------
703            ! Get shape function values
704            !-------
705            SELECT CASE (order)
706            CASE (2)
707               DO i = 1,order
708                  xval(i) = shape2(zeta,i)
709               ENDDO
710            CASE (3)
711               DO i = 1,order
712                  xval(i) = shape3(zeta,i,dx)
713               ENDDO
714            CASE (4)
715               DO i = 1,order
716                  xval(i) = shape4(zeta,i,dx)
717               ENDDO
718            CASE (5)
719               DO i = 1,order
720                  xval(i) = shape5(zeta,i,dx)
721               ENDDO
722            CASE (6)
723               DO i = 1,order
724                  xval(i) = shape6(zeta,i,dx)
725               ENDDO
726            END SELECT
727         CASE('csi')
728            ! order = SIZE(coor,1)
729            iorig = (order+1)/2
730            !-------
731            ! Find out center cell widths
732            !-------
733            ALLOCATE(zetacsi(order))
734            !Zetacsi as defined in Yueng and Pope hence the name
735            !The defintions for zetacsi are true only for a structured grid
736            !if (order.eq.4) then
737     
738            !end if
739            !-------
740            ! Get shape function values
741            !-------
742            SELECT CASE (order)
743            CASE (4)
744               DO i = 1, order
745                  zetacsi(i) = (-ppos + coor(i))/dx(1)
746               END DO
747               DO i = 1,order
748                  xval(i) = shape4csi(zetacsi(i),i,dx,1)
749               ENDDO
750     
751            end SELECT
752     
753            DEALLOCATE(zetacsi)
754         end SELECT !Scheme
755     
756         !-------
757         ! Calculate interpolated value
758         !-------
759         interp_scl = 0.0
760         DO i = 1,order
761            interp_scl = interp_scl + scalar(i)*xval(i)
762         ENDDO
763     
764         !-------
765         ! Return the weights (optional)
766         !-------
767         IF (PRESENT(weight_pointer)) THEN
768            weight_pointer => xval
769         ENDIF
770     
771       END SUBROUTINE interp_oned_scalar
772     
773     
774     
775     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
776       SUBROUTINE interp_oned_vector(coor,vector,ppos,interp_vec,order,&
777           fun, weight_pointer)
778     
779     ! Interpolate an arbitrary sized array in one space dimension.
780     ! Uses the scalar interpolation to obtain the weights.
781     
782     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
783     
784     !-----------------------------------------------
785     ! Local variables
786     !-----------------------------------------------
787         REAL(prcn), DIMENSION(:), INTENT(in):: coor
788         REAL(prcn), DIMENSION(:,:), INTENT(in):: vector
789         REAL(prcn), INTENT(in):: ppos
790         REAL(prcn), DIMENSION(:), INTENT(out):: interp_vec
791         INTEGER, INTENT(in) :: order
792         CHARACTER(LEN=5) :: fun
793         REAL(prcn), DIMENSION(:), POINTER, OPTIONAL:: weight_pointer
794     
795         INTEGER:: vec_size, nv, i
796         REAL(prcn), DIMENSION(:), POINTER:: weights_scalar
797     !-----------------------------------------------
798     
799         !order    = SIZE(coor)
800         !print*,'In Interp_onedvector:order = ',order
801         vec_size = SIZE(vector,2)
802     
803         !-------
804         ! Interpolate first component and get weights
805         !-------
806         CALL interp_oned_scalar(coor,vector(:,1),ppos             &
807              ,interp_vec(1),order, fun, weights_scalar)
808     
809         !-------
810         ! Interpolate remaining components
811         !-------
812         DO nv = 2,vec_size
813            interp_vec(nv) = 0.0
814            DO i = 1,order
815               interp_vec(nv) =  interp_vec(nv)  &
816                    + vector(i,nv)*weights_scalar(i)
817            ENDDO
818         ENDDO
819     
820         !-------
821         ! Return the weights (optional)
822         !-------
823         IF (PRESENT(weight_pointer)) THEN
824            weight_pointer => weights_scalar
825         ENDIF
826     
827       END SUBROUTINE interp_oned_vector
828     
829     
830     
831     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
832       SUBROUTINE interp_twod_scalar(coor,scalar,ppos,interp_scl,order,&
833           isch,weight_pointer)
834     
835     ! Interpolate a scalar quantity in two dimensions.
836     
837     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
838     
839     !-----------------------------------------------
840     ! Local variables
841     !-----------------------------------------------
842         REAL(prcn), DIMENSION(:,:,:), INTENT(in):: coor
843         REAL(prcn), DIMENSION(:,:), INTENT(in):: scalar
844         REAL(prcn), DIMENSION(2), INTENT(in):: ppos
845         REAL(prcn), INTENT(out):: interp_scl
846         INTEGER, INTENT(in):: order
847         CHARACTER(LEN=5), INTENT(in) :: isch
848         REAL(prcn), DIMENSION(:,:,:), POINTER, OPTIONAL:: weight_pointer
849     
850         REAL(prcn), DIMENSION(:,:), ALLOCATABLE:: zetacsi
851         INTEGER:: i, j
852         INTEGER:: iorig
853         REAL(prcn), DIMENSION(2):: zeta
854     !-----------------------------------------------
855     
856         !-------
857         ! Get order of interpolation
858         !-------
859         !
860         weights  = zero
861         DO i = 1,order-1
862            dx(i) = coor(i+1,1,1)-coor(i,1,1)
863            dy(i) = coor(1,i+1,2)-coor(1,i,2)
864            !dz(i) = coor(1,1,i+1,3)-coor(1,1,i,3)
865         ENDDO
866     
867         SELECT CASE(isch)
868     
869         CASE('lpi')
870            !order = SIZE(coor,1)
871            iorig = order/2
872     
873            !-------
874            ! Find out center cell widths
875            !-------
876     
877            !-------
878            ! Zeta as defined in Bala/Maxey
879            !-------
880            zeta(1:2) = ppos(1:2) - coor(iorig,iorig,1:2)
881            zeta(1) = zeta(1)/dx(iorig)
882            zeta(2) = zeta(2)/dy(iorig)
883     
884            !-------
885            ! Get shape function values
886            !-------
887            SELECT CASE (order)
888            CASE (2)
889               DO i = 1,order
890                  xval(i) = shape2(zeta(1),i)
891                  yval(i) = shape2(zeta(2),i)
892                  !zval(i) = shape2(zeta(3),i)
893               ENDDO
894            CASE (3)
895               DO i = 1,order
896                  xval(i) = shape3(zeta(1),i,dx)
897                  yval(i) = shape3(zeta(2),i,dy)
898                  !zval(i) = shape3(zeta(3),i,dz)
899               ENDDO
900            CASE (4)
901               DO i = 1,order
902                  ! print*, 'in interp....zetayp(3,i) = ', zetayp(3,i),zval(i),i
903                  xval(i) = shape4(zeta(1),i,dx)
904                  yval(i) = shape4(zeta(2),i,dy)
905                  !zval(i) = shape4(zeta(3),i,dz)
906                  !Print*, ppos(1:3)
907                  !Print*,'int',i,xval(i),yval(i),zval(i)
908     !!$          xval(i) = shape4new(ppos(1),coor(1:order,1,1,1),i)
909     !!$          yval(i) = shape4new(ppos(2),coor(1,1:order,1,2),i)
910     !!$          zval(i) = shape4new(ppos(3),coor(1,1,1:order,3),i)
911               ENDDO
912            CASE (5)
913               DO i = 1,order
914                  xval(i) = shape5(zeta(1),i,dx)
915                  yval(i) = shape5(zeta(2),i,dy)
916                  !zval(i) = shape5(zeta(3),i,dz)
917               ENDDO
918            CASE (6)
919               DO i = 1,order
920                  xval(i) = shape6(zeta(1),i,dx)
921                  yval(i) = shape6(zeta(2),i,dy)
922                  !zval(i) = shape6(zeta(3),i,dz)
923               ENDDO
924            END SELECT
925     
926         CASE('csi')
927            ! order = SIZE(coor,1)
928            iorig = (order+1)/2
929     
930            !-------
931            ! Find out center cell widths
932            !-------
933            ALLOCATE(zetacsi(2,order))
934            !Zetacsi as defined in Yueng and Pope hence the name
935            !The defintions for zetacsi are true only for a structured grid
936     
937            !-------
938            ! Get shape function values
939            !-------
940            SELECT CASE (order)
941            CASE (4)
942               DO i = 1, order
943                  zetacsi(1,i) = (-ppos(1) + coor(i,1,1))/dx(1)
944                  zetacsi(2,i) = (-ppos(2) + coor(1,i,2))/dy(1)
945                  !zetacsi(3,i) = (-ppos(3) + coor(1,1,i,3))/dz(1)
946               END DO
947               DO i = 1,order
948                  xval(i) = shape4csi(zetacsi(1,i),i,dx,1)
949                  yval(i) = shape4csi(zetacsi(2,i),i,dy,2)
950                  !zval(i) = shape4csi(zetacsi(3,i),i,dz,3)
951     
952                  !Print*,'zetacsi  = ', zetacsi(3,i), coor(1,1,i,3), zval(i)
953               ENDDO
954            CASE(3)
955               DO i = 1, order
956     
957                  zetacsi(1,i) = ((-ppos(1) + coor(i,1,1))/dx(1))
958                  zetacsi(2,i) =((-ppos(2) + coor(1,i,2))/dy(1))
959                  !zetacsi(3,i) = ((-ppos(3) + coor(1,1,i,3))/dz(1))
960     !!$             zetacsi(1,i) = (ppos(1) - coor(1,1,1,1))/(coor(order,1,1&
961     !!$                  &,1)-coor(1,1,1,1))
962     !!$             zetacsi(2,i) = (ppos(2) - coor(1,1,1,2))/(coor(1,order,1&
963     !!$                  &,2)-coor(1,1,1,2))
964     !!$             zetacsi(3,i) = (ppos(3) - coor(1,1,1,3))/(coor(1,1,order&
965     !!$                  &,3)-coor(1,1,1,3))
966               END DO
967               DO i = 1,order
968                  if((xval(1)-coor(1,1,1)).lt.dx(1)) then
969                     xval(i) = shape3csileft(zetacsi(1,i),i,dx,1)
970                  else
971                     xval(i) = shape3csiright(zetacsi(1,i),i,dx,1)
972                  endif
973                  if((yval(1)-coor(1,1,2)).lt.dy(1)) then
974     
975                     yval(i) = shape3csileft(zetacsi(2,i),i,dy,2)
976                  else
977     
978                     yval(i) = shape3csiright(zetacsi(2,i),i,dy,2)
979                  end if
980     
981                  print*,'zeta = ',zetacsi(1,i), xval(i),i
982               ENDDO
983            END SELECT
984     
985            DEALLOCATE(zetacsi)
986         END SELECT !SCHEME
987         !-------
988         ! Calculate weights for the different nodes
989         !-------
990         ! DO 10 k = 1,order
991     !!$    DO 10 j = 1,order
992     !!$    DO 10 i = 1,order
993     !!$
994     !!$    10 CONTINUE
995         !If(order.eq.3) Print*,'in interpo...sum wt=,',sum(weights),order
996     
997         !-------
998         ! Calculate the interpolated value
999         !-------
1000         interp_scl = 0.0
1001         !DO 20 k = 1,order
1002         DO  j = 1,order
1003            DO  i = 1,order
1004               weights(i,j,1) = xval(i)*yval(j)
1005               interp_scl = interp_scl + scalar(i,j)*weights(i,j,1)
1006     
1007            end DO
1008         end DO
1009     
1010         !-------
1011         ! Return the weights for the force distribution (optional)
1012         !-------
1013         IF (PRESENT(weight_pointer)) THEN
1014            weight_pointer => weights
1015         ENDIF
1016     
1017       END SUBROUTINE interp_twod_scalar
1018     
1019     
1020     
1021     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1022       SUBROUTINE interp_twod_vector(coor,vector,ppos,interp_vec,order,&
1023           fun,weight_pointer)
1024     
1025     ! Interpolate an arbitrary sized array in two dimensions.
1026     ! Uses the scalar interpolation to obtain the weights.
1027     
1028     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1029     
1030     !-----------------------------------------------
1031     ! Local variables
1032     !-----------------------------------------------
1033         REAL(prcn), DIMENSION(:,:,:), INTENT(in):: coor, vector
1034         REAL(prcn), DIMENSION(2), INTENT(in):: ppos
1035         REAL(prcn), DIMENSION(:), INTENT(out):: interp_vec
1036         INTEGER, INTENT(in):: order
1037         CHARACTER(LEN=5) :: fun
1038         REAL(prcn), DIMENSION(:,:,:), POINTER, OPTIONAL:: weight_pointer
1039     
1040         INTEGER :: vec_size
1041         INTEGER:: i, j, nv
1042         REAL(prcn), DIMENSION(:,:,:), POINTER:: weights_scalar
1043     !-----------------------------------------------
1044     
1045         !-------
1046         ! Get order of interpolation
1047         !-------
1048         !order    = SIZE(coor,1)
1049         vec_size = SIZE(vector,3)
1050         !print*,'In Interp_threedvector:order = ',order,vec_size
1051     
1052         CALL interp_twod_scalar(coor,vector(:,:,1),ppos  &
1053              ,interp_vec(1),order,fun,weights_scalar)
1054     
1055         !-------
1056         ! Calculate the interpolated velocities
1057         !-------
1058         DO nv = 2,vec_size
1059            interp_vec(nv) = 0.0
1060            !DO  k = 1,order
1061            DO j = 1,order
1062               DO  i = 1,order
1063                  interp_vec(nv) =  interp_vec(nv)  &
1064                       + vector(i,j,nv)*weights_scalar(i,j,1)
1065               ENDDO
1066            ENDDO
1067         ENDDO
1068     
1069     
1070         !-------
1071         ! Return the weights for the force distribution (optional)
1072         !-------
1073         IF (PRESENT(weight_pointer)) THEN
1074            weight_pointer => weights_scalar
1075         ENDIF
1076     
1077       END SUBROUTINE interp_twod_vector
1078     
1079     
1080     
1081     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1082       SUBROUTINE interp_threed_scalar(coor,scalar,ppos,interp_scl,order,&
1083           isch,weight_pointer)
1084     
1085     ! Interpolate a scalar quantity in three dimensions.
1086     
1087     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1088     
1089     !-----------------------------------------------
1090     ! Local variables
1091     !-----------------------------------------------
1092         REAL(prcn), DIMENSION(:,:,:,:), INTENT(in):: coor
1093         REAL(prcn), DIMENSION(:,:,:), INTENT(in):: scalar
1094         REAL(prcn), DIMENSION(3), INTENT(in):: ppos
1095         REAL(prcn), INTENT(out):: interp_scl
1096         INTEGER, INTENT(in):: order
1097         CHARACTER(LEN=5), INTENT(in) :: isch
1098         REAL(prcn), DIMENSION(:,:,:), POINTER, OPTIONAL:: weight_pointer
1099     
1100         REAL(prcn), DIMENSION(:,:), ALLOCATABLE:: zetacsi
1101         INTEGER:: i, j, k
1102         INTEGER:: iorig
1103         REAL(prcn) :: zeta(3), zetasph, sigma, bandwidth
1104         LOGICAL :: calcwts
1105     !-----------------------------------------------
1106     
1107         !-------
1108         ! Get order of interpolation
1109         !-------
1110         !
1111         weights  = zero
1112         DO i = 1,order-1
1113            dx(i) = coor(i+1,1,1,1)-coor(i,1,1,1)
1114            dy(i) = coor(1,i+1,1,2)-coor(1,i,1,2)
1115            dz(i) = coor(1,1,i+1,3)-coor(1,1,i,3)
1116         ENDDO
1117     
1118         !Print*,'In interpolator', isch
1119         SELECT CASE(isch)
1120     
1121         CASE('lpi')
1122            calcwts = .true.
1123            !order = SIZE(coor,1)
1124            iorig = order/2
1125     
1126            !-------
1127            ! Find out center cell widths
1128            !-------
1129     
1130            !-------
1131            ! Zeta as defined in Bala/Maxey
1132            !-------
1133            zeta(1:3) = ppos(1:3) - coor(iorig,iorig,iorig,1:3)
1134            zeta(1) = zeta(1)/dx(iorig)
1135            zeta(2) = zeta(2)/dy(iorig)
1136            zeta(3) = zeta(3)/dz(iorig)
1137     
1138            !-------
1139            ! Get shape function values
1140            !-------
1141            SELECT CASE (order)
1142            CASE (2)
1143               DO i = 1,order
1144                  xval(i) = shape2(zeta(1),i)
1145                  yval(i) = shape2(zeta(2),i)
1146                  zval(i) = shape2(zeta(3),i)
1147               ENDDO
1148            CASE (3)
1149               DO i = 1,order
1150                  xval(i) = shape3(zeta(1),i,dx)
1151                  yval(i) = shape3(zeta(2),i,dy)
1152                  zval(i) = shape3(zeta(3),i,dz)
1153               ENDDO
1154            CASE (4)
1155               DO i = 1,order
1156                  ! print*, 'in interp....zetayp(3,i) = ', zetayp(3,i),zval(i),i
1157                  xval(i) = shape4(zeta(1),i,dx)
1158                  yval(i) = shape4(zeta(2),i,dy)
1159                  zval(i) = shape4(zeta(3),i,dz)
1160                  !Print*, ppos(1:3)
1161                  !Print*,'int',i,xval(i),yval(i),zval(i)
1162     !!$          xval(i) = shape4new(ppos(1),coor(1:order,1,1,1),i)
1163     !!$          yval(i) = shape4new(ppos(2),coor(1,1:order,1,2),i)
1164     !!$          zval(i) = shape4new(ppos(3),coor(1,1,1:order,3),i)
1165               ENDDO
1166            CASE (5)
1167               DO i = 1,order
1168                  xval(i) = shape5(zeta(1),i,dx)
1169                  yval(i) = shape5(zeta(2),i,dy)
1170                  zval(i) = shape5(zeta(3),i,dz)
1171               ENDDO
1172            CASE (6)
1173               DO i = 1,order
1174                  xval(i) = shape6(zeta(1),i,dx)
1175                  yval(i) = shape6(zeta(2),i,dy)
1176                  zval(i) = shape6(zeta(3),i,dz)
1177               ENDDO
1178            END SELECT
1179     
1180         CASE('csi')
1181     
1182            calcwts = .true.
1183            ! order = SIZE(coor,1)
1184            iorig = (order+1)/2
1185     
1186            !-------
1187            ! Find out center cell widths
1188            !-------
1189            ALLOCATE(zetacsi(3,order))
1190            !Zetacsi as defined in Yueng and Pope hence the name
1191            !The defintions for zetacsi are true only for a structured grid
1192     
1193            !-------
1194            ! Get shape function values
1195            !-------
1196            SELECT CASE (order)
1197            CASE (4)
1198               DO i = 1, order
1199                  zetacsi(1,i) = (-ppos(1) + coor(i,1,1,1))/dx(1)
1200                  zetacsi(2,i) = (-ppos(2) + coor(1,i,1,2))/dy(1)
1201                  zetacsi(3,i) = (-ppos(3) + coor(1,1,i,3))/dz(1)
1202               END DO
1203               DO i = 1,order
1204                  xval(i) = shape4csi(zetacsi(1,i),i,dx,1)
1205                  yval(i) = shape4csi(zetacsi(2,i),i,dy,2)
1206                  zval(i) = shape4csi(zetacsi(3,i),i,dz,3)
1207     
1208                  !Print*,'zetacsi  = ', zetacsi(3,i), coor(1,1,i,3), zval(i)
1209               ENDDO
1210            CASE(3)
1211               DO i = 1, order
1212     
1213                  zetacsi(1,i) = ((-ppos(1) + coor(i,1,1,1))/dx(1))
1214                  zetacsi(2,i) =((-ppos(2) + coor(1,i,1,2))/dy(1))
1215                  zetacsi(3,i) = ((-ppos(3) + coor(1,1,i,3))/dz(1))
1216     !!$             zetacsi(1,i) = (ppos(1) - coor(1,1,1,1))/(coor(order,1,1&
1217     !!$                  &,1)-coor(1,1,1,1))
1218     !!$             zetacsi(2,i) = (ppos(2) - coor(1,1,1,2))/(coor(1,order,1&
1219     !!$                  &,2)-coor(1,1,1,2))
1220     !!$             zetacsi(3,i) = (ppos(3) - coor(1,1,1,3))/(coor(1,1,order&
1221     !!$                  &,3)-coor(1,1,1,3))
1222               END DO
1223               DO i = 1,order
1224                  if((xval(1)-coor(1,1,1,1)).lt.dx(1)) then
1225                     xval(i) = shape3csileft(zetacsi(1,i),i,dx,1)
1226                  else
1227                     xval(i) = shape3csiright(zetacsi(1,i),i,dx,1)
1228                  endif
1229                  if((yval(1)-coor(1,1,1,2)).lt.dy(1)) then
1230     
1231                     yval(i) = shape3csileft(zetacsi(2,i),i,dy,2)
1232                  else
1233     
1234                     yval(i) = shape3csiright(zetacsi(2,i),i,dy,2)
1235                  end if
1236                  if((zval(1)-coor(1,1,1,3)).lt.dz(1)) then
1237                     zval(i) = shape3csileft(zetacsi(3,i),i,dz,3)
1238                  else
1239     
1240                     zval(i) = shape3csiright(zetacsi(3,i),i,dz,3)
1241                  endif
1242     
1243                  print*,'zeta = ',zetacsi(1,i), xval(i),i
1244               ENDDO
1245            END SELECT
1246     
1247            DEALLOCATE(zetacsi)
1248     
1249         CASE('sph')
1250     
1251            iorig = (order+1)/2
1252            calcwts = .false.
1253            !-------
1254            SELECT CASE (order)
1255            CASE (4)
1256               bandwidth  = one*(dx(1)*dy(1))**(half)
1257               sigma = one/pi
1258               do k = 1, order
1259                  do j = 1, order
1260                     DO i = 1, order
1261                        zetasph = (-ppos(1) + coor(i,j,k,1))**2.0
1262                        zetasph = zetasph + (-ppos(2) + coor(i,j,k,2))**2.0
1263                        !zetasph = zetasph + (-ppos(3) + coor(i,j,k,3))**2.0
1264                        zetasph = sqrt(zetasph)
1265     
1266                        zetasph = zetasph/bandwidth
1267     
1268                        if(zetasph.ge.zero.and.zetasph.lt.one) then
1269                           weights(i,j,k) = one - (three/two)*zetasph**two&
1270                                & + (three/four)*zetasph**three
1271                        elseif(zetasph.ge.one.and.zetasph.lt.two) then
1272                           weights(i,j,k) = fourth*(two-zetasph)**three
1273                        else
1274                           weights(i,j,k) = zero
1275                        end if
1276                        weights(i,j,k) = (sigma*weights(i,j,k))!/bandwidth&
1277                             !&**three
1278     
1279                        !Print*,weights(i,j,k), i,j,k
1280                     END DO
1281                  end do
1282               end DO
1283            END SELECT
1284     
1285         END SELECT !SCHEME
1286         !-------
1287         ! Calculate weights for the different nodes
1288         !-------
1289         if(calcwts) then
1290            DO  k = 1,order
1291               DO  j = 1,order
1292                  DO  i = 1,order
1293                     weights(i,j,k) = xval(i)*yval(j)*zval(k)
1294                  end DO
1295               end DO
1296            end DO
1297         end if
1298     
1299         !If(order.eq.3) Print*,'in interpo...sum wt=,',sum(weights),order
1300     
1301         !-------
1302         ! Calculate the interpolated value
1303         !-------
1304         interp_scl = 0.0
1305         DO  k = 1,order
1306            DO  j = 1,order
1307               DO  i = 1,order
1308                  interp_scl = interp_scl + scalar(i,j,k)*weights(i,j,k)
1309               end DO
1310            end DO
1311         end DO
1312     
1313     
1314         !-------
1315         ! Return the weights for the force distribution (optional)
1316         !-------
1317         IF (PRESENT(weight_pointer)) THEN
1318            weight_pointer => weights
1319         ENDIF
1320       END SUBROUTINE interp_threed_scalar
1321     
1322     
1323     
1324     
1325     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1326       SUBROUTINE interp_threed_vector(coor,vector,ppos,interp_vec,order,&
1327           fun,weight_pointer)
1328     
1329     ! Interpolate an arbitrary sized array in three dimensions.
1330     ! Uses the scalar interpolation to obtain the weights.
1331     
1332     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1333     
1334     !-----------------------------------------------
1335     ! Local variables
1336     !-----------------------------------------------
1337         REAL(prcn), DIMENSION(:,:,:,:), INTENT(in):: coor, vector
1338         REAL(prcn), DIMENSION(3), INTENT(in):: ppos
1339         REAL(prcn), DIMENSION(:), INTENT(out):: interp_vec
1340         INTEGER, INTENT(in):: order
1341         CHARACTER(LEN=5) :: fun
1342         REAL(prcn), DIMENSION(:,:,:), POINTER, OPTIONAL:: weight_pointer
1343     
1344         INTEGER :: vec_size
1345         INTEGER:: i, j, k, nv
1346         REAL(prcn), DIMENSION(:,:,:), POINTER:: weights_scalar
1347     !-----------------------------------------------
1348     
1349         !-------
1350         ! Get order of interpolation
1351         !-------
1352         !order    = SIZE(coor,1)
1353         vec_size = SIZE(vector,4)
1354         !print*,'In Interp_threedvector:order = ',order,vec_size
1355     
1356         CALL interp_threed_scalar(coor,vector(:,:,:,1),ppos  &
1357              ,interp_vec(1),order,fun,weights_scalar)
1358     
1359         !-------
1360         ! Calculate the interpolated velocities
1361         !-------
1362         DO nv = 2,vec_size
1363            interp_vec(nv) = 0.0
1364            DO  k = 1,order
1365               DO j = 1,order
1366                  DO  i = 1,order
1367                     interp_vec(nv) =  interp_vec(nv)  &
1368                          + vector(i,j,k,nv)*weights_scalar(i,j,k)
1369                  ENDDO
1370               ENDDO
1371            ENDDO
1372         ENDDO
1373     
1374         !-------
1375         ! Return the weights for the force distribution (optional)
1376         !-------
1377         IF (PRESENT(weight_pointer)) THEN
1378            weight_pointer => weights_scalar
1379         ENDIF
1380     
1381       END SUBROUTINE interp_threed_vector
1382     
1383     
1384     
1385     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1386       SUBROUTINE calc_weightderiv_threed(coor,ppos,weight_pointer,order, isch)
1387     
1388     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1389     
1390     !-----------------------------------------------
1391     ! Local variables
1392     !-----------------------------------------------
1393         REAL(prcn), DIMENSION(:,:,:,:), INTENT(in):: coor
1394         REAL(prcn), DIMENSION(3), INTENT(in):: ppos
1395         REAL(prcn), DIMENSION(:,:,:,:) :: weight_pointer
1396         INTEGER, INTENT(in) :: order
1397         CHARACTER(len=5), INTENT(in) :: isch
1398     
1399         INTEGER:: i, j, k, kk
1400         Real(prcn) :: dx1,dy1,dz1
1401         INTEGER :: iorig
1402         REAL(prcn):: zeta(3), zetasph, bandwidth, sigma, tmp
1403         REAL(prcn), DIMENSION(3,order) :: zetacsi !right now only
1404         ! for cubic spline interp rg 03/14/05
1405     !-----------------------------------------------
1406     
1407         weight_pointer = zero
1408         !-------
1409         ! Get order of interpolation
1410         !-------
1411         !order = SIZE(coor,1)
1412         iorig = order/2
1413         ! print*,'in interp...Iorig = ',iorig !Debug
1414     
1415         !-------
1416         ! Find out center cell widths
1417         !-------
1418         DO i = 1,order-1
1419            dx(i) = coor(i+1,1,1,1)-coor(i,1,1,1)
1420            dy(i) = coor(1,i+1,1,2)-coor(1,i,1,2)
1421            dz(i) = coor(1,1,i+1,3)-coor(1,1,i,3)
1422         ENDDO
1423     
1424         !Zetacsi as defined in Yueng and Pope hence the name
1425         !The defintions for zetacsi are true only for a structured grid
1426         !if (order.eq.4) then
1427     
1428         !end if
1429         !-------
1430         ! Zeta as defined in Bala/Maxey
1431         !-------
1432     
1433     
1434         !-------
1435         ! Get shape function values
1436         !-------
1437         SELECT CASE (isch)
1438         CASE ('csi')
1439            !Allocate(zetacsi(3,order))
1440     
1441            DO k = 1,3
1442               SELECT CASE(order)
1443               CASE(4)
1444                  DO i = 1, order
1445                     zetacsi(1,i) = (-ppos(1) + coor(i,1,1,1))/dx(1)
1446                     zetacsi(2,i) = (-ppos(2) + coor(1,i,1,2))/dy(1)
1447                     zetacsi(3,i) = (-ppos(3) + coor(1,1,i,3))/dz(1)
1448                  END DO
1449                  DO i = 1,order
1450                     IF(k.EQ.1) THEN
1451                        xval(i) = (shape4deriv_csi(zetacsi(1,i),i,dx))/dx(1)
1452                        yval(i) = shape4csi(zetacsi(2,i),i,dy,2)
1453                        zval(i) = shape4csi(zetacsi(3,i),i,dz,3)
1454                     ELSEIF(k.EQ.2) THEN
1455                        xval(i) = shape4csi(zetacsi(1,i),i,dx,1)
1456                        yval(i) = (shape4deriv_csi(zetacsi(2,i),i,dy))/(dy(1))
1457                        zval(i) = shape4csi(zetacsi(3,i),i,dz,3)
1458                     ELSEIF(k.EQ.3) THEN
1459                        xval(i) = shape4csi(zetacsi(1,i),i,dx,1)
1460                        yval(i) = shape4csi(zetacsi(2,i),i,dy,2)
1461                        zval(i) = (shape4deriv_csi(zetacsi(3,i),i,dz))/(dz(1))
1462                     ENDIF
1463                  ENDDO
1464               CASE(3)
1465                  dx1 = coor(order,1,1,1)-coor(1,1,1,1)
1466                  dy1 = coor(1,order,1,2)-coor(1,1,1,2)
1467                  dz1 = coor(1,1,order,3)-coor(1,1,1,3)
1468                  zetacsi(1,i) = (ppos(1) - coor(1,1,1,1))/dx1
1469                  zetacsi(2,i) = (ppos(2) - coor(1,1,1,2))/dy1
1470                  zetacsi(3,i) = (ppos(3) - coor(1,1,1,3))/dz1
1471                  DO i = 1,order
1472                     IF(k.EQ.1) THEN
1473                        xval(i) = (shape3deriv_csi(zetacsi(1,i),i,dx))/dx1
1474                        yval(i) = shape3csileft(zetacsi(2,i),i,dy,2)
1475                        zval(i) = shape3csileft(zetacsi(3,i),i,dz,3)
1476                     ELSEIF(k.EQ.2) THEN
1477                        xval(i) = shape3csileft(zetacsi(1,i),i,dx,1)
1478                        yval(i) = (shape3deriv_csi(zetacsi(2,i),i,dy))/(dy1)
1479                        zval(i) = shape3csileft(zetacsi(3,i),i,dz,3)
1480                     ELSEIF(k.EQ.3) THEN
1481                        xval(i) = shape3csileft(zetacsi(1,i),i,dx,1)
1482                        yval(i) = shape3csileft(zetacsi(2,i),i,dy,2)
1483                        zval(i) = (shape3deriv_csi(zetacsi(3,i),i,dz))/(dz1)
1484                     ENDIF
1485                  ENDDO
1486               END SELECT!order
1487               DO kk = 1,order
1488                  DO j = 1,order
1489                     DO i = 1,order
1490                        weight_pointer(i,j,kk,k) = xval(i)*yval(j)*zval(kk)
1491                     ENDDO
1492                  ENDDO
1493               ENDDO
1494            ENDDO!end loop over the coordinate directions
1495            !deallocate(zetacsi)
1496         CASE('lpi')
1497            zeta(1:3) = ppos(1:3) - coor(iorig,iorig,iorig,1:3)
1498            zeta(1) = zeta(1)/dx(iorig)
1499            zeta(2) = zeta(2)/dy(iorig)
1500            zeta(3) = zeta(3)/dz(iorig)
1501            DO k = 1,3
1502               SELECT CASE(order)
1503               CASE(4)
1504                  DO i = 1,order
1505                     IF(k.EQ.1) THEN
1506                        xval(i) = (shape4deriv(zeta(1),i,dx))/dx(iorig)
1507                        yval(i) = shape4(zeta(2),i,dy)
1508                        zval(i) = shape4(zeta(3),i,dz)
1509                     ELSEIF(k.EQ.2) THEN
1510                        xval(i) = shape4(zeta(1),i,dx)
1511                        yval(i) = (shape4deriv(zeta(2),i,dy))/(dy(iorig))
1512                        zval(i) = shape4(zeta(3),i,dz)
1513                     ELSEIF(k.EQ.3) THEN
1514                        xval(i) = shape4(zeta(1),i,dx)
1515                        yval(i) = shape4(zeta(2),i,dy)
1516                        zval(i) = (shape4deriv(zeta(3),i,dz))/(dz(iorig))
1517                     ENDIF
1518                  ENDDO
1519               CASE(2)
1520                  DO i = 1,order
1521                     IF(k.EQ.1) THEN
1522                        xval(i) = (shape2deriv(zeta(1),i))/dx(iorig)
1523                        yval(i) = shape2(zeta(2),i)
1524                        zval(i) = shape2(zeta(3),i)
1525                     ELSEIF(k.EQ.2) THEN
1526                        xval(i) = shape2(zeta(1),i)
1527                        yval(i) = (shape2deriv(zeta(2),i))/(dy(iorig))
1528                        zval(i) = shape2(zeta(3),i)
1529                     ELSEIF(k.EQ.3) THEN
1530                        xval(i) = shape2(zeta(1),i)
1531                        yval(i) = shape2(zeta(2),i)
1532                        zval(i) = (shape2deriv(zeta(3),i))/(dz(iorig))
1533                     ENDIF
1534                  ENDDO
1535               end SELECT
1536     
1537               DO kk = 1,order
1538                  DO j = 1,order
1539                     DO i = 1,order
1540                        weight_pointer(i,j,kk,k) = -xval(i)*yval(j)*zval(kk)
1541                     ENDDO
1542                  ENDDO
1543               ENDDO
1544            ENDDO!end loop over the coordinate directions
1545     
1546         CASE('sph')
1547            SELECT CASE (order)
1548     
1549            CASE (4)
1550                         bandwidth  = one*(dx(1)*dy(1))**(half)
1551               sigma = one/pi
1552     
1553               do k = 1, order
1554                  do j = 1, order
1555                     DO i = 1, order
1556                        zetasph = (-ppos(1) + coor(i,j,k,1))**2.0
1557                        zetasph = zetasph + (-ppos(2) + coor(i,j,k,2))**2.0
1558                        !zetasph = zetasph + (-ppos(3) + coor(i,j,k,3))**2.0
1559                        zetasph = sqrt(zetasph)
1560     
1561                        zetasph = zetasph/bandwidth
1562     
1563                        if(zetasph.ge.zero.and.zetasph.lt.one) then
1564                           tmp = -two*(three/two)*zetasph &
1565                                & + three*(three/four)*zetasph**two
1566                        elseif(zetasph.ge.one.and.zetasph.lt.two) then
1567                           tmp = -three*fourth*(two-zetasph)**two
1568                        else
1569                           tmp = zero
1570                        end if
1571     
1572                        weight_pointer(i,j,k,1) = (tmp/zetasph)*(-ppos(1) &
1573                             &+ coor(i,j,k,1))
1574                        weight_pointer(i,j,k,2) = (tmp/zetasph)*(-ppos(2) &
1575                             &+ coor(i,j,k,2))
1576                        weight_pointer(i,j,k,3) = (tmp/zetasph)*(-ppos(3) &
1577                             &+ coor(i,j,k,3))
1578                        weight_pointer(i,j,k,:) = (sigma*weight_pointer(i,j&
1579                             &,k,:))/bandwidth&
1580                             &**two
1581     
1582                        !Print*,weights(i,j,k), i,j,k
1583                     END DO
1584                  end do
1585               end DO
1586            END SELECT
1587     
1588         END SELECT
1589       END SUBROUTINE calc_weightderiv_threed
1590     
1591     
1592     
1593     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1594       SUBROUTINE calc_weightderiv_twod(coor,ppos,weight_pointer,order, isch)
1595     
1596     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1597     
1598     !-----------------------------------------------
1599     ! Local variables
1600     !-----------------------------------------------
1601         REAL(prcn), DIMENSION(:,:,:), INTENT(in):: coor
1602         REAL(prcn), DIMENSION(2), INTENT(in):: ppos
1603         REAL(prcn), DIMENSION(:,:,:,:) :: weight_pointer
1604         INTEGER, INTENT(in) :: order
1605         CHARACTER(len=5), INTENT(in) :: isch
1606     
1607         INTEGER:: i, j, k
1608         REAL(prcn) :: dx1,dy1
1609         INTEGER :: iorig
1610         REAL(prcn):: zeta(3), zetasph, bandwidth, sigma, tmp
1611         REAL(prcn), DIMENSION(2,order) :: zetacsi !right now only
1612         ! for cubic spline interp rg 03/14/05
1613     !-----------------------------------------------
1614     
1615         weight_pointer = zero
1616         !-------
1617         ! Get order of interpolation
1618         !-------
1619         !order = SIZE(coor,1)
1620         iorig = order/2
1621         ! print*,'in interp...Iorig = ',iorig !Debug
1622     
1623         !-------
1624         ! Find out center cell widths
1625         !-------
1626         DO i = 1,order-1
1627            dx(i) = coor(i+1,1,1)-coor(i,1,1)
1628            dy(i) = coor(1,i+1,2)-coor(1,i,2)
1629         ENDDO
1630     
1631         !Zetacsi as defined in Yueng and Pope hence the name
1632         !The defintions for zetacsi are true only for a structured grid
1633         !if (order.eq.4) then
1634     
1635         !end if
1636         !-------
1637         ! Zeta as defined in Bala/Maxey
1638         !-------
1639     
1640     
1641         !-------
1642         ! Get shape function values
1643         !-------
1644         SELECT CASE (isch)
1645         CASE ('csi')
1646            !Allocate(zetacsi(3,order))
1647     
1648            DO k = 1,2
1649               SELECT CASE(order)
1650               CASE(4)
1651                  DO i = 1, order
1652                     zetacsi(1,i) = (-ppos(1) + coor(i,1,1))/dx(1)
1653                     zetacsi(2,i) = (-ppos(2) + coor(1,i,2))/dy(1)
1654     
1655                  END DO
1656                  DO i = 1,order
1657                     IF(k.EQ.1) THEN
1658                        xval(i) = (shape4deriv_csi(zetacsi(1,i),i,dx))/dx(1)
1659                        yval(i) = shape4csi(zetacsi(2,i),i,dy,2)
1660                     ELSEIF(k.EQ.2) THEN
1661                        xval(i) = shape4csi(zetacsi(1,i),i,dx,1)
1662                        yval(i) = (shape4deriv_csi(zetacsi(2,i),i,dy))/(dy(1))
1663                     ENDIF
1664                  ENDDO
1665               CASE(3)
1666                  dx1 = coor(order,1,1)-coor(1,1,1)
1667                  dy1 = coor(1,order,2)-coor(1,1,2)
1668                  zetacsi(1,i) = (ppos(1) - coor(1,1,1))/dx1
1669                  zetacsi(2,i) = (ppos(2) - coor(1,1,2))/dy1
1670                  DO i = 1,order
1671                     IF(k.EQ.1) THEN
1672                        xval(i) = (shape3deriv_csi(zetacsi(1,i),i,dx))/dx1
1673                        yval(i) = shape3csileft(zetacsi(2,i),i,dy,2)
1674                     ELSEIF(k.EQ.2) THEN
1675                        xval(i) = shape3csileft(zetacsi(1,i),i,dx,1)
1676                        yval(i) = (shape3deriv_csi(zetacsi(2,i),i,dy))/(dy1)
1677                     ENDIF
1678                  ENDDO
1679               END SELECT!order
1680                  DO j = 1,order
1681                     DO i = 1,order
1682                        weight_pointer(i,j,1,k) = xval(i)*yval(j)
1683                     ENDDO
1684                  ENDDO
1685               ENDDO!end loop over the coordinate directions
1686            !deallocate(zetacsi)
1687         CASE('lpi')
1688            zeta(1:2) = ppos(1:2) - coor(iorig,iorig,1:2)
1689            zeta(1) = zeta(1)/dx(iorig)
1690            zeta(2) = zeta(2)/dy(iorig)
1691            DO k = 1,2
1692               SELECT CASE(order)
1693               CASE(4)
1694                  DO i = 1,order
1695                     IF(k.EQ.1) THEN
1696                        xval(i) = (shape4deriv(zeta(1),i,dx))/dx(iorig)
1697                        yval(i) = shape4(zeta(2),i,dy)
1698                     ELSEIF(k.EQ.2) THEN
1699                        xval(i) = shape4(zeta(1),i,dx)
1700                        yval(i) = (shape4deriv(zeta(2),i,dy))/(dy(iorig))
1701                     ENDIF
1702                  ENDDO
1703               CASE(2)
1704                  DO i = 1,order
1705                     IF(k.EQ.1) THEN
1706                        xval(i) = (shape2deriv(zeta(1),i))/dx(iorig)
1707                        yval(i) = shape2(zeta(2),i)
1708                     ELSEIF(k.EQ.2) THEN
1709                        xval(i) = shape2(zeta(1),i)
1710                        yval(i) = (shape2deriv(zeta(2),i))/(dy(iorig))
1711                     ENDIF
1712                  ENDDO
1713               end SELECT
1714     
1715               DO j = 1,order
1716                     DO i = 1,order
1717                        weight_pointer(i,j,1,k) = -xval(i)*yval(j)
1718                     ENDDO
1719                  ENDDO
1720               ENDDO!end loop over the coordinate directions
1721     
1722         CASE('sph')
1723            SELECT CASE (order)
1724     
1725            CASE (4)
1726               bandwidth  = one*(dx(1)*dy(1))**(half)
1727               sigma = one/pi
1728     
1729               do j = 1, order
1730                  DO i = 1, order
1731                     zetasph = (-ppos(1) + coor(i,j,1))**2.0
1732                     zetasph = zetasph + (-ppos(2) + coor(i,j,2))**2.0
1733                     !zetasph = zetasph + (-ppos(3) + coor(i,j,k,3))**2.0
1734                     zetasph = sqrt(zetasph)
1735     
1736                     zetasph = zetasph/bandwidth
1737     
1738                     if(zetasph.ge.zero.and.zetasph.lt.one) then
1739                        tmp = -two*(three/two)*zetasph &
1740                             & + three*(three/four)*zetasph**two
1741                     elseif(zetasph.ge.one.and.zetasph.lt.two) then
1742                        tmp = -three*fourth*(two-zetasph)**two
1743                     else
1744                        tmp = zero
1745                     end if
1746     
1747                     weight_pointer(i,j,1,1) = (tmp/zetasph)*(-ppos(1) &
1748                             &+ coor(i,j,1))
1749                        weight_pointer(i,j,1,2) = (tmp/zetasph)*(-ppos(2) &
1750                             &+ coor(i,j,2))
1751                        !weight_pointer(i,j,1,3) = (tmp/zetasph)*(-ppos(3) &
1752                        !     &+ coor(i,j,3))
1753                        weight_pointer(i,j,1,:) = (sigma*weight_pointer(i,j&
1754                             &,1,:))/bandwidth&
1755                             &**two
1756     
1757                        !Print*,weights(i,j,k), i,j,k
1758                     END DO
1759                  end do
1760               END SELECT
1761     
1762            END SELECT
1763     
1764          END SUBROUTINE calc_weightderiv_twod
1765     
1766     
1767     
1768     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1769       FUNCTION justweights(coor,ppos)
1770     ! To calculate just the weights using trilinear interpolation
1771     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1772     
1773     !-----------------------------------------------
1774     ! Local variables
1775     !-----------------------------------------------
1776         REAL(prcn), DIMENSION(:,:,:), POINTER:: justweights
1777         REAL(prcn), DIMENSION(:,:,:,:), INTENT(IN):: coor
1778         REAL(prcn), DIMENSION(3), INTENT(IN):: ppos
1779     
1780         INTEGER, PARAMETER:: order=2
1781         INTEGER:: i, j, k, iorig
1782         REAL(prcn):: dxl, dyl, dzl
1783         REAL(prcn), DIMENSION(3):: zeta
1784     !-----------------------------------------------
1785     
1786         iorig = order/2
1787     
1788         dxl = coor(order,1,1,1)-coor(order-1,1,1,1)
1789         dyl = coor(1,order,1,2)-coor(1,order-1,1,2)
1790         dzl = coor(1,1,order,3)-coor(1,1,order-1,3)
1791     
1792         zeta(1:3) = ppos(1:3) - coor(iorig,iorig,iorig,1:3)
1793         zeta(1) = zeta(1)/dxl
1794         zeta(2) = zeta(2)/dyl
1795         zeta(3) = zeta(3)/dzl
1796     
1797         DO i = 1,order
1798            xval(i) = shape2(zeta(1),i)
1799            yval(i) = shape2(zeta(2),i)
1800            zval(i) = shape2(zeta(3),i)
1801         ENDDO
1802     
1803         !-------
1804         ! Calculate weights for the different nodes
1805         !-------
1806         DO  k = 1,order
1807            DO  j = 1,order
1808               DO  i = 1,order
1809                  weights(i,j,k) = xval(i)*yval(j)*zval(k)
1810               end DO
1811            end DO
1812         end DO
1813     
1814     
1815         justweights => weights
1816       END FUNCTION justweights
1817     
1818     
1819     
1820     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1821       FUNCTION shape2(zeta,i)
1822     ! Second-order (linear) shape functions
1823     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1824     
1825     !-----------------------------------------------
1826     ! Local variables
1827     !-----------------------------------------------
1828         REAL(prcn):: shape2
1829         REAL(prcn):: zeta
1830         INTEGER:: i
1831     !-----------------------------------------------
1832         SELECT CASE (i)
1833         CASE (1)
1834            shape2 = 1 - zeta
1835         CASE (2)
1836            shape2 = zeta
1837         END SELECT
1838       END FUNCTION shape2
1839     
1840     
1841     
1842     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1843       FUNCTION shape3(zeta,i,dx)
1844     ! Third-order (quadratic) shape functions
1845     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1846     
1847     !-----------------------------------------------
1848     ! Local variables
1849     !-----------------------------------------------
1850         REAL(prcn):: shape3
1851         REAL(prcn):: zeta
1852         REAL(prcn), DIMENSION(:):: dx
1853         INTEGER:: i
1854         REAL(prcn):: zh, num, denom
1855     !-----------------------------------------------
1856     
1857         SELECT CASE (i)
1858         CASE (1)
1859            zh    = dx(1)*zeta
1860            denom = dx(1)*(dx(1)+dx(2))
1861            num   = (zh - dx(1))*(zh - dx(1) - dx(2))
1862            shape3 = num/denom
1863         CASE (2)
1864            zh    = dx(1)*zeta
1865            denom = -dx(1)*dx(2)
1866            num   = zh*(zh - dx(1) - dx(2))
1867            shape3 = num/denom
1868         CASE (3)
1869            zh    = dx(1)*zeta
1870            denom = dx(2)*(dx(1)+dx(2))
1871            num   = zh*(zh - dx(1))
1872            shape3 = num/denom
1873         END SELECT
1874       END FUNCTION shape3
1875     
1876     
1877     
1878     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1879       FUNCTION shape4(zeta,i,dx)
1880     ! Fourth-order (cubic) shape functions
1881     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1882     
1883     !-----------------------------------------------
1884     ! Local variables
1885     !-----------------------------------------------
1886         REAL(prcn):: shape4
1887         REAL(prcn):: zeta
1888         REAL(prcn), DIMENSION(:):: dx
1889         INTEGER:: i
1890         REAL(prcn):: r1, r2, c1, c2, c3
1891     !-----------------------------------------------
1892     
1893         r1 = dx(2)/dx(1)
1894         r2 = dx(3)/dx(1)
1895     
1896         SELECT CASE (i)
1897         CASE (1)
1898            c1 = 1.0/(1.0 + r1)
1899            c3 = 1.0/(1.0 + r1 + r1*r2)
1900            shape4 = -c1*c3*(r1**3.0)*(zeta)*(zeta-1.0)*(zeta-(1.0+r2))
1901            !shape4 = -(one/6.)*(zeta)*(zeta-one)*(zeta-two)
1902         CASE (2)
1903            c2 = 1.0/(1.0 + r2)
1904            shape4 = c2*(zeta-1.0)*(zeta*r1+1.0)*(zeta-(1.0+r2))
1905            ! shape4 = half*(zeta-one)*(zeta+one)*(zeta-two)
1906         CASE (3)
1907            c1 = 1.0/(1.0 + r1)
1908            shape4 = -(c1/r2)*(zeta)*(zeta*r1+1.0)*(zeta-(1.0+r2))
1909            ! shape4 = -half*(zeta)*(zeta+1)*(zeta-two)
1910         CASE (4)
1911            c2 = 1.0/(1.0 + r2)
1912            c3 = 1.0/(1.0 + r1 + r1*r2)
1913            shape4 = (c3*c2/r2)*(zeta)*(zeta*r1+1.0)*(zeta-1.0)
1914            ! shape4 = (one/6.)*(zeta)*(zeta+one)*(zeta-one)
1915         END SELECT
1916       END FUNCTION shape4
1917     
1918     
1919     
1920     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1921       FUNCTION shape5(zeta,i,dx)
1922     ! Fifth-order (power of 4) shape functions
1923     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1924     
1925     !-----------------------------------------------
1926     ! Local variables
1927     !-----------------------------------------------
1928         REAL(prcn):: shape5
1929         REAL(prcn):: zeta
1930         REAL(prcn), DIMENSION(:):: dx
1931         INTEGER:: i
1932         REAL(prcn):: d1, d2, d3, d4
1933         REAL(prcn):: num, denom, zh
1934     !-----------------------------------------------
1935     
1936         SELECT CASE (i)
1937         CASE (1)
1938            d1 = -dx(1)
1939            d2 = d1 - dx(2)
1940            d3 = d2 - dx(3)
1941            d4 = d3 - dx(4)
1942            denom = d1*d2*d3*d4
1943            zh = zeta*dx(2)
1944            num = zh*(zh -dx(2))*(zh -dx(2) -dx(3)) &
1945                 *(zh -dx(2) -dx(3) -dx(4))
1946            shape5 = num/denom
1947         CASE (2)
1948            d1 =  dx(1)
1949            d2 = -dx(2)
1950            d3 =  d2 - dx(3)
1951            d4 =  d3 - dx(4)
1952            denom = d1*d2*d3*d4
1953            zh = zeta*dx(2)
1954            num = (zh +dx(1))*(zh -dx(2))*(zh -dx(2) -dx(3)) &
1955                 *(zh -dx(2) -dx(3) -dx(4))
1956            shape5 = num/denom
1957         CASE (3)
1958            d1 =  dx(1) + dx(2)
1959            d2 =  dx(2)
1960            d3 = -dx(3)
1961            d4 =  d3 - dx(4)
1962            denom = d1*d2*d3*d4
1963            zh = zeta*dx(2)
1964            num = (zh +dx(1))*(zh)*(zh -dx(2) -dx(3)) &
1965                 *(zh -dx(2) -dx(3) -dx(4))
1966            shape5 = num/denom
1967         CASE (4)
1968            d1 =  dx(1) + dx(2) + dx(3)
1969            d2 =  d1 - dx(1)
1970            d3 =  d2 - dx(2)
1971            d4 = -dx(4)
1972            denom = d1*d2*d3*d4
1973            zh = zeta*dx(2)
1974            num = (zh +dx(1))*(zh)*(zh -dx(2)) &
1975                 *(zh -dx(2) -dx(3) -dx(4))
1976            shape5 = num/denom
1977         CASE (5)
1978            d1 =  dx(1) + dx(2) + dx(3) + dx(4)
1979            d2 =  d1 - dx(1)
1980            d3 =  d2 - dx(2)
1981            d4 =  d3 - dx(3)
1982            denom = d1*d2*d3*d4
1983            zh = zeta*dx(2)
1984            num = (zh +dx(1))*(zh)*(zh -dx(2)) &
1985                 *(zh -dx(2) -dx(3))
1986            shape5 = num/denom
1987         END SELECT
1988       END FUNCTION shape5
1989     
1990     
1991     
1992     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1993       FUNCTION shape6(zeta,i,dx)
1994     ! Sixth-order (power of 5) shape functions
1995     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1996     !-----------------------------------------------
1997     ! Local variables
1998     !-----------------------------------------------
1999         REAL(prcn):: shape6
2000         REAL(prcn):: zeta
2001         REAL(prcn), DIMENSION(:):: dx
2002         INTEGER:: i
2003         REAL(prcn):: d1, d2, d3, d4, d5
2004         REAL(prcn):: num, denom, zh
2005     !-----------------------------------------------
2006     
2007         SELECT CASE (i)
2008         CASE (1)
2009            d1 = -dx(1)
2010            d2 = d1 - dx(2)
2011            d3 = d2 - dx(3)
2012            d4 = d3 - dx(4)
2013            d5 = d4 - dx(5)
2014            denom = d1*d2*d3*d4*d5
2015            zh = zeta*dx(3)
2016            num = (zh + dx(2))*(zh)*(zh - dx(3)) &
2017                 *(zh -dx(3) -dx(4))*(zh - dx(3) -dx(4) -dx(5))
2018            shape6 = num/denom
2019         CASE (2)
2020            d1 =  dx(1)
2021            d2 = -dx(2)
2022            d3 =  d2 - dx(3)
2023            d4 =  d3 - dx(4)
2024            d5 =  d4 - dx(5)
2025            denom = d1*d2*d3*d4*d5
2026            zh = zeta*dx(3)
2027            num = (zh +dx(1) +dx(2))*(zh)*(zh -dx(3)) &
2028                 *(zh -dx(3) -dx(4))*(zh - dx(3) -dx(4) -dx(5))
2029            shape6 = num/denom
2030         CASE (3)
2031            d1 =  dx(1) + dx(2)
2032            d2 =  dx(2)
2033            d3 = -dx(3)
2034            d4 =  d3 - dx(4)
2035            d5 =  d4 - dx(5)
2036            denom = d1*d2*d3*d4*d5
2037            zh = zeta*dx(3)
2038            num = (zh +dx(1) +dx(2))*(zh +dx(2))*(zh -dx(3)) &
2039                 *(zh -dx(3) -dx(4))*(zh - dx(3) -dx(4) -dx(5))
2040            shape6 = num/denom
2041         CASE (4)
2042            d1 =  dx(1) + dx(2) + dx(3)
2043            d2 =  d1 - dx(1)
2044            d3 =  d2 - dx(2)
2045            d4 = -dx(4)
2046            d5 =  d4 - dx(5)
2047            denom = d1*d2*d3*d4*d5
2048            zh = zeta*dx(3)
2049            num = (zh +dx(1) +dx(2))*(zh +dx(2))*(zh) &
2050                 *(zh -dx(3) -dx(4))*(zh - dx(3) -dx(4) -dx(5))
2051            shape6 = num/denom
2052         CASE (5)
2053            d1 =  dx(1) + dx(2) + dx(3) + dx(4)
2054            d2 =  d1 - dx(1)
2055            d3 =  d2 - dx(2)
2056            d4 =  d3 - dx(3)
2057            d5 = -dx(5)
2058            denom = d1*d2*d3*d4*d5
2059            zh = zeta*dx(3)
2060            num = (zh +dx(1) +dx(2))*(zh +dx(2))*(zh) &
2061                 *(zh -dx(3))*(zh - dx(3) -dx(4) -dx(5))
2062            shape6 = num/denom
2063         CASE (6)
2064            d1 =  dx(1) + dx(2) + dx(3) + dx(4) + dx(5)
2065            d2 =  d1 - dx(1)
2066            d3 =  d2 - dx(2)
2067            d4 =  d3 - dx(3)
2068            d5 =  d4 - dx(4)
2069            denom = d1*d2*d3*d4*d5
2070            zh = zeta*dx(3)
2071            num = (zh +dx(1) +dx(2))*(zh +dx(2))*(zh) &
2072                 *(zh -dx(3))*(zh - dx(3) -dx(4))
2073            shape6 = num/denom
2074         END SELECT
2075       END FUNCTION shape6
2076     
2077     
2078     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2079       FUNCTION shape3deriv_csi(zeta,i,dx)
2080     
2081     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2082     !-----------------------------------------------
2083     ! Local variables
2084     !-----------------------------------------------
2085         REAL(prcn):: shape3deriv_csi
2086         REAL(prcn):: zeta
2087         REAL(prcn), DIMENSION(:):: dx
2088         INTEGER:: i
2089     !-----------------------------------------------
2090     
2091     !!$    IF (zeta.GE.-two.AND.zeta.LE.-one) THEN
2092     !!$       shape3deriv_csi = (half)*(two+zeta)**2.0
2093     !!$    ELSEIF(zeta.GT.-one.AND.zeta.LE.zero) THEN
2094     !!$       shape3deriv_csi = (one/six)*(-9.0*zeta**2.0-12.0*zeta)
2095     !!$    ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2096     !!$       shape3deriv_csi = (one/six)*(9.0*zeta**2.0-12.0*zeta)
2097     !!$    ELSEIF(zeta.GT.one.AND.zeta.LE.two) THEN
2098     !!$       shape3deriv_csi = (-half)*(two-zeta)**2.0
2099     !!$    ELSE
2100     !!$       shape3deriv_csi = zero
2101     !!$    ENDIF
2102         SELECT CASE (i)
2103         CASE (1)
2104            shape3deriv_csi = -two*(1-zeta)
2105         CASE (2)
2106            shape3deriv_csi = -two*zeta+two*(1-zeta)
2107         CASE (3)
2108            shape3deriv_csi = two*zeta
2109         END SELECT
2110       END FUNCTION shape3deriv_csi
2111     
2112     
2113     
2114     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2115       FUNCTION shape4deriv_csi(zeta,i,dx)
2116     
2117     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2118     !-----------------------------------------------
2119     ! Local variables
2120     !-----------------------------------------------
2121         REAL(prcn):: shape4deriv_csi
2122         REAL(prcn):: zeta
2123         REAL(prcn), DIMENSION(:):: dx
2124         INTEGER:: i
2125     !-----------------------------------------------
2126     
2127         IF (zeta.GE.-two.AND.zeta.LE.-one) THEN
2128            shape4deriv_csi = (half)*(two+zeta)**2.0
2129         ELSEIF(zeta.GT.-one.AND.zeta.LE.zero) THEN
2130            shape4deriv_csi = (one/six)*(-9.0*zeta**2.0-12.0*zeta)
2131         ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2132            shape4deriv_csi = (one/six)*(9.0*zeta**2.0-12.0*zeta)
2133         ELSEIF(zeta.GT.one.AND.zeta.LE.two) THEN
2134            shape4deriv_csi = (-half)*(two-zeta)**2.0
2135         ELSE
2136            shape4deriv_csi = zero
2137         ENDIF
2138       END FUNCTION shape4deriv_csi
2139     
2140     
2141     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2142       FUNCTION shape4deriv(zeta,i,dx)
2143     
2144     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2145     !-----------------------------------------------
2146     ! Local variables
2147     !-----------------------------------------------
2148         REAL(prcn):: shape4deriv
2149         REAL(prcn):: zeta
2150         REAL(prcn), DIMENSION(:):: dx
2151         INTEGER:: i
2152         REAL(prcn) :: tmp
2153         REAL(prcn):: r1, r2, c1, c2, c3
2154     !-----------------------------------------------
2155     
2156         r1 = dx(2)/dx(1)
2157         r2 = dx(3)/dx(1)
2158     
2159         SELECT CASE (i)
2160         CASE (1)
2161            c1 = 1.0/(1.0 + r1)
2162            c3 = 1.0/(1.0 + r1 + r1*r2)
2163            tmp = -c1*c3*(r1**3.0)
2164            !shape4deriv = -c1*c3*(r1**3.0)*(zeta)*(zeta-1.0)*(zeta-(1.0+r2))
2165            shape4deriv = (1.0)*(zeta-1.0)*(zeta-(1.0+r2))&
2166                 + (zeta)*(1.0)*(zeta-(1.0+r2))&
2167                 +(zeta)*(zeta-1.0)*(1.0)
2168            shape4deriv = tmp*shape4deriv
2169         CASE (2)
2170            c2 = 1.0/(1.0 + r2)
2171     
2172            shape4deriv = (1.0)*(zeta*r1+1.0)*(zeta-(1.0+r2)) &
2173                 +(zeta-1.0)*(r1)*(zeta-(1.0+r2)) &
2174                 +(zeta-1.0)*(zeta*r1+1.0)*(1.0)
2175            shape4deriv = c2*shape4deriv
2176         CASE (3)
2177            c1 = 1.0/(1.0 + r1)
2178            tmp = -c1/r2
2179            shape4deriv = (1.0)*(zeta*r1+1.0)*(zeta-(1.0+r2))&
2180                 +(zeta)*(r1)*(zeta-(1.0+r2))&
2181                 +(zeta)*(zeta*r1+1.0)*(1.0)
2182            shape4deriv = tmp*shape4deriv
2183         CASE (4)
2184            c2 = 1.0/(1.0 + r2)
2185            c3 = 1.0/(1.0 + r1 + r1*r2)
2186            tmp = (c3*c2/r2)
2187            shape4deriv = (1.0)*(zeta*r1+1.0)*(zeta-1.0)&
2188                 +(zeta)*(r1)*(zeta-1.0)&
2189                 +(zeta)*(zeta*r1+1.0)*(1.0)
2190            shape4deriv = tmp*shape4deriv
2191         END SELECT
2192       END FUNCTION shape4deriv
2193     
2194     
2195     
2196     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2197       FUNCTION shape2deriv(zeta,i)
2198     
2199     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2200     !-----------------------------------------------
2201     ! Local variables
2202     !-----------------------------------------------
2203         REAL(prcn):: shape2deriv
2204         REAL(prcn):: zeta
2205         INTEGER:: i
2206     !-----------------------------------------------
2207     
2208         SELECT CASE (i)
2209         CASE (1)
2210            shape2deriv = -one!zeta
2211         CASE (2)
2212            shape2deriv = one
2213         END SELECT
2214       END FUNCTION shape2deriv
2215     
2216     
2217     
2218     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2219       FUNCTION shape4csi(zeta,i,dx,dim)
2220     
2221     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2222     !-----------------------------------------------
2223     ! Local variables
2224     !-----------------------------------------------
2225         REAL(prcn):: shape4csi
2226         REAL(prcn),INTENT(in):: zeta
2227         REAL(prcn), DIMENSION(:),INTENT(in):: dx
2228         INTEGER,INTENT(in):: i,dim
2229     !-----------------------------------------------
2230     
2231         IF (zeta.GE.-two.AND.zeta.LE.-one) THEN
2232            shape4csi = (one/six)*(two+zeta)**3.0
2233         ELSEIF(zeta.GT.-one.AND.zeta.LE.zero) THEN
2234            shape4csi = (one/six)*(-3.0*zeta**3.0-six*zeta**2.0+4.0)
2235         ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2236            shape4csi = (one/six)*(3.0*zeta**3.0-six*zeta**2.0+4.0)
2237         ELSEIF(zeta.GT.one.AND.zeta.LE.two) THEN
2238            shape4csi = (one/six)*(two-zeta)**3.0
2239         ELSE
2240            shape4csi = zero
2241            !print*,'shape4rg .... zeta=',zeta,dim
2242         ENDIF
2243     
2244       END FUNCTION shape4csi
2245     
2246     
2247     
2248     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2249       FUNCTION shape3csileft(zeta,i,dx,dim)
2250     
2251     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2252     !-----------------------------------------------
2253     ! Local variables
2254     !-----------------------------------------------
2255         REAL(prcn):: shape3csileft
2256         REAL(prcn),INTENT(in):: zeta
2257         REAL(prcn), DIMENSION(:),INTENT(in):: dx
2258         INTEGER,INTENT(in):: i,dim
2259     !-----------------------------------------------
2260     
2261         IF (zeta.GE.-one.AND.zeta.LE.zero) THEN
2262            shape3csileft = (half)*(zeta+one)**2.0
2263         ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2264            shape3csileft = (half)*(-two*zeta**2.+2.0*zeta+1.0)
2265         ELSEIF(zeta.GT.one.AND.zeta.LE.two) THEN
2266            shape3csileft = (half)*(zeta-two)**2.0
2267         ELSE
2268            shape3csileft = zero
2269            !print*,'shape4rg .... zeta=',zeta,dim
2270         ENDIF
2271     
2272     !!$    IF (zeta.GE.-two.AND.zeta.LE.-one) THEN
2273     !!$       shape3csileft = (half)*(two+zeta)**2.0
2274     !!$    ELSEIF(zeta.GT.-one.AND.zeta.LE.zero) THEN
2275     !!$       shape3csileft = (half)*(-two*zeta**2.-two*zeta+one)
2276     !!$    ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2277     !!$       shape3csileft = (half)*(one-zeta)**2.0
2278     !!$    ELSEIF(zeta.GT.one.AND.zeta.LE.two) THEN
2279     !!$       shape3csileft = (half)*(two-zeta)**2.0
2280     !!$    ELSE
2281     !!$       shape3csileft = zero
2282     !!$       !print*,'shape4rg .... zeta=',zeta,dim
2283     !!$    ENDIF
2284     !!$    SELECT CASE (i)
2285     !!$    CASE (1)
2286     !!$       shape3csileft = (1-zeta)**2.
2287     !!$       !shape4 = -(one/6.)*(zeta)*(zeta-one)*(zeta-two)
2288     !!$    CASE (2)
2289     !!$        shape3csileft = two*zeta*(1-zeta)
2290     !!$     CASE (3)
2291     !!$        shape3csileft = zeta**2.
2292     !!$     END SELECT
2293       END FUNCTION shape3csileft
2294     
2295     
2296     
2297     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2298       FUNCTION shape3csiright(zeta,i,dx,dim)
2299     
2300     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2301     !-----------------------------------------------
2302     ! Local variables
2303     !-----------------------------------------------
2304         REAL(prcn):: shape3csiright
2305         REAL(prcn),INTENT(in):: zeta
2306         REAL(prcn), DIMENSION(:),INTENT(in):: dx
2307         INTEGER,INTENT(in):: i,dim
2308     !-----------------------------------------------
2309     
2310         IF (zeta.GE.-two.AND.zeta.LE.-one) THEN
2311            shape3csiright = (half)*(two+zeta)**2.0
2312         ELSEIF(zeta.GT.-one.AND.zeta.LE.zero) THEN
2313            shape3csiright = (half)*(-two*zeta**2.-two*zeta+one)
2314         ELSEIF(zeta.GT.zero.AND.zeta.LE.one) THEN
2315            shape3csiright = (half)*(one-zeta)**2.0
2316         ELSE
2317            shape3csiright = zero
2318            !print*,'shape4rg .... zeta=',zeta,dim
2319         ENDIF
2320       END FUNCTION shape3csiright
2321     
2322     
2323     
2324     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2325       FUNCTION shape4new(pos,x,i)
2326     
2327     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2328     !-----------------------------------------------
2329     ! Local variables
2330     !-----------------------------------------------
2331         REAL(prcn):: shape4new
2332         REAL(prcn):: pos
2333         REAL(prcn), DIMENSION(:):: x
2334         INTEGER:: i
2335         REAL(prcn):: num,den
2336     !-----------------------------------------------
2337     
2338         SELECT CASE (i)
2339         CASE (1)
2340            num = (pos-x(2))*(pos - x(3))*(pos - x(4))
2341            den = (x(1) - x(2))*(x(1) - x(3))*(x(1) - x(4))
2342            shape4new = num/den
2343         CASE (2)
2344            num = (pos-x(1))*(pos - x(3))*(pos - x(4))
2345            den = (x(2) - x(1))*(x(2) - x(3))*(x(2) - x(4))
2346            shape4new = num/den
2347         CASE (3)
2348            num = (pos-x(1))*(pos - x(2))*(pos - x(4))
2349            den = (x(3) - x(1))*(x(3) - x(2))*(x(3) - x(4))
2350            shape4new = num/den
2351         CASE (4)
2352            num = (pos-x(1))*(pos - x(2))*(pos - x(3))
2353            den = (x(4) - x(1))*(x(4) - x(2))*(x(4) - x(3))
2354            shape4new = num/den
2355         END SELECT
2356       END FUNCTION shape4new
2357     
2358     
2359     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2360       SUBROUTINE set_interpolation_scheme(choice)
2361     
2362     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2363         !USE discretelement, ONLY : scheme, interp_scheme, order,ob2l,ob2r,&
2364         !     & gstencil, vstencil, sstencil, wtderivp
2365           USE param
2366         IMPLICIT NONE
2367     !-----------------------------------------------
2368     ! Local variables
2369     !-----------------------------------------------
2370         INTEGER, INTENT(in) :: choice
2371         INTEGER :: order_orig
2372     !-----------------------------------------------
2373     
2374         order_orig = order
2375         IF(choice.EQ.1) THEN
2376            interp_scheme = 'lpi'
2377            scheme = '4-order'
2378         ELSE IF(choice.EQ.2) THEN
2379            interp_scheme = 'lpi'
2380            scheme = '2-order'
2381         ELSE IF(choice.EQ.3) THEN
2382            interp_scheme = 'csi'
2383            scheme = '4-order'
2384         ENDIF
2385         SELECT CASE(scheme)
2386         CASE("2-order")
2387            order = 2
2388         CASE("3-order")
2389            order = 3
2390         CASE("4-order")
2391            order = 4
2392         CASE("5-order")
2393            order = 5
2394         CASE("6-order")
2395            order = 6
2396         END SELECT
2397     
2398     ! if ob2l is even then ob2l will equal ob2r since results after decimal are
2399     ! discarded (truncated)
2400         ob2l = (order+1)/2
2401         ob2r = order/2
2402     
2403         IF(.not.allocated(gstencil)) THEN
2404     ! max(1*(3-dimn), order*(dimn-2)) =order (in 3D) or =1 (in 2D)
2405            ALLOCATE(gstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2406         ELSEIF(ALLOCATED(gstencil).AND.order_orig.NE.order) THEN
2407            DEALLOCATE(gstencil)
2408            ALLOCATE(gstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2409         ENDIF
2410     
2411         IF(.not.allocated(vstencil)) THEN
2412            ALLOCATE(vstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2413         ELSEIF(ALLOCATED(vstencil).AND.order_orig.NE.order) THEN
2414            DEALLOCATE(vstencil)
2415            ALLOCATE(vstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2416         ENDIF
2417     
2418         IF(.not.allocated(pgradstencil)) THEN
2419            ALLOCATE(pgradstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2420         ELSEIF(ALLOCATED(pgradstencil).AND.order_orig.NE.order) THEN
2421            DEALLOCATE(pgradstencil)
2422            ALLOCATE(pgradstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2423         ENDIF
2424     
2425         IF(.not.allocated(psgradstencil)) THEN
2426            ALLOCATE(psgradstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2427         ELSEIF(ALLOCATED(psgradstencil).AND.order_orig.NE.order) THEN
2428            DEALLOCATE(psgradstencil)
2429            ALLOCATE(psgradstencil  (order,order,max(1*(3-dimn), order*(dimn-2)),3))
2430         ENDIF
2431     
2432         IF(.not.allocated(VEL_SOL_STENCIL)) THEN
2433            ALLOCATE(VEL_SOL_STENCIL  (order,order,max(1*(3-dimn), order*(dimn-2)),3, DIM_M))
2434         ELSEIF(ALLOCATED(VEL_SOL_STENCIL).AND.order_orig.NE.order) THEN
2435            DEALLOCATE(VEL_SOL_STENCIL)
2436            ALLOCATE(VEL_SOL_STENCIL  (order,order,max(1*(3-dimn), order*(dimn-2)),3, DIM_M))
2437         ENDIF
2438     
2439         IF(.not.allocated(sstencil)) THEN
2440            ALLOCATE(sstencil  (order,order,max(1*(3-dimn), order*(dimn-2))))
2441         ELSEIF(ALLOCATED(sstencil).AND.order_orig.NE.order) THEN
2442            DEALLOCATE(sstencil)
2443            ALLOCATE(sstencil  (order,order,max(1*(3-dimn), order*(dimn-2))))
2444         ENDIF
2445     
2446       END SUBROUTINE set_interpolation_scheme
2447     
2448     
2449     
2450     END MODULE interpolation
2451