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

#include "grid_size_M.f"

      parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) )
*-----------------------------------------------------
      integer ibr(nbr3)
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
      parameter ( PI = 3.141592D0 ) 

      real cdp_uw,cdp_pt,cdp_k,cdp_ptd ! カラートーン刻み
      real tlev1,tlev2          ! カラーレベル
      integer ipat              ! トーン番号
      integer cmax
      parameter ( cmax = 60 )   ! カラー数

      real wmax, wmin, umax, umin, ptmax, ptmin
      real ptdmax, ptdmin, kmax, kmin

      parameter ( wmax = 25., wmin = -25., umax = 25., umin = -25. )
C      parameter ( wmax = 15., wmin = -15., umax = 15., umin = -15. )
      parameter ( ptmax = 280., ptmin = 180. )
      parameter ( ptdmax = 4.0, ptdmin = -4.0 ) 
C      parameter ( ptdmax = 2.5, ptdmin = -2.5 ) 
      parameter ( kmax = 50., kmin = 1.e-5 )

      real*4 pti_u(2,cmax+1), pti_w(2,cmax+1)
      real*4 pti_pt(2,cmax+1), pti_ptd(2,cmax+1)
      real*4 pti_k(2,cmax+1)
*-----------------------------------------------------
      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
     \     u ( -2:nxb, -2:nzb ),
     \     w ( -2:nxb, -2:nzb ),
     \     tpot ( -2:nxb, -2:nzb ),
     \     kturb ( -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 
     \      up    ( nx+1,nz+1 ),
     \      wp    ( nx+1,nz+1 ),
     \      tpotp ( nx+1, nz+1 ),
     \      tpall ( nx+1, nz+1 ),
     \      tpave ( nz+1 ),
     \      tpdev ( nx+1, nz+1 ),
     \      kturbp( nx+1, nz+1 )

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

      character*100 dumy
      character CSGI*1

*-----------------------------------------------------
      
      open (11,file='exparam')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read (11,*) ylat, als

      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)/1000.0
         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
      write(6,*) 'SELECT TEMP. MODE (I)' 
      write(6,*) '1; RAW PTEMP 2; DEVIATION ?' 
      read (5,*) IDN

      call sgpwsn
      read(*,*) iws

      CALL SWISET('IWIDTH',  900)
      CALL SWISET('IHEIGHT', 450)
      call swlset( 'LDUMP', .TRUE. )
      call swlset( 'LWAIT', .FALSE. )

      call gropn(iws)
      call sldiv( 'T', 2, 2 )
c      call gropn(-iws)
c      call sldiv( 'T', 1, 3 )
      call swlset( 'LALT',  .TRUE. )
      call sglset( 'LFULL', .TRUE. )
      call sglset( 'LCORNER', .FALSE. )
      call uzfact(0.7)
      call udlset( 'LMSG', .FALSE. )
      call udiset( 'INDXMJ', 2 )


      cdp_uw = (wmax - wmin)/cmax
      cdp_pt = (ptmax - ptmin)/cmax
      cdp_k = (kmax - kmin)/cmax
      cdp_ptd = (ptdmax - ptdmin)/cmax

      do k=1,cmax+1
        pti_w(1,k) = wmin + (k-1)*cdp_uw
        pti_w(2,k) = wmin + (k-1)*cdp_uw
        pti_pt(1,k) = ptmin + (k-1)*cdp_pt
        pti_pt(2,k) = ptmin + (k-1)*cdp_pt
        pti_k(1,k) = kmin + (k-1)*cdp_k
        pti_k(2,k) = kmin + (k-1)*cdp_k
        pti_ptd(1,k) = ptdmin + (k-1)*cdp_ptd
        pti_ptd(2,k) = ptdmin + (k-1)*cdp_ptd
      end do

      do it = 1, ite

      read(31) time, ((u(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(33) time, ((w(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb)
      read(38) time, ((kturb(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 へ************************************
            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**************************************************************
c**************************************************************

            do iz = 1,nz+1            
               do ix = 1, nx+1
                  up(ix,iz)    = u(ix-1,iz-1)

                  wp(ix,iz)    = ( w(ix-2,iz-1) +
     \                 2.0*w(ix-1,iz-1) + w(ix,iz-1))/4.0

                  kturbp(ix,iz) = kturb(ix-1,iz-1)

                  tpotp(ix,iz) = ( tpot(ix-2,iz-1) +
     \                 2.0*tpot(ix-1,iz-1) + tpot(ix,iz-1))/4.0

                  tpall(ix,iz) = tpotp(ix,iz) + 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


            if ( IDN .eq. 1 ) then 

C*********************************************************************
C pt
C*********************************************************************

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

            call udiset( 'ICYCLE', 5 )
            call grfrm
            call grfig
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call uwsgya ( z0p, nz+1 )
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( CSGI(135)//'[K]', 'X[km]', '', 'Z[km]' )
            call uetonf ( tpall, nx+1, nx+1, nz+1 )
            call usdaxs
            call uxsttl ( 'T', header, 1. )
C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call uwsgyz( .FALSE. )
            call grswnd( 0.0, 1.0,  ptmin, ptmax )
            call grsvpt ( 0.95, 0.98, 0.1, 0.42)
            call grstrn( 1 )
            call grstrf
            call uetone( pti_pt, 2, 2, cmax+1 )
            call slpvpr( 3 )
            call uyaxdv( 'L', 2., 10. )
            call uzinit
            call ueitlv

            else 
C*********************************************************************
C ptemp_dev
C*********************************************************************

            call udiset( 'ICYCLE', 10 )

            TLEV2 = ptdmin + cdp_ptd 
            TLEV1 = -10.0
            call uestlv(TLEV1, TLEV2, 30999)

            do i = 2,cmax-1
               TLEV2 =  ptdmin + i*cdp_ptd 
               TLEV1 =  TLEV2 - cdp_ptd
               ipat  = (30+(i-1))*1000 + 999
               call uestlv(TLEV1, TLEV2, IPAT)
            end do

            TLEV2 = 10.0
            TLEV1 = ptdmax - cdp_ptd 
            call uestlv(TLEV1, TLEV2, 89999)

            call grfrm
            call grfig
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call uwsgya ( z0p, nz+1 )
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( CSGI(131)
     \           //CSGI(135)//'[K]', 'X[km]', '', 'Z[km]' )
            call uetonf ( tpdev, nx+1, nx+1, nz+1 )
            call usdaxs
            call uxsttl ( 'T', header, 1. )

c            call udgcla ( -10.0, 10.0, 2.0 )
c            call uddclv ( 2.0 )
c            call uddclv ( -2.0 )
c            call udsclv ( 0.0, 1, 1, '', '0.0' )
c            call udcntz ( tpdev, nx+1, nx+1, nz+1, ibr, nbr3)



C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call uwsgyz( .FALSE. )
            call grswnd( 0.0, 1.0,  ptdmin, ptdmax )
            call grsvpt ( 0.95, 0.98, 0.1, 0.42)
            call grstrn( 1 )
            call grstrf
            call uetone( pti_ptd, 2, 2, cmax+1 )
            call slpvpr( 3 )
            call uyaxdv( 'L', 0.2, 1.0 )
            call uzinit
            call ueitlv

            end if
*
*
C*********************************************************************
C kturb
C*********************************************************************

            do i = 1,cmax
               TLEV2 =  kmin + i*cdp_k 
               TLEV1 =  TLEV2 - cdp_k
               ipat  = (30+(i-1))*1000 + 999
               call uestlv(TLEV1, TLEV2, IPAT)
            end do

            call grfrm
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call uwsgya ( z0p, nz+1 )
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( 'kturb[m^2/s]', 'X[km]', '', 'Z[km]' )
            call uetonf ( kturbp, nx+1, nx+1, nz+1 )
            call usdaxs
            call uxsttl ( 'T', header, 1. )

C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call uwsgyz( .FALSE. )
            call grswnd( 0.0, 1.0,  kmin, kmax )
            call grsvpt( 0.95, 0.98,  0.10,  0.42 )
            call grstrn( 1 )
            call grstrf
            call uetone( pti_k, 2, 2, cmax+1 )
            call slpvpr( 3 )
            call uyaxdv( 'L', 5., 10. )
            call uzinit
            call ueitlv

*
*
C*********************************************************************
C w
C*********************************************************************
            call udiset( 'ICYCLE', 4 )

            TLEV2 =  wmin
            TLEV1 =  -50.0
            call uestlv(TLEV1, TLEV2, 30999)

            cdp_uw = (wmax - wmin)/(cmax-2)
            do i = 1,cmax-2
               TLEV2 =  wmin + i*cdp_uw
               TLEV1 =  TLEV2 - cdp_uw
               ipat  = (30+i)*1000 + 999
               call uestlv(TLEV1, TLEV2, IPAT)
            end do

            TLEV2 =  50.0
            TLEV1 =  wmax
            call uestlv(TLEV1, TLEV2, 89999)

            call grfrm
            call grfig
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call uwsgya ( z0p, nz+1 )
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( 'W[m/sec]', 'X[km]', '', 'Z[km]' )
            call uetonf ( wp, nx+1, nx+1, nz+1 )
            call usdaxs
            call uxsttl ( 'T', header, 1. )

            call udgcla ( -50.0, 50.0, 5.0 )
            call udsclv ( 0.0, 1, 1, '', '0.0' )
            call udcntz ( wp, nx+1, nx+1, nz+1, ibr, nbr3)

C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call uwsgyz( .FALSE. )
            call grswnd( 0.0, 1.0,  wmin, wmax )
            call grsvpt ( 0.95, 0.98, 0.1, 0.42 )
            call grstrn( 1 )
            call grstrf
            call uetone( pti_w, 2, 2, cmax+1 )
            call slpvpr( 3 )
            call uyaxdv( 'L', 5., 10. )
            call uzinit

C*********************************************************************
C u
C*********************************************************************

            call grfrm
            call grfig
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
            call grstrn (1)
            call uspfit
            call grstrf
            call uwsgya ( z0p, nz+1 )
            call udpset ( 'rsizet', 0.01 )
            call ussttl ( 'U[m/sec]', 'X[km]', '', 'Z[km]' )
            call uetonf ( up, nx+1, nx+1, nz+1 )
            call usdaxs
            call uxsttl ( 'T', header, 1. )

            call udgcla ( -50.0, 50.0, 5.0 )
            call udsclv ( 0.0, 1, 1, '', '0.0' )
            call udcntz ( up, nx+1, nx+1, nz+1, ibr, nbr3)

C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call uwsgyz( .FALSE. )
            call grswnd( 0.0, 1.0,  umin, umax )
            call grsvpt ( 0.95, 0.98, 0.1, 0.42)
            call grstrn( 1 )
            call grstrf
            call uetone( pti_w, 2, 2, cmax+1 )
            call slpvpr( 3 )
            call uyaxdv( 'L', 5., 10. )
            call uzinit
            call ueitlv
*
*


         end if
      end if

      end do

      call grcls

      call FCLOSE

      end
