!= deepconv/arare 湿潤大気対流計算用主プログラム (三次元版)
!
!= deepconv/arare main program for moist atmospheric convection (3D)
!
! Authors::   杉山耕一朗(SUGIYAMA Ko-ichiro), ODAKA Masatsugu
! Version::   $Id: arare.f90,v 1.32 2012-06-01 05:29:15 odakker Exp $
! Tag Name::  $Name: arare5-20120828 $
! 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
  use gtool_history, only: HistoryPut
  use gtool_historyauto, only: HistoryAutoPut

  ! 初期設定モジュール
  ! Initialize module
  !
  use mpi_wrapper, only : MPIWrapperInit, MPIWrapperFinalize, myrank
  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, TimeB, NstepShort, NstepOutput, &
    &                 EndTime, RestartTime, &
    &                 DelTimeLong, DelTimeShort, DelTimeOutput, 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
  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

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

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

  ! 放射強制計算用モジュール
  ! Radiative forceing module
  !
  use Radiation_simple,  only: Radiation_simple_init, xyz_PTempRadConst, xyz_PTempRadVary, &
    &                          xyz_ExnerRadConst, xyz_ExnerRadVary
  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, only: Surfaceflux_Bulk_init, Surfaceflux_Bulk_forcing
  use Surfaceflux_Const, only: Surfaceflux_Const_init, Surfaceflux_Const_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           :: IDRadMethod = 0
  integer, parameter:: IDRadMethodHeatConst = 1
  integer, parameter:: IDRadMethodHeatVary  = 2
  integer, parameter:: IDRadMethodHeatBalance = 3
  integer           :: IDSurfaceMethod = 0
  integer, parameter:: IDSurfaceMethodDiff = 1
  integer, parameter:: IDSurfaceMethodBulk = 2
  integer, parameter:: IDSurfaceMethodConst = 3
  integer           :: IDCloudMethod = 0
  integer, parameter:: IDCloudMethodK1969 = 1
  integer, parameter:: IDCloudMethodMarsCond = 2

  !------------------------------------------
  ! 初期化手続き ; 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 )

    !-------------------------------
    ! 物理過程: 乱流
    !
    select case ( IDTurbMethod )

    case ( IDTurbMethodKW1978 )

      ! Km の移流
      !
      call Dynamics_Km_forcing(               &
        & pyz_VelXNl, xqz_VelYNl, xyr_VelZNl, & ! (in)
        & xyz_KmBl,    xyz_KmNl,              & ! (in)
        & xyz_DKmDtNl                         & ! (inout)
        & )

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

    !-------------------------------
    ! 物理過程: 放射
    !
    select case (IDRadMethod)

    case (IDRadMethodHeatConst)

      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 (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)
      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)
        & )
    case (IDSurfaceMethodConst)
      call Surfaceflux_Const_forcing(       &
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,        &!(inout)
        &   xyz_DPTempDtNl, xyz_DExnerDtNl,      &!(inout)
        &   xyzf_DQMixDtNl                       &!(inout)
        & )

    end select

    !-----------------------------------------
    ! 凝結過程
    !
    !
    select case (IDCloudMethod)
    case (IDCloudMethodK1969)
      call CloudPhys_K1969_FallRain( &
        &   xyzf_QMixNl,       & !(in)
        &   xyz_DExnerDtNl,    & !(inout)
        &   xyzf_DQMixDtNl     & !(inout)
        & )
    end select

    !------------------------------------------
    ! スポンジ層; sponge layer
    !
    !!!!!!!!!!!!!!!!!!!!!!!
    !
    ! test
    !
    !!!!!!!!!!!!!!!!!!!!!!!

    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)
      & )

    !-----------------------------------------
    ! 移流拡散.
    ! 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)
      &   xyzf_QMixBl, xyzf_QMixNl,   & ! (in) mixing ratio  [kg/kg]
      &   xyz_CDensBl, xyz_CDensNl,   & ! (in) Cloud density [kg/m^-3]
      &   pyz_DVelXDtNl,              & ! (inout)
      &   xqz_DVelYDtNl,              & ! (inout)
      &   xyr_DVelZDtNl,              & ! (inout)
      &   xyz_DPTempDtNl,             & ! (inout)
      &   xyzf_DQMixDtNl,             & ! (inout)
      &   xyz_DCDensDtNl,             & ! (inout)
      &   xyz_PTempAl,                & ! (out)
      &   xyzf_QMixAl                 & ! (out)
      & )

    !------------------------------------------
    ! 凝結過程
    !
    select case (IDCloudMethod)
    case (IDCloudMethodK1969)
      call Cloudphys_K1969_forcing(   &
        &   xyz_ExnerNl,              &!(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

    ! 短い時間ステップの時間積分. オイラー法を利用.
    ! Time integration with short time step.
    !

    ! 圧力関数への強制は 0 にする
    !
!    xyz_DExnerDtNl = 0.0d0

    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

      ! 陽解法: 速度 u, v の計算.
      ! Time integration horizontal velocity (u).
      !
      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

    !
    ! 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

    !------------------------------------------
    ! 移流に対する CFL 条件のチェック
    ! CFL condtion check for advection
    !
!    call CFLCheckTimeLongVelX( pyz_VelXNl ) !(in)
!    call CFLCheckTimeLongVelY( xqz_VelYNl ) !(in)
!    call CFLCheckTimeLongVelZ( xyr_VelZNl ) !(in)

    ! ヒストリファイル出力. 時間フィルタをかけた後の数値を出力.
    ! 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

    !----------------------------------------------------
    ! リスタートファイルの作成
    ! Make restartfile.
    !
    if (mod(t, NstepOutput) == 0) then

      ! ファイル出力
      !
      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 Time
  !
  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

    ! 時間ループの初期化
    !
    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

    real(DP) :: aaa_VarNl(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: aaaa_VarNl(imin:imax,jmin:jmax,kmin:kmax,ncmax)

    aaa_VarNl  = pyz_VelXNl  + tfil * (pyz_VelXAl - 2.0d0 * pyz_VelXNl + pyz_VelXBl)
    pyz_VelXNl = aaa_VarNl

    aaa_VarNl  = xqz_VelYNl  + tfil * (xqz_VelYAl - 2.0d0 * xqz_VelYNl + xqz_VelYBl)
    xqz_VelYNl = aaa_VarNl

    aaa_VarNl  = xyr_VelZNl  + tfil * (xyr_VelZAl - 2.0d0 * xyr_VelZNl + xyr_VelZBl)
    xyr_VelZNl = aaa_VarNl

    aaa_VarNl  = xyz_PTempNl + tfil * (xyz_PTempAl - 2.0d0 * xyz_PTempNl + xyz_PTempBl)
    xyz_PTempNl= aaa_VarNl

    aaa_VarNl  = xyz_ExnerNl + tfil * (xyz_ExnerAl - 2.0d0 * xyz_ExnerNl + xyz_ExnerBl)
    xyz_ExnerNl= aaa_VarNl

    aaa_VarNl  = xyz_KmNl    + tfil * (xyz_KmAl - 2.0d0 * xyz_KmNl + xyz_KmBl)
    xyz_KmNl   = aaa_VarNl

    aaa_VarNl  = xyz_KhNl    + tfil * (xyz_KhAl - 2.0d0 * xyz_KhNl + xyz_KhBl)
    xyz_KhNl   = aaa_VarNl

    aaa_VarNl  = xyz_CDensNl + tfil * (xyz_CDensAl - 2.0d0 * xyz_CDensNl + xyz_CDensBl)
    xyz_CDensNl = aaa_VarNl

    aaaa_VarNl  = xyzf_QMixNl + tfil * (xyzf_QMixAl - 2.0d0 * xyzf_QMixNl + xyzf_QMixBl)
    xyzf_QMixNl = aaaa_VarNl

  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 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( '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)
      & )

    ! CPU time
    !
    call ClocksetPreStop

    ! 最初の一回目の場合はオイラー法で回す.
    !
    if ( FlagInitialRun ) call TimesetDelTimeHalf

  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 = ""

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

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

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

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

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

    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 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 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 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

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

  end subroutine deepconv_main

end program deepconv_arare