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