source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinterpbasic.F @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 5 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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