!************************************************************
!***Nonlinear Function Optimization (Minimization) Routine
!***This operates on a subset of the total model parameters 
!***(PARMMODEL,NPARM) that will be optimized (OPTPARM,NOPTPARM), 
!***the other parameters being fixed because at min/max limits 
!***or fixed as part of searches for their confidence limits 
!************************************************************
	SUBROUTINE OPTIMIZE()
	USE SHAREDVARS
	REAL*8 FUNCCALC,PARMLIM(4)
	INTEGER IPOS,IOPTPARM,JOPTPARM,IERR,ITER
!***
!***Reset Variables
!***
	DO IPOS=1,4
		FUNC(IPOS)=0.0
		STEP(IPOS)=0.0
	ENDDO
	DO IPARM=1,4
		DFUNC(IPARM)=0.0
		VECTOR(IPARM)=0.0
		DO JPARM=1,4
			DFUNC2(IPARM,JPARM)=0.0
			MATRIX(IPARM,JPARM)=0.0
		ENDDO
	ENDDO
	DO IERR=0,8
		OPTIMIZERROR(IERR)=0
	ENDDO
!***
!***Specify initial function value and derivatives, and store in history for IPHASE=0 
!***If NOPTPARM=NPARM, search is for all parameters with specified initial guesses
!***If NOPTPARM<NPARM, search is for CL of a fixed parameter, so optimization is just for other parameters with their guesses.  
!***
	FUNC(0)=FUNCCALC(OPTPARMSEARCH(1,0))
	CALL DFUNCCALC(OPTPARMSEARCH(1,0))
	IF(IPHASE.EQ.0) THEN
          	HISTFUNC(0)=FUNC(0)
          	DO IOPTPARM=1,NOPTPARM
          		HISTPARM(0,IOPTPARM)=OPTPARMSEARCH(IOPTPARM,0)
          		HISTDFUN(0,IOPTPARM)=DFUNC(IOPTPARM)
          	ENDDO
	ENDIF
	IF(PRINTFLAG) TYPE*,'*****IN OPTIMIZATION - INITIAL CONDITIONS*****'
	IF(PRINTFLAG) TYPE"('PARMSEARCH  ',4G15.6)",(OPTPARMSEARCH(IOPTPARM,0),IOPTPARM=1,NOPTPARM)
	IF(PRINTFLAG) TYPE"('INITIAL FUNC',G15.6)",FUNC(0)
	IF(PRINTFLAG) TYPE"('FIRSTDERIV  ',4G15.6)",(DFUNC(IOPTPARM),IOPTPARM=1,NOPTPARM)
	DO IOPTPARM=1,NOPTPARM
		IF(PRINTFLAG) TYPE"('SECONDDERIV ',4G15.6)",(DFUNC2(IOPTPARM,JOPTPARM),JOPTPARM=1,NOPTPARM)
	ENDDO
!***
!***Main Optimization Iteration Loop
!***
	DO ITER=1,MAXIT
		IF(PRINTFLAG2) TYPE"('*****IN OPTIMIZE - LOOP',I3,' *****',G15.6)",ITER,FUNC(0)
!***Determine search direction
		IF(PRINTFLAG) TYPE"('CALLING NEWTON:')"
		CALL NEWTON
		IF(PRINTFLAG) TYPE"('EXITING NEWTON, VECTOR IS:',3F15.6)",(VECTOR(IOPTPARM),IOPTPARM=1,NOPTPARM)
!***If parameter 1 already at max and vector would increase it, recalculate vector for remaining parameters and set VECTOR(1) to 0
		IF(OPTPARMSEARCH(1,0).GE.0.9998D0.AND.VECTOR(1).GE.0D0) THEN
			IF(PRINTFLAG2) TYPE"('Y0 AT MAX - REDUCE OPTIMIZATION PARAMETERS')"
			NOPTPARM=NOPTPARM-1
			FIXPARMLIMID=1
			FIXPARMLIMVAL=OPTPARMSEARCH(1,0)
			DO IOPTPARM=1,NOPTPARM  !Shift parameter ids and values downward to form reduced set
				OPTPARMID(IOPTPARM)=OPTPARMID(IOPTPARM+1)
				PARMLIM(IOPTPARM)=OPTPARMSEARCH(IOPTPARM+1,0)
			ENDDO
			IF(PRINTFLAG) TYPE"('MODIFIED,REDUCED OPTPARM IS',3F15.6)",(PARMLIM(IOPTPARM),IOPTPARM=1,NOPTPARM)
			CALL DFUNCCALC(PARMLIM)
			CALL NEWTON
			IF(PRINTFLAG) TYPE"('MODIFIED REDUCED VECTOR IS:',3F15.6)",(VECTOR(IOPTPARM),IOPTPARM=1,NOPTPARM)
			NOPTPARM=NOPTPARM+1
			DO IOPTPARM=NOPTPARM,2,-1   !Shift parameter ids and values back up and add back fixed parameter as vector(1)=0   
				VECTOR(IOPTPARM)=VECTOR(IOPTPARM-1)  
				OPTPARMID(IOPTPARM)=OPTPARMID(IOPTPARM-1)
			ENDDO
			VECTOR(1)=0D0
			OPTPARMID(1)=1
			FIXPARMLIMID=0
			IF(PRINTFLAG) TYPE"('MODIFIED FULL SIZED VECTOR IS:',3F15.6)",(VECTOR(IOPTPARM),IOPTPARM=1,NOPTPARM)
		END IF
!***Find minimum along search direction
		IF(PRINTFLAG) TYPE"('CALLING LINEAR SEARCH')"		
		CALL LINEAR
		IF(PRINTFLAG) TYPE"('LEAVING LINEAR SEARCH')"
		IF(PRINTFLAG) TYPE"('STEP',6F15.6)",STEP
		DO IOPTPARM=1,NOPTPARM
			IF(PRINTFLAG) TYPE"('PARM',6F15.6)",(OPTPARMSEARCH(IOPTPARM,IPOS),IPOS=0,5)
		ENDDO
		IF(PRINTFLAG) TYPE"('FUNC',5F15.6)",(FUNC(IPOS),IPOS=0,4)
!***Determine First and Second Derivatives at Current Parameter Values
		CALL DFUNCCALC(OPTPARMSEARCH(1,BEST))
		IF(PRINTFLAG) TYPE"('NEW FIRSTDERIV ',4G15.6)",(DFUNC(IOPTPARM),IOPTPARM=1,NOPTPARM)
		DO IOPTPARM=1,NOPTPARM
			IF(PRINTFLAG) TYPE"('NEW SECONDDERIV',4G15.6)",(DFUNC2(IOPTPARM,JOPTPARM),JOPTPARM=1,NOPTPARM)
		ENDDO
!***Save iphase=0 iteration history for output
		IF(IPHASE.EQ.0) THEN
          		HISTFUNC(ITER)=FUNC(BEST)
          		DO IOPTPARM=1,NOPTPARM
          			HISTPARM(ITER,IOPTPARM)=OPTPARMSEARCH(IOPTPARM,BEST)
          			HISTDFUN(ITER,IOPTPARM)=DFUNC(IOPTPARM)
          		ENDDO
		ENDIF	
!***Terminate if convergence requirement met and store best estimates for next iteration or terminatio
		CONVERGE=.TRUE.
		IF(FUNC(BEST).GT.1D-20) THEN
			IF(ABS((FUNC(BEST)-FUNC(0))/FUNC(BEST)).GT.CONVRGREQ) CONVERGE=.FALSE.
     			IF(OPTPARMSEARCH(1,BEST).LT.0.9998.AND.ABS(DFUNC(1)/FUNC(BEST)).GT.CONVRGREQ) CONVERGE=.FALSE.
			DO IOPTPARM=2,NOPTPARM 
				IF(ABS(DFUNC(IOPTPARM)/FUNC(BEST)).GT.CONVRGREQ) CONVERGE=.FALSE.
     			ENDDO
		END IF
		CALL SHIFT(BEST,0)
		IF(CONVERGE) THEN
			IF(IPHASE.EQ.0) HISTITER=ITER
			EXIT
		ENDIF
		IF(PRINTFLAG) THEN
			TYPE "('*****END OF OPTIM LOOP',I5,' *****',G15.6)",ITER,FUNC(BEST)
			!ACCEPT *,IGO
		ENDIF
	END DO
	IF(IPHASE.EQ.0) THEN
		IF(ITER.GT.MAXIT) HISTITER=50
          	FUNCSOLUT=FUNC(BEST)
          	DO IOPTPARM=1,NOPTPARM
          		PARMSOLUT(IOPTPARM)=OPTPARMSEARCH(IOPTPARM,BEST)
          	ENDDO
	ENDIF
!***
!***Set Error Code if Maximum Iterations Reached or Parameter Value at Maximum or Minimum
!***
	DO IOPTPARM=1,NOPTPARM
		GAP=1D-6*(PARMMAX(IOPTPARM)-PARMMIN(IOPTPARM))
		IF(OPTPARMSEARCH(IOPTPARM,0).GE.(PARMMAX(OPTPARMID(IOPTPARM))-GAP).OR.OPTPARMSEARCH(IOPTPARM,0).LE.(PARMMIN(OPTPARMID(IOPTPARM))+GAP)) OPTIMIZERROR(IOPTPARM)=1
	END DO
	IF(.NOT.CONVERGE) OPTIMIZERROR(0)=1
!***
!***End Routine and Return to Calling Routine
!***
	RETURN
	END SUBROUTINE OPTIMIZE

