!=================================================================
!  鉛直 1 次元火星大気放射対流平衡モデル: メインプログラム
!  拡散係数パラメタリゼーション(Crank-Nikolson 法)
!
!  1999/10/17 小高正嗣
!  1999/11/06 小高正嗣
!  2000/02/14 小高正嗣
!  2000/02/22 小高正嗣; 大気上端からの下向き放射を考慮
!  2000/04/02 小高正嗣; 拡散係数を KW78 で評価, 
!=================================================================
      program mars1D_CN

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

      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
     \   DTPRAD(-2:NZB),
     \   DTPDI0(-2:NZB)

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

      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)

      dimension EQWDTH(0:NZ+1,0:NZ+1)
      dimension work(-2:NZB)

      logical SWITCH

      parameter ( SIGMA = 5.67D-8 )

      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 SETGRD
     O   ( ZGRID, FDZ, FDZM, DX,
     D     NZ, NZB                )

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

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

!-----------------------------------------------------------------      
! 時間積分パラメータの設定
!-----------------------------------------------------------------      

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

!-----------------------------------------------------------------      
! 初期条件の設定
!-----------------------------------------------------------------      

      IF ( TIME0  .EQ.  0.0D0 ) THEN

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

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

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

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

         TIME = TIME0

      ELSE

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

      END IF

      IF ( IRUN .EQ. 1 ) then

         CALL FLWRIT
     I      (  TIME ,
     I         TPOT0, TEMP0, CDTURB, DTPRAD, DTPDI0, 
     D         NZ   , NZB                   )

      END IF

      do IZ = -2, NZB
         WRITE(16,*) IZ, ZGRID(IZ), TEMP0(IZ), CDTURB(IZ)
      end do

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

      do INL1 = 1, NLOOP1
         do INL2 = 1, NLOOP2

         TIME = TIME + DTIME

!         do IX = -2, NXB
            TSFC (IX) = 270.0D0
!            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 CLRAD             ! 放射加熱率の計算
     I      ( TPOT0 , EQWDTH, PSFC,
     O        DTPRAD, FXRADU, FXRADD, FXRADN, BBRAD,
     B        TEMP0 , PRES0 , DENS0, FDZ  ,
     C        TSFC  , CPDRY , RDRY , DX   , 
     C        SIGMA ,
     D        NX    , NXB   , NZ   , NZB     )

!! 高度 NZ-5 層以上は放射強制固定

         do IZ = NZ-4, NZ
            DTPRAD( IZ ) = DTPRAD(NZ-5)
         end do

! 地表風速の計算(2次元計算結果をフィッテイング)

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

         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     )

c         do IX = -2, NXB
c            CDRAGV( IX ) = CD2( 0 )
c         end do

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

! 1 ステップ時間積分

         do IZ = 0, NZ-1
            TPOT0(IZ) = TPOT0(IZ) +
     \           DTIME  *  DTPRAD(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
            call CLDIFEQ
     I         ( TPOT0, CD1,
     B           DENS0, 
     C           DTIME,
     G           DX, FDZ, FDZM,
     D           NZ   , NZB     )
         end if

         call BOUND ( TPOT0, NZ, NZB )

! 温位から温度を計算

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

         call BOUND( TEMP0, NZ, NZB )

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

c         call CLTSFC
c     I   ( FHEAT , FXRADD, FXRADS,
c     O     TSFC  ,
c     C     SIGMA ,
c     D     NX    , NXB   , 
c     D     NZ    , NZB            )

!-----------------------------------------------------------------      
! 乱流拡散係数の計算
! 乱流エネルギーの格子点が拡散係数のそれとは異ることに注意
!-----------------------------------------------------------------      

         call ADDIF0            ! 拡散加熱率の計算
     I      ( TKENRG,  CD2 ,
     O        DTQDI0,
     B        DENS0,
     G        DX, FDZ, FDZM,
     D        NZ, NZB )

         call CLTURB            ! 乱流生成と散逸の計算
     O      ( SSTURB,
     I        TKENRG, CDTURB,
     B        TPOT0 , 
     G        ZGRID , FDZ   , FDZM  , DX    ,
     C        GRAV  , 
     W        STABZ , 
     D        NX    , NZ    , NXB   , NZB    )

         do IZ = 0, NZ-1
            TKENRG( IZ ) = TKENRG( IZ ) 
     \           + DTI * ( DTQDI0(IZ) + SSTURB(IZ) )
         end do

         do IZ = 0, NZ-1

            MIXLNG = MIN( ZGRID(IZ), FDZ(IZ) )

            CDTURB( IZ ) = 
     \           MAX( 0.0D0, 0.2D0 * sqrt( TKENRG(IZ) ) *  MIXLNG )

         end do

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

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

         CALL FLWRIT
     I     (  TIME ,
     I        TPOT0, TEMP0, CDTURB, DTPRAD, DTPDI0, 
     D        NZ   , NZB              )

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

         call dclanim4
     I     (  ZGRID , TIME,
!     I        FXRADU, FXRADD,  FXRADN, TEMP0,
     I        TPOT0 , CDTURB,  TEMP0, DTPRAD, 
     D        NZ    , NZB          )

      end do

      call grcls

      CALL FCLOSE

      stop
      end

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

