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/dev_r12970_AGRIF_CMEMS/AGRIF_FILES – NEMO

source: vendors/AGRIF/dev_r12970_AGRIF_CMEMS/AGRIF_FILES/modinterpbasic.F90 @ 13144

Last change on this file since 13144 was 13027, checked in by rblod, 4 years ago

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