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 @ 13924

Last change on this file since 13924 was 13924, checked in by jchanut, 3 years ago

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