!= Module HeatFlux
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: heatflux.f90,v 1.1 2009-03-05 05:39:42 yamasita Exp $
! Tag Name::  $Name: arare4-20090306 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
! 火星湿潤対流計算用. 
!
! 下部境界からのフラックスによる温度と凝結成分混合比の変化率をバルク法に
! 基づいて計算するモジュール
!
! 熱フラックスを Fh と凝結物質のフラックス Fq とすると, 温度と凝結物質
! の変化率 H, Q は
!
!   H = Fh/Δz_1  
!   Q = Fq/Δz_1  
!
! と表される. ここで Δz_1 は最下層の格子間隔である.  
!
! 熱フラックス Fh と凝結物質のフラックス Fq は以下の式にしたがって
! 計算する.
!
!   Fh = - Cdρ|V| * π_1θ_1 
!   Fh = - Cdρ|V| * Q_1 
!
! ここで θ_1, Q_1, π_1 は最下層の温度と凝結物質の混合比および無次元
! 圧力関数, T_sfc は下部境界の温度, Q*(T_sfc) は T_sfc で決まる飽和混
! 合比である. 
!
! ここではモデル最下層の温度を初期条件に戻すようにフラックスを与える. 
!
! 風速値 v とバルク係数 Cd は一定とする. 無次元圧力関数 π_1 は基本場
! の値を用いる.
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module HeatFlux
  !
  !下部境界でのフラックスの計算モジュール
  !
  
  !モジュール読み込み
  use gridset,  only: DimXMin,         & !x 方向の配列の下限
    &                 DimXMax,         & !x 方向の配列の上限
    &                 DimZMin,         & !z 方向の配列の下限
    &                 DimZMax,         & !z 方向の配列の上限
    &                 RegZMin,         & !z 方向の物理領域の下限
    &                 SpcNum,          & !
    &                 DelZ               !z 方向の格子点間隔
  use basicset, only: xz_ExnerBasicZ     !エクスナー関数の基本場
  use StorePotTemp,  only: StorePotTempFlux
  use StoreMixRt,    only: StoreMixRtFlux

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public xz_HeatFluxBulk
  public xz_MixRtFluxBulk
  public xz_HeatFluxDiff
  public xz_MixRtFluxDiff
  public xza_MixRtFluxDiff

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

  !値を保存
  save Bulk
  save Kappa

contains


!!!------------------------------------------------------------------------!!!
  function xz_HeatFluxBulk( xz_PotTemp )
    !
    ! 下部境界からのフラックスによる温度の変化率をバルク法に
    ! 基づいて計算する.  
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8)               :: xz_HeatFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    real(8)  :: VelX = 100.0d0             !下層での水平速度を決めうち

    !初期化
    !  * 全ての値をゼロに固定
    xz_HeatFluxBulk = 0.0d0
    
    !地表面熱フラックスを計算
    !  * 単位は K/s
    !  * 格子点 xz では, 物理領域の最下端の添え字は RegZMin+1
    !
    xz_HeatFluxBulk(:,RegZMin+1) =                        &
      & max(                                              &
      &         0.0d0,                                    &
      &       - Bulk * VelX * xz_ExnerBasicZ(:,RegZMin+1) &
      &          * xz_PotTemp(:,RegZMin+1) / DelZ         &
      &    )
    
  end function xz_HeatFluxBulk
  
!!!------------------------------------------------------------------------!!!
  function xz_MixRtFluxBulk( xz_MixRt )
    !
    ! 下部境界からのフラックスによる凝結成分混合比の変化率をバルク法に
    ! 基づいて計算する.  
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xz_MixRt(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8)               :: xz_MixRtFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面熱フラックス
    real(8)  :: VelX = 100.0d0             !下層での水平速度を決めうち

    !初期化
    !  * 全ての値をゼロに固定
    xz_MixRtFluxBulk = 0.0d0

    
    xz_MixRtFluxBulk(:,RegZMin+1) =                        &
      & max(                                               &         
      &         0.0d0,                                     &
      &       - Bulk * VelX * xz_MixRt(:,RegZMin+1) / DelZ &
      &    )
    
  end function xz_MixRtFluxBulk


!!!------------------------------------------------------------------------!!!
  function xz_HeatFluxDiff( xz_PotTemp )

    implicit none

    real(8), intent(in)   :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8)               :: xz_HeatFluxDiff(DimXMin:DimXMax,DimZMin:DimZMax)

    xz_HeatFluxDiff = 0.0d0

    xz_HeatFluxDiff(:,RegZMin+1) =                                        &
      &  - Kappa * xz_PotTemp(:,RegZMin+1) * xz_ExnerBasicZ(:, RegZMin+1) &
      &    / ( ( DelZ * 5.0d-1 ) ** 2.0d0 )

    call StorePotTempFlux( xz_HeatFluxDiff )

  end function xz_HeatFluxDiff


!!!------------------------------------------------------------------------!!!
  function xz_MixRtFluxDiff( xz_MixRt )

    implicit none

    real(8), intent(in)   :: xz_MixRt(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8)               :: xz_MixRtFluxDiff(DimXMin:DimXMax,DimZMin:DimZMax)


    xz_MixRtFluxDiff = 0.0d0

    xz_MixRtFluxDiff(:,RegZMin+1) =                 &
      & max(                                        &
      &       0.0d0,                                &
      &     - Kappa * xz_MixRt(:,RegZMin+1)         &
      &        / ( ( DelZ * 5.0d-1 ) ** 2.0d0 )     &
      &    )

!    xz_MixRtFluxDiff(:,RegZMin+1) =               &
!      &  - Kappa * xz_MixRt(:,RegZMin+1)          &
!      &    / ( ( DelZ * 5.0d-1 ) ** 2.0d0 )       

  end function xz_MixRtFluxDiff


!!!------------------------------------------------------------------------!!!
  function xza_MixRtFluxDiff( xza_MixRt )

    implicit none

    real(8), intent(in)   :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)
    real(8)               :: xza_MixRtFluxDiff(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)

    xza_MixRtFluxDiff = 0.0d0

    xza_MixRtFluxDiff(:,RegZMin+1,:) =         &
      &  - Kappa * xza_MixRt(:,RegZMin+1,:)    &
      &    / ( ( DelZ * 5.0d-1 ) ** 2.0d0 )       

!    xza_MixRtFluxDiff(:,RegZMin+1,:) =                 &
!      & max(                                           &
!      &      0.0d0,                                    &
!      &    - Kappa * xza_MixRt(:,RegZMin+1,:)          &
!      &      / ( ( DelZ * 5.0d-1 ) ** 2.0d0 )          &
!      &  )

    call StoreMixRtFlux( xza_MixRtFluxDiff )

  end function xza_MixRtFluxDiff


end module HeatFlux



!= Module HeatFlux_N1994
!
! Authors::   ODAKA Masatsugu 
! Version::   $Id: heatflux.f90,v 1.1 2009-03-05 05:39:42 yamasita Exp $
! Tag Name::  $Name: arare4-20090306 $
! 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_N1994
  !
  !下部境界でのフラックスの計算モジュール
  !
  
  !モジュール読み込み
  use gridset,  only: DimXMin,         & !x 方向の配列の下限
    &                 DimXMax,         & !x 方向の配列の上限
    &                 DimZMin,         & !z 方向の配列の下限
    &                 DimZMax,         & !z 方向の配列の上限
    &                 RegZMin,         & !z 方向の物理領域の下限
    &                 DelZ               !z 方向の格子点間隔
  use basicset, only: xz_ExnerBasicZ  ,& !エクスナー関数の基本場
    &                 xz_PotTempBasicZ,& !温位の基本場
    &                 TempSfc         ,& !地表面温度(実際は最下層の温度)
    &                 PressSfc        ,& !地表面圧力 [Pa]
    &                 CpDry           ,& !乾燥成分の定圧比熱 [J/K kg]
    &                 GasRDry         ,& !乾燥成分の気体定数 [J/K kg]
    &                 TempGnd            !地表面温度(地表面フラックス計算用)

  use average, only:  xz_avr_pz          
  use chemcalc, only: SvapPress

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public xz_HeatFluxBulk
  public xz_MixRtFluxBulk

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

  !値を保存
  save Bulk

contains


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

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

    real(8), allocatable  :: xz_VelX(:,:)  !水平風速 (xz 格子)


    !作業配列の割り付け
    allocate(xz_VelX(DimXMin:DimXMax,DimZMin:DimZMax))


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

    kz = RegZMin+1
    xz_VelX = xz_avr_pz(pz_VelX)

    xz_HeatFluxBulk(:,kz) =                                       &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk                                                &
      &       * SQRT(xz_VelX(:,kz)**2.0 + VelX**2.0)              &
      &       * (xz_ExnerBasicZ(:,kz)                             &
      &            * (xz_PotTemp(:,kz)  + xz_PotTempBasicZ(:,kz)) &
      &          - TempSfc                                        &
      &         )/DelZ                                            &
      &    )

    !作業配列の解放
    deallocate(xz_VelX)
    
  end function xz_HeatFluxBulk
  
!!!------------------------------------------------------------------------!!!
  function xz_MixRtFluxBulk( xz_MixRt, pz_VelX, xz_Exner, ID )
    !
    ! 下部境界からのフラックスによる凝結成分混合比の変化率を,
    ! バルク方法に基づいて計算する.
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)   :: xz_MixRt(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !温位の擾乱成分    
    real(8), intent(in)   :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !水平風速
    real(8), intent(in)   :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !圧力の擾乱成分    
    integer, intent(in)   :: ID            !凝結成分の ID 番号
    real(8)               :: xz_MixRtFluxBulk(DimXMin:DimXMax,DimZMin:DimZMax)
                                           !地表面混合比フラックス
    real(8)               :: VelX = 3.0d0  !下層での水平速度嵩上げ値
    integer               :: kz            !配列添字

    real(8), allocatable  :: xz_VelX(:,:)  !水平風速 (xz 格子)


    !作業配列の割り付け
    allocate(xz_VelX(DimXMin:DimXMax,DimZMin:DimZMax))


    !初期化
    !  * 全ての値をゼロに固定
    xz_MixRtFluxBulk = 0.0d0

    kz = RegZMin+1
    xz_VelX = xz_avr_pz(pz_VelX)
    
    xz_MixRtFluxBulk(:,kz) =                                     &
      & max(                                                     &         
      &       0.0d0,                                             &
      &     - Bulk                                               &
      &       * SQRT(xz_VelX(:,kz)**2.0 + VelX**2.0)             &
      &       * (xz_MixRt(:,kz)                                  & 
      &            - SvapPress( ID, TempSfc ) / (                &
      &              ( xz_Exner(:,kz) + xz_ExnerBasicZ(:,kz) )   &
      &              **( CpDry / GasRDry )*PressSfc )            &
      &         )/DelZ                                           &
      &    )

    !作業配列の解放
    deallocate(xz_VelX)

    
  end function xz_MixRtFluxBulk
  
end module HeatFlux_N1994
