!= Module HeatFlux
!
! Authors::   ODAKA Masatsugu 
! Version::   $Id: surfaceflux_bulk.f90,v 1.7 2011-06-03 06:11:22 sugiyama Exp $
! Tag Name::  $Name:  $
! 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 Surfaceflux_bulk
  !
  !下部境界でのフラックスの計算モジュール
  !
  
  !モジュール読み込み
  use dc_types, only: DP, STRING
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

  use mpi_wrapper,only: myrank
  use gridset,  only: imin,         & !x 方向の配列の下限
    &                 imax,         & !x 方向の配列の上限
    &                 jmin,         & !y 方向の配列の下限
    &                 jmax,         & !y 方向の配列の上限
    &                 kmin,         & !z 方向の配列の下限
    &                 kmax,         & !z 方向の配列の上限
    &                 nx, ny, nz, ncmax
  use axesset, only:  z_dz,         & !z 方向の格子点間隔
    &                 xyz_avr_pyz,  &
    &                 xyz_avr_xqz
  use basicset, only: xyz_ExnerBZ,  & !エクスナー関数の基本場
    &                 xyz_PressBZ,  & !
    &                 xyz_PTempBZ,  & !温位の基本場
    &                 xyz_TempBZ,   & !
    &                 xyzf_QMixBZ     !温位の基本場
  use constants,only: MolWtDry
  use composition,only : IdxCC, IdxCG, SpcWetID, CondNum, MolWtWet, SpcWetSymbol
  use chemcalc,only : SvapPress


  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public surfaceflux_bulk_init
  public surfaceflux_bulk_forcing

  !変数定義
  real(DP), save  :: Bulk = 1.5d-3    !熱・運動量フラックスのバルク係数
  real(DP), save  :: VelX0 = 0.0d0    !下層での水平速度嵩上げ値
  real(DP), save  :: VelY0 = 0.0d0    !下層での水平速度嵩上げ値


contains
!!!------------------------------------------------------------------------!!!
  subroutine Surfaceflux_Bulk_init( cfgfile )
    !
    !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
    !

    !暗黙の型宣言禁止
    implicit none

    !入力変数
    character(STRING), intent(in)      :: cfgfile

    !内部変数
    integer    :: l, unit

    !---------------------------------------------------------------    
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_bulk_nml/ Bulk, VelX0, VelY0
    
    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=surfaceflux_bulk_nml)
    close(unit)  

    if (myrank == 0) then 
      call MessageNotify( "M", "SurfaceSfc", "Bulk = %f", d=(/Bulk/))
      call MessageNotify( "M", "SurfaceSfc", "VelX0 = %f", d=(/VelX0/))
      call MessageNotify( "M", "SurfaceSfc", "VelY0 = %f", d=(/VelY0/))
    end if

    call HistoryAutoAddVariable(  &
      & varname='PTempSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='surface flux of potential temperature', &
      & units='kg.kg-1.s-1',            &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='flux of heating from surface', &
      & units='kg.m.s-1',            &
      & xtype='double')

    call HistoryAutoAddVariable(  &
      & varname='VelYSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='flux of momentum from surface', &
      & units='kg.m.s-1',            &
      & xtype='double')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Sfc', & 
        & dims=(/'x','y','z','t'/),     &
        & longname='Surface Flux term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='double')
    end do

  end subroutine Surfaceflux_Bulk_init


!!!------------------------------------------------------------------------!!!
  subroutine Surfaceflux_Bulk_forcing( &
    &   Time, pyz_VelX, xqz_VelY, xyz_PTemp, xyz_Exner, xyzf_QMix, &
    &   pyz_DVelXDt, xqz_DVelYDt, xyz_DPTempDt, xyzf_DQMixDt &
    & )
    ! 
    ! 下部境界からのフラックスによる温度の変化率を,
    ! バルク方法に基づいて計算する.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)   :: Time         !時刻
    real(DP), intent(in)   :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速
    real(DP), intent(in)   :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !温位の擾乱成分    
    real(DP), intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: pyz_Momflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !運動量フラックス
    real(DP)               :: xqz_Momflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !運動量フラックス
    real(DP)               :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !地表面熱フラックス
    real(DP)               :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !物質的フラックス
    real(DP)               :: pyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xqz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyz_VelXSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_VelYSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: TempSfc 
    real(DP)               :: PressSfc 
    real(DP)               :: xyz_VelSfc(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    integer                :: kz            !配列添字
    integer                :: s             !ループ変数

    ! 初期化
    !
    kz = 1
    xyz_Heatflux = 0.0d0
    xyzf_QMixflux = 0.0d0
    pyz_Momflux = 0.0d0
    xqz_Momflux = 0.0d0

    pyz_DVelXDt0  = pyz_DVelXDt
    xqz_DVelYDt0  = xqz_DVelYDt
    xyz_DPTempDt0 = xyz_DPTempDt
    xyzf_DQMixDt0 = xyzf_DQMixDt

    xyz_TempAll  = (xyz_Exner + xyz_ExnerBZ) * (xyz_PTemp + xyz_PTempBZ)
    xyzf_QMixAll = xyzf_QMix + xyzf_QMixBZ
    
    pyz_VelXSfc = pyz_VelX + VelX0
    xqz_VelYSfc = xqz_VelY + VelY0    

    ! 地表面での値
    ! 
    TempSfc = xyz_TempBZ(1,1,kz)  !水平一様なので, i=j=1 で代表させる
    PressSfc= xyz_PressBZ(1,1,kz) !水平一様なので, i=j=1 で代表させる
    
    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1
    xyz_VelXSfc = xyz_avr_pyz(pyz_VelXSfc)
    xyz_VelYSfc = xyz_avr_xqz(xqz_VelYSfc)    
    
    xyz_VelSfc = SQRT(  xyz_VelXSfc ** 2.0 + xyz_VelYSfc ** 2.0 )

    xyz_Heatflux(:,:,kz) =                                        &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk * xyz_VelSfc(:,:,kz)                           &
      &       * ( xyz_TempAll(:,:,kz) - TempSfc )                 &
      &       / ( z_dz(kz) * 5.0d-1 )                             &
      &    )

    do s = 1, CondNum
      xyzf_QMixflux(:,:,kz,IdxCG(s)) =                            &
        & MAX(                                                    &
        &       0.0d0,                                            & 
        &     - Bulk * xyz_VelSfc(:,:,kz)                         & 
        &       * (  xyzf_QMixAll(:,:,kz,s)                       &  
        &          - SvapPress( SpcWetID(IdxCC(s)), TempSfc )     &
        &             / PressSfc * (MolWtWet(IdxCG(s)) / MolWtDry)&
        &         ) / ( z_dz(kz) * 5.0d-1 )                       &
        &    )
    end do

    !地表面運動量フラックスによる変化率を計算
    !  * 単位は m/s^2
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1
    pyz_MomFlux(:,:,kz) =                                         &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk * abs(pyz_VelXSfc(:,:,kz)) * pyz_VelX(:,:,kz)  &
      &       / ( z_dz(kz) * 5.0d-1 )                             &
      &   )

    xqz_MomFlux(:,:,kz) =                                         &
      & MAX(                                                      &
      &       0.0d0,                                              & 
      &     - Bulk * abs(xyz_VelYSfc(:,:,kz)) * xqz_VelY(:,:,kz) &
      &       / ( z_dz(kz) * 5.0d-1 )                             &
      &   )

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux
    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux
    pyz_DVelXDt  = pyz_DVelXDt0  + pyz_MomFlux
    xqz_DVelYDt  = xqz_DVelYDt0  + xqz_MomFlux

    call HistoryAutoPut(Time, 'PTempSfc', xyz_HeatFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(Time, 'VelXSfc',  pyz_MomFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(Time, 'VelYSfc',  xqz_MomFlux(1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(Time, trim(SpcWetSymbol(s))//'_Sfc', xyzf_Qmixflux(1:nx,1:ny,1:nz,s))
    end do
  end subroutine Surfaceflux_Bulk_forcing
  
end module Surfaceflux_bulk
