!= Module HeatFluxBS1998
!
! Authors::   KAWABATA Takuya
! Version::   $Id: surfaceflux_bs1998.f90,v 1.1 2012-05-30 01:33:03 odakker Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2012. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
! ä¸Šéƒ¨, ä¸‹éƒ¨å¢ƒç•Œã‹ã‚‰ã®ä¸€å®šã®ç†±ãƒ»é‹å‹•é‡ãƒ»ç‰©è³ªãƒ•ãƒ©ãƒƒã‚¯ã‚¹ã‚’ä¸ŽãˆãŸå ´åˆã«ã‚ˆã‚‹
! é‹å‹•é‡, æ¸©åº¦, å‡çµæˆåˆ†ã®å¤‰åŒ–çŽ‡ã‚’è¨ˆç®—ã™ã‚‹ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«.
! Module HeatFluxConst ã‚’å…ƒã«ä½œæˆ
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!
module Surfaceflux_BS1998
  !
  !ä¸Šä¸‹å¢ƒç•Œã§ã®å®šãƒ•ãƒ©ãƒƒã‚¯ã‚¹è¨ˆç®—ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«
  !

  !ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«èªã¿è¾¼ã¿
  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:  xyz_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 DExnerDt, only: xyz_DExnerDt_xyz_xyzf
  use Radiation_Simple, only: xyz_PTempRadBS1998

  !æš—é»™ã®åž‹å®£è¨€ç¦æ¢
  implicit none

  !å±žæ€§ã®æŒ‡å®š
  private

  !é–¢æ•°ã‚’ public ã«è¨å®š
  public surfaceflux_bs1998_init
  public surfaceflux_bs1998_forcing

  !å¤‰æ•°å®šç¾©
  real(DP), save :: SfcXMomFluxBtm = 0.0d0
  real(DP), save :: SfcXMomFluxTop = 0.0d0
  real(DP), save :: SfcYMomFluxBtm = 0.0d0
  real(DP), save :: SfcYMomFluxTop = 0.0d0
  real(DP), save :: SfcHeatFluxBtm = 0.0d0
  real(DP), save :: SfcHeatFluxTop = 0.0d0
  real(DP), save :: SfcQmixFluxBtm = 0.0d0
  real(DP), save :: SfcQmixFluxTop = 0.0d0

  character(*), parameter:: module_name = 'surfaceflux_bs1998'
                              ! ãƒ¢ã‚¸ãƒ¥ãƒ¼ãƒ«ã®åç§°.
                              ! Module name

contains
!!!------------------------------------------------------------------------!!!
  subroutine Surfaceflux_BS1998_init
    !
    !NAMELIST ã‹ã‚‰å¿…è¦ãªæƒ…å ±ã‚’èªã¿å–ã‚Š, å‡ºåŠ›å¤‰æ•°ã®è¨å®šã‚’è¡Œã†.
    !

    !æš—é»™ã®åž‹å®£è¨€ç¦æ¢
    implicit none

    !å†…éƒ¨å¤‰æ•°
    integer    :: l, unit

    !---------------------------------------------------------------
    ! NAMELIST ã‹ã‚‰æƒ…å ±ã‚’å–å¾—
    !
    NAMELIST /surfaceflux_bs1998_nml/ &
      & SfcXMomFluxBtm, SfcXMomFluxTop, SfcYMomFluxBtm, SfcYMomFluxTop, &
      & SfcHeatFluxBtm, SfcHeatFluxTop, SfcQmixFluxBtm, SfcQmixFluxTop

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

    if (myrank == 0) then
      call MessageNotify( "M", module_name, "SfcXMomFluxBtm = %f", &
        &                  d=(/SfcXMomFluxBtm/) )
      call MessageNotify( "M", module_name, "SfcXMomFluxTop = %f", &
        &                  d=(/SfcXMomFluxTop/) )
      call MessageNotify( "M", module_name, "SfcYMomFluxBtm = %f", &
        &                  d=(/SfcYMomFluxBtm/) )
      call MessageNotify( "M", module_name, "SfcYMomFluxTop = %f", &
        &                  d=(/SfcYMomFluxTop/) )
      call MessageNotify( "M", module_name, "SfcHeatFluxBtm = %f", &
        &                  d=(/SfcHeatFluxBtm/) )
      call MessageNotify( "M", module_name, "SfcHeatFluxTop = %f", &
        &                  d=(/SfcHeatFluxTop/) )
      call MessageNotify( "M", module_name, "SfcQmixFluxBtm = %f", &
        &                  d=(/SfcQmixFluxBtm/) )
      call MessageNotify( "M", module_name, "SfcQmixFluxTop = %f", &
        &                  d=(/SfcQmixFluxTop/) )
    end if

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

    call HistoryAutoAddVariable(      &
      & varname='ExnerSfcFlux',       &
      & dims=(/'x','y','z','t'/),         &
      & longname='surface exner function flux (heat flux divided by density and specific heat)', &
      & units='s-1',             &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfcFlux',    &
      & dims=(/'x','y','z','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','z','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','z','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='ExnerSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='exner function 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','z','t'/),          &
      & longname='surface heat flux',  &
      & units='W.m-2',                 &
      & xtype='float')

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

    call HistoryAutoAddVariable(                      &
      & varname='SfcYMomFlux',                        &
      & dims=(/'x','y','z','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','z','t'/),                              &
        & longname=trim(SpcWetSymbol(l))//' surface mass flux',  &
        & units='kg.m-2.s-1',                                    &
        & xtype='float')
    end do
  end subroutine Surfaceflux_BS1998_init

!!!------------------------------------------------------------------------!!!
  subroutine Surfaceflux_BS1998_forcing( &
    &   pyz_DVelXDt, xqz_DVelYDt, xyz_DPTempDt, xyz_DExnerDt, xyzf_DQMixDt &
    & )
    !

    !æš—é»™ã®åž‹å®£è¨€ç¦æ¢
    implicit none

    !å¤‰æ•°å®šç¾©
    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):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP) :: pyr_VelXflux (imin:imax,jmin:jmax,kmin:kmax) !é€Ÿåº¦ãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: xqr_VelYflux (imin:imax,jmin:jmax,kmin:kmax) !é€Ÿåº¦ãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: xyr_PTempFlux(imin:imax,jmin:jmax,kmin:kmax) !æ¸©åº¦ãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: xyr_ExnerFlux(imin:imax,jmin:jmax,kmin:kmax) !Exner é–¢æ•°å¼·åˆ¶é …
    real(DP) :: xyrf_QMixFlux(imin:imax,jmin:jmax,kmin:kmax,ncmax)
                                                             !ç‰©è³ªãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: pyr_SfcXMomFlux(imin:imax,jmin:jmax,kmin:kmax) !é‹å‹•é‡ãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: xqr_SfcYMomFlux(imin:imax,jmin:jmax,kmin:kmax) !é‹å‹•é‡ãƒ•ãƒ©ãƒƒã‚¯ã‚¹
    real(DP) :: xyr_SfcHeatFlux(imin:imax,jmin:jmax,kmin:kmax) !ç†±ãƒ•ãƒ©ãƒƒã‚¯ã‚¹

    real(DP) :: pyz_DVelXDtFlux (imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xqz_DVelYDtFlux (imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_DPTempDtFlux(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_DExnerDtFlux(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyzf_DQMixDtFlux(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    integer  :: s             !ãƒ«ãƒ¼ãƒ—å¤‰æ•°

pyr_VelXflux = 0.0d0
xqr_VelYflux = 0.0d0
xyr_PTempFlux = 0.0d0
xyr_ExnerFlux = 0.0d0
xyrf_QMixFlux = 0.0d0

pyr_SfcXMomFlux = 0.0d0
xqr_SfcYMomFlux = 0.0d0
xyr_SfcHeatFlux = 0.0d0

pyz_DVelXDtFlux = 0.0d0
xqz_DVelYDtFlux = 0.0d0
xyz_DPTempDtFlux = 0.0d0
xyz_DExnerDtFlux = 0.0d0
xyzf_DQMixDtFlux = 0.0d0

!-----

    ! NAMLIST ã§è¨å®šã—ãŸå€¤ã‚’ä»£å…¥
    !
!    xyr_SfcHeatFlux(:,:,1)  = - xyz_PTempRadBS1998(:,:,1) * CpDry * xyz_DensBZ(:,:,1) * xyz_dz(:,:,1)
!    xyr_SfcHeatFlux(:,:,nz) = - xyz_PTempRadBS1998(:,:,nz) * CpDry * xyz_DensBZ(:,:,nz) * xyz_dz(:,:,nz)
!    write(*,*) 'xyr_SfcHeatFluxBtm: ', xyr_SfcHeatFlux(1,1,1)
!    write(*,*) 'xyr_SfcHeatFluxTop: ', xyr_SfcHeatFlux(1,1,nz)

    xyr_SfcHeatFlux(:,:,1)  = SfcHeatFluxBtm
    xyr_SfcHeatFlux(:,:,nz) = SfcHeatFluxTop
!    write(*,*) 'xyr_SfcHeatFluxBtmNamelist: ', xyr_SfcHeatFlux(1,1,1)
!    write(*,*) 'xyr_SfcHeatFluxTopNamelist: ', xyr_SfcHeatFlux(1,1,nz)

    pyr_SfcXMomFlux(:,:,1)  = SfcXMomFluxBtm
    pyr_SfcXMomFlux(:,:,nz) = SfcXMomFluxTop

    xqr_SfcYMomFlux(:,:,1)  = SfcYMomFluxBtm
    xqr_SfcYMomFlux(:,:,nz) = SfcYMomFluxTop

    xyrf_QMixFlux(:,:,1,:)  = SfcQmixFluxBtm
    xyrf_QMixFlux(:,:,nz,:) = SfcQmixFluxTop

    ! é€Ÿåº¦ãƒ•ãƒ©ãƒƒã‚¯ã‚¹, æ¸©åº¦ãƒ•ãƒ©ãƒƒã‚¯ã‚¹ã¸ã®å¤‰æ›
    !
    pyr_VelXFlux  = pyr_SfcXMomFlux / xyz_DensBz
    xqr_VelYFlux  = xqr_SfcYMomFlux / xyz_DensBz
    xyr_PtempFlux = xyr_SfcHeatFlux / xyz_DensBz / CpDry

    ! Exner é–¢æ•°äºˆå ±å¼ã®å¼·åˆ¶é …ã®è¨ˆç®—
    !
    xyr_ExnerFlux = xyz_DExnerDt_xyz_xyzf(xyr_PTempFlux, xyrf_QMixFlux)

    ! æ™‚é–“å¤‰åŒ–çŽ‡ã®è¨ˆç®—
    !
    xyz_DPTempDtFlux = xyr_PTempFlux / xyz_dz
    xyz_DExnerDtFlux = xyr_ExnerFlux / xyz_dz
    do s = 1, ncmax
      xyzf_DQMixDtFlux(:,:,:,s) = xyrf_QMixFlux(:,:,:,s) / xyz_dz
    end do
    pyz_DVelXDtFlux  = pyr_VelXFlux / xyz_dz
    xqz_DVelYDtFlux  = xqr_VelYFlux / xyz_dz

    ! ä»®å¼•æ•°é…åˆ—ã¸ã®åŠ ç®—
    !
    xyz_DPTempDt = xyz_DPTempDt + xyz_DPTempDtFlux
    xyz_DExnerDt = xyz_DExnerDt + xyz_DExnerDtFlux
    do s = 1, ncmax
      xyzf_DQMixDt(:,:,:,s) = xyzf_DQMixDt(:,:,:,s) + xyzf_DQMixDtFlux(:,:,:,s)
    end do
    pyz_DVelXDt = pyz_DVelXDt + pyz_DVelXDtFlux
    xqz_DVelYDt = xqz_DVelYDt + xqz_DVelYDtFlux

    ! å‡ºåŠ›
    !
    call HistoryAutoPut(TimeN, 'PTempSfc', xyz_DPTempDtFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerSfc', xyz_DExnerDtFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfc',  pyz_DVelXDtFlux (1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfc',  xqz_DVelYDtFlux (1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', &
        & xyzf_DQMixDtFlux(1:nx,1:ny,1:nz,s))
    end do
    call HistoryAutoPut(TimeN, 'PTempSfcFlux', xyr_PTempFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerSfcFlux', xyr_ExnerFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfcFlux',  pyr_VelXFlux (1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfcFlux',  xqr_VelYFlux (1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcFlux', &
        & xyrf_QMixFlux(1:nx,1:ny,1:nz,s))
    end do

    call HistoryAutoPut(TimeN, 'SfcHeatFlux', xyr_SfcHeatFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'SfcXMomFlux', pyr_SfcXMomFlux(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'SfcXMomFlux', xqr_SfcYMomFlux(1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcMassFlux', &
        & xyz_DensBZ(1:nx,1:ny,1:nz) * xyrf_QMixFlux(1:nx,1:ny,1:nz,s))
    end do

  end subroutine Surfaceflux_BS1998_forcing

end module Surfaceflux_BS1998