C***************************************************************

      SUBROUTINE TYPE66(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)

C     Version from: 12/17/89
C     SAMPLE INPUT AND OUTPUT ARE LOCATED AFTER THE VARIABLE LIST
C--------------------------------------------------------------
C     This subroutine represents a model for a separately 
C     exited DC motor (permanent magnet or any other const.
C     current field source) and a series DC motor. It is
C     an analytical model based on fundamental motor equations.
C     For a direct coupled PV system MODE should be set equal to
C     one and current and speed will be the significant inputs.
C     For a system with battery storage MODE should be set equal
C     to zero, i.e. voltage and speed will be significant inputs.
C     To support convergence between the motor component and
C     any motor load component a convergence promotor is included.
C****************************************************************


C****************************************************************
C     Variables:
C     CVISC -- viscous friction [Nm*sec]
C     CHANGE -- auxiliary variables
C     FAIL -- to output error message
C              =1: max. speed is exceeded
C              =2: rated current is exceeded	
C              =3: max. current is exceeded	
C     ETOL,TOL -- tolerances used in the convergence promoter
C     EFLAG -- variable to keep track of significant inputs 
C     FLAG -- determines whether bisection or secant method
C             is being continued
C     G -- objective function
C     GLAST,GNEG,GPOS -- objective functions at the limits
C                        and from the previous timestep
C     I  -- Current through motor [amps]
C     IMAX -- max. current [amps]
C     IRAT -- rated current [amps]
C     KF -- motor constant,relates the permanent magnet flux
C           through the armature current and the emf [Nm/amp]
C     MAF -- mutual inductance [henry]
C     MEMO -- keeps track of the operating mode
C     MODE -- determines whether I or V and respectively N or T
C             are significant inputs
C	      =0: V and N are significant
C	      =1: I and N are significant
C     N  -- motor speed [1/sec]
C     NMAX -- max. speed [1/sec]
C     NN -- relative addresse used in convergence promoter
C     RA -- motor armature resistance [ohms]
C     RF -- field resistance [ohms]
C     TEM -- motor torque (electromagnetical) [Nm]
C     TLOSS -- loss torque
C     TSHAFT -- shaft torque
C     TYPE -- determines type of DC motor
C	      =1: separately excited
C	      =2: series
C     V  -- terminal voltage [volts]
C     X,XLAST,XNEG,XPOS,XNEW -- converging variable - torque 
C***************************************************************
C
C	SAMPLE INPUTS AND OUTPUTS
C
C     INPUTS             SAMPLE VALUES
C     XIN(1)  I              4.5
C     XIN(2)  V            125.0
C     XIN(3)  N             13.0
C     XIN(4)  TSHAFT        21.0
C     XIN(5)  MODE           1.0
C
C     PARAMETERS         SAMPLE VALUES
C     PAR(1)  TYPE           2.0
C     PAR(2)  NMAX          70.0
C     PAR(3)  IRAT          20.0
C     PAR(4)  IMAX          30.0
C     PAR(5)  MAF             .59
C     PAR(6)  CSTAT           .5
C     PAR(7)  CVISC           .02
C     PAR(8)  RA             1.44
C     PAR(9)  RF              .70
C
C     OUTPUTS            SAMPLE VALUES
C     OUT(1)  I              4.5
C     OUT(2)  V            226.5
C     OUT(3)  N             13.0
C     OUT(4)  TSHAFT        11.19
C     OUT(5)  PEL         1019.0
C     OUT(6)  PMECH        913.8
C     OUT(7)  EFF             .8966
C     OUT(8)  FAIL           0.0
C
C*********************************************************************
      IMPLICIT NONE

	DOUBLE PRECISION XIN,OUT
      INTEGER INFO,TYPE,FAIL
      INTEGER NSTORE,IAV,ISTORE	
      INTEGER NN,MEMO,FLAG	
      REAL TIME,T,DTDT,ICNTRL,PAR,S
      REAL V,I,N,TEM,RA,KF,PI,IRAT,TLOSS,TSHAFT
      REAL MAF,RF,IMAX,NMAX,CSTAT,CVISC,MODE
      REAL PEL,PMECH,EFF,TOL,ETOL,CHANGE,EFLAG
      REAL GLAST,GNEG,GPOS,XLAST,X,XNEW,XPOS,XNEG,G	
      DIMENSION XIN(5), OUT(8), PAR(9), INFO(15)

      COMMON /STORE/ NSTORE,IAV,S(5000)	
C**************************************************************

C--------------------------------------------------------------
C**** initial call in simulation
      IF(INFO(7).LT.0)THEN

C****   accurate expression for pi
        PI=ACOS(-1.)

        INFO(6)=8
        INFO(9)=0



C****   SET PARAMETERS
        TYPE=PAR(1)

        IF(TYPE.EQ.1)THEN
          CALL TYPECK(1,INFO,5,8,0)
	ELSE
          CALL TYPECK(1,INFO,5,9,0)
	ENDIF

	NMAX=PAR(2)
	IRAT=PAR(3)
	IMAX=PAR(4)

	IF(TYPE.EQ.1)THEN
          KF=PAR(5)
	ELSE
	  MAF=PAR(5)
	ENDIF

	CSTAT=PAR(6)
	CVISC=PAR(7)
        RA=PAR(8)

	IF(TYPE.GT.1)THEN
	  RF=PAR(9)
	ENDIF

	INFO(10)=11
	IF(TYPE.GT.1)THEN
	  CALL TYPECK(1,INFO,5,9,0)
	ELSE
	  CALL TYPECK(1,INFO,5,8,0)
	ENDIF
	ISTORE=INFO(10)

C****   initialization of storage values
	S(ISTORE)=0.
	S(ISTORE+1)=0.
	S(ISTORE+2)=0.
	S(ISTORE+3)=0.
	S(ISTORE+4)=0.
	S(ISTORE+5)=0.
	S(ISTORE+6)=0.
	S(ISTORE+7)=0.
	S(ISTORE+8)=0.
	S(ISTORE+9)=0.
	S(ISTORE+10)=0.

	TOL=0.001
	ETOL=0.01

	NN=0
	MEMO=0
C---------------------------------------------------------
C**** first and following calls in timestep
      ELSE

	NN=S(ISTORE+7)
	MEMO=S(ISTORE+8)
	FLAG=S(ISTORE+9)
	EFLAG=S(ISTORE+10)

        FAIL=0

C****   set inputs
        I=XIN(1)
        V=XIN(2)
        N=XIN(3)
        TSHAFT=XIN(4)
      	MODE=XIN(5)

	IF(INFO(7).EQ.0)THEN
	  MEMO=INT(MODE)
	  FLAG=0
          NN=0
          IF(MODE.EQ.1.)THEN
      	    EFLAG=I
          ELSE
            EFLAG=V
      	  ENDIF
	ENDIF

C****   if mode changes during timestep, the convergence
C****   promoter is reset to the first step 
	IF(MEMO.NE.INT(MODE))THEN
	  NN=INFO(7)
      	  MEMO=INT(MODE)
	ENDIF
      	
C****   if significant inputs change during timestep, the convergence
C****   promoter is reset to the first step 
      	IF((MODE.EQ.1.).AND.(I.NE.EFLAG))THEN
      	  NN=INFO(7)
          EFLAG=I
      	ELSEIF((MODE.EQ.0.).AND.(V.NE.EFLAG))THEN
      	  NN=INFO(7)
      	  EFLAG=V
      	ENDIF

C---------------------------------------------------------
C**** 	characteristic motor equations

	TLOSS=CSTAT+CVISC*N

      	IF(MODE.EQ.1) THEN
C****     current and speed serve as input, voltage and 
C****     torque are computed
		
	  IF(TYPE.EQ.1)THEN
C****       separately excited motor
            V=KF*2.*PI*N+I*RA
            TEM=KF*I

	  ELSEIF(TYPE.EQ.2)THEN
C****       series motor
	    V=I*(2.*PI*N*MAF+RA+RF)
	    TEM=MAF*I**2
	  ENDIF

        ELSE
C****     speed and voltage serve as input, torque and current
C****     are computed.
	
	  IF(TYPE.EQ.1)THEN
C****       separately excited motor
            I=(V-KF*2.*PI*N)/RA
            TEM=KF*I

	  ELSEIF(TYPE.EQ.2)THEN
C****       series motor
	    I=V/(2.*PI*N*MAF+RA+RF)
	    TEM=MAF*I**2

	  ENDIF

        ENDIF

        TSHAFT=TEM-TLOSS

C****   check on TSHAFT
        IF(TSHAFT.LT.0.)THEN
          TSHAFT=0.
        ENDIF

C---------------------------------------------------------
C****   Convergence promotion on torque

	X=S(ISTORE)
	GPOS=S(ISTORE+1)
	XPOS=S(ISTORE+2)
	GNEG=S(ISTORE+3)
	XNEG=S(ISTORE+4)
	GLAST=S(ISTORE+5)
	XLAST=S(ISTORE+6)


	IF(INFO(7).EQ.NN)THEN
C****     successive substitution

	  XNEW=TSHAFT
          FLAG=0

	ELSEIF((INFO(7).EQ.NN+1).OR.(INFO(7).EQ.NN+2))THEN
C****     successive substitution

	  XNEW=TSHAFT
	  G=XNEW-X
	  IF(G.GT.0.)THEN
	    GPOS=G
	    XPOS=X
	  ELSE
	    GNEG=G
	    XNEG=X
	  ENDIF 

	ELSE

	  XNEW=TSHAFT
	  G=XNEW-X

C****     Convergence check
      	  
          IF(ABS(G).GT.TOL)THEN

	    IF((ABS(ABS(G)-ABS(GNEG)).LT.ETOL).OR.(ABS(ABS(G)
     &        -ABS(GPOS)).LT.ETOL).OR.(FLAG.EQ.1))THEN

C****         use bisection method
	      IF(G.GT.0.)THEN
	        XPOS=X
	      ELSE
	        XNEG=X
	      ENDIF
	
	      XNEW=(XPOS+XNEG)/2.
	      FLAG=1

	    ELSE
C****         use secant method unless residuals increase
      	      IF(INFO(7).EQ.NN+3)THEN
C****           two most reasonable previous guesses are used
      			
	        IF((G.GT.GNEG).AND.(G.GT.GPOS))THEN	
	          G=GPOS
	          X=XPOS
	          GLAST=GNEG
	          XLAST=XNEG
	        ELSE
	          IF(G.LT.0.)THEN
	            IF(ABS(G).GT.ABS(GNEG))THEN
		      GLAST=GNEG
		      XLAST=XNEG
		    ENDIF
	          ELSE
	            IF(ABS(G).GT.GPOS)THEN
		      GLAST=GPOS
		      XLAST=XPOS
	            ENDIF
	          ENDIF		             	    
 	        ENDIF
      	      ELSE

      		IF(G.GT.GLAST)THEN
C****      	  switch to bisection method
      	          XNEW=(XPOS+XNEG)/2.
      		  FLAG=1
	        ENDIF                

      	      ENDIF		

	      CHANGE=G*(X-XLAST)/(G-GLAST)

	      XNEW=X-CHANGE

              IF((XNEW.GT.XNEG).OR.(XNEW.LT.XPOS))THEN
C****           switch to bisection
      	        XNEW=(XPOS+XNEG)/2.
                FLAG=1
              ENDIF 
	    ENDIF
            TSHAFT=XNEW
          ELSE
C****       solution is found
            TSHAFT=X
          ENDIF
	ENDIF

	XLAST=X
	GLAST=G
	X=XNEW
	S(ISTORE)=X
	S(ISTORE+1)=GPOS
	S(ISTORE+2)=XPOS
	S(ISTORE+3)=GNEG
	S(ISTORE+4)=XNEG
	S(ISTORE+5)=GLAST
	S(ISTORE+6)=XLAST
	S(ISTORE+7)=NN
	S(ISTORE+8)=MEMO
	S(ISTORE+9)=FLAG
	S(ISTORE+10)=EFLAG

      ENDIF			


      IF(N.GT.NMAX)THEN
	FAIL=1
      ENDIF
      IF(I.GT.IRAT)THEN
	IF(I.GT.IMAX)THEN
	  FAIL=3
	ELSE 
	  FAIL=2
	ENDIF
      ENDIF	

      PEL=I*V
      PMECH=2.*PI*N*TSHAFT
      IF(PEL.NE.0.)THEN
	EFF=PMECH/PEL
      ELSE
	EFF=0.
      ENDIF					
      
C**** SET OUTPUTS
      OUT(1)=I
      OUT(2)=V
      OUT(3)=N
      OUT(4)=TSHAFT
      OUT(5)=PEL	
      OUT(6)=PMECH	
      OUT(7)=EFF	
      OUT(8)=FAIL

	RETURN 1
      END