!= deepconv/arare 湿潤大気対流計算用主プログラム (三次元版)
!
!= deepconv/arare main program for moist atmospheric convection (3D)
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare.f90,v 1.19 2011-06-14 02:14:22 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
  use gtool_history, only: HistoryPut
  use gtool_historyauto, only: HistoryAutoPut
  use dc_clock, only : CLOCK, DCClockCreate, &
    & DCClockClose, DCClockStart, DCClockStop, &
    & DCClockResult, DCClockPredict, &
    & operator(+)                             ! Access module (モジュール指定)

  ! 初期設定モジュール
  ! Initialize module
  !
  use mpi_wrapper,   only : MPIWrapperInit, MPIWrapperFinalize, myrank
  use argset,   only: argset_init
  use gridset,  only: gridset_init, imin, imax, jmin, jmax, kmin, kmax, nx, ny, nz, ncmax
  use timeset,  only: timeset_init, DelTimeLong, DelTimeShort, EndTime, RestartTime, &
    &                 NstepLong, NstepShort, NstepRStat
  use axesset,  only: axesset_init, xyz_avr_pqz, 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

  ! 下請けモジュール
  ! Utility modules
  !
  use cflcheck, only : CFLCheckTimeShort, &
    &                  CFLCheckTimeLongVelX, &
    &                  CFLCheckTimeLongVelY, &
    &                  CFLCheckTimeLongVelZ
  use timefilter, only : AsselinFilter
  use damping, only: Damping_init, DampSponge_xyz, DampSponge_xyr, &
    &                DampSponge_pyz,DampSponge_xqz

  ! 力学過程計算用関数モジュール
  ! Dynamical processes module
  !
  use DynamicsHEVI, only: Dynamics_Init,  &
    &                      Dynamics_Long_forcing, Dynamics_Short_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_RadHeatConst, xyz_RadHeatVary
  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

  ! 湿潤過程計算用モジュール
  ! 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(STRING) :: cfgfile
                             ! NAMELIST ファイル名 ; NAMELIST fine 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) :: TimeN, TimeA   ! 時刻 ; Time 
  real(DP), allocatable :: DelTimeLFrog(:)      
                             ! リープフロッグスキーム用時間格子間隔
                             ! Time interval for Leap-frog scheme
  real(DP) :: DelTimeEuler   ! オイラースキーム用時間格子
                             ! Time interval for Eular scheme
  integer, allocatable :: NStepEuler(:) 
                             ! オイラースキーム用時間ステップ数
                             ! The number of time step for Eular scheme
  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 :: t, tau, f   ! do ループ変数 ; do loop variable 

  character(STRING) :: FlagTurbulence = ""
  character(STRING) :: FlagRadiation  = ""
  character(STRING) :: FlagCloudMicroPhys = ""
  character(STRING) :: FlagSurfaceHeating = ""
  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           :: IDCloudMethod = 0
  integer, parameter:: IDCloudMethodK1969 = 1
  integer, parameter:: IDCloudMethodMarsCond = 2
  type(CLOCK)       :: clock_init, clock_loop  ! Variables for CPU time counting 
                                               ! CPU 時間計測用変数

  !------------------------------------------
  ! 初期化手続き ; Initialize procedure 
  !
  call MainInit
  
  !------------------------------------------
  ! 時間積分 time integration 
  !
  call MessageNotify( "M", "main", "Time Integration Start" )

  ! 時刻の初期化
  !
  TimeN = RestartTime
  TimeA = RestartTime + DelTimeLong
  t = 1

  ! 時間発展ループのスタート
  !
  do while (TimeA <= EndTime + DelTimeLong) 

    ! CPU Time
    call DCClockStart ( clk = clock_loop ) ! (inout) ! Start CPU time counting 
                                                     ! (CPU 時間計測開始)
    
    !-------------------------------
    ! 物理過程: 乱流
    !
    select case ( IDTurbMethod )
      
    case ( IDTurbMethodKW1978 )
      call turbulence_KW1978_forcing(                      &
        &   TimeN, DelTimeLFrog(t),                         &!(in)
        &   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_RadHeatConst
      call HistoryAutoPut(TimeN, 'PTempRad', xyz_RadHeatConst(1:nx,1:ny,1:nz))

    case (IDRadMethodHeatVary)

      xyz_DPTempDtNl = xyz_DPTempDtNl + xyz_RadHeatVary
      call HistoryAutoPut(TimeN, 'PTempRad', xyz_RadHeatVary(1:nx,1:ny,1:nz))

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

    end select

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

    case (IDSurfaceMethodDiff)
      call Surfaceflux_Diff_forcing(   &
        &   TimeN, xyz_PTempNl, xyzf_QMixNl, & !(in)
        &   xyz_DPTempDtNl, xyzf_DQMixDtNl   & !(out)
        & )
     
    case (IDSurfaceMethodBulk)
      call Surfaceflux_Bulk_forcing( &
        &   TimeN,                               &!(in)
        &   pyz_VelXNl, xqz_VelYNl, xyz_PTempNl, &!(in)
        &   xyz_ExnerNl, xyzf_QMixNl,            &!(in)
        &   pyz_DVelXDtNl, xqz_DVelYDtNl,        &!(inout)
        &   xyz_DPTempDtNl, xyzf_DQMixDtNl       &!(inout)
        & )
      
    end select


    !-----------------------------------------
    ! 凝結過程
    ! 
    !
    select case (IDCloudMethod)
    case (IDCloudMethodK1969)
      call CloudPhys_K1969_FallRain( & 
        &   TimeN, xyzf_QMixNl,       & !(in)
        &   xyzf_DQMixDtNl           & !(inout)
        & )
    end select
    
    !-----------------------------------------
    ! 移流拡散.
    ! Advection and diffusion
    !
    call Dynamics_Long_forcing(       &
      &   TimeN, DelTimeLFrog(t),     & ! (in)
      &   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_KmBl,    xyz_KmNl,      & ! (in)
      &   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_DKmDtNl,                & ! (inout)
      &   xyz_DCDensDtNl,             & ! (inout)
      &   xyz_PTempAl,                & ! (out)
      &   xyzf_QMixAl                 & ! (out)
      & )

    !------------------------------------------
    ! 凝結過程
    ! 
    select case (IDCloudMethod)
    case (IDCloudMethodK1969)
      call Cloudphys_K1969_forcing(   &
        &   TimeN, DelTimeLFrog(t),    &!(in)
        &   xyz_ExnerNl,              &!(in)
        &   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.
    !
    Euler: do tau = 1, NstepEuler(t)

      ! 火星計算の場合. 凝結量の評価はここで行う. 
      ! 
      ! * 確認事項
      !  * 凝結量と潜熱加熱は Ns, As のどちらを使って見積もるのか?
      !
      select case (IDCloudMethod)
      case (IDCloudMethodMarsCond)

        call cloudphys_marscond_forcing(  &
          &   TimeN,               &  !(in) 時刻
          &   DelTimeEuler,        &  !(in) 時間間隔
          &   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(         &
        &   TimeN, DelTimeEuler,  & ! (in)
        &   pyz_VelXNs,          & ! (in)
        &   xqz_VelYNs,          & ! (in)
        &   xyr_VelZNs,          & ! (in)
        &   xyz_ExnerNs,         & ! (in)
        &   pyz_DVelXDtNl,       & ! (in)
        &   xqz_DVelYDtNl,       & ! (in)
        &   xyr_DVelZDtNl,       & ! (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
      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. 
    !
    if (RestartTime == 0.0d0 .AND. t == 1) then 
      write(*,*) "skip: ", RestartTime, t    
    else
      call AsselinFilter( &
        & pyz_VelXAl,     & ! (in)
        & pyz_VelXNl,     & ! (inout)
        & pyz_VelXBl      & ! (in)
        & )
      call AsselinFilter( &
        & xqz_VelYAl,     & ! (in)
        & xqz_VelYNl,     & ! (inout)
        & xqz_VelYBl      & ! (in)
        & )
      call AsselinFilter( &
        & xyr_VelZAl,     & ! (in)
        & xyr_VelZNl,     & ! (inout)
        & xyr_VelZBl      & ! (in)
        & )
      call AsselinFilter( &
        & xyz_PTempAl,  & ! (in)
        & xyz_PTempNl,  & ! (inout)
        & xyz_PTempBl   & ! (in)
        & )
      call AsselinFilter( &
        &  xyz_ExnerAl,   & ! (in)
        &  xyz_ExnerNl,   & ! (inout)
        &  xyz_ExnerBl    & ! (in)
        & )         
      call AsselinFilter( &
        & xyz_KmAl,       & ! (in)
        & xyz_KmNl,       & ! (inout)
        & xyz_KmBl        & ! (in)
        & )
      call AsselinFilter( &
        & xyz_CDensAl,       & ! (in)
        & xyz_CDensNl,       & ! (inout)
        & xyz_CDensBl        & ! (in)
        & )
      do f = 1, ncmax
        call AsselinFilter(         &
          &   xyzf_QMixAl(:,:,:,f), &
          &   xyzf_QMixNl(:,:,:,f), &
          &   xyzf_QMixBl(:,:,:,f)  &
          & )
      end do
    end if

    ! スポンジ層
    ! Numerical dumping.
    !
    call DampSponge_pyz( pyz_VelXAl,  pyz_VelXBl,  DelTimeLFrog(t) )
    call DampSponge_xqz( xqz_VelYAl,  xqz_VelYBl,  DelTimeLFrog(t) )
    call DampSponge_xyr( xyr_VelZAl,  xyr_VelZBl,  DelTimeLFrog(t) )
    call DampSponge_xyz( xyz_PTempAl, xyz_PTempBl, DelTimeLFrog(t) )
    call DampSponge_xyz( xyz_ExnerAl, xyz_ExnerBl, DelTimeLFrog(t) )
    call DampSponge_xyz( xyz_KmAl,    xyz_KmBl,    DelTimeLFrog(t) )
    
    ! 移流に対する 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_pqz(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 f = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f)), xyzf_QMixNl(1:nx, 1:ny, 1:nz, f))
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'All', xyzf_QMixNl(1:nx, 1:ny, 1:nz, f) + xyzf_QMixBZ(1:nx, 1:ny, 1:nz, f))
    end do
    
    !----------------------------------------------------
    ! リスタートファイルの作成
    ! Make restartfile.
    !    
    if (t /= 1 .AND. mod(t, NStepRStat) == 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) 

      ! CPU Time
      ! 
      if (myrank == 0) then 
        write(*,*) ""
        call MessageNotify( "M", "main", "Time = %f", d=(/TimeA/) )
        
        call DCClockStop( clk = clock_loop ) ! (inout)  ! Stop CPU time counting 
                                                        ! (CPU 時間計測終了)
        call DCClockPredict( &                       ! Estimate remaining time (残り時間の予測)
          &   clk = clock_init + clock_loop,       & ! (in)
          &   progress = real(t) / real(NStepLong) & ! (in) 
          &  )
      end if
    end if

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

    t = t + 1
    TimeN = TimeN + DelTimeLong
    TimeA = TimeA + DelTimeLong

  end do

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

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

  !----------------------------------------------------
  ! CPU Time
  ! 
  call DCClockResult( &                    ! Display total CPU time (全 CPU 時間の表示)
    & clks = (/clock_init, clock_loop/), & ! (in)
    & total_auto = .true. )                ! (in)
  call DCClockClose( clk = clock_init )    ! (inout)       ! Finalize (後処理)
  call DCClockClose( clk = clock_loop )    ! (inout)       ! Finalize (後処理)

  
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( DelTimeLFrog(NstepLong+1) )
    allocate( NStepEuler(NStepLong+1) )

    DelTimeEuler = 0.0d0
    DelTimeLFrog = 0.0d0 
    NStepEuler = 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

  end subroutine VariableAllocate
  
  !-----------------------------------------------------------------------
  subroutine MainInit
    implicit none

    ! CPU Time
    !
    call DCClockCreate( &           ! Initialize (初期化)
      & clk = clock_init, &         ! (out)
      & name = 'initialization' )   ! (in)
    call DCClockCreate( &           ! Initialize (初期化)
      & clk = clock_loop, &         ! (out)
      & name = 'time-integration' ) ! (in)

    call DCClockStart( clk = clock_init ) ! (inout)   ! Start CPU time counting 
                                                      ! (CPU 時間計測開始)    

    ! MPI
    !
    call MPIWrapperInit

    ! NAMELIST ファイル名の読み込み
    ! Loading NAMELIST file.
    !
    call argset_init( cfgfile ) !(out)
    
    ! 時間積分の初期化
    ! Initialization of time integration.
    !
    call timeset_init( cfgfile ) !(in)
    
    ! 格子点情報の初期化
    ! Initialization of grid arrangement.
    !
    call gridset_init( cfgfile ) !(in)
    
    ! 化学計算ルーチンの初期化
    ! Initialization of chemical routines.
    !
    call chemcalc_init( )
    
    ! 定数の情報の初期化
    ! Initialization of constant variables.
    !
    call constants_init( cfgfile ) !(in)

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

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

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

    
    ! 内部変数の初期化
    ! 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.
    !
    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( cfgfile ) !(in)  
    
    ! 負の湿潤量の補填計算の初期化
    ! Initialization of negative moist value correction.
    !
    call FillNegative_Init() 
   

    ! 物理過程の設定
    !
    !
    call deepconv_main( cfgfile ) !(in) 

       
    ! 時刻とループ回数の初期化
    ! Initialization of time integration.
    !
    NstepEuler   = NstepShort 
    DelTimeLFrog = DelTimeLong * 2.0d0 
    DelTimeEuler = DelTimeShort
    
    if ( RestartTime == 0.0d0 ) then    
      NstepEuler(1)   = NstepShort /2       
      DelTimeLFrog(1) = DelTimeLong         
    end if
    
    call MessageNotify( "M", "main", "NstepLong= %d", i=(/NstepLong/) )
    call MessageNotify( "M", "main", "NstepEuler= %d", i=(/NstepEuler/) )
    call MessageNotify( "M", "main", "DelTimeLFrog= %f", d=(/DelTimeLFrog/) )
    call MessageNotify( "M", "main", "DelTimeEuler= %f", d=(/DelTimeEuler/) )
    
    !-------------------------------------------------------------
    ! 基本場のファイル出力
    !
    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 DCClockStop( clk = clock_init ) ! (inout)    ! Stop CPU time counting 
                                                      ! (CPU 時間計測終了)
    
  end subroutine MainInit
  

  subroutine deepconv_main( cfgfile  )

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

    implicit none

    character(STRING), intent(in) :: cfgfile  ! nml file
    integer                       :: unit     !装置番号

    NAMELIST /deepconv_main_nml/  &
      & FlagTurbulence,            &
      & FlagRadiation,             &
      & FlagCloudMicroPhys,        &
      & FlagSurfaceHeating

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


    if (FlagTurbulence == "KW1978") then
      IDTurbMethod = IDTurbMethodKW1978 
      call turbulence_kw1978_init( cfgfile )
    end if

    if (FlagRadiation == "HeatConst") then 
      IDRadMethod = IDRadMethodHeatConst 
      call Radiation_Simple_init(cfgfile)

    elseif (FlagRadiation == "HeatVary") then 
      IDRadMethod = IDRadMethodHeatVary
      call Radiation_Simple_init(cfgfile)

    elseif (FlagRadiation == "HeatBalance") then 
      IDRadMethod = IDRadMethodHeatBalance
      call radiation_heatbalance_init(cfgfile)
    end if

    if (FlagSurfaceHeating == "Diff") then 
      IDSurfaceMethod = IDSurfaceMethodDiff
      call Surfaceflux_Diff_init( cfgfile )

    elseif (FlagSurfaceHeating == "Bulk") then 
      IDSurfaceMethod = IDSurfaceMethodBulk
      call Surfaceflux_Bulk_init( cfgfile )

    end if

    if (FlagCloudMicroPhys == "K1969") then 
      IDCloudMethod = IDCloudMethodK1969
      call cloudphys_K1969_init(cfgfile)

    elseif (FlagCloudMicroPhys == "MarsCond") then 
      IDCloudMethod = IDCloudMethodMarsCond
      call cloudphys_marscond_init(cfgfile)
    end if

    call MessageNotify( "M", "main", "FlagTurbulence     = %c", c1=trim(FlagTurbulence)) 
    call MessageNotify( "M", "main", "IDTurbMethod       = %d", i=(/IDTurbMethod/))
    call MessageNotify( "M", "main", "FlagRadiation      = %c", c1=trim(FlagRadiation)) 
    call MessageNotify( "M", "main", "IDRadMethod        = %d", i=(/IDRadMethod/))
    call MessageNotify( "M", "main", "FlagSurfaceHeating = %c", c1=trim(FlagSurfaceHeating)) 
    call MessageNotify( "M", "main", "IDSurfaceMethod    = %d", i=(/IDSurfaceMethod/))
    call MessageNotify( "M", "main", "FlagCloudMicroPhys = %c", c1=trim(FlagCloudMicroPhys)) 
    call MessageNotify( "M", "main", "IDCloudMethod      = %d", i=(/IDCloudMethod/))
    
  end subroutine deepconv_main

      
end program deepconv_arare
