!**************************************************************************************************************
!***Routine for Search Direction by Damped Newton-Raphson Method (Dahlquist and Bjorck "Numerical Methods")***
!**************************************************************************************************************
	SUBROUTINE NEWTON
	USE SHAREDVARS
	REAL*8 SUM1,SUM2,SUM3
	INTEGER ICODE,ITYPE,IOPTPARM,JOPTPARM,MPARM
	LOGICAL OK_FLAG
!***Invert second derivative matrix
	IF(PRINTFLAG) TYPE *,'CALLING INVERT'
     	CALL INVERT(NOPTPARM,DFUNC2,MATRIX,ok_flag)
	IF(PRINTFLAG) TYPE *,OK_FLAG
	DO IOPTPARM=1,NOPTPARM
		IF(PRINTFLAG) TYPE "(<NOPTPARM>G15.8,5X,<NOPTPARM>G15.8)",(DFUNC2(IOPTPARM,JOPTPARM),JOPTPARM=1,NOPTPARM),(MATRIX(IOPTPARM,JOPTPARM),JOPTPARM=1,NOPTPARM)
	ENDDO
!***Set Search Vector to Newton Direction if Descent; Otherwise to Steepest Descent Direction 
	SUM1=0
	SUM2=0
	SUM3=0
	DO IOPTPARM=1,NOPTPARM
		SUM1=SUM1+DFUNC(IOPTPARM)**2
		DO JOPTPARM=1,NOPTPARM
			SUM2=SUM2+DFUNC(IOPTPARM)*DFUNC(JOPTPARM)*DFUNC2(IOPTPARM,JOPTPARM)
			SUM3=SUM3+DFUNC(IOPTPARM)*DFUNC(JOPTPARM)*MATRIX(IOPTPARM,JOPTPARM)
		END DO
	END DO
	IF(SUM3.GT.MAX(SUM1**2/SUM2,0.0)) THEN
		IF(PRINTFLAG) TYPE *,'NEWTON DIRECTION'
		DO IOPTPARM=1,NOPTPARM
			VECTOR(IOPTPARM)=0.0
			DO JOPTPARM=1,NOPTPARM
				VECTOR(IOPTPARM)=VECTOR(IOPTPARM)-MATRIX(IOPTPARM,JOPTPARM)*DFUNC(JOPTPARM)
			END DO
		END DO
		STEP(3)=1.0
	ELSE
		IF(PRINTFLAG) TYPE *,'STEEPEST DESCENT'
		DO IOPTPARM=1,NOPTPARM
			VECTOR(IOPTPARM)=-DFUNC(IOPTPARM)
		END DO
		IF(SUM3.GT.0) THEN
			STEP(3)=SUM1/SUM2
		ELSE
			STEP(3)=1
		ENDIF
	END IF
	END SUBROUTINE NEWTON
