SUBROUTINE INITAL C... C... M1 - BLASIUS PROBLEM C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ ETA, NSTOP, NORUN + /Y/ G1, G2, G3 + /F/ G1E, G2E, G3E + /S/ G3S, DG3S, FLAGM, FLAGP + /I/ IP C... C... SET THE INITIAL ESTIMATE OF G3, INCREMENT IN G3, C... CONTROL VARIABLES FOR ITERATION ON G3 IF(NORUN.EQ.1)THEN G3S =0.3D0 DG3S=0.2D0 FLAGM=-1.0 FLAGP=-1.0 END IF C... C... INITIAL CONDITIONS G1=0.0D0 G2=0.0D0 G3=G3S C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ ETA, NSTOP, NORUN + /Y/ G1, G2, G3 + /F/ G1E, G2E, G3E + /S/ G3S, DG3S, FLAGM, FLAGP + /I/ IP C... C... ODES G1E=G2 G2E=G3 G3E=-0.5D0*G1*G3 RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ ETA, NSTOP, NORUN + /Y/ G1, G2, G3 + /F/ G1E, G2E, G3E + /S/ G3S, DG3S, FLAGM, FLAGP + /I/ IP C... C... INCREMENT THE COUNTER TO DETERMINE THE END OF EACH RUN IP=IP+1 C... C... PRINT THE SOLUTION FOR THE FIRST AND LAST RUNS IF((NORUN.EQ.1).OR.(NORUN.EQ.20))THEN C... C... PRINT A HEADING IF(IP.EQ.1)WRITE(NO,1)NORUN 1 FORMAT(//,' NORUN = ',I3,//, + 7X,'eta',8X,'g1',8X,'g2',8X,'g3') C... C... PRINT THE SOLUTION WRITE(NO,2)ETA,G1,G2,G3 2 FORMAT(F10.2,3F10.5) END IF C... C... AT THE END OF EACH RUN, INCREMENT G3 DEPENDING ON THE SIGN C... DEVIATION FROM THE PRESCRIBED BOUNDARY CONDITION G2(INF) C... = 1 IF(IP.EQ.11)THEN CALL SEARCH C... C... MONITOR THE ITERATION WRITE(*,3)NORUN,G2,G3S,DG3S,FLAGP,FLAGM 3 FORMAT(I5,3F12.6,2F10.3) END IF RETURN END SUBROUTINE SEARCH C... C... SUBROUTINE SEARCH SEARCHES FOR A ROOT OF A FUNCTION BY C... INTERVAL HALVING C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ ETA, NSTOP, NORUN + /Y/ G1, G2, G3 + /F/ G1E, G2E, G3E + /S/ G3S, DG3S, FLAGM, FLAGP + /I/ IP C... C... IF (1 - G2(INF)) IS NONNEGATIVE, INCREASE G3(0) IF(FUNC(G2).GE.0.0D0)THEN C... C... IF A SMALLER INCREMENT IN G3(0) IS REQUIRED, HALVE C... DG3S IF(FLAGM.GE.0.0D0)DG3S=DG3S/2.0D0 C... C... INCREMENT G3(0) POSITIVELY FOR THE NEXT RUN G3S=G3S+DG3S FLAGP=1.0D0 C... C... IF (1 - G2(INF)) IS NEGATIVE, DECREASE G3(0) ELSE IF(FUNC(G2).LT.0.0D0)THEN C... C... IF A SMALLER INCREMENT IN G3(0) IS REQUIRED, HALVE C... DG3S IF(FLAGP.GE.0.0D0)DG3S=DG3S/2.0D0 C... C... INCREMENT G3(0) NEGATIVELY FOR THE NEXT RUN G3S=G3S-DG3S FLAGM=1.0D0 END IF RETURN END DOUBLE PRECISION FUNCTION FUNC(G2) C... C... FUNCTION FUNC DEFINES THE FUNCTION FOR WHICH A ZERO IS TO BE C... COMPUTED C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) FUNC=1.0D0-G2 RETURN END