!= 
!
!= Utility module for advection test
!
! Authors::   Hiroki Kashimura, Yoshiyuki O. Takahashi
! Version::   $Id: dynamics_1d_utils.f90,v 1.1 2015/01/31 06:16:26 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module adv_test
  !
  != 
  !
  != Utility module for advection test
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

  !== References
  !
!!$  !  Chou, M.-D.,
!!$  !    Atmospheric solar heating rate in the water vapor bands,
!!$  !    J. Climate Appl. Meteor., 25, 1532-1542, 1986.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
  ! NAMELIST#set_1d_profile_nml
  !

  ! USE statements
  !

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

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

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only:                              &
    &                    ncmax
                               ! 成分の数
                               ! Number of composition

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

  use constants0, only : PI

  implicit none

  private

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

  integer, save      :: IDVelDist
  integer, parameter :: IDVelDistNCARASPSummerCol = 1
  integer, parameter :: IDVelDistSymHadley        = 2
  integer, parameter :: IDVelDistAsymHadley       = 3

  real(DP), save :: Alpha
  real(DP), save :: Tau
  real(DP), save :: Press0
  real(DP), save :: Omega0
  real(DP), save :: HCSigmaTop
  integer , save :: HCNumCell    ! Number of Hadley cell in each hemisphere
  real(DP), save :: HCV0


  public :: AdvTestSetHorVels
  public :: AdvTestSetVerVel
  public :: AdvTestSetICs
  public :: AdvTestInit


  character(*), parameter:: module_name = 'adv_test'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: dynamics_1d_utils.f90,v 1.1 2015/01/31 06:16:26 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

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

  subroutine AdvTestSetHorVels(   &
    & xyz_UN, xyz_VN              & ! (out)
    & )
    !-------セミラグのテスト用の流速分布を与える--------
    !Gives a velocity for Test. Only used for debug.

    use timeset, only: &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 
    use axesset   , only : x_Lon, y_Lat, z_Sigma, r_Sigma
    use constants , only: RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      &                   Grav

    real(DP), intent(out) :: xyz_UN     (0:imax-1, 1:jmax, 1:kmax)
                              ! 東西風速
                              ! Zonal Wind    
    real(DP), intent(out) :: xyz_VN     (0:imax-1, 1:jmax, 1:kmax)
                              ! 南北風速
                              ! Meridional Wind    

    ! 作業変数
    ! Work variables
    !
    real(DP) :: U0
    real(DP) :: Time
    real(DP) :: StreamFunc

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

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


    select case ( IDVelDist )
    case ( IDVelDistNCARASPSummerCol )
      ! 水平流速分布を与える
      U0 = 2.0_DP * PI * RPlanet / ( 86400.0_DP * 12.0_DP )
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyz_UN(i,j,k) = &
              &   U0 * ( cos( y_Lat(j) ) * cos( Alpha ) + cos( x_Lon(i) ) * sin( y_Lat(j) ) * sin( Alpha ) )
            xyz_VN(i,j,k) = &
              & - U0 * ( sin( x_Lon(i) ) * sin( Alpha ) )
          end do
        end do
      end do
    case ( IDVelDistSymHadley )
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1

            if ( z_Sigma(k) < HCSigmaTop ) then
              StreamFunc = 0.0_DP

              xyz_UN(i,j,k) = 0.0_DP
              xyz_VN(i,j,k) = 0.0_DP
            else

!!$                StreamFunc = &
!!$                  & StreamFunc0 &
!!$                  & * sin( ( y_Lat(j) - 0.0_DP ) / ( HCLatMax - 0.0_DP ) * PI ) &
!!$                  & * sin( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )
!!$
!!$                xyz_UN(i,j,k) = 0.0_DP
!!$                xyz_VN(i,j,k) =                                                 &
!!$                  & StreamFunc0 * Grav / Press0                                 &
!!$                  & * sin( ( y_Lat(j) - 0.0_DP ) / ( HCLatMax - 0.0_DP ) * PI ) &
!!$                  & * PI/(r_Sigma(0)-HCSigmaTop)                                &
!!$                  & * cos( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )

              Time = TimeN

!!$                StreamFunc = &
!!$                  & StreamFunc0 &
!!$                  & * sin( HCNumCell * y_Lat(j) ) &
!!$                  & * sin( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )
!!$
!!$                xyz_UN(i,j,k) = 0.0_DP
!!$                xyz_VN(i,j,k) =                                               &
!!$                  & StreamFunc0 * Grav / Press0                               &
!!$                  & * sin( HCNumCell * y_Lat(j) )                             &
!!$                  & * PI/(r_Sigma(0)-HCSigmaTop)                              &
!!$                  & * cos( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                  & * cos( PI * Time / Tau )

!!$              StreamFunc = &
!!$                &   PI * RPlanet**2 * Press0 / Grav * SigDot0         &
!!$                &   * (   sin( dble( 2 * HCNumCell + 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell + 1 )               &
!!$                &       + sin( dble( 2 * HCNumCell - 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell - 1 ) )             &
!!$                &   * sin( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                &   * cos( PI * Time / Tau )
!!$
!!$              xyz_UN(i,j,k) = 0.0_DP
!!$              xyz_VN(i,j,k) =                                              &
!!$                & - RPlanet / ( 2.0_DP * cos( y_Lat(j) ) ) * SigDot0       &
!!$                &   * PI/(r_Sigma(0)-HCSigmaTop)                           &
!!$                &   * (   sin( dble( 2 * HCNumCell + 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell + 1 )               &
!!$                &       + sin( dble( 2 * HCNumCell - 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell - 1 ) )             &
!!$                & * cos( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                & * cos( PI * Time / Tau )

              StreamFunc =                                            &
                & - 2.0_DP * PI * RPlanet                             &
                &   * (r_Sigma(0)-HCSigmaTop) / PI                    &
                &   * Press0 / Grav * HCV0                            &
                &   * sin( dble( 2 * HCNumCell ) * y_Lat(j) )         &
                &   * cos( y_Lat(j) )**2                              &
                &   * sin( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
                &   * cos( PI * Time / Tau )

              xyz_UN(i,j,k) = 0.0_DP
              xyz_VN(i,j,k) =                                              &
                &   HCV0                                                   &
                &   * sin( dble( 2 * HCNumCell ) * y_Lat(j) )              &
                &   * cos( y_Lat(j) )                                      &
                &   * cos( (z_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
                &   * cos( PI * Time / Tau )

            end if

          end do
        end do
      end do
    case ( IDVelDistAsymHadley )
    end select


  end subroutine AdvTestSetHorVels

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

  subroutine AdvTestSetVerVel( &
    & xyr_SigDotN,             & ! (out)
    & xyr_StreamFunc           & ! (out) optional
    & )
    !-------セミラグのテスト用の流速分布を与える--------
    !Gives a velocity for Test. Only used for debug.

    use axesset, only : y_Lat, r_Sigma

    use timeset, only: &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    use constants , only: RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      &                   Grav

    real(DP), intent(out) :: xyr_SigDotN(0:imax-1, 1:jmax, 0:kmax)
                              ! 鉛直流速（SigmaDot）
    real(DP), intent(out), optional :: xyr_StreamFunc(0:imax-1, 1:jmax, 0:kmax)

    ! 作業変数
    ! Work variables
    !
    real(DP) :: Time
    real(DP) :: Shape
    real(DP) :: StreamFunc

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

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


    select case ( IDVelDist )
    case ( IDVelDistNCARASPSummerCol )
      !鉛直流速分布を与える
      Time = TimeN
      do k = 0, kmax
        Shape = min( 1.0_DP, &
          &          2.0_DP * sqrt( sin( (r_Sigma(k)-r_Sigma(kmax))/(r_Sigma(0)-r_Sigma(kmax)) * PI ) ) &
          &        )
        do j = 1, jmax
          do i = 0, imax-1
            xyr_SigDotN(i,j,k) = Omega0 / Press0 * cos( 2.0_DP*PI/Tau*Time )*sin( Shape*PI/2.0_DP )
          end do
        end do
      end do

      if ( present( xyr_StreamFunc ) ) then
        xyr_StreamFunc = -999.0_DP
      end if

    case ( IDVelDistSymHadley )

      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( r_Sigma(k) < HCSigmaTop ) then
              StreamFunc = 0.0_DP

              xyr_SigDotN(i,j,k) = 0.0_DP
            else

!!$                StreamFunc = &
!!$                  & StreamFunc0 &
!!$                  & * sin( ( y_Lat(j) - 0.0_DP ) / ( HCLatMax - 0.0_DP ) * PI ) &
!!$                  & * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )
!!$
!!$                xyr_SigDotN(i,j,k) = &
!!$                  & - StreamFunc0 * Grav / Press0 / ( RPlanet * cos( y_Lat(j) ) )     &
!!$                  & * (   PI / ( HCLatMax - 0.0_DP ) &
!!$                  &       * cos( ( y_Lat(j) - 0.0_DP ) / ( HCLatMax - 0.0_DP ) * PI ) &
!!$                  &       * cos( y_Lat(j) )   &
!!$                  &     -   sin( ( y_Lat(j) - 0.0_DP ) / ( HCLatMax - 0.0_DP ) * PI ) &
!!$                  &       * sin( y_Lat(j) ) ) &
!!$                  & * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )

              Time = TimeN

!!$              StreamFunc = &
!!$                & StreamFunc0 &
!!$                & * sin( HCNumCell * y_Lat(j) ) &
!!$                & * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI )
!!$
!!$              xyr_SigDotN(i,j,k) = &
!!$                & - StreamFunc0 * Grav / Press0 / ( RPlanet * cos( y_Lat(j) ) ) &
!!$                & * (   HCNumCell &
!!$                &       * cos( HCNumCell * y_Lat(j) ) * cos( y_Lat(j) )   &
!!$                &     -   sin( HCNumCell * y_Lat(j) ) * sin( y_Lat(j) ) ) &
!!$                & * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                & * cos( PI * Time / Tau )

!!$              StreamFunc = &
!!$                &   PI * RPlanet**2 * Press0 / Grav * SigDot0         &
!!$                &   * (   sin( dble( 2 * HCNumCell + 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell + 1 )               &
!!$                &       + sin( dble( 2 * HCNumCell - 1 ) * y_Lat(j) ) &
!!$                &           / dble( 2 * HCNumCell - 1 ) )             &
!!$                &   * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                &   * cos( PI * Time / Tau )
!!$
!!$              xyr_SigDotN(i,j,k) = &
!!$                &   SigDot0                                                       &
!!$                &   * cos( dble( 2 * HCNumCell ) * y_Lat(j) )                     &
!!$                &   * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
!!$                &   * cos( PI * Time / Tau )

              StreamFunc =                                            &
                & - 2.0_DP * PI * RPlanet                             &
                &   * (r_Sigma(0)-HCSigmaTop) / PI                    &
                &   * Press0 / Grav * HCV0                            &
                &   * sin( dble( 2 * HCNumCell ) * y_Lat(j) )         &
                &   * cos( y_Lat(j) )**2                              &
                &   * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
                &   * cos( PI * Time / Tau )

              xyr_SigDotN(i,j,k) = &
                & - (r_Sigma(0)-HCSigmaTop) / ( PI * RPlanet )             &
                &   * HCV0                                                 &
                &   * (   dble( 2 * HCNumCell )                            &
                &         * cos( dble( 2 * HCNumCell ) * y_Lat(j) )        &
                &         * cos( y_Lat(j) )                                &
                &       - 2.0_DP &
                &         * sin( dble( 2 * HCNumCell ) * y_Lat(j) )        &
                &         * sin( y_Lat(j) ) )                              &
                &   * sin( (r_Sigma(k)-HCSigmaTop)/(r_Sigma(0)-HCSigmaTop) * PI ) &
                &   * cos( PI * Time / Tau )

            end if

            if ( present( xyr_StreamFunc ) ) then
              xyr_StreamFunc(i,j,k) = StreamFunc
            end if

          end do
        end do
      end do

      xyr_SigDotN(:,:,0   ) = 0.0_DP
      xyr_SigDotN(:,:,kmax) = 0.0_DP


    case ( IDVelDistAsymHadley )

      if ( present( xyr_StreamFunc ) ) then
        xyr_StreamFunc = -999.0_DP
      end if

    end select


  end subroutine AdvTestSetVerVel

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

  subroutine AdvTestSetICs(     &
    & xy_Ps, xyz_Temp,          & ! (in)
    & xyz_U, xyz_V, xyzf_QMix   & ! (out)
    & )
    !-------セミラグのテスト用の流速分布を与える--------
    !Gives a velocity for Test. Only used for debug.

    use axesset   , only : x_Lon, y_Lat
    use constants , only: RPlanet 
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet

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


    real(DP), intent(in ) :: xy_Ps    (0:imax-1, 1:jmax)
                              ! 東西風速
                              ! Zonal Wind
    real(DP), intent(in ) :: xyz_Temp (0:imax-1, 1:jmax, 1:kmax)
                              ! 温度
                              ! Temperature
    real(DP), intent(out) :: xyz_U    (0:imax-1, 1:jmax, 1:kmax)
                              ! 東西風速
                              ! Zonal Wind
    real(DP), intent(out) :: xyz_V    (0:imax-1, 1:jmax, 1:kmax)
                              ! 南北風速
                              ! Meridional Wind
    real(DP), intent(out) :: xyzf_QMix(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! 
                              ! Mass mixing ratio

    ! 作業変数
    ! Work variables
    !
    real(DP) :: Lon0
    real(DP) :: Lat0
    real(DP) :: Height0
    real(DP) :: Height1
    real(DP) :: Height2
    real(DP) :: RScale
    real(DP) :: ZScale

    real(DP) :: xyz_QVap     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xy_SurfHeight(0:imax-1, 1:jmax)
    real(DP) :: xyz_Height   (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: xyz_HorDist(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_VerDist(0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: D1
    real(DP) :: D2

    integer:: i               ! 東西方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in zonal direction
    integer:: j               ! 南北方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in meridional direction
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: n


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


    xyz_QVap      = 0.0_DP
    xy_SurfHeight = 0.0_DP

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    call AuxVars( &
      & xy_Ps, xyz_Temp, xyz_QVap,     & ! (in )
      & xy_SurfHeight = xy_SurfHeight, & ! (in ) optional
      & xyz_Height    = xyz_Height     & ! (out) optional
      & )


    !
    ! Utility module for advection test
    !
    call AdvTestSetHorVels(   &
      & xyz_U, xyz_V          & ! (out)
      & )

    select case ( IDVelDist )
    case ( IDVelDistNCARASPSummerCol )

      Lon0    = 3.0_DP * PI / 2.0_DP
      Lat0    = 0.0_DP

      Height0 = 4.5e3_DP

      RScale = RPlanet / 3.0_DP
      ZScale = 1.0e3_DP

      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyz_HorDist(i,j,k) = &
              &   RPlanet &
              & * acos(   sin( Lat0 ) * sin( y_Lat(j) ) &
              &         + cos( Lat0 ) * cos( y_Lat(j) ) * cos( x_Lon(i) - Lon0 ) )
          end do
        end do
      end do
      xyz_VerDist = xyz_Height - Height0

      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              D2 =   ( xyz_HorDist(i,j,k) / RScale )**2 &
                &  + ( xyz_VerDist(i,j,k) / ZScale )**2
              D1 = min( 1.0_DP, D2 )
              if ( mod(n,2) == 1 ) then
                ! q5
                xyzf_QMix(i,j,k,n) = 1.0_DP / 2.0_DP * ( 1.0_DP + cos( PI * D1 ) )
              else if ( mod(n,2) == 0 ) then
                ! q6
                if ( D2 <= 1.0_DP ) then
                  xyzf_QMix(i,j,k,n) = 1.0_DP
                else
                  xyzf_QMix(i,j,k,n) = 0.0_DP
                end if
                if ( ( xyz_Height(i,j,k) > Height0 ) .and. &
                  &  ( ( Lat0 - 1.0_DP / 8.0_DP < y_Lat(j) ) .and. &
                  &    ( y_Lat(j) < Lat0 + 1.0_DP / 8.0_DP ) ) ) then
                  xyzf_QMix(i,j,k,n) = 0.0_DP
                end if
              else
                xyzf_QMix(i,j,k,n) = 0.0_DP
              end if
            end do
          end do
        end do
      end do

    case ( IDVelDistSymHadley )

      Height1 = 2.0e3_DP
      Height2 = 5.0e3_DP
      Height0 = ( Height1 + Height2 ) / 2.0_DP

      xyz_VerDist = xyz_Height - Height0

      do n = 1, ncmax
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              if ( ( Height1 < xyz_Height(i,j,k) ) .and. &
                &  ( xyz_Height(i,j,k) < Height2 ) ) then
                xyzf_QMix(i,j,k,n) = &
                  & ( 1.0_DP + cos( 2.0_DP * PI * xyz_VerDist(i,j,k) / ( Height2 - Height1 ) ) ) / 2.0_DP
              else
                xyzf_QMix(i,j,k,n) = 0.0_DP
              end if
            end do
          end do
        end do
      end do

    case ( IDVelDistAsymHadley )
    end select



  end subroutine AdvTestSetICs

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

  subroutine AdvTestInit

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

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

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

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

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

    use constants , only: &
      & RPlanet, &
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
      & Grav                  ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration

    ! 宣言文 ; Declaration statements
    !
    character(STRING) :: VelDist
    real(DP) :: AlphaInDeg
!!$    real(DP) :: HCW0
!!$    real(DP) :: Temp0

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


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /adv_test_nml/ &
      & VelDist,       &
      & AlphaInDeg,    &
      & Tau,           &
      & Press0,        &
      & Omega0,        &
      & HCSigmaTop,    &
      & HCNumCell,     &
      & HCV0
          !
          ! デフォルト値については初期化手続 "set_GATE_profile#SetGATEProfileInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "set_GATE_profile#SetGATEProfileInit" for the default values.
          !

    ! デフォルト値の設定
    ! Default values settings
    !
    VelDist       = "NCARASPSummerCol"
    AlphaInDeg    = 0.0_DP
    Tau           = 345600.0_DP
    Press0        = 1000.0e2_DP
    Omega0        = PI * 4.0e4_DP / Tau
    HCSigmaTop    = 0.1_DP
    HCNumCell     = 3
    HCV0          = 10.0_DP


    ! 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 = adv_test_nml,             & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

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


    select case ( VelDist )
    case ( "NCARASPSummerCol" )
      IDVelDist = IDVelDistNCARASPSummerCol
    case ( "SymHadley" )
      IDVelDist = IDVelDistSymHadley
    case ( "AsymHadley" )
      IDVelDist = IDVelDistAsymHadley
    case default
      call MessageNotify( 'E', module_name, 'VelDist of %c is not supported.', c1 = trim( VelDist ) )
    end select

    Alpha    = AlphaInDeg * PI / 180.0_DP


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


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'VelDist       = %c', c1 = trim( VelDist ) )
    call MessageNotify( 'M', module_name, 'Alpha         = %f', d = (/ Alpha /) )
    call MessageNotify( 'M', module_name, 'Tau           = %f', d = (/ Tau /) )
    call MessageNotify( 'M', module_name, 'Press0        = %f', d = (/ Press0 /) )
    call MessageNotify( 'M', module_name, 'Omega0        = %f', d = (/ Omega0 /) )
    call MessageNotify( 'M', module_name, 'HCSigmaTop    = %f', d = (/ HCSigmaTop /) )
    call MessageNotify( 'M', module_name, 'HCNumCell     = %d', i = (/ HCNumCell /) )
!!$    call MessageNotify( 'M', module_name, 'InFileNameForcing  = %c', c1 = trim(InFileNameForcing ) )
    call MessageNotify( 'M', module_name, 'HCV0          = %f', d = (/ HCV0 /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    adv_test_inited = .true.


  end subroutine AdvTestInit

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

end module adv_test