!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(地中/地面温度)
!
!  2000/09/22  小高正嗣
!=================================================================
      program tgrand
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = NZG )
      parameter ( NT  = 100 )
      parameter ( DAY = 86400.0D0 )
      parameter ( HOUR= 3600.0D0 )
!-----------------------------------------------------------------      
      dimension TEMP0(-2:NZB)

      dimension
     \     TEMPG (NZG),         ! 地中温度
     \     GGRID (-1:NZGB),      ! 鉛直座標
     \     GDZ   (0:NZGB),     ! 鉛直格子間隔(無次元)
     \     GDZM  (0:NZG)       ! 鉛直格子間隔(半整数格子上)

      real*4
     \   GGRID_p(NZP),
     \   TEMPG_p(NZP) 

      real*4
     \   ZMIN , ZMAX,  TMIN,  TMAX
  
      real*4
     \   TSFC(0:NT), TIMEL(0:NT), TEMPA(0:NT)

      character*14 dheader
      character*18 pheader
      character*100 model

      character*5 CLAT
      character*5 CLS
      character*2 timell
      character*2 date
      character*2 CSGI


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

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

      write(CLS,  '(F5.1)') ALS
      write(CLAT, '(F5.1)') YLAT
      pheader = 'Ls=' // CLS // ' ' // CSGI(172) // '=' // CLAT

      call FOPEN ( irun+1 )

      DATA ( GGRID(I), I=-1,NZGB )
     1     / -6.0D0 ,
     2       -6.0D0 , -4.0D0 , -2.7D0 , -1.8D0, -1.2D0, 
     3       -0.79D0, -0.53D0, -0.35D0, -0.2D0, -0.1D0, 
     4        0.0D0 ,  0.0D0 /

      call SETCSG 
     C   ( ALBEDO, DIFG, CPGRND, DENSG, CDTURBG, SKDPS,
     G     GGRID , GDZ , GDZM, 
     D     NZG   , NZGB ) 

      do iz = 1, nzp
         GGRID_p(iz)= GGRID(iz)
      end do

      ZMIN    = GGRID_p(1)
      ZMAX    = 0.0
      TMIN    = 160.0
      TMAX    = 290.0

      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 sglset( 'LCHAR', .TRUE. )
      call sldiv ('Y', 2, 1 )

      call grfrm
      CALL sglset( 'LCHAR', .TRUE. )
      call grswnd( TMIN, TMAX, ZMIN, ZMAX )   
      call grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Temperature', '[K]', 'Depth', 'Z[km]' )
      call usdaxs

      do it = 1, ite


         read(36) TIME, ( TEMPG(IZ), IZ=1,NZG )
         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
               write(date, '(I2.2)') INT( time/DAY ) + 1

c*** 時刻変換: 秒から LT へ************************************
               time = MOD( time, DAY ) + TIMES
               time = time / HOUR

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

               write(timell,'(I2.2)') INT(time)
               dheader = 'day' // date // ' LT=' // timell // ':00'

               do iz = 1, nzp
                  TEMPG_p(iz)= TEMPG(iz)
               end do

               CALL sgrset( 'FFCT', 0.05 )
               CALL sgsplc( timell )
               call uulinz( NZP, TEMPG_p, GGRID_p, 1, 3  )

            end if
         end if

      end do

!      Call uxsttl( 'T', dheader, 1. )      
      CALL uxsttl( 'T', pheader, 1. )      
      CALL uxsttl( 'T', model, 1. )      

      call FCLOSE

      call FOPEN ( irun+1 )

      ilate = INT(TIMES / DTIME / NLOOP2 )
     
      do it = 1, nloop1
         read(32) TIME, ( TEMP0(IZ), IZ=-2,NZB )
         read(36) TIME, ( TEMPG(IZ), IZ=1,NZG )

         TEMPA(it+ilate) = TEMP0(0) 
         TSFC (it+ilate) = TEMPG(NZG) 
         TIMEL(it+ilate) = (MOD ( TIME, DAY ) + TIMES ) / HOUR
      end do

      TEMPA(ilate) = TEMPA(nloop1+ilate)
      TSFC(ilate)  = TSFC(nloop1+ilate)
      TIMEL(ilate) = TIMEL(nloop1+ilate)


      il = 0
      do it = 1, nloop1+ilate
         if ( TIMEL(it) .ge. 24.0 ) then
            TIMEL(il) = TIMEL(it) - 24.0
            TSFC (il) = TSFC(it)
            TEMPA(il) = TEMPA(it)
            il = il + 1
         end if
      end do

      call grfrm
      CALL sglset( 'LCHAR', .FALSE. )
      call grswnd( 0.0, 24.0, TMIN, TMAX )   
      call grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Local Time', '[Hour]', 'Temperature', '[K]' )
      call usdaxs
      call uulinz(nloop1+1, TIMEL, TSFC, 1, 3  )
      call uulinz(nloop1+1, TIMEL, TEMPA, 3, 3  )
               
!      Call uxsttl( 'T', dheader, 1. )      
      CALL uxsttl( 'T', pheader, 1. )      
      CALL uxsttl( 'T', model, 1. )      

      call grcls                                                        

      call FCLOSE


      end

