!=================================================================
!  正味放射フラックスの描画
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!  2000/03/20  小高正嗣
!=================================================================
      program rad
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( SIGMA = 5.67D-8 )

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

      dimension
     \  TPOT ( -2:NXB, -2:NZB )

      dimension
     \  TEMP ( -2:NXB, -2:NZRB ) 

      DIMENSION
     \   FDZ  ( -2:NZB ), FDZM ( -2:NZB )

      dimension 
     \   OPL   (-2:NZRB),
     \   RGRID (-2:NZRB),
     \   RDZ   (-2:NZRB), 
     \   RDZM  (-2:NZRB)

      dimension 
     \   EQWDTH(0:NZR+1,0:NZR+1)

      dimension
     \   BBRAD ( -2:NXB, -2:NZRB ),
     \   FXRADU( -2:NXB, -2:NZRB ),
     \   FXRADD( -2:NXB, -2:NZRB ),
     \   FXRADN( -2:NXB, -2:NZRB ),
     \   DTPRAD( -2:NXB, -2:NZRB )

      dimension
     \  FXRADN_mean(-2:NZRB),
     \  DTPRAD_mean(-2:NZRB),
     \  TEMP_mean(-2:NZRB)

      dimension TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB )

      real*4
     \   ZGRID_p(0:NZB),
     \   QVAL1_p(0:NZB) 

      real*4
     \   ZMIN , ZMAX

      character*14 header
      character*2 timel
      character*2 day

      character*100 dumy

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

!-----------------------------------------------------------------      
! 物理定数の設定
!-----------------------------------------------------------------      

      call SETCST
     O   ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT )

      call SETCS2
     O   ( CDRAG, VSFC0, TSFC, QVSFC,
     O     FCOLI,
     O     PSFC,
     C     RMVAP,
     C     FXSOL,
     D     NX, NXB                      )

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

      open(10,file='exparam')
      read(10,*) IRUN
      read(10,*) TIME0, NLOOP1, NLOOP2, DTIME

      write(day,'(I2.2)') irun

      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

      do iz = 0, nzb
         ZGRID_p(iz)= real( ZGRID(iz) )/ 1000.0
      end do

      ZMIN    = ZGRID_p(0)
      ZMAX    = ZGRID_p(NZB)

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

      call SETGRD
     O   ( ZGRID, FDZ, FDZM, DX,
     O     RGRID, RDZ, RDZM,
     D     NZ, NZB, NZR, NZRB             )

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

      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( -10.0, 10.0, ZMIN, ZMAX )   
c      CALL grswnd( 100.0, 350.0, ZMIN, ZMAX )   
c      CALL grswnd( 250.0, 350.0, ZMIN, ZMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'heating rate', '[K/day]', '', 'Z[km]' )
c      CALL ussttl( 'netflux', '[W/m2]', '', 'Z[km]' )
      CALL usdaxs

      il = 1

      do it = 1, ite

         read(34) TIME, ((TPOT(IX,IZ),IX=-2,nxb), IZ=-2,NZB )
         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
              
         do IX = -2, NXB
            TSFC (IX) = 180.0D0 + MAX( 0.0D0, 
     \           100.0D0 * sin( 2.0D0 * 3.14D0 * TIME / 86400.0D0 ) )
         end do

         call RBASIC
     I   ( RGRID, PSFC  , RDZ  , RDZM  ,
     O     OPL  , EQWDTH,
     C     GRAV , RDRY  ,
     D     NZ   , NZB   , NZR  , NZRB  , DX  )

!-----------------------------------------------------------------      
! 温度場の計算

      do iz = 0, NZ
         do ix = -2, NXB
            TEMP( ix, iz ) = ( TPOT( ix, iz ) + TPOT0( iz ) )  
     \           * ( PRES0( iz ) / PSFC ) ** ( RDRY / CPDRY )
         end do
      end do

      do iz = NZ+1, NZR
         do ix = -2, NXB
          TEMP( ix, iz ) = TEMP( ix, NZ )
         end do
      end do

      do ix = -2, NXB
         TEMP( ix, -1 )    = TEMP( ix, 0 )
         TEMP( ix, -2 )    = TEMP( ix, 0 )
         TEMP( ix, NZR+1 ) = TEMP( ix, NZ )
         TEMP( ix, NZR+2 ) = TEMP( ix, NZ )
      end do

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

         call CLBBRAD
     I   ( TEMP , 
     O     BBRAD,
     D     NX   , NXB, 
     D     NZR  , NZRB  )

         call CLFXRDU
     I   ( BBRAD , EQWDTH, TEMP,
     O     FXRADU,
     C     SIGMA , TSFC  ,
     D     NZR   , NZRB  , NX    , NXB    )      

         call CLFXRDD
     I   ( BBRAD , EQWDTH, TEMP,
     O     FXRADD,
     D     NZR   , NZRB  , NX    , NXB    )      

         call CLFXRDN
     I   ( FXRADU, FXRADD,
     O     FXRADN,
     D     NX    , NXB   , 
     D     NZR   , NZRB     )

         call CLQRAD
     I   ( FXRADN, DENS0, TEMP0 , TPOT0,
     I     FDZ   , 
     O     DTPRAD,
     C     CPDRY , DX   ,
     D     NX    , NXB  , NZ    , NZB  ,
     D     NZR   , NZRB                  )

         
c*** 時刻変換: 秒から LT へ************************************
               time = time - 86400.* ( irun -1 )
               time = time / 3600. + 6

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

               write(timel,'(I2.2)') INT(time)
               header = 'day' // day // ' LT=' // timel // ':00'

               do iz = -2, NZB
                  FXRADN_mean(iz) = 0.0D0
                  DTPRAD_mean(iz) = 0.0D0
                  TEMP_mean(iz) = 0.0D0
               end do

               do iz = 0, NZB
               do ix = 0, NX-1
                  FXRADN_mean(iz) = FXRADN_mean(iz) + FXRADN(ix,iz)/NX
                  DTPRAD_mean(iz) = DTPRAD_mean(iz) + DTPRAD(ix,iz)/NX
               end do
               end do

               do iz = 0, NZB
c                  QVAL1_p(iz) = real( FXRADN_mean(iz) )
                  QVAL1_p(iz) = real( DTPRAD_mean(iz) ) * 86400.0
               end do

!               CALL uuslni( 2 )
               CALL uulinz( NZB+1, QVAL1_p, ZGRID_p, 1, 2  )
!               CALL uxsttl( 'T', header, 1. )      

!               il = il + 1

            end if
         end if

      end do

      CALL grcls                                                        

      stop
      end

