!= 譎る俣繝輔ぅ繝ｫ繧ｿ繝ｼ (Williams, 2009)
!
!= Time filter (Williams, 2009)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: timefilter_williams2009.f90,v 1.4 2015/01/29 12:09:04 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. 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 2015/01/29 12:09:04 yot 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

    ! 螳｣險譁� ; 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 ) 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