!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(地表からの熱フラックス)
!
!  2000/02/29  小高正嗣
!=================================================================
      program hflux
!-----------------------------------------------------------------      
      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), QVAL1(-2:NZB)

      DIMENSION
     \     CDRAGV( -2:NXB ),
     \     RICHD( -2:NXB ),
     \     TSFC ( -2:NXB )

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

      real*4
     \   ZMIN , ZMAX

      character*14 header
      character*2 timel
      character*2 day
      
      parameter(PI = 3.141592D0)

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

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

      call SETCST( grav, psd, amdry, rdry, cpdry, rmvap, alatent )

      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., 24.0, -5., 20. )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'LT', '[hour]', 'heat flux', '[W/m^2]' )
      CALL usdaxs

      do it = 1, ite

         read(32) TIME, ( TEMP0(IZ), IZ=-2,NZB )
         read(31) TIME, ( TPOT0(IZ), IZ=-2,NZB )
c         write(0,*) 'reading time = ', time
c         write(0,*) 'writing time = ', time

         do IX = -2, NX+2
            TSFC(IX) = 180.0D0 + 1.0D+2 * 
     \           MAX( sin( 2.0d0 * PI * TIME / 8.64D+4 ), 0.0D0 )
         end do
        
c*** 時刻変換: 秒から LT へ************************************
         time = time - 86400.* ( irun -1 )
         time = time / 3600. + 6

         time = it + 6
         if ( time .gt. 24 ) then
            time = time - 24
         end if

         TIME_p(it) = real( time )
c**************************************************************
         write(timel,'(I2.2)') INT(time)
         header = 'day' // day // ' LT=' // timel // ':00'

         VSFC0 = 1.0D0

         call CLDRAG
     I   ( TPOT0, 
     O     RICHD, CDRAGV, 
     G     ZGRID, 
     C     GRAV , TSFC , VSFC0,
     D     NX   , NXB  , NZB   )

         dens0 = psd / RDRY / TEMP0(0)
         FLUX(it) = - CDRAGV(0) * VSFC0 * CPDRY *dens0
     \        * ( TEMP0(0) - TSFC(0) )

      end do

      CALL uuslni( 2 )
      CALL uumrkz( NZP, TIME_p, FLUX, 2, 1, 0.01  )

      CALL grcls                                                        

      stop
      end

