!=================================================================
!  鉛直 1 次元火星大気放射対流平衡モデル: メインプログラム
!  拡散係数パラメタリゼーション(Crank-Nikolson 法)
!
!  1999/10/17 小高正嗣
!  1999/11/06 小高正嗣
!  2000/02/14 小高正嗣
!  2000/02/22 小高正嗣; 大気上端からの下向き放射を考慮
!  2000/05/11 小高正嗣; バンドモデルへ
!  2000/05/19 小高正嗣; 高速化
!  2000/09/22 小高正嗣; 地表面温度を計算
!  2000/09/25 小高正嗣; 近赤外放射を計算
!  2000/10/19 小高正嗣; ダスト放射を計算
!  2000/11/01 小高正嗣; 透過率計算のバグを修正
!=================================================================
      program mars1D_CLEAR

#include "grid_size_M.f"
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      

      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.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), TPOT1(-2:NZB)

      dimension
     \   TEMPG(NZG)

      dimension
     \   DTPRAD(-2:NZB), DTSRAD(-2:NZB),  DTPDI0(-2:NZB),
     \   DTPQRAD(-2:NZB), DTQSRAD(-2:NZB)

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

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

      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 
     \   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
     \   FXQRADU (-2:NZB), FXQRADD (-2:NZB), FXQRADN (-2:NZB),
     \   FXQRADU1(-2:NZB), FXQRADD1(-2:NZB), FXQRADN1(-2:NZB), 
     \   FXQRADU2(-2:NZB), FXQRADD2(-2:NZB), FXQRADN2(-2:NZB),
     \   QRAD1   (-2:NZB), QRAD2   (-2:NZB)

      dimension
     \   SXQRADU(-2:NZB), SXQRADD(-2:NZB), SXQRADN(-2:NZB), 
     \   SXSOL  (-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 
     \   QTAU   (-2:NZB), QDTAU   (-2:NZB),
     \   QTAUIR1(-2:NZB), QDTAUIR1(-2:NZB),
     \   QTAUIR2(-2:NZB), QDTAUIR2(-2:NZB)

      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 
     \     ASR(0:NZ,0:NZ), BSR(0:NZ,0:NZ), 
     \     CSR(0:NZ,0:NZ), DSR(0:NZ,0:NZ), 
     \     RSR(0:NZ), SSR(0:NZ)

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

      dimension 
     \   TKENRG(-2:NZB),
     \   CDTURB(-2:NZB), 
     \   CD1   (-2:NZB), 
     \   CD2   (-2:NZB),
     \   STABZ (-2:NZB),
     \   SSTURB(-2:NZB),
     \   DTQDI0(-2:NZB)

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

      dimension
     \   FHEAT ( 2 )

      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)

      dimension
     \   GGRID(-1:NZGB), GDZ(0:NZGB), GDZM(0:NZG)

      logical SWITCH

      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                      )

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

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

      call SETDUST
     O   ( QEXT_SR, QEXT_IR1, QEXT_IR2, 
     O     QALB_SR, QALB_IR1, QALB_IR2, 
     O     QASM_SR, QASM_IR1, QASM_IR2,
     O     QTAU_ALL, QDENS  , 
     O     QREFF   , QRMOD              )

!-----------------------------------------------------------------      
! 鉛直格子の設定
!-----------------------------------------------------------------      

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

      DATA ( GGRID(IZ), IZ=-1,NZGB )
     1     / -6.0D0 ,
     2       -6.0D0 , -4.0D0 , -2.7D0 , -1.8D0, -1.2D0, 
     3       -0.79D0, -0.53D0, -0.35D0, -0.2D0, -0.1D0, 
     4        0.0D0 ,  0.0D0 /

      call SETGGD
     I   ( SKDPS, GGRID, 
     O     GDZ  , GDZM, 
     D     NZG  , NZGB       ) 

!-----------------------------------------------------------------      
! 実験番号の設定と入出力ファイル設定
!-----------------------------------------------------------------      

      READ(5,*) IRUN
      CALL FOPEN ( IRUN )

!-----------------------------------------------------------------      
! 時間積分パラメータの設定
!-----------------------------------------------------------------      
! TIME0   計算開始時刻(sec)
! NLOOP1  1出力毎の計算ステップ数(I)
! NLOOP2  出力回数(I)
! DTIME   時間刻(sec)
! TIMES   計算開始時刻の LT=0:00 からのずれ(sec)

      WRITE(6,*)'TIME0, NLOOP1, NLOOP2, DTIME, TIMES'
      READ (5,*) TIME0, NLOOP1, NLOOP2, DTIME, TIMES
      WRITE(6,*) TIME0, NLOOP1, NLOOP2, DTIME, TIMES

! YLAT   緯度(deg)
! ALS    火星中心黄径(deg)

      WRITE(6,*)'YLAT, ALS'
      READ (5,*) YLAT, ALS
      WRITE(6,*) YLAT, ALS

      YLAT = YLAT * PI/180.0D0  ! rad へ変換
      ALS  = ALS  * PI/180.0D0  ! rad へ変換
!-----------------------------------------------------------------      
! 初期条件の設定
!-----------------------------------------------------------------      

      IF ( TIME0  .EQ.  0.0D0 ) THEN

         do IZ = 0, NZ
            TEMP0( IZ ) = TEMPZ( ZGRID(IZ) )
         end do

         do IZ = 0, NZ
            CDTURB( IZ ) = 0.0D0
         end do

         do IZ = 1, NZG
            TEMPG( IZ ) = 210.0D0
         end do

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

         call BOUND ( TEMP0, NZ, NZB )
         call BOUND ( CDTURB, NZ, NZB )

         TIME = TIME0
      
      ELSE

         CALL FLREAD
     I     (  TIME0, 
     O        TIME ,
     O        TPOT0, TEMP0, CDTURB, DTPRAD, DTPDI0, 
     O        TEMPG, 
     D        NZ   , NZB  , NZG              )

      END IF

      IF ( IRUN .EQ. 1 ) then

         CALL FLWRIT
     I     (  TIME ,
     I        TPOT0, TEMP0, CDTURB, DTPRAD, DTPDI0, 
     I        TEMPG, DTSRAD, DTPQRAD, DTQSRAD,
     D        NZ   , NZB  , NZG                   )

      END IF

!-----------------------------------------------------------------      
! 放射場の初期条件設定; バンドパラメターの設定
! (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 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 QBASIC               ! ダスト分布の計算
     I   ( DENS0  , QTAU_ALL, 
     O     QVAP0  , BETA ,  
     C     QEXT_SR, 
     C     PI, QREFF, QDENS,
     G     FDZ, DX, 
     D     NZ , NZB )

      do IZ = -2, NZB
         WRITE(16,*) IZ, ZGRID(IZ), PRES0(IZ), DENS0(IZ),
     \                   TEMP0(IZ), TPOT0(IZ), PPAI0(IZ), 
     \                   QVAP0(IZ), RVAP0(IZ)
      end do

      call CLFXSOL
     I   ( S00  , ALS   , YLAT , TIME, TIMES,
     O     S0   , COSZ  , 
     C     EXCTL, OBLIC , OMGA0, DAY , PI    )

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

      call CLTAUN               ! CO2 の等価関数の計算
     I   ( OPL , SNU, RNU, PRES0,
     O     TAUN,
     C     DELNU, pi, 
     D     NZ  , NZB, NBND  )

      if ( COSZ .gt. 0.0D0 ) then

         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 )

      end if

      call SETQTAU              ! ダストの光学的厚さの計算
     I   ( DENS0, QVAP0, BETA,
     O     QTAU , QDTAU, 
     C     QALB_SR, QASM_SR, 
     G     FDZ  , DX,
     D     NZ   , NZB   )

! 赤外に対するダストの光学的厚さ

      do iz = -2, NZB
         QTAUIR1 (iz) = QTAU (iz) * 0.5
         QDTAUIR1(iz) = QDTAU(iz) * 0.5
      end do

      do iz = -2, NZB
         QTAUIR2 (iz) = QTAUIR1 (iz) * QEXT_IR2 / QEXT_IR1
         QDTAUIR2(iz) = QDTAUIR1(iz) * QEXT_IR2 / QEXT_IR1
      end do

! δ-Eddington 近似

      QALB_SRE = QALB_DE( QALB_SR, QASM_SR ) 
      QASM_SRE = QASM_DE( QASM_SR )

      QALB_IR1 = QALB_DE( QALB_IR1, QASM_IR1 )
      QASM_IR1 = QASM_DE( QASM_IR1 )

      QALB_IR2 = QALB_DE( QALB_IR2, QASM_IR2 )
      QASM_IR2 = QASM_DE( QASM_IR2 )

!-----------------------------------------------------------------      
! ループ計算開始
!-----------------------------------------------------------------      

      do INL1 = 1, NLOOP1

! 1 ステップ前の値を保存
         do IZ = -2, NZB
            TPOT1(iz)=TPOT0(iz)
         end do

         do INL2 = 1, NLOOP2
!-----------------------------------------------------------------      

         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 QBASIC            ! ダスト分布の計算
     I      ( DENS0  , QTAU_ALL, 
     O        QVAP0  , BETA ,  
     C        QEXT_SR, 
     C        PI, QREFF, QDENS,
     G        FDZ, DX, 
     D        NZ , NZB )

!-----------------------------------------------------------------      
! 放射計算
! 放射計算は 1 分に 1 回計算する(2000/11/01).

         if ( TIME .eq. TIME0 .or. MOD(TIME, 60.0) .eq. 0.0D0 ) then

!-----------------------------------------------------------------      
! 赤外放射計算
! 光路長, 透過関数は HOUR ごとに 1 回計算
         if ( MOD(TIME, HOUR) == 0.0D0 ) then
            call SETOPL         ! CO2 光路長の計算
     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 SETQTAU        ! ダストの光学的厚さの計算
     I         ( DENS0, QVAP0, BETA,
     O           QTAU , QDTAU, 
     C           QALB_SR, QASM_SR, 
     G           FDZ  , DX,
     D           NZ   , NZB   )

            do iz = -2, NZB
               QTAUIR1 (iz) = QTAU (iz) * 0.5
               QDTAUIR1(iz) = QDTAU(iz) * 0.5
            end do

            do iz = -2, NZB
               QTAUIR2 (iz) = QTAUIR1 (iz) * QEXT_IR2 / QEXT_IR1
               QDTAUIR2(iz) = QDTAUIR1(iz) * QEXT_IR2 / QEXT_IR1
            end do

         end if

!-----------------------------------------------------------------      
! CO2 赤外放射加熱率の計算

         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  )

! ダストによる赤外放射加熱率の計算; (赤外 5-11.6 μm)

         DINU1  = 900.0D0
         DINU2  = 2000.0D0
         DELNUD = DELNU * 4.0 

         call CLDIRAD
     I      ( QTAUIR1, QDTAUIR1, PRES0 , DENS0, TEMP0, TPOT0, 
     I        DINU1  , DINU2   , DELNUD,
     O        FXQRADU1, FXQRADD1, FXQRADN1, QRAD1 , 
     C        QALB_IR1, QASM_IR1, 
     C        CPDRY , TSFC  , GRAV  , PI   , 
     W        ASR   , BSR   , CSR   , DSR  ,
     W        RSR   , SSR   ,
     G        FDZ   , FDZM  , DX    ,
     D        NX    , NXB   , NZ    , NZB          )


! ダストによる赤外放射加熱率の計算; (赤外 20-200 μm)

         DINU1  = 50.0D0
         DINU2  = 500.0D0
         DELNUD = DELNU * 2.0 

         call CLDIRAD
     I      ( QTAUIR2, QDTAUIR2 , PRES0 , DENS0, TEMP0, TPOT0, 
     I        DINU1  , DINU2    , DELNUD,
     O        FXQRADU2, FXQRADD2, FXQRADN2, QRAD2 , 
     C        QALB_IR2, QASM_IR2, 
     C        CPDRY , TSFC  , GRAV  , PI   , 
     W        ASR   , BSR   , CSR   , DSR  ,
     W        RSR   , SSR   ,
     G        FDZ   , FDZM  , DX    ,
     D        NX    , NXB   , NZ    , NZB          )

         do iz = -2, NZB
            DTPQRAD(iz) = QRAD1(iz) + QRAD2(iz)
            FXQRADU(iz) = FXQRADU1(iz) + FXQRADU2(iz)
            FXQRADD(iz) = FXQRADD1(iz) + FXQRADD2(iz)
            FXQRADN(iz) = FXQRADN1(iz) + FXQRADN2(iz)
         end do

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

            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 )
         end if 
!-----------------------------------------------------------------
! 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

! ダスト

         call CLDSRAD
     I      ( QTAU  , QDTAU  , PRES0  , DENS0, TPOT0, 
     I        S0    , COSZ   , 
     O        SXQRADU, SXQRADD, SXQRADN, DTQSRAD,  SXSOL, 
     C        ALBEDO, QALB_SRE, QASM_SRE,
     C        DELNU , CPDRY  , TSOL   , SIGMA, 
     C        GRAV  , PI     , 
     W        ASR   , BSR    , CSR    , DSR  ,
     W        RSR   , SSR    ,
     G        FDZ   , FDZM  , DX    ,
     D        NZ    , NZB                      )

!-----------------------------------------------------------------      
         else 
         do iz = -2, NZB
            SXRADU(iz) = 0.0D0
            SXRADD(iz) = 0.0D0
            SXRADN(iz) = 0.0D0
            SXQRADU(iz)= 0.0D0
            SXQRADD(iz)= 0.0D0
            SXQRADN(iz)= 0.0D0
            SXSOL(iz)  = 0.0D0
            SQRAD1(iz) = 0.0D0
            SQRAD2(iz) = 0.0D0 
            SQRAD3(iz) = 0.0D0
            DTSRAD(iz) = 0.0D0
            DTQSRAD(iz)= 0.0D0
         end do
!-----------------------------------------------------------------      
         end if
!-----------------------------------------------------------------      
         end if
!-----------------------------------------------------------------      
! 地表風速の計算(2次元計算結果をフィッテイング)

         TIMEL = MOD( TIME, 86400.0 ) + 21600.0D0
         VSFC0 = 0.5D0 + 10.0D0 * 
     \        exp ( - ( ( TIMEL - 54000.0D0 ) / 10800.0D0 )**2 )

         call CLDRAG            ! 抵抗係数の計算
     I      ( TPOT0, 
     O        RICHD, CDRAGV,
     G        ZGRID, 
     C        GRAV , TSFC , VSFC0,
     D        NX   , NXB  , NZB   )

         call CDCOEF            ! 半整数格子上の拡散係数の計算(含リミッター)
     I      ( CDTURB,
     O        CD1   , CD2   ,
     B        ZGRID ,
     D        NZ    , NZB     )

         call FXSURF            ! 地表熱フラックスの計算
     I      ( TPOT0, CDRAGV, VSFC0, 
     O        FHEAT, 
     B        DENS0, PPAI0,
     C        CPDRY, GRAV ,
     S        TSFC , 
     D        NX   , NXB  , NZ   , NZB )

!-----------------------------------------------------------------      
! 1 ステップ時間積分
!-----------------------------------------------------------------      

         do IZ = 0, NZ
            TPOT0(IZ) = TPOT0(IZ) + DTIME * 
     \           (  DTPRAD(IZ)  + DTSRAD(IZ)  
     \          +   DTPQRAD(iz) + DTQSRAD(iz) ) 
     \           * ( TPOT0( iz ) / TEMP0( iz ) )
         end do

! 地表からの熱フラックスを加算

         TPOT0(0) = TPOT0(0)
     \     + DTIME * FHEAT(1) / DX / FDZ(0) / DENS0(0) / CPDRY
!         TPOT0(NZ-1) = TPOT0(NZ-1)
!     \     + DTIME * FHEAT(2) / DX / FDZ(NZ-1) / DENS0(NZ-1) / CPDRY

!-----------------------------------------------------------------      
! 拡散項の計算

         SWITCH = .FALSE.
         do iz = NZB, -2, -1
            if ( CD1(iz) .ne. 0.0D0 ) then
               SWITCH = .TRUE.
            end if
         end do

         if ( SWITCH ) then
!            write(6,*) TIME
!            write(6,*) 'call CLDIFEQ'
            call CLDIFEQ
     I         ( TPOT0, CD2,
     B           DENS0, 
     C           DTIME,
     G           DX, FDZ, FDZM,
     D           NZ   , NZB     )
!            write(6,*) 'finished CLDIFEQ'
         end if

!         TPOT0(   -1 ) = TPOT0(    0 )
!         TPOT0(   -2 ) = TPOT0(    1 )
         call BOUND ( TPOT0, NZ, NZB )

! 温位から温度を計算

         do IZ = 0, NZ
            TEMP0( IZ ) = TPOT0( IZ ) * PPAI0( IZ )
         end do

         TEMP0(   -1 ) = TEMP0(    0 )
         TEMP0(   -2 ) = TEMP0(    1 )
!         call BOUND( TPOT0, NZ, NZB )
!         call BOUND( TEMP0, NZ, NZB )

!-----------------------------------------------------------------      
! 地表面温度の計算
!-----------------------------------------------------------------      

         F0 = - SXQRADN(0) - ( SXRADD(NZ+1) - SXRADD(0) )
!         F0 = S0 * COSZ - ( SXRADD(NZ+1) - SXRADD(0) )

         do iz = -2, NZB
            FXRADD(iz) = FXRADD(iz) + FXQRADD(iz)
         end do

         call CLTSFC
     I   ( FHEAT , FXRADD, F0,
     O     TEMPG , TSFC  ,
     C     SIGMA , ALBEDO, CPGRND, 
     C     CDTURBG,DTIME , DENSG ,
     G     GDZ   , GDZM  , 
     D     NX    , NXB   , NZ   , NZB   ,
     D     NZG   , NZGB                   )
         
!-----------------------------------------------------------------      
! 乱流拡散係数の計算
!-----------------------------------------------------------------      

         call CLKDIF 
     I   ( TPOT0 , TSFC , 
     O     CDTURB,
     C     GRAV  ,
     G     ZGRID , FDZ  , FDZM , DX,
     W     STABZ ,
     D     NX    , NXB  , NZ   , NZB )

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

         TIME = TIME + DTIME

         end do
!-----------------------------------------------------------------      
! 出力
!-----------------------------------------------------------------      

! 2-grid ノイズの除去
         do iz = 1, NZ-1
            CD2(iz) = ( CD2(iz+1) + CD2(iz)*2.0D0 + CD2(iz-1) ) *0.25D0
         end do

         TPOT_DIFF = 0.0
         do IZ = 0, NZ
            TPOT_DIFF = MAX( TPOT_DIFF, ABS(TPOT0(iz) - TPOT1(iz)) )
         end do

         write(*,*) TIME, 'TPOT_DIFF=',TPOT_DIFF

c         do IZ = -2, NZB
c            DTPRAD(IZ) = DTPRAD(IZ) + DTSRAD(IZ) 
c     \           + DTPQRAD(IZ) + DTQSRAD(IZ) 
c         end do

         CALL FLWRIT
     I     (  TIME ,
     I        TPOT0, TEMP0, CDTURB, DTPRAD, DTPDI0, 
     I        TEMPG, DTSRAD, DTPQRAD, DTQSRAD,
     D        NZ    , NZB , NZG )

         write(16,*) 'TIME=',TIME
         write(72) TIME, TSFC(1), FHEAT(1)

      end do

      CALL FCLOSE

      stop
      end

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