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: DynamicsHsplVAS83

  ! ����������¿ã���ç®�����������å­����
  ! A dynamics for calculation with physical processes only
  !
  use dynamics_physicsonly, only: DynamicsPhysicsOnly

  !
  ! Dynamical process for TWP-ICE experiment
  !
  use dynamics_twpice_scm_exp, only : DynamicsTWPICESCMExp

  ! 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 rad_DennouAGCM, only: RadDennouAGCMFlux

  ! �¾å��¢é�£ã���¼ã����
  ! Routines for radiation calculation
  !
  use rad_utils, only : RadDTempDt, RadFluxOutput

  ! �°ç��大æ������¾å��¢ã���� Ver. 2
  ! radiation model for the Earth's atmosphere Ver. 2
  !
  use rad_Earth_V2, only: RadEarthV2Flux

  ! wrapper of RRTMG
  ! wrapper of RRTMG
  !
  use rad_rrtmg_wrapper, only: RadRRTMGWrapperFlux

  ! ����大æ������¾å��¢ã���� Ver. 1
  ! radiation model for the Mars' atmosphere Ver. 1
  !
  use rad_Mars_V1, only: RadMarsV1Flux

!!$  ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$  ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$  !
!!$  use rad_15m_NLTE, only: rad15mNLTEMergeHR

  ! ����è¨�ç®����赤å����±ç�è¨�ç®�
  ! Calculation of near infrared heating rate in the case of Mars
  !
  use rad_Mars_NIR, only : RadMarsNIRINOUT

  ! Schneider and Liu (2009) ���¾å��¢ã����
  ! Radiation model by Schneider and Liu (2009)
  !
  use rad_SL09, only : RadSL09Flux

  ! ç°¡å���¾å��¢ã����
  ! Simple radiation model
  !
  use rad_simple, only : RadSimpleFlux

  ! ä½����������¾å��¢ã����
  ! radiation model with no absorption and no scattering
  !
  use rad_none, only : RadNoneFlux

  ! ���´æ�¡æ�£ã����������
  ! Vertical diffusion flux
  !
  use vdiffusion_my, only: VDiffusionMY25, VDiffusion, VDiffusionOutput
  use vdiffusion_my, only: VDiffusionMY251DWrapper3D
  use vdiffusion_my, only: VDiffusionMY25GBT94

  ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
  ! JMA turbulent mixing module
  !
  use vdiffusion_jma_my_wrapper, only : VDiffusionJMAMYWrapper3D

  ! �°ä����������±ã�����´æ�¡æ��
  ! Vertical diffusion of heat under the ground
  !
  use subsurface_diffusion_heat, only: SubsurfaceDiffusion

  ! Gravity wave drag by McFarlane (1987)
  ! Gravity wave drag by McFarlane (1987)
  !
  use gwd_m1987, only : GWDM1987

  ! ç©��²ã�����¡ã�¿ã���¼ã�¼ã�·ã�§ã�� (対æ�調ç�)
  ! Cumulus parameterization (convection adjust)
  !
  use moist_conv_adjust, only: MoistConvAdjust

  ! Relaxed Arakawa-Schubert scheme
  ! Relaxed Arakawa-Schubert scheme
  !
  use relaxed_arakawa_schubert, only : RAS
  use relaxed_arakawa_schubert, only : RAS1DWrapper3D
  use relaxed_arakawa_schubert, only : RASWithIce1DWrapper3DWrapper

  ! 大è�模å��çµ� (��対æ��§å��çµ�)
  ! Large scale condensation
  !
  use lscond, only: LScaleCond
  use lscond, only: LScaleCond1D3DWrapper
  use lscond, only: LScaleCondLL911D3DWrapper

  ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
  ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
  !
!!$  use lscond_LL91, only : LScaleCondLL91

  !
  != Saturation adjustment
  !
  use saturation_adjust, only : SaturationAdjust

  ! �°è¡¨�¢ã����������
  ! Surface flux
  use surface_flux_bulk, only: SurfaceFlux, SurfaceFluxOutput

  !
  ! set dust flux
  !
!!$  use set_dust_flux, only : SetDustFlux

  ! ����������������
  ! Lower boundary flux
  !
  use lb_flux_simple, only : LBFluxSimple

  ! ä¹¾ç�¥å�æµ�調ç�
  ! Dry convective adjustment
  !
  use dry_conv_adjust, only: DryConvAdjust

  ! ������
  ! Mass fixer
  !
  use mass_fixer, only: MassFixerColumn

  ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
  ! 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:              &
!!$    & PhyImplSDHSetMethodFromMatthews,     &
!!$    & PhyImplSDHTendency,                  &
!!$    & PhyImplSDHCorSOTempBySnowMelt

!!$  ! �°è§£æ³���������������� (ä¸��������������)
!!$  ! Matrices handling for implicit scheme (for a part of physical processes)
!!$  !
!!$  use phy_implicit_sdh_V2, only:       &
!!$    & PhyImplSDHV2SetMethodMatthews,   &
!!$    & PhyImplSDHV2Tendency,            &
!!$    & PhyImplSDHV2CorSOTempBySnowMelt

  ! �°è§£æ³���������������� (ä¸��������������)
  ! Matrices handling for implicit scheme (for a part of physical processes)
  !
  use phy_implicit_sdh_V3, only: PhyImplSDHV3SetMethodMatthews, PhyImplSDHV3Tendency, PhyImplSDHV3CorSOTempBySnowMelt

  ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
  ! 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 : BucketSetFlagOceanFromMatthews, BucketModEvapAndLatentHeatFlux, BucketIntegration, BucketPRCPAdjust

  ! �¿ã�¤ã�����£ã���¿ã�� (Asselin, 1972)
  ! Time filter (Asselin, 1972)
  !
  use timefilter_asselin1972, only: TimeFilter, TimeFilterSurfVars

  ! �������£ã���¿ã�� (Williams, 2009)
  ! Time filter (Williams, 2009)
  !
  use timefilter_williams2009, only: TimeFilterWilliams2009, TimeFilterWilliams2009SurfVars

  ! ���»ç���
  ! 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

  ! ��, æ°·ã���²å��
  ! snow/ice fraction
  !
  use snowice_frac, only : CalcSnowFrac

  ! ���¹ã�������¼ã�¿å�ºå��
  ! 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      , LatentHeatFusion, Grav
                              ! $ g $ [m s-2]. 
                              ! ���������. 
                              ! Gravitational acceleration

  !
  !
  use dc_message, only: MessageNotify

  ! ç¨��¥å�������¡ã��
  ! Kind type parameter
  !
  use dc_types, only: DP, STRING, TOKEN      ! �­ã�¼ã���¼ã��.   Keywords. 

  ! �²ã�����¢ã����
  ! No cloud model
  !
  use cloud_none, only : CloudNone, CloudNoneWithIce

  ! ç°¡å���²ã�¢ã����
  ! Simple cloud
  !
  use cloud_simple, only: CloudSimple, CloudSimpleWithIce, CloudSimpleCalcCloudCover, CloudSimpleDivideWatAndIce, CloudSimpleCalcPRCPKeyLLTemp3D

  ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
  ! Cloud model based on Tiedtke (1993)
  !
  use cloud_T1993base, only : CloudT1993baseWithIce

  ! ���� H2O �²ã�¢ã����
  ! Mars H2O cloud model
  !
  use cloud_mars_h2o, only : CloudMarsH2O

  ! ����������
  ! Gravitational sedimentation process
  !
  use grav_sed, only : GravSed

  ! 主æ�����¸å���
  ! Phase change of atmospheric major component
  !
  use major_comp_phase_change, only : MajorCompPhaseChangeInAtm

  ! äº��±å��°ã���¤ã��¢ºèª�
  ! Check values of prognostic variables
  !
  use check_prog_vars, only: CheckProgVars

  ! Output frequently used variables
  ! Output frequently used variables
  !
  use output_freq_used_vars, only : OutputFreqUsedVars


  ! 宣�� ; 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_SurfMajCompIceB (:,:)
                              ! $ M_mcs (t-\Delta t) $ . (kg m-2)
                              ! Surface major component ice amount (kg m-2)
  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_SurfMajCompIceN (:,:)
                              ! $ M_mcs (t) $ . (kg m-2)
                              ! Surface major component ice 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_SurfMajCompIceA (:,:)
                              ! $ M_mcs (t+\Delta t) $ . (kg m-2)
                              ! Surface major component ice 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:: xyz_DTurKinEneDt(:,:,:)
                              ! 
                              ! Turbulent kinetic energy tendency (m2 s-3)

  real(DP), allocatable:: xy_SurfHeight (:,:)
                              ! $ z_s $ . �°è¡¨�¢é�åº� (m)
                              ! Surface height (m)
  real(DP), allocatable:: xy_SurfHeightStd (:,:)
                              ! 
                              ! Surface height standard deviation(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_SurfRoughLenMom (:,:)
                              ! �°è¡¨ç²�åº��� (m)
                              ! Surface rough length for momentum (m)
  real(DP), allocatable:: xy_SurfRoughLenHeat(:,:)
                              ! �°è¡¨ç²�åº��� (m)
                              ! Surface rough length for heat (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)

  real(DP), allocatable:: xy_SnowFrac(:,:)
                              ! 
                              ! Snow fraction (1)

!!$  logical , allocatable:: xy_FlagMatthewsLand(:,:)
!!$                              !
!!$                              ! Flag for land grid point based on Matthews
  integer , allocatable:: xy_PhyImplSDHIndexCalcMethod(:,:)
                              !
                              ! Index for calculation method used in PhyImplSDHTendency
  logical , allocatable:: xy_BucketFlagOceanGrid(:,:)
                              !
                              ! Flag for ocean grid point used in bucket model

  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_VelDiffCoef (:,:,:)
                              ! �¡æ�£ä��°ï�������. 
                              ! Diffusion coefficient: velocity
  real(DP), allocatable:: xyr_TempDiffCoef (:,:,:)
                              ! �¡æ�£ä��°ï�æ¸�º¦. 
                              ! Transfer coefficient: temperature
  real(DP), allocatable:: xyr_QMixDiffCoef (:,:,:)
                              ! �¡æ�£ä��°ï�æ¯�æ¹�. 
                              ! Diffusion coefficient: specific humidity


  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_SurfMOLength(:,:)
                              ! 
                              ! Monin-Obukov length

  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_DTempDtRadL (:,:,:)
                              ! �·æ³¢���±ç�. 
                              ! Temperature tendency with longwave
  real(DP), allocatable:: xyz_DTempDtRadS (:,:,:)
                              ! �­æ³¢���±ç�. 
                              ! Temperature tendency with shortwave

  real(DP), allocatable:: xyz_DUDtGWD (:,:,:)
                              ! $ \DP{u}{t} $ . �±è¥¿é¢���å¤��� (m s-2)
                              ! Eastward wind tendency by gravity wave drag (m s-2)
  real(DP), allocatable:: xyz_DVDtGWD (:,:,:)
                              ! $ \DP{v}{t} $ . ���������� (m s-2)
                              ! Northward wind tendency by gravity wave drag (m s-2)

  real(DP), allocatable:: xyz_OMG (:,:,:)
                              ! Vertical velocity in pressure coordinate

  real(DP), allocatable:: xy_Rain (:,:)
                              ! ��水�. 
                              ! Precipitation
  real(DP), allocatable:: xy_RainCumulus  (:,:)
                              ! 
                              ! Rain due to moist convection
  real(DP), allocatable:: xy_RainLsc      (:,:)
                              ! 
                              ! Rain due to non-convective condensation

  real(DP), allocatable:: xy_Snow (:,:)
                              ! 
                              ! Snow fall
  real(DP), allocatable:: xy_SnowCumulus  (:,:)
                              ! 
                              ! Snow fall due to moist convection
  real(DP), allocatable:: xy_SnowLsc      (:,:)
                              ! 
                              ! Snow fall due to non-convective condensation

  real(DP), allocatable:: xyz_DTempDtCum(:,:,:)
  real(DP), allocatable:: xyz_DQVapDtCum(:,:,:)
  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_DQH2OSolDtCum(:,:,:)
                              ! Production rate of cloud ice in the layer 
                              ! due to condensation in cumulus convection 
                              ! parameterization (kg kg-1)
  real(DP), allocatable:: xyz_DUDtCum(:,:,:)
  real(DP), allocatable:: xyz_DVDtCum(:,:,:)

  real(DP), allocatable:: xyz_MoistConvDetTend       (:,:,:)
  real(DP), allocatable:: xyz_MoistConvSubsidMassFlux(:,:,:)

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

  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)
  real(DP), allocatable:: xyz_CloudCoverforRad(:,:,:)
                              ! Cloud cover (cloud fraction)

  real(DP), allocatable:: xy_SurfDustGravSedFlux(:,:)


  ! �業��
  ! Work variables
  !
  integer           :: IDDynMode                 ! 使�����������
                                                 ! Dynamics used for an experiment
  !
  integer, parameter:: IDDynModeHSPLVAS83       = 0
  integer, parameter:: IDDynModeNoHorAdv        = 1
  integer, parameter:: IDDynModeTWPICE          = 2

  integer           :: IDPhysMode                 ! 使������������
                                                 ! Physics used for an experiment
  !
  integer, parameter:: IDPhysModeNoPhysics       = 0
  integer, parameter:: IDPhysModeFullPhysics     = 1
  integer, parameter:: IDPhysModeHS94            = 2
  integer, parameter:: IDPhysModeVenusSimple     = 3
  integer, parameter:: IDPhysModeJupiterSimple   = 4
  integer, parameter:: IDPhysModeJupiterSimpleV2 = 5

  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:: IDRadMethodEarthV2     = 21
  integer, parameter:: IDRadMethodMarsV1      = 22
  integer, parameter:: IDRadMethodSL09        = 23
!!$  integer, parameter:: IDRadMethodVenusSimple = 24
  integer, parameter:: IDRadMethodSimple      = 25
  integer, parameter:: IDRadMethodRRTMG       = 26
  integer, parameter:: IDRadMethodNone        = 27


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


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

  integer           :: IDMCMethod                ! 湿潤対æ����ç®��¹æ�
                                                 ! Method for moist convection

  integer, parameter:: IDMCMethodNone        = 30
  integer, parameter:: IDMCMethodMCA         = 31
  integer, parameter:: IDMCMethodRAS         = 32
  integer, parameter:: IDMCMethodRASWithIce  = 33

  integer           :: IDLSCMethod               ! 大è�模å��çµ� (��対æ��§å��çµ�) ���ç®��¹æ�
                                                 ! Method for large scale condensation
                                                 ! (non-convective condensation)
  !
  integer, parameter:: IDLSCMethodNone        = 40
  integer, parameter:: IDLSCMethodM65         = 41
  integer, parameter:: IDLSCMethodLL91        = 42
  integer, parameter:: IDLSCMethodSatAdjM65   = 43
  integer, parameter:: IDLSCMethodM65WithIce  = 44
  integer, parameter:: IDLSCMethodLL91WithIce = 45

  integer           :: IDCloudMethod             ! �²ã���ç®��¹æ�
                                                 ! Method for cloud
  !
  integer, parameter:: IDCloudMethodNone          = 50
  integer, parameter:: IDCloudMethodSimple        = 51
  integer, parameter:: IDCloudMethodSimpleWithIce = 52
  integer, parameter:: IDCloudMethodT1993WithIce  = 53
  integer, parameter:: IDCloudMethodMarsH2OCloud  = 54

  integer           :: IDSfcMoistMethod          ! ����è¡��¢æ°´�����ç®��¹æ�
                                                 ! Method for surface moisture calculation
  !
  integer, parameter:: IDSfcMoistMethodNone   = 60
  integer, parameter:: IDSfcMoistMethodBucket = 61

  integer           :: IDDCMethod                ! ä¹¾ç�¥å�æµ����ç®��¹æ�
                                                 ! Method for dry convection
  !
  integer, parameter:: IDDCMethodNone = 70
  integer, parameter:: IDDCMethodDCA  = 71

  integer           :: IDGWDMethod               ! 
                                                 ! Method for gravity wave drag
  !
  integer, parameter:: IDGWDMethodNone  = 80
  integer, parameter:: IDGWDMethodM1987 = 81


  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

  integer:: i
  integer:: j
  integer:: k


  ! ���� ; 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 ( IDPhysMode )
    case ( IDPhysModeNoPhysics )

      xyz_DUDt     = 0.0_DP
      xyz_DVDt     = 0.0_DP
      xyz_DTempDt  = 0.0_DP
      xy_DPsDt     = 0.0_DP
      xyzf_DQMixDt = 0.0_DP

    case ( IDPhysModeHS94 )

      ! 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)
      xy_DPsDt     = 0.0_DP
      xyzf_DQMixDt = 0.0_DP

    case ( IDPhysModeVenusSimple )

      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_Temp, xyz_VirTemp, xyr_VirTemp, 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, xyz_VirTemp, xyr_VirTemp, xy_PsB, xyz_Press, xyr_Press, xyr_Temp, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyz_DUDt, xyz_DVDt, xyz_DTempDt )

      xy_DPsDt     = 0.0_DP
      xyzf_DQMixDt = 0.0_DP


    case ( IDPhysModeJupiterSimple )

      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_VirTemp   = xyr_VirTemp, xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xyz_Exner     = xyz_Exner, xy_SurfHeight = xy_SurfHeight, xyz_Height    = xyz_Height )


      ! Schneider and Liu (2009) ���������´æ··��課ç�
      ! Vertical diffusion by Schneider and Liu (2009)
      !
      call SL09Diffusion( xy_SurfHeight, xyz_Height, xyz_UB, xyz_VB, xyzf_QMixB, xyr_Press, xyr_VirTemp, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt )

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Flux( xyr_Press, xyz_Press, xyz_TempB, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      ! Net flux is calculated.
      !
      xyr_RadSFlux     = xyr_RadSUwFlux     - xyr_RadSDwFlux
      xyr_RadLFlux     = xyr_RadLUwFlux     - xyr_RadLDwFlux
      xyra_DelRadLFlux = xyra_DelRadLUwFlux - xyra_DelRadLDwFlux


      ! �¾å�������æ¸�º¦å¤�����
      ! Temperature tendency with radiation
      !
      call RadDTempDt( xyr_RadLFlux, xyr_RadSFlux, xyr_Press, xyz_DTempDtRadL, xyz_DTempDtRadS )          ! (out)
      xyr_RadLDwFlux     = 0.0_DP
      xyra_DelRadLDwFlux = 0.0_DP

      ! ���­ç�±å���±ç���������ç®�
      ! Sum all diabatic heating rates
      !
      xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS

      xy_DPsDt     = 0.0_DP


    case ( IDPhysModeJupiterSimpleV2 )

      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level,
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_VirTemp   = xyr_VirTemp, xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xy_SurfHeight = xy_SurfHeight, xyz_Height    = xyz_Height, xyz_Exner     = xyz_Exner,  xyr_Exner  = xyr_Exner )


      xyr_MomFluxX (:,:,1:kmax)   = 0.0_DP
      xyr_MomFluxY (:,:,1:kmax)   = 0.0_DP
      xyr_HeatFlux (:,:,1:kmax)   = 0.0_DP
      xyrf_QMixFlux(:,:,1:kmax,:) = 0.0_DP

      xyr_VelDiffCoef  = 0.0_DP
      xyr_TempDiffCoef = 0.0_DP
      xyr_QMixDiffCoef = 0.0_DP

      ! ����������������
      ! Lower boundary flux
      !
      call LBFluxSimple( xyz_UB, xyz_VB, xyz_TempB, xyr_VirTemp, xyzf_QMixB, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xyr_MomFluxX(:,:,0:0), xyr_MomFluxY(:,:,0:0), xyr_HeatFlux(:,:,0:0), xyrf_QMixFlux(:,:,0:0,:), xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef )

      call PhyImplAtmOnlyTendency( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt )

      xy_DPsDt             = 0.0_DP

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Flux( xyr_Press, xyz_Press, xyz_TempB, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      ! Net flux is calculated.
      !
      xyr_RadSFlux     = xyr_RadSUwFlux     - xyr_RadSDwFlux
      xyr_RadLFlux     = xyr_RadLUwFlux     - xyr_RadLDwFlux
      xyra_DelRadLFlux = xyra_DelRadLUwFlux - xyra_DelRadLDwFlux


      ! �¾å�������æ¸�º¦å¤�����
      ! Temperature tendency with radiation
      !
      call RadDTempDt( xyr_RadLFlux, xyr_RadSFlux, xyr_Press, xyz_DTempDtRadL, xyz_DTempDtRadS )          ! (out)
      xyr_RadLDwFlux     = 0.0_DP
      xyra_DelRadLDwFlux = 0.0_DP

      ! ���­ç�±å���±ç���������ç®�
      ! Sum all diabatic heating rates
      !
      xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS


    case ( IDPhysModeFullPhysics )

      ! �°è¡¨�¢æ�¡ä»¶��¨­å®�
      ! Configure surface conditions
      !
      call SetSurfaceProperties( xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xy_SurfTemp            = xy_SurfTemp, xy_SurfAlbedo          = xy_SurfAlbedo, xy_SurfHumidCoef       = xy_SurfHumidCoef, xy_SurfRoughLenMom     = xy_SurfRoughLenMom, xy_SurfRoughLenHeat    = xy_SurfRoughLenHeat, xy_SurfHeatCapacity    = xy_SurfHeatCapacity, xy_DeepSubSurfHeatFlux = xy_DeepSubSurfHeatFlux, xy_SurfCond            = xy_SurfCond, xy_SurfType            = xy_SurfType, xy_SeaIceConc          = xy_SeaIceConc, xy_SoilHeatCap         = xy_SoilHeatCap, xy_SoilHeatDiffCoef    = xy_SoilHeatDiffCoef, xy_SurfHeightStd       = xy_SurfHeightStd )

      ! ��, æ°·ã���²å��
      ! snow/ice fraction
      !
      call CalcSnowFrac( xy_SurfSnowB, xy_SnowFrac )

      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsB, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xy_SurfHeight = xy_SurfHeight, xy_SurfTemp = xy_SurfTemp, 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 RadDennouAGCMFlux( xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_Press, xy_SurfTemp, xy_SurfAlbedo, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )


      case ( IDRadMethodEarthV2 )

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

        ! 
        ! Cloud model
        !
        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = 0.0_DP
        case ( IDCloudMethodSimple )
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = 0.0_DP
        case ( IDCloudMethodSimpleWithIce )
          if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
            call MessageNotify( 'E', prog_name, 'H2OSol is not found.' )
          end if
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OSol' ))
        case ( IDCloudMethodT1993WithIce )
          if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
            call MessageNotify( 'E', prog_name, 'H2OSol is not found.' )
          end if
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OSol' ))
        end select

!!$        call CloudSimpleDivideWatAndIce(                           &
!!$          & xyz_TempB,                                             & ! (in )
!!$          & xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' )),     & ! (in )
!!$          & xyz_QH2OLiqforRad, xyz_QH2OSolforRad                   & ! (out)
!!$          & )

        ! Cloud cover is calculated.
        !
        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          xyz_CloudCoverforRad = 1.0_DP
        case ( IDCloudMethodSimple )
          call CloudSimpleCalcCloudCover( xyz_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_CloudCoverforRad )
        case ( IDCloudMethodSimpleWithIce )
          call CloudSimpleCalcCloudCover( xyz_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_CloudCoverforRad )
        case ( IDCloudMethodT1993WithIce )
          xyz_CloudCoverforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'CloudCover' ))
        end select

        call RadEarthV2Flux( xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_QH2OLiqforRad, xyz_QH2OSolforRad, xyz_CloudCoverforRad, xy_SurfTemp, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      case ( IDRadMethodRRTMG )

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

        ! 
        ! Cloud model
        !
        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = 0.0_DP
        case ( IDCloudMethodSimple )
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = 0.0_DP
        case ( IDCloudMethodSimpleWithIce )
          if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
            call MessageNotify( 'E', prog_name, 'H2OSol is not found.' )
          end if
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OSol' ))
        case ( IDCloudMethodT1993WithIce )
          if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
            call MessageNotify( 'E', prog_name, 'H2OSol is not found.' )
          end if
          xyz_QH2OLiqforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OLiq' ))
          xyz_QH2OSolforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'H2OSol' ))
        end select

        ! Cloud cover is calculated.
        !
        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          xyz_CloudCoverforRad = 1.0_DP
        case ( IDCloudMethodSimple )
          call CloudSimpleCalcCloudCover( xyz_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_CloudCoverforRad )
        case ( IDCloudMethodSimpleWithIce )
          call CloudSimpleCalcCloudCover( xyz_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_CloudCoverforRad )
        case ( IDCloudMethodT1993WithIce )
          xyz_CloudCoverforRad = xyzf_QMixB(:,:,:,CompositionInqIndex( 'CloudCover' ))
        end select


        ! wrapper of RRTMG
        ! wrapper of RRTMG
        !
        call RadRRTMGWrapperFlux( xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_TempB, xyr_Temp, xyzf_QMixB(:,:,:,IndexH2OVap), xyz_QH2OLiqforRad, xyz_QH2OSolforRad, xyz_CloudCoverforRad, xy_SurfTemp, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      case ( IDRadMethodMarsV1 )

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

        call RadMarsV1Flux( xy_SurfType, xy_SurfMajCompIceB, xy_SurfAlbedo, xyz_Press, xyr_Press, xyz_TempB, xyr_Temp, xy_SurfTemp, xyzf_QMixB(:,:,:,CompositionInqIndex( 'Dust' )), xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      case ( IDRadMethodSL09 )

        ! Schneider and Liu (2009) ���¾å��¢ã����
        ! Radiation model by Schneider and Liu (2009)
        !
        call RadSL09Flux( xyr_Press, xyz_Press, xyz_TempB, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      case ( IDRadMethodSimple )

        ! ç°¡å���¾å��¢ã����
        ! Simple radiation model
        !
        call RadSimpleFlux( xy_SurfAlbedo, xy_SurfTemp, xyr_Press, xyz_Press, xyz_TempB, xyzf_QMixB(:,:,:,IndexH2OVap), xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      case ( IDRadMethodNone )

        ! ä½����������¾å��¢ã����
        ! radiation model with no absorption and no scattering
        !
        call RadNoneFlux( xy_SurfAlbedo, xy_SurfTemp, xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux )

      end select

      ! Net flux is calculated.
      !
      xyr_RadSFlux     = xyr_RadSUwFlux     - xyr_RadSDwFlux
      xyr_RadLFlux     = xyr_RadLUwFlux     - xyr_RadLDwFlux
      xyra_DelRadLFlux = xyra_DelRadLUwFlux - xyra_DelRadLDwFlux


      ! �°è¡¨�¢ã����������
      ! Surface flux
      !
      select case ( IDSfcFluxMethod )
      case ( IDSfcFluxMethodL82 )

        call SurfaceFlux( 'L82', xyz_UB, xyz_VB, xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, xyzf_QMixB, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfRoughLenMom, xy_SurfRoughLenHeat, xy_SnowFrac, xy_SurfMomFluxX, xy_SurfMomFluxY, xy_SurfHeatFlux, xyf_SurfQMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xy_SurfMOLength )

      case ( IDSfcFluxMethodBH91B94 )

        call SurfaceFlux( 'BH91B94', xyz_UB, xyz_VB, xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xy_SurfVirTemp, xyzf_QMixB, xyr_Press, xy_SurfHeight, xyz_Height, xyz_Exner, xyr_Exner, xy_SurfTemp, xy_SurfHumidCoef, xy_SurfRoughLenMom, xy_SurfRoughLenHeat, xy_SnowFrac, xy_SurfMomFluxX, xy_SurfMomFluxY, xy_SurfHeatFlux, xyf_SurfQMixFlux, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xy_SurfMOLength )

      end select
      !
      ! set dust flux
      !   This is ad hoc treatment now (yot, 2013/09/28)
      !
!!$      if ( CompositionInqIndex('Dust') > 0 ) then
!!$        call SetDustFlux(                                       &
!!$          & xyf_SurfQMixFlux(:,:,CompositionInqIndex('Dust'))   & ! (out)
!!$          & )
!!$      end if


      ! ���´æ�¡æ�£ã����������
      ! Vertical diffusion flux
      !
      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY2 )

        call VDiffusion( xyz_UB,     xyz_VB,     xyzf_QMixB, xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner,    xyr_Exner, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )

      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, xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xyzf_QMixB(:,:,:,CompositionInqIndex( 'TKE' )), xy_SurfMomFluxX, xy_SurfMomFluxY, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyz_DTurKinEneDt )

      case ( IDVDiffMethodJMA )

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

        ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
        ! JMA turbulent mixing module
        !
        call VDiffusionJMAMYWrapper3D( xyz_UB,     xyz_VB,     xyzf_QMixB, xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyz_Press, xyr_Press, xy_SurfHeight, xyz_Height, xyr_Height, xyz_Exner, xyr_Exner, xy_SurfMOLength, xyzf_QMixB(:,:,:,CompositionInqIndex( 'TKE' )), xy_SurfMomFluxX, xy_SurfMomFluxY, xy_SurfHeatFlux, xyf_SurfQMixFlux, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xyz_DTurKinEneDt )

      end select

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


      ! ä¸������������������å¤��������ç®� (�°è§£æ³�)
      ! 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_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt, xy_DSurfTempDt )

        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 SubsurfaceDiffusionFlagLandMatthews(  &
!!$          & xy_SurfType,                           & ! (in)
!!$          & xy_FlagMatthewsLand                    & ! (out)
!!$          & )
        call SubsurfaceDiffusion( xy_DeepSubSurfHeatFlux, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xy_SurfTemp, xyz_SoilTemp, xy_SurfSnowB, xyr_SoilTempTransCoef, xyr_SoilHeatFlux )


        ! This value is not correct, if snow evaporates. 
        ! If a bucket model is used, this surface latent heat flux is 
        ! corrected in BucketModEvapAndLatentHeatFlux by the use of surface 
        ! moisture and surface snow amount. 
        !
!!$        xy_SurfLatentHeatFluxA = LatentHeat * xyrf_QMixFlux(:,:,0,IndexH2OVap)


!!$        select case ( IDSfcMoistMethod )
!!$        case ( IDSfcMoistMethodBucket )
!!$          ! ���±ã���¢ã�������������°è¡¨�¢ã���������¹ä¿®æ­�
!!$          ! Modification of surface flux for bucket model
!!$          !
!!$          call BucketSetFlagOceanFromMatthews( &
!!$            & xy_SurfType,                     & ! (in)
!!$            & xy_BucketFlagOceanGrid           & ! (out)
!!$            & )
!!$          call BucketModEvapAndLatentHeatFlux(                         &
!!$            & xy_BucketFlagOceanGrid, xy_SoilMoistB, xy_SurfSnowB,     & ! (in   )
!!$            & xyrf_QMixFlux(:,:,0,IndexH2OVap), xy_SurfLatentHeatFluxA & ! (inout)
!!$            & )
!!$        end select



!!$        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)
!!$          & )
!!$        select case ( IDSfcMoistMethod )
!!$        case ( IDSfcMoistMethodBucket )
!!$          call BucketSetFlagOceanFromMatthews( &
!!$            & xy_SurfType,                     & ! (in)
!!$            & xy_BucketFlagOceanGrid           & ! (out)
!!$            & )
!!$        case default
!!$          xy_BucketFlagOceanGrid = .true.
!!$        end select
!!$        call PhyImplSDHV2Tendency(                                   &
!!$          & xy_PhyImplSDHIndexCalcMethod, xy_BucketFlagOceanGrid,    & ! (in)
!!$          & xy_SnowFrac,                                             & ! (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)
!!$          & xyzf_QMixB,                                              & ! (in)
!!$          & xy_SurfHumidCoef,                                        & ! (in)
!!$          & xy_SurfHeatCapacity,                                     & ! (in)
!!$          & xy_SoilHeatCap, xy_SoilHeatDiffCoef,                     & ! (in)
!!$          & xyra_DelRadLFlux,                                        & ! (in)
!!$          & xyr_Press, xyz_Exner, xyr_Exner,                         & ! (in)
!!$          & xyr_VirTemp, xyz_Height,                                 & ! (in)
!!$          & xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef,     & ! (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)
!!$          & )

        call PhyImplSDHV3SetMethodMatthews( xy_SurfCond, xy_SurfType, xy_SeaIceConc, xy_PhyImplSDHIndexCalcMethod )
        select case ( IDSfcMoistMethod )
        case ( IDSfcMoistMethodBucket )
          call BucketSetFlagOceanFromMatthews( xy_SurfType, xy_BucketFlagOceanGrid )
        case default
          xy_BucketFlagOceanGrid = .true.
        end select
        call PhyImplSDHV3Tendency( xy_PhyImplSDHIndexCalcMethod, xy_BucketFlagOceanGrid, xy_SnowFrac, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA, xyr_SoilHeatFlux, xyr_RadSFlux, xyr_RadLFlux, xy_DeepSubSurfHeatFlux, xyz_TempB, xy_SurfTemp, xyz_SoilTemp, xyzf_QMixB, xy_SurfHumidCoef, xy_SoilHeatCap, xy_SoilHeatDiffCoef, xyra_DelRadLFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, xy_SurfVelTransCoef, xy_SurfTempTransCoef, xy_SurfQVapTransCoef, xyr_SoilTempTransCoef, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xyz_DUDt, xyz_DVDt, xyz_DTempDtVDiff, xyzf_DQMixDt, xy_DSurfTempDt, xyz_DSoilTempDt, xy_DPsDt, xy_DSurfMajCompIceDt, xy_DSoilMoistDt, xy_DSurfSnowDt )

      case ( IDPhyTendMethodImpAtmOnly )

        call PhyImplAtmOnlyTendency( xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef, 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_DPsDt             = 0.0_DP
        xy_DSurfMajCompIceDt = 0.0_DP
        xy_DSoilMoistDt      = 0.0_DP
        xy_DSurfSnowDt       = 0.0_DP

        xy_SurfH2OVapFluxA     = 0.0_DP
        xy_SurfLatentHeatFluxA = 0.0_DP

      end select


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


      ! Gravity wave drag
      ! Gravity wave drag
      !
      select case ( IDGWDMethod )
      case ( IDGWDMethodNone )

        xyz_DUDtGWD = 0.0_DP
        xyz_DVDtGWD = 0.0_DP

      case ( IDGWDMethodM1987 )

        ! Gravity wave drag by McFarlane (1987)
        ! Gravity wave drag by McFarlane (1987)
        !
        call GWDM1987( xyz_UB, xyz_VB, xyz_TempB, xyz_Press, xyr_Press, xyz_Exner, xyz_Height, xy_SurfHeight, xy_SurfHeightStd, xyz_DUDtGWD, xyz_DVDtGWD )
        xyz_DUDt = xyz_DUDt + xyz_DUDtGWD
        xyz_DVDt = xyz_DVDt + xyz_DVDtGWD

      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 RadDTempDt( xyr_RadLFluxA, xyr_RadSFlux, xyr_Press, xyz_DTempDtRadL, xyz_DTempDtRadS )          ! (out)

!!$      select case ( IDRadMethod )
!!$      case ( IDRadMethodMarsV1 )
!!$        ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$        ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$        !
!!$        call rad15mNLTEMergeHR(                &
!!$          & xyz_Press, xyz_TempB, xyz_VirTemp, &
!!$          & xyz_DTempDtRadL                    &
!!$          & )
!!$      end select

      ! ���­ç�±å���±ç���������ç®�
      ! Sum all diabatic heating rates
      !
      xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS


      select case ( IDRadMethod )
      case ( IDRadMethodMarsV1 )
        ! ����è¨�ç®����赤å����±ç�è¨�ç®�
        ! Calculation of near infrared heating rate in the case of Mars
        !
        call RadMarsNIRINOUT( xyz_Press, xyz_DTempDt )
      end select


      ! ���´æ�¡æ�£ã���������¹ã���ºå�� 
      !   * �ºå�����¿ã���µã�����¼ã���³ã�§ã����, è¨�ç®�����½±�¿ã������
      ! 
      ! 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, xyr_Press, xyz_Exner, xyr_Exner, xyr_VirTemp, xyz_Height, xyr_VelDiffCoef, xyr_TempDiffCoef, xyr_QMixDiffCoef )

      ! �°è¡¨�¢ã���������¹ã���ºå�� 
      !   * �ºå�����¿ã���µã�����¼ã���³ã�§ã����, è¨�ç®�����½±�¿ã������
      ! 
      ! Output surface fluxes
      !   * This subroutine works for output only, 
      !     so it does not influence a calculation.
      !
      call SurfaceFluxOutput( xy_SnowFrac, xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, xy_SurfH2OVapFluxA, xy_SurfLatentHeatFluxA, 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 RadFluxOutput( xyr_RadSUwFlux, xyr_RadSDwFlux, xyr_RadLUwFlux, xyr_RadLDwFlux, xyra_DelRadLUwFlux, xyra_DelRadLDwFlux, xy_DSurfTempDt, xyz_DTempDtVDiff )

    end select


    ! �°é�¢æ¸©åº��»å��å£�æ¸�º¦�»å��å£�æ°´å���»ç���������
    ! Time integration of surface temperature, soil temperature, soil 
    ! moisture, and surface snow amount
    !
    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      ! �°é�¢æ¸©åº��»å��å£�æ¸�º¦������ç©���
      ! Time integration of surface temperature and soil temperature
      !
      call IntegralSurfTemp( xy_DSurfTempDt, xyz_DSoilTempDt, xy_SurfTemp   , xyz_SoilTemp )

      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodNone )
        xy_SoilMoistA = xy_SoilMoistB
        xy_SurfSnowA  = xy_SurfSnowB
      case ( IDSfcMoistMethodBucket )
        ! ��å£�æ°´å���»å�°é�¢ç����������ç©���
        ! Time integration of soil moisture and snow amount
        !
        call BucketSetFlagOceanFromMatthews( xy_SurfType, xy_BucketFlagOceanGrid )
        call BucketIntegration( xy_BucketFlagOceanGrid, xy_DSoilMoistDt, xy_DSurfSnowDt, xy_SoilMoistB, xy_SurfSnowB, xy_SoilMoistA, xy_SurfSnowA )
      end select

      xy_SurfMajCompIceA = xy_SurfMajCompIceB + xy_DSurfMajCompIceDt * ( 2.0_DP * DelTime )

    end select


    ! ������
    ! Dynamical core
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      call DynamicsHSplVAS83( 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, xyz_ArgOMG = xyz_OMG )
    case ( IDDynModeNoHorAdv )
      call DynamicsPhysicsOnly( xyz_Exner, xy_SurfHeight, xyz_Height, xyz_DUDt, xyz_DVDt, xyz_DTempDt, xyzf_DQMixDt, xy_PsB, xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsN, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsA, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA )
      xyz_OMG = 0.0_DP
    case ( IDDynModeTWPICE )
      !
      ! Dynamical process for TWP-ICE experiment
      !
      call DynamicsTWPICESCMExp( TimeN, xyz_Press, xyz_Exner, xyz_Height, 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 )
      xyz_OMG = 0.0_DP
    end select


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


    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimple )
      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press = xyz_Press, xyr_Press = xyr_Press )

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjust( xyz_TempA, xyz_UA, xyz_VA, xyzf_QMixA, xyz_Press, xyr_Press )
      end select

    case ( IDPhysModeJupiterSimpleV2 )
      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press = xyz_Press, xyr_Press = xyr_Press )

      ! 湿潤対æ�調ç�
      ! Moist convective adjustment
      !
      call MoistConvAdjust( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press, xyr_Press, xyz_DQH2OLiqDtCum )

      ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
      ! Large scale condensation (non-convective condensation) 
      ! (Manabe, 1965)
      !
      call LScaleCond( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press, xyr_Press, xyz_DQH2OLiqDtLSC )

      call CloudSimpleCalcPRCPKeyLLTemp3D( xyr_Press, xyz_Press, xyz_DQH2OLiqDtCum, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_RainCumulus, xy_SnowCumulus )
      call CloudSimpleCalcPRCPKeyLLTemp3D( xyr_Press, xyz_Press, xyz_DQH2OLiqDtLsc, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_RainLsc, xy_SnowLsc )

      xy_Rain = xy_RainCumulus + xy_RainLsc
      xy_Snow = xy_SnowCumulus + xy_SnowLsc

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjust( xyz_TempA, xyz_UA, xyz_VA, xyzf_QMixA, xyz_Press, xyr_Press )
      end select

      ! ������³ª����æ­�
      ! Fix masses of constituents
      !
      call MassFixerColumn( xyr_Press, xyzf_QMixA )

    case ( IDPhysModeFullPhysics )

      ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
      ! Interpolate temperature on half sigma level, 
      ! and calculate pressure and height
      !
      call AuxVars( xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xyz_Exner     = xyz_Exner, xyr_Exner = xyr_Exner, xy_SurfHeight = xy_SurfHeight, xyr_Height    = xyr_Height, xyz_Height    = xyz_Height )


      ! ç©��²ã�����¡ã�¿ã���¼ã�¼ã�·ã�§ã��
      ! Cumulus parameterization
      !
      select case ( IDMCMethod )
      case ( IDMCMethodNone )

        xyz_DQH2OLiqDtCum = 0.0_DP
        xyz_DQH2OSolDtCum = 0.0_DP

        xyz_MoistConvDetTend        = 0.0_DP
        xyz_MoistConvSubsidMassFlux = 0.0_DP

        xy_RainCumulus = 0.0_DP
        xy_SnowCumulus = 0.0_DP

      case ( IDMCMethodMCA )

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

        ! 湿潤対æ�調ç�
        ! Moist convective adjustment
        !
        call MoistConvAdjust( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press, xyr_Press, xyz_DQH2OLiqDtCum )

        ! It would be better that lines below are included in 
        ! MoistConvAdjust subroutine.
        xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) = xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) + xyz_DQH2OLiqDtCum * 2.0_DP * DelTime
        xyz_DQH2OLiqDtCum = 0.0_DP
        xyz_DQH2OSolDtCum = 0.0_DP

        xyz_MoistConvDetTend        = 0.0_DP
        xyz_MoistConvSubsidMassFlux = 0.0_DP

        xy_RainCumulus = 0.0_DP
        xy_SnowCumulus = 0.0_DP

      case ( IDMCMethodRAS )

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

        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        call RAS( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_DQH2OLiqDtCum, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux )

        ! It would be better that lines below are included in 
        ! RAS subroutine.
        xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) = xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) + xyz_DQH2OLiqDtCum * 2.0_DP * DelTime
        xyz_DQH2OLiqDtCum = 0.0_DP
        xyz_DQH2OSolDtCum = 0.0_DP

        xy_RainCumulus = 0.0_DP
        xy_SnowCumulus = 0.0_DP

      case ( IDMCMethodRASWithIce )
        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OSol is not found.' )
        end if

        ! It should be noted that H2OLiq and H2OSol have updated. 
        ! But, DQH2OLiqDtCum and DQH2OSolDtCum are not zero, so those
        ! will be set zero below.
        call RASWithIce1DWrapper3DWrapper( xy_SurfTemp, xyz_Press, xyr_Press, xyz_Exner, xyr_Exner, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')), xyzf_QMixA(:,:,:,CompositionInqIndex('H2OSol')), xyz_UA, xyz_VA, xyz_DTempDtCum, xyz_DQVapDtCum, xyz_DQH2OLiqDtCum, xyz_DQH2OSolDtCum, xyz_DUDtCum, xyz_DVDtCum, xy_RainCumulus, xy_SnowCumulus, xyz_MoistConvDetTend, xyz_MoistConvSubsidMassFlux )

        xyz_TempA = xyz_TempA + xyz_DTempDtCum * ( 2.0_DP * DelTime )
        xyzf_QMixA(:,:,:,IndexH2OVap) = xyzf_QMixA(:,:,:,IndexH2OVap) + xyz_DQVapDtCum * ( 2.0_DP * DelTime )
        xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')) = xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')) + xyz_DQH2OLiqDtCum * ( 2.0_DP * DelTime )
        xyzf_QMixA(:,:,:,CompositionInqIndex('H2OSol')) = xyzf_QMixA(:,:,:,CompositionInqIndex('H2OSol')) + xyz_DQH2OSolDtCum * ( 2.0_DP * DelTime )
        xyz_UA = xyz_UA + xyz_DUDtCum * ( 2.0_DP * DelTime )
        xyz_VA = xyz_VA + xyz_DVDtCum * ( 2.0_DP * DelTime )

        ! ������³ª����æ­�
        ! Fix masses of constituents
        !
        call MassFixerColumn( xyr_Press, xyzf_QMixA )

        xyz_DTempDtCum    = 0.0_DP
        xyz_DQVapDtCum    = 0.0_DP
        xyz_DQH2OLiqDtCum = 0.0_DP
        xyz_DQH2OSolDtCum = 0.0_DP
        xyz_DUDtCum       = 0.0_DP
        xyz_DVDtCum       = 0.0_DP


      end select


      ! 大è�模å��çµ� (��対æ��§å��çµ�)
      ! Large scale condensation
      !
      select case ( IDLSCMethod )
      case ( IDLSCMethodNone )

        xyz_DQH2OLiqDtLSC = 0.0_DP
        xyz_DQH2OSolDtLSC = 0.0_DP

      case ( IDLSCMethodM65 )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) 
        ! (Manabe, 1965)
        !
        call LScaleCond( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_Press, xyr_Press, xyz_DQH2OLiqDtLSC )

        ! It would be better that lines below are included in LScaleCond
        ! subroutine.
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) = xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) + xyz_DQH2OLiqDtLSC * 2.0_DP * DelTime
        xyz_DQH2OLiqDtLSC = 0.0_DP
        xyz_DQH2OSolDtLSC = 0.0_DP

!!$      case ( IDLSCMethodLL91 )
!!$          ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
!!$          ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
!!$          !
!!$          call LScaleCondLL91(                                 &
!!$            & xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap),        &  ! (inout)
!!$            & xyz_DTempDtCond, xyz_DQVapDtCond,                &  ! (inout)
!!$            & xyz_Press, xyr_Press,                            &  ! (in)
!!$            & xyz_DQH2OLiqDtLSC                                &  ! (out)
!!$            & )
      case ( IDLSCMethodSatAdjM65 )
        !
        != Saturation adjustment
        !
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        call SaturationAdjust( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )), xyz_Press, xyr_Press, xyz_DQH2OLiqDtLSC )

      case ( IDLSCMethodM65WithIce )
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if

        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) 
        ! (Manabe, 1965)
        !
        ! It should be noted that H2OLiq and H2OSol have updated in above 
        ! subroutine.

        call LScaleCond1D3DWrapper( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OSol' )), xyz_Press, xyr_Press, xyz_DQH2OLiqDtLSC, xyz_DQH2OSolDtLSC )

        ! Temporal treatment
        xyz_DQH2OLiqDtLSC = 0.0_DP
        xyz_DQH2OSolDtLSC = 0.0_DP

      case ( IDLSCMethodLL91WithIce )
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if

        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        ! It should be noted that H2OLiq and H2OSol have updated in above 
        ! subroutine.

        call LScaleCondLL911D3DWrapper( xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OSol' )), xyz_Press, xyr_Press, xyz_DQH2OLiqDtLSC, xyz_DQH2OSolDtLSC )

        ! Temporal treatment
        xyz_DQH2OLiqDtLSC = 0.0_DP
        xyz_DQH2OSolDtLSC = 0.0_DP

      end select


      ! 
      ! Cloud model
      !
      select case ( IDCloudMethod )
      case ( IDCloudMethodNone )

        ! �²ã�����¢ã����
        ! No cloud model
        !
        if ( CompositionInqIndex( 'H2OSol' ) <= 0 ) then
          if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
            xy_RainLsc = 0.0_DP
            xy_SnowLsc = 0.0_DP
          else
            call CloudNone( xyr_Press, xyz_Press, xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) / ( 2.0_DP * DelTime ), xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_RainLsc )
            xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) = 0.0_DP
            xy_SnowLsc = 0.0_DP
          end if
        else
          call CloudNoneWithIce( xyr_Press, xyz_Press, xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) / ( 2.0_DP * DelTime ), xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OSol' )) / ( 2.0_DP * DelTime ), xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xy_RainLsc, xy_SnowLsc )
          xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OLiq' )) = 0.0_DP
          xyzf_QMixA(:,:,:,CompositionInqIndex( 'H2OSol' )) = 0.0_DP
        end if

      case ( IDCloudMethodSimple )

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

        ! Update cloud water
        !
        call CloudSimple( xyr_Press, xyz_Press, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')), xy_RainLsc, xy_SnowLsc )

      case ( IDCloudMethodSimpleWithIce )

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

        ! Update cloud water
        !
        call CloudSimpleWithIce( xyr_Press, xyz_Press, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap)                  , xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')), xyzf_QMixA(:,:,:,CompositionInqIndex('H2OSol')), xy_RainLsc, xy_SnowLsc )

!!$        call PhyImplSDHV2SetMethodMatthews(  &
!!$          & xy_SurfType, xy_SeaIceConc,      & ! (in)
!!$          & xy_PhyImplSDHIndexCalcMethod     & ! (out)
!!$          & )
!!$        call PhyImplSDHV2CorSOTempBySnowMelt(  &
!!$          & xy_PhyImplSDHIndexCalcMethod,    & ! (in   )
!!$          & xy_Snow,                         & ! (in   )
!!$          & xy_SurfTemp                      & ! (inout)
!!$          & )

      case ( IDCloudMethodT1993WithIce )

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

        ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
        ! Interpolate temperature on half sigma level, 
        ! and calculate pressure and height
        !
        call AuxVars( xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap), xyr_Press     = xyr_Press, xyz_Press     = xyz_Press, xyz_VirTemp   = xyz_VirTemp )

        ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
        ! Cloud model based on Tiedtke (1993)
        !
        call CloudT1993baseWithIce( xyz_Press, xyr_Press, xyz_VirTemp, xyz_DQH2OLiqDtCum, xyz_DQH2OSolDtCum, xyz_MoistConvDetTend, xyz_OMG, xyz_MoistConvSubsidMassFlux, xyz_DTempDt, xyz_TempA, xyzf_QMixA(:,:,:,IndexH2OVap)                  , xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')), xyzf_QMixA(:,:,:,CompositionInqIndex('H2OSol')), xyzf_QMixA(:,:,:,CompositionInqIndex('CloudCover')), xy_RainLsc, xy_SnowLsc )

      case ( IDCloudMethodMarsH2OCloud )

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

        ! æ¸�º¦�����´æ�°Ï������������, æ°��§ã���åº������
        ! Interpolate temperature on half sigma level, 
        ! and calculate pressure and height
        !
        call AuxVars( xy_PsB, xyz_TempB, xyzf_QMixA(:,:,:,IndexH2OVap), xyz_VirTemp = xyz_VirTemp, xyr_Press     = xyr_Press, xy_SurfHeight = xy_SurfHeight, xyr_Height = xyr_Height )

        ! ���� H2O �²ã�¢ã����
        ! Mars H2O cloud model
        !
        call CloudMarsH2O( xyr_Press, xyz_VirTemp, xyr_Height, xyz_DQH2OLiqDtCum, xyz_DQH2OLiqDtLSC, xyzf_QMixA(:,:,:,CompositionInqIndex('H2OLiq')), xy_RainLsc, xy_SnowLsc )

      end select


      ! Sum of cumulus and non-convective (large scale) condensation
      xy_Rain = xy_RainCumulus + xy_RainLsc
      xy_Snow = xy_SnowCumulus + xy_SnowLsc


      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodBucket )
        ! ���±ã���¢ã����, ��æ°´ã��¼´���°è¡¨æ°´é�鎀������ç®�
        ! bucket model, calculation of change of surface water due to 
        ! precipitation
        !
        call BucketSetFlagOceanFromMatthews( xy_SurfType, xy_BucketFlagOceanGrid )
        call BucketPRCPAdjust( xy_BucketFlagOceanGrid, xy_Rain, xy_Snow, xy_SoilMoistA, xy_SurfSnowA )
      end select


      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjust( xyz_TempA, xyz_UA, xyz_VA, xyzf_QMixA, xyz_Press, xyr_Press )
      end select



      ! ����������
      ! Gravitational sedimentation process
      !
      if ( CompositionInqIndex('Dust') > 0 ) then
        call GravSed( 'MarsDust', xyr_Press, xyz_VirTemp, xyr_Height, xyzf_QMixA(:,:,:,CompositionInqIndex('Dust')), xy_SurfDustGravSedFlux )
      end if

!!$      if ( IDPhysMode == IDPhysModeFullPhysics ) then

        ! 主æ�����¸å���
        ! Phase change of atmospheric major component
        !
        call MajorCompPhaseChangeInAtm( xyr_Press, xyz_Press, xy_PsA, xyz_TempA, xy_SurfMajCompIceA )

!!$      end if

      ! ������³ª����æ­�
      ! Fix masses of constituents
      !
      call MassFixerColumn( xyr_Press, xyzf_QMixA )

    end select


    ! �������£ã���¿ã�� (Asselin, 1972)
    ! Time filter (Asselin, 1972)
    !
!!$    if ( .not. flag_initial .or. .not. firstloop ) then
!!$      call TimeFilter( &
!!$        & 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  )   ! (in)
!!$
!!$      select case ( IDPhysMode )
!!$      case ( IDPhysModeFullPhysics )
!!$        call TimeFilterSurfVars(                             &
!!$          & xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, &   ! (in)
!!$          & xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN, &   ! (inout)
!!$          & xy_SurfMajCompIceA, xy_SoilMoistA, xy_SurfSnowA  &   ! (in)
!!$          & )
!!$      end select
!!$    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, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA )

      select case ( IDPhysMode )
      case ( IDPhysModeFullPhysics )
        call TimeFilterWilliams2009SurfVars( xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN, xy_SurfMajCompIceA, xy_SoilMoistA, xy_SurfSnowA, xy_PsA )
      end select
    end if


    ! äº��±å��°ã���¤ã��¢ºèª�
    ! 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 )

    ! Output frequently used variables
    ! Output frequently used variables
    !
    call OutputFreqUsedVars( xy_PsA, xyz_TempA, xyzf_QMixA, xy_SurfHeight )

    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimpleV2 )
      call HistoryAutoPut( TimeN, 'Rain', xy_Rain )
      call HistoryAutoPut( TimeN, 'Snow', xy_Snow )
      call HistoryAutoPut( TimeN, 'PRCP', xy_Rain + xy_Snow )

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

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

!!$      call HistoryAutoPut( TimeN, 'Rain', xy_Rain * LatentHeat )
!!$      call HistoryAutoPut( TimeN, 'Rain', ( xy_Rain + xy_Snow ) * LatentHeat )
      call HistoryAutoPut( TimeN, 'Rain', xy_Rain )
      call HistoryAutoPut( TimeN, 'Snow', xy_Snow )
      call HistoryAutoPut( TimeN, 'PRCP', xy_Rain + xy_Snow )
      call HistoryAutoPut( TimeN, 'RainCum', xy_RainCumulus )
      call HistoryAutoPut( TimeN, 'SnowCum', xy_SnowCumulus )
      call HistoryAutoPut( TimeN, 'PRCPCum', xy_RainCumulus + xy_SnowCumulus )
      call HistoryAutoPut( TimeN, 'RainLsc', xy_RainLsc )
      call HistoryAutoPut( TimeN, 'SnowLsc', xy_SnowLsc )
      call HistoryAutoPut( TimeN, 'PRCPLsc', xy_RainLsc + xy_SnowLsc )

      call HistoryAutoPut( TimeN, 'SeaIceConc'      , xy_SeaIceConc       )
      call HistoryAutoPut( TimeN, 'SurfAlbedo'      , xy_SurfAlbedo       )
      call HistoryAutoPut( TimeN, 'SurfRoughLenMom' , xy_SurfRoughLenMom  )
      call HistoryAutoPut( TimeN, 'SurfRoughLenHeat', xy_SurfRoughLenHeat )

      if ( CompositionInqIndex('Dust') > 0 ) then
        call HistoryAutoPut( TimeN, 'SurfDustGravSedFlux', xy_SurfDustGravSedFlux )
      end if
    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 ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      xy_SurfMajCompIceB = xy_SurfMajCompIceN
      xy_SurfMajCompIceN = xy_SurfMajCompIceA
      xy_SurfMajCompIceA = 0.0_DP

      xy_SoilMoistB = xy_SoilMoistN
      xy_SoilMoistN = xy_SoilMoistA
      xy_SoilMoistA = 0.0_DP

      xy_SurfSnowB  = xy_SurfSnowN
      xy_SurfSnowN  = xy_SurfSnowA
      xy_SurfSnowA  = 0.0_DP
    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, xy_PsB, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN )

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿å�ºå��
      ! Restart data of surface temperature output
      !
      call RestartSurfTempOutput( xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN )
    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

    ! ���»ç���
    ! 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: RestartFileInit, RestartFileOpen, RestartFileGet

    ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempInit, 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

    ! ��ç³»å�����¼ã�¿ã����¿è¾¼��
    ! 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

    ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only : HS94Init

    ! Yamamoto and Takahashi (2003) ����£ã��ç°¡å������è¨�ç®���������¼·��
    ! forcing for simple Venus calculation following Yamamoto and Takahashi (2003)
    !
    use yt2003_forcing, only : YT2003ForcingInit

    ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
    ! Setting planetary surface properties
    !
    use surface_properties, only : SurfacePropertiesInit

    ! ��, æ°·ã���²å��
    ! snow/ice fraction
    !
    use snowice_frac, only : SnowIceFracInit

    ! �°ä����������±ã�����´æ�¡æ��
    ! 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
    ! �°ä��±ä�å°��¢ã�������������´å�����°è§£æ³�����������ç©���
    !
    ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
    use phy_implicit_sdh_V3, only : PhyImplSDHV3Init

    ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
    ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
    !
    use phy_implicit_atmonly, only : PhyImplAtmOnlyInit

    ! �°è§£æ³�����������ç©��������������¼ã����
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplUtilsInit

    ! ���±ã���¢ã����
    ! Bucket model
    !
    use Bucket_Model, only : BucketModelInit

    ! �¾å����������� (GFD �»è�³å�¶æ¥½�����ºã���¾å��¢ã����)
    ! Radiation flux (radiation model developed by GFD Dennou Club)
    !
    use rad_DennouAGCM, only : RadDennouAGCMInit

    ! �°ç��大æ������¾å��¢ã���� Ver. 2
    ! radiation model for the Earth's atmosphere Ver. 2
    !
    use rad_Earth_V2, only : RadEarthV2Init

    ! wrapper of RRTMG
    ! wrapper of RRTMG
    !
    use rad_rrtmg_wrapper, only : RadRRTMGWrapperInit

    ! ����大æ������¾å��¢ã���� Ver. 1
    ! radiation model for the Mars' atmosphere Ver. 1
    !
    use rad_Mars_V1, only : RadMarsV1Init

!!$    ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$    ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$    !
!!$    use rad_15m_NLTE, only: Rad15mNLTEInit

    ! ����è¨�ç®����赤å����±ç�è¨�ç®�
    ! Calculation of near infrared heating rate in the case of Mars
    !
    use rad_Mars_NIR, only : RadMarsNIRInit

    ! Schneider and Liu (2009) ���¾å��¢ã����
    ! Radiation model by Schneider and Liu (2009)
    !
    use rad_SL09, only : RadSL09Init

    ! ç°¡å���¾å��¢ã����
    ! Simple radiation model
    !
    use rad_simple, only : RadSimpleInit

    ! ä½����������¾å��¢ã����
    ! radiation model with no absorption and no scattering
    !
    use rad_none, only : RadNoneInit

    ! �¾å��¢é�£ã���¼ã����
    ! Routines for radiation calculation
    !
    use rad_utils, only : RadUtilsInit

    ! ���´æ�¡æ�£ã����������
    ! Vertical diffusion flux
    !
    use vdiffusion_my, only : VDiffusionInit

    ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
    ! JMA turbulent mixing module
    !
    use vdiffusion_jma_my_wrapper, only : VDiffusionJMAInit

    ! Schneider and Liu (2009) ���������´æ··��課ç�
    ! Vertical diffusion by Schneider and Liu (2009)
    !
    use sl09_diffusion, only : SL09DiffusionInit

    ! Gravity wave drag by McFarlane (1987)
    ! Gravity wave drag by McFarlane (1987)
    !
    use gwd_m1987, only : GWDM1987Init

    ! �²ã�����¢ã����
    ! No cloud model
    !
    use cloud_none, only : CloudNoneInit

    ! ç°¡å���²ã�¢ã����
    ! Simple cloud
    !
    use cloud_simple, only: CloudSimpleInit

    ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
    ! Cloud model based on Tiedtke (1993)
    !
    use cloud_T1993base, only : CloudT1993baseInit

    ! ���� H2O �²ã�¢ã����
    ! Mars H2O cloud model
    !
    use cloud_mars_h2o, only : CloudMarsH2OInit

    ! �°è¡¨�¢ã���������� (�������)
    ! Surface flux (Bulk method)
    !
    use surface_flux_bulk, only : SurfaceFluxInit

    !
    ! set dust flux
    !
!!$    use set_dust_flux, only : SetDustFluxInit

    ! ����������������
    ! Lower boundary flux
    !
    use lb_flux_simple, only : LBFluxSimpleInit

    ! �°é�¢æ¸©åº�, ��å£�æ¸�º¦������ç©���
    ! Time integration of surface temperature and soil temperature
    !
    use intg_surftemp, only : IntgSurfTempInit

    ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only : DynamicsHSplVAS83Init

    ! ����������¿ã���ç®�����������å­����
    ! A dynamics for calculation with physical processes only
    !
    use dynamics_physicsonly, only : DynamicsPhysicsOnlyInit

    !
    ! Dynamical process for TWP-ICE experiment
    !
    use dynamics_twpice_scm_exp, only : DynamicsTWPICESCMExpInit

    ! 湿潤対æ�調ç�
    ! Moist convective adjustment
    !
    use moist_conv_adjust, only : MoistConvAdjustInit

    ! Relaxed Arakawa-Schubert scheme
    ! Relaxed Arakawa-Schubert scheme
    !
    use relaxed_arakawa_schubert, only : RASInit

    ! 大è�模å��çµ� (��対æ��§å��çµ�)
    ! Large scale condensation (non-convective condensation)
    !
    use lscond, only : LScaleCondInit

    ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
    ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
    !
!!$    use lscond_LL91, only : LScaleCondLL91Init

    !
    != Saturation adjustment
    !
    use saturation_adjust, only : SaturationAdjustInit

    ! ä¹¾ç�¥å�æµ�調ç�
    ! Dry convective adjustment
    !
    use dry_conv_adjust, only : DryConvAdjustInit

    ! ����������
    ! Gravitational sedimentation process
    !
    use grav_sed, only : GravSedInit

    ! 主æ�����¸å���
    ! Phase change of atmospheric major component
    !
    use major_comp_phase_change, only : MajorCompPhaseChangeInit

    ! �¿ã�¤ã�����£ã���¿ã�� (Asselin, 1972)
    ! Time filter (Asselin, 1972)
    !
    use timefilter_asselin1972, only : TimeFiltInit

    ! �������£ã���¿ã�� (Williams, 2009)
    ! Time filter (Williams, 2009)
    !
    use timefilter_williams2009, only: TimeFilterWilliams2009Init

    ! äº��±å��°ã���¤ã��¢ºèª�
    ! Check values of prognostic variables
    !
    use check_prog_vars, only : CheckProgVarsInit

    ! ������
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit

    ! Output frequently used variables
    ! Output frequently used variables
    !
    use output_freq_used_vars, only : OutputFreqUsedVarsInit

    ! 宣�� ; Declaration statements
    !
    implicit none

    character(*), parameter:: version = '$Name:  $' // '$Id: dcpam_main.f90,v 1.67 2015/03/11 04:54:54 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)      :: RadModel
                                ! Radiation model 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
    character(STRING)      :: MCMethod
                                ! Moist convection parameterization
    character(STRING)      :: LSCMethod
                                ! Large scale condensation parameterization
    character(STRING)      :: CloudMethod
                                ! Cloud model
    character(STRING)      :: SfcMoistMethod
                                ! Surface moist model
    character(STRING)      :: GWDMethod
                                ! Gravity wave drag
    character(STRING)      :: DCMethod
                                ! Dry convection

    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):: briefexplrad
                              ! å®�è¡����¡ã�¤ã����°¡æ½������ (�¾å����)
                              ! Brief account of executable file (radiation)
    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)

    logical:: FlagBucketModel

    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, RadModel, SfcFluxMethod, VDiffMethod, PhysImpMode, MCMethod, LSCMethod, CloudMethod, SfcMoistMethod, GWDMethod, DCMethod, 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, prog_name )

    ! NAMELIST ���¡ã�¤ã�����¥å��
    ! Input NAMELIST file name
    !
    call NmlutilInit( namelist_filename )

    ! �����������¤ã��¨­å®�
    ! Default values settings
    !
    DynMode                 = 'HSPLVAS83'
!!$    DynMode                 = 'NoHorAdv'

    PhysMode                = 'FullPhysics'
!!$    PhysMode                = 'HS94'
!!$    PhysMode                = 'VenusSimple'
!!$    PhysMode                = 'JupiterSimple'
!!$    PhysMode                = 'NoPhysics'

    RadModel                = 'DennouAGCM'
!!$    RadModel                = 'Earth'
!!$    RadModel                = 'Mars'
!!$    RadModel                = 'SL09'

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

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

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

!!$    MCMethod                = 'None'
    MCMethod                = 'MCA'
!!$    MCMethod                = 'RAS'

!!$    LSCMethod               = 'None'
    LSCMethod               = 'M65'
!!$    LSCMethod               = 'M65WithIce'
!!$    LSCMethod               = 'LL91'
!!$    LSCMethod               = 'LL91WithIce'

    CloudMethod             = 'None'
!!$    CloudMethod             = 'Simple'
!!$    CloudMethod             = 'T1993WithIce'

    SfcMoistMethod          = 'None'
!!$    SfcMoistMethod          = 'Bucket'

    GWDMethod               = 'None'
!!$    GWDMethod               = 'M1987'

!!$    DCMethod                = 'None'
    DCMethod                = 'DCA'


    FlagSnow                = .false.

    FlagMajCompPhaseChange  = .false.

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


    ! è¨�ç®��¢ã�¼ã����¨­å®�
    ! 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, 'DynMode=<%c>.', c1 = trim(DynMode) )
    !
    select case ( DynMode )
    case ( 'HSPLVAS83' )
      IDDynMode = IDDynModeHSPLVAS83
    case ( 'NoHorAdv' )
      IDDynMode = IDDynModeNoHorAdv
    case ( 'TWPICE' )
      IDDynMode = IDDynModeTWPICE
    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 ( 'HS94' )
      IDPhysMode = IDPhysModeHS94
    case ( 'VenusSimple' )
      IDPhysMode = IDPhysModeVenusSimple
    case ( 'JupiterSimple' )
      IDPhysMode = IDPhysModeJupiterSimple
    case ( 'JupiterSimpleV2' )
      IDPhysMode = IDPhysModeJupiterSimpleV2
    case ( 'NoPhysics' )
      IDPhysMode = IDPhysModeNoPhysics
    case default
      call MessageNotify( 'E', prog_name, 'PhysMode=<%c> is not supported.', c1 = trim(PhysMode) )
    end select


    ! Identification of calculation method for radiation
    !
    call MessageNotify( 'M', prog_name, 'RadModel=<%c>.', c1 = trim(RadModel) )
    !
    select case ( RadModel )
    case ( 'DennouAGCM' )
      IDRadMethod = IDRadMethodDennouAGCM
    case ( 'Earth' )
      IDRadMethod = IDRadMethodEarthV2
    case ( 'Mars' )
      IDRadMethod = IDRadMethodMarsV1
    case ( 'SL09' )
      IDRadMethod = IDRadMethodSL09
!!$    case ( 'VenusSimple' )
!!$      IDRadMethod = IDRadMethodVenusSimple
    case ( 'Simple' )
      IDRadMethod = IDRadMethodSimple
    case ( 'RRTMG' )
      IDRadMethod = IDRadMethodRRTMG
    case ( 'none' )
      IDRadMethod = IDRadMethodNone
    case default
      call MessageNotify( 'E', prog_name, 'RadModel=<%c> is not supported.', c1 = trim(RadModel) )
    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 ( 'JMA' )
      IDVDiffMethod = IDVDiffMethodJMA
    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 ( 'AtmOnly' )
      IDPhyTendMethod = IDPhyTendMethodImpAtmOnly
      FlagPhysImpSoilModelSO = .false.
    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
    if ( ( IDPhysMode /= IDPhysModeFullPhysics ) .and. ( IDPhyTendMethod == IDPhyTendMethodImpAtmOnly ) ) then
      call MessageNotify( 'E', prog_name, 'PhysMode has to be "FullPhysics" true, when PhyImpMode is "AtmOnly".' )
    end if


    ! Identification of calculation method for moist convection
    !
    call MessageNotify( 'M', prog_name, 'MCMethod=<%c>.', c1 = trim(MCMethod) )
    !
    select case ( MCMethod )
    case ( 'None' )
      IDMCMethod = IDMCMethodNone
    case ( 'MCA' )
      IDMCMethod = IDMCMethodMCA
    case ( 'RAS' )
      IDMCMethod = IDMCMethodRAS
    case ( 'RASWithIce' )
      IDMCMethod = IDMCMethodRASWithIce
    case default
      call MessageNotify( 'E', prog_name, 'MCMethod=<%c> is not supported.', c1 = trim(MCMethod) )
    end select


    ! Identification of calculation method for large scale condensation
    !
    call MessageNotify( 'M', prog_name, 'LSCMethod=<%c>.', c1 = trim(LSCMethod) )
    !
    select case ( LSCMethod )
    case ( 'None' )
      IDLSCMethod = IDLSCMethodNone
    case ( 'M65' )
      IDLSCMethod = IDLSCMethodM65
!!$    case ( 'LL91' )
!!$      IDLSCMethod = IDLSCMethodLL91
    case ( 'SatAdjM65' )
      IDLSCMethod = IDLSCMethodSatAdjM65
    case ( 'M65WithIce' )
      IDLSCMethod = IDLSCMethodM65WithIce
    case ( 'LL91WithIce' )
      IDLSCMethod = IDLSCMethodLL91WithIce
    case default
      call MessageNotify( 'E', prog_name, 'LSCMethod=<%c> is not supported.', c1 = trim(LSCMethod) )
    end select


    ! Identification of calculation method for cloud
    !
    call MessageNotify( 'M', prog_name, 'CloudMethod=<%c>.', c1 = trim(CloudMethod) )
    !
    select case ( CloudMethod )
    case ( 'None' )
      IDCloudMethod = IDCloudMethodNone
    case ( 'Simple' )
      IDCloudMethod = IDCloudMethodSimple
    case ( 'SimpleWithIce' )
      IDCloudMethod = IDCloudMethodSimpleWithIce
    case ( 'T1993WithIce' )
      IDCloudMethod = IDCloudMethodT1993WithIce
    case ( 'MarsH2OCloud' )
      IDCloudMethod = IDCloudMethodMarsH2OCloud
    case default
      call MessageNotify( 'E', prog_name, 'LSCloudMethod=<%c> is not supported.', c1 = trim(CloudMethod) )
    end select
    ! check a non-convective condensation
    if ( IDCloudMethod == IDCloudMethodT1993WithIce ) then
      if ( IDLSCMethod /= IDLSCMethodNone ) then
        call MessageNotify( 'E', prog_name, 'If CloudMethod="T1993WithIce", LscMethod has to be "None".' )
      end if
    end if

    ! Identification of calculation method for surface moisture
    !
    call MessageNotify( 'M', prog_name, 'SfcMoistMethod=<%c>.', c1 = trim(SfcMoistMethod) )
    !
    select case ( SfcMoistMethod )
    case ( 'None' )
      IDSfcMoistMethod = IDSfcMoistMethodNone
      FlagBucketModel  = .false.
    case ( 'Bucket' )
      IDSfcMoistMethod = IDSfcMoistMethodBucket
      FlagBucketModel  = .true.
    case default
      call MessageNotify( 'E', prog_name, 'SfcMoistMethod=<%c> is not supported.', c1 = trim(SfcMoistMethod) )
    end select


    ! Identification of calculation method for gravity wave drag
    !
    call MessageNotify( 'M', prog_name, 'GWDMethod=<%c>.', c1 = trim(GWDMethod) )
    !
    select case ( GWDMethod )
    case ( 'None' )
      IDGWDMethod = IDGWDMethodNone
    case ( 'M1987' )
      IDGWDMethod = IDGWDMethodM1987
    case default
      call MessageNotify( 'E', prog_name, 'GWDMethod=<%c> is not supported.', c1 = trim(GWDMethod) )
    end select


    ! Identification of calculation method for dry convection
    !
    call MessageNotify( 'M', prog_name, 'DCMethod=<%c>.', c1 = trim(DCMethod) )
    !
    select case ( DCMethod )
    case ( 'None' )
      IDDCMethod = IDDCMethodNone
    case ( 'DCA' )
      IDDCMethod = IDDCMethodDCA
    case default
      call MessageNotify( 'E', prog_name, 'DCMethod=<%c> is not supported.', c1 = trim(DCMethod) )
    end select



    ! è¨�ç®��¢ã�¼ã����¡¨ç¤�
    ! Display calculation mode
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      briefexpldyn = 'used'
    case ( IDDynModeNoHorAdv )
      briefexpldyn = 'not used'
    case ( IDDynModeTWPICE )
      briefexpldyn = 'TWP-ICE'
    end select

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      briefexplphy = 'parameterization suite is used'
    case ( IDPhysModeHS94 )
      briefexplphy = 'forcing for Held and Suarez (1994) dynamical core test'
    case ( IDPhysModeVenusSimple )
      briefexplphy = 'simple forcing for a Venus-like planet'
    case ( IDPhysModeJupiterSimple )
      briefexplphy = 'simple forcing for a Jupiter-like planet'
    case ( IDPhysModeJupiterSimpleV2 )
      briefexplphy = 'simple forcing for a Jupiter-like planet Ver. 2'
    case ( IDPhysModeNoPhysics )
      briefexplphy = 'not used'
    end select

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )
        briefexplrad = 'dennou AGCM5 default'
      case ( IDRadMethodEarthV2 )
        briefexplrad = 'Earth'
      case ( IDRadMethodMarsV1 )
        briefexplrad = 'Mars'
      case ( IDRadMethodSL09 )
        briefexplrad = 'Schneider and Liu (2009) (Jupiter-like planet)'
!!$      case ( IDRadMethodVenusSimple )
!!$        briefexplrad = 'Venus simple radiation model (tentative)'
      case ( IDRadMethodSimple )
        briefexplrad = 'Simple radiation model (tentative)'
      case ( IDRadMethodRRTMG )
        briefexplrad = 'RRTMG'
      case ( IDRadMethodNone )
        briefexplrad = 'none'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplrad', '' )
      end select

      select case ( IDSfcFluxMethod )
      case ( IDSfcFluxMethodL82 )
        briefexplsfcflux = 'surface flux by the method of Louis et al. (1982)'
      case ( IDSfcFluxMethodBH91B94 )
        briefexplsfcflux = 'surface flux by the 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 = 'vertical diffusion with the method of Mellor and Yamada level 2'
      case ( IDVDiffMethodMY25 )
        briefexplvdiff = 'vertical diffusion with the method of Mellor and Yamada level 2.5'
      case ( IDVDiffMethodJMA )
        briefexplvdiff = 'vertical diffusion with the method of JMA physical library of 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 ( 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 ( IDPhysMode == IDPhysModeFullPhysics ) then
      call MessageNotify( 'M', prog_name, '|    Radiation model : %c', c1 = trim(briefexplrad) )
      call MessageNotify( 'M', prog_name, '|    Surface flux method : %c', c1 = trim(briefexplsfcflux) )
      call MessageNotify( 'M', prog_name, '|    Vertical diffusion method : %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, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, 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_SurfMajCompIceB(0:imax-1, 1:jmax) )
      allocate( xy_SoilMoistB     (0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowB      (0:imax-1, 1:jmax) )

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

      allocate( xy_SurfMajCompIceA(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)
      call RestartSurfTempInit
      call RestartSurfTempGet( xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN )

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

!!$      xy_SurfMajCompIceN = 0.0_DP
!!$      xy_SurfMajCompIceB = xy_SurfMajCompIceN

    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 )

    ! Output frequently used variables
    ! Output frequently used variables
    !
    call OutputFreqUsedVarsInit


    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimpleV2 )

      ! ���¹ã�������¼ã�¿å�ºå�����������¸ã����°ç�»é��
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'Rain', (/ 'lon ', 'lat ', 'time' /), 'precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'Snow', (/ 'lon ', 'lat ', 'time' /), 'precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCP', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation', 'kg m-2 s-1' )

    case ( IDPhysModeFullPhysics )

      ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿ã���¡ã�¤ã����������
      ! 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( 'SurfMajCompIce' , (/ 'lon ', 'lat ', 'time' /), 'surface major component ice', 'kg m-2' )     ! (in)

      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 (rain)', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'Snow', (/ 'lon ', 'lat ', 'time' /), 'precipitation (snow)', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'PRCP', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'RainCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation', 'kg m-2 s-1' )

      ! ���¹ã�������¼ã�¿å�ºå�� (�¹ã�¿ã�¼ã������)
      ! History data output (Start time)
      !
      call HistoryAutoPut( TimeN, 'SurfMajCompIce', xy_SurfMajCompIceN )
      call HistoryAutoPut( TimeN, 'SoilMoist'     , xy_SoilMoistN      )
      call HistoryAutoPut( TimeN, 'SurfSnow'      , xy_SurfSnowN       )
    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( xyz_OMG(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) )
    allocate( xyz_Exner    (0:imax-1, 1:jmax, 1:kmax) )

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      allocate( xy_SurfAlbedo         (0:imax-1, 1:jmax) )
      allocate( xy_SurfHumidCoef      (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLenMom    (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLenHeat   (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_SurfHeightStd      (0:imax-1, 1:jmax) )
      allocate( xy_SnowFrac           (0:imax-1, 1:jmax) )

!!$      allocate( xy_FlagMatthewsLand         (0:imax-1, 1:jmax) )
      allocate( xy_PhyImplSDHIndexCalcMethod(0:imax-1, 1:jmax) )
      allocate( xy_BucketFlagOceanGrid      (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( 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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixDiffCoef(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( xy_SurfMOLength     (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_DUDtGWD(0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DVDtGWD(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( xy_RainCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_RainLsc      (0:imax-1, 1:jmax) )
      allocate( xy_Snow         (0:imax-1, 1:jmax) )
      allocate( xy_SnowCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_SnowLsc      (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQVapDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OSolDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DUDtCum       (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DVDtCum       (0:imax-1,1:jmax,1:kmax) )

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

      allocate( xyz_MoistConvDetTend       (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_MoistConvSubsidMassFlux(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) )
      allocate( xyz_CloudCoverforRad(0:imax-1,1:jmax,1:kmax) )

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

      ! ���¹ã�������¼ã�¿å�ºå�����������¸ã����°ç�»é��
      ! 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' )

      call HistoryAutoAddVariable( 'SurfRoughLenMom' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length for momentum', 'm' )

      call HistoryAutoAddVariable( 'SurfRoughLenHeat' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length for heat', 'm' )

      call HistoryAutoAddVariable( 'SurfDustGravSedFlux' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length', 'm' )

    case ( IDPhysModeVenusSimple )

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

    case ( IDPhysModeJupiterSimple )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_VirTemp(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_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( 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) )
      allocate( xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax,  0:1) )
      allocate( xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax,  0:1) )


    case ( IDPhysModeJupiterSimpleV2 )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_VirTemp(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( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixDiffCoef(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_RadLFlux      (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( 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) )
      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( xy_Rain         (0:imax-1, 1:jmax) )
      allocate( xy_RainCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_RainLsc      (0:imax-1, 1:jmax) )
      allocate( xy_Snow         (0:imax-1, 1:jmax) )
      allocate( xy_SnowCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_SnowLsc      (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQVapDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OSolDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DUDtCum       (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DVDtCum       (0:imax-1,1:jmax,1:kmax) )

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

    end select







    ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
    ! Setting planetary surface properties
    !
    select case ( IDSfcMoistMethod )
    case ( IDSfcMoistMethodNone )
      call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .false., FlagSnow )
    case ( IDSfcMoistMethodBucket )
      call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .true. , FlagSnow )
    end select


    select case ( IDPhysMode )
    case ( IDPhysModeNoPhysics )

      ! Nothing to do

    case ( IDPhysModeHS94 )

      ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
      ! Forcing and dissipation suggested by Held and Suarez (1994)
      !
      call HS94Init

    case ( IDPhysModeVenusSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! Yamamoto and Takahashi (2003) ����£ã��ç°¡å������è¨�ç®���������¼·��
      ! forcing for simple Venus calculation following Yamamoto and Takahashi (2003)
      !
      call YT2003ForcingInit

    case ( IDPhysModeJupiterSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! Schneider and Liu (2009) ���������´æ··��課ç�
      ! Vertical diffusion by Schneider and Liu (2009)
      !
      call SL09DiffusionInit

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Init

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    case ( IDPhysModeJupiterSimpleV2 )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! ����������������
      ! Lower boundary flux
      !
      call LBFluxSimpleInit

      ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
      ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
      !
      call PhyImplAtmOnlyInit

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Init

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    case ( IDPhysModeFullPhysics )

      ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
      ! Setting planetary surface properties
      !
      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodNone )
        call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .false., FlagSnow )
      case ( IDSfcMoistMethodBucket )
        call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .true. , FlagSnow )
      end select

      ! ��, æ°·ã���²å��
      ! snow/ice fraction
      !
      call SnowIceFracInit

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )

        ! �¾å����������� (GFD �»è�³å�¶æ¥½�����ºã���¾å��¢ã����)
        ! Radiation flux (radiation model developed by GFD Dennou Club)
        !
        call RadDennouAGCMInit( flag_rst = .not. flag_initial )

      case ( IDRadMethodEarthV2 )

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

        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          ! �²ã�����¢ã����
          ! No cloud model
          !
          call CloudNoneInit( FlagSnow )
        case ( IDCloudMethodSimple )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodSimpleWithIce )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodT1993WithIce )
          ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
          ! Cloud model based on Tiedtke (1993)
          !
          call CloudT1993baseInit( FlagSnow )
        end select

        ! �°ç��大æ������¾å��¢ã���� Ver. 2
        ! radiation model for the Earth's atmosphere Ver. 2
        !
        call RadEarthV2Init( FlagSnow )

      case ( IDRadMethodRRTMG )

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

        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          ! �²ã�����¢ã����
          ! No cloud model
          !
          call CloudNoneInit( FlagSnow )
        case ( IDCloudMethodSimple )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodSimpleWithIce )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodT1993WithIce )
          ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
          ! Cloud model based on Tiedtke (1993)
          !
          call CloudT1993baseInit( FlagSnow )
        end select

        ! wrapper of RRTMG
        ! wrapper of RRTMG
        !
        call RadRRTMGWrapperInit( FlagSnow )

      case ( IDRadMethodMarsV1 )

        ! ����大æ������¾å��¢ã���� Ver. 1
        ! radiation model for the Mars' atmosphere Ver. 1
        !
        call RadMarsV1Init

      case ( IDRadMethodSL09 )

        ! Schneider and Liu (2009) ���¾å��¢ã����
        ! Radiation model by Schneider and Liu (2009)
        !
        call RadSL09Init

      case ( IDRadMethodSimple )

        ! ç°¡å���¾å��¢ã����
        ! Simple radiation model
        !
        call RadSimpleInit

      case ( IDRadMethodNone )

        ! ä½����������¾å��¢ã����
        ! radiation model with no absorption and no scattering
        !
        call RadNoneInit

      end select


      ! ���´æ�¡æ�£ã����������
      ! Vertical diffusion flux
      !
      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY2, IDVDiffMethodMY25 )

        ! ���´æ�¡æ�£ã���������� (Mellor and Yamada, 1974)
        ! Vertical diffusion flux (Mellor and Yamada, 1974)
        !
        call VDiffusionInit

      case ( IDVDiffMethodJMA )
        ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
        ! JMA turbulent mixing module
        !
        call VDiffusionJMAInit

      end select

      ! �°è¡¨�¢ã���������� (�������)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

      !
      ! set dust flux
      !
!!$      call SetDustFluxInit

      ! ä¸������������������å¤��������ç®� (�°è§£æ³�)
      ! Calculate tendency by a part of physical processes (implicit)
      !
      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )

        ! �°è§£æ³�����������ç©���
        ! Time integration with implicit scheme
        !
        call PhyImplInit( FlagBucketModel, FlagSnow )

      case ( IDPhyTendMethodImpSoilModel )

        select case ( IDSfcMoistMethod )
        case ( IDSfcMoistMethodBucket )
          ! ���±ã���¢ã����
          ! Bucket model
          !
          call BucketModelInit( FlagSnow )

        end select

        ! �°ä����������±ã�����´æ�¡æ��
        ! Vertical diffusion of heat under the ground
        !
        call SubsurfaceDiffusionInit

        ! �°ä��±ä�å°��¢ã�������������´å�����°è§£æ³�����������ç©���
        !
        ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
!!$        call PhyImplSDHInit(                                                 &
!!$          & FlagBucketModel, FlagSnow,                                       &
!!$          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
!!$          & )

!!$        call PhyImplSDHV2Init(                                                &
!!$          & FlagBucketModel, FlagSnow,                                       &
!!$          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
!!$          & )
        call PhyImplSDHV3Init( FlagBucketModel, FlagSnow, FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName )

      case ( IDPhyTendMethodImpAtmOnly )

        ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
        ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
        !
        call PhyImplAtmOnlyInit

      end select


      ! Gravity wave drag by McFarlane (1987)
      ! Gravity wave drag by McFarlane (1987)
      !
      call GWDM1987Init

      ! �°è§£æ³�����������ç©��������������¼ã����
      ! Routines for time integration with implicit scheme
      !
      call PhyImplUtilsInit

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit


!!$      select case ( IDRadMethod )
!!$      case ( IDRadMethodMarsV1 )
!!$        ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$        ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$        !
!!$        call Rad15mNLTEInit
!!$      end select

      select case ( IDRadMethod )
      case ( IDRadMethodMarsV1 )
        ! ����è¨�ç®����赤å����±ç�è¨�ç®�
        ! Calculation of near infrared heating rate in the case of Mars
        !
        call RadMarsNIRInit
      end select


      ! ���´æ�¡æ�£ã���������� (Mellor and Yamada, 1974)
      ! Vertical diffusion flux (Mellor and Yamada, 1974)
      !
      call VDiffusionInit

      ! �°è¡¨�¢ã���������� (�������)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    end select



    ! �°é�¢æ¸©åº��»å��å£�æ¸�º¦�»å��å£�æ°´å���»ç���������
    ! Time integration of surface temperature, soil temperature, soil moisture, 
    ! and surface snow amount
    !
    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      ! �°é�¢æ¸©åº�, ��å£�æ¸�º¦������ç©���
      ! Time integration of surface temperature and soil temperature
      !
      call IntgSurfTempInit

      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodBucket )
        ! ���±ã���¢ã����
        ! Bucket model
        !
        call BucketModelInit( FlagSnow )
      end select

    end select




    ! ������
    ! Dynamical core
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
      ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
      !
      call DynamicsHSplVAS83Init
    case ( IDDynModeNoHorAdv )
      ! ����������¿ã���ç®�����������å­����
      ! A dynamics for calculation with physical processes only
      !
      call DynamicsPhysicsOnlyInit
    case ( IDDynModeTWPICE )
      !
      ! Dynamical process for TWP-ICE experiment
      !
      call DynamicsTWPICESCMExpInit
    end select


    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

    case ( IDPhysModeJupiterSimpleV2 )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! 湿潤対æ�調ç�
      ! Moist convective adjustment
      !
      call MoistConvAdjustInit

      ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
      ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
      !
      call LScaleCondInit( FlagSnow )

      ! ç°¡å���²ã�¢ã����
      ! Simple cloud
      !
      call CloudSimpleInit( FlagSnow )

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

      ! ������
      ! Mass fixer
      !
      call MassFixerInit

    case ( IDPhysModeFullPhysics )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDMCMethod )
      case ( IDMCMethodMCA )
        ! 湿潤対æ�調ç�
        ! Moist convective adjustment
        !
        call MoistConvAdjustInit
      case ( IDMCMethodRAS )
        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        call RASInit
      case ( IDMCMethodRASWithIce )
        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        call RASInit
      end select

      select case ( IDLSCMethod )
      case ( IDLSCMethodM65 )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
!!$      case ( IDLSCMethodLL91 )
!!$        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
!!$        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
!!$        !
!!$        call LScaleCondLL91Init
      case ( IDLSCMethodSatAdjM65 )
        !
        != Saturation adjustment
        !
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        call SaturationAdjustInit( FlagSnow )
      case ( IDLSCMethodM65WithIce )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
      case ( IDLSCMethodLL91WithIce )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
      end select


      ! 
      ! Cloud model
      !
      select case ( IDCloudMethod )
      case ( IDCloudMethodNone )
        ! �²ã�����¢ã����
        ! No cloud model
        !
        call CloudNoneInit( FlagSnow )
      case ( IDCloudMethodSimple )
        ! ç°¡å���²ã�¢ã����
        ! Simple cloud
        !
        call CloudSimpleInit( FlagSnow )
      case ( IDCloudMethodSimpleWithIce )
        ! ç°¡å���²ã�¢ã����
        ! Simple cloud
        !
        call CloudSimpleInit( FlagSnow )
      case ( IDCloudMethodT1993WithIce )
        ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
        ! Cloud model based on Tiedtke (1993)
        !
        call CloudT1993baseInit( FlagSnow )
      case ( IDCloudMethodMarsH2OCloud )
        ! ���� H2O �²ã�¢ã����
        ! Mars H2O cloud model
        !
        call CloudMarsH2OInit
      end select


      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodBucket )
        ! ���±ã���¢ã����
        ! Bucket model
        !
        call BucketModelInit( FlagSnow )
      end select

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

      ! ����������
      ! Gravitational sedimentation process
      !
      call GravSedInit

      ! 主æ�����¸å���
      ! Phase change of atmospheric major component
      !
      call MajorCompPhaseChangeInit( FlagMajCompPhaseChange, CondMajCompName )

      ! ������
      ! Mass fixer
      !
      call MassFixerInit

    end select


    ! �¿ã�¤ã�����£ã���¿ã�� (Asselin, 1972)
    ! Time filter (Asselin, 1972)
    !
!!$    call TimeFiltInit

    ! �������£ã���¿ã�� (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

    ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only : DynamicsHSplVAS83Finalize

    !
    ! Dynamical process for TWP-ICE experiment
    !
    use dynamics_twpice_scm_exp, only : DynamicsTWPICESCMExpFinalize

    ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only: Hs94Finalize

    ! �¾å����������� (���³ã���¢ã����)
    ! Radiation flux (band model)
    !
    use rad_DennouAGCM, only: RadDennouAGCMFinalize

    ! 座æ����¼ã�¿è¨­å®�
    ! 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

    ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! 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 ( IDPhysMode == IDPhysModeFullPhysics ) 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( xyz_DTurKinEneDt )

    deallocate( xyz_OMG )

    deallocate( xy_SurfHeight )
    deallocate( xyz_Height )

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      deallocate( xy_SurfAlbedo          )
      deallocate( xy_SurfHumidCoef       )
      deallocate( xy_SurfRoughLenMom     )
      deallocate( xy_SurfRoughLenHeat    )
      deallocate( xy_SurfHeatCapacity    )
      deallocate( xy_SeaIceConc          )
      deallocate( xy_SurfCond            )
      deallocate( xy_SurfType            )
      deallocate( xy_DeepSubSurfHeatFlux )
      deallocate( xy_SoilHeatCap         )
      deallocate( xy_SoilHeatDiffCoef    )
      deallocate( xy_SurfHeightStd       )
      deallocate( xy_SnowFrac            )

!!$      deallocate( xy_FlagMatthewsLand          )
      deallocate( xy_PhyImplSDHIndexCalcMethod )
      deallocate( xy_BucketFlagOceanGrid       )

      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_VelDiffCoef  )
      deallocate( xyr_TempDiffCoef )
      deallocate( xyr_QMixDiffCoef )

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

      deallocate( xy_SurfMOLength )

      deallocate( xy_DSurfTempDt )
      deallocate( xyz_DSoilTempDt )

      deallocate( xy_DPsDt )

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

      deallocate( xyz_DTempDtVDiff )

      deallocate( xyz_DUDtGWD )
      deallocate( xyz_DVDtGWD )

      deallocate( xyz_DTempDtRadL )
      deallocate( xyz_DTempDtRadS )


      deallocate( xy_Rain          )
      deallocate( xy_RainCumulus   )
      deallocate( xy_RainLsc       )
      deallocate( xy_SnowCumulus   )
      deallocate( xy_SnowLsc       )

      deallocate( xyz_DTempDtCum     )
      deallocate( xyz_DQVapDtCum     )
      deallocate( xyz_DQH2OLiqDtCum  )
      deallocate( xyz_DQH2OSolDtCum  )
      deallocate( xyz_DUDtCum        )
      deallocate( xyz_DVDtCum        )

      deallocate( xyz_DQH2OLiqDtLSC  )
      deallocate( xyz_DQH2OSolDtLSC  )

      deallocate( xyz_MoistConvDetTend        )
      deallocate( xyz_MoistConvSubsidMassFlux )

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



    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      ! ���¢ã�¸ã�¥ã�¼ã��������°ã���²ä�解é��
      ! Dellocation of variables in modules
      !
      call DynamicsHSplVAS83Finalize
    case ( IDDynModeNoHorAdv )
    case ( IDDynModeTWPICE )
      !
      ! Dynamical process for TWP-ICE experiment
      !
      call DynamicsTWPICESCMExpFinalize
    end select


    call AuxVarsFinalize

    select case ( IDPhysMode )
    case ( IDPhysModeHS94 )
      call Hs94Finalize
    end select

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      ! �²ä�解é�¤ã�����¹ã�¿ã�¼ã�����¡ã�¤ã�����äº�����
      ! Dellocation and close a restart file
      !
      call RadDennouAGCMFinalize
    end select

    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

    ! ���»ç���
    ! 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: RestartFileInit, RestartFileOpen, RestartFileGet

    ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! Restart data of surface temperature input/output
    !
    use restart_surftemp_io, only: RestartSurfTempInit, 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

    ! ��ç³»å�����¼ã�¿ã����¿è¾¼��
    ! 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

    ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only : HS94Init

    ! Yamamoto and Takahashi (2003) ����£ã��ç°¡å������è¨�ç®���������¼·��
    ! forcing for simple Venus calculation following Yamamoto and Takahashi (2003)
    !
    use yt2003_forcing, only : YT2003ForcingInit

    ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
    ! Setting planetary surface properties
    !
    use surface_properties, only : SurfacePropertiesInit

    ! ��, æ°·ã���²å��
    ! snow/ice fraction
    !
    use snowice_frac, only : SnowIceFracInit

    ! �°ä����������±ã�����´æ�¡æ��
    ! 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
    ! �°ä��±ä�å°��¢ã�������������´å�����°è§£æ³�����������ç©���
    !
    ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
    use phy_implicit_sdh_V3, only : PhyImplSDHV3Init

    ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
    ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
    !
    use phy_implicit_atmonly, only : PhyImplAtmOnlyInit

    ! �°è§£æ³�����������ç©��������������¼ã����
    ! Routines for time integration with implicit scheme
    !
    use phy_implicit_utils, only : PhyImplUtilsInit

    ! ���±ã���¢ã����
    ! Bucket model
    !
    use Bucket_Model, only : BucketModelInit

    ! �¾å����������� (GFD �»è�³å�¶æ¥½�����ºã���¾å��¢ã����)
    ! Radiation flux (radiation model developed by GFD Dennou Club)
    !
    use rad_DennouAGCM, only : RadDennouAGCMInit

    ! �°ç��大æ������¾å��¢ã���� Ver. 2
    ! radiation model for the Earth's atmosphere Ver. 2
    !
    use rad_Earth_V2, only : RadEarthV2Init

    ! wrapper of RRTMG
    ! wrapper of RRTMG
    !
    use rad_rrtmg_wrapper, only : RadRRTMGWrapperInit

    ! ����大æ������¾å��¢ã���� Ver. 1
    ! radiation model for the Mars' atmosphere Ver. 1
    !
    use rad_Mars_V1, only : RadMarsV1Init

!!$    ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$    ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$    !
!!$    use rad_15m_NLTE, only: Rad15mNLTEInit

    ! ����è¨�ç®����赤å����±ç�è¨�ç®�
    ! Calculation of near infrared heating rate in the case of Mars
    !
    use rad_Mars_NIR, only : RadMarsNIRInit

    ! Schneider and Liu (2009) ���¾å��¢ã����
    ! Radiation model by Schneider and Liu (2009)
    !
    use rad_SL09, only : RadSL09Init

    ! ç°¡å���¾å��¢ã����
    ! Simple radiation model
    !
    use rad_simple, only : RadSimpleInit

    ! ä½����������¾å��¢ã����
    ! radiation model with no absorption and no scattering
    !
    use rad_none, only : RadNoneInit

    ! �¾å��¢é�£ã���¼ã����
    ! Routines for radiation calculation
    !
    use rad_utils, only : RadUtilsInit

    ! ���´æ�¡æ�£ã����������
    ! Vertical diffusion flux
    !
    use vdiffusion_my, only : VDiffusionInit

    ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
    ! JMA turbulent mixing module
    !
    use vdiffusion_jma_my_wrapper, only : VDiffusionJMAInit

    ! Schneider and Liu (2009) ���������´æ··��課ç�
    ! Vertical diffusion by Schneider and Liu (2009)
    !
    use sl09_diffusion, only : SL09DiffusionInit

    ! Gravity wave drag by McFarlane (1987)
    ! Gravity wave drag by McFarlane (1987)
    !
    use gwd_m1987, only : GWDM1987Init

    ! �²ã�����¢ã����
    ! No cloud model
    !
    use cloud_none, only : CloudNoneInit

    ! ç°¡å���²ã�¢ã����
    ! Simple cloud
    !
    use cloud_simple, only: CloudSimpleInit

    ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
    ! Cloud model based on Tiedtke (1993)
    !
    use cloud_T1993base, only : CloudT1993baseInit

    ! ���� H2O �²ã�¢ã����
    ! Mars H2O cloud model
    !
    use cloud_mars_h2o, only : CloudMarsH2OInit

    ! �°è¡¨�¢ã���������� (�������)
    ! Surface flux (Bulk method)
    !
    use surface_flux_bulk, only : SurfaceFluxInit

    !
    ! set dust flux
    !
!!$    use set_dust_flux, only : SetDustFluxInit

    ! ����������������
    ! Lower boundary flux
    !
    use lb_flux_simple, only : LBFluxSimpleInit

    ! �°é�¢æ¸©åº�, ��å£�æ¸�º¦������ç©���
    ! Time integration of surface temperature and soil temperature
    !
    use intg_surftemp, only : IntgSurfTempInit

    ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only : DynamicsHSplVAS83Init

    ! ����������¿ã���ç®�����������å­����
    ! A dynamics for calculation with physical processes only
    !
    use dynamics_physicsonly, only : DynamicsPhysicsOnlyInit

    !
    ! Dynamical process for TWP-ICE experiment
    !
    use dynamics_twpice_scm_exp, only : DynamicsTWPICESCMExpInit

    ! 湿潤対æ�調ç�
    ! Moist convective adjustment
    !
    use moist_conv_adjust, only : MoistConvAdjustInit

    ! Relaxed Arakawa-Schubert scheme
    ! Relaxed Arakawa-Schubert scheme
    !
    use relaxed_arakawa_schubert, only : RASInit

    ! 大è�模å��çµ� (��対æ��§å��çµ�)
    ! Large scale condensation (non-convective condensation)
    !
    use lscond, only : LScaleCondInit

    ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
    ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
    !
!!$    use lscond_LL91, only : LScaleCondLL91Init

    !
    != Saturation adjustment
    !
    use saturation_adjust, only : SaturationAdjustInit

    ! ä¹¾ç�¥å�æµ�調ç�
    ! Dry convective adjustment
    !
    use dry_conv_adjust, only : DryConvAdjustInit

    ! ����������
    ! Gravitational sedimentation process
    !
    use grav_sed, only : GravSedInit

    ! 主æ�����¸å���
    ! Phase change of atmospheric major component
    !
    use major_comp_phase_change, only : MajorCompPhaseChangeInit

    ! �¿ã�¤ã�����£ã���¿ã�� (Asselin, 1972)
    ! Time filter (Asselin, 1972)
    !
    use timefilter_asselin1972, only : TimeFiltInit

    ! �������£ã���¿ã�� (Williams, 2009)
    ! Time filter (Williams, 2009)
    !
    use timefilter_williams2009, only: TimeFilterWilliams2009Init

    ! äº��±å��°ã���¤ã��¢ºèª�
    ! Check values of prognostic variables
    !
    use check_prog_vars, only : CheckProgVarsInit

    ! ������
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit

    ! Output frequently used variables
    ! Output frequently used variables
    !
    use output_freq_used_vars, only : OutputFreqUsedVarsInit

    ! 宣�� ; Declaration statements
    !
    implicit none

    character(*), parameter:: version = '$Name:  $' // '$Id: dcpam_main.f90,v 1.67 2015/03/11 04:54:54 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)      :: RadModel
                                ! Radiation model 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
    character(STRING)      :: MCMethod
                                ! Moist convection parameterization
    character(STRING)      :: LSCMethod
                                ! Large scale condensation parameterization
    character(STRING)      :: CloudMethod
                                ! Cloud model
    character(STRING)      :: SfcMoistMethod
                                ! Surface moist model
    character(STRING)      :: GWDMethod
                                ! Gravity wave drag
    character(STRING)      :: DCMethod
                                ! Dry convection

    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):: briefexplrad
                              ! å®�è¡����¡ã�¤ã����°¡æ½������ (�¾å����)
                              ! Brief account of executable file (radiation)
    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)

    logical:: FlagBucketModel

    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, RadModel, SfcFluxMethod, VDiffMethod, PhysImpMode, MCMethod, LSCMethod, CloudMethod, SfcMoistMethod, GWDMethod, DCMethod, 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, prog_name )

    ! NAMELIST ���¡ã�¤ã�����¥å��
    ! Input NAMELIST file name
    !
    call NmlutilInit( namelist_filename )

    ! �����������¤ã��¨­å®�
    ! Default values settings
    !
    DynMode                 = 'HSPLVAS83'
!!$    DynMode                 = 'NoHorAdv'

    PhysMode                = 'FullPhysics'
!!$    PhysMode                = 'HS94'
!!$    PhysMode                = 'VenusSimple'
!!$    PhysMode                = 'JupiterSimple'
!!$    PhysMode                = 'NoPhysics'

    RadModel                = 'DennouAGCM'
!!$    RadModel                = 'Earth'
!!$    RadModel                = 'Mars'
!!$    RadModel                = 'SL09'

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

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

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

!!$    MCMethod                = 'None'
    MCMethod                = 'MCA'
!!$    MCMethod                = 'RAS'

!!$    LSCMethod               = 'None'
    LSCMethod               = 'M65'
!!$    LSCMethod               = 'M65WithIce'
!!$    LSCMethod               = 'LL91'
!!$    LSCMethod               = 'LL91WithIce'

    CloudMethod             = 'None'
!!$    CloudMethod             = 'Simple'
!!$    CloudMethod             = 'T1993WithIce'

    SfcMoistMethod          = 'None'
!!$    SfcMoistMethod          = 'Bucket'

    GWDMethod               = 'None'
!!$    GWDMethod               = 'M1987'

!!$    DCMethod                = 'None'
    DCMethod                = 'DCA'


    FlagSnow                = .false.

    FlagMajCompPhaseChange  = .false.

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


    ! è¨�ç®��¢ã�¼ã����¨­å®�
    ! 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, 'DynMode=<%c>.', c1 = trim(DynMode) )
    !
    select case ( DynMode )
    case ( 'HSPLVAS83' )
      IDDynMode = IDDynModeHSPLVAS83
    case ( 'NoHorAdv' )
      IDDynMode = IDDynModeNoHorAdv
    case ( 'TWPICE' )
      IDDynMode = IDDynModeTWPICE
    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 ( 'HS94' )
      IDPhysMode = IDPhysModeHS94
    case ( 'VenusSimple' )
      IDPhysMode = IDPhysModeVenusSimple
    case ( 'JupiterSimple' )
      IDPhysMode = IDPhysModeJupiterSimple
    case ( 'JupiterSimpleV2' )
      IDPhysMode = IDPhysModeJupiterSimpleV2
    case ( 'NoPhysics' )
      IDPhysMode = IDPhysModeNoPhysics
    case default
      call MessageNotify( 'E', prog_name, 'PhysMode=<%c> is not supported.', c1 = trim(PhysMode) )
    end select


    ! Identification of calculation method for radiation
    !
    call MessageNotify( 'M', prog_name, 'RadModel=<%c>.', c1 = trim(RadModel) )
    !
    select case ( RadModel )
    case ( 'DennouAGCM' )
      IDRadMethod = IDRadMethodDennouAGCM
    case ( 'Earth' )
      IDRadMethod = IDRadMethodEarthV2
    case ( 'Mars' )
      IDRadMethod = IDRadMethodMarsV1
    case ( 'SL09' )
      IDRadMethod = IDRadMethodSL09
!!$    case ( 'VenusSimple' )
!!$      IDRadMethod = IDRadMethodVenusSimple
    case ( 'Simple' )
      IDRadMethod = IDRadMethodSimple
    case ( 'RRTMG' )
      IDRadMethod = IDRadMethodRRTMG
    case ( 'none' )
      IDRadMethod = IDRadMethodNone
    case default
      call MessageNotify( 'E', prog_name, 'RadModel=<%c> is not supported.', c1 = trim(RadModel) )
    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 ( 'JMA' )
      IDVDiffMethod = IDVDiffMethodJMA
    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 ( 'AtmOnly' )
      IDPhyTendMethod = IDPhyTendMethodImpAtmOnly
      FlagPhysImpSoilModelSO = .false.
    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
    if ( ( IDPhysMode /= IDPhysModeFullPhysics ) .and. ( IDPhyTendMethod == IDPhyTendMethodImpAtmOnly ) ) then
      call MessageNotify( 'E', prog_name, 'PhysMode has to be "FullPhysics" true, when PhyImpMode is "AtmOnly".' )
    end if


    ! Identification of calculation method for moist convection
    !
    call MessageNotify( 'M', prog_name, 'MCMethod=<%c>.', c1 = trim(MCMethod) )
    !
    select case ( MCMethod )
    case ( 'None' )
      IDMCMethod = IDMCMethodNone
    case ( 'MCA' )
      IDMCMethod = IDMCMethodMCA
    case ( 'RAS' )
      IDMCMethod = IDMCMethodRAS
    case ( 'RASWithIce' )
      IDMCMethod = IDMCMethodRASWithIce
    case default
      call MessageNotify( 'E', prog_name, 'MCMethod=<%c> is not supported.', c1 = trim(MCMethod) )
    end select


    ! Identification of calculation method for large scale condensation
    !
    call MessageNotify( 'M', prog_name, 'LSCMethod=<%c>.', c1 = trim(LSCMethod) )
    !
    select case ( LSCMethod )
    case ( 'None' )
      IDLSCMethod = IDLSCMethodNone
    case ( 'M65' )
      IDLSCMethod = IDLSCMethodM65
!!$    case ( 'LL91' )
!!$      IDLSCMethod = IDLSCMethodLL91
    case ( 'SatAdjM65' )
      IDLSCMethod = IDLSCMethodSatAdjM65
    case ( 'M65WithIce' )
      IDLSCMethod = IDLSCMethodM65WithIce
    case ( 'LL91WithIce' )
      IDLSCMethod = IDLSCMethodLL91WithIce
    case default
      call MessageNotify( 'E', prog_name, 'LSCMethod=<%c> is not supported.', c1 = trim(LSCMethod) )
    end select


    ! Identification of calculation method for cloud
    !
    call MessageNotify( 'M', prog_name, 'CloudMethod=<%c>.', c1 = trim(CloudMethod) )
    !
    select case ( CloudMethod )
    case ( 'None' )
      IDCloudMethod = IDCloudMethodNone
    case ( 'Simple' )
      IDCloudMethod = IDCloudMethodSimple
    case ( 'SimpleWithIce' )
      IDCloudMethod = IDCloudMethodSimpleWithIce
    case ( 'T1993WithIce' )
      IDCloudMethod = IDCloudMethodT1993WithIce
    case ( 'MarsH2OCloud' )
      IDCloudMethod = IDCloudMethodMarsH2OCloud
    case default
      call MessageNotify( 'E', prog_name, 'LSCloudMethod=<%c> is not supported.', c1 = trim(CloudMethod) )
    end select
    ! check a non-convective condensation
    if ( IDCloudMethod == IDCloudMethodT1993WithIce ) then
      if ( IDLSCMethod /= IDLSCMethodNone ) then
        call MessageNotify( 'E', prog_name, 'If CloudMethod="T1993WithIce", LscMethod has to be "None".' )
      end if
    end if

    ! Identification of calculation method for surface moisture
    !
    call MessageNotify( 'M', prog_name, 'SfcMoistMethod=<%c>.', c1 = trim(SfcMoistMethod) )
    !
    select case ( SfcMoistMethod )
    case ( 'None' )
      IDSfcMoistMethod = IDSfcMoistMethodNone
      FlagBucketModel  = .false.
    case ( 'Bucket' )
      IDSfcMoistMethod = IDSfcMoistMethodBucket
      FlagBucketModel  = .true.
    case default
      call MessageNotify( 'E', prog_name, 'SfcMoistMethod=<%c> is not supported.', c1 = trim(SfcMoistMethod) )
    end select


    ! Identification of calculation method for gravity wave drag
    !
    call MessageNotify( 'M', prog_name, 'GWDMethod=<%c>.', c1 = trim(GWDMethod) )
    !
    select case ( GWDMethod )
    case ( 'None' )
      IDGWDMethod = IDGWDMethodNone
    case ( 'M1987' )
      IDGWDMethod = IDGWDMethodM1987
    case default
      call MessageNotify( 'E', prog_name, 'GWDMethod=<%c> is not supported.', c1 = trim(GWDMethod) )
    end select


    ! Identification of calculation method for dry convection
    !
    call MessageNotify( 'M', prog_name, 'DCMethod=<%c>.', c1 = trim(DCMethod) )
    !
    select case ( DCMethod )
    case ( 'None' )
      IDDCMethod = IDDCMethodNone
    case ( 'DCA' )
      IDDCMethod = IDDCMethodDCA
    case default
      call MessageNotify( 'E', prog_name, 'DCMethod=<%c> is not supported.', c1 = trim(DCMethod) )
    end select



    ! è¨�ç®��¢ã�¼ã����¡¨ç¤�
    ! Display calculation mode
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      briefexpldyn = 'used'
    case ( IDDynModeNoHorAdv )
      briefexpldyn = 'not used'
    case ( IDDynModeTWPICE )
      briefexpldyn = 'TWP-ICE'
    end select

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      briefexplphy = 'parameterization suite is used'
    case ( IDPhysModeHS94 )
      briefexplphy = 'forcing for Held and Suarez (1994) dynamical core test'
    case ( IDPhysModeVenusSimple )
      briefexplphy = 'simple forcing for a Venus-like planet'
    case ( IDPhysModeJupiterSimple )
      briefexplphy = 'simple forcing for a Jupiter-like planet'
    case ( IDPhysModeJupiterSimpleV2 )
      briefexplphy = 'simple forcing for a Jupiter-like planet Ver. 2'
    case ( IDPhysModeNoPhysics )
      briefexplphy = 'not used'
    end select

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )
        briefexplrad = 'dennou AGCM5 default'
      case ( IDRadMethodEarthV2 )
        briefexplrad = 'Earth'
      case ( IDRadMethodMarsV1 )
        briefexplrad = 'Mars'
      case ( IDRadMethodSL09 )
        briefexplrad = 'Schneider and Liu (2009) (Jupiter-like planet)'
!!$      case ( IDRadMethodVenusSimple )
!!$        briefexplrad = 'Venus simple radiation model (tentative)'
      case ( IDRadMethodSimple )
        briefexplrad = 'Simple radiation model (tentative)'
      case ( IDRadMethodRRTMG )
        briefexplrad = 'RRTMG'
      case ( IDRadMethodNone )
        briefexplrad = 'none'
      case default
        call MessageNotify( 'E', 'Unexpected error in setting briefexplrad', '' )
      end select

      select case ( IDSfcFluxMethod )
      case ( IDSfcFluxMethodL82 )
        briefexplsfcflux = 'surface flux by the method of Louis et al. (1982)'
      case ( IDSfcFluxMethodBH91B94 )
        briefexplsfcflux = 'surface flux by the 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 = 'vertical diffusion with the method of Mellor and Yamada level 2'
      case ( IDVDiffMethodMY25 )
        briefexplvdiff = 'vertical diffusion with the method of Mellor and Yamada level 2.5'
      case ( IDVDiffMethodJMA )
        briefexplvdiff = 'vertical diffusion with the method of JMA physical library of 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 ( 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 ( IDPhysMode == IDPhysModeFullPhysics ) then
      call MessageNotify( 'M', prog_name, '|    Radiation model : %c', c1 = trim(briefexplrad) )
      call MessageNotify( 'M', prog_name, '|    Surface flux method : %c', c1 = trim(briefexplsfcflux) )
      call MessageNotify( 'M', prog_name, '|    Vertical diffusion method : %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, xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, 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_SurfMajCompIceB(0:imax-1, 1:jmax) )
      allocate( xy_SoilMoistB     (0:imax-1, 1:jmax) )
      allocate( xy_SurfSnowB      (0:imax-1, 1:jmax) )

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

      allocate( xy_SurfMajCompIceA(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)
      call RestartSurfTempInit
      call RestartSurfTempGet( xy_SurfTemp, xyz_SoilTemp, xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN )

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

!!$      xy_SurfMajCompIceN = 0.0_DP
!!$      xy_SurfMajCompIceB = xy_SurfMajCompIceN

    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 )

    ! Output frequently used variables
    ! Output frequently used variables
    !
    call OutputFreqUsedVarsInit


    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimpleV2 )

      ! ���¹ã�������¼ã�¿å�ºå�����������¸ã����°ç�»é��
      ! Register of variables for history data output
      !
      call HistoryAutoAddVariable( 'Rain', (/ 'lon ', 'lat ', 'time' /), 'precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'Snow', (/ 'lon ', 'lat ', 'time' /), 'precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCP', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation', 'kg m-2 s-1' )

    case ( IDPhysModeFullPhysics )

      ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿ã���¡ã�¤ã����������
      ! 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( 'SurfMajCompIce' , (/ 'lon ', 'lat ', 'time' /), 'surface major component ice', 'kg m-2' )     ! (in)

      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 (rain)', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'Snow', (/ 'lon ', 'lat ', 'time' /), 'precipitation (snow)', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'PRCP', (/ 'lon ', 'lat ', 'time' /), 'precipitation', 'kg m-2 s-1' )

      call HistoryAutoAddVariable( 'RainCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPCum', (/ 'lon ', 'lat ', 'time' /), 'convective precipitation', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'RainLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (rain)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'SnowLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation (snow)', 'kg m-2 s-1' )
      call HistoryAutoAddVariable( 'PRCPLsc', (/ 'lon ', 'lat ', 'time' /), 'non-convective precipitation', 'kg m-2 s-1' )

      ! ���¹ã�������¼ã�¿å�ºå�� (�¹ã�¿ã�¼ã������)
      ! History data output (Start time)
      !
      call HistoryAutoPut( TimeN, 'SurfMajCompIce', xy_SurfMajCompIceN )
      call HistoryAutoPut( TimeN, 'SoilMoist'     , xy_SoilMoistN      )
      call HistoryAutoPut( TimeN, 'SurfSnow'      , xy_SurfSnowN       )
    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( xyz_OMG(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) )
    allocate( xyz_Exner    (0:imax-1, 1:jmax, 1:kmax) )

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      allocate( xy_SurfAlbedo         (0:imax-1, 1:jmax) )
      allocate( xy_SurfHumidCoef      (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLenMom    (0:imax-1, 1:jmax) )
      allocate( xy_SurfRoughLenHeat   (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_SurfHeightStd      (0:imax-1, 1:jmax) )
      allocate( xy_SnowFrac           (0:imax-1, 1:jmax) )

!!$      allocate( xy_FlagMatthewsLand         (0:imax-1, 1:jmax) )
      allocate( xy_PhyImplSDHIndexCalcMethod(0:imax-1, 1:jmax) )
      allocate( xy_BucketFlagOceanGrid      (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( 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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixDiffCoef(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( xy_SurfMOLength     (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_DUDtGWD(0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_DVDtGWD(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( xy_RainCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_RainLsc      (0:imax-1, 1:jmax) )
      allocate( xy_Snow         (0:imax-1, 1:jmax) )
      allocate( xy_SnowCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_SnowLsc      (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQVapDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OSolDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DUDtCum       (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DVDtCum       (0:imax-1,1:jmax,1:kmax) )

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

      allocate( xyz_MoistConvDetTend       (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyz_MoistConvSubsidMassFlux(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) )
      allocate( xyz_CloudCoverforRad(0:imax-1,1:jmax,1:kmax) )

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

      ! ���¹ã�������¼ã�¿å�ºå�����������¸ã����°ç�»é��
      ! 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' )

      call HistoryAutoAddVariable( 'SurfRoughLenMom' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length for momentum', 'm' )

      call HistoryAutoAddVariable( 'SurfRoughLenHeat' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length for heat', 'm' )

      call HistoryAutoAddVariable( 'SurfDustGravSedFlux' , (/ 'lon ', 'lat ', 'time' /), 'surface roughness length', 'm' )

    case ( IDPhysModeVenusSimple )

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

    case ( IDPhysModeJupiterSimple )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_VirTemp(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_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( 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) )
      allocate( xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax,  0:1) )
      allocate( xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax,  0:1) )


    case ( IDPhysModeJupiterSimpleV2 )

      allocate( xyz_Press  (0:imax-1, 1:jmax, 1:kmax) )
      allocate( xyr_Press  (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_VirTemp(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( xyz_DTempDtVDiff(0:imax-1, 1:jmax, 1:kmax) )

      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_VelDiffCoef (0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_TempDiffCoef(0:imax-1, 1:jmax, 0:kmax) )
      allocate( xyr_QMixDiffCoef(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_RadLFlux      (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( 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) )
      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( xy_Rain         (0:imax-1, 1:jmax) )
      allocate( xy_RainCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_RainLsc      (0:imax-1, 1:jmax) )
      allocate( xy_Snow         (0:imax-1, 1:jmax) )
      allocate( xy_SnowCumulus  (0:imax-1, 1:jmax) )
      allocate( xy_SnowLsc      (0:imax-1, 1:jmax) )

      allocate( xyz_DTempDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQVapDtCum    (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OLiqDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DQH2OSolDtCum (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DUDtCum       (0:imax-1,1:jmax,1:kmax) )
      allocate( xyz_DVDtCum       (0:imax-1,1:jmax,1:kmax) )

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

    end select







    ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
    ! Setting planetary surface properties
    !
    select case ( IDSfcMoistMethod )
    case ( IDSfcMoistMethodNone )
      call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .false., FlagSnow )
    case ( IDSfcMoistMethodBucket )
      call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .true. , FlagSnow )
    end select


    select case ( IDPhysMode )
    case ( IDPhysModeNoPhysics )

      ! Nothing to do

    case ( IDPhysModeHS94 )

      ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
      ! Forcing and dissipation suggested by Held and Suarez (1994)
      !
      call HS94Init

    case ( IDPhysModeVenusSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! Yamamoto and Takahashi (2003) ����£ã��ç°¡å������è¨�ç®���������¼·��
      ! forcing for simple Venus calculation following Yamamoto and Takahashi (2003)
      !
      call YT2003ForcingInit

    case ( IDPhysModeJupiterSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! Schneider and Liu (2009) ���������´æ··��課ç�
      ! Vertical diffusion by Schneider and Liu (2009)
      !
      call SL09DiffusionInit

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Init

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    case ( IDPhysModeJupiterSimpleV2 )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! ����������������
      ! Lower boundary flux
      !
      call LBFluxSimpleInit

      ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
      ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
      !
      call PhyImplAtmOnlyInit

      ! Schneider and Liu (2009) ���¾å��¢ã����
      ! Radiation model by Schneider and Liu (2009)
      !
      call RadSL09Init

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    case ( IDPhysModeFullPhysics )

      ! ����è¡��¢ã���¼ã�¿ã��¨­å®�
      ! Setting planetary surface properties
      !
      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodNone )
        call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .false., FlagSnow )
      case ( IDSfcMoistMethodBucket )
        call SurfacePropertiesInit( FlagPhysImpSoilModelSO, .true. , FlagSnow )
      end select

      ! ��, æ°·ã���²å��
      ! snow/ice fraction
      !
      call SnowIceFracInit

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDRadMethod )
      case ( IDRadMethodDennouAGCM )

        ! �¾å����������� (GFD �»è�³å�¶æ¥½�����ºã���¾å��¢ã����)
        ! Radiation flux (radiation model developed by GFD Dennou Club)
        !
        call RadDennouAGCMInit( flag_rst = .not. flag_initial )

      case ( IDRadMethodEarthV2 )

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

        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          ! �²ã�����¢ã����
          ! No cloud model
          !
          call CloudNoneInit( FlagSnow )
        case ( IDCloudMethodSimple )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodSimpleWithIce )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodT1993WithIce )
          ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
          ! Cloud model based on Tiedtke (1993)
          !
          call CloudT1993baseInit( FlagSnow )
        end select

        ! �°ç��大æ������¾å��¢ã���� Ver. 2
        ! radiation model for the Earth's atmosphere Ver. 2
        !
        call RadEarthV2Init( FlagSnow )

      case ( IDRadMethodRRTMG )

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

        select case ( IDCloudMethod )
        case ( IDCloudMethodNone )
          ! �²ã�����¢ã����
          ! No cloud model
          !
          call CloudNoneInit( FlagSnow )
        case ( IDCloudMethodSimple )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodSimpleWithIce )
          ! ç°¡å���²ã�¢ã����
          ! Simple cloud
          !
          call CloudSimpleInit( FlagSnow )
        case ( IDCloudMethodT1993WithIce )
          ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
          ! Cloud model based on Tiedtke (1993)
          !
          call CloudT1993baseInit( FlagSnow )
        end select

        ! wrapper of RRTMG
        ! wrapper of RRTMG
        !
        call RadRRTMGWrapperInit( FlagSnow )

      case ( IDRadMethodMarsV1 )

        ! ����大æ������¾å��¢ã���� Ver. 1
        ! radiation model for the Mars' atmosphere Ver. 1
        !
        call RadMarsV1Init

      case ( IDRadMethodSL09 )

        ! Schneider and Liu (2009) ���¾å��¢ã����
        ! Radiation model by Schneider and Liu (2009)
        !
        call RadSL09Init

      case ( IDRadMethodSimple )

        ! ç°¡å���¾å��¢ã����
        ! Simple radiation model
        !
        call RadSimpleInit

      case ( IDRadMethodNone )

        ! ä½����������¾å��¢ã����
        ! radiation model with no absorption and no scattering
        !
        call RadNoneInit

      end select


      ! ���´æ�¡æ�£ã����������
      ! Vertical diffusion flux
      !
      select case ( IDVDiffMethod )
      case ( IDVDiffMethodMY2, IDVDiffMethodMY25 )

        ! ���´æ�¡æ�£ã���������� (Mellor and Yamada, 1974)
        ! Vertical diffusion flux (Mellor and Yamada, 1974)
        !
        call VDiffusionInit

      case ( IDVDiffMethodJMA )
        ! JMA ä¹±æ�æ··å���¢ã�¸ã�¥ã�¼ã��
        ! JMA turbulent mixing module
        !
        call VDiffusionJMAInit

      end select

      ! �°è¡¨�¢ã���������� (�������)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

      !
      ! set dust flux
      !
!!$      call SetDustFluxInit

      ! ä¸������������������å¤��������ç®� (�°è§£æ³�)
      ! Calculate tendency by a part of physical processes (implicit)
      !
      select case ( IDPhyTendMethod )
      case ( IDPhyTendMethodImp1LayModel )

        ! �°è§£æ³�����������ç©���
        ! Time integration with implicit scheme
        !
        call PhyImplInit( FlagBucketModel, FlagSnow )

      case ( IDPhyTendMethodImpSoilModel )

        select case ( IDSfcMoistMethod )
        case ( IDSfcMoistMethodBucket )
          ! ���±ã���¢ã����
          ! Bucket model
          !
          call BucketModelInit( FlagSnow )

        end select

        ! �°ä����������±ã�����´æ�¡æ��
        ! Vertical diffusion of heat under the ground
        !
        call SubsurfaceDiffusionInit

        ! �°ä��±ä�å°��¢ã�������������´å�����°è§£æ³�����������ç©���
        !
        ! Time integration by using implicit scheme in case using subsurface thermal diffusion model
!!$        call PhyImplSDHInit(                                                 &
!!$          & FlagBucketModel, FlagSnow,                                       &
!!$          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
!!$          & )

!!$        call PhyImplSDHV2Init(                                                &
!!$          & FlagBucketModel, FlagSnow,                                       &
!!$          & FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName  &
!!$          & )
        call PhyImplSDHV3Init( FlagBucketModel, FlagSnow, FlagPhysImpSoilModelSO, FlagMajCompPhaseChange, CondMajCompName )

      case ( IDPhyTendMethodImpAtmOnly )

        ! �°è§£æ³�����������ç©��� (大æ����� / ����è¡��¢æ¸©åº��»å��å£�æ¸�º¦è¨�ç®�����)
        ! Time integration by using implicit scheme in case without calculation of surface and soil temperature
        !
        call PhyImplAtmOnlyInit

      end select


      ! Gravity wave drag by McFarlane (1987)
      ! Gravity wave drag by McFarlane (1987)
      !
      call GWDM1987Init

      ! �°è§£æ³�����������ç©��������������¼ã����
      ! Routines for time integration with implicit scheme
      !
      call PhyImplUtilsInit

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit


!!$      select case ( IDRadMethod )
!!$      case ( IDRadMethodMarsV1 )
!!$        ! (����大æ�����) Non-LTE �¾å��¢ã����
!!$        ! Non-NLTE radiation model (for the Mars' atmosphere)
!!$        !
!!$        call Rad15mNLTEInit
!!$      end select

      select case ( IDRadMethod )
      case ( IDRadMethodMarsV1 )
        ! ����è¨�ç®����赤å����±ç�è¨�ç®�
        ! Calculation of near infrared heating rate in the case of Mars
        !
        call RadMarsNIRInit
      end select


      ! ���´æ�¡æ�£ã���������� (Mellor and Yamada, 1974)
      ! Vertical diffusion flux (Mellor and Yamada, 1974)
      !
      call VDiffusionInit

      ! �°è¡¨�¢ã���������� (�������)
      ! Surface flux (Bulk method)
      !
      call SurfaceFluxInit

      ! �¾å��¢é�£ã���¼ã����
      ! Routines for radiation calculation
      !
      call RadUtilsInit

    end select



    ! �°é�¢æ¸©åº��»å��å£�æ¸�º¦�»å��å£�æ°´å���»ç���������
    ! Time integration of surface temperature, soil temperature, soil moisture, 
    ! and surface snow amount
    !
    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )

      ! �°é�¢æ¸©åº�, ��å£�æ¸�º¦������ç©���
      ! Time integration of surface temperature and soil temperature
      !
      call IntgSurfTempInit

      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodBucket )
        ! ���±ã���¢ã����
        ! Bucket model
        !
        call BucketModelInit( FlagSnow )
      end select

    end select




    ! ������
    ! Dynamical core
    !
    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
      ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
      !
      call DynamicsHSplVAS83Init
    case ( IDDynModeNoHorAdv )
      ! ����������¿ã���ç®�����������å­����
      ! A dynamics for calculation with physical processes only
      !
      call DynamicsPhysicsOnlyInit
    case ( IDDynModeTWPICE )
      !
      ! Dynamical process for TWP-ICE experiment
      !
      call DynamicsTWPICESCMExpInit
    end select


    select case ( IDPhysMode )
    case ( IDPhysModeJupiterSimple )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

    case ( IDPhysModeJupiterSimpleV2 )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      ! 湿潤対æ�調ç�
      ! Moist convective adjustment
      !
      call MoistConvAdjustInit

      ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
      ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
      !
      call LScaleCondInit( FlagSnow )

      ! ç°¡å���²ã�¢ã����
      ! Simple cloud
      !
      call CloudSimpleInit( FlagSnow )

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

      ! ������
      ! Mass fixer
      !
      call MassFixerInit

    case ( IDPhysModeFullPhysics )

      ! è£���������°ã��è¨�ç®������µã�����¼ã���³ã�»é�¢æ�°ç¾¤
      ! Subroutines and functions for calculating auxiliary variables
      !
      call AuxVarsInit

      select case ( IDMCMethod )
      case ( IDMCMethodMCA )
        ! 湿潤対æ�調ç�
        ! Moist convective adjustment
        !
        call MoistConvAdjustInit
      case ( IDMCMethodRAS )
        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        call RASInit
      case ( IDMCMethodRASWithIce )
        ! Relaxed Arakawa-Schubert scheme
        ! Relaxed Arakawa-Schubert scheme
        !
        call RASInit
      end select

      select case ( IDLSCMethod )
      case ( IDLSCMethodM65 )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
!!$      case ( IDLSCMethodLL91 )
!!$        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
!!$        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
!!$        !
!!$        call LScaleCondLL91Init
      case ( IDLSCMethodSatAdjM65 )
        !
        != Saturation adjustment
        !
        if ( CompositionInqIndex( 'H2OLiq' ) <= 0 ) then
          call MessageNotify( 'E', prog_name, 'H2OLiq is not found.' )
        end if
        call SaturationAdjustInit( FlagSnow )
      case ( IDLSCMethodM65WithIce )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Manabe, 1965)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
      case ( IDLSCMethodLL91WithIce )
        ! 大è�模å��çµ� (��対æ��§å��çµ�) (Le Treut and Li, 1991)
        ! Large scale condensation (non-convective condensation) (Le Treut and Li, 1991)
        !
        call LScaleCondInit( FlagSnow )
      end select


      ! 
      ! Cloud model
      !
      select case ( IDCloudMethod )
      case ( IDCloudMethodNone )
        ! �²ã�����¢ã����
        ! No cloud model
        !
        call CloudNoneInit( FlagSnow )
      case ( IDCloudMethodSimple )
        ! ç°¡å���²ã�¢ã����
        ! Simple cloud
        !
        call CloudSimpleInit( FlagSnow )
      case ( IDCloudMethodSimpleWithIce )
        ! ç°¡å���²ã�¢ã����
        ! Simple cloud
        !
        call CloudSimpleInit( FlagSnow )
      case ( IDCloudMethodT1993WithIce )
        ! Tiedtke (1993) ���ºã�¥ã��²ã�¢ã����
        ! Cloud model based on Tiedtke (1993)
        !
        call CloudT1993baseInit( FlagSnow )
      case ( IDCloudMethodMarsH2OCloud )
        ! ���� H2O �²ã�¢ã����
        ! Mars H2O cloud model
        !
        call CloudMarsH2OInit
      end select


      select case ( IDSfcMoistMethod )
      case ( IDSfcMoistMethodBucket )
        ! ���±ã���¢ã����
        ! Bucket model
        !
        call BucketModelInit( FlagSnow )
      end select

      select case ( IDDCMethod )
      case ( IDDCMethodDCA )
        ! ä¹¾ç�¥å�æµ�調ç�
        ! Dry convective adjustment
        !
        call DryConvAdjustInit
      end select

      ! ����������
      ! Gravitational sedimentation process
      !
      call GravSedInit

      ! 主æ�����¸å���
      ! Phase change of atmospheric major component
      !
      call MajorCompPhaseChangeInit( FlagMajCompPhaseChange, CondMajCompName )

      ! ������
      ! Mass fixer
      !
      call MassFixerInit

    end select


    ! �¿ã�¤ã�����£ã���¿ã�� (Asselin, 1972)
    ! Time filter (Asselin, 1972)
    !
!!$    call TimeFiltInit

    ! �������£ã���¿ã�� (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

    ! ��å­���� (�¹ã���������, Arakawa and Suarez (1983))
    ! Dynamical process (Spectral method, Arakawa and Suarez (1983))
    !
    use dynamics_hspl_vas83, only : DynamicsHSplVAS83Finalize

    !
    ! Dynamical process for TWP-ICE experiment
    !
    use dynamics_twpice_scm_exp, only : DynamicsTWPICESCMExpFinalize

    ! Held and Suarez (1994) ������å¼·å�¶ã���£é��
    ! Forcing and dissipation suggested by Held and Suarez (1994)
    !
    use held_suarez_1994, only: Hs94Finalize

    ! �¾å����������� (���³ã���¢ã����)
    ! Radiation flux (band model)
    !
    use rad_DennouAGCM, only: RadDennouAGCMFinalize

    ! 座æ����¼ã�¿è¨­å®�
    ! 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

    ! �°è¡¨�¢æ¸©åº����¹ã�¿ã�¼ã�����¼ã�¿å�¥å�ºå��
    ! 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 ( IDPhysMode == IDPhysModeFullPhysics ) 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( xyz_DTurKinEneDt )

    deallocate( xyz_OMG )

    deallocate( xy_SurfHeight )
    deallocate( xyz_Height )

    if ( IDPhysMode == IDPhysModeFullPhysics ) then
      deallocate( xy_SurfAlbedo          )
      deallocate( xy_SurfHumidCoef       )
      deallocate( xy_SurfRoughLenMom     )
      deallocate( xy_SurfRoughLenHeat    )
      deallocate( xy_SurfHeatCapacity    )
      deallocate( xy_SeaIceConc          )
      deallocate( xy_SurfCond            )
      deallocate( xy_SurfType            )
      deallocate( xy_DeepSubSurfHeatFlux )
      deallocate( xy_SoilHeatCap         )
      deallocate( xy_SoilHeatDiffCoef    )
      deallocate( xy_SurfHeightStd       )
      deallocate( xy_SnowFrac            )

!!$      deallocate( xy_FlagMatthewsLand          )
      deallocate( xy_PhyImplSDHIndexCalcMethod )
      deallocate( xy_BucketFlagOceanGrid       )

      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_VelDiffCoef  )
      deallocate( xyr_TempDiffCoef )
      deallocate( xyr_QMixDiffCoef )

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

      deallocate( xy_SurfMOLength )

      deallocate( xy_DSurfTempDt )
      deallocate( xyz_DSoilTempDt )

      deallocate( xy_DPsDt )

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

      deallocate( xyz_DTempDtVDiff )

      deallocate( xyz_DUDtGWD )
      deallocate( xyz_DVDtGWD )

      deallocate( xyz_DTempDtRadL )
      deallocate( xyz_DTempDtRadS )


      deallocate( xy_Rain          )
      deallocate( xy_RainCumulus   )
      deallocate( xy_RainLsc       )
      deallocate( xy_SnowCumulus   )
      deallocate( xy_SnowLsc       )

      deallocate( xyz_DTempDtCum     )
      deallocate( xyz_DQVapDtCum     )
      deallocate( xyz_DQH2OLiqDtCum  )
      deallocate( xyz_DQH2OSolDtCum  )
      deallocate( xyz_DUDtCum        )
      deallocate( xyz_DVDtCum        )

      deallocate( xyz_DQH2OLiqDtLSC  )
      deallocate( xyz_DQH2OSolDtLSC  )

      deallocate( xyz_MoistConvDetTend        )
      deallocate( xyz_MoistConvSubsidMassFlux )

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



    select case ( IDDynMode )
    case ( IDDynModeHSPLVAS83 )
      ! ���¢ã�¸ã�¥ã�¼ã��������°ã���²ä�解é��
      ! Dellocation of variables in modules
      !
      call DynamicsHSplVAS83Finalize
    case ( IDDynModeNoHorAdv )
    case ( IDDynModeTWPICE )
      !
      ! Dynamical process for TWP-ICE experiment
      !
      call DynamicsTWPICESCMExpFinalize
    end select


    call AuxVarsFinalize

    select case ( IDPhysMode )
    case ( IDPhysModeHS94 )
      call Hs94Finalize
    end select

    select case ( IDPhysMode )
    case ( IDPhysModeFullPhysics )
      ! �²ä�解é�¤ã�����¹ã�¿ã�¼ã�����¡ã�¤ã�����äº�����
      ! Dellocation and close a restart file
      !
      call RadDennouAGCMFinalize
    end select

    call AxessetFinalize

    ! ���»ç���çµ�äº�����
    ! Termination of time control
    !
    call TimesetClose

    ! Finalize MPI
    !
    call MPIWrapperFinalize


  end subroutine MainTerminate