!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2006. All rights reserved.
!---------------------------------------------------------------------

module phys_vfilter_conserve_mod
  !
  != 物理過程 鉛直フィルター(保存バージョン)モジュール
  !
  !== 概要
  !
  ! 鉛直方向にフィルターけける.
  ! $\sum T_{k} \Delta p_{k}$
  ! が保存するように補正を行う.
  !
  use type_mod,    only : REKIND, DBKIND, INTKIND, TOKEN, STRING

  implicit none

  private
  public :: phys_vfilter_conserve
contains

  subroutine phys_vfilter_conserve( & 
    & xyz_Temp   , & !(inout)
    & xyz_VelLon , & !(inout)
    & xyz_VelLat , & !(inout)
    & xyz_DTempDtVertFiltCons   , & !(out)
    & xyz_DVelLonDtVertFiltCons , & !(out) 
    & xyz_DVelLatDtVertFilstCons, & !(out)
    & xyr_Temp   , & !(in) 
    & xyr_Press  , & !(in)
    & DelTimePhy )   !(in)
    !
    != 鉛直フィルター 3 点トリオで保存させる版
    !
    !== 概要
    !
    ! 鉛直方向のギザギザを消去するため鉛直方向にフィルターを
    ! かける.
    ! 温度場に対してフィルターをかけることが目的であるが
    ! 東西風, 南北風にフィルターをかけることも実装されている.
    ! 全層を調節する.
    ! 下から順に 3 点トリオで調節していき, 誤差はその 3 点にばらまく.
    !
    !== HISTORY
    ! * 2006-12-1  M.Ishiwatari AGCM5 石渡バージョンの p2grstU.F より移植
    !
    use type_mod,      only: REKIND, DBKIND, INTKIND, TOKEN, STRING
    use grid_3d_mod,   only: im, jm, km

    real(DBKIND), intent(inout) :: xyz_Temp   ( im*jm, km )   ! 温度Ｔ
    real(DBKIND), intent(inout) :: xyz_VelLon   ( im*jm, km ) ! Ｕ
    real(DBKIND), intent(inout) :: xyz_VelLat   ( im*jm, km ) ! Ｖ

    real(DBKIND), intent(in) :: xyr_Press ( im*jm, km+1 ) ! 気圧
    real(DBKIND), intent(in) :: xyr_Temp  ( im*jm, km+1 ) ! 温度
    real(DBKIND), intent(in) :: DelTimePhy                ! 2 Δt

    real(DBKIND), intent(out) :: xyz_DTempDtVertFiltCons (im*jm,km) ! 温度変化率(鉛直フィルター)
    real(DBKIND), intent(out) :: xyz_DVelLonDtVertFiltCons (im*jm,km) ! Ｕ変化率(鉛直フィルター)
    real(DBKIND), intent(out) :: xyz_DVelLatDtVertFilstCons   ( im*jm, km )   ! Ｖ変化率(鉛直フィルター)

    real(DBKIND) :: FGS                        ! フィルター係数 (T)
    real(DBKIND) :: FGSU                       ! 調節する度合(u,v)
    integer(INTKIND) :: KminVfil               ! 調節開始レベル
    integer(INTKIND) :: KmaxVfil               ! 調節終了レベル

    DATA       FGS    / 0.1 /
    DATA       FGSU   / 0.1 /
    DATA       KminVfil / 2   /
    ! 注: KmaxVfil のデフォルト値は下で決める.

    ! とりあえず NAMELIST によるパラメータ入力の実装はサボる
    !  NAMELIST /NMGRST/ FGS, FGSU, KminVfil, KmaxVfil
    !  LOGICAL    OFIRST
    !  DATA       OFIRST / .TRUE. /
    !  SAVE

    real(DBKIND) :: xyz_DelPress  ( im*jm, km )   !" ΔＰ
    real(DBKIND) :: xyz_TempBasic    ( im*jm, km )   !" T 基本場
    real(DBKIND) :: xyz_VelLonBasic    ( im*jm, km )   !" T 基本場
    real(DBKIND) :: xyz_VelLatBasic    ( im*jm, km )   !" T 基本場
    real(DBKIND) :: xyz_TempOrg  ( im*jm, km )   !" Ｔ(調節前)
    real(DBKIND) :: xyz_VelLonOrg  ( im*jm, km )   !" Ｕ(調節前)
    real(DBKIND) :: xyz_VelLatOrg  ( im*jm, km )   !" Ｖ(調節前)
    real(DBKIND) :: SumTempDelPress
    real(DBKIND) :: SumVelLonDelPress
    real(DBKIND) :: SumVelLatDelPress
    real(DBKIND) :: SumDelPress
    character(STRING),  parameter:: subname = "phys_vfilter_conserve"
    integer(INTKIND) :: ij
    integer(INTKIND) :: k

    ! 開始処理
    call BeginSub(subname)

    ! KmaxVfil デフォルト値
    KmaxVfil = km-1

    ! 本来ならここで, NAMELIST によるパラメータ入力
    ! をおこなうべき.

    ! 1. ファクター計算 >
    do k = 1, km
      do ij = 1, im*jm
        xyz_DelPress  (ij,k) = xyr_Press(ij,k) - xyr_Press(ij,k+1)
        xyz_TempBasic (ij,k) = ( xyr_Temp(ij,k) + xyr_Temp(ij,k+1) )/2.0d0
      enddo  
    enddo

    do ij = 1, im*jm
      xyz_VelLonBasic (ij,1) &
        & = ( 3.0d0 * xyz_VelLon(ij,1) + xyz_VelLon(ij,2) )/4.0d0
      xyz_VelLatBasic (ij,1) &
        & = ( 3.0d0 * xyz_VelLat(ij,1) + xyz_VelLat(ij,2) )/4.0d0
      do k = 2, km-1
        xyz_VelLonBasic (ij,k) &
          & = (   xyz_VelLon(ij,k-1) + 2.0d0 * xyz_VelLon(ij,k) &
          &     + xyz_VelLon(ij,k+1) )/4.0d0
        xyz_VelLatBasic (ij,k) &
          & = (   xyz_VelLat(ij,k-1) + 2.0 * xyz_VelLat(ij,k) &
          &     + xyz_VelLat(ij,k+1) )/4.0d0
      enddo
      xyz_VelLonBasic (ij,km) &
        & = ( xyz_VelLon(ij,km-1) + 3.0d0 * xyz_VelLon(ij,km) )/4.0d0
      xyz_VelLatBasic (ij,km) &
        & = ( xyz_VelLat(ij,km-1) + 3.0d0 * xyz_VelLat(ij,km) )/4.0d0
    enddo

    ! 調節前の値の保存
    xyz_TempOrg = xyz_Temp
    xyz_VelLonOrg = xyz_VelLon
    xyz_VelLatOrg = xyz_VelLat

    !  2. 調節
    xyz_DTempDtVertFiltCons = 0.0d0
    xyz_DVelLonDtVertFiltCons = 0.0d0
    xyz_DVelLatDtVertFilstCons = 0.0d0

    SumTempDelPress=0.0
    SumDelPress=0.0

    do k = KminVfil, KmaxVfil
      do ij = 1, im*jm
        SumTempDelPress = (xyz_TempBasic(ij,k-1) - xyz_Temp(ij,k-1) ) &
          &      * xyz_DelPress(ij,k-1) &
          & + (xyz_TempBasic(ij,k) - xyz_Temp(ij,k) ) * xyz_DelPress(ij,k) &
          & + (xyz_TempBasic(ij,k+1) - xyz_Temp(ij,k+1) ) & 
          &   * xyz_DelPress(ij,k+1)
        SumVelLonDelPress = (xyz_VelLonBasic(ij,k-1) - xyz_VelLon(ij,k-1) ) &
          &   * xyz_DelPress(ij,k-1) &
          & + (xyz_VelLonBasic(ij,k) - xyz_VelLon(ij,k) ) &
          &   * xyz_DelPress(ij,k) &
          & + (xyz_VelLonBasic(ij,k+1) - xyz_VelLon(ij,k+1) ) &
          &   * xyz_DelPress(ij,k+1)
        SumVelLatDelPress = (xyz_VelLatBasic(ij,k-1) - xyz_VelLat(ij,k-1) ) &
          &   * xyz_DelPress(ij,k-1) &
          & + (xyz_VelLatBasic(ij,k) - xyz_VelLat(ij,k) ) &
          &   * xyz_DelPress(ij,k) &
          & + (xyz_VelLatBasic(ij,k+1) - xyz_VelLat(ij,k+1) ) &
          &   * xyz_DelPress(ij,k+1)
        SumDelPress =  xyz_DelPress(ij,k-1) + xyz_DelPress(ij,k) &
          & + xyz_DelPress(IJ,K+1)

        xyz_Temp(ij,k-1) = xyz_Temp(ij,k-1) &
          & + FGS*(   xyz_TempBasic(ij,k-1) - xyz_Temp(ij,k-1) &
          &         - SumTempDelPress/SumDelPress )
        xyz_VelLon(ij,k-1) = xyz_VelLon(ij,k-1) &
          & + FGSU*(   xyz_VelLonBasic(ij,k-1) - xyz_VelLon(ij,k-1) &
          &          - SumVelLonDelPress/SumDelPress )
        xyz_VelLat(ij,k-1) = xyz_VelLat(ij,k-1) &
          & + FGSU*(   xyz_VelLatBasic(ij,k-1) - xyz_VelLat(ij,k-1) &
          &          - SumVelLatDelPress/SumDelPress )

        xyz_Temp(ij,k) = xyz_Temp(ij,k) &
          & + FGS*(   xyz_TempBasic(ij,k) - xyz_Temp(ij,k) &
          &         - SumTempDelPress/SumDelPress )
        xyz_VelLon(ij,k) = xyz_VelLon(ij,k) &
          & + FGSU*(   xyz_VelLonBasic(ij,k) - xyz_VelLon(ij,k) &
          &          - SumVelLonDelPress/SumDelPress )
        xyz_VelLat(ij,k) = xyz_VelLat(ij,k) &
          & + FGSU*(   xyz_VelLatBasic(ij,k) - xyz_VelLat(ij,k) &
          &          - SumVelLatDelPress/SumDelPress )

        xyz_Temp(ij,k+1) = xyz_Temp(ij,k+1) &
          & + FGS*(    xyz_TempBasic(ij,k+1) - xyz_Temp(ij,k+1) &
          &          - SumTempDelPress/SumDelPress )
        xyz_VelLon(ij,k+1) = xyz_VelLon(ij,k+1) &
          & + FGSU*(   xyz_VelLonBasic(ij,k+1) - xyz_VelLon(ij,k+1) &
          &          - SumVelLonDelPress/SumDelPress )
        xyz_VelLat(ij,k+1) = xyz_VelLat(ij,k+1) &
          & + FGSU*(   xyz_VelLatBasic(ij,k+1) - xyz_VelLat(ij,k+1) &
          &          - SumVelLatDelPress/SumDelPress )
      enddo
    enddo

    do k =1,km
      do ij = 1,im*jm
        xyz_DTempDtVertFiltCons(ij,k) = ( xyz_Temp(ij,k) - xyz_TempOrg(ij,k) )/DelTimePhy
        xyz_DVelLonDtVertFiltCons(ij,k) = ( xyz_VelLon(ij,k) - xyz_VelLonOrg(ij,k) )/DelTimePhy
        xyz_DVelLatDtVertFilstCons(ij,k) = ( xyz_VelLat(ij,k) - xyz_VelLatOrg(ij,k) )/DelTimePhy
      enddo
    enddo

    ! 終了処理
    call EndSub(subname)

  end subroutine phys_vfilter_conserve
end module phys_vfilter_conserve_mod

