SUBROUTINE INCHS90
!use consts
use global
implicit none  
    
REAL(DP) :: CRIN(NSTLAY),kloc,irmin,irmax,jjj,jzmin,jzmax,ir,rbh,rbhs,ra,volij,a1,cstor 
integer :: i,j,jz,iz,izloc,izzz,iloc,ist,jst
save

   96 FORMAT(10I10)
   98 FORMAT(10G10.4)
      DO 14 I=1,IDIM
      DO 14 J=1,JDIM
   14 GVOHS(I,J)=0.
!-define the flow path: if ITOPFL=0, the first cell is set at the bottom of
!-the store. If ITOPFL<>0, it is set at the top. The default value is 0 -DP
      KLOC=0
      DO 180 I=1,NRLOC
       IRMIN=IRLOC(I)
       IRMAX=IRLOC(I+1)-1
       DO 181 J=1,NZLOC
        IF (ITOPFL.EQ.0) THEN
           JJJ=NZLOC+1-J
!----------the first cell is set at the bottom of the store
          ELSE
           JJJ=J
!----------the first cell is set at the top of the store
          ENDIF
        JZMIN=JZLOC(JJJ)
        JZMAX=JZLOC(JJJ+1)-1
        KLOC=KLOC+1
         DO 182 IR=IRMIN,IRMAX
          DO 183 JZ=JZMIN,JZMAX
           GVOHS(IR,JZ)=KLOC
  183     CONTINUE
  182    CONTINUE
  181   CONTINUE
  180 CONTINUE
!-Initialize the thermal properties of the ground within the store volume.-DP
!-DP OCT 1996
      DO 170 I=2,IRR
        DO 171 IZ=1,NLAYST
          DO 172 J=JTZS(IZ),JBZS(IZ)
            RLAM(I,J)=RLAMR(IZ)
            CIN(I,J)=CR(IZ)
  172     CONTINUE
  171   CONTINUE
  170 CONTINUE
!-Fill the ILAYST(MAXZLO) array; it will give the IZ ground layer in the store
!-knowing the IZLOC vertical location.
      IZLOC=0
      DO 400 IZ=1,NLAYST
        IF (ITOPFL.EQ.0) THEN
          IZZZ=NLAYST+1-IZ
        ELSE
          IZZZ=IZ
        END IF
        DO 401 JZ=1,NZLOCS(IZZZ)
          IZLOC=IZLOC+1
          ILAYST(IZLOC)=IZZZ
  401   CONTINUE
  400 CONTINUE
!-END DP OCT 1996
!-end DP95-------------------------------------------------------------------
      NLOC=0
      DO 15 I=1,IMAX1
      DO 15 J=1,JMAX1
      INVOL(I,J)=dint(GVOHS(I,J))
   15 NLOC=MAX(INVOL(I,J),NLOC)
      NLOC1=NLOC+1
      IMAXR1=IMAXR+1

!-Check if the maximum number of local solutions is exceeded---------------DP
      IF(NLOC.GE.IDIMLO) THEN
!-------DP OCT 1996
        WRITE(*,23) NLOC,IDIMLO,NRLOC,NZLOC
        RETURN
      END IF
   23 FORMAT(' ***ERROR***IN TYPE 60',/, &
     ' THE NUMBER OF LOCAL SOLUTIONS (NLOC=',I3,')',/, & 
     ' EXCEEDS OR IS EQUAL TO IDIMLO=',I3,/, &
     ' NLOC IS EQUAL TO NRLOC x NZLOC (',I3,' x ',I3,').',/, &
     ' IF YOU STILL WANT TO RUN YOUR CASE, INCREASE THE VALUE OF',/, &
     ' "IDIMLO" IN THE "TRNVDST.INC" INCLUDE FILE AND RECOMPILE.')
!-----END DP OCT 1996

!.............................. Calculate Borehole Termal Resistance   6/93/LM
      IF(NPIPES.GT.0) THEN
        CALL BHRES90(HEIGHT,Vref,DPIPES,NPIPES,Rb,Rpo,Rpi,L,Lb,Lpipe,Lgap, &
                      Gap,Rful0,Tref,Tflg,Pflg,Rbh,Rbhs,Ra,Reynol )
 
        RMPIPE=Rbh
        RMCONS=RMPIPE
        RAPIPE=Ra
      ENDIF
!.....................................................................end 6/93/LM
      RMPOLD=RMPIPE
      T0=10.
      R0=RR(2)
      R1=RR(IMAXR1)
!-----DP OCT 1996
      DO 410 IZ=1,NLAYST
        D1(IZ)=R1*R1*(dLOG(R1/R0)-0.75+2.*PI*RLAMR(IZ)*RMPIPE)/2.
        D1(IZ)=SQRT(D1(IZ))
  410 CONTINUE
!-----END DP OCT 1996
      IVOHSN=IMAX
      IVOHSX=2
      JVOHSN=JMAX
      JVOHSX=2
      GVOHST=0.
      VOLST=0.
      DO 103 ILOC=1,NLOC
  103 VOLOC(ILOC)=0.
      DO 100 I=2,IMAX
      DO 100 J=2,JMAX
      IF(GVOHS(I,J).LT..000001) GO TO 100
      IST=I
      JST=J
      IVOHSN=MIN0(IVOHSN,I)
      IVOHSX=MAX0(IVOHSX,I)
      JVOHSN=MIN0(JVOHSN,J)
      JVOHSX=MAX0(JVOHSX,J)
      VOLIJ=PI*(R(I+1)**2-R(I)**2)*(Z(J+1)-Z(J))
!-----DP OCT 1996
      VOL(I,J)=VOLIJ
!-----END DP OCT 1996
      ILOC=INVOL(I,J)
      VOLOC(ILOC)=VOLOC(ILOC)+VOLIJ
      VOLST=VOLST+VOLIJ
!-----DP OCT 1996
      IZLOC=MOD(ILOC-1,NZLOC)+1
      IZ=ILAYST(IZLOC)
      GVOHS(I,J)=RLAM(I,J)*VOLIJ/(D1(IZ)*D1(IZ))
!-----END DP OCT 1996
      GVOHST=GVOHST+GVOHS(I,J)
  100 CONTINUE

!-----DP OCT 1996
      DO 420 IZ=1,NLAYST
        TSF1(IZ)=R1*R1/(2.*D1(IZ)*D1(IZ))
 420  CONTINUE
!-----END DP OCT 1996
      PLST=VOLST/(PI*R1*R1)
!-----CPST=CW*PLST*PI*R0*R0
      DO 115 ILOC=1,NLOC
      DO 112 I=1,IMAXR
      TR(ILOC,I)=0.
  112 CONTINUE
  115 CONTINUE
      DO 130 I=1,IMAX
      DO 122 J=1,JMAX
      TSF(I,J)=0.
  122 CONTINUE
  130 CONTINUE
      TSFH(1)=0.
      DO 140 I=2,IMAXR
      A1=RMR(I)/R1
      TSFH(I)=0.5*A1*A1-dLOG(A1)-0.75
  140 CONTINUE
!-----CST=CR*PI*R1*R1
      DO 260 ILOC=1,NLOC
      VOLOC(ILOC)=VOLOC(ILOC)/VOLST
  260 PLSTL(ILOC)=VOLOC(ILOC)*PLST
!-----DP OCT 1996
      DO 430 IZ=1,NLAYST
        CRIN(IZ)=CR(IZ)*(1.-R0*R0/(R1*R1))
  430 CONTINUE
      DO 270 I=IVOHSN,IVOHSX
      DO 270 J=JVOHSN,JVOHSX
      IF(INVOL(I,J).EQ.0) GO TO 270
      ILOC=INVOL(I,J)
      IZLOC=MOD(ILOC-1,NZLOC)+1
      IZ=ILAYST(IZLOC)
      CIN(I,J)=CRIN(IZ)
  270 CONTINUE
!-----DP OCT 1996
!-----Calculate CAPSTO, the mean vol. heat cap. of the ground in stor. region
      CAPSTO=0.
      DO 440 IZ=1,NLAYST
        CAPSTO=CAPSTO+CR(IZ)*THLST(IZ)
  440 CONTINUE
      CAPSTO=CAPSTO/HEIGHT
      CALL TSTAV90

      CSTOR=CAPSTO*VOLST
!-----END DP OCT 1996
      IF(LFPRT.NE.0) WRITE(LUW_DST,300)VOLST,CSTOR,R0,DISPIP,PLST
  300 FORMAT(' D U C T   S T O R A G E',/,' TOTAL VOLUME  ',1PE12.5, &
     '   TOTAL HEAT CAPACITY  ',E12.5,/,' PIPES : RADIUS ',E12.5, &
     '   SPACING ',E12.5,'   TOTAL LENGTH ',E12.5)     
!write(516,*) rbh,rmpipe,ra,capsto,cstor

RETURN
END subroutine inchs90