************************************************************************
*****                                                              *****
*****                       PORCYL                                 *****
*****                       ______                                 *****
*****                                                              *****
************************************************************************
*****                                                              *****
*****  THIS PROGRAM SOLVES FOR BOUNDARY LAYER FLOW OVER A CYLINDER *****
*****                                                              *****
*****  IN A POROUS MEDIUM USING THE FINITE DIFFERENCE TECHNIQUE.   *****
*****                                                              *****
*****     A CONSTANT (I.E. UNIFORM) Y-GRID SPACING IS USED.        *****
*****                                                              *****
*****      THE WALL TEMPERATURE IS ASSUMED TO BE UNIFORM           *****
*****                                                              *****
************************************************************************
C
      DIMENSION V(500),T(2,500),E(500),F(500)
      DIMENSION G(500),H(500),Y(500),A(500)
      REAL ND,NDM
*
**********************************************************************
*
      OPEN(1,FILE='PORCYPRT.DAT')
      OPEN(2,FILE='PORCYPLT.DAT')
*
**********************************************************************
*
      WRITE (1,1000)
      WRITE (1,1001)
      WRITE(6,6000)
      WRITE(6,1000)
      WRITE(6,1001)
      WRITE(6,6000)
*
      SND=0.0
      X = 0.0
      N = 50
      N1=N-1
      Y(1) = 0.0
      DY = 0.050
*
      DO 1 J = 2,N
        Y(J) = Y(J-1) + DY
 1    CONTINUE
      DX = 0.001
      DXMAX=0.005
      XMAX=1.5707
*
*************************** SOLUTION BEGINS **************************
*
 100  CONTINUE
*
      WRITE (1,2000) X
*
*
      CALL VELDV (X,U,DUX)
*
****************** SOLVE ENERGY EQUATION TO GET "T" ********************
*
      DO 12 J = 1,N
        V(J)= -Y(J)*DUX
        E(J) = (2.0/(DY*DY))+(U/DX)
        F(J) = V(J)/(2.0*DY)-1.0/(DY*DY)
        G(J) = -V(J)/(2.0*DY)-1.0/(DY*DY)
        H(J) = U*T(1,J)/DX
 12   CONTINUE
*
      E(1)=1.0
      F(1)=0.0
      G(1)=0.0
      H(1)=1.0
      E(N)=1.0
      F(N)=0.0
      G(N)=0.0
      H(N)=0.0
*
      CALL TRISOL(N,E,F,G,H,A)
*
*
      DO 13 J = 1,N
        T(2,J) = A(J)
 13   CONTINUE
*
************* FIND THE BOUNDARY LAYER THICKNESSES ********************
*
      TT = 0.01*T(2,1)
      DO 22 J=1,N1
        IF (T(2,J).LE.TT) GO TO 23
  22  CONTINUE
  23  DE = Y(J)-(Y(J)-Y(J-1))*(TT-T(2,J))/(T(2,J-1)-T(2,J))
*
************* CHECK POSITION OF OUTER GRID POINTS ********************
*
      IF (DE.LE.Y(N-5)) GO TO 24
      NMIN = N+1
      NMAX = N+10
      IF (NMAX.GT.500) GO TO 101
      DO 25 I=NMIN,NMAX
        Y(I) = Y(I-1) + DY
        T(2,I) = 0.0
  25  CONTINUE
      WRITE (1,3000)
      WRITE (6,3000)
      N = NMAX
  24  CONTINUE
*
* *************************** RETURN VALUES ****************************
*
      DO 30 J=1,N
      T(1,J)=T(2,J)
  30  CONTINUE
*
************************** FIND HEAT TRANSFER *************************
*
      ND = (T(2,1)-T(2,2))/(Y(2)-Y(1))
      IF( X.GT.0.0) THEN
       SND=SND+(ND+NDM)*DX/2.0
      ENDIF
      NDM=ND
*
************************** WRITE THE RESULTS **************************
*
      WRITE(1,4000) ND,DE
      WRITE(6,4200) X,ND
      WRITE(2,4100) X,ND,DE
      WRITE(1,4001)
      DO 40 J=1,N
        WRITE(1,4002)J,Y(J),T(2,J)
 40   CONTINUE
*
************************** CHANGE DX **********************************
*
      IF (DX.LT.DXMAX) THEN
        DX = 1.05*DX
      ELSE
        DX = DXMAX
      ENDIF
*
      X = X + DX
*
************************** CHECK FOR MAXIMUM X ************************
*
      IF(X.GT.XMAX) GO TO 102
      GO TO 100
 101  WRITE (1,5000)
      WRITE (6,5000)
      GO TO 103
 102  WRITE(1,5001)
      WRITE(6,5001)
      WRITE(1,5500) SND/XMAX
      WRITE(6,5500) SND/XMAX
 103  CONTINUE
      CLOSE(1)
      CLOSE(2)
      STOP
*
************************** FORMAT STATEMENTS **************************
*
 1000 FORMAT (10X,'BOUNDARY LAYER FLOW OVER AN ISOTHERMAL CYLINDER')
 1001 FORMAT (10X,'-----------------------------------------',//)
 2000 FORMAT (/,5X,'SOLUTION FOR X =',F12.4,/
     &         ,5X,'___________________________',/)
 3000 FORMAT (5X,'****** NUMBER OF GRID POINTS INCREASED ******')
 4000 FORMAT (5X,'ND/PD 1/2 =',E12.4,' DELT =',E12.4,/)
 4001 FORMAT (2X,'J',7X,'Y',12X,'T')
 4002 FORMAT (I4,2E12.4)
 4100 FORMAT (F10.5,',',F10.5,',',F10.5)
 4200 FORMAT (5X,' X =',F10.4,' ND/PD 1/2 =',F10.4)
 5000 FORMAT (5X,'N EXCEEDS 500')
 5001 FORMAT (/,5X,'********** XMAX REACHED **********')
 5500 FORMAT (5X,' MEAN ND/PD 1/2 =',F10.4)
 6000 FORMAT(////////)
      END
*
*
********************************************************************
*
      SUBROUTINE TRISOL(NN,A,B,C,D,H)
*
*********** THIS IS A TRI-DIAGONAL MATRIX SOLVER *******************
*                     ---                 ---
*
*   THIS TRIDIAGONAL MATRIX SOLVER USES THE THOMAS ALGORITHM
*
      DIMENSION A(500),B(500),C(500),D(500),H(500),W(500),Q(500),G(500)
      W(1)=A(1)
      G(1)=D(1)/W(1)
      DO 2 K=2,NN
        K1=K-1
        Q(K1)=B(K1)/W(K1)
        W(K)=A(K)-C(K)*Q(K1)
        G(K)=(D(K)-C(K)*G(K1))/W(K)
    2 CONTINUE
      H(NN)=G(NN)
      N1=NN-1
      DO 3 K=1,N1
        KK=NN-K
        H(KK)=G(KK)-Q(KK)*H(KK+1)
    3 CONTINUE
      RETURN
      END
*
********************************************************************
*
      SUBROUTINE VELDV (X,U,DUX)
*
****** THIS DETERMINES X-WISE VELOCITY AND VELOCITY GRADIENT ******
*
      U = 2.0*SIN(2.0*X)
      DUX = 4.0*COS(2.0*X)
      RETURN
      END