C*********************************************************************
C
C     PROGRAM xzcont_pt
C
C     x-z平面等値線図作成プログラム(時刻, 条件指定)
C
C     履歴 1997/07/24 中島健介
C          1997/07/25 小高正嗣
C          1997/07/26 小高正嗣
C          1997/07/26 小高正嗣
C          1998/09/17 小高正嗣
C          1998/10/02 小高正嗣
C          1999/01/10 小高正嗣
C          1999/02/04 小高正嗣
C          1999/02/07 小高正嗣
C
C*********************************************************************
      program xzcont_pt
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------

#include "grid_size_M.f"

      parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) )
*-----------------------------------------------------
      integer ibr(nbr3)
*-----------------------------------------------------
      real cdp                  ! カラートーン刻み
      real tlev1,tlev2          ! カラーレベル
      integer ipat              ! トーン番号
      integer cmax
      parameter ( cmax = 60 )
      real ptmax, ptmin
      parameter ( ptmax = 3.0, ptmin = -3.0 )

      real*4 pti(2,cmax)
*-----------------------------------------------------
      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*8
     \     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 
     \      tpotp ( nx+1, nz+1 ),
     \      tpall ( nx+1, nz+1 ),
     \      tpave ( nz+1 ),
     \      tpdev ( nx+1, nz+1 )

      character*6 stime
c      character*2 stime
      character*11 ttime
      character*100 dumy
      character*8 ltime
*-----------------------------------------------------
      real*4 xvpt
      real*4 yvpt
      real*4 tpdev_p(nz+1, nz+1)
*-----------------------------------------------------
      xvpt = 0.05
      yvpt = 0.4

      open (11,file='exparam.1')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu

      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)   = zgrid(iz)
         p0p(iz)   = pres0(iz)
         d0p(iz)   = dens0(iz)
         temp0p(iz)= temp0(iz)
         tpot0p(iz)= tpot0(iz)
         pai0p(iz) = ppai0(iz)
         rvap0p(iz)= rvap0(iz)
         qvap0p(iz)= qvap0(iz)
      end do

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

      cdp = (ptmax - ptmin)/cmax
      do i = 1,cmax
         TLEV2 =  ptmin + i*cdp 
         TLEV1 =  TLEV2 - cdp
         ipat  = (30+(i-1))*1000 + 999
         call uestlv(TLEV1, TLEV2, IPAT)
      end do

      do k=1,cmax+1
        pti(1,k) = ptmin + (k-1)*cdp
        pti(2,k) = ptmin + (k-1)*cdp
      end do

      iflag = 0

      call sgpwsn
      read(*,*) iws

      call gropn(iws)
      call swlset( 'LALT',  .TRUE. )
      call sglset( 'LFULL', .TRUE. )
      call sglset( 'LCORNER', .FALSE. )
      call uzfact(0.5)
      call UDLSET( 'LMSG', .FALSE. )
*
      do it = 1, ite

      read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)

      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 へ************************************
c            time = time - 86400.*2
c            time = time / 3600. + 6
c**************************************************************

            write(stime,100) INT(time)
 100        format(I6.6)
c 100        format(I2.2)
c            ltime = 'LT=' // stime // ':00'
            ttime = 't=' // stime // 'sec'

            do iz = 1,nz+1
               do ix = 1, nx+1
                  tpotp(ix,iz) = tpot(ix-1,iz-1)
                  tpall(ix,iz) = tpot(ix-1,iz-1) + tpot0(iz-1)
               end do
            end do
*
            do iz = 1, nz+1
               tpave(iz) = 0.
               do ix = 1,nx+1
                  tpave(iz) = tpave(iz) + tpall(ix,iz) /(nx+1)
               end do
            end do

            do iz = 1,nz+1
               do ix = 1, nx+1
                  tpdev(ix,iz) = tpall(ix,iz) - tpave(iz) 
               end do
            end do

            do iz = 1,nz+1
               do ix = 1, nz+1
                  tpdev_p(ix,iz) = tpdev(ix,iz)
               end do
            end do

            if (iflag .eq. 4) then
               iflag = 1
               xvpt = 0.05
            else 
               iflag = iflag + 1
            end if

            if ( iflag .eq. 1 ) then
               call grfrm
            else
               call grfig
            end if

            call grswnd ( 0.0, 10.0, 0.0, 10.0 )
c            call grsvpt ( 0.15, 0.95, 0.25, 0.406 )
            call grsvpt ( xvpt, xvpt+0.2, 0.25, 0.45 )
            call grstrn (1)
            call uspfit
            call grstrf
            call udpset ( 'rsizet', 0.01 )
c            call ussttl ( 'anomaly THETA[K]', 'X[km]', '', 'Z[km]' )
            call udpset ( 'ICYCLE', 2 )

c            call uestlv ( -10., 0., 603)
c            call uetone ( tpdev_p, nz+1, nz+1, nz+1 )
            call uetonf ( tpdev_p, nz+1, nz+1, nz+1 )
c            call ueitlv

            if ( iflag .eq. 1 ) then
               call ussttl ( '', '', '', 'Z[km]' )
               call usdaxs
            else if ( iflag .eq. 4 ) then
c               call slpvpr( 3 )
c               call uxaxdv( 'B', 1., 2. )
c               call uxaxdv( 'T', 1., 2. )
c               call uyaxdv( 'L', 1., 2. )
c               call uyaxdv( 'R', 1., 2. )
               call ussttl ( '', 'X[km]', '', '' )
               call usdaxs
            else
               call slpvpr( 3 )
               call uxaxdv( 'B', 1., 2. )
               call uxaxdv( 'T', 1., 2. )
               call uyaxdv( 'L', 1., 2. )
               call uyaxdv( 'R', 1., 2. )
            end if

c            if ( it .eq. 2 ) then
c               call udgcla ( -0.01, 0.01, 0.005 )
c               call uddclv(0.)
c            else

            call udgcla ( -4., 4., 0.5 )

c            end if

c            call udcntz ( tpdev_p, nz+1, nz+1, nz+1, ibr, nbr3)
            call uxsttl ( 'T', ttime, 1. )
c            call uxsttl ( 'T', ltime, 1. )
            call uzinit
*
*
            xvpt = xvpt + 0.23

         end if
      end if

      end do

      call grcls

      call FCLOSE

      end
