New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
modinterpbasic.F90 in vendors/AGRIF/CMEMS_2020/AGRIF_FILES – NEMO

source: vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90 @ 10087

Last change on this file since 10087 was 10087, checked in by rblod, 6 years ago

update AGRIF library
  • Property svn:keywords set to Id
File size: 47.7 KB
Line 
1!
2! $Id$
3!
4!     AGRIF (Adaptive Grid Refinement In Fortran)
5!
6!     Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
7!                        Christophe Vouland (Christophe.Vouland@imag.fr)
8!
9!     This program is free software; you can redistribute it and/or modify
10!     it under the terms of the GNU General Public License as published by
11!     the Free Software Foundation; either version 2 of the License, or
12!     (at your option) any later version.
13!
14!     This program is distributed in the hope that it will be useful,
15!     but WITHOUT ANY WARRANTY; without even the implied warranty of
16!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17!     GNU General Public License for more details.
18!
19!     You should have received a copy of the GNU General Public License
20!     along with this program; if not, write to the Free Software
21!     Foundation, Inc., 59 Temple Place- Suite 330, Boston, MA 02111-1307, USA.
22!
23!
24!> Module Agrif_InterpBasic
25!>
26!> Contains different procedures of interpolation (linear,lagrange, spline,...) used in
27!! the Agrif_Interpolation module.
28!
29module Agrif_InterpBasic
30!
31    use Agrif_Types
32!
33    implicit none
34!
35    real, dimension(5,Agrif_MaxRaff,3)      :: tabppm
36    real, dimension(Agrif_MaxRaff)          :: tabdiff2, tabdiff3
37    real, dimension(:),   allocatable       :: tabtest4
38    real, dimension(:,:), allocatable       :: coeffparent
39    integer, dimension(:,:), allocatable    :: indparent
40    integer, dimension(:,:), allocatable    :: indparentppm, indchildppm
41    integer, dimension(:), allocatable      :: indparentppm_1d, indchildppm_1d
42!
43
44    private :: Agrif_limiter_vanleer
45!
46contains
47!
48!===================================================================================================
49!  subroutine Agrif_basicinterp_linear1D
50!
51!> Linear 1D interpolation on a child grid (vector y) from its parent grid (vector x).
52!---------------------------------------------------------------------------------------------------
53subroutine Agrif_basicinterp_linear1D ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
54!---------------------------------------------------------------------------------------------------
55    real, dimension(np), intent(in)     :: x            !< Coarse input data from parent
56    real, dimension(nc), intent(out)    :: y            !< Fine output data to child
57    integer,             intent(in)     :: np           !< Length of input array
58    integer,             intent(in)     :: nc           !< Length of output array
59    real(kind=8),                intent(in)     :: s_parent     !< Parent grid position (s_root = 0)
60    real(kind=8),                intent(in)     :: s_child      !< Child  grid position (s_root = 0)
61    real(kind=8),                intent(in)     :: ds_parent    !< Parent grid dx (ds_root = 1)
62    real(kind=8),                intent(in)     :: ds_child     !< Child  grid dx (ds_root = 1)
63!
64    integer :: i, coeffraf, locind_parent_left
65    real(kind=8)    :: globind_parent_left, globind_parent_right
66    real(kind=8)    :: invds, invds2, ypos, ypos2, diff
67!
68    coeffraf = nint(ds_parent/ds_child)
69!
70    if ( coeffraf == 1 ) then
71        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
72        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
73        return
74    endif
75!
76    ypos = s_child
77    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
78    globind_parent_left  = s_parent + (locind_parent_left - 1)*ds_parent
79    globind_parent_right = globind_parent_left + ds_parent
80!
81    invds  = 1./ds_parent
82    invds2 = ds_child/ds_parent
83    ypos2 = ypos*invds
84    globind_parent_right = globind_parent_right*invds
85!
86    do i = 1,nc-1
87!
88        if (ypos2 > globind_parent_right) then
89            locind_parent_left = locind_parent_left + 1
90            globind_parent_right = globind_parent_right + 1.
91            ypos2 = ypos*invds+(i-1)*invds2
92        endif
93!
94        diff = globind_parent_right - ypos2
95! quick fix for roundoff error
96        diff=nint(diff*coeffraf)/real(coeffraf)
97
98        y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1))
99
100        ypos2 = ypos2 + invds2
101!
102    enddo
103!
104    ypos = s_child + (nc-1)*ds_child
105    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
106!
107    if (locind_parent_left == np) then
108        y(nc) = x(np)
109    else
110        globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent
111        diff=(globind_parent_left + ds_parent - ypos)*invds
112
113! quick fix for roundoff error
114        diff=nint(diff*coeffraf)/real(coeffraf)
115!        y(nc) = ((globind_parent_left + ds_parent - ypos)*x(locind_parent_left)  &
116!                           + (ypos - globind_parent_left)*x(locind_parent_left+1))*invds
117        y(nc) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1))
118    endif
119
120!---------------------------------------------------------------------------------------------------
121end subroutine Agrif_basicinterp_linear1D
122!===================================================================================================
123!
124!===================================================================================================
125!  subroutine Linear1dPrecompute2d
126!
127!> Computes 2D coefficients and index for a linear 1D interpolation on a child grid (vector y)
128!! from its parent grid (vector x).
129!---------------------------------------------------------------------------------------------------
130subroutine Linear1dPrecompute2d ( np2, np, nc, s_parent, s_child, ds_parent, ds_child, dir )
131!---------------------------------------------------------------------------------------------------
132    integer, intent(in) :: np,nc,np2
133    real(kind=8),    intent(in) :: s_parent, s_child
134    real(kind=8),    intent(in) :: ds_parent, ds_child
135    integer, intent(in) :: dir
136!
137    integer :: i,coeffraf,locind_parent_left,inc,inc1,inc2
138    integer, dimension(:,:), allocatable :: indparent_tmp
139    real, dimension(:,:), allocatable :: coeffparent_tmp
140    real(kind=8)    :: ypos,globind_parent_left,globind_parent_right
141    real(kind=8)    :: invds, invds2, invds3
142    real(kind=8) :: ypos2,diff
143!
144    coeffraf = nint(ds_parent/ds_child)
145!
146    ypos = s_child
147    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
148    globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent
149    globind_parent_right = globind_parent_left + ds_parent
150!
151    invds = 1./ds_parent
152    invds2 = ds_child/ds_parent
153    invds3 = 0.5/real(coeffraf)
154    ypos2 = ypos*invds
155    globind_parent_right=globind_parent_right*invds
156!
157    if (.not.allocated(indparent)) then
158        allocate(indparent(nc*np2,3),coeffparent(nc*np2,3))
159    else
160        if ( size(indparent,1) < nc*np2 ) then
161            allocate(coeffparent_tmp(size(indparent,1),size(indparent,2)))
162            allocate(  indparent_tmp(size(indparent,1),size(indparent,2)))
163            coeffparent_tmp = coeffparent
164            indparent_tmp   = indparent
165            deallocate(indparent,coeffparent)
166            allocate(indparent(nc*np2,3),coeffparent(nc*np2,3))
167            coeffparent(1:size(coeffparent_tmp,1),1:size(coeffparent_tmp,2)) = coeffparent_tmp
168            indparent(  1:size(indparent_tmp,  1),1:size(indparent_tmp,  2)) = indparent_tmp
169            deallocate(indparent_tmp,coeffparent_tmp)
170        endif
171    endif
172!
173    do i = 1,nc-1
174!
175        if (ypos2 > globind_parent_right) then
176            locind_parent_left = locind_parent_left + 1
177            globind_parent_right = globind_parent_right + 1.d0
178            ypos2 = ypos*invds+(i-1)*invds2
179        endif
180!
181        diff = globind_parent_right - ypos2
182        diff = invds3*nint(2*coeffraf*diff)
183        indparent(i,dir) = locind_parent_left
184        coeffparent(i,dir) = diff
185        ypos2 = ypos2 + invds2
186!
187    enddo
188!
189    ypos = s_child + (nc-1)*ds_child
190    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
191
192    if (locind_parent_left == np) then
193        indparent(nc,dir) = locind_parent_left-1
194        coeffparent(nc,dir) = 0.
195    else
196        globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent
197        indparent(nc,dir) = locind_parent_left
198        diff = (globind_parent_left + ds_parent - ypos) * invds
199        diff = invds3*nint(2*coeffraf*diff)
200        coeffparent(nc,dir) = diff
201    endif
202
203    do i=2, np2
204        inc  =  i*nc
205        inc1 = (i-1)*nc
206        inc2 = (i-2)*nc
207!CDIR ALTCODE
208        indparent(1+inc1:inc,dir) = indparent(1+inc2:inc1,dir)+np
209!CDIR ALTCODE
210        coeffparent(1+inc1:inc,dir) =coeffparent(1:nc,dir)
211    enddo
212!---------------------------------------------------------------------------------------------------
213end subroutine Linear1dPrecompute2d
214!===================================================================================================
215!
216!===================================================================================================
217!  subroutine Linear1dAfterCompute
218!
219!> Carries out a linear 1D interpolation on a child grid (vector y) from its parent grid (vector x)
220!! using precomputed coefficient and index.
221!---------------------------------------------------------------------------------------------------
222subroutine Linear1dAfterCompute ( x, y, np, nc, dir )
223!---------------------------------------------------------------------------------------------------
224    integer,             intent(in)     :: np, nc
225    real, dimension(np), intent(in)     :: x
226    real, dimension(nc), intent(out)    :: y
227    integer,             intent(in)     :: dir
228!
229    integer :: i
230!
231!CDIR ALTCODE
232!CDIR NODEP
233    do i = 1,nc
234        y(i) = coeffparent(i,dir)  * x(MAX(indparent(i,dir),1)) + &
235           (1.-coeffparent(i,dir)) * x(indparent(i,dir)+1)
236    enddo
237!---------------------------------------------------------------------------------------------------
238end subroutine Linear1dAfterCompute
239!===================================================================================================
240!
241!===================================================================================================
242!  subroutine Lagrange1d
243!
244!> Carries out a lagrange 1D interpolation on a child grid (vector y) from its parent grid
245!! (vector x).
246!---------------------------------------------------------------------------------------------------
247subroutine Lagrange1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
248!---------------------------------------------------------------------------------------------------
249    integer,             intent(in)     :: np, nc
250    real, dimension(np), intent(in)     :: x
251    real, dimension(nc), intent(out)    :: y
252    real(kind=8),                intent(in)     :: s_parent, s_child
253    real(kind=8),                intent(in)     :: ds_parent, ds_child
254!
255    integer :: i, coeffraf, locind_parent_left
256    real(kind=8)    :: ypos,globind_parent_left
257    real(kind=8)    :: deltax, invdsparent
258    real    :: t2,t3,t4,t5,t6,t7,t8
259!
260    if (np <= 2) then
261        call Agrif_basicinterp_linear1D(x,y,np,nc,s_parent,s_child,ds_parent,ds_child)
262        return
263    endif
264!
265    coeffraf = nint(ds_parent/ds_child)
266!
267    if (coeffraf == 1) then
268        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
269        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
270        return
271    endif
272!
273    invdsparent = 1./ds_parent
274    ypos = s_child
275!
276    do i = 1,nc
277!
278        locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
279        globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent
280
281        deltax = invdsparent*(ypos-globind_parent_left)
282        deltax = nint(coeffraf*deltax)/real(coeffraf)
283
284        ypos = ypos + ds_child
285        if (abs(deltax) <= 0.0001) then
286            y(i)=x(locind_parent_left)
287            cycle
288        endif
289!
290        t2 = deltax - 2.
291        t3 = deltax - 1.
292        t4 = deltax + 1.
293
294        t5 = -(1./6.)*deltax*t2*t3
295        t6 = 0.5*t2*t3*t4
296        t7 = -0.5*deltax*t2*t4
297        t8 = (1./6.)*deltax*t3*t4
298
299        y(i) = t5*x(locind_parent_left-1) + t6*x(locind_parent_left)    &
300              +t7*x(locind_parent_left+1) + t8*x(locind_parent_left+2)
301!
302    enddo
303!---------------------------------------------------------------------------------------------------
304end subroutine Lagrange1d
305!===================================================================================================
306!
307!===================================================================================================
308!  subroutine Constant1d
309!
310!> Carries out a constant 1D interpolation on a child grid (vector y) from its parent grid (vector x).
311!---------------------------------------------------------------------------------------------------
312subroutine Constant1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
313!---------------------------------------------------------------------------------------------------
314    integer,             intent(in)     :: np, nc
315    real, dimension(np), intent(in)     :: x
316    real, dimension(nc), intent(out)    :: y
317    real(kind=8),                intent(in)     :: s_parent, s_child
318    real(kind=8),                intent(in)     :: ds_parent, ds_child
319!
320    integer :: i, coeffraf, locind_parent
321    real(kind=8)    :: ypos
322!
323    coeffraf = nint(ds_parent/ds_child)
324!
325    if (coeffraf == 1) then
326        locind_parent = 1 + nint((s_child - s_parent)/ds_parent)
327        y(1:nc) = x(locind_parent:locind_parent+nc-1)
328        return
329    endif
330!
331    ypos = s_child
332!
333    do i = 1,nc
334!
335        locind_parent = 1 + nint((ypos - s_parent)/ds_parent)
336        y(i) = x(locind_parent)
337        ypos = ypos + ds_child
338!
339    enddo
340!---------------------------------------------------------------------------------------------------
341end subroutine Constant1d
342!===================================================================================================
343!
344!===================================================================================================
345!  subroutine Linear1dConserv
346!
347!> Carries out a conservative linear 1D interpolation on a child grid (vector y) from its parent
348!! grid (vector x).
349!---------------------------------------------------------------------------------------------------
350subroutine Linear1dConserv ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
351!---------------------------------------------------------------------------------------------------
352    integer,             intent(in)     :: np, nc
353    real, dimension(np), intent(in)     :: x
354    real, dimension(nc), intent(out)    :: y
355    real(kind=8),                intent(in)     :: s_parent, s_child
356    real(kind=8),                intent(in)     :: ds_parent, ds_child
357!
358    real, dimension(:), allocatable :: ytemp
359    integer :: i,coeffraf,locind_parent_left,locind_parent_last
360    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope
361    integer :: i1,i2,ii
362    integer :: diffmod
363!
364    coeffraf = nint(ds_parent/ds_child)
365!
366    if (coeffraf == 1) then
367        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
368        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
369        return
370    endif
371!
372    diffmod = 0
373    if (mod(coeffraf,2) == 0) diffmod = 1
374
375    xdiffmod = real(diffmod)/2.
376
377    allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
378!
379    ypos = s_child
380!
381    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
382    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
383
384    xpmin = s_parent + (locind_parent_left-1)*ds_parent
385    xpmax = s_parent + (locind_parent_last-1)*ds_parent
386
387    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
388    i2 = 1+agrif_int((xpmax-s_child)/ds_child)
389
390    i = i1
391
392    if (locind_parent_left == 1) then
393        slope = (x(locind_parent_left+1)-x(locind_parent_left))/(coeffraf)
394    else
395        slope = (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
396    endif
397
398    do ii = i-coeffraf/2+diffmod,i+coeffraf/2
399        ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
400    enddo
401
402    locind_parent_left = locind_parent_left + 1
403
404    do i = i1+coeffraf, i2-coeffraf,coeffraf
405        slope = (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
406        do ii = i-coeffraf/2+diffmod,i+coeffraf/2
407            ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
408        enddo
409        locind_parent_left = locind_parent_left + 1
410    enddo
411
412    i = i2
413
414    if (locind_parent_left == np) then
415        slope = (x(locind_parent_left)-x(locind_parent_left-1))/(coeffraf)
416    else
417        slope = (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
418    endif
419
420    do ii = i-coeffraf/2+diffmod,nc
421        ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
422    enddo
423!
424    y(1:nc)=ytemp(1:nc)
425!
426    deallocate(ytemp)
427!---------------------------------------------------------------------------------------------------
428end subroutine Linear1dConserv
429!===================================================================================================
430!
431!===================================================================================================
432!  subroutine Linear1dConservLim
433!
434!> Carries out a limited conservative linear 1D interpolation on a child grid (vector y) from
435!! its parent grid (vector x).
436!---------------------------------------------------------------------------------------------------
437subroutine Linear1dConservLim ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
438!---------------------------------------------------------------------------------------------------
439    integer,             intent(in)     :: np, nc
440    real, dimension(np), intent(in)     :: x
441    real, dimension(nc), intent(out)    :: y
442    real(kind=8),                intent(in)     :: s_parent, s_child
443    real(kind=8),                intent(in)     :: ds_parent, ds_child
444!
445    real, dimension(:), allocatable :: ytemp
446    integer :: i,coeffraf,locind_parent_left,locind_parent_last
447    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope
448    integer :: i1,i2,ii
449    integer :: diffmod
450!
451    coeffraf = nint(ds_parent/ds_child)
452!
453    if (coeffraf == 1) then
454        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
455        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
456        return
457    endif
458!
459    if (coeffraf /= 3) then
460        print *,'Linear1dConservLim not ready for refinement ratio = ', coeffraf
461        stop
462    endif
463!
464    diffmod = 0
465    if (mod(coeffraf,2) == 0) diffmod = 1
466
467    xdiffmod = real(diffmod)/2.
468
469    allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
470!
471    ypos = s_child
472!
473    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
474    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
475
476    xpmin = s_parent + (locind_parent_left-1)*ds_parent
477    xpmax = s_parent + (locind_parent_last-1)*ds_parent
478
479    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
480    i2 = 1+agrif_int((xpmax-s_child)/ds_child)
481
482    i = i1
483
484    if (locind_parent_left == 1) then
485        slope=0.
486    else
487        slope = Agrif_limiter_vanleer(x(locind_parent_left-1:locind_parent_left+1))
488        slope = slope / coeffraf
489    endif
490
491    do ii = i-coeffraf/2+diffmod,i+coeffraf/2
492        ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
493    enddo
494
495    locind_parent_left = locind_parent_left + 1
496
497    do i = i1+coeffraf, i2-coeffraf,coeffraf
498        slope = Agrif_limiter_vanleer(x(locind_parent_left-1:locind_parent_left+1))
499        slope = slope / coeffraf
500        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
501            ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
502        enddo
503        locind_parent_left = locind_parent_left + 1
504    enddo
505
506    i = i2
507
508    if (locind_parent_left == np) then
509        slope=0.
510    else
511        slope = Agrif_limiter_vanleer(x(locind_parent_left-1:locind_parent_left+1))
512        slope = slope / coeffraf
513    endif
514
515    do ii=i-coeffraf/2+diffmod,nc
516        ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
517    enddo
518!
519    y(1:nc) = ytemp(1:nc)
520!
521    deallocate(ytemp)
522!---------------------------------------------------------------------------------------------------
523end subroutine Linear1dConservLim
524!===================================================================================================
525!
526!===================================================================================================
527!  subroutine PPM1d
528!
529!> Carries out a 1D interpolation and apply monotonicity constraints using piecewise parabolic
530!! method (PPM) on a child grid (vector y) from its parent grid (vector x).
531!---------------------------------------------------------------------------------------------------
532subroutine PPM1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
533!---------------------------------------------------------------------------------------------------
534    integer,             intent(in)     :: np, nc
535    real, dimension(np), intent(in)     :: x
536    real, dimension(nc), intent(out)    :: y
537    real(kind=8),                intent(in)     :: s_parent, s_child
538    real(kind=8),                intent(in)     :: ds_parent, ds_child
539!
540    integer :: i,coeffraf,locind_parent_left,locind_parent_last
541    integer :: iparent,ipos,pos,nmin,nmax
542    real(kind=8)    :: ypos
543    integer :: i1,jj
544    real(kind=8) :: xpmin
545    real :: a
546!
547    real, dimension(np) :: xl,delta,a6,slope
548    integer :: diffmod
549    real :: invcoeffraf
550!
551    coeffraf = nint(ds_parent/ds_child)
552!
553    if (coeffraf == 1) then
554        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
555!CDIR ALTCODE
556!CDIR SHORTLOOP
557        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
558        return
559    endif
560!
561    invcoeffraf = ds_child/ds_parent
562!
563    if( .not. allocated(tabtest4) ) then
564        allocate(tabtest4(-2*coeffraf:nc+2*coeffraf))
565    else
566        if (size(tabtest4) < nc+4*coeffraf+1) then
567            deallocate( tabtest4 )
568            allocate(tabtest4(-2*coeffraf:nc+2*coeffraf))
569        endif
570    endif
571!
572    ypos = s_child
573!
574    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
575    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1)*ds_child - s_parent)/ds_parent)
576!
577    xpmin = s_parent + (locind_parent_left-1)*ds_parent
578    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
579!
580!CDIR NOVECTOR
581    do i=1,coeffraf
582        tabdiff2(i) = (real(i)-0.5)*invcoeffraf
583    enddo
584!
585    a = invcoeffraf**2
586    tabdiff3(1) = (1./3.)*a
587    a = 2.*a
588!CDIR NOVECTOR
589    do i=2,coeffraf
590        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a
591    enddo
592!
593    if      ( locind_parent_last+2 <= np ) then
594        nmax = locind_parent_last+2
595    else if ( locind_parent_last+1 <= np ) then
596        nmax = locind_parent_last+1
597    else
598        nmax = locind_parent_last
599    endif
600!
601    if (locind_parent_left-1 >= 1) then
602        nmin = locind_parent_left-1
603    else
604        nmin = locind_parent_left
605    endif
606!
607!CDIR ALTCODE
608!CDIR SHORTLOOP
609    do i = nmin,nmax
610        slope(i) = x(i) - x(i-1)
611    enddo
612!
613!CDIR ALTCODE
614!CDIR SHORTLOOP
615    do i = nmin+1,nmax-1
616        xl(i)= 0.5*(x(i-1)+x(i))-0.08333333333333*(slope(i+1)-slope(i-1))
617    enddo
618!
619! apply parabolic monotonicity
620!CDIR ALTCODE
621!CDIR SHORTLOOP
622    do i = locind_parent_left,locind_parent_last
623        delta(i) = xl(i+1) - xl(i)
624        a6(i) = 6.*x(i)-3.*(xl(i) +xl(i+1))
625    enddo
626!
627    diffmod = 0
628    if (mod(coeffraf,2) == 0) diffmod = 1
629!
630    ipos = i1
631!
632    do iparent = locind_parent_left,locind_parent_last
633        pos=1
634!CDIR ALTCODE
635!CDIR SHORTLOOP
636        do jj = ipos-coeffraf/2+diffmod,ipos+coeffraf/2
637!
638            tabtest4(jj) = xl(iparent) + tabdiff2(pos) *  (delta(iparent)+a6(iparent)) &
639                                       - tabdiff3(pos) *  a6(iparent)
640            pos = pos+1
641        enddo
642        ipos = ipos + coeffraf
643    enddo
644!
645!CDIR ALTCODE
646!CDIR SHORTLOOP
647    y(1:nc) = tabtest4(1:nc)
648!---------------------------------------------------------------------------------------------------
649end subroutine PPM1d
650!===================================================================================================
651!
652!===================================================================================================
653!  subroutine PPM1dPrecompute2d
654!
655!> Computes 2D coefficients and index for a 1D interpolation using piecewise parabolic method
656!---------------------------------------------------------------------------------------------------
657subroutine PPM1dPrecompute2d ( np2, np, nc, s_parent, s_child, ds_parent, ds_child, dir )
658!---------------------------------------------------------------------------------------------------
659    integer,             intent(in)     :: np2, np, nc
660    real(kind=8),                intent(in)     :: s_parent, s_child
661    real(kind=8),                intent(in)     :: ds_parent, ds_child
662    integer,             intent(in)     :: dir
663!
664    integer, dimension(:,:), allocatable :: indparent_tmp
665    integer, dimension(:,:), allocatable :: indchild_tmp
666    integer :: i,coeffraf,locind_parent_left,locind_parent_last
667    integer :: iparent,ipos,pos
668    real    :: ypos
669    integer :: i1,jj,k,l,j
670    real(kind=8) :: xpmin
671    real :: a
672!
673    integer :: diffmod
674    real :: invcoeffraf
675!
676    coeffraf = nint(ds_parent/ds_child)
677!
678    invcoeffraf = ds_child/ds_parent
679!
680    if (.not.allocated(indparentppm)) then
681        allocate(indparentppm(np2*nc,3),indchildppm(np2*nc,3))
682    else
683        if (size(indparentppm,1) < np2*nc) then
684            allocate(   &
685                indparent_tmp(size(indparentppm,1),size(indparentppm,2)), &
686                indchild_tmp( size(indparentppm,1),size(indparentppm,2)))
687            indparent_tmp = indparentppm
688            indchild_tmp  = indchildppm
689            deallocate(indparentppm,indchildppm)
690            allocate(indparentppm(np2*nc,3),indchildppm(np2*nc,3))
691            indparentppm(1:size(indparent_tmp,1),1:size(indparent_tmp,2)) = indparent_tmp
692            indchildppm( 1:size(indparent_tmp,1),1:size(indparent_tmp,2)) = indchild_tmp
693            deallocate(indparent_tmp,indchild_tmp)
694        endif
695    endif
696
697    if (.not.allocated(indparentppm_1d)) then
698        allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf),    &
699                 indchildppm_1d( -2*coeffraf:nc+2*coeffraf))
700    else
701        if (size(indparentppm_1d) < nc+4*coeffraf+1) then
702            deallocate(indparentppm_1d,indchildppm_1d)
703            allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf),&
704                     indchildppm_1d( -2*coeffraf:nc+2*coeffraf))
705        endif
706    endif
707!
708    ypos = s_child
709!
710    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
711    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1)*ds_child - s_parent)/ds_parent)
712!
713    xpmin = s_parent + (locind_parent_left-1)*ds_parent
714    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
715!
716    do i = 1,coeffraf
717        tabdiff2(i)=(real(i)-0.5)*invcoeffraf
718    enddo
719!
720    a = invcoeffraf**2
721    tabdiff3(1) = (1./3.)*a
722    a = 2.*a
723!CDIR ALTCODE
724    do i = 2,coeffraf
725        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a
726    enddo
727
728!CDIR ALTCODE
729    do i = 1,coeffraf
730        tabppm(1,i,dir) = 0.08333333333333*(-1.+4*tabdiff2(i)-3*tabdiff3(i))
731        tabppm(2,i,dir) = 0.08333333333333*(7.-26.*tabdiff2(i)+18.*tabdiff3(i))
732        tabppm(3,i,dir) = 0.08333333333333*(7.+30*tabdiff2(i)-30*tabdiff3(i))
733        tabppm(4,i,dir) = 0.08333333333333*(-1.-10.*tabdiff2(i)+18.*tabdiff3(i))
734        tabppm(5,i,dir) = 0.08333333333333*(2*tabdiff2(i)-3*tabdiff3(i))
735    enddo
736!
737    diffmod = 0
738    if (mod(coeffraf,2) == 0) diffmod = 1
739!
740    ipos = i1
741!
742    do iparent = locind_parent_left,locind_parent_last
743        pos=1
744!CDIR ALTCODE
745        do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
746            indparentppm_1d(jj) = iparent-2
747            indchildppm_1d(jj)  = pos
748            pos = pos+1
749        enddo
750        ipos = ipos + coeffraf
751    enddo
752!
753
754    k=1
755    l=0
756    do i=1,np2
757      do j=1,nc
758        indchildppm(k,dir) = indchildppm_1d(j)
759        indparentppm(k,dir) = indparentppm_1d(j) + l
760        k=k+1
761      enddo
762      l=l+np
763    enddo
764!---------------------------------------------------------------------------------------------------
765end subroutine PPM1dPrecompute2d
766!===================================================================================================
767!
768
769!===================================================================================================
770!
771!===================================================================================================
772!  subroutine PPM1dAfterCompute
773!
774! Carries out a 1D interpolation and apply monotonicity constraints using piecewise parabolic
775! method (PPM) on a child grid (vector y) from its parent grid (vector x).
776! Use precomputed coefficient and index.
777!---------------------------------------------------------------------------------------------------
778subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc )
779!---------------------------------------------------------------------------------------------------
780    integer,             intent(in)     :: np, nc
781    real, dimension(np), intent(in)     :: x
782    real, dimension(nc), intent(out)    :: y
783    integer,             intent(in)     :: dir
784    integer, dimension(1:),intent(in) :: indchildppmloc
785    real, dimension(1:,1:),intent(in) :: tabppmloc
786!
787    integer :: i,j,jp1,jp2,jp3,jp4,k
788
789   
790    do i = 1,nc
791        j = indparentppm(i,dir)
792        jp1=j+1
793        jp2=jp1+1
794        jp3=jp2+1
795        jp4=jp3+1
796        y(i) = tabppmloc(1,indchildppmloc(i)) * x(j) + &
797               tabppmloc(2,indchildppmloc(i)) * x(jp1) + &
798               tabppmloc(3,indchildppmloc(i)) * x(jp2) + &
799               tabppmloc(4,indchildppmloc(i)) * x(jp3) + &
800               tabppmloc(5,indchildppmloc(i)) * x(jp4)
801    enddo
802
803!---------------------------------------------------------------------------------------------------
804end subroutine PPM1dAfterCompute
805
806!===================================================================================================
807!
808!===================================================================================================
809!  subroutine weno1d
810!
811! Carries out a 1D interpolation and apply monotonicity constraints
812! using WENO method on a child grid (vector y) from its parent grid (vector x).
813!---------------------------------------------------------------------------------------------------
814!subroutine weno1dnew ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
815!---------------------------------------------------------------------------------------------------
816!
817!CC   Description:
818!CC   subroutine to do a 1D interpolation and apply monotonicity constraints
819!CC   using piecewise parabolic method
820!CC   on a child grid (vector y) from its parent grid (vector x).
821!C    Method:
822!
823!     Declarations:
824!
825!      Implicit none
826!
827!     Arguments
828!      Integer             :: np,nc
829!      Real, Dimension(np) :: x
830!      Real, Dimension(nc) :: y
831!      Real, Dimension(:),Allocatable :: ytemp
832!      Real                :: s_parent,s_child,ds_parent,ds_child
833!
834!     Local scalars
835!      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
836!      Integer :: iparent,ipos,pos,nmin,nmax
837!      Real    :: ypos
838!      integer :: i1,jj
839!      Real :: xpmin,cavg,a,b
840!
841!      Real :: xrmin,xrmax,am3,s2,s1
842!      Real, Dimension(np) :: xr,xl,delta,a6,slope,slope2,smooth
843!      Real, Dimension(:),Allocatable  :: diff,diff2,diff3
844!      INTEGER :: diffmod
845!      REAL :: invcoeffraf
846!      integer :: s,l,k
847!      integer :: etan, etap
848!      real :: delta0, delta1, delta2
849!      real :: epsilon
850!      parameter (epsilon = 1.D-8)
851!      real, dimension(:,:), allocatable :: ak, ck
852!
853!      coeffraf = nint(ds_parent/ds_child)
854!
855!      If (coeffraf == 1) Then
856!          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
857!          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
858!          return
859!      End If
860!      invcoeffraf = ds_child/ds_parent
861!      Allocate(ak(0:1,coeffraf))
862!      Allocate(ck(0:1,coeffraf))
863!
864!
865!      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
866!      ypos = s_child
867!
868!      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
869!      locind_parent_last = 1 +&
870!      agrif_ceiling((ypos +(nc - 1)&
871!      *ds_child - s_parent)/ds_parent)
872!
873!      xpmin = s_parent + (locind_parent_left-1)*ds_parent
874!      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
875!
876!      Allocate( diff(coeffraf),diff2(coeffraf),diff3(coeffraf) )
877!
878!      diff(1)=0.5*invcoeffraf
879!      do i=2,coeffraf
880!         diff(i) = diff(i-1)+invcoeffraf
881!      enddo
882!
883!      ak = 0.
884!      ck = 0.
885!
886!      do i=1,coeffraf
887!         do k=0,1
888!         do s=0,2
889!         do l=0,2
890!           if (l /= s) then
891!            ak(k,i) = ak(k,i)+(diff(i)-(k-l+1.))
892!           endif
893!         enddo
894!         enddo
895!         enddo
896!
897!         etap = 0
898!         etan = 0
899!         do k=0,1
900!          if (ak(k,i) > 0) then
901!            etap = etap+1
902!          else if (ak(k,i) < 0) then
903!            etan = etan + 1
904!          endif
905!         enddo
906!
907!         do k=0,1
908!           if (ak(k,i) == 0) then
909!            Ck(k,i) = 1.
910!           else if (ak(k,i) > 0) then
911!            Ck(k,i) = 1./(etap * ak(k,i))
912!           else
913!            Ck(k,i) = -1./(etan * ak(k,i))
914!           endif
915!         enddo
916!      enddo
917!
918!
919!      a = 0.
920!      b = invcoeffraf
921!      Do i=1,coeffraf
922!         diff2(i) = 0.5*(b*b - a*a)
923!         diff3(i) = (1./3.)*(b*b*b - a*a*a)
924!         a = a + invcoeffraf
925!         b = b + invcoeffraf
926!      End do
927!
928!      if( locind_parent_last+2 <= np ) then
929!           nmax = locind_parent_last+2
930!      elseif( locind_parent_last+1 <= np ) then
931!           nmax = locind_parent_last+1
932!      else
933!           nmax = locind_parent_last
934!      endif
935!
936!      if(locind_parent_left-2 >= 1) then
937!          nmin = locind_parent_left-2
938!      elseif(locind_parent_left-1 >= 1) then
939!          nmin = locind_parent_left-1
940!      else
941!          nmin = locind_parent_left
942!      endif
943!
944!      Do i = nmin+1,nmax
945!         slope(i) = (x(i) - x(i-1))
946!      Enddo
947!      DO i=nmin+2,nmax
948!        smooth(i) = 0.5*(slope(i)**2+slope(i-1)**2)&
949!       +(slope(i)-slope(i-1))**2
950!      enddo
951!
952!        diffmod = 0
953!        IF (mod(coeffraf,2) == 0) diffmod = 1
954!
955!        ipos = i1
956!
957!        Do iparent = locind_parent_left,locind_parent_last
958!             pos=1
959!
960!            delta0=1./(epsilon+smooth(iparent  ))**3
961!            delta1=1./(epsilon+smooth(iparent+1))**3
962!            delta2=1./(epsilon+smooth(iparent+2))**3
963!
964!             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
965!
966!               pos = pos+1
967!             End do
968!             ipos = ipos + coeffraf
969!
970!        End do
971!
972!
973!        y(1:nc)=ytemp(1:nc)
974!        deallocate(ytemp)
975!        deallocate(diff, diff2, diff3)
976!
977!        deallocate(ak,ck)
978!
979!      Return
980!      End subroutine weno1dnew
981!===================================================================================================
982!
983!===================================================================================================
984!  subroutine WENO1d
985!
986!> Carries out a a 1D interpolation using WENO method on a child grid (vector y) from its parent
987!! grid (vector x).
988!---------------------------------------------------------------------------------------------------
989subroutine WENO1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
990!---------------------------------------------------------------------------------------------------
991    integer,             intent(in)     :: np, nc
992    real, dimension(np), intent(in)     :: x
993    real, dimension(nc), intent(out)    :: y
994    real(kind=8),                intent(in)     :: s_parent, s_child
995    real(kind=8),                intent(in)     :: ds_parent, ds_child
996!
997    real, dimension(:), allocatable :: ytemp
998    integer :: i,coeffraf,locind_parent_left,locind_parent_last
999    integer :: iparent,ipos,pos,nmin,nmax
1000    real(kind=8)    :: ypos
1001    integer :: i1,jj
1002    real(kind=8) :: xpmin
1003!
1004    real, dimension(np) :: slope
1005    real, dimension(:), allocatable  :: diff
1006    integer :: diffmod
1007    real :: invcoeffraf
1008    real :: delta0, delta1, sumdelta
1009    real, parameter :: epsilon = 1.d-8
1010!
1011    coeffraf = nint(ds_parent/ds_child)
1012!
1013    if (coeffraf == 1) then
1014        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
1015        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
1016        return
1017    endif
1018!
1019    invcoeffraf = ds_child/ds_parent
1020!
1021    allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
1022    ypos = s_child
1023!
1024    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
1025    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
1026!
1027    xpmin = s_parent + (locind_parent_left-1)*ds_parent
1028    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
1029!
1030    allocate(diff(coeffraf))
1031    diff(1) = 0.5*invcoeffraf
1032    do i = 2,coeffraf
1033        diff(i) = diff(i-1)+invcoeffraf
1034    enddo
1035!
1036    if ( locind_parent_last+2 <= np ) then
1037        nmax = locind_parent_last+2
1038    else if ( locind_parent_last+1 <= np ) then
1039        nmax = locind_parent_last+1
1040      else
1041        nmax = locind_parent_last
1042    endif
1043!
1044    if(locind_parent_left-1 >= 1) then
1045        nmin = locind_parent_left-1
1046    else
1047        nmin = locind_parent_left
1048    endif
1049!
1050    do i = nmin+1,nmax
1051        slope(i) = x(i) - x(i-1)
1052    enddo
1053!
1054    diffmod = 0
1055    if (mod(coeffraf,2) == 0) diffmod = 1
1056!
1057    ipos = i1
1058!
1059    do iparent = locind_parent_left,locind_parent_last
1060        pos=1
1061        delta0 = 1./(epsilon+slope(iparent  )**2)**2
1062        delta1 = 1./(epsilon+slope(iparent+1)**2)**2
1063        sumdelta = 1./(delta0+delta1)
1064        do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
1065            ytemp(jj) = x(iparent)+(diff(pos)-0.5)*( delta0*slope(iparent) + &
1066                                                     delta1*slope(iparent+1))*sumdelta
1067            pos = pos+1
1068        enddo
1069        ipos = ipos + coeffraf
1070    enddo
1071!
1072    y(1:nc) = ytemp(1:nc)
1073    deallocate(ytemp)
1074    deallocate(diff)
1075!---------------------------------------------------------------------------------------------------
1076end subroutine WENO1d
1077!===================================================================================================
1078!
1079!===================================================================================================
1080!  subroutine ENO1d
1081!
1082!> Carries out a 1D interpolation using piecewise polynomial ENO reconstruction technique
1083!! on a child grid (vector y) from its parent grid (vector x).
1084!! \see ---- p 163-164 Computational gasdynamics ----
1085!---------------------------------------------------------------------------------------------------
1086subroutine ENO1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
1087!---------------------------------------------------------------------------------------------------
1088    integer,             intent(in)     :: np, nc
1089    real, dimension(np), intent(in)     :: x
1090    real, dimension(nc), intent(out)    :: y
1091    real(kind=8),                intent(in)     :: s_parent, s_child
1092    real(kind=8),                intent(in)     :: ds_parent, ds_child
1093!
1094    integer :: i,coeffraf,locind_parent_left,locind_parent_last
1095    integer :: ipos, pos
1096    real(kind=8)    :: ypos,xi
1097    integer :: i1,jj
1098    real(kind=8) :: xpmin
1099!
1100    real, dimension(:),   allocatable  :: ytemp
1101    real, dimension(:,:), allocatable  :: xbar
1102    real, dimension(1:np+1)            :: xhalf
1103    real, dimension(3,np)              :: dd, c
1104    integer :: diffmod, left
1105!
1106    coeffraf = nint(ds_parent/ds_child)
1107!
1108    if (coeffraf == 1) then
1109        locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
1110        y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
1111        return
1112    end if
1113
1114    diffmod = 0
1115    if (mod(coeffraf,2) == 0) diffmod = 1
1116!
1117    allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
1118    ypos = s_child
1119!
1120    locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
1121    locind_parent_last = 1 + agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
1122!
1123    xpmin = s_parent + (locind_parent_left-1)*ds_parent
1124    i1 = 1+agrif_int((xpmin-s_child)/ds_child)
1125!
1126    do i = 1,np+1
1127        xhalf(i) = i - 0.5
1128    enddo
1129!
1130! Compute divided differences
1131!
1132    dd(1,1:np)   = x(1:np)
1133    dd(2,1:np-1) = 0.5*( dd(1,2:np) - dd(1,1:np-1) )
1134    dd(3,1:np-2) = (1./3.)*( dd(2,2:np-1) - dd(2,1:np-2) )
1135!
1136    allocate( xbar( coeffraf,2 ) )
1137    xi = 0.5
1138    do i = 1,coeffraf
1139        xbar(i,1) = (i-1)*ds_child/ds_parent - xi
1140        xbar(i,2) =  i   *ds_child/ds_parent - xi
1141    enddo
1142!
1143    ipos = i1
1144!
1145    do i = locind_parent_left,locind_parent_last
1146        left = i
1147        do jj = 2,3
1148            if(abs(dd(jj,left)) > abs(dd(jj,left-1))) left = left-1
1149        enddo
1150!
1151!       convert to Taylor series form
1152        call taylor(i,xhalf(left:left+2),dd(1:3,left),c(1:3,i))
1153    enddo
1154!
1155! Evaluate the reconstruction on each cell
1156!
1157    do i = locind_parent_left,locind_parent_last
1158!
1159        pos = 1
1160!
1161        do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
1162            ytemp(jj) = ( c(1,i)*(xbar(pos,2)-xbar(pos,1))              &
1163                        + c(2,i)*(xbar(pos,2)*xbar(pos,2) -             &
1164                                  xbar(pos,1)*xbar(pos,1))              &
1165                        + c(3,i)*(xbar(pos,2)*xbar(pos,2)*xbar(pos,2) - &
1166                                  xbar(pos,1)*xbar(pos,1)*xbar(pos,1))  &
1167                        ) * coeffraf
1168            pos = pos+1
1169        enddo
1170        ipos = ipos + coeffraf
1171!
1172    enddo
1173!
1174    y(1:nc) = ytemp(1:nc)
1175    deallocate(ytemp,xbar)
1176!---------------------------------------------------------------------------------------------------
1177end subroutine ENO1d
1178!===================================================================================================
1179!
1180!     ************************************************************************** 
1181!CC   Subroutine ppm1d_lim
1182!     **************************************************************************
1183!
1184      Subroutine ppm1d_lim(x,y,np,nc,s_parent,s_child,ds_parent,ds_child) 
1185!
1186!CC   Description:
1187!CC   Subroutine to do a 1D interpolation and apply monotonicity constraints
1188!CC   using piecewise parabolic method 
1189!CC   on a child grid (vector y) from its parent grid (vector x).
1190!C    Method:
1191!
1192!     Declarations:
1193!
1194      Implicit none
1195!       
1196!     Arguments
1197      Integer             :: np,nc     
1198      Real, Dimension(np) :: x     
1199      Real, Dimension(nc) :: y
1200      Real, Dimension(:),Allocatable :: ytemp
1201      Real(kind=8)        :: s_parent,s_child,ds_parent,ds_child
1202!
1203!     Local scalars
1204      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
1205      Integer :: iparent,ipos,pos,nmin,nmax
1206      Real(kind=8)    :: ypos
1207      integer :: i1,jj
1208      Real(kind=8) :: xpmin
1209      real :: cavg,a,b
1210!     
1211      Real :: xrmin,xrmax,am3,s2,s1 
1212      Real, Dimension(np) :: dela,xr,xl,delta,a6,slope,slope2
1213      Real, Dimension(:),Allocatable  :: diff,diff2,diff3   
1214      INTEGER :: diffmod
1215!     
1216      coeffraf = nint(ds_parent/ds_child)
1217!
1218      If (coeffraf == 1) Then
1219          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
1220          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
1221          return
1222      End If
1223!     
1224      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
1225      ypos = s_child 
1226!
1227      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
1228      locind_parent_last = 1 + &
1229           agrif_ceiling((ypos +(nc - 1)  &
1230           *ds_child - s_parent)/ds_parent) 
1231!
1232      xpmin = s_parent + (locind_parent_left-1)*ds_parent       
1233      i1 = 1+agrif_int((xpmin-s_child)/ds_child)       
1234!     
1235      Allocate( diff(coeffraf),diff2(coeffraf),diff3(coeffraf) )
1236!     
1237         diff(:) = ds_child/ds_parent
1238!     
1239      Do i=1,coeffraf
1240         a = real(i-1)*ds_child/ds_parent
1241         b = real(i)*ds_child/ds_parent
1242         diff2(i) = 0.5*(b*b - a*a) 
1243         diff3(i) = (1./3.)*(b*b*b - a*a*a)
1244      End do
1245!
1246      if( locind_parent_last+2 <= np ) then
1247           nmax = locind_parent_last+2   
1248      else if( locind_parent_last+1 <= np ) then
1249           nmax = locind_parent_last+1
1250      else
1251           nmax = locind_parent_last 
1252      endif     
1253!     
1254      if(locind_parent_left-1 >= 1) then
1255          nmin = locind_parent_left-1
1256      else
1257          nmin = locind_parent_left
1258      endif   
1259!
1260      Do i = nmin,nmax
1261         slope(i) = x(i) - x(i-1)
1262         slope2(i) = 2.*abs(slope(i))
1263      Enddo
1264!
1265      Do i = nmin,nmax-1
1266         dela(i) = 0.5 * ( slope(i) + slope(i+1) )
1267! Van Leer slope limiter
1268         dela(i) = min( abs(dela(i)),slope2(i), &
1269                       slope2(i+1) )*sign(1.,dela(i))
1270         IF( slope(i)*slope(i+1) <= 0. ) dela(i) = 0.
1271      Enddo
1272!
1273      Do i = nmin,nmax-2
1274         xr(i) = x(i) + (1./2.)*slope(i+1) + (-1./6.)*dela(i+1) &
1275                                          + ( 1./6. )*dela(i)
1276      Enddo
1277!
1278      Do i = nmin,nmax-2
1279         xrmin = min(x(i),x(i+1))
1280         xrmax = max(x(i),x(i+1))
1281         xr(i) = min(xr(i),xrmax)
1282         xr(i) = max(xr(i),xrmin)
1283         xl(i+1) = xr(i)         
1284      Enddo
1285! apply parabolic monotonicity
1286       Do i = locind_parent_left,locind_parent_last
1287          If( ( (xr(i)-x(i))* (x(i)-xl(i)) ) .le. 0. ) then
1288             xl(i) = x(i) 
1289             xr(i) = x(i)
1290          Endif         
1291          delta(i) = xr(i) - xl(i)
1292          am3 = 3. * x(i)
1293          s1  = am3 - 2. * xr(i)
1294          s2  = am3 - 2. * xl(i)
1295          IF( delta(i) * (xl(i) - s1) .le. 0. ) xl(i) = s1
1296          IF( delta(i) * (s2 - xr(i)) .le. 0. ) xr(i) = s2
1297          delta(i) = xr(i) - xl(i)
1298          a6(i) = 6.*x(i)-3.*(xl(i) +xr(i))
1299!
1300       End do   
1301!
1302        diffmod = 0
1303       IF (mod(coeffraf,2) == 0) diffmod = 1           
1304!
1305        ipos = i1
1306!               
1307        Do iparent = locind_parent_left,locind_parent_last       
1308             pos=1
1309             cavg = 0.
1310             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
1311!
1312               ytemp(jj) = (diff(pos)*xl(iparent)   &
1313                  + diff2(pos) &
1314                  *  (delta(iparent)+a6(iparent)) &
1315                  - diff3(pos)*a6(iparent))*coeffraf
1316                             
1317               cavg = cavg + ytemp(jj)
1318               pos = pos+1 
1319             End do
1320             ipos = ipos + coeffraf
1321!
1322        End do     
1323!
1324!
1325        y(1:nc)=ytemp(1:nc)                                 
1326        deallocate(ytemp)               
1327        deallocate(diff, diff2, diff3)
1328      Return
1329      End Subroutine ppm1d_lim
1330!===================================================================================================
1331!  subroutine taylor
1332!---------------------------------------------------------------------------------------------------
1333subroutine taylor ( ind, xhalf, dd, c )
1334!---------------------------------------------------------------------------------------------------
1335    integer,            intent(in)  :: ind
1336    real, dimension(3), intent(in)  :: xhalf
1337    real, dimension(3), intent(in)  :: dd
1338    real, dimension(3), intent(out) :: c
1339!
1340    real, dimension(0:3,0:3) :: d
1341    integer                  :: i, j
1342!
1343    d(0,0:3) = 1.0
1344    do i = 1,3
1345        d(i,0) = (ind-xhalf(i))*d(i-1,0)
1346    enddo
1347!
1348    do i = 1,3
1349        do j = 1,3-i
1350            d(i,j) = d(i,j-1) + (ind-xhalf(i+j))*d(i-1,j)
1351        enddo
1352    enddo
1353!
1354    do j = 1,3
1355        c(j) = 0.
1356        do i=0,3-j
1357            c(j) = c(j) + d(i,j)*dd(i+j)
1358        enddo
1359    enddo
1360!---------------------------------------------------------------------------------------------------
1361end subroutine taylor
1362!===================================================================================================
1363!
1364!===================================================================================================
1365!  function Agrif_limiter_vanleer
1366!---------------------------------------------------------------------------------------------------
1367real function Agrif_limiter_vanleer ( tab ) result(res)
1368!---------------------------------------------------------------------------------------------------
1369    real, dimension(3), intent(in) :: tab
1370!
1371    real :: p1, p2, p3
1372
1373    p1 =    (tab(3)-tab(1))/2.
1374    p2 = 2.*(tab(2)-tab(1))
1375    p3 = 2.*(tab(3)-tab(2))
1376
1377    if ((p1>0.).AND.(p2>0.).AND.(p3>0)) then
1378        res = minval((/p1,p2,p3/))
1379    elseif ((p1<0.).AND.(p2<0.).AND.(p3<0)) then
1380        res = maxval((/p1,p2,p3/))
1381    else
1382        res=0.
1383    endif
1384!---------------------------------------------------------------------------------------------------
1385end function Agrif_limiter_vanleer
1386!===================================================================================================
1387!
1388end module Agrif_InterpBasic
Note: See TracBrowser for help on using the repository browser.