source: branches/UKMO/r6232_tracer_advection/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinterpbasic.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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