!= deepconv/arare 木星大気対流計算用主プログラム
!
!= deepconv/arare main program for Jovian atmospheric convection
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: arare_jupiter.f90,v 1.7 2010-08-11 07:23:05 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. 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, DimZMin, DimZMax, &
    &                 SpcNum, FileXMin, FileXMax, 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, only: xza_MixRtBasicZ, xz_DensBasicZ, xz_EffMolWtBasicZ, &
    &                 xz_PotTempBasicZ, xz_TempBasicZ, xz_PressBasicZ, &
    &                 xz_VelSoundBasicZ, xz_ExnerBasicZ
  use ChemCalc, only: ChemCalc_init
  use namelist_util, only: NmlutilInit, namelist_filename

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

  ! 力学過程計算用関数モジュール
  ! Dynamical processes module
  !
  use DynFunc,       only: xz_AdvScalar, xz_AdvKm, xza_AdvScalar, &
    &                      pz_AdvVelX, xr_Buoy, xr_AdvVelZ, pz_GradPi
  use DynImpFunc,    only: xz_Exner_init, xz_Exner, xr_GradPi

  ! 乱流拡散計算用モジュール
  ! Turbulent diffusion module
  !
  use Turbulence,    only: Turbulence_Init, xz_TurbScalar, xza_TurbScalar,  &
    &                      pz_TurbVelX, xr_TurbVelZ, xz_ShearKm, xz_DispKm, &
    &                      xz_DispHeat, xz_TurbKm, xz_BuoyKm
  
  ! 境界からのフラックス計算用モジュール
  ! Surface flux module
  !
  use HeatFlux,      only: xz_HeatFluxDiff, xza_MixRtFluxDiff

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

  ! 湿潤過程計算用モジュール
  ! Moist processes modules
  !
!  use moistset,      only: moistset_init
  use MoistAdjust,   only: MoistAdjustSvapPress, MoistAdjustNH4SH
  use WarmRainPrm,   only: WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, &
    &                      xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH,         &
    &                      xza_Cloud2Rain, xza_FallRain
  use MoistBuoyancy, only: MoistBuoy_Init, xz_BuoyMoistKm, xr_BuoyMolWt, &
    &                      xr_BuoyDrag
  use fillnegative,  only: FillNegative_init, xza_FillNegative_xza

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

  ! 数値拡散/摩擦計算用モジュール
  ! Numerical diffussion /dumping module
  !
  use NumDiffusion4th,only: NumDiffusion_Init, xz_NumDiffScalar, &
    &                       xz_NumDiffKm, xza_NumDiffScalar,     &
    &                       pz_NumDiffVelX, xr_NumDiffVelZ
  use damping,        only: damping_init, xz_Sponge, pz_Sponge, &
    &                       xr_Sponge

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

  ! 下請けモジュール
  ! Utility modules
  !
  use cflcheck,      only: CFLCheckTimeShort, CFLCheckTimeLongVelX, &
    &                      CFLCheckTimeLongVelZ
  use boundary,      only: BoundaryXCyc_xz, BoundaryZSym_xz, &
    &                      BoundaryXCyc_xza, BoundaryZSym_xza, BoundaryXCyc_pz, &
    &                      BoundaryZSym_pz, BoundaryXCyc_xr, BoundaryZAntiSym_xr
  use setmargin, only: SetMargin_init, SetMargin_xyzf, SetMargin_xyz, &
    &                  SetMargin_pyz, SetMargin_xqz, SetMargin_xyr, &
    &                  SetMargin_az, SetMargin_ar, SetMargin_xza

  use average, only: xz_avr_pz, xz_avr_xr


  ! 変数宣言 Variables definition 
  !
  implicit none

  ! 内部変数
  ! Internal variables
  !
!  character(80) :: cfgfile  ! NAMELIST ファイル名 NAMELIST file name
  real(8), allocatable :: xz_VelXNl(:,:) 
                            ! $ u(t-\Delta t) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXBl(:,:) 
                            ! $ u(t-\Delta t) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXNl(:,:) 
                            ! $ u(t) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXAl(:,:) 
                            ! $ u (t+\Delta t)a $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXNs(:,:) 
                            ! $ u (\tau) $ 水平風 Horizontal wind
  real(8), allocatable :: pz_VelXAs(:,:) 
                            ! $ u (\tau +\Delta \tau) 水平風 Horizontal wind
  real(8), allocatable :: xz_VelZNl(:,:) 
                            ! $ w (t-\Delta t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZBl(:,:) 
                            ! $ w (t-\Delta t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZNl(:,:) 
                            ! $ w (t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZAl(:,:) 
                            ! $ w (t+\Delta t) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZNs(:,:) 
                            ! $ w (\tau) 鉛直風 Vertical wind
  real(8), allocatable :: xr_VelZAs(:,:) 
                            ! $ w (\tau +\Delta \tau) $ 鉛直風 Vertical wind
  real(8), allocatable :: xz_ExnerBl(:,:)
                            ! $ \pi (t-\Delta t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerNl(:,:)
                            ! $ \pi (t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerAl(:,:)
                            ! $ \pi (t+\Delta t) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerNs(:,:)
                            ! $ \pi (\tau) $ 圧力関数 Exner function
  real(8), allocatable :: xz_ExnerAs(:,:)
                            ! $ \pi (\tau +\Delta \tau) $ 圧力関数 Exner function
  real(8), allocatable :: xz_PotTempBl(:,:) 
                            ! $ \theta (t-\Delta t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempNl(:,:) 
                            ! $ \theta (t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempAl(:,:) 
                            ! $ \theta (t+\Delta t) $ 温位 Potential temp.
  real(8), allocatable :: xz_PotTempWork(:,:)
                            ! 温位の作業配列 Work array for potential temp.
  real(8), allocatable :: xz_PotTempOrg(:,:)
                            ! 温位の作業配列 Work array for potential temp.
  real(8), allocatable :: xz_KmBl(:,:)
                            ! $ Km (t-\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KmNl(:,:)         
                            ! $ K_m (t) 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KmAl(:,:)
                            ! $ K_m (t+\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhBl(:,:)
                            ! $ K_h (t-\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhNl(:,:)
                            ! $ K_h (t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff. 
  real(8), allocatable :: xz_KhAl(:,:)
                            ! $ K_h (t+\Delta t) $ 乱流拡散係数 
                            ! Turbulent diffusion coeff.
  real(8), allocatable :: xza_MixRtBl(:,:,:)
                            ! $ q (t-\Delta t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtNl(:,:,:) 
                            ! $ q (t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtAl(:,:,:)
                            ! $ q (t+\Delta t) $ 湿潤量の混合比 
                            ! Mixing ratio of moist variables.
  real(8), allocatable :: xza_MixRtWork(:,:,:) 
                            ! 湿潤量の作業配列
                            ! Work array for mixing ratio.
  real(8), allocatable :: xza_MixRtOrg(:,:,:) 
                            ! 湿潤量の作業配列
                            ! Work array for mixing ratio.
  real(8), allocatable :: pz_AccelVelXNl(:,:)  
                            ! 気圧傾度力を除いた $u$ の変化率
                            ! Tendency of $u$ except for pressure gradient term
  real(8), allocatable :: xr_AccelVelZNl(:,:)  
                            ! 気圧傾度力を除いた $w$ の変化率
                            ! Tendency of $w$ except for pressure gradient term
  real(8), allocatable :: xza_DelMixRt(:,:,:)  
                            ! 湿潤量の混合比 $ q $ の増分
                            ! Mixing ratio variation.
  real(8) :: DelTimeLFrog
                            ! リープフロッグスキーム用時間格子間隔
                            ! Time interval for Leap-frog scheme

  real(8), allocatable :: Km1(:,:), Km2(:,:), Km3(:,:), Km4(:,:), Km5(:,:), Km6(:,:), Km7(:,:)
  real(8), allocatable :: Tnd1(:,:,:), Tnd2(:,:,:), Tnd3(:,:,:), Tnd4(:,:,:), Tnd5(:,:,:)

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

  logical :: CloudFlag = .true.
  logical :: TurbFlag  = .true.
  logical :: BuoyMoistFlag = .true.
  
  !------------------------------------------
  ! 初期化手続き Initialize procedure 
  !
  call MainInit

  !------------------------------------------
  ! 時間積分 time integration 
  !
!!  do t1 = 1, NstepLFrog / NstepDisp
!!    do t2 = 1, NstepDisp
 
!  do t = 1, NstepLFrog + 1
  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.
      !
      if (TurbFlag) then 
!      xz_KmAl =                                                       &
!        & xz_KmBl                                                     &
!        & + DelTimeLFrog                                           &
!        &   * (                                                       &
!        &     + xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)               &
!        &     + xz_NumDiffKm(xz_KmBl)                                 &
!        &     + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) &
!        &     + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)             &
!        &     + xz_TurbKm(xz_KmBl)                                    &
!        &     + xz_DispKm(xz_KmBl)                                    &
!        &    )

      Km1 = xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl)
      Km2 = xz_NumDiffKm(xz_KmBl)
      Km3 = xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl)
      Km4 = xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)
      Km5 = xz_TurbKm(xz_KmBl)
      Km6 = xz_DispKm(xz_KmBl)
      Km7 = xz_BuoyKm(xz_PotTempBl)

      xz_KmAl =                                          &
        & xz_KmBl                                        &
        & + DelTimeLFrog                                 &
        &   * (                                          &
        &      Km1 + Km2 + Km3 + Km4 + Km5 + Km6         &
        &    )

      call HistoryAutoPut(TimeN,'DKmDtAdv',  Km1(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtDiff', Km2(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtBuoy', Km3(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtShear',Km4(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtTurb', Km5(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtDisp', Km6(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN,'DKmDtBuoyT',Km7(FileXMin:FileXMax,FileZMin:FileZMax))

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

      ! 境界条件 Boundary condition
      call BoundaryXCyc_xz( xz_KmAl )  ! (inout)
      call BoundaryZSym_xz( xz_KmAl )  ! (inout)
!      call SetMargin_aaz( xz_KmAl )  ! (inout)

      else

        xz_KmBl = 0.0d0
        xz_KmNl = 0.0d0
        xz_KmAl = 0.0d0

      end if

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

      ! 温位の移流拡散の計算 
      ! Advection and diffusion of potential temperature.
      !
!      xz_PotTempAl =                                                  &
!        &   xz_PotTempBl                                              &
!        & + DelTimeLFrog                                              &
!        &   * (                                                       &
!        &     + xz_AdvScalar( xz_PotTempNl,     pz_VelXNl, xr_VelZNl) &
!        &     + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl) &
!        &     + xz_NumDiffScalar(xz_PotTempBl)                        &
!        &     + xz_RadHeatConst( xz_ExnerBl )                         &
!        &     + xz_Sponge( xz_PotTempBl )                             &
!        &     + xz_TurbScalar(xz_PotTempBl,     xz_KhBl)              &
!        &     + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)              &
!        &     + xz_DispHeat( xz_KmBl )                                &
!        &     + xz_HeatFluxDiff( xz_PotTempNl )                       &
!        &      )

      Km1 = + xz_AdvScalar( xz_PotTempNl,     pz_VelXNl, xr_VelZNl) &
        &   + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl)
      Km2 = + xz_NumDiffScalar(xz_PotTempBl)
      Km3 = + xz_RadHeatConst( xz_ExnerBl )
      Km4 = + xz_Sponge( xz_PotTempBl )
      Km5 = + xz_TurbScalar(xz_PotTempBl, xz_KhBl) &
        &   + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)
      Km6 = + xz_DispHeat( xz_KmBl )
      Km7 = xz_HeatFluxDiff( xz_PotTempNl )

      xz_PotTempAl = xz_PotTempBl + DelTimeLFrog * ( Km1 + Km2 + Km3 + Km4 + Km5 + Km6 + Km7 )
      
      call HistoryAutoPut(TimeN, 'PotTempAdv',  Km1(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'PotTempDiff', Km2(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'PotTempRad',  Km3(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'PotTempSpng', Km4(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'PotTempTurb', Km5(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'PotTempDisp', Km6(FileXMin:FileXMax,FileZMin:FileZMax))        
      call HistoryAutoPut(TimeN, 'PotTempFlux', Km7(FileXMin:FileXMax,FileZMin:FileZMax))        

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xz( xz_PotTempAl ) ! (inout)
      call BoundaryZSym_xz( xz_PotTempAl ) ! (inout)


      ! 凝縮成分混合比の移流拡散 
      ! Advection and diffusion of vapor, cloud and rain mixing ratios.
      !
!        xza_MixRtAl =                                                 & 
!          &   xza_MixRtBl                                             &
!          & + DelTimeLFrog                                            &
!          &   * (                                                     &
!          &    + xza_AdvScalar(xza_MixRtNl,     pz_VelXNl, xr_VelZNl) &
!          &    + xza_AdvScalar(xza_MixRtBasicZ, pz_VelXNl, xr_VelZNl) & 
!          &    + xza_NumDiffScalar(xza_MixRtBl)                       &
!          &    + xza_FallRain(xza_MixRtBl)                            &
!          &    + xza_TurbScalar(xza_MixRtBl,    xz_KhBl)              &
!          &    + xza_TurbScalar(xza_MixRtBasicZ,xz_KhBl)              &
!          &    + xza_MixRtFluxDiff(xza_MixRtNl)                       &
!          &   )
        
        Tnd1 = + xza_AdvScalar(xza_MixRtNl,     pz_VelXNl, xr_VelZNl) &
          &    + xza_AdvScalar(xza_MixRtBasicZ, pz_VelXNl, xr_VelZNl)
        Tnd2 = + xza_NumDiffScalar(xza_MixRtBl)
        Tnd3 = + xza_FallRain(xza_MixRtBl)
        Tnd4 = + xza_TurbScalar(xza_MixRtBl,    xz_KhBl)           &
          &    + xza_TurbScalar(xza_MixRtBasicZ,xz_KhBl)
        Tnd5 = + xza_MixRtFluxDiff(xza_MixRtNl)

        xza_MixRtAl = xza_MixRtBl + DelTimeLFrog * ( Tnd1 + Tnd2 + Tnd3 + Tnd4 + Tnd5 )

        do s = 1, SpcNum
          call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtAdv',  &
            & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s))
          call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtDiff', &
            & Tnd2(FileXMin:FileXMax, FileZMin:FileZMax, s))
          call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtFall', &
            & Tnd3(FileXMin:FileXMax, FileZMin:FileZMax, s))
          call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtTurb',  &
            & Tnd4(FileXMin:FileXMax, FileZMin:FileZMax, s))
          call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtFlux',  &
            & Tnd5(FileXMin:FileXMax, FileZMin:FileZMax, s))
        end do

      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)
!      call SetMargin_xyzf( xza_MixRtAl )  ! (inout)

      ! 移流によって負になった部分を埋める
      ! Negative values due to advection are corrected.
      !
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = xza_FillNegative_xza( xza_MixRtWork ) 
      
      ! 埋めた/削った量を保管
      ! Correction value is stored.
      !
!      call StoreMixRtFill1( &
!        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )    
      Tnd1 = (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtFill1',  &
          & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do      

      ! 埋め切れなかった部分をゼロにする. Fill2 に保管
      ! Negative values mixing ratios are corrected.
      xza_MixRtWork = xza_MixRtAl
      xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
!      call StoreMixRtFill2( &
!        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )
      Tnd1 = (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog 
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtFill2',  &
          & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do
       
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)
!      call SetMargin_xyzf( xza_MixRtAl )  ! (inout)


      !!!------------------------------------------------------
      !!! 凝結過程
      !!!
      if (CloudFlag) then 
        call ArareCloudPhys
      end if
      
      ! 速度の移流拡散.
      ! Advection and diffusion of velocity components.
      !
!      pz_AccelVelXNl =                                    &
!        & + pz_AdvVelX(pz_VelXNl, xr_VelZNl)              &
!        & + pz_TurbVelX(xz_KmBl, pz_VelXBl, xr_VelZBl)   &
!        & + pz_Sponge(pz_VelXBl)                          &
!        & + pz_NumDiffVelX(pz_VelXBl)

      Km1 = pz_AdvVelX(pz_VelXNl, xr_VelZNl)
      Km2 = pz_Sponge(pz_VelXBl)
      Km3 = pz_NumDiffVelX(pz_VelXBl)
      Km4 = pz_TurbVelX(xz_KmBl, pz_VelXBl, xr_VelZBl)

      pz_AccelVelXNl = Km1 + Km2 + Km3 + Km4

      call HistoryAutoPut(TimeN, 'DVelXDtAdv',  Km1(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelXDtSpng', Km2(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelXDtDiff', Km3(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelXDtTurb', Km4(FileXMin:FileXMax,FileZMin:FileZMax))


!      xr_AccelVelZNl =                                    &
!        & + xr_AdvVelZ(xr_VelZNl, pz_VelXNl)              &
!        & + xr_Buoy(xz_PotTempNl)                         &    
!        & + xr_Sponge(xr_VelZBl)                          &
!        & + xr_NumDiffVelZ(xr_VelZBl)                     &
!        & + xr_TurbVelZ(xz_KmBl, pz_VelXBl, xr_VelZBl)    &

      Km1 = xr_AdvVelZ(xr_VelZNl, pz_VelXNl)
      Km2 = xr_Buoy(xz_PotTempNl)
      Km3 = xr_Sponge(xr_VelZBl) 
      Km4 = xr_NumDiffVelZ(xr_VelZBl)
      Km5 = xr_TurbVelZ(xz_KmBl, pz_VelXBl, xr_VelZBl)
      Km6 = xr_BuoyMolWt(xza_MixRtNl)
      Km7 = xr_BuoyDrag(xza_MixRtNl)

      if ( BuoyMoistFlag ) then 
        xr_AccelVelZNl = Km1 + Km2 + Km3 + Km4 + Km5 + Km6 + Km7
      else
        xr_AccelVelZNl = Km1 + Km2 + Km3 + Km4 + Km5
      end if

      call HistoryAutoPut(TimeN, 'DVelZDtAdv',  Km1(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtBuoyT',Km2(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtSpng', Km3(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtDiff', Km4(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtTurb', Km5(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtBuoyM',Km6(FileXMin:FileXMax,FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'DVelZDtBuoyD',Km7(FileXMin:FileXMax,FileZMin:FileZMax))

      ! 短い時間ステップの初期値作成.
      ! Initial values set up for time integration with short time step.
      !
      pz_VelXNs  = pz_VelXBl
      xr_VelZNs  = xr_VelZBl
      xz_ExnerNs = xz_ExnerBl
      
      ! 短い時間ステップの時間積分. オイラー法を利用.o
      ! Time integration with short time step.
      !
      Eular: do tau = 1, NstepShort
        
        ! 速度 u の計算.
        ! Time integration horizontal velocity.
        !
!        pz_VelXAs = &
!          & pz_VelXNs                                          &
!          &  + DelTimeShort                                    &
!          &    * (                                             &
!          &     - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)  &
!          &     + pz_AccelVelXNl                               &
!          &     )

        Km1 = - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)

        pz_VelXAs = pz_VelXNs + DelTimeShort * ( Km1 + pz_AccelVelXNl )

        call HistoryAutoPut(TimeN, 'DVelXDtPGrad',Km1(FileXMin:FileXMax,FileZMin:FileZMax))

        ! 境界条件 Boundary conditions.
        ! 
        call BoundaryXCyc_pz( pz_VelXAs ) ! (inout)
        call BoundaryZSym_pz( pz_VelXAs ) ! (inout)
!        call SetMargin_aaz( pz_VelXAs )  ! (inout)
       

        ! エクスナー関数の計算.
        ! Time integration exner function.
        !
        xz_ExnerAs = xz_Exner( &
          & xr_AccelVelZNl,    &
          & pz_VelXNs,         &
          & pz_VelXAs,         &
          & xr_VelZNs,         &
          & xz_ExnerNs)
        
        ! 境界条件  Boundary condition
        !
        call BoundaryXCyc_xz( xz_ExnerAs ) ! (inout)
        call BoundaryZSym_xz( xz_ExnerAs ) ! (inout)
!        call SetMargin_aaz( xz_ExnerAs )  ! (inout)


        ! 速度 w の計算
        ! Time integration vertical velocity.
        !
!        xr_VelZAs =                                                    &
!          & xr_VelZNs                                                  &
!          &  + DelTimeShort                                            &
!          &    * (                                                     &
!          &     - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
!          &     + xr_AccelVelZNl                                       &
!          &     )

        Km1 = - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs)

        xr_VelZAs = xr_VelZNs + DelTimeShort * ( Km1 + xr_AccelVelZNl )

        call HistoryAutoPut(TimeN, 'DVelZDtPGrad',Km1(FileXMin:FileXMax,FileZMin:FileZMax))        
        ! 境界条件 ; Boundary condition
        !
        call BoundaryXCyc_xr( xr_VelZAs )     ! (inout)
        call BoundaryZAntiSym_xr( xr_VelZAs ) ! (inout)
!        call SetMargin_aar( xr_VelZAs )  ! (inout)

        
        ! 短い時間ステップのループを回すための処置
        ! Renew prognostic variables for next short time step integration.
        !
        xz_ExnerNs  = xz_ExnerAs
        pz_VelXNs   = pz_VelXAs
        xr_VelZNs   = xr_VelZAs
        
      end do Eular
      
      ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
      ! Renew prognostic variables for next long time step integration.
      !
      xz_ExnerAl  = xz_ExnerAs
      pz_VelXAl   = pz_VelXAs
      xr_VelZAl   = xr_VelZAs

!DEBUG 
      call HistoryAutoPut(TimeN, 'PotTempBl', xz_PotTempBl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'ExnerBl', xz_ExnerBl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelXBl',  pz_VelXBl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelZBl',  xr_VelZBl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'KmBl',    xz_KmBl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'KhBl',    xz_KhBl(FileXMin:FileXMax, FileZMin:FileZMax))
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'Bl', xza_MixRtBl(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do

      call HistoryAutoPut(TimeN, 'PotTempNl', xz_PotTempNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'ExnerNl', xz_ExnerNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelXNl',  pz_VelXNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelZNl',  xr_VelZNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'KmNl',    xz_KmNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'KhNl',    xz_KhNl(FileXMin:FileXMax, FileZMin:FileZMax))
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'Nl', xza_MixRtNl(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do


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

      ! ヒストリファイル出力. 時間フィルタをかけた後の数値を出力. 
      ! Out put to history file.
      !
      xz_VelXNl = xz_avr_pz(pz_VelXNl)
      xz_VelZNl = xz_avr_xr(xr_VelZNl)
      
      call HistoryAutoPut(TimeN, 'PotTemp', xz_PotTempNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Exner', xz_ExnerNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelX',  xz_VelXNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'VelZ',  xz_VelZNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Km',    xz_KmNl(FileXMin:FileXMax, FileZMin:FileZMax))
      call HistoryAutoPut(TimeN, 'Kh',    xz_KhNl(FileXMin:FileXMax, FileZMin:FileZMax))
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s)), xza_MixRtNl(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do


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

      ! ファイル出力
      !
      call HistoryPut( 't',       TimeN,        rstat)
      call HistoryPut( 'VelX',    pz_VelXNl,    rstat)
      call HistoryPut( 'VelZ',    xr_VelZNl,    rstat)
      call HistoryPut( 'Exner',   xz_ExnerNl,   rstat)
      call HistoryPut( 'PotTemp', xz_PotTempNl, rstat)
      call HistoryPut( 'Km',      xz_KmNl,      rstat)
      call HistoryPut( 'Kh',      xz_KhNl,      rstat)
      call HistoryPut( 'MixRt',   xza_MixRtNl,  rstat)    

      call HistoryPut( 't',       TimeA,        rstat)
      call HistoryPut( 'VelX',    pz_VelXAl,    rstat)
      call HistoryPut( 'VelZ',    xr_VelZAl,    rstat)
      call HistoryPut( 'Exner',   xz_ExnerAl,   rstat)
      call HistoryPut( 'PotTemp', xz_PotTempAl, rstat)
      call HistoryPut( 'Km',      xz_KmAl,      rstat)
      call HistoryPut( 'Kh',      xz_KhAl,      rstat)
      call HistoryPut( 'MixRt',   xza_MixRtAl,  rstat) 

      ! Show CPU time 
      ! 
      call ClocksetPredict
      
    end if

      ! 長い時間ステップのループを回すための処置.
      ! Renew prognostic variables for next long time step integration.
      !
      xz_PotTempBl = xz_PotTempNl
      xza_MixRtBl  = xza_MixRtNl
      xz_ExnerBl   = xz_ExnerNl
      pz_VelXBl    = pz_VelXNl
      xr_VelZBl    = xr_VelZNl
      xz_KmBl      = xz_KmNl
      xz_KhBl      = xz_KhNl
      
      xz_PotTempNl = xz_PotTempAl
      xza_MixRtNl  = xza_MixRtAl
      xz_ExnerNl   = xz_ExnerAl
      pz_VelXNl    = pz_VelXAl
      xr_VelZNl    = xr_VelZAl
      xz_KmNl      = xz_KmAl
      xz_KhNl      = xz_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 ArareCloudPhys

      ! 値の保存. 
      xza_MixRtOrg  = xza_MixRtAl
      xz_PotTempOrg = xz_PotTempAl

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

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xza_MixRtWork = xza_MixRtAl
      
      ! 雨への変化量を計算
      ! Conversion values are calculated.
      !
      xza_MixRtAl = &
        & xza_MixRtWork + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog )

      ! 雲から雨への変換量を保管
      ! Conversion values are sotred.
      !
!      call StoreMixRtCond( &
!        & (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog & ! (in)
!        & )
      Tnd1 = (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog 
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtC2R',  &
          & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do
      
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)
!      call SetMargin_xyzf( xza_MixRtAl )  ! (inout) 
     
      ! 暖かい雨のパラメタリゼーション.
      ! * 蒸気<-->雨 の変換を行う
      !
      ! Warm rain parameterization.
      ! * Conversion from rain to vapor.

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork = xza_MixRtAl
      
      ! 雨から蒸気への混合比変化を求める
      ! * 温位の計算において, 混合比変化が必要となるため, 
      !   混合比変化を 1 つの配列として用意する.
      ! Conversion values are calculated.
      !
      xza_DelMixRt = 0.0d0
      xza_DelMixRt =                                                   &
        & (                                                            &
        &   + xza_Rain2Gas(                                            &
        &       xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog &
        &       )                                                      &
        &   + xza_Rain2GasNH4SH(                                       &
        &       xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog &
        &       )                                                      &
        &  )    

      ! 温位の計算. 雨から蒸気への変換に伴う潜熱・反応熱を追加.
      !
      xz_PotTempAl =                                                    &   
        & xz_PotTempWork                                                &
        & + (                                                           &
        &     xz_Rain2GasHeat( xz_PotTempAl, xz_ExnerNl, xza_DelMixRt ) & 
        &   + xz_Rain2GasHeatNH4SH( xz_ExnerNl, xza_DelMixRt )          &
        &    )
      
      ! 混合比の計算. 雨から蒸気への変換分を追加
      !
      xza_MixRtAl   = xza_MixRtWork + xza_DelMixRt
      
      ! 雨から蒸気の値の保管
      !
!      call StorePotTempCond( &
!        & (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog & ! (in)
!        & ) 
!      call StoreMixRtCond( &
!        & xza_DelMixRt / DelTimeLFrog & ! (in)
!        & ) 
      Km1  = (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog 
      Tnd1 = xza_DelMixRt / DelTimeLFrog
      call HistoryAutoPut(TimeN, 'PotTempR2G',  &
        & Km1(FileXMin:FileXMax, FileZMin:FileZMax) )
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtR2G',  &
          & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s) )
      end do
      
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xz( xz_PotTempAl ) ! (inout)
      call BoundaryZSym_xz( xz_PotTempAl ) ! (inout)
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)
!      call SetMargin_aaz( xz_PotTempAl )  ! (inout)
!      call SetMargin_xyzf( xza_MixRtAl )  ! (inout)

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

      ! これまでの値を作業配列に保管
      ! Previous values are stored to work area.
      !
      xz_PotTempWork = xz_PotTempAl
      xza_MixRtWork  = xza_MixRtAl
      
      ! 湿潤調節法を適用
      ! Moist adjustment is applied.
      !
      call MoistAdjustSvapPress( &
        & xz_ExnerNl,            & ! (in)
        & xz_PotTempAl,          & ! (inout)
        & xza_MixRtAl            & ! (inout)
        & )
      call MoistAdjustNH4SH(     &
        & xz_ExnerNl,            & ! (in)
        & xz_PotTempAl,          & ! (inout)
        & xza_MixRtAl            & ! (inout)
        & )

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

      Km1 = (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog
      Tnd1= (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog
      call HistoryAutoPut(TimeN, 'PotTempG2C', Km1(FileXMin:FileXMax,FileZMin:FileZMax))
      do s = 1, SpcNum
        call HistoryAutoPut(TimeN, 'D'//trim(SpcWetSymbol(s))//'DtG2C', &
          & Tnd1(FileXMin:FileXMax, FileZMin:FileZMax, s))
      end do
      
      ! 境界条件 Boundary condition
      !
      call BoundaryXCyc_xz( xz_PotTempAl ) ! (inout)
      call BoundaryZSym_xz( xz_PotTempAl ) ! (inout)
      call BoundaryXCyc_xza( xza_MixRtAl ) ! (inout)
      call BoundaryZSym_xza( xza_MixRtAl ) ! (inout)

  end subroutine ArareCloudPhys

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

!    pz_VelXNl  =  pz_VelXNl   + tfil * (pz_VelXBl  - 2.0d0 * pz_VelXNl  + pz_VelXAl)
!    xr_VelZNl  =  xr_VelZNl   + tfil * (xr_VelZBl  - 2.0d0 * xr_VelZNl  + xr_VelZAl)
!    xz_PotTempNl =  xz_PotTempNl  + tfil * (xz_PotTempBl - 2.0d0 * xz_PotTempNl + xz_PotTempAl)
!    xz_ExnerNl =  xz_ExnerNl  + tfil * (xz_ExnerBl - 2.0d0 * xz_ExnerNl + xz_ExnerAl)
!    xz_KmNl    =  xz_KmNl     + tfil * (xz_KmBl    - 2.0d0 * xz_KmNl    + xz_KmAl)
!    xz_KhNl    =  xz_KhNl     + tfil * (xz_KhBl    - 2.0d0 * xz_KhNl    + xz_KhAl)
!    xza_MixRtNl = xza_MixRtNl + tfil * (xza_MixRtBl - 2.0d0 * xza_MixRtNl + xza_MixRtAl)

    pz_VelXNl  =  tfil * ( pz_VelXBl + pz_VelXAl )       + tfil2 * pz_VelXNl
    xr_VelZNl  =  tfil * ( xr_VelZBl + xr_VelZAl )       + tfil2 * xr_VelZNl
    xz_PotTempNl= tfil * ( xz_PotTempBl + xz_PotTempAl ) + tfil2 * xz_PotTempNl
    xz_ExnerNl =  tfil * ( xz_ExnerBl + xz_ExnerAl )     + tfil2 * xz_ExnerNl 
    xz_KmNl    =  tfil * ( xz_KmBl + xz_KmAl )           + tfil2 * xz_KmNl
    xz_KhNl    =  tfil * ( xz_KhBl + xz_KhAl )           + tfil2 * xz_KhNl
    xza_MixRtNl = tfil * ( xza_MixRtBl + xza_MixRtAl )   + tfil2 * xza_MixRtNl
    
  end subroutine AsselinTimeFilter

  !-----------------------------------------------------------------------
  subroutine MainInit
    ! 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 ファイルから取得する.
    ! 
    call ReStartFileio_BasicZ_Get
    call ReStartFileio_Var_Get(   &
      & pz_VelXBl,  pz_VelXNl,    & ! (out)
      & xr_VelZBl,  xr_VelZNl,    & ! (out)
      & xz_PotTempBl, xz_PotTempNl,   & ! (out)
      & xz_ExnerBl, xz_ExnerNl,   & ! (out)
      & xza_MixRtBl,xza_MixRtNl,  & ! (out)
      & xz_KmBl,    xz_KmNl,      & ! (out)
      & xz_KhBl,    xz_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( &
      & xza_MixRtBasicZ,    & ! (in)
      & xz_DensBasicZ       & ! (in)
      & )

    ! 放射強制の初期化
    !  Initialization of radiative forcing.
    !
    call Radiation_Init( &
      & cfgfile          & ! (in)    
      & )
    
    ! 湿潤気塊の浮力計算の初期化
    ! Initialization of moist buoyancy calculation.
    !
    call MoistBuoy_Init( ) 

    ! 圧力計算用係数行列の初期化 
    ! Initialization of coefficient matrix for exner function calculation.
    !
    call xz_Exner_Init( )              
    
    !-------------------------------------------------------------
    ! 基本場のファイル出力
    !
    call HistoryPut( 'DensBasicZ',   xz_DensBasicZ     , rstat)
    call HistoryPut( 'ExnerBasicZ',  xz_ExnerBasicZ    , rstat)
    call HistoryPut( 'PotTempBasicZ',  xz_PotTempBasicZ  , rstat)
    call HistoryPut( 'VelSoundBasicZ', xz_VelSoundBasicZ , rstat)
    call HistoryPut( 'TempBasicZ',    xz_TempBasicZ     , rstat)
    call HistoryPut( 'PressBasicZ',   xz_PressBasicZ    , rstat)
    call HistoryPut( 'MixRtBasicZ',    xza_MixRtBasicZ  , rstat)
    call HistoryPut( 'EffMolWtBasicZ', xz_EffMolWtBasicZ, rstat)
    
    ! 音波に対する CFL 条件のチェック
    ! CFL condtion check for sound wave.
    !
    call CFLCheckTimeShort( &
      & xz_VelSoundBasicZ   & ! (in)
      & )

    ! CPU time
    !
    call ClocksetPreStop

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

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

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

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate(                                                  &
      & xz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xz_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xz_ExnerBl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax),          &
!
      & xz_PotTempOrg(DimXMin:DimXMax, DimZMin:DimZMax),       &
      & xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax),      &
      & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAl(DimXMin:DimXMax, DimZMin:DimZMax),        &
!
      & xz_KmBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xz_KhBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xza_MixRtOrg(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),  &
      & xza_MixRtWork(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
      & xza_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
      & xza_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),   &
!
      & pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax),        & 
!
      & xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)   &
      & )
    
    pz_VelXNs = 0.0d0;     pz_VelXAs = 0.0d0    
    xr_VelZNs = 0.0d0;     xr_VelZAs = 0.0d0
    xz_ExnerNs = 0.0d0;    xz_ExnerAs = 0.0d0

    pz_VelXBl = 0.0d0;     pz_VelXNl = 0.0d0;     pz_VelXAl = 0.0d0
    xr_VelZBl = 0.0d0;     xr_VelZNl = 0.0d0;     xr_VelZAl = 0.0d0
    xz_ExnerBl = 0.0d0;    xz_ExnerNl = 0.0d0;    xz_ExnerAl = 0.0d0
    xz_KmBl = 0.0d0;       xz_KmNl = 0.0d0;       xz_KmAl = 0.0d0
    xz_KhBl = 0.0d0;       xz_KhNl = 0.0d0;       xz_KhAl = 0.0d0
    xz_PotTempBl = 0.0d0;  xz_PotTempNl = 0.0d0;  xz_PotTempAl = 0.0d0
    xza_MixRtBl = 0.0d0;    xza_MixRtNl = 0.0d0;    xza_MixRtAl = 0.0d0

    pz_AccelVelXNl = 0.0d0
    xr_AccelVelZNl = 0.0d0

    xza_DelMixRt = 0.0d0

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

    allocate(                                              &
         & Tnd1(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
         & Tnd2(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
         & Tnd3(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
         & Tnd4(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum),  &
         & Tnd5(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)  &
         &) 
    Tnd1 = 0.0d0
    Tnd2 = 0.0d0
    Tnd3 = 0.0d0
    Tnd4 = 0.0d0
    Tnd5 = 0.0d0

    ! 時間ループの初期化
    !
!    t = 1
    
  end subroutine ArareAlloc
!-----------------------------------------------------------------------
end program arare
