!= deepconv/arare 湿潤大気対流計算用主プログラム 
!
!= deepconv/arare main program for moist atmospheric convection 
!
! Authors::   杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu
! Version::   $Id: arare.f90,v 1.52 2015/02/25 04:25:17 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
! ソースコード書き換え済み (村橋 究理基)

program deepconv_arare
  !
  ! 非静力学モデル deepconv/arare 湿潤大気対流計算用主プログラム (三次元版)
  !

  ! モジュール引用  use statement 
  !

  ! gtool5 関連 
  ! gtool5 modules
  !
  use dc_types,      only: STRING, DP
  use dc_message,    only: MessageNotify
  ! dc_string を変数読み込みに利用するため追加 (mkuriki)
  use dc_string, only: toChar
  ! use gtool_history, only: HistoryPut
  ! gtool の nc データ読み込み用の関数を追加 (mkuriki)
  use gtool_history, only: HistoryPut, HistoryGet
  use gtool_historyauto, only: HistoryAutoPut

  ! 初期設定モジュール
  ! Initialize module
  !
  use mpi_wrapper, only : myrank, MPIWrapperInit, MPIWrapperFinalize
  use argset,   only: argset_init
  use clockset, only: ClocksetInit, ClocksetClose, ClocksetPredict, &
    &                 ClockSetPreStart, ClockSetLoopStart, ClockSetPreStop, ClockSetLoopStop
  use gridset,  only: gridset_init, imin, imax, jmin, jmax, kmin, kmax, nx, ny, nz, ncmax
  use timeset,  only: timeset_init, TimesetDelTimeHalf, TimesetProgress, &
    &                 TimeA, TimeN, NstepShort, NstepOutput, &
    &                 EndTime, FlagInitialRun
  use axesset,  only: axesset_init, xyz_avr_pyz, xyz_avr_xqz, xyz_avr_xyr
  use constants,only: constants_init
  use composition, only: composition_init, SpcWetSymbol
  use fileset,  only: fileset_init
  use basicset, only: xyzf_QMixBZ, xyz_DensBZ, xyz_EffMolWtBZ, &
    &                 xyz_PTempBZ, xyz_TempBZ, xyz_PressBZ, &
    &                 xyz_VelSoundBZ, xyz_ExnerBZ, xyz_VPTempBZ
  use ChemCalc, only: ChemCalc_init, ChemCalc_init2
!  use ChemCalc, only: ChemCalc_init
  use namelist_util, only: NmlutilInit, namelist_filename

  ! 下請けモジュール
  ! Utility modules
  !
  use cflcheck, only : CFLCheckTimeShort, &
    &                  CFLCheckTimeLongVelX, &
    &                  CFLCheckTimeLongVelY, &
    &                  CFLCheckTimeLongVelZ
  use damping, only: Damping_init, SpongeLayer_forcing
  use energymonit, only: EnergyMonit_init, EnergyMonit_exec
  use setmargin, only: SetMargin_init

  ! 力学過程計算用関数モジュール
  ! Dynamical processes module
  !
  use DynamicsHEVI, only: Dynamics_Init,           &
    &                     Dynamics_Long_forcing,   &
    &                     Dynamics_Short_forcing,  &
    &                     Dynamics_Long_integrate
  use fillnegative,only: FillNegative_init

  ! 乱流拡散計算用モジュール
  ! Turbulent diffusion module
  !
  use Turbulence_kw1978,  only: Turbulence_kw1978_Init, Turbulence_KW1978_forcing
  use Turbulence_constKm, only: Turbulence_constKm_Init, Turbulence_constKm_forcing

  ! 放射強制計算用モジュール
  ! Radiative forceing module
  !
  use Radiation_simple,  only: Radiation_simple_init, xyz_PTempRadConst, xyz_PTempRadVary, &
    &                          xyz_ExnerRadConst, xyz_ExnerRadVary
  use Radiation_sounding,  only: Radiation_sounding_init, xyz_PTempRadSndg,  &
    &                           xyz_ExnerRadSndg
  use Radiation_bs1998,   only: Radiation_bs1998_init, xyz_PTempRadBS1998,  &
    &                           xyz_ExnerRadBS1998
  use radiation_heatbalance, only: Radiation_heatbalance_init, Radiation_heatbalance_forcing

  ! 境界からのフラックス計算用モジュール
  ! Surface flux module
  !
  use Surfaceflux_Diff, only: Surfaceflux_Diff_init, Surfaceflux_Diff_forcing
!  use Surfaceflux_Bulk_l1982, only: Surfaceflux_Bulk_init, Surfaceflux_Bulk_forcing
  ! --------------------------------------
  ! surfaceflux_Bulk を surfaceflux_Bulk_o2001 を利用するように変更 (mkuriki)
  ! use Surfaceflux_Bulk, only: Surfaceflux_Bulk_init, Surfaceflux_Bulk_forcing
  use Surfaceflux_Bulk_o2001, only: Surfaceflux_Bulk_init, Surfaceflux_Bulk_forcing
  ! use Surfaceflux_Bulk_l1979, only: Surfaceflux_Bulk_init, Surfaceflux_Bulk_forcing
  use Surfaceflux_Const, only: Surfaceflux_Const_init, Surfaceflux_Const_forcing
  use Surfaceflux_bs1998, only: Surfaceflux_bs1998_init, Surfaceflux_bs1998_forcing

  ! 湿潤過程計算用モジュール
  ! Moist processes modules
  !
  use Cloudphys_K1969, only: Cloudphys_K1969_Init, Cloudphys_K1969_forcing !, &
!    &                        Cloudphys_K1969_FallRain
  use Cloudphys_marscond, only: cloudphys_marscond_Init, cloudphys_marscond_forcing


  ! ファイル入出力モジュール
  ! File I/O module
  !
  use RestartFileIO, only : ReStartFileio_init, ReStartFileio_Finalize, &
    &                       ReStartFileio_BZ_Get, ReStartFileio_Var_Get, rstat
  use HistoryFileIO, only: HistoryFileio_init, HistoryFileio_Finalize


  ! 暗黙の型宣言禁止
  !
  implicit none

  ! 内部変数
  ! Internal variables
  !
  character(*), parameter:: prog_name = 'arare'
                            ! 主プログラム名.
                            ! Main program name

  real(DP), allocatable :: pyz_VelXBl(:,:,:)    
                             ! $ u (t-\Delta t) $ 東西風 ; zonal wind
  real(DP), allocatable :: pyz_VelXNl(:,:,:)    
                             ! $ u (t) $          東西風 ; zonal wind
  real(DP), allocatable :: xyz_VelXNl(:,:,:)    
                             ! $ u (t) $          東西風 ; zonal wind
  real(DP), allocatable :: pyz_VelXAl(:,:,:)    
                             ! $ u (t+\Delta t) $ 東西風 ; zonal wind
  real(DP), allocatable :: pyz_VelXNs(:,:,:)    
                             ! $ u (\tau) $ 東西風 ; zonal wind
  real(DP), allocatable :: pyz_VelXAs(:,:,:)    
                             ! $ u (\tau +\Delta \tau) $ 東西風 ; zonal wind
  real(DP), allocatable :: xqz_VelYBl(:,:,:)    
                             ! $ v (t-\Delta t) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xqz_VelYNl(:,:,:)    
                             ! $ v (t) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xyz_VelYNl(:,:,:)    
                             ! $ v (t) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xqz_VelYAl(:,:,:)    
                             ! $ v (t+\Delta t) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xqz_VelYNs(:,:,:)   
                             ! $ v (\tau -\tau) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xqz_VelYAs(:,:,:)
                             ! $ v (\tau) $ 南北風 ; meridonal wind
  real(DP), allocatable :: xyr_VelZBl(:,:,:)    
                             ! $ w (t-\Delta t) $ 鉛直風 ; vertical wind
  real(DP), allocatable :: xyr_VelZNl(:,:,:)    
                             ! $ w (t) $ 鉛直風 ; vertical wind
  real(DP), allocatable :: xyz_VelZNl(:,:,:)    
                             ! $ w (t) $ 鉛直風 ; vertical wind
  real(DP), allocatable :: xyr_VelZAl(:,:,:)    
                             ! $ w (t+\Delta t) $ 鉛直風 ; vertical wind
  real(DP), allocatable :: xyr_VelZNs(:,:,:)    
                             ! $ w (\tau) $ 鉛直風 ; vertical wind
  real(DP), allocatable :: xyr_VelZAs(:,:,:) 
                             ! $ w (\tau +\Delta \tau)  鉛直風 ; vertical wind
  real(DP), allocatable :: xyz_ExnerBl(:,:,:)   
                             ! $ \pi (t-\Delta t) $ 圧力関数 ; Exner function
  real(DP), allocatable :: xyz_ExnerNl(:,:,:)   
                             ! $ \pi (t) $ 圧力関数 ; Exner function
  real(DP), allocatable :: xyz_ExnerAl(:,:,:)
                             ! $ \pi (t+\Delta t) $ 圧力関数 ; Exner function
  real(DP), allocatable :: xyz_ExnerNs(:,:,:)   
                             ! $ \pi (\tau -\Delta \tau) $ 圧力関数 ; Exner function
  real(DP), allocatable :: xyz_ExnerAs(:,:,:)   
                             ! $ \pi (\tau) $ 圧力関数 ; Exner function
  real(DP), allocatable :: xyz_PTempBl(:,:,:) 
                             ! $ \theta (t-\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PTempNl(:,:,:) 
                             ! $ \theta (t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PTempAl(:,:,:) 
                             ! $ \theta (t+\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PTempNs(:,:,:) 
                             ! $ \theta (t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PTempAs(:,:,:) 
                             ! $ \theta (t+\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_CDensBl(:,:,:) 
                             ! $ \theta (t-\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_CDensNl(:,:,:) 
                             ! $ \theta (t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_CDensAl(:,:,:) 
                             ! $ \theta (t+\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_CDensNs(:,:,:) 
                             ! $ \theta (t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_CDensAs(:,:,:) 
                             ! $ \theta (t+\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_KmBl(:,:,:)
                             ! $ Km (t-\Delta t) $ 乱流拡散係数 
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyz_KmNl(:,:,:)
                             ! $ K_m (t) $ 乱流拡散係数 
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyz_KmAl(:,:,:)
                             ! $ K_m (t+\Delta t) $ 乱流拡散係数 
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyz_KhBl(:,:,:)      
                             ! $ K_h (t-\Delta t) $ 乱流拡散係数
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyz_KhNl(:,:,:)
                             ! $ K_h (t) $ 乱流拡散係数 
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyz_KhAl(:,:,:)
                             ! $ K_h (t+\Delta t) $ 乱流拡散係数
                             ! Turbulent diffusion coeff. 
  real(DP), allocatable :: xyzf_QMixBl(:,:,:,:) 
                             ! $ q (t-\Delta t) $ 湿潤量の混合比
                             ! Mixing ratio of moist variables.
  real(DP), allocatable :: xyzf_QMixNl(:,:,:,:) 
                             ! $ q (t) $ 湿潤量の混合比 
                             ! Mixing ratio of moist variables
  real(DP), allocatable :: xyzf_QMixAl(:,:,:,:) ! 
                             ! $ q (t+\Delta t) $ 湿潤量の混合比 
                             !Mixing ratio of moist variables

  real(DP), allocatable :: pyz_DVelXDtNl(:,:,:)
  real(DP), allocatable :: xqz_DVelYDtNl(:,:,:)
  real(DP), allocatable :: xyr_DVelZDtNl(:,:,:)
  real(DP), allocatable :: xyz_DPTempDtNl(:,:,:)
  real(DP), allocatable :: xyz_DExnerDtNl(:,:,:)
  real(DP), allocatable :: xyz_DExnerDtNs(:,:,:)
  real(DP), allocatable :: xyzf_DQMixDtNl(:,:,:,:)
  real(DP), allocatable :: xyz_DKmDtNl(:,:,:)
  real(DP), allocatable :: xyz_DCDensDtNl(:,:,:)
  
  integer :: s, t, tau  ! do ループ変数 ; do loop variable 

  integer           :: IDTurbMethod           = 0
  integer, parameter:: IDTurbMethodKW1978     = 2
  integer, parameter:: IDTurbMethodConstKm    = 3
  integer           :: IDRadMethod            = 0
  integer, parameter:: IDRadMethodHeatConst   = 1
  integer, parameter:: IDRadMethodHeatVary    = 2
  integer, parameter:: IDRadMethodHeatBalance = 3
  integer, parameter:: IDRadMethodSounding    = 4
  integer, parameter:: IDRadMethodBaker1998   = 5
  integer           :: IDSurfaceMethod        = 0
  integer, parameter:: IDSurfaceMethodDiff    = 1
  integer, parameter:: IDSurfaceMethodBulk    = 2
  integer, parameter:: IDSurfaceMethodConst   = 3
  integer, parameter:: IDSurfaceMethodBaker1998 = 4
  integer           :: IDCloudMethod          = 0
  integer, parameter:: IDCloudMethodK1969     = 1
  integer, parameter:: IDCloudMethodMarsCond  = 2
  integer           :: IDDebugMethod          = 0
  integer, parameter:: IDDebugMethodNoTendencyLong= 1
  integer, parameter:: IDDebugMethodWindConst     = 2

  !------------------------------------------
  ! 自分で追加した変数 (mkuriki)
  ! HistoryGet 用の時間指定 文字型変数
  character(STRING) :: gtioTime
  ! HistoryGet 用の読み込み配列 配列長は読み込む数と一致しなくてはならないらしい
  real(DP) :: myz_Rad(1:205) = 0.0d0
  ! 読み込んだ値が新しいかどうか判断するための変数
  ! real(DP) :: myz_Rad_prevVal = 0.0d0
  ! 地表面温度 を読み込むための単精度変数 / 2 次元配列変数
  real(DP) :: my_sfcTemp = 0.0d0
  ! real(DP), allocatable :: my_xy_sfcTemp(:,:)
  ! 読み込んだ値が新しいかどうか判断するための変数
  ! real(DP) :: my_sfcTemp_prevVal = 0.0d0
  ! ローカルタイムを記録する変数
  ! integer           :: localTime = 7
  ! 配列操作用のループ変数
  integer           :: myxx = 1
  integer           :: myyy = 1
  integer           :: myzz = 1
  ! 時間変数比較のため integer にするための変数
  ! integer           :: myInt = 0
  ! 切り出す時間を指定する変数
  real(DP) :: myTime = 0.0d0;

  !------------------------------------------
  ! 初期化手続き ; Initialize procedure 
  !
  call MainInit

  !------------------------------------------
  ! 時間積分 time integration 
  !
  if (myrank == 0) call MessageNotify( "M", "main", "Time Integration Start" )

  ! CPU 時間計測開始
  ! Start CPU time counting 
  !
  call ClocksetLoopStart
  
  ! 時間発展ループのスタート
  !
  do while (TimeN <= EndTime ) 

    !-----------------------------------------
    ! 移流拡散.
    ! Advection and diffusion
    !
    call Dynamics_Long_forcing(       &
      &   pyz_VelXBl,  pyz_VelXNl,    & ! (in)
      &   xqz_VelYBl,  xqz_VelYNl,    & ! (in)
      &   xyr_VelZBl,  xyr_VelZNl,    & ! (in)
      &   xyz_PTempBl, xyz_PTempNl,   & ! (in) 
      &   xyz_ExnerBl, xyz_ExnerNl,   & ! (in)
      &   xyzf_QMixBl, xyzf_QMixNl,   & ! (in)
      &   xyz_KmBl,    xyz_KmNl,      & ! (in) 
      &   pyz_DVelXDtNl,              & ! (inout)
      &   xqz_DVelYDtNl,              & ! (inout)
      &   xyr_DVelZDtNl,              & ! (inout)
      &   xyz_DPTempDtNl,             & ! (inout)
      &   xyz_DExnerDtNl,             & ! (inout)
      &   xyzf_DQMixDtNl,             & ! (inout)
      &   xyz_DKmDtNl                 & ! (inout)
      )
    
    !-------------------------------
    ! 物理過程: 乱流
    !
    select case ( IDTurbMethod )
      
    case ( IDTurbMethodKW1978 )

      call turbulence_KW1978_forcing(                      &
        &   pyz_VelXBl,  xqz_VelYBl,  xyr_VelZBl,          &!(in)
        &   xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl,         &!(in)
        &   xyz_KmBl,    xyz_KhBl,                         &!(in)
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,  xyr_DVelZDtNl,  &!(inout)
        &   xyz_DPTempDtNl,xyz_DExnerDtNl, xyzf_DQMixDtNl, &!(inout)
        &   xyz_DKmDtNl,                                   &!(inout)
        &   xyz_KmAl, xyz_KhAl                             &!(out)
        & )

    case ( IDTurbMethodConstKm )

      call turbulence_ConstKm_forcing(                     &
        &   pyz_VelXBl,  xqz_VelYBl,  xyr_VelZBl,          &!(in)
        &   xyz_PTempBl,                                   &!(in)
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,  xyr_DVelZDtNl,  &!(inout)
        &   xyz_DPTempDtNl,xyz_DExnerDtNl,                 &!(inout)
        &   xyz_KmAl, xyz_KhAl                             &!(out)
        & )
    end select
    
    !-------------------------------
    ! 物理過程: 放射
    !
    select case (IDRadMethod)
      
    case (IDRadMethodHeatConst)
      
      !-------------------------------
      ! 1 時間ごとに温度分布を与えるように変更 (mkuriki)
      ! 時刻 7:00 から計算するとする. (データがローカルタイムで記録されている)
      ! localTime の初期値は 7, 最大 30
      !

      ! 時間のループ変数を比較のため, 整数型に変換
      ! myInt = nint(TimeN)

      ! 整数型にしてから剰余で 1 時間を判定
      ! if (mod(myInt, 3600) == 0) then
      ! 整数型にしてから剰余で 5 分を判定
      ! if (mod(myInt, 300) == 0) then
      ! 一つ前の値を保存
      ! myz_Rad_prevVal = myz_Rad(1)
      ! qprof.nc の読み込み
      ! gtioTime ='t='//trim(toChar(localTime))
      ! 6 時から始めるとして, 現在の積分時刻をローカルタイムに変換して読み込む値とする
      ! gtool の仕様上, 指定した値にもっとも近い値を使うようなので, 
      ! 読み込む値が 150 秒ずれこんでいるため, 補正をかけている. 
      ! さらにさらに出力間隔と読み込み間隔が微妙にずれ込むことによってグラフがガタつくので
      ! 5 秒程度の微妙な時間差をつけて補正する. 
      ! gtioTime ='t='//trim(toChar(6.0d0 + (TimeN + 150.0d0)/3600.d0))
      myTime = (TimeN + 145.0d0)/3600.d0
      ! 時刻が6 < myTime <= 30 となるように調整
      myTime = myTime - 24.0d0*(int(myTime/24)) + 6.0d0
      ! 0 時スタートに変更
      ! if (myTime < 6) then
      !   myTime = myTime - 24.0d0*(int(myTime/24)) + 24.0d0
      ! else
      !   myTime = myTime - 24.0d0*(int(myTime/24))
      ! endif
      gtioTime ='t='//trim(toChar(myTime))
      ! call HistoryGet('ncSrc/qprof.nc', 'dptrad', myz_Rad, range=gtioTime) !   1 時間用
      call HistoryGet('ncSrc/qprof_05m.nc', 'dptrad', myz_Rad, range=gtioTime) !   5 分用
      ! ローカルタイムのカウントアップ
      ! このあと地表面温度においてもローカルタイムを増やすので, ここでは処理しない
      ! localTime = localTime + 1

      ! 元データの下層における格子点は 1.5625, 4.6875, 9.375, 18.75, 37.5, 75, 150, 250, 350 ... m のため
      ! 計算に合わせて, 50(75), 150, 250 ... のみを取り出すためのループ
      ! 使う配列要素は (6:200)
      ! do myxx = imin, imax
      !   do myyy = jmin, jmax
      !     do myzz = kmin, kmax
      ! 前回の値と異なる場合にのみ配列に代入し直す
      ! if (myz_Rad_prevVal /= myz_Rad(1)) then
        do myxx = 1, nx
          do myyy = 1, ny
            do myzz = 1, nz
              xyz_PTempRadConst(myxx, myyy, myzz) = myz_Rad(myzz + 5)
            end do
          end do
        end do

      ! endif ! 放射分布の読み込み終わり
      
      xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempRadConst
      xyz_DExnerDtNl = xyz_DExnerDtNl + xyz_ExnerRadConst
      
      call HistoryAutoPut(TimeN, 'PTempRad', xyz_PTempRadConst(1:nx,1:ny,1:nz))
      call HistoryAutoPut(TimeN, 'ExnerRad', xyz_ExnerRadConst(1:nx,1:ny,1:nz))
      
    case (IDRadMethodHeatVary)

      xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempRadVary
      xyz_DExnerDtNl = xyz_DExnerDtNl + xyz_ExnerRadVary

      call HistoryAutoPut(TimeN, 'PTempRad', xyz_PTempRadVary(1:nx,1:ny,1:nz))
      call HistoryAutoPut(TimeN, 'ExnerRad', xyz_ExnerRadVary(1:nx,1:ny,1:nz))

    case (IDRadMethodSounding)

      xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempRadSndg
      xyz_DExnerDtNl = xyz_DExnerDtNl + xyz_ExnerRadSndg

      call HistoryAutoPut(TimeN, 'PTempRad', xyz_PTempRadSndg(1:nx,1:ny,1:nz))
      call HistoryAutoPut(TimeN, 'ExnerRad', xyz_ExnerRadSndg(1:nx,1:ny,1:nz))

    case (IDRadMethodBaker1998)

      xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_PTempRadBS1998
      xyz_DExnerDtNl = xyz_DExnerDtNl + xyz_ExnerRadBS1998

      call HistoryAutoPut(TimeN, 'PTempRad', xyz_PTempRadBS1998(1:nx,1:ny,1:nz))
      call HistoryAutoPut(TimeN, 'ExnerRad', xyz_ExnerRadBS1998(1:nx,1:ny,1:nz))

    case (IDRadMethodHeatBalance)
      call radiation_heatbalance_forcing(   &
        &   xyz_ExnerNl, xyz_PTempNl,       & !(in)
        &   xyz_DPTempDtNl, xyz_DExnerDtNl  & !(inout)
        & )

    end select

    !--------------------------------
    ! 境界からの熱・運動量輸送
    !
    select case (IDSurfaceMethod)

    case (IDSurfaceMethodDiff)
      call Surfaceflux_Diff_forcing(        &
        &   xyz_PTempNl, xyzf_QMixNl,       & !(in)
        &   xyz_DPTempDtNl, xyz_DExnerDtNl, & !(inout)
        &   xyzf_DQMixDtNl                  & !(inout)
        & )
     
    case (IDSurfaceMethodBulk)
      !-------------------------------
      ! 1 時間ごとに温度分布を与えるように変更 (mkuriki)
      ! 時刻 7:00 から計算するとする. (データがローカルタイムで記録されている)
      ! localTime の初期値は 7, 最大 30
      ! 放射を読み込んでいるところで, 計算タイムステップは整数型 myInt に変換済み
      !
      ! 剰余で 1 時間を判定
      ! if (mod(myInt, 3600) == 0) then
      ! 剰余で 5 分を判定
      ! if (mod(myInt, 300) == 0) then
      ! tsfc.nc の読み込み
      ! 一つ前の値を保存
      ! my_sfcTemp_prevVal = my_sfcTemp
      ! gtioTime ='t='//trim(toChar(localTime))
      ! call HistoryGet('ncSrc/tsfc.nc', 'tsfc', my_sfcTemp, range=gtioTime) !   1 時間用
      call HistoryGet('ncSrc/tsfc_05m.nc', 'tsfc', my_sfcTemp, range=gtioTime) !   5 分用
      ! ローカルタイムのカウントアップ
      ! localTime = localTime + 1

      ! 読み込んだ地面温度を 2次元配列に代入するためのループ
      ! 前回の値と異なる場合にのみ配列に代入し直す
      ! if (my_sfcTemp_prevVal /= my_sfcTemp) then
      !   do myxx = 1, nx
      !     do myyy = 1, ny
      !         my_xy_sfcTemp(myxx, myyy) = my_sfcTemp
      !     end do
      !   end do

      ! endif ! 地面温度の読み込み終わり

      !--------------------------------
      ! 地表面温度を与える関数に変更 (mkuriki)
      ! 地表面温度を表す 2次元配列 my_xy_sfcTemp を (in) な引数として追加
      ! 地表面温度を表す 1次元配列 my_sfcTemp を (in) な引数として追加
      ! 
      ! call Surfaceflux_Bulk_forcing(             &
      !   &   pyz_VelXNl, xqz_VelYNl, xyz_PTempNl, &!(in)
      !   &   xyz_ExnerNl, xyzf_QMixNl,            &!(in)
      !   &   pyz_DVelXDtNl, xqz_DVelYDtNl,        &!(inout)
      !   &   xyz_DPTempDtNl, xyz_DExnerDtNl,      &!(inout)
      !   &   xyzf_DQMixDtNl                       &!(inout)
      !   & )
      call Surfaceflux_Bulk_forcing(             &
        &   pyz_VelXNl, xqz_VelYNl, xyz_PTempNl, &!(in)
        &   xyz_ExnerNl, xyzf_QMixNl,            &!(in)
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,        &!(inout)
        &   xyz_DPTempDtNl, xyz_DExnerDtNl,      &!(inout)
        &   xyzf_DQMixDtNl,                       &!(inout)
        ! &   my_xy_sfcTemp &!(in)
        &   my_sfcTemp &!(in)
        & )

    case (IDSurfaceMethodConst)
      call Surfaceflux_Const_forcing(            &
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,        &!(inout)
        &   xyz_DPTempDtNl, xyz_DExnerDtNl,      &!(inout)
        &   xyzf_DQMixDtNl                       &!(inout)
        & )

    case (IDSurfaceMethodBaker1998)
      call Surfaceflux_bs1998_forcing(           &
        &   xyz_DPTempDtNl, xyz_DExnerDtNl       &!(inout)
        & )
      
    end select
    
    !------------------------------------------
    ! スポンジ層; sponge layer
    !
    call SpongeLayer_forcing(                          &
      &   pyz_VelXBl, xqz_VelYBl, xyr_VelZBl,          & !(in)
      &   xyz_PTempBl, xyz_ExnerBl,                    & !(in)
      &   pyz_DVelXDtNl, xqz_DVelYDtNl, xyr_DVelZDtNl, & !(inout)
      &   xyz_DPTempDtNl, xyz_DExnerDtNl               & !(inout)
      & )
    
    !-----------------------------------------
    ! 長い時間ステップでの積分
    ! Integration
    !
    call Dynamics_Long_integrate(           &
      & xyz_PTempBl, xyz_DPTempDtNl,        & !(IN)
      & xyzf_QMixBl, xyzf_DQMixDtNl,        & !(IN)
      & xyz_PTempAl,                        & !(OUT)
      & xyzf_QMixAl                         & !(OUT)
      )
    
    !------------------------------------------
    ! 凝結過程. Al な値を入れ替え.
    ! 
    select case (IDCloudMethod)
    case (IDCloudMethodK1969)
      call Cloudphys_K1969_forcing(   &
        &   xyz_ExnerNl,              &!(in)
        &   xyzf_QMixNl,              &!(in)
        &   xyz_DExnerDtNl,           &!(inout)
        &   xyz_PTempAl, xyzf_QMixAl  &!(inout)
        & )
    end select

    ! 短い時間ステップの初期値作成.
    ! Initial values set up for time integration with short time step.
    !
    pyz_VelXNs  = pyz_VelXBl
    xqz_VelYNs  = xqz_VelYBl
    xyr_VelZNs  = xyr_VelZBl
    xyz_ExnerNs = xyz_ExnerBl
    xyz_PTempNs = xyz_PTempBl
    xyz_CDensNs = xyz_CDensBl

    ! DEBUG: 長い時間ステップでの tendency を無視する場合
    !
    select case (IDDebugMethod)
    case (IDDebugMethodNoTendencyLong)
      pyz_DVelXDtNl  = 0.0d0
      xqz_DVelYDtNl  = 0.0d0
      xyr_DVelZDtNl  = 0.0d0
      xyz_DPTempDtNl = 0.0d0
      xyz_DExnerDtNl = 0.0d0
      xyz_DExnerDtNs = 0.0d0
      xyzf_DQMixDtNl = 0.0d0
      xyz_DKmDtNl    = 0.0d0
      xyz_DCDensDtNl = 0.0d0
    end select
    
    ! 短い時間ステップの時間積分. オイラー法を利用.
    ! Time integration with short time step.
    !
    Euler: do tau = 1, NstepShort

      ! 火星計算の場合. 凝結量の評価はここで行う. 
      ! 
      select case (IDCloudMethod)
      case (IDCloudMethodMarsCond)

        call cloudphys_marscond_forcing(  &
          &   xyz_PTempNs,         &  !(in) 温位
          &   xyz_ExnerNs,         &  !(in) エクスナー関数
          &   xyz_CDensNs,         &  !(in) 
          &   xyz_DPTempDtNl,      &  !(in)    
          &   xyz_DExnerDtNl,      &  !(in)    
          &   xyz_DCDensDtNl,      &  !(in)    
          &   xyz_PTempAs,         &  !(out) 温位
          &   xyz_CDensAs,         &  !(out) 雲密度
          &   xyz_DExnerDtNs       &  !(out) 
          & )

      end select

      ! HE-VI : 速度 u, v の計算.
      !
      call Dynamics_Short_forcing( &
        &   pyz_VelXNs,          & ! (in)
        &   xqz_VelYNs,          & ! (in)
        &   xyr_VelZNs,          & ! (in)
        &   xyz_ExnerNs,         & ! (in)
        &   pyz_DVelXDtNl,       & ! (in)
        &   xqz_DVelYDtNl,       & ! (in)
        &   xyr_DVelZDtNl,       & ! (in)
        &   xyz_DExnerDtNl,      & ! (in)
        &   xyz_DExnerDtNs,      & ! (in)
        &   pyz_VelXAs,          & ! (out)
        &   xqz_VelYAs,          & ! (out)
        &   xyr_VelZAs,          & ! (out)
        &   xyz_ExnerAs          & ! (out)
        & )
      
      ! 短い時間ステップのループを回すための処置
      ! Renew prognostic variables for next short time step integration.
      !
      xyz_ExnerNs  = xyz_ExnerAs
      pyz_VelXNs   = pyz_VelXAs
      xqz_VelYNs   = xqz_VelYAs
      xyr_VelZNs   = xyr_VelZAs
      xyz_PTempNs  = xyz_PTempAs
      xyz_CDensNs  = xyz_CDensAs

      !
      ! clear tendency
      !
      xyz_DExnerDtNs = 0.0d0

    end do Euler
    
    ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
    ! Renew prognostic variables for next long time step integration.
    !
    xyz_ExnerAl  = xyz_ExnerAs
    pyz_VelXAl   = pyz_VelXAs
    xqz_VelYAl   = xqz_VelYAs
    xyr_VelZAl   = xyr_VelZAs
    select case (IDCloudMethod)
    case (IDCloudMethodMarsCond)
      xyz_PTempAl = xyz_PTempAs
      xyz_CDensAl = xyz_CDensAs
    end select

    ! 速度場を初期速度で固定する場合
    !
    select case (IDDebugMethod)
    case (IDDebugMethodWindConst)
      pyz_VelXAl = pyz_VelXBl
      xqz_VelYAl = xqz_VelYBl
      xyr_VelZAl = xyr_VelZBl
    end select
   
    !
    ! clear tendency
    !
    pyz_DVelXDtNl = 0.0d0
    xqz_DVelYDtNl = 0.0d0
    xyr_DVelZDtNl = 0.0d0
    xyz_DKmDtNl   = 0.0d0
    xyz_DPTempDtNl = 0.0d0
    xyz_DExnerDtNl = 0.0d0
    xyz_DCDensDtNl = 0.0d0
    xyzf_DQMixDtNl = 0.0d0

    ! 時間フィルタ. 
    ! Time filter. 
    !
    call AsselinTimeFilter 

    ! ヒストリファイル出力. 時間フィルタをかけた後の数値を出力. 
    ! Out put to history file.
    !
    xyz_VelXNl = xyz_avr_pyz(pyz_VelXNl)
    xyz_VelYNl = xyz_avr_xqz(xqz_VelYNl)
    xyz_VelZNl = xyz_avr_xyr(xyr_VelZNl)
    
    call HistoryAutoPut(TimeN, 'PTemp', xyz_PTempNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'PTempAll', xyz_PTempNl(1:nx, 1:ny, 1:nz) + xyz_PTempBZ(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'Exner', xyz_ExnerNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'ExnerAll', xyz_ExnerNl(1:nx, 1:ny, 1:nz) + xyz_ExnerBZ(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'VelX',  xyz_VelXNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'VelY',  xyz_VelYNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'VelZ',  xyz_VelZNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'Km',    xyz_KmNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'Kh',    xyz_KhNl(1:nx, 1:ny, 1:nz))
    call HistoryAutoPut(TimeN, 'CDens', xyz_CDensNl(1:nx, 1:ny, 1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s)), xyzf_QMixNl(1:nx, 1:ny, 1:nz, s))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'All', xyzf_QMixNl(1:nx, 1:ny, 1:nz, s) + xyzf_QMixBZ(1:nx, 1:ny, 1:nz, s))
    end do

    !------------------------------------------
    ! 保存量の出力
    ! Out put conservation variables
    !
    call EnergyMonit_exec(                    &
      &  pyz_VelXNl,  xqz_VelYNL, xyr_VelZNl, & !(in)
      &  xyz_ExnerNl, xyz_PTempNl             & !(in)
      & )
    
    !----------------------------------------------------
    ! リスタートファイルの作成
    ! Make restartfile.
    !    
    if (mod(t, NstepOutput) == 0) then 

      ! 移流に対する CFL 条件のチェック 
      ! CFL condtion check for advection
      !
      call CFLCheckTimeLongVelX( pyz_VelXNl ) !(in)
      call CFLCheckTimeLongVelY( xqz_VelYNl ) !(in)
      call CFLCheckTimeLongVelZ( xyr_VelZNl ) !(in)
      
      ! リスタートファイルに変数を出力
      !
      call HistoryPut( 't',     TimeN,       rstat)
      call HistoryPut( 'VelX',  pyz_VelXNl,  rstat)
      call HistoryPut( 'VelY',  xqz_VelYNl,  rstat)
      call HistoryPut( 'VelZ',  xyr_VelZNl,  rstat)
      call HistoryPut( 'Exner', xyz_ExnerNl, rstat)
      call HistoryPut( 'PTemp', xyz_PTempNl, rstat)
      call HistoryPut( 'Km',    xyz_KmNl,    rstat)
      call HistoryPut( 'Kh',    xyz_KhNl,    rstat)
      call HistoryPut( 'CDens', xyz_CDensNl, rstat)
      call HistoryPut( 'QMix',  xyzf_QMixNl, rstat)    

      call HistoryPut( 't',     TimeA,       rstat)
      call HistoryPut( 'VelX',  pyz_VelXAl,  rstat)
      call HistoryPut( 'VelY',  xqz_VelYAl,  rstat)
      call HistoryPut( 'VelZ',  xyr_VelZAl,  rstat)
      call HistoryPut( 'Exner', xyz_ExnerAl, rstat)
      call HistoryPut( 'PTemp', xyz_PTempAl, rstat)
      call HistoryPut( 'Km',    xyz_KmAl,    rstat)
      call HistoryPut( 'Kh',    xyz_KhAl,    rstat)
      call HistoryPut( 'CDens', xyz_CDensAl, rstat)
      call HistoryPut( 'QMix',  xyzf_QMixAl, rstat) 

      ! Show CPU time 
      ! 
      call ClocksetPredict
      
    end if

    ! 長い時間ステップのループを回すための処置.
    ! Renew prognostic variables for next long time step integration.
    !
    pyz_VelXBl  = pyz_VelXNl;  pyz_VelXNl  = pyz_VelXAl
    xqz_VelYBl  = xqz_VelYNl;  xqz_VelYNl  = xqz_VelYAl
    xyr_VelZBl  = xyr_VelZNl;  xyr_VelZNl  = xyr_VelZAl
    xyz_PTempBl = xyz_PTempNl; xyz_PTempNl = xyz_PTempAl
    xyz_ExnerBl = xyz_ExnerNl; xyz_ExnerNl = xyz_ExnerAl
    xyz_KmBl    = xyz_KmNl;    xyz_KmNl    = xyz_KmAl
    xyz_KhBl    = xyz_KhNl;    xyz_KhNl    = xyz_KhAl
    xyz_CDensBl = xyz_CDensNl; xyz_CDensNl = xyz_CDensAl
    xyzf_QMixBl = xyzf_QMixNl; xyzf_QMixNl = xyzf_QMixAl
    t = t + 1

    ! 時刻の進行
    ! Progress time
    !
    call TimesetProgress

  end do

  !----------------------------------------------------
  ! 出力ファイルのクローズ
  ! Close out put files.
  !
  call HistoryFileio_finalize
  call ReStartFileio_finalize

  !----------------------------------------------------
  ! MPI END
  !
  call MPIWrapperFinalize

  !----------------------------------------------------
  ! CPU 時間の計測終了
  ! 
  call ClocksetLoopStop
  call ClocksetClose
  
contains
!-----------------------------------------------------------------------
  subroutine VariableAllocate
    !
    !初期化として, 配列を定義し, 値としてゼロを代入する.
    !

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate( pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( pyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_VelXNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( pyz_VelXAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( pyz_VelXNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( pyz_VelXAs(imin:imax,jmin:jmax,kmin:kmax) )

    pyz_VelXBl  = 0.0d0;    pyz_VelXNl  = 0.0d0;    pyz_VelXAl  = 0.0d0
    pyz_VelXNs  = 0.0d0;    pyz_VelXAs = 0.0d0    
    xyz_VelXNl  = 0.0d0

    allocate( xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xqz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_VelYNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xqz_VelYAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xqz_VelYNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xqz_VelYAs(imin:imax,jmin:jmax,kmin:kmax) )

    xqz_VelYBl  = 0.0d0;    xqz_VelYNl  = 0.0d0;    xqz_VelYAl  = 0.0d0
    xqz_VelYNs  = 0.0d0;    xqz_VelYAs = 0.0d0    
    xyz_VelYNl  = 0.0d0

    allocate( xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_VelZNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_VelZNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_VelZAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_VelZNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_VelZAs(imin:imax,jmin:jmax,kmin:kmax) )

    xyr_VelZBl  = 0.0d0;    xyr_VelZNl  = 0.0d0;    xyr_VelZAl  = 0.0d0
    xyr_VelZNs  = 0.0d0;    xyr_VelZAs = 0.0d0
    xyz_VelZNl  = 0.0d0

    allocate( xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_ExnerNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_ExnerAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_ExnerNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_ExnerAs(imin:imax,jmin:jmax,kmin:kmax) )

    xyz_ExnerBl = 0.0d0;    xyz_ExnerNl = 0.0d0;    xyz_ExnerAl = 0.0d0
    xyz_ExnerNs = 0.0d0;    xyz_ExnerAs = 0.0d0

    allocate( xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_PTempNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_PTempAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_PTempNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_PTempAs(imin:imax,jmin:jmax,kmin:kmax) )

    xyz_PTempBl = 0.0d0;  xyz_PTempNl = 0.0d0;  xyz_PTempAl = 0.0d0
    xyz_PTempNs = 0.0d0;  xyz_PTempAs = 0.0d0

    allocate( xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_CDensNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_CDensAl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_CDensNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_CDensAs(imin:imax,jmin:jmax,kmin:kmax) )

    xyz_CDensBl = 0.0d0;  xyz_CDensNl = 0.0d0;  xyz_CDensAl = 0.0d0
    xyz_CDensNs = 0.0d0;  xyz_CDensAs = 0.0d0

    allocate( xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_KmNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax) )

    xyz_KmBl    = 0.0d0;    xyz_KmNl    = 0.0d0;    xyz_KmAl    = 0.0d0

    allocate( xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_KhNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax) )

    xyz_KhBl    = 0.0d0;    xyz_KhNl    = 0.0d0;    xyz_KhAl    = 0.0d0

    allocate( xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax,ncmax) )
    allocate( xyzf_QMixNl(imin:imax,jmin:jmax,kmin:kmax,ncmax) )
    allocate( xyzf_QMixAl(imin:imax,jmin:jmax,kmin:kmax,ncmax) )

    xyzf_QMixBl = 0.0d0;   xyzf_QMixNl = 0.0d0;   xyzf_QMixAl = 0.0d0

    allocate( pyz_DVelXDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xqz_DVelYDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyr_DVelZDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_DKmDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_DPTempDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_DExnerDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_DExnerDtNs(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyz_DCDensDtNl(imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xyzf_DQMixDtNl(imin:imax,jmin:jmax,kmin:kmax,ncmax) )

    pyz_DVelXDtNl = 0.0d0
    xqz_DVelYDtNl = 0.0d0
    xyr_DVelZDtNl = 0.0d0
    xyz_DKmDtNl   = 0.0d0
    xyz_DPTempDtNl = 0.0d0
    xyz_DExnerDtNl = 0.0d0
    xyz_DExnerDtNs = 0.0d0
    xyz_DCDensDtNl = 0.0d0
    xyzf_DQMixDtNl = 0.0d0

    ! ----------------------------------
    ! 地表面温度を読み込むための 2次元配列を初期化 (mkuriki)
    ! allocate( my_xy_sfcTemp(imin:imax,jmin:jmax) )
    ! my_xy_sfcTemp = 0.0d0

    ! 時間ループの初期化
    !
    t = 1

  end subroutine VariableAllocate


  !-----------------------------------------------------------------------
  subroutine AsselinTimeFilter
    !
    ! 時間フィルター; Asselin のタイムフィルターを利用
    !   t = 0.0 の場合には tfil = 0.0d0, それ以外は tfil = 1.0d-1
    !   (t = 0 の時はオイラー法で, それ以外はリープフロッグ法で積分するため)
    !
    use TimeSet, only: tfil

    ! 暗黙の型宣言禁止    
    !
    implicit none

    pyz_VelXNl  =  pyz_VelXNl   + tfil * (pyz_VelXAl  - 2.0d0 * pyz_VelXNl  + pyz_VelXBl)
    xqz_VelYNl  =  xqz_VelYNl   + tfil * (xqz_VelYAl  - 2.0d0 * xqz_VelYNl  + xqz_VelYBl)
    xyr_VelZNl  =  xyr_VelZNl   + tfil * (xyr_VelZAl  - 2.0d0 * xyr_VelZNl  + xyr_VelZBl)
    xyz_PTempNl =  xyz_PTempNl  + tfil * (xyz_PTempAl - 2.0d0 * xyz_PTempNl + xyz_PTempBl)
    xyz_ExnerNl =  xyz_ExnerNl  + tfil * (xyz_ExnerAl - 2.0d0 * xyz_ExnerNl + xyz_ExnerBl)
    xyz_KmNl    =  xyz_KmNl     + tfil * (xyz_KmAl    - 2.0d0 * xyz_KmNl    + xyz_KmBl)
    xyz_KhNl    =  xyz_KhNl     + tfil * (xyz_KhAl    - 2.0d0 * xyz_KhNl    + xyz_KhBl)
    xyz_CDensNl =  xyz_CDensNl  + tfil * (xyz_CDensAl - 2.0d0 * xyz_CDensNl + xyz_CDensBl)
    xyzf_QMixNl =  xyzf_QMixNl  + tfil * (xyzf_QMixAl - 2.0d0 * xyzf_QMixNl + xyzf_QMixBl)
    
  end subroutine AsselinTimeFilter


  !-----------------------------------------------------------------------
  subroutine MainInit
    
    ! 暗黙の型宣言禁止    
    !
    implicit none

    character(STRING) :: cfgfile
                             ! NAMELIST ファイル名 ; NAMELIST fine name

    ! MPI
    !
    call MPIWrapperInit

    ! 時間計測
    !
    call ClocksetInit      !初期化
    call ClocksetPreStart  !初期化ルーチン計測開始
    
    ! NAMELIST ファイル名の読み込み
    ! Loading NAMELIST file.
    !
    call argset_init( cfgfile ) !(out)

    ! NAMELIST ファイル名のモジュールへの登録
    ! Loading NAMELIST file.
    !
    call NmlutilInit( cfgfile ) !(in)
    
    ! 時間積分の初期化
    ! Initialization of time integration.
    !
    call timeset_init

    ! 格子点情報の初期化
    ! Initialization of grid arrangement.
    !
    call gridset_init
    
    ! 化学計算ルーチンの初期化
    ! Initialization of chemical routines.
    !
    call chemcalc_init
    
    ! 定数の情報の初期化
    ! Initialization of constant variables.
    !
    call constants_init

    ! 軸の情報の初期化
    ! Initialization of axis variables.
    !
    call axesset_init
    
    ! I/O ファイル名の初期化
    ! Initialization of output file name. 
    !
    call fileset_init
    
    ! 湿潤過程共有変数の初期化
    ! Initialization of common variables for moist process.
    !
    call composition_init

    ! ヒストリーファイル・リスタートファイルの初期化
    ! Initialize restart & history files.
    !
    call HistoryFileio_init
    call ReStartFileio_init

    ! 数値摩擦係数の初期化
    ! Initialization of numerical friction coefficient.
    !
    call Damping_Init

    ! マージンの設定の初期化
    ! Initialization of margin
    !
    call SetMargin_Init
   
    ! 内部変数の初期化
    ! Initialization of internal variables.
    !
    call VariableAllocate

    ! 初期値の代入 
    ! * ReStartFile が設定されている場合にはファイルを読み込む. 
    !   設定されていない場合にはデフォルトの基本場と擾乱場を作る. 
    !
    ! Initial value set up.
    ! * Read restartfile if it is specified. If not, make default basic
    !   state and disturbance.
    !
    if (myrank == 0) call MessageNotify( "M", "main", "Initial value setup." )

    ! 基本場, 擾乱場の初期値を netCDF ファイルから取得する.
    ! 
    call ReStartFileio_BZ_Get
    call ReStartFileio_Var_Get(     &
      & pyz_VelXBl,  pyz_VelXNl,    & ! (out)
      & xqz_VelYBl,  xqz_VelYNl,    & ! (out)
      & xyr_VelZBl,  xyr_VelZNl,    & ! (out)
      & xyz_PTempBl, xyz_PTempNl,   & ! (out)
      & xyz_ExnerBl, xyz_ExnerNl,   & ! (out)
      & xyzf_QMixBl, xyzf_QMixNl,   & ! (out)
      & xyz_KmBl,    xyz_KmNl,      & ! (out)
      & xyz_KhBl,    xyz_KhNl,      & ! (out)
      & xyz_CDensBl, xyz_CDensNl    & ! (out)
      & )
    
    ! 力学モジュールの初期化
    ! Initialization of dynamical modules
    !
    call Dynamics_init
    
    ! 負の湿潤量の補填計算の初期化
    ! Initialization of negative moist value correction.
    !
    call FillNegative_Init
   
    ! 物理過程の設定
    !
    call deepconv_main

    !-------------------------------------------------------------
    ! 基本場のファイル出力
    !
    call HistoryPut( 'DensBZ',     xyz_DensBZ,     rstat)
    call HistoryPut( 'ExnerBZ',    xyz_ExnerBZ,    rstat)
    call HistoryPut( 'PTempBZ',    xyz_PTempBZ,    rstat)
    call HistoryPut( 'VPTempBZ',   xyz_VPTempBZ,   rstat)
    call HistoryPut( 'VelSoundBZ', xyz_VelSoundBZ, rstat)
    call HistoryPut( 'TempBZ',     xyz_TempBZ,     rstat)
    call HistoryPut( 'PressBZ',    xyz_PressBZ,    rstat)
    call HistoryPut( 'QMixBZ',     xyzf_QMixBZ,    rstat)
    call HistoryPut( 'EffMolWtBZ', xyz_EffMolWtBZ, rstat)
    
    ! 音波に対する CFL 条件のチェック
    ! CFL condtion check for sound wave.
    !
    call CFLCheckTimeShort( &
      & xyz_VelSoundBZ   & ! (in)
      & )
    
    ! 保存量のモニタリング
    !
    call EnergyMonit_init

    ! 化学計算ルーチンの初期化(2)
    ! 鉛直方向に計算をさぼる場合に利用する初期化ルーチン. 
    ! 基本場が決まらないと初期化出来ない. 
    !
    call chemcalc_init2

    ! t=0 から数値積分を行う場合は, 最初の一回目はオイラー法を利用. 
    !
    if ( FlagInitialRun ) call TimesetDelTimeHalf

    ! 初期化ルーチンの CPU 時間の計測終了
    !
    call ClocksetPreStop

    
  end subroutine MainInit
  

  subroutine deepconv_main

    use dc_message,    only: MessageNotify
    use dc_iounit,     only : FileOpen    

    implicit none

    integer            :: unit                    !装置番号
    character(STRING)  :: FlagTurbMethod    = ""  !乱流拡散に関する設定
    character(STRING)  :: FlagRadMethod     = ""  !放射に関する設定
    character(STRING)  :: FlagCloudMethod   = ""  !雲微物理に関する設定
    character(STRING)  :: FlagSurfaceMethod = ""  !地表面過程に関する設定
    character(STRING)  :: FlagWindMethod    = ""  !速度場に関する設定
    character(STRING)  :: FlagDebugMethod   = ""  !デバッグ用のフラグ

    NAMELIST /deepconv_main_nml / &
      & FlagTurbMethod,           &
      & FlagRadMethod,            &
      & FlagCloudMethod,          &
      & FlagSurfaceMethod,        &
      & FlagWindMethod,           &
      & FlagDebugMethod

    ! デフォルト値の設定
    ! Default values settings
    !
    FlagTurbMethod       = "Nothing"
!!$    FlagTurbMethod    = "KW1978"
!!$    FlagTurbMethod    = "ConstKm"

    FlagRadMethod        = "Nothing"
!!$    FlagRadMethod     = "HeatConst"
!!$    FlagRadMethod     = "HeatVary"
!!$    FlagRadMethod     = "HeatBalance"

    FlagSurfaceMethod    = "Nothing"
!!$    FlagSurfaceMethod = "Diff"
!!$    FlagSurfaceMethod = "Bulk"
!!$    FlagSurfaceMethod = "Const"
!!$    FlagSurfaceMethod = "Baker1998"

    FlagCloudMethod      = "Nothing"
!!$    FlagCloudMethod   = "K1969"
!!$    FlagCloudMethod   = "MarsCond"

    FlagWindMethod       = "Nothing"
!!$    FlagWindMethod    = "Const"

    FlagDebugMethod      = "Nothing"
!!$    FlagDebugMethod   = "NoTendencyLong"


    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=deepconv_main_nml)
    close(unit)

    select case ( FlagTurbMethod )
    case ( "Nothing" )
    case ( "KW1978" )
      IDTurbMethod = IDTurbMethodKW1978 
      call turbulence_kw1978_init
    case ( "ConstKm" )
      IDTurbMethod = IDTurbMethodConstKm
      call turbulence_constKm_init
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagTurbMethod=<%c> is not supported.', &
        & c1 = trim(FlagTurbMethod) )
    end select

    select case ( FlagRadMethod )
    case ( "Nothing" )
    case ( "HeatConst" )
      IDRadMethod = IDRadMethodHeatConst 
      call Radiation_Simple_init
    case ( "HeatVary" )
      IDRadMethod = IDRadMethodHeatVary
      call Radiation_Simple_init
    case ( "HeatBalance" )
      IDRadMethod = IDRadMethodHeatBalance
      call radiation_heatbalance_init
    case ( "Sounding" )
      IDRadMethod = IDRadMethodSounding
      call radiation_sounding_init
    case ( "Baker1998" )
      IDRadMethod = IDRadMethodBaker1998
      call Radiation_BS1998_init
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagRadMethod=<%c> is not supported.', &
        & c1 = trim(FlagRadMethod) )
    end select

    select case ( FlagSurfaceMethod )
    case ( "Nothing" )
    case ( "Diff" )
      IDSurfaceMethod = IDSurfaceMethodDiff
      call Surfaceflux_Diff_init
    case ( "Bulk" )
      IDSurfaceMethod = IDSurfaceMethodBulk
      call Surfaceflux_Bulk_init
    case ( "Const" )
      IDSurfaceMethod = IDSurfaceMethodConst
      call Surfaceflux_Const_init
    case ( "Baker1998" )
      IDSurfaceMethod = IDSurfaceMethodBaker1998
      call Surfaceflux_bs1998_init
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagSurfaceMethod=<%c> is not supported.', &
        & c1 = trim(FlagSurfaceMethod) )
    end select

    select case ( FlagCloudMethod )
    case ( "Nothing" )
    case ( "K1969" )
      IDCloudMethod = IDCloudMethodK1969
      call cloudphys_K1969_init
    case ( "MarsCond" )
      IDCloudMethod = IDCloudMethodMarsCond
      call cloudphys_marscond_init
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagCloudMethod=<%c> is not supported.', &
        & c1 = trim(FlagCloudMethod) )
    end select

    select case ( FlagWindMethod )
    case ( "Nothing" )
    case ( "Const" )
      IDDebugMethod = IDDebugMethodWindConst
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagWindMethod=<%c> is not supported.', &
        & c1 = trim(FlagWindMethod) )
    end select

    select case ( FlagDebugMethod )
    case ( "Nothing" )
    case ( "WindConst" )
      IDDebugMethod = IDDebugMethodWindConst
    case ( "NoTendencyLong" )
      IDDebugMethod = IDDebugMethodNoTendencyLong
    case default
      call MessageNotify( 'E', prog_name, &
        & 'FlagDebugdMethod=<%c> is not supported.', &
        & c1 = trim(FlagWindMethod) )
    end select

    if (myrank == 0) then 
      call MessageNotify( "M", "main", "FlagTurbMethod    = %c", c1=trim(FlagTurbMethod))
      call MessageNotify( "M", "main", "IDTurbMethod      = %d", i=(/IDTurbMethod/))
      call MessageNotify( "M", "main", "FlagRadMethod     = %c", c1=trim(FlagRadMethod))
      call MessageNotify( "M", "main", "IDRadMethod       = %d", i=(/IDRadMethod/))
      call MessageNotify( "M", "main", "FlagSurfaceMethod = %c", c1=trim(FlagSurfaceMethod))
      call MessageNotify( "M", "main", "IDSurfaceMethod   = %d", i=(/IDSurfaceMethod/))
      call MessageNotify( "M", "main", "FlagCloudMethod   = %c", c1=trim(FlagCloudMethod))
      call MessageNotify( "M", "main", "IDCloudMethod     = %d", i=(/IDCloudMethod/))
      call MessageNotify( "M", "main", "FlagDebugMethod   = %c", c1=trim(FlagDebugMethod))
      call MessageNotify( "M", "main", "IDDebugMethod     = %d", i=(/IDDebugMethod/))
    end if

  end subroutine deepconv_main

      
end program deepconv_arare
