!= 時間フィルター (Williams, 2009)
!
!= Time filter (Williams, 2009)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: timefilter_williams2009.f90,v 1.4 2025/09/17 22:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module timefilter_williams2009
  !
  != 時間フィルター (Williams, 2009)
  !
  != Time filter (Williams, 2009)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 計算モードを抑えるための時間フィルター (Williams, 2009)
  ! を適用するためのモジュールです.
  ! Time filter by Williams (2009) is the modified version of that 
  ! by Asselin (1972).
  !
!!$  ! すなわち
!!$  ! \[
!!$  !   \bar{A}^{t}
!!$  !   = {A}^{t}
!!$  !     + \epsilon_f
!!$  !        \left( 
!!$  !               \bar{A}^{t-\Delta t} - 2{A}^{t} + {A}^{t+\Delta t}
!!$  !        \right) \] \[
!!$  !   \qquad
!!$  !   = ( 1-2 \epsilon_f ) {A}^{t}
!!$  !     + \epsilon_f 
!!$  !        \left( \bar{A}^{t-\Delta t} + {A}^{t+\Delta t} \right)
!!$  ! \]
!!$  ! として $ \bar{A}^{t} $ を求めます. 添え字
!!$  ! $ t-\Delta t, t, t+\Delta t $ はそれぞれ時間ステップを表し,
!!$  ! タイムフィルターが適用された項には上線 ( $ \bar{\quad} $ ) 
!!$  ! が付いています.
!!$  !
!!$  ! タイムフィルターの係数 $\epsilon_f$ は Create の際に設定します.
!!$  !
!!$  !--
!!$  ! \[
!!$  !   \bar{\mathscr A}^{t}
!!$  !     = ( 1-2 \epsilon_f ) {\mathscr A}^{t}
!!$  !     +  \epsilon_f 
!!$  !         \left( \bar{\mathscr A}^{t-\Delta t} + {\mathscr A}^{t+\Delta t} \right)
!!$  ! \]
!!$  ! として$\bar{\mathscr A}^{t}$を求めます. 添え字
!!$  !++
!!$  !
!!$  ! Time filter (Asselin, 1972) for suppression of computational mode 
!!$  ! is applied.
!!$  !
!!$  ! Concretely, $ \bar{A}^{t} $ is derived as follows.
!!$  ! \[
!!$  !   \bar{A}^{t}
!!$  !   = {A}^{t}
!!$  !     + \epsilon_f
!!$  !        \left( 
!!$  !               \bar{A}^{t-\Delta t} - 2{A}^{t} + {A}^{t+\Delta t}
!!$  !        \right) \] \[
!!$  !   \qquad
!!$  !   = ( 1-2 \epsilon_f ) {A}^{t}
!!$  !     + \epsilon_f 
!!$  !        \left( \bar{A}^{t-\Delta t} + {A}^{t+\Delta t} \right),
!!$  ! \]
!!$  ! where suffices $ t-\Delta t, t, t+\Delta t $ represent time step, 
!!$  ! and over-bar $ \bar{\quad} $ represents the terms applied 
!!$  ! time filter.
!!$  !
!!$  ! Time filter coefficient $ \epsilon_f $ is specified by "Create".
  !
  !== Procedures List
  ! 
  ! TimeFilterWilliams2009    :: 時間フィルターの適用
  ! ------------------------  :: ------------
  ! TimeFilterWilliams2009    :: Apply time filter
  !
  !== NAMELIST
  !
  ! NAMELIST#timefilter_williams2009_nml
  !
  !== References
  !
  ! * Williams, P. D., 2009: 
  !   A proposed modification to the Robert-Asselin time filter, 
  !   <i>Mon. Wea. Rev.</i>, <b>137</b>, 2538--2546. 
  !
  ! * Asselin, R. A., 1972: 
  !   Frequency filter for time integrations.
  !   <i>Mon. Wea. Rev.</i>, <b>100</b>, 487--490.
  !

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

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

  use composition, only: ncmax  ! 微量成分の数
                                ! Number of composition

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


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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: TimeFilterWilliams2009
  public :: TimeFilterWilliams2009SurfVars
  public :: TimeFilterWilliams2009Init

  ! 公開変数
  ! Public variables
  !

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

  real(DP), save :: FilterParameterRA
                              ! Robert-Asselin 時間フィルターの係数. 
                              ! Time filter coefficient for Robert-Asselin filter
                              ! This is equal to nu/2 in Williams (2009).
  real(DP), save :: FilterParameterW
                              ! Williams 時間フィルターの係数. 
                              ! Time filter coefficient for Williams filter
                              ! This is equal to alpha in Williams (2009).

  character(*), parameter:: module_name = 'timefilter_williams2009'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: timefilter_williams2009.f90,v 1.4 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

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

  subroutine TimeFilterWilliams2009(                 &
    & xyz_UB, xyz_VB, xyz_TempB, xyzf_QMixB, xy_PsB, &   ! (in)
    & xyz_UN, xyz_VN, xyz_TempN, xyzf_QMixN, xy_PsN, &   ! (inout)
    & xyz_UA, xyz_VA, xyz_TempA, xyzf_QMixA, xy_PsA  &   ! (inout)
    & )
    !
    ! 予報変数にタイムフィルターを適用します. 
    !
    ! Time filter is applied to predictional variables.
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimesetClockStart, TimesetClockStop

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

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only: MassFixer

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_UB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t-\Delta t) $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_VB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t-\Delta t) $ .   南北風速. Northward wind
    real(DP), intent(in):: xyz_TempB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t-\Delta t) $ .   温度. Temperature
    real(DP), intent(in):: xyzf_QMixB (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t-\Delta t) $ .   混合比. Mass mixing ratio of constituents
    real(DP), intent(in):: xy_PsB (0:imax-1, 1:jmax)
                              ! $ p_s (t-\Delta t) $ . 地表面気圧. Surface pressure
    real(DP), intent(inout):: xyz_UN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t) $ .     東西風速. Eastward wind
    real(DP), intent(inout):: xyz_VN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t) $ .     南北風速. Northward wind
    real(DP), intent(inout):: xyz_TempN (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t) $ .     温度. Temperature
    real(DP), intent(inout):: xyzf_QMixN (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t) $ .     混合比. Mass mixing ratio of constituents
    real(DP), intent(inout):: xy_PsN (0:imax-1, 1:jmax)
                              ! $ p_s (t) $ .   地表面気圧. Surface pressure
    real(DP), intent(inout):: xyz_UA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t+\Delta t) $ .   東西風速. Eastward wind
    real(DP), intent(inout):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t+\Delta t) $ .   南北風速. Northward wind
    real(DP), intent(inout):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t+\Delta t) $ .   温度. Temperature
    real(DP), intent(inout):: xyzf_QMixA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio of constituents
    real(DP), intent(inout):: xy_PsA (0:imax-1, 1:jmax)
                              ! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_Delta  (0:imax-1, 1:jmax)
    real(DP) :: xyz_Delta (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyzf_Delta(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)

    !   Variables for mass fix
    real(DP) :: xy_PsNSave    (0:imax-1, 1:jmax)
    real(DP) :: xy_PsASave    (0:imax-1, 1:jmax)
    real(DP) :: xyzf_QMixNSave(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP) :: xyzf_QMixASave(0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
    real(DP) :: xyr_PressNSave(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_PressASave(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_PressN    (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_PressA    (0:imax-1, 1:jmax, 0:kmax)

    ! 実行文 ; Executable statement
    !

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

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! FilterParameterRA の有効性のチェック
    ! Check validity of "FilterParameterRA"
    !
    if ( FilterParameterRA == 0.0_DP ) return


    ! Variables are saved
    xy_PsNSave     = xy_PsN
    xy_PsASave     = xy_PsA
    xyzf_QMixNSave = xyzf_QMixN
    xyzf_QMixASave = xyzf_QMixA


    ! タイムフィルターを適用
    ! Apply time filter
    !
    xyz_Delta  =  &
      & + FilterParameterRA * ( xyz_UB     - 2.0_DP * xyz_UN     + xyz_UA     )
    xyz_UN     =  xyz_UN    &
      & + FilterParameterW              * xyz_Delta
!!$    xyz_UA     =  xyz_UA    &
!!$      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xyz_Delta
    xyz_UA     =  xyz_UA    &
      & - ( 1.0_DP - FilterParameterW ) * xyz_Delta

    xyz_Delta  =  &
      & + FilterParameterRA * ( xyz_VB     - 2.0_DP * xyz_VN     + xyz_VA     )
    xyz_VN     =  xyz_VN    &
      & + FilterParameterW              * xyz_Delta
!!$    xyz_VA     =  xyz_VA    &
!!$      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xyz_Delta
    xyz_VA     =  xyz_VA    &
      & - ( 1.0_DP - FilterParameterW ) * xyz_Delta

    xyz_Delta  =  &
      & + FilterParameterRA * ( xyz_TempB  - 2.0_DP * xyz_TempN  + xyz_TempA  )
    xyz_TempN  =  xyz_TempN &
      & + FilterParameterW              * xyz_Delta
!!$    xyz_TempA  =  xyz_TempA &
!!$      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xyz_Delta
    xyz_TempA  =  xyz_TempA &
      & - ( 1.0_DP - FilterParameterW ) * xyz_Delta

    xy_Delta   =  &
      & + FilterParameterRA * ( xy_PsB     - 2.0_DP * xy_PsN     + xy_PsA     )
    xy_PsN     =  xy_PsN    &
      & + FilterParameterW              * xy_Delta
!!$    xy_PsA     =  xy_PsA    &
!!$      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xy_Delta
    xy_PsA     =  xy_PsA    &
      & - ( 1.0_DP - FilterParameterW ) * xy_Delta

    xyzf_Delta =  &
      & + FilterParameterRA * ( xyzf_QMixB - 2.0_DP * xyzf_QMixN + xyzf_QMixA )
!!$    xyzf_QMixN = xyzf_QMixN + xyzf_Delta
    xyzf_QMixN = xyzf_QMixN &
      & + FilterParameterW              * xyzf_Delta
!    xyzf_QMixA = xyzf_QMixA &
!      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xyzf_Delta
    xyzf_QMixA = xyzf_QMixA &
      & - ( 1.0_DP - FilterParameterW ) * xyzf_Delta
    !
    ! Mass fix
    !   MEMO (1/2):
    !     If mass fixing below is commented out, filter behaves like 
    !     Robert-Asselin filter.
    !   -----
    call AuxVars( &
      & xy_PsNSave, xyz_TempN, xyzf_QMixNSave(:,:,:,IndexH2OVap), & ! (in )
      & xyr_Press = xyr_PressNSave                                & ! (out) optiona
      & )
    call AuxVars( &
      & xy_PsN    , xyz_TempN, xyzf_QMixN    (:,:,:,IndexH2OVap), & ! (in )
      & xyr_Press = xyr_PressN                                    & ! (out) optiona
      & )
    call MassFixer(                    &
      & xyr_PressN,                    & ! (in)
      & xyzf_QMixN,                    & ! (inout)
      & xyr_PressRef = xyr_PressNSave, & ! (in) optional
      & xyzf_QMixRef = xyzf_QMixNSave  & ! (in) optional
      & )
    call AuxVars( &
      & xy_PsASave, xyz_TempA, xyzf_QMixASave(:,:,:,IndexH2OVap), & ! (in )
      & xyr_Press = xyr_PressASave                                & ! (out) optiona
      & )
    call AuxVars( &
      & xy_PsA    , xyz_TempA, xyzf_QMixA    (:,:,:,IndexH2OVap), & ! (in )
      & xyr_Press = xyr_PressA                                    & ! (out) optiona
      & )
    call MassFixer(                    &
      & xyr_PressA,                    & ! (in)
      & xyzf_QMixA,                    & ! (inout)
      & xyr_PressRef = xyr_PressASave, & ! (in) optional
      & xyzf_QMixRef = xyzf_QMixASave  & ! (in) optional
      & )
    !   -----


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine TimeFilterWilliams2009

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

  subroutine TimeFilterWilliams2009SurfVars(           &
    & xy_SurfMajCompIceB, xy_SoilMoistB, xy_SurfSnowB, & ! (in)
    & xy_SurfMajCompIceN, xy_SoilMoistN, xy_SurfSnowN, & ! (inout)
    & xy_SurfMajCompIceA, xy_SoilMoistA, xy_SurfSnowA, & ! (inout)
    & xy_PsA                                           & ! (inout)
    & )
    !
    ! 予報変数にタイムフィルターを適用します. 
    !
    ! Time filter is applied to predictional variables.
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimesetClockStart, TimesetClockStop

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only:  &
      & Grav
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in   ):: xy_SurfMajCompIceB(0:imax-1, 1:jmax)
                              ! $ M_mcs (t-\Delta t) $ . 
                              ! Surface major component ice amount
    real(DP), intent(in   ):: xy_SoilMoistB (0:imax-1, 1:jmax)
                              ! $ M_ws (t-\Delta t) $ . 地表水質量. 
                              ! Soil moisture
    real(DP), intent(in   ):: xy_SurfSnowB  (0:imax-1, 1:jmax)
                              ! $ M_ss (t-\Delta t) $ . 積雪量.
                              ! Surface snow amount
!!$    real(DP), intent(in   ):: xy_SOSeaIceMassB  (0:imax-1, 1:jmax)
!!$                              ! $ M_si (t-\Delta t) $ . 
!!$                              ! Slab sea ice mass
    real(DP), intent(inout):: xy_SurfMajCompIceN(0:imax-1, 1:jmax)
                              ! $ M_mcs (t) $ . 
                              ! Surface major component ice amount
    real(DP), intent(inout):: xy_SoilMoistN (0:imax-1, 1:jmax)
                              ! $ M_ws (t-\Delta t) $ . 地表水質量. 
                              ! Soil moisture
    real(DP), intent(inout):: xy_SurfSnowN  (0:imax-1, 1:jmax)
                              ! $ M_ss (t) $ . 積雪量.
                              ! Surface snow amount
!!$    real(DP), intent(inout):: xy_SOSeaIceMassN  (0:imax-1, 1:jmax)
!!$                              ! $ M_si (t) $ . 
!!$                              ! Slab sea ice mass
    real(DP), intent(inout):: xy_SurfMajCompIceA(0:imax-1, 1:jmax)
                              ! $ M_mcs (t+\Delta t) $ . 
                              ! Surface major component ice amount
    real(DP), intent(inout):: xy_SoilMoistA (0:imax-1, 1:jmax)
                              ! $ M_ws (t-\Delta t) $ . 地表水質量. 
                              ! Soil moisture
    real(DP), intent(inout):: xy_SurfSnowA  (0:imax-1, 1:jmax)
                              ! $ M_ss (t+\Delta t) $ . 積雪量.
                              ! Surface snow amount
!!$    real(DP), intent(inout):: xy_SOSeaIceMassA  (0:imax-1, 1:jmax)
!!$                              ! $ M_si (t+\Delta t) $ . 
!!$                              ! Slab sea ice mass
    real(DP), intent(inout):: xy_PsA (0:imax-1, 1:jmax)
                              ! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xy_Delta  (0:imax-1, 1:jmax)
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude


    ! 実行文 ; Executable statement
    !

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


    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )


    ! FilterParameterRA の有効性のチェック
    ! Check validity of "FilterParameterRA"
    !
    if ( FilterParameterRA == 0.0_DP ) return

    ! タイムフィルターを適用
    ! Apply time filter
    !
    xy_Delta   =  &
      & + FilterParameterRA                                         &
      & * ( xy_SurfMajCompIceB - 2.0_DP * xy_SurfMajCompIceN + xy_SurfMajCompIceA )
    xy_SurfMajCompIceN = xy_SurfMajCompIceN                         &
      & + FilterParameterW              * xy_Delta
!!$    xy_SurfMajCompIceA = xy_SurfMajCompIceA                         &
!!$      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xy_Delta
    xy_SurfMajCompIceA = xy_SurfMajCompIceA                         &
      & - ( 1.0_DP - FilterParameterW ) * xy_Delta
!!$    xy_SurfMajCompIceN = max( xy_SurfMajCompIceN, 0.0_DP )
!!$    xy_SurfMajCompIceA = max( xy_SurfMajCompIceA, 0.0_DP )
    ! Mass fix
    !   MEMO (2/2):
    !     If mass fixing below is commented out, filter behaves like 
    !     Robert-Asselin filter.
    !   -----
    do j = 1, jmax
      do i = 0, imax-1
        if ( xy_SurfMajCompIceA(i,j) < 0.0_DP ) then
          xy_PsA(i,j) = xy_PsA(i,j) + xy_SurfMajCompIceA(i,j) * Grav
          xy_SurfMajCompIceA(i,j) = 0.0_DP
        end if
      end do
    end do
    !   -----

    xy_Delta   =  &
      & + FilterParameterRA                                         &
      & * ( xy_SoilMoistB - 2.0_DP * xy_SoilMoistN + xy_SoilMoistA )
    xy_SoilMoistN = xy_SoilMoistN + xy_Delta
!!$    xy_SoilMoistN = xy_SoilMoistN                                   &
!!$      & + FilterParameterW              * xy_Delta
!!$!    xy_SoilMoistA = xy_SoilMoistA                                   &
!!$!      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xy_Delta
!!$    xy_SoilMoistA = xy_SoilMoistA                                   &
!!$      & - ( 1.0_DP - FilterParameterW ) * xy_Delta
!!$    xy_SoilMoistN = max( xy_SoilMoistN, 0.0_DP )
!!$    xy_SoilMoistA = max( xy_SoilMoistA, 0.0_DP )

    xy_Delta   =  &
      & + FilterParameterRA                                         &
      & * ( xy_SurfSnowB  - 2.0_DP * xy_SurfSnowN  + xy_SurfSnowA  )
    xy_SurfSnowN  = xy_SurfSnowN + xy_Delta
!!$    xy_SurfSnowN  = xy_SurfSnowN                                    &
!!$      & + FilterParameterW              * xy_Delta
!!$!    xy_SurfSnowA  = xy_SurfSnowA                                    &
!!$!      & - ( 1.0_DP - FilterParameterW ) * FilterParameterW * xy_Delta
!!$    xy_SurfSnowA  = xy_SurfSnowA                                    &
!!$      & - ( 1.0_DP - FilterParameterW ) * xy_Delta
!!$    xy_SurfSnowN  = max( xy_SurfSnowN, 0.0_DP )
!!$    xy_SurfSnowA  = max( xy_SurfSnowA, 0.0_DP )

!!$    xy_Delta   =  &
!!$      & + FilterParameterRA                                         &
!!$      & * ( xy_SOSeaIceMassB  - 2.0_DP * xy_SOSeaIceMassN  + xy_SOSeaIceMassA  )
!!$    xy_SOSeaIceMassN  = xy_SOSeaIceMassN + xy_Delta


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine TimeFilterWilliams2009SurfVars

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

  subroutine TimeFilterWilliams2009Init
    !
    ! timefilter_williams2009 モジュールの初期化を行います. 
    ! NAMELIST#timefilter_williams2009_nml の読み込みはこの手続きで行われます. 
    !
    ! "timefilter_williams2009" module is initialized. 
    ! "NAMELIST#timefilter_williams2009_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; 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

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 質量の補正
    ! Mass fixer
    !
    use mass_fixer, only : MassFixerInit

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage

    
    ! 宣言文 ; 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 /timefilter_williams2009_nml/ &
      & FilterParameterRA, &
      & FilterParameterW
          !
          ! デフォルト値については初期化手続 "timefilter_asselin1972#Hs94Init" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "timefilter_asselin1972#Hs94Init" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( timefilter_williams2009_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FilterParameterRA = 0.05_DP
    FilterParameterW  = 0.53_DP
                    ! recommended value for alpha in Williams (2009)

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

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = timefilter_williams2009_nml )
    end if


    ! 質量の補正
    ! Mass fixer
    !
    call MassFixerInit


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  FilterParameterRA = %f', d = (/ FilterParameterRA /) )
    call MessageNotify( 'M', module_name, '  FilterParameterW  = %f', d = (/ FilterParameterW  /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    timefilter_williams2009_inited = .true.

  end subroutine TimeFilterWilliams2009Init

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

end module timefilter_williams2009