!subroutine GHX(VOLUME, DEEP, DEPTH, RBORE, R0, NSERIE, NRLOC, NZLOC, RLAMRin, CRin, &    
!               NPIPES, Rpo, Rpi, DPIPES, Lb, Lpipe, Lgap, Gap, Vref, Tref, &
!               Par21, CW, ROW, ISO, FRISO, THISO, RISLAM, TIMO3, TSTMAX, TSTARG, &   
!               TGRAD, IPRE, TCMAX, TCMIN, TCPH, TAIRM, ATAIR, TAPH, ILAY, RLAMLin, &   
!               CLin, THLin, LFPRT, LPRT, TFCHS, QWCHS, TACHSH, TACHS, CIRCUL,time,delt,out)
               
subroutine GHX(Par21, par,TFCHS, QWCHS, TACHSH, TACHS, CIRCUL,time2,out)
!,Tij,zdist,rdist,FRRij)
               
!-----------------------------------------------------------------------------------------------------------------------
!***********************************************************
!**********           LUND-DST            ******************
!**********      DUCT STORAGE MODEL       ******************
!**********     CYLINDRICAL SYMMETRY      ******************
!********** CLAESSON  EFTRING  HELLSTROM  ******************
!**********          JOHANSSON            ******************
!********** DEP. OF MATHEMATICAL PHYSICS  ******************
!********** LUND INSTITUTE OF TECHNOLOGY  ******************
!********** BOX 118,S-221 00 LUND SWEDEN  ******************
!**********           JULY 1981           ******************
!***********************************************************
!=========     TRNSYS  VERSION  BY        ==================
!=========       L. MAZZARELLA            ==================
!=========  INSTITUT FUER THERMODYNAMIC   ==================
!=========       UND WAERMETECHNIK        ==================
!=========  UNIVERSITY OF STUTTGART - FRG ==================
!=========              or                ==================
!=========  DIPARTIMENTO DI ENERGETICA    ==================
!=========    POLITECNICO DI MILANO       ==================
!=========    P.ZA L. DA VINCI  32        ==================
!=========    20133 MILANO - ITALY        ==================
!=========         OCTOBER 1989           ==================
!===========================================================     
!---------         REVISION R.4           ------------------
!---------           D. PAHUD             ------------------
!---------  LABORATORY OF ENERGY SYSTEMS  ------------------ 
!---------   SWISS FEDERAL INSTITUTE OF   ------------------
!---------     TECHNOLOGY IN LAUSANNE     ------------------ 
!---------         (LASEN, EPFL)          ------------------
!---------  CH-1015 LAUSANNE, SWITZERLAND ------------------    
!---------         OCTOBER 1996           ------------------
!-----------------------------------------------------------
!
!     LAST MODIFIED:
!       2/03 -- JWT - MODIFIED FOR TRNSYS VERSION 16
!-----------------------------------------------------------------------------------------------------------------------
! Copyright © 2004 Thermal Energy System Specialists, LLC. All rights reserved.
!use consts
use global
implicit none  
       
real(DP) :: Par21, RLAMLin,clin, &   
               THLin, TFCHS, QWCHS, TACHSH, TACHS, CIRCUL,rlamrin,crin
real(DP) :: qchs,par(15)
real(DP) :: tmold
integer :: j

real(DP) :: touts,out(10)
real(DP) :: ip1,ip2,ip3,ip4,dt1,dt2,dt3,dt4,thltot,depthz
real(DP) :: ti1, ti2, ti3, ti4, ti5, ti6
integer :: iflag
real(DP) :: Tij(nzloc,Nloc),zdist(nzloc),rdist(nloc),FRRij(nzloc,Nloc) 
save
CHARACTER(LEN=500)  :: MESSAGE1,MESSAGE2,MESSAGE3,MESSAGE4

!integer :: IDIMR,NSTLAY,IDIMLO,MAXTR1,MRLOM1,NLAY,MAXSPE,MAXRLO,MAXZLO
!
!    DATA STATEMENTS
DATA MESSAGE1/'The total number of ground layers specified exceeds &
the maximum number allowed. If you would like to run at this numb & 
er of ground layers you will have to increase the NLAY variable in &
the include file "TRNVDST.INC", recompile the source code and bui &
ld a new DLL.'/
DATA MESSAGE2/'The total number of ground layers in the storage vo &
lume exceeds the maximum number allowed. If you would like to run &
at this number of ground layers you will have to increase the NSTL &
AY variable in the include file "TRNVDST.INC", recompile the sourc &
e code and build a new DLL.'/
DATA MESSAGE3/'The minimum number of vertical subregions exceeds t &
he maximum number allowed due to the number of ground layers in th &
e storage volume. If you would like to run at this number of groun &
d layers you will have to increase the MAXZLO variable in the incl &
ude file "TRNVDST.INC", recompile the source code and build a new &
DLL.'/
DATA MESSAGE4/'The number of specified output temperatures of the &
ground exceeds the maximum number allowed. If you would like to ru &
n at this number of output temperatures you will have to increase &
the MAXSPE variable in the include file "TRNVDST.INC", recompile t &
he source code and build a new DLL.'/

integer :: i
real(DP) :: times,time2,delts,tighx
integer :: t1,t2,clock_rate,clock_max,temp
logical :: first_run
save
namelist /dst_params/ DEPTH,R0,NPIPES,Rpo,Rpi,DPIPES,Lb,Lpipe,Lgap,Gap,Tref,Par21,ISO,FRISO, &
                      THISO,RISLAM,TIMO3,TSTMAX,TGRAD,IPRE,TCMAX,TCMIN,TCPH,TAIRM,ATAIR,TAPH, &
                      ILAY,LFPRT,LPRT
                         

delts = delt*3600
times = time2*3600
nlayst = 1

!    DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME
!IF (.not.first_run) THEN
if (time2 < delt) then
    open(unit = 3, status = "old", file = "..\..\namelist_files\op_set.txt")
    read(3,nml = dst_params)
    close(3)
    ti1 = 0
    ti2 = 0
    ti3 = 0
    ti4 = 0
    ti5 = 0
    tighx = 0
    !       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
    volume = par(1)
    Deep = par(2) + depth
    rbore = par(3)
    Vref = REAL(PAR(4))/1000./3600. 
    NZLOC = jfix(PAR(9)+0.5D0)
    NSERIE = JFIX(PAR(7)+0.5D0)
    NRLOC = JFIX(PAR(8)+0.5D0)
    Tstarg = par(11)
    rlamr(1) = par(12)/3.6 !thermal conductivity
    cr(1) = par(13)*1000.   !heat capacity
    rlaml(1) = par(14)/3.6  !thermal conductivity
    cl(1) = par(15)*1000.   !heat capacity
    
    !open(unit = 519,file = "test.txt",recl = 2140000)  
    LUW_DST=50 
    IF(NRLOC > NSERIE)   NRLOC=NSERIE
    IF(NRLOC <= 0)       NRLOC=1
    IF(NZLOC <= 0)        NZLOC=1

    !       ITOPFL: FLAG TO INDICATE WHERE THE FIRST LOCATION IS SET
    !       ITOPFL<>0 THE FIRST LOCATION IS SET AT THE TOP OF THE STORE
    !       ITOPFL= 0 THE FIRST LOCATION IS SET AT THE BOTTOM OF THE STORE
    !       THE DEFAULT VALUE IS ITOPFL=0
    ITOPFL=0
    LURRB=-1
    !write(10,*) volume,deep,depth,rbore,nserie,nrloc,nzloc 
    IF(npipes >= 0.) THEN
        RMPIPE = REAL(npipes)*3.6
        RMCONS = RMPIPE
        RAPIPE = ABS(REAL(Rpo))*3.6
        NPIPES = 0
        J = 12
    ELSE IF (npipes >= -11.01) THEN
        NPIPES = -npipes
        Rb = R0                                                       !m
        Rpo = REAL(Rpo)                                            !m
        Rpi = REAL(Rpi)                                            !m
        DPIPES = REAL(DPIPES)                                          !m
        L = REAL(RLAMR(1))                                            !W/mK
        Lb = REAL(Lb)                                        !W/mK
        Lpipe = REAL(Lpipe)                                      !W/mK
        Lgap = REAL(Lgap)                                      !W/mK
        Gap = REAL(Gap)                                           !m
        Vref = REAL(Vref)                               !m3/s
        Tref = REAL(Tref)                                           !C
        Rful0 = .0002                                                 !W/m2K

        Tflg = 1    !Tflg=0 constant temperature, =1 constant flux heat trasf. cor.
        Pflg = 0    !Pflg=0 emisimmetric , =1 simmetric U pipes arrangement
        J = 21
    ELSE
        !          PAR(11)>=12: logical unit through which Rb(flow,Temp) values are provided
        LURRB = -JFIX(npipes - 0.01)
        NPIPES = 0
        J = 12
    ENDIF

    !-----  IRAFLG=0 Ra effects not calculated; IRAFLG=1 Ra effects calculated
    IRAFLG = 0
    IF (ABS(par21) > 0.0000010)    IRAFLG = 1

    !-----  IRAMET=0 const. heat flux for Ra effects; IRAMET=1 const. wall temp. (U-pipe)
    IRAMET = 0
    IF (par21 < (-0.0000010))      IRAMET =1 

    CW = REAL(par(5))
    ROW = REAL(par(6))
    CW = CW*ROW*1000.
    FRISO = REAL(FRISO)
    THISO = REAL(THISO)
    RISLAM = REAL(RISLAM)
    TIMO3 = REAL(TIMO3)
    TSTMAX = REAL(TSTMAX)
    TSTARG = REAL(TSTARG)
    TGRAD = REAL(TGRAD)
    TCMAX = REAL(TCMAX)
    TCMIN = REAL(TCMIN)
    TCPH = REAL(TCPH)
    TAIRM = REAL(TAIRM)
    ATAIR = REAL(ATAIR)
    TAPH = REAL(TAPH)

    THL(ILAY) = par(10)+1000.

    THLST(NLAYST) = DEEP-DEPTH

    !QWCHS = QWCHS/ROW/3600. 

    IF (rlamr(1) < (-1.D-6)) THEN
        !          The thermal properties inside the storage volume may vary from
        !          layer to layer. They are set as those outside the storage volume.
        !          DEPTH: depth of the top of the store.
        !          DEEP:  depth of the bottom of the store.
        IFLAG=1
        THLTOT=0.
        DEPTHZ=DEPTH
        DO 15 I=1,ILAY
           THLTOT=THLTOT+THL(I)
           IF ((THLTOT > DEPTHZ).AND.(IFLAG == 1)) THEN
              RLAMR(NLAYST)=RLAML(I)
              CR(NLAYST)   =CL(I)
              THLST(NLAYST)=THLTOT-DEPTHZ
              IF (THLST(NLAYST) >= (DEEP-DEPTHZ))  THLST(NLAYST)=DEEP-DEPTHZ
              DEPTHZ=DEPTHZ+THLST(NLAYST)

        !                IFLAG: FLAG TO STOP PROCESS WHEN DEPTHZ.GE.DEEP
              IF (DEPTHZ >= DEEP) IFLAG=0
              NLAYST=NLAYST+1
           ENDIF
        15          CONTINUE
        NLAYST=NLAYST-1

        !          REDEFINE L=RLAMR(1) AS THE ARITHMETIC MEAN OF STORAGE THERMAL CONDUCT.
        L=0.
        DO 16 I=1,NLAYST
           L=L+RLAMR(I)
        16          CONTINUE
        L=L/NLAYST

        !          NUMBER OF VERTICAL LOCAL SOLUTIONS GREATER OR EQUAL TO NLAYST
        IF (NZLOC.LT.NLAYST) NZLOC=NLAYST

    ENDIF

    IF(LPRT /= 0 .OR. LFPRT /= 0) THEN
        !LUW_DST=JFIX(PAR(J+2)+0.5D0)
        !IP1=JFIX(PAR(J+3)+0.5D0)
        !IP1=1
        DT1=DT1*3600.
        !IP2=JFIX(PAR(J+5)+0.5D0)
        DT2=DT2*3600.
        !IP3=JFIX(PAR(J+7)+0.5D0)
        DT3=DT3*3600.
        !IP4=JFIX(PAR(J+9)+0.5D0)
        !IP4=1
        DT4=DT4*3600.
        !IBLOSS=JFIX(PAR(J+11)+0.5D0)
        !NSPEC=JFIX(PAR(J+12)+0.5D0)


        !    DO 20 I=1,NSPEC
        !       IGSPEC(I)=JFIX(PAR(J)+0.5D0)
        !       JGSPEC(I)=JFIX(PAR(J+1)+0.5D0)
        !       ILSPEC(I)=JFIX(PAR(J+2)+0.5D0)
        !       J=J+3
        !20       CONTINUE
    ENDIF

    IF((LPRT /= 0 .OR. LFPRT /= 0).AND. (LUW_DST /= LUW1)) &
      OPEN(UNIT=LUW_DST,FILE='DST.DAT',STATUS='UNKNOWN')

    !       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
    !call system_clock(t1,clock_rate,clock_max)
    CALL CHSIN90(TMSTOR,CAPSTO)
    !call system_clock(t2,clock_rate,clock_max)
    !ti1 = ti1 + real(t2 - t1)/real(clock_rate) 
    
    OUT(1)=TMSTOR
    OUT(2)=0.D0
    OUT(3)=TMSTOR
       
    first_run = .true.
endif 
!else

QWCHS = QWCHS/ROW/3600.
!
IF (CIRCUL > 0.) THEN
    QWCHS = ABS(QWCHS)
!--------The fluid circulate from SUBREGION 1 to SUBREGION NLOC=NRLOC*NZLOC
!--------NZLOC=2>1: from the bottom to the top of the store     I  2  I 4 I6I
!--------NRLOC=3>1: from the center to the border of the store  I  1  I 3 I5I
ELSE
    QWCHS = -ABS(QWCHS)
!--------The fluid circulate from SUBREGION NLOC downto SUBREGION 1
!--------NZLOC=2>1: from the top to the bottom of the store     I  2  I 4 I6I
!--------NRLOC=3>1: from the border to the center of the store  I  1  I 3 I5I
ENDIF
!
IF(QWCHS == 0.) GO TO 40
IF(TMSTOR > TSTMAX) GO TO 40
    LOMODE = 3
    GO TO 45
40 CONTINUE
    QWCHS = 0.
    TFCHS = 0.
    LOMODE = 0
45 CONTINUE
QCHS=0.
!-TOUT: INSTANTANEOUS OUTLET TEMPERATURE AT THE END OF DELTS, THE TRNSYS
!-SIMULATION TIMESTEP. IF NO FLOW, TOUT IS SET TO THE INSTANTANEOUS
!-TEMPERATURE OF THE GROUND HEAT EXCHANGER AT THE END OF DELTS.
!call system_clock(t1,clock_rate,clock_max)
CALL CHS90(TIMEs,DELTs,QCHS,TFCHS,TACHSH,TACHS,QWCHS,LOMODE,TOUT,TMSTOR,TMOLD)
!,Tij,zdist,rdist,FRRij)
!call system_clock(t2,clock_rate,clock_max)

!ti2 = ti2 + real(t2 - t1)/real(clock_rate)

iF(QWCHS /= 0.) THEN
!---------TOUTS: MEAN OUTLET TEMPERATURE RETURNED BY DST
  TOUTS = TFCHS-(DEPLUS-DENEG)/(ABS(QWCHS)*CW*DELTS)
ELSE
!---------IF NO FLOW, THE MEAN OUTLET TEMPERATURE RETURNED BY DST IS SET TO THE
!---------MEAN TEMPERATURE OF THE GROUND HEAT EXCHANGER DURING THE TIMESTEP DELTS
  TOUTS = TUTMED(0)
ENDIF

OUT(1) = TOUTS
OUT(2) = ABS(QWCHS)*ROW*3600.
OUT(3) = TMSTOR
OUT(4) = (DEPLUS-DENEG)/DELTS*3.6
OUT(5) = Rmpipe  !DETOP/DELTS*3.6
OUT(6) =  tat  !DESIDE/DELTS*3.6
OUT(7) = tout   !DEBOT/DELTS*3.6 


!-EFFECTIVE VOL. HEAT CAP. OF STORAGE IS "CAPSTO*(1.-R0*R0/(R1*R1))" 
      OUT(8) = (TMSTOR-TMOLD)*CAPSTO*(1.-R0*R0/(R1*R1))*VOLUME/DELTS*3.6      
!-END DP OCT 1996
!-THE NEXT OUTPUTS ARE MEAN TEMPERATURES DURING THE SIMULATION TIMESTEP DELTS
!-OUT(9) =TUTMED(0): AVERAGE SOIL TEMPERATURE NEAR THE BOREHOLES
!-OUT(10)=TUTMED(1): AV. SOIL TEMP. NEAR THE BORE. IN RADIAL LOCATION 1 (CENTER)
!-OUT(11)=TUTMED(NRLOC):AV. SOIL TEMP. NEAR THE BORE. IN RAD. LOC. NRLOC(BORDER)
!-OUT(12)=TUTMED(NRLOC-1):AV. SOIL TEMP. NEAR THE BORE. IN RAD. LOC. (NRLOC-1)
!-   ...
!-OUT(9+NRLOC)=TUTMED(2):AV. SOIL TEMP. NEAR THE BORE. IN RAD. LOC. 2
!-IF NRLOC>11, THE LAST OUTPUTS ARE NOT RETURNED (MAX 20 OUTPUTS)
OUT(9) = TUTMED(0)               
OUT(10) =  TUTMED(1)

!write(18,*) time,taht,tat,tout

DO 47 I = 1,NRLOC-1
IF (I <= 10) OUT(10+I) = TUTMED(NRLOC+1-I)
47 CONTINUE

IF(LPRT /= 0) CALL CHSUT90(TIME2,IP1,IP2,IP3,IP4,DT1,DT2,DT3,DT4)

!if(time2 == 8760*15)then
!!    write(10,110) "chsin90","chs90"
!!    110 format(2x,a,2x,a,2x)
!    write(*,*) ti1,ti2
!endif

!endif

return
end subroutine GHX