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 branches/UKMO/dev_r5518_GO6_package/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinterpbasic.F90 @ 7993

Last change on this file since 7993 was 7993, checked in by frrh, 7 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.
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.