!= TWP-ICE SCM 実験用力学過程
!
!= dynamical process for TWP-ICE SCM Experiment
!
! Authors::   Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dynamics_twpice_scm_exp.f90,v 1.1 2025/12/17 00:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module dynamics_twpice_scm_exp

  !
  != dcpam 地球大気向け短波放射モデル Ver. 2
  !
  != dcpam short wave radiation model for the Earth's atmosphere Ver. 1
  !
  ! <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.

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

  implicit none

  private


  character(len=STRING), save :: InFileNameSounding
  character(len=STRING), save :: InFileNameForcing

  integer              , save :: Inkmax
  integer              , save :: Intmax
  real(DP), allocatable, save :: z_InPress(:)
  real(DP), allocatable, save :: a_InTime (:)

  real(DP), save :: TimeAtDataStart

  logical              , save :: FlagDynExp

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

  public :: DynamicsTWPICESCMExp
  public :: DynamicsTWPICESCMExpInit
  public :: DynamicsTWPICESCMExpFinalize

  character(*), parameter:: module_name = 'dynamics_twpice_scm_exp'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: dynamics_twpice_scm_exp.f90,v 1.1 2015/02/11 11:53:16 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

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

  subroutine DynamicsTWPICESCMExp( &
    & Time,                                                      &
    & xyz_Press, xyz_Exner, xyz_Height,                          &
    & xyz_DUDtPhy, xyz_DVDtPhy, xyz_DTempDtPhy, xyzf_DQMixDtPhy, &
    & xy_PsB, xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB,             &
    & xy_PsA, xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA              &
    & )

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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Omega
                              ! $ \Omega $ [s-1].
                              ! 回転角速度.
                              ! Angular velocity

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

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

    ! 格子点設定
    ! 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
      &                    IndexH2OVap

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime            ! $ \Delta t $ [s]

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

    ! 1 次元計算用力学過程ユーティリティモジュール
    ! Utility module for dynamics for 1-D calculation
    !
    use dynamics_1d_utils, only : Dynamics1DUtilsVerAdv

    real(DP), intent(in ) :: Time
    real(DP), intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_Exner(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_Height(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_DUDtPhy(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_DVDtPhy(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_DTempDtPhy (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyzf_DQMixDtPhy(0:imax-1,1:jmax,1:kmax,1:ncmax)
    real(DP), intent(in ) :: xy_PsB    (0:imax-1,1:jmax)
    real(DP), intent(in ) :: xyz_UB    (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_VB    (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_TempB (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyzf_QMixB(0:imax-1,1:jmax,1:kmax,1:ncmax)
    real(DP), intent(out) :: xy_PsA    (0:imax-1,1:jmax)
    real(DP), intent(out) :: xyz_UA    (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_VA    (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_TempA (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyzf_QMixA(0:imax-1,1:jmax,1:kmax,1:ncmax)


    !
    ! local variables
    !
    character(len=STRING) :: InFileName
    character(len=STRING) :: InVarName

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

    real(DP) :: xyz_PTempB  (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_PTempA  (0:imax-1,1:jmax,1:kmax)

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

    real(DP) :: xyz_UObs    (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_VObs    (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_WObs    (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DUDtVAdv(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DVDtVAdv(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DPTempDtHAdv(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DQVapDtHAdv (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DPTempDtVAdv(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DQVapDtVAdv (0:imax-1,1:jmax,1:kmax)

    real(DP) :: xyz_DUDt    (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DVDt    (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DPTempDt(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DQVapDt (0:imax-1,1:jmax,1:kmax)

    real(DP), parameter :: UVTimeConst = 2.0_DP * 60.0_DP * 60.0_DP

    real(DP) :: DelTimeX2

    integer :: j
    integer :: k
    integer :: n


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


    DelTimeX2 = 2.0_DP * DelTime


    do k = 1, kmax
      do j = 1, jmax
        xyz_CorPar(:,j,k) = 2.0_DP * Omega * sin( y_Lat(j) )
      end do
    end do


    InFileName = InFileNameSounding
    !
    InVarName  = 'u'
    call DynamicsTWPICESCMExpPTInterp( &
      & InFileName, InVarName, Time,         &
      & xyz_Press,                           &
      & xyz_UObs                             &
      & )
    InVarName  = 'v'
    call DynamicsTWPICESCMExpPTInterp( &
      & InFileName, InVarName, Time,         &
      & xyz_Press,                           &
      & xyz_VObs                             &
      & )

    InFileName = InFileNameForcing
    !
    InVarName  = 'w'
    call DynamicsTWPICESCMExpPTInterp( &
      & InFileName, InVarName, Time,         &
      & xyz_Press,                           &
      & xyz_WObs                             &
      & )
    InVarName  = 'dptdtadvh'
    call DynamicsTWPICESCMExpPTInterp( &
      & InFileName, InVarName, Time,         &
      & xyz_Press,                           &
      & xyz_DPTempDtHAdv                     &
      & )
    InVarName  = 'dqdtadvh'
    call DynamicsTWPICESCMExpPTInterp( &
      & InFileName, InVarName, Time,         &
      & xyz_Press,                           &
      & xyz_DQVapDtHAdv                      &
      & )

    ! 1 次元計算用力学過程ユーティリティモジュール
    ! Utility module for dynamics for 1-D calculation
    !
    call Dynamics1DUtilsVerAdv(  &
      & xyz_WObs, xyz_Height, xyz_UB,  &
      & xyz_DUDtVAdv                   &
      & )
    call Dynamics1DUtilsVerAdv(  &
      & xyz_WObs, xyz_Height, xyz_VB,  &
      & xyz_DVDtVAdv                   &
      & )
    xyz_PTempB = xyz_TempB / xyz_Exner
    call Dynamics1DUtilsVerAdv(  &
      & xyz_WObs, xyz_Height, xyz_PTempB,  &
      & xyz_DPTempDtVAdv                   &
      & )
    call Dynamics1DUtilsVerAdv(  &
      & xyz_WObs, xyz_Height, xyzf_QMixB(:,:,:,IndexH2OVap), &
      & xyz_DQVapDtVAdv                                      &
      & )


    if ( FlagDynExp ) then
!!$      xyz_DUDt = xyz_DUDtVAdv + xyz_CorPar * xyz_VB &
!!$        & + xyz_DUDtPhy - ( xyz_UB - xyz_UObs ) / UVTimeConst
!!$      xyz_DVDt = xyz_DVDtVAdv - xyz_CorPar * xyz_UB &
!!$        & + xyz_DVDtPhy - ( xyz_VB - xyz_VObs ) / UVTimeConst

      ! Horizontal wind velocities are relaxed to observed ones. 
      xyz_DUDt = - ( xyz_UB - xyz_UObs ) / UVTimeConst
      xyz_DVDt = - ( xyz_VB - xyz_VObs ) / UVTimeConst

      xyz_UA = xyz_UB + xyz_DUDt * DelTimeX2
      xyz_VA = xyz_VB + xyz_DVDt * DelTimeX2

    else
!!$      xyz_UT = xyz_UB + ( xyz_DUDtVAdv + xyz_DUDtPhy ) * DelTimeX2
!!$      xyz_VT = xyz_VB + ( xyz_DVDtVAdv + xyz_DVDtPhy ) * DelTimeX2
!!$      !
!!$      xyz_UA = &
!!$        &   (   ( 1.0_DP + DelTimeX2 / UVTimeConst )                  &
!!$        &         * ( xyz_UT + DelTimeX2 / UVTimeConst * xyz_UObs )   &
!!$        &     + ( xyz_CorPar * DelTimeX2 )                            &
!!$        &         * ( xyz_VT + DelTimeX2 / UVTimeConst * xyz_VObs ) ) &
!!$        & / (   ( 1.0_DP + DelTimeX2 / UVTimeConst )**2               &
!!$        &     + ( xyz_CorPar * DelTimeX2 )**2                       )
!!$      xyz_VA = &
!!$        & - (   ( xyz_CorPar * DelTimeX2 )                            &
!!$        &         * ( xyz_UT + DelTimeX2 / UVTimeConst * xyz_UObs )   &
!!$        &     - ( 1.0_DP + DelTimeX2 / UVTimeConst )                  &
!!$        &         * ( xyz_VT + DelTimeX2 / UVTimeConst * xyz_VObs ) ) &
!!$        & / (   ( 1.0_DP + DelTimeX2 / UVTimeConst )**2               &
!!$        &     + ( DelTimeX2 * xyz_CorPar )**2                       )

      xyz_UA = &
        &   ( xyz_UB + xyz_UObs * DelTimeX2 / UVTimeConst ) &
        & / ( 1.0_DP + DelTimeX2 / UVTimeConst )
      xyz_VA = &
        &   ( xyz_VB + xyz_VObs * DelTimeX2 / UVTimeConst ) &
        & / ( 1.0_DP + DelTimeX2 / UVTimeConst )
    end if



    xyz_DPTempDt = xyz_DPTempDtHAdv + xyz_DPTempDtVAdv &
      & + xyz_DTempDtPhy / xyz_Exner
    xyz_DQVapDt  = xyz_DQVapDtHAdv  + xyz_DQVapDtVAdv  &
      & + xyzf_DQMixDtPhy(:,:,:,IndexH2OVap)

    call DynamicsTWPICESCMExpPsTInterp( &
      & Time,                               &
      & xy_PsA                              &
      & )

    xyz_PTempA = xyz_PTempB + xyz_DPTempDt * DelTimeX2
    do n = 1, IndexH2OVap-1
      xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) &
        & + xyzf_DQMixDtPhy(:,:,:,n) * DelTimeX2
    end do
    n = IndexH2OVap
    xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) &
      & + xyz_DQVapDt * DelTimeX2
    do n = IndexH2OVap+1, ncmax
      xyzf_QMixA(:,:,:,n) = xyzf_QMixB(:,:,:,n) &
        & + xyzf_DQMixDtPhy(:,:,:,n) * DelTimeX2
    end do
    xyzf_QMixA = max( xyzf_QMixA, 0.0_DP )

    n = IndexH2OVap
    call AuxVars(                               &
      & xy_PsA, xyz_TempA, xyzf_QMixA(:,:,:,n), &  ! (in )
      & xyz_Exner = xyz_ExnerA                  &  ! (out) optional
      & )
    xyz_TempA = xyz_PTempA * xyz_ExnerA


    call HistoryAutoPut( Time, 'TWPICEWObs' , xyz_WObs )
    call HistoryAutoPut( Time, 'TWPICEDUDtVAdv' , xyz_DUDtVAdv )
    call HistoryAutoPut( Time, 'TWPICEDVDtVAdv' , xyz_DVDtVAdv )
    call HistoryAutoPut( Time, 'TWPICEDPTempDtVAdv' , xyz_DPTempDtVAdv )
    call HistoryAutoPut( Time, 'TWPICEDPTempDtHAdv' , xyz_DPTempDtHAdv )
    call HistoryAutoPut( Time, 'TWPICEDPTempDtPhy' , xyz_DTempDtPhy / xyz_Exner )
    call HistoryAutoPut( Time, 'TWPICEDQVapDtVAdv' , xyz_DQVapDtVAdv )
    call HistoryAutoPut( Time, 'TWPICEDQVapDtHAdv' , xyz_DQVapDtHAdv )
    call HistoryAutoPut( Time, 'TWPICEDQVapDtPhy' , xyzf_DQMixDtPhy(:,:,:,IndexH2OVap) )

  end subroutine DynamicsTWPICESCMExp

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

  subroutine DynamicsTWPICESCMExpPsTInterp( &
    & Time,                                     &
    & xy_Ps                                     &
    & )

    ! MPI
    !
    use mpi_wrapper, only : flag_mpi_init

    ! gtool データ入力
    ! Gtool data input
    !
    use gtool_history, only: HistoryGet

    ! 格子点設定
    ! 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

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

!!$    character(*), intent(in ) :: InFileName
!!$    character(*), intent(in ) :: InVarName
    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(out) :: xy_Ps(0:imax-1,1:jmax)


    !
    ! local variables
    !
    real(DP) :: xyz_Array(0:imax-1,1:jmax,1:Inkmax)
    real(DP) :: z_InArray(1:Inkmax)
    character(len=STRING) :: InFileName
    character(len=STRING) :: InVarName
    real(DP) :: Ps1
    real(DP) :: Ps2
    integer :: k
    integer :: t
    integer :: TIndex

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


    if ( ( Time < a_InTime(1) ) .or. ( Time > a_InTime(Intmax) ) ) then
      call MessageNotify( 'E', module_name, &
        & 'Time of %f is out of range, %f < Time < %f.', &
        & d = (/Time, a_InTime(1), a_InTime(Intmax)/) )
    end if

    do t = 2, Intmax
      if ( Time <= a_InTime(t) ) exit
    end do


    InFileName = InFileNameSounding
    InVarName  = "plev"

    TIndex = t-1
    call HistoryGet(                   &
      & InFileName, InVarName,         & ! (in)
      & xyz_Array,                     & ! (out)
      & range = 'time=^'//toChar(t-1), & ! (in)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )
    do k = 1, Inkmax
      z_InArray(k) = xyz_Array(0,1,k)
    end do
    Ps1 = z_InArray(1)

    TIndex = t
    call HistoryGet(                   &
      & InFileName, InVarName,         & ! (in)
      & xyz_Array,                     & ! (out)
      & range = 'time=^'//toChar(t-1), & ! (in)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )
    do k = 1, Inkmax
      z_InArray(k) = xyz_Array(0,1,k)
    end do
    Ps2 = z_InArray(1)

    xy_Ps = ( Ps2 - Ps1 ) &
      & / ( a_InTime(t) - a_InTime(t-1) ) &
      & * ( Time        - a_InTime(t-1) ) &
      & + Ps1


  end subroutine DynamicsTWPICESCMExpPsTInterp

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

  subroutine DynamicsTWPICESCMExpPTInterp( &
    & InFileName, InVarName, Time,               &
    & xyz_Press,                                 &
    & xyz_Array                                  &
    & )

    ! MPI
    !
    use mpi_wrapper, only : flag_mpi_init

    ! gtool データ入力
    ! Gtool data input
    !
    use gtool_history, only: HistoryGet

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

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

    ! 格子点設定
    ! 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

    ! 1 次元計算用力学過程ユーティリティモジュール
    ! Utility module for dynamics for 1-D calculation
    !
    use dynamics_1d_utils, only : Dynamics1DUtilsVerInterp


    character(*), intent(in ) :: InFileName
    character(*), intent(in ) :: InVarName
    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP)    , intent(out) :: xyz_Array(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP) :: xyz_InArray(0:imax-1,1:jmax,1:Inkmax)
    real(DP) :: z_InArray(1:Inkmax)
    real(DP) :: xyz_Array1(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_Array2(0:imax-1,1:jmax,1:kmax)
    integer :: k
    integer :: t
    integer :: TIndex

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


    if ( ( Time < a_InTime(1) ) .or. ( Time > a_InTime(Intmax) ) ) then
      call MessageNotify( 'E', module_name, &
        & 'Time of %f is out of range, %f < Time < %f.', &
        & d = (/Time, a_InTime(1), a_InTime(Intmax)/) )
    end if

    do t = 2, Intmax
      if ( Time <= a_InTime(t) ) exit
    end do


    TIndex = t-1
    call HistoryGet(                   &
      & InFileName, InVarName,         & ! (in)
      & xyz_InArray,                   & ! (out)
      & range = 'time=^'//toChar(t-1), & ! (in)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )
    do k = 1, Inkmax
      z_InArray(k) = xyz_InArray(0,1,k)
    end do
    call Dynamics1DUtilsVerInterp( &
      & Inkmax, z_InPress, z_InArray,    &
      & xyz_Press,                       &
      & xyz_Array1                       &
      & )

    TIndex = t
    call HistoryGet(                   &
      & InFileName, InVarName,         & ! (in)
      & xyz_InArray,                   & ! (out)
      & range = 'time=^'//toChar(t),   & ! (in)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )
    do k = 1, Inkmax
      z_InArray(k) = xyz_InArray(0,1,k)
    end do
    call Dynamics1DUtilsVerInterp( &
      & Inkmax, z_InPress, z_InArray,    &
      & xyz_Press,                       &
      & xyz_Array2                       &
      & )

    xyz_Array = ( xyz_Array2 - xyz_Array1 ) &
      & / ( a_InTime(t) - a_InTime(t-1) ) &
      & * ( Time        - a_InTime(t-1) ) &
      & + xyz_Array1


  end subroutine DynamicsTWPICESCMExpPTInterp

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

  subroutine DynamicsTWPICESCMExpInit

    ! MPI
    !
    use mpi_wrapper, only : flag_mpi_init

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

    ! gtool データ入力
    ! Gtool data input
    !
    use gtool_history, only: HistoryGet, HistoryGetAttr

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

    ! NetCDF のラッパープログラム
    ! NetCDF wrapper
    !
    use netcdf_wrapper, only : NWInqDimLen, NWGetAtt

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & InitialDate             ! 計算開始日時.
                                ! Start date of calculation

    ! 暦と日時の取り扱い
    ! Calendar and Date handler
    !
    use dc_calendar, only: &
      & DC_CAL_DATE, &          ! 日時を表現するデータ型.
                                ! Data type for date and time
      & DCCalDateCreate, DCCalDateDifference

    ! 1 次元計算用力学過程ユーティリティモジュール
    ! Utility module for dynamics for 1-D calculation
    !
    use dynamics_1d_utils, only : Dynamics1DUtilsInit


    ! 宣言文 ; Declaration statements
    !
    character(len=STRING) :: InFileName

    character(len=STRING) :: Units

    integer  :: TWPICEDataStartYear
    integer  :: TWPICEDataStartMonth
    integer  :: TWPICEDataStartDay
    integer  :: TWPICEDataStartHour
    integer  :: TWPICEDataStartMin
    real(DP) :: TWPICEDataStartSec

    type(DC_CAL_DATE):: TWPICEDataStartDate

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


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

    ! デフォルト値の設定
    ! Default values settings
    !
    InFileNameSounding = 'sounding.nc'
    InFileNameForcing  = 'forcing.nc'

    FlagDynExp = .false.

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

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


    InFileName = InFileNameSounding

    call NWInqDimLen( &
      & InFileName,   & ! (in )
      & 'levelnum',   & ! (in )
      & Inkmax        & ! (out)
      & )
    call NWInqDimLen( &
      & InFileName,   & ! (in )
      & 'time',       & ! (in )
      & Intmax        & ! (out)
      & )

    allocate( z_InPress( Inkmax ) )
    allocate( a_InTime ( Intmax ) )

    call HistoryGet(                   &
      & InFileName, 'plev',            & ! (in)
      & z_InPress,                     & ! (out)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )
    call HistoryGet(                   &
      & InFileName, 'time',            & ! (in)
      & a_InTime,                      & ! (out)
      & flag_mpi_split = flag_mpi_init & ! (in) optional
      & )

    ! Check unit
    call NWGetAtt(  &
      & InFileName, 'time', 'units', & ! (in )
      & Units       & ! (out)
      & )
    if ( ( Units(1:4) /= "hour" ) .and. &
      &  ( Units(1:4) /= "Hour" ) .and. &
      &  ( Units(1:4) /= "HOUR" ) ) then
      call MessageNotify( 'E', module_name, 'Unit of time, %c is inappropriate.', c1 = trim(Units) )
    end if
    ! Unit conversion from hour to second
    a_InTime = a_InTime * ( 60.0_DP * 60.0_DP )

    ! Offset
    TWPICEDataStartYear  = 2006
    TWPICEDataStartMonth =    1
    TWPICEDataStartDay   =   18
    TWPICEDataStartHour  =    0
    TWPICEDataStartMin   =    0
    TWPICEDataStartSec   =    0.0_DP
    call DCCalDateCreate(    &
      & year  = TWPICEDataStartYear,   & ! (in)
      & month = TWPICEDataStartMonth,  & ! (in)
      & day   = TWPICEDataStartDay,    & ! (in)
      & hour  = TWPICEDataStartHour,   & ! (in)
      & min   = TWPICEDataStartMin,    & ! (in)
      & sec   = TWPICEDataStartSec,    & ! (in)
      & date  = TWPICEDataStartDate )    ! (out) optional
    TimeAtDataStart = DCCalDateDifference( &
      &                 start_date = InitialDate,        & ! (in)
      &                 end_date   = TWPICEDataStartDate & ! (in)
      &                 )
    a_InTime = a_InTime + TimeAtDataStart


    ! 1 次元計算用力学過程ユーティリティモジュール
    ! Utility module for dynamics for 1-D calculation
    !
    call Dynamics1DUtilsInit


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'TWPICEWObs', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'vertical wind for TWP-ICE experiment', 'm s-1' )
    call HistoryAutoAddVariable( 'TWPICEDUDtVAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'zonal wind tendency due to vertical advection', 'm s-2' )
    call HistoryAutoAddVariable( 'TWPICEDVDtVAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'meridional wind tendency due to vertical advection', 'm s-2' )
    call HistoryAutoAddVariable( 'TWPICEDPTempDtVAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'potential temperature tendency due to vertical advection', 'K s-1' )
    call HistoryAutoAddVariable( 'TWPICEDPTempDtHAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'potential temperature tendency due to horizontal advection', 'K s-1' )
    call HistoryAutoAddVariable( 'TWPICEDPTempDtPhy', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'potential temperature tendency due to physical processes', 'K s-1' )
    call HistoryAutoAddVariable( 'TWPICEDQVapDtVAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'specific humidity tendency due to vertical advection', 's-1' )
    call HistoryAutoAddVariable( 'TWPICEDQVapDtHAdv', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'specific humidity tendency due to horizontal advection', 's-1' )
    call HistoryAutoAddVariable( 'TWPICEDQVapDtPhy', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'specific humidity tendency due to physical processes', 's-1' )


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'InFileNameSounding = %c', c1 = trim(InFileNameSounding) )
    call MessageNotify( 'M', module_name, 'InFileNameForcing  = %c', c1 = trim(InFileNameForcing ) )
    call MessageNotify( 'M', module_name, 'FlagDynExp         = %b', l = (/ FlagDynExp /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    dynamics_twpice_scm_exp_inited = .true.


  end subroutine DynamicsTWPICESCMExpInit

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

  subroutine DynamicsTWPICESCMExpFinalize


    ! 宣言文 ; Declaration statements
    !


    deallocate( z_InPress )
    deallocate( a_InTime  )


    dynamics_twpice_scm_exp_inited = .false.


  end subroutine DynamicsTWPICESCMExpFinalize

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

end module dynamics_TWPICE_SCM_exp