!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(水平平均温位)
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!=================================================================
      program plank

!-----------------------------------------------------------------      
      implicit real*4 ( a-h,o-z )
!-----------------------------------------------------------------      

!-----------------------------------------------------------------      
      parameter ( NZ    = 100     )
      parameter ( SIGMA = 5.67e-8 )
      parameter ( NU    = 667.0   )

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

      TSFC_MAX = 280.0
      TSFC_MIN = 180.0
      DTSFC = ( TSFC_MAX - TSFC_MIN )/NZ 

      TIME_MAX = 86400.0
      TIME_MIN = 0.0
      DTIME = ( TIME_MAX - TIME_MIN )/NZ 

      TIME = TIME_MIN  
      BSFC_MEAN = 0.0

      do ix = 1, NZ+1
         TSFC = 180.0D0 + 1.0D+2 * 
     \              MAX( sin( 2.0d0 * PI * TIME / 8.64D+4 ), 0.0D0 )
         BSFC = 1.91e-8 * nu**3 / 
     \        ( exp( 1.438 * nu / TSFC ) - 1.0 )

         BSFC_MEAN = BSFC_MEAN + BSFC / (NZ+1)

         TIME = TIME + DTIME
      end do


      BT = 0.5 * BSFC_MEAN

      TEMP = 1.438 * nu / 
     \     log( 1.0 + ( 1.91e-8 * NU**3 ) / BT )

      WRITE(*,*) TEMP

      stop

!      TEMP_MIN = TEMP(1)
!      TEMP_MAX = TEMP(NZ+1)

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

      CALL grfrm
      CALL grswnd( TSFC_MIN, TSFC_MAX, TEMP_MIN, TEMP_MAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'Tsfc', '[K]', 'Tco2', '[K]' )
      CALL usdaxs

      CALL uuslni( 2 )
!      CALL uulin( NZ+1, TSFC, TEMP )
      
      CALL grcls                                                        


      end

