PROGRAM CREATE_ZONE_AVERAGES_TANH
c
c   This program first creates a computational grid or mesh of N
c   zones .The grid begins at 0.0 and has a total length of 1.0.
c   It consists of N+1 lefthand interface locations.
c
c   Zone averages of a discontinuous-like distribution are created
c   by finding  the average value of tanh(a(x-x0)) in each zone.  the 
c   offset x0 is chosen such to ensure that the function is centered, 
c   and a is chosen such that 90% of the change in  tanh occurs over
c   a specified number of zones.
c
c   A file "OUTPUT.DAT" is created which may be used as input for
c   subsequent example programs.

          ALLOCATABLE ::  XL(:), A_AVG(:)

c
c   The number of grid zones is selected interactively
c
      PRINT *,' Please enter desired number of grid zones (integer):'
      READ *,N
c
c   The width of the discontinuity in zones
c
      PRINT *,' Kindly enter the approximate width of the '
      PRINT *,' discontinuity in zones'
      READ *,NWIDTH
c
c   The type of grid (uniform, non-uniform, discontinuous) is
c                                              selected interactively
      PRINT *,' Please enter -1 for a discontinuous grid'
      PRINT *,'               0 for a non-uniform grid or'
      PRINT *,'              +1 for a uniform grid:'
      READ *,IGRID


      PRINT *,' Thank you. '
      PRINT *,'   You have requested a discontinuity of approx.',
     ^        NWIDTH,' zones. The total number of grid zones is',N

          ALLOCATE (XL(N+1), A_AVG(N))

          X_BEGIN = 0.0
          X_LENGTH = 1.0


      IF (IGRID .EQ. 1 ) THEN
          PRINT *,' on a uniform grid'
          CALL CREATE_UNIFORM_GRID( XL(1), X_BEGIN, X_LENGTH, N )
      ELSE
         IF (IGRID .EQ. 0 ) THEN
             PRINT *,' on a non-uniform grid'
             CALL CREATE_NONUNIFORM_GRID( XL(1), X_BEGIN, X_LENGTH, N )
         ELSE
            IF (IGRID .EQ. -1 ) THEN
                PRINT *,' on a discontinuous-uniform grid'
               CALL CREATE_DISCONTINUOUS_GRID( XL(1), X_BEGIN,
     ^                                                  X_LENGTH, N )
            ELSE
c
c   unknown grid type
               PRINT *,' UNKNOWN GRID TYPE:',IGRID
               CALL ABORT()
            ENDIF
         ENDIF
      ENDIF



c
c   The functional form is tanh(alpha*(x-x0))
c
c   The average of this function on the interval (x1,x2) is
c
c              [ln(cosh(alpha*(x2-x0))) - 
c                                 ln(cosh(alpha*(x2-x0)))]/[alpha*(x2-x1)]
c
      X0 = 0.5
c
c    90% of the variation of tanh occurs on the range (-1.5:1.5)
c
      NHALF = N/2
      NZL = NWIDTH/2
      NZR = NWIDTH-NZL
      DELTAX = XL(NHALF+NZL) - XL(NHALF-NZR) 
      ALPHA = 3.0 / DELTAX
          DO I = 1,N
             XLEFT = ALOG( COSH( ALPHA*(XL(I)-X0) ) )
             XRIGHT = ALOG( COSH( ALPHA*(XL(I+1)-X0) ) )
             DX = XL(I+1)-XL(I)

             A_AVG(I) =  (XRIGHT - XLEFT)/(ALPHA*DX)
          ENDDO
          OPEN (UNIT=11, FILE='OUTPUT.AVG')
	  write(11,50) N
50    FORMAT(1X,I10)

          WRITE (6,100)
100   FORMAT(6X,'I',6X,'XMID',10X,'A_AVG',/,40('_'))

          DO I=1,N
            XMID = 0.5 * (XL(I+1)+XL(I))
            WRITE (6,200) I, XMID, A_AVG(I)
200     FORMAT (4X,I4,t10,1P,2(E12.5,3x))
            WRITE (11,300)  XL(I), A_AVG(I)
300     FORMAT (4X,1P,2(E12.5,3x))
          ENDDO
            WRITE (11,300)  XL(N+1), A_AVG(N)
            CLOSE(11)

          DEALLOCATE (XL, A_AVG)

          END ! PROGRAM CREATE_ZONE_AVERAGES_TANH