!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi
! Version::   $Id: dcpam_main_gabls2.f90,v 1.1 2013/09/21 14:58:39 yot Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2010. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

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

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

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

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

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

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

  ! 陰解法のための行列処理 (一部の物理過程用)
  ! 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:              &
    & PhyImplSDHSetMethodFromMatthews,     &
    & PhyImplSDHTendency,                  &
    & PhyImplSDHCorSOTempBySnowMelt

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

  ! 時間フィルター (Williams, 2009)
  ! Time filter (Williams, 2009)
  !
  use timefilter_williams2009, only: &
    & TimeFilterWilliams2009, TimeFilterWilliams2009SurfVars

  ! 時刻管理
  ! Time control
  !
  use timeset, only: TimesetProgress, &
    & TimeB, &                ! ステップ $ t - \Delta t $ の時刻. 
                              ! Time of step $ t - \Delta t $. 
    & TimeN, &                ! ステップ $ t $ の時刻. 
                              ! Time of step $ t $. 
    & TimeA, &                ! ステップ $ t + \Delta t $ の時刻. 
                              ! Time of step $ t + \Delta t $. 
    & EndTime, &              ! 計算終了時刻. 
                              ! End time of calculation
    & DelTime                 ! $ \Delta t $ [s]

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

  ! 惑星表面特性の設定
  ! 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, &
                              ! 成分の数
                              ! Number of composition
    &                    a_QMixName, &
                              ! 成分の変数名
                              ! Name of variables for composition
    &                    a_QMixLongName, &
                              ! 成分の長い変数名
                              ! Long name of variables for composition
    &                    IndexH2OVap, &
    &                    CompositionInqIndex


  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax     ! 鉛直層数. 
                              ! Number of vertical level

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only:  &
    & LatentHeat      , &
    & LatentHeatFusion, &
    & Grav
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

  !
  !
  use dc_message, only: MessageNotify

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

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

  !
  ! Routines for GABLS tests
  !
  use gabls, only : SetGabls2SurfTemp

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

  ! 宣言文 ; Declaration statements
  !
  implicit none

  character(*), parameter:: prog_name = 'dcpam_main'
                            ! 主プログラム名. 
                            ! Main program name

  ! 予報変数 (ステップ $ 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_SurfSnowB (:,:)
                              ! $ M_ss (t-\Delta t) $ . 積雪量 (kg m-2)
                              ! Surface snow amount (kg m-2)

  real(DP), allocatable:: xy_SurfMajCompIceB (:,:)
                              ! $ M_mcs (t-\Delta t) $ . (kg m-2)
                              ! Surface major component ice amount (kg m-2)
  real(DP), allocatable:: xy_SurfMajCompIceN (:,:)
                              ! $ M_mcs (t) $ . (kg m-2)
                              ! Surface major component ice amount (kg m-2)
  real(DP), allocatable:: xy_SurfMajCompIceA (:,:)
                              ! $ M_mcs (t+\Delta t) $ . (kg m-2)
                              ! Surface major component ice 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:: xyz_DTurKinEneDt(:,:,:)
                              ! 
                              ! Turbulent kinetic energy tendency (m2 s-3)

  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_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)
  integer , allocatable:: xy_SurfType (:,:)
                              ! 惑星表面タイプ (土地利用, Matthews 分布) (1)
                              ! Surface type (land use type classified by Matthews) (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)

  integer , allocatable:: xy_PhyImplSDHIndexCalcMethod(:,:)
                              !
                              ! Index for calculation method used in PhyImplSDHTendency

  real(DP), allocatable:: xyr_Temp (:,:,:)
                              ! $ \hat{T} $ . 温度 (半整数レベル) (K)
                              ! Temperature (half level) (K)
  real(DP), allocatable:: xyz_VirTemp (:,:,:)
                              ! $ T_v $ . 仮温度 (K)
                              ! Virtual temperature (K)
  real(DP), allocatable:: xyr_VirTemp (:,:,:)
                              ! $ \hat{T}_v $ . 仮温度 (半整数レベル) (K)
                              ! Virtual temperature (half level) (K)
  real(DP), allocatable:: xy_SurfVirTemp (:,:)
                              ! $ \hat{T}_{v,s} $ . 仮温度 (惑星表面) (K)
                              ! Virtual temperature (surface) (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_RadLUwFlux (:,:,:)
                              ! 長波フラックス (W m-2)
                              ! Upward longwave flux (W m-2)
  real(DP), allocatable:: xyr_RadLDwFlux (:,:,:)
                              ! 長波フラックス (W m-2)
                              ! Downward longwave flux (W m-2)
  real(DP), allocatable:: xyr_RadLFlux  (:,:,:)
  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 (:,:,:)
  real(DP), allocatable:: xyr_RadSUwFlux (:,:,:)
                              ! 短波 (日射) フラックス (W m-2)
                              ! Upward shortwave flux (W m-2)
  real(DP), allocatable:: xyr_RadSDwFlux (:,:,:)
                              ! 短波 (日射) フラックス (W m-2)
                              ! Downward shortwave flux (W m-2)
  real(DP), allocatable:: xyra_DelRadLFlux   (:,:,:,:)
                              ! 長波地表温度変化 (W m-2)
  real(DP), allocatable:: xyra_DelRadLUwFlux (:,:,:,:)
                              ! 長波地表温度変化 (W m-2)
                              ! 
  real(DP), allocatable:: xyra_DelRadLDwFlux (:,:,:,:)
                              ! 長波地表温度変化 (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:: xy_SurfMomFluxX (:,:)
                              ! 惑星表面東西方向運動量フラックス
                              ! Eastward momentum flux at surface
  real(DP), allocatable:: xy_SurfMomFluxY (:,:)
                              ! 惑星表面南北方向運動量フラックス. 
                              ! Northward momentum flux at surface
  real(DP), allocatable:: xy_SurfHeatFlux (:,:)
                              ! 惑星表面熱フラックス. 
                              ! Heat flux at surface
  real(DP), allocatable:: xyf_SurfQMixFlux(:,:,:)
                              ! 惑星表面成分質量フラックス. 
                              ! Mass flux of compositions at surface

  real(DP), allocatable:: xy_SurfH2OVapFluxA(:,:)
                              ! 惑星表面水蒸気フラックス.
                              ! Water vapor flux at the surface
  real(DP), allocatable:: xy_SurfLatentHeatFluxA(:,:)
                              ! 惑星表面潜熱フラックス.
                              ! Latent heat flux at the surface
        ! NOTE:
        ! Only if the evaporation of liquid water is considered, a variable, 
        ! xy_SurfLatentHeatFlux is not required, since a latent heat flux
        ! at the surface is equal to water mass flux times latent heat. 
        ! But, if the evaporation of snow is considered, that is not the case 
        ! and a variable for the latent heat flux is required in addition to 
        ! that for the water mass flux.
        !

  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_DPsDt (:,:)
                              !  (Pa s-1)
                              ! Surface pressure tendency (Pa s-1)

  real(DP), allocatable:: xy_DSurfMajCompIceDt (:,:)
                              !  (kg m-2 s-1)
                              ! Surface major component ice tendency (kg m-2 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_DTempDtCond (:,:,:)
                              ! 凝結加熱率. 
                              ! Condensation heating
  real(DP), allocatable:: xyz_DQVapDtCond (:,:,:)
                              ! 凝結比湿変化. 
                              ! Condensation specific humidity tendency

  real(DP), allocatable:: xyz_QH2OLiqforRad(:,:,:)
                              ! Array for liquid water for radiation calculation
                              ! (kg kg-1)
  real(DP), allocatable:: xyz_QH2OSolforRad(:,:,:)
                              ! Array for solid water (ice) for radiation calculation
                              ! (kg kg-1)


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

  integer           :: IDPhysMode                 ! 使用する物理過程
                                                 ! Physics used for an experiment
  !
  integer, parameter:: IDPhysModeFullPhysics     = 1

  integer           :: IDPhyTendMethod           ! 物理過程による変化率の計算方法
                                                 ! Method calculating physics tendency
  !
  integer, parameter:: IDPhyTendMethodImp1LayModel = 10
  integer, parameter:: IDPhyTendMethodImpSoilModel = 11

  integer           :: IDVDiffMethod            ! 
                                                ! Method for vertical diffusion
  !
  integer, parameter:: IDVDiffMethodMY2   = 80
  integer, parameter:: IDVDiffMethodMY25  = 81


  integer           :: IDSfcFluxMethod          ! 
                                                ! Method for surface flux evaluation
  !
  integer, parameter:: IDSfcFluxMethodL82     = 90
  integer, parameter:: IDSfcFluxMethodBH91B94 = 91


  logical :: FlagDCA


  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  & ! (inout) optional
!!$      & )
    xy_SurfHeight = 0.0_DP


    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      call SetGabls2SurfTemp(     &
        & xy_SurfTemp                   & ! (out)
        & )

      ! 地表面条件の設定
      ! Configure surface conditions
      !
      call SetSurfaceProperties(                           &
        & xy_SurfMajCompIceB, & ! (in)    optional
!!$        & xy_SurfTemp            = xy_SurfTemp,            & ! (inout) optional
        & xy_SurfHumidCoef       = xy_SurfHumidCoef,       & ! (inout) optional
        & xy_SurfRoughLength     = xy_SurfRoughLength      & ! (inout) optional
!!$        & xy_SurfHeatCapacity    = xy_SurfHeatCapacity,    & ! (inout) optional
!!$        & xy_DeepSubSurfHeatFlux = xy_DeepSubSurfHeatFlux, & ! (inout) optional
!!$        & xy_SurfCond            = xy_SurfCond,            & ! (inout) optional
!!$        & xy_SurfType            = xy_SurfType,            & ! (inout) optional
!!$        & xy_SeaIceConc          = xy_SeaIceConc,          & ! (inout) optional
!!$        & xy_SoilHeatCap         = xy_SoilHeatCap,         & ! (inout) optional
!!$        & xy_SoilHeatDiffCoef    = xy_SoilHeatDiffCoef     & ! (inout) optional
        & )

      xy_SurfHeatCapacity    = 0.0_DP
      xy_DeepSubSurfHeatFlux = 0.0_DP
      xy_SurfCond            = 0
      xy_SurfType            = 0
      xy_SeaIceConc          = 0.0_DP
      xy_SoilHeatCap         = 0.0_DP
      xy_SoilHeatDiffCoef    = 0.0_DP

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


      ! Net flux is calculated.
      !
      xyr_RadSFlux     = 0.0_DP
      xyr_RadLFlux     = 0.0_DP
      xyra_DelRadLFlux = 0.0_DP


      ! 地表面フラックス
      ! Surface flux
      !
      select case ( IDSfcFluxMethod )
      case ( IDSfcFluxMethodL82 )
        call SurfaceFlux(                                                 &
          & 'L82',                                                        & ! (in)
          & xyz_UB, xyz_VB,                                               & ! (in)
          & xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, & ! (in)
          & xyzf_QMixB,                                                    & ! (in)
          & xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner,   & ! (in)
          & xy_SurfTemp, xy_SurfHumidCoef,                                & ! (in)
          & xy_SurfRoughLength, xy_SurfRoughLength / 10.0_DP,             & ! (in)
          & xy_SurfMomFluxX, xy_SurfMomFluxY, xy_SurfHeatFlux, xyf_SurfQMixFlux,      & ! (out)
          & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                    & ! (out)
          & xy_SurfQVapTransCoef                                          & ! (out)
          & )
      case ( IDSfcFluxMethodBH91B94 )
        call SurfaceFlux(                                                 &
          & 'BH91B94',                                                    & ! (in)
          & xyz_UB, xyz_VB,                                               & ! (in)
          & xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, & ! (in)
          & xyzf_QMixB,                                                    & ! (in)
          & xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner,   & ! (in)
          & xy_SurfTemp, xy_SurfHumidCoef,                                & ! (in)
          & xy_SurfRoughLength, xy_SurfRoughLength / 10.0_DP,             & ! (in)
          & xy_SurfMomFluxX, xy_SurfMomFluxY, xy_SurfHeatFlux, xyf_SurfQMixFlux,      & ! (out)
          & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                    & ! (out)
          & xy_SurfQVapTransCoef                                          & ! (out)
          & )
      end select


      ! 鉛直拡散フラックス
      ! Vertical diffusion flux
      !
      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY2 )

        call VDiffusion(                                              &
          & xyz_UB,     xyz_VB,     xyzf_QMixB,                       & ! (in)
          & xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in)
          & xy_SurfHeight,                                            & ! (in)
          & xyz_Height, xyr_Height, xyz_Exner,    xyr_Exner,          & ! (in)
          & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux,  & ! (out)
          & xyr_VelTransCoef, xyr_TempTransCoef,                      & ! (out)
          & xyr_QMixTransCoef                                         & ! (out)
          & )

      case ( IDVDiffMethodMY25 )

        if ( CompositionInqIndex( 'TKE' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'TKE is not found.' )
        end if

        call VDiffusionMY25(                                          &
          & xyz_UB,     xyz_VB,     xyzf_QMixB,                       & ! (in)
          & xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, & ! (in)
          & xy_SurfHeight,                                            & ! (in)
          & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,             & ! (in)
          & xyzf_QMixB(:,:,:,CompositionInqIndex( 'TKE' )),           & ! (in )
          & xy_SurfMomFluxX, xy_SurfMomFluxY,                         & ! (in)
          & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux,  & ! (out)
          & xyr_VelTransCoef, xyr_TempTransCoef,                      & ! (out)
          & xyr_QMixTransCoef,                                        & ! (out)
          & xyz_DTurKinEneDt                                          & ! (out)
          & )

      end select

      xyr_MomFluxX (:,:,0)   = xy_SurfMomFluxX
      xyr_MomFluxY (:,:,0)   = xy_SurfMomFluxY
      xyr_HeatFlux (:,:,0)   = xy_SurfHeatFlux
      xyrf_QMixFlux(:,:,0,:) = xyf_SurfQMixFlux

      xy_SurfLatentHeatFluxA = LatentHeat * xyf_SurfQMixFlux(:,:,IndexH2OVap)



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

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

        xy_SurfH2OVapFluxA     = 0.0_DP
        xy_SurfLatentHeatFluxA = 0.0_DP

        xyz_DSoilTempDt      = 0.0_DP
        xy_DPsDt             = 0.0_DP
        xy_DSurfMajCompIceDt = 0.0_DP
        xy_DSoilMoistDt      = 0.0_DP
        xy_DSurfSnowDt       = 0.0_DP

      case ( IDPhyTendMethodImpSoilModel )


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


        call PhyImplSDHSetMethodFromMatthews(  &
          & xy_SurfType, xy_SeaIceConc,        & ! (in)
          & xy_PhyImplSDHIndexCalcMethod       & ! (out)
          & )
        call PhyImplSDHTendency(                                      &
          & xy_PhyImplSDHIndexCalcMethod,                             & ! (in)
          & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux,  & ! (in)
          & xy_SurfLatentHeatFluxA,                                   & ! (in)
          & xyr_SoilHeatFlux,                                         & ! (in)
          & xyr_RadSFlux, xyr_RadLFlux,                               & ! (in)
          & xy_DeepSubSurfHeatFlux,                                   & ! (in)
          & xy_SurfTemp, xyz_SoilTemp,                                & ! (in)
          & xy_SurfHumidCoef,                                         & ! (in)
          & xy_SurfHeatCapacity,                                      & ! (in)
          & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                      & ! (in)
          & xyra_DelRadLFlux,                                         & ! (in)
          & xyr_Press, xyz_Exner, xyr_Exner,                          & ! (in)
          & xyr_VelTransCoef, xyr_TempTransCoef,                      & ! (in)
          & xyr_QMixTransCoef,                                        & ! (in)
          & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                & ! (in)
          & xy_SurfQVapTransCoef,                                     & ! (in)
          & xyr_SoilTempTransCoef,                                    & ! (in)
          & xy_SurfMajCompIceB,                                       & ! (in)
          & xy_SurfSnowB,                                             & ! (in)
          & xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt,       & ! (out)
          & xy_DSurfTempDt,                                           & ! (out)
          & xyz_DSoilTempDt,                                          & ! (out)
          & xy_DPsDt, xy_DSurfMajCompIceDt,                           & ! (out)
          & xy_DSoilMoistDt,                                          & ! (out)
          & xy_DSurfSnowDt                                            & ! (out)
          & )


!!$        ! This is temporal treatment.
!!$        xy_SurfH2OVapFluxA = xyrf_QMixFlux(:,:,0,IndexH2OVap)

!!$        call PhyImplSDHV2SetMethodMatthews(  &
!!$          & xy_SurfType, xy_SeaIceConc,      & ! (in)
!!$          & xy_PhyImplSDHIndexCalcMethod     & ! (out)
!!$          & )
!!$        call PhyImplSDHV2Tendency(                                    &
!!$          & xy_SurfType,                                              & ! (in)
!!$          & xy_PhyImplSDHIndexCalcMethod,                             & ! (in)
!!$          & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux,  & ! (in)
!!$          & xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA,               & ! (out)
!!$          & xyr_SoilHeatFlux,                                         & ! (in)
!!$          & xyr_RadSFlux, xyr_RadLFlux,                               & ! (in)
!!$          & xy_DeepSubSurfHeatFlux,                                   & ! (in)
!!$          & xyz_TempB, xy_SurfTemp, xyz_SoilTemp,                     & ! (in)
!!$          & xy_SurfHumidCoef,                                         & ! (in)
!!$          & xy_SurfHeatCapacity,                                      & ! (in)
!!$          & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                      & ! (in)
!!$          & xyra_DelRadLFlux,                                         & ! (in)
!!$          & xyr_Press, xyz_Exner, xyr_Exner,                          & ! (in)
!!$          & xyr_VelTransCoef, xyr_TempTransCoef,                      & ! (in)
!!$          & xyr_QMixTransCoef,                                        & ! (in)
!!$          & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                & ! (in)
!!$          & xy_SurfQVapTransCoef,                                     & ! (in)
!!$          & xyr_SoilTempTransCoef,                                    & ! (in)
!!$          & xy_SurfMajCompIceB,                                       & ! (in)
!!$          & xy_SoilMoistB, xy_SurfSnowB,                              & ! (in)
!!$          & xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt,       & ! (out)
!!$          & xy_DSurfTempDt,                                           & ! (out)
!!$          & xyz_DSoilTempDt,                                          & ! (out)
!!$          & xy_DPsDt, xy_DSurfMajCompIceDt,                           & ! (out)
!!$          & xy_DSoilMoistDt,                                          & ! (out)
!!$          & xy_DSurfSnowDt                                            & ! (out)
!!$          & )

      end select


      ! Overwrite tendency of turbulent kinetic energy
      !
      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY25 )
        if ( CompositionInqIndex( 'TKE' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'TKE is not found.' )
        end if
        xyzf_DQMixDt(:,:,:,CompositionInqIndex( 'TKE' )) = &
          & xyz_DTurKinEneDt
      end select


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



      ! 鉛直拡散フラックスの出力 
      !   * 出力のみのサブルーチンであり, 計算には影響しない
      ! 
      ! 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, & ! (in)
        & xyz_DUDt,  xyz_DVDt,  xyz_DTempDtVDiff,  xyzf_DQMixDt,   & ! (in)
        & xyz_Exner, xyr_Exner,                                    & ! (in)
        & xyr_VelTransCoef, xyr_TempTransCoef,                     & ! (in)
        & xyr_QMixTransCoef                                        & ! (in)
        & )

      ! 地表面フラックスの出力 
      !   * 出力のみのサブルーチンであり, 計算には影響しない
      ! 
      ! 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,  & ! (in)
        & xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA,               & ! (in)
        & xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt,       & ! (in)
        & xy_SurfTemp, xy_DSurfTempDt,                              & ! (in)
        & xyr_Press, xyz_Exner, xyr_Exner, xy_SurfHumidCoef,        & ! (in)
        & xy_SurfVelTransCoef, xy_SurfTempTransCoef,                & ! (in)
        & xy_SurfQVapTransCoef                                      & ! (in)
        & )

    end select


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

    ! Surface pressure is temporarily adjusted, here.
    !
    xy_PsA = xy_PsA + xy_DPsDt * ( 2.0_DP * DelTime )


    if ( FlagDCA ) then

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

      ! 乾燥対流調節
      ! Dry convective adjustment
      !
      call DryConvAdjust(                        &
        & xyz_TempA, xyz_UA, xyz_VA, xyzf_QMixA, &  ! (inout)
        & xyz_Press, xyr_Press                   &  ! (in)
        & )
    end if


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


    ! Messages for debug run
    !
!!$    select case ( IDPhysMode )
!!$    case ( IDPhysModeFullPhysics )
!!$      write( 6, * ) xyz_TempA(0,jmax/2+1,1), &
!!$        & xyz_VirTemp(0,jmax/2+1,1) - xyz_TempB(0,jmax/2+1,1), &
!!$        & xyr_VirTemp(0,jmax/2+1,1) - xyr_Temp (0,jmax/2+1,1)
!!$    case default
!!$      write( 6, * ) xyz_TempA(0,jmax/2+1,1)
!!$    end select


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

    ! ヒストリデータ出力
    ! 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 )

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    call AuxVars( &
      & xy_PsN, xyz_TempN, xyzf_QMixN(:,:,:,IndexH2OVap),    & ! (in )
      & xy_SurfHeight = xy_SurfHeight,                       & ! (in ) optional
      & xyz_Height    = xyz_Height                           & ! (out) optional
      & )
    call HistoryAutoPut( TimeN, 'Height', xyz_Height )

    call AuxVars( &
      & xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap),    & ! (in )
      & xyz_Exner     = xyz_Exner                            & ! (out) optional
      & )
    call HistoryAutoPut( TimeN, 'PotTemp', xyz_TempA / xyz_Exner )

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

      call HistoryAutoPut( TimeN, 'DTempDtCond', xyz_DTempDtCond )
      call HistoryAutoPut( TimeN, 'DQVapDtCond', xyz_DQVapDtCond )

      call HistoryAutoPut( TimeN, 'SeaIceConc'     , xy_SeaIceConc      )
      call HistoryAutoPut( TimeN, 'SurfRoughLength', xy_SurfRoughLength )
    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.


    ! 時刻の進行
    ! 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, xy_PsB, &  ! (in)
      & xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN  &  ! (in)
      & )

    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

    ! 時刻管理
    ! 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, & ! 経度格子点数. 
                                 ! Number of grid points in longitude
      &                  jmax, & ! 緯度格子点数. 
                                 ! Number of grid points in latitude
      &                  kmax, & ! 鉛直層数. 
                                 ! Number of vertical level
      &                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: RestartFileInit, RestartFileOpen, RestartFileGet

    ! ヒストリデータ出力
    ! 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

    ! 時系列データの読み込み
    ! Reading time series
    !
    use read_time_series, only : ReadTimeSeriesInit

    ! 地表面データ提供
    ! Prepare surface data
    !
    use surface_data, only : SurfDataInit

    ! ファイルから 1 次元プロファイルを読んで設定する. 
    ! read 1-D profile from a file and set it 
    !
    use set_1d_profile, only : Set1DProfileInit


    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    use auxiliary, only : AuxVarsInit, AuxVars

    ! 惑星表面データの設定
    ! Setting planetary surface properties
    !
    use surface_properties, only : SurfacePropertiesInit

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

    ! 陰解法による時間積分
    ! Time integration with implicit scheme
    !
    use phy_implicit, only : PhyImplInit

    ! 地下熱伝導モデルを用いた場合の陰解法による時間積分
    !
    ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
    use phy_implicit_sdh, only : PhyImplSDHInit

    ! 地下熱伝導モデルを用いた場合の陰解法による時間積分
    !
    ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
    use phy_implicit_sdh_V2, only : PhyImplSDHV2Init

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

    ! 鉛直拡散フラックス (Mellor and Yamada, 1974)
    ! Vertical diffusion flux (Mellor and Yamada, 1974)
    !
    use vdiffusion_my, only : VDiffusionInit

    ! 地表面フラックス (バルク法)
    ! Surface flux (Bulk method)
    !
    use surface_flux_bulk, only : SurfaceFluxInit

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

    ! 時間フィルター (Williams, 2009)
    ! Time filter (Williams, 2009)
    !
    use timefilter_williams2009, only: TimeFilterWilliams2009Init

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

    !
    ! Routines for GABLS tests
    !
    use gabls, only : GablsInit

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    character(*), parameter:: version = &
      & '$Name:  $' // &
      & '$Id: dcpam_main_gabls2.f90,v 1.1 2013/09/21 14:58:39 yot Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version

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

    character(STRING)      :: DynMode
                                ! Dynamics used in calculation
    character(STRING)      :: PhysMode
                                ! Physics used in calculation
    character(STRING)      :: SfcFluxMethod
                                ! Method for surface flux evaluation used in calculation
    character(STRING)      :: VDiffMethod
                                ! Method for vertical diffusion evaluation used in calculation
    character(STRING)      :: PhysImpMode
                                ! Mode for implicit method used in calculation

    logical                :: FlagPhysImpSoilModelSO
                                ! flag for use of slab ocean
    logical                :: FlagSnow
                                ! flag for treating snow
    logical                :: FlagMajCompPhaseChange
                                ! flag for use of major component phase change

    character(STRING):: CondMajCompName
                                ! name of condensable major component

    character(STRING):: briefexpldyn
                              ! 実行ファイルの簡潔な説明 (力学過程)
                              ! Brief account of executable file (dynamics)
    character(STRING):: briefexplphy
                              ! 実行ファイルの簡潔な説明 (物理過程)
                              ! Brief account of executable file (physics)
    character(STRING):: briefexplsfcflux
                              ! 実行ファイルの簡潔な説明 (惑星表面フラックス)
                              ! Brief account of executable file (surface flux)
    character(STRING):: briefexplvdiff
                              ! 実行ファイルの簡潔な説明 (鉛直拡散過程)
                              ! Brief account of executable file (vertical diffusion)
    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:: n               ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /dcpam_main_nml/                                       &
      & DynMode, PhysMode, SfcFluxMethod, VDiffMethod, PhysImpMode,  &
      & FlagDCA,                                                    &
      & FlagSnow, FlagMajCompPhaseChange, CondMajCompName
          !
          ! デフォルト値については初期化手続 "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(       &
      & namelist_filename,    & ! (out)
      & prog_name            & ! (in )
      & )

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

    ! デフォルト値の設定
    ! Default values settings
    !
    DynMode                 = 'NoHorAdv'

    PhysMode                = 'FullPhysics'
!!$    PhysMode                = 'NoPhysics'

    SfcFluxMethod           = 'L82'
!!$    SfcFluxMethod           = 'BH91B94'

    VDiffMethod             = 'MY2'
!!$    VDiffMethod             = 'MY2.5'

    PhysImpMode             = '1LayModel'
!!$    PhysImpMode             = 'SoilModel'
!!$    PhysImpMode             = 'SoilModelSO'

    FlagDCA = .false.

    FlagSnow                = .false.

    FlagMajCompPhaseChange  = .false.

    CondMajCompName         = ''
!!$    CondMajCompName         = 'CO2'


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

      rewind( unit_nml )
      read( unit_nml, &            ! (in)
        & nml = dcpam_main_nml, &  ! (out)
        & 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, &
      & 'DynMode=<%c>.', &
      & c1 = trim(DynMode) )
    !
    select case ( DynMode )
    case ( 'NoHorAdv' )
      IDDynMode = IDDynModeNoHorAdv
    case default
      call MessageNotify( 'E', prog_name, &
        & 'DynMode=<%c> is not supported.', &
        & c1 = trim(DynMode) )
    end select


    ! Identification of calculation method for physics
    !
    call MessageNotify( 'M', prog_name, &
      & 'PhysMode=<%c>.', &
      & c1 = trim(PhysMode) )
    !
    select case ( PhysMode )
    case ( 'FullPhysics' )
      IDPhysMode = IDPhysModeFullPhysics
    case default
      call MessageNotify( 'E', prog_name, &
        & 'PhysMode=<%c> is not supported.', &
        & c1 = trim(PhysMode) )
    end select


    ! Identification of calculation method for surface flux
    !
    call MessageNotify( 'M', prog_name, &
      & 'SfcFluxMethod=<%c>.', &
      & c1 = trim(SfcFluxMethod) )
    !
    select case ( SfcFluxMethod )
    case ( 'L82' )
      IDSfcFluxMethod = IDSfcFluxMethodL82
    case ( 'BH91B94' )
      IDSfcFluxMethod = IDSfcFluxMethodBH91B94
    case default
      call MessageNotify( 'E', prog_name, &
        & 'SfcFluxMethod=<%c> is not supported.', &
        & c1 = trim(SfcFluxMethod) )
    end select


    ! Identification of calculation method for vertical diffusion
    !
    call MessageNotify( 'M', prog_name, &
      & 'VDiffMethod=<%c>.', &
      & c1 = trim(VDiffMethod) )
    !
    select case ( VDiffMethod )
    case ( 'MY2' )
      IDVDiffMethod = IDVDiffMethodMY2
    case ( 'MY2.5' )
      IDVDiffMethod = IDVDiffMethodMY25
    case default
      call MessageNotify( 'E', prog_name, &
        & 'VDiffMethod=<%c> is not supported.', &
        & c1 = trim(VDiffMethod) )
    end select


    ! Identification of calculation method for solving simultaneous linear equations 
    ! of physics
    !
    call MessageNotify( 'M', prog_name, &
      & 'PhysImpMode=<%c>.', &
      & c1 = trim(PhysImpMode) )
    !
    select case ( PhysImpMode )
    case ( '1LayModel' )
      IDPhyTendMethod = IDPhyTendMethodImp1LayModel
      FlagPhysImpSoilModelSO = .false.
    case ( 'SoilModel' )
      IDPhyTendMethod = IDPhyTendMethodImpSoilModel
      FlagPhysImpSoilModelSO = .false.
    case ( 'SoilModelSO' )
      IDPhyTendMethod = IDPhyTendMethodImpSoilModel
      FlagPhysImpSoilModelSO = .true.
    case default
      call MessageNotify( 'E', prog_name, &
        & 'PhysImpMode=<%c> is not supported.', &
        & c1 = trim(PhysImpMode) )
    end select
    !
    !   Check for value of FlagFullPhysics
    !
    if ( ( IDPhysMode /= IDPhysModeFullPhysics ) .and. &
      &  ( IDPhyTendMethod == IDPhyTendMethodImpSoilModel ) ) then
      call MessageNotify( 'E', prog_name, &
        & 'PhysMode has to be "FullPhysics" true, when PhyImpMode is "SoilModel" or "SoilModelSO".' )
    end if


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

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      briefexplphy = 'parameterization suite is used'
    end select

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      select case ( IDSfcFluxMethod )
      case ( IDSfcFluxMethodL82 )
        briefexplsfcflux = 'surface flux is estimated with method of Louis et al. (1982)'
      case ( IDSfcFluxMethodBH91B94 )
        briefexplsfcflux = 'surface flux is estimated with method of Beljaars and Holtslag (1991), Beljaars (1994)'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplsfcflux', '' )
      end select

      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY2 )
        briefexplvdiff = 'system with Mellor and Yamada level 2'
      case ( IDVDiffMethodMY25 )
        briefexplvdiff = 'system with Mellor and Yamada level 2.5'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplvdiff', '' )
      end select

      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )
        briefexplimp = 'system with surface 1 layer model'
      case ( IDPhyTendMethodImpSoilModel )
        briefexplimp = 'system with thermal diffusion soil model'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplimp', '' )
      end select
    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 ( IDPhysMode == IDPhysModeFullPhysics ) then
      call MessageNotify( 'M', prog_name, '|    Surface flux model : %c', c1 = trim(briefexplsfcflux) )
      call MessageNotify( 'M', prog_name, '|    Vertical diffusion model : %c', c1 = trim(briefexplvdiff) )
      call MessageNotify( 'M', prog_name, '|    Implicit method : %c', c1 = trim(briefexplimp) )
      call MessageNotify( 'M', prog_name, '|    Major component phase change : %b', l = (/ FlagMajCompPhaseChange /) )
    end if
    call MessageNotify( 'M', prog_name,   '|  -- version = %c', c1 = trim(version) )
    call MessageNotify( 'M', prog_name,   '+-------------------------------------' )
    call MessageNotify( 'M', prog_name, '' )


    ! Initialization of modules used in this module
    !


    ! 時刻管理
    ! 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



    ! 時系列データの読み込み
    ! Reading time series
    !
    call ReadTimeSeriesInit

    ! 地表面データ提供
    ! Prepare surface data
    !
    call SurfDataInit

    ! ファイルから 1 次元プロファイルを読んで設定する. 
    ! read 1-D profile from a file and set it 
    !
    call Set1DProfileInit





    ! 予報変数の割付
    ! 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 RestartFileInit
    call RestartFileGet( &
      & xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, & ! (out)
      & xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, & ! (out)
      & flag_initial )                                   ! (out) optional


    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      ! 地表面温度, 土壌温度の割付
      ! 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_SurfSnowB      (0:imax-1, 1:jmax) )

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


      xy_SurfSnowB       = 0.0_DP
      xy_SurfMajCompIceB = 0.0_DP


!!$      xy_SurfSnowN = 1.0d1
!!$      xy_SurfSnowB = xy_SurfSnowN

      ! 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' , &         ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'eastward wind', 'm s-1' )               ! (in)

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

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

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

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

    call HistoryAutoAddVariable( 'Height' , &    ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'geopotential height', 'm' )             ! (in)

    call HistoryAutoAddVariable( 'PotTemp' ,  &  ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'potential temperature', 'K' )           ! (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 ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      ! ヒストリデータ出力のためのへの変数登録
      ! 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( '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( xy_DPsDt    (0:imax-1, 1:jmax) )
    allocate( xyzf_DQMixDt(0:imax-1, 1:jmax, 1:kmax, 1:ncmax) )

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

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

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      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_SurfType           (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( xy_PhyImplSDHIndexCalcMethod(0:imax-1, 1:jmax) )

      allocate( xyr_Temp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyz_VirTemp   (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_VirTemp   (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xy_SurfVirTemp(0:imax-1, 1:jmax) )
      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0: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_RadLUwFlux    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadLDwFlux    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSFlux      (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSUwFlux    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_RadSDwFlux    (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyra_DelRadLFlux  (0:imax-1, 1:jmax, 0:kmax,  0:1) )
      allocate( xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax,  0:1) )
      allocate( xyra_DelRadLDwFlux(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( xy_SurfMomFluxX  (0:imax-1, 1:jmax) )
      allocate( xy_SurfMomFluxY  (0:imax-1, 1:jmax) )
      allocate( xy_SurfHeatFlux  (0:imax-1, 1:jmax) )
      allocate( xyf_SurfQMixFlux (0:imax-1, 1:jmax, 1:ncmax) )

      allocate( xy_SurfH2OVapFluxA    (0:imax-1, 1:jmax) )
      allocate( xy_SurfLatentHeatFluxA(0:imax-1, 1:jmax) )

      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_DSurfMajCompIceDt(0:imax-1, 1:jmax) )
      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_DTempDtCond (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DQVapDtCond (0:imax-1, 1:jmax, 1:kmax) )

      allocate( xyz_QH2OLiqforRad   (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_QH2OSolforRad   (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( 'SurfRoughLength' , &
        & (/ 'lon ', 'lat ', 'time' /), &
        & 'surface roughness length', 'm' )

    end select




    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      call GablsInit

      ! 惑星表面データの設定
      ! Setting planetary surface properties
      !
      call SurfacePropertiesInit(                    &
        & FlagPhysImpSoilModelSO, .false., FlagSnow  & ! (in)
        & )


      ! 補助的な変数を計算するサブルーチン・関数群
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! 鉛直拡散フラックス (Mellor and Yamada, 1974)
      ! Vertical diffusion flux (Mellor and Yamada, 1974)
      !
      call VDiffusionInit

      ! 地表面フラックス (バルク法)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

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


      case ( IDPhyTendMethodImpSoilModel )

        ! 地下における熱の鉛直拡散
        ! Vertical diffusion of heat under the ground
        !
        call SubsurfaceDiffusionInit

        ! 地下熱伝導モデルを用いた場合の陰解法による時間積分
        !
        ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
        call PhyImplSDHInit(                                                 &
          & .false., FlagSnow,                                               &
          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
          & )

!!$        call PhyImplSDHV2Init(                                                &
!!$          & .false., FlagSnow,                                       &
!!$          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
!!$          & )

      end select


      ! 陰解法による時間積分のためのルーチン
      ! Routines for time integration with implicit scheme
      !
      call PhyImplUtilsInit

      ! 鉛直拡散フラックス (Mellor and Yamada, 1974)
      ! Vertical diffusion flux (Mellor and Yamada, 1974)
      !
      call VDiffusionInit

      ! 地表面フラックス (バルク法)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

    end select


    ! 力学過程
    ! Dynamical core
    !
    select case ( IDDynMode )
    case ( IDDynModeNoHorAdv )
      ! 物理過程のみの計算のための力学過程
      ! A dynamics for calculation with physical processes only
      !
      call DynamicsPhysicsOnlyInit
    end select


    !
    ! Dry convective adjustment
    !
    if ( FlagDCA ) then
      call DryConvAdjustInit
    end if


    ! 時間フィルター (Williams, 2009)
    ! Time filter (Williams, 2009)
    !
    call TimeFilterWilliams2009Init

    ! 予報変数の値の確認
    ! Check values of prognostic variables
    !
    call CheckProgVarsInit

    ! 補助的な変数を計算するサブルーチン・関数群
    ! Subroutines and functions for calculating auxiliary variables
    !
    call AuxVarsInit

    !
    ! End of initialization
    !


    ! 初回だけはオイラー法を用いるため, Δ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

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

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

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

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

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

    ! ヒストリデータファイルクローズ
    ! 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( xyz_DTurKinEneDt )

    deallocate( xy_SurfHeight )
    deallocate( xyz_Height )

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      deallocate( xy_SurfHumidCoef       )
      deallocate( xy_SurfRoughLength     )
      deallocate( xy_SurfHeatCapacity    )
      deallocate( xy_SeaIceConc          )
      deallocate( xy_SurfCond            )
      deallocate( xy_SurfType            )
      deallocate( xy_DeepSubSurfHeatFlux )
      deallocate( xy_SoilHeatCap         )
      deallocate( xy_SoilHeatDiffCoef    )

      deallocate( xy_PhyImplSDHIndexCalcMethod )

      deallocate( xyr_Temp   )
      deallocate( xyz_VirTemp    )
      deallocate( xyr_VirTemp    )
      deallocate( xy_SurfVirTemp )
      deallocate( xyz_Press  )
      deallocate( xyr_Press  )
      deallocate( xyr_Height )
      deallocate( xyz_Exner  )
      deallocate( xyr_Exner  )

      deallocate( xyr_RadLFlux     )
      deallocate( xyr_RadLFluxA    )
      deallocate( xyr_RadLUwFlux   )
      deallocate( xyr_RadLDwFlux   )
      deallocate( xyr_RadSFlux     )
      deallocate( xyr_RadSUwFlux     )
      deallocate( xyr_RadSDwFlux     )
      deallocate( xyra_DelRadLUwFlux )
      deallocate( xyra_DelRadLDwFlux )

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

      deallocate( xy_SurfMomFluxX  )
      deallocate( xy_SurfMomFluxY  )
      deallocate( xy_SurfHeatFlux  )
      deallocate( xyf_SurfQMixFlux )

      deallocate( xy_SurfH2OVapFluxA     )
      deallocate( xy_SurfLatentHeatFluxA )

      deallocate( xyr_SoilHeatFlux )

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

      deallocate( xy_DSurfMajCompIceDt )
      deallocate( xy_DSoilMoistDt      )
      deallocate( xy_DSurfSnowDt       )

      deallocate( xyz_DTempDtVDiff )

      deallocate( xyz_DTempDtCond  )
      deallocate( xyz_DQVapDtCond  )

      deallocate( xyz_QH2OLiqforRad    )
      deallocate( xyz_QH2OSolforRad    )
    end if ! FlagFullPhysics



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

    call AxessetFinalize

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

    ! Finalize MPI
    !
    call MPIWrapperFinalize


  end subroutine MainTerminate

end program dcpam_main_gabls2