!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  DCL による描画(成分ごとの鉛直加熱率分布)
!
!  2000/10/18  小高正嗣
!=================================================================
      program qrate
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NBND   = 16       )
      parameter ( NBNDSR1= 10 ) ! 近赤外バンド数1(4.3μm)
      parameter ( NBNDSR2= 10 ) ! 近赤外バンド数2(2.7μm)
      parameter ( NBNDSR3=  8 ) ! 近赤外バンド数3(2.0μm)
      parameter ( DELNU  = 25.0D0   )
      parameter ( TSOL   = 5760.0D0 )! 太陽表面温度
      parameter ( SIGMA  = 5.67D-8  )
      parameter ( PI  = 3.141592D0 )
      parameter ( NT  = 100 )
      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), TPOT1(-2:NZB)

      dimension
     \     TEMPG (NZG)         ! 地中温度

      dimension
     \   DTPRAD(-2:NZB),
     \   DTSRAD(-2:NZB),
     \   DTPDI0(-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),
     \   SXRADU_NU1(-2:NZB,NBNDSR1), SXRADD_NU1(-2:NZB,NBNDSR1),
     \   SXRADU_NU3(-2:NZB,NBNDSR3), SXRADD_NU3(-2:NZB,NBNDSR3)

      dimension
     \   SXRADU (-2:NZB), SXRADD (-2:NZB), SXRADN (-2:NZB), 
     \   SXRADU1(-2:NZB), SXRADD1(-2:NZB), SXRADN1(-2:NZB), 
     \   SXRADU2(-2:NZB), SXRADD2(-2:NZB), SXRADN2(-2:NZB), 
     \   SXRADU3(-2:NZB), SXRADD3(-2:NZB), SXRADN3(-2:NZB), 
     \   SQRAD1 (-2:NZB), SQRAD2 (-2:NZB), SQRAD3 (-2:NZB)

      dimension 
     \   TAUN    (0:NZ+1,0:NZ+1,NBND),
     \   TAUN_SR1(0:NZ+1,0:NZ+1,NBNDSR1),
     \   TAUN_SR2(0:NZ+1,0:NZ+1,NBNDSR2),
     \   TAUN_SR3(0:NZ+1,0:NZ+1,NBNDSR3)

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

      dimension
     \   WNU_SR1(NBNDSR1), SNU_SR1(NBNDSR1), RNU_SR1(NBNDSR1),
     \   WNU_SR2(NBNDSR2), SNU_SR2(NBNDSR2), RNU_SR2(NBNDSR2),
     \   WNU_SR3(NBNDSR3), SNU_SR3(NBNDSR3), RNU_SR3(NBNDSR3)

      dimension 
     \   CDTURB(-2:NZB)

      dimension
     \   TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB ),
     \   CDRAGV( -2:NXB ),
     \   RICHD( -2:NXB )

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

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

      real*4
     \   QRAD_p(0:NZ), SQRAD_p(0:NZ), QCONV(0:NZ), ZGRID_p(0:NZ),
     \   QMIN, QMAX, ZMIN, ZMAX

      character*14 dheader
      character*18 pheader
      character*100 model

      character*2 date
      character*2 timel
      character*5 CLAT
      character*5 CLS
      character*2 CSGI
      
!-----------------------------------------------------------------      

      call SETCST( grav, psd, amdry, rdry, cpdry, rmvap, alatent )

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

      write(CLS,  '(F5.1)') ALS
      write(CLAT, '(F5.1)') YLAT

      pheader = 'Ls=' // CLS // ' ' // CSGI(172) // '=' // CLAT

      YLAT = YLAT * PI/180.0D0  ! rad へ変換
      ALS  = ALS  * PI/180.0D0  ! rad へ変換

      QMIN = - 100.0
      QMAX =   100.0

      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                      )

      call SETCSG 
     I   ( DAY   ,
     O     ALBEDO, DIFG, CPGRND, DENSG, CDTURBG, SKDPS )

      call SETCSP
     I   ( PI ,
     C     S00, EXCTL, OBLIC, OMGA0 )

      call SETGRD
     O   ( ZGRID, FDZ, FDZM, DX,
     D     NZ, NZB                )

!-----------------------------------------------------------------      
! 放射場の初期条件設定; バンドパラメターの設定
! (Houghton, 1986;The physics of atmosphere, 2nd ed., 
!  Cambridge Univ. Press)

      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  /

! CO2 4.3 μm

      DATA ( WNU_SR1(I), I=1,NBNDSR1 )
     1     /  2212.5D0, 2237.5D0, 2262.5D0, 2287.5D0,
     2        2312.5D0, 2337.5D0, 2362.5D0, 2387.5D0,  
     3        2412.5D0, 2437.5D0                         /

      DATA ( SNU_SR1(I), I=1,NBNDSR1 )
     1     /  9.504D-1,  2.217D+2,  4.566D+3,  7.965D+3, 
     2        1.055D+5,  5.587D+5,  6.819D+5,  1.256D+4,
     3        7.065D-2,  8.522D-2                        /

      DATA ( RNU_SR1(I), I=1,NBNDSR1 )
     1     /  2.866D0 ,  3.000D+1,  1.134D+2,  2.011D+2, 
     2        5.880D+2,  1.206D+3,  1.182D+3,  8.873D+1,
     3        3.404D-1,  4.236D-1                        /

! CO2 2.7 μm

      DATA ( WNU_SR2(I), I=1,NBNDSR2 )
     1     /  3150.0D0, 3250.0D0, 3350.0D0, 3450.0D0,
     2        3550.0D0, 3650.0D0, 3750.0D0, 3850.0D0,  
     3        3950.0D0, 4050.0D0                         /

      DATA ( SNU_SR2(I), I=1,NBNDSR2 )
     1     /  1.324D-1,  7.731D-2,  1.232D0 ,  5.159D0, 
     2        4.299D+3,  1.543D+4,  1.649D+4,  1.180D-1,
     3        1.464D-2,  1.251D-2                        /

      DATA ( RNU_SR2(I), I=1,NBNDSR2 )
     1     /  9.836D-1,  4.900D-1,  2.952D0 ,  7.639D0 , 
     2        1.914D+2,  3.245D+2,  2.722D+2,  9.535D-1,
     3        2.601D-1,  2.021D-1                        /


! CO2 2.0 μm

      DATA ( WNU_SR3(I), I=1,NBNDSR3 )
     1     /  4650.0D0, 4750.0D0, 4850.0D0, 4950.0D0,
     2        5050.0D0, 5150.0D0, 5250.0D0, 5350.0D0     /

      DATA ( SNU_SR3(I), I=1,NBNDSR3 )
     1     /  2.185D-1,  2.040D0 ,  1.197D+2,  4.829D+2, 
     2        8.778D+1,  8.346D+1,  8.518D-2,  4.951D-1  /

      DATA ( RNU_SR3(I), I=1,NBNDSR3 )
     1     /  1.916D0,   6.475D0 ,  3.112D+1,  5.759D+1, 
     2        2.012D+1,  1.804D+1,  8.474D-1,  1.597D0   / 

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

      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

      write(6,*) 'data number = ', nloop1
      write(6,*) 'INPUT DATA BEGINING, END, & INTERVAL' 
      write(6,*) 'itb ite interval ?'
      read (*,*) itb, ite, interv

      do IZ = 0, NZ
         ZGRID_p(iz) = REAL( ZGRID(iz) ) / 1000.0
      end do
      
      ZMIN = 0.0
      ZMAX = ZGRID_p(NZ)

      CALL sgpwsn
      read(*,*) iws
      CALL gropn(iws)
      CALL sglset( 'LCORNER', .FALSE. )
      CALL sglset( 'LCHAR', .TRUE. )
      CALL sldiv ( 'Y', 2, 3 )
      CALL sglset( 'LFULL', .TRUE. )

      do it = 1, ite

         read(31) TIME, ( TPOT0(IZ), IZ=-2,NZB )
         read(32) TIME, ( TEMP0(IZ), IZ=-2,NZB )
         read(33) TIME, ( CDTURB(IZ), IZ=-2,NZB )
         read(36) TIME, ( TEMPG(IZ), IZ=1,NZG )
         read(72) TIME, TSFC(1), HEAT

         do ix = -2, NXB
            TSFC(ix) = TEMPG(NZG) 
         end do

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

         write(0,*) 'writing time = ', time
         write(date, '(I2.2)') INT( time/DAY ) + 1

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)
         dheader = 'day' // date // ' LT=' // timel // ':00'
!=======================================================================
         call CLFXSOL
     I      ( S00  , ALS   , YLAT , TIME, TIMES,
     O        S0   , COSZ  , 
     C        EXCTL, OBLIC , OMGA0, DAY , PI    )

         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  )
         
         call CLRAD             
     I      ( TPOT0 , TAUN  , PSFC  , OPL   ,
     I        WNU   , SNU   , RNU   ,
     O        DTPRAD, FXRADU, FXRADD, FXRADN, BBRAD,
     B        TEMP0 , PRES0 , DENS0 , FDZ   ,
     C        TSFC  , CPDRY , RDRY  , DX    , SIGMA,
     C        DELNU , GRAV  , PI    ,
     W        FXRADU_NU, FXRADD_NU,
     D        NX    , NXB   , NZ    , NZB   , NBND  )

!-----------------------------------------------------------------      
! 近赤外放射加熱率の計算
! 天頂角パラメータ COSZ が 0<COS<1 のときのみ計算
! COSZ=0 の場合, 加熱率は 0 
!-----------------------------------------------------------------      
         if ( COSZ .gt. 0.0D0 ) then
!-----------------------------------------------------------------      
! 光路長, 透過関数は HOUR ごとに 1 回計算

            do IZ = -2, NZB     ! 光路長を天頂角方向に射影
               OPL_S(IZ) = OPL(IZ) / COSZ
            end do

            call CLTAUN
     I         ( OPL_S, SNU_SR1, RNU_SR1, PRES0,
     O           TAUN_SR1, 
     C           DELNU, PI, 
     D           NZ   , NZB, NBNDSR1 )

            call CLTAUN
     I         ( OPL_S, SNU_SR2, RNU_SR2, PRES0,
     O           TAUN_SR2, 
     C           DELNU*4.0, PI, 
     D           NZ   , NZB, NBNDSR2 )

            call CLTAUN
     I         ( OPL_S, SNU_SR3, RNU_SR3, PRES0,
     O           TAUN_SR3, 
     C           DELNU*4.0, PI, 
     D           NZ   , NZB, NBNDSR3 )

!-----------------------------------------------------------------
! CO2 4.3 μm

            call CLSRAD
     I         ( TPOT0,  TAUN_SR1  , 
     I           WNU_SR1, SNU_SR1, RNU_SR1, 
     I           S0    , COSZ  ,
     O           SXRADU1, SXRADD1, SXRADN1, SQRAD1,
     B           TEMP0 , PRES0 , DENS0 , 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU , GRAV  , PI    ,
     W           SXRADU_NU1, SXRADD_NU1,
     G           FDZ   , 
     D           NZ    , NZB   , NBNDSR1 )

! CO2 2.7 μm

            call CLSRAD
     I         ( TPOT0,  TAUN_SR2  , 
     I           WNU_SR2, SNU_SR2, RNU_SR2, 
     I           S0    , COSZ  ,
     O           SXRADU2, SXRADD2, SXRADN2, SQRAD2,
     B           TEMP0 , PRES0 , DENS0 , 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU*4.0 , GRAV  , PI    ,
     W           SXRADU_NU1, SXRADD_NU1,
     G           FDZ   , 
     D           NZ    , NZB   , NBNDSR2 )

! CO2 2.0 μm

            call CLSRAD
     I         ( TPOT0,  TAUN_SR3, 
     I           WNU_SR3, SNU_SR3, RNU_SR3, 
     I           S0    , COSZ  ,
     O           SXRADU3, SXRADD3, SXRADN3, SQRAD3,
     B           TEMP0 , PRES0 , DENS0 , 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU*4.0 , GRAV  , PI    ,
     W           SXRADU_NU3, SXRADD_NU3,
     G           FDZ   , 
     D           NZ    , NZB   , NBNDSR3 )

            do iz = -2, NZB
               DTSRAD(iz) = SQRAD1(iz) + SQRAD2(iz) + SQRAD3(iz)
               SXRADD(iz) = SXRADD1(iz) + SXRADD2(iz) + SXRADD3(iz)
            end do
!-----------------------------------------------------------------      
         else 
            do iz = -2, NZB
               SXRADU(iz) = 0.0D0
               SXRADD(iz) = 0.0D0
               SXRADN(iz) = 0.0D0
               SQRAD1(iz) = 0.0D0
               SQRAD2(iz) = 0.0D0 
               SQRAD3(iz) = 0.0D0
               DTSRAD(iz) = 0.0D0
            end do
!-----------------------------------------------------------------      
         end if
!-----------------------------------------------------------------      

         do IZ = 0, NZ
            QRAD_p (iz) = REAL( DTPRAD(iz)*( TPOT0( iz )/TEMP0( iz ) ))
     \           * DAY
            SQRAD_p(iz) = REAL( DTSRAD(iz)*( TPOT0( iz )/TEMP0( iz ) ))
     \           * DAY

            QRAD_p (iz) = QRAD_p (iz) + SQRAD_p(iz)

            if ( QRAD_p (iz) > QMAX ) then
               QRAD_p (iz) = QMAX
            else if ( QRAD_p (iz) < QMIN) then
               QRAD_p (iz) = QMIN
            end if

         end do

         QCONV(0) = 1.0 / ( DENS0(0) * FDZ(0) * DX ) * (
     \           ( DENS0 (1) + DENS0 (0) ) * 0.5 *
     \           ( CDTURB(1) + CDTURB(0) ) * 0.5 *
     \           ( TPOT0 (1) - TPOT0 (0) ) / ( FDZM(1) * DX )  
     \            - HEAT / CPDRY ) * DAY 

         do IZ = 1, NZ
            QCONV(iz) = 1.0 / ( DENS0(iz) * FDZ(iz) * DX ) * (
     \           ( DENS0 (iz+1) + DENS0 (iz) ) * 0.5 *
     \           ( CDTURB(iz+1) + CDTURB(iz) ) * 0.5 *
     \           ( TPOT0 (iz+1) - TPOT0 (iz) ) / (FDZM(iz+1)+DX) -
     \           ( DENS0 (iz) + DENS0 (iz-1) ) * 0.5 *
     \           ( CDTURB(iz) + CDTURB(iz-1) ) * 0.5 *
     \           ( TPOT0 (iz) - TPOT0 (iz-1) ) / (FDZM(iz)+DX) )* DAY
         end do
         
         do iz = NZ, 0, -1
            if ( CDTURB(iz) > 0.0D0 ) then
               QCONV(iz+1) = QCONV(iz+2)
               QCONV(iz)   = QCONV(iz-1)
               exit
            end if
         end do

         do iz = 0, NZ
            if ( QCONV(iz) > QMAX ) then
               QCONV (iz) = QMAX
            else if ( QCONV(iz) < QMIN) then
               QCONV (iz) = QMIN
            end if
         end do

         CALL grfrm
         CALL grswnd( QMIN, QMAX, 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]' )
         CALL usdaxs

         call uulinz(NZ, QRAD_p , ZGRID_p, 1, 3  )
!         call uulinz(NZ, SQRAD_p, ZGRID_p, 2, 3  )
         call uulinz(NZ, QCONV  , ZGRID_p, 3, 3  )

         CALL uxsttl( 'T', dheader, 1. )      
         CALL uxsttl( 'T', pheader, 1. )      
         CALL uxsttl( 'T', model, 1. )      

         end if
         end if

      end do

      CALL grcls                                                        

      call FCLOSE

      stop
      end

