C*********************************************************************
C
C     PROGRAM hfluxprf
C
C     水平平均鉛直分布図作成プログラム(熱フラックス)
C
C     1997/07/24 中島健介
C     1999/04/15 小高正嗣
C     2000/11/21 小高正嗣; D 論用に書き直し
C
C*********************************************************************
      program hfluxprf
*-----------------------------------------------------
      implicit real*8 ( a-h, o-z )
*-----------------------------------------------------
#include "grid_size_M.f"
*-----------------------------------------------------
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )

      real*8 
     \     wread(-2:nx+2,-2:nz+2), tread(-2:nx+2,-2:nz+2)

      real*4 
     \     w(nx+1,nz+1), t(nx+1,nz+1), tt(nx+1,nz+1), 
     \     ave_w(nz+1) , ave_t(nz+1) , 
     \     wd(nx+1,nz+1), td(nx+1,nz+1), wdtd(nx+1,nz+1), 
     \     ave_wt(nz+1), ave_wdtd(nz+1), wt_mean(nz+1)

      character*2 mint
      character*2 lt
      character*2 date
      character*17 header
      character*100 dumy

*-----------------------------------------------------
      real*8 z0(-2:nzb), p0(-2:nzb), d0(-2:nzb), 
     &       temp0(-2:nzb), tpot0(-2:nzb), ppai0(-2:nzb)
      real*4 zgrid(nz+1)
*-----------------------------------------------------
      do iz = 1, nz+1
         wt_mean(iz) = 0.0
      end do

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

      open (11,file='exparam')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, times
      read (11,*) ylat, als

      call FOPEN ( irun+1 )

      do iz = 1,9
         read(16,'(a100)') dumy
      end do

      do iz = -2, nzb
         read(16,*) ii, z0(iz), p0(iz), d0(iz),
     &              temp0(iz), tpot0(iz), ppai0(iz), dum1, dum2
      end do

      do iz = 1, nz+1
         zgrid(iz) = z0(iz-1)/1000.
      end do

      write(6,*) 'data number = ', nloop1
      write(0,*) 'ibegin(2-', nloop1, 
     &           ') iend(2-', nloop1,  
     &           ') inter(1-', nloop1, ')'
      read (5,*) itb, ite, interv

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

      call grfrm
      call grswnd( -25., 25., 0., 20. )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn(1)
      call grstrf
      call ussttl ( 'Heat Flux ', 'W/m^2', '', 'Z[km]' )
      call usdaxs

      il = 1
      do it = 1, ite

      read(33) time, ((wread(ix,iz),ix=-2,nx+2), iz=-2,nz+2)
      read(34) time, ((tread(ix,iz),ix=-2,nx+2), iz=-2,nz+2)

      write(0,*) 'reading time = ', time

      if ( it .ge. itb   .and.  it .le. ite ) then

      write(0,*) 'writing time = ', time


c*** 時刻変換: 秒から LT へ************************************
      time = time + times
      write(date,'(I2.2)') INT( time / day ) + 1

      time = ( time - INT( time / day ) * day ) / hour 

      if ( INT(time) .ge. 24 ) then
         time = time - 24
      end if

      write(lt, '(I2.2)') INT(time)
      write(mint, '(I2.2)') INT( MOD(time,1)*60.0 )

      header = 'day=' // date 
c**************************************************************

      do ix = 1, nx+1
         do iz = 1,nz+1
            w(ix,iz) = wread(ix-1,iz-1)
            t(ix,iz) = tread(ix-1,iz-1)
            tt(ix,iz) = tread(ix-1,iz-1) + tpot0(iz-1)
         end do
      end do


      do iz = 1,nz+1
         ave_w(iz) = 0.0
         ave_t(iz) = 0.0
         do ix = 1, nx+1
            ave_w(iz) = ave_w(iz) + w(ix,iz) / (nx+1)
            ave_t(iz) = ave_t(iz) + tt(ix,iz) / (nx+1)
         end do
      end do

      do iz = 1,nz+1
         do ix = 1, nx+1
            wd(ix,iz) = w(ix,iz) - ave_w(iz)
            td(ix,iz) = tt(ix,iz) - ave_t(iz)
         end do
      end do

      do iz = 1,nz+1
         ave_wt(iz) = cpdry * d0(iz) * ave_w(iz)* ave_t(iz)
         do ix = 1, nx+1
            wdtd(ix,iz) = wd(ix,iz)*td(ix,iz)
         end do
      end do

      do iz = 1,nz+1
         ave_wdtd(iz) = 0.0
         do ix = 1, nx
            ave_wdtd(iz) = ave_wdtd(iz) + wdtd(ix,iz) / (nx+1)
         end do
         ave_wdtd(iz) = cpdry * d0(iz) * ave_wdtd(iz)
      end do

      do iz = 1,nz+1
            wt_mean(iz) = wt_mean(iz) + ave_wdtd(iz) / (ite-itb+1)
      end do

      if ( mod( (it-itb), interv ) .eq. 0 ) then
         CALL sgsplc( lt )
         call uulinz(nz+1, ave_wdtd, zgrid, il, 1)
         il = il + 1
      end if

      end if

      end do

      CALL sglset( 'LCHAR', .FALSE. )
      call uulinz(nz+1, wt_mean, zgrid, 1, 3)
      call uxsttl ( 'T', header, 1. )

      call grcls
      call FCLOSE ( irun+1 )


      end



