subroutine chs90(time2,DTDEM,QCHS,TFCHS,TACHSH,TACHS,QWCHS,LMODE,TUTCHS,ASTOR,TMOLD)
!,Tij,zdist,rdist,FRRij)
!use consts
use global
implicit none  
      real(DP) :: lmode
      integer :: clock_rate,clock_max,tstart,tend
      real(DP) :: time2,dtdem,ti3    
      real(DP) :: Tij(nzloc,Nloc),zdist(nzloc),rdist(nloc),FRRij(nzloc,Nloc) 
      real(DP) :: oldtim,tmold,old1,old2,oldaht,oldat,tutchs,tutold,qchs1,tachsh,tachs
      real(DP) :: tfchs,qwchs,qchs,dtdem2,reftem,eta,cotan,refflu,astor,ebchs
      integer :: iloc,i,it,j,iz,iflag,izloc
      SAVE
!
     LOMODE=LMODE
      OLDTIM=TIME
      TIME=Time2
      IF(TIME.NE.OLDTIM) THEN
!  ............................SAVE THE FIELD VARIABLES
        CALL TMSTO90
        TMOLD=TMSTOR
        OLD1=QWGMS
        OLD2=TIMGLO
        OLDAHT=TAHT
        OLDAT=TAT
        TUTCHS=TOUT
        TUTOLD=TOUT
        QCHS1=QDEM1
        ETOP=ETOP+DETOP
        ESIDE=ESIDE+DESIDE
        EBOT=EBOT+DEBOT
        EPLUS=EPLUS+DEPLUS
        ENEG=ENEG+DENEG
        DO 4 ILOC=1,NLOC
        OLDER(ILOC)=ESUMR(ILOC)
        DO 3 I=1,IMAXR
    3   TROLD(ILOC,I)=TR(ILOC,I)
    4   CONTINUE
        IF(IFLAGG.EQ.0)GO TO 17
        DO 5 ILOC=1,NLOC
        DO 5 IT=1,ITR1
    5   OLDE(ILOC,IT)=ESUM(ILOC,IT)
        DO 6 I=1,IMAX
        DO 6 J=1,JMAX
    6   TOLD(I,J)=T(I,J)
   17   CONTINUE

      ELSE
!............REBUILT THE FIELD VARIABLES TO THE INITIAL VALUE

        QWGMS=OLD1
        TIMGLO=OLD2
        TAHT=OLDAHT
        TAT=OLDAT
        TOUT=TUTOLD
        DO 14 ILOC=1,NLOC
        ESUMR(ILOC)=OLDER(ILOC)
        DO 13 I=1,IMAXR
   13   TR(ILOC,I)=TROLD(ILOC,I)
   14   CONTINUE
        IF(IFLAGG.EQ.0)GO TO 18
        DO 15 ILOC=1,NLOC
        DO 15 IT=1,ITR1
   15   ESUM(ILOC,IT)=OLDE(ILOC,IT)
        DO 16 I=1,IMAX
        DO 16 J=1,JMAX
   16   T(I,J)=TOLD(I,J)
   18   CONTINUE
!
      ENDIF
!
      DO 19 I=0,NRLOC
        TUTME1(I)=0.
   19 CONTINUE
      IFLAGG=0
      DEPLUS=0.
      DENEG=0.
      DEBOT=0.
      DESIDE=0.
      DETOP=0.
      TAIRH=TACHSH
      TAIR=TACHS
      TDEM=TFCHS
      QWDEM=QWCHS
      QDEM=QCHS
      DTDEM1=DTDEM
      DTDEM2=DTDEM
!-----Calculate new value of borehole thermal resistance if flow rate <> 0.
      IF (QWDEM.NE.0.) THEN
          RMPOLD=RMPIPE
!---------Flow rate per borehole
          XRB(1)=ABS(QWDEM)*NSERIE/RBORE
!---------Extrapolate Rb and Ra from data if asked for
!-------------reset the constant value to RMPIPE (case of constant RMPIPE)
          RMPIPE=RMCONS
          IF (IRAFLG.EQ.1) THEN
!-------------Calculate the effective borehole thermal resistance.
              IF (RMPIPE.LE.0.) THEN
                  WRITE(*,112)
  112             FORMAT(//,1X,'** Warning : ', &
     !        /,15X,'RMPIPE <= 0. RA EFFECTS CAN NOT BE CALCULATED!'
             /,15X,'SET Ra TO ZERO TO IGNORE THEM')
	            RETURN
                ENDIF
              IF (RAPIPE.LE.0.) THEN
                  WRITE(*,113)
  113             FORMAT(//,1X,'** Warning : ', &
     !        /,15X,'RAPIPE <= 0. RA EFFECTS CAN NOT BE CALCULATED!'
             /,15X,'THE RAPIPE VALUES MUST BE LARGER THAN ZERO')
                  RETURN
                ENDIF
              IF (IRAMET.EQ.1) THEN
!-----------------U-PIPE INSTALLATION: uniform wall temperature along the
!-----------------borehole; symmetrical case: R1=R2 => R12=4RbRa/(4Rb-Ra)
                  REFTEM=RMPIPE
                  eta= HEIGHT / (CW*XRB(1)*SQRT(RAPIPE*RMPIPE))
                  IF (eta.NE.0.) THEN
                      call coth90(eta,cotan)
                      REFTEM=RMPIPE * eta * cotan
!-------------------ELSE eta*COTH(eta)=1 and REFTEM=RMPIPE
                    ENDIF
                  RMPIPE=REFTEM
                ELSE
!-----------------COAXIAL PIPE INSTALLATION: Rbeff: minimum between uniform
!-----------------heat flux and uniform wall temperature along the borehole.
!-----------------Coaxial pipe: R12 = Ra
                  REFTEM=RMPIPE
                  eta= HEIGHT / (CW*XRB(1)*2.*RMPIPE) &
                      * SQRT(1.+4.*RMPIPE/RAPIPE)
                  IF (eta.NE.0.) THEN
                      CALL coth90(eta,cotan)
                      REFTEM=RMPIPE * eta * cotan
!-------------------ELSE eta*COTH(eta)=1 and REFTEM=RMPIPE
                    ENDIF
                  REFFLU=RMPIPE + &
                    (1./(3.*RAPIPE)) * ( HEIGHT/(CW*XRB(1)) )**2
                  IF (REFTEM.LT.REFFLU) THEN
                      RMPIPE=REFTEM
                    ELSE
                      RMPIPE=REFFLU
                    ENDIF
                ENDIF
            ENDIF
          IF (RMPOLD.NE.RMPIPE) THEN
!-------------Recalculate the first heat transfer coefficient
!-------------DP OCT 1996
              DO 201 IZ=1,NLAYST
              GFRR(2,IZ)=2.*PI/(dLOG(RMR(2)/RR(2))/RLAMR(IZ)+ &
                        RMPIPE*2.*PI)
  201         CONTINUE
!-------------END DP OCT 1996
            ENDIF
        ENDIF

   10 IFLAG=0
      IF((TIMGLO+DTDEM1).LT.DT) GO TO 20
      DTDEM1=DT-TIMGLO
      IFLAG=1
      IFLAGG=1
   20 CONTINUE
      QWGMS=QWGMS+QWDEM*DTDEM1

!********** LOCAL HEAT FLOW AROUND PIPES
      call system_clock(tstart,clock_rate,clock_max)
      CALL LOCAL90(time2)
      !,Tij,zdist,rdist,FRRij)
      call system_clock(tend,clock_rate,clock_max)

      ti3 = ti3 + real(tend - tstart)/real(clock_rate)
     
      DO 21 I=0,NRLOC
        TUTME1(I)=TUTMED(I)*DTDEM1+TUTME1(I)
   21 CONTINUE
      TAHT=TAHT + TAIRH*DTDEM1
      TAT =TAT  + TAIR *DTDEM1
      IF(IFLAG.EQ.0) GO TO 40
      DTDEM2=DTDEM2-DTDEM1
      DTDEM1=DTDEM2
      DO 543 ILOC=1,NLOC
!-------DP OCT 1996
        IZLOC=MOD(ILOC-1,NZLOC)+1
        IZ=ILAYST(IZLOC)
        DTEMP(ILOC)=ESUM(ILOC,1)/CRTOT(IZ)
!-------END DP OCT 1996
        DO 533 I=1,ITR1M1
          ESUM(ILOC,I)=ESUM(ILOC,I+1)
  533   CONTINUE
        ESUM(ILOC,ITR1)=ESUMR(ILOC)
        ESUMR(ILOC)=0.
        DO 100 I=1,IMAXR
          TR(ILOC,I)=TR(ILOC,I)-DTEMP(ILOC)
  100   CONTINUE
  543 CONTINUE

!********** MEAN FLOW RATE * HEAT CAPACITIVITY
      QWGM=QWGMS/DT
      CWQWGM=ABS(CW*QWGM)
!********** HEAT FLOWS
!-----MEAN AIR TEMPERATURES DURING THE GLOBAL TIMESTEP DT
      TAIRH= TAHT/DT
      TAIR = TAT /DT
      TAHT = 0.
      TAT  = 0.
      CALL FLOW290

!********* NEW TEMPERATURES IN INTERNAL POINTS
      CALL NTCHS90

      CALL NTLOC90

      CALL NEWT290

      CALL TSTAV90

!========== CHECK ON ADIABATIC BOUNDARY CONDICTIONS CONSISTENCE

      IF(ABS(T(IMAX-1,JREF)-T(IMAX,JREF))/DRBREF.GT.DTBREF.OR. &
         ABS(T(2,JMAX-1)-T(2,JMAX))/DZBREF.GT.DTBREF) THEN

        WRITE(*,111)
  111   FORMAT(//,1X,'** Warning : ', &
     !       /,15X,'THE GLOBAL FIELD MESH BOUNDARIES ARE TOO NEAR TO',
            /,15X,'THE STORAGE BOUNDARIES !!!', &
            /,15X,'THE GOVERNING PARAMETER FOR THE MESH GENERATION', &
            /,15X,'MUST BE INCREASED (SY = SIMULATION YEARS).')
        RETURN
      ENDIF

      IBAL=1
      QWGMS=0.
      TIMGLO=0.
      GO TO 10

   40 CONTINUE
      CALL TMSTO90

      TIMGLO=TIMGLO+DTDEM1
      ASTOR=TMSTOR
      EBCHS=DETOP+DESIDE+DEBOT
      DO 41 I=0,NRLOC
        TUTMED(I)=TUTME1(I)/DTDEM
   41 CONTINUE

!if(time2/3600 == 8760*5)then
!    write(10,110) "local"
!    110 format(2x,a,2x)
!    write(*,*) ti3
!endif

return 
end subroutine chs90