program dcpam_main
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! モデルの使い方については {チュートリアル}[link:../../../doc/tutorial/rakuraku/] を
  ! 参照してください.
  !
  ! See {Tutorial}[link:../../../doc/tutorial/rakuraku/index.htm.en] for usage of the 
  ! model. 
  !

  ! モジュール引用 ; USE statements
  !

  ! 力学過程 (スペクトル法, Arakawa and Suarez (1983))
  ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
  !
  use dynamics_hspl_vas83, only: Dynamics

  ! 物理過程のみの計算のための力学過程
  ! A dynamics for calculation with physical processes only
  !
  use dynamics_physicsonly, only: DynamicsPhysicsOnly

  ! Held and Suarez (1994) による強制と散逸
  ! Forcing and dissipation suggested by Held and Suarez (1994)
  !
  use held_suarez_1994, only: HS94Forcing

  ! 簡単金星計算のための強制
  ! forcing for simple Venus calculation
  !
  use yt2003_forcing, only: YT2003Forcing

  ! Schneider and Liu (2009) による鉛直混合課程
  ! Vertical diffusion by Schneider and Liu (2009)
  !
  use sl09_diffusion, only : SL09Diffusion

  ! 放射フラックス (GFD 電脳倶楽部開発の放射モデル)
  ! Radiation flux (radiation model developed by GFD Dennou Club)
  !
  use radiation_DennouAGCM, only: RadiationDennouAGCMFlux

  ! 放射関連ルーチン
  ! Routines for radiation calculation
  !
  use radiation_utils, only : RadiationDTempDt, RadiationFluxOutput

  use radiation_dcpam_E_V2, only: RadiationDcpamEV2Flux

  ! dcpam 火星大気向け放射モデル Ver. 1
  ! dcpam radiation model for the Mars' atmosphere Ver. 1
  !
  use radiation_dcpam_M_V1, only: RadiationDcpamMV1Flux

  ! 火星計算用近赤外加熱率計算
  ! Calculation of near infrared heating rate in the case of Mars
  !
  use radiation_dcpam_M_NIR, only : RadiationDcpamMNIRINOUT

  ! Schneider and Liu (2009) の放射モデル
  ! Radiation model by Schneider and Liu (2009)
  !
  use radiation_SL09, only : RadiationSL09Flux

  ! 鉛直拡散フラックス (Mellor and Yamada, 1974, レベル 2)
  ! Vertical diffusion flux (Mellor and Yamada, 1974, Level 2)
  !
  use vdiffusion_my1974, only: VDiffusion, VDiffusionOutput

  ! 地下における熱の鉛直拡散
  ! Vertical diffusion of heat under the ground
  !
  use subsurface_diffusion_heat, only: SubsurfaceDiffusion

  ! 積雲パラメタリゼーション (対流調節)
  ! Cumulus parameterization (convection adjust)
  !
  use moist_conv_adjust, only: MoistConvAdjust

  ! 大規模凝結
  ! Large scale condensation
  !
  use lscond, only: LScaleCond

  ! 地表面フラックス
  ! Surface flux
  use surface_flux_bulk, only: SurfaceFlux, SurfaceFluxOutput

  ! 乾燥対流調節
  ! Dry convective adjustment
  !
  use dry_conv_adjust, only: DryConvAdjust

  ! 質量の補正
  ! Mass fixer
  !
  use mass_fixer, only: MassFixer

  ! 温度の半整数σレベルの補間, 気圧と高度の算出
  ! Interpolate temperature on half sigma level, 
  ! and calculate pressure and height
  !
  use auxiliary, only: AuxVars

  ! 陰解法による時間積分のためのルーチン
  ! Routines for time integration with implicit scheme
  !
  use phy_implicit_utils, only : PhyImplEvalRadLFluxA

  ! 陰解法のための行列処理 (一部の物理過程用)
  ! Matrices handling for implicit scheme (for a part of physical processes)
  !
  use phy_implicit, only: PhyImplTendency

  ! 陰解法のための行列処理 (一部の物理過程用)
  ! Matrices handling for implicit scheme (for a part of physical processes)
  !
  use phy_implicit_sdh, only: PhyImplSDHTendency

  ! 陰解法による時間積分 (大気のみ / 惑星表面温度・土壌温度計算なし)
  ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
  !
  use phy_implicit_atmonly, only : PhyImplAtmOnlyTendency

  ! 地面温度の時間積分・地表面放射補正
  ! Time integration of surface temperature, correction of flux on surface
  !
  use intg_surftemp, only: IntegralSurfTemp

  ! バケツモデル
  ! Bucket model
  !
  use Bucket_Model, only : BucketModQvapFlux, BucketIntegration, BucketPRCPAdjust

  ! 鉛直フィルター (Ishiwatari et al., 2002)
  ! Vertical filter (Ishiwatari et al., 2002)
  !
  use vertical_filter, only: VerticalFilter

  ! タイムフィルター (Asselin, 1972)
  ! Time filter (Asselin, 1972)
  !
  use timefilter_asselin1972, only: TimeFilter, TimeFilterSurfVars

  ! 時刻管理
  ! Time control
  !
  use timeset, only: TimesetProgress, TimeB, TimeN, TimeA, EndTime, DelTime                 ! $ \Delta t $ [s]

  ! リスタートデータ入出力
  ! Restart data input/output
  !
  use restart_file_io, only: RestartFileOutPut

  ! 地表面温度リスタートデータ入出力
  ! Restart data of surface temperature input/output
  !
  use restart_surftemp_io, only: RestartSurfTempOutPut

  ! 惑星表面特性の設定
  ! Setting of surface properties
  !
  use surface_properties, only: SetSurfaceProperties

  ! ヒストリデータ出力
  ! History data output
  !
  use gtool_historyauto, only: HistoryAutoPut, HistoryAutoAllVarFix

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only: ncmax, a_QMixName, a_QMixLongName, IndexH2OVap, CompositionInqIndex


  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, jmax, kmax     ! 鉛直層数. 
                              ! Number of vertical level

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only: LatentHeat

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, STRING, TOKEN      ! キーワード.   Keywords. 

  use set_cloud, only: SetCloudCloudWaterLossRateInOut, SetCloudRegDQCloudWaterDt


  ! CO2 相変化 (火星計算用)
  ! Phase change of CO2 (for the use in Mars calculation)
  !
  use co2_phase_change, only : CO2PhaseChangeLimitTemp

  ! 予報変数の値の確認
  ! Check values of prognostic variables
  !
  use check_prog_vars, only: CheckProgVars

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 予報変数 (ステップ $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  ! Prediction variables  (Step $ t-\Delta t $ , $ t $ , $ t+\Delta t $ )
  !
  real(DP), allocatable:: xyz_UB (:,:,:)
                              ! $ u (t-\Delta t) $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyz_VB (:,:,:)
                              ! $ v (t-\Delta t) $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable:: xyz_TempB (:,:,:)
                              ! $ T (t-\Delta t) $ .   温度. Temperature (K)
  real(DP), allocatable:: xyzf_QMixB(:,:,:,:)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xy_PsB (:,:)
                              ! $ p_s (t-\Delta t) $ . 地表面気圧. Surface pressure (Pa)
  real(DP), allocatable:: xyz_UN (:,:,:)
                              ! $ u (t) $ .     東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyz_VN (:,:,:)
                              ! $ v (t) $ .     南北風速. Northward wind (m s-1)
  real(DP), allocatable:: xyz_TempN (:,:,:)
                              ! $ T (t) $ .     温度. Temperature (K)
  real(DP), allocatable:: xyzf_QMixN(:,:,:,:)
                              ! $ q (t) $ .     混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xy_PsN (:,:)
                              ! $ p_s (t) $ .   地表面気圧. Surface pressure (Pa)
  real(DP), allocatable:: xyz_UA (:,:,:)
                              ! $ u (t+\Delta t) $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable:: xyz_VA (:,:,:)
                              ! $ v (t+\Delta t) $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable:: xyz_TempA (:,:,:)
                              ! $ T (t+\Delta t) $ .   温度. Temperature (K)
  real(DP), allocatable:: xyzf_QMixA(:,:,:,:)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable:: xy_PsA (:,:)
                              ! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure (Pa)

  real(DP), allocatable:: xy_SoilMoistB (:,:)
                              ! $ M_ws (t-\Delta t) $ . 土壌水分 (kg m-2)
                              ! Soil moisture (kg m-2)
  real(DP), allocatable:: xy_SurfSnowB (:,:)
                              ! $ M_ss (t-\Delta t) $ . 積雪量 (kg m-2)
                              ! Surface snow amount (kg m-2)
  real(DP), allocatable:: xy_SoilMoistN (:,:)
                              ! $ M_ws (t) $          . 土壌水分 (kg m-2)
                              ! Soil moisture (kg m-2)
  real(DP), allocatable:: xy_SurfSnowN (:,:)
                              ! $ M_ss (t) $ . 積雪量 (kg m-2)
                              ! Surface snow amount (kg m-2)
  real(DP), allocatable:: xy_SoilMoistA (:,:)
                              ! $ M_ws (t+\Delta t) $ . 土壌水分 (kg m-2)
                              ! Soil moisture (kg m-2)
  real(DP), allocatable:: xy_SurfSnowA (:,:)
                              ! $ M_ss (t+\Delta t) $ . 積雪量 (kg m-2)
                              ! Surface snow amount (kg m-2)


  ! 診断変数, 他
  ! Diagnostic variables, etc.
  !
  real(DP), allocatable:: xyz_DUDt (:,:,:)
                              ! $ \DP{u}{t} $ . 東西風速変化 (m s-2)
                              ! Eastward wind tendency (m s-2)
  real(DP), allocatable:: xyz_DVDt (:,:,:)
                              ! $ \DP{v}{t} $ . 南北風速変化 (m s-2)
                              ! Northward wind tendency (m s-2)
  real(DP), allocatable:: xyz_DTempDt (:,:,:)
                              ! $ \DP{T}{t} $ . 温度変化 (K s-1)
                              ! Temperature tendency (K s-1)
  real(DP), allocatable:: xyzf_DQMixDt(:,:,:,:)
                              ! $ \DP{q}{t} $ . 混合比変化 (s-1)
                              ! Mass mixing ratio tendency (s-1)

  real(DP), allocatable:: xy_SurfHeight (:,:)
                              ! $ z_s $ . 地表面高度 (m)
                              ! Surface height (m)

  real(DP), allocatable:: xy_SurfTemp (:,:)
                              ! 地表面温度 (K)
                              ! Surface temperature (K)
  real(DP), allocatable:: xyz_SoilTemp(:,:,:)
                              ! 土壌温度 (K)
                              ! Soil temperature (K)


  real(DP), allocatable:: xy_SurfAlbedo (:,:)
                              ! 地表アルベド (1)
                              ! Surface albedo (1)
  real(DP), allocatable:: xy_SurfHumidCoef (:,:)
                              ! 地表湿潤度 (1)
                              ! Surface humidity coefficient (1)
  real(DP), allocatable:: xy_SurfRoughLength (:,:)
                              ! 地表粗度長 (m)
                              ! Surface rough length (m)
  real(DP), allocatable:: xy_SurfHeatCapacity (:,:)
                              ! 地表熱容量 (J m-2 K-1)
                              ! Surface heat capacity (J m-2 K-1)
  real(DP), allocatable:: xy_SeaIceConc(:,:)
                              ! 海氷密度 (0 <= xy_SeaIceConc <= 1) (1)
                              ! Sea ice concentration (0 <= xy_SeaIceConc <= 1) (1)
  integer , allocatable:: xy_SurfCond (:,:)
                              ! 地表状態 (0: 固定, 1: 可変) (1)
                              ! Surface condition (0: fixed, 1: variable) (1)
  real(DP), allocatable:: xy_DeepSubSurfHeatFlux (:,:)
                              ! 地中熱フラックス (W m-2)
                              ! "Deep subsurface heat flux" (W m-2)
                              ! Heat flux at the bottom of surface/soil layer.
  real(DP), allocatable:: xy_SoilHeatCap (:,:)
                              ! 土壌熱容量 (J K-1 kg-1)
                              ! Specific heat of soil (J K-1 kg-1)
  real(DP), allocatable:: xy_SoilHeatDiffCoef (:,:)
                              ! 土壌熱伝導係数 (J m-3 K-1)
                              ! Heat conduction coefficient of soil (J m-3 K-1)

  real(DP), allocatable:: xyr_Temp (:,:,:)
                              ! $ \hat{T} $ . 温度 (半整数レベル) (K)
                              ! Temperature (half level) (K)
  real(DP), allocatable:: xyz_Press (:,:,:)
                              ! $ p $ . 気圧 (整数レベル) (Pa)
                              ! Air pressure (full level) (Pa)
  real(DP), allocatable:: xyr_Press (:,:,:)
                              ! $ \hat{p} $ . 気圧 (半整数レベル) (Pa)
                              ! Air pressure (half level) (Pa)
  real(DP), allocatable:: xyz_Height (:,:,:)
                              ! 高度 (整数レベル) (m)
                              ! Height (full level) (m)
  real(DP), allocatable:: xyr_Height (:,:,:)
                              ! 高度 (半整数レベル) (m)
                              ! Height (half level) (m)
  real(DP), allocatable:: xyz_Exner (:,:,:)
                              ! Exner 関数 (整数レベル) (1)
                              ! Exner function (full level) (1)
  real(DP), allocatable:: xyr_Exner (:,:,:)
                              ! Exner 関数 (半整数レベル) (1)
                              ! Exner function (half level) (1)

  real(DP), allocatable:: xyr_RadLFlux (:,:,:)
                              ! 長波フラックス (W m-2)
                              ! Longwave flux (W m-2)
  real(DP), allocatable:: xyr_RadLFluxA (:,:,:)
                              ! 陰解法で解いた地表面熱収支と整合的な $ t+\Delta t $ に
                              ! おける長波フラックスの計算
                              !   * ここで計算された値が直接次のステップ $ t $ における
                              !     長波フラックスとして用いられるわけではない
                              !   * 現在の時間ステップにおける長波放射加熱率の計算に使
                              !     われる
                              !
                              ! Evaluate longwave flux at $ t+\Delta t $ consistent 
                              ! with surface energy balance solved with implicit method
                              !   * The evaluated value is not used directly as Longwave
                              !     flux at next step $ t $.
                              !   * The evaluated value is used to calculate long wave 
                              !     radiative heating rate in the current time step.
  real(DP), allocatable:: xyr_RadSFlux (:,:,:)
                              ! 短波 (日射) フラックス (W m-2)
                              ! Shortwave (insolation) flux (W m-2)
  real(DP), allocatable:: xyra_DelRadLFlux (:,:,:,:)
                              ! 長波地表温度変化 (W m-2)
                              ! Surface temperature tendency with longwave (W m-2)

  real(DP), allocatable:: xyr_MomFluxX (:,:,:)
                              ! 東西方向運動量フラックス
                              ! Eastward momentum flux
  real(DP), allocatable:: xyr_MomFluxY (:,:,:)
                              ! 南北方向運動量フラックス. 
                              ! Northward momentum flux
  real(DP), allocatable:: xyr_HeatFlux (:,:,:)
                              ! 熱フラックス. 
                              ! Heat flux
  real(DP), allocatable:: xyrf_QMixFlux(:,:,:,:)
                              ! 成分質量フラックス. 
                              ! Mass flux of compositions

  real(DP), allocatable:: xyr_SoilHeatFlux (:,:,:)
                              ! 土壌中の熱フラックス (W m-2)
                              ! Heat flux in sub-surface soil (W m-2)

  real(DP), allocatable:: xyr_VelTransCoef (:,:,:)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
  real(DP), allocatable:: xyr_TempTransCoef (:,:,:)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
  real(DP), allocatable:: xyr_QMixTransCoef(:,:,:)
                              ! 輸送係数：混合比.
                              ! Transfer coefficient: mixing ratio

  real(DP), allocatable:: xy_SurfVelTransCoef (:,:)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
  real(DP), allocatable:: xy_SurfTempTransCoef (:,:)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
  real(DP), allocatable:: xy_SurfQVapTransCoef (:,:)
                              ! 輸送係数：水蒸気
                              ! Transfer coefficient: water vapor
  real(DP), allocatable:: xyr_SoilTempTransCoef (:,:,:)
                              ! 輸送係数：土壌温度.
                              ! Transfer coefficient: soil temperature

  real(DP), allocatable:: xy_DSurfTempDt (:,:)
                              ! 地表面温度変化率. 
                              ! Surface temperature tendency
  real(DP), allocatable:: xyz_DSoilTempDt (:,:,:)
                              ! $ \DP{Tg}{t} $ . 土壌温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

  real(DP), allocatable:: xy_DSoilMoistDt (:,:)
                              ! 土壌温度時間変化率 (kg m-2 s-1)
                              ! Soil temperature tendency (kg m-2 s-1)
  real(DP), allocatable:: xy_DSurfSnowDt (:,:)
                              ! 積雪率時間変化率 (kg m-2 s-1)
                              ! Surface snow amount tendency (kg m-2 s-1)

  real(DP), allocatable:: xyz_DTempDtVDiff(:,:,:)
                              ! 鉛直拡散による加熱率 (K s-1)
                              ! Temperature tendency due to vertical diffusion (K s-1)

  real(DP), allocatable:: xyz_DTempDtRadL (:,:,:)
                              ! 長波加熱率. 
                              ! Temperature tendency with longwave
  real(DP), allocatable:: xyz_DTempDtRadS (:,:,:)
                              ! 短波加熱率. 
                              ! Temperature tendency with shortwave

  real(DP), allocatable:: xy_Rain (:,:)
                              ! 降水量. 
                              ! Precipitation
  real(DP), allocatable:: xyz_DTempDtCond (:,:,:)
                              ! 凝結加熱率. 
                              ! Condensation heating
  real(DP), allocatable:: xyz_DQVapDtCond (:,:,:)
                              ! 凝結比湿変化. 
                              ! Condensation specific humidity tendency

  real(DP), allocatable:: xyz_DDelLWDtCum(:,:,:)
                              ! Production rate of liquid water in the layer 
                              ! due to condensation in cumulus convection 
                              ! parameterization (kg m-2 s-1)
  real(DP), allocatable:: xyz_DDelLWDtLSC(:,:,:)
                              ! Production rate of liquid water in the layer 
                              ! due to condensation in large scale condensation 
                              ! (kg m-2 s-1)

  real(DP), allocatable:: xyz_DQH2OLiqDtCum(:,:,:)
                              ! Production rate of cloud water in the layer 
                              ! due to condensation in cumulus convection 
                              ! parameterization (kg kg-1)
  real(DP), allocatable:: xyz_DQH2OLiqDtLSC(:,:,:)
                              ! Production rate of cloud water in the layer 
                              ! due to condensation in large scale condensation 
                              ! (kg kg-1)


  ! 作業変数
  ! Work variables
  !
  integer           :: IDDynMode                 ! 使用する力学過程
                                                 ! Dynamics used for an experiment

  integer, parameter:: IDDynModeHSPLVAS83       = 0
  integer, parameter:: IDDynModeNoHorAdv        = 1

  integer           :: IDPhyMode                 ! 使用する物理過程
                                                 ! Physics used for an experiment

  integer, parameter:: IDPhyModeNoPhysics     = 0
  integer, parameter:: IDPhyModeFullPhysics   = 1
  integer, parameter:: IDPhyModeHS94          = 2
  integer, parameter:: IDPhyModeVenusSimple   = 3
  integer, parameter:: IDPhyModeJupiterSimple = 4

  integer           :: IDPhyTendMethod           ! 物理過程による変化率の計算方法
                                                 ! Method calculating physics tendency

  integer, parameter:: IDPhyTendMethodImp1LayModel = 10
  integer, parameter:: IDPhyTendMethodImpSoilModel = 11
  integer, parameter:: IDPhyTendMethodImpAtmOnly   = 12

  integer           :: IDRadMethod               ! 放射過程の計算方法
                                                 ! Method for radiation

  integer, parameter:: IDRadMethodDennouAGCM = 20
  integer, parameter:: IDRadMethodDcpamEV1   = 21
  integer, parameter:: IDRadMethodDcpamMV1   = 22
  integer, parameter:: IDRadMethodSL09       = 23

  logical:: FlagPhyImpSoilModelSO   ! flag for use of slab ocean
  logical:: FlagVerticalFilter      ! flag for use of vertical filter


  logical:: firstloop = .true.
                              ! 初回のループであることを示すフラグ. 
                              ! Flag implying first loop

  logical:: flag_initial
                              ! 内部サブルーチン MainInit で設定されます. 
                              ! リスタートデータを読み込む場合には, 
                              ! .false. が, 初期値データを読み込む場合には
                              ! .true. が設定されます. 
                              ! 
                              ! This variable is set in an internal 
                              ! subroutine "MainInit". 
                              ! If restart data is loaded, .false. is set. 
                              ! On the other hand, if initial data is loaded, 
                              ! .true. is set.

  integer:: n                 ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

  ! 実行文 ; Executable statement
  !

  ! 主プログラムの初期化 (内部サブルーチン)
  ! Initialization for the main program (Internal subroutine)
  !
  call MainInit


  ! 時間積分
  ! Time integration
  !
  loop_time : do while ( TimeB < EndTime )

    ! 地表面高度の設定
    ! Set surface height
    !
    call SetSurfaceProperties( xy_SurfHeight = xy_SurfHeight )


    select case ( IDPhyMode )
    case ( IDPhyModeNoPhysics )

      xyz_DUDt     = 0.0d0
      xyz_DVDt     = 0.0d0
      xyz_DTempDt  = 0.0d0
      xyzf_DQMixDt = 0.0d0

    case ( IDPhyModeHS94 )

      ! Held and Suarez (1994) による強制と散逸
      ! Forcing and dissipation suggested by Held and Suarez (1994)
      !
      call HS94Forcing( xyz_UB, xyz_VB, xyz_TempB, xy_PsB, xyz_DUDt, xyz_DVDt, xyz_DTempDt )    ! (out)
      xyzf_DQMixDt = 0.

    case ( IDPhyModeVenusSimple )

      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB,     xyz_TempB, xyr_Temp, xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xy_SurfHeight = xy_SurfHeight, xyz_Height    = xyz_Height, xyr_Height    = xyr_Height, xyz_Exner     = xyz_Exner,  xyr_Exner  = xyr_Exner )

      ! 簡単金星計算のための強制
      ! forcing for simple Venus calculation
      !
      call YT2003Forcing( xy_SurfHeight, xyz_UB, xyz_VB, xyz_TempB, xy_PsB, xyz_Press, xyr_Press, xyr_Temp, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_DUDt, xyz_DVDt, xyz_DTempDt )

      xyzf_DQMixDt = 0.0d0


    case ( IDPhyModeJupiterSimple )

      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB,     xyz_TempB, xyr_Press     = xyr_Press )


      ! Schneider and Liu (2009) による鉛直混合課程
      ! Vertical diffusion by Schneider and Liu (2009)
      !
      call SL09Diffusion( xyz_UB, xyz_VB, xyr_Press, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff )

      ! Schneider and Liu (2009) の放射モデル
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadiationSL09Flux( xyr_Press, xyz_TempB, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux )

      ! 放射による温度変化率
      ! Temperature tendency with radiation
      !
      call RadiationDTempDt( xyr_RadLFlux, xyr_RadSFlux, xyr_Press, xyz_DTempDtRadL, xyz_DTempDtRadS )          ! (out)

      ! 非断熱加熱率の総和の計算
      ! Sum all diabatic heating rates
      !
      xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS

      xyzf_DQMixDt = 0.0_DP


    case ( IDPhyModeFullPhysics )

      ! 地表面条件の設定
      ! Configure surface conditions
      !
      call SetSurfaceProperties( FlagPhyImpSoilModelSO, xy_SoilMoistB, xy_SurfSnowB, xy_SurfTemp            = xy_SurfTemp, xy_SurfAlbedo          = xy_SurfAlbedo, xy_SurfHumidCoef       = xy_SurfHumidCoef, xy_SurfRoughLength     = xy_SurfRoughLength, xy_SurfHeatCapacity    = xy_SurfHeatCapacity, xy_DeepSubSurfHeatFlux = xy_DeepSubSurfHeatFlux, xy_SurfCond            = xy_SurfCond, xy_SeaIceConc          = xy_SeaIceConc, xy_SoilHeatCap         = xy_SoilHeatCap, xy_SoilHeatDiffCoef    = xy_SoilHeatDiffCoef )

      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB,     xyz_TempB, xyr_Temp, xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xy_SurfHeight = xy_SurfHeight, xyz_Height    = xyz_Height, xyr_Height = xyr_Height, xyz_Exner     = xyz_Exner,  xyr_Exner  = xyr_Exner )


      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )

        ! 放射フラックス (GFD 電脳倶楽部開発の放射モデル)
        ! Radiation flux (radiation model developed by GFD Dennou Club)
        !
        call RadiationDennouAGCMFlux( xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_Press, xy_SurfTemp, xy_SurfAlbedo, xyr_RadLFlux, xyr_RadSFlux, xyra_DelRadLFlux, flag_rst = .not. flag_initial )

      case ( IDRadMethodDcpamEV1 )

        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          stop 'H2OLiq is not present.'
        end if

        call RadiationDcpamEV2Flux( xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' )), xyz_Height, xy_SurfTemp, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux )

      case ( IDRadMethodDcpamMV1 )

        call RadiationDcpamMV1Flux( xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_TempB, xy_SurfTemp, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux )

      case ( IDRadMethodSL09 )

        ! Schneider and Liu (2009) の放射モデル
        ! Radiation model by Schneider and Liu (2009)
        !
        call RadiationSL09Flux( xyr_Press, xyz_TempB, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux )

        ! This is done no to correct long wave radiation flux by PhyImplEvalRadLFluxA.
        ! (yot, 2010/09/28)
        !
        xyra_DelRadLFlux = 0.0_DP

      end select


      ! 鉛直拡散フラックス
      ! Vertical diffusion flux
      !
      call VDiffusion( xyz_UB,     xyz_VB,     xyzf_QMixB, xyz_TempB,  xyr_Temp,   xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner,    xyr_Exner, xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )

      ! 地表面フラックス
      ! Surface flux
      !
      call SurfaceFlux( xyz_UB, xyz_VB, xyz_TempB, xyr_Temp, xyzf_QMixB, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfRoughLength, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )

      ! 一部の物理過程の時間変化率の計算 (陰解法)
      ! Calculate tendency by a part of physical processes (implicit)
      !
      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )

        call PhyImplTendency( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfCond, xy_SurfHeatCapacity, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt, xy_DSurfTempDt )

        xyz_DSoilTempDt = 0.0_DP
        xy_DSoilMoistDt = 0.0_DP
        xy_DSurfSnowDt  = 0.0_DP

      case ( IDPhyTendMethodImpSoilModel )

        ! バケツモデルのための地表面フラックス修正
        ! Modification of surface flux for bucket model
        !
        call BucketModQvapFlux( xy_SurfCond, xy_SoilMoistB, xy_SurfSnowB, xyrf_QMixFlux(:,:,0,IndexH2OVap) )

        ! 地下における熱の鉛直拡散
        ! Vertical diffusion of heat under the ground
        !
        call SubsurfaceDiffusion( xy_DeepSubSurfHeatFlux, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_SurfTemp, xyz_SoilTemp, xyr_SoilTempTransCoef, xyr_SoilHeatFlux )

        call PhyImplSDHTendency( FlagPhyImpSoilModelSO, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xy_SurfTemp, xyz_SoilTemp, xy_SurfHumidCoef, xy_SurfCond, xy_SurfHeatCapacity, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_SeaIceConc, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfSnowB, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DSoilMoistDt, xy_DSurfSnowDt )

      case ( IDPhyTendMethodImpAtmOnly )

        call PhyImplAtmOnlyTendency( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt )

        xy_DSurfTempDt  = 0.0_DP
        xyz_DSoilTempDt = 0.0_DP
        xy_DSoilMoistDt = 0.0_DP
        xy_DSurfSnowDt  = 0.0_DP

      end select


      ! 陰解法で解いた地表面熱収支と整合的な $ t+\Delta t $ における長波フラックスの
      ! 計算
      !   * ここで計算された値が直接次のステップ $ t $ における長波フラックスとして
      !   * 用いられるわけではない
      !   * 現在の時間ステップにおける長波放射加熱率の計算に使われる
      !
      ! Evaluate longwave flux at $ t+\Delta t $ consistent with surface energy balance
      ! solved with implicit method
      !   * The evaluated value is not used directly as Longwave flux at next step $t$.
      !   * The evaluated value is used to calculate long wave radiative heating rate 
      !     in the current time step.
      !
      call PhyImplEvalRadLFluxA( xyr_RadLFlux, xyz_DTempDtVDiff, xy_DSurfTempDt, xyra_DelRadLFlux, xyr_RadLFluxA )                                  ! (out)

      ! 放射による温度変化率
      ! Temperature tendency with radiation
      !
      call RadiationDTempDt( xyr_RadLFluxA, xyr_RadSFlux, xyr_Press, xyz_DTempDtRadL, xyz_DTempDtRadS )          ! (out)


      ! 非断熱加熱率の総和の計算
      ! Sum all diabatic heating rates
      !
      xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS


      select case ( IDRadMethod )
      case ( IDRadMethodDcpamMV1 )
        ! 火星計算用近赤外加熱率計算
        ! Calculation of near infrared heating rate in the case of Mars
        !
        call RadiationDcpamMNIRINOUT( xyz_Press, xyz_DTempDt )
      end select


      if ( CompositionInqIndex( 'H2OLiq' ) > 0 ) then

        ! Add loss rate of cloud water
        !
        call SetCloudCloudWaterLossRateInOut( xyzf_QMixB(:,:,:,CompositionInqIndex('H2OLiq')), xyzf_DQMixDt(:,:,:,CompositionInqIndex('H2OLiq')) )

      end if


      ! 鉛直拡散フラックスの出力 
      !   * 出力のみのサブルーチンであり, 計算には影響しない
      ! 
      ! Output Vertical diffusion fluxes
      !   * This subroutine works for output only, 
      !     so it does not influence a calculation.
      !
      call VDiffusionOutput( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt,  xyz_DVDt,  xyz_DTempDtVDiff,  xyzf_DQMixDt, xyz_Exner, xyr_Exner, xyr_VelTransCoef, xyr_TempTransCoef, xyr_QMixTransCoef )

      ! 地表面フラックスの出力 
      !   * 出力のみのサブルーチンであり, 計算には影響しない
      ! 
      ! Output surface fluxes
      !   * This subroutine works for output only, 
      !     so it does not influence a calculation.
      !
      call SurfaceFluxOutput( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt, xy_SurfTemp, xy_DSurfTempDt, xyr_Press, xyz_Exner, xyr_Exner, xy_SurfHumidCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )

      ! 放射フラックスの出力 
      !   * 出力のみのサブルーチンであり, 計算には影響しない
      ! 
      ! Output radiation fluxes
      !   * This subroutine works for output only, 
      !     so it does not influence a calculation.
      !
      call RadiationFluxOutput( xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux, xy_DSurfTempDt, xyz_DTempDtVDiff )

    end select


    ! 地面温度・土壌温度・土壌水分・積雪量の積分
    ! Time integration of surface temperature, soil temperature, soil moisture, 
    ! and surface snow amount
    !
    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )

      ! 地面温度・土壌温度の時間積分
      ! Time integration of surface temperature and soil temperature
      !
      call IntegralSurfTemp( xy_DSurfTempDt, xyz_DSoilTempDt, xy_SurfTemp   , xyz_SoilTemp )

      ! 土壌水分・地面積雪量の時間積分
      ! Time integration of soil moisture and snow amount
      !
      call BucketIntegration( xy_SurfCond, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_SoilMoistB, xy_SurfSnowB, xy_SoilMoistA, xy_SurfSnowA )

    end select


    ! 力学過程
    ! Dynamical core
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      call Dynamics( xyz_UB,   xyz_VB,   xyz_TempB,   xyzf_QMixB,   xy_PsB, xyz_UN,   xyz_VN,   xyz_TempN,   xyzf_QMixN,   xy_PsN, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_SurfHeight, xyz_UA,   xyz_VA,   xyz_TempA,   xyzf_QMixA,   xy_PsA )
    case ( IDDynModeNoHorAdv )
      call DynamicsPhysicsOnly( xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_PsB, xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsA, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA )
    end select


    select case ( IDPhyMode )
    case ( IDPhyModeJupiterSimple )
      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA,     xyz_TempA, xyz_Press = xyz_Press, xyr_Press = xyr_Press )  ! (out) optional

      ! 乾燥対流調節
      ! Dry convective adjustment
      !
      call DryConvAdjust( xyz_TempA, xyz_Press, xyr_Press )

    case ( IDPhyModeFullPhysics )
      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA,     xyz_TempA, xyz_Press = xyz_Press, xyr_Press = xyr_Press )  ! (out) optional

      ! 積雲パラメタリゼーション
      ! Cumulus parameterization
      !
      xy_Rain         = 0.
      xyz_DTempDtCond = 0.
      xyz_DQVapDtCond = 0.
      ! new routine
      call MoistConvAdjust( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_Rain, xyz_DTempDtCond, xyz_DQVapDtCond, xyz_Press, xyr_Press, xyz_DDelLWDtCum, xyz_DQH2OLiqDtCum )

      ! 大規模凝結
      ! Large scale condensation
      !
      call LScaleCond( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_Rain, xyz_DTempDtCond, xyz_DQVapDtCond, xyz_Press, xyr_Press, xyz_DDelLWDtLSC, xyz_DQH2OLiqDtLSC )


      call SetCloudRegDQCloudWaterDt( 'CUM', xyz_DQH2OLiqDtCum )
      call SetCloudRegDQCloudWaterDt( 'LSC', xyz_DQH2OLiqDtLSC )


      ! バケツモデル, 降水に伴う地表水量変化の計算
      ! bucket model, calculation of change of surface water due to precipitation
      !
      call BucketPRCPAdjust( xy_SurfCond, xy_PsA, xy_Rain, xyz_TempA, xy_SoilMoistA, xy_SurfSnowA )


      ! 乾燥対流調節
      ! Dry convective adjustment
      !
      call DryConvAdjust( xyz_TempA, xyz_Press, xyr_Press )


      ! CO2 相変化 (火星計算用)
      ! Phase change of CO2 (for the use in Mars calculation)
      !
      call CO2PhaseChangeLimitTemp( xyr_Press, xyz_Press, xy_SurfTemp, xyz_TempA )

      ! 成分の質量の補正
      ! Fix masses of constituents
      !
      call MassFixer( xyr_Press, xyzf_QMixA )

    end select

    ! Vertical filter
    ! 鉛直フィルタ
    ! 
    if ( FlagVerticalFilter ) then

      ! 温度の半整数σレベルの補間, 気圧と高度の算出
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA,     xyz_TempA, xyr_Temp  = xyr_Temp, xyr_Press = xyr_Press )  ! (out) optional

      call VerticalFilter( xyz_UA, xyz_VA, xyz_TempA, xyr_Press, xyr_Temp )

    end if

    ! 時間フィルター
    ! Time filter
    !
    if ( .not. flag_initial .or. .not. firstloop ) then
      call TimeFilter( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA  )   ! (in)

      select case ( IDPhyMode )
      case ( IDPhyModeFullPhysics )
        call TimeFilterSurfVars( xy_SoilMoistB, xy_SurfSnowB, xy_SoilMoistN, xy_SurfSnowN, xy_SoilMoistA, xy_SurfSnowA )
      end select
    end if

    write( 6, * ) xyz_TempA(0,jmax/2+1,1)

    ! 予報変数の値の確認
    ! Check values of prognostic variables
    !
    call CheckProgVars( xy_PsA, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA )

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeA, 'U',    xyz_UA )
    call HistoryAutoPut( TimeA, 'V',    xyz_VA )
    call HistoryAutoPut( TimeA, 'Temp', xyz_TempA )
    do n = 1, ncmax
      call HistoryAutoPut( TimeA, a_QMixName(n), xyzf_QMixA(:,:,:,n) )
    end do
    call HistoryAutoPut( TimeA, 'Ps',   xy_PsA )

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      call HistoryAutoPut( TimeN, 'SurfTemp', xy_SurfTemp )
      if ( size( xyz_SoilTemp ) /= 0 ) then
        call HistoryAutoPut( TimeN, 'SoilTemp', xyz_SoilTemp )
      end if

      call HistoryAutoPut( TimeA, 'SoilMoist', xy_SoilMoistA )
      call HistoryAutoPut( TimeA, 'SurfSnow' , xy_SurfSnowA  )

      call HistoryAutoPut( TimeN, 'Rain', xy_Rain * LatentHeat )
      call HistoryAutoPut( TimeN, 'DTempDtCond', xyz_DTempDtCond )
      call HistoryAutoPut( TimeN, 'DQVapDtCond', xyz_DQVapDtCond )

      call HistoryAutoPut( TimeN, 'SeaIceConc', xy_SeaIceConc )
      call HistoryAutoPut( TimeN, 'SurfAlbedo', xy_SurfAlbedo )
    end select

    ! 次の時間ステップに向けて予報変数の入れ替え
    ! Exchange prediction variables for the next time step
    !
    xyz_UB     = xyz_UN     
     xyz_UN     = xyz_UA     
     xyz_UA     = 0.
    xyz_VB     = xyz_VN     
     xyz_VN     = xyz_VA     
     xyz_VA     = 0.
    xyz_TempB  = xyz_TempN  
     xyz_TempN  = xyz_TempA  
     xyz_TempA  = 0.
    xyzf_QMixB = xyzf_QMixN 
     xyzf_QMixN = xyzf_QMixA 
     xyzf_QMixA = 0.
    xy_PsB     = xy_PsN     
     xy_PsN     = xy_PsA     
     xy_PsA     = 0.


    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      xy_SoilMoistB = xy_SoilMoistN
      xy_SoilMoistN = xy_SoilMoistA
      xy_SoilMoistA = 0.0d0

      xy_SurfSnowB  = xy_SurfSnowN
      xy_SurfSnowN  = xy_SurfSnowA
      xy_SurfSnowA  = 0.0d0
    end select

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

    ! NAMELIST から読み込んだ変数名に無効なものが存在したかどうかをチェック
    ! HistoryAutoAddVariable で登録した変数名を印字
    !
    ! Check that invalid variable names are loaded from NAMELIST or not
    ! Print registered variable names by "HistoryAutoAddVariable"
    !
    !!! if ( firstloop ) call HistoryAutoAllVarFix

    ! リスタートデータ出力
    ! Restart data output
    !
    call RestartFileOutput( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN(:,:,:,IndexH2OVap), xy_PsN  )  ! (in)

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      ! 地表面温度リスタートデータ出力
      ! Restart data of surface temperature output
      !
      call RestartSurfTempOutput( xy_SurfTemp  )             ! (in)
    end select

    firstloop = .false.

  ! 時間積分終了
  ! Time integration is finished
  !
  end do loop_time

  ! 主プログラムの終了処理 (内部サブルーチン)
  ! Termination for the main program (Internal subroutine)
  !
  call MainTerminate



contains

  !-------------------------------------------------------------------

  subroutine MainInit
    !
    ! 主プログラムの初期化手続き. 
    !
    ! Initialization procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperInit

    use dc_message, only: MessageNotify

    ! コマンドライン引数処理
    ! Command line option parser
    !
    use option_parser, only: OptParseInit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: NmlutilInit, NmlutilMsg, namelist_filename

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetInit, TimesetDelTimeHalf, TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    use fileset, only: FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: GridsetInit, imax, jmax, kmax, kslmax    ! 地下の鉛直層数. 
                                 ! Number of subsurface vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: CompositionInit

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: ConstantsInit

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: ConstantsSnowSeaIceInit

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetInit

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileOpen, RestartFileGet

    ! 地表面温度リスタートデータ入出力
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempOpen, RestartSurfTempGet

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileOpen
    use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT               ! 標準出力の装置番号. Unit number of standard output

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 宣言文 ; Declaration statements
    !
    implicit none
    character(*), parameter:: prog_name = 'dcpam_main'
                              ! 主プログラム名. 
                              ! Main program name
    character(*), parameter:: version = '$Name: dcpam5-20110221-2 $' // '$Id: dcpam_main.f90,v 1.18 2011-02-18 04:45:17 yot Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version

    character(STRING)      :: opt_nml_filename
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name

    logical:: FlagDynamics      ! 力学過程計算用フラグ
                                ! Flag for dynamics

    logical:: FlagFullPhysics   ! 全物理過程計算用フラグ
                                ! Flag for full physics
    logical:: FlagHS94          ! Held and Suarez (1994) 強制オン／オフ. 
                                ! Held and Suarez (1994) forcing on/off. 
    logical:: FlagVenusSimple   ! 金星簡単強制用フラグ
                                ! Flag for simple forcing for Venus
    logical:: FlagJupiterSimple !
                                ! Flag for simple forcing for Jupiter planets

    logical:: FlagRadiationDennouAGCM
                                ! flag for use of radiation model used in DennouAGCM
    logical:: FlagRadiationDcpamEV1
                                ! flag for use of radiation model used in DcpamEV1
    logical:: FlagRadiationDcpamMV1
                                ! flag for use of radiation model used in DcpamMV1
    logical:: FlagRadiationSL09
                                ! flag for use of radiation model used in DcpamMV1

    logical:: FlagPhyImp1LayModel
                             ! flag for use of implicit solver with surface 1 layer model
    logical:: FlagPhyImpSoilModel
                             ! flag for use of implicit solver with soil model
    logical:: FlagPhyImpAtmOnly
                             ! flag for use of implicit solver of atmospheric only

    character(STRING):: briefexpldyn
                              ! 実行ファイルの簡潔な説明 (力学過程)
                              ! Brief account of executable file (dynamics)
    character(STRING):: briefexplphy
                              ! 実行ファイルの簡潔な説明 (物理過程)
                              ! Brief account of executable file (physics)
    character(STRING):: briefexplrad
                              ! 実行ファイルの簡潔な説明 (放射過程)
                              ! Brief account of executable file (radiation)
    character(STRING):: briefexplimp
                              ! 実行ファイルの簡潔な説明 (陰解法方程式構築方法)
                              ! Brief account of executable file (implicit method)

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    integer:: k
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /dcpam_main_nml/ FlagDynamics, FlagFullPhysics, FlagHS94, FlagVenusSimple, FlagJupiterSimple, FlagRadiationDennouAGCM, FlagRadiationDcpamEV1, FlagRadiationDcpamMV1, FlagRadiationSL09, FlagPhyImp1LayModel, FlagPhyImpSoilModel, FlagPhyImpSoilModelSO, FlagPhyImpAtmOnly, FlagVerticalFilter
          !
          ! デフォルト値については初期化手続 "main/dcpam_main.F90#MainInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "main/dcpam_main.F90#MainInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    ! Initialize MPI
    !
    call MPIWrapperInit

    ! コマンドライン引数処理
    ! Command line option parser
    !
    call OptParseInit( opt_nml_filename, prog_name )

    ! NAMELIST ファイル名入力
    ! Input NAMELIST file name
    !
    call NmlutilInit( opt_nml_filename )

    ! デフォルト値の設定
    ! Default values settings
    !
    FlagDynamics            = .true.
    FlagFullPhysics         = .true.
    FlagHS94                = .false.
    FlagVenusSimple         = .false.
    FlagJupiterSimple       = .false.

    FlagRadiationDennouAGCM = .true.
    FlagRadiationDcpamEV1   = .false.
    FlagRadiationDcpamMV1   = .false.
    FlagRadiationSL09       = .false.
    FlagPhyImp1LayModel     = .true.
    FlagPhyImpSoilModel     = .false.
    FlagPhyImpSoilModelSO   = .false.
    FlagPhyImpAtmOnly       = .false.
    FlagVerticalFilter      = .false.

    ! 計算モードの設定
    ! Configure calculation mode
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = dcpam_main_nml, iostat = iostat_nml )    ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, prog_name ) ! (in)
      if ( iostat_nml == 0 ) write( STDOUT, nml = dcpam_main_nml )
    end if


    ! Identification of calculation method for dynamics
    !
    call MessageNotify( 'M', prog_name, 'FlagDynamics=<%b>.', l = (/FlagDynamics/) )

    if ( FlagDynamics ) then
      IDDynMode = IDDynModeHSPLVAS83
    else
      IDDynMode = IDDynModeNoHorAdv
    end if


    ! Identification of calculation method for physics
    !
    if ( FlagFullPhysics .and. FlagHS94 ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagHS94=<%b>.', l = (/FlagFullPhysics, FlagHS94/) )
    else if ( FlagFullPhysics .and. FlagVenusSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagVenusSimple=<%b>.', l = (/FlagFullPhysics, FlagVenusSimple/) )
    else if ( FlagFullPhysics .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagFullPhysics, FlagJupiterSimple/) )
    else if ( FlagHS94 .and. FlagVenusSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagHS94=<%b> conflicts with FlagVenusSimple=<%b>.', l = (/FlagHS94, FlagVenusSimple/) )
    else if ( FlagHS94 .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagHS94=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagHS94, FlagJupiterSimple/) )
    else if ( FlagVenusSimple .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagVenusSimple=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagVenusSimple, FlagJupiterSimple/) )
    end if
    !
    if ( FlagFullPhysics ) then
      IDPhyMode = IDPhyModeFullPhysics
    else if ( FlagHS94 ) then
      IDPhyMode = IDPhyModeHS94
    else if ( FlagVenusSimple ) then
      IDPhyMode = IDPhyModeVenusSimple
    else if ( FlagJupiterSimple ) then
      IDPhyMode = IDPhyModeJupiterSimple
    else
      IDPhyMode = IDPhyModeNoPhysics
    end if


    ! Identification of calculation method for radiation
    !
    call MessageNotify( 'M', prog_name, 'FlagRadiationDennouAGCM=<%b>.', l = (/FlagRadiationDennouAGCM/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationDcpamEV1  =<%b>.', l = (/FlagRadiationDcpamEV1/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationDcpamMV1  =<%b>.', l = (/FlagRadiationDcpamMV1/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationSL09      =<%b>.', l = (/FlagRadiationSL09/) )
    !
    !   Check for conflict
    !
    if ( FlagRadiationDennouAGCM .and. FlagRadiationDcpamEV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationDcpamEV1=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationDcpamEV1/) )
    else if ( FlagRadiationDennouAGCM .and. FlagRadiationDcpamMV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationDcpamMV1=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationDcpamMV1/) )
    else if ( FlagRadiationDennouAGCM .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationSL09/) )
    else if ( FlagRadiationDcpamEV1 .and. FlagRadiationDcpamMV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamEV1=<%b> conflicts with FlagRadiationDcpamMV1=<%b>.', l = (/FlagRadiationDcpamEV1, FlagRadiationDcpamMV1/) )
    else if ( FlagRadiationDcpamEV1 .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamEV1=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDcpamEV1, FlagRadiationSL09/) )
    else if ( FlagRadiationDcpamMV1 .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamMV1=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDcpamMV1, FlagRadiationSL09/) )
    end if
    !
    if ( FlagRadiationDennouAGCM ) then
      IDRadMethod = IDRadMethodDennouAGCM
    else if ( FlagRadiationDcpamEV1 ) then
      IDRadMethod = IDRadMethodDcpamEV1
    else if ( FlagRadiationDcpamMV1 ) then
      IDRadMethod = IDRadMethodDcpamMV1
    else if ( FlagRadiationSL09 ) then
      IDRadMethod = IDRadMethodSL09
    else
      if ( FlagFullPhysics ) then
        call MessageNotify( 'E', prog_name, 'None of radiation model is selected.' )
      end if
    end if


    ! Identification of calculation method for solving simultaneous linear equations 
    ! of physics
    !
    call MessageNotify( 'M', prog_name, 'FlagPhyImp1LayModel=<%b>.', l = (/FlagPhyImp1LayModel/) )
    call MessageNotify( 'M', prog_name, 'FlagPhyImpSoilModel=<%b>.', l = (/FlagPhyImpSoilModel/) )
    call MessageNotify( 'M', prog_name, 'FlagPhyImpAtmOnly  =<%b>.', l = (/FlagPhyImpAtmOnly/) )
    !
    !   Check for conflict
    !
    if ( FlagPhyImp1LayModel .and. FlagPhyImpSoilModel ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImp1LayModel=<%b> conflicts with FlagPhyImpSoilModel=<%b>.', l = (/FlagPhyImp1LayModel, FlagPhyImpSoilModel/) )
    else if ( FlagPhyImpSoilModel .and. FlagPhyImpAtmOnly ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpSoilModel=<%b> conflicts with FlagPhyImpAtmOnly=<%b>.', l = (/FlagPhyImpSoilModel, FlagPhyImpAtmOnly/) )
    else if ( FlagPhyImpAtmOnly .and. FlagPhyImp1LayModel ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpAtmOnly=<%b> conflicts with FlagPhyImp1LayModel=<%b>.', l = (/FlagPhyImpAtmOnly, FlagPhyImp1LayModel/) )
    end if
    !
    !   Check for value of FlagFullPhysics
    !
    if ( ( .not. FlagFullPhysics ) .and. FlagPhyImpSoilModel ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics has to be true, when FlagPhyImpSoilModel is true.' )
    end if
    if ( ( .not. FlagFullPhysics ) .and. FlagPhyImpAtmOnly ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics has to be true, when FlagPhyImpAtmOnly is true.' )
    end if
    !
    !   Check for slab ocean
    !
    if ( ( .not. FlagPhyImpSoilModel ) .and. FlagPhyImpSoilModelSO ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpSoilModel has to be true, when FlagPhyImpSoilModelSO is true.' )
    end if
    !
    if ( FlagPhyImp1LayModel ) then
      IDPhyTendMethod = IDPhyTendMethodImp1LayModel
    else if ( FlagPhyImpSoilModel ) then
      IDPhyTendMethod = IDPhyTendMethodImpSoilModel
    else if ( FlagPhyImpAtmOnly ) then
      IDPhyTendMethod = IDPhyTendMethodImpAtmOnly
    else
      if ( .not. FlagFullPhysics ) then
        call MessageNotify( 'E', prog_name, 'None of calculation method for simultaneous linear equations of physics is selected.' )
      end if
    end if


    ! 計算モードの表示
    ! Display calculation mode
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      briefexpldyn = 'used'
    case ( IDDynModeNoHorAdv )
      briefexpldyn = 'not used'
    end select

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      briefexplphy = 'parameterization suite is used'
    case ( IDPhyModeHS94 )
      briefexplphy = 'forcing for Held and Suarez (1994) dynamical core test'
    case ( IDPhyModeVenusSimple )
      briefexplphy = 'simple forcing for a Venus-like planet'
    case ( IDPhyModeJupiterSimple )
      briefexplphy = 'simple forcing for a Jupiter-like planet'
    case ( IDPhyModeNoPhysics )
      briefexplphy = 'not used'
    end select

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )
        briefexplrad = 'dennou AGCM5 default'
      case ( IDRadMethodDcpamEV1 )
        briefexplrad = 'Earth-like planet'
      case ( IDRadMethodDcpamMV1 )
        briefexplrad = 'Mars'
      case ( IDRadMethodSL09 )
        briefexplrad = 'Schneider and Liu (2009) (Jupiter-like planet)'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplrad', '' )
      end select

      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )
        briefexplimp = 'system with surface 1 layer model'
      case ( IDPhyTendMethodImpSoilModel )
        briefexplimp = 'system with thermal diffusion soil model'
      case ( IDPhyTendMethodImpAtmOnly )
        briefexplimp = 'system only with atmosphere'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplimp', '' )
      end select
    else
      briefexplrad = ''
    end if

    call MessageNotify( 'M', prog_name, '' )
    call MessageNotify( 'M', prog_name,   '+-------------------------------------' )
    call MessageNotify( 'M', prog_name,   '|  Dynamics: %c', c1 = trim(briefexpldyn) )
    call MessageNotify( 'M', prog_name,   '|  Physics : %c', c1 = trim(briefexplphy) )
    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      call MessageNotify( 'M', prog_name, '|    Radiation model : %c', c1 = trim(briefexplrad) )
      call MessageNotify( 'M', prog_name, '|    Implicit method : %c', c1 = trim(briefexplimp) )
    end if
    call MessageNotify( 'M', prog_name,   '|  -- version = %c', c1 = trim(version) )
    call MessageNotify( 'M', prog_name,   '+-------------------------------------' )
    call MessageNotify( 'M', prog_name, '' )

    ! 時刻管理
    ! Time control
    !
    call TimesetInit

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    call FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    call GridsetInit

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    call CompositionInit

    ! 物理定数設定
    ! Physical constants settings
    !
    call ConstantsInit

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    call ConstantsSnowSeaIceInit

    ! 座標データ設定
    ! Axes data settings
    !
    call AxessetInit

    ! 予報変数の割付
    ! Allocation of prediction variables
    !
    allocate( xyz_UB    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VB    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempB (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsB    (0:imax-1, 1:jmax) )

    allocate( xyz_UN    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VN    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempN (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsN    (0:imax-1, 1:jmax) )

    allocate( xyz_UA    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VA    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempA (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsA    (0:imax-1, 1:jmax) )


    ! リスタートデータ入力
    ! Restart data input
    !
    call RestartFileGet( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, flag_initial )                                  ! (out) optional


    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )

      ! 地表面温度, 土壌温度の割付
      ! Allocation of surface temperature and soil temperature
      !
      allocate( xy_SurfTemp (0:imax-1, 1:jmax) )
      allocate( xyz_SoilTemp(0:imax-1, 1:jmax, 1:kslmax) )

      allocate( xy_SoilMoistB(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowB (0:imax-1, 1:jmax) )

      allocate( xy_SoilMoistN(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowN (0:imax-1, 1:jmax) )

      allocate( xy_SoilMoistA(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowA (0:imax-1, 1:jmax) )


      ! 地表面温度リスタートデータ入力
      ! Restart data of surface temperature input
      !
      call RestartSurfTempGet( xy_SurfTemp  )          ! (out)


      ! THIS IS A TEMPORARY LINE.
      ! 土壌温度, ..., の初期値設定, いずれリスタートファイルから読むようにする
      ! Setting of initial values of soil temperature, ..., these values are input from restart file in near future
      !
      do k = 1, kslmax
        xyz_SoilTemp(:,:,k) = xy_SurfTemp
      end do

      xy_SoilMoistN = 0.0_DP
      xy_SurfSnowN  = 0.0_DP

      xy_SoilMoistB = xy_SoilMoistN
      xy_SurfSnowB  = xy_SurfSnowN

      xy_SoilMoistA = 0.0_DP
      xy_SurfSnowA  = 0.0_DP

    end select


    ! リスタートデータファイルの初期化
    ! Initialization of restart data file
    !
    call RestartFileOpen

    ! ヒストリデータファイルの初期化
    ! Initialization of history data files
    !
    call HistoryFileOpen

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'U' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'eastward wind', 'm s-1' )               ! (in)

    call HistoryAutoAddVariable( 'V' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'northward wind', 'm s-1' )              ! (in)

    call HistoryAutoAddVariable( 'Temp' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'temperature', 'K' )                     ! (in)

    do n = 1, ncmax
      call HistoryAutoAddVariable( a_QMixName(n) , (/ 'lon ', 'lat ', 'sig ', 'time' /), a_QMixLongName(n), 'kg kg-1' )             ! (in)
    end do

    call HistoryAutoAddVariable( 'Ps' , (/ 'lon ', 'lat ', 'time' /), 'surface pressure', 'Pa' )               ! (in)


    ! ヒストリデータ出力 (スタート時刻)
    ! History data output (Start time)
    !
    call HistoryAutoPut( TimeN, 'U', xyz_UN )
    call HistoryAutoPut( TimeN, 'V', xyz_VN )
    call HistoryAutoPut( TimeN, 'Temp', xyz_TempN )
    do n = 1, ncmax
      call HistoryAutoPut( TimeN, a_QMixName(n), xyzf_QMixN(:,:,:,n) )
    end do
    call HistoryAutoPut( TimeN, 'Ps', xy_PsN )


    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )

      ! 地表面温度リスタートデータファイルの初期化
      ! Initialization of restart data file of surface temperature
      !
      call RestartSurfTempOpen

      ! ヒストリデータ出力のためのへの変数登録
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'SurfTemp' , (/ 'lon ', 'lat ', 'time' /), 'surface temperature', 'K' )

      call HistoryAutoAddVariable( 'SoilTemp', (/ 'lon ', 'lat ', 'ssz ', 'time' /), 'soil temperature', 'K' )

      call HistoryAutoAddVariable( 'SoilMoist' , (/ 'lon ', 'lat ', 'time' /), 'soil moisture', 'kg m-2' )              ! (in)

      call HistoryAutoAddVariable( 'SurfSnow'  , (/ 'lon ', 'lat ', 'time' /), 'surface snow amount', 'kg m-2' )        ! (in)

      call HistoryAutoAddVariable( 'Rain', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'W m-2' )

      call HistoryAutoAddVariable( 'DTempDtCond' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'condensation heating', 'K s-1' )

      call HistoryAutoAddVariable( 'DQVapDtCond' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'condensation moistening', 'kg kg-1 s-1' )

    end select


    ! 診断変数の割付
    ! Allocation of diagnostic variables
    !
    allocate( xyz_DUDt    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_DVDt    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

    allocate( xy_SurfHeight (0:imax-1, 1:jmax) )

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      allocate( xy_SurfAlbedo         (0:imax-1, 1:jmax) )
      allocate( xy_SurfHumidCoef      (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLength    (0:imax-1, 1:jmax) )
      allocate( xy_SurfHeatCapacity   (0:imax-1, 1:jmax) )
      allocate( xy_SeaIceConc         (0:imax-1, 1:jmax) )
      allocate( xy_SurfCond           (0:imax-1, 1:jmax) )
      allocate( xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) )
      allocate( xy_SoilHeatCap        (0:imax-1, 1:jmax) )
      allocate( xy_SoilHeatDiffCoef   (0:imax-1, 1:jmax) )

      allocate( xyr_Temp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Height (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Height (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Exner  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Exner  (0:imax-1, 1:jmax, 0:kmax) )

      allocate( xyr_RadLFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadLFluxA    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax,  0:1) )

      allocate( xyr_MomFluxX  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_MomFluxY  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_HeatFlux  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyrf_QMixFlux (0:imax-1, 1:jmax, 0:kmax, 1:ncmax) )

      allocate( xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) )

      allocate( xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) )

      allocate( xy_SurfVelTransCoef (0:imax-1, 1:jmax) )
      allocate( xy_SurfTempTransCoef(0:imax-1, 1:jmax) )
      allocate( xy_SurfQVapTransCoef(0:imax-1, 1:jmax) )

      allocate( xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax) )

      allocate( xy_DSurfTempDt (0:imax-1, 1:jmax) )
      allocate( xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) )

      allocate( xy_DSoilMoistDt(0:imax-1, 1:jmax) )
      allocate( xy_DSurfSnowDt (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax) )


      allocate( xy_Rain         (0:imax-1, 1:jmax) )
      allocate( xyz_DTempDtCond (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DQVapDtCond (0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyz_DDelLWDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DDelLWDtLSC (0:imax-1,1:jmax,1:kmax) )

      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtLSC (0:imax-1,1:jmax,1:kmax) )


      ! ヒストリデータ出力のためのへの変数登録
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'SeaIceConc' , (/ 'lon ', 'lat ', 'time' /), 'sea ice concentration', '1' )

      call HistoryAutoAddVariable( 'SurfAlbedo' , (/ 'lon ', 'lat ', 'time' /), 'surface albedo', '1' )

    case ( IDPhyModeVenusSimple )

      allocate( xyr_Temp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Height (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Height (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Exner  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Exner  (0:imax-1, 1:jmax, 0:kmax) )

    case ( IDPhyModeJupiterSimple )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )

      allocate( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyr_HeatFlux  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadLFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax,  0:1) )

    end select

    ! 初回だけはオイラー法を用いるため, Δt を半分に
    ! Delta t is reduced to half in order to use Euler method at initial step
    !
    if ( flag_initial ) then
      call TimesetDelTimeHalf
    end if

  end subroutine MainInit

  !-------------------------------------------------------------------

  subroutine MainTerminate
    !
    ! 主プログラムの終了処理手続き. 
    !
    ! Termination procedure for the main program. 
    !

    ! モジュール引用 ; USE statements
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperFinalize

    ! 力学過程 (スペクトル法, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only: DynamicsFinalize

    ! Held and Suarez (1994) による強制と散逸
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only: Hs94Finalize

    ! 放射フラックス (バンドモデル)
    ! Radiation flux (band model)
    !
    use radiation_DennouAGCM, only: RadiationDennouAGCMFinalize

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetFinalize

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxFinalize

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClose

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileClose

    ! 地表面温度リスタートデータ入出力
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempClose

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileClose

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

    ! リスタートデータファイルクローズ
    ! Close restart data file
    !
    call RestartFileClose

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      ! 地表面温度リスタートデータファイルクローズ
      ! Close restart data file of surface temperature
      !
      call RestartSurfTempClose
    end if

    ! ヒストリデータファイルクローズ
    ! Close history data files
    !
    call HistoryFileClose

    ! 予報変数の割付解除
    ! Deallocation of prediction variables
    !
    deallocate( xyz_UB     )
    deallocate( xyz_VB     )
    deallocate( xyz_TempB  )
    deallocate( xyzf_QMixB )
    deallocate( xy_PsB     )

    deallocate( xyz_UN     )
    deallocate( xyz_VN     )
    deallocate( xyz_TempN  )
    deallocate( xyzf_QMixN )
    deallocate( xy_PsN     )

    deallocate( xyz_UA     )
    deallocate( xyz_VA     )
    deallocate( xyz_TempA  )
    deallocate( xyzf_QMixA )
    deallocate( xy_PsA     )

    ! 診断変数の割付解除
    ! Dellocation of diagnostic variables
    !
    deallocate( xyz_DUDt     )
    deallocate( xyz_DVDt     )
    deallocate( xyz_DTempDt  )
    deallocate( xyzf_DQMixDt )

    deallocate( xy_SurfHeight )

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      deallocate( xy_SurfAlbedo          )
      deallocate( xy_SurfHumidCoef       )
      deallocate( xy_SurfRoughLength     )
      deallocate( xy_SurfHeatCapacity    )
      deallocate( xy_SeaIceConc          )
      deallocate( xy_SurfCond            )
      deallocate( xy_DeepSubSurfHeatFlux )
      deallocate( xy_SoilHeatCap         )
      deallocate( xy_SoilHeatDiffCoef    )

      deallocate( xyr_Temp   )
      deallocate( xyz_Press  )
      deallocate( xyr_Press  )
      deallocate( xyz_Height )
      deallocate( xyr_Height )
      deallocate( xyz_Exner  )
      deallocate( xyr_Exner  )

      deallocate( xyr_RadLFlux     )
      deallocate( xyr_RadLFluxA    )
      deallocate( xyr_RadSFlux     )
      deallocate( xyra_DelRadLFlux )

      deallocate( xyr_MomFluxX  )
      deallocate( xyr_MomFluxY  )
      deallocate( xyr_HeatFlux  )
      deallocate( xyrf_QMixFlux )

      deallocate( xyr_VelTransCoef  )
      deallocate( xyr_TempTransCoef )
      deallocate( xyr_QMixTransCoef )

      deallocate( xy_SurfVelTransCoef  )
      deallocate( xy_SurfTempTransCoef )
      deallocate( xy_SurfQVapTransCoef )

      deallocate( xy_DSurfTempDt )
      deallocate( xyz_DSoilTempDt )

      deallocate( xy_DSoilMoistDt )
      deallocate( xy_DSurfSnowDt )

      deallocate( xyz_DTempDtVDiff )

      deallocate( xyz_DTempDtRadL )
      deallocate( xyz_DTempDtRadS )


      deallocate( xy_Rain          )
      deallocate( xyz_DTempDtCond  )
      deallocate( xyz_DQVapDtCond  )

      deallocate( xyz_DDelLWDtCum  )
      deallocate( xyz_DDelLWDtLSC  )

      deallocate( xyz_DQH2OLiqDtCum  )
      deallocate( xyz_DQH2OLiqDtLSC  )
    end if ! FlagFullPhysics



    ! 各モジュール内の変数の割付解除
    ! Dellocation of variables in modules
    !
    call DynamicsFinalize
    call AuxFinalize

    if ( IDPhyMode == IDPhyModeHS94 ) then
      call Hs94Finalize
    end if

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      ! 割付解除とリスタートファイルの終了処理
      ! Dellocation and close a restart file
      !
      call RadiationDennouAGCMFinalize
    end if

    call AxessetFinalize

    ! 時刻管理終了処理
    ! Termination of time control
    !
    call TimesetClose

    ! Finalize MPI
    !
    call MPIWrapperFinalize


  end subroutine MainTerminate

end program dcpam_main

  subroutine MainInit
    !
    ! 主プログラムの初期化手続き. 
    !
    ! Initialization procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperInit

    use dc_message, only: MessageNotify

    ! コマンドライン引数処理
    ! Command line option parser
    !
    use option_parser, only: OptParseInit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: NmlutilInit, NmlutilMsg, namelist_filename

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetInit, TimesetDelTimeHalf, TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    use fileset, only: FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: GridsetInit, imax, jmax, kmax, kslmax    ! 地下の鉛直層数. 
                                 ! Number of subsurface vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: CompositionInit

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: ConstantsInit

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    use constants_snowseaice, only: ConstantsSnowSeaIceInit

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetInit

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileOpen, RestartFileGet

    ! 地表面温度リスタートデータ入出力
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempOpen, RestartSurfTempGet

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileOpen
    use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT               ! 標準出力の装置番号. Unit number of standard output

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 宣言文 ; Declaration statements
    !
    implicit none
    character(*), parameter:: prog_name = 'dcpam_main'
                              ! 主プログラム名. 
                              ! Main program name
    character(*), parameter:: version = '$Name: dcpam5-20110221-2 $' // '$Id: dcpam_main.f90,v 1.18 2011-02-18 04:45:17 yot Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version

    character(STRING)      :: opt_nml_filename
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name

    logical:: FlagDynamics      ! 力学過程計算用フラグ
                                ! Flag for dynamics

    logical:: FlagFullPhysics   ! 全物理過程計算用フラグ
                                ! Flag for full physics
    logical:: FlagHS94          ! Held and Suarez (1994) 強制オン／オフ. 
                                ! Held and Suarez (1994) forcing on/off. 
    logical:: FlagVenusSimple   ! 金星簡単強制用フラグ
                                ! Flag for simple forcing for Venus
    logical:: FlagJupiterSimple !
                                ! Flag for simple forcing for Jupiter planets

    logical:: FlagRadiationDennouAGCM
                                ! flag for use of radiation model used in DennouAGCM
    logical:: FlagRadiationDcpamEV1
                                ! flag for use of radiation model used in DcpamEV1
    logical:: FlagRadiationDcpamMV1
                                ! flag for use of radiation model used in DcpamMV1
    logical:: FlagRadiationSL09
                                ! flag for use of radiation model used in DcpamMV1

    logical:: FlagPhyImp1LayModel
                             ! flag for use of implicit solver with surface 1 layer model
    logical:: FlagPhyImpSoilModel
                             ! flag for use of implicit solver with soil model
    logical:: FlagPhyImpAtmOnly
                             ! flag for use of implicit solver of atmospheric only

    character(STRING):: briefexpldyn
                              ! 実行ファイルの簡潔な説明 (力学過程)
                              ! Brief account of executable file (dynamics)
    character(STRING):: briefexplphy
                              ! 実行ファイルの簡潔な説明 (物理過程)
                              ! Brief account of executable file (physics)
    character(STRING):: briefexplrad
                              ! 実行ファイルの簡潔な説明 (放射過程)
                              ! Brief account of executable file (radiation)
    character(STRING):: briefexplimp
                              ! 実行ファイルの簡潔な説明 (陰解法方程式構築方法)
                              ! Brief account of executable file (implicit method)

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    integer:: k
    integer:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /dcpam_main_nml/ FlagDynamics, FlagFullPhysics, FlagHS94, FlagVenusSimple, FlagJupiterSimple, FlagRadiationDennouAGCM, FlagRadiationDcpamEV1, FlagRadiationDcpamMV1, FlagRadiationSL09, FlagPhyImp1LayModel, FlagPhyImpSoilModel, FlagPhyImpSoilModelSO, FlagPhyImpAtmOnly, FlagVerticalFilter
          !
          ! デフォルト値については初期化手続 "main/dcpam_main.F90#MainInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "main/dcpam_main.F90#MainInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    ! Initialize MPI
    !
    call MPIWrapperInit

    ! コマンドライン引数処理
    ! Command line option parser
    !
    call OptParseInit( opt_nml_filename, prog_name )

    ! NAMELIST ファイル名入力
    ! Input NAMELIST file name
    !
    call NmlutilInit( opt_nml_filename )

    ! デフォルト値の設定
    ! Default values settings
    !
    FlagDynamics            = .true.
    FlagFullPhysics         = .true.
    FlagHS94                = .false.
    FlagVenusSimple         = .false.
    FlagJupiterSimple       = .false.

    FlagRadiationDennouAGCM = .true.
    FlagRadiationDcpamEV1   = .false.
    FlagRadiationDcpamMV1   = .false.
    FlagRadiationSL09       = .false.
    FlagPhyImp1LayModel     = .true.
    FlagPhyImpSoilModel     = .false.
    FlagPhyImpSoilModelSO   = .false.
    FlagPhyImpAtmOnly       = .false.
    FlagVerticalFilter      = .false.

    ! 計算モードの設定
    ! Configure calculation mode
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = dcpam_main_nml, iostat = iostat_nml )    ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, prog_name ) ! (in)
      if ( iostat_nml == 0 ) write( STDOUT, nml = dcpam_main_nml )
    end if


    ! Identification of calculation method for dynamics
    !
    call MessageNotify( 'M', prog_name, 'FlagDynamics=<%b>.', l = (/FlagDynamics/) )

    if ( FlagDynamics ) then
      IDDynMode = IDDynModeHSPLVAS83
    else
      IDDynMode = IDDynModeNoHorAdv
    end if


    ! Identification of calculation method for physics
    !
    if ( FlagFullPhysics .and. FlagHS94 ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagHS94=<%b>.', l = (/FlagFullPhysics, FlagHS94/) )
    else if ( FlagFullPhysics .and. FlagVenusSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagVenusSimple=<%b>.', l = (/FlagFullPhysics, FlagVenusSimple/) )
    else if ( FlagFullPhysics .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagFullPhysics, FlagJupiterSimple/) )
    else if ( FlagHS94 .and. FlagVenusSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagHS94=<%b> conflicts with FlagVenusSimple=<%b>.', l = (/FlagHS94, FlagVenusSimple/) )
    else if ( FlagHS94 .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagHS94=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagHS94, FlagJupiterSimple/) )
    else if ( FlagVenusSimple .and. FlagJupiterSimple ) then
      call MessageNotify( 'E', prog_name, 'FlagVenusSimple=<%b> conflicts with FlagJupiterSimple=<%b>.', l = (/FlagVenusSimple, FlagJupiterSimple/) )
    end if
    !
    if ( FlagFullPhysics ) then
      IDPhyMode = IDPhyModeFullPhysics
    else if ( FlagHS94 ) then
      IDPhyMode = IDPhyModeHS94
    else if ( FlagVenusSimple ) then
      IDPhyMode = IDPhyModeVenusSimple
    else if ( FlagJupiterSimple ) then
      IDPhyMode = IDPhyModeJupiterSimple
    else
      IDPhyMode = IDPhyModeNoPhysics
    end if


    ! Identification of calculation method for radiation
    !
    call MessageNotify( 'M', prog_name, 'FlagRadiationDennouAGCM=<%b>.', l = (/FlagRadiationDennouAGCM/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationDcpamEV1  =<%b>.', l = (/FlagRadiationDcpamEV1/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationDcpamMV1  =<%b>.', l = (/FlagRadiationDcpamMV1/) )
    call MessageNotify( 'M', prog_name, 'FlagRadiationSL09      =<%b>.', l = (/FlagRadiationSL09/) )
    !
    !   Check for conflict
    !
    if ( FlagRadiationDennouAGCM .and. FlagRadiationDcpamEV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationDcpamEV1=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationDcpamEV1/) )
    else if ( FlagRadiationDennouAGCM .and. FlagRadiationDcpamMV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationDcpamMV1=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationDcpamMV1/) )
    else if ( FlagRadiationDennouAGCM .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDennouAGCM=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDennouAGCM, FlagRadiationSL09/) )
    else if ( FlagRadiationDcpamEV1 .and. FlagRadiationDcpamMV1 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamEV1=<%b> conflicts with FlagRadiationDcpamMV1=<%b>.', l = (/FlagRadiationDcpamEV1, FlagRadiationDcpamMV1/) )
    else if ( FlagRadiationDcpamEV1 .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamEV1=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDcpamEV1, FlagRadiationSL09/) )
    else if ( FlagRadiationDcpamMV1 .and. FlagRadiationSL09 ) then
      call MessageNotify( 'E', prog_name, 'FlagRadiationDcpamMV1=<%b> conflicts with FlagRadiationSL09=<%b>.', l = (/FlagRadiationDcpamMV1, FlagRadiationSL09/) )
    end if
    !
    if ( FlagRadiationDennouAGCM ) then
      IDRadMethod = IDRadMethodDennouAGCM
    else if ( FlagRadiationDcpamEV1 ) then
      IDRadMethod = IDRadMethodDcpamEV1
    else if ( FlagRadiationDcpamMV1 ) then
      IDRadMethod = IDRadMethodDcpamMV1
    else if ( FlagRadiationSL09 ) then
      IDRadMethod = IDRadMethodSL09
    else
      if ( FlagFullPhysics ) then
        call MessageNotify( 'E', prog_name, 'None of radiation model is selected.' )
      end if
    end if


    ! Identification of calculation method for solving simultaneous linear equations 
    ! of physics
    !
    call MessageNotify( 'M', prog_name, 'FlagPhyImp1LayModel=<%b>.', l = (/FlagPhyImp1LayModel/) )
    call MessageNotify( 'M', prog_name, 'FlagPhyImpSoilModel=<%b>.', l = (/FlagPhyImpSoilModel/) )
    call MessageNotify( 'M', prog_name, 'FlagPhyImpAtmOnly  =<%b>.', l = (/FlagPhyImpAtmOnly/) )
    !
    !   Check for conflict
    !
    if ( FlagPhyImp1LayModel .and. FlagPhyImpSoilModel ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImp1LayModel=<%b> conflicts with FlagPhyImpSoilModel=<%b>.', l = (/FlagPhyImp1LayModel, FlagPhyImpSoilModel/) )
    else if ( FlagPhyImpSoilModel .and. FlagPhyImpAtmOnly ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpSoilModel=<%b> conflicts with FlagPhyImpAtmOnly=<%b>.', l = (/FlagPhyImpSoilModel, FlagPhyImpAtmOnly/) )
    else if ( FlagPhyImpAtmOnly .and. FlagPhyImp1LayModel ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpAtmOnly=<%b> conflicts with FlagPhyImp1LayModel=<%b>.', l = (/FlagPhyImpAtmOnly, FlagPhyImp1LayModel/) )
    end if
    !
    !   Check for value of FlagFullPhysics
    !
    if ( ( .not. FlagFullPhysics ) .and. FlagPhyImpSoilModel ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics has to be true, when FlagPhyImpSoilModel is true.' )
    end if
    if ( ( .not. FlagFullPhysics ) .and. FlagPhyImpAtmOnly ) then
      call MessageNotify( 'E', prog_name, 'FlagFullPhysics has to be true, when FlagPhyImpAtmOnly is true.' )
    end if
    !
    !   Check for slab ocean
    !
    if ( ( .not. FlagPhyImpSoilModel ) .and. FlagPhyImpSoilModelSO ) then
      call MessageNotify( 'E', prog_name, 'FlagPhyImpSoilModel has to be true, when FlagPhyImpSoilModelSO is true.' )
    end if
    !
    if ( FlagPhyImp1LayModel ) then
      IDPhyTendMethod = IDPhyTendMethodImp1LayModel
    else if ( FlagPhyImpSoilModel ) then
      IDPhyTendMethod = IDPhyTendMethodImpSoilModel
    else if ( FlagPhyImpAtmOnly ) then
      IDPhyTendMethod = IDPhyTendMethodImpAtmOnly
    else
      if ( .not. FlagFullPhysics ) then
        call MessageNotify( 'E', prog_name, 'None of calculation method for simultaneous linear equations of physics is selected.' )
      end if
    end if


    ! 計算モードの表示
    ! Display calculation mode
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      briefexpldyn = 'used'
    case ( IDDynModeNoHorAdv )
      briefexpldyn = 'not used'
    end select

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      briefexplphy = 'parameterization suite is used'
    case ( IDPhyModeHS94 )
      briefexplphy = 'forcing for Held and Suarez (1994) dynamical core test'
    case ( IDPhyModeVenusSimple )
      briefexplphy = 'simple forcing for a Venus-like planet'
    case ( IDPhyModeJupiterSimple )
      briefexplphy = 'simple forcing for a Jupiter-like planet'
    case ( IDPhyModeNoPhysics )
      briefexplphy = 'not used'
    end select

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )
        briefexplrad = 'dennou AGCM5 default'
      case ( IDRadMethodDcpamEV1 )
        briefexplrad = 'Earth-like planet'
      case ( IDRadMethodDcpamMV1 )
        briefexplrad = 'Mars'
      case ( IDRadMethodSL09 )
        briefexplrad = 'Schneider and Liu (2009) (Jupiter-like planet)'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplrad', '' )
      end select

      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )
        briefexplimp = 'system with surface 1 layer model'
      case ( IDPhyTendMethodImpSoilModel )
        briefexplimp = 'system with thermal diffusion soil model'
      case ( IDPhyTendMethodImpAtmOnly )
        briefexplimp = 'system only with atmosphere'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplimp', '' )
      end select
    else
      briefexplrad = ''
    end if

    call MessageNotify( 'M', prog_name, '' )
    call MessageNotify( 'M', prog_name,   '+-------------------------------------' )
    call MessageNotify( 'M', prog_name,   '|  Dynamics: %c', c1 = trim(briefexpldyn) )
    call MessageNotify( 'M', prog_name,   '|  Physics : %c', c1 = trim(briefexplphy) )
    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      call MessageNotify( 'M', prog_name, '|    Radiation model : %c', c1 = trim(briefexplrad) )
      call MessageNotify( 'M', prog_name, '|    Implicit method : %c', c1 = trim(briefexplimp) )
    end if
    call MessageNotify( 'M', prog_name,   '|  -- version = %c', c1 = trim(version) )
    call MessageNotify( 'M', prog_name,   '+-------------------------------------' )
    call MessageNotify( 'M', prog_name, '' )

    ! 時刻管理
    ! Time control
    !
    call TimesetInit

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    call FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    call GridsetInit

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    call CompositionInit

    ! 物理定数設定
    ! Physical constants settings
    !
    call ConstantsInit

    ! 雪と海氷の定数の設定
    ! Setting constants of snow and sea ice
    !
    call ConstantsSnowSeaIceInit

    ! 座標データ設定
    ! Axes data settings
    !
    call AxessetInit

    ! 予報変数の割付
    ! Allocation of prediction variables
    !
    allocate( xyz_UB    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VB    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempB (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixB(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsB    (0:imax-1, 1:jmax) )

    allocate( xyz_UN    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VN    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempN (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixN(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsN    (0:imax-1, 1:jmax) )

    allocate( xyz_UA    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_VA    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_TempA (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_QMixA(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )
    allocate( xy_PsA    (0:imax-1, 1:jmax) )


    ! リスタートデータ入力
    ! Restart data input
    !
    call RestartFileGet( xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, flag_initial )                                  ! (out) optional


    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )

      ! 地表面温度, 土壌温度の割付
      ! Allocation of surface temperature and soil temperature
      !
      allocate( xy_SurfTemp (0:imax-1, 1:jmax) )
      allocate( xyz_SoilTemp(0:imax-1, 1:jmax, 1:kslmax) )

      allocate( xy_SoilMoistB(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowB (0:imax-1, 1:jmax) )

      allocate( xy_SoilMoistN(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowN (0:imax-1, 1:jmax) )

      allocate( xy_SoilMoistA(0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowA (0:imax-1, 1:jmax) )


      ! 地表面温度リスタートデータ入力
      ! Restart data of surface temperature input
      !
      call RestartSurfTempGet( xy_SurfTemp  )          ! (out)


      ! THIS IS A TEMPORARY LINE.
      ! 土壌温度, ..., の初期値設定, いずれリスタートファイルから読むようにする
      ! Setting of initial values of soil temperature, ..., these values are input from restart file in near future
      !
      do k = 1, kslmax
        xyz_SoilTemp(:,:,k) = xy_SurfTemp
      end do

      xy_SoilMoistN = 0.0_DP
      xy_SurfSnowN  = 0.0_DP

      xy_SoilMoistB = xy_SoilMoistN
      xy_SurfSnowB  = xy_SurfSnowN

      xy_SoilMoistA = 0.0_DP
      xy_SurfSnowA  = 0.0_DP

    end select


    ! リスタートデータファイルの初期化
    ! Initialization of restart data file
    !
    call RestartFileOpen

    ! ヒストリデータファイルの初期化
    ! Initialization of history data files
    !
    call HistoryFileOpen

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'U' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'eastward wind', 'm s-1' )               ! (in)

    call HistoryAutoAddVariable( 'V' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'northward wind', 'm s-1' )              ! (in)

    call HistoryAutoAddVariable( 'Temp' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'temperature', 'K' )                     ! (in)

    do n = 1, ncmax
      call HistoryAutoAddVariable( a_QMixName(n) , (/ 'lon ', 'lat ', 'sig ', 'time' /), a_QMixLongName(n), 'kg kg-1' )             ! (in)
    end do

    call HistoryAutoAddVariable( 'Ps' , (/ 'lon ', 'lat ', 'time' /), 'surface pressure', 'Pa' )               ! (in)


    ! ヒストリデータ出力 (スタート時刻)
    ! History data output (Start time)
    !
    call HistoryAutoPut( TimeN, 'U', xyz_UN )
    call HistoryAutoPut( TimeN, 'V', xyz_VN )
    call HistoryAutoPut( TimeN, 'Temp', xyz_TempN )
    do n = 1, ncmax
      call HistoryAutoPut( TimeN, a_QMixName(n), xyzf_QMixN(:,:,:,n) )
    end do
    call HistoryAutoPut( TimeN, 'Ps', xy_PsN )


    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )

      ! 地表面温度リスタートデータファイルの初期化
      ! Initialization of restart data file of surface temperature
      !
      call RestartSurfTempOpen

      ! ヒストリデータ出力のためのへの変数登録
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'SurfTemp' , (/ 'lon ', 'lat ', 'time' /), 'surface temperature', 'K' )

      call HistoryAutoAddVariable( 'SoilTemp', (/ 'lon ', 'lat ', 'ssz ', 'time' /), 'soil temperature', 'K' )

      call HistoryAutoAddVariable( 'SoilMoist' , (/ 'lon ', 'lat ', 'time' /), 'soil moisture', 'kg m-2' )              ! (in)

      call HistoryAutoAddVariable( 'SurfSnow'  , (/ 'lon ', 'lat ', 'time' /), 'surface snow amount', 'kg m-2' )        ! (in)

      call HistoryAutoAddVariable( 'Rain', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'W m-2' )

      call HistoryAutoAddVariable( 'DTempDtCond' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'condensation heating', 'K s-1' )

      call HistoryAutoAddVariable( 'DQVapDtCond' , (/ 'lon ', 'lat ', 'sig ', 'time' /), 'condensation moistening', 'kg kg-1 s-1' )

    end select


    ! 診断変数の割付
    ! Allocation of diagnostic variables
    !
    allocate( xyz_DUDt    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_DVDt    (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_DTempDt (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

    allocate( xy_SurfHeight (0:imax-1, 1:jmax) )

    select case ( IDPhyMode )
    case ( IDPhyModeFullPhysics )
      allocate( xy_SurfAlbedo         (0:imax-1, 1:jmax) )
      allocate( xy_SurfHumidCoef      (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLength    (0:imax-1, 1:jmax) )
      allocate( xy_SurfHeatCapacity   (0:imax-1, 1:jmax) )
      allocate( xy_SeaIceConc         (0:imax-1, 1:jmax) )
      allocate( xy_SurfCond           (0:imax-1, 1:jmax) )
      allocate( xy_DeepSubSurfHeatFlux(0:imax-1, 1:jmax) )
      allocate( xy_SoilHeatCap        (0:imax-1, 1:jmax) )
      allocate( xy_SoilHeatDiffCoef   (0:imax-1, 1:jmax) )

      allocate( xyr_Temp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Height (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Height (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Exner  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Exner  (0:imax-1, 1:jmax, 0:kmax) )

      allocate( xyr_RadLFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadLFluxA    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax,  0:1) )

      allocate( xyr_MomFluxX  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_MomFluxY  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_HeatFlux  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyrf_QMixFlux (0:imax-1, 1:jmax, 0:kmax, 1:ncmax) )

      allocate( xyr_SoilHeatFlux(0:imax-1, 1:jmax, 0:kslmax) )

      allocate( xyr_VelTransCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempTransCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixTransCoef(0:imax-1, 1:jmax, 0:kmax) )

      allocate( xy_SurfVelTransCoef (0:imax-1, 1:jmax) )
      allocate( xy_SurfTempTransCoef(0:imax-1, 1:jmax) )
      allocate( xy_SurfQVapTransCoef(0:imax-1, 1:jmax) )

      allocate( xyr_SoilTempTransCoef (0:imax-1, 1:jmax, 0:kslmax) )

      allocate( xy_DSurfTempDt (0:imax-1, 1:jmax) )
      allocate( xyz_DSoilTempDt(0:imax-1, 1:jmax, 1:kslmax) )

      allocate( xy_DSoilMoistDt(0:imax-1, 1:jmax) )
      allocate( xy_DSurfSnowDt (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax) )


      allocate( xy_Rain         (0:imax-1, 1:jmax) )
      allocate( xyz_DTempDtCond (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DQVapDtCond (0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyz_DDelLWDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DDelLWDtLSC (0:imax-1,1:jmax,1:kmax) )

      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtLSC (0:imax-1,1:jmax,1:kmax) )


      ! ヒストリデータ出力のためのへの変数登録
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'SeaIceConc' , (/ 'lon ', 'lat ', 'time' /), 'sea ice concentration', '1' )

      call HistoryAutoAddVariable( 'SurfAlbedo' , (/ 'lon ', 'lat ', 'time' /), 'surface albedo', '1' )

    case ( IDPhyModeVenusSimple )

      allocate( xyr_Temp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Height (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Height (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_Exner  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Exner  (0:imax-1, 1:jmax, 0:kmax) )

    case ( IDPhyModeJupiterSimple )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )

      allocate( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyr_HeatFlux  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadLFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSFlux     (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax,  0:1) )

    end select

    ! 初回だけはオイラー法を用いるため, Δt を半分に
    ! Delta t is reduced to half in order to use Euler method at initial step
    !
    if ( flag_initial ) then
      call TimesetDelTimeHalf
    end if

  end subroutine MainInit

  subroutine MainTerminate
    !
    ! 主プログラムの終了処理手続き. 
    !
    ! Termination procedure for the main program. 
    !

    ! モジュール引用 ; USE statements
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperFinalize

    ! 力学過程 (スペクトル法, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only: DynamicsFinalize

    ! Held and Suarez (1994) による強制と散逸
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only: Hs94Finalize

    ! 放射フラックス (バンドモデル)
    ! Radiation flux (band model)
    !
    use radiation_DennouAGCM, only: RadiationDennouAGCMFinalize

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetFinalize

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level, 
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxFinalize

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClose

    ! リスタートデータ入出力
    ! Restart data input/output
    !
    use restart_file_io, only: RestartFileClose

    ! 地表面温度リスタートデータ入出力
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempClose

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileClose

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

    ! リスタートデータファイルクローズ
    ! Close restart data file
    !
    call RestartFileClose

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      ! 地表面温度リスタートデータファイルクローズ
      ! Close restart data file of surface temperature
      !
      call RestartSurfTempClose
    end if

    ! ヒストリデータファイルクローズ
    ! Close history data files
    !
    call HistoryFileClose

    ! 予報変数の割付解除
    ! Deallocation of prediction variables
    !
    deallocate( xyz_UB     )
    deallocate( xyz_VB     )
    deallocate( xyz_TempB  )
    deallocate( xyzf_QMixB )
    deallocate( xy_PsB     )

    deallocate( xyz_UN     )
    deallocate( xyz_VN     )
    deallocate( xyz_TempN  )
    deallocate( xyzf_QMixN )
    deallocate( xy_PsN     )

    deallocate( xyz_UA     )
    deallocate( xyz_VA     )
    deallocate( xyz_TempA  )
    deallocate( xyzf_QMixA )
    deallocate( xy_PsA     )

    ! 診断変数の割付解除
    ! Dellocation of diagnostic variables
    !
    deallocate( xyz_DUDt     )
    deallocate( xyz_DVDt     )
    deallocate( xyz_DTempDt  )
    deallocate( xyzf_DQMixDt )

    deallocate( xy_SurfHeight )

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      deallocate( xy_SurfAlbedo          )
      deallocate( xy_SurfHumidCoef       )
      deallocate( xy_SurfRoughLength     )
      deallocate( xy_SurfHeatCapacity    )
      deallocate( xy_SeaIceConc          )
      deallocate( xy_SurfCond            )
      deallocate( xy_DeepSubSurfHeatFlux )
      deallocate( xy_SoilHeatCap         )
      deallocate( xy_SoilHeatDiffCoef    )

      deallocate( xyr_Temp   )
      deallocate( xyz_Press  )
      deallocate( xyr_Press  )
      deallocate( xyz_Height )
      deallocate( xyr_Height )
      deallocate( xyz_Exner  )
      deallocate( xyr_Exner  )

      deallocate( xyr_RadLFlux     )
      deallocate( xyr_RadLFluxA    )
      deallocate( xyr_RadSFlux     )
      deallocate( xyra_DelRadLFlux )

      deallocate( xyr_MomFluxX  )
      deallocate( xyr_MomFluxY  )
      deallocate( xyr_HeatFlux  )
      deallocate( xyrf_QMixFlux )

      deallocate( xyr_VelTransCoef  )
      deallocate( xyr_TempTransCoef )
      deallocate( xyr_QMixTransCoef )

      deallocate( xy_SurfVelTransCoef  )
      deallocate( xy_SurfTempTransCoef )
      deallocate( xy_SurfQVapTransCoef )

      deallocate( xy_DSurfTempDt )
      deallocate( xyz_DSoilTempDt )

      deallocate( xy_DSoilMoistDt )
      deallocate( xy_DSurfSnowDt )

      deallocate( xyz_DTempDtVDiff )

      deallocate( xyz_DTempDtRadL )
      deallocate( xyz_DTempDtRadS )


      deallocate( xy_Rain          )
      deallocate( xyz_DTempDtCond  )
      deallocate( xyz_DQVapDtCond  )

      deallocate( xyz_DDelLWDtCum  )
      deallocate( xyz_DDelLWDtLSC  )

      deallocate( xyz_DQH2OLiqDtCum  )
      deallocate( xyz_DQH2OLiqDtLSC  )
    end if ! FlagFullPhysics



    ! 各モジュール内の変数の割付解除
    ! Dellocation of variables in modules
    !
    call DynamicsFinalize
    call AuxFinalize

    if ( IDPhyMode == IDPhyModeHS94 ) then
      call Hs94Finalize
    end if

    if ( IDPhyMode == IDPhyModeFullPhysics ) then
      ! 割付解除とリスタートファイルの終了処理
      ! Dellocation and close a restart file
      !
      call RadiationDennouAGCMFinalize
    end if

    call AxessetFinalize

    ! 時刻管理終了処理
    ! Termination of time control
    !
    call TimesetClose

    ! Finalize MPI
    !
    call MPIWrapperFinalize


  end subroutine MainTerminate