!= Module HeatFlux
!
! Authors::   ODAKA Masatsugu, TAKAHASHI Yoshiyuki
! Version::   $Id: surfaceflux_bulk.f90,v 1.13 2011-10-10 15:44:24 yot Exp $
! Tag Name::  $Name: arare5-20111010 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!
! *** Explanation below are obsolete. (YOT, 2011/09/01) ***
!
!
! 下部境界からのフラックスによる温度と凝結成分の変化率をバルク方法に
! 基づいて計算するモジュール. これは中島 (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,  &
    &                 pyz_avr_xyz,  &
    &                 xqz_avr_xyz
  use basicset, only: xyz_ExnerBZ,  & !エクスナー関数の基本場
    &                 xyz_PressBZ,  & !
    &                 xyz_PTempBZ,  & !温位の基本場
    &                 xyz_TempBZ,   & !
    &                 xyzf_QMixBZ,  & !温位の基本場
    &                 xyz_DensBZ      !基本場の密度
  use constants,only: MolWtDry, PressBasis, TempSfc, PressSfc, CpDry, GasRDry
  use composition,only : IdxCC, IdxCG, SpcWetID, CondNum, MolWtWet, SpcWetSymbol
  use chemcalc,only : SvapPress
  use namelist_util, only: namelist_filename
  use timeset, only:  TimeN
  use setmargin,only: SetMargin_xyz, SetMargin_xyzf, SetMargin_pyz, SetMargin_xqz


  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

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

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


  real(DP), save  :: Vel0 = 0.0d0    !下層での水平速度嵩上げ値


  character(*), parameter:: module_name = 'surfaceflux_bulk'
                              ! モジュールの名称.
                              ! Module name

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

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer    :: l, unit

    !---------------------------------------------------------------
    ! NAMELIST から情報を取得
    !
    NAMELIST /surfaceflux_bulk_nml/ Bulk, Vel0

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=surfaceflux_bulk_nml)
    close(unit)  

    if (myrank == 0) then 
      call MessageNotify( "M", module_name, "Bulk = %f", d=(/Bulk/) )
      call MessageNotify( "M", module_name, "Vel0 = %f", d=(/Vel0/))
    end if


    call HistoryAutoAddVariable(      &
      & varname='PTempSfcFlux',       &
      & dims=(/'x','y','t'/),         &
      & longname='surface potential temperature flux (heat flux divided by density and specific heat)', &
      & units='K.m.s-1',             &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfcFlux',    &
      & dims=(/'x','y','t'/),     &
      & longname='surface flux of x-component of velocity (momentum flux divided by density)', &
      & units='m2.s-2',           &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelYSfcFlux',    &
      & dims=(/'x','y','t'/),     &
      & longname='surface flux of y-component of velocity (momentum flux divided by density)', &
      & units='m2.s-2',           &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_SfcFlux', & 
        & dims=(/'x','y','t'/),     &
        & longname='surface flux of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio (mass flux divided by density)',  &
        & units='m.s-1',    &
        & xtype='float')
    end do


    call HistoryAutoAddVariable(  &
      & varname='PTempSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='potential temperature tendency by surface flux', &
      & units='K.s-1',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='x-component velocity tendency by surface flux', &
      & units='m.s-2',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelYSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='y-component velocity tendency by surface flux', &
      & units='m.s-2',            &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Sfc', & 
        & dims=(/'x','y','z','t'/),     &
        & longname=trim(SpcWetSymbol(l))//' mixing ratio tendency by surface flux',  &
        & units='s-1',    &
        & xtype='float')
    end do


    call HistoryAutoAddVariable(       &
      & varname='SfcHeatFlux',         &
      & dims=(/'x','y','t'/),          &
      & longname='surface heat flux',  &
      & units='W.m-2',                 &
      & xtype='float')

    call HistoryAutoAddVariable(                      &
      & varname='SfcXMomFlux',                        &
      & dims=(/'x','y','t'/),                         &
      & longname='surface x-component momentum flux', &
      & units='kg.m-2.s-1',                           &
      & xtype='float')

    call HistoryAutoAddVariable(                      &
      & varname='SfcYMomFlux',                        &
      & dims=(/'x','y','t'/),                         &
      & longname='surface y-component momentum flux', &
      & units='kg.m-2.s-1',                           &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(                               &
        & varname=trim(SpcWetSymbol(l))//'_SfcMassFlux',         &
        & dims=(/'x','y','t'/),                              &
        & longname=trim(SpcWetSymbol(l))//' surface mass flux',  &
        & units='kg.m-2.s-1',                                    &
        & xtype='float')
    end do

  end subroutine Surfaceflux_Bulk_init


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

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    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)               :: py_VelXflux (imin:imax,jmin:jmax)
                                           !運動量フラックス
                                           !(strictly speaking, this value is not 
                                           !momentum flux, but is it divided by density)
    real(DP)               :: xq_VelYflux (imin:imax,jmin:jmax)
                                           !運動量フラックス
                                           !(strictly speaking, this value is not 
                                           !momentum flux, but is it divided by density)
    real(DP)               :: xy_PTempFlux(imin:imax,jmin:jmax)
                                           !地表面熱フラックス
                                           !(strictly speaking, this value is not 
                                           !heat flux, but is it divided by density and 
                                           !specific heat)
    real(DP)               :: xyf_QMixFlux(imin:imax,jmin:jmax,ncmax)
                                           !物質的フラックス
    real(DP)               :: xyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: xyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xyz 格子)
    real(DP)               :: pyz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (pyz 格子)
    real(DP)               :: xqz_AbsVel(imin:imax,jmin:jmax,kmin:kmax)
                                           !水平風速 (xqz 格子)
    real(DP)               :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
                                           !Total value of potential temperature
    real(DP)               :: xyz_TempAll (imin:imax,jmin:jmax,kmin:kmax)
                                           !Total value of temperature
    real(DP)               :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                           !Total value of mixing ratios
    real(DP)               :: xy_DPTempDtBulk(imin:imax,jmin:jmax)
                                        !potential temperature tendency by surface flux
    real(DP)               :: xyf_DQMixDtBulk(imin:imax,jmin:jmax, ncmax)
                                        !mixing ratio tendency by surface flux
    real(DP)               :: py_DVelXDtBulk (imin:imax,jmin:jmax)
                                        !x-component velocity tendency by surface flux
    real(DP)               :: xq_DVelYDtBulk (imin:imax,jmin:jmax)
                                        !y-component velocity tendency by surface flux

    real(DP)               :: xyz_DPTempDtBulk(imin:imax,jmin:jmax,kmin:kmax)
                                        ! variable for output
    real(DP)               :: xyzf_DQMixDtBulk(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                        ! variable for output
    real(DP)               :: pyz_DVelXDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                                        ! variable for output
    real(DP)               :: xqz_DVelYDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                                        ! variable for output

    real(DP)               :: ExnerBZSfc
                                        ! Basic state Exner function at the surface
    real(DP)               :: xy_PressSfc(imin:imax,jmin:jmax)
                                        ! Total pressure at the surface

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

    ! 初期化
    !
    kz = 1

    xyz_PTempAll  = xyz_PTemp + xyz_PTempBZ
    xyz_TempAll   = (xyz_Exner + xyz_ExnerBZ) * xyz_PTempAll
    xyzf_QMixAll  = xyzf_QMix + xyzf_QMixBZ

    ExnerBZSfc    = (PressSfc / PressBasis) ** (GasRDry / CpDry)
    xy_PressSfc   = PressBasis * ( ExnerBZSfc + xyz_Exner(:,:,kz) )**(CpDry / GasRDry)
                    ! Perturbation component of Exner function at the surface is assumed 
                    ! to be same as that at the lowest layer. (YOT, 2011/09/03)


    ! Velocities at xyz grid points are calculated.
    xyz_VelX = xyz_avr_pyz(pyz_VelX)
    xyz_VelY = xyz_avr_xqz(xqz_VelY)

    xyz_AbsVel = SQRT( xyz_VelX**2 + xyz_VelY**2 + Vel0**2 )
    pyz_AbsVel = pyz_avr_xyz(xyz_AbsVel)
    xqz_AbsVel = xqz_avr_xyz(xyz_AbsVel)


    ! Something like heat, mass, and momentum fluxes are calculated.
    ! The values below are not heat, mass, and momentum fluxes, but are those divided by
    ! by density and specific heat, density, and density, respectively.
    !
    xy_PTempFlux = - Bulk * xyz_AbsVel(:,:,kz)                                    &
      & * ( xyz_PTempAll(:,:,kz) - TempSfc / ( ExnerBZSfc + xyz_Exner(:,:,kz) ) )
    !
    xyf_QMixFlux = 0.0d0
    do s = 1, CondNum
      xyf_QMixFlux(:,:,IdxCG(s)) =                                 &
        &     - Bulk * xyz_AbsVel(:,:,kz)                          &
        &       * (                                                &
        &            xyzf_QMixAll(:,:,kz,s)                        &
        &          - SvapPress( SpcWetID(IdxCC(s)), TempSfc )      &
        &             / ( xy_PressSfc )                            &
        &             * (MolWtWet(IdxCG(s)) / MolWtDry)            &
        &         )
    end do
    !
    py_VelXFlux = - Bulk * pyz_AbsVel(:,:,kz) * pyz_VelX(:,:,kz)
    !
    xq_VelYFlux = - Bulk * xqz_AbsVel(:,:,kz) * xqz_VelY(:,:,kz)


    ! Something like heat flux and mass flux are restricted.
    ! This would be arbitrary treatment. 
    xy_PTempFlux = max( 0.0d0, xy_PTempFlux )
    xyf_QMixFlux = max( 0.0d0, xyf_QMixFlux )


    ! Tendencies by surface fluxes (convergences of fluxes) are calculated.
    !
    xy_DPTempDtBulk = - ( 0.0d0 - xy_PTempFlux ) / z_dz(kz)
    !
    xyf_DQMixDtBulk = - ( 0.0d0 - xyf_QMixFlux ) / z_dz(kz)
    !
    py_DVelXDtBulk = - ( 0.0d0 - py_VelXFlux ) / z_dz(kz)
    !
    xq_DVelYDtBulk = - ( 0.0d0 - xq_VelYFlux ) / z_dz(kz)


    ! Add tendency by surface flux convergence
    !
    xyz_DPTempDt(:,:,kz) = xyz_DPTempDt(:,:,kz) + xy_DPTempDtBulk
    do s = 1, ncmax
      xyzf_DQMixDt(:,:,kz,s) = xyzf_DQMixDt(:,:,kz,s) + xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDt (:,:,kz) = pyz_DVelXDt (:,:,kz) + py_DVelXDtBulk
    xqz_DVelYDt (:,:,kz) = xqz_DVelYDt (:,:,kz) + xq_DVelYDtBulk



    ! Output
    !
    xyz_DPTempDtBulk = 0.0d0
    xyzf_DQMixDtBulk = 0.0d0
    pyz_DVelXDtBulk  = 0.0d0
    xqz_DVelYDtBulk  = 0.0d0

    xyz_DPTempDtBulk(:,:,kz) = xy_DPTempDtBulk
    do s = 1, ncmax
      xyzf_DQMixDtBulk(:,:,kz,s) = xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDtBulk (:,:,kz) = py_DVelXDtBulk
    xqz_DVelYDtBulk (:,:,kz) = xq_DVelYDtBulk
    !
    call HistoryAutoPut(TimeN, 'PTempSfc', xyz_DPTempDtBulk(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfc',  pyz_DVelXDtBulk (1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfc',  xqz_DVelYDtBulk (1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', &
        & xyzf_DQMixDtBulk(1:nx,1:ny,1:nz,s))
    end do


    call HistoryAutoPut(TimeN, 'PTempSfcFlux', xy_PTempFlux(1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelXSfcFlux',  py_VelXFlux (1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelYSfcFlux',  xq_VelYFlux (1:nx,1:ny))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcFlux', &
        & xyf_QMixFlux(1:nx,1:ny,s))
    end do

    call HistoryAutoPut(TimeN, 'SfcHeatFlux', &
      & CpDry * xyz_DensBZ(1:nx,1:ny,1) * xy_PTempFlux(1:nx,1:ny) &
      &   * ( ExnerBZSfc + xyz_Exner(1:nx,1:ny,1) ) )
    call HistoryAutoPut(TimeN, 'SfcXMomFlux', &
      & xyz_DensBZ(1:nx,1:ny,1) * py_VelXFlux (1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'SfcYMomFlux', &
      & xyz_DensBZ(1:nx,1:ny,1) * xq_VelYFlux (1:nx,1:ny))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcMassFlux', &
        & xyz_DensBZ(1:nx,1:ny,1) * xyf_QMixFlux(1:nx,1:ny,s))
    end do

    ! Set Margin
    !
    call SetMargin_xyz( xyz_DPTempDt )
    call SetMargin_pyz( pyz_DVelXDt )
    call SetMargin_xqz( xqz_DVelYDt )
    call SetMargin_xyzf(xyzf_DQMixDt )

  end subroutine Surfaceflux_Bulk_forcing
  
end module Surfaceflux_bulk
