!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(水平平均温位)
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!=================================================================
      program ttend
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = NZ )

      dimension
     \   TEMP0(-2:NZB) , TPOT0(-2:NZB),
     \   CDTURB(-2:NZB), ZGRID(-2:NZB),
     \   DTPRAD(-2:NZB), DTPDI0(-2:NZB),
     \   PRES0(-2:NZB) , DENS0(-2:NZB),
     \   PPAI0(-2:NZB) , QVAP0(-2:NZB), RVAP0(-2:NZB)

      real*4
     \   ZGRID_p(0:NZP),
     \   QVAL1_p(0:NZP) ,
     \   QVAL2_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), PRES0(IZ), DENS0(IZ),
     \                   TEMP0(IZ), TPOT0(IZ), PPAI0(IZ), 
     \                   QVAP0(IZ), RVAP0(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 select'
      write(6,*) '1: IR heating rate'
      write(6,*) '2: SR heating rate'
      write(6,*) '3: dust IR heating rate'      
      write(6,*) '4: dust SR heating rate'      
      read(*,*) idata


      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
c      CALL grswnd( -10.0, 10.0, ZMIN, ZMAX )   
      CALL grswnd( -200.0, 200.0, ZMIN, ZMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'THETA TENDENCY', '[K/day]', '', 'Z[km]' )
      CALL usdaxs


      do it = 1, ite

         if ( idata == 1 ) then
            read(34) TIME, ( DTPRAD(IZ), IZ=-2,NZB )
         else if ( idata == 2 ) then
            read(37) TIME, ( DTPRAD(IZ), IZ=-2,NZB )
         else if ( idata == 3 ) then
            read(38) TIME, ( DTPRAD(IZ), IZ=-2,NZB )
         else if ( idata == 4 ) then
            read(39) TIME, ( DTPRAD(IZ), IZ=-2,NZB )
         else 
            write(*,*) 'data number error: stop'
            stop
         end if

c         read(35) TIME, ( DTPDI0(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
         
!*** 時刻変換: 秒から LT へ************************************
               time = time - 86400.* ( irun -1 )
               time = time / 3600. + 6

               if ( time .gt. 24 ) then
                  time = time - 24
               end if
!**************************************************************

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

!               CALL uuslnt( 3 )
!               CALL uuslni( 2 )
               CALL uulinz( NZP, QVAL1_p, ZGRID_p, 1, 2 )
!               CALL uulinz( NZP, QVAL1_p, ZGRID_p, 1, 0.4 )
!               CALL uuslnt( 1 )
!               CALL uulinz( NZP, QVAL2_p, ZGRID_p, 1, 0.45 )

!               CALL uxsttl( 'T', header, 1. )      

            end if
         end if

      end do

      header = '1D model'
      CALL uxsttl( 'T', header, 1. )      

      CALL grcls                                                        


      stop
      end

