C********************************************************************
C
C     PROGRAM QVTEND
C     
C     水平平均混合非増加率鉛直分布図作成プログラム
C
C     履歴 1999/06/15 小高正嗣
C          2000/11/21 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      program qvt
C---------------------------------------------------------------------  
      implicit real*8 ( a-h, o-z )
C---------------------------------------------------------------------  
#include "grid_size_M.f"
      parameter ( nzprof = 100 )
      parameter ( nzp = nz+1   )
C---------------------------------------------------------------------  
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
C                                                                       
      character*5  nmrun
      character*4  from
      character*1  chara
      character*1  pz

      character*2  start
      character*2  finish
      character*2  date      
      character*3  qrad
      character*20 header

      real*8  ekvis(nzp), eknlv(nzp), ekbuoy(nzp)

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

      real*8  qvtend(nzp),qvmean(nzp),dqvadv(nzp),
     &        dqvdif(nzp),dqvdnl(nzp),
     &        dqvdi0(nzp),dqvnfl(nzp),qverror(nzp)

      real*4  qvtend_p(nzp),qvmean_p(nzp),dqvadv_p(nzp),
     &        dqvdif_p(nzp),dqvdnl_p(nzp),zgrid_p(nzp),
     &        dqvdi0_p(nzp),dqvnfl_p(nzp),qverror_p(nzp)

      real*4  QVMAX, QVMIN

      integer begin,end,interv,fnum
      integer nloop1

      character CSGI*1
      character title*5
*-----------------------------------------------------

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

*-----------------------------------------------------
      QVMAX = 5.0e-10
      QVMIN = - QVMAX

      title = 'dq/dt'

      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, z0(iz), p0(iz), d0(iz), temp0(iz),
     \                     tpot0(iz), ppai0(iz), dum1, dum2
      end do

      call FCLOSE

      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 )
      call sldiv ('Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )

      do it = 1, ite

      fnum = 50 + it
      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,*) zgrid(iz), 
     &                        tperror(iz),
     &                        tptend(iz),
     &                        tpmean(iz),
     &                        dtpadv(iz),
     &                        dtpdif(iz),
     &                        dtpdnl(iz),
     &                        dtpad0(iz),
     &                        dtpdi0(iz),
     &                        dtprad(iz),
     &                        dtirad(iz),
     &                        dtsrad(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,*) zgrid(iz),
     &                        qvtend(iz),
     &                        qvmean(iz),
     &                        dqvadv(iz),
     &                        dqvdif(iz),
     &                        dqvdnl(iz),
     &                        dqvad0,
     &                        dqvdi0(iz),
     &                        dqvnfl(iz),
     &                        qverror(iz)
      end do

      do iz = 1, nzp-1
         zgrid_p(iz) = REAL(zgrid(iz)) / 1000.

         qvtend_p(iz) = REAL(qvtend(iz))
         qvmean_p(iz) = REAL(qvmean(iz))
         dqvadv_p(iz) = REAL(dqvadv(iz) + dqvdnl(iz) )
         dqvdif_p(iz) = REAL(dqvdif(iz) + dqvdi0_p(iz) )
         qverror_p(iz) = REAL(qverror(iz))
      end do


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

      write(date, '(I2.2)') INT( t_begin / day + 1 )

      begin = INT( t_begin -  INT( t_begin / day ) * day )
      end   = INT( t_end -  INT( t_begin / day ) * day )

      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(start,'(I2.2)') INT(t_begin)
      write(finish,'(I2.2)') INT(t_end)

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

      header = 'day' // date // 
     \     ' LT=' // start // ':00-' // finish // ':00'
   
      close(fnum)

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

      call grfrm
      call grswnd( QVMIN, QVMAX, 0., 20. )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn(1)
      call grstrf
      call ussttl ( '', '', '', 'Z[km]' )
      call usdaxs
      call uxsttl ( 'T', header, 1. )
      call uxsttl ( 'B', title, 0. )
      call sglset( 'LCHAR', .TRUE. )

C      call sgsplc( 'tend' ) 
C      call uulinz(nzp-1, qvtend_p, zgrid_p, 3, 3)
C      call sgsplc( 'adv' ) 
C      call uulinz(nzp-1, dqvadv_p, zgrid_p, 1, 3)
C      call sgsplc( 'dif' ) 
C      call uulinz(nzp-1, dqvdif_p, zgrid_p, 2, 3)
c      call sgsplc( 'nfl' ) 
c      call uulinz(nzp-1, dqvnfl_p, zgrid_p, 4, 3)

      call sgsplc( 'tend' ) 
      call SGSPLI( 23 )
      call sgplu(nzp-1, qvtend_p, zgrid_p)
      call sgsplc( 'adv' ) 
      call SGSPLI( 33 )
      call sgplu(nzp-1, dqvadv_p, zgrid_p)
      call sgsplc( 'dif' ) 
      call SGSPLI( 43 )
      call sgplu(nzp-1, dqvdif_p, zgrid_p)
      call sgsplc( 'nfl' ) 
      call SGSPLI( 53 )
      call sgplu(nzp-1, dqvnfl_p, zgrid_p)

      call sglset( 'LCHAR', .FALSE. )

      call grfrm
      call grswnd( QVMIN, QVMAX, 0., 20. )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn(1)
      call grstrf
      call ussttl ( '', '', '', 'Z[km]' )
      call usdaxs
      call uxsttl ( 'T', header, 1. )
      call uxsttl ( 'B', title, 0. )
      call sglset( 'LCHAR', .TRUE. )

C      call sgsplc( 'tend' ) 
C      call uulinz(nzp-1, qvtend_p, zgrid_p, 3, 3)
      call sgsplc( 'adv' ) 
      call uulinz(nzp-1, dqvadv_p, zgrid_p, 1, 3)
      call sgsplc( 'dif' ) 
      call uulinz(nzp-1, dqvdif_p, zgrid_p, 2, 3)

      call sgsplc( 'nfl' ) 
      call uulinz(nzp-1, dqvnfl_p, zgrid_p, 4, 3)
      call sglset( 'LCHAR', .FALSE. )



      end if
      end if

      end do

      call grcls

      end



