C*********************************************************************
C
C     PROGRAM xzcont_q
C
C     x-z平面等値線図作成プログラム(時刻, 条件指定)
C
C     履歴 1997/07/24 中島健介
C          1999/02/04 小高正嗣
C          2000/11/21 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      program xzcont_q
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------
      real*4 cdp                  ! カラートーン刻み
      real*4 tlev1,tlev2          ! カラーレベル
      integer ipat              ! トーン番号
      integer cmax
      parameter ( cmax = 60 )
      real*4 qvmax, qvmin
C      parameter ( qvmax = 1.0e-4, qvmin = 0.0 )
C      parameter ( qvmax = 1.0e-4, qvmin = 1.0e-5 )

      real*4 pti(2,cmax)
*-----------------------------------------------------

#include "grid_size_M.f"

      parameter ( nbr3 = 3*((nx+6)*(nz+6)/16+1) )
C      parameter ( nxp = nx, nzp= 100)
      parameter ( nxp = 256, nzp= 102)
C      parameter ( nzdiff = 105 )
*-----------------------------------------------------
      integer ibr(nbr3)
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
      parameter ( PI = 3.141592D0 ) 

      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
     \     w ( -2:nxb, -2:nzb ),
     \     qvap ( -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
     \      wp ( nx+1,nz+1 ),
     \      qvapp ( nxp+1, nzp+1 )

      real*4 z0pp( 0:nzp )


      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)   = REAL(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 DATA (I)' 
      write(6,*) '1; QVAP 2; QCLW 3; QRAI ?'
      read (*,*) IDN

      call sgpwsn
      read(*,*) iws
      CALL SWISET('IWIDTH',  640)
      CALL SWISET('IHEIGHT', 300)
C      call swlset( 'LDUMP', .TRUE. )
C      call swlset( 'LWAIT', .FALSE. )

      call gropn(iws)
C      call sldiv( 'T', 2, 3 )
C      call sldiv( 'T', 1, 3 )
      call swlset( 'LALT',  .TRUE. )
      call sglset( 'LFULL', .TRUE. )
      call sglset( 'LCORNER', .FALSE. )
      call udlset( 'LMSG', .FALSE. )
      call uzfact(0.6)
*
      qvmax = 1.0e-4
      qvmin = 1.0e-8
      cdp = (qvmax - qvmin)/cmax

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

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

      do it = 1, ite

         if (IDN .eq. 1) then
            read(35) time, ((qvap(ix,iz),ix=-2,nxb), iz=-2,nzb)
         else if (IDN .eq. 2 ) then
            read(36) time, ((qvap(ix,iz),ix=-2,nxb), iz=-2,nzb)
         else if (IDN .eq. 3 ) then
            read(37) time, ((qvap(ix,iz),ix=-2,nxb), iz=-2,nzb)
         else
            write(*,*) 'IDN is not correct; Stop.'
            stop
         end if

      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** 2-grid noize の除去 ***************************************
            
            do iz = 1,nzp+1
               do ix = 1, nxp+1
                  qvapp(ix,iz) = ( qvap(2*(ix-1)-1,2*(iz-1))
     \                 + 2.0 * qvap(2*(ix-1),2*(iz-1))
     \                 + qvap(2*(ix-1)+1,2*(iz-1)) ) / 4.0 
               end do
               z0pp(iz-1) = z0p( 2*(iz-1) ) 
            end do
            z0pp(nzp) = z0p(nz)
*
c**************************************************************

            call grfrm
            call grswnd ( 0.0, 51.2, 0.0, 20.0 )
            call grsvpt ( 0.1, 0.9, 0.1, 0.42 )
C            call grsvpt ( 0.15, 0.95, 0.1, 0.42 )
            call grstrn(1)
            call uspfit
            call grstrf

            call udpset ( 'rsizet', 0.01 )
            call ussttl ( 'Dust Mixing Ratio[kg/kg]', 
     \           'X[km]', '', 'Z[km]' )
C            call udpset ( 'ICYCLE', 5 )
            call uwsgya ( z0pp, nzp+1 )
C            call uetone ( qvapp, nxp+1, nxp+1, nzp+1 )
            call uetonf ( qvapp, nxp+1, nxp+1, nzp+1 )

C            call udgcla (  0.0, 5.0e-5, 1.0e-5 )
C            call udcntz ( qvapp, nxp+1, nxp+1, nzp+1, ibr, nbr3)

            call usdaxs
            call uxsttl ( 'T', header, 1. )

C*********************************************************************
C トーンバー
C*********************************************************************
            call grfig
            call grswnd( 0.0, 1.0,  qvmin, qvmax )
            call grsvpt( 0.95, 0.98,  0.1,  0.42 )
            call grstrn( 1 )
            call grstrf
            call uwsgyb( qvmin, qvmax, cmax )
            call uetone( pti, 2, 2, cmax )

            call slpvpr( 3 )
            call uyaxdv( 'L', 1.0e-5, 2.0e-5 )

*
*
         end if
      end if
 
      end do

      call grcls

      call FCLOSE

      end