!=======================================================================
!  エネルギー収支解析プログラム
!
!  2001/05/31  小高正嗣
!
!=======================================================================
      program energy
!-----------------------------------------------------------------------
      implicit real*8 (a-h, o-z)
!-----------------------------------------------------------------------
      parameter (ndim = 300)

      real*4 time(ndim), ETHERM(ndim), EKINET(ndim), ETOTAL(ndim),
     \     CTHKIN(ndim), DEKVIS(ndim), DEKNLV(ndim), EKPROG(ndim)

      real*4 hour, TMAX, ETMAX, ETMIN, EKMAX, EGMAX, EGMIN
      real*4 xvpmin, xvpmax, yvpmin, yvpmax

!-----------------------------------------------------------------------
      xvpmin = 0.2
      xvpmax = 0.8
      yvpmin = 0.2
      yvpmax = 0.6

      hour =3600.0

      open(11)
      read(11,*) idn
      close(11)

      open(10)

      do id = 1, idn
         read(10,*) time(id)
         read(10,*) ETHERM(id)
         read(10,*) EKINET(id)
         read(10,*) ETOTAL(id)
         read(10,*) CTHKIN(id)
         read(10,*) DEKVIS(id)
         read(10,*) DEKNLV(id)
      end do

      EKPROG(1) = CTHKIN(1) + DEKVIS(1) +DEKNLV(1)
      do id = 2, idn
         EKPROG(id) = EKPROG(id-1) + CTHKIN(id) + DEKVIS(id) +DEKNLV(id)
      end do

      close(10)

      do id = 1, idn
         time(id) = time(id) / hour
      end do

      ETMAX = 0. ; ETMIN = 0. ; EGMAX = 0. ; EGMIN = 0. ; EKMAX = 0. 

      do id = 1, idn
         ETMAX = max( ETMAX, ETOTAL(id) )
         ETMIN = min( ETMIN, ETOTAL(id) )
         EGMAX = max( EGMAX, CTHKIN(id), DEKVIS(id), DEKNLV(id) )
         EGMIN = min( EGMIN, CTHKIN(id), DEKVIS(id), DEKNLV(id) )
         EKMAX = max( EKMAX, EKINET(id) )
      end do

      ETMAX = 6.0E+10
      EKMAX = 6.0E+8

      call sgpwsn
      read(*,*) iws

      CALL SWISET('IHEIGHT', 280)
      CALL SWISET('IWIDTH', 770)
      call swlset( 'LDUMP', .TRUE. )

      call gropn ( iws )
C      call sldiv ( 'Y', 1, 2 )
      call sldiv ( 'Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )
      call sglset( 'LFULL', .TRUE. )
      call uzfact(1.2)

      call grfrm
      call grswnd( time(1), time(idn), ETMIN, ETMAX )
      call grsvpt( xvpmin, xvpmax, yvpmin, yvpmax )
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Time', 'hour', 'Energy', 'J' ) 
      call usdaxs

      call sgspli( 43 )
      call sgplu ( idn, time, ETOTAL )

C      call sgspli( 23 )
C      call sgplu ( idn, time, EKINET )

      call grfrm
      call grswnd( time(1), time(idn), 0.0, EKMAX )
      call grsvpt( xvpmin, xvpmax, yvpmin, yvpmax )
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Time', 'hour', 'Kinetic Energy', 'J' ) 
      call usdaxs

      call sgspli( 23 )
      call sgsplc( 'kinetic' )
      call sgplu ( idn, time, EKINET )

      call grcls

      end