!= タイムフィルター (Asselin, 1972)
!
!= Time filter (Asselin, 1972)
!
! Authors::   Yasuhiro Morikawa, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: timefilter_asselin1972.f90,v 1.13 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_asselin1972
  !
  != タイムフィルター (Asselin, 1972)
  !
  != Time filter (Asselin, 1972)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 計算モードを抑えるためのタイムフィルター (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
  ! 
  ! TimeFilter    :: タイムフィルターの適用
  ! ------------  :: ------------
  ! TimeFilter    :: Apply time filter
  !
  !== NAMELIST
  !
  ! NAMELIST#timefilter_asselin1972_nml
  !
  !== References
  !
  ! * 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 :: TimeFilter
  public :: TimeFilterSurfVars
  public :: TimeFiltInit

  ! 公開変数
  ! Public variables
  !

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

  real(DP), save:: FilterParameter
                              ! タイムフィルターの係数. 
                              ! Time filter coefficient
!!$  integer, save:: IntStep     ! タイムフィルターを適用するステップ間隔. 
!!$                              ! Step interval of time filter application
!!$  integer, save:: CurStep     ! TimeFilter が呼ばれた回数. 
!!$                              ! Number of times that "TimeFilter" is called
!!$  integer, save:: PrevStep    ! 前回タイムフィルターを適用したステップ数. 
!!$                              ! Step number that is filtered previously

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

contains

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

  subroutine TimeFilter( &
    & 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  &   ! (in)
    & )
    !
    ! 予報変数にタイムフィルターを適用します. 
    !
    ! Time filter is applied to predictional variables.
    !

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

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

    ! 宣言文 ; 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(in):: xyz_UA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ u (t+\Delta t) $ .   東西風速. Eastward wind
    real(DP), intent(in):: xyz_VA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ v (t+\Delta t) $ .   南北風速. Northward wind
    real(DP), intent(in):: xyz_TempA (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T (t+\Delta t) $ .   温度. Temperature
    real(DP), intent(in):: xyzf_QMixA (0:imax-1, 1:jmax, 1:kmax, 1:ncmax)
                              ! $ q (t+\Delta t) $ .   混合比. Mass mixing ratio of constituents
    real(DP), intent(in):: xy_PsA (0:imax-1, 1:jmax)
                              ! $ p_s (t+\Delta t) $ . 地表面気圧. Surface pressure

    ! 作業変数
    ! Work variables
    !

    ! 実行文 ; Executable statement
    !

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

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


!!$    ! ステップのチェック
!!$    ! Check Step
!!$    !
!!$    CurStep = CurStep + 1
!!$    if ( CurStep - PrevStep < IntStep ) return
!!$    PrevStep = CurStep

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

    ! タイムフィルターを適用
    ! Apply time filter
    !
    xyz_UN     =  xyz_UN    &
      & + FilterParameter * ( xyz_UB     - 2.0d0 * xyz_UN     + xyz_UA     )

    xyz_VN     =  xyz_VN    &
      & + FilterParameter * ( xyz_VB     - 2.0d0 * xyz_VN     + xyz_VA     )

    xyz_TempN  =  xyz_TempN &
      & + FilterParameter * ( xyz_TempB  - 2.0d0 * xyz_TempN  + xyz_TempA  )

    xy_PsN     =  xy_PsN    &
      & + FilterParameter * ( xy_PsB     - 2.0d0 * xy_PsN     + xy_PsA     )

    xyzf_QMixN = xyzf_QMixN &
      & + FilterParameter * ( xyzf_QMixB - 2.0d0 * xyzf_QMixN + xyzf_QMixA )


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

  end subroutine TimeFilter

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

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

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

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

    ! 宣言文 ; 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(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(in   ):: xy_SurfMajCompIceA(0:imax-1, 1:jmax)
                              ! $ M_mcs (t+\Delta t) $ . 
                              ! Surface major component ice amount
    real(DP), intent(in   ):: xy_SoilMoistA (0:imax-1, 1:jmax)
                              ! $ M_ws (t-\Delta t) $ . 地表水質量. 
                              ! Soil moisture
    real(DP), intent(in   ):: xy_SurfSnowA  (0:imax-1, 1:jmax)
                              ! $ M_ss (t+\Delta t) $ . 積雪量.
                              ! Surface snow amount

    ! 作業変数
    ! Work variables
    !
!!$    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_asselin1972_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


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


!!$    ! ステップのチェック
!!$    ! Check Step
!!$    !
!!$    CurStep = CurStep + 1
!!$    if ( CurStep - PrevStep < IntStep ) return
!!$    PrevStep = CurStep

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

    ! タイムフィルターを適用
    ! Apply time filter
    !
    xy_SurfMajCompIceN = xy_SurfMajCompIceN                         &
      & + FilterParameter                                           &
      & * ( xy_SurfMajCompIceB - 2.0_DP * xy_SurfMajCompIceN + xy_SurfMajCompIceA )

    xy_SoilMoistN = xy_SoilMoistN                                   &
      & + FilterParameter                                           &
      & * ( xy_SoilMoistB - 2.0_DP * xy_SoilMoistN + xy_SoilMoistA )

    xy_SurfSnowN  = xy_SurfSnowN                                    &
      & + FilterParameter                                           &
      & * ( xy_SurfSnowB  - 2.0_DP * xy_SurfSnowN  + xy_SurfSnowA  )

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

  end subroutine TimeFilterSurfVars

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

  subroutine TimeFiltInit
    !
    ! timefilter_asselin1972 モジュールの初期化を行います. 
    ! NAMELIST#timefilter_asselin1972_nml の読み込みはこの手続きで行われます. 
    !
    ! "timefilter_asselin1972" module is initialized. 
    ! "NAMELIST#timefilter_asselin1972_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

    ! メッセージ制御
    ! 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_asselin1972_nml/ &
      & FilterParameter   !, IntStep
          !
          ! デフォルト値については初期化手続 "timefilter_asselin1972#Hs94Init" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "timefilter_asselin1972#Hs94Init" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( timefilter_asselin1972_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FilterParameter = 0.05_DP
!!$    IntStep         = 1
!!$    CurStep         = 0

    ! 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_asselin1972_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_asselin1972_nml )
    end if

!!$    ! IntStep のチェック
!!$    ! Check "IntStep"
!!$    !
!!$    if ( IntStep < 1 ) then
!!$      call MessageNotify( 'E', module_name, &
!!$        & 'IntStep=<%d> must be greater than 0', &
!!$        & i = (/ IntStep /) )
!!$    end if
!!$
!!$    ! PrevStep の設定
!!$    ! Configure "PrevStep"
!!$    !
!!$    PrevStep = - IntStep

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

    timefilter_asselin1972_inited = .true.

  end subroutine TimeFiltInit

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

end module timefilter_asselin1972