!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(水平平均温位)
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!=================================================================
      program cdprof

!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = 50 )

      dimension
     \   TEMP0(-2:NZB) , TPOT0(-2:NZB),
     \   CDTURB(-2:NZB), ZGRID(-2:NZB)

      real*4
     \   ZGRID_p(0:NZP),
     \   QVAL1_p(0:NZP) 

      real*4
     \   ZMIN , ZMAX

      character*14 header
      character*2 timel
      character*2 day

      NAMELIST /CONST1/ GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT
      NAMELIST /CONST2/ CDRAG, VSFC0, PSFC

!-----------------------------------------------------------------      

      open(10,file='exparam')
      read(10,*) IRUN
      read(10,*) TIME0, NLOOP1, NLOOP2, DTIME

      write(day,'(I2.2)') irun

      call FOPEN ( irun+1 )

      do IZ = -2, NZB
         READ(16,*) ITW, ZGRID(IZ), TEMP0(IZ), CDTURB(IZ)
      end do

      do iz = 0, nzp
         ZGRID_p(iz)= real( ZGRID(iz) )/ 1000.0
      end do

      ZMIN    = ZGRID_p(0)
      ZMAX    = ZGRID_p(NZP)

      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)
      CALL sglset( 'LCORNER', .FALSE. )

      CALL grfrm
      CALL grswnd( 0.0, 1.2e+5, ZMIN, ZMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'Kturb', '[]', '', 'Z[km]' )
      CALL usdaxs

      do it = 1, ite

         read(33) TIME, ( CDTURB(IZ), 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 - 86400.* ( irun -1 )
               time = time / 3600. + 6
               write(*,*) time
               if ( time .gt. 24 ) then
                  time = time - 24
               end if
c**************************************************************

               write(timel,'(I2.2)') INT(time)
               header = 'day' // day // ' LT=' // timel // ':00'
    
               do iz = 0, NZP
                  QVAL1_p(iz) = real( CDTURB(iz) )
               end do

               CALL uuslni( 2 )
               CALL uulin( NZP+1, QVAL1_p, ZGRID_p )
c               CALL uxsttl( 'T', header, 1. )      

            end if
         end if

      end do

      CALL grcls                                                        


      stop
      end

