!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(地表からの熱フラックス)
!
!  2000/02/29  小高正嗣
!=================================================================
      program hflux
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter( NZP = 50 )
      parameter( PI  = 3.141592D0 )
      parameter ( NT  = 100 )
      parameter ( DAY = 86400.0D0 )
      parameter ( HOUR= 3600.0D0 )
!-----------------------------------------------------------------      

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

      dimension
     \     TEMPG (NZG)         ! 地中温度

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


      real*4
     \   TGRD(0:NT), TIMEL(0:NT), FLUX(0:NT),
     \   FMIN, FMAX

      character*18 pheader
      character*100 model

      character*5 CLAT
      character*5 CLS
      character*2 CSGI
      
!-----------------------------------------------------------------      

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

      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

      FMIN = -5.0
      FMAX = 20.0

      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

c      write(6,*) 'data number = ', nloop1
c      write(6,*) 'INPUT DATA BEGINING, END, & INTERVAL' 
c      write(6,*) 'itb ite interval ?'
c      read (*,*) itb, ite, interv


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

         TGRD (it+ilate) = TEMPG(NZG) 
         TIMEL(it+ilate) = ( MOD( TIME, DAY ) + TIMES ) / HOUR

         do ix = -2, NXB
            TSFC(ix) = TEMPG(NZG) 
         end do

         VSFC0 = 0.5D0 + 10.0D0 * 
     \        exp ( - ( ( TIMEL(it+ilate) * HOUR - 54000.0D0 ) 
     \                   / 10800.0D0 )**2 )

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

c         FLUX(it+ilate) = - CDRAGV(0) * VSFC0 * CPDRY * DENS0(0)
c     \        * ( TPOT0(0)*PPAI0(0) - TEMPG(NZG) )

         FLUX(it+ilate) = real( HEAT )
         
c         write(*,*) CDRAGV(0), VSFC0, CPDRY, DENS0(0), 
c     \        TEMP0(0) - TEMPG(NZG), FLUX(it+ilate)
         
      end do

      TIMEL(ilate) = TIMEL(nloop1+ilate)
      TGRD(ilate)  = TGRD(nloop1+ilate)
      FLUX(ilate)  = FLUX(nloop1+ilate)

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

      CALL sgpwsn
      read(*,*) iws
      CALL gropn(iws)
      CALL sglset( 'LCORNER', .FALSE. )

      CALL grfrm
c      CALL grswnd( 0., 24.0, 160., 300. )   
      CALL grswnd( 0., 24.0, FMIN, FMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )  
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'Local Time', '[Hour]', 'Heat flux', '[W/m^2]' )
      CALL usdaxs
      call uulinz(nloop1+1, TIMEL, FLUX, 1, 3  )
c      call uulinz(nloop1+1, TIMEL, TGRD, 1, 3  )


      CALL uxsttl( 'T', pheader, 1. )      
      CALL uxsttl( 'T', model, 1. )      

      CALL grcls                                                        

      call FCLOSE

      stop
      end

