!= Module HeatFlux
!
! Authors::   ODAKA Masatsugu 
! Version::   $Id: surfaceflux_diff.f90,v 1.8 2011-10-04 05:16:27 sugiyama Exp $
! Tag Name::  $Name: arare5-20111010 $
! 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_diff
  !
  !下部境界でのフラックスの計算モジュール
  !
  
  !モジュール読み込み
  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 方向の格子点間隔
  use basicset, only: xyz_ExnerBZ     !エクスナー関数の基本場
  use composition,only : GasNum, SpcWetSymbol
  use namelist_util, only: namelist_filename
  use timeset, only:  TimeN
  use setmargin,only: SetMargin_xyz, SetMargin_xyzf

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public に設定
  public surfaceflux_diff_init
  public surfaceflux_diff_forcing

  !変数定義
  real(DP), save  :: Kappa = 800.0d0

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

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer    :: l, unit

    !---------------------------------------------------------------    
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_diff_nml/ Kappa
    
    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=surfaceflux_diff_nml)
    close(unit)  

    if (myrank == 0) then 
      call MessageNotify( "M", "SurfaceFlux", "Kappa = %f", d=(/Kappa/))
    end if

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

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

  end subroutine Surfaceflux_Diff_init


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

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                                           !温位の擾乱成分    
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !温位の擾乱成分    
    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)               :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)               :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)
    real(DP)               :: xyz_Heatflux(imin:imax,jmin:jmax,kmin:kmax)
                                           !地表面熱フラックス
    real(DP)               :: xyzf_QMixflux(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    integer                :: kz            !配列添字
    integer                :: l             !ループ変数

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

    xyz_DPTempDt0 = xyz_DPTempDt
    xyzf_DQMixDt0 = xyzf_DQMixDt

    !地表面熱フラックスによる加熱率を計算
    !  * 単位は K/s
    !  * エクスナー関数は基本場の値で代表させる.     
    !  * 格子点 xz では, 物理領域の最下端の添え字は kz = 1

    xyz_HeatFlux(:,:,kz) =                                  &
      &  - Kappa * xyz_PTemp(:,:,kz) * xyz_ExnerBZ(:,:,kz)  &
      &    / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) 

    do l = 1, GasNum
      xyzf_QMixFlux(:,:,kz,l) =                       &
        & max(                                        &
        &       0.0d0,                                &
        &     - Kappa * xyzf_QMix(:,:,kz,l)           &
        &        / ( ( z_dz(kz) * 5.0d-1 ) ** 2.0d0 ) &
        &    )
    end do

    xyz_DPTempDt = xyz_DPTempDt0 + xyz_Heatflux
    xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Qmixflux

    call HistoryAutoPut(TimeN, 'PTempFlux', xyz_HeatFlux(1:nx,1:ny,1:nz))
    do l = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Flux', xyzf_Qmixflux(1:nx,1:ny,1:nz,l))
    end do    

    ! Set margin
    !
    call SetMargin_xyz(xyz_DPTempDt)
    call SetMargin_xyzf(xyzf_DQMixDt)

  end subroutine Surfaceflux_Diff_forcing
  
end module Surfaceflux_diff