SUBROUTINE CHSUT90(time2,IP1,IP2,IP3,IP4,DT1,DT2,DT3,DT4)
!use consts
use global
implicit none      
integer :: iday,ih,i,j
real(DP) :: time2
real(DP) ::  T1 ,T2,T3,T4,ctime,iy,rs,ip1,dt1,ig,jg,iloc,il,tspec
real(DP) :: rmrad,ip2,dt2,ip3,dt3,ip4,dt4
save

T1 = 0.
T2=0.
T3=0.
T4=0.
CTIME = time2
IY = int(CTIME/31536000.)
IDAY = int(CTIME/86400.)
IDAY = MOD(IDAY,365)
IH = int(CTIME/3600.)
IH = MOD(IH,24)
RS = MOD(CTIME,3600.)
69 FORMAT(' *DST*',20(1H*),1X,I3,'Y ',I3,'D ',I2,'H ',F5.0, &
 'S ',28(1H*))

!********** PRINT OF RESULTS

IF(IP1 == 0) GO TO 100
IF(time2-T1 < DT1) GOTO 100
T1 = time2
WRITE(LUW_DST,69) IY,IDAY,IH,RS
WRITE(LUW_DST,61) TAIRH,TAIR,TMSTOR,TIN,TOUT,TDEM,EPLUS,ENEG, &
 QDEM1,QDEM,LOMODE,QWDEM
61    FORMAT(' *DST1* TAIRH=',F6.2,' TAIR=',F6.2,' TMSTOR=',F6.2, &
 ' TIN=',F6.2,' TOUT=',F6.2,' TDEM=',F6.2,/, &
 ' *DST1* E+=',G10.5,'    E-=',G10.5,'   Q=',G11.5,'  QDEM=', &
 G11.5,/,' *DST1* LOMODE=',I2,8X,'QWDEM=',G11.5)
IF(NSPEC == 0) GO TO 68
DO 65 I = 1,NSPEC
IG = IGSPEC(I)
JG = JGSPEC(I)
ILOC = INVOL(IG,JG)
IL = ILSPEC(I)
IF(IL > IMAXR) GO TO 55
IF(ILOC == 0) THEN
TSPEC = T(IG,JG)
ELSE
TSPEC = T(IG,JG)+TR(ILOC,IL)+TSF(IG,JG)*TSFH(IL)
ENDIF
RMRAD = RMR(IL)
GO TO 56
55 TSPEC = T(IG,JG)
RMRAD = 1.E+10
56 CONTINUE
WRITE(LUW_DST,62) I,RM(IG),ZM(JG),RMRAD,TSPEC
62 FORMAT(' *DST1* POINT ',I2,' RADIUS=',F6.2,' DEPTH=', &
 F6.2,' LOCAL RADIUS=',F6.3,' TEMP=',F7.2)
65 CONTINUE
68 CONTINUE
IF(IBLOSS == 0) GO TO 79
WRITE(LUW_DST,71) ETOP,ESIDE,EBOT
71 FORMAT(' *DST1* ACCUMULATED HEAT LOSS THROUGH : TOP=', &
 G11.5,' SIDE= ',G11.5,/,' *DST1* BOTTOM=',G11.5)
79 CONTINUE
WRITE(LUW_DST,720)

100 CONTINUE

IF(IP2 == 0) GO TO 200
IF(time2-T2 < DT2) GOTO 200
T2=time2
WRITE(LUW_DST,69) IY,IDAY,IH,RS
WRITE(LUW_DST,718)
718   FORMAT(/,' *DST2* GLOBAL TEMPERATURE FIELD')
CALL PRDST390(T,IDIM,JDIM,IMAX1,JMAX1,RM,ZM,LUW_DST)

200 IF(IP3 == 0) GO TO 300
IF(time2-T3 < DT3) GOTO 300
T3 = time2
WRITE(LUW_DST,69) IY,IDAY,IH,RS
WRITE(LUW_DST,128)
128 FORMAT(' *DST3* GLOBAL HEAT FLOW FIELD')
WRITE(LUW_DST,130)
130   FORMAT(/' *DST3* HEAT FLOW IN POSITIVE R-DIRECTION  (W/M2)')
DO 131 J = 2,JMAX
DO 131 I = 2,IMAX
131   RLAM(I,J) = FR(I,J)/AREAR(I,J)
CALL PRDST390(RLAM,IDIM,JDIM,IMAX,JMAX,R,ZM,LUW_DST)

WRITE(LUW_DST,132)
132   FORMAT(/' *DST3* HEAT FLOW IN POSITIVE Z-DIRECTION  (W/M2)')
DO 133 J = 2,JMAX
DO 133 I = 2,IMAX
133   RLAM(I,J) = FZ(I,J)/AREAZ(I,J)
CALL PRDST390(RLAM,IDIM,JDIM,IMAX,JMAX,RM,Z,LUW_DST)

719   WRITE(LUW_DST,720)
720   FORMAT(1X,75(1H*))

300 CONTINUE
IF(IP4 == 0) GO TO 400
IF(time2-T4 < DT4) GOTO 400
T4 = time2
WRITE(LUW_DST,69) IY,IDAY,IH,RS
WRITE(LUW_DST,201)
201 FORMAT(/' *DST4* LOCAL TEMPERATURE FIELD',/, &
'        THE TEMPERATURE AROUND A PIPE IS GIVEN BY :',/, &
 '        T(GLOBAL)+TR(LOCAL)+TSF(AMPLITUDE)*TSFH')

CALL PRLOC90

WRITE(LUW_DST,202)
202 FORMAT(' *DST4* STEADY FLUX SOLUTION - FUNCTION H')

CALL PRDST190(TSFH,RMR,IMAXR,LUW_DST)

WRITE(LUW_DST,203)
203 FORMAT(/' *DST4* STEADY FLUX SOLUTION - AMPLITUDES')

CALL PRDST390(TSF,IDIM,JDIM,IVOHSX,JVOHSX,RM,ZM,LUW_DST)

WRITE(LUW_DST,720)

400 CONTINUE

RETURN
END subroutine chsut90