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.F in trunk/AGRIF/AGRIF_FILES – NEMO

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

Last change on this file since 779 was 779, checked in by rblod, 16 years ago

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