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

C**************************************************************************
C  LIQUID COLLECTOR-STORAGE SUBSYSTEM - THIS COMPONENT MODELS THE COMMON
C  PIECES OF A SOLAR DOMESTIC HOT WATER SYSTEM; THE COLLECTOR, PUMPS, HEAT
C  EXCHANGER AND TANK. THIS ROUTINE WAS WRITTEN TO GREATLY ENHANCE THE 
C  SIMULATION SPEED COMPARED TO USING TYPES 1 TO 5.
C
C  Modifications
C
C  JWT - 05/93 This subroutine was last seen in TRNSYS version 14.2 . I have brought it
C              back with no changes to be used in TRNSYS version 15.
C
C  GMB - 07/06 Translated to version 16 standard.
C   
C**************************************************************************
       
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: TYPE852

C-----------------------------------------------------------------------------------------------------------------------
C    ACCESS TRNSYS FUNCTIONS
	USE TrnsysFunctions
	USE TrnsysConstants

C-----------------------------------------------------------------------------------------------------------------------
C DECLARATIONS

      IMPLICIT NONE						!REQUIRES THE USER TO DEFINE ALL VARIABLES BEFORE USING THEM
	DOUBLE PRECISION XIN				!THE ARRAY FROM WHICH THE INPUTS TO THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION OUT				!THE ARRAY WHICH WILL BE USED TO STORE THE OUTPUTS FROM THIS TYPE
	DOUBLE PRECISION TIME				!THE CURRENT SIMULATION TIME - YOU MAY USE THIS VARIABLE BUT DO NOT SET IT!
	DOUBLE PRECISION PAR				!THE ARRAY FROM WHICH THE PARAMETERS FOR THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION T					!AN ARRAY CONTAINING THE RESULTS FROM THE DIFFERENTIAL EQUATION SOLVER
	DOUBLE PRECISION DTDT				!AN ARRAY CONTAINING THE DERIVATIVES TO BE PASSED TO THE DIFF.EQ. SOLVER
	DOUBLE PRECISION TIME0,TFINAL,DELT	!SIMULATION START TIME, STOP TIME, AND TIMESTEP
	INTEGER*4 ICNTRL					!AN ARRAY FOR HOLDING VALUES OF CONTROL FUNCTIONS WITH THE NEW SOLVER
      INTEGER*4 INFO(15)					!THE INFO ARRAY STORES AND PASSES VALUABLE INFORMATION TO AND FROM THIS TYPE
	INTEGER*4 NP,NI,ND                  !Number of Parameters, Inputs, and Derivatives
	INTEGER*4 NoutToCheck               !Number of Outputs for checking connections using
	INTEGER*4 NoutActual                !Number of Outputs in the OUT() array
      CHARACTER*3 OCHECK					!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE OUTPUTS
	CHARACTER*3 YCHECK					!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE INPUTS

      DOUBLE PRECISION AA,BB,T0,TBAR      !Variables used in calling Differential_Eqn() subroutine
      DOUBLE PRECISION TL     ! XIN(1)
      DOUBLE PRECISION FLOWL  ! XIN(2)
      DOUBLE PRECISION TA     ! XIN(3)
      DOUBLE PRECISION TENV   ! XIN(4)
      DOUBLE PRECISION GT     ! XIN(5)
      DOUBLE PRECISION GH     ! XIN(6)
      DOUBLE PRECISION GD     ! XIN(7)
      DOUBLE PRECISION RHO    ! XIN(8)
      DOUBLE PRECISION THETA  ! XIN(9)
      DOUBLE PRECISION SLP    ! XIN(10)

      DOUBLE PRECISION A      ! PAR(1)
      DOUBLE PRECISION XNS    ! PAR(2) = Number of collectors in series
      DOUBLE PRECISION FLOWC  ! PAR(3)
      DOUBLE PRECISION CPC    ! PAR(4)
      DOUBLE PRECISION FLOWH  ! PAR(5)
      DOUBLE PRECISION CP     ! PAR(6)
      DOUBLE PRECISION GTEST  ! PAR(7)
      DOUBLE PRECISION FRTAN  ! PAR(8)
      DOUBLE PRECISION FRUL   ! PAR(9)
      DOUBLE PRECISION B0     ! PAR(10)
      DOUBLE PRECISION TMAX   ! PAR(11)
      DOUBLE PRECISION EPS    ! PAR(12)
      DOUBLE PRECISION STOR   ! PAR(13)
      DOUBLE PRECISION RHOT   ! PAR(14)
      DOUBLE PRECISION SLOSS  ! PAR(15)
      DOUBLE PRECISION TRATIO ! PAR(16)
      DOUBLE PRECISION TI     ! PAR(17)

      DOUBLE PRECISION NS     ! Number of collectors in series
      DOUBLE PRECISION RDCONV ! Radians per degree
      DOUBLE PRECISION PI     ! Value of pi
      DOUBLE PRECISION STICK  ! Number of flip-flops of controller before sticking
      DOUBLE PRECISION CAPC
	DOUBLE PRECISION CAPH
	DOUBLE PRECISION CAPSTR
	DOUBLE PRECISION FPUL
C	DOUBLE PRECISION FRUL
	DOUBLE PRECISION RTEST
	DOUBLE PRECISION DIA
	DOUBLE PRECISION UAS
	DOUBLE PRECISION R1
	DOUBLE PRECISION R2
	DOUBLE PRECISION R3
	DOUBLE PRECISION XK
	DOUBLE PRECISION RATIO
	DOUBLE PRECISION ACAPH
	DOUBLE PRECISION CAPL
	DOUBLE PRECISION XKAT
	DOUBLE PRECISION EFFSKY
	DOUBLE PRECISION EFFGND
	DOUBLE PRECISION COSSLP
	DOUBLE PRECISION FSKY
	DOUBLE PRECISION FGND
	DOUBLE PRECISION GDSKY
	DOUBLE PRECISION GDGND
	DOUBLE PRECISION XKATB
	DOUBLE PRECISION XKATDS
	DOUBLE PRECISION XKATDG
	DOUBLE PRECISION TC
	DOUBLE PRECISION XON
	DOUBLE PRECISION QDUMP
	DOUBLE PRECISION QU
	DOUBLE PRECISION QLOSS
	DOUBLE PRECISION QL
	DOUBLE PRECISION DELTU

C Functions:
      DOUBLE PRECISION Type852_IAM

      INTEGER  IUNIT

      LOGICAL ON

      PARAMETER (NP=17,NI=10,NoutActual=20,NoutToCheck=7,ND=0)
       
      DIMENSION XIN(NI),OUT(NoutActual),PAR(NP)
	
	DIMENSION YCHECK(NI),OCHECK(NoutToCheck)
       
      DATA IUNIT/0/,RDCONV/0.017453/

      DATA PI/3.1415927/,STICK/5./

      DATA YCHECK /'TE1','MF1','TE1','TE1','IR1',
     &             'IR1','IR1','DM1','DG1','DG1'/

      DATA OCHECK /'TE1','MF1','PW1','PW1','PW1',
     &             'EN1','PW1'/

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

C-----------------------------------------------------------------------------------------------------------------------
C Post-Convergence Operations:
      IF(INFO(13).GT.0) THEN
	  RETURN 1
	ENDIF

C-----------------------------------------------------------------------------------------------------------------------
C Final Call after simulation is done:
      IF(INFO(8).EQ.-1) THEN
	  RETURN 1
	ENDIF
C-----------------------------------------------------------------------------------------------------------------------
 
C     FIRST CALL OF SIMULATION
       
      IF (INFO(7).EQ.-1) THEN

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

        CALL RCHECK(INFO,YCHECK,OCHECK)
       
        OUT(1)=TI
        OUT(14)=TI
        OUT(15)=TI
       
        INFO(6)=NoutActual
        INFO(9)=1
        INFO(10)=0

C     RETURN TO THE CALLING PROGRAM
        RETURN 1

	ENDIF

C   SET LOCAL VARIABLES TO PARAMETER VALUES. PARAMETERS WILL NOT BE RESET 
C   IF THE UNIT NUMBER HAS NOT CHANGED.
      IF (INFO(1).NE.IUNIT) THEN

        IUNIT=INFO(1)

        A=PAR(1)
        XNS=PAR(2)
        NS=INT(XNS+0.1)
        FLOWC=PAR(3)
        CPC=PAR(4)
        FLOWH=PAR(5)
        CP=PAR(6)
        GTEST=PAR(7)
        FRTAN=PAR(8)
        FRUL=PAR(9)
        B0=PAR(10)
        TMAX=PAR(11)
        EPS=PAR(12)
        STOR=PAR(13)
        RHOT=PAR(14)
        SLOSS=PAR(15)
        TRATIO=PAR(16)
        TI=PAR(17)
       
C     NO HEAT EXCHANGER PRESENT
        IF(EPS.LE.0.) THEN
          FLOWH=FLOWC
          CP=CPC
	  ENDIF

	ENDIF

C     MAIN CALCULATIONS

      CAPC=FLOWC*CPC
      CAPH=FLOWH*CP
      CAPSTR=STOR*CP*RHOT
      FPUL=-GTEST*CPC*DLOG(1.-FRUL/GTEST/CPC)
      RTEST=GTEST*CPC*(1.-EXP(-FPUL/GTEST/CPC))
      DIA=(4.*STOR/PI/TRATIO)**(1.0/3.0)
      UAS=PI*DIA*DIA*(TRATIO+0.5)*SLOSS
       
C     MODIFYING FACTORS       
      R1=XNS*CAPC/A*(1.-EXP(-FPUL*A/XNS/CAPC))/RTEST
      XK=R1*A*FRUL/CAPC/XNS
      R2=(1.-(1.-XK)**NS)/XNS/XK
      R3=1.
      
	IF(EPS.GT.0) THEN 
	  R3=1./(1.+(A*R1*R2*FRUL)*(1./EPS/DMIN1(CAPC,CAPH)-1./CAPC))
	ENDIF

      RATIO=R1*R2*R3

C     SET INPUT VARIABLES
10     TL=XIN(1)
       FLOWL=XIN(2)
       TA=XIN(3)
       TENV=XIN(4)
       GT=XIN(5)
       GH=XIN(6)
       GD=XIN(7)
       RHO=XIN(8)
       THETA=XIN(9)
       SLP=XIN(10)
       
       ACAPH=CAPH
       CAPL=FLOWL*CP
       
       IF(GT.GT.0. .AND. THETA.LT.90.) GO TO 15
C    NO RADIATION
       XKAT=0.
       GO TO 16
 
C     USE RELATIONS OF BRANDEMUEHL FOR EFFECTIVE INCIDENCE ANGLES
C     FOR DIFFUSE
15     EFFSKY=59.98-0.1388*SLP+0.001497*SLP*SLP
       EFFGND=90.-0.5788*SLP+0.002693*SLP*SLP
       
       COSSLP=COS(SLP*RDCONV)
       FSKY=(1.+COSSLP)/2.
       FGND=(1.-COSSLP)/2.
       
       GDSKY=FSKY*GD
       GDGND=RHO*FGND*GH
       
       XKATB=Type852_IAM(B0,THETA)   ! Beam incidence angle modifier
       XKATDS=Type852_IAM(B0,EFFSKY) ! Sky Diffuse incidence angle modifier
       XKATDG=Type852_IAM(B0,EFFGND) ! Ground Diffuse incidence angle modifier

       XKAT=(XKATB*(GT-GDSKY-GDGND)+XKATDS*GDSKY+XKATDG*GDGND)/GT
 
C     OUT(14) IS THE TANK TEMPERATURE FROM THE LAST CALL TO THIS
C     UNIT. OUT(15) IS THE TANK TEMPERATURE FROM THE LAST TIMESTEP.
C     OUT(16) IS THE PUMP CONTROL FUNCTION FROM THE LAST ITERATION
C     AND OUT(17) IS THE NUMBER OF TIMES THE CONTROL FUNCTION HAS
C     BEEN CHANGED IN THE CURRENT TIMESTEP.
16     IF (INFO(7) .EQ. 0) OUT(15) = OUT(14)
       IF (INFO(7) .LT. 1) OUT(17) = 0.
       T0 = OUT(15)
       TBAR = T0
       IF(INFO(7) .GT. 0.) TBAR = OUT(1)

C     SET COLLECTOR OUTLET TEMPERATURE FOR PUMP CONTROL
       TC = RATIO*A*(GT*FRTAN*XKAT-FRUL*(TBAR-TA))/CAPH+TBAR

C     SET PUMP CONTROL FUNCTION
       ON = OUT(16) .GT. 0.5
       IF (OUT(17) .GT. STICK) GO TO 17
       ON = TC .GT. TBAR
       XON = 0.
       IF (ON) XON = 1.
       IF (XON .NE. OUT(16)) OUT(17) = OUT(17) + 1.
       OUT(16) = XON

17     IF (ON) GO TO 20

C     PUMP OFF
       TC = TBAR
       ACAPH=0.
       GO TO 60

C     RELIEF VALVE
20     IF (TC .LE. TMAX) GO TO 50
       TC=TMAX
       GO TO 60

C    SET COEFFICIENTS FOR ANALYTICAL SOLUTION TO DIFFEQ
C     COLLECTOR PUMP ON
50     AA = -(RATIO*A*FRUL+CAPL+UAS)/CAPSTR
       BB = (RATIO*A*(FRTAN*XKAT*GT+FRUL*TA)+CAPL*TL+UAS*TENV)/CAPSTR
       GO TO 80

C     COLLECTOR OFF OR RELIEF VALVE OPEN
60     AA = -(ACAPH+CAPL+UAS)/CAPSTR
       BB=(ACAPH*TC+CAPL*TL+UAS*TENV)/CAPSTR

C     SOLVE THE DIFFEQ AND DETERMINE RATES
80     CALL Differential_Eqn(TIME,AA,BB,T0,T,TBAR)
       TC = RATIO*A*(GT*FRTAN*XKAT-FRUL*(TBAR-TA))/CAPH+TBAR
       QDUMP = DMAX1(ACAPH*(TC-TMAX),0.0)
       QU = ACAPH*(TC-TBAR) - QDUMP
       QLOSS = UAS*(TBAR-TENV)
       QL = CAPL*(TBAR-TL)
       DELTU = CAPSTR*(T-TI)
       
       OUT(14) = T
       OUT(1)=TBAR
       OUT(2)=FLOWL
       OUT(3)=QU
       OUT(4)=QLOSS
       OUT(5)=QL
       OUT(6)=DELTU
       OUT(7)=QDUMP
       
       RETURN 1
       END

C-----------------------------------------------------------------------------------------
C Subroutines and Functions
C-----------------------------------------------------------------------------------------


C-----------------------------------------------------------------------------------------
      FUNCTION Type852_IAM(B0,IncAngle)
C Returns IAM (Incidence Angle Modifier) as a function of B0 and IncAngle

C Arguments:
C B0 : ASHRAE incidence angle modifier constant
C IncAngle : incidence angle

	   IMPLICIT NONE

         DOUBLE PRECISION Type852_IAM
	   DOUBLE PRECISION B0,IncAngle

	   Type852_IAM=1.-B0*(1./DMAX1(0.5,COS(IncAngle*0.017453))-1.)
     &              - (1.-B0)*(DMAX1(60.,IncAngle)-60.)/30.

	END FUNCTION
C-----------------------------------------------------------------------------------------
