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

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

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

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

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(99) :: InFileFlux
                          ! 入力 ファイル名
                          ! input file name
  character(99) :: OutFileVFlux
                          ! 出力ファイル名
                          ! output file name
  real(8) :: ResWN
                          ! 波数分解能 [m-1]
                          ! resolution of wavenumber [m-1]
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  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]

  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: kmax
                          ! 鉛直レベル最大値
                          ! vertical level maximum
  integer :: k_w, k_r
                          ! 作業変数
                          ! work variables
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: r_RadFluxUp(:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: r_RadFluxDn(:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: r_RadNetFlux(:)
                          ! 惑星大気の射出に由来する正味上向きフラックス [W m-2]
                          ! net upward flux originated from emittance of planetary atmospheire [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]


  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /vflux_nml/ &
    & InFileFlux,      &
    & OutFileVFlux

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

  print *, InFileFlux
  print *, OutFileVFlux

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

  kmax = size( a_Axis ) - 1

  Nullify( a_Axis )

  !print *, kmax

  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( r_Press       (0:kmax        ) )
  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( 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        ) )


  ! 鉛直プロファイルの読み込み
  ! read profile data
  !
  call HistoryGet(''//InFileFlux//'', 'r_Press', r_Press)
  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

  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                    &
    & file=''//OutFileVFlux//'', &
    & title='vertical flux; '//OutFileVFlux//'', &
    & source='sum_verticalprof.f90;'//InFileFlux//'', &
    & institution='DCRTM',       &
    & dims=(/'r_Press'/), dimsizes=(/kmax+1/),      &
    & longnames=(/'pressure'/),       &
    & units=(/'Pa'/)                    )

  call HistoryPut('r_Press',r_Press)

  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)

  call HistoryClose

end program sum_verticalflux
