[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[TRNSYS-users] Fwd: new component



Dear TRNSYS users


I tried to modify a trnsys component (type 160) by introducing a new electrochemical model. I have followed the procedure and I recived this message:

 

“Fatal Error at time   :      4001.000000

    Generated by Unit     :     8

    Generated by Type     :   208

    TRNSYS Message      6 : The TRNSYS Program attempted to evaluate a logarithm (base 10) on a value that is less than or equal to zero. Please check the EQUATION formulation and re-run the simulation

    Reported information: Not available”

 

Then I tried to create a new component that have the same program as Type 160 (the original one) and I received the same message.

 

Attached is the program of type206 the modified version of type160.


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

C**********************************************************************************************
C                                                                             ØYSTEIN ULLEBERG 
C                                                                                          IFE 
C                                                                      
C   TYPE160: ADVANCED ALKALINE ELECTROLYZER - Version 4.5              
C                                                                      
C   Versions:                                                          
C   1.0  1998       - Originally developed for ØU's PhD-thesis [2]     
C   2.0  2000-08-07 - Modified for TRNSYS 14.2 W/IISiBat, ØU           
C   3.0  2001-03-14 - Modified for TRNSYS 15.0 w/IISiBat3, ØU          
C   4.0  2001-10-10 - Modifications for TYPE160 TRNSYS15.0 w/IISiBat3, 
C                     modularization of TYPE160, new ELYFARAD equation, 
C                     improved thermal models in ELYOFF and ELYTHERM,  
C                     improved source code readability, ØU             
C   4.1  2002-01-23 - Bugs in FORMAT statements fixed                  
C   4.2  2002-04-14 - NUMSTR removed from declaration                  
C   4.3  2003-10-30 - Checking voltage & current limits are skipped    
C   4.4  2004-08-17 - Updated to TRNSYS 16.0 coding standard, DB
C   4.5  2005-03-21 - fixed a bug. This component resets values in its own XIN() array. I
C                     added an internal array called XIN_INT(), which is filled from the
C                     normal XIN() array at each timestep. Internal modifications are now
C                     made to the XIN_INT() array and not to XIN() itself.       
C   4.6  2006-07-27 - DAA: Added a check so Type160 does not issue a warning if it is turned
C                     off and the voltage is low.
C                                                                      
C                                                                      *
C   MAIN EQUATIONS:                                                    *
C   ---------------                                                    *
C                                                                      *
C   THE MAIN REACTION - DECOMPOSITION OF WATER IS :                    *
C                                                                      *
C       H2O -> H2 + 1/2 O2                                             *
C                                                                      *
C   THE MODEL CONSISTS OF THREE MAIN PARTS:                            *
C                                                                      *
C   1. THE I-U CHARACTERISTICS USED IS OF THE FORM:                    *
C                                                                      *
C       U = Urev(T) + r(T)*I/A + s*log10[(t(T)*I/A+1]                  *
C                                                                      *
C       where:                                                         *
C                                                                      *
C                   DG(T)                                              *
C         Urev(T) = -----              reversible voltage              *
C                   n*F                                                *
C                                                                      *
C         r(T) = r1 + r2*T             resistance in electrolyte       *
C         s    = s1                    overvoltage on electrodes       *
C         t(T) = t1 + t2/T + t3/T^2    overvoltage on electrodes       *
C                                                                      *
C         P = constant                                                 *
C                                                                      *
C                                                                      *
C   2. THERMAL ENERGY BALANCE:                                         *
C                                                                      *
C       Qstore = Qgen - Qloss - Qcool                                  *
C                                                                      *
C       where:                                                         *
C                                                                      *
C         Qstore = Ct*dT/dt                        stored heat         *
C         Qgen   = ncells*U*I*[1-eta_e]            generated heat      *
C         Qloss  = 1/Rt*(T-Troom)                  heat losses         *
C         Qcool  = n_H2O*cp_hH2O*(TH2Oin-TH2Oout)  cooling             *
C                                                                      *
C                                                                      *
C   3. OXYGEN & HYDROGEN PRODUCTION:                                   *
C                                                                      *
C                (I/A)^2                                               *
C      Eta_f = ------------ * a2       Faraday efficiency              *
C              a1 + (I/A)^2                                            *
C                                                                      *
C                               I                                      *
C      n_H2 = Eta_f * ncells * ---     molar H2 production rate        *
C                              n*F                                     *
C                                                                      *
C      n_O2 = 1/2 * n_H2               molar O2 production rate        *
C                                                                      *
C                                                                      *
C   SUBROUTINES                                                        *
C                                                                      *
C   TYPE160 calls the following 6 subroutines:                         *
C                                                                      *
C   1. ELYOFF:                                                         *
C      - Simplified calculations when the electrolyzer is switched OFF *
C                                                                      *
C   2. ELYGIBBS:                                                       *
C      - Thermodynamics in a water splitting reaction (per cell basis) *
C                                                                      *
C   3. ELYELEC:                                                        *
C      - Electrochemical model calculations (per cell basis)           *
C                                                                      *
C   4. ELYFARAD:                                                       *
C      - Faraday efficiency calculations                               *
C                                                                      *
C   5. ELYSTACK                                                        *
C      - Overall stack calculations                                    *
C                                                                      *
C   6. ELYTHERM                                                        *
C      - Thermal energy equation calculaltions (per stack basis)       *
C                                                                      *
C                                                                      *
C   ASSUMPTIONS:                                                       *
C   See reference [1].                                                 *
C                                                                      *
C   REFERENCES:                                                        *
C   1. Ulleberg Ø (2001) Modeling of advanced alkaline electrolyzers:  *
C      A systems Simulation Approach.  Paper to be published in 2001.  *
C   2. Ulleberg Ø (1998) Stand-Alone Power Systems for the Future:     *
C      Optimal Design, Operation & Control of Solar-Hydrogen Energy    *
C      Systems, PhD thesis.  Norwegian University of Science and       *
C      Technology, Trondheim.  ISBN 82-471-0344-3.                     *
C                                                                      *
C                                                                      *
C   PARAMETERS - TRNSYS DECK:                                          *
C   1.  TMODE    [1-3]      TEMPERATURE MODE                           *
C                           MODE 1 = T IS GIVEN AS INPUT XIN(7)        *
C                           MODE 2 = T FROM SIMPLE THERMAL MODEL       *
C                           MODE 3 = T FROM COMPLEX THERMAL MODEL      *
C   2.  AREA     [m^2]      AREA OF ELECTRODES                         *
C   3.  NCELLS   [#]        NUMBER OF CELLS IN SERIES                  *
C   4.  NSTACKS  [#]        NUMBER OF STACKS IN PARALLEL PER UNIT      *
C   5.  IDMAX    [mA/cm^2]  MAXIMUM ALLOWABLE CURRENT DENSITY          *
C   6.  TMAX     [°C]       MAXIMUM ALLOWABLE OPERATING TEMPERATURE    *
C   7.  UCMIN    [V/cell]   MINIMUM ALLOWABLE OPERATING CELL VOLTAGE   *
C   8.  R_T      [K/W]      THERMAL RESISTANCE                         *
C   9.  TAU_T    [hr]       THERMAL TIME CONSTANT                      *
C   10. ELYTYPE  [#]        TYPE OF ELECTOLYZER                        *
C   11. LUD      [#]        LOGICAL UNIT OF PARAMETER FILE             *
C                                                                      *
C   PARAMETERS - EXTERNAL FILE:                                        *
C   1.  R1    [ohm m^2]     OHMIC RESISTANCE                           *
C   2.  R2    [ohm m^2/C]   OHMIC RESISTANCE                           *
C   3.  S1    [V]           OVERVOLTAGE ON ELECTRODES                  *
C   4.  T1    [m^2/A]       OVERVOLTAGE ON ELECTRODES                  *
C   5.  T2    [m^2 °C/A]    OVERVOLTAGE ON ELECTRODES                  *
C   6.  T3    [m^2 °C^2/A]  OVERVOLTAGE ON ELECTRODES                  *
C   7.  A1    [mA/cm2]      FARADAY EFFICIENCY                         *
C   8.  A2    [0...1]       FARADAY EFFICIENCY                         *
C   9.  HX1   [W/°C]        HEAT EXCHANGER COEFFICIENT                 *
C   10. HX2   [W/°C per A]  HEAT EXCHANGER COEFFICIENT                 *
C                                                                      *
C   INPUTS:                                                            *
C   1. SWITCH     [0 or 1]  ELECTROLYZER OPERATING SWITCH              *
C   2. IELY       [A]       CURRENT THROUGH ELECTROLYZER               *
C   3. PELY       [bar]     PRESSURE                                   *
C   4. TROOM      [°C]      AMBIENT (ROOM) TEMPERATURE                 *
C   5. TCOOLIN    [°C]      TEMPERATURE - COOLING WATER INLET          *
C   6. VCOOL      [Nm3/h]   COOLING WATER FLOW RATE                    *
C   7. TELY/TINI  [°C]      TEMERATURE - ELECTROLYZER                  *
C                           (TELY->TCMODE=1, TINI->TCMODE=2 OR 3)      *
C                                                                      *
C   VARIABLES:                                                         *
C   T             [°C]      TEMPERATURE                                *
C   P             [bar]     PRESSURE                                   *
C                                                                      *
C   SUBSCRIPTS:                                                        *
C   COOL       COOLING WATER                                           *
C   ELY        ELECTROLYZER                                            *
C   0,STD      STANDARD CONDITIONS                                     *
C   REF        REFERENCE                                               *
C   TOT        TOTAL/OVERALL                                           *
C   H2         HYDROGEN                                                *
C   O2         OXYGEN                                                  *
C   H2O        WATER                                                   *
C   GAS        IDEAL GAS                                               *
C   FLAG       FOR ITERATION PURPOSES                                  *
C   DIFF       FOR ITERATION PURPOSES                                  *
C   OLD        FOR ITERATION PURPOSES                                  *
C   INI        VALUE AT BEGINNING OF TIMESTEP                          *
C                                                                      *
C   OUTPUTS:                                                           *
C   1.  IELY      [A]       CURRENT DRAWN BY ELECTROLYZER              *
C   2.  UELY      [V]       TOTAL VOLTAGE ACROSS ELECTROLYZER          *
C   3.  PTOT      [W]       TOTAL POWER DRAWN BY ELECTROLYZER          *
C   4.  VDOT_H2   [Nm3/hr]  HYDROGEN PRODUCTION RATE                   *
C   5.  VDOT_O2   [Nm3/hr]  OXYGEN PRODUCTION RATE                     *
C   6.  ETA_TOT   [0...1]   OVERALL EFFICIENCY                         *
C   7.  ETA_E     [0...1]   ENERGY EFFICIENCY                          *
C   8.  ETA_F     [0...1]   FARADAY EFFICIENCY                         *
C   9.  QGEN      [W]       HEAT GENERATED BY ELECTROLYZER             *
C   10. QLOSS     [W]       THERMAL HEAT LOSS                          *
C   11. QCOOL     [W]       AUXILARY COOLING                           *
C   12. QSTORE    [W]       ENERGY STORAGE RATE                        *
C   13. TELY      [°C]      TEMERATURE - ELECTROLYZER                  *
C   14. TCOOLOUT  [°C]      TEMPERATURE - COOLING WATER OUTLET         *
C   15. IDENSITY  [mA/cm2]  CURRENT DENSITY                            *
C   16. UCELL     [V/CELL]  VOLTAGE ACROSS A SINGLE CELL               *
C   17. UREV      [V/CELL]  REVERSIBLE VOLTAGE                         *
C   18. UTN       [V/CELL]  THERMONEUTRAL VOLTAGE                      *
C   19. void                for testing purposes                       *
C   20. void                for testing purposes                       *
C                                                                      *
C**********************************************************************************************
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: TYPE206

C    USE STATEMENTS
      USE TrnsysFunctions
	USE TrnsysConstants

      IMPLICIT NONE !force explicit declaration of local variables

C    TRNSYS DECLARATIONS
	DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,STORED
      INTEGER*4 INFO(15),NPMAX,NI,NO,ND,NS,IUNIT,ITYPE,ICNTRL
      CHARACTER*3 OCHECK,YCHECK

C    SET THE MAXIMUM NUMBER OF PARAMETERS(NP),INPUTS(NIMAX),OUTPUTS(NO),AND DERIVATIVES(ND)
C    THAT MAY BE SUPPLIED FOR THIS TYPE
      PARAMETER (NPMAX=20,NI=7,NO=18,ND=0,NS=2)

C    REQUIRED TRNSYS DIMENSIONS
      DIMENSION XIN(NI),OUT(NO),PAR(NPMAX),YCHECK(NI),OCHECK(NO),
     & STORED(NS)

C    COMMON BLOCKS
      !none required by this Type
	      
C    LOCAL VARIABLE DECLARATIONS
      DOUBLE PRECISION ELYPAR(10)
      INTEGER TMODE,NSTACKS,NELY,ELYTYPE,LUD,I,NP
	CHARACTER (len=60) ELYTEXT
	CHARACTER (len=30) IDENSITYStr,IDMAXStr,UCELLStr,UCMINStr,
	1TELYStr,TMAXStr
	CHARACTER (len=maxMessageLength) T160Msg
	DOUBLE PRECISION IDMAX,TMAX,UCMIN
      DOUBLE PRECISION R1,R2,S1,T1,T2,T3,A1,A2,HX1,HX2
	DOUBLE PRECISION IDENSITY,UCELL,IELY,TELY,TINI,TNEW
	DOUBLE PRECISION BETA,DIF,TIME0,DELT,DIFF,XIN_INT(NI)
      INTEGER*4 IS,TDEBUG1,TDEBUG2,LUW
	LOGICAL DeckPars,LUPars


C    DATA STATEMENTS
	DATA IUNIT/0/
	DATA TDEBUG1/-2/,TDEBUG2/-1/
      DATA YCHECK/'DM1','EL1','PR1','TE1','TE1','SF1','TE1'/
	DATA OCHECK/'EL1','EL2','PW2','SF1','SF1','DM1','DM1','DM1','PW2',
     &            'PW2','PW2','PW2','TE1','TE1','EL3','EL4','EL4','EL4'/ 

C--------------------------------------------------------------------------------------------------
C    GET GLOBAL TRNSYS SIMULATION VARIABLES
      TIME0   = getSimulationStartTime()
      DELT    = getSimulationTimeStep()
	LUW     = getListingFileLogicalUnit()

C--------------------------------------------------------------------------------------------------
C    SET THE VERSION INFORMATION FOR TRNSYS
      IF(INFO(7).EQ.-2) THEN
	  INFO(12)=16
	  RETURN 1
	ENDIF

C--------------------------------------------------------------------------------------------------
C    PERFORM LAST CALL MANIPULATIONS 
      IF (INFO(8).EQ.-1) THEN
	  RETURN 1
	ENDIF

C--------------------------------------------------------------------------------------------------
C    PERFORM POST CONVERGENCE MANIPULATIONS 
      IF(INFO(13).GT.0) THEN
	  CALL getStorageVars(STORED,NS,INFO)
	  STORED(1)=STORED(2)
	  CALL setStorageVars(STORED,NS,INFO)
	  RETURN 1
	ENDIF

C--------------------------------------------------------------------------------------------------
C    PERFORM FIRST CALL MANIPULATIONS
      IF (INFO(7).EQ.-1) THEN
        !retrieve unit and type number for this component from the INFO array
        IUNIT=INFO(1)
	  ITYPE=INFO(2)
        !reserve space in the OUT array using INFO(6)
	  INFO(6)  = NO
	  !set the way in which this component is to be called
        INFO(9)  = 1
	  !set required number of spots in the single precision storage structure    
	  INFO(10) = 0 
	  !reserve required number of spots in the double precision storage structure  
	  CALL setStorageSize(NS,INFO)
        !check for the correct number of PARAMETERS.
	  IF (INFO(4).EQ.11) THEN     !some parameters are contained in the input file, others in an external file.
	    DeckPars = .FALSE.
	    LUPars   = .TRUE.
	    NP=11
	  ELSEIF (INFO(4).EQ.19) THEN !all parameters are contained in the input file.
	    DeckPars = .TRUE.
	    LUPars   = .FALSE.
	    NP=19
	  ELSE 
          CALL TYPECK(4,INFO,0,0,0)
        ENDIF
	  !call the type check routine to compare what this type requires to what is in the input file
    	  CALL TYPECK(1,INFO,NI,NP,ND)
	  !call the input-output check subroutine to set the correct input and output units.
        CALL RCHECK(INFO,YCHECK,OCHECK)
        !RETURN TO THE CALLING PROGRAM
        RETURN 1
      ENDIF
	
C--------------------------------------------------------------------------------------------------
C    PERFORM INITIAL TIMESTEP MANIPULATIONS
      IF (TIME.LT.(TIME0+DELT/2.)) THEN
        !set the UNIT number for future calls
        IUNIT=INFO(1)
	  ITYPE=INFO(2)
        !determine whether parameters have been provided in the deck or in the external file.
	  IF (INFO(4).EQ.11) THEN
	    DeckPars = .FALSE.
	    LUPars   = .TRUE.
	  ELSEIF (INFO(4).EQ.19) THEN
	    DeckPars = .TRUE.
	    LUPars   = .FALSE.
	  ENDIF
        !read the parameter values that are always contained in the deck.
	  !NOTE: not all parameters are required by the main routine.
        TMODE   = JFIX(PAR(1)+0.1d0)
	  NSTACKS = JFIX(PAR(4)+0.1d0)
        IDMAX   = PAR(5)                
        TMAX    = PAR(6)                
        UCMIN   = PAR(7)                
        !read the parameters that may be in the deck and may be in the external file.
        IF (DeckPars.EQ..FALSE.) THEN   !the remaining parameters are in an external file
          ELYTYPE  = JFIX(PAR(10)+0.1)      
          LUD      = JFIX(PAR(11)+0.1)
	    NELY     = 10
          !get the remaining parameters from the external file.
          CALL PARREAD(NELY,LUD,ELYTYPE,ELYPAR,ELYTEXT,INFO)
          R1       = ELYPAR(1)
          R2       = ELYPAR(2)
          S1       = ELYPAR(3)
          T1       = ELYPAR(4)
          T2       = ELYPAR(5)
          T3       = ELYPAR(6)
          A1       = ELYPAR(7)
          A2       = ELYPAR(8)
          HX1      = ELYPAR(9)
          HX2      = ELYPAR(10)
        ELSE                            !the remaining parameters are in the deck
          R1       = PAR(10) ; ELYPAR(1)  = R1
          R2       = PAR(11) ; ELYPAR(2)  = R2
          S1       = PAR(12) ; ELYPAR(3)  = S1
          T1       = PAR(13) ; ELYPAR(4)  = T1
          T2       = PAR(14) ; ELYPAR(5)  = T2
          T3       = PAR(15) ; ELYPAR(6)  = T3
          A1       = PAR(16) ; ELYPAR(7)  = A1
          A2       = PAR(17) ; ELYPAR(8)  = A2
          HX1      = PAR(18) ; ELYPAR(9)  = HX1
          HX2      = PAR(19) ; ELYPAR(10) = HX2
	    ELYTEXT  = 'Electrolyzer description unavailable'
        ENDIF
        !write electrolyzer information to the list file
        WRITE(LUW,110) INFO(1),ELYTEXT,R1,R2,S1,T1,T2,T3,A1,A2,HX1,HX2
        !check the parameters for problems and RETURN if any are found
        	IF ((TMODE.LT.1).OR.(TMODE.GT.3)) CALL TYPECK(-4,INFO,0,1,0)
        !set initial values of the global storage array variables
	  STORED(1:NS) = XIN(7)
	  CALL setStorageVars(STORED,NS,INFO)
	  !set initial output values
	  OUT(1:NO)=0.d0
        RETURN 1
      ENDIF

C--------------------------------------------------------------------
C    THIS IS AN ITERATIVE CALL TO THIS COMPONENT ***
C--------------------------------------------------------------------

C    RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED SINCE THE LAST 
C    TIME THEY WERE READ IN
      IF(INFO(1).NE.IUNIT) THEN
        IUNIT = INFO(1)
        ITYPE = INFO(2)
        !reread parameter values (those contained in TRNSYS input file)
        !determine whether parameters have been provided in the deck or in the external file.
	  IF (INFO(4).EQ.11) THEN
	    DeckPars = .FALSE.
	    LUPars   = .TRUE.
	  ELSEIF (INFO(4).EQ.19) THEN
	    DeckPars = .TRUE.
	    LUPars   = .FALSE.
	  ENDIF
        !read the parameter values that are always contained in the deck.
	  !NOTE: not all parameters are required by the main routine.
        TMODE   = JFIX(PAR(1)+0.1d0)
	  NSTACKS = JFIX(PAR(4)+0.1d0)
        IDMAX   = PAR(5)                
        TMAX    = PAR(6)                
        UCMIN   = 0                
        !read the parameters that may be in the deck and may be in the external file.
        IF (DeckPars.EQ..FALSE.) THEN   !the remaining parameters are in an external file
          ELYTYPE  = JFIX(PAR(10)+0.1)      
          LUD      = JFIX(PAR(11)+0.1)
	    NELY     = 10
          !get the remaining parameters from the external file.
          CALL PARREAD(NELY,LUD,ELYTYPE,ELYPAR,ELYTEXT,INFO)
          R1       = ELYPAR(1)
          R2       = ELYPAR(2)
          S1       = ELYPAR(3)
          T1       = ELYPAR(4)
          T2       = ELYPAR(5)
          T3       = ELYPAR(6)
          A1       = ELYPAR(7)
          A2       = ELYPAR(8)
          HX1      = ELYPAR(9)
          HX2      = ELYPAR(10)
        ELSE                            !the remaining parameters are in the deck
          R1       = PAR(10) ; ELYPAR(1)  = R1
          R2       = PAR(11) ; ELYPAR(2)  = R2
          S1       = PAR(12) ; ELYPAR(3)  = S1
          T1       = PAR(13) ; ELYPAR(4)  = T1
          T2       = PAR(14) ; ELYPAR(5)  = T2
          T3       = PAR(15) ; ELYPAR(6)  = T3
          A1       = PAR(16) ; ELYPAR(7)  = A1
          A2       = PAR(17) ; ELYPAR(8)  = A2
          HX1      = PAR(18) ; ELYPAR(9)  = HX1
          HX2      = PAR(19) ; ELYPAR(10) = HX2
	    ELYTEXT  = 'Electrolyzer description unavailable'
        ENDIF
	ENDIF


C---START OF MAIN PROGRAM----------------------------------------------

C    This component resets values in its own XIN() array for internal calculations. So that the 
C     global values in the XIN() array are not modified, the contents of the XIN() array need to
C     be read into an internal array for those internal calulations.
      XIN_INT(1:NI) = XIN(1:NI)

C    Retrieve stored values for this component
      CALL getStorageVars(STORED,NS,INFO)


C   Get the Values of the Inputs
      IELY   = XIN(2)/NSTACKS		!Current per stack
	XIN_INT(2) = IELY				!Resetting current for internal calculations

      TINI   = STORED(1)			!Initial electrolyzer temperature
	TELY   = TINI				!Temperature (first guess)


C---ELECTROLYZER IS SWITCHED 'OFF'---
      IF ((JFIX(XIN_INT(1)+0.01) .EQ. 0).OR.(IELY.LT.1E-6)) THEN
	  CALL ELYOFF(ELYPAR,PAR,XIN_INT,OUT,DELT,TELY,INFO)

C---ELECTROLYZER IS SWITCHED 'ON'---
	ELSE

C   Thermal steady state operation (TMODE=1)
        IF (TMODE.EQ.1) THEN
	    TINI = XIN_INT(7)									 !T is given
		TELY = TINI										 !constant T
	    CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)		 !Thermodynamics
	    IF (ErrorFound () ) RETURN 1
	    CALL ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)		 !Electrochemistry
	    IF (ErrorFound () ) RETURN 1
		CALL ELYFARAD(ELYPAR,PAR,XIN_INT,OUT)				 !Efficiency calcs.
  	    CALL ELYSTACK(ELYPAR,PAR,XIN_INT,OUT)				 !Stack calculations
	    CALL ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELY) !Thermal energy calcs

C   Thermal non-steady state operation (TMODE=2 or 3)
	  ELSEIF ((TMODE.EQ.2).or.(TMODE.EQ.3)) THEN
		BETA = 0.6
		DIFF = 10.
		I = 0				
		DO WHILE ((DIFF.gt.0.0001).and.(I.lt.100))
		  CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)		!Thermodynamics
		  IF (ErrorFound () ) RETURN 1		
		  CALL ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)		!Electrochemistry
		  IF (ErrorFound () ) RETURN 1		
		  CALL ELYFARAD(ELYPAR,PAR,XIN_INT,OUT)					!Efficiency			
  		  CALL ELYSTACK(ELYPAR,PAR,XIN_INT,OUT)					!Stack calcs			
		  CALL ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELY)	!Thermal energy
		  TNEW = OUT(13)									!T from ELYTHERM				
		  DIFF = DABS(TELY-TNEW)							!T difference
		  TELY = BETA*TNEW+(1-BETA)*TELY					!Resetting T
		  I = I+1	
		END DO
	  ENDIF
	ENDIF
   

C   Warnings
	IDENSITY = OUT(15)
	UCELL = OUT(16)
      
	IF (IDENSITY.GT.IDMAX) THEN 
	  WRITE(IDENSITYStr,*) JFIX(IDENSITY+0.1d0)
	  WRITE(IDMAXStr,*) JFIX(IDMAX+0.1d0)
        T160Msg = 'The electrolyzer current ('//TRIM(ADJUSTL(IDENSITYStr
	1))//' [mA/cm^2]) is too high. The maximum allowable electrolyzer c
     1urrent is '//TRIM(ADJUSTL(IDMAXStr))//' [mA/cm^2].'
	  CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
      ENDIF 

      IF ((UCELL.LT.UCMIN).and.
     &	(JFIX(XIN_INT(1)+0.01).NE.0)) THEN 
	  WRITE(UCELLStr,*) UCELL
	  WRITE(UCMINStr,*) UCMIN
        T160Msg = 'The electrolyzer voltage ('//TRIM(ADJUSTL(UCELLStr))
	1//' [V]) is too low. The minimum allowable electrolyzer voltage i
     1s '//TRIM(ADJUSTL(UCMINStr)) //' [V].'
	  CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
      ENDIF 

      IF (TELY.GT.TMAX) THEN 
	  WRITE(TELYStr,*) TELY
	  WRITE(TMAXStr,*) TMAX
        T160Msg = 'The electrolyzer temperature ('//TRIM(ADJUSTL(TELYStr
	1))//' [C]) is too high. The maximum allowable electrolyzer tempera
     1ture is '//TRIM(ADJUSTL(TMAXStr)) //' [C].'
	  CALL MESSAGES(-1,T160Msg,'warning',IUNIT,ITYPE)
      ENDIF 


C    WRITE DEBUGGING INFORMATION 
	IF((TIME.GE.TDEBUG1).AND.(TIME.LE.TDEBUG2)) THEN
	  WRITE(LUW,*) 'UNIT        = ',INFO(1)
	  WRITE(LUW,*) 'TYPE        = ',INFO(2)
	  WRITE(LUW,*) 'TIME        = ',TIME 
	  WRITE(LUW,*) 'ITER        = ',INFO(7) 
	  WRITE(LUW,*) 'PELYTOT     = ',OUT(3)*NSTACKS	
	  WRITE(LUW,*) ' '	
	ENDIF

C    STORE THE FINAL ELECTROLYZER TEMPERATURE AS THE NEXT TIMESTEP'S INITIAL TEMPERATURE.
      STORED(2) = OUT(13)	!TELY from ELYTHERM
	CALL setStorageVars(STORED,NS,INFO)

C---OUTPUTS---
	OUT(1)  = IELY* NSTACKS	! IELY
	OUT(2)  = OUT(2)				! USTACK
	OUT(3)  = OUT(3)* NSTACKS	! PTOT
	OUT(4)  = OUT(4)* NSTACKS	! VH2
	OUT(5)  = OUT(5)* NSTACKS	! VO2
	OUT(6)  = OUT(6)			! ETA_TOT
	OUT(7)  = OUT(7)			! ETA_E
	OUT(8)  = OUT(8)			! ETA_F
	OUT(9)  = OUT(9)* NSTACKS	! QGEN
	OUT(10) = OUT(10)* NSTACKS	! QLOSS
	OUT(11) = OUT(11)* NSTACKS	! QCOOL
	OUT(12) = OUT(12)* NSTACKS 	! QSTORE
	OUT(13) = OUT(13)			! TELY; final temperature
	OUT(14) = OUT(14)			! TCOOLOUT 
	OUT(15) = OUT(15)  			! IDENSITY
	OUT(16) = OUT(16)			! UCELL
	OUT(17) = OUT(17)			! UREV
	OUT(18) = OUT(18)			! UTN


110   FORMAT(//2X,72H***************************************************
     .*********************,//30X,13HS U M M A R Y,//24X,23HELECTROLYZER 
     . PARAMETERS,7H - UNIT,I3,//4X,A50,//4X,5HR1 = ,E10.4,3X,9H[ohm m^2
     .],/4X,5HR2 =  ,E10.4,3X,11H[ohm m^2/C],/4X,5HS1 = ,E10.4,3X,3H[V],
     ./4X,5HT1 = ,E10.4,3X,7H[m^2/A],/4X,5HT2 = ,E10.4,3X,10H[m^2 °C/A],
     ./4X,5HT3 = ,E10.4,3X,12H[m^2 °C^2/A],/4X,5HA1 = ,F12.4,1X,3H[-],/4
     .X,5HA2 = ,F12.4,1X,7H[0...1],/4X,6HHX1 = ,F11.4,1X,6H[W/°C],/4X,6H
     .HX2 = ,F11.4,1X,12H[W/°C per A],//2X,72H**************************
     .**********************************************,/)

150   FORMAT(/16X,25H*********WARNING*********,/16X,4HUNIT,I3,2X,6HTIME
     .=, F9.3,/16X,25HINCOMPLETE CONVERGENCE!!!,/16X,25H****************
     .*********,/)

      RETURN 1
      END

C---END OF ELECTROLYZER SUBROUTINE--------------------------------------



C***********************************************************************
C    EVERYTHING BELOW THIS LINE BELONGS TO THE ELECTROLYZER MODEL      *
C***********************************************************************
 	
	

	SUBROUTINE ELYOFF(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,INFO)
C***********************************************************************
C                                                                      *
C   ELECTROLYZER SWITCHED OFF                                          *
C                                                                      *
C   NOTE! THERMAL CALCULATIONS ON A PER STACK BASIS                    *
C                                                                      *
C   PARAMETERS:                                                        *
C   TMODE      TEMPERATURE MODE (1=T GIVEN, 2 & 3 = T CALCULATED)      *
C   NCELLS     NUMBER OF CELLS IN SERIES [#]                           *
C   R_T        THERMAL RESISTANCE [K/W]                                *
C   TAU_T      THERMAL TIME CONSTANT [hr]                              *
C                                                                      *
C   INPUTS:                                                            *
C   TROOM      AMBIENT (ROOM) TEMPERATURE [°C]                         *
C   DELT       TIME STEP [hr]                                          *
C   TINI       INITIAL ELECTROLYZER TEMPERATURE [°C]                   *
C                                                                      *
C   VARIABLES:                                                         *
C   C_T        THERMAL CAPACITY OF ELECTROLYZER [J/K]                  *
C   UREST      RESTING VOLTAGE [V]                                     *
C   TINI       INITIAL ELECTROLYZER TEMPERATURE [°C]                   *
C                                                                      *
C   OUTPUTS:                                                           *
C   1.  IELY       CURRENT DRAWN BY ELECTROLYZER [A]                   *
C   2.  UELY       TOTAL VOLTAGE ACROSS ELECTROLYZER [V]               *
C   3.  PTOT       TOTAL POWER DRAWN BY ELECTROLYZER [W]               *
C   4.  VDOT_H2    HYDROGEN PRODUCTION RATE [Nm3/hr]                   *
C   5.  VDOT_O2    OXYGEN PRODUCTION RATE [Nm3/hr]                     *
C   6.  ETA_TOT    OVERALL EFFICIENCY [-]                              *
C   7.  ETA_E      ENERGY EFFICIENCY [-]                               *
C   8.  ETA_F      FARADAY EFFICIENCY [-]                              *
C   9.  QGEN       HEAT GENERATED BY ELECTROLYZER [W]                  *
C   10. QLOSS      THERMAL HEAT LOSS [W]                               *
C   11. QCOOL      AUXILARY COOLING [W]                                *
C   12. QSTORE     ENERGY STORAGE RATE [W]                             *
C   13. TELY       ELECTROLYZER TEMPERATURE [°C]                       *
C   14. TCOOLOUT   TEMPERATURE - COOLING WATER OUTLET [°C]             *
C   15. IDENSITY   CURRENT DENSITY [mA/cm2]                            *
C   16. UCELL      VOLTAGE ACROSS A SINGLE CELL [V/CELL]               *
C   17. UREV       REVERSIBLE VOLTAGE [V/CELL]                         *
C   18. UTN        THERMONEUTRAL VOLTAGE [V/CELL]                      *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
	DOUBLE PRECISION NCELLS,R_T,TAU_T,C_T,AA,BB
	DOUBLE PRECISION QGEN,QCOOL,QLOSS,QSTORE,UREV,UTN,UREST
	DOUBLE PRECISION TROOM,TINI,DELT,TELY
	DOUBLE PRECISION BETA,DIFF,TOLD,TNEW
	INTEGER TMODE,I,INFO(15)


C   PARAMETERS
      TMODE  = JFIX(PAR(1)+0.1)      
	NCELLS = PAR(3)
      R_T    = PAR(8)                
      TAU_T  = PAR(9)           

C   INPUTS
      TROOM   = XIN_INT(4)

C   THERMAL MODEL
C   (On a per stack basis)
      C_T   = TAU_T*3600/R_T			!3600 s = 1 hr
      QGEN  = 0.
      QCOOL = 0.

C   TEMPERATURE (FINAL)
	IF (TMODE.EQ.1) THEN
        TELY = TINI
	  QLOSS = 1/R_T*(TELY-TROOM)
      ELSEIF (TMODE.EQ.2) THEN
        TELY = TINI
	  BETA = 0.6
	  DIFF = 10.
	  I = 0				
	  DO WHILE ((DIFF.gt.0.0001).and.(I.lt.100))
	    TOLD = TELY										!intermediate T
	    QLOSS = 1/R_T*(TELY-TROOM)
          TELY = TINI+(DELT*3600)/C_T*(QGEN-QCOOL-QLOSS)	!1hr=3600s
		TNEW = TELY																
		DIFF = DABS(TNEW-TOLD)							!T difference
		TELY = BETA*TNEW+(1-BETA)*TELY					!Resetting T
		I = I+1	
	  END DO
      ELSEIF (TMODE.EQ.3) THEN
        AA = 1/C_T*(1/R_T)
        BB = 1/C_T*(QGEN+TROOM/R_T)
        TELY = (TINI-BB/AA)*DEXP(-AA*DELT*3600)+BB/AA      !1hr=3600s
	  QLOSS = 1/R_T*(TELY-TROOM)
      ENDIF

      QSTORE = QGEN-QCOOL-QLOSS

C   RESTING VOLTAGE
	CALL ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO)
	UREV  = OUT(17)
	UTN   = OUT(18)
	UREST = UREV

C   OUTPUTS
      OUT(1)  = 0.				!IELY
      OUT(2)  = UREST*NCELLS		!UELY
      OUT(3)  = 0.				!PTOT
      OUT(4)  = 0.				!VDOT_H2
      OUT(5)  = 0.				!VDOT_O2
	OUT(6)  = 0.				!ETA_TOT
      OUT(7)  = 0.				!ETA_E
      OUT(8)  = 0.				!ETA_F
	OUT(9)  = QGEN				
	OUT(10) = QLOSS				
	OUT(11) = QCOOL				
	OUT(12) = QSTORE			
	OUT(13) = TELY				!Final temperature				
	OUT(14) = XIN_INT(5)			!TCOOLOUT=TCOOLIN
	OUT(15) = 0.				!IDENSITY
      OUT(16) = UREST				!UCELL 
      OUT(17) = UREV				
      OUT(18) = UTN				

	RETURN
	END



                
	SUBROUTINE ELYGIBBS(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) 
C***********************************************************************
C                                                                      *
C   THIS SUBROUTINE CALCULATES THE THERMODYNAMICS IN A                 *
C   WATER SPLITTING REACTON (WATER ELECTROLYSIS):                      *
C                                                                      *
C                H2O --> 1/2 O2 + H2                                   *
C                                                                      *
C                                                                      *
C   DATA:                                                              *
C   TREF       REFERENCE TEMPERATURE [C]                               *
C   PREF       REFERENCE PRESSURE [bar]                                *
C   FARAD      FARADAY CONSTANT [C/mol] or [As/mol]                    *
C   RGAS       UNIVERSAL GAS CONSTANT [J/K-mol]                        *
C   NELEC      NUMBER OF ELECTRONS PER MOLECULE OF H2 [#]              *
C                                                                      *
C   INPUTS:                                                            *
C   TELY       TEMPERATURE  [°C]                                       *
C   PELY       PRESSURE [bar]                                          *
C                                                                      *
C   VARIABLES:                                                         *
C   DH         ENTHALPY CHANGE FOR DECOMPOSITION OF WATER [J/mol]      *
C   DS         ENTROPY CHANGE FOR DECOMPOSITION OF WATER [J/mol-K]     *
C   DG         CHANGE IN GIBBS ENERGY - DECOMPOSITION OF WATER [J/mol] *
C                                                                      *
C   OUTPUTS:                                                           *
C   UTN	     THERMONEUTRAL VOLTAGE [V/cell ]                         *
C   UREV       THERMODYNAMIC REVERSIBLE VOLTAGE [V/cell]               *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
	DOUBLE PRECISION TELY,PELY,UTN,UREV
	DOUBLE PRECISION TREF,PREF,FARAD,RGAS,T_H2,T_O2,T_H2O,DH_H2,
	1DH_O2,DH_H2O
	DOUBLE PRECISION DH0F_H2O,DH0F_H2,DH0F_O2,S0F_H2O,S0F_H2,S0F_O2
	DOUBLE PRECISION S_H2,S_O2,S_H2O,CP0_H2O,CP0_H2,CP0_O2
	DOUBLE PRECISION DH,DS,DG
	INTEGER NELEC,INFO(15)
      DATA TREF/25.0/		!C
      DATA PREF/1.0/		!bar
      DATA FARAD/96485/   !As/mol
      DATA RGAS/8.3145/	!J/K-mol
	DATA NELEC/2/		!mol electrons/mol H2
      DATA DH0F_H2O/-286E3/,  DH0F_H2/0/,     DH0F_O2/0/       !J/mol
      DATA S0F_H2O/70/,       S0F_H2/131/,    S0F_O2/205/      !J/K-mol
      DATA CP0_H2O/75/,       CP0_H2/29/,     CP0_O2/29/       !J/K-mol

C   INPUTS
	PELY = XIN_INT(3)

C   TEMPERATURE
	T_H2  = TELY
	T_O2  = TELY
	T_H2O = TELY

C   ENTHALPY
      DH_H2  = CP0_H2*(T_H2-TREF)+DH0F_H2
      DH_O2  = CP0_O2*(T_O2-TREF)+DH0F_O2
      DH_H2O = CP0_H2O*(T_H2O-TREF)+DH0F_H2O
      DH = DH_H2+0.5*DH_O2-DH_H2O

C   ENTROPY - CHECK LOGARITHM ARGUMENTS BEFORE TRYING THEM.
      IF (((T_H2+273.15)/(TREF+273.15)).LE.0) THEN
	  CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF

      IF ((PELY/PREF).LE.0) THEN
	  CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF

      S_H2 = CP0_H2*DLOG((T_H2+273.15)/(TREF+273.15))
     &     -RGAS*DLOG(PELY/PREF)+S0F_H2
 

      IF (((T_O2+273.15)/(TREF+273.15)).LE.0) THEN
	  CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF
 
      IF ((PELY/PREF).LE.0) THEN
	  CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF

      S_O2 = CP0_O2*DLOG((T_O2+273.15)/(TREF+273.15))
     &     -RGAS*DLOG(PELY/PREF)+S0F_O2


      IF (((T_H2O+273.15)/(TREF+273.15)).LE.0) THEN
	  CALL MESSAGES(8,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF

      S_H2O = CP0_H2O*DLOG((T_H2O+273.15)/(TREF+273.15))+S0F_H2O
      DS = S_H2+0.5*S_O2-S_H2O

C   GIBBS FREE ENERGY   
      DG = DH-(TELY+273.15)*DS
 
C   THERMONEUTRAL VOLTAGE (PER CELL)
      UTN = DH/(NELEC*FARAD)
 
C   REVERSIBLE VOLTAGE (PER CELL)
      UREV = DG/(NELEC*FARAD)

C   OUTPUTS
	OUT(17) = UREV
	OUT(18) = UTN
     
	RETURN
      END




 	SUBROUTINE ELYELEC(ELYPAR,PAR,XIN_INT,OUT,TELY,INFO) 
C***********************************************************************
C                                                                      *
C   ELECTROCHEMICAL MODEL                                              *
C                                                                      *
C   PARAMETERS:                                                        *
C   AREA  AREA OF ELECTRODES                          [m^2]            *
C   R1    OHMIC RESISTANCE							[ohm-m^2]        *
C   R2           "									[ohm-m^2/C]      *
C   R3           "									[ohm-m^2/C^2]    *
C   S1    PARAMETER FOR OVERVOLTAGE ON ELECTRODES		[V]              *
C   S2           "               "					[V/C]            *
C   S3           "               "					[V/C^2]          *
C   T1    PARAMETER FOR OVERVOLTAGE ON ELECTRODES		[m^2/A]          *
C   T2             "        "							[m^2/A-C]        *
C   T3             "        "							[m^2/A-C^2]      *
C                                                                      *
C   INPUTS:                                                            *
C   TELY       TEMPERATURE  [C]                                        *
C   IELY       CURRENT DRAWN BY ELECTROLYZER [A]                       *
C   UREV       THERMODYNAMIC REVERSIBLE VOLTAGE [V/cell]               *
C                                                                      *
C   OUTPUTS:                                                           *
C   IDENSITY   CURRENT DENSITY [mA/cm2]                                *
C   UCELL      VOLTAGE ACROSS A SINGLE CELL [V/CELL]                   *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20),TELY
	DOUBLE PRECISION IELY,UREV
	DOUBLE PRECISION R1,R2,R3,S1,S2,S3,T1,T2,T3,AREA
	DOUBLE PRECISION UOHMIC,UOVERPOT,UCELL,IDENSITY
	INTEGER INFO(15)
 
C   PARAMETERS
	R1     = ELYPAR(1)
	R2     = ELYPAR(2)
	S1     = ELYPAR(3)
	T1     = ELYPAR(4)
	T2     = ELYPAR(5)
	T3     = ELYPAR(6)
	AREA   = PAR(2)

C   INPUTS
	IELY  = XIN_INT(2)
	UREV  = OUT(17)

C   CURRENT-VOLTAGE CHARACTERISTIC (PER CELL)
      IF (((T1+T2/TELY+T3/(TELY**2))*IELY/AREA+1).LE.0) THEN
	  CALL MESSAGES(6,'','fatal',INFO(1),INFO(2))
	  RETURN
	ENDIF

	IF (IELY.EQ.0) THEN 

	 UOHMIC = (R1+R2*TELY)*IELY/AREA
      UOVERPOT = S1*DLOG10((T1+T2/TELY+T3/(TELY**2))*IELY/AREA+1)
      UCELL = UREV+UOHMIC+UOVERPOT
      ELSE
      UOHMIC = 8.3145*TELY*LOG((IELY+0.1)/0.1)/(2*0.9*96484)
	T1= exp(1267*(1/303-1/TELY))
	UOHMIC = UOHMIC + IELY*0.2/(AREA*10000*(0.005139*14+0.00326)*T1)
      UOVERPOT = (8.3145*TELY/(2*96484))*LOG(1+IELY)
      UCELL = UREV+UOHMIC+UOVERPOT
	ENDIF

C   CURRENT DENSITY
	IDENSITY = IELY/AREA/10		!10 mA/cm2 = 1 A/m2

C   OUTPUTS
	OUT(15) = IDENSITY
	OUT(16) = UCELL

	RETURN
	END



 	SUBROUTINE ELYFARAD(ELYPAR,PAR,XIN_INT,OUT) 
C***********************************************************************
C                                                                      *
C   THIS SUBROUTINE CALCULATES THE FARADAY EFFICIENCY                  *
C                                                                      *
C   PARAMETERS:                                                        *
C   A1			 PARAMETER FOR FARADAY EFFICIENCY  [mA/cm2]    *
C                                                                      *
C   INPUTS:                                                            *
C   IDENSITY       CURRENT DENSITY [mA/cm2]                            *
C                                                                      *
C   OUTPUTS:                                                           *
C   ETA_F          FARADAY EFFICIENCY [0-1]                            *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
	DOUBLE PRECISION A1,A2,IDENSITY,ETA_F

C   PARAMETERS
	A1 = ELYPAR(7)
	A2 = ELYPAR(8) 

C   INPUTS
	IDENSITY = OUT(15)

C   FARADAY EFFICIENCY
	ETA_F = IDENSITY**2/(A1+IDENSITY**2)*A2 

C   OUTPUTS
	OUT(8) = ETA_F
	
      RETURN
      END




 	SUBROUTINE ELYSTACK(ELYPAR,PAR,XIN_INT,OUT) 
C***********************************************************************
C                                                                      *
C   STACK EFFICIENCY, FLOW RATES, VOLTAGE, AND POWER                   *
C                                                                      *
C   DATA:                                                              *
C   TREF       REFERENCE TEMPERATURE [C]                               *
C   PREF       REFERENCE PRESSURE [bar]                                *
C   FARAD      FARADAY CONSTANT [C/mol] or [As/mol]                    *
C   RGAS       UNIVERSAL GAS CONSTANT [J/K-mol]                        *
C   NELEC      NUMBER OF ELECTRONS PER MOLECULE OF H2 [#]              *
C                                                                      *
C   PARAMETERS:                                                        *
C   NCELLS     NUMBER OF CELLS IN SERIES [#]                           *
C                                                                      *
C   INPUTS:                                                            *
C   IELY       CURRENT THROUGH ELECTROLYZER [A]                        *
C   ETA_F      FARADAY EFFICIENCY [0-1]                                *
C   UCELL      VOLTAGE ACROSS A SINGLE CELL [V/CELL]                   *
C   UTN	     THERMONEUTRAL VOLTAGE [V/cell ]                         *
C                                                                      *
C   VARIABLES:                                                         *
C   RHO_STD    DENSITY OF GAS AT STANDARD CONDITIONS [mol/m3]          *
C   NDOT_H2    HYDROGEN PRODUCTION RATE [mol/s]                        *
C   NDOT_O2    OXYGEN PRODUCTION RATE [mol/s]                          *
C                                                                      *
C   OUTPUTS:                                                           *
C   UELY       TOTAL VOLTAGE ACROSS ELECTROLYZER [V]                   *
C   PTOT       TOTAL POWER DRAWN BY ELECTROLYZER [W]                   *
C   VDOT_H2    HYDROGEN PRODUCTION RATE [Nm3/hr]                       *
C   VDOT_O2    OXYGEN PRODUCTION RATE [Nm3/hr]                         *
C   ETA_TOT    OVERALL EFFICIENCY [-]                                  *
C   ETA_E      ENERGY EFFICIENCY [-]                                   *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
	DOUBLE PRECISION IELY,UTN,UCELL,ETA_F
	DOUBLE PRECISION FARAD,TSTD,PSTD,RGAS
	DOUBLE PRECISION RHO_STD,NDOT_H2,NDOT_O2
	DOUBLE PRECISION ETA_E,ETA_TOT,VDOT_H2,VDOT_O2,UELY,PTOT
	INTEGER NELEC,NCELLS
	DATA NELEC/2/		!mol electrons/mol H2
	DATA FARAD/96485/   !As/mol
      DATA RGAS/8.3145/	!J/K-mol
	DATA TSTD/0/		!C
      DATA PSTD/1.0/		!bar

C   PARAMETERS
	NCELLS = JFIX(PAR(3)+0.1d0)

C   INPUTS
	IELY  = XIN_INT(2)
	ETA_F = OUT(8)
	UCELL = OUT(16)
	UTN   = OUT(18)

C   ENERGY EFFICIENCY
      ETA_E = UTN/UCELL

C   OVERALL EFFICIENCY
      ETA_TOT = ETA_E*ETA_F

C   HYDROGEN & OXYGEN PRODUCTION---
      NDOT_H2 = ETA_F*NCELLS*IELY/(NELEC*FARAD)
      NDOT_O2 = 0.5*NDOT_H2

C   PRODUCTION FLOW RATES (STANDARD CUBIC METERS)
      RHO_STD = PSTD*1E5/(RGAS*(TSTD+273.15))		!1 bar = 1E5 Pa
      VDOT_H2 = NDOT_H2/RHO_STD*3600				!Nm3/hr
      VDOT_O2 = NDOT_O2/RHO_STD*3600				!Nm3/hr

C   STACK VOLTAGE
      UELY = NCELLS*UCELL

C   STACK POWER
	PTOT = UELY*IELY

C   OUTPUTS
	OUT(2)  = UELY
	OUT(3)  = PTOT
      OUT(4)  = VDOT_H2
      OUT(5)  = VDOT_O2
	OUT(6)  = ETA_TOT
	OUT(7)  = ETA_E
	
      RETURN
      END



 	SUBROUTINE ELYTHERM(ELYPAR,PAR,XIN_INT,OUT,DELT,TINI,TELYINI) 
C***********************************************************************
C                                                                      *
C   THERMAL ENERGY MODEL (ON A PER STACK BASIS)                        *
C                                                                      *
C   DATA:                                                              *
C   TREF       REFERENCE TEMPERATURE [C]                               *
C   PREF       REFERENCE PRESSURE [bar]                                *
C   FARAD      FARADAY CONSTANT [C/mol] or [As/mol]                    *
C   RGAS       UNIVERSAL GAS CONSTANT [J/K·mol]                        *
C   NELEC      NUMBER OF ELECTRONS PER MOLECULE OF H2 [#]              *
C                                                                      *
C   PARAMETERS:                                                        *
C   TMODE      TEMPERATURE MODE (1=T GIVEN, 2 & 3 = T CALCULATED)      *
C   NCELLS     NUMBER OF CELLS IN SERIES [#]                           *
C   R_T        THERMAL RESISTANCE [K/W]                                *
C   TAU_T      THERMAL TIME CONSTANT [hr]                              *
C   HX1        PARAMETER FOR HEAT EXCHANGER COEFFICIENT [W/°C]         *
C   HX2        PARAMETER FOR HEAT EXCHANGER COEFFICIENT [W/°C per A]   *
C                                                                      *
C   INPUTS:                                                            *
C   IELY       CURRENT THROUGH ELECTROLYZER [A]                        *
C   TROOM      AMBIENT (ROOM) TEMPERATURE [C]                          *
C   TCOOLIN    TEMPERATURE - COOLING WATER INLET [C]                   *
C   VCOOL      COOLING WATER FLOW RATE [Nm3/h]                         *
C   UCELL      VOLTAGE ACROSS A SINGLE CELL [V/CELL]                   *
C   UTN	     THERMONEUTRAL VOLTAGE [V/cell ]                         *
C   DELT       TIME STEP [hr]                                          *
C   TINI       INITIAL ELECTROLYZER TEMPERATURE [°C]                   *
C   TELYINI    ELECTROLYZER TEMERATURE - AT ITERATION N [°C]           *
C                                                                      *
C   VARIABLES:                                                         *
C   UA_HX      OVERALL HEAT TRANSFER COEFFICIENT f(IELY) [W/°C]        *
C   CP_H2O     HEAT CAPACITY OF COOLING WATER [W/K]                    *
C   C_T        THERMAL CAPACITY OF ELECTROLYZER [J/K]                  *
C                                                                      *
C   OUTPUTS:                                                           *
C   QGEN       HEAT GENERATED BY ELECTROLYZER [W]                      *
C   QLOSS      THERMAL HEAT LOSS [W]                                   *
C   QCOOL      AUXILARY COOLING [W]                                    *
C   QSTORE     ENERGY STORAGE RATE [W]                                 *
C   TELYFIN    ELECTROLYZER TEMERATURE - AT ITERATION N+1  [°C]        *
C   TCOOLOUT   TEMPERATURE - COOLING WATER OUTLET [°C]                 *
C                                                                      *
C***********************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION ELYPAR(10),PAR(11),XIN_INT(7),OUT(20)
	DOUBLE PRECISION R_T,TAU_T,HX1,HX2,RHO_H2O,M_H2O,CP0_H2O
	DOUBLE PRECISION IELY,TROOM,VCOOL,UCELL,UTN
	DOUBLE PRECISION QGEN,QLOSS,QCOOL,QSTORE,UA_HX,CP_H2O,C_T,AA,BB
	DOUBLE PRECISION TCOOLIN,TCOOLOUT,DELT,TINI,TELYINI,TELYFIN,TELY
	INTEGER TMODE,NCELLS
      DATA M_H2O/18.016/		!g/mol
      DATA RHO_H2O/1000/		!kg/m3
      DATA CP0_H2O/75/		!J/K-mol


C   PARAMETERS
      TMODE  = JFIX(PAR(1)+0.1)      
	NCELLS = PAR(3)
      R_T    = PAR(8)                
      TAU_T  = PAR(9)
	HX1    = ELYPAR(9)
	HX2    = ELYPAR(10)           

C   INPUTS
	IELY    = XIN_INT(2)
	TROOM   = XIN_INT(4)
      TCOOLIN = XIN_INT(5)
      VCOOL   = XIN_INT(6)
	UCELL   = OUT(16)
	UTN     = OUT(18)

C   TEMPERATURE (INITIAL)
	TELY = TELYINI

C   GENERATED THERMAL ENERGY
      QGEN  = NCELLS*IELY*(UCELL-UTN)

C   HEAT LOSSES TO AMBIENT
	QLOSS = 1/R_T*(TELY-TROOM)

C   AUXILIARY COOLING DEMAND
	CP_H2O = VCOOL*RHO_H2O/M_H2O*CP0_H2O*1000/3600    !1kg=1000g, 1hr=3600s  
	IF (TMODE.EQ.1) THEN
	  QCOOL = QGEN-QLOSS
	  TCOOLOUT = TCOOLIN+QCOOL/CP_H2O
	ELSE 
	  UA_HX = HX1+HX2*IELY									
        TCOOLOUT = TCOOLIN+(TELY-TCOOLIN)*(1-DEXP(-UA_HX/CP_H2O))             
        QCOOL = CP_H2O*(TCOOLOUT-TCOOLIN)
	ENDIF

C   THERMAL CAPACITANCE
	C_T   = TAU_T*3600/R_T		!3600 s = 1 hr

C   TEMPERATURE (FINAL)
      IF (TMODE.EQ.1) THEN
	  TELYFIN = TINI
      ELSEIF (TMODE.EQ.2) THEN
        TELYFIN = TINI+(DELT*3600)/C_T*(QGEN-QCOOL-QLOSS)	!1hr=3600s
      ELSEIF (TMODE.EQ.3) THEN
        AA = 1/C_T*(1/R_T+CP_H2O*(1-DEXP(-UA_HX/CP_H2O)))
        BB = 1/C_T*(QGEN+TROOM/R_T
     &      +CP_H2O*(1-DEXP(-UA_HX/CP_H2O))*TCOOLIN)
        TELYFIN = (TINI-BB/AA)*DEXP(-AA*DELT*3600)+BB/AA		 !1hr=3600s
      ENDIF

C   TOTAL THERMAL ENERGY BALANCE
      QSTORE = QGEN-QCOOL-QLOSS

C   OUTPUTS
      OUT(9)  = QGEN          
      OUT(10) = QLOSS         
      OUT(11) = QCOOL         
      OUT(12) = QSTORE        
      OUT(13) = TELYFIN
      OUT(14) = TCOOLOUT

	RETURN
	END