!= 地表面データ提供
!
!= Prepare surface data
!
! Authors::   Yoshiyuki O. Takahashi, Yasuhiro MORIKAWA, Shin-ichi Takehiro
! Version::   $Id: surface_data.f90,v 1.18 2026/03/10 09:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module surface_data
  !
  != 地表面データ提供
  !
  != Prepare surface data
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! GCM で用いる地表面データを生成します.
  ! 現在は暫定的に Hosaka et al. (1998) の SST 分布を与えます.
  !
  ! Surface data for GCM is generated.
  ! Now, SST profile in Hosaka et al. (1998) is provided tentatively.
  !
  !== Procedures List
  !
  ! SetSurfData   :: 地表面データの取得
  ! ------------  :: ------------
  ! SetSurfData   :: Set surface data
  !
  !== NAMELIST
  !
  ! NAMELIST#surface_data_nml
  !

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

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

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

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP     ! 倍精度実数型. Double precision. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public:: SetSurfData
  public:: SurfDataInit

  ! 公開変数
  ! Public variables
  !
  logical, save :: surface_data_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  ! 非公開変数
  ! Private variables
  !
  character(STRING), save:: Pattern
                              ! 地表面データのパターン. 
                              ! 以下のパターンを選択可能. 
                              ! 
                              ! Surface data pattern. 
                              ! Available patterns are as follows.
                              ! 
                              ! * "Hosaka et al. (1998)"
                              ! * "Homogeneous"
                              ! 
  real(DP), save:: SurfTemp
                              ! 地表面温度の基準値. 
                              ! Standard value of surface temperature
  real(DP), save:: Albedo
                              ! 地表アルベド. 
                              ! Surface albedo
  real(DP), save:: HumidCoef
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
  real(DP), save:: RoughLength
                              ! 地表粗度長. 
                              ! Surface rough length
  real(DP), save:: HeatCapacity
                              ! 地表熱容量. 
                              ! Surface heat capacity
  real(DP), save:: TempFlux
                              ! 地中熱フラックス. 
                              ! Ground temperature flux
  integer, save:: SurfType
                              ! 土地利用.
                              ! Surface index
  integer, save:: SurfCond
                              ! 地表状態 (0: 固定, 1: 可変). 
                              ! Surface condition (0: fixed, 1: variable)
  real(DP), save:: SeaIceConc
                              ! 海氷面密度
                              ! Sea ice concentration
  real(DP), save:: SoilHeatCap
                              ! 土壌熱容量 (J K-1 m-3)
                              ! Specific heat of soil (J K-1 m-3)
  real(DP), save:: SoilHeatDiffCoef
                              ! 土壌熱伝導係数 (W m-1 K-1)
                              ! Heat conduction coefficient of soil (W m-1 K-1)
  real(DP), save:: SurfHeightStd
                              ! 
                              ! Standard deviation of surface height

  character(*), parameter:: module_name = 'surface_data'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: surface_data.f90,v 1.17 2015/02/14 07:26:43 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


contains

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

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

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

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

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $.
                              ! 円周率. Circular constant

    ! ファイルから 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_SurfType (0:imax-1, 1:jmax)
                              ! 土地利用
                              ! Surface index
    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)
    real(DP), intent(out), optional:: xy_SurfHeightStd(0:imax-1, 1:jmax)
                              ! 
                              ! Standard deviation of surface height (m)

    ! 作業変数
    ! 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


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. surface_data_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    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  & ! (out)
          & )
      end if

    case ( 'nh00_control' )
      ! Neale and Hoskins (2000) の control experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & 'control',     & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_peaked' )
      ! Neale and Hoskins (2000) の peaked experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & 'peaked',    & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_flat' )
      ! Neale and Hoskins (2000) の flat experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & 'flat',      & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_control-5n' )
      ! Neale and Hoskins (2000) の control-5n experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & 'control-5n',& ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_qobs' )
      ! Neale and Hoskins (2000) の qobs experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & 'qobs',      & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_1keq' )
      ! Neale and Hoskins (2000) の 1keq experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & '1keq',      & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_3keq' )
      ! Neale and Hoskins (2000) の 3keq experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & '3keq',      & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      end if

    case ( 'nh00_3kw1' )
      ! Neale and Hoskins (2000) の 3kw1 experiment において用いられた SST
      ! SST used for Control experiment by Neale and Hoskins (2000)
      !
      if ( present( xy_SurfTemp ) ) then
        call NH00SST( &
          & '3kw1',      & ! (in   )
          & xy_SurfTemp  & ! (inout)
          & )
      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_SurfType            ) ) xy_SurfType            = SurfType
    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
    if ( present(xy_SurfHeightStd       ) ) xy_SurfHeightStd       = SurfHeightStd


  end subroutine SetSurfData

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

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

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

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

#ifdef LIB_MPI
    use w_mpi_module_base, only: y_Lat_global => y_Lat, js => j1, je => j2
#else
    use w_module_base, only: y_Lat_global => y_Lat
#endif
    
    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, 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) :: y_SurfTempGlobal (1:jmax_global)
                              ! 地表面温度. 
                              ! Surface temperature

    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

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. surface_data_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! 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_global
      Phi = abs( y_Lat_global(j) - LatCenter * PI / 180.0_DP )
      y_SurfTempGlobal (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_global
      if ( y_Lat_global(j) <= LatAPlus ) then
        jp = j
        if ( j == jmax_global ) jp = jp - 1
      end if
      if ( y_Lat_global(j) < LatAMinus ) then
        jm = j
        if ( j == jmax_global ) jm = jm - 1
      end if
    end do

    if ( jmax_global /= 1 ) then
      SurfTempMx = &
        & (   y_SurfTempGlobal(jm) * ( y_Lat_global(jm+1) - LatAMinus ) &
        &   + y_SurfTempGlobal(jm+1) * ( LatAMinus - y_Lat_global(jm) ) &
        & ) &
        & / ( y_Lat_global(jm+1) - y_Lat_global(jm) )

      y_SurfTempGlobal(jm+1:jp) = SurfTempMx
    end if

#ifdef LIB_MPI
    xy_SurfTemp = spread(y_SurfTempGlobal(js:je),1,imax)
#else
    xy_SurfTemp = spread(y_SurfTempGlobal,1,imax)
#endif

  end subroutine Hosakaetal98SST

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

  recursive subroutine NH00SST( &
    & SSTType,     & ! (in   )
    & xy_SurfTemp  & ! (inout)
    & )
    !
    ! Set SST described by Neale and Hoskins (2000)
    !

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

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

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, 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)

    real(DP) :: TAmp
    real(DP) :: LonCen
    real(DP) :: LonWid
    real(DP) :: LatWid

    ! 作業変数
    ! 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.15_DP

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

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


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

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

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

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

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

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

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

      call NH00SST( &
        & 'control',     & ! (in   )
        & xy_SurfTempTmp1  & ! (inout)
        & )
      call NH00SST( &
        & 'flat',     & ! (in   )
        & xy_SurfTempTmp2  & ! (inout)
        & )
      xy_SurfTemp = ( xy_SurfTempTmp1 + xy_SurfTempTmp2 ) * 0.5_DP

    else if ( SSTType == '1keq' ) then
      ! Neale and Hoskins (2000) の 1KEQ experiment において用いられた SST
      ! SST used for 1KEQ experiment by Neale and Hoskins (2000)
      !

      call NH00SST( &
        & 'control',     & ! (in   )
        & xy_SurfTemp  & ! (inout)
        & )

      TAmp   =   1.0_DP
      LonCen = 180.0_DP * PI / 180.0_DP
      LonWid =  30.0_DP * PI / 180.0_DP
      LatWid =  15.0_DP * PI / 180.0_DP
      do j = 1, jmax
        do i = 0, imax-1
          if ( ( abs( x_Lon(i) - LonCen ) < LonWid ) .and. &
            &  ( abs( y_Lat(j)          ) < LatWid ) ) then
            xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                              &
              & + TAmp * cos( PI/2.0_DP * ( x_Lon(i) - LonCen ) / LonWid )**2 &
              &        * cos( PI/2.0_DP * y_Lat(j)              / LatWid )**2
          end if
        end do
      end do

    else if ( SSTType == '3keq' ) then
      ! Neale and Hoskins (2000) の 3KEQ experiment において用いられた SST
      ! SST used for 1KEQ experiment by Neale and Hoskins (2000)
      !

      call NH00SST( &
        & 'control',     & ! (in   )
        & xy_SurfTemp  & ! (inout)
        & )

      TAmp   =   3.0_DP
      LonCen = 180.0_DP * PI / 180.0_DP
      LonWid =  30.0_DP * PI / 180.0_DP
      LatWid =  15.0_DP * PI / 180.0_DP
      do j = 1, jmax
        do i = 0, imax-1
          if ( ( abs( x_Lon(i) - LonCen ) < LonWid ) .and. &
            &  ( abs( y_Lat(j)          ) < LatWid ) ) then
            xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                              &
              & + TAmp * cos( PI/2.0_DP * ( x_Lon(i) - LonCen ) / LonWid )**2 &
              &        * cos( PI/2.0_DP * y_Lat(j)              / LatWid )**2
          end if
        end do
      end do

    else if ( SSTType == '3kw1' ) then
      ! Neale and Hoskins (2000) の 3KW1 experiment において用いられた SST
      ! SST used for 1KEQ experiment by Neale and Hoskins (2000)
      !

      call NH00SST( &
        & 'control',     & ! (in   )
        & xy_SurfTemp  & ! (inout)
        & )

      TAmp   =   3.0_DP
      LonCen = 180.0_DP * PI / 180.0_DP
      LatWid =  30.0_DP * PI / 180.0_DP
      do j = 1, jmax
        do i = 0, imax-1
          if ( abs( y_Lat(j) ) < LatWid ) then
            xy_SurfTemp(i,j) = xy_SurfTemp(i,j)                  &
              & + TAmp * cos( x_Lon(i) - LonCen )                &
              &        * cos( PI/2.0_DP * y_Lat(j) / LatWid )**2
          end if
        end do
      end do

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


  end subroutine NH00SST

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

  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

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


    ! 宣言文 ; 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, SurfType, SurfCond, &
      & SeaIceConc,                                 &
      & SoilHeatCap, SoilHeatDiffCoef,              &
      & SurfHeightStd
          !
          ! デフォルト値については初期化手続 "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


    ! デフォルト値の設定 (まずは Pattern のみ)
    ! Default values settings (At first, "Pattern" only)
    !
!!$    Pattern      = 'Hosaka et al. (1998)'
    Pattern      = 'homogeneous'
    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
    SurfType     = 0
    SurfCond     = 0
    SeaIceConc   = 0.0_DP

    SoilHeatCap      = 2.1e6_DP
    !   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.2e0_DP
    !   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


    SurfHeightStd = 0.0_DP


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

      rewind( unit_nml )
      read( unit_nml, &             ! (in)
        & nml = surface_data_nml, & ! (out)
        & 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, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

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

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


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


    ! 印字 ; 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, '  SurfType         = %d', i = (/ SurfType         /) )
    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, '  SurfHeightStd    = %f', d = (/ SurfHeightStd /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    surface_data_inited = .true.

  end subroutine SurfDataInit

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

end module surface_data