C*********************************************************************
C
C     PROGRAM turbprf
C
C     水平平均鉛直分布図作成プログラム(拡散係数)
C
C     履歴 1997/07/24  中島健介
C          1999/06/15  小高正嗣
C          2000/11/30  小高正嗣
C
C*********************************************************************
      program turbprf
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------
#include "grid_size_M.f"
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )

      real*8 time, kturbread(-2:nxb,-2:nzb)

      real*4 kturb(nx+1,nz+1), kturbz(nz+1), zgrid(nz+1)

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

      real*4 Kmix, Kmax, zmax
      parameter ( Kmin =0.0, Kmax =40.0 )

*-----------------------------------------------------

      real*8 z0(-2:nzb), p0(-2:nzb), d0(-2:nzb), 
     \       temp0(-2:nzb), tpot0(-2:nzb), ppai0(-2:nzb)

*-----------------------------------------------------

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

      nloop1 = nloop1 + 1

      call FOPEN ( irun+1 )

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

      do iz = -2, nzb
         read(16,*) ii, z0(iz), p0(iz), d0(iz), temp0(iz),
     \                     tpot0(iz), ppai0(iz), dum1, dum2
      end do

      write(6,*) 'data number = ', nloop1
      write(0,*) 'ibegin(2-', nloop1, 
     &           ') iend(2-', nloop1,  
     &           ') inter(1-', nloop1, ')'
      read (5,*) itb, ite, interv

      write(6,*) 'z-hight number ?(I;1-',nz,')'
      read(5,*) izz

      do iz = 1, izz
         zgrid(iz) = z0(iz-1) / 1000.0
      end do

      zmax = zgrid(izz)

      call sgpwsn
      read(*,*) iws

      CALL SWISET('IWIDTH', 420)
      CALL SWISET('IHEIGHT', 420)
      call sglset( 'LFULL', .TRUE. )
      call swlset( 'LDUMP', .TRUE. )

      call gropn(iws)
      call sglset( 'LCORNER', .FALSE. )
      CALL sglset( 'LCHAR', .TRUE. )

      call grfrm
      call grswnd( Kmin, Kmax, 0., zmax )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn(1)
      call grstrf
      call ussttl ( 'Kturb(z)', 'm^2/s', '', 'Z[km]' )
      call usdaxs

      li = 23

      do it = 1, ite

      read(38) time, ((kturbread(ix,iz),ix=-2,nx+2), iz=-2,nz+2)

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

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

            write(0,*) 'writing 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 ix = 1, nx+1
               do iz = 1,nz+1
                  kturb(ix,iz) = kturbread(ix-1,iz-1)
               end do
            end do

            do iz = 1,nz+1
               kturbz(iz) = 0.
               do ix = 1, nx+1
                  kturbz(iz) = kturbz(iz) + kturb(ix,iz)
               end do
               kturbz(iz) = kturbz(iz)/(nx+1)
            end do

            call sglset( 'LCHAR', .TRUE. )
            call sgsplc( lt ) 
            call sglset ( 'LROT', .TRUE. )
C            call sgiset ( 'IROT', 90 )

            call SGSPLI( li )
            call sgplu(izz, kturbz, zgrid)

            call sglset( 'LCHAR', .FALSE. )
            call sglset ( 'LROT', .FALSE. )

            li = li + 10
            
         end if
      end if
      
      end do

      header='day=' // date
      call uxsttl ( 'T', header, 1. )

      call grcls
      call FCLOSE ( irun+1 )

      end