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
      real*8 kturbread(-2:nxb,-2:nzb)

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

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

      character*4 xpos

      real*4 para1
      integer para2
      integer para3
      integer para4

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

      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

      do iz = 1, nz+1
         zgrid(iz) = z0(iz-1) / 1000.
      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,*) 'snapshot location = ?(1-',nx,')'
      write(6,*) 'ix1, ix2, ix3'
      read (5,*) ix1, ix2, ix3

      call sgpwsn
      read(*,*) iws

      call gropn(iws)
      call sldiv( 'Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )

      call grfrm
      call grswnd( 0., 35., 0., 1. )
C      call grswnd( 0., 25., 0., 20. )
      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

      para1 = 0.01
      para2 = 1
      para3 = 1
      para4 = 1

      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. )
            call sgiset ( 'IROT', 90 )
            call sgrset ( 'FFCT', para1 )
            
C            call uulinz(nz+1, kturbz, zgrid, para3, para4)
            call uulinz(15, kturbz, zgrid, para3, para4)

            call sglset( 'LCHAR', .FALSE. )
            call sglset ( 'LROT', .FALSE. )
            para1 = para1 + 0.03
            para2 = para2 + 1
            para3 = para2 - INT((para2-1)/4)*4
            para4 = para4 + INT((para2-1)/4.)
            
            if ( it .eq. ite ) then

               call grfrm
               call grswnd( 0., 50., 0., 20. )
               call grsvpt( 0.2, 0.8, 0.2, 0.8 )
               call grstrn(1)
               call grstrf
               call ussttl ( 'Kturb', '[m^2/s]', '', 'Z[km]' )
               call usdaxs
               call uxsttl ( 'T', header, 1. )

               call uulinz(nz+1, kturbz, zgrid, 1, 2)

               do iz = 1,nz+1
                  kturbz2(iz) = kturb(ix1,iz)
                  kturbz3(iz) = kturb(ix2,iz)
                  kturbz4(iz) = kturb(ix3,iz)
               end do
               
               call sglset( 'LCHAR', .TRUE. )
               call sgrset ( 'FFCT', 0.5 )

               write(xpos,400) REAL((ix1-1)/10.)
 400           format(F4.1)
               call sgsplc( xpos ) 
               call uulinz(nz+1, kturbz2, zgrid, 2, 1)
               
               write(xpos,400) REAL((ix2-1)/10.)
               call sgsplc( xpos ) 
               call uulinz(nz+1, kturbz3, zgrid, 3, 1)
         
               write(xpos,400) REAL((ix3-1)/10.)
               call sgsplc( xpos ) 
               call uulinz(nz+1, kturbz4, zgrid, 4, 1)

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

            
         end if
      end if
      
      end do

      call grcls
      call FCLOSE ( irun+1 )

      end