PROGRAM MAP_3
c
c Givn two uniform grids defined by the locations of the lefthand
c zone interfaces XL and XNUL. The grids are offset by an amount
c DX, that is XNUL(I) = XL(I) + DX.
c This program creates piecewise parabolic interpolations
c from zone averaged data, defined on the XL grid.
c Monotonicity constraints are applied. For those zones
c considered to be in a discontinuity, an interpolation more
c appropriate for discontinuous functions is used, producing
c a better representation of unresolved shocks.
c Regions of over lap between the two grid are determined and
c advective fluxes found by integrating the PPM interpolations
c over these regions. The zone averages on the grid XNUL are found
c by differeing these fluxes.
c
c Note, the interpolations require 3 boundary zones at each end, and
c the integrations will require an additional zone.
c
PARAMETER ( NBDY = 4 )
ALLOCATABLE :: XL(:), A_AVG(:), AL(:), AR(:), DA(:), A6(:),
^ SMALLDA(:), UNSMOOTH(:), DISCONT(:), DX(:), DA_AVL(:),
^ DVOLL(:), XF(:), XF1(:), XLNU(:), A_AVGNU(:)
c
OPEN (unit = 11, file ='INPUT.AVG')
c The data consists of N+2 lines. The first line
c contains the integer N, the number to zone averages in the data set.
c This is followed by N lines, each line having two numbers, a value for
c the location of the left edge of the zone XL, and the zone average for
c that zone A_AVG. The lines are assumed to be in sequential order according
c to XL. An additional line has only one value, the final left zone edge. A
c grid of N zones requires N+1 left zone edge locations. An appropriate data
c set can be created from the program CREATE_ZONE_AVERAGES. The XL values
c are used for plotting, not in the interpolation process.
c Read in the number of zones in the file
READ (11,*) N
c
c Next, allocate the appropriate space
ALLOCATE (XL(1-NBDY:N+1+NBDY))
ALLOCATE (AL(1-NBDY:N+1+NBDY))
ALLOCATE (AR(1-NBDY:N+1+NBDY))
ALLOCATE (DA(1-NBDY:N+NBDY))
ALLOCATE (A6(1-NBDY:N+NBDY))
ALLOCATE (A_AVG(1-NBDY:N+NBDY))
ALLOCATE (DA_AVL(1-NBDY:N+NBDY))
ALLOCATE (XLNU(1-NBDY:N+1+NBDY))
ALLOCATE (A_AVGNU(1-NBDY:N+NBDY))
ALLOCATE (SMALLDA(1-NBDY:N+NBDY))
ALLOCATE (UNSMOOTH(1-NBDY:N+NBDY))
ALLOCATE (DISCONT(1-NBDY:N+NBDY))
ALLOCATE (DVOLL(1-NBDY:N+NBDY))
ALLOCATE (XF(1-NBDY:N+NBDY))
ALLOCATE (XF1(1-NBDY:N+NBDY))
c The following arrays will hold various geometrical factors
c same grid:
ALLOCATE (DX(1-NBDY:N+NBDY))
DO I=1,N
READ (11,*) XL(I),A_AVG(I)
ENDDO
READ (11,*) XL(N+1)
c
c
c Query the user for desired boundary conditions
c
PRINT *,' Please enter 0 for reflecting boundary conditions'
PRINT *,' or 1 for periodic boundary conditions'
READ *, IBDY_TYPE
IF ( IBDY_TYPE .EQ. 0 ) THEN
PRINT *,' Selecting reflecting boundary conditions'
ELSE
IF ( IBDY_TYPE .EQ. 1 ) THEN
PRINT *,' Selecting periodic boundary conditions'
ELSE
PRINT *,' Unknown boundary type:',IBDY_TYPE
call abort()
ENDIF
ENDIF
c
c Initialize the fake zones:
CALL BOUNDARY( XL(1-NBDY), A_AVG(1-NBDY), IBDY_TYPE, N,NBDY )
c
c
c While a uniform grid is assumed, it is convient to calculate
c the individual zone widths.
DO I=1-NBDY,N+NBDY
DX(I) = XL(I+1) - XL(I)
ENDDO
c
c Now determine the offset between the two grids. For stability, we require
c that the offset be no greater than one zone. That is
c XL(I-1) <= XLNU(I) <= XL(I+1), when express as a fraction COURNO
c
c COURNO = ABS(XL(I) - XLNU(I))/DX
c this requirement becores
c
c 0 <= COURNO(I) <= 1.0,
c
c A familiar Courant condition or stability condition.
c
c
1 PRINT *,' Enter the fractional offset (between -1.0 and 1.0)'
PRINT *,' between the two grids: DXF. '
READ *,DXF
IF( DXF .LT. -1.0 .OR. DXF .GT. 1.0 ) THEN
PRINT *,' You entered',DXF
GOTO 1
ENDIF
c
c Determine the locations of the lefthand zone interface of the new grid
c (The first and last interface locations will not be needed)
c
DO I=2-NBDY,N+NBDY
IF (DXF .LE. 0.0) THEN
XLNU(I) = XL(I) + DXF * DX(I-1)
ELSE
XLNU(I) = XL(I) + DXF * DX(I)
ENDIF
ENDDO
c Fill an array, SMALLDA, which contains what is to be considered a trivial difference of A_AVG.
c In this Program, A_AVG is not positive definite, so a small value is appropriate.
DO I=1-NBDY,N+NBDY
SMALLDA(I) = 0.0005
ENDDO
c
c DA_AVL(I), is defined for 5-NBDY<= I <= N+NBDY-4
nnbdy = NBDY-3
DO I=1-nnbdy,N+nnbdy
DVOLL(i) = XL(I) - XLNU(I)
IF (DVOLL(I) .GE. 0.0) THEN
c The data will come from the
c zone (i-1) at the left,
c
XF(I) = DVOLL(I)/DX(i-1)
c
ELSE
c
c The data will come from zone (i)
XF(I) = -DVOLL(I)/DX(I)
c
ENDIF
ENDDO
mbdy = 1
c
c Determinve fluxes
c
CALL PPM98_PASSIVE_FLUX0
^ ( XF(1-NBDY), DVOLL(1-NBDY), A_AVG(1-NBDY),
^ SMALLDA(1-NBDY), UNSMOOTH(1-NBDY), DISCONT(1-NBDY),
^ DA_AVL(1-NBDY), N, NBDY )
c
c Now use the fluxes DA_AVL to find the new zone averages
c
DO I=5-NBDY,N+NBDY-4
DXNU = XLNU(I+1)-XLNU(I)
A_AVGNU(I) = ( A_AVG(I)*DX(I) + DA_AVL(I)-DA_AVL(i+1))/DXNU
ENDDO
c
c Display old and new zone averages
PRINT *,' I XMID A_AVG XMIDNU A_AVGNU'
DO I=5-NBDY,N+NBDY-4
XMID = 0.5*(XL(I+1)+XL(I))
XMIDNU = 0.5*(XLNU(I+1)+XLNU(I))
WRITE(6,100) I, XMID, A_AVG(I), XMIDNU, A_AVGNU(I)
100 FORMAT(3x,i4,1p,4(1x,e12.5))
ENDDO
DEALLOCATE (XL, A_AVG, AL, AR, DA, A6, SMALLDA, UNSMOOTH, DISCONT)
DEALLOCATE (DX, A_AVGNU, XLNU, DA_AVL, DVOLL, XF, XF1)
END ! PROGRAM MAP_3