!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(放射フラックス)
!
!  2000/02/29  小高正嗣
!=================================================================
      program rflux
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = NZ )
      parameter ( NBND   = 16       )
      parameter ( DELNU  = 25.0D0   )
      parameter ( SIGMA  = 5.67D-8  )
      parameter ( PI = 3.14159263D0 )
      parameter ( DAY    = 86400.0D0)
      parameter ( HOUR   = 3600.0D0 ) 
!-----------------------------------------------------------------      
      dimension
     \   PRES0(-2:NZB), TEMP0(-2:NZB), TPOT0(-2:NZB),
     \   PPAI0(-2:NZB), DENS0(-2:NZB), 
     \   RVAP0(-2:NZB), QVAP0(-2:NZB), ZGRID(-2:NZB)

      dimension
     \     TEMPG(NZG)

      dimension
     \   FDZ(-2:NZB), FDZM(-2:NZB), OPL(-2:NZB)

      dimension
     \   BBRAD ( -2:NZB ), FXRADU( -2:NZB ), FXRADD( -2:NZB ),
     \   FXRADN( -2:NZB )

      dimension 
     \   FXRADU_NU(-2:NZB,NBND), FXRADD_NU(-2:NZB,NBND)

      dimension 
     \   TAUN (0:NZ+1,0:NZ+1,NBND)

      dimension
     \   WNU(NBND), SNU(NBND), RNU(NBND)

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

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

      real*4
     \   ZMIN , ZMAX

      character*2 timel
      character*5 header

      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                      )

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

      DATA ( WNU(I), I=1,NBND )
     1     /  512.5D0,  537.5D0,  562.5D0,  587.5D0,
     2        612.5D0,  637.5D0,  662.5D0,  687.5D0,  
     3        712.5D0,  737.5D0,  762.5D0,  787.5D0,
     4        812.5D0,  837.5D0,  862.5D0,  887.5D0  / 

      DATA ( SNU(I), I=1,NBND )
     1     /  1.952D-2,  2.785D-1,  5.495D-1,  5.331D0, 
     2        5.196D+2,  7.778D+3,  8.746D+4,  2.600D+4,
     3        1.232D+3,  2.042D+2,  7.278D0 ,  1.337D0,
     4        3.974D-1,  1.280D-2,  2.501D-3,  3.937D-3  /

      DATA ( RNU(I), I=1,NBND )
     1     /  2.870D-1,  1.215D0 ,  2.404D0 ,  1.958D+1, 
     2        5.804D+1,  2.084D+2,  7.594D+2,  2.635D+2,
     3        8.387D+1,  2.852D+1,  6.239D0 ,  2.765D0 ,
     4        8.897D-1,  3.198D-1,  1.506D-1,  1.446D-1  /

!-----------------------------------------------------------------      
      open(10,file='exparam')
      read(10,*) IRUN
      READ(10,*) TIME0, NLOOP1, NLOOP2, DTIME, TIMES
      read(10,*) YLAT, ALS

      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

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

      ZMIN    = ZGRID_p(0)
      ZMAX    = ZGRID_p(NZP)
      DX = 100.0D0
!-----------------------------------------------------------------      

      FDZ(0) = ZGRID(0) * 2.0D0 / DX

      do IZ = 1, NZ-1
         FDZ(IZ) = ( ZGRID(IZ) - 
     \               ( ZGRID(IZ-1) + FDZ(IZ-1)*0.5D0*DX ) ) * 2.0D0 / DX
      end do

      FDZ(-1)   = FDZ(0)
      FDZ(-2)   = FDZ(1)
      FDZ(NZ)   = FDZ(NZ-1)
      FDZ(NZ+1) = FDZ(NZ-2)
      FDZ(NZ+2) = FDZ(NZ-3)

      do IZ = 0, NZ
         FDZM(IZ) = ( FDZ(IZ-1) + FDZ(IZ) ) / 2.0D0
      end do
C
      FDZM(-2)   = FDZM(2)
      FDZM(-1)   = FDZM(1)
      FDZM(NZ+1) = FDZM(NZ-1)
      FDZM(NZ+2) = FDZM(NZ-2)


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

      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.0, 300.0, ZMIN, ZMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'IR flux', '[W/m2]', '', 'Z[km]' )
      CALL usdaxs

      do it = 1, ite

         read(32) TIME, ( TEMP0(IZ), IZ=-2,NZB )
         read(36) TIME, ( TEMPG(IZ), IZ= 1,NZG )

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

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

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

            TSFC(0) = TEMPG(NZG)
              
            call VPROF          ! 大気構造の計算
     I         ( ZGRID, PSFC,
     O           PRES0, DENS0, TEMP0, TPOT0, PPAI0,
     O           QVAP0, RVAP0,
     C           GRAV , RDRY , CPDRY, PSD,  RMVAP,
     G           FDZ  , FDZM , DX   ,
     D           NZ   , NZB                         )

            call SETOPL         ! 光路長の計算
     I         ( PRES0, PSFC  ,
     O           OPL  , 
     C           GRAV , 
     D           NZ   , NZB     )

            call CLTAUN         ! 透過関数の計算(15 μm バンド)
     I         ( OPL , SNU, RNU, PRES0,
     O           TAUN,
     C           DELNU, pi, 
     D           NZ  , NZB, NBND  )
            
            do inb = 1, NBND

               call CLBBRAD 
     I            ( TEMP0, WNU(inb),
     O              BBRAD, BBRAD_S ,
     C              TSFC , pi      ,
     D              NZ   , NZB     , NXB )     

               call CLFXRDU
     I            ( BBRAD , BBRAD_S, TAUN, inb,
     O              FXRADU, 
     D              NZ    , NZB    , NBND  )
               
               call CLFXRDD
     I            ( BBRAD , TAUN, inb,
     O              FXRADD, 
     D              NZ    , NZB , NBND  )

               do iz = -2, NZB
                  FXRADU_NU(iz, inb) = FXRADU(iz)
                  FXRADD_NU(iz, inb) = FXRADD(iz)
               end do

               do iz = -2, NZB
                  FXRADU(iz) = 0.0D0
                  FXRADD(iz) = 0.0D0
               end do

            end do
            
            do iz = -2, NZB
               do inb = 1, NBND/2
                  FXRADU(iz) = FXRADU(iz) + DELNU * 
     \            ( FXRADU_NU(iz, inb) + FXRADU_NU(iz, NBND+1-inb) )
                  FXRADD(iz) = FXRADD(iz) + DELNU * 
     \            ( FXRADD_NU(iz, inb) + FXRADD_NU(iz, NBND+1-inb) )
               end do
            end do

            do iz = -2, NZB
               FXRADU(iz) = sigma * TSFC(0)**4 + FXRADU(iz)
            end do

            call CLFXRDN
     I         ( FXRADU, FXRADD,
     O           FXRADN,
     D           NZ    , NZB     )


c*** 時刻変換: 秒から LT へ************************************
            time = MOD( time, DAY ) + TIMES
            time = time / HOUR

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

            write(timel,'(I2.2)') INT(time)
            header = 'LT=' // timel
            
            do iz = 0, NZP
               QVAL1_p(iz) = real( FXRADU(iz) ) 
            end do

            do iz = 35, 40
               write(*,*) iz, TEMP0(iz), PRES0(iz), 
     \              FXRADU(iz), FXRADD(iz)
            end do
            write(*,*) NZ+1, TEMP0(NZ+1), PRES0(NZ+1), 
     \           FXRADU(NZ+1), FXRADD(NZ+1)
               
            CALL uulinz( NZP+1, QVAL1_p, ZGRID_p, 1, 2  )

            do iz = 0, NZP
               QVAL1_p(iz) = real( FXRADD(iz) ) 
            end do

            CALL uulinz( NZP+1, QVAL1_p, ZGRID_p, 2, 2  )

            CALL uxsttl( 'T', header, 1. )      


            end if
         end if

      end do

      CALL grcls                                                        


      stop
      end

