C*********************************************************************
C
C     PROGRAM tgrand
C
C     地中・地表温度日変化図作成プログラム(時刻, 条件指定)
C
C     履歴 2000/09/22 小高正嗣
C          2000/11/21 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      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
     \     TPOT  (-2:NXB,-2:NZB),
     \     TEMPG (-2:NXB,NZG),  ! 地中温度
     \     GGRID (-1:NZGB),     ! 鉛直座標
     \     GDZ   (0:NZGB),      ! 鉛直格子間隔(無次元)
     \     GDZM  (0:NZG)        ! 鉛直格子間隔(半整数格子上)

      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*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*5 CLAT
      character*5 CLS
      character*2 mint
      character*2 lt
      character*2 date
      character*1 CSGI
      character*100 dumy

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

      open(10,file='exparam')
      read(10,*) IRUN
      read(10,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read(10,*) YLAT, ALS

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

c*** 時刻変換: 日付設定 ************************************
      time = time0 + times
      write(date, '(I2.2)') INT( time/day ) + 1
      dheader = 'day=' // date
c**************************************************************

      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

      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 /

      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 SWISET('IWIDTH', 840)
      CALL SWISET('IHEIGHT', 420)
      call sglset( 'LFULL', .TRUE. )
      call swlset( 'LDUMP', .TRUE. )

      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( 'Grand Temperature', 'K', 'Depth', 'Z' )
      call usdaxs

      li = 23

      do it = 1, ite

         read(40) TIME, ((TEMPG(IX,IZ), IX=-2,NXB),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

c*** 時刻変換: 秒から LT へ************************************
            time = time + times
            time = ( time - INT( time / day ) * day ) / hour

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

            write(lt, '(I2.2)') INT(time)
c**************************************************************

            do IZ = 1, NZG
               TEMPG_p(IZ) = 0.0
               end do

               do IZ = 1, NZG
                  do IX = 1, NX
                     TEMPG_p(IZ) = TEMPG_p(IZ) + TEMPG(IX,IZ) / NX
                  end do
               end do

               CALL sgrset( 'FFCT', 0.05 )
               CALL sgsplc( lt )
               call SGSPLI( li )
               call sgplu(nzp, tempg_p, ggrid_p)
               li = li + 10

            end if
         end if

      end do

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

      call FCLOSE

      call FOPEN ( irun+1 )

      ilate = INT(TIMES / HOUR )

      do it = 1, nloop1

         read(34) TIME, ((TPOT(IX,IZ), IX=-2,NXB),IZ=-2,NZB)
         read(40) TIME, ((TEMPG(IX,IZ),IX=-2,NXB),IZ= 1,NZG)

         TEMPA(it+ilate) = 0.0
         TSFC (it+ilate) = 0.0

         do IX = 1, NX
            TEMPA(it+ilate) = TEMPA(it+ilate) +
     \           ( TPOT(IX,0) + TPOT0(0) ) / NX
            TSFC (it+ilate) = TSFC (it+ilate) +
     \           TEMPG(IX,NZG) / NX
         end do

         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',
     \     'Grand Surface Temp.', 'K' )
      call usdaxs

      call SGSPLI( 22 )
      call sgplu(nloop1+1, TIMEL, TSFC)
      call SGSPLI( 42 )
      call sgplu(nloop1+1, TIMEL, TEMPA)

      call uxsttl( 'T', dheader, 1. )
c      call uxsttl( 'T', pheader, 1. )

      call grcls

      call FCLOSE

      end