!= dcrtm 鉛直フラックス, 加熱率 計算 プログラム
!
!= dcrtm vertical flux calculation main program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20160720 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2015. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

program calc_heatrate
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

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

  ! netcdf データ入出力
  ! netcdf data input or output
  !
  use gtool_history

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

  ! 加熱率 計算
  ! heating rate calculation
  !
  use mod_HR, only: SumHeatingRate, &
    &               SumHeatingRateOD, &
    &               HeatingRate,     &
    &               HeatingRateOD


  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(99) :: InFileProf
                          ! 入力 ファイル名 (圧力, 温度, 混合比)
                          ! input file name (pressure , temperature, mixing ratio)
  character(99) :: InFileOptDep
                          ! 入力 ファイル名 (optical depth)
                          ! input file name (optical depth)
  character(99) :: InFileFlux
                          ! 入力 ファイル名 (フラックス)
                          ! input file name (flux)
  character(99) :: OutFileHeatRate
                          ! 出力ファイル名
                          ! output file name
  real(8) :: ResWN
                          ! 波数分解能 [m-1]
                          ! resolution of wavenumber [m-1]
  real(8), allocatable :: z_Press(:)
                          ! 気圧 (整数レベル) [Pa]
                          ! pressure (full level) [Pa]
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: r_Temp(:)
                          ! 気温 (半整数レベル) [K]
                          ! temperature (half level) [K]
  real(8), allocatable :: m_MolWt(:)
                          ! 分子量 [kg mol-1]
                          ! molecular weight [kg mol-1]
  real(8), allocatable :: rm_MixRatio(:,:)
                          ! 体積混合比 (半整数レベル) [1]
                          ! volume mixing ratio (half level) [1]
  real(8), allocatable :: m_CpMol(:)
                          ! 定圧モル比熱 [J mol-1 K-1]
                          ! molar heat capacity [J mol-1 K-1]
  real(8), allocatable :: r_CpMix(:)
                          ! 混合大気の定圧比熱 (半整数レベル) [J kg-1 K-1]
                          ! heat capacity of mix atmosphere (half level) [J kg-1 K-1]
  real(8), allocatable :: z_CpMix(:)
                          ! 混合大気の定圧比熱 (整数レベル) [J kg-1 K-1]
                          ! heat capacity of mix atmosphere (full level) [J kg-1 K-1]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-1]
  real(8), allocatable :: rw_OptDepTOA(:,:)
                          ! 大気上端からの光学的厚さ (半整数レベル) [1]
                          ! optical depth from top of atmosphere (half level) [1]
  real(8), allocatable :: rw_RadFluxUp(:,:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2 m]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2 m]
  real(8), allocatable :: rw_RadFluxDn(:,:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2 m]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2 m]
  real(8), allocatable :: rw_SolFluxUp(:,:)
                          ! 中心星入射に由来する上向きフラックス [W m-2 m]
                          ! upward flux originated from emittance of central star [W m-2 m]
  real(8), allocatable :: rw_SolFluxDn(:,:)
                          ! 中心星入射に由来する下向きフラックス [W m-2 m]
                          ! downward flux originated from emittance of central star [W m-2 m]
  real(8) :: MolWt0
                          ! 非吸収分子の平均分子量 [kg mol-1]
                          ! molecular weight of non-absorption gas [kg mol-1]
  real(8) :: CpMol0
                          ! 非吸収分子の定圧モル比熱 [J mol-1 K-1]
                          ! molar heat capacity of non-absorption gas [J mol-1 K-1]
  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: fi   = 13
                          ! 装置番号(Cp file)
                          ! device number : Cp file                          
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: kmax
                          ! 鉛直レベル最大値
                          ! vertical level maximum
  integer :: mmax
                          ! 分子数最大値
                          ! molecular number maximum
  integer :: k_w, k_r, k_m
                          ! 作業変数
                          ! work variables
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array
  character(99) :: charDMMY
                          ! dammy 配列
                          ! dammy array

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: r_RadFluxUp(:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of planetary atmosphere [W m-2]
  real(8), allocatable :: r_RadFluxDn(:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of planetary atmosphere [W m-2]
  real(8), allocatable :: r_RadNetFlux(:)
                          ! 惑星大気の射出に由来する正味上向きフラックス [W m-2]
                          ! net upward flux originated from emittance of planetary atmosphere [W m-2]
  real(8), allocatable :: r_SolFluxUp(:)
                          ! 中心星入射に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: r_SolFluxDn(:)
                          ! 中心星入射に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: r_SolNetFlux(:)
                          ! 中心星入射に由来する正味上向きフラックス [W m-2]
                          ! net upward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: z_RadHeatRate(:)
                          ! 惑星大気の射出に由来する加熱率 [K s-1]
                          ! heating ratio originated from emittance of planetary atmosphere [K s-1]
  real(8), allocatable :: z_SolHeatRate(:)
                          ! 中心星入射に由来する加熱率 [K s-1]
                          ! heating ratio originated from emittance of central star [K s-1]
  real(8), allocatable :: zw_RadHeatRate(:,:)
                          ! 惑星大気の射出に由来する加熱率 [K s-1]
                          ! heating ratio originated from emittance of planetary atmosphere [K s-1]
  real(8), allocatable :: zw_SolHeatRate(:,:)
                          ! 中心星入射に由来する加熱率 [K s-1]
                          ! heating ratio originated from emittance of central star [K s-1]
  real(8), allocatable :: r_MolAve(:)
                          ! 平均分子量 (半整数レベル) [kg mol-1]
                          ! average molerular weight [kg mol-1]
  real(8), allocatable :: m_SumMixRatio(:)
                          ! 混合比の和
                          ! summation of volume mixing ratio
  real(8) :: MaxMixRatio
                          ! 最大混合比
                          ! maximum volume mixing ratio

  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /heatrate_nml/ &
    & InFileProf,         &
    & InFileOptDep,       &
    & InFileFlux,         &
    & MolWt0,             &
    & CpMol0,             &
    & OutFileHeatRate

  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=heatrate_nml)
  close(fnml)
  !print *, MinWN, MaxWN, ResWN

  ! 配列サイズを調べる
  ! check array size
  !
  call HistoryGetPointer(''//InFileProf//'', 'r_Press', a_Axis )

  kmax = size( a_Axis ) - 1

  Nullify( a_Axis )

  !print *, kmax

  call HistoryGetPointer(''//InFileProf//'', 'm_MolWt', a_Axis)

  mmax = size( a_Axis )

  Nullify( a_Axis )


  call HistoryGetPointer(''//InFileFlux//'', 'w_WaveNum', a_Axis )

  wmax = size( a_Axis ) - 1

  ResWN = a_Axis(2) - a_Axis(1)

  Nullify( a_Axis )

  !print *, wmax
  !print *, ResWN

  !wmax = abs((MaxWN - MinWN)/ResWN)

  allocate( z_Press       (1:kmax                ) )
  allocate( r_Press       (0:kmax                ) )
  allocate( r_Temp        (0:kmax                ) )
  allocate( r_MolAve      (0:kmax                ) )
  allocate( m_MolWt       (        0:mmax        ) )
  allocate( rm_MixRatio   (0:kmax, 0:mmax        ) )
  allocate( m_CpMol       (        0:mmax        ) )
  allocate( m_SumMixRatio (        0:mmax        ) )
  allocate( r_CpMix       (0:kmax                ) )
  allocate( z_CpMix       (1:kmax                ) )
  allocate( w_WaveNum     (                0:wmax) )
  allocate( rw_OptDepTOA  (0:kmax,         0:wmax) )
  allocate( rw_RadFluxUp  (0:kmax,         0:wmax) )
  allocate( rw_RadFluxDn  (0:kmax,         0:wmax) )
  allocate( rw_SolFluxUp  (0:kmax,         0:wmax) )
  allocate( rw_SolFluxDn  (0:kmax,         0:wmax) )
  allocate( zw_RadHeatRate(1:kmax,         0:wmax) )
  allocate( zw_SolHeatRate(1:kmax,         0:wmax) )
  allocate( r_RadFluxUp   (0:kmax                ) )
  allocate( r_RadFluxDn   (0:kmax                ) )
  allocate( r_SolFluxUp   (0:kmax                ) )
  allocate( r_SolFluxDn   (0:kmax                ) )
  allocate( r_RadNetFlux  (0:kmax                ) )
  allocate( r_SolNetFlux  (0:kmax                ) )
  allocate( z_RadHeatRate (1:kmax                ) )
  allocate( z_SolHeatRate (1:kmax                ) )

  ! 鉛直プロファイルの読み込み
  ! read profile data
  !
  call HistoryGet(''//InFileProf//'', 'r_Press', r_Press)
  call HistoryGet(''//InFileProf//'', 'r_Temp', r_Temp)
  call HistoryGet(''//InFileProf//'', 'm_MolWt', m_MolWt(1:mmax))
  call Historyget(''//InFileProf//'', 'rm_MixRatio', rm_MixRatio(:,1:mmax))

  call HistoryGet(''//InFileOptDep//'', 'w_WaveNum', w_WaveNum)
  call HistoryGet(''//InFileOptDep//'', 'rw_OptDep', rw_OptDepTOA)

  call HistoryGet(''//InFileFlux//'', 'rw_RadFluxUp', rw_RadFluxUp)
  call HistoryGet(''//InFileFlux//'', 'rw_RadFluxDn', rw_RadFluxDn)
  call HistoryGet(''//InFileFlux//'', 'rw_SolFluxUp', rw_SolFluxUp)
  call HistoryGet(''//InFileFlux//'', 'rw_SolFluxDn', rw_SolFluxDn)

  ! 波数積分
  ! summation flux
  !
  r_RadFluxUp(:) = (rw_RadFluxUp(:,0) + rw_RadFluxUp(:,wmax))*0.5_8
  r_RadFluxDn(:) = (rw_RadFluxDn(:,0) + rw_RadFluxDn(:,wmax))*0.5_8
  r_SolFluxUp(:) = (rw_SolFluxUp(:,0) + rw_SolFluxUp(:,wmax))*0.5_8
  r_SolFluxDn(:) = (rw_SolFluxDn(:,0) + rw_SolFluxDn(:,wmax))*0.5_8

  do k_w = 1, wmax-1
    r_RadFluxUp(:) = r_RadFluxUp(:) + rw_RadFluxUp(:,k_w)
    r_RadFluxDn(:) = r_RadFluxDn(:) + rw_RadFluxDn(:,k_w)
    r_SolFluxUp(:) = r_SolFluxUp(:) + rw_SolFluxUp(:,k_w)
    r_SolFluxDn(:) = r_SolFluxDn(:) + rw_SolFluxDn(:,k_w)
  end do

  r_RadFluxUp(:) = r_RadFluxUp(:)*ResWN
  r_RadFluxDn(:) = r_RadFluxDn(:)*ResWN
  r_SolFluxUp(:) = r_SolFluxUp(:)*ResWN
  r_SolFluxDn(:) = r_SolFluxDn(:)*ResWN

  r_RadNetFlux = r_RadFluxUp - r_RadFluxDn
  r_SolNetFlux = r_SolFluxUp - r_SolFluxDn

  ! 非吸収分子の体積混合比の計算
  ! calculate volume mixing ratio of non-absrption gas
  !
  m_MolWt(0) = MolWt0

  !if(m_MolWt(0) > 0.0_8) then

  rm_MixRatio(:,0) = 1.0_8
  do k_r = 0, kmax

    do k_m = 1, mmax
      rm_MixRatio(k_r,0) = rm_MixRatio(k_r,0) - rm_Mixratio(k_r,k_m)
    end do

    if(rm_MixRatio(k_r,0) < 0.0_8) then
      print *, 'sum of vlume mixing ratio over 1.0:', 1.0_8 - rm_MixRatio(k_r,0), k_r

      MaxMixRatio = rm_MixRatio(k_r,0)
      do k_m = 1, mmax
        if(rm_MixRatio(k_r,k_m) > MaxMixRatio) then
          MaxMixRatio = rm_MixRatio(k_r,k_m)
        end if
      end do

      do k_m = 1, mmax
        if(rm_MixRatio(k_r,k_m) == MaxMixRatio) then
          exit
        end if
      end do
      rm_MixRatio(k_r,k_m) = rm_MixRatio(k_r,k_m) + rm_MixRatio(k_r,0)
      rm_MixRatio(k_r,0) = 0.0_8

      print *, 'modify volume mixing ratio, k_r=', k_r
      do k_m = 1, mmax
        print *, k_m, rm_MixRatio(k_r,k_m)
      end do
    end if

  end do


  m_SumMixRatio = 0.0_8
  do k_r = 0, kmax
    m_SumMixRatio(0:mmax) = m_SumMixRatio(0:mmax) + rm_MixRatio(k_r,0:mmax)
  end do


  ! 平均分子量の計算
  ! calculate average molecular weight
  !
  do k_r = 0, kmax
    r_MolAve(k_r) = 0.0_8
    do k_m = 0, mmax
      r_MolAve(k_r) = r_MolAve(k_r) + m_MolWt(k_m)*rm_MixRatio(k_r,k_m)
    end do
  end do

  ! 定圧比熱の計算
  ! calculate heat capacity
  !
  m_CpMol = 0.0_8
  m_CpMol(0) = CpMol0

  open(fi, file='./src/data/CpMol')
  read(fi, * ) CharDMMY

  do k_m = 1, mmax
    read(fi, * ) CharDMMY
    read(fi, * ) m_CpMol(k_m)
    print *, 'CpMol=', m_CpMol(k_m)
  end do

  close(fi)  

  r_CpMix = 0.0_8

  do k_r = 0, kmax
    do k_m = 0, mmax
      r_CpMix(k_r) = r_CpMix(k_r) + m_CpMol(k_m) * rm_MixRatio(k_r,k_m)
    end do ! k_m in k_r
  end do ! k_r

  r_CpMix(:) = r_CpMix(:)/r_MolAve(:) ! [J kg-1 K-1]

  do k_r = 1, kmax
    z_CpMix(k_r) = ( &
      & r_CpMix(k_r-1) * r_Press(k_r-1)/r_Temp(k_r-1) + &
      & r_CpMix(k_r  ) * r_Press(k_r  )/r_Temp(k_r  ) ) &
      & /(r_Press(k_r-1)/r_Temp(k_r-1) + r_Press(k_r)/r_Temp(k_r))
    !print *, z_CpMix(k_r)
  end do !k_r

  do k_r = 1, kmax

    z_Press(k_r) = ( &
      & r_Press(k_r-1) * r_Press(k_r-1)/r_Temp(k_r-1) + &
      & r_Press(k_r  ) * r_Press(k_r  )/r_Temp(k_r  ) ) &
      & /(r_Press(k_r-1)/r_Temp(k_r-1) + r_Press(k_r)/r_Temp(k_r))

  end do


  do k_w = 0, wmax

    call HeatingRate(kmax,         & !(in )
      &                 r_Press,      & !(in )
      &                 z_CpMix,      & !(in )
      &                 rw_SolFluxUp(:,k_w), & !(in )
      &                 rw_SolFluxDn(:,k_w), & !(in )
      &                 zw_SolHeatRate(:,k_w) ) !(out)

    call HeatingRateOD(kmax,         & !(in )
      &                 r_Press,      & !(in )
      &                 z_CpMix,      & !(in )
      &                 rw_OptDepTOA(:,k_w), & !(in )
      &                 rw_RadFluxUp(:,k_w), & !(in )
      &                 rw_RadFluxDn(:,k_w), & !(in )
      &                 zw_RadHeatRate(:,k_w) ) !(out)
    
  end do

!  zw_RadHeatRate = abs(zw_RadHeatRate)


  call SumHeatingRate(kmax,         & !(in )
    &                 wmax,         & !(in )
    &                 ResWN,        & !(in )
    &                 r_Press,      & !(in )
    &                 z_CpMix,      & !(in )
    &                 rw_SolFluxUp, & !(in )
    &                 rw_SolFluxDn, & !(in )
    &                 z_SolHeatRate ) !(out)

  !do k_r = 1, kmax
  !  print *, z_SolHeatRate(k_r)
  !end do

  call SumHeatingRateOD(kmax,         & !(in )
    &                 wmax,         & !(in )
    &                 ResWN,        & !(in )
    &                 r_Press,      & !(in )
    &                 z_CpMix,      & !(in )
    &                 rw_OptDepTOA, & !(in )
    &                 rw_RadFluxUp, & !(in )
    &                 rw_RadFluxDn, & !(in )
    &                 z_RadHeatRate ) !(out)

  !do k_r = 1, kmax
  !  print *, z_RadHeatRate(k_r)
  !end do

  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                    &
    & file=''//OutFileHeatRate//'', &
    & title='vertical flux; '//OutFileHeatRate//'', &
    & source='sum_verticalprofHR.f90;'//InFileProf//'', &
    & institution='DCRTM',       &
    & dims=(/'r_Press  ','z_Press  ','w_WaveNum'/), dimsizes=(/kmax+1,kmax,wmax+1/),      &
    & longnames=(/'pressure  ', 'pressure  ','wavenumber'/),       &
    & units=(/'Pa ','Pa ','m-1'/)                    )

  call HistoryPut('r_Press  ',r_Press)
  call HistoryPut('z_Press  ',z_Press)
  call HistoryPut('w_WaveNum', w_WaveNum)

  call HistoryAddVariable( &
    & varname='zw_RadHeatRate', dims=(/'z_Press  ','w_WaveNum'/), &
    & longname='heating rate (planet)', units='K s-1 m', xtype='double')
  call HistoryPut('zw_RadHeatRate',zw_RadHeatRate)

  call HistoryAddVariable( &
    & varname='zw_SolHeatRate', dims=(/'z_Press  ','w_WaveNum'/), &
    & longname='heating rate (solar)', units='K s-1 m', xtype='double')
  call HistoryPut('zw_SolHeatRate',zw_SolHeatRate)

  call HistoryAddVariable( &
    & varname='z_RadHeatRate', dims=(/'z_Press  '/), &
    & longname='heating rate (planet)', units='K s-1 m', xtype='double')
  call HistoryPut('z_RadHeatRate',z_RadHeatRate)

  call HistoryAddVariable( &
    & varname='z_SolHeatRate', dims=(/'z_Press  '/), &
    & longname='heating rate (solar)', units='K s-1 m', xtype='double')
  call HistoryPut('z_SolHeatRate',z_SolHeatRate)

if(.false.) then
  call HistoryAddVariable( &
    & varname='r_RadFluxUp', dims=(/'r_Press'/), &
    & longname='upward flux(planet)', units='W m-2', xtype='double')
  call HistoryPut('r_RadFluxUp',r_RadFluxUp)

  call HistoryAddVariable( &
    & varname='r_RadFluxDn', dims=(/'r_Press'/), &
    & longname='upward flux(planet)', units='W m-2', xtype='double')
  call HistoryPut('r_RadFluxDn',r_RadFluxDn)

  call HistoryAddVariable( &
    & varname='r_SolFluxUp', dims=(/'r_Press'/), &
    & longname='upward flux(solar)', units='W m-2', xtype='double')
  call HistoryPut('r_SolFluxUp',r_SolFluxUp)

  call HistoryAddVariable( &
    & varname='r_SolFluxDn', dims=(/'r_Press'/), &
    & longname='upward flux(solar)', units='W m-2', xtype='double')
  call HistoryPut('r_SolFluxDn',r_SolFluxDn)

  call HistoryAddVariable( &
    & varname='r_RadNetFlux', dims=(/'r_Press'/), &
    & longname='net upward flux(planet)', units='W m-2', xtype='double')
  call HistoryPut('r_RadNetFlux',r_RadNetFlux)

  call HistoryAddVariable( &
    & varname='r_SolNetFlux', dims=(/'r_Press'/), &
    & longname='net upward flux(solar)', units='W m-2', xtype='double')
  call HistoryPut('r_SolNetFlux',r_SolNetFlux)
end if

  call HistoryClose

end program calc_heatrate
