	PROGRAM CIP01

	INTEGER :: NX,NT,I,K,ISG,IUP 
	PARAMETER(NX=100,NT=100) 
	REAL :: ZZ,XX,YU(NX),DT,DX,FDIF,XAM1,XBM1,U 
	REAL :: F(NX),G(NX),FN(NX),GN(NX) 
	PARAMETER (DT=0.002,DX=0.01) 
	COMMON /BL1/U 

	DO 50 I=1,NX
		YU(I)=1.0
		IF ((I.GE.34).AND.(I.LE.46)) THEN
		   F(I)=1.0
		ELSE IF ((I.EQ.33).OR.(I.EQ.47)) THEN
		   F(I)=0.75
		ELSE IF ((I.EQ.32).OR.(I.EQ.48)) THEN
		   F(I)=0.5
		ELSE IF ((I.EQ.31).OR.(I.EQ.49)) THEN
		   F(I)=0.25
		ELSE
		   F(I)=0.0
		ENDIF	
	IF(I.EQ.1) THEN
	   G(I)=0.0
	ELSE
  	   G(I)=(F(I)-F(I-1))/DX
	ENDIF	
	WRITE(*,*) 'F(',I,')=',F(I)
50	CONTINUE

	OPEN(UNIT=10,FILE='pulse_fi.dat',STATUS='UNKNOWN')

	DO 300 I=1,NX
		WRITE(10,310) I,F(I)
300	CONTINUE

310	FORMAT(1X,I6,3X,E12.5)

	CLOSE(UNIT=10,STATUS='KEEP')


	DO 200 K=1,NT
		WRITE(*,*) 'K=',K
	DO 100 I=2,NX
		ZZ=1.0
		U=YU(I)
		CALL SIGN(ZZ)
		ISG=-ZZ
		IUP=I+ISG
		   XX=-YU(I)*DT
		   FDIF=(F(IUP)-F(I))/DX*ISG
		   XAM1=(G(I)+G(IUP)-2.0*FDIF)/(DX*DX)
		   XBM1=(3.0*FDIF-2.0*G(I)-G(IUP))/DX*ISG
		FN(I)=((XAM1*XX+XBM1)*XX+G(I))*XX+F(I)
		GN(I)=(3.0*XAM1*XX+2.0*XBM1)*XX+G(I)
		F(I)=FN(I)
		G(I)=GN(I)
100	CONTINUE
200	CONTINUE

	OPEN(UNIT=10,FILE='pulse_fo.dat',STATUS='UNKNOWN')

	DO 320 I=1,NX
		WRITE(10,330) I,F(I)
320	CONTINUE

330	FORMAT(1X,I6,3X,E12.5)

	CLOSE(UNIT=10,STATUS='KEEP')

	STOP
	END

	SUBROUTINE SIGN(ZZ)
	REAL :: ZZ,U
	
	IF (U.GE.0.0) THEN
	   ZZ=ZZ
	ELSE IF (U.LT.0.0) THEN
	   ZZ=-ZZ
	ENDIF

	RETURN 
	END


