!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(水平平均温位)
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!=================================================================
      program tprof
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NZP = 50 )

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

      dimension
     \   EQWDTH(0:NZ+1,0:NZ+1),
     \   FDZ(-2:NZB), FDZM(-2:NZB),OPL(-2:NZB)

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

      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*14 header
      character*2 timel
      character*2 day

      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 = -2, NZB
         READ(16,*) ITW, ZGRID(IZ), TEMP0(IZ), CDTURB(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)

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

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

      do IZ = 1, NZ-1
         FDZ(IZ) = ( ZGRID(IZ) - 
     \               ( ZGRID(IZ-1) + FDZ(IZ-1)*0.5D0 ) ) * 2.0D0
      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)

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

      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( 100.0, 240.0, ZMIN, ZMAX )   
      CALL grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      CALL grstrn( 1 )
      CALL grstrf
      CALL ussttl( 'netflux', '[W/m2]', '', 'Z[km]' )
      CALL usdaxs

      il = 1

      do it = 1, ite

         read(31) TIME, ( TEMP0(IZ), 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 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 CLEQWD            ! 等価幅の計算
     I      ( TEMP0 , PRES0, OPL , 
     I        ZGRID , FDZ  , FDZM, DX, 
     O        EQWDTH,                       
     D        NZ    , NZB      )
         
         call CLBBRAD 
     I      ( TEMP0,
     O        BBRAD, 
     D        NZ   , NZB  )

         call CLFXRDU
     I      ( BBRAD , EQWDTH, TEMP0,
     O        FXRADU,
     C        SIGMA , TSFC  ,
     D        NZ    , NZB   ,  NXB   )      


         call CLFXRDD           ! 放射場の計算; 下向き放射束密度
     I      ( BBRAD , EQWDTH, TEMP0,
     O        FXRADD,
     D        NZ    , NZB           )


         call CLFXRDN           ! 放射場の計算; 正味の放射束密度
     I      ( FXRADU, FXRADD,
     O        FXRADN,
     D        NZ    , NZB     )
         
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 = 0, NZP
                  QVAL1_p(iz) = real( FXRADN(iz) )
               end do
               
!               CALL uuslni( 2 )
               CALL uulinz( NZP+1, QVAL1_p, ZGRID_p, il, 2  )
!               CALL uxsttl( 'T', header, 1. )      

               il = il + 1

            end if
         end if

      end do

      CALL grcls                                                        


      stop
      end

