!  This component has been assigned Type Number 206. If that number conflicts with 
!  another user Type number, you will need to change it and recompile the appropriate
!  dll.

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

!DEC$ATTRIBUTES DLLEXPORT :: TYPE206


C
C***********************************************************************
C    4PORT store Model                                                 *
C                                                                      *
C    This subroutine models a hot water store with internal or mantle  *
C    heat exchanger.                                                   *
C                                                                      *
C    The Model was developed by R.H. Marshall and C.L.W. Li, University*
C    of Wales College of Cardiff in the framework of the Solar Storage *
C    Testing Group between 1986 and 1989. The model was extensively    *
C    validated against experimental data by the five participating     *
C    institutes.                                                       *
C                                                                      *
C    The TYPE74 TRNSYS-module was adapted by R. Kubler and F. Muller,  *
C    ITW, University of Stuttgart with support from the authors of the *
C    model in 1990                                                     *
C                                                                      *
C    Version 1.1                                                       *
C    26. sept.1991                                                     *
C                                                                      *
C    MODIFIED BY N. BLAIR 8-13-92 : NODEMX ADDED TO THE CALL           *
C    STATEMENTS FOR MIX						       *
C                                                                      *
C***********************************************************************
C
C     SAMPLE INPUTS AND OUTPUTS
C
C     PARAMETERS     DESCRIPTION		            SAMPLE VALUES
C     PAR(1) VS      Storage volume(m3)             	         6.28
C     PAR(2) Cps     Storage material specific heat cap.          .98
C     PAR(3) RHOs    Storage material density                  102.4
C     PAR(4) (UA)ss  Storage material eff. diffusivity          20.0
C     PAR(5) Vf      Heat exchanger volume                       0.3
C     PAR(6) Cpf     Heat transfer fluid specific heat cap.      0.98
C     PAR(7) RHOf    Heat transfer fluid density                99.2
C     PAR(8) Hs      Storage height                              2.0
C     PAR(9) Lh      Length of auxiliary heater                   .4
C                     if installed in the top of storage
C     PAR(10) HTRMODE   0 : if none                              2.0
C                       1 : if installed from the top
C                       2 : if installed from the side 
C     PAR(11) NF     Number of heat exchanger nodes              3.0
C     PAR(12) NS     Number of storage nodes                     8.0
C     PAR(13) NDEAD  Number of nodes in the dead zone            1.0
C     PAR(14) NH     number of the node, in which the aux.       4.0
C                    heater is installed (only of HTRMODE=SIDE)
C     PAR(15) NT     Number of the node, in which the thermostat 4.0
C                    for the aux. heater is located
C     PAR(16) (UA)sa Heat loss capacity rate from store material 5.4
C                    to ambient                                
C     PAR(17) (UA)fa Heat loss capacity rate from heat transfer  7.2
C                    fluid to ambient
C     PAR(18) (UA)fs Heat transfer capacity rate of the heat    10.7
C                    exchanger
C     PAR(19) SIGMAs Weight factor for explicit/implicit          .25
C                    solution in the storage 
C     PAR(20) SIGMAf Weight factor for explicit/implicit          .25
C                    solution in the heat exchanger
C     PAR(21) Tf,ini Initial temperature in the heat exchanger  25.0
C     PAR(22) Ts,ini Initial temperature in the heat exchanger  21.0
C
C     INPUTS
C     XIN(1) Tf,i    Heat transfer fluid inlet temperature      24.7
C     XIN(2) dM f    Heat transfer fluid mass flow rate          4.8
C     XIN(3) Ts,i    Storage medium inlet temperature           21.0
C     XIN(4) dM s    Storage medium flow mass rate               8.6
C     XIN(5) T a     Ambient Temperature                        18.0
C     XIN(6) dQ aux  Auxiliary heater input                    500.0
C
C     OUTPUTS
C     OUT(1) T f,o    Heat transfer fluid outlet temperature    23.27
C     OUT(2) dM f     Heat transfer fluid flow rate              4.8
C     OUT(3) T s.o    Storage medium outlet temperature         21.0
C     OUT(4) dM s     Storage medium flow rate                   8.6
C     OUT(5) dQ l     Heat loss rate of the storage             18.0
C     OUT(6) dQ s     Heating power supplied to the load       500.0
C     OUT(7) T s      Mean storage temperature                  21.52
C     OUT(8) dQ f     Heating power supplied by the heat         6.725
C                       transfer fluid
C     OUT(9) dQ aux   Power supplied to auxiliary heater       500.0
C     OUT(10) Tf      Mean heat transfer fluid temperature      22.7
C     OUT(11) Ts,1    Storage temperature of node 1             21.87
C     OUT(12) Ts,2    Storage temperature of node 2             21.87
C     OUT(13) Ts,3    Storage temperature of node 3             21.87
C     OUT(14) Ts,4    Storage temperature of node 4             21.87
C     OUT(15) Ts,5    Storage temperature of node 5             21.39
C     OUT(16) Ts,6    Storage temperature of node 6             21.19
C     OUT(17) Ts,7    Storage temperature of node 7             21.06
C     OUT(18) Ts,8    Storage temperature of node 8             21.01
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      PARAMETER (NODEMX=50, PI=3.1417)
C
      DIMENSION XIN(10),OUT(70),PAR(22),INFO(15)
C
      COMMON/SIM/ TIME0,TFINAL,DELT
      COMMON/STORE/ NSTORE,IAV,S(5000)
	COMMON /LUNITS/LUR,LUW,IFORM,LUK
C
C     Number of INPUTS and PARAMETERS
C
      INTEGER NI,NP
C
C     Implicit or explicit method (1:implicit, 2:explicit)
C
      INTEGER METHOD
C
C     Weight factors for upwind/downwind differencing methods
C
      REAL BETAF,BETAS
      REAL MDOTS,MDOTF,CPMDS,CPMDF
C
C     Weight factors for implicit/explicit methods
C
C     0 = fully explicit, 1 = fully implicit,
C     between 0 and 1 partly implicit
C
      REAL SIGF,SIGS

C    Parameter values
C
C      VS    : Store volume                             [m3]
C      CPS   : Specific Heat Capacity of storage medium [kJ/kgK]
C      RHOS  : Density of storage medium                [kg/m3]
C      CONDS : Thermal conductivity storage medium      [kJ/hK]
C      VF    : Heat exchanger volume                    [m3]
C      CPF   : Specific Heat Capacity of heat transfer  [kJ/kgK]
C              medium
C      RHOF  : Density of heat transfer medium          [kg/m3]
C      STOHEI: Storage height                           [m]
C      HRTLEN: Length of auxiliary heater (top)         [m]
C      HTRPOS: auxiliary heater position (top, side, none)
C      NF    : number of nodes in the heat exchanger
C      NS    : number of nodes in the store
C      NDEAD : number of nodes in the dead zone
C      NFSTA : Nodal position of the heat exchanger inlet
C      NH    : number of immersion heater nodes in the model
C      UASA  : heat loss capacity rate store-ambient    [kJ/hK]
C      UAFA  : heat loss capacity rate fluid-ambient    [kJ/hK]
C      UA    : heat transfer capacity rate heat exch.  [kJ/hK]
C      UASTJ : turbulent mixing  coefficient in the store
C            : NOTE: this is set to 0
C      SIGF  : weight factor fluid  flow
C      SIGS  : weight factor store  flow
C
C derived quantities
C
C      FMCP  : heat capacity of the heat exchanger      [kJ/K]
C      SMCP  : heat capacity of the storag medium       [kJ/K]
C
C
C      Heat capacity of storage nodes
C
C
C      NFS   : Minimum of NF and NS
C
      REAL HTRLEN,STOHEI,CONDS,VS,VF,FMCP,SMCP
      INTEGER NF,NS,NDEAD,NFSTA,NH,NFS
C
C      coupling of nodes
C
C      IPTCFT :  Array for holding the postion node of the
C                heat exchanger fluid
C      IPTCST :  Array for holding the position node of the
C                storage medium
C      HSONOF :  Array for holding the on-off switch between the
C                heater nodes and the storage nodes
C      SAONOF :  Array for holding the on-off switch between the
C                storage nodes and the ambient temperature
C      SFONOF :  Array fo holding the on-off switch between the
C                storage nodes and the heat exchanger nodes
C
      INTEGER IPTCFT,IPTCST,HSONOF,SAONOF,SFONOF
      DIMENSION IPTCFT(NODEMX),IPTCST(NODEMX),HSONOF(NODEMX)
      DIMENSION SAONOF(NODEMX),SFONOF(NODEMX)
C
      DIMENSION CPMF(NODEMX),CPMS(NODEMX)
      DIMENSION OMF(NODEMX),OMS(NODEMX)
C
C      Nodal temperatures
C
      DIMENSION TS(NODEMX),TF(NODEMX),TDUM(NODEMX)
C
C      storage temperatures before mixing
C
      DIMENSION  TSNEW(NODEMX), QFS(NODEMX)
C
C      Temperatures
C
C      TFINI : initial heat exchanger temperature      [øC]
C      TSINI : initial storage temperature             [øC]
C      TSI   : storage inlet temperature               [øC]
C      TSE   : storage exit temperature                [øC]
C      TFI   : heat exchanger inlet temperature        [øC]
C      TFE   : heat exchanger exit temperature         [øC]
C      TAMB  : storage ambient temperature             [øC]
C
      REAL TFINI,TSINI,TSI,TSIN,TAMB,TFI,TFIN,TFE,TSE
C
C      capacity flows
C
C      CPMDS : capacity flow through storage
C      CPMDF : capacity flow through heat exchanger
C
C
C      Position of auxiliary heater (top, side, none)
C
      INTEGER HTRMODE
      CHARACTER HTRPOS*4

!    Set the version information for TRNSYS
      IF (INFO(7).EQ.-2) THEN
	  INFO(12) = 15
	  RETURN 1
	ENDIF

C
      UASTJ = 0.0
      INO=INFO(15)
      IF (INFO(7).NE.-1) GO TO 10
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     SET UP Parameter values
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C  FIRST CALL OF SIMULATION
      NI = 5
      IF (INFO(3).eq.6) NI=6
       NP = 21
      INFO(4)=NP
      INFO(6)=INT(PAR(12))+9
C
C Allocate storage in the S-array to save variables
C
       INFO(10)=(INT(PAR(12))+INT(PAR(11))+10)*2
C
C      Set unit no. to 0 to allocate parameter values to local variables
C
      IUNIT = 0
C
      CALL TYPECK(1,INFO,NI,NP,0)
       INO=INFO(10)
C
C
C  SET LOCAL VARIABLES TO PARAMETER VALUES.
C
10    IF (INFO(1).EQ.IUNIT) GO TO 30
      IUNIT=INFO(1)
C
C      Initial temperatures
C
      TFINI  =  PAR(21)
      TSINI  =  PAR(22)
      NS     =  INT(PAR(12))
      NF     =  INT(PAR(11))
C
C----Allocate initial temperatures
C
       S(INO+11)=TIME0
       S(INO+12)=TFINI
       S(INO+13)=TSINI
       S(INO+14)=XIN(5)
       S(INO+15)=0.0
       S(INO+16)=TFINI
       S(INO+17)=TSINI
       S(INO+18)=TSINI
       S(INO+19)=TSINI
       S(INO+20)=TFINI
C
C
C Initialize storage and fluid temperatures
C
      DO 380 J=1,NS
       IS   =INO+20+NS+J
       S(IS)=TSINI
380   CONTINUE
      DO 390 I=1,NF
       IS   =INO+20+2*NS+NF+I
       S(IS)=TFINI
390   CONTINUE
C
C----Set up initial conditions
C
      RETURN 1
C
30    VS     =  PAR(1)
      CPS    =  PAR(2)
      RHOS   =  PAR(3)
      CONDS  =  PAR(4)
      VF     =  PAR(5)
      CPF    =  PAR(6)
      RHOF   =  PAR(7)
      STOHEI =  PAR(8)
      HTRLEN =  PAR(9)
      HTRMODE=  INT(PAR(10))
      NF     =  INT(PAR(11))
      NS     =  INT(PAR(12))
      NDEAD  =  INT(PAR(13))
      NH     =  INT(PAR(14))
      NT     =  PAR(15)
      UASA   =  PAR(16)
      UAFA   =  PAR(17)
      UA     =  PAR(18)
      SIGS   =  PAR(19)
      SIGF   =  PAR(20)
C
      IF(SIGF.GT.0.99) METHOD=1
C
C
C Update information of n-1 timelevel if iteration was successful
C The data of the n-1 timelevel are stored in the s-array in the
C following manner:
C S(INO) -    S(INO+10)            variables needed for n-timelevel
C S(INO+21)-S(INO+20+NS)           nodal storage tempereratures
C S(INO+21+2*NS)-S(INO+20+2*NS+NF) nodal fluid temperatures
C
C The results from the last call are stored in the following manner
C S(INO+11)-S(INO+20)              variables needed for n-timelevel
C S(INO+21+NS)-S(INO+20+2*NS)      nodal storage temperatures
C S(INO+21+2*NS+NF)-S(INO+20+2*NS+2*NF) nodal fluid temperatures
C
C On the first call in a timestep, the results of the last call
C become the data of n-1 timelevel
C
       IF(INFO(7).EQ.0) THEN
      DO 116 IN=1,NS
       IS   =INO+20+IN
       IS1  =IS+NS
       S(IS)=S(IS1)
116   CONTINUE
      DO 117 I=1,NF
       IS   =INO+20+2*NS+I
       IS1  =IS+NF
       S(IS)=S(IS1)
117   CONTINUE
      DO 118 IP=1,10
       IS   =INO+IP
       IS1  =IS+10
       S(IS)=S(IS1)
118   CONTINUE
C
      ENDIF
C
       DO 381 J=1,NS
        IS   =INO+20+J
        TS(J)=S(IS)
381    CONTINUE
       DO 391 I=1,NF
        IS   =INO+20+2*NS+I
        TF(I)=S(IS)
391    CONTINUE
C
C derived quantities
C
      FMCP=  RHOF*VF*CPF
      SMCP=  RHOS*VS*CPS
C
      IF (HTRMODE.eq.0) HTRPOS = 'NONE'
      IF (HTRMODE.eq.1) HTRPOS = 'TOP'
      IF (HTRMODE.eq.2) HTRPOS = 'SIDE'
C
C      calculation using upwind differencing
C
      BETAF=0.
      BETAS=0.
C
C----Check for correct setup of NF, NS and NDEAD.
      IF((NF+NDEAD).GT.NS)THEN
       STOP' Error! NF+NDEAD must not greater than NS.'
      ENDIF
C
      NFS=MIN0(NF,NS)
C----Set up on-off contact between fluid, store and immersion
C    heater.
      IJ=NS-NDEAD-NF
      DO 310 I=1,NF
       IJ=IJ+1
       IPTCST(I)=IJ
310   CONTINUE
      DO 320 I=1,NS
       HSONOF(I)=0.
       SFONOF(I)=0.
       IPTCFT(I)=1
320   CONTINUE
      IJ=0
      NFSTA=NS-NDEAD-NF+1
      NFEND=NS-NDEAD
      DO 330 I=NFSTA,NFEND
       IJ=IJ+1
       SFONOF(I)=1.
       IPTCFT(I)=IJ
330   CONTINUE
      DO 340 I=1,NS
       IF(UAFA .GT. 0.)THEN
        SAONOF(I)=ABS(-1.+SFONOF(I))
       ELSE
        SAONOF(I)=1.
       ENDIF
340   CONTINUE
      IF(HTRPOS .EQ. 'TOP ') THEN
       NH=NINT(HTRLEN/STOHEI*FLOAT(NS))
       IF(NH .LE. 0)NH=1
       IF(NH .GT. NS)STOP' (NH) must not > (NS). Check HTRLEN !'
       DO 350 I=1,NH
        HSONOF(I)=1.
350    CONTINUE
      ENDIF
      IF(HTRPOS .EQ. 'SIDE')THEN
        HSONOF(NH)=1.
      ENDIF
C----Set up nodal heat capacities.
      DO 360 I=1,NF
       CPMF(I)=FMCP/FLOAT(NF)
       OMF(I)=0.0
360   CONTINUE
      DO 370 J=1,NS
       CPMS(J)=SMCP/FLOAT(NS)
       OMS(J)=0.0
370   CONTINUE
C
C+++++++++++++++++++++++++++++++++++++++
C
C End of set up
C
C++++++++++++++++++++++++++++++++++++++
C
C  SET INPUTS
C
      TFIN   =XIN(1)
      MDOTF  =XIN(2)
      TSIN   =XIN(3)
      MDOTS  =XIN(4)
      TAMB   =XIN(5)
      PAUX   =XIN(6)
      TIMN   =TIME
C
C
       TIMNM1=S(INO+1)
       TFINM1=S(INO+2)
       TSINM1=S(INO+3)
       TANM1 =S(INO+4)
       PAUNM1=S(INO+5)
       TFAV2 =S(INO+6)
       TSAV2 =S(INO+7)
       TSIJ2 =S(INO+8)
       TSEND2=S(INO+9)
       TFEND2=S(INO+10)
C
C
C  Return if
C  (second call from TRNSYS)
C
      IF(INFO(8).EQ.2) RETURN	1
C
C Start simulation
C
C
C Use UA=f(Re,Pr,Gr) is ISWICH=0, else use UA=const
C     IF(ISWICH(IDSET).EQ.0) CALL FUNCUA(MDOTF,IST(N),A,UA)
C
C --- Set counter for each simulation time step
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C $$  SETDT automates the calculation of each simulation time step    $$
C $$        size and sets the counting number(s) for each             $$
C $$        time step size.                                           $$
C $$  (note : This subroutine always set the simulation time step     $$
C $$          size less than or equal to the simulation time step     $$
C $$          size specified for TRNSYS.)                             $$
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ZERO=0.1E-8
      IF(UA .LT. 0.)PAUSE' UAFS<= 0! Check !!'
C----Compute storage nodal conductivity
      UASKJ=FLOAT(NS)*CONDS*VS/STOHEI**2
C---- Compute the simulation time step size for the heat store
      CPMDF=0.
      CPMDS=0.
      CPMDF=MDOTF*CPF
      CPMDS=MDOTS*CPS
      FNF=FLOAT(NF)
      FNS=FLOAT(NS)
      FNFS=FLOAT(NFS)
      DTS=SMCP/FNS/((CPMDS+UA/FNFS+UASA/FNS+UASKJ+UASTJ+ZERO)*
     :(1.0-SIGS+ZERO))
      DTF=FMCP/FNF/((CPMDF+UA/FNFS+UAFA/FNF+ZERO)*
     :(1.0-SIGF+ZERO))
C--- Pick up the smallest time step size
      DELTIM=AMIN1(DTF,DTS)*0.9999
C--- Set counter (NCT) to "1" if DELTIM is GE DELT
      IF(DELTIM.GE.DELT) THEN
       DELTIM=DELT
       NCT=1
      ENDIF
      IF(DTF.GE.DTS) GOTO 2
      IF((DELTIM.LT.DELT).AND.(METHOD.NE.1)) THEN
C
C----Using implicit method for the heat exchanger equation to allow
C    a large delta time if it is in discharge or stand-by mode.
C
      IF(((CPMDS.LT.0.1).AND.(CPMDF.LE.0.1)).OR.(CPMDF.LE.0.1))
     : THEN
1      SIGF=SIGF+0.1
       IF(SIGF.GE.1.0) SIGF=1.0
      DTF=FMCP/FNF/((CPMDF+UA/FNFS+UAFA/FNF+ZERO)*
     :(1.0-SIGF+ZERO))
      IF((DTF.LT.DTS).AND.(SIGF.LT.0.35)) GOTO 1
      DELTIM = AMIN1(DTF,DTS)*0.9999
       IF(DELTIM.GT.DELT) THEN
        DELTIM=DELT
        NCT=1
        GOTO 35
       ENDIF
      ENDIF
      ENDIF
2     IF(DELTIM.LE.0.) THEN
         WRITE(LUW,'(A)') ' DTS,DTF,FMCP,SMCP '
         WRITE(LUW,'(4E12.3)')  DTS,DTF,FMCP,SMCP
       STOP'DELTIM negative??, please check'
      ENDIF
C---- Set up the correct counter for DELTIM<DELT
      NCT=0
12    NCT=NCT+1
      IF((DELTIM*FLOAT(NCT)).LT.DELT) GOTO 12
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C END of SETDT
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
35    CONTINUE
C
      DELBAK = DELTIM
      TIMX   = TIMNM1
C
C Reset sums
C
      TFES = 0.
      TSES = 0.
      TFSIJS=0.
      TFIS  =0.
      TSIS  =0.
      TSAVS =0.
      TFAVS =0.
C
      DO 20 k=1,NCT
         IF(NCT.NE.1.AND.K.EQ.1) THEN
          DELTIM=DELT-(DELBAK*FLOAT(NCT-1))
         ELSE
          DELTIM=DELBAK
         ENDIF
C
C---- Set boundary conditions
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C $$  SETBC sets the boundary conditions.
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      TIMX=TIMX+DELTIM
       TFI=TFIN
       TSI=TSIN
C----Set boundary conditions.
      IF((CPMDF.GT.0.1).AND.(CPMDS.LT.0.1)) THEN
       UASKJJ=0.
      ELSE
       UASKJJ=UASKJ
      ENDIF
       IF (HTRPOS.EQ.'TOP')  PAUXJ=PAUX/FLOAT(NH)
       IF (HTRPOS.EQ.'SIDE') PAUXJ=PAUX
       IF (HTRPOS.EQ.'NONE') PAUXJ=0.
C
C---- Calculate Fluid and Store temperatures
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C $$  MODEL calculates the temperature                                $$
C $$  profile at the end of each simulation time step size.           $$
C $$                                                                  $$
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C----Set up heat transfer coefficients.
      UAFSIJ=UA/FLOAT(NFS)
      UAFAI=UAFA/FLOAT(NF)
      SDIA=SQRT(4.*VS/(PI*STOHEI))
      SASIDE=PI*SDIA*STOHEI
      TOPBOT=PI*SDIA*SDIA/4.
      SATOL=SASIDE+2.*TOPBOT
      UASATB=UASA*TOPBOT/SATOL
      IF(UAFA .GT. 0.)THEN
       FNS=FLOAT(NS-NF)
       IF(FNS .LE. 0.)FNS=1.0
       SIDEJ=SASIDE/FNS
C       UASAJJ=UASA/FNS
      ELSE
C       UASAJJ=UASA/FLOAT(NS)
       SIDEJ=SASIDE/FLOAT(NS)
      ENDIF
      UASAS=UASA*SIDEJ/SATOL
      SIGFDT=SIGF*DELTIM
      SIGSDT=SIGS*DELTIM
C----Calculate TS based on old time level of TF.
      DO 5 J=1,NS
       IF((J.EQ.1).OR.(J.EQ.NS))THEN
        UASAJJ=UASATB+UASAS
       ELSE
        UASAJJ=UASAS
       ENDIF
       UASKJ1=UASKJ
       UASKJ3=UASKJ
       UASTJ1=UASTJ
       UASTJ3=UASTJ
       QELEC=PAUXJ*HSONOF(J)
       UAFS=UAFSIJ*SFONOF(J)
       UASAJ=UASAJJ*SAONOF(J)
       TSJ=TS(J)
       JM1=J-1
       JP1=J+1
       TFII=TF(IPTCFT(J))
       IF(J .EQ. 1)THEN
        TSJM1=TSJ
        TSKJM1=TSJ
        UASKJ1=0.
        UASTJ1=0.
       ELSE
        TSKJM1=TS(JM1)
        TSJM1=TS(JM1)
       ENDIF
C----The lowermost storage node is fed from the cold water source
C    at temperature TSIN.
       IF(J.EQ.NS) THEN
        TSJP1=TSI
        TSKJP1=TSJ
        UASKJ3=0.
        UASTJ3=0.
       ELSE
        TSJP1=TS(JP1)
        TSKJP1=TS(JP1)
       ENDIF
       DS=BETAS*CPMDS*(2.0*TSJ-TSJM1-TSJP1)+CPMDS*
     :    (TSJP1-TSJ)+UAFS*(TFII-TSJ)+UASAJ*(TAMB-TSJ)+
     :    (UASKJ3+UASTJ3)*(TSKJP1-TSJ)+(UASKJ1+UASTJ1)*
     :    (TSKJM1-TSJ)+QELEC
       BS=CPMS(J)+SIGSDT*(CPMDS+UAFS+UASAJ+UASKJ3+UASTJ3
     :    +UASKJ1+UASTJ1-2.0*BETAS*CPMDS)
       OMS(J)=DS/BS
       TSNEW(J)=TS(J)
5     CONTINUE
C----Update to get new time level of TS
      DO 7 J=1,NS
       TSNEW(J)=TSNEW(J)+DELTIM*OMS(J)
7     CONTINUE
      CALL MIX(NS,TSNEW,NODEMX)
C----Calculate new time level of TF.
C    Use 0.5*(TS+TSNEW) for TSJ.
      DO 101 I=1,NF
       TFII=TF(I)
       IM1=I-1
       IP1=I+1
       TSJ=0.5*(TS(IPTCST(I))+TSNEW(IPTCST(I)))
C----The uppermost (heat exchanger) fluid node is fed from the
C    heating source at temperature TFI.
       IF(I.EQ.1) THEN
        TFIM1=TFI
       ELSE
        TFIM1=TF(IM1)
       ENDIF
       IF(I.EQ.NF)THEN
        TFIP1=TFII
       ELSE
        TFIP1=TF(IP1)
       ENDIF
       DF=BETAF*CPMDF*(2.0*TFII-TFIM1-TFIP1)+CPMDF*
     : (TFIM1-TFII)+UAFSIJ*(TSJ-TFII)+UAFAI*(TAMB-TFII)
       BF=CPMF(I)+SIGFDT*(CPMDF+UAFSIJ+UAFAI-2.*BETAF*CPMDF)
       QFS(I)=UAFSIJ*(TSJ-TFII)
       OMF(I)=DF/BF
101    CONTINUE
C----Calculate new time level of TS.
C    Use the same amount of QFS for TFII.
      DO 201 J=1,NS
       IF((J.EQ.1).OR.(J.EQ.NS))THEN
        UASAJJ=UASATB+UASAS
       ELSE
        UASAJJ=UASAS
       ENDIF
       UASKJ1=UASKJ
       UASKJ3=UASKJ
       UASTJ1=UASTJ
       UASTJ3=UASTJ
       QELEC=PAUXJ*HSONOF(J)
       UAFS=UAFSIJ*SFONOF(J)
       UASAJ=UASAJJ*SAONOF(J)
       TSJ=TS(J)
       JM1=J-1
       JP1=J+1
       IF(J .EQ. 1)THEN
        TSJM1=TSJ
        TSKJM1=TSJ
        UASKJ1=0.
        UASTJ1=0.
       ELSE
        TSKJM1=TS(JM1)
        TSJM1=TS(JM1)
       ENDIF
C----The lowermost storage node is fed from the cold water source
C    at temperature TSI.
       IF(J.EQ.NS) THEN
        TSJP1=TSI
        TSKJP1=TSJ
        UASKJ3=0.
        UASTJ3=0.
       ELSE
        TSJP1=TS(JP1)
        TSKJP1=TS(JP1)
       ENDIF
       QSF=-QFS(IPTCFT(J))*SFONOF(J)
       DS=BETAS*CPMDS*(2.0*TSJ-TSJM1-TSJP1)+CPMDS*
     :    (TSJP1-TSJ)+QSF+UASAJ*(TAMB-TSJ)+
     :    (UASKJ3+UASTJ3)*(TSKJP1-TSJ)+(UASKJ1+UASTJ1)*
     :    (TSKJM1-TSJ)+QELEC
C       DS=BETAS*CPMDS*(2.0*TSJ-TSJM1-TSJP1)+CPMDS*
C     :    (TSJP1-TSJ)+UAFS*(TFII-TSJ)+UASAJ*(TAMB-TSJ)+
C     :    (UASKJ3+UASTJ3)*(TSKJP1-TSJ)+(UASKJ1+UASTJ1)*
C     :    (TSKJM1-TSJ)+QELEC
       BS=CPMS(J)+SIGSDT*(CPMDS+UAFS+UASAJ+UASKJ3+UASTJ3
     :    +UASKJ1+UASTJ1-2.0*BETAS*CPMDS)
       OMS(J)=DS/BS
201    CONTINUE
C---- for direct charge/discharge storage (no HX) use 
C     flowrate weighted OMF and OMS
C
      SMDOT = MDOTF+MDOTS
       RMD1=0.5
       RMD2=0.5
      IF (SMDOT.GT.ZERO) RMD1=MDOTF/SMDOT
      IF (SMDOT.GT.ZERO) RMD2=MDOTS/SMDOT
C
C----Update
      DO 301 I=1,NF
      IF(UA.LT.ZERO) OMF(I)=RMD1*OMF(I)+RMD2*OMS(I)
      TF(I)=TF(I)+DELTIM*OMF(I)
301    CONTINUE
      DO 401 J=1,NS
      IF(UA.LT.ZERO) OMS(J)=RMD1*OMF(J)+RMD2*OMS(J)
      TS(J)=TS(J)+DELTIM*OMS(J)
401    CONTINUE
C----Mix any storage node which has a hotter fluid layer
C    beneath it.
      CALL MIX(NS,TS,NODEMX)
C
C----Sum up data
C
      TFES = TFES+TF(NF)*DELTIM
      TSES = TSES+TS(1)*DELTIM
      TFIS = TFIS+TFI*DELTIM
      TSIS = TSIS+TSI*DELTIM
      TFAVS= TFAVS+TMEAN(TF,NF)*DELTIM
      TSAVS= TSAVS+TMEAN(TS,NS)*DELTIM
20    CONTINUE
C
C
C----Locate the temperatures of the storage nodes which are in contact
C    with the heat exchanger fluid nodes.
        DO 31 I=1,NF
         TDUM(I)=TS(IPTCST(I))
31      CONTINUE
        TFSIJS=TFSIJS+TMEAN(TDUM,NF)
C----Find the mean storage nodal temperatures which are in contact with
C    the ambient temperature. This is necessary for calculating the
C    store heat loss for those types of mantle heat exchanger.
        IF(UAFA .GT. 0.)THEN
         J=0.
         DO 40 I=1,NS
          TSAJ=TS(I)*SAONOF(I)
          IF(TSAJ .GT. 0.)THEN
           J=J+1
           TDUM(J)=TSAJ
          ENDIF
40       CONTINUE
         IF((NS-NF) .EQ. 0)THEN
          TSAJS=0.
         ELSE
          TSAJS=TSAJS+TMEAN(TDUM,(NS-NF))
         ENDIF
        ENDIF
C----Average and store data.

       TFE   =TFES/DELT
       TSE   =TSES/DELT
       TFAV  =TFAVS/DELT
       TSAV  =TSAVS/DELT
       TFEND =TMEAN(TF,NF)
       TSEND =TMEAN(TS,NS)
       TSIJAV=TFSIJS/DELT
       TFI   =TFIS/DELT
       TSI   =TSIS/DELT
C
C Check for temperature limits
C
       IF((TSEND.GE.120.).OR.(TSEND.LE.0.)) THEN
         WRITE(LUW,'(A,E12.2)')' Storage Temperature exceeds limits',TSEND
       ENDIF
C
C Calculate heat loss during time step
C
       QLOSS = UASA*(TSAV-TAMB)*DELT
C
C Calculate heat added to the store
C
       QCHAR = CPMDF*(TFI-TFE)*DELT
C
C Calculate heat extracted from store
C
       QDCAR = CPMDS*(TSE-TSI)*DELT
C
C Change in internal energy content
C
       DQS   = SMCP*(TSEND-TSEND2) + FMCP*(TFEND-TFEND2)
C
C Check heat balance
C
       DQSERR=(QCHAR-QDCAR-QLOSS-DQS)/(SMCP+FMCP)
C
       IF(DQSERR.GT.0.1) THEN
         WRITE(LUW,'(A,F6.2)') ' Heat balance error in K ',DQSERR
       ENDIF
C
C----Update information.
C
       S(INO+11)=TIME
       S(INO+12)=TFIN
       S(INO+13)=TSIN
       S(INO+14)=TAMB
       S(INO+15)=PAUX
       S(INO+16)=TFAV
       S(INO+17)=TSAV
       S(INO+18)=TSIJAV
       S(INO+19)=TSEND
       S(INO+20)=TFEND
C
C
C----Allocate local variables to output array
C
C
      OUT(1)=TFE
      OUT(2)=MDOTF
      OUT(3)=TSE
      OUT(4)=MDOTS
      OUT(5)=QLOSS/DELT
      OUT(6)=QDCAR/DELT
      OUT(7)=TSEND
      OUT(8)=QCHAR/DELT
      OUT(9)=PAUX
      OUT(10)=TFEND
      DO 110 IN=11,NS+10
       INODE=IN-10
       OUT(IN)=TS(INODE)
       IS=INO+20+NS+INODE
       S(IS)=TS(INODE)
110   CONTINUE
      DO 111 I=1,NF
       IS=INO+20+2*NS+NF+I
       S(IS)=TF(I)
111   CONTINUE
C
C
      RETURN 1
      END
C-----------------------------------------------------------------------
      SUBROUTINE INTPOL(XM1,X,YM1,Y,XI,YI)
      A=(X-XI)/(X-XM1)
      B=1.-A
      YI=A*YM1+B*Y
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE MIX(NODE,T,NODEMX)
      DIMENSION T(NODEMX)
      IF(NODE .EQ. 1)RETURN
1     NMAX=1
C----Find the lowest unstable node, NMAX.
      DO 1000 J=2,NODE
      TJ=T(J)
      TJM1=T(J-1)
      IF(TJ.GT.TJM1) THEN
       NMAX=J
      ENDIF
 1000 CONTINUE
      IF(NMAX.EQ.1) RETURN
C----Find the highest unstable node, NMIN.
      NMIN=NMAX
      DO 2000 I=2,NMAX
      J=NMAX-I+1
      IF(T(J).LT.T(NMAX)) THEN
       NMIN=J
      ENDIF
 2000 CONTINUE
C----Compute the average.
      TAV=T(NMIN)
      NMINP1=NMIN+1
      DO 3000 J=NMINP1,NMAX
      TAV=TAV+T(J)
 3000 CONTINUE
      NDIFF=NMAX-NMIN+1
      IF((NMAX-NMIN).LT.1) THEN
       STOP ' Some error in mix subroutine. Check NMIN and NMAX.'
      ENDIF
      TAV=TAV/FLOAT(NDIFF)
      DO 4000 J=NMIN,NMAX
      T(J)=TAV
 4000 CONTINUE
      GOTO 1
      END
C-----------------------------------------------------------------------
      REAL FUNCTION TMEAN(X,N)
      DIMENSION X(N)
      AV=0.
      DO 10 I=1,N
      AV=AV+X(I)
10    CONTINUE
      TMEAN=AV/FLOAT(N)
      END