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