!= 鉛直フィルター (Ishiwatari et al., 2002)
!
!= Vertical filter (Ishiwatari et al., 2002)
!
! Authors::   Satoshi Noda
! Version::   $Id: vertical_filter_wt_VF-UVT.f90,v 1.2 2010/05/08 06:53:35 noda Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2010. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module vertical_filter
  !
  != 鉛直フィルター (Ishiwatari et al., 2002)
  !
  != Vertical filter (Ishiwatari et al., 2002)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 鉛直モードを抑えるための鉛直フィルター (Ishiwatari et al., 2002)
  ! を適用するためのモジュールです.
  !
  ! すなわち
  ! \begin{align}
  !   T_{B,k} &\equiv \frac{T_{k+1/2} + T_{k-1/2}}{2} \\
  ! \end{align}
  ! と定義した上で, 
  ! \begin{align}
  !   \bar{T}_{k-1} &= T_{k-1} + F_T
  !     ( T_{B, k-1} - T_{k-1} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}, \\
  !   \bar{T}_{k} &= T_{k} + F_T
  !     ( T_{B, k} - T_{k} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}, \\
  !   \bar{T}_{k+1} &= T_{k+1} + F_T
  !     ( T_{B, k+1} - T_{k+1} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}
  ! \end{align}
  ! という操作を下の層から順に適用することで $ \bar{T} $ を求めます. 
  ! ここで添え字 $ k $ は層の番号です.
  ! 鉛直フィルターが適用された項には上線 ( $ \bar{\quad} $ )
  ! が付いています.
  !
  ! 鉛直フィルターの係数 $F_T$ は Create の際に設定します.
  !
  !++
  !
  ! Vertical filter (Ishiwatari et al., 2002) for suppression of vertical mode 
  ! is applied.
  !
  ! Concretely, $ \bar{A}^{t} $ is derived as follows.
  ! Reference Temperature is defined as
  ! \begin{align}
  !   T_{B,k} &\equiv \frac{T_{k+1/2} + T_{k-1/2}}{2}.
  ! \end{align}
  ! Operations as follows are applied from lower layer to higher layer.
  ! \begin{align}
  !   \bar{T}_{k-1} &= T_{k-1} + F_T
  !     ( T_{B, k-1} - T_{k-1} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}, \\
  !   \bar{T}_{k} &= T_{k} + F_T
  !     ( T_{B, k} - T_{k} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}, \\
  !   \bar{T}_{k+1} &= T_{k+1} + F_T
  !     ( T_{B, k+1} - T_{k+1} )
  !     - \frac{ \Sigma_{k'=k-1}^{k+1} F_T (T_{B,k'} - T_{k'}) \delta p_{k'} }{\Sigma_{k'=k-1}^{k+1}  \delta p_{k'}}
  ! \end{align}
  ! where a suffix $ k $ is the number of layer.
  ! Over-bar ( $ \bar{\quad} $ ) represents the terms
  ! applied vertical filter.
  !
  ! Vertical filter coefficient $ F_T $ is specified by "Create".
  !
  !== Procedures List
  ! 
  ! VerticalFilter    :: 鉛直フィルターの適用
  ! ----------------  :: ---------------------
  ! VerticalFilter    :: Apply vertical filter
  !
  !== NAMELIST
  !
  ! NAMELIST#vertical_filter_nml
  !
  !== References
  !
  ! * Ishiwatari, M., Takehiro S., Nakajima, K. and Hayashi, Y. -Y. 2002: 
  !   A Numerical Study on Appearance of the Runaway Greenhouse State of a
  !   Three-Dimensional Gray Atmosphere.
  !   <i>J. Atmos. Sci.</i>, <b>59</b>, 3223--3238.
  !

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

  use timeset, only: DelTime

  ! 種別型パラメタ
  ! 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:: VerticalFilter

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

  ! 非公開変数
  ! Private variables
  !
  real(DP), save:: FilterParameter
                              ! 鉛直フィルターの係数. 
                              ! Vertical filter coefficient
!!$  integer, save:: IntStep     ! 鉛直フィルターを適用するステップ間隔. 
!!$                              ! Step interval of vertical filter application
!!$  integer, save:: CurStep     ! VerticalFilter が呼ばれた回数. 
!!$                              ! Number of times that "VerticalFilter" is called
!!$  integer, save:: PrevStep    ! 前回鉛直フィルターを適用したステップ数. 
!!$                              ! Step number that is filtered previously


  character(*), parameter:: module_name = 'vertical_filter'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: vertical_filter_wt_VF-UVT.f90,v 1.2 2010/05/08 06:53:35 noda Exp $'
                              ! モジュールのバージョン
                              ! Module version

  ! INTERFACE 文 ; INTERFACE statements
  !
  interface VerticalFilter
    module procedure VerticalFilter
  end interface

contains

  subroutine VerticalFilter( &
    & xyz_UA, xyz_VA, xyz_TempA, &
    & xyr_Press, xyr_Temp &
    & )
    !
    ! 予報変数に鉛直フィルターを適用します. 
    !
    ! Vertical filter is applied to predictional variables.
    !

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    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(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ p_{1/2} (t+\Delta t) $ .   気圧. Pressure
    real(DP), intent(in):: xyr_Temp (0:imax-1, 1:jmax, 0:kmax)
                              ! $ T_{1/2} (t+\Delta t) $ .   温度. Temperature


    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: k1              ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    real(DP):: xyz_TempRef(0:imax-1, 1:jmax, 1:kmax), &
      & Temp_corr(0:imax-1, 1:jmax)

    integer  :: NLDump, NLNewtonZonal, NLRayleigh
    real(DP) :: CoefDump(1:kmax), CoefNewtonZonal(1:kmax), CoefRayleigh(1:kmax)
    real(DP) :: TimeScaleDump


    !
    ! testing (鉛直フィルタ (Ishiwatari et al., 2002))
    !
    ! Ishiwatari et al. で使われていた定数 (F_T = 0.1)
    
    NLDump = 7
    TimeScaleDump = 10800.0_DP
    
    do k= 1, kmax - NLDump
      CoefDump(k) = 0.0_DP
      print *, '*** DEBUG ***', k, CoefDump(k)
    end do
    
    do k= kmax - NLDump + 1, kmax
      CoefDump(k) = DelTime / (TimeScaleDump * ( kmax - k + 1 ))
      print *, '*** DEBUG ***', k, CoefDump(k)
    end do
    
    NLNewtonZonal = NLDump
    CoefNewtonZonal = CoefDump
   
    NLRayleigh = NLDump
    CoefRayleigh = CoefDump


    ! 実行文 ; Executable statement
    !

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

    ! 初期化
    ! Initialization
    !
    if ( .not. vertical_filter_inited ) call VerticalFilterInit

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

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

    ! 鉛直フィルターを適用
    ! Apply vertical filter
    !


!    call VerticalFilterEachVar(Var, FilterParameter)

    ! Temp
    do k=1,kmax
      xyz_TempRef(:,:,k)  = 0.5_DP * ( xyr_Temp(:,:,k-1) + xyr_Temp(:,:,k) )
    end do

    do k=2,kmax-1
      Temp_corr = 0.0_DP
      do k1=-1,1
        Temp_corr = Temp_corr + FilterParameter &
          & * ( xyz_TempRef(:,:,k+k1) - xyz_TempA(:,:,k+k1) ) &
          & * ( xyr_Press(:,:,k+k1-1) - xyr_Press(:,:,k+k1) )
      end do
      Temp_corr = Temp_corr / ( xyr_Press(:,:,k-2) - xyr_Press(:,:,k+1) )
      do k1=-1,1
        xyz_TempA(:,:,k+k1) = xyz_TempA(:,:,k+k1) &
          & + FilterParameter * ( xyz_TempRef(:,:,k+k1) - xyz_TempA(:,:,k+k1) ) - Temp_corr
      end do
    end do


    ! U
    do k=2,kmax-1
      xyz_UA(:,:,k) = xyz_UA(:,:,k) &
        & + 0.5 * FilterParameter   &
        & * ( xyz_UA(:,:,k-1) + xyz_UA(:,:,k+1) - 2 * xyz_UA(:,:,k))
    end do

    ! V
    do k=2,kmax-1
      xyz_VA(:,:,k) = xyz_VA(:,:,k) &
        & + 0.5 * FilterParameter   &
        & * ( xyz_VA(:,:,k-1) + xyz_VA(:,:,k+1) - 2 * xyz_VA(:,:,k))
    end do



    !
    ! testing (ニュートン冷却(帯状方向への平滑化, Ishiwatari et al. 2002))
    !

    xyz_TempRef = 0.0_DP
    ! 参照温度の計算
    do k = kmax-NLNewtonZonal+1, kmax
!      do j=1,jmax
        do i=0,imax-1
          xyz_TempRef(0,:,k) = xyz_TempRef(0,:,k) + xyz_TempA(i,:,k)
        end do
        xyz_TempRef(0,:,k) = xyz_TempRef(0,:,k) / imax
        do i=1,imax-1
          xyz_TempRef(i,:,k) = xyz_TempRef(0,:,k)
        end do
!      end do
    end do

    ! 温度の補正
    do k = kmax-NLNewtonZonal+1, kmax
      xyz_TempA(:,:,k) = xyz_TempA(:,:,k) &
        & - CoefNewtonZonal(k) * ( xyz_TempA(:,:,k) - xyz_TempRef(:,:,k) )
    end do

    !
    ! testing (レイリー摩擦)
    !
    
    do k= kmax-NLRayleigh+1,kmax
      xyz_UA(:,:,k) = (1.0_DP - CoefRayleigh(k)) * xyz_UA(:,:,k)
      xyz_VA(:,:,k) = (1.0_DP - CoefRayleigh(k)) * xyz_VA(:,:,k)
    end do







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

  end subroutine VerticalFilter

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

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

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

    ! 実行文 ; Executable statement
    !

    if ( vertical_filter_inited ) return
    call InitCheck

    ! デフォルト値の設定
    ! Default values settings
    !
    FilterParameter = 0.1_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 = vertical_filter_nml, &  ! (out)
        & iostat = iostat_nml )        ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 ) write( STDOUT, nml = vertical_filter_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) )

    vertical_filter_inited = .true.
  end subroutine VerticalFilterInit

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

  subroutine InitCheck
    !
    ! 依存モジュールの初期化チェック
    !
    ! Check initialization of dependency modules

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

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_inited

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: constants_inited

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: timeset_inited

    ! 実行文 ; Executable statement
    !

    if ( .not. namelist_util_inited ) &
      & call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) &
      & call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) &
      & call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) &
      & call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

    if ( .not. timeset_inited ) &
      & call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

  end subroutine InitCheck

end module vertical_filter
