source: trunk/AGRIF/AGRIF_FILES/modinterpbasic.F @ 396

Last change on this file since 396 was 396, checked in by opalod, 15 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.9 KB
Line 
1!
2! $Id$
3!
4C     AGRIF (Adaptive Grid Refinement In Fortran)
5C
6C     Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
7C                        Christophe Vouland (Christophe.Vouland@imag.fr)   
8C
9C     This program is free software; you can redistribute it and/or modify
10C     it under the terms of the GNU General Public License as published by
11C     the Free Software Foundation; either version 2 of the License, or
12C     (at your option) any later version.
13C
14C     This program is distributed in the hope that it will be useful,
15C     but WITHOUT ANY WARRANTY; without even the implied warranty of
16C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17C     GNU General Public License for more details.
18C
19C     You should have received a copy of the GNU General Public License
20C     along with this program; if not, write to the Free Software
21C     Foundation, Inc., 59 Temple Place- Suite 330, Boston, MA 02111-1307, USA.
22C
23C
24C
25CCC   Module Agrif_Interpbasic
26C
27      Module Agrif_Interpbasic
28C
29CCC   Description:
30CCC   Module containing different procedures of interpolation (linear,lagrange,
31CCC   spline,...) used in the Agrif_Interpolation module.
32C
33C     Modules used:
34      USE Agrif_types
35C
36      IMPLICIT NONE
37C             
38      CONTAINS
39C     Define procedures contained in this module
40C 
41C
42C     ************************************************************************** 
43CCC   Subroutine Linear1d 
44C     ************************************************************************** 
45C 
46      Subroutine Linear1d(x,y,np,nc,
47     &                    s_parent,s_child,ds_parent,ds_child) 
48C
49CCC   Description:
50CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
51CCC   its parent grid (vector x). 
52C
53CC    Method:
54C
55C     Declarations:
56C
57     
58C       
59C     Arguments
60      INTEGER             :: np,nc     
61      REAL, DIMENSION(np) :: x     
62      REAL, DIMENSION(nc) :: y 
63      REAL                :: s_parent,s_child,ds_parent,ds_child
64C
65C     Local scalars
66      INTEGER :: i,coeffraf,locind_parent_left
67      REAL    :: ypos,globind_parent_left
68C
69C
70
71      coeffraf = nint(ds_parent/ds_child)
72C
73      if (coeffraf == 1) then
74C
75          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
76C       
77          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
78C
79          return
80C
81      endif                         
82C
83      ypos = s_child     
84C
85      do i = 1,nc-1
86C
87        locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
88C
89        globind_parent_left = s_parent 
90     &                        + (locind_parent_left - 1)*ds_parent
91C       
92        y(i) = ((globind_parent_left + ds_parent - ypos)
93     &          *x(locind_parent_left)
94     &        + (ypos - globind_parent_left)
95     &          *x(locind_parent_left+1))
96     &         / ds_parent               
97C
98        ypos = ypos + ds_child
99C
100      enddo
101C
102      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
103C
104      if (locind_parent_left == np) then
105C
106          y(nc) = x(np)
107C
108        else
109C
110          globind_parent_left = s_parent 
111     &                        + (locind_parent_left - 1)*ds_parent 
112C     
113          y(nc) = ((globind_parent_left + ds_parent - ypos)
114     &            *x(locind_parent_left)
115     &          + (ypos - globind_parent_left)
116     &            *x(locind_parent_left+1))
117     &           / ds_parent     
118C
119      endif                                         
120C           
121      Return
122C
123C       
124      End Subroutine Linear1d   
125C
126C
127C
128C     ************************************************************************** 
129CCC   Subroutine Lagrange1d 
130C     **************************************************************************
131C
132      Subroutine Lagrange1d(x,y,np,nc,
133     &                      s_parent,s_child,ds_parent,ds_child)
134C
135CCC   Description:
136CCC   Subroutine to do a lagrange 1D interpolation on a child grid (vector y) 
137CCC   from its parent grid (vector x). 
138C
139CC    Method:
140C
141C     Declarations:
142C
143     
144C             
145C     Arguments
146      INTEGER             :: np,nc     
147      REAL, DIMENSION(np) :: x     
148      REAL, DIMENSION(nc) :: y 
149      REAL                :: s_parent,s_child,ds_parent,ds_child 
150C
151C     Local scalars
152      INTEGER :: i,coeffraf,locind_parent_left
153      REAL    :: ypos,globind_parent_left
154      REAL    :: X1,X2,X3 
155C
156C 
157      if (np <= 2) then
158C     
159          Call Linear1D(x,y,np,nc,
160     &                  s_parent,s_child,ds_parent,ds_child)
161C         
162         Return
163C 
164      endif
165C
166      coeffraf = nint(ds_parent/ds_child)
167C
168      if (coeffraf == 1) then
169C
170          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
171C       
172          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
173C
174          return
175C
176      endif
177C
178      ypos = s_child     
179C
180      do i = 1,nc
181C
182        locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
183C
184        globind_parent_left = s_parent 
185     &                        + (locind_parent_left - 1)*ds_parent 
186C
187        if (locind_parent_left+2 <= np) then           
188C
189            X1 = (x(locind_parent_left+1)-x(locind_parent_left))
190     &           /ds_parent
191C 
192            X2 = (x(locind_parent_left+2)-x(locind_parent_left+1))
193     &           /ds_parent
194C
195            X3 = (X2 - X1)/(2.*ds_parent)                 
196C
197            y(i) = x(locind_parent_left) + 
198     &             (ypos - globind_parent_left)*X1 +
199     &             (ypos - globind_parent_left)*
200     &             (ypos - globind_parent_left - ds_parent)*X3 
201C
202          elseif (locind_parent_left+1 <= np) then
203C
204            X1 = (x(locind_parent_left)-x(locind_parent_left-1))
205     &           /ds_parent
206C     
207            X2 = (x(locind_parent_left+1)-x(locind_parent_left))
208     &           /ds_parent
209C 
210            X3 = (X2 - X1)/(2.*ds_parent)                 
211C
212            y(i) = x(locind_parent_left-1) + 
213     &             (ypos - globind_parent_left - ds_parent)*X1 +
214     &             (ypos - globind_parent_left - ds_parent)*
215     &             (ypos - globind_parent_left)*X3
216C
217          else
218C
219            X1 = (x(locind_parent_left-1)-x(locind_parent_left-2))
220     &           /ds_parent
221C     
222            X2 = (x(locind_parent_left)-x(locind_parent_left-1))
223     &           /ds_parent
224C 
225            X3 = (X2 - X1)/(2.*ds_parent)                 
226C
227            y(i) = x(locind_parent_left-2) + 
228     &             (ypos - globind_parent_left - 2.*ds_parent)*X1 +
229     &             (ypos - globind_parent_left - 2.*ds_parent)*
230     &             (ypos - globind_parent_left - ds_parent)*X3
231C
232        endif   
233C       
234        ypos = ypos + ds_child 
235C
236      enddo
237C
238      return
239C
240C       
241      End Subroutine Lagrange1d 
242C
243C
244C     ************************************************************************** 
245CCC   Subroutine Constant1d 
246C     ************************************************************************** 
247C 
248      Subroutine constant1d(x,y,np,nc,
249     &                    s_parent,s_child,ds_parent,ds_child) 
250C
251CCC   Description:
252CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
253CCC   its parent grid (vector x). 
254C
255CC    Method:
256C
257C     Declarations:
258C
259     
260C       
261C     Arguments
262      INTEGER             :: np,nc     
263      REAL, DIMENSION(np) :: x     
264      REAL, DIMENSION(nc) :: y 
265      REAL                :: s_parent,s_child,ds_parent,ds_child
266C
267C     Local scalars
268      INTEGER :: i,coeffraf,locind_parent
269      REAL    :: ypos
270C
271C
272
273      coeffraf = nint(ds_parent/ds_child)
274C
275      if (coeffraf == 1) then
276C
277          locind_parent = 1 + nint((s_child - s_parent)/ds_parent)
278C       
279          y(1:nc) = x(locind_parent:locind_parent+nc-1)
280C
281          return
282C
283      endif                         
284C
285      ypos = s_child     
286C
287      do i = 1,nc
288C
289        locind_parent = 1 + nint((ypos - s_parent)/ds_parent)
290C
291        y(i) = x(locind_parent)
292C
293        ypos = ypos + ds_child
294C
295      enddo
296C           
297      Return
298C
299C       
300      End Subroutine constant1d   
301C
302C     ************************************************************************** 
303CCC   Subroutine Linear1dconserv
304C     ************************************************************************** 
305C 
306      Subroutine Linear1dconserv(x,y,np,nc,
307     &                    s_parent,s_child,ds_parent,ds_child) 
308C
309CCC   Description:
310CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
311CCC   its parent grid (vector x). 
312C
313CC    Method:
314C
315C     Declarations:
316C
317      Implicit none
318C       
319C     Arguments
320      Integer             :: np,nc     
321      Real, Dimension(np) :: x     
322      Real, Dimension(nc) :: y
323      Real, Dimension(:),Allocatable :: ytemp
324      Real                :: s_parent,s_child,ds_parent,ds_child
325C
326C     Local scalars
327      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
328      Real    :: ypos
329      integer :: i1,i2,ii
330      real :: xpmin,xpmax,slope
331      INTEGER :: diffmod
332      REAL :: xdiffmod
333
334C
335C
336
337      coeffraf = nint(ds_parent/ds_child)
338C
339      If (coeffraf == 1) Then
340C
341          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
342C       
343          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
344C
345          return
346C
347      End If
348C             
349      diffmod = 0
350      IF (mod(coeffraf,2) == 0) diffmod = 1 
351
352      xdiffmod = real(diffmod)/2.
353                         
354      allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
355C
356      ypos = s_child 
357C     
358      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
359
360      locind_parent_last = 1 +
361     &  agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
362     
363      xpmin = s_parent + (locind_parent_left-1)*ds_parent 
364      xpmax = s_parent + (locind_parent_last-1)*ds_parent         
365     
366      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
367      i2 = 1+agrif_int((xpmax-s_child)/ds_child)
368
369      i = i1
370     
371      if (locind_parent_left == 1) then
372        slope=
373     &   (x(locind_parent_left+1)-x(locind_parent_left))/(coeffraf)
374      else
375         slope=
376     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
377      endif
378     
379        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
380          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
381        enddo 
382
383        locind_parent_left = locind_parent_left + 1
384                     
385      do i=i1 +  coeffraf, i2 - coeffraf,coeffraf
386        slope=
387     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
388        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
389          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
390        enddo
391        locind_parent_left = locind_parent_left + 1
392      enddo
393     
394      i = i2
395     
396      if (locind_parent_left == np) then
397        slope=
398     &   (x(locind_parent_left)-x(locind_parent_left-1))/(coeffraf)
399      else
400         slope=
401     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
402      endif
403     
404        do ii=i-coeffraf/2+diffmod,nc
405          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
406        enddo     
407C
408      y(1:nc)=ytemp(1:nc)                                   
409C           
410      deallocate(ytemp)
411      Return
412C       
413      End Subroutine Linear1dconserv
414     
415C
416C     ************************************************************************** 
417CCC   Subroutine Linear1dconservlim
418C     ************************************************************************** 
419C 
420      Subroutine Linear1dconservlim(x,y,np,nc,
421     &                    s_parent,s_child,ds_parent,ds_child) 
422C
423CCC   Description:
424CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
425CCC   its parent grid (vector x). 
426C
427CC    Method:
428C
429C     Declarations:
430C
431      Implicit none
432C       
433C     Arguments
434      Integer             :: np,nc     
435      Real, Dimension(np) :: x     
436      Real, Dimension(nc) :: y
437      Real, Dimension(:),Allocatable :: ytemp
438      Real                :: s_parent,s_child,ds_parent,ds_child
439C
440C     Local scalars
441      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
442      Real    :: ypos
443      integer :: i1,i2,ii
444      real :: xpmin,xpmax,slope
445      INTEGER :: diffmod
446      real :: xdiffmod
447C
448C
449
450      coeffraf = nint(ds_parent/ds_child)
451C
452      If (coeffraf == 1) Then
453C
454          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
455C       
456          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
457C
458          return
459C
460      End If
461C
462      IF (coeffraf .NE.3) THEN
463      print *,'LINEARCONSERVLIM not ready for refinement ratio = ',
464     &   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))
474C
475      ypos = s_child 
476C     
477      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
478
479      locind_parent_last = 1 +
480     &  agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
481     
482      xpmin = s_parent + (locind_parent_left-1)*ds_parent 
483      xpmax = s_parent + (locind_parent_last-1)*ds_parent         
484     
485      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
486      i2 = 1+agrif_int((xpmax-s_child)/ds_child)
487
488      i = i1
489     
490      if (locind_parent_left == 1) then
491        slope=0.       
492      else
493        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
494        slope = slope / coeffraf   
495      endif
496     
497        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
498          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
499        enddo 
500
501        locind_parent_left = locind_parent_left + 1
502                     
503      do i=i1 +  coeffraf, i2 - coeffraf,coeffraf     
504        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
505        slope = slope / coeffraf
506
507        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
508          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
509        enddo
510        locind_parent_left = locind_parent_left + 1
511      enddo
512     
513      i = i2
514     
515      if (locind_parent_left == np) then
516        slope=0.     
517      else
518        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
519        slope = slope / coeffraf     
520      endif
521     
522        do ii=i-coeffraf/2+diffmod,nc
523          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
524        enddo     
525C
526      y(1:nc)=ytemp(1:nc)                                   
527C           
528      deallocate(ytemp)
529      Return
530C       
531      End Subroutine Linear1dconservlim     
532C         
533
534C     ************************************************************************** 
535CCC   Subroutine ppm1d
536C     ************************************************************************** 
537C 
538      Subroutine ppm1d(x,y,np,nc,
539     &                    s_parent,s_child,ds_parent,ds_child) 
540C
541CCC   Description:
542CCC   Subroutine to do a 1D interpolation and apply monotonicity constraints
543CCC   using piecewise parabolic method 
544CCC   on a child grid (vector y) from its parent grid (vector x).
545CC    Method:
546C
547C     Declarations:
548C
549      Implicit none
550C       
551C     Arguments
552      Integer             :: np,nc     
553      Real, Dimension(np) :: x     
554      Real, Dimension(nc) :: y
555      Real, Dimension(:),Allocatable :: ytemp
556      Real                :: s_parent,s_child,ds_parent,ds_child
557C
558C     Local scalars
559      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
560      Integer :: iparent,ipos,pos,nmin,nmax
561      Real    :: ypos
562      integer :: i1,jj
563      Real :: xpmin,cavg,a,b
564C     
565      Real :: xrmin,xrmax,am3,s2,s1 
566      Real, Dimension(np) :: dela,xr,xl,delta,a6,slope,slope2
567      Real, Dimension(:),Allocatable  :: diff,diff2,diff3   
568      INTEGER :: diffmod
569C     
570      coeffraf = nint(ds_parent/ds_child)
571C
572      If (coeffraf == 1) Then
573          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
574          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
575          return
576      End If
577C     
578      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
579      ypos = s_child 
580C
581      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
582      locind_parent_last = 1 +
583     &      agrif_ceiling((ypos +(nc - 1) 
584     &      *ds_child - s_parent)/ds_parent) 
585C
586      xpmin = s_parent + (locind_parent_left-1)*ds_parent       
587      i1 = 1+agrif_int((xpmin-s_child)/ds_child)       
588C     
589      Allocate( diff(coeffraf),diff2(coeffraf),diff3(coeffraf) )
590C     
591         diff(:) = ds_child/ds_parent
592C     
593      Do i=1,coeffraf
594         a = real(i-1)*ds_child/ds_parent
595         b = real(i)*ds_child/ds_parent
596         diff2(i) = 0.5*(b*b - a*a) 
597         diff3(i) = (1./3.)*(b*b*b - a*a*a)
598      End do
599C
600      if( locind_parent_last+2 <= np ) then
601           nmax = locind_parent_last+2   
602      else if( locind_parent_last+1 <= np ) then
603           nmax = locind_parent_last+1
604      else
605           nmax = locind_parent_last 
606      endif     
607C     
608      if(locind_parent_left-1 >= 1) then
609          nmin = locind_parent_left-1
610      else
611          nmin = locind_parent_left
612      endif   
613C 
614      Do i = nmin,nmax
615         slope(i) = x(i) - x(i-1)
616         slope2(i) = 2.*abs(slope(i))
617      Enddo
618C
619      Do i = nmin,nmax-1
620         dela(i) = 0.5 * ( slope(i) + slope(i+1) )
621C Van Leer slope limiter
622         dela(i) = min( abs(dela(i)),slope2(i),
623     &                  slope2(i+1) )*sign(1.,dela(i))
624         IF( slope(i)*slope(i+1) <= 0. ) dela(i) = 0.
625      Enddo
626C
627      Do i = nmin,nmax-2
628         xr(i) = x(i) + (1./2.)*slope(i+1) + (-1./6.)*dela(i+1)
629     &                                     + ( 1./6. )*dela(i)
630      Enddo
631C
632      Do i = nmin,nmax-2
633         xrmin = min(x(i),x(i+1))
634         xrmax = max(x(i),x(i+1))
635         xr(i) = min(xr(i),xrmax)
636         xr(i) = max(xr(i),xrmin)
637         xl(i+1) = xr(i)         
638      Enddo
639C apply parabolic monotonicity
640       Do i = locind_parent_left,locind_parent_last
641          If( ( (xr(i)-x(i))* (x(i)-xl(i)) ) .le. 0. ) then
642             xl(i) = x(i) 
643             xr(i) = x(i)
644          Endif         
645          delta(i) = xr(i) - xl(i)
646          am3 = 3. * x(i)
647          s1  = am3 - 2. * xr(i)
648          s2  = am3 - 2. * xl(i)
649          IF( delta(i) * (xl(i) - s1) .le. 0. ) xl(i) = s1
650          IF( delta(i) * (s2 - xr(i)) .le. 0. ) xr(i) = s2
651          delta(i) = xr(i) - xl(i)
652          a6(i) = 6.*x(i)-3.*(xl(i) +xr(i))
653C
654       End do   
655C
656        diffmod = 0
657   IF (mod(coeffraf,2) == 0) diffmod = 1           
658C
659        ipos = i1
660C               
661        Do iparent = locind_parent_left,locind_parent_last       
662             pos=1
663             cavg = 0.
664             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
665C
666               ytemp(jj) = (diff(pos)*xl(iparent)   
667     &             + diff2(pos)
668     &             *  (delta(iparent)+a6(iparent))
669     &             - diff3(pos)*a6(iparent))*coeffraf
670                             
671               cavg = cavg + ytemp(jj)
672               pos = pos+1 
673             End do
674             ipos = ipos + coeffraf
675C
676        End do     
677C
678C
679        y(1:nc)=ytemp(1:nc)                                 
680        deallocate(ytemp)               
681        deallocate(diff, diff2, diff3)
682      Return
683      End Subroutine ppm1d
684C                             
685C     ************************************************************************** 
686CCC   Subroutine eno1d
687C     ************************************************************************** 
688C 
689      Subroutine eno1d(x,y,np,nc,
690     &                    s_parent,s_child,ds_parent,ds_child) 
691C
692CCC   Description:
693CCC   ---- p 163-164 Computational gasdynamics ----
694CCC   Subroutine to do a 1D interpolation 
695CCC   using piecewise polynomial ENO reconstruction technique 
696CCC   on a child grid (vector y) from its parent grid (vector x).
697CC    Method:
698C
699C     Declarations:
700C
701      Implicit none
702C       
703C     Arguments
704      Integer             :: np,nc     
705      Real, Dimension(np) :: x     
706      Real, Dimension(nc) :: y
707      Real, Dimension(:),Allocatable :: ytemp
708      Real                :: s_parent,s_child,ds_parent,ds_child
709C
710C     Local scalars
711      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
712      Integer :: ipos,pos
713      Real    :: ypos,xi
714      integer :: i1,jj
715      Real :: xpmin,cavg
716C 
717      Real, Dimension(3,np) :: dd,c
718      Integer :: left
719C           
720      Real, DImension(1:np+1) :: xhalf
721      Real, Dimension(:,:),Allocatable  :: Xbar 
722      INTEGER :: diffmod     
723C
724      coeffraf = nint(ds_parent/ds_child)
725C     
726      If (coeffraf == 1) Then
727          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
728          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
729          return
730      End If     
731     
732      diffmod = 0
733      IF (mod(coeffraf,2) == 0) diffmod = 1 
734C     
735      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
736      ypos = s_child 
737      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
738      locind_parent_last = 1 +
739     &      agrif_ceiling((ypos +(nc - 1) *ds_child - 
740     &      s_parent)/ds_parent)       
741      xpmin = s_parent + (locind_parent_left-1)*ds_parent     
742      i1 = 1+agrif_int((xpmin-s_child)/ds_child)           
743C     
744      xhalf(np+1) = np + 0.5
745      Do i = 1,np
746          xhalf(i) = i - 0.5
747      Enddo
748C
749C compute divided differences
750C
751      dd(1,1:np) = x(1:np)
752      dd(2,1:np-1) = 0.5*( dd(1,2:np) - dd(1,1:np-1) )
753      dd(3,1:np-2) = (1./3.)*( dd(2,2:np-1) - dd(2,1:np-2) )
754C
755      Allocate( Xbar( coeffraf,2 ) )
756      xi = 0.5
757      Do i = 1,coeffraf
758        Xbar(i,1) = (i-1)*ds_child/ds_parent - xi
759        Xbar(i,2) = i*ds_child/ds_parent - xi
760      Enddo
761C
762      ipos = i1
763C           
764      DO i = locind_parent_left,locind_parent_last           
765         left = i           
766         do jj = 2,3
767             If(abs(dd(jj,left)) .gt. abs(dd(jj,left-1)))
768     &               left = left-1           
769         enddo
770C           
771C  convert to Taylor series form
772C
773         Call Taylor(i,xhalf(left:left+2),dd(1:3,left),c(1:3,i))
774      ENDDO     
775C
776C evaluate the reconstruction on each cell
777C
778       DO i = locind_parent_left,locind_parent_last 
779C
780         cavg = 0.
781         pos = 1.           
782C         
783         Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
784           ytemp(jj) =(c(1,i)*(Xbar(pos,2)-Xbar(pos,1))
785     &                +c(2,i)*(Xbar(pos,2)*Xbar(pos,2)-
786     &                         Xbar(pos,1)*Xbar(pos,1))
787     &                +c(3,i)*(Xbar(pos,2)*Xbar(pos,2)*Xbar(pos,2)-
788     &                         Xbar(pos,1)*Xbar(pos,1)*Xbar(pos,1)))
789     &                         *coeffraf
790           cavg = cavg + ytemp(jj)
791           pos = pos+1
792         Enddo               
793         ipos = ipos + coeffraf                 
794      ENDDO
795C
796      y(1:nc)=ytemp(1:nc)                                           
797      deallocate(ytemp,Xbar)                 
798C     
799      Return       
800      End Subroutine eno1d
801C 
802C     
803C     ************************************************************************** 
804CCC   Subroutine taylor
805C     ************************************************************************** 
806C       
807      subroutine taylor(ind,xhalf,dd,c)     
808C     
809      Integer :: ind
810      real,dimension(3) :: dd,c     
811      real,dimension(0:3,0:3) :: d 
812      real,dimension(3) :: xhalf   
813      integer ::i,j
814C     
815C
816      d(0,0:3)=1.
817      do i = 1,3
818         d(i,0)=(ind-xhalf(i))*d(i-1,0)
819      enddo 
820C     
821      do i = 1,3
822         do j = 1,3-i
823           d(i,j) = d(i,j-1) + (ind-xhalf(i+j))*d(i-1,j)
824         enddo
825      enddo
826C         
827      do j = 1,3       
828         c(j) = 0.
829         do i=0,3-j
830            c(j) = c(j) + d(i,j)*dd(i+j)         
831         enddo
832      enddo
833C     
834      end subroutine taylor   
835     
836     
837      REAL FUNCTION vanleer(tab)
838      REAL, DIMENSION(3) :: tab
839       real res1
840       real p1,p2,p3
841       
842       p1=(tab(3)-tab(1))/2.
843       p2=2.*(tab(2)-tab(1))
844       p3=2.*(tab(3)-tab(2))
845       
846       if ((p1>0.).AND.(p2>0.).AND.(p3>0)) then
847          res1=minval((/p1,p2,p3/))
848       elseif ((p1<0.).AND.(p2<0.).AND.(p3<0)) then
849          res1=maxval((/p1,p2,p3/))
850       else
851          res1=0.
852       endif
853         
854          vanleer = res1   
855     
856     
857      END FUNCTION vanleer 
858
859C     
860      End Module Agrif_Interpbasic
Note: See TracBrowser for help on using the repository browser.