C********************************************************************
C
C     PROGRAM EKPROF
C     
C     水平平均熱収支鉛直分布図作成プログラム
C
C     20001/05/23  小高正嗣
C
C*********************************************************************
      program ekprof
      implicit real*8 ( a-h, o-z )
C---------------------------------------------------------------------  
#include "grid_size_M.f"
      parameter ( nzprof = 100 )
      parameter ( nzp = nz+1   )
C---------------------------------------------------------------------  
      parameter ( day = 86400.0 )
      parameter ( hour = 3600.0 )
C                                                                       
      character*5  nmrun
      character*4  from
      character*1  chara
      character*1  pz

      character*2  start
      character*2  finish
      character*2  date      
      character*20 header

      real*8  u ( -2:nxb, -2:nzb ), w ( -2:nxb, -2:nzb )
      real*8  kturb ( -2:nxb, -2:nzb ), tpot( -2:nxb, -2:nzb )
      real*4  tkturb( -2:nxb, -2:nzb ), tkbuoy( -2:nxb, -2:nzb )
      real*4  fdz(-2:nzb), fdzm(-2:nzb)

      real*4 zgrid(nzp), zgrid_p(nzp)
      real*4 ekvis(nzp), eknlv(nzp), ekbuoy(nzp)
      real*4 tkturz(nzp), tkbuoyz(nzp)

      real*4 RMAX1, RMIN1, RMAX2, RMIN2

      integer fnum, begin, end, interv
      integer nloop1
      character CSGI*1
      character*100 dumy

      real*4 ekvisall, eknlvall, ekbuoyall, tkturball, tkbyouall
      real*4 ekgen, eklostvsi, eklostnlv, eklost

*-----------------------------------------------------

      real*8 z0(-2:nzb), p0(-2:nzb), d0(-2:nzb), 
     \       temp0(-2:nzb), tpot0(-2:nzb), ppai0(-2:nzb)

*-----------------------------------------------------

      CHARACTER*36 CNUM
      CHARACTER*4  FNIN
      CHARACTER*2  SNAME
      CHARACTER*20 FNAME(100)

*-----------------------------------------------------

      RMAX1 = 5.0e-4
      RMIN1 = - RMAX1
      RMAX2 = 5.0e-3
      RMIN2 = - RMAX2

      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
      close(11)
      
      tinterv = nloop2 * ( dtime / 2. )

      call FOPEN ( irun+1 )


      CNUM = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
      SNAME = 'T1'
      FNIN  = '@' // SNAME // CNUM(IRUN+1  :IRUN+1  )
      FNAME(34) = FNIN // '.TPOT.DATA'
      FNAME(38) = FNIN // '.TURB.DATA'

      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

      call FCLOSE

      fdz(-2)= z0(0)
      fdz(-1)= z0(0)
      fdz(0) = z0(0)*2
      do iz = 1, nz-1
         fdz(iz) = ( z0(iz) - z0(iz-1) ) * 2.0 - fdz(iz-1)
      end do
      fdz(nz) = fdz(nz-1)
      fdz(nz+1) = fdz(nz-1)
      fdz(nz+2) = fdz(nz-1)

      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)
C      call sldiv ('Y', 3, 2 )
      call sldiv ('Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )


      OPEN (34, FILE = FNAME(34), STATUS = 'OLD', FORM = 'UNFORMATTED' )
      OPEN (38, FILE = FNAME(38), STATUS = 'OLD', FORM = 'UNFORMATTED' )
C      OPEN (10, FILE = 'ekall' ) 

      do it = 1, ite

      read(34) time, ((tpot(ix,iz),ix=-2,nxb), iz=-2,nzb )
      read(38) time, ((kturb(ix,iz),ix=-2,nxb), iz=-2,nzb)


      fnum = 50 + it
      open(fnum)

!********************************************************************
! * kinetic energy production & dissipation *

      read(fnum,'(a2,a15,g15.7,a6,g15.7,a4)') 
     &     chara, nmrun, t_begin, from, t_end,
     &     pz

      do iz = 1, nzp
         read(fnum,*) zgrid(iz), 
     &           ekvis(iz),
     &           eknlv(iz),
     &           ekbuoy(iz)
         end do

!********************************************************************
      if ( t_begin .eq. 0.0 ) then
         t_begin = t_end - tinterv
      end if

      begin = INT( t_begin -  INT( t_begin / day ) * day )
      end   = INT( t_end -  INT( t_begin / day ) * day )

      t_begin = t_begin + times
      t_end   = t_end   + times

      write(date, '(I2.2)') INT( t_begin / day + 1 )

      t_begin = ( t_begin - INT( t_begin / day ) * day  ) / hour 
      t_end   = ( t_end   - INT( t_end   / day ) * day  ) / hour 
      
      if ( INT(t_begin) .ge. 24 ) then
         t_begin = t_begin - 24
      end if

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

      write(start,'(I2.2)') INT(t_begin)
      write(finish,'(I2.2)') INT(t_end)

!********************************************************************

      header = 'day' // date // 
     \     ' LT=' // start // ':00-' // finish // ':00'

      close(fnum)

      if ( it .ge. itb   .and.  it .le. ite ) then
      if ( mod( (it-itb), interv ) .eq. 0 ) then


!**********************************************************************
!* turbulent keinetic energey dissipation ver turbulent diffusion

         do iz = -2, nzb
            do ix = -2, nxb
               tkturb(ix,iz) = d0(iz) * kturb(ix,iz)**3.0 /
     \              ( 0.2**2 * MIN( fdz(iz), z0(iz) )**4.0 ) 
     \              * fdz(iz) * 100.0
            end do
         end do


!* turbulent keinetic energey generation by byoyency force

         do iz = -2, nzb
            do ix = 02, nzb
               tkbuoy(ix,iz) =  
     \            + kturb(ix,iz) * d0(iz) * MAX ( 0.0, 
     \              ( ( tpot(ix,iz+1) + tpot0(iz+1) ) -
     \                ( tpot(ix,iz-1) + tpot0(iz-1) ) ) ) * grav /
     \              ( 2.0 * tpot0(iz) * fdz(iz) ) * fdz(iz) * 100.0
            end do
         end do

         do iz = 1, nzp
            tkturz(iz) = 0.0
            do ix = 1, nx
               tkturz(iz) = tkturz(iz) + tkturb(ix,iz-1) / REAL(nx)
            end do
         end do

         do iz = 1, nzp
            tkbuoyz(iz) = 0.0
            do ix = 1, nx
               tkbuoyz(iz) = tkbuoyz(iz) + tkbuoy(ix,iz-1) / REAL(nx)
            end do
         end do

!**********************************************************************
! time & volume integration

         ekvisall = 0.0
         eknlvall = 0.0
         ekbuoyall = 0.0
         tkturball = 0.0
         tkbuoyall = 0.0

         do iz = 1, nzp
            ekvisall  = ekvisall  + ekvis (iz) 
            eknlvall  = eknlvall  + eknlv (iz) 
            ekbuoyall = ekbuoyall + ekbuoy(iz) 
            tkturball = tkturball + tkturz(iz) 
            tkbuoyall = tkbuoyall + tkbuoyz(iz) 
         end do

         ekgen  = ekgen + ekbuoyall
         eklostvis = eklostvis + ekvisall
         eklostnlv = eklostnlv + eknlvall
         eklost = eklost + ekvisall + eknlvall
         tkturbsum = tkturbsum + tkturball
         tkbuoysum = tkbuoybsum + tkbuoyall

C         write(10,*) 'Time                  ',
C     \     ' LT=' // start // ':00-' // finish // ':00'
C         write(10,*) 'Thermal to Keinetic   ', 
C     \        ekbuoyall*tinterv*REAL(nx)
C         write(10,*) 'Viscos dissipation    ', 
C     \        ekvisall*tinterv*REAL(nx)
C         write(10,*) 'Non Linear dissipation', 
C     \        eknlvall*tinterv*REAL(nx)
C         write(10,*) 'Turbulent Energy Gen.  ', 
C     \        tkbuoyall*tinterv*REAL(nx)
C         write(10,*) 'Turbulent dissipation ', 
C     \        tkturball*tinterv*REAL(nx)
C         write(10,*) '**********************'

!**********************************************************************
! 
         do iz = 1, nzp
            ekvis  (iz) = ekvis (iz) / 100.0 / fdz(iz)
            eknlv  (iz) = eknlv (iz) / 100.0 / fdz(iz)
            ekbuoy (iz) = ekbuoy(iz) / 100.0 / fdz(iz)
            tkturz (iz) = tkturz(iz) / 100.0 / fdz(iz)
            tkbuoyz(iz) = tkbuoyz(iz)/ 100.0 / fdz(iz)
         end do

!**********************************************************************


      do iz = 1, nzp
         zgrid_p(iz) = zgrid(iz-1) / 1000.
      end do

      call grfrm
      call grswnd( RMIN1, RMAX1, 0., 20. )
      call grsvpt( 0.2, 0.8, 0.35, 0.85 )
      call grstrn(1)
      call grstrf
      call ussttl ( '', '', '', 'Z[km]' )
      call usdaxs
      call uxsttl ( 'T', header, 1. )
      call sglset( 'LCHAR', .TRUE. )

c      call sgsplc( 'gen' )
c      call uulinz(nzp, ekbuoy, zgrid_p, 1, 3)
c      call sgsplc( 'vis' )
c      call uulinz(nzp, ekvis , zgrid_p, 2, 3)
c      call sgsplc( 'nlv' )
c      call uulinz(nzp, eknlv , zgrid_p, 3, 3)
c      call sglset( 'LCHAR', .FALSE. )
c      call uulinz(nzp, tkturz , zgrid_p, 4, 3)

c      call sgsplc( 'buoy' )
c      call uulinz(nzp, tkbuoyz, zgrid_p, 1, 1)

      call sgsplc( 'gen' )
      call SGSPLI( 24 )
      call sgplu(nzp, ekbuoy, zgrid_p)
      call sgsplc( 'vis' )
      call SGSPLI( 34 )
      call sgplu(nzp, ekvis, zgrid_p)
      call sgsplc( 'nlv' )
      call SGSPLI( 44 )
      call sgplu(nzp, eknlv, zgrid_p)
      call sgsplc( 'turb' )
      call SGSPLI( 94 )
      call sgplu(nzp, tkturz , zgrid_p)
      call sglset( 'LCHAR', .FALSE. )

      call grfig
      call grswnd( RMIN2, RMAX2, 0., 1. )
      call grsvpt( 0.2, 0.8, 0.2, 0.3 )
      call grstrn(1)
      call grstrf
      call uxaxdv ( 'T', RMAX2*0.1, RMAX2*0.2 )
      call uxaxdv ( 'B', RMAX2*0.1, RMAX2*0.5 )
      call uyaxdv ( 'L', 0.2, 1. )
      call uyaxdv ( 'R', 0.2, 1. )
      call uzfact(0.8)
C      call uxsttl( 'B', '[m2/sec3]', 1.2 )
      call uxsttl( 'B', '[J/m^3/sec]', 1.2 )
      call uzfact(1.25)
      call uxsttl( 'B', 'dEk/dt', 0. )

      call sglset( 'LCHAR', .TRUE. )

c      call sgsplc( 'gen' )
c      call uulinz(15, ekbuoy, zgrid_p, 1, 3)
c      call sgsplc( 'vis' )
c      call uulinz(15, ekvis , zgrid_p, 2, 3)
c      call sgsplc( 'nlv' )
c      call uulinz(15, eknlv , zgrid_p, 3, 3)

c      call sglset( 'LCHAR', .FALSE. )
c      call uulinz(15, tkturz , zgrid_p, 4, 3)
c      call sgsplc( 'buoy' )
c      call uulinz(15, tkbuoyz, zgrid_p, 1, 1)

      call sgsplc( 'gen' )
      call SGSPLI( 24 )
      call sgplu(15, ekbuoy, zgrid_p)
      call SGSPLI( 34 )
      call sgsplc( 'vis' )
      call sgplu(15, ekvis, zgrid_p)
      call sgsplc( 'nlv' )
      call SGSPLI( 44 )
      call sgplu(15, eknlv, zgrid_p)
      call SGSPLI( 94 )
      call sgsplc( 'turb' )
      call sgplu(15, tkturz , zgrid_p)

      call sglset( 'LCHAR', .FALSE. )

      end if
      end if

      end do

      call grcls

C      write(10,*) 'Total Energy Generation', 
C     \        ekgen*tinterv*REAL(nx)
C      write(10,*) 'Total Energy Lost by VIS', 
C     \        eklostvis*tinterv*REAL(nx)
C      write(10,*) 'Total Energy Lost by NLV', 
C     \        eklostnlv*tinterv*REAL(nx)
C      write(10,*) 'Total Energy Lost       ', 
C     \        eklost*tinterv*REAL(nx)
C      write(10,*) 'Turbulent dissipation   ', 
C     \        (tkturbsum+tkbuoysum)*tinterv*REAL(nx)
C      write(10,*) 'Generation - Lost       ', 
C     \        (ekgen + eklost)*tinterv*REAL(nx)

C      close(10)
      close(34)
      close(38)

      end



