************************************************************************
*****                                                              *****
*****                       SIMVART                                *****
*****                       _______                                *****
*****                                                              *****
************************************************************************
*****                                                              *****
*****  THIS PROGRAM GIVES SIMILARITY SOLUTION RESULTS FOR          *****
*****                                                              *****
*****  LAMINAR BOUNDARY LAYER FLOW OVER A FLAT PLATE WHOSE         *****
*****                                                              *****
*****  SURFACE TEMPERATURE IS ( T1 + C x**n )                      *****
*****                                                              *****
************************************************************************
*
***********************************************************************
*
      DIMENSION T(1500),F(1500),G(1500),H(1500)
     & ,HGUESS(1500),FPIN(1500),P(1500)
      REAL KF0,KG0,KH0,KF1,KG1,KH1,KF2,KG2,KH2,KF3,KG3,KH3
     & ,KT0,KT1,KT2,KT3
      COMMON DETA,XN,PR
*
***********************************************************************
*
      OPEN(UNIT=1,FILE='SIMVARPR.DAT')
      OPEN(UNIT=2,FILE='SIMVARPL.DAT')
*
***********************************************************************
*
      WRITE(6,4250)
      WRITE(6,4290)
      WRITE(6,4260)
      READ(5,*) PR
      WRITE(6,4250)
      WRITE(6,4265)
      READ(5,*) XN
      WRITE(6,4250)
      DETA=0.01
      N=1500
      WRITE(1,4290)
      WRITE(1,4270) PR
*
***********************************************************************
*
*     N=NUMBER OF GRID POINTS
*     PR=PRANDTL NUMBER
*     DETA=ETA STEP SIZE
*
***********************************************************************
*
      INTR=1
      HGUESS(INTR)=1.0
  100 F(1)=0.0
      G(1)=0.0
      H(1)=HGUESS(INTR)
      DO 1000 I=2,N
      F0=F(I-1)
      G0=G(I-1)
      H0=H(I-1)
      FC=F0
      GC=G0
      HC=H0
      CALL RUTKUN(FC,GC,HC,KF0,KG0,KH0)
      FC=F0+0.5*KF0
      GC=G0+0.5*KG0
      HC=H0+0.5*KH0
      CALL RUTKUN(FC,GC,HC,KF1,KG1,KH1)
      FC=F0+0.5*KF1
      GC=G0+0.5*KG1
      HC=H0+0.5*KH1
      CALL RUTKUN(FC,GC,HC,KF2,KG2,KH2)
      FC=F0+KF2
      GC=G0+KG2
      HC=H0+KH2
      CALL RUTKUN(FC,GC,HC,KF3,KG3,KH3)
      F(I)=F(I-1)+(KF0+2.0*KF1+2.0*KF2+KF3)/6.0
      G(I)=G(I-1)+(KG0+2.0*KG1+2.0*KG2+KG3)/6.0
      H(I)=H(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0
 1000 CONTINUE
      FPIN(INTR)=G(N)
      IF (INTR.EQ.1) THEN
        INTR=2
        HGUESS(INTR)=HGUESS(1)-0.0001
        GO TO 100
      ELSE
        IF(INTR.GT.200) GO TO 200
        IF(ABS(1.0-FPIN(INTR)).LT.0.0000005) GO TO 300
        INTR=INTR+1
        DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)-
     &  HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1))
        HGUESS(INTR)=HGUESS(INTR-1)+DGUESS
        GO TO 100
      ENDIF
  200 WRITE(6,2000)
      GO TO 7777
  300 WRITE(6,3000) INTR
 7777 CONTINUE
      HU=H(1)
      INTR=1
      HGUESS(INTR)=G(1)
  101 T(1)=0.0
      H(1)=HGUESS(INTR)
      DO 1001 I=2,N
      FC=(F(I)+F(I-1))/2.0
      GC=(G(I)+G(I-1))/2.0
      T0=T(I-1)
      H0=H(I-1)
      TC=T0
      HC=H0
      CALL RUTTUN(TC,HC,KT0,KH0,FC,GC)
      TC=T0+0.5*KT0
      HC=H0+0.5*KH0
      CALL RUTTUN(TC,HC,KT1,KH1,FC,GC)
      TC=T0+0.5*KT1
      HC=H0+0.5*KH1
      CALL RUTTUN(TC,HC,KT2,KH2,FC,GC)
      TC=T0+KT2
      HC=H0+KH2
      CALL RUTTUN(TC,HC,KT3,KH3,FC,GC)
      T(I)=T(I-1)+(KT0+2.0*KT1+2.0*KT2+KT3)/6.0
      H(I)=H(I-1)+(KH0+2.0*KH1+2.0*KH2+KH3)/6.0
 1001 CONTINUE
      FPIN(INTR)=T(N)
      IF (INTR.EQ.1) THEN
        INTR=2
        HGUESS(INTR)=HGUESS(1)-0.0001
        GO TO 101
      ELSE
        IF(INTR.GT.200) GO TO 201
        IF(ABS(1.0-FPIN(INTR)).LT.0.00001) GO TO 301
        INTR=INTR+1
        DGUESS=+(FPIN(INTR-1)-1.0)*(HGUESS(INTR-1)-
     &  HGUESS(INTR-2))/(FPIN(INTR-2)-FPIN(INTR-1))
        HGUESS(INTR)=HGUESS(INTR-1)+DGUESS
        GO TO 101
      ENDIF
  201 WRITE(6,2000)
      GO TO 7778
  301 WRITE(6,3001) INTR
 7778 CONTINUE
      WRITE(1,4500)
      DO 7000 I=1,N
      ETA=(I-1)*DETA
      WRITE(1,4000) ETA,F(I),G(I),H(I),T(I)
      WRITE(2,5500) ETA,G(I),T(I)
 7000 CONTINUE
      TGRAD=(T(2)-T(1))/DETA
      WRITE(1,4280) TGRAD
      WRITE(6,4280) TGRAD
      WRITE(1,4281) HU
      WRITE(6,4281) HU
      STOP
      CLOSE(1)
      CLOSE(2)
 2000 FORMAT(' ******** FAILURE TO CONVERGE *********')
 3000 FORMAT(' VELOCITY FUNCTION CONVERGENCE IN ',I5,' ITERATIONS')
 3001 FORMAT(' TEMPERATURE FUNCTION CONVERGENCE IN ',I5,' ITERATIONS')
 4000 FORMAT(5F10.6)
 4250 FORMAT(////)
 4260 FORMAT(' INPUT THE VALUE OF THE PRANDTL NUMBER THEN ')
 4265 FORMAT(' INPUT THE VALUE OF n THEN ')
 4270 FORMAT(//,' PRANDTL NUMBER = ',F12.3,//)
 4280 FORMAT(//,' THETA GRADIENT AT WALL = ',F9.3)
 4281 FORMAT(//,' F PRIME GRADIENT AT WALL = ',F9.3)
 4290 FORMAT(/,' SIMILARITY SOLUTION FOR FLOW OVER ISOTHERMAL PLATE',/,
     $ ' __________________________________________________',//)
 4500 FORMAT('     ETA        F     dF/DETA   D2F/DETA2   THETA')
 5500 FORMAT(F10.6,',',F10.6,',',F10.6)
      END
*
***********************************************************************
*
*     RUTTA-KUNGE SOLUTION ROUTINE FOR VELOCITY
*
***********************************************************************
*
      SUBROUTINE RUTKUN(F,G,H,KF,KG,KH)
      REAL KF,KG,KH
      COMMON DETA,XN,PR
      KF=DETA*G
      KG=DETA*H
      KH=-0.5*F*H*DETA
      RETURN
      END
*
***********************************************************************
*
*     RUTTA-KUNGE SOLUTION ROUTINE FOR TEMPERATURE
*
***********************************************************************
*
      SUBROUTINE RUTTUN(T,H,KT,KH,F,G)
      REAL KT,KH
      COMMON DETA,XN,PR
      KT=DETA*H
      KH=-0.5*PR*F*H*DETA+XN*G*PR*(T-1.0)*DETA
      RETURN
      END