您所在位置：主题：各位大神帮我解答下FORTAN
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-----------------------------------------------------------------------------------------------------------------------