qq:478672328(有偿服务)  告诉下我什么意思  需要做个外部的TXT数据作为输入

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

 

C-----------------------------------------------------------------------------------------------------------------------

C    DESCRIPTION:

C     THIS SUBROUTINE MODELS A VERTICAL CYLINDRICAL STORAGE TANK WITH AN IMMERSED

C     HEAT EXCHANGER. THIS ROUTINE SOLVES THE COUPLED DIFFERENTIAL EQUATIONS

C     IMPOSED BY CONSIDERING THE MASS IN THE STORAGE TANK AND THE MASS IN THE

C     HEAT EXCHANGER.

C  

C     NEW FEATURES:

C        - USERS CAN NOW SPECIFY THE FRACTION OF INLET FLOW THAT GOES TO EACH NODE

C        - LOSSES TO A GAS FLUE CAN NOW BE SPECIFIED (THESE LOSSES ARE SET TO ZERO

C            IF AUX ENERGY IS BEING ADDED TO THE TANK

C        - USERS CAN ADD/REMOVE ENERGY FROM ANY NODE BY THE USE OF "MISCELLANEOUS"

C            ENERGY FLOWS.  THIS FEATURE CAN BE USED TO SIMULATE A PILOT LIGHT.

C        - TEMPERATURE SEEKING MODE FOR THE FLUID INLETS HAS BEEN ADDED

C        - NEW NODAL EDGE LOSSES SO USERS CAN PARTIALLY BURY THEIR TANK ETC...

C        - MULTIPLE HEAT EXCHANGERS NOW ADDED

C

C     CARE SHOULD BE TAKEN WHEN CONSIDERING HIGH MASS FLOW RATES OF MIXING BETWEEN

C     NODES AS THIS CAN CAUSE AN ENERGY BALANCE PROBLEM WITH THE MODEL IF THE "CONVERGED"

C     PARAMETER IS TOO LARGE.

C     

C     UPDATES:

C        JULY 11TH, 2004 - FIXED A BUG WITH 1-NODE TANKS (FLOW_OUT_LOAD)

C        JULY 20TH, 2004 - FIXED A BUG WITH TANK VOLUME CALCULATION FOR HX REDUCTIONS

C        JULY 30TH, 2004 - REPLACED A FEW EXP( WITH DEXP(

C        SEPTEMBER 28TH, 2004 - ADDED THE ABILITY TO READ THE PARAMETERS FROM A DATA FILE SO WE CAN HAVE

C           A NICE APPLICATION PROGRAM THROUGH THE STUDIO

C        NOVEMBER 2, 2004 - FIXED A BUG IN THE TNK VOLUME CALCULATION FOR IHXGEOM=1

C        MARCH 25TH, 2005 - RESET THE NUMBER OF HX NODES TO ZERO FOR SYSTEMS WITHOUT A HX

C-----------------------------------------------------------------------------------------------------------------------

! Copyright ?2004 Thermal Energy System Specialists, LLC. All rights reserved.

 

C-----------------------------------------------------------------------------------------------------------------------

      !Export this subroutine for its use in external DLLs.

      !DEC$ATTRIBUTES DLLEXPORT :: TYPE534

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    ACCESS TRNSYS FUNCTIONS

USE TrnsysConstants

USE TrnsysFunctions

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    TRNSYS DECLARATIONS

      IMPLICIT NONE

DOUBLE PRECISION XIN,OUT,TIME,PAR,T,DTDT,STORED,TIME0,TFINAL,DELT

      INTEGER*4 INFO(15),NP,NI,NOUT,ND,IUNIT,ITYPE,ICNTRL,NPAR,NIN,

1   N_HX_MAX

      CHARACTER*3 YCHECK,OCHECK

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    USER DECLARATIONS

INTEGER MAXNODES,MAXHXNODES,MISCMAX,MAXPORTS,NSTORED,NITEMS


      PARAMETER (MISCMAX=10,MAXPORTS=10,MAXNODES=20,N_HX_MAX=5)

PARAMETER (NIN=4+2*MAXNODES+MISCMAX+2*MAXPORTS+2*N_HX_MAX)

PARAMETER (MAXHXNODES=20)

      PARAMETER (NSTORED=4*MAXNODES+2*N_HX_MAX*MAXHXNODES*MAXNODES)

      PARAMETER (NPAR=21+MAXPORTS*(2+MAXNODES)+2*MAXNODES+MISCMAX+

1   N_HX_MAX*17+MAX(MAXNODES,2*MAXHXNODES)*N_HX_MAX)

      PARAMETER (NOUT=11+3*MAXPORTS+MAXNODES+7*N_HX_MAX)

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    REQUIRED TRNSYS DIMENSIONS

      DIMENSION XIN(NIN),OUT(NOUT),PAR(NPAR),YCHECK(NIN),OCHECK(NOUT),

1   STORED(NSTORED),T(MAXNODES),DTDT(MAXNODES)

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES

INTEGER ITS_DT(MAXNODES),ITERS_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

1   N_HX,ITERS_TANK(MAXNODES),N_HXFLOW(MAXHXNODES,N_HX_MAX),NODES,

     1   IPASSES(MAXNODES,N_HX_MAX),MAXITERS,INLET_MODE(MAXPORTS),

     1   IHXGEOM(N_HX_MAX),K,N_HXNODES(N_HX_MAX),N_TNKDIV(N_HX_MAX),J,

     1   I,ITFLUID,N_ENTRY(MAXPORTS),NSINK,N_EXIT(MAXPORTS),LU_DATA,

     1   IHXFLUID(N_HX_MAX),ITERS_TANK_MAX,ITERS_HX_MAX,

     1   ITERS_SURF,NODE_CRIT,ICASE,JJ,LL,KK,N_MISC,NODE_MISC(MISCMAX),

1   N_PORTS,N_SPOT,NODE_MIN,N_STORED

      DOUBLE PRECISION K_WALL(N_HX_MAX),L_TUBE(N_HX_MAX),VAL(NPAR),

1   N_TUBES(N_HX_MAX),PI,T_ELOSS(MAXNODES),Q_HX(N_HX_MAX),

1   XNODES,VOLTANK,TANKHEIGHT,U_EDGE(MAXNODES),U_BOTTOM,T_DIFF,

1   U_TOP,FRAC_IN(MAXPORTS,MAXNODES),FRAC_IN_TOT,BETA_TANK,

1   C(N_HX_MAX),RA_EXP(N_HX_MAX),GF(N_HX_MAX),GF_EXP(N_HX_MAX),

     1   HEADERVOL(N_HX_MAX),ACROSS_HX(N_HX_MAX),DIA_COIL(N_HX_MAX),

1   TOTFRAC,VOLTUBE_I(N_HX_MAX),VOLTUBE_O(N_HX_MAX),RADIUS,

     1   TFLOW_IN(MAXPORTS),FLOW_IN(MAXPORTS),THXFLOW_IN(N_HX_MAX),

     1   FLOWHX_IN(N_HX_MAX),T_TLOSS,T_BLOSS,DELT_USED,DELT_NOW,

     1   QHX_TOT,QLOSS_E,QLOSS_B,QLOSS_T,QAUX_TOT,QSTORE_HX(N_HX_MAX),

     1   QSTORE_T,QDELIV,QDELIV_HX(N_HX_MAX),AVETANK,DT_OLD,FX,

     1   TOT_VOL,CP_TANK,RHO_TANK,K_TANK,MU_TANK,AVEHX,TSURF_OLD,GRAV,

     1   T_PROPS,T_REF,CP_REF,RHO_REF,MU_REF,K_REF,CP_NODE,K_NODE,

     1   MU_NODE,RHO_NODE,BETA,ALPHA_NODE,RA_NODE,H_OUT,FLOW_TUBE,

1   RE_HX,PR_HX,RE_CRIT,HE,H_IN,Q_TEMP,TSURF_I,TSURF_O,

     1   FLOW_MAINS(MAXPORTS,MAXNODES),UA_TOP,UA_EDGE,QMISC_TOT,

     1   T_IN_LOAD(MAXPORTS,MAXNODES),T_FLUE,QLOSS_F,QDEL(MAXPORTS),

     1   FLOW_IN_LOAD(MAXPORTS,MAXNODES),AVETEMP_HX(N_HX_MAX),

     1   FLOW_OUT_LOAD(MAXPORTS,MAXNODES),UA_BOT,TINIT_HX(N_HX_MAX),

1   TEMPHX_O(N_HX_MAX),AVETEMP_T,TOTCAP_T,EB_TOP_T,EB_BOT_T,

1   EB_ERROR_T,EB_TOP_HX,EB_BOT_HX,EB_ERROR_HX(N_HX_MAX),DT_CRIT,

     1   FPX,DT_NEW,TMIX,CMIX,U_T,V_T,QTMIX,UA_FLUE(MAXNODES),

     1   OA_TUBE(N_HX_MAX),IA_TUBE(N_HX_MAX),OA_NODE,IA_NODE,L_NODE,

     1   TANKVOL(MAXNODES),HXVOL_O(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   PERVOL_T,PERVOL_HX(N_HX_MAX),

     1   TAVE_HX(MAXNODES,MAXHXNODES,N_HX_MAX),IN_DIA(N_HX_MAX),

     1   OUT_DIA(N_HX_MAX),OUT_RAD(N_HX_MAX),

     1   IN_RAD(N_HX_MAX),

     1   TI_TANK(MAXNODES),TF_TANK(MAXNODES),TAVE_TANK(MAXNODES),

     1   TI_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   TF_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

1   HXVOL_I(MAXNODES,MAXHXNODES,N_HX_MAX),CAP_T(MAXNODES),

     1   Q_MISC(MISCMAX),

     1   CAP_HX(MAXNODES,MAXHXNODES,N_HX_MAX),CP_HX(N_HX_MAX),

     1   K_HX(N_HX_MAX),MU_HX(N_HX_MAX),RHO_HX(N_HX_MAX),

     1   AA(MAXNODES),BB(MAXNODES),NU_NODE,

     1   UA_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   NUSSELT,CONVERGED,FLOW_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   TIN_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   AA_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   BB_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   TF_HX_NEW(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   TAVE_HX_NEW(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   TAVE_TANK_NEW(MAXNODES),

     1   TF_TANK_NEW(MAXNODES),K_ADDCOND,ATOP(MAXNODES),AEDGE(MAXNODES),

     1   ABOT(MAXNODES),ACOND(MAXNODES),FLOWI_MIX(MAXNODES),

     1   FLOWF_MIX(MAXNODES),MIXFLOW,NUSSELT_HX,TANKAVE(MAXNODES),

     1   TSTART_TANK(MAXNODES),TSTART_HX(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   WALL_RESIST(MAXNODES,MAXHXNODES,N_HX_MAX),L_COND(MAXNODES),

     1   PITCH_COIL(N_HX_MAX),

     1   RESIST_OUT(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   RESIST_IN(MAXNODES,MAXHXNODES,N_HX_MAX),

     1   DELT_MIN,HXAVE(MAXNODES,MAXHXNODES,N_HX_MAX),Q_AUX(MAXNODES),

1   FRACHX(MAX(MAXHXNODES,MAXNODES),N_HX_MAX)

      CHARACTER (LEN=MaxmessageLength) MESSAGE1,MESSAGE2

CHARACTER*1 ADUM

LOGICAL NOTCONV_HX,CHECKED(MAXNODES),LAMINAR_HX,HX_CONVERGE,

1   NOTCONV_TANK,TANK_CONVERGE,DONT_MIX,DONT_COND,DONT_LOSS,

     1   CHECKING_TEMPS,HAVE_AUX

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    DATA STATEMENTS

      DATA PI/3.14159265358979/,CONVERGED/0.0001/,DELT_MIN/0.001/

DATA DONT_MIX/.FALSE./,DONT_LOSS/.FALSE./

MESSAGE1='The current TRNSYS timestep is less than the minimum tim

1estep set in this model.  Either increase the TRNSYS timestep or d

     1ecrease the minimum timestep allowed by this model (DELT_MIN).'

MESSAGE2='The storage tank model was unable to read the parameter

1values from the supplied data file.  Please check the location and

     1 format of this data file carefully.'

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    GET GLOBAL TRNSYS SIMULATION VARIABLES

      TIME0=getSimulationStartTime()

      TFINAL=getSimulationStopTime()

      DELT=getSimulationTimeStep()

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    SET THE VERSION INFORMATION FOR TRNSYS

      IF(INFO(7).EQ.-2) THEN

   INFO(12)=16

   RETURN 1

ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE

      IF (INFO(8).EQ.-1) THEN

   RETURN 1

ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    UPDATE THE VAL ARRAY WHICH IS USED AS A SUBSTITUTE FOR THE PAR ARRAY

      DO J=1,INFO(4)

   VAL(J)=PAR(J)

ENDDO

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    PERFORM ANY "AFTER-ITERATION" MANIPULATIONS THAT ARE REQUIRED

      IF(INFO(13).GT.0) THEN

 

   IUNIT=0

 

C       RETRIEVE THE CRITICAL PARAMETERS TO KNOW HOW TO SET THE STORAGE ARRAY

         LU_DATA=JFIX(VAL(1)+0.1)

         NODES=JFIX(VAL(2)+0.1)

         N_PORTS=JFIX(VAL(3)+0.1)

   N_HX=JFIX(VAL(4)+0.1)

         N_MISC=JFIX(VAL(5)+0.1)

 

         IF(LU_DATA.GE.10) THEN

      REWIND(LU_DATA)

      READ(LU_DATA,*,END=10,ERR=9000) ADUM

      READ(LU_DATA,*,END=10,ERR=9000) ADUM

      DO J=1,NPAR

         READ(LU_DATA,*,END=10,ERR=10) VAL(5+J)

      ENDDO

   ENDIF

 

10    ITFLUID=JFIX(VAL(8)+0.1)      

         IF(ITFLUID.EQ.0) THEN

      N_SPOT=16+NODES

         ELSE IF(ITFLUID.EQ.1) THEN

      N_SPOT=11+NODES

         ELSE IF(ITFLUID.EQ.2) THEN

      N_SPOT=12+NODES

         ELSE IF(ITFLUID.EQ.3) THEN

      N_SPOT=12+NODES

   ENDIF

 

   DO J=1,N_PORTS

            INLET_MODE(J)=JFIX(VAL(N_SPOT+1)+0.1)

 

            IF(INLET_MODE(J).EQ.1) THEN

         N_SPOT=N_SPOT+3

            ELSE IF(INLET_MODE(J).EQ.2) THEN

         N_SPOT=N_SPOT+2+NODES

      ELSE

         N_SPOT=N_SPOT+2

      ENDIF

   ENDDO

         N_SPOT=N_SPOT+NODES+N_MISC

 

   N_HXNODES=0

   DO K=1,N_HX

      IHXGEOM(K)=JFIX(VAL(N_SPOT+1)+0.1)

            N_HXNODES(K)=JFIX(VAL(N_SPOT+2)+0.1)

            IHXFLUID(K)=JFIX(VAL(N_SPOT+3)+0.1)                      

 

            IF(IHXFLUID(K).EQ.0) THEN

               N_SPOT=N_SPOT+18

            ELSE IF(IHXFLUID(K).EQ.1) THEN

               N_SPOT=N_SPOT+14

            ELSE IF(IHXFLUID(K).EQ.2) THEN

               N_SPOT=N_SPOT+15

            ELSE IF(IHXFLUID(K).EQ.3) THEN

               N_SPOT=N_SPOT+15

      ENDIF

 

      IF(IHXGEOM(K).EQ.3) THEN

         N_SPOT=N_SPOT+2

      ENDIF

      

  IF(IHXGEOM(K).EQ.1) THEN

         N_TNKDIV(K)=NODES

         N_SPOT=N_SPOT+NODES

      ELSE

         N_TNKDIV(K)=1

         N_SPOT=N_SPOT+2*N_HXNODES(K)

      ENDIF

   ENDDO

 

C       CALCULATE THE AMOUNT OF STORAGE SPOTS NEEDED

         NITEMS=0

         DO K=1,N_HX

            NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

   ENDDO

   NITEMS=2*NITEMS+4*NODES

 

C       RETRIEVE THE STORAGE ARRAY

   CALL getStorageVars(STORED,NITEMS,INFO)

   IF(ErrorFound()) RETURN 1

 

C       SET THE INITIAL TANK TEMPS TO THE FINAL TANK TEMPS FROM THE PREVIOUS TIMESTEP

   DO J=1,NODES

      TF_TANK(J)=STORED(NODES+J)

      FLOWF_MIX(J)=STORED(3*NODES+J)

      STORED(J)=TF_TANK(J)

      STORED(2*NODES+J)=FLOWF_MIX(J)

         ENDDO

 

C       GET THE RIGHT SPOT IN THE STORAGE ARRAY FOR THE FINAL TEMPS

   N_STORED=4*NODES

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

C       RETRIEVE THE FINAL TEMPERATURES

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           TF_HX(I,J,K)=STORED(N_STORED+1)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

C       SET THE NEW INITIAL TEMPERATURES

   N_STORED=4*NODES

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           STORED(N_STORED+1)=TF_HX(I,J,K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

C       SET THE STORAGE VARIABLES AND GET OUT

   CALL setStorageVars(STORED,NITEMS,INFO)

 

   RETURN 1

ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE

      IF (INFO(7).EQ.-1) THEN

 

C       RETRIEVE THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT FROM THE INFO ARRAY

         IUNIT=INFO(1)

   ITYPE=INFO(2)

 

C       SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK

         LU_DATA=JFIX(VAL(1)+0.1)

         NODES=JFIX(VAL(2)+0.1)

         N_PORTS=JFIX(VAL(3)+0.1)

   N_HX=JFIX(VAL(4)+0.1)

         N_MISC=JFIX(VAL(5)+0.1)

 

   IF((NODES.LT.1).OR.(NODES.GT.MAXNODES))

1      CALL TYPECK(-4,INFO,0,2,0)

   IF((N_PORTS.LT.0).OR.(N_PORTS.GT.MAXPORTS))

1      CALL TYPECK(-4,INFO,0,3,0)

   IF((N_HX.LT.0).OR.(N_HX.GT.N_HX_MAX))

1      CALL TYPECK(-4,INFO,0,4,0)

   IF((N_MISC.LT.0).OR.(N_MISC.GT.MISCMAX))

1      CALL TYPECK(-4,INFO,0,5,0)

 

         IF(LU_DATA.GE.10) THEN

      REWIND(LU_DATA)

      READ(LU_DATA,*,END=40,ERR=9000) ADUM

      READ(LU_DATA,*,END=40,ERR=9000) ADUM

      DO J=1,NPAR-5

         READ(LU_DATA,*,END=40,ERR=40) VAL(5+J)

      ENDDO

   ENDIF

 

40    ITFLUID=JFIX(VAL(8)+0.1)      

         IF(ITFLUID.EQ.0) THEN

      N_SPOT=16+NODES

         ELSE IF(ITFLUID.EQ.1) THEN

      N_SPOT=11+NODES

         ELSE IF(ITFLUID.EQ.2) THEN

      N_SPOT=12+NODES

         ELSE IF(ITFLUID.EQ.3) THEN

      N_SPOT=12+NODES

   ENDIF

 

   DO J=1,N_PORTS

            INLET_MODE(J)=JFIX(VAL(N_SPOT+1)+0.1)

 

            IF(INLET_MODE(J).EQ.1) THEN

         N_SPOT=N_SPOT+3

            ELSE IF(INLET_MODE(J).EQ.2) THEN

         N_SPOT=N_SPOT+2+NODES

      ELSE

         N_SPOT=N_SPOT+2

      ENDIF

   ENDDO

         N_SPOT=N_SPOT+NODES+N_MISC

 

   N_HXNODES=0

   DO K=1,N_HX

      IHXGEOM(K)=JFIX(VAL(N_SPOT+1)+0.1)

            N_HXNODES(K)=JFIX(VAL(N_SPOT+2)+0.1)

            IHXFLUID(K)=JFIX(VAL(N_SPOT+3)+0.1)                      

 

      IF((IHXGEOM(K).LT.1).OR.(IHXGEOM(K).GT.4))

1         CALL TYPECK(-4,INFO,0,N_SPOT+1,0)

      IF(N_HXNODES(K).GT.MAXHXNODES)

1         CALL TYPECK(-4,INFO,0,N_SPOT+2,0)

      IF((IHXFLUID(K).LT.0).OR.(IHXFLUID(K).GT.3))

1         CALL TYPECK(-4,INFO,0,N_SPOT+3,0)

 

            IF(IHXFLUID(K).EQ.0) THEN

               N_SPOT=N_SPOT+18

            ELSE IF(IHXFLUID(K).EQ.1) THEN

               N_SPOT=N_SPOT+14

            ELSE IF(IHXFLUID(K).EQ.2) THEN

               N_SPOT=N_SPOT+15

            ELSE IF(IHXFLUID(K).EQ.3) THEN

               N_SPOT=N_SPOT+15

      ENDIF

 

      IF(IHXGEOM(K).EQ.3) THEN

         N_SPOT=N_SPOT+2

      ENDIF

 

      IF(IHXGEOM(K).EQ.1) THEN

         N_TNKDIV(K)=NODES

         N_SPOT=N_SPOT+NODES

      ELSE

         N_TNKDIV(K)=1

         N_SPOT=N_SPOT+N_HXNODES(K)*2

      ENDIF

   ENDDO

 

         IF(LU_DATA.GE.10) THEN

      NP=5

   ELSE

      NP=N_SPOT

   ENDIF

   NI=2*N_PORTS+2*N_HX+4+2*NODES+N_MISC

   ND=NODES

 

         INFO(6)=11+3*N_PORTS+NODES+7*N_HX

         INFO(9)=1

         INFO(10)=0

       

C       CALL THE TYPE CHECK SUBROUTINE TO COMPARE WHAT THIS COMPONENT REQUIRES TO WHAT IS SUPPLIED IN THE INPUT FILE

   CALL TYPECK(1,INFO,NI,NP,ND)

 

C       SET THE YCHECK AND OCHECK ARRAYS TO CONTAIN THE CORRECT VARIABLE TYPES FOR THE INPUTS AND OUTPUTS

         DO J=1,N_PORTS

      YCHECK(2*J-1)='TE1'

      YCHECK(2*J)='MF1'

   ENDDO

 

         DO J=1,N_HX

      YCHECK(2*N_PORTS+2*J-1)='TE1'

      YCHECK(2*N_PORTS+2*J)='MF1'

   ENDDO

 

   N_SPOT=2*N_PORTS+2*N_HX

 

         YCHECK(N_SPOT+1)='TE1'

 

 

   DO J=1,NODES

            YCHECK(N_SPOT+1+J)='TE1'

         ENDDO

   N_SPOT=N_SPOT+1+NODES

 

         YCHECK(N_SPOT+1)='TE1'

   YCHECK(N_SPOT+2)='TE1'

         YCHECK(N_SPOT+3)='MF1'

 

   DO J=1,NODES

      YCHECK(N_SPOT+3+J)='PW1'

   ENDDO

   

   DO J=1,N_MISC

     YCHECK(N_SPOT+3+NODES+J)='PW1'

   ENDDO

 

         DO I=1,N_PORTS

            OCHECK(2*I-1)='TE1'

            OCHECK(2*I)='MF1'

   ENDDO

 

   OCHECK(2*N_PORTS+1)='TE1'

   OCHECK(2*N_PORTS+2)='PW1'

 

         DO J=1,N_PORTS

      OCHECK(2*N_PORTS+2+J)='PW1'

   ENDDO

 

   OCHECK(3*N_PORTS+3)='PW1'

   OCHECK(3*N_PORTS+4)='PW1'

   OCHECK(3*N_PORTS+5)='PW1'

   OCHECK(3*N_PORTS+6)='PW1'

   OCHECK(3*N_PORTS+7)='PW1'

   OCHECK(3*N_PORTS+8)='PW1'

   OCHECK(3*N_PORTS+9)='PW1'

   OCHECK(3*N_PORTS+10)='PW1'

   OCHECK(3*N_PORTS+11)='PC1'

 

   DO J=1,NODES

      OCHECK(3*N_PORTS+11+J)='TE1'

         ENDDO

 

   DO J=1,N_HX

      OCHECK(3*N_PORTS+NODES+11+2*J-1)='TE1'

      OCHECK(3*N_PORTS+NODES+11+2*J)='MF1'

   ENDDO

 

   DO J=1,N_HX

      OCHECK(3*N_PORTS+11+NODES+2*N_HX+J)='TE1'

   ENDDO

 

   DO J=1,N_HX

      OCHECK(3*N_PORTS+11+NODES+3*N_HX+1+4*(J-1))='PW1'

      OCHECK(3*N_PORTS+11+NODES+3*N_HX+2+4*(J-1))='PW1'

      OCHECK(3*N_PORTS+11+NODES+3*N_HX+3+4*(J-1))='PW1'

      OCHECK(3*N_PORTS+11+NODES+3*N_HX+4+4*(J-1))='PC1'

   ENDDO

 

C       CALL THE RCHECK SUBROUTINE TO SET THE CORRECT INPUT AND OUTPUT TYPES FOR THIS COMPONENT

         CALL RCHECK(INFO,YCHECK,OCHECK)

 

C       SET THE NUMBER OF STORAGE SPOTS NEEDED FOR THIS COMPONENT

         NITEMS=0

         DO K=1,N_HX

            NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

   ENDDO

   NITEMS=2*NITEMS+4*NODES

         CALL setStorageSize(NITEMS,INFO)

   IF(ErrorFound()) RETURN 1

 

C       CHECK THE MINIMUM TIMESTEP

         IF (DELT.LT.DELT_MIN) THEN

      CALL MESSAGES(-1,MESSAGE1,'FATAL',IUNIT,ITYPE)

   ENDIF

 

C       RETURN TO THE CALLING PROGRAM

         RETURN 1

 

      ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INITIAL TIME

      IF (TIME.LT.(TIME0+DELT/2.D0)) THEN

 

C       SET THE UNIT NUMBER FOR FUTURE CALLS

         IUNIT=INFO(1)

 

C       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER

         LU_DATA=JFIX(VAL(1)+0.5)                  !DIMENSIONLESS,>10 = PARAMETERS READ FROM A FILE

         NODES=JFIX(VAL(2)+0.5)                    !DIMENSIONLESS

         XNODES=DBLE(NODES)                        !DIMENSIONLESS

         N_PORTS=JFIX(VAL(3)+0.5)                  !DIMENSIONLESS

         N_HX=JFIX(VAL(4)+0.1)  !DIMENSIONLESS

         N_MISC=JFIX(VAL(5)+0.5)                   !DIMENSIONLESS

   

C       READ THE PARAMETERS FROM THE DATA FILE

         IF(LU_DATA.GE.10) THEN

      REWIND(LU_DATA)

      READ(LU_DATA,*,END=100,ERR=9000) ADUM

      READ(LU_DATA,*,END=100,ERR=9000) ADUM

      DO J=1,NPAR

         READ(LU_DATA,*,END=100,ERR=100) VAL(5+J)

      ENDDO

   ENDIF

     

100    VOLTANK=VAL(6)                            !M

         IF(VOLTANK.LE.0.) CALL TYPECK(-4,INFO,0,6,0)

 

   TANKHEIGHT=VAL(7)                         !M

         IF(TANKHEIGHT.LE.0.) CALL TYPECK(-4,INFO,0,7,0)

 

   ITFLUID=JFIX(VAL(8)+0.1)                  !DIMENSIONLESS

C                                                  0=PROVIDED PARAMETERS

C                                                      1=WATER-VARIABLE

C                                                      2=PROPYLENE GLYCOL

C                                                      3=ETHYLENE GLYCOL

         IF((ITFLUID.LT.0).OR.(ITFLUID.GT.3)) CALL TYPECK(-4,INFO,0,8,0)

 

         IF(ITFLUID.EQ.0) THEN

      CP_TANK=VAL(9)                         !KJ/KG.K

      IF(CP_TANK.LE.0.) CALL TYPECK(-4,INFO,0,9,0)

 

      RHO_TANK=VAL(10)                        !KJ/KG.K

      IF(RHO_TANK.LE.0.) CALL TYPECK(-4,INFO,0,10,0)

 

      K_TANK=VAL(11)                          !KJ/KG.K

      IF(K_TANK.LE.0.) CALL TYPECK(-4,INFO,0,11,0)

 

      MU_TANK=VAL(12)                         !KJ/KG.K

      IF(MU_TANK.LE.0.) CALL TYPECK(-4,INFO,0,12,0)

 

      BETA_TANK=VAL(13)                       !KJ/KG.K

      IF(BETA_TANK.LE.0.) CALL TYPECK(-4,INFO,0,13,0)

 

            N_SPOT=13

 

         ELSE IF(ITFLUID.EQ.1) THEN

            N_SPOT=8

 

         ELSE IF(ITFLUID.EQ.2) THEN

      PERVOL_T=VAL(9)                        !%

            N_SPOT=9

 

            IF((PERVOL_T.LT.0.).OR.(PERVOL_T.GT.100.))

     1         CALL TYPECK(-4,INFO,0,9,0)

 

         ELSE IF(ITFLUID.EQ.3) THEN

      PERVOL_T=VAL(9)                        !%

            N_SPOT=9

 

            IF((PERVOL_T.LT.0.).OR.(PERVOL_T.GT.100.))

     1         CALL TYPECK(-4,INFO,0,9,0)

 

   ELSE

      CONTINUE

   ENDIF

 

         U_TOP=VAL(N_SPOT+1)                        !KJ/H.M2.K

         IF(U_TOP.LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+1,0)

         

   DO J=1,NODES

      U_EDGE(J)=VAL(N_SPOT+1+J)               !KJ/H.M2.K

            IF(U_EDGE(J).LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+1+J,0)

   ENDDO

 

         U_BOTTOM=VAL(N_SPOT+2+NODES)               !KJ/H.M2.K

         IF(U_BOTTOM.LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+3,0)

 

   K_ADDCOND=VAL(N_SPOT+3+NODES)              !KJ/H.M.K

 

         N_SPOT=N_SPOT+3+NODES

   IF(ErrorFound()) RETURN 1

 

         DO J=1,N_PORTS

            INLET_MODE(J)=JFIX(VAL(N_SPOT+1)+0.1)      !DIMENSIONLESS

C                                                      1=PROVIDE INLET AND OUTLET LOCATIONS

C                                                      2=PROVIDE FRACTION OF INLET FLOW TO NODES AND OUTLET LOCATIONS

C                                                      3=TEMPERATURE SEEKING MODE

 

 

            IF(INLET_MODE(J).EQ.1) THEN

 

C             READ THE LOAD FLOW PARAMETERS FROM THE INPUT FILE

         N_ENTRY(J)=JFIX(VAL(N_SPOT+2)+0.1)                  !DIMENSIONLESS

         N_EXIT(J)=JFIX(VAL(N_SPOT+3)+0.1)                   !DIMENSIONLESS

 

C             CHECK THE ENTRY LOCATIONS FOR PROBLEMS

         IF((N_ENTRY(J).LT.1).OR.(N_ENTRY(J).GT.NODES)) THEN

                  CALL TYPECK(-4,INFO,0,N_SPOT+2,0)

            RETURN 1

         ENDIF

 

         IF((N_EXIT(J).LT.1).OR.(N_EXIT(J).GT.NODES)) THEN

                  CALL TYPECK(-4,INFO,0,N_SPOT+3,0)

            RETURN 1

         ENDIF

 

C             SET THE FRACTION OF INLET FLOW TO 1 FOR THOSE NODES CONTAINING THE INLET

               DO K=1,NODES

            IF(K.EQ.N_ENTRY(J)) THEN

               FRAC_IN(J,K)=1.

            ELSE

               FRAC_IN(J,K)=0.

            ENDIF

         ENDDO

         N_SPOT=N_SPOT+3

 

      ELSE IF(INLET_MODE(J).EQ.2) THEN

         N_EXIT(J)=JFIX(VAL(N_SPOT+2)+0.1)                   !DIMENSIONLESS

 

         IF((N_EXIT(J).LT.1).OR.(N_EXIT(J).GT.NODES)) THEN

                  CALL TYPECK(-4,INFO,0,N_SPOT+2,0)

            RETURN 1

         ENDIF

 

C             GET THE FRACTION OF INLET FLOW TO EACH NODE

               FRAC_IN_TOT=0.

         DO K=1,NODES

            FRAC_IN(J,K)=VAL(N_SPOT+2+K)

            FRAC_IN_TOT=FRAC_IN_TOT+FRAC_IN(J,K)

            IF((FRAC_IN(J,K).LT.0.).OR.(FRAC_IN(J,K).GT.1.))

1               CALL TYPECK(-4,INFO,0,N_SPOT+2+K,0)

         ENDDO

         IF(FRAC_IN_TOT.LE.0.)

1            CALL TYPECK(-4,INFO,0,N_SPOT+2+NODES,0)

 

C             NORMALIZE THE FRACTIONS BACK TO A TOTAL OF 1

         DO K=1,NODES

            FRAC_IN(J,K)=FRAC_IN(J,K)/FRAC_IN_TOT

         ENDDO

 

         N_SPOT=N_SPOT+2+NODES

      

  ELSE

         N_EXIT(J)=JFIX(VAL(N_SPOT+2)+0.1)                   !DIMENSIONLESS

 

         IF((N_EXIT(J).LT.1).OR.(N_EXIT(J).GT.NODES)) THEN

                  CALL TYPECK(-4,INFO,0,N_SPOT+2,0)

            RETURN 1

         ENDIF

         N_SPOT=N_SPOT+2

      ENDIF

   ENDDO

   IF(ErrorFound()) RETURN 1

 

C       GET THE FLUE LOSS COEFFICIENT

         DO J=1,NODES

      UA_FLUE(J)=VAL(N_SPOT+J)

      IF(UA_FLUE(J).LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+J,0)

   ENDDO

   N_SPOT=N_SPOT+NODES

 

C       GET THE NUMBER OF MISCELLANEOUS HEAT FLOWS INTO/OUT OF THE TANK (IN=POSITIVE)

         DO J=1,N_MISC

      NODE_MISC(J)=JFIX(VAL(N_SPOT+J))

      IF((NODE_MISC(J).LT.1).OR.(NODE_MISC(J).GT.NODES)) THEN

         CALL TYPECK(-4,INFO,0,N_SPOT+J,0)

      ENDIF

   ENDDO

         N_SPOT=N_SPOT+N_MISC

 

C       RESET THE NUMBER OF HX NODES

   N_HXNODES=0

 

C       LOOP THROUGH EACH HEAT EXCHANGER AND RETRIEVE THE PARAMETERS

         DO K=1,N_HX

 

C          DETERMINE WHICH TYPE OF HEAT EXCHANGER IT IS

            IHXGEOM(K)=JFIX(VAL(N_SPOT+1)+0.1)           !DIMENSIONLESS,

C                                                         1=HORIZONTAL TUBE BANK

C                                                         2=VERTICAL TUBE BANK

C                                                         3=COILED HEAT EXCHANGER

C                                                         4=SERPENTINE TUBE

            IF((IHXGEOM(K).LT.1).OR.(IHXGEOM(K).GT.4)) THEN

         CALL TYPECK(-4,INFO,0,N_SPOT+1,0)

         RETURN 1

      ENDIF

 

C          RESET THE NUMBER OF TANK DIVISIONS BASED ON THE HX GEOMETRY

      IF(IHXGEOM(K).EQ.1) THEN

         N_TNKDIV(K)=NODES

      ELSE

         N_TNKDIV(K)=1

      ENDIF

 

      N_HXNODES(K)=JFIX(VAL(N_SPOT+2)+0.1)                      !DIMENSIONLESS

      IF((N_HXNODES(K).LT.1).OR.(N_HXNODES(K).GT.MAXHXNODES))

     1         CALL TYPECK(-4,INFO,0,N_SPOT+2,0)

 

C          GET THE REST OF THE HEAT EXCHANGER PARAMETERS

            IHXFLUID(K)=JFIX(VAL(N_SPOT+3)+0.1)                       !DIMENSIONLESS,

C                                                  0=PROVIDED PARAMETERS

C                                                      1=WATER-VARIABLE

C                                                      2=PROPYLENE GLYCOL

C                                                      3=ETHYLENE GLYCOL

            IF((IHXFLUID(K).LT.0).OR.(IHXFLUID(K).GT.3)) THEN

               CALL TYPECK(-4,INFO,0,N_SPOT+3,0)

         RETURN 1

      ENDIF

 

            IF(IHXFLUID(K).EQ.0) THEN

         CP_HX(K)=VAL(N_SPOT+4)                         !KJ/KG.K

         RHO_HX(K)=VAL(N_SPOT+5)                        !KJ/KG.K

         K_HX(K)=VAL(N_SPOT+6)                          !KJ/KG.K

         MU_HX(K)=VAL(N_SPOT+7)                         !KJ/KG.K

      

         IF(CP_HX(K).LE.0.) CALL TYPECK(-4,INFO,0,N_SPOT+4,0)

         IF(RHO_HX(K).LE.0.) CALL TYPECK(-4,INFO,0,N_SPOT+5,0)

         IF(K_HX(K).LE.0.) CALL TYPECK(-4,INFO,0,N_SPOT+6,0)

         IF(MU_HX(K).LE.0.) CALL TYPECK(-4,INFO,0,N_SPOT+7,0)

 

               N_SPOT=N_SPOT+7

 

            ELSE IF(IHXFLUID(K).EQ.1) THEN

               N_SPOT=N_SPOT+3

 

            ELSE IF(IHXFLUID(K).EQ.2) THEN

         PERVOL_HX(K)=VAL(N_SPOT+4)                        !%

               IF((PERVOL_HX(K).LT.0.).OR.(PERVOL_HX(K).GT.100.))

     1            CALL TYPECK(-4,INFO,0,N_SPOT+4,0)

               N_SPOT=N_SPOT+4

 

            ELSE IF(IHXFLUID(K).EQ.3) THEN

         PERVOL_HX(K)=VAL(N_SPOT+4)                        !%

               IF((PERVOL_HX(K).LT.0.).OR.(PERVOL_HX(K).GT.100.))

     1            CALL TYPECK(-4,INFO,0,N_SPOT+4,0)

               N_SPOT=N_SPOT+4

 

      ELSE

         CONTINUE

      ENDIF

 

      C(K)=VAL(N_SPOT+1)                                         !DIMENSIONLESS

      RA_EXP(K)=VAL(N_SPOT+2)                                    !DIMENSIONLESS

      GF(K)=VAL(N_SPOT+3)                                        !DIMENSIONLESS

      GF_EXP(K)=VAL(N_SPOT+4)                                    !DIMENSIONLESS

            IN_DIA(K)=VAL(N_SPOT+5)                                    !M

      OUT_DIA(K)=VAL(N_SPOT+6)                                  !M

      K_WALL(K)=VAL(N_SPOT+7)                                   !KJ/H.M.K

      L_TUBE(K)=VAL(N_SPOT+8)                                   !M

      N_TUBES(K)=DBLE(JFIX(VAL(N_SPOT+9)+0.1))                  !DIMENSIONLESS

      HEADERVOL(K)=VAL(N_SPOT+10)                                !M3

            ACROSS_HX(K)=VAL(N_SPOT+11)                                !M2

 

            IF(C(K).LT.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+1,0)

            IF(GF(K).LT.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+3,0)

            IF(IN_DIA(K).LE.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+5,0)

            IF(OUT_DIA(K).LE.IN_DIA(K))

     1         CALL TYPECK(-4,INFO,0,N_SPOT+6,0)

            IF(K_WALL(K).LE.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+7,0)

            IF(L_TUBE(K).LE.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+8,0)

            IF(N_TUBES(K).LT.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+9,0)

            IF(HEADERVOL(K).LT.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+10,0)

            IF(ACROSS_HX(K).LT.0.)

     1         CALL TYPECK(-4,INFO,0,N_SPOT+11,0)

 

C          READ IN THE COIL DIAMETER AND PITCH IF ITS A COILED HEAT EXCHANGER

            IF(IHXGEOM(K).EQ.3) THEN

        DIA_COIL(K)=VAL(N_SPOT+12)

        PITCH_COIL(K)=VAL(N_SPOT+13)

            

              IF(DIA_COIL(K).LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+12,0)

              IF(PITCH_COIL(K).LT.0.) CALL TYPECK(-4,INFO,0,N_SPOT+13,0)

 

        N_SPOT=N_SPOT+13

      ELSE

        N_SPOT=N_SPOT+11

      ENDIF

 

C          GET THE FRACTION OF THE TUBES THAT RESIDE IN EACH OF THE TANK NODES

            IF(IHXGEOM(K).EQ.1) THEN

               TOTFRAC=0.

               DO I=1,NODES

            FRACHX(I,K)=VAL(N_SPOT+I)                  !DIMENSIONLESS

                  TOTFRAC=TOTFRAC+FRACHX(I,K)

            IF((FRACHX(I,K).LT.0.).OR.(TOTFRAC.GT.1.00001))

     1               CALL TYPECK(-4,INFO,0,N_SPOT+I,0)

               ENDDO

         N_SPOT=N_SPOT+NODES

 

C          GET THE ADJACENT TANK NODE AND FRACTION OF HEAT EXCHANGER LENGTH IN THAT TANK NODE

            ELSE

               TOTFRAC=0.

               DO J=1,N_HXNODES(K)

            N_HXFLOW(J,K)=JFIX(VAL(N_SPOT+(2*J)-1)+0.1)    !DIMENSIONLESS

            FRACHX(J,K)=VAL(N_SPOT+(2*J))                  !DIMENSIONLESS

                  TOTFRAC=TOTFRAC+FRACHX(J,K)

            IF((N_HXFLOW(J,K).LT.1).OR.(N_HXFLOW(J,K).GT.NODES))

     1               CALL TYPECK(-4,INFO,0,N_SPOT+(2*J)-1,0)

            IF((FRACHX(J,K).LT.0.).OR.(TOTFRAC.GT.1.00001))

     1               CALL TYPECK(-4,INFO,0,N_SPOT+(2*J),0)

            IF(J.GT.1) THEN

               IF(ABS(N_HXFLOW(J,K)-N_HXFLOW(J-1,K)).GT.1)

     1                  CALL TYPECK(-4,INFO,0,N_SPOT+(2*J)-1,0)

            ENDIF

               ENDDO

         N_SPOT=N_SPOT+2*N_HXNODES(K)

            ENDIF

 

      IF(ERRORFOUND()) RETURN 1

 

C          DETERMINE THE NUMBER OF TIMES THE INTERFACE BETWEEN TANK NODES IS BROKEN BY A HX PASS

            DO J=1,NODES

              IPASSES(J,K)=0

            ENDDO

 

            IF(IHXGEOM(K).EQ.1) THEN

               DO J=1,NODES-1

            IF((FRACHX(J,K).GT.0.).AND.(FRACHX(J+1,K).GT.0.)) THEN

               IPASSES(J,K)=IPASSES(J,K)+1

            ENDIF

               ENDDO

            ELSE

               DO J=1,N_HXNODES(K)-1

            IF(N_HXFLOW(J,K).LT.N_HXFLOW(J+1,K)) THEN

               IPASSES(N_HXFLOW(J,K),K)=IPASSES(N_HXFLOW(J,K),K)+1

            ELSE IF(N_HXFLOW(J,K).GT.N_HXFLOW(J+1,K)) THEN

               IPASSES(N_HXFLOW(J,K)-1,K)=

1                  IPASSES(N_HXFLOW(J,K)-1,K)+1

            ELSE

               IPASSES(N_HXFLOW(J,K),K)=IPASSES(N_HXFLOW(J,K),K)

            ENDIF

               ENDDO

            ENDIF

 

C          SET SOME INITIAL CONDITIONS

      OA_TUBE(K)=PI*OUT_DIA(K)*L_TUBE(K)

      IA_TUBE(K)=PI*IN_DIA(K)*L_TUBE(K)

      VOLTUBE_O(K)=PI*OUT_DIA(K)*OUT_DIA(K)*L_TUBE(K)/4.

      VOLTUBE_I(K)=PI*IN_DIA(K)*IN_DIA(K)*L_TUBE(K)/4.

      OUT_RAD(K)=OUT_DIA(K)/2.

      IN_RAD(K)=IN_DIA(K)/2.

   ENDDO

 

C       SEE IF CONDUCTION EFFECTS SHOULD BE INCLUDED

   IF(K_ADDCOND.GE.0.) THEN

            DONT_COND=.FALSE.

   ELSE

            DONT_COND=.TRUE.

   ENDIF

 

C       TURN OFF CONDUCTION AND MIXING IF THERE IS ONLY 1 NODE (FULLY MIXED)

         IF(NODES.EQ.1) THEN

      DONT_COND=.TRUE.

      DONT_MIX=.TRUE.

         ENDIF

 

C       CALCULATE THE NOMINAL VOLUME OF EACH TANK AND HEAT EXCHANGER NODE

         DO J=1,NODES

      TANKVOL(J)=VOLTANK/DBLE(NODES)

         ENDDO

 

         DO K=1,N_HX

      IF(IHXGEOM(K).EQ.1) THEN

         DO I=1,N_TNKDIV(K)

                  DO J=1,N_HXNODES(K)

               IF((J.EQ.1).OR.(J.EQ.N_HXNODES(K))) THEN

                        HXVOL_O(I,J,K)=HEADERVOL(K)*FRACHX(I,K)+

     1                     VOLTUBE_O(K)/DBLE(N_HXNODES(K))*N_TUBES(K)*

     1                     FRACHX(I,K)

                        HXVOL_I(I,J,K)=HEADERVOL(K)*FRACHX(I,K)+

     1                     VOLTUBE_I(K)/DBLE(N_HXNODES(K))*N_TUBES(K)*

     1                     FRACHX(I,K)

               ELSE

                        HXVOL_O(I,J,K)=VOLTUBE_O(K)/DBLE(N_HXNODES(K))*

1                     N_TUBES(K)*FRACHX(I,K)

                        HXVOL_I(I,J,K)=VOLTUBE_I(K)/DBLE(N_HXNODES(K))*

1                     N_TUBES(K)*FRACHX(I,K)

               ENDIF

               TANKVOL(I)=TANKVOL(I)-HXVOL_O(I,J,K)

                  ENDDO

               ENDDO

            ELSE IF((IHXGEOM(K).GT.1).AND.(IHXGEOM(K).LE.4)) THEN

         DO I=1,N_TNKDIV(K)

                  DO J=1,N_HXNODES(K)

               IF((J.EQ.1).OR.(J.EQ.N_HXNODES(K))) THEN

                        HXVOL_O(I,J,K)=HEADERVOL(K)+VOLTUBE_O(K)*

1                     N_TUBES(K)*FRACHX(J,K)

                        HXVOL_I(I,J,K)=HEADERVOL(K)+VOLTUBE_I(K)*

1                     N_TUBES(K)*FRACHX(J,K)

               ELSE

                        HXVOL_O(I,J,K)=VOLTUBE_O(K)*N_TUBES(K)*

1                     FRACHX(J,K)

                        HXVOL_I(I,J,K)=VOLTUBE_I(K)*N_TUBES(K)*

1                     FRACHX(J,K)

               ENDIF

               TANKVOL(N_HXFLOW(J,K))=TANKVOL(N_HXFLOW(J,K))-

1                  HXVOL_O(I,J,K)

                  ENDDO

               ENDDO

            ENDIF

   ENDDO

 

C       CHECK TO MAKE SURE THAT THE TANK VOLUME IS GREATER THAN ZERO

         DO J=1,NODES

      IF(TANKVOL(J).LE.0.) CALL TYPECK(-4,INFO,0,2,0)

         ENDDO

 

C       CALCULATE THE EDGE, TOP, BOTTOM, AND COMMON AREAS FOR EACH OF THE TANK NODES

C       THIS MODEL ASSUMES THAT THE NODES ARE OF EQUAL VOLUME (BEFORE THE

C       HEAT EXCHANGER VOLUME REDUCTION).

 

C       CALCULATE THE TANK RADIUS

         RADIUS=(VOLTANK/PI/TANKHEIGHT)**0.5

 

C       IF THERE IS ONLY ONE NODE, THE AREA CALCULATIONS ARE SIMPLE

         IF(NODES.EQ.1) THEN

      ATOP(1)=PI*RADIUS*RADIUS

      ABOT(1)=PI*RADIUS*RADIUS

      AEDGE(1)=2.*PI*RADIUS*TANKHEIGHT

            ACOND(1)=PI*RADIUS*RADIUS

            L_COND(1)=TANKHEIGHT/2.

      GOTO 126

   ENDIF

 

C       CALCULATE THE AREAS FOR A MULTI-NODE TANK

         DO J=1,NODES

      IF(J.EQ.1) THEN

         ATOP(J)=PI*RADIUS*RADIUS

         ABOT(J)=0.

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ELSE IF(J.EQ.NODES) THEN

         ATOP(J)=0.

         ABOT(J)=PI*RADIUS*RADIUS

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ELSE

         ATOP(J)=0.

         ABOT(J)=0.

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ENDIF

         ENDDO

 

C       MODIFY THE CONDUCTION AREA TO ACCOUNT FOR HX PASSES

   DO K=1,N_HX

      IF(IHXGEOM(K).GT.0) THEN

         DO J=1,NODES

            ACOND(J)=ACOND(J)-DBLE(IPASSES(J,K))*ACROSS_HX(K)

            IF(ACOND(J).LT.0.) CALL TYPECK(-4,INFO,NI,2,0)

               ENDDO

      ENDIF

   ENDDO

 

126      CONTINUE

 

C       INITIALIZE THE TEMPERATURES

         DO J=1,NODES

      TI_TANK(J)=T(J)

      TF_TANK(J)=T(J)

      TAVE_TANK(J)=T(J)

      FLOWF_MIX(J)=0.

      FLOWI_MIX(J)=0.

         ENDDO

 

         DO K=1,N_HX

      DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

            IF(IHXGEOM(K).EQ.1) THEN

               TINIT_HX(K)=T(I)

            ELSE

               TINIT_HX(K)=T(N_HXFLOW(J,K))

            ENDIF

            TI_HX(I,J,K)=TINIT_HX(K)

            TF_HX(I,J,K)=TINIT_HX(K)

            TAVE_HX(I,J,K)=TINIT_HX(K)

            WALL_RESIST(I,J,K)=9.9999D+20

            RESIST_OUT(I,J,K)=9.9999D+20

            RESIST_IN(I,J,K)=9.9999D+20

               ENDDO

            ENDDO

         ENDDO

 

C       SET THE AVERAGE TEMPERATURE OF THE STORAGE TANK

         AVETEMP_T=0.

         DO J=1,NODES

      AVETEMP_T=AVETEMP_T+TAVE_TANK(J)

         ENDDO

         AVETEMP_T=AVETEMP_T/XNODES

 

C        SET THE OUTLET AND AVERAGE TEMPERATURES FOR THE HX

   AVETEMP_HX=0.

   TEMPHX_O=0.

   DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

            IF(IHXGEOM(K).EQ.1) THEN

               AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*

1                  FRACHX(I,K)/DBLE(N_HXNODES(K))

            ELSE

               AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*

1                  FRACHX(J,K)

            ENDIF

               ENDDO

               

         IF(IHXGEOM(K).EQ.1) THEN

            TEMPHX_O(K)=TEMPHX_O(K)+TAVE_HX(I,N_HXNODES(K),K)*

1               FRACHX(I,K)

         ELSE

            TEMPHX_O(K)=TAVE_HX(I,N_HXNODES(K),K)

         ENDIF

            ENDDO

         ENDDO

 

C       SET THE INITIAL VALUES OF THE OUTPUTS HERE

         DO J=1,N_PORTS

            OUT(2*J-1)=TAVE_TANK(N_EXIT(J))

            OUT(2*J)=0.

   ENDDO

 

   OUT(2*N_PORTS+1)=AVETEMP_T

   OUT(2*N_PORTS+2)=0.

 

         DO J=1,N_PORTS

      OUT(2*N_PORTS+2+J)=0.

   ENDDO

 

   OUT(3*N_PORTS+3)=0.

   OUT(3*N_PORTS+4)=0.

   OUT(3*N_PORTS+5)=0.

   OUT(3*N_PORTS+6)=0.

   OUT(3*N_PORTS+7)=0.

   OUT(3*N_PORTS+8)=0.

   OUT(3*N_PORTS+9)=0.

   OUT(3*N_PORTS+10)=0.

   OUT(3*N_PORTS+11)=0.

 

   DO J=1,NODES

      OUT(3*N_PORTS+11+J)=TAVE_TANK(J)

         ENDDO

 

   DO J=1,N_HX

      OUT(3*N_PORTS+NODES+11+2*J-1)=TEMPHX_O(J)

      OUT(3*N_PORTS+NODES+11+2*J)=0.

   ENDDO

 

   DO J=1,N_HX

      OUT(3*N_PORTS+NODES+11+2*N_HX+J)=AVETEMP_HX(J)

   ENDDO

 

   DO J=1,N_HX

      OUT(3*N_PORTS+11+NODES+3*N_HX+1+4*(J-1))=0.

      OUT(3*N_PORTS+11+NODES+3*N_HX+2+4*(J-1))=0.

      OUT(3*N_PORTS+11+NODES+3*N_HX+3+4*(J-1))=0.

      OUT(3*N_PORTS+11+NODES+3*N_HX+4+4*(J-1))=0.

   ENDDO

 

C       SET THE INITIAL STORAGE VARIABLES HERE

         DO J=1,NODES

      STORED(J)=TI_TANK(J)

      STORED(NODES+J)=TF_TANK(J)

      STORED(2*NODES+J)=FLOWI_MIX(J)

      STORED(3*NODES+J)=FLOWF_MIX(J)

         ENDDO

 

C       STORE THE HX TEMPERATURES IN THE STORAGE ARRAY AT THE FIRST ITERATION

         NITEMS=0

         DO K=1,N_HX

            NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

   ENDDO

   NITEMS=2*NITEMS+4*NODES

 

   N_STORED=4*NODES

 

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           STORED(N_STORED+1)=TI_HX(I,J,K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           STORED(N_STORED+1)=TF_HX(I,J,K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

C       PUT THE STORED ARRAY IN THE GLOBAL STORED ARRAY

   CALL setStorageVars(STORED,NITEMS,INFO)

 

C       RETURN TO THE CALLING PROGRAM

         RETURN 1

 

      ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    *** ITS AN ITERATIVE CALL TO THIS COMPONENT ***

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    RE-READ THE PARAMETERS IF ANOTHER UNIT OF THIS TYPE HAS BEEN CALLED

      IF(INFO(1).NE.IUNIT) THEN

 

C       RESET THE UNIT NUMBER

   IUNIT=INFO(1)

   ITYPE=INFO(2)

    

C       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER

         LU_DATA=JFIX(VAL(1)+0.5)                  !DIMENSIONLESS,>10 = PARAMETERS READ FROM A FILE

         NODES=JFIX(VAL(2)+0.5)                    !DIMENSIONLESS

         XNODES=DBLE(NODES)                        !DIMENSIONLESS

         N_PORTS=JFIX(VAL(3)+0.5)                  !DIMENSIONLESS

         N_HX=JFIX(VAL(4)+0.1)  !DIMENSIONLESS

         N_MISC=JFIX(VAL(5)+0.5)                   !DIMENSIONLESS

   

C       READ THE PARAMETERS FROM THE DATA FILE

         IF(LU_DATA.GE.10) THEN

      REWIND(LU_DATA)

      READ(LU_DATA,*,END=30,ERR=9000) ADUM

      READ(LU_DATA,*,END=30,ERR=9000) ADUM

      DO J=1,NPAR

         READ(LU_DATA,*,END=30,ERR=30) VAL(5+J)

      ENDDO

   ENDIF

     

30    VOLTANK=VAL(6)                             !M

   TANKHEIGHT=VAL(7)                          !M

   ITFLUID=JFIX(VAL(8)+0.1)                   !DIMENSIONLESS

         IF(ITFLUID.EQ.0) THEN

      CP_TANK=VAL(9)                          !KJ/KG.K

      RHO_TANK=VAL(10)                        !KJ/KG.K

      K_TANK=VAL(11)                          !KJ/KG.K

      MU_TANK=VAL(12)                         !KJ/KG.K

      BETA_TANK=VAL(13)                       !KJ/KG.K

            N_SPOT=13

         ELSE IF(ITFLUID.EQ.1) THEN

            N_SPOT=8

         ELSE IF(ITFLUID.EQ.2) THEN

      PERVOL_T=VAL(9)                         !%

            N_SPOT=9

         ELSE IF(ITFLUID.EQ.3) THEN

      PERVOL_T=VAL(9)                         !%

            N_SPOT=9

   ELSE

      CONTINUE

   ENDIF

         U_TOP=VAL(N_SPOT+1)                        !KJ/H.M2.K

   DO J=1,NODES

      U_EDGE(J)=VAL(N_SPOT+1+J)               !KJ/H.M2.K

   ENDDO

         U_BOTTOM=VAL(N_SPOT+2+NODES)               !KJ/H.M2.K

   K_ADDCOND=VAL(N_SPOT+3+NODES)              !KJ/H.M.K

         N_SPOT=N_SPOT+3+NODES

 

         DO J=1,N_PORTS

            INLET_MODE(J)=JFIX(VAL(N_SPOT+1)+0.1)

            IF(INLET_MODE(J).EQ.1) THEN

         N_ENTRY(J)=JFIX(VAL(N_SPOT+2)+0.1)                  !DIMENSIONLESS

         N_EXIT(J)=JFIX(VAL(N_SPOT+3)+0.1)                   !DIMENSIONLESS

               DO K=1,NODES

            IF(K.EQ.N_ENTRY(J)) THEN

               FRAC_IN(J,K)=1.

            ELSE

               FRAC_IN(J,K)=0.

            ENDIF

         ENDDO

         N_SPOT=N_SPOT+3

      ELSE IF(INLET_MODE(J).EQ.2) THEN

         N_EXIT(J)=JFIX(VAL(N_SPOT+2)+0.1)                   !DIMENSIONLESS

               FRAC_IN_TOT=0.

         DO K=1,NODES

            FRAC_IN(J,K)=VAL(N_SPOT+2+K)

            FRAC_IN_TOT=FRAC_IN_TOT+FRAC_IN(J,K)

         ENDDO

         DO K=1,NODES

            FRAC_IN(J,K)=FRAC_IN(J,K)/FRAC_IN_TOT

         ENDDO

         N_SPOT=N_SPOT+2+NODES

  ELSE

         N_EXIT(J)=JFIX(VAL(N_SPOT+2)+0.1)                   !DIMENSIONLESS

         N_SPOT=N_SPOT+2

      ENDIF

   ENDDO

 

C       GET THE FLUE LOSS COEFFICIENT

         DO J=1,NODES

      UA_FLUE(J)=VAL(N_SPOT+J)

   ENDDO

   N_SPOT=N_SPOT+NODES

 

C       GET THE NUMBER OF MISCELLANEOUS HEAT FLOWS INTO/OUT OF THE TANK (IN=POSITIVE)

         DO J=1,N_MISC

      NODE_MISC(J)=JFIX(VAL(N_SPOT+J))

   ENDDO

         N_SPOT=N_SPOT+N_MISC

 

C       RESET THE NUMBER OF HX NODES

   N_HXNODES=0

 

C       LOOP THROUGH EACH HEAT EXCHANGER AND RETRIEVE THE PARAMETERS

         DO K=1,N_HX

 

C          DETERMINE WHICH TYPE OF HEAT EXCHANGER IT IS

            IHXGEOM(K)=JFIX(VAL(N_SPOT+1)+0.1)           !DIMENSIONLESS,

 

C          RESET THE NUMBER OF TANK DIVISIONS BASED ON THE HX GEOMETRY

      IF(IHXGEOM(K).EQ.1) THEN

         N_TNKDIV(K)=NODES

      ELSE

         N_TNKDIV(K)=1

      ENDIF

 

      N_HXNODES(K)=JFIX(VAL(N_SPOT+2)+0.1)                      !DIMENSIONLESS

 

C          GET THE REST OF THE HEAT EXCHANGER PARAMETERS

            IHXFLUID(K)=JFIX(VAL(N_SPOT+3)+0.1)                       !DIMENSIONLESS,

            IF(IHXFLUID(K).EQ.0) THEN

         CP_HX(K)=VAL(N_SPOT+4)                         !KJ/KG.K

         RHO_HX(K)=VAL(N_SPOT+5)                        !KJ/KG.K

         K_HX(K)=VAL(N_SPOT+6)                          !KJ/KG.K

         MU_HX(K)=VAL(N_SPOT+7)                         !KJ/KG.K

               N_SPOT=N_SPOT+7

            ELSE IF(IHXFLUID(K).EQ.1) THEN

               N_SPOT=N_SPOT+3

            ELSE IF(IHXFLUID(K).EQ.2) THEN

         PERVOL_HX(K)=VAL(N_SPOT+4)                        !%

               N_SPOT=N_SPOT+4

            ELSE IF(IHXFLUID(K).EQ.3) THEN

         PERVOL_HX(K)=VAL(N_SPOT+4)                        !%

               N_SPOT=N_SPOT+4

      ENDIF

 

      C(K)=VAL(N_SPOT+1)                                         !DIMENSIONLESS

      RA_EXP(K)=VAL(N_SPOT+2)                                    !DIMENSIONLESS

      GF(K)=VAL(N_SPOT+3)                                        !DIMENSIONLESS

      GF_EXP(K)=VAL(N_SPOT+4)                                    !DIMENSIONLESS

            IN_DIA(K)=VAL(N_SPOT+5)                                    !M

      OUT_DIA(K)=VAL(N_SPOT+6)                                  !M

      K_WALL(K)=VAL(N_SPOT+7)                                   !KJ/H.M.K

      L_TUBE(K)=VAL(N_SPOT+8)                                   !M

      N_TUBES(K)=DBLE(JFIX(VAL(N_SPOT+9)+0.1))                  !DIMENSIONLESS

      HEADERVOL(K)=VAL(N_SPOT+10)                                !M3

            ACROSS_HX(K)=VAL(N_SPOT+11)                                !M2

 

C          READ IN THE COIL DIAMETER AND PITCH IF ITS A COILED HEAT EXCHANGER

            IF(IHXGEOM(K).EQ.3) THEN

        DIA_COIL(K)=VAL(N_SPOT+12)

        PITCH_COIL(K)=VAL(N_SPOT+13)

        N_SPOT=N_SPOT+13

      ELSE

        N_SPOT=N_SPOT+11

      ENDIF

 

C          GET THE FRACTION OF THE TUBES THAT RESIDE IN EACH OF THE TANK NODES

            IF(IHXGEOM(K).EQ.1) THEN

               TOTFRAC=0.

               DO I=1,NODES

            FRACHX(I,K)=VAL(N_SPOT+I)                  !DIMENSIONLESS

                  TOTFRAC=TOTFRAC+FRACHX(I,K)

               ENDDO

         N_SPOT=N_SPOT+NODES

 

C          GET THE ADJACENT TANK NODE AND FRACTION OF HEAT EXCHANGER LENGTH IN THAT TANK NODE

            ELSE

               TOTFRAC=0.

               DO J=1,N_HXNODES(K)

            N_HXFLOW(J,K)=JFIX(VAL(N_SPOT+(2*J)-1)+0.1)    !DIMENSIONLESS

            FRACHX(J,K)=VAL(N_SPOT+(2*J))                  !DIMENSIONLESS

                  TOTFRAC=TOTFRAC+FRACHX(J,K)

               ENDDO

         N_SPOT=N_SPOT+2*N_HXNODES(K)

            ENDIF

 

C          DETERMINE THE NUMBER OF TIMES THE INTERFACE BETWEEN TANK NODES IS BROKEN BY A HX PASS

            DO J=1,NODES

              IPASSES(J,K)=0

            ENDDO

 

            IF(IHXGEOM(K).EQ.1) THEN

               DO J=1,NODES-1

            IF((FRACHX(J,K).GT.0.).AND.(FRACHX(J+1,K).GT.0.)) THEN

               IPASSES(J,K)=IPASSES(J,K)+1

            ENDIF

               ENDDO

            ELSE

               DO J=1,N_HXNODES(K)-1

            IF(N_HXFLOW(J,K).LT.N_HXFLOW(J+1,K)) THEN

               IPASSES(N_HXFLOW(J,K),K)=IPASSES(N_HXFLOW(J,K),K)+1

            ELSE IF(N_HXFLOW(J,K).GT.N_HXFLOW(J+1,K)) THEN

               IPASSES(N_HXFLOW(J,K)-1,K)=

1                  IPASSES(N_HXFLOW(J,K)-1,K)+1

            ELSE

               IPASSES(N_HXFLOW(J,K),K)=IPASSES(N_HXFLOW(J,K),K)

            ENDIF

               ENDDO

            ENDIF

 

C          SET SOME INITIAL CONDITIONS

      OA_TUBE(K)=PI*OUT_DIA(K)*L_TUBE(K)

      IA_TUBE(K)=PI*IN_DIA(K)*L_TUBE(K)

      VOLTUBE_O(K)=PI*OUT_DIA(K)*OUT_DIA(K)*L_TUBE(K)/4.

      VOLTUBE_I(K)=PI*IN_DIA(K)*IN_DIA(K)*L_TUBE(K)/4.

      OUT_RAD(K)=OUT_DIA(K)/2.

      IN_RAD(K)=IN_DIA(K)/2.

   ENDDO

 

C       SEE IF CONDUCTION EFFECTS SHOULD BE INCLUDED

   IF(K_ADDCOND.GE.0.) THEN

            DONT_COND=.FALSE.

   ELSE

            DONT_COND=.TRUE.

   ENDIF

 

C       TURN OFF CONDUCTION AND MIXING IF THERE IS ONLY 1 NODE (FULLY MIXED)

         IF(NODES.EQ.1) THEN

      DONT_COND=.TRUE.

      DONT_MIX=.TRUE.

         ENDIF

 

C       CALCULATE THE NOMINAL VOLUME OF EACH TANK AND HEAT EXCHANGER NODE

         DO J=1,NODES

      TANKVOL(J)=VOLTANK/DBLE(NODES)

         ENDDO

 

         DO K=1,N_HX

      IF(IHXGEOM(K).EQ.1) THEN

         DO I=1,N_TNKDIV(K)

                  DO J=1,N_HXNODES(K)

               IF((J.EQ.1).OR.(J.EQ.N_HXNODES(K))) THEN

                        HXVOL_O(I,J,K)=HEADERVOL(K)*FRACHX(I,K)+

     1                     VOLTUBE_O(K)/DBLE(N_HXNODES(K))*N_TUBES(K)*

     1                     FRACHX(I,K)

                        HXVOL_I(I,J,K)=HEADERVOL(K)*FRACHX(I,K)+

     1                     VOLTUBE_I(K)/DBLE(N_HXNODES(K))*N_TUBES(K)*

     1                     FRACHX(I,K)

               ELSE

                        HXVOL_O(I,J,K)=VOLTUBE_O(K)/DBLE(N_HXNODES(K))*

1                     N_TUBES(K)*FRACHX(I,K)

                        HXVOL_I(I,J,K)=VOLTUBE_I(K)/DBLE(N_HXNODES(K))*

1                     N_TUBES(K)*FRACHX(I,K)

               ENDIF

               TANKVOL(I)=TANKVOL(I)-HXVOL_O(I,J,K)

                  ENDDO

               ENDDO

            ELSE IF((IHXGEOM(K).GT.1).AND.(IHXGEOM(K).LE.4)) THEN

         DO I=1,N_TNKDIV(K)

                  DO J=1,N_HXNODES(K)

               IF((J.EQ.1).OR.(J.EQ.N_HXNODES(K))) THEN

                        HXVOL_O(I,J,K)=HEADERVOL(K)+VOLTUBE_O(K)*

1                     N_TUBES(K)*FRACHX(J,K)

                        HXVOL_I(I,J,K)=HEADERVOL(K)+VOLTUBE_I(K)*

1                     N_TUBES(K)*FRACHX(J,K)

               ELSE

                        HXVOL_O(I,J,K)=VOLTUBE_O(K)*N_TUBES(K)*

1                     FRACHX(J,K)

                        HXVOL_I(I,J,K)=VOLTUBE_I(K)*N_TUBES(K)*

1                     FRACHX(J,K)

               ENDIF

               TANKVOL(N_HXFLOW(J,K))=TANKVOL(N_HXFLOW(J,K))-

1                  HXVOL_O(I,J,K)

                  ENDDO

               ENDDO

            ENDIF

   ENDDO

 

C       CHECK TO MAKE SURE THAT THE TANK VOLUME IS GREATER THAN ZERO

         DO J=1,NODES

      IF(TANKVOL(J).LE.0.) CALL TYPECK(-4,INFO,0,2,0)

         ENDDO

 

C       CALCULATE THE EDGE, TOP, BOTTOM, AND COMMON AREAS FOR EACH OF THE TANK NODES

C       THIS MODEL ASSUMES THAT THE NODES ARE OF EQUAL VOLUME (BEFORE THE

C       HEAT EXCHANGER VOLUME REDUCTION).

 

C       CALCULATE THE TANK RADIUS

         RADIUS=(VOLTANK/PI/TANKHEIGHT)**0.5

 

C       IF THERE IS ONLY ONE NODE, THE AREA CALCULATIONS ARE SIMPLE

         IF(NODES.EQ.1) THEN

      ATOP(1)=PI*RADIUS*RADIUS

      ABOT(1)=PI*RADIUS*RADIUS

      AEDGE(1)=2.*PI*RADIUS*TANKHEIGHT

            ACOND(1)=PI*RADIUS*RADIUS

            L_COND(1)=TANKHEIGHT/2.

      GOTO 136

   ENDIF

 

C       CALCULATE THE AREAS FOR A MULTI-NODE TANK

         DO J=1,NODES

      IF(J.EQ.1) THEN

         ATOP(J)=PI*RADIUS*RADIUS

         ABOT(J)=0.

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ELSE IF(J.EQ.NODES) THEN

         ATOP(J)=0.

         ABOT(J)=PI*RADIUS*RADIUS

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ELSE

         ATOP(J)=0.

         ABOT(J)=0.

         AEDGE(J)=2.*PI*RADIUS*TANKHEIGHT/DBLE(NODES)

         ACOND(J)=PI*RADIUS*RADIUS

         L_COND(J)=TANKHEIGHT/DBLE(NODES)

      ENDIF

         ENDDO

 

C       MODIFY THE CONDUCTION AREA TO ACCOUNT FOR HX PASSES

   DO K=1,N_HX

      IF(IHXGEOM(K).GT.0) THEN

         DO J=1,NODES

            ACOND(J)=ACOND(J)-DBLE(IPASSES(J,K))*ACROSS_HX(K)

               ENDDO

      ENDIF

   ENDDO

 

136      CONTINUE

 

      ENDIF

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    RETRIEVE THE VALUES IN THE STORAGE ARRAY FOR THIS ITERATION

      NITEMS=0

      DO K=1,N_HX

         NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

ENDDO

NITEMS=2*NITEMS+4*NODES

 

      CALL getStorageVars(STORED,NITEMS,INFO)

IF(ErrorFound()) RETURN 1

 

      DO J=1,NODES

   TI_TANK(J)=STORED(J)

   TSTART_TANK(J)=STORED(J)

   TF_TANK(J)=STORED(NODES+J)

   FLOWI_MIX(J)=STORED(2*NODES+J)

   FLOWF_MIX(J)=STORED(3*NODES+J)

      ENDDO

 

      N_STORED=4*NODES

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

        TI_HX(I,J,K)=STORED(N_STORED+1)

        TSTART_HX(I,J,K)=STORED(N_STORED+1)

        N_STORED=N_STORED+1

            ENDDO

         ENDDO

      ENDDO

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

        TF_HX(I,J,K)=STORED(N_STORED+1)

        N_STORED=N_STORED+1

            ENDDO

         ENDDO

      ENDDO

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    RETRIEVE THE CURRENT VALUES OF THE INPUTS TO THIS MODEL FROM THE XIN ARRAY IN SEQUENTIAL ORDER

      DO J=1,N_PORTS

   TFLOW_IN(J)=XIN(2*J-1)                     !C

         FLOW_IN(J)=XIN(2*J)                        !KG/HR

         IF(FLOW_IN(J).LT.0.) CALL TYPECK(-3,INFO,2*J,0,0)

ENDDO

 

      DO K=1,N_HX

         THXFLOW_IN(K)=XIN(2*N_PORTS+2*(K-1)+1)     !C

         FLOWHX_IN(K)=XIN(2*N_PORTS+2*(K-1)+2)      !KG/HR

         IF(FLOWHX_IN(K).LT.0.)

1      CALL TYPECK(-3,INFO,2*N_PORTS+2*(K-1)+2,0,0)

      ENDDO


T_TLOSS=XIN(2*N_PORTS+2*N_HX+1)               !C

      

DO J=1,NODES

   T_ELOSS(J)=XIN(2*N_PORTS+2*N_HX+1+J)           !C

ENDDO

N_SPOT=2*N_PORTS+2*N_HX+1+NODES

      

T_BLOSS=XIN(N_SPOT+1)                         !C

      T_FLUE=XIN(N_SPOT+2)                          !KG/HR

      MIXFLOW=XIN(N_SPOT+3)                         !KG/HR

      DO J=1,NODES

   Q_AUX(J)=XIN(N_SPOT+3+J)                   !KJ/HR

      ENDDO

      DO J=1,N_MISC

   Q_MISC(J)=XIN(N_SPOT+3+NODES+J)            !KJ/HR

      ENDDO

      

C    TURN OFF MIXING IF THE MIXING FLOW IS ZER0 OR LESS

      IF(MIXFLOW.LE.0.) DONT_MIX=.TRUE.

IF((MIXFLOW.GT.0.).AND.(NODES.GT.1)) DONT_MIX=.FALSE.

 

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    CHECK THE INPUTS FOR PROBLEMS

IF(ErrorFound()) RETURN 1

C-----------------------------------------------------------------------------------------------------------------------

 

C-----------------------------------------------------------------------------------------------------------------------

C    PERFORM ALL THE CALCULATION HERE FOR THIS MODEL.

MAXITERS=20

 

C    SET THE INITIAL MIXING FLOW TO WHAT IT WAS AT THE END OF THE PREVIOUS TIMESTEP IF ITS STILL NEEDED, OTHERWISE TURN IT OFF

      DO J=1,NODES-1

         IF((TSTART_TANK(J)+CONVERGED).GE.TSTART_TANK(J+1)) THEN

      FLOWF_MIX(J)=0.

   ELSE

      FLOWF_MIX(J)=FLOWI_MIX(J)

   ENDIF

      ENDDO

 

C    SET SOME AVERAGE TEMPERATURES FOR THE FIRST INTERNAL ITERATION

      DO J=1,NODES

         TAVE_TANK(J)=TI_TANK(J)

      ENDDO

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

               TAVE_HX(I,J,K)=TI_HX(I,J,K)

            ENDDO

         ENDDO

      ENDDO

 

C    SET THE DEFAULT VALUES FOR THE CHECKED ARRAY FOR THIS ITERATION

      DO J=1,NODES

   CHECKED(J)=.FALSE.

      ENDDO

 

C    SET SOME INITIAL CONDITIONS FOR THE ITERATIVE CALCULATIONS

      DO J=1,NODES

         ITERS_TANK(J)=1

   TANKAVE(J)=0.

      ENDDO

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

         HXAVE(I,J,K)=0.

            ENDDO

         ENDDO

ENDDO

 

      ITERS_TANK_MAX=1

DELT_USED=0.

DELT_NOW=DELT

CHECKING_TEMPS=.TRUE.

      QHX_TOT=0.

QLOSS_E=0.

QLOSS_B=0.

QLOSS_T=0.

QLOSS_F=0.

QAUX_TOT=0.

QMISC_TOT=0.

QSTORE_HX=0.

QSTORE_T=0.

QDELIV=0.

QDELIV_HX=0.

Q_HX=0.

QDEL=0.

 

C    START THE ITERATIVE CALCULATION HERE - EVERYTHING ABOVE THIS POINT IS INDEPENDENT OF THE CONVERGENCE SCHEME

200   CONTINUE

 

C    DETERMINE THE THERMAL PROPERTIES OF THE FLUIDS

      AVETANK=0.

TOT_VOL=0.

DO J=1,NODES

   AVETANK=AVETANK+TAVE_TANK(J)*TANKVOL(J)

   TOT_VOL=TOT_VOL+TANKVOL(J)

      ENDDO

AVETANK=AVETANK/TOT_VOL


IF(ITFLUID.GT.0) THEN

   CALL FLUIDPROPS(ITFLUID,AVETANK,PERVOL_T,CP_TANK,RHO_TANK,

1      MU_TANK,K_TANK)

ENDIF

 

C    DETERMINE THE CAPACITANCE OF EACH OF THE NODES FOR THE TANK AND HEAT EXCHANGER, AND THE

C    THERMAL PROPERTIES OF THE FLUIDS

      DO J=1,NODES

         CAP_T(J)=TANKVOL(J)*RHO_TANK*CP_TANK

      ENDDO

 

C    THE DIFFERENTIAL EQUATIONS WILL BE SOLVED ANALYTICALLY BY PUTTING THEM INTO THE

C    FORM dT/dt=aT+b

      DO J=1,NODES

   AA(J)=0.

   BB(J)=0.

      ENDDO

 

C    SEE IF THERE IS AN INTERNAL HEAT EXCHANGER IN THIS TANK AND SKIP THE NEXT SECTION IF NOT

      DO K=1,N_HX

 

C       EVALUATE THE HEAT EXCHANGER CONTRIBUTION AND DETERMINE NEW a AND b VALUES

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

               ITERS_HX(I,J,K)=1

            ENDDO

         ENDDO

         ITERS_HX_MAX=1

 

205      CONTINUE

 

C       DETERMINE THE THERMAL PROPERTIES OF THE HX FLUID

         AVEHX=0.

   TOT_VOL=0.

   DO I=1,N_TNKDIV(K)

      DO J=1,N_HXNODES(K)

         AVEHX=AVEHX+TAVE_HX(I,J,K)*HXVOL_I(I,J,K)

         TOT_VOL=TOT_VOL+HXVOL_I(I,J,K)

      ENDDO

         ENDDO

   AVEHX=AVEHX/TOT_VOL

   IF(IHXFLUID(K).GT.0) THEN

      CALL FLUIDPROPS(IHXFLUID(K),AVEHX,PERVOL_HX(K),CP_HX(K),

1         RHO_HX(K),MU_HX(K),K_HX(K))

   ENDIF

 

C       PERFORM THE CALCULATION FOR EACH TANK NODE AND EACH HEAT EXCHANGER NODE

         DO 209 I=1,N_TNKDIV(K)

 

C       SKIP THIS SECTION IF THE TANK NODE DOESN'T HAVE A HEAT EXCHANGER

   IF(FRACHX(I,K).LE.0.) GOTO 209

 

         DO J=1,N_HXNODES(K)

 

C          SET THE TANK NODE IN WHICH THIS HEAT EXCHANGER NODES RESIDES

            IF(IHXGEOM(K).EQ.1) THEN

         NSINK=I

      ELSE

         NSINK=N_HXFLOW(J,K)

      ENDIF

 

C          STEP 1, GUESS A SURFACE TEMPERATURE FOR THE OUTSIDE OF THE HEAT EXCHANGER

            TSURF_OLD=(TAVE_TANK(NSINK)+TAVE_HX(I,J,K))/2.

      ITERS_SURF=1

 

C          EVALUATE THE RAYLEIGH NUMBER FOR THIS NODE

208         GRAV=9.8*3600.*3600.                                   !M/H2

            T_PROPS=(TSURF_OLD+TAVE_TANK(NSINK))/2.                !C

            T_REF=T_PROPS-1.

      IF(ITFLUID.GT.0) THEN

         CALL FLUIDPROPS(ITFLUID,T_REF,PERVOL_T,CP_REF,

     1                  RHO_REF,MU_REF,K_REF)

         CALL FLUIDPROPS(ITFLUID,T_PROPS,PERVOL_T,CP_NODE,

     1                  RHO_NODE,MU_NODE,K_NODE)

         BETA=-1./RHO_NODE*(RHO_REF-RHO_NODE)/(T_REF-T_PROPS)

      ELSE

         CP_REF=CP_TANK

         K_REF=K_TANK

         MU_REF=MU_TANK

         RHO_REF=RHO_TANK

         BETA=BETA_TANK

 

         CP_NODE=CP_TANK

         K_NODE=K_TANK

         MU_NODE=MU_TANK

         RHO_NODE=RHO_TANK

            ENDIF

 

      ALPHA_NODE=K_NODE/RHO_NODE/CP_NODE                     !M2/H

      NU_NODE=MU_NODE/RHO_NODE                               !M2/H

            RA_NODE=GRAV*BETA*DABS((TSURF_OLD-TAVE_TANK(NSINK)))

     1              *(OUT_DIA(K)**3.)/ALPHA_NODE/NU_NODE           !DIMENSIONLESS

 

 

C          CALCULATE THE NUSSELT NUMBER FOR THIS NODE

            IF(RA_NODE.LE.0.) THEN

         NUSSELT=0.

            ELSE

               NUSSELT=C(K)*(RA_NODE**RA_EXP(K))*(GF(K)**GF_EXP(K))               !DIMENSIONLESS

            ENDIF

 

C          RESISTANCES ARE CALCULATED FOR ONE TUBE FOR THE LENGTH OF THE HEAT EXCHANGER NODE IN THE FLOW DIRECTION,

C          WE'LL DEFINE L_NODE AS BEING THE LENGTH OF ONE NODE OF ONE TUBE OF THE HEAT EXCHANGER

            IF(IHXGEOM(K).EQ.1) THEN

               L_NODE=L_TUBE(K)/DBLE(N_HXNODES(K))

      ELSE

         L_NODE=L_TUBE(K)*FRACHX(J,K)

      ENDIF

 

C          DEFINE THE OUTER AND INNER SURFACE AREAS FOR ONE TUBE OF ONE NODE

            OA_NODE=PI*OUT_DIA(K)*L_NODE

            IA_NODE=PI*IN_DIA(K)*L_NODE

 

C          CALCULATE THE NATURAL CONVECTION HEAT TRANSFER COEFFICIENT FOR THIS NODE

            H_OUT=NUSSELT*K_NODE/OUT_DIA(K)                              !KJ/H.M2.K

      IF(H_OUT.LE.0.) THEN

         RESIST_OUT(I,J,K)=9.9999D+20

      ELSE

               RESIST_OUT(I,J,K)=1./H_OUT/OA_NODE                        !H.K/KJ

            ENDIF

 

C          CALCULATE THE HEAT EXCHANGER PIPE WALL RESISTANCE FOR THIS NODE

            WALL_RESIST(I,J,K)=DLOG(OUT_RAD(K)/IN_RAD(K))/2./PI/

1         L_NODE/K_WALL(K)    !H.K/KJ

 

C          START THE CALCULATIONS FOR DETERMINING THE HEAT TRANSFER COEFFICIENT FOR INTERNAL

C          FLOW IN THE HEAT EXCHANGER TUBES BY CALCULATING THE DIMENSIONLESS PARAMETERS

            FLOW_TUBE=FLOWHX_IN(K)/N_TUBES(K)

      RE_HX=4.*FLOW_TUBE/PI/IN_DIA(K)/MU_HX(K)

      PR_HX=CP_HX(K)*MU_HX(K)/K_HX(K)

 

C          IF THERE IS NO FLOW, ASSUME THAT THE INNER RESISTANCE TO HEAT TRANSFER IS NEGLIGIBLE

            IF(RE_HX.LE.0.) THEN

         NUSSELT_HX=1.

         GOTO 207

      ENDIF

 

C          SET THE CRITICAL REYNOLDS NUMBER FOR LAMINAR FLOW

            IF(IHXGEOM(K).NE.3) THEN

         RE_CRIT=2300.

      ELSE

               RE_CRIT=20000.*((IN_DIA(K)/DIA_COIL(K))**0.32)

            ENDIF

 

C          DETERMINE IF THE FLOW IS LAMINAR OR TURBULENT

            IF(RE_HX.LE.RE_CRIT) THEN

         LAMINAR_HX=.TRUE.

      ELSE

         LAMINAR_HX=.FALSE.

            ENDIF

 

C          CALCULATE THE NUSSELT NUMBER BASED ON THE TYPE OF HEAT EXCHANGER AND THE

C          FLOW CHARACTERIZATION.  COILED TUBE HX'S HAVE DIFFERENT CORRELATIONS THAN

C          TUBE BANKS AND SERPENTINE TUBE HEAT EXCHANGERS - THE ASSUMPTION MADE HERE

C          IS THAT THE INTERNAL CONVECTION COEFFICIENT REMAINS CONSTANT OVER THE LENGTH

C          OF THE TUBE

            IF(LAMINAR_HX) THEN

         IF(IHXGEOM(K).NE.3) THEN

            IF((L_TUBE(K)/IN_DIA(K)).LE.(0.0425*RE_HX*PR_HX)) THEN

               NUSSELT_HX=((3.66**3.)+(1.61**3.)*RE_HX*PR_HX*

     1                           (IN_DIA(K)/L_TUBE(K)))**(1./3.)

            ELSE

               NUSSELT_HX=4.364

            ENDIF

         ELSE

            HE=(RE_HX*((IN_DIA(K)/DIA_COIL(K))**0.5))/

     1               (1.+((PITCH_COIL(K)/PI/IN_DIA(K))**2.))

                  NUSSELT_HX=((((48./11.)+(51./11./

     1                    (1.+((1342./PR_HX/HE/HE)**2.))))**3.)+

     1                    (1.816*((HE/(1.+1.15/PR_HX))**1.5))

     1                    )**(1./3.)

         ENDIF

            ELSE

         IF(IHXGEOM(K).NE.3) THEN

                  IF(PR_HX.LE.1.5) THEN

               NUSSELT_HX=0.0214*((RE_HX**0.8)-100.)*(PR_HX**0.4)

            ELSE

               NUSSELT_HX=0.012*((RE_HX**0.87)-280.)*(PR_HX**0.4)

            ENDIF

         ELSE

                  NUSSELT_HX=0.023*(RE_HX**0.85)*(PR_HX**0.4)*

     1                    ((IN_DIA(K)/DIA_COIL(K))**0.1)

         ENDIF

            ENDIF

 

C          NOW CALCULATE THE INTERNAL HEAT TRANSFER COEFFICIENT

207         H_IN=NUSSELT_HX*K_HX(K)/IN_DIA(K)

      IF(H_IN.LE.0.) THEN

         RESIST_IN(I,J,K)=9.9999D+20

      ELSE

               RESIST_IN(I,J,K)=1./H_IN/IA_NODE                           !H.K/KJ

      ENDIF

 

C          NOW CALCULATE THE OVERALL HEAT TRANSFER COEFFICIENT FOR THIS NODE

            UA_HX(I,J,K)=1./(RESIST_OUT(I,J,K)+WALL_RESIST(I,J,K)+

1         RESIST_IN(I,J,K))                          !KJ/H.K

      IF(IHXGEOM(K).EQ.1) THEN

         UA_HX(I,J,K)=UA_HX(I,J,K)*N_TUBES(K)*FRACHX(I,K)

      ELSE

         UA_HX(I,J,K)=UA_HX(I,J,K)*N_TUBES(K)

      ENDIF

 

      IF((RESIST_OUT(I,J,K).GT.9.99D+20).OR.

     1         (RESIST_IN(I,J,K).GT.9.99D+20).OR.

     1         (RESIST_OUT(I,J,K).GT.9.99D+20)) UA_HX(I,J,K)=0.

 

C          WITH THE OVERALL HEAT TRANSFER COEFFICIENT KNOWN, CALCULATE THE SURFACE

C          TEMPERATURES OF THE PIPE BASED ON THE LAST VALUES OF STORAGE AND HEAT EXCHANGER

C          TEMPERATURES FOR THIS NODE

            Q_TEMP=UA_HX(I,J,K)*(TAVE_HX(I,J,K)-TAVE_TANK(NSINK))

            IF(IHXGEOM(K).EQ.1) THEN

         TSURF_I=TAVE_HX(I,J,K)-Q_TEMP/N_TUBES(K)/FRACHX(I,K)*

1                 RESIST_IN(I,J,K)

         TSURF_O=TSURF_I-Q_TEMP/N_TUBES(K)/FRACHX(I,K)*

1            WALL_RESIST(I,J,K)

      ELSE

         TSURF_I=TAVE_HX(I,J,K)-Q_TEMP/N_TUBES(K)*RESIST_IN(I,J,K)

         TSURF_O=TSURF_I-Q_TEMP/N_TUBES(K)*WALL_RESIST(I,J,K)

      ENDIF

 

C          CHECK TO SEE IF THE SURFACE TEMPERATURE HAS CONVERGED - IF NOT, REDO THE CALCS

            IF((DABS(TSURF_O-TSURF_OLD).GT.CONVERGED).AND.

     1         (ITERS_SURF.LE.50)) THEN

         TSURF_OLD=TSURF_O

         ITERS_SURF=ITERS_SURF+1

         GOTO 208

      ENDIF

         ENDDO

209      CONTINUE

 

C       TIME TO EVALUATE THE HEAT EXCHANGER - START BY SETTING THE INLET TEMPERATURE AND

C       FLOWRATE TO EACH OF THE HEAT EXCHANGER NODES

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

 

C             THIS HX NODE CONTAINS THE ENTRY - SET TO INLET TEMP AND FLOW

         IF(J.EQ.1) THEN

            TIN_HX(I,J,K)=THXFLOW_IN(K)

C             THIS HX NODE RECEIVES FLOW FROM ANOTHER NODE

         ELSE

            TIN_HX(I,J,K)=TAVE_HX(I,J-1,K)

               ENDIF

 

         IF(IHXGEOM(K).EQ.1) THEN

            FLOW_HX(I,J,K)=FLOWHX_IN(K)*FRACHX(I,K)

         ELSE

            FLOW_HX(I,J,K)=FLOWHX_IN(K)

               ENDIF             

            ENDDO

         ENDDO

 

C       EVALUATE THE EFFECTS OF THE HEAT EXCHANGER AND INLET FLOW AND CALCULATE a AND b

C       FOR THE HEAT EXCHANGER NODES

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

               IF(IHXGEOM(K).EQ.1) THEN

            NSINK=I

         ELSE

            NSINK=N_HXFLOW(J,K)

         ENDIF

               CAP_HX(I,J,K)=HXVOL_I(I,J,K)*RHO_HX(K)*CP_HX(K)

 

         IF(CAP_HX(I,J,K).LE.0.) THEN

            AA_HX(I,J,K)=0.

            BB_HX(I,J,K)=0.

         ELSE

              AA_HX(I,J,K)=(-FLOW_HX(I,J,K)*CP_HX(K)-UA_HX(I,J,K))/

     1                 CAP_HX(I,J,K)

            BB_HX(I,J,K)=(+FLOW_HX(I,J,K)*CP_HX(K)*TIN_HX(I,J,K)

     1                 +UA_HX(I,J,K)*TAVE_TANK(NSINK))/CAP_HX(I,J,K)

         ENDIF

            ENDDO

         ENDDO

 

C       DETERMINE THE FINAL AND AVERAGE TEMPERATURES FOR THE HEAT EXCHANGER NODES

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

         IF((AA_HX(I,J,K).EQ.0.).OR.(DELT_NOW.LE.0.)) THEN

            TF_HX_NEW(I,J,K)=TI_HX(I,J,K)+BB_HX(I,J,K)*DELT_NOW

            TAVE_HX_NEW(I,J,K)=TI_HX(I,J,K)+BB_HX(I,J,K)*

1               DELT_NOW/2.

         ELSE

                  TF_HX_NEW(I,J,K)=(TI_HX(I,J,K)+BB_HX(I,J,K)/

1               AA_HX(I,J,K))*DEXP(AA_HX(I,J,K)*DELT_NOW)-

     1               BB_HX(I,J,K)/AA_HX(I,J,K)

                  TAVE_HX_NEW(I,J,K)=(TI_HX(I,J,K)+BB_HX(I,J,K)/

1               AA_HX(I,J,K))*(DEXP(AA_HX(I,J,K)*DELT_NOW)-1.)/

     1               AA_HX(I,J,K)/DELT_NOW-BB_HX(I,J,K)/AA_HX(I,J,K)

               ENDIF

            ENDDO

         ENDDO

 

C       SEE IF THE HEAT EXCHANGER NODE TEMPERATURES HAVE CONVERGED - IF SO, MOVE ON

C      AND DO THE TANK NODES, IF NOT, GO BACK TO THE HEAT EXCHANGER START AND TRY AGAIN

         HX_CONVERGE=.TRUE.

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

         IF((DABS(TAVE_HX_NEW(I,J,K)-TAVE_HX(I,J,K)).GT.CONVERGED)

1            .AND.(ITERS_HX(I,J,K).LE.50)) THEN

                  HX_CONVERGE=.FALSE.

            ITERS_HX(I,J,K)=ITERS_HX(I,J,K)+1

            ITERS_HX_MAX=MAX(ITERS_HX(I,J,K),ITERS_HX_MAX)

        ENDIF

            ENDDO

         ENDDO

      

        IF(HX_CONVERGE) THEN

      IF(ITERS_HX_MAX.GE.50) THEN

         NOTCONV_HX=.TRUE.

      ELSE

         NOTCONV_HX=.FALSE.

      ENDIF

 

      DO I=1,N_TNKDIV(K)

         DO J=1,N_HXNODES(K)

            TAVE_HX(I,J,K)=TAVE_HX_NEW(I,J,K)

            TF_HX(I,J,K)=TF_HX_NEW(I,J,K)

               ENDDO

            ENDDO

         ELSE

      DO I=1,N_TNKDIV(K)

         DO J=1,N_HXNODES(K)

            TAVE_HX(I,J,K)=TAVE_HX_NEW(I,J,K)

            TF_HX(I,J,K)=TF_HX_NEW(I,J,K)

               ENDDO

            ENDDO

            GOTO 205

   ENDIF

 

C       NOW SET THE TANK AA AND BB VALUES BASED ON THE HEAT EXCHANGER CONTRIBUTION

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

               IF(IHXGEOM(K).EQ.1) THEN

            NSINK=I

         ELSE

            NSINK=N_HXFLOW(J,K)

         ENDIF

         AA(NSINK)=AA(NSINK)-UA_HX(I,J,K)/CAP_T(NSINK)

         BB(NSINK)=BB(NSINK)+UA_HX(I,J,K)*TAVE_HX(I,J,K)/

     1            CAP_T(NSINK)

            ENDDO

         ENDDO

 

C    END OF HEAT EXCHANGER SECTION

353   ENDDO

 

C    SET THE INLET TEMPERATURES AND FLOW RATES TO EACH OF THE NODES

DO I=1,N_PORTS

         IF(INLET_MODE(I).EQ.3) THEN

            T_DIFF=1.D+20

      NODE_MIN=NODES

      DO J=1,NODES

         IF(DABS(TSTART_TANK(J)-TFLOW_IN(I)).LE.T_DIFF) THEN

            NODE_MIN=J

            T_DIFF=DABS(TAVE_TANK(J)-TFLOW_IN(I))

         ENDIF

      ENDDO

      

  DO J=1,NODES

         IF(J.EQ.NODE_MIN) THEN

           FRAC_IN(I,J)=1.0

         ELSE

           FRAC_IN(I,J)=0.0

         ENDIF

      ENDDO

         ENDIF

ENDDO

 

      DO I=1,N_PORTS

   DO J=1,NODES

      FLOW_MAINS(I,J)=FLOW_IN(I)*FRAC_IN(I,J)

   ENDDO

ENDDO

 

C    START AT THE BOTTOM AND WORK TOWARDS THE OUTLET NODE

      DO I=1,N_PORTS

         DO J=NODES,1,-1

      IF(J.EQ.N_EXIT(I)) EXIT

 

      IF(J.EQ.NODES) THEN

         T_IN_LOAD(I,J)=TAVE_TANK(J)

         FLOW_IN_LOAD(I,J)=0.

         FLOW_OUT_LOAD(I,J)=FLOW_MAINS(I,J)

      ELSE

         T_IN_LOAD(I,J)=TAVE_TANK(J+1)

         FLOW_IN_LOAD(I,J)=FLOW_OUT_LOAD(I,J+1)

         FLOW_OUT_LOAD(I,J)=FLOW_MAINS(I,J)+FLOW_IN_LOAD(I,J)

      ENDIF

   ENDDO

      

C       NOW START AT THE TOP AND WORK TOWARDS THE OUTLET NODE

         DO J=1,NODES

      IF(J.EQ.N_EXIT(I)) EXIT

 

      IF(J.EQ.1) THEN

         T_IN_LOAD(I,J)=TAVE_TANK(J)

         FLOW_IN_LOAD(I,J)=0.

         FLOW_OUT_LOAD(I,J)=FLOW_MAINS(I,J)

      ELSE

         T_IN_LOAD(I,J)=TAVE_TANK(J-1)

         FLOW_IN_LOAD(I,J)=FLOW_OUT_LOAD(I,J-1)

         FLOW_OUT_LOAD(I,J)=FLOW_MAINS(I,J)+FLOW_IN_LOAD(I,J)

      ENDIF

   ENDDO

ENDDO

 

IF(NODES.EQ.1) THEN

         DO I=1,N_PORTS

      FLOW_OUT_LOAD(I,1)=FLOW_MAINS(I,1)

   ENDDO

ENDIF

 

C    NOW SET THE CORRECT AA AND BB TERMS FOR INLET FLOW 1

      DO I=1,N_PORTS

   DO J=1,NODES

      IF(NODES.EQ.1) THEN

               AA(J)=AA(J)-FLOW_MAINS(I,J)*CP_TANK/CAP_T(J)

         BB(J)=BB(J)+FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)/CAP_T(J)

      ELSE IF(J.NE.N_EXIT(I)) THEN

               AA(J)=AA(J)-FLOW_OUT_LOAD(I,J)*CP_TANK/CAP_T(J)

         BB(J)=BB(J)+(FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)+

1            FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J))/CAP_T(J)

      ELSE

     IF(J.EQ.1) THEN

            FLOW_OUT_LOAD(I,J)=FLOW_IN(I)

            FLOW_IN_LOAD(I,J)=FLOW_OUT_LOAD(I,J+1)

            T_IN_LOAD(I,J)=TAVE_TANK(J+1)

 

                  AA(J)=AA(J)-FLOW_OUT_LOAD(I,J)*CP_TANK/CAP_T(J)

            BB(J)=BB(J)+(FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)+

1               FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J))/CAP_T(J)

      

     ELSE IF(J.EQ.NODES) THEN

            FLOW_OUT_LOAD(I,J)=FLOW_IN(I)

            FLOW_IN_LOAD(I,J)=FLOW_OUT_LOAD(I,J-1)

            T_IN_LOAD(I,J)=TAVE_TANK(J-1)

 

                  AA(J)=AA(J)-FLOW_OUT_LOAD(I,J)*CP_TANK/CAP_T(J)

            BB(J)=BB(J)+(FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)+

1               FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J))/CAP_T(J)

      

     ELSE

            FLOW_OUT_LOAD(I,J)=FLOW_IN(I)

            FLOW_IN_LOAD(I,J)=FLOW_OUT_LOAD(I,J-1)+

1               FLOW_OUT_LOAD(I,J+1)

            IF(FLOW_IN(I).GT.0.) THEN

               T_IN_LOAD(I,J)=(TAVE_TANK(J-1)*FLOW_OUT_LOAD(I,J-1)

1                  +TAVE_TANK(J+1)*FLOW_OUT_LOAD(I,J+1))/

     1                  (FLOW_OUT_LOAD(I,J-1)+FLOW_OUT_LOAD(I,J+1))

            ELSE

               T_IN_LOAD(I,J)=TAVE_TANK(J)

            ENDIF

 

                  AA(J)=AA(J)-FLOW_OUT_LOAD(I,J)*CP_TANK/CAP_T(J)

            BB(J)=BB(J)+(FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)+

1               FLOW_OUT_LOAD(I,J-1)*CP_TANK*TAVE_TANK(J-1)+

1               FLOW_OUT_LOAD(I,J+1)*CP_TANK*TAVE_TANK(J+1))

     1               /CAP_T(J)

     ENDIF

      ENDIF

         ENDDO

      ENDDO

 

C    DETERMINE THE EFFECTS OF THE AUXILIARY HEATERS

      HAVE_AUX=.FALSE.

DO J=1,NODES

         BB(J)=BB(J)+Q_AUX(J)/CAP_T(J)

   IF(Q_AUX(J).GT.0) HAVE_AUX=.TRUE.

      ENDDO

 

C    EVALUATE THE EFFECTS OF GAS FLUE LOSSES AND DETERMINE NEW a AND b VALUES - TURN OFF THE LOSSES

C    TO THE GAS FLUE IF THERE IS AUX ENERGY BEING ADDED TO THE TANK

      IF(.NOT.HAVE_AUX) THEN

   DO J=1,NODES

      AA(J)=AA(J)-UA_FLUE(J)/CAP_T(J)

      BB(J)=BB(J)+UA_FLUE(J)*T_FLUE/CAP_T(J)

         ENDDO

ENDIF

 

C    EVALUATE THE EFFECTS OF TOP LOSSES AND DETERMINE NEW a AND b VALUES

      IF(DONT_LOSS) GOTO 211

      DO J=1,NODES

   UA_TOP=U_TOP*ATOP(J)

   AA(J)=AA(J)-UA_TOP/CAP_T(J)

   BB(J)=BB(J)+UA_TOP*T_TLOSS/CAP_T(J)

      ENDDO

 

C    EVALUATE THE EFFECTS OF EDGE LOSSES AND DETERMINE NEW a AND b VALUES

      DO J=1,NODES

   UA_EDGE=U_EDGE(J)*AEDGE(J)

   AA(J)=AA(J)-UA_EDGE/CAP_T(J)

   BB(J)=BB(J)+UA_EDGE*T_ELOSS(J)/CAP_T(J)

      ENDDO

 

C    EVALUATE THE EFFECTS OF BOTTOM LOSSES AND DETERMINE NEW a AND b VALUES

      DO J=1,NODES

   UA_BOT=U_BOTTOM*ABOT(J)

   AA(J)=AA(J)-UA_BOT/CAP_T(J)

   BB(J)=BB(J)+UA_BOT*T_BLOSS/CAP_T(J)

      ENDDO

 

C    EVALUATE THE EFFECTS OF MIXING BETWEEN NODES AND DETERMINE NEW a AND b VALUES.

211   IF(DONT_MIX) GOTO 215

      DO J=1,NODES

   IF(J.EQ.1) THEN

      AA(J)=AA(J)-FLOWF_MIX(J)*CP_TANK/CAP_T(J)

      BB(J)=BB(J)+FLOWF_MIX(J)*CP_TANK*TAVE_TANK(J+1)/CAP_T(J)

   ELSE IF(J.EQ.NODES) THEN

      AA(J)=AA(J)-FLOWF_MIX(J-1)*CP_TANK/CAP_T(J)

      BB(J)=BB(J)+FLOWF_MIX(J-1)*CP_TANK*TAVE_TANK(J-1)/CAP_T(J)

   ELSE

      AA(J)=AA(J)-FLOWF_MIX(J)*CP_TANK/CAP_T(J)

     1                 -FLOWF_MIX(J-1)*CP_TANK/CAP_T(J)

      BB(J)=BB(J)+FLOWF_MIX(J)*CP_TANK*TAVE_TANK(J+1)/CAP_T(J)

     1                 +FLOWF_MIX(J-1)*CP_TANK*TAVE_TANK(J-1)/CAP_T(J)

   ENDIF

      ENDDO

 

C    EVALUATE THE EFFECTS OF CONDUCTIVE EXCHANGE BETWEEN NODES AND DETERMINE NEW a AND b VALUES.

215   IF(DONT_COND) GOTO 217

 

      DO J=1,NODES

   IF(J.EQ.1) THEN

      AA(J)=AA(J)-(K_TANK+K_ADDCOND)*ACOND(J)/L_COND(J)/CAP_T(J)

      BB(J)=BB(J)+(K_TANK+K_ADDCOND)*ACOND(J)/L_COND(J)*

1         TAVE_TANK(J+1)/CAP_T(J)

   ELSE IF(J.EQ.NODES) THEN

      AA(J)=AA(J)-(K_TANK+K_ADDCOND)*ACOND(J-1)/L_COND(J-1)/

1         CAP_T(J)

      BB(J)=BB(J)+(K_TANK+K_ADDCOND)*ACOND(J-1)/L_COND(J-1)*

     1         TAVE_TANK(J-1)/CAP_T(J)

   ELSE

      AA(J)=AA(J)-(K_TANK+K_ADDCOND)*ACOND(J)/L_COND(J)/CAP_T(J)

     1                 -(K_TANK+K_ADDCOND)*ACOND(J-1)/L_COND(J-1)/

     1                  CAP_T(J)

      BB(J)=BB(J)+(K_TANK+K_ADDCOND)*ACOND(J)/L_COND(J)*

1         TAVE_TANK(J+1)/CAP_T(J)+(K_TANK+K_ADDCOND)*ACOND(J-1)/

     1         L_COND(J-1)*TAVE_TANK(J-1)/CAP_T(J)

   ENDIF

      ENDDO

 

C    DETERMINE THE EFFECTS OF THE MISCELLANEOUS ENERGY FLOWS

217   DO J=1,N_MISC

         BB(NODE_MISC(J))=BB(NODE_MISC(J))+Q_MISC(J)/CAP_T(NODE_MISC(J))

      ENDDO

 

C    DETERMINE THE FINAL AND AVERAGE TEMPERATURES FOR THE TANK NODES

      DO J=1,NODES

   IF((AA(J).EQ.0.).OR.(DELT_NOW.LE.0.)) THEN

      TF_TANK_NEW(J)=TI_TANK(J)+BB(J)*DELT_NOW

      TAVE_TANK_NEW(J)=TI_TANK(J)+BB(J)*DELT_NOW/2.

   ELSE

            TF_TANK_NEW(J)=(TI_TANK(J)+BB(J)/AA(J))*

     1              DEXP(AA(J)*DELT_NOW)-BB(J)/AA(J)

            TAVE_TANK_NEW(J)=(TI_TANK(J)+BB(J)/AA(J))*

     1             (DEXP(AA(J)*DELT_NOW)-1.)/AA(J)/DELT_NOW-BB(J)/AA(J)

         ENDIF

      ENDDO

 

C    SEE IF THE TANK NODE TEMPERATURES HAVE CONVERGED - IF SO, MOVE ON

C    IF NOT, GO BACK TO THE TANK START AND TRY AGAIN

      TANK_CONVERGE=.TRUE.

      DO J=1,NODES

      IF((DABS(TAVE_TANK_NEW(J)-TAVE_TANK(J)).GT.CONVERGED).AND.

     1         (ITERS_TANK(J).LE.MAXITERS)) THEN

               TANK_CONVERGE=.FALSE.

         ITERS_TANK(J)=ITERS_TANK(J)+1

         ITERS_TANK_MAX=MAX(ITERS_TANK_MAX,ITERS_TANK(J))

      ENDIF

      ENDDO

      

      IF(TANK_CONVERGE) THEN

   IF(ITERS_TANK_MAX.GE.MAXITERS) THEN

      NOTCONV_TANK=.TRUE.

   ELSE

      NOTCONV_TANK=.FALSE.

   ENDIF

 

   DO J=1,NODES

      TAVE_TANK(J)=TAVE_TANK_NEW(J)

      TF_TANK(J)=TF_TANK_NEW(J)

         ENDDO

      ELSE

   DO J=1,NODES

      TAVE_TANK(J)=TAVE_TANK_NEW(J)

      TF_TANK(J)=TF_TANK_NEW(J)

         ENDDO

         GOTO 200

ENDIF

 

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

C     NO MIXING SECTION - CALCULATE TEMPS AND HEAT FLOWS AND GET OUT

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

299   IF(DONT_MIX) THEN

 

C       SINCE WE ARE NOW DONE AT THIS ITERATION, STORE THE TEMPERATURES IN THE S-ARRAY

         DO J=1,NODES

      STORED(J)=TSTART_TANK(J)

      STORED(NODES+J)=TF_TANK(J)

      STORED(2*NODES+J)=FLOWI_MIX(J)

      STORED(3*NODES+J)=FLOWF_MIX(J)

         ENDDO

 

C       STORE THE HX TEMPERATURES IN THE STORAGE ARRAY AT THE FIRST ITERATION

         NITEMS=0

         DO K=1,N_HX

            NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

   ENDDO

   NITEMS=2*NITEMS+4*NODES

 

   N_STORED=4*NODES

 

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           STORED(N_STORED+1)=TSTART_HX(I,J,K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

         DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

           STORED(N_STORED+1)=TF_HX(I,J,K)

           N_STORED=N_STORED+1

               ENDDO

            ENDDO

         ENDDO

 

C       PUT THE STORED ARRAY IN THE GLOBAL STORED ARRAY

   CALL setStorageVars(STORED,NITEMS,INFO)

   IF(ErrorFound()) RETURN 1

 

C       SET THE OUTLET TEMPS

   TEMPHX_O=0.

   AVETEMP_HX=0.

 

   DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

            IF(IHXGEOM(K).EQ.1) THEN

               AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*

1                  FRACHX(I,K)/DBLE(N_HXNODES(K))

            ELSE

               AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*

1                  FRACHX(J,K)

            ENDIF

               ENDDO

               

         IF(IHXGEOM(K).EQ.1) THEN

            TEMPHX_O(K)=TEMPHX_O(K)+TAVE_HX(I,N_HXNODES(K),K)*

1               FRACHX(I,K)

         ELSE

            TEMPHX_O(K)=TAVE_HX(I,N_HXNODES(K),K)

         ENDIF

            ENDDO

         ENDDO

 

C       SET THE AVERAGE TEMPERATURE OF THE STORAGE TANK

306      AVETEMP_T=0.

   TOTCAP_T=0.

         DO J=1,NODES

      AVETEMP_T=AVETEMP_T+TAVE_TANK(J)*CAP_T(J)

            TOTCAP_T=TOTCAP_T+CAP_T(J)

         ENDDO

         AVETEMP_T=AVETEMP_T/TOTCAP_T

 

C       CALCULATE THE HEAT TRANSFER FLOWS AND AVERAGE TOP AND BOTTOM TEMPERATURES

   DO J=1,NODES

      UA_EDGE=U_EDGE(J)*AEDGE(J)

      QLOSS_E=QLOSS_E+UA_EDGE*(TAVE_TANK(J)-T_ELOSS(J))

  UA_BOT=U_BOTTOM*ABOT(J)

      QLOSS_B=QLOSS_B+UA_BOT*(TAVE_TANK(J)-T_BLOSS)

      UA_TOP=U_TOP*ATOP(J)

      QLOSS_T=QLOSS_T+UA_TOP*(TAVE_TANK(J)-T_TLOSS)

      QAUX_TOT=QAUX_TOT+Q_AUX(J)

            QSTORE_T=QSTORE_T+CAP_T(J)*(TF_TANK(J)-TSTART_TANK(J))/DELT

      QLOSS_F=QLOSS_F+UA_FLUE(J)*(TAVE_TANK(J)-T_FLUE)

     

  DO I=1,N_PORTS

     QDEL(I)=QDEL(I)+FLOW_OUT_LOAD(I,J)*CP_TANK*TAVE_TANK(J)-

     1            FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)-

     1            FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J)

            ENDDO

         ENDDO

 

   DO I=1,N_PORTS

      QDELIV=QDELIV+QDEL(I)

   ENDDO

 

   DO J=1,N_MISC

     QMISC_TOT=QMISC_TOT+Q_MISC(J)

   ENDDO

 

         DO K=1,N_HX

      DO I=1,N_TNKDIV(K)

         DO J=1,N_HXNODES(K)

            IF(IHXGEOM(K).EQ.1) THEN

               NSINK=I

            ELSE

               NSINK=N_HXFLOW(J,K)

            ENDIF

 

            Q_HX(K)=Q_HX(K)+UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

     1               TAVE_HX(I,J,K))

            QHX_TOT=QHX_TOT+UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

     1               TAVE_HX(I,J,K))

                  QSTORE_HX(K)=QSTORE_HX(K)+CAP_HX(I,J,K)*(TF_HX(I,J,K)-

     1                      TSTART_HX(I,J,K))/DELT

                  QDELIV_HX(K)=QDELIV_HX(K)+FLOW_HX(I,J,K)*CP_HX(K)*

     1                      (TAVE_HX(I,J,K)-TIN_HX(I,J,K))

               ENDDO

            ENDDO

         ENDDO

 

C       CALCULATE THE ENERGY BALANCE ERROR FOR THE TANK AND HX

         EB_TOP_T=QSTORE_T+QLOSS_E+QLOSS_B+QHX_TOT+QLOSS_T+QDELIV-

1            QAUX_TOT+QMISC_TOT+QLOSS_F

   EB_BOT_T=DABS(QLOSS_T)+DABS(QDELIV)+DABS(QLOSS_E)+

1            DABS(QLOSS_B)+DABS(QHX_TOT)+DABS(QAUX_TOT)+

     1            DABS(QMISC_TOT)+DABS(QLOSS_F)

   EB_ERROR_T=DABS(EB_TOP_T)/DMAX1(1.0,EB_BOT_T)*100.

 

         DO K=1,N_HX

            EB_TOP_HX=QSTORE_HX(K)-Q_HX(K)+QDELIV_HX(K)

      EB_BOT_HX=DABS(Q_HX(K))+DABS(QDELIV_HX(K))

      EB_ERROR_HX(K)=DABS(EB_TOP_HX)/DMAX1(1.0,EB_BOT_HX)*100.

         ENDDO

 

         GOTO 50

      ENDIF

 

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

C     ITERATION AT THE CURRENT TIMESTEP - SEE IF ANY OF THE NODES ARE UNSTABLE.

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

      IF(CHECKING_TEMPS) THEN

         DT_CRIT=DELT_NOW

   NODE_CRIT=0

 

C       CHECK ALL BUT THE BOTTOM NODES

         DO 262 J=1,NODES-1

      ITS_DT(J)=1

      DT_OLD=0.

            IF(CHECKED(J)) GOTO 263

 

C          CASE 1 NODE WASN'T UNSTABLE BEFORE, BUT IS NOW

      IF(((TI_TANK(J)+CONVERGED).GE.TI_TANK(J+1)).AND.

     1         ((TF_TANK(J)+CONVERGED).LT.TF_TANK(J+1))) THEN

 

C             USE AN ITERATIVE NEWTON'S METHOD TO SOLVE FOR THE TIME AT WHICH THE

C             NODE BECOMES UNSTABLE

264          IF(AA(J+1).NE.0.) THEN

            FX=(TI_TANK(J+1)+BB(J+1)/AA(J+1))*

     1               DEXP(AA(J+1)*DT_OLD)-BB(J+1)/AA(J+1)

            FPX=AA(J+1)*(TI_TANK(J+1)+BB(J+1)/AA(J+1))*

     1               DEXP(AA(J+1)*DT_OLD)

               ELSE

            FX=TI_TANK(J+1)+BB(J+1)*DT_OLD

            FPX=BB(J+1)

         ENDIF

         IF(AA(J).NE.0.) THEN

            FX=FX-((TI_TANK(J)+BB(J)/AA(J))*

     1               DEXP(AA(J)*DT_OLD)-BB(J)/AA(J))

            FPX=FPX-(AA(J)*(TI_TANK(J)+BB(J)/AA(J))*

     1               DEXP(AA(J)*DT_OLD))

               ELSE

            FX=FX-(TI_TANK(J)+BB(J)*DT_OLD)

            FPX=FPX-BB(J)

         ENDIF

 

C             WITH F(X) AND F'(X) CALCULATED, DETERMINE THE CRITICAL TIMESTEP AT WHICH THE NODES

C             BECOME UNSTABLE AND SEE IF THE TIME CALCULATION HAS CONVERGED

               DT_NEW=DMIN1(DELT_NOW,(DMAX1(0.,(DT_OLD-FX/FPX))))

         DT_NEW=DMAX1(DT_NEW,DELT_MIN)

         IF((DABS(DT_NEW-DT_OLD).LE.0.0001).OR.(ITS_DT(J).GT.50))

     1            THEN

            IF(DT_NEW.LT.DT_CRIT) THEN

                     DT_CRIT=DMIN1(DELT_NOW,DT_NEW)

               NODE_CRIT=J

               ICASE=1

            ENDIF

         ELSE

            DT_OLD=DT_NEW

            ITS_DT(J)=ITS_DT(J)+1

            GOTO 264

         ENDIF

 

C          CASE 2 NODE WAS UNSTABLE BEFORE, BUT ISN'T NOW

      ELSE IF(((TI_TANK(J)+CONVERGED).LT.TI_TANK(J+1)).AND.

     1         (TF_TANK(J).GE.TF_TANK(J+1))) THEN

C             USE AN ITERATIVE NEWTON'S METHOD TO SOLVE FOR THE TIME AT WHICH THE

C             NODE BECOMES UNSTABLE

266          IF(AA(J+1).NE.0.) THEN

            FX=(TI_TANK(J+1)+BB(J+1)/AA(J+1))*

     1               DEXP(AA(J+1)*DT_OLD)-BB(J+1)/AA(J+1)

            FPX=AA(J+1)*(TI_TANK(J+1)+BB(J+1)/AA(J+1))*

     1               DEXP(AA(J+1)*DT_OLD)

               ELSE

            FX=TI_TANK(J+1)+BB(J+1)*DT_OLD

            FPX=BB(J+1)

         ENDIF

         IF(AA(J).NE.0.) THEN

            FX=FX-((TI_TANK(J)+BB(J)/AA(J))*

     1               DEXP(AA(J)*DT_OLD)-BB(J)/AA(J))

            FPX=FPX-(AA(J)*(TI_TANK(J)+BB(J)/AA(J))*

     1               DEXP(AA(J)*DT_OLD))

               ELSE

            FX=FX-(TI_TANK(J)+BB(J)*DT_OLD)

            FPX=FPX-BB(J)

         ENDIF

 

C             WITH F(X) AND F'(X) CALCULATED, DETERMINE THE CRITICAL TIMESTEP AT WHICH THE NODES

C             BECOME UNSTABLE AND SEE IF THE TIME CALCULATION HAS CONVERGED

               DT_NEW=DMIN1(DELT_NOW,(DMAX1(0.,(DT_OLD-FX/FPX))))

         DT_NEW=DMAX1(DT_NEW,DELT_MIN)

         IF((DABS(DT_NEW-DT_OLD).LE.0.0001).OR.(ITS_DT(J).GT.50))

     1            THEN

            IF(DT_NEW.LT.DT_CRIT) THEN

                     DT_CRIT=DMIN1(DELT_NOW,DT_NEW)

               NODE_CRIT=J

               ICASE=2

            ENDIF

         ELSE

            DT_OLD=DT_NEW

            ITS_DT(J)=ITS_DT(J)+1

            GOTO 266

         ENDIF

 

C          CASE 3 - ANOTHER CASE OF NODE WASN'T UNSTABLE BEFORE, BUT IS NOW.  THIS ONE IS

C          CAUSED BY AN EARLIER NODAL MIXING THAT RESETS THE INITIAL TEMPERATURE

      ELSE IF((FLOWF_MIX(J).LE.0.).AND.

     1         ((TF_TANK(J)+CONVERGED).LT.TF_TANK(J+1))) THEN

               DT_NEW=DELT_MIN

         IF(DT_NEW.LT.DT_CRIT) THEN

                  DT_CRIT=DMIN1(DELT_NOW,DT_NEW)

            NODE_CRIT=J

            ICASE=1

         ENDIF

      ELSE

               NODE_CRIT=NODE_CRIT

      ENDIF

263         CONTINUE

262      CONTINUE

 

C       CHECK TO SEE WHICH NODE IS CRITICAL

         IF(NODE_CRIT.GT.0) THEN

      IF(ICASE.EQ.1) THEN

         FLOWF_MIX(NODE_CRIT)=MIXFLOW

         CHECKED(NODE_CRIT)=.TRUE.

      ELSE

         FLOWF_MIX(NODE_CRIT)=0.

         CHECKED(NODE_CRIT)=.TRUE.

      ENDIF

 

C          UPDATE THE TIMESTEP VARIABLES

      DELT_NOW=DT_CRIT

 

            CHECKING_TEMPS=.FALSE.

      GOTO 200

 

C       NOTHING IS CRITICAL, CALCULATE THE TEMPERATURES AND ENERGY FLOWS

   ELSE

      DELT_NOW=DELT-DELT_USED

      DO J=1,NODES

         TANKAVE(J)=TANKAVE(J)+TAVE_TANK(J)*DELT_NOW/DELT

 

         UA_EDGE=U_EDGE(J)*AEDGE(J)

 UA_BOT=U_BOTTOM*ABOT(J)

         UA_TOP=U_TOP*ATOP(J)

         QLOSS_E=QLOSS_E+(UA_EDGE*(TAVE_TANK(J)-T_ELOSS(J)))*

     1                 DELT_NOW/DELT

         QLOSS_B=QLOSS_B+(UA_BOT*(TAVE_TANK(J)-T_BLOSS))*

     1                 DELT_NOW/DELT

         QLOSS_T=QLOSS_T+UA_TOP*(TAVE_TANK(J)-T_TLOSS)*

1                 DELT_NOW/DELT

         QLOSS_F=QLOSS_F+(UA_FLUE(J)*(TAVE_TANK(J)-T_FLUE))*

     1                 DELT_NOW/DELT

         QAUX_TOT=QAUX_TOT+Q_AUX(J)*DELT_NOW/DELT

               QSTORE_T=QSTORE_T+(CAP_T(J)*(TF_TANK(J)-TI_TANK(J))

     1                  /DELT_NOW)*DELT_NOW/DELT

               

         DO I=1,N_PORTS

        QDEL(I)=QDEL(I)+(FLOW_OUT_LOAD(I,J)*CP_TANK*

1               TAVE_TANK(J)-FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)-

     1               FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J))

     1               *DELT_NOW/DELT

         ENDDO

            ENDDO

 

      DO I=1,N_PORTS

         QDELIV=QDELIV+QDEL(I)

      ENDDO

 

      DO J=1,N_MISC

         QMISC_TOT=QMISC_TOT+Q_MISC(J)*DELT_NOW/DELT

      ENDDO

 

            DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

              IF(IHXGEOM(K).EQ.1) THEN

                 NSINK=I

              ELSE

                 NSINK=N_HXFLOW(J,K)

              ENDIF

                    QHX_TOT=QHX_TOT+(UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

1                 TAVE_HX(I,J,K)))*DELT_NOW/DELT

                    Q_HX(K)=Q_HX(K)+(UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

1                 TAVE_HX(I,J,K)))*DELT_NOW/DELT

              HXAVE(I,J,K)=HXAVE(I,J,K)+TAVE_HX(I,J,K)*

1                 DELT_NOW/DELT

                    QSTORE_HX(K)=QSTORE_HX(K)+(CAP_HX(I,J,K)*

1                 (TF_HX(I,J,K)-TI_HX(I,J,K))/DELT_NOW)*

     1                  DELT_NOW/DELT

              QDELIV_HX(K)=QDELIV_HX(K)+(FLOW_HX(I,J,K)*CP_HX(K)*

     1                       (TAVE_HX(I,J,K)-TIN_HX(I,J,K)))*

     1                        DELT_NOW/DELT

                  ENDDO

               ENDDO

            ENDDO

   ENDIF

 

      ELSE

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

C     JUST FINISHED A RE-RUN AT SOME CRITICAL TIMESTEP AND NEED TO CALCULATE THE TEMPERATURES

C     AND ENERGY FLOWS, THEN SET THE TIMESTEP TO THAT REMAINING, AND RUN AGAIN

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

   DO J=1,NODES

      TANKAVE(J)=TANKAVE(J)+TAVE_TANK(J)*DELT_NOW/DELT

      UA_EDGE=U_EDGE(J)*AEDGE(J)

  UA_BOT=U_BOTTOM*ABOT(J)

      UA_TOP=U_TOP*ATOP(J)

      QLOSS_E=QLOSS_E+(UA_EDGE*(TAVE_TANK(J)-T_ELOSS(J)))*

1              DELT_NOW/DELT

      QLOSS_B=QLOSS_B+(UA_BOT*(TAVE_TANK(J)-T_BLOSS))*

1              DELT_NOW/DELT

      QLOSS_T=QLOSS_T+UA_TOP*(TAVE_TANK(J)-T_TLOSS)*DELT_NOW/DELT

      QLOSS_F=QLOSS_F+(UA_FLUE(J)*(TAVE_TANK(J)-T_FLUE))*

     1              DELT_NOW/DELT

      QAUX_TOT=QAUX_TOT+Q_AUX(J)*DELT_NOW/DELT

            QSTORE_T=QSTORE_T+(CAP_T(J)*(TF_TANK(J)-TI_TANK(J))

     1               /DELT_NOW)*DELT_NOW/DELT

      DO I=1,N_PORTS

     QDEL(I)=QDEL(I)+(FLOW_OUT_LOAD(I,J)*CP_TANK*

1            TAVE_TANK(J)-FLOW_MAINS(I,J)*CP_TANK*TFLOW_IN(I)-

     1            FLOW_IN_LOAD(I,J)*CP_TANK*T_IN_LOAD(I,J))

     1            *DELT_NOW/DELT

      ENDDO

         ENDDO

 

   DO I=1,N_PORTS

      QDELIV=QDELIV+QDEL(I)

   ENDDO

 

   DO J=1,N_MISC

      QMISC_TOT=QMISC_TOT+Q_MISC(J)*DELT_NOW/DELT

   ENDDO

 

   DO K=1,N_HX

            DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

            IF(IHXGEOM(K).EQ.1) THEN

               NSINK=I

            ELSE

               NSINK=N_HXFLOW(J,K)

            ENDIF

            HXAVE(I,J,K)=HXAVE(I,J,K)+TAVE_HX(I,J,K)*DELT_NOW/DELT

                  QHX_TOT=QHX_TOT+(UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

1               TAVE_HX(I,J,K)))*DELT_NOW/DELT

                  Q_HX(K)=Q_HX(K)+(UA_HX(I,J,K)*(TAVE_TANK(NSINK)-

1               TAVE_HX(I,J,K)))*DELT_NOW/DELT

                  QSTORE_HX(K)=QSTORE_HX(K)+(CAP_HX(I,J,K)*(TF_HX(I,J,K)

1               -TI_HX(I,J,K))/DELT_NOW)*DELT_NOW/DELT

            QDELIV_HX(K)=QDELIV_HX(K)+(FLOW_HX(I,J,K)*CP_HX(K)*

1                  (TAVE_HX(I,J,K)-TIN_HX(I,J,K)))*DELT_NOW/DELT

               ENDDO

            ENDDO

         ENDDO

 

C       UPDATE THE TIMESTEP VARIABLES

   DELT_USED=DELT_USED+DELT_NOW

   DELT_NOW=DELT-DELT_USED

 

C       SET THE INTIAL TEMPERATURES TO THE FINAL TEMPERATURES

         DO J=1,NODES

      TI_TANK(J)=TF_TANK(J)

         ENDDO

         DO K=1,N_HX

      DO I=1,N_TNKDIV(K)

               DO J=1,N_HXNODES(K)

            TI_HX(I,J,K)=TF_HX(I,J,K)

               ENDDO

            ENDDO

   ENDDO

 

C       NOW RE-RUN THE CALCULATIONS AT THE REMAINING TIMESTEP

         CHECKING_TEMPS=.TRUE.

   GOTO 200

      

      ENDIF

 

 

****************************************************************************************

C    SINCE WE ARE NOW DONE AT THIS ITERATION, STORE THE TEMPERATURES IN THE S-ARRAY

      DO J=1,NODES

   STORED(J)=TSTART_TANK(J)

   STORED(NODES+J)=TF_TANK(J)

   STORED(2*NODES+J)=FLOWI_MIX(J)

   STORED(3*NODES+J)=FLOWF_MIX(J)

      ENDDO

 

C    STORE THE HX TEMPERATURES IN THE STORAGE ARRAY AT THE FIRST ITERATION

      NITEMS=0

      DO K=1,N_HX

         NITEMS=NITEMS+N_TNKDIV(K)*N_HXNODES(K)

ENDDO

NITEMS=2*NITEMS+4*NODES

 

N_STORED=4*NODES

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

        STORED(N_STORED+1)=TSTART_HX(I,J,K)

        N_STORED=N_STORED+1

            ENDDO

         ENDDO

      ENDDO

 

      DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

        STORED(N_STORED+1)=TF_HX(I,J,K)

        N_STORED=N_STORED+1

            ENDDO

         ENDDO

      ENDDO

 

C    PUT THE STORED ARRAY IN THE GLOBAL STORED ARRAY

CALL setStorageVars(STORED,NITEMS,INFO)

IF(ErrorFound()) RETURN 1

 

C    IF THERE IS NO HEAT EXCHANGER, THE OUTLET TEMP IS THE NODE TEMP CONTAINING THE EXIT

AVETEMP_HX=0.

TEMPHX_O=0.

DO K=1,N_HX

         DO I=1,N_TNKDIV(K)

            DO J=1,N_HXNODES(K)

         IF(IHXGEOM(K).EQ.1) THEN

            AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*FRACHX(I,K)

1               /DBLE(N_HXNODES(K))

         ELSE

            AVETEMP_HX(K)=AVETEMP_HX(K)+TAVE_HX(I,J,K)*FRACHX(J,K)

         ENDIF

            ENDDO

               

      IF(IHXGEOM(K).EQ.1) THEN

         TEMPHX_O(K)=TEMPHX_O(K)+TAVE_HX(I,N_HXNODES(K),K)*

1            FRACHX(I,K)

      ELSE

         TEMPHX_O(K)=TAVE_HX(I,N_HXNODES(K),K)

      ENDIF

         ENDDO

      ENDDO

 

C    SET THE AVERAGE TEMPERATURE OF THE STORAGE TANK

605   AVETEMP_T=0.

TOTCAP_T=0.

      DO J=1,NODES

   AVETEMP_T=AVETEMP_T+TANKAVE(J)*CAP_T(J)

         TOTCAP_T=TOTCAP_T+CAP_T(J)

      ENDDO

      AVETEMP_T=AVETEMP_T/TOTCAP_T

 

      EB_TOP_T=QSTORE_T+QLOSS_T+QDELIV+QLOSS_E+QLOSS_B+QHX_TOT-

1    QAUX_TOT-QMISC_TOT+QLOSS_F

EB_BOT_T=DABS(QLOSS_T)+DABS(QDELIV)+DABS(QLOSS_E)+

1    DABS(QLOSS_B)+DABS(QHX_TOT)+DABS(QAUX_TOT)+DABS(QMISC_TOT)

1   +DABS(QLOSS_F)

EB_ERROR_T=DABS(EB_TOP_T)/DMAX1(1.0,EB_BOT_T)*100.

 

      DO K=1,N_HX

         EB_TOP_HX=QSTORE_HX(K)-Q_HX(K)+QDELIV_HX(K)

   EB_BOT_HX=DABS(Q_HX(K))+DABS(QDELIV_HX(K))

   EB_ERROR_HX(K)=DABS(EB_TOP_HX)/DMAX1(1.0,EB_BOT_HX)*100.

      ENDDO

 

C    PERFORM AN ADIABATIC MIXING IF TEMPERATURE INVERSIONS HAVE OCCURRED AND THE MIXING ALGORITHM

C    IS TURNED OFF

50    IF(MIXFLOW.GE.0.) GOTO 5500

 

      NITEMS=4*NODES+2*N_HXNODES(K)*N_TNKDIV(K)         

CALL getStorageVars(STORED,NITEMS,INFO)

IF(ErrorFound()) RETURN 1

 

5150  DO JJ=1,NODES-1

          IF(TF_TANK(JJ).LT.TF_TANK(JJ+1)) GO TO 5165

      ENDDO

      GO TO 5500

5165  TMIX = 0.

      CMIX = 0.

      DO LL=JJ,NODES

         TMIX = TMIX + TF_TANK(LL)*CAP_T(LL)

         CMIX = CMIX + CAP_T(LL)

         IF(LL.EQ.NODES) GO TO 5185

         IF(TMIX/CMIX.GT.TF_TANK(LL+1)) GO TO 5185

      ENDDO

5185  TMIX = TMIX/CMIX

      DO KK=JJ,LL

        IF(AA(KK).EQ.0.) THEN

          U_T = DELT/CAP_T(KK)

          V_T = U_T/2.

  ELSE

          U_T=(DEXP(AA(KK)*DELT)-1.)/AA(KK)/CAP_T(KK)

          V_T=((DEXP(AA(KK)*DELT)-1.)/AA(KK)/DELT-1.)/AA(KK)/CAP_T(KK)

        ENDIF

         QTMIX=(TMIX-TF_TANK(KK))/U_T

         TF_TANK(KK)=TMIX

   TAVE_TANK(KK)=TAVE_TANK(KK)+QTMIX*V_T

   STORED(NODES+KK)=TF_TANK(KK)

   CALL setStorageVars(STORED,NITEMS,INFO)

   IF(ErrorFound()) RETURN 1

      ENDDO

GO TO 5150

C-----------------------------------------------------------------------------------------------------------------------

 

 

C-----------------------------------------------------------------------------------------------------------------------

C    SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT

C    SET OUTPUTS

5500  CONTINUE

      DO I=1,N_PORTS

         OUT(2*I-1)=TAVE_TANK(N_EXIT(I))

         OUT(2*I)=FLOW_IN(I)

ENDDO

 

OUT(2*N_PORTS+1)=AVETEMP_T

OUT(2*N_PORTS+2)=QDELIV

 

      DO J=1,N_PORTS

   OUT(2*N_PORTS+2+J)=QDEL(J)

ENDDO

 

OUT(3*N_PORTS+3)=QLOSS_T

OUT(3*N_PORTS+4)=QLOSS_E

OUT(3*N_PORTS+5)=QLOSS_B

OUT(3*N_PORTS+6)=QLOSS_F

OUT(3*N_PORTS+7)=QAUX_TOT

OUT(3*N_PORTS+8)=QMISC_TOT

OUT(3*N_PORTS+9)=QSTORE_T

OUT(3*N_PORTS+10)=QHX_TOT

OUT(3*N_PORTS+11)=EB_ERROR_T

 

DO J=1,NODES

   OUT(3*N_PORTS+11+J)=TAVE_TANK(J)

      ENDDO

 

DO J=1,N_HX

   OUT(3*N_PORTS+11+NODES+2*J-1)=TEMPHX_O(J)

   OUT(3*N_PORTS+11+NODES+2*J)=FLOWHX_IN(J)

ENDDO

 

DO J=1,N_HX

   OUT(3*N_PORTS+11+NODES+2*N_HX+J)=AVETEMP_HX(J)

ENDDO

 

DO J=1,N_HX

   OUT(3*N_PORTS+11+NODES+3*N_HX+1+4*(J-1))=QDELIV_HX(J)

   OUT(3*N_PORTS+11+NODES+3*N_HX+2+4*(J-1))=QSTORE_HX(J)

   OUT(3*N_PORTS+11+NODES+3*N_HX+3+4*(J-1))=Q_HX(J)

   OUT(3*N_PORTS+11+NODES+3*N_HX+4+4*(J-1))=EB_ERROR_HX(J)

ENDDO

 

RETURN 1

C-----------------------------------------------------------------------------------------------------------------------

9000  CALL MESSAGES(-1,MESSAGE2,'FATAL',IUNIT,ITYPE)

      END

C-----------------------------------------------------------------------------------------------------------------------