C*********************************************************************
C
C     PROGRAM PTDVLOP_BL
C
C     水平平均温度・温位鉛直分布作成プログラム(時間発展)
C
C     1999/06/15 小高正嗣
C     2000/11/21 小高正嗣; D 論用に書き直し
C     2000/11/29 小高正嗣; 境界層用
C
C*********************************************************************
      program ptdvlop_bl
C---------------------------------------------------------------------  
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C---------------------------------------------------------------------  
#include "grid_size_M.f"
      parameter ( nzprof = 100 )                              
      parameter ( nzp = nz+1   )
C---------------------------------------------------------------------  
C                                                                       
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )

      integer fnum
      character*5  nmrun
      character*4  from
      character*1  chara
      character*1  pz
      character*2  end

      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 )

      real*4  dtpad0,dqvad0
      real*4  tptend(nzp),tpmean(nzp),dtpadv(nzp),
     &        dtpdif(nzp),dtpdnl(nzp),
     &        dtpdi0(nzp),dtprad(nzp),tperror(nzp),
     &        dtirad(nzp),dtsrad(nzp)

      real*4 tmean(nzp), zgridp(nzp), zmax, Tmin, Tmax
      real*4 ekvis(nzp), eknlv(nzp), ekbuoy(nzp)

      real*4 line(2)
      real*4 temp(2)

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

      character*100 dumy
      character CSGI*1
C
C=====================================================================  

      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

      call FCLOSE ( irun+1 )

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

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

      do iz = 1, izz
         zgridp(iz) = zgrid(iz-1) / 1000.
      end do

      zmax = zgridp(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
C      call grswnd( 200., 250., 0., zmax )
C      call grswnd( 190., 290., 0., zmax )
      call grswnd( Tmin, Tmax, 0., zmax )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn(1)
      call grstrf
      call ussttl( CSGI(135), 'K', '', 'Z[km]' ) 
      call usdaxs
      CALL sglset( 'LCHAR', .TRUE. )

      do i = ibegin, iend, inter

         type = type + 0.1

         fnum = 50 + i
         open(fnum)

!********************************************************************
! * kinetic energy production & dissipation *

         read(fnum,'(a2,a15,g15.7,a6,g15.7,a4)') 
     &        chara, nmrun, t_begin, from, t_end,
     &        pz

         do iz = 1, nzp
            read(fnum,*) zgrid(iz), 
     &           ekvis(iz),
     &           eknlv(iz),
     &           ekbuoy(iz)
         end do

!********************************************************************

         read(fnum,'(a2,a15,g15.7,a6,g15.7,a4)') 
     &        chara, nmrun, t_begin, from, t_end,
     &        pz

         do iz = 1, nzp
            read(fnum,'(10g15.7)') zgrid(iz), 
     &           tperror(iz),
     &           tptend(iz),
     &           tpmean(iz),
     &           dtpadv(iz),
     &           dtpdif(iz),
     &           dtpdnl(iz),
     &           dtpad0,
     &           dtpdi0(iz),
     &           dtprad(iz),
     &           dtirad(iz),
     &           dtsrad(iz)
         end do

!********************************************************************
         write(date, '(I2.2)') INT( time0 / day + 1 )

         if ( t_begin .eq. 0.0 ) then
            t_begin = t_end - tinterv
         end if

         t_begin = t_begin + times
         t_end   = t_end   + times

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

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

         header = 'day=' // date // ' ' // 'LT=' // lt // ':' // mint
!********************************************************************
C         CALL sgrset( 'FFCT', 0.8 )
C         CALL sgsplc( lt )
         call SGSPLI( 43 )
         call sgplu(izz, tpmean, zgridp)

         line(1) = 0.05
         line(2) = 0.05
         temp(1) = Tmin
         temp(2) = Tmax

         call SGSPLI( 23 )
         call sgplu(2, temp, line)

         line(1) = 0.4
         line(2) = 0.4
         temp(1) = Tmin
         temp(2) = Tmax

         call SGSPLI( 153 )
         call sgplu(2, temp, line)

         close(fnum)

      end do

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




    
      call grcls

      end