C*********************************************************************
C
C     PROGRAM explot
C
C     運動エネルギー図作成プログラム(時刻, 条件指定)
C
C     履歴 2001/06/04 小高正嗣
C
C*********************************************************************
      program ekplot
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------

#include "grid_size_M.f"

      parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) )
      parameter ( nzp = nz+1   )
      parameter ( NT = 100 ) 
*-----------------------------------------------------
      integer ibr(nbr3)
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
      parameter ( PI = 3.141592D0 ) 

      real*8
     \   zgrid ( -2:nzb ), temp0 ( -2:nzb ), tpot0 ( -2:nzb ),
     \   dens0 ( -2:nzb ), ppai0 ( -2:nzb ), pres0 ( -2:nzb ),
     \   rvap0 ( -2:nzb ), qvap0 ( -2:nzb ),
     \   FDZ(-2:NZB)     , FDZM(-2:NZB)    

      real*8
     \     u ( -2:nxb, -2:nzb ), w ( -2:nxb, -2:nzb ),
     \     tpot ( -2:nxb, -2:nzb )

      real*4
     \    z0p ( 0:nz ), p0p ( 0:nz )   , d0p ( 0:nz ),
     \ temp0p ( 0:nz ), tpot0p ( 0:nz ), pai0p ( 0:nz ),
     \ rvap0p ( 0:nz ), qvap0p ( 0:nz )

      real*4
     \     ETHERM(NZ,NT), EKINET(NZ,NT), ETOTAL(NZ,NT), zgrid_p(nz),
     \     RMIN, RMAX

      real*4
     \     DTHERM(NZ), DKINET(NZ), DTOTAL(NZ)

      character*2 mint
      character*2 lt
      character*2 date
      character*17 header

      character*100 dumy
      character CSGI*1
*-----------------------------------------------------

      RMAX = 5.0E-4
      RMIN = - RMAX

      open (11,file='exparam')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read (11,*) ylat, als
      close(11)

      tinterv = nloop2 * ( dtime / 2. )

      call FOPEN ( irun+1 )

      do iz = 1,9
         read(16,'(a100)') dumy
      end do

      do iz = -2, nzb
         read(16,*) ii, zgrid(iz), pres0(iz), dens0(iz), temp0(iz),
     \                     tpot0(iz), ppai0(iz), qvap0(iz), rvap0(iz)
      end do

      do iz = 0, nz
         z0p(iz)   = REAL(zgrid(iz)) / 1000.0
         p0p(iz)   = REAL(pres0(iz))
         d0p(iz)   = REAL(dens0(iz))
         temp0p(iz)= REAL(temp0(iz))
         tpot0p(iz)= REAL(tpot0(iz))
         pai0p(iz) = REAL(ppai0(iz))
         rvap0p(iz)= REAL(rvap0(iz))
         qvap0p(iz)= REAL(qvap0(iz))
      end do

      do iz = 1, nz
         zgrid_p(iz) = zgrid(iz-1) / 1000.
      end do

      call setcst
     O   ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT )

      call setgrd
     O   ( ZGRID, FDZ, FDZM, DX,
     D     NZ, NZB                )

      write(6,*) 'data number = ', nloop1
      write(6,*) 'INPUT DATA BEGINING, END, & INTERVAL' 
      write(6,*) 'itb ite interval ?'
      read (*,*) itb, ite, interv

      call sgpwsn
      read(*,*) iws
      call gropn(iws)
C      call sldiv ('Y', 3, 2 )
C      call sldiv ('Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )

      do it = 1, ite

      read(31) time, ((u(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(33) time, ((w(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)

      write(0,*) 'reading time = ', time

c*** 時刻変換: 秒から LT へ************************************
      time = time + times
      write(date,'(I2.2)') INT( time / day ) + 1

      time = ( time - INT( time / day ) * day ) / hour 
      
      if ( INT(time) .ge. 24 ) then
         time = time - 24
      end if

      write(lt, '(I2.2)') INT(time)
      write(mint, '(I2.2)') INT( MOD(time,1)*60.0 )

      header = 'day=' // date // ' ' // 'LT=' // lt // ':' // mint
c**************************************************************

      do IZ = 1, NZ
         ETHERM(iz,it) = 0.0
         do IX = 0, NX-1
            ETHERM(iz,it) = ETHERM(iz,it) +
     \           CPDRY * DENS0(iz) * 
     \           tpot(ix, iz) * PPAI0(IZ)
         end do
         ETHERM(iz,it) = ETHERM(iz,it)/ REAL(NX)
      end do
      

      DO IZ = 1, NZ
         EKINET(iz,it) = 0.0
         DO IX = 1, NX
            EKINET(iz,it) = EKINET(iz,it) +
     \           DENS0(iz) * ( U(IX,IZ) ** 2 ) / 2.0D0
         end do
      end do

      DO IZ = 1, NZ-1
         DO IX = 1, NX
            EKINET(iz,it) = EKINET(IZ,it)
     \           +   ( DENS0(IZ-1) + DENS0(IZ) ) / 2.0D0
     \           * W(IX,IZ) ** 2 / 2.0D0
         end do
         EKINET(iz,it) = EKINET(iz,it) / REAL(NX)
      end do

      do ix = 1, NZ
         ETOTAL(iz,it) = EKINET(iz,it) + ETHERM(iz,it)
      end do

c**************************************************************

      if ( it .ge. itb   .and.  it .le. ite ) then
         if ( mod( (it-itb), interv ) .eq. 0 ) then

            write(0,*) 'writing time = ', time
            
            do iz = 1, NZ
               DKINET(iz) = 
     \              ( EKINET(iz,it) - EKINET(iz,it-1) ) / tinterv

               DTHERM(iz) = 
     \              ( ETHERM(iz,it) - ETHERM(iz,it-1) ) / tinterv

               DTOTAL(iz) = DKINET(iz) + DTHERM(iz)
            end do

            call grfrm
            call grswnd( RMIN, RMAX, 0., 20. )
            call grsvpt( 0.2, 0.8, 0.2, 0.8 )
            call grstrn(1)
            call grstrf
            call ussttl ( 'energy', 'J/mm', '', 'Z[km]' )
            call usdaxs
            call uxsttl ( 'T', header, 1. )
            call sglset( 'LCHAR', .TRUE. )

            call sgsplc( 'kinet' )
            call SGSPLI( 24 )
            call sgplu(nz-1, DKINET, zgrid_p)
c            call sgsplc( 'therm' )
c            call SGSPLI( 34 )
c            call sgplu(nz, DTHERM, zgrid_p)
c            call sgsplc( 'total' )
c            call SGSPLI( 44 )
c            call sgplu(nz, ETOTAL, zgrid_p)
            call sglset( 'LCHAR', .FALSE. )
*
         end if
      end if

      end do

      call grcls

      call FCLOSE

      end
