!= 加熱率計算 サブルーチン
!
!= heating rate calculation sub-program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20160720 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

module mod_HR
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

  ! モジュール引用 ; USE statements
  !

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: Grav


  ! 宣言文 ; Declaration statements
  !
  implicit none

  private

  public :: SumHeatingRate
  public :: SumHeatingRateOD
  public :: HeatingRate
  public :: HeatingRateOD

contains

  subroutine SumHeatingRate(kmax,      & !(in ) 
    &                       wmax,      & !(in )
    &                       ResWN,     & !(in )
    &                       r_Press,   & !(in )
    &                       z_Cp,      & !(in )
    &                       rw_FluxUp, & !(in )
    &                       rw_FluxDn, & !(in )
    &                       z_HeatRate ) !(out)

    integer, intent(in ) :: kmax
    integer, intent(in ) :: wmax
    real(8), intent(in ) :: ResWN
    real(8), intent(in ) :: r_Press    (0:kmax)
    real(8), intent(in ) :: z_Cp       (1:kmax)
    real(8), intent(in ) :: rw_FluxUp  (0:kmax, 0:wmax)
    real(8), intent(in ) :: rw_FluxDn  (0:kmax, 0:wmax)
    real(8), intent(out) :: z_HeatRate (1:kmax)
    real(8) :: z_HeatRateW (1:kmax)
    integer :: k_r, k_w

    z_HeatRate = 0.0_8

    k_w = 0
    call HeatingRate(kmax,      & !(in )
      &              r_Press,   & !(in )
      &              z_Cp,      & !(in )
      &              rw_FluxUp(:,k_w),  & !(in )
      &              rw_FluxDn(:,k_w),  & !(in )
      &              z_HeatRateW ) !(out)

    z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)*0.5_8

    k_w = wmax
    call HeatingRate(kmax,      & !(in )
      &              r_Press,   & !(in )
      &              z_Cp,      & !(in )
      &              rw_FluxUp(:,k_w),  & !(in )
      &              rw_FluxDn(:,k_w),  & !(in )
      &              z_HeatRateW ) !(out)

    z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)*0.5_8

    do k_w = 1, wmax-1

      call HeatingRate(kmax,      & !(in )
        &              r_Press,   & !(in )
        &              z_Cp,      & !(in )
        &              rw_FluxUp(:,k_w),  & !(in )
        &              rw_FluxDn(:,k_w),  & !(in )
        &              z_HeatRateW ) !(out)

      z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)

    end do

    z_HeatRate(:) = z_HeatRate(:) * ResWN

  end subroutine SumHeatingRate

  subroutine HeatingRate(kmax,      & !(in )
    &                    r_Press,   & !(in )
    &                    z_Cp,      & !(in )
    &                    r_FluxUp,  & !(in )
    &                    r_FluxDn,  & !(in )
    &                    z_HeatRate ) !(out)

    integer, intent(in ) :: kmax
    real(8), intent(in ) :: r_Press    (0:kmax)
    real(8), intent(in ) :: z_Cp       (1:kmax)
    real(8), intent(in ) :: r_FluxUp   (0:kmax)
    real(8), intent(in ) :: r_FluxDn   (0:kmax)
    real(8), intent(out) :: z_HeatRate (1:kmax)
    integer :: k_r
    real(8) :: DelNetFlux
    real(8) :: EPS


    EPS = 1.0d-16

    do k_r = 1, kmax
      DelNetFlux = (r_FluxUp(k_r-1)-r_FluxDn(k_r-1)) - &
        &          (r_FluxUp(k_r  )-r_FluxDn(k_r  ))
      if(abs(DelNetFlux) < max(r_FluxUp(k_r-1), &
        &                      r_FluxUp(k_r  ), &
        &                      r_FluxDn(k_r-1), &
        &                      r_FluxDn(k_r  )  )*EPS ) then

        z_HeatRate(k_r) = 0.0_8
      else
        z_HeatRate(k_r) = DelNetFlux * Grav/&
          & (z_Cp(k_r) * (r_Press(k_r-1)-r_Press(k_r))) !* ResWN
      end if
    end do !k_r 

  end subroutine HeatingRate

  subroutine SumHeatingRateOD(kmax,      & !(in ) 
    &                       wmax,      & !(in )
    &                       ResWN,     & !(in )
    &                       r_Press,   & !(in )
    &                       z_Cp,      & !(in )
    &                       rw_OptDepTOA,& !(in )
    &                       rw_FluxUp, & !(in )
    &                       rw_FluxDn, & !(in )
    &                       z_HeatRate ) !(out)

    integer, intent(in ) :: kmax
    integer, intent(in ) :: wmax
    real(8), intent(in ) :: ResWN
    real(8), intent(in ) :: r_Press      (0:kmax)
    real(8), intent(in ) :: z_Cp         (1:kmax)
    real(8), intent(in ) :: rw_OptDepTOA (0:kmax, 0:wmax)
    real(8), intent(in ) :: rw_FluxUp    (0:kmax, 0:wmax)
    real(8), intent(in ) :: rw_FluxDn    (0:kmax, 0:wmax)
    real(8), intent(out) :: z_HeatRate   (1:kmax)
    real(8) :: z_HeatRateW (1:kmax)
    integer :: k_r, k_w

    z_HeatRate = 0.0_8

    k_w = 0
    call HeatingRateOD(kmax,      & !(in )
      &              r_Press,   & !(in )
      &              z_Cp,      & !(in )
      &              rw_OptDepTOA(:,k_w), & !(in )
      &              rw_FluxUp(:,k_w),  & !(in )
      &              rw_FluxDn(:,k_w),  & !(in )
      &              z_HeatRateW ) !(out)

    z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)*0.5_8

    k_w = wmax
    call HeatingRateOD(kmax,      & !(in )
      &              r_Press,   & !(in )
      &              z_Cp,      & !(in )
      &              rw_OptDepTOA(:,k_w), & !(in )
      &              rw_FluxUp(:,k_w),  & !(in )
      &              rw_FluxDn(:,k_w),  & !(in )
      &              z_HeatRateW ) !(out)

    z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)*0.5_8

    do k_w = 1, wmax-1

      call HeatingRateOD(kmax,      & !(in )
        &              r_Press,   & !(in )
        &              z_Cp,      & !(in )
        &              rw_OptDepTOA(:,k_w), & !(in )
        &              rw_FluxUp(:,k_w),  & !(in )
        &              rw_FluxDn(:,k_w),  & !(in )
        &              z_HeatRateW ) !(out)

      z_HeatRate(:) = z_HeatRate(:) + z_HeatRateW(:)

    end do

    z_HeatRate(:) = z_HeatRate(:) * ResWN

  end subroutine SumHeatingRateOD

  subroutine HeatingRateOD(kmax,       & !(in )
    &                      r_Press,    & !(in )
    &                      z_Cp,       & !(in )
    &                      r_OptDepTOA,& !(in )
    &                      r_FluxUp,   & !(in )
    &                      r_FluxDn,   & !(in )
    &                      z_HeatRate  ) !(out)

    integer, intent(in ) :: kmax
    real(8), intent(in ) :: r_Press    (0:kmax)
    real(8), intent(in ) :: z_Cp       (1:kmax)
    real(8), intent(in ) :: r_OptDepTOA(0:kmax)
    real(8), intent(in ) :: r_FluxUp   (0:kmax)
    real(8), intent(in ) :: r_FluxDn   (0:kmax)
    real(8), intent(out) :: z_HeatRate (1:kmax)
    integer :: k_r
    real(8) :: DelNetFlux
    real(8) :: EPS


    EPS = 1.0d-16

    do k_r = 1, kmax
      DelNetFlux = (r_FluxUp(k_r-1)-r_FluxDn(k_r-1)) - &
        &          (r_FluxUp(k_r  )-r_FluxDn(k_r  ))
      if(abs(DelNetFlux) < max(r_FluxUp(k_r-1), &
        &                      r_FluxUp(k_r  ), &
        &                      r_FluxDn(k_r-1), &
        &                      r_FluxDn(k_r  )  )*EPS ) then

        z_HeatRate(k_r) = 0.0_8
      else

        if( (r_OptDepTOA(k_r-1)-r_OptDepTOA(k_r)) < 1d-10 ) then
          z_HeatRate(k_r) = 0.0_8
        else

          z_HeatRate(k_r) = DelNetFlux * Grav/&
            & (z_Cp(k_r) * (r_Press(k_r-1)-r_Press(k_r))) !* ResWN
        end if

      end if
    end do !k_r 

  end subroutine HeatingRateOD

end module mod_HR
