!= deepconv/arare 湿潤大気対流計算用主プログラム (三次元版)
!
!= deepconv/arare main program for moist atmospheric convection (3D)
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare_3d.f90,v 1.10 2009-02-28 08:19:10 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

program 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, &
    &                 DimXMin, DimXMax, DimYMin, DimYMax, DimZMin, DimZMax, &
    &                 SpcNum, FileXMin, FileXMax, FileYMin, FileYMax, FileZMin, FileZMax
  use timeset,  only: timeset_init, TimesetDelTimeHalf, TimesetProgress, &
    &                 TimeA, TimeN, TimeB, NstepShort, NstepOutput, &
    &                 EndTime, RestartTime, tfil, &
    &                 DelTimeLong, DelTimeShort, DelTimeOutput, FlagInitialRun
  use axesset,  only: axesset_init
  use constants,only: constants_init
  use composition, only: composition_init, SpcWetSymbol
  use fileset,  only: fileset_init
  use basicset_3d, only: xyza_MixRtBasicZ, xyz_DensBasicZ, xyz_EffMolWtBasicZ, &
    &                 xyz_PotTempBasicZ, xyz_TempBasicZ, xyz_PressBasicZ, &
    &                 xyz_VelSoundBasicZ, xyz_ExnerBasicZ
  use ChemCalc_3d, only: ChemCalc_init
  use namelist_util, only: NmlutilInit, namelist_filename

  ! 化学量計算モジュール
  ! Chemical calculation modules
  !
  use chemdata,      only: chemdata_init

  ! 力学過程計算用関数モジュール
  ! Dynamical processes module
  !
  use DynFunc_3d,    only: xyz_AdvScalar, xyza_AdvScalar, xyz_AdvKm, &
    &                         pyz_AdvVelX, xqz_AdvVelY,              &
    &                         xyr_Buoy, xyr_AdvVelZ, pyz_GradPi, xqz_GradPi
  use DynImpFunc_3d, only: xyz_Exner_init, xyz_Exner, xyr_GradPi

  ! 乱流拡散計算用モジュール
  ! Turbulent diffusion module
  !
  use Turbulence_3d,  only: Turbulence_Init, xyz_BuoyKm, xyz_TurbKm, &
    &                       xyz_TurbScalar, xyza_TurbScalar, pyz_TurbVelX, &
    &                       xyr_TurbVelZ  , xyz_ShearKm    , xyz_DispKm,   &
    &                       xyz_DispHeat  , xqz_TurbVelY   , EddyViscosity

!  ! 境界からのフラックス計算用モジュール
!  ! Surface flux module
!  !
!  use HeatFlux_3d,    only: xyz_HeatFluxBulk, pyz_MomFluxBulk, xqz_MomFluxBulk

  ! 境界からのフラックス計算用モジュール
  ! Surface flux module
  !
  use HeatFlux_3d,   only: xyz_HeatFluxDiff, xyza_MixRtFluxDiff

  ! 放射強制計算用モジュール
  ! Radiative forceing module
  !
  use Radiation_3d,   only: Radiation_init, xyz_RadHeatConst, xyz_NewtonCool

  ! 湿潤過程計算用モジュール
  ! Moist processes modules
  !
!  use moistset,       only: moistset_init
  use MoistAdjust_3d, only: MoistAdjustSvapPress, MoistAdjustNH4SH
  use WarmRainPrm_3d, only: WarmRainPrm_Init, xyz_Rain2GasHeat, xyza_Rain2Gas, &
    &                       xyza_Rain2GasNH4SH, xyz_Rain2GasHeatNH4SH,         &
    &                       xyza_Cloud2Rain, xyza_FallRain
  use MoistBuoyancy_3d,only: MoistBuoy_Init, xyz_BuoyMoistKm, xyr_BuoyMolWt, &
    &                        xyr_BuoyDrag
  use fillnegative_3d,only: FillNegative_init, xyza_FillNegative_xyza

  ! 安定度の計算
  ! Static stability calculation module
  !
!  use ECCM_3D,        only: ECCM_Stab

  ! 数値拡散/摩擦計算用モジュール
  ! Numerical diffussion /dumping module
  !
  use NumDiffusion4th_3d,only: NumDiffusion_Init, xyz_NumDiffScalar, &
    &                       xyz_NumDiffKm, xyza_NumDiffScalar, &
    &                       pyz_NumDiffVelX, xqz_NumDiffVelY,  &
    &                       xyr_NumDiffVelZ
  use damping_3d,     only: damping_init, xyz_Sponge, pyz_Sponge, &
    &                       xqz_Sponge, xyr_Sponge

  ! ファイル入出力モジュール
  ! File I/O module
  !
  use RestartFileIO_3d, only : ReStartFileio_init, ReStartFileio_Finalize, &
    &                          ReStartFileio_BasicZ_Get, ReStartFileio_Var_Get, rstat
  use HistoryFileIO_3d, only: HistoryFileIO_init, HistoryFileIO_finalize

  ! 下請けモジュール
  ! Utility modules
  !
  use cflcheck,       only : CFLCheckTimeShort, &
    &                        CFLCheckTimeLongVelX, &
    &                        CFLCheckTimeLongVelY, &
    &                        CFLCheckTimeLongVelZ
  use xyz_bc_module,  only : BoundaryXCyc_xyz, BoundaryYCyc_xyz, &
    &                        BoundaryZSym_xyz, BoundaryZCyc_xyz, &
    &                        BoundaryXCyc_pyz, BoundaryYCyc_pyz, &
    &                        BoundaryZSym_pyz, BoundaryZCyc_pyz, &
    &                        BoundaryXCyc_xqz, BoundaryYCyc_xqz, &
    &                        BoundaryZSym_xqz, BoundaryZCyc_xqz, &
    &                        BoundaryXCyc_xyr, BoundaryYCyc_xyr, &
    &                        BoundaryZSym_xyr, BoundaryZCyc_xyr, &
    &                        BoundaryZAsym_xyr
  use xyz_base_module, only: xyz_avr_pyz, xyz_avr_xyr, xyz_avr_xqz

  implicit none

  ! 内部変数
  ! Internal variables
  !
  real(DP), allocatable :: xyz_VelXNl(:,:,:)
  real(DP), allocatable :: xyz_VelYNl(:,:,:)
  real(DP), allocatable :: xyz_VelZNl(:,:,:)

  real(DP), allocatable :: pyz_VelXBl(:,:,:)
                             ! $ u (t-\Delta t) $ 東西風 ; zonal wind
  real(DP), allocatable :: pyz_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 :: 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 :: 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_PotTempBl(:,:,:)
                             ! $ \theta (t-\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PotTempNl(:,:,:)
                             ! $ \theta (t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PotTempAl(:,:,:)
                             ! $ \theta (t+\Delta t) $ 温位 ; Potential temp.
  real(DP), allocatable :: xyz_PotTempWork(:,:,:)
                             ! 温位 $ \theta $ の作業配列 ; Work array
  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 :: xyza_MixRtBl(:,:,:,:)
                             ! $ q (t-\Delta t) $ 湿潤量の混合比
                             ! Mixing ratio of moist variables.
  real(DP), allocatable :: xyza_MixRtNl(:,:,:,:)
                             ! $ q (t) $ 湿潤量の混合比
                             ! Mixing ratio of moist variables
  real(DP), allocatable :: xyza_MixRtAl(:,:,:,:) !
                             ! $ q (t+\Delta t) $ 湿潤量の混合比
                             !Mixing ratio of moist variables
  real(DP), allocatable :: xyza_MixRtWork(:,:,:,:)
                             ! 湿潤量の作業配列
                             ! Work array for mixing ratio.
  real(DP), allocatable :: pyz_AccelVelXNl(:,:,:)
                             ! 気圧傾度力を除いた $u$ の変化率
                             ! Tendency of $u$ except for pressure gradient term
  real(DP), allocatable :: xqz_AccelVelYNl(:,:,:)
                             ! 気圧傾度力を除いた $v$ の変化率
                             ! Tendency of $v$ except for pressure gradient term
  real(DP), allocatable :: xyr_AccelVelZNl(:,:,:)
                             ! 気圧傾度力を除いた $w$ の変化率
                             ! Tendency of $w$ except for pressure gradient term
  real(DP), allocatable :: xyza_DelMixRt(:,:,:,:)
                             ! 湿潤量の混合比 $ q $ の増分
                             ! Mixing ratio variation.

  real(DP)              :: DelTimeLFrog
                             ! リープフロッグスキーム用時間格子間隔
                             ! Time interval for Leap-frog scheme

  integer :: &
    & t,     &
    & tau,   & ! do ループ変数 ; do loop variable
    & t1,    & ! do ループ変数 ; do loop variable
    & t2,    & ! do ループ変数 ; do loop variable
    & s        ! do ループ変数 ; do loop variable

  real(8), allocatable :: Km1(:,:,:), Km2(:,:,:), Km3(:,:,:), Km4(:,:,:), Km5(:,:,:), Km6(:,:,:), Km7(:,:,:)

  logical :: CloudFlag = .false.

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

  !------------------------------------------
  ! 時間積分 time integration
  !
!  do t1 = 1, NstepLFrog / NstepDisp
!    do t2 = 1, NstepDisp
!
  t = 1
  do while (TimeN <= EndTime )

      ! 時刻の設定
      ! Time setting.
      !
!      Time = Time + DelTimeLong
!      t = (t1 - 1) * NstepDisp + t2
      DelTimeLFrog = DelTimeLong * 2.0d0

      ! 渦拡散係数の移流拡散
      ! Advection and diffusion of turbulent diffusion coefficient.
      !
      xyz_KmAl =                                                          &
        & xyz_KmBl                                                        &
        & + DelTimeLFrog                                                  &
        &   * (                                                           &
        &     + xyz_AdvKm(xyz_KmNl, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)   &
        &     + xyz_NumDiffKm(xyz_KmBl)                                   &
        &     + xyz_BuoyMoistKm(xyz_PotTempBl, xyz_ExnerBl, xyza_MixRtBl)  &
        &     + xyz_ShearKm(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl) &
        &     + xyz_TurbKm(xyz_KmBl)                                      &
        &     + xyz_DispKm(xyz_KmBl)                                      &
        &    )

      ! 値の上限下限の設定
      !  * 値は正になることを保証する
      !  * 値の上限は 800 とする. 中島健介(1994, 学位論文)参照
      !
      ! Upper and lower bound value are specified.
      !
      xyz_KmAl = max( 0.0d0, min( xyz_KmAl, 800.0d0 ) )

      ! 境界条件 ; Boundary condition
      !
      call BoundaryXCyc_xyz( xyz_KmAl ) ! (inout)
      call BoundaryYCyc_xyz( xyz_KmAl ) ! (inout)
      call BoundaryZSym_xyz( xyz_KmAl ) ! (inout)

      ! スカラーに対する渦拡散係数の計算
      ! Specify turbulent diffusion coefficient for scalar variables.
      !
      xyz_KhAl = 3.0d0 * xyz_KmAl

!      write(*,*) "Km", minval( xyz_KmAl ), maxval( xyz_KmAl )

      ! 温位の移流拡散の計算
      ! Advection and diffusion of potential temperature.
      !
      xyz_PotTempAl =                                              &
        &   xyz_PotTempBl                                          &
        & + DelTimeLFrog                                           &
        &   * (                                                    &
        &     + xyz_AdvScalar( xyz_PotTempNl,                      &
        &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) &
        &     + xyz_AdvScalar( xyz_PotTempBasicZ,                  &
        &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) &
        &     + xyz_NumDiffScalar(xyz_PotTempBl)                   &
        &     + xyz_RadHeatConst( xyz_ExnerBl )                    &
        &     + xyz_Sponge( xyz_PotTempBl )                        &
        &     + xyz_TurbScalar(xyz_PotTempBl,     xyz_KhBl)        &
        &     + xyz_TurbScalar(xyz_PotTempBasicZ, xyz_KhBl)        &
        &     + xyz_DispHeat( xyz_KmBl )                           &
        &     + xyz_HeatFluxDiff( xyz_PotTempNl )                  &
!!      &     + xyz_HeatFluxBulk( xyz_PotTempBl, pyz_VelXBl, xqz_VelYBl ) &
!!      &     + xyz_NewtonCool( xyz_PotTempBl )                    &
      &      )

        Km1 = + xyz_AdvScalar( xyz_PotTempNl,                      &
        &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl) &
        &     + xyz_AdvScalar( xyz_PotTempBasicZ,                  &
        &                      pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
        Km2 = + xyz_NumDiffScalar(xyz_PotTempBl)
        Km3 = + xyz_RadHeatConst( xyz_ExnerBl )
        Km4 = + xyz_Sponge( xyz_PotTempBl )
        Km5 = + xyz_TurbScalar(xyz_PotTempBl,     xyz_KhBl)       &
        &     + xyz_TurbScalar(xyz_PotTempBasicZ, xyz_KhBl)
        Km6 = + xyz_DispHeat( xyz_KmBl )
        Km7 = + xyz_HeatFluxDiff( xyz_PotTempNl )

      ! 境界条件 ; Boundary condition
      !
      call BoundaryXCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryYCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryZSym_xyz( xyz_PotTempAl ) ! (inout)

      write(*,*) "PTemp_1", minval( Km1(1:nx,1:ny,1:nz) ), maxval( Km1(1:nx,1:ny,1:nz) ), Km1(45:51,1,1)
      write(*,*) "PTemp_2", minval( Km2(1:nx,1:ny,1:nz) ), maxval( Km2(1:nx,1:ny,1:nz) ), Km2(45:51,1,1)
      write(*,*) "PTemp_3", minval( Km3(1:nx,1:ny,1:nz) ), maxval( Km3(1:nx,1:ny,1:nz) ), Km3(45:51,1,1)
      write(*,*) "PTemp_4", minval( Km4(1:nx,1:ny,1:nz) ), maxval( Km4(1:nx,1:ny,1:nz) ), Km4(45:51,1,1)
      write(*,*) "PTemp_5", minval( Km5(1:nx,1:ny,1:nz) ), maxval( Km5(1:nx,1:ny,1:nz) ), Km5(45:51,1,1)
      write(*,*) "PTemp_6", minval( Km6(1:nx,1:ny,1:nz) ), maxval( Km6(1:nx,1:ny,1:nz) ), Km6(45:51,1,1)
      write(*,*) "PTemp_7", minval( Km7(1:nx,1:ny,1:nz) ), maxval( Km7(1:nx,1:ny,1:nz) ), Km7(45:51,1,1)

!      write(*,*) "PotTemp", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl )

      ! 凝縮成分混合比の移流拡散
      ! Advection and diffusion of vapor, cloud and rain mixing ratios.
      !
      xyza_MixRtAl =                                               &
        &   xyza_MixRtBl                                           &
        & + DelTimeLFrog                                           &
        &   * (                                                    &
        &    + xyza_AdvScalar (xyza_MixRtNl,     pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)&
        &    + xyza_AdvScalar (xyza_MixRtBasicZ, pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)&
        &    + xyza_NumDiffScalar(xyza_MixRtBl)                    &
        &    + xyza_FallRain(xyza_MixRtBl)                         &
        &    + xyza_TurbScalar(xyza_MixRtBl,    xyz_KhBl)          &
        &    + xyza_TurbScalar(xyza_MixRtBasicZ,xyz_KhBl)          &
        &    + xyza_MixRtFluxDiff(xyza_MixRtNl)                      &
!!      &    + xyza_MixRtFluxBulk(xyza_MixRtNl)                      &
        &   )

      ! 境界条件 Boundary condition
      !
      do s = 1, SpcNum
        call BoundaryXCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryYCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryZSym_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
      end do

!      write(*,*) "MixRt", minval( xyza_MixRtAl ), maxval( xyza_MixRtAl )

      ! 移流によって負になった部分を埋める
      ! Negative values due to advection are corrected.
      !
      xyza_MixRtWork = xyza_MixRtAl
      xyza_MixRtAl   = xyza_FillNegative_xyza( xyza_MixRtWork )

      ! 埋めた/削った量を保管
      ! Correction value is stored.
      !
!      call StoreMixRtFill1( &
!        & (xyza_MixRtAl - xyza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )

      ! 埋め切れなかった部分をゼロにする. Fill2 に保管
      ! Negative values mixing ratios are corrected.
      xyza_MixRtWork = xyza_MixRtAl
      xyza_MixRtAl = max( - xyza_MixRtBasicZ, xyza_MixRtWork )

      ! 境界条件 Boundary condition
      !
      do s = 1, SpcNum
        call BoundaryXCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryYCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryZSym_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
      end do

!      write(*,*) "9"
!      write(*,*) "PotTemp", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl )
!      write(*,*) "MixRt", minval( xyza_MixRtAl ), maxval( xyza_MixRtAl )
!      write(*,*) "ExnerN", minval( xyz_ExnerNl ), maxval( xyz_ExnerNl )
!      write(*,*) "ExnerB", minval( xyz_ExnerBl ), maxval( xyz_ExnerBl )

      !!!------------------------------------------------------
      !!! 凝結過程
      !!!
      if (CloudFlag) then

      ! 暖かい雨のパラメタリゼーション.
      ! * 雲<-->雨 の変換を行う.
      !
      ! Warm rain parameterization.
      ! * Conversion from cloud to rain.

      !これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xyza_MixRtWork = xyza_MixRtAl

      !雨への変化量を計算
      ! Conversion values are calculated.
      !
      xyza_MixRtAl   = xyza_MixRtWork &
        &             + xyza_Cloud2Rain( xyza_MixRtAl, DelTimeLFrog )

      ! 雲から雨への変換量を保管
      ! Conversion values are sotred.
      !
!      call StoreMixRtCond( &
!        & (xyza_MixRtAl - xyza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )

      ! 境界条件 Boundary condition
      !
      do s = 1, SpcNum
        call BoundaryXCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryYCyc_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
        call BoundaryZSym_xyz( xyza_MixRtAl(:,:,:,s) ) ! (inout)
      end do

!      write(*,*) "10"
!      write(*,*) "PotTemp", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl )
!      write(*,*) "MixRt", minval( xyza_MixRtAl ), maxval( xyza_MixRtAl )
!      write(*,*) "ExnerN", minval( xyz_ExnerNl ), maxval( xyz_ExnerNl )
!      write(*,*) "ExnerB", minval( xyz_ExnerBl ), maxval( xyz_ExnerBl )

      ! 暖かい雨のパラメタリゼーション.
      ! * 蒸気<-->雨 の変換を行う
      !
      ! Warm rain parameterization.
      ! * Conversion from rain to vapor.

      !これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xyz_PotTempWork = xyz_PotTempAl
      xyza_MixRtWork  = xyza_MixRtAl

      ! 雨から蒸気への混合比変化を求める
      ! * 温位の計算において, 混合比変化が必要となるため,
      !   混合比変化を 1 つの配列として用意する.
      !
      ! Conversion values are calculated.
      !
      xyza_DelMixRt = 0.0d0
      xyza_DelMixRt =                                                    &
        & (                                                              &
        &   + xyza_Rain2Gas(                                             &
        &       xyz_ExnerNl, xyz_PotTempAl, xyza_MixRtAl, DelTimeLFrog   &
        &     )                                                          &
        &   + xyza_Rain2GasNH4SH(                                        &
        &       xyz_ExnerNl, xyz_PotTempAl, xyza_MixRtAl, DelTimeLFrog   &
        &       )                                                        &
        &  )

      ! 温位の計算. 雨から蒸気への変換に伴う潜熱・反応熱を追加.
      !
      xyz_PotTempAl =                                                         &
        & xyz_PotTempWork                                                     &
        & + (                                                                 &
        &     + xyz_Rain2GasHeat( xyz_PotTempAl, xyz_ExnerNl, xyza_DelMixRt ) &
        &     + xyz_Rain2GasHeatNH4SH( xyz_ExnerNl, xyza_DelMixRt )           &
        &    )

      ! 混合比の計算. 雨から蒸気への変換分を追加
      !
      xyza_MixRtAl   = xyza_MixRtWork + xyza_DelMixRt

      !雨から蒸気の値の保管
      !
!      call StorePotTempCond( &
!        & (xyz_PotTempAl - xyz_PotTempWork) / DelTimeLFrog & ! (in)
!        & )
!      call StoreMixRtCond( &
!        & xyza_DelMixRt / DelTimeLFrog & ! (in)
!        & )

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryYCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryZSym_xyz( xyz_PotTempAl ) ! (inout)

      do s = 1, SpcNum
        call BoundaryXCyc_xyz( xyza_MixRtAl(:,:,:,s) )
        call BoundaryYCyc_xyz( xyza_MixRtAl(:,:,:,s) )
        call BoundaryZSym_xyz( xyza_MixRtAl(:,:,:,s) )
      end do

!      write(*,*) "11"
!      write(*,*) "PotTemp", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl )
!      write(*,*) "MixRt", minval( xyza_MixRtAl ), maxval( xyza_MixRtAl )
!      write(*,*) "ExnerN", minval( xyz_ExnerNl ), maxval( xyz_ExnerNl )
!      write(*,*) "ExnerB", minval( xyz_ExnerBl ), maxval( xyz_ExnerBl )

      ! 湿潤飽和調節
      ! * 蒸気<-->雲の変換を行う.
      !
      ! Moist adjustment.
      ! * Conversion from vapor to cloud.

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xyz_PotTempWork = xyz_PotTempAl
      xyza_MixRtWork  = xyza_MixRtAl

      ! 湿潤調節法を適用
      ! Moist adjustment is applied.
      !
      call MoistAdjustSvapPress( &
        & xyz_ExnerNl,           & ! (in)
        & xyz_PotTempAl,         & ! (inout)
        & xyza_MixRtAl           &
        & )
      call MoistAdjustNH4SH(     &
        & xyz_ExnerNl,           & ! (in)
        & xyz_PotTempAl,         & ! (inout)
        & xyza_MixRtAl           & ! (inout)
        & )

      ! 湿潤調節法による温位と混合比の変化量を保管
      ! Adjustment values are stored.
      !
!      call StorePotTempCond( &
!        & (xyz_PotTempAl - xyz_PotTempWork) / DelTimeLFrog & ! (in)
!        & )
!      call StoreMixRtCond( &
!        & (xyza_MixRtAl - xyza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryYCyc_xyz( xyz_PotTempAl ) ! (inout)
      call BoundaryZSym_xyz( xyz_PotTempAl ) ! (inout)

      do s = 1, SpcNum
        call BoundaryXCyc_xyz( xyza_MixRtAl(:,:,:,s) )
        call BoundaryYCyc_xyz( xyza_MixRtAl(:,:,:,s) )
        call BoundaryZSym_xyz( xyza_MixRtAl(:,:,:,s) )
      end do

!      write(*,*) "12"
!      write(*,*) "PotTemp", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl )
!      write(*,*) "MixRt", minval( xyza_MixRtAl ), maxval( xyza_MixRtAl )
!      write(*,*) "ExnerN", minval( xyz_ExnerNl ), maxval( xyz_ExnerNl )
!      write(*,*) "ExnerB", minval( xyz_ExnerBl ), maxval( xyz_ExnerBl )

      end if

      ! 速度の移流拡散.
      ! Advection and diffusion of velocity components.
      !
      pyz_AccelVelXNl =                                      &
        & + pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)  &
        & + pyz_Sponge(pyz_VelXBl)                           &
        & + pyz_NumDiffVelX(pyz_VelXBl)                      &
        & + pyz_TurbVelX(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl)
!        & + pyz_MomFluxBulk( pyz_VelXBl )

      Km1 = pyz_AdvVelX(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
      Km2 = pyz_Sponge(pyz_VelXBl)
      Km3 = pyz_NumDiffVelX(pyz_VelXBl)
      Km4 = pyz_TurbVelX(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl)

      write(*,*) "DVelXDt_long_1", minval( Km1(1:nx,1:ny,1:nz) ), maxval( Km1(1:nx,1:ny,1:nz) ), Km1(45:51,1,1)
      write(*,*) "DVelXDt_long_2", minval( Km2(1:nx,1:ny,1:nz) ), maxval( Km2(1:nx,1:ny,1:nz) ), Km2(45:51,1,1)
      write(*,*) "DVelXDt_long_3", minval( Km3(1:nx,1:ny,1:nz) ), maxval( Km3(1:nx,1:ny,1:nz) ), Km3(45:51,1,1)
      write(*,*) "DVelXDt_long_4", minval( Km4(1:nx,1:ny,1:nz) ), maxval( Km4(1:nx,1:ny,1:nz) ), Km4(45:51,1,1)

      xqz_AccelVelYNl =                                      &
        & + xqz_AdvVelY(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)  &
        & + xqz_Sponge(xqz_VelYBl)                           &
        & + xqz_NumDiffVelY(xqz_VelYBl)                      &
        & + xqz_TurbVelY(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl)
!        & + xqz_MomFluxBulk( xqz_VelYBl )

      xyr_AccelVelZNl =                                      &
        & + xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)  &
        & + xyr_Buoy(xyz_PotTempNl)                          &
        & + xyr_Sponge(xyr_VelZBl)                           &
        & + xyr_NumDiffVelZ(xyr_VelZBl)                      &
        & + xyr_TurbVelZ(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl)   &
        & + xyr_BuoyMolWt(xyza_MixRtNl)                                  &
        & + xyr_BuoyDrag(xyza_MixRtNl)

      Km1 = xyr_AdvVelZ(pyz_VelXNl, xqz_VelYNl, xyr_VelZNl)
      Km2 = xyr_Buoy(xyz_PotTempNl)
      Km3 = xyr_Sponge(xyr_VelZBl)
      Km4 =  xyr_NumDiffVelZ(xyr_VelZBl)
      Km5 = xyr_TurbVelZ(xyz_KmBl, pyz_VelXBl, xqz_VelYBl, xyr_VelZBl)
      Km6 = xyr_BuoyMolWt(xyza_MixRtNl)
      Km7 = xyr_BuoyDrag(xyza_MixRtNl)

      write(*,*) "DVelZDt_long_1", minval( Km1(1:nx,1:ny,0:nz) ), maxval( Km1(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_1 0  ", Km1(45:51,1,0)
      write(*,*) "DVelZDt_long_1 nz ", Km1(45:51,1,nz)
      write(*,*) "DVelZDt_long_2", minval( Km2(1:nx,1:ny,0:nz) ), maxval( Km2(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_2 0  ", Km2(45:51,1,0)
      write(*,*) "DVelZDt_long_2 nz ", Km2(45:51,1,nz)
      write(*,*) "DVelZDt_long_3", minval( Km3(1:nx,1:ny,0:nz) ), maxval( Km3(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_3 0  ", Km3(45:51,1,0)
      write(*,*) "DVelZDt_long_3 nz ", Km3(45:51,1,nz)
      write(*,*) "DVelZDt_long_4", minval( Km4(1:nx,1:ny,0:nz) ), maxval( Km4(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_4 0  ", Km4(45:51,1,0)
      write(*,*) "DVelZDt_long_4 nz ", Km4(45:51,1,nz)
      write(*,*) "DVelZDt_long_5", minval( Km5(1:nx,1:ny,0:nz) ), maxval( Km5(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_5 0  ", Km5(45:51,1,0)
      write(*,*) "DVelZDt_long_5 nz ", Km5(45:51,1,nz)
      write(*,*) "DVelZDt_long_6", minval( Km6(1:nx,1:ny,0:nz) ), maxval( Km6(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_6 0  ", Km6(45:51,1,0)
      write(*,*) "DVelZDt_long_6 nz ", Km6(45:51,1,nz)
      write(*,*) "DVelZDt_long_7", minval( Km7(1:nx,1:ny,0:nz) ), maxval( Km7(1:nx,1:ny,0:nz) )
      write(*,*) "DVelZDt_long_7 0  ", Km7(45:51,1,0)
      write(*,*) "DVelZDt_long_7 1  ", Km7(45:51,1,1)
      write(*,*) "DVelZDt_long_7 nz ", Km7(45:51,1,nz)

      ! 短い時間ステップの初期値作成.
      ! 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

      write(*,*) "DVelXDt_long", minval( pyz_AccelVelXNl(1:nx,1:ny,1:nz) ), maxval( pyz_AccelVelXNl(1:nx,1:ny,1:nz) )
      write(*,*) "DVelZDt_long", minval( xyr_AccelVelZNl(1:nx,1:ny,0:nz) ), maxval( xyr_AccelVelZNl(1:nx,1:ny,0:nz) )

      ! 短い時間ステップの時間積分. オイラー法を利用.
      ! Time integration with short time step.
      !
      Euler: do tau = 1, NstepShort

        ! 速度 u の計算.
        ! Time integration horizontal velocity (u).
        !
        pyz_VelXAs = &
          & pyz_VelXNs                                          &
          &  + DelTimeShort                                     &
          &    * (                                              &
          &     - pyz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs) &
          &     + pyz_AccelVelXNl                               &
          &     )

        Km5 = - pyz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs)
        write(*,*) "DVelXDt_short", minval( Km5(1:nx,1:ny,1:nz) ), maxval( Km5(1:nx,1:ny,1:nz) )

        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_pyz( pyz_VelXAs ) ! (inout)
        call BoundaryYCyc_pyz( pyz_VelXAs ) ! (inout)
        call BoundaryZSym_pyz( pyz_VelXAs ) ! (inout)

        write(*,*) "VelX", minval( pyz_VelXAs ), maxval( pyz_VelXAs )

        ! 速度 u の計算.
        ! Time integration horizontal velocity (v).
        !
        xqz_VelYAs = &
          & xqz_VelYNs                                          &
          &  + DelTimeShort                                     &
          &    * (                                              &
          &     - xqz_GradPi(xyz_ExnerNs, pyz_VelXNs, xqz_VelYNs, xyr_VelZNs)&
          &     + xqz_AccelVelYNl                               &
          &     )

        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_xqz( xqz_VelYAs ) ! (inout)
        call BoundaryYCyc_xqz( xqz_VelYAs ) ! (inout)
        call BoundaryZSym_xqz( xqz_VelYAs ) ! (inout)

        write(*,*) "VelY", minval( xqz_VelYAs ), maxval( xqz_VelYAs )

        ! エクスナー関数の計算.
        ! Time integration exner function.
        !
        xyz_ExnerAs = xyz_Exner( &
          & xyr_AccelVelZNl,     &
          & pyz_VelXNs,          &
          & pyz_VelXAs,          &
          & xqz_VelYNs,          &
          & xqz_VelYAs,          &
          & xyr_VelZNs,          &
          & xyz_ExnerNs          &
          & )

        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_xyz( xyz_ExnerAs ) ! (inout)
        call BoundaryYCyc_xyz( xyz_ExnerAs ) ! (inout)
        call BoundaryZSym_xyz( xyz_ExnerAs ) ! (inout)

        write(*,*) "Exner", minval( xyz_ExnerAs ), maxval( xyz_ExnerAs )

        ! 速度 w の計算
        ! Time integration vertical velocity.
        !
        xyr_VelZAs =                                                   &
          & xyr_VelZNs                                                 &
          &  + DelTimeShort                                            &
          &    * (                                                     &
          &     - xyr_GradPi(xyz_ExnerAs,xyz_ExnerNs,                  &
          &                  pyz_VelXNs,xqz_VelYNs, xyr_VelZNs)        &
          &     + xyr_AccelVelZNl                                      &
          &     )

        ! 境界条件 Boundary condition
        !
        call BoundaryXCyc_xyr( xyr_VelZAs ) ! (inout)
        call BoundaryYCyc_xyr( xyr_VelZAs ) ! (inout)
        call BoundaryZAsym_xyr( xyr_VelZAs ) ! (inout)

        write(*,*) "VelZ", minval( xyr_VelZAs ), maxval( xyr_VelZAs )

        ! 短い時間ステップのループを回すための処置
        ! 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

      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

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

      write(*,*) "3D ", TimeN
      write(*,*) "VelX  ", minval( pyz_VelXAl ), maxval( pyz_VelXAl ), pyz_VelXAl(45:51,1,1)
      write(*,*) "VelY  ", minval( xqz_VelYAl ), maxval( xqz_VelYAl ), xqz_VelYAl(45:51,1,1)
      write(*,*) "VelZ  ", minval( xyr_VelZAl ), maxval( xyr_VelZAl ), xyr_VelZAl(45:51,1,1)
      write(*,*) "Exner ", minval( xyz_ExnerAl ), maxval( xyz_ExnerAl ), xyz_ExnerAl(45:51,1,1)
      write(*,*) "PTemp ", minval( xyz_PotTempAl ), maxval( xyz_PotTempAl ), xyz_PotTempAl(45:51,1,1)
      write(*,*) "Km    ", minval( xyz_KmAl ), maxval( xyz_KmAl ), xyz_KmAl(45:51,1,1)

!      ! 混合比がゼロ以下にならないための処置.
!      ! Negative values mixing ratios are corrected.
!      !
!      xyza_MixRtWork = xyza_MixRtAl
!      xyza_MixRtAl = max( - xyza_MixRtBasicZ, xyza_MixRtWork )
!!     call StoreMixRtFill2( &
!!       & (xyza_MixRtAl - xyza_MixRtWork) / DelTimeLFrog & ! (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)

!      write(*,*) "Exner", maxval( xyz_ExnerNl ), minval( xyz_ExnerNl )

      call HistoryAutoPut(TimeN, 'PotTemp', xyz_PotTempNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Exner', xyz_ExnerNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelX',  xyz_VelXNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelY',  xyz_VelYNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelZ',  xyz_VelZNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Km',    xyz_KmNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Kh',    xyz_KhNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax))
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s)), xyza_MixRtNl(FileXMin:FileXMax, FileYMin:FileYMax, FileZMin:FileZMax, 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( 'PotTemp', xyz_PotTempNl, rstat)
      call HistoryPut( 'Km',      xyz_KmNl,      rstat)
      call HistoryPut( 'Kh',      xyz_KhNl,      rstat)
      call HistoryPut( 'MixRt',   xyza_MixRtNl,  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( 'PotTemp', xyz_PotTempAl, rstat)
      call HistoryPut( 'Km',      xyz_KmAl,      rstat)
      call HistoryPut( 'Kh',      xyz_KhAl,      rstat)
      call HistoryPut( 'MixRt',   xyza_MixRtAl, rstat)

      ! Show CPU time
      !
      call ClocksetPredict

    end if

      ! 長い時間ステップのループを回すための処置.
      ! Renew prognostic variables for next long time step integration.
      !
      xyz_PotTempBl = xyz_PotTempNl
      xyza_MixRtBl  = xyza_MixRtNl
      xyz_ExnerBl   = xyz_ExnerNl
      pyz_VelXBl    = pyz_VelXNl
      xqz_VelYBl    = xqz_VelYNl
      xyr_VelZBl    = xyr_VelZNl
      xyz_KmBl      = xyz_KmNl
      xyz_KhBl      = xyz_KhNl

      xyz_PotTempNl = xyz_PotTempAl
      xyza_MixRtNl  = xyza_MixRtAl
      xyz_ExnerNl   = xyz_ExnerAl
      pyz_VelXNl    = pyz_VelXAl
      xqz_VelYNl    = xqz_VelYAl
      xyr_VelZNl    = xyr_VelZAl
      xyz_KmNl      = xyz_KmAl
      xyz_KhNl      = xyz_KhAl

      ! 時刻の進行
      ! Progress time
      !
      call TimesetProgress
      t = t + 1
  end do

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

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

  !----------------------------------------------------
  ! CPU Time
  !
  call ClocksetLoopStop
  call ClocksetClose

contains

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

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

    real(8) :: tfil2

    tfil2 = 1.0d0 - 2.0d0 * tfil

    pyz_VelXNl  =  pyz_VelXNl   + tfil * (pyz_VelXBl  - 2.0d0 * pyz_VelXNl  + pyz_VelXAl)
    xqz_VelYNl  =  xqz_VelYNl   + tfil * (xqz_VelYBl  - 2.0d0 * xqz_VelYNl  + xqz_VelYAl)
    xyr_VelZNl  =  xyr_VelZNl   + tfil * (xyr_VelZBl  - 2.0d0 * xyr_VelZNl  + xyr_VelZAl)
    xyz_PotTempNl = xyz_PotTempNl  + tfil * (xyz_PotTempBl - 2.0d0 * xyz_PotTempNl + xyz_PotTempAl)
    xyz_ExnerNl =  xyz_ExnerNl  + tfil * (xyz_ExnerBl - 2.0d0 * xyz_ExnerNl + xyz_ExnerAl)
    xyz_KmNl    =  xyz_KmNl     + tfil * (xyz_KmBl    - 2.0d0 * xyz_KmNl    + xyz_KmAl)
    xyz_KhNl    =  xyz_KhNl     + tfil * (xyz_KhBl    - 2.0d0 * xyz_KhNl    + xyz_KhAl)
    xyza_MixRtNl = xyza_MixRtNl + tfil * (xyza_MixRtBl- 2.0d0 * xyza_MixRtNl + xyza_MixRtAl)

!    pz_VelXNl  =  tfil * ( pz_VelXBl + pz_VelXAl )       + tfil2 * pz_VelXNl
!    xr_VelZNl  =  tfil * ( xr_VelZBl + xr_VelZAl )       + tfil2 * xr_VelZNl
!    xyz_PotTempNl= tfil * ( xyz_PotTempBl + xyz_PotTempAl ) + tfil2 * xyz_PotTempNl
!    xyz_ExnerNl =  tfil * ( xyz_ExnerBl + xyz_ExnerAl )     + tfil2 * xyz_ExnerNl
!    xyz_KmNl    =  tfil * ( xyz_KmBl + xyz_KmAl )           + tfil2 * xyz_KmNl
!    xyz_KhNl    =  tfil * ( xyz_KhBl + xyz_KhAl )           + tfil2 * xyz_KhNl
!    xza_MixRtNl = tfil * ( xza_MixRtBl + xza_MixRtAl )   + tfil2 * xza_MixRtNl

  end subroutine AsselinTimeFilter

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

    ! NAMELIST ファイル名の読み込み
    ! Loading NAMELIST file.
    !
    implicit none

    ! 変数の設定
    !
    character(STRING) :: cfgfile  ! NAMELIST ファイル名

    ! 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 constatns.
    !
    call chemdata_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

    ! 湿潤過程共有変数の初期化
    ! Initialization of common variables for moist process.
    !
!    call moistset_init( )

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

    ! マージンの設定の初期化
    ! Initialization of margin
    !
!    call SetMargin_Init

    ! 内部変数の初期化
    ! Initialization of internal variables.
    !
    call ArareAlloc

    ! 初期値の代入
    ! * 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 ファイルから取得する.
    !
!    write(*,*) "start"
    call ReStartFileio_BasicZ_Get
!    write(*,*) "end"
    call ReStartFileio_Var_Get(     &
      & pyz_VelXBl,  pyz_VelXNl,    & ! (out)
      & xqz_VelYBl,  xqz_VelYNl,    & ! (out)
      & xyr_VelZBl,  xyr_VelZNl,    & ! (out)
      & xyz_PotTempBl, xyz_PotTempNl,   & ! (out)
      & xyz_ExnerBl, xyz_ExnerNl,   & ! (out)
      & xyza_MixRtBl,xyza_MixRtNl,  & ! (out)
      & xyz_KmBl,    xyz_KmNl,      & ! (out)
      & xyz_KhBl,    xyz_KhNl       & ! (out)
      & )

    ! 化学計算ルーチンの初期化
    ! Initialization of chemical routines.
    !
    call chemcalc_init

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

    ! 数値拡散項の初期化
    ! Initialization of numerical diffusion term.
    !
    call NumDiffusion_Init( &
      & cfgfile             & ! (in)
      & )

    ! 乱流拡散項の初期化
    ! Initialization of turbulent diffusion term.
    !
    call Turbulence_Init( )

    ! 暖かい雨のパラメタリゼーションの初期化
    ! Initialization of warmrain parameterization.
    !
    call WarmRainPrm_Init( &
      & cfgfile            & ! (in)
      & )

    ! 負の湿潤量の補填計算の初期化
    ! Initialization of negative moist value correction.
    !
    call FillNegative_Init( &
      & xyza_MixRtBasicZ,    & ! (in)
      & xyz_DensBasicZ       & ! (in)
      & )

    ! 放射強制の初期化
    !  Initialization of radiative forcing.
    !
    call Radiation_Init( &
      & cfgfile          & ! (in)
      & )

!    write(*,*) "21"

    ! 湿潤気塊の浮力計算の初期化
    ! Initialization of moist buoyancy calculation.
    !
    call MoistBuoy_Init( )

!    write(*,*) "22"

    ! 圧力計算用係数行列の初期化
    ! Initialization of coefficient matrix for exner function calculation.
    !
    call xyz_Exner_Init( )

!    write(*,*) "23"
!    write(*,*) "24"

    !-------------------------------------------------------------
    ! 基本場のファイル出力
    !
!    call HistoryPut( 'DensBasicZ',     xyz_DensBasicZ     , rstat)
!    call HistoryPut( 'ExnerBasicZ',    xyz_ExnerBasicZ    , rstat)
!    call HistoryPut( 'PotTempBasicZ',  xyz_PotTempBasicZ  , rstat)
!    call HistoryPut( 'VelSoundBasicZ', xyz_VelSoundBasicZ , rstat)
!    call HistoryPut( 'TempBasicZ',     xyz_TempBasicZ     , rstat)
!    call HistoryPut( 'PressBasicZ',    xyz_PressBasicZ    , rstat)
!    call HistoryPut( 'MixRtBasicZ',    xyza_MixRtBasicZ   , rstat)
!    call HistoryPut( 'EffMolWtBasicZ', xyz_EffMolWtBasicZ , rstat)

!    write(*,*) "24"

    ! 音波に対する CFL 条件のチェック
    ! CFL condtion check for sound wave.
    !
    call CFLCheckTimeShort( &
      & xyz_VelSoundBasicZ   & ! (in)
      & )

!    write(*,*) "25"

    ! CPU time
    !
    call ClocksetPreStop

!    ! 保存量のモニタリング
!    !
!    call EnergyMonit_init

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

!    write(*,*) "26"

  end subroutine MainInit

  !-----------------------------------------------------------------------
  subroutine ArareAlloc
    !
    !初期化として, 配列を定義し, 値としてゼロを代入する.
    !

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate(                                                       &
      & xyz_VelXNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_VelYNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_VelZNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & pyz_VelXBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & pyz_VelXAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xqz_VelYBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_VelYAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyr_VelZBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_VelZAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_ExnerBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerNs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_ExnerAs(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_PotTempWork(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_PotTempAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_KmBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KmNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KmAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyz_KhBl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KhNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyz_KhAl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyza_MixRtWork(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum), &
      & xyza_MixRtBl  (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum), &
      & xyza_MixRtNl  (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum), &
      & xyza_MixRtAl  (DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum), &
!
      & pyz_AccelVelXNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xqz_AccelVelYNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
      & xyr_AccelVelZNl(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax),&
!
      & xyza_DelMixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum) )

    pyz_VelXNs  = 0.0d0;    pyz_VelXAs = 0.0d0
    xqz_VelYNs  = 0.0d0;    xqz_VelYAs = 0.0d0
    xyr_VelZNs  = 0.0d0;    xyr_VelZAs = 0.0d0
    xyz_ExnerNs = 0.0d0;    xyz_ExnerAs = 0.0d0

    pyz_VelXBl  = 0.0d0;    pyz_VelXNl  = 0.0d0;    pyz_VelXAl  = 0.0d0
    xqz_VelYBl  = 0.0d0;    xqz_VelYNl  = 0.0d0;    xqz_VelYAl  = 0.0d0
    xyr_VelZBl  = 0.0d0;    xyr_VelZNl  = 0.0d0;    xyr_VelZAl  = 0.0d0
    xyz_ExnerBl = 0.0d0;    xyz_ExnerNl = 0.0d0;    xyz_ExnerAl = 0.0d0
    xyz_KmBl    = 0.0d0;    xyz_KmNl    = 0.0d0;    xyz_KmAl    = 0.0d0
    xyz_KhBl    = 0.0d0;    xyz_KhNl    = 0.0d0;    xyz_KhAl    = 0.0d0
    xyz_PotTempBl = 0.0d0;  xyz_PotTempNl = 0.0d0;  xyz_PotTempAl = 0.0d0
    xyza_MixRtBl = 0.0d0;   xyza_MixRtNl = 0.0d0;   xyza_MixRtAl = 0.0d0

    xyz_PotTempWork = 0.0d0
    xyza_MixRtWork = 0.0d0

    pyz_AccelVelXNl = 0.0d0
    xqz_AccelVelYNl = 0.0d0
    xyr_AccelVelZNl = 0.0d0

    xyza_DelMixRt = 0.0d0

    DelTimeLFrog = 0.0d0

    allocate(                                            &
         & Km1(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km2(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km3(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km4(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km5(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km6(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax),        &
         & Km7(DimXMin:DimXMax, DimYMin:DimYMax, DimZMin:DimZMax)         &
         &)
    Km1 = 0.0d0
    Km2 = 0.0d0
    Km3 = 0.0d0
    Km4 = 0.0d0
    Km5 = 0.0d0
    Km6 = 0.0d0
    Km7 = 0.0d0

  end subroutine ArareAlloc
!-----------------------------------------------------------------------
end program arare