subroutine SetSurfData( xy_SurfTemp, xy_SurfAlbedo, xy_SurfHumidCoef, xy_SurfRoughLength, xy_SurfHeatCapacity, xy_DeepSubSurfHeatFlux, xy_SurfCond, xy_SeaIceConc, xy_SoilHeatCap, xy_SoilHeatDiffCoef )
    !
    ! GCM 用の地表面データを返します.
    !
    ! Return surface data for GCM.
    !

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

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: LChar

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat                 ! $ \varphi $ [rad.] . 緯度. Latitude

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


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

    ! 作業変数
    ! Work variables
    !
!!$    integer:: i               ! 経度方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
!!$    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement

    if ( .not. surface_data_inited ) call SurfDataInit

    select case ( LChar( trim(Pattern) ) )

    case ( 'homogeneous' )
      ! SST 一様
      ! SST is homogeneous
      !
      if ( present(xy_SurfTemp) ) xy_SurfTemp = SurfTemp

    case ( 'hosaka et al. (1998)' )
      ! Hosaka et al. (1998) において用いられた SST
      ! SST used in Hosaka et al. (1998)
      !
      if ( present( xy_SurfTemp ) ) then
        call Hosakaetal98SST( xy_SurfTemp )
      end if

    case ( 'nh01_control' )
      ! Neale and Hoskins (2001) の Control experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2001)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH01SST( 'control', xy_SurfTemp )
      end if

    case ( '1-d profile' )

      if ( present( xy_SurfTemp ) ) then
        call Set1DProfileSurfTemp( xy_SurfTemp )
      end if

    case default
      call MessageNotify( 'E', module_name, 'Pattern=<%c> is invalid.', c1 = trim(Pattern) )
    end select

    if ( present(xy_SurfAlbedo          ) ) xy_SurfAlbedo          = Albedo
    if ( present(xy_SurfHumidCoef       ) ) xy_SurfHumidCoef       = HumidCoef
    if ( present(xy_SurfRoughLength     ) ) xy_SurfRoughLength     = RoughLength
    if ( present(xy_SurfHeatCapacity    ) ) xy_SurfHeatCapacity    = HeatCapacity
    if ( present(xy_DeepSubSurfHeatFlux ) ) xy_DeepSubSurfHeatFlux = TempFlux
    if ( present(xy_SurfCond            ) ) xy_SurfCond            = SurfCond
    if ( present(xy_SeaIceConc          ) ) xy_SeaIceConc          = SeaIceConc
    if ( present(xy_SoilHeatCap         ) ) xy_SoilHeatCap         = SoilHeatCap
    if ( present(xy_SoilHeatDiffCoef    ) ) xy_SoilHeatDiffCoef    = SoilHeatDiffCoef


  end subroutine SetSurfData

  subroutine Hosakaetal98SST( xy_SurfTemp )
    !
    ! GCM 用の地表面データを返します.
    !
    ! Return surface data for GCM.
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat                 ! $ \varphi $ [rad.] . 緯度. Latitude

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: LChar

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(out) :: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature

    ! 作業変数 (Hosaka et al. (1998))
    ! Work variables (Hosaka et al. (1998))
    !
    real(DP):: TempEq         ! 赤道上 (正確には LatCenter 上) での温度.
                              ! Temperature on the equator 
                              ! (on LatCenter, to be exact)
    real(DP):: LatCenter      ! 温度最高の緯度. 
                              ! Latitude on which temperature is maximum.
    real(DP):: LatFlatWidth   ! 温度が平坦化される緯度幅. 
                              ! Latitude width in which temperature is flattened
    integer:: jp
    integer:: jm

    real(DP):: LatA, Alpha, Beta, Gamma

    real(DP):: Phi1, AlphaBeta4, Phi, LatAPlus, LatAMinus
    real(DP):: SurfTempMx

    ! 作業変数
    ! Work variables
    !
!!$    integer:: i               ! 経度方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
!!$    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement

    if ( .not. surface_data_inited ) call SurfDataInit

    ! Hosaka et al. (1998) において用いられた SST
    ! SST used in Hosaka et al. (1998)
    !

!!$    TempEq       = SurfTemp
    TempEq       = 302.0_DP
    LatCenter    =   0.0_DP
    LatFlatWidth =   7.0_DP

    LatA         =  30.0_DP
    Alpha        =  60.0_DP
    Beta         =  32.0_DP
    Gamma        =   0.0_DP

    Phi1 = abs( LatA * PI / 180.0_DP )
    AlphaBeta4  = 2.0_DP *( Phi1**3 ) * ( Beta / Alpha )

    do j = 1, jmax
      Phi = abs( y_Lat(j) - LatCenter * PI / 180.0_DP )
      xy_SurfTemp (:,j) = TempEq - Alpha / 2.0_DP * ( Phi - max( sqrt( Phi1**2 + AlphaBeta4 ) - sqrt( ( Phi - Phi1 )**2 + AlphaBeta4 ), 0.0_DP ) ) + Gamma * ( Phi**3 )
    end do

    ! 中心 LatCenter +/- LatFlatWidth の間を平坦に
    ! Flatten between LatCenter +/- LatFlatWidth
    !
    if ( LatFlatWidth < 0.0_DP ) then
      LatFlatWidth = - LatFlatWidth
    end if
    LatAPlus = ( LatCenter + LatFlatWidth ) * PI / 180.0_DP
    LatAMinus = ( LatCenter - LatFlatWidth ) * PI / 180.0_DP

    jp = 1
    jm = jmax
    do j = 1, jmax
      if ( y_Lat(j) <= LatAPlus ) then
        jp = j
        if ( j == jmax ) jp = jp - 1
      end if
      if ( y_Lat(j) < LatAMinus ) then
        jm = j
        if ( j == jmax ) jm = jm - 1
      end if
    end do

    if ( jmax /= 1 ) then
      SurfTempMx = (   xy_SurfTemp(0,jm) * ( y_Lat(jm+1) - LatAMinus ) + xy_SurfTemp(0,jm+1) * ( LatAMinus - y_Lat(jm) ) ) / ( y_Lat(jm+1) - y_Lat(jm) )

      xy_SurfTemp(:,jm+1:jp) = SurfTempMx
    end if


  end subroutine Hosakaetal98SST

  subroutine InitCheck
    !
    ! 依存モジュールの初期化チェック
    !
    ! Check initialization of dependency modules

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

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_inited

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

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


    ! 実行文 ; Executable statement

    if ( .not. namelist_util_inited ) call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

  end subroutine InitCheck

  recursive subroutine NH01SST( SSTType, xy_SurfTemp )
    !
    ! Set SST described by Neale and Hoskins (2001)
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: x_Lon, y_Lat                 ! $ \varphi $ [rad.] . 緯度. Latitude

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: LChar

    ! 宣言文 ; Declaration statements
    !
    implicit none
    character(len=*), intent(in   ) :: SSTType
    real(DP)        , intent(inout) :: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度.
                              ! Surface temperature

    ! 作業変数
    ! Work variables
    !
    real(DP) :: Temp0         ! Zero degree Celsius
                              ! Latitude width in which temperature is flattened
    real(DP) :: xy_SurfTempTmp1 (0:imax-1, 1:jmax)
    real(DP) :: xy_SurfTempTmp2 (0:imax-1, 1:jmax)

    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
!!$    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement

    Temp0 = 273.15d0

    if ( SSTType == 'control' ) then
      ! Neale and Hoskins (2001) の control experiment において用いられた SST
      ! SST used for control experiment by Neale and Hoskins (2001)
      !

      do j = 1, jmax
        if ( abs( y_Lat(j) ) < PI / 3.0d0 ) then
          xy_SurfTemp(:,j) = 27.0d0 * ( 1.0d0 - sin( 3.0d0 * y_Lat(j) / 2.0d0 )**2 )
        else
          xy_SurfTemp(:,j) = 0.0d0
        end if
      end do
      xy_SurfTemp = xy_SurfTemp + Temp0


    else if ( SSTType == 'peaked' ) then
      ! Neale and Hoskins (2001) の Peaked experiment において用いられた SST
      ! SST used for Peaked experiment by Neale and Hoskins (2001)
      !

      do j = 1, jmax
        if ( abs( y_Lat(j) ) < PI / 3.0d0 ) then
          xy_SurfTemp(:,j) = 27.0d0 * ( 1.0d0 - 3.0d0 * abs( y_Lat(j) ) / PI )
        else
          xy_SurfTemp(:,j) = 0.0d0
        end if
      end do
      xy_SurfTemp = xy_SurfTemp + Temp0

    else if ( SSTType == 'flat' ) then
      ! Neale and Hoskins (2001) の Flat experiment において用いられた SST
      ! SST used for Flat experiment by Neale and Hoskins (2001)
      !

      do j = 1, jmax
        if ( abs( y_Lat(j) ) < PI / 3.0d0 ) then
          xy_SurfTemp(:,j) = 27.0d0 * ( 1.0d0 - sin( 3.0d0 * y_Lat(j) / 2.0d0 )**4 )
        else
          xy_SurfTemp(:,j) = 0.0d0
        end if
      end do
      xy_SurfTemp = xy_SurfTemp + Temp0

    else if ( SSTType == 'control-5n' ) then
      ! Neale and Hoskins (2001) の Control-5N experiment において用いられた SST
      ! SST used for Control-5N experiment by Neale and Hoskins (2001)
      !

      do j = 1, jmax
        if ( y_Lat(j) < - PI / 3.0d0 ) then
          xy_SurfTemp(:,j) = 0.0d0
        else if ( y_Lat(j) < PI / 36.0d0 ) then
          xy_SurfTemp(:,j) = 27.0d0 * ( 1.0d0 - sin( 90.0d0/65.0d0 * ( y_Lat(j) - PI/36.0d0 ) )**2 )
        else if ( y_Lat(j) < PI / 3.0d0 ) then
          xy_SurfTemp(:,j) = 27.0d0 * ( 1.0d0 - sin( 90.0d0/55.0d0 * ( y_Lat(j) - PI/36.0d0 ) )**2 )
        else
          xy_SurfTemp(:,j) = 0.0d0
        end if
      end do
      xy_SurfTemp = xy_SurfTemp + Temp0

    else if ( SSTType == 'qobs' ) then
      ! Neale and Hoskins (2001) の Qobs experiment において用いられた SST
      ! SST used for Qobs experiment by Neale and Hoskins (2001)
      !

      call NH01SST( 'control', xy_SurfTempTmp1 )
      call NH01SST( 'flat', xy_SurfTempTmp2 )
      xy_SurfTemp = ( xy_SurfTempTmp1 + xy_SurfTempTmp2 ) * 0.5d0

!!$    else if ( SSTType == '1keq' ) then
!!$      ! Neale and Hoskins (2001) の 1KEQ experiment において用いられた SST
!!$      ! SST used for 1KEQ experiment by Neale and Hoskins (2001)
!!$      !
!!$
!!$      call NH01SST( &
!!$        & 'control',     & ! (in   )
!!$        & xy_SurfTemp  & ! (inout)
!!$        & )
!!$
!!$      do j = 1, jmax
!!$        do i = 1,
!!$          if ( ( abs( x_Lon(i) - PI ) < PI /  3.0d0 ) .and. &
!!$            &  ( abs( y_Lat(j)      ) < PI / 12.0d0 ) ) then
!!$            xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                                 &
!!$              & + 1.0d0 * cos( PI/2.0d0 * ( x_Lon(i) - PI ) / ( PI ) )**2       &
!!$              &         * cos( PI/2.0d0 * y_Lat(j)          / ( PI ) )**2 )
!!$          end if
!!$        end do
!!$      end do

    else
      call MessageNotify( 'E', module_name, 'SSTType=<%c> is invalid.', c1 = trim(SSTType) )
    end if


  end subroutine NH01SST

  subroutine SurfDataInit

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

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

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

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: LChar

    ! 宣言文 ; Declaration statements
    !
    implicit none

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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /surface_data_nml/ Pattern, SurfTemp, Albedo, HumidCoef, RoughLength, HeatCapacity, TempFlux, SurfCond, SeaIceConc, SoilHeatCap, SoilHeatDiffCoef
          !
          ! デフォルト値については初期化手続 "surface_data#SurfDataInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "surface_data#SurfDataInit" for the default values. 
          !


    ! 実行文 ; Executable statement

    if ( surface_data_inited ) return
    call InitCheck

    ! デフォルト値の設定 (まずは Pattern のみ)
    ! Default values settings (At first, "Pattern" only)
    !
    Pattern      = 'Hosaka et al. (1998)'
    SurfTemp     = 273.15_DP
    Albedo       = 0.15_DP
    HumidCoef    = 1.0_DP
    RoughLength  = 1.0e-4_DP
    HeatCapacity = 0.0_DP
    TempFlux     = 0.0_DP
    SurfCond     = 0
    SeaIceConc   = 0.0_DP

    SoilHeatCap      = 2.1d6
    !   volumetric heat capacity (J m-3 K-1)
    !   Value of Clay for porosity f=0.4, volumetric wetness theta=0.2 in Table 12.3 by
    !   Hillel (2004).
    !   Note that the unit of Table 12.3 of Hillel (2004) would be wrong. Although the
    !   unit in the Table is wrong, the volumetric heat capacity of 2.1d6 J m-3 K-1 is
    !   within the range of typical value of it.

    SoilHeatDiffCoef = 1.2d0
    !   thermal conductivity (W m-1 K-1)
    !   Value of Clay for porosity f=0.4, volumetric wetness theta=0.2 in Table 12.3 by
    !   Hillel (2004).

    !   Reference
    !
    !   Hillet, D.,
    !     Introduction to Environmental Soil Physics,
    !     Elsevier Academic Press, pp494, 2004.


    !   Sample values for Mars
    !     These values were obtained from Kiefer (1976) and Kieffer et al. (1977).
    !   Reference
    !     Kieffer, Science, 194, 1344-1346, 1976.
    !     Kieffer et al., JGR, 82, 4249-4291, 1977.
    !
    !   Standard model: see Kieffer et al. (1977) p. 4286,
    !   albedo,                 A   = 0.25
    !     (Kieffer et al., 1977)
    !   thermal inertia,        I   = 6.5e-3 cal cm-2 s-1/2 K-1 = 272 J m-2 s-1/2 K-1
    !     (Kieffer et al., 1977)
    !   density,                rho = 1.65 g cm-3 = 1650 kg m-3
    !     (Kieffer, 1976)
    !   specific heat,          cp  = 0.14 cal g-1 K-1 = 586 J kg-1 K-1
    !     (Kieffer, 1976)
    !
    !   heat capacity,          cp*rho = 0.97e6 J m-3 K-1
    !   conduction coefficient, k = I**2 / (cp*rho) = 7.6e-2   J m-1 s-1 K-1
    !
!!$    SoilHeatCap      = 0.97d6
!!$    SoilHeatDiffCoef = 0.076d0

    ! NOTE:
    !   Values by Kieffer (1976) and Kieffer et al. (1977) would be appropriate for GCM
    !   experiment.


    !   Sample values for Mars
    !     These values were obtained from Savijarvi (1995).
    !   Reference
    !     Savijarvi, H., Mars boundary layer modeling: Diurnal moisture cycle and soil
    !       properties at the Viking lander 1 site, Icarus, 117, 120-127, 1995.
    !
!!$    SoilHeatCap      = 0.8d6
!!$    SoilHeatDiffCoef = 0.18d0


    ! NAMELIST の読み込み (まずは Pattern のみ)
    ! NAMELIST is input (At first, "Pattern" only)
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

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

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

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

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  Pattern          = %c', c1 = trim(Pattern) )
    call MessageNotify( 'M', module_name, '  SurfTemp         = %f', d = (/ SurfTemp         /) )
    call MessageNotify( 'M', module_name, '  Albedo           = %f', d = (/ Albedo           /) )
    call MessageNotify( 'M', module_name, '  HumidCoef        = %f', d = (/ HumidCoef        /) )
    call MessageNotify( 'M', module_name, '  RoughLength      = %f', d = (/ RoughLength      /) )
    call MessageNotify( 'M', module_name, '  HeatCapacity     = %f', d = (/ HeatCapacity     /) )
    call MessageNotify( 'M', module_name, '  TempFlux         = %f', d = (/ TempFlux         /) )
    call MessageNotify( 'M', module_name, '  SurfCond         = %d', i = (/ SurfCond         /) )
    call MessageNotify( 'M', module_name, '  SeaIceConc       = %f', d = (/ SeaIceConc       /) )
    call MessageNotify( 'M', module_name, '  SoilHeatCap      = %f', d = (/ SoilHeatCap      /) )
    call MessageNotify( 'M', module_name, '  SoilHeatDiffCoef = %f', d = (/ SoilHeatDiffCoef /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    surface_data_inited = .true.
  end subroutine SurfDataInit