!= Module HeatFlux_3d
!
! Authors::   ODAKA Masatsugu 
! Version::   $Id: heatflux_3d.f90,v 1.2 2007/10/01 08:57:59 odakker Exp $
! Tag Name::  $Name: arare4-20071012 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
! 下部境界からのフラックスによる温度と凝結成分の変化率をバルク方法に
! 基づいて計算するモジュール. これは中島 (1994) で用いられた方法である.
!
! 熱フラックスを Fh と凝結物質のフラックス Fq とすると, 温度と凝結物質
! の変化率 H, Q は
!
!   H = Fh/Δz_1  
!   Q = Fq/Δz_1  
!
! と表される. ここで Δz_1 は最下層の格子間隔である.  
!
! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって
! 計算する.
!
!   Fh = - Cdρ|V| * (π_1θ_1 - T_sfc)
!   Fh = - Cdρ|V| * (Q_1 - Q*(T_sfc))
!
! ここで θ_1, Q_1, π_1 は最下層の温度と凝結物質の混合比および無次元
! 圧力関数, T_sfc は下部境界の温度, Q*(T_sfc) は T_sfc で決まる飽和混
! 合比である. 
!
! バルク係数 Cd は一定とする. 無次元圧力関数 π_1 は基本場の値を用いる.
! 風速値 |V| は
!
!   V = ( V^2 + V_0^2 )^(1/2)
!
! と計算する. 
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!

module HeatFlux_3d
  !
  !下部境界でのフラックスの計算モジュール
  !
  
  !モジュール読み込み
  use gridset_3d,  only: DimXMin,         & !x 方向の配列の下限
    &                 DimXMax,         & !x 方向の配列の上限
    &                 DimYMin,         & !y 方向の配列の下限
    &                 DimYMax,         & !y 方向の配列の上限
    &                 DimZMin,         & !z 方向の配列の下限
    &                 DimZMax,         & !z 方向の配列の上限
    &                 RegZMin,         & !z 方向の物理領域の下限
    &                 z_dz               !z 方向の格子点間隔
  use basicset_3d, only: xyz_ExnerBasicZ  ,& !エクスナー関数の基本場
    &                 xyz_PotTempBasicZ,& !温位の基本場
    &                 TempSfc         ,& !地表面温度
    &                 PressSfc        ,& !地表面圧力 [Pa]
    &                 CpDry           ,& !乾燥成分の定圧比熱 [J/K kg]
    &                 GasRDry            !乾燥成分の気体定数 [J/K kg]

  use xyz_module, only:  xyz_avr_pyz, xyz_avr_xqz
  use chemcalc, only: SvapPress
  use StorePotTemp_3d,  only: StorePotTempFlux


  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public xyz_HeatFluxBulk
  public pyz_MomFluxBulk
  public xqz_MomFluxBulk

  !変数定義
!  real(8)  :: Bulk = 1.5d-3        !熱・運動量フラックスのバルク係数
  real(8)  :: Bulk = 1.0d-2        !熱・運動量フラックスのバルク係数

  !値を保存
  save Bulk

  interface pyz_MomFluxBulk
    module procedure aaz_MomFluxBulk
  end interface

  interface xqz_MomFluxBulk
    module procedure aaz_MomFluxBulk
  end interface


contains


!!!------------------------------------------------------------------------!!!
  function aaz_MomFluxBulk( Vel )
    ! 
    ! 下部境界からのフラックスによる運動量の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: Vel(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !水平風速
    real(8)               :: aaz_MomFluxBulk(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    integer               :: kz            !配列添字


    !初期化
    !  * 全ての値をゼロに固定
    aaz_MomFluxBulk = 0.0d0
    
    !地表面運動量フラックスによる変化率を計算
    !  * 単位は m/s^2
    !  * 格子点 xz では, 物理領域の最下端の添え字は RegZMin+1

    kz = RegZMin

    aaz_MomFluxBulk(:,:,kz) =                                     &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk * abs(Vel(:,:,kz)) * Vel(:,:,kz) /z_dz(kz) )

  end function aaz_MomFluxBulk
  
!!!------------------------------------------------------------------------!!!
  function xyz_HeatFluxBulk( xyz_PotTemp, pyz_VelX, xqz_VelY )
    ! 
    ! 下部境界からのフラックスによる温度の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8), intent(in)   :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !水平風速
    real(8), intent(in)   :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !水平風速
    real(8)               :: xyz_HeatFluxBulk(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    real(8)               :: VelX = 0.0d0  !下層での水平速度嵩上げ値
    real(8)               :: TempGrndSfc = 270.0d0  ! 
    integer               :: kz            !配列添字

    real(8), allocatable  :: xyz_VelX(:,:,:)  !水平風速 (xyz 格子)
    real(8), allocatable  :: xyz_VelY(:,:,:)  !水平風速 (xyz 格子)


    !作業配列の割り付け
    allocate(xyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax))
    allocate(xyz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax))


    !初期化
    !  * 全ての値をゼロに固定
    xyz_HeatFluxBulk = 0.0d0
    
    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は RegZMin+1

    kz = RegZMin
    xyz_VelX = xyz_avr_pyz(pyz_VelX)
    xyz_VelY = xyz_avr_xqz(xqz_VelY)

    xyz_HeatFluxBulk(:,:,kz) =                                    &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk                                                &
      &       * SQRT( (xyz_VelX(:,:,kz) + VelX)**2.0              &
      &               + xyz_VelY(:,:,kz)**2.0 )                   &
      &       * (xyz_ExnerBasicZ(:,:,kz)                          &
      &            * (xyz_PotTemp(:,:,kz)  + xyz_PotTempBasicZ(:,:,kz)) &
      &          - TempGrndSfc                                    &
      &         )/z_dz(kz)                                        &
      &    )

    call StorePotTempFlux( xyz_HeatFluxBulk )

    !作業配列の解放
    deallocate(xyz_VelX)
    deallocate(xyz_VelY)
   
  end function xyz_HeatFluxBulk

end module HeatFlux_3d