!20140514
!20140604

program main_flux_icrccm_dcpam

  use gtool_history
  use gridset, only: GridsetSet, GridsetPrint, imax, jmax, kmax
  use rad_rte_two_stream_app, only:RadRTETwoStreamAppInit, RadRTETwoStreamAppLW
  use mod_flux
  implicit none

!  real(8) :: Pai = 3.1415926535897932_8
!  real(8) :: GasRUniv = 8.314_8      !Gas constant
!  real(8) :: LatentHeat = 43655.0_8  !Latent heat
  real(8) :: Grav = 9.8_8              !Acceleration of gravity
  real(8) :: AvogadroNum = 6.02214129e+23_8 !CODATA 2010
                                            ! $ \N_A $ [mol-1]
  integer :: NumMinWaveNum
  integer :: NumRangeWaveNum
  real(8), allocatable :: r_Press(:)
  !real(8), allocatable :: z_Press(:)
  real(8), allocatable :: r_Temp(:)
  !real(8), allocatable :: z_Temp(:)
  integer, allocatable :: WaveNum(:)     !2014/02/28 real(8) -> integer
  !real(8), allocatable :: z_AbsorpCo(:,:)
  !real(8), allocatable :: ScatCS(:)
  !real(8), allocatable :: z_DelOptDep(:,:)
  real(8), allocatable :: r_RadUwFlux(:,:)
  real(8), allocatable :: r_RadDwFlux(:,:)
  real(8), allocatable :: r_FluxUp(:)
  real(8), allocatable :: r_FluxDn(:)
  real(8), allocatable :: z_FluxConv(:)
  real(8), allocatable :: r_OptDepTOA(:,:)
  !real(8) :: ScatCS_0
  integer :: k, kk !, kkk
  integer :: fo = 10
  integer :: mfc = 11
  character( 4) :: RunCode
  character(10) :: DataBase
  character( 3) :: charData
  integer :: k_w
  integer :: Flag_LandF = 1

  integer :: fot = 13

  ! for dcpam5 radiation code
  !integer :: imax
  !integer :: jmax
  !integer :: kmax
  real(8), allocatable :: xyz_SSA          (:, :, :) !single scahharing albedo =0.0
  real(8), allocatable :: xyz_AF           (:, :, :) !asynmetory factor = 0
  real(8), allocatable :: xyr_OptDep       (:, :, :) !opt depth from TOA
  real(8), allocatable :: xy_SurfAlbedo    (:, :) !albedo =0
  real(8), allocatable :: xyr_PFInted      (:, :, :) !planck function integral kv
  real(8), allocatable :: xy_SurfPFInted   (:, :) !planck function in land surface
  real(8), allocatable :: xy_SurfDPFDTInted(:, :) !=0 planck function の温度微分
  real(8), allocatable :: xyr_RadUwFlux    (:, :, :) !FluxUp
  real(8), allocatable :: xyr_RadDwFlux    (:, :, :) !FluxDn
  real(8), allocatable :: xyra_DelRadUwFlux(:, :, :, :) !
  real(8), allocatable :: xyra_DelRadDwFlux(:, :, :, :)

  !==read condition file
  namelist /flux_icrccm_nml/ &
    & RunCode,         &
    & DataBase,        &
    & kmax,            &
    & NumMinWaveNum,   &
    & NumRangeWaveNum

  open(mfc, file='Flux_ICRCCM.nml')
  read(mfc, nml=flux_icrccm_nml)
  close(mfc)
  print *, RunCode, DataBase, kmax, NumMinWaveNum, NumRangeWaveNum

  write(charData, '(A1, A2)') 'H', DataBase(9:10)




  call RadRTETwoStreamAppInit

  !== set1: grid
  imax = 1
  jmax = 1
  !kmax = 32

  call GridsetSet(imax, jmax, kmax)
  !call GridsetPrint

  !== set0: filename
  !Infilename='IAPWS_k00600R-1T2000' !F1 が点く場合は注意!'
  !NumMinWaveNum = 1
  !NumRangeWaveNum = 300000
  !RunCode = 'C35x'

  allocate (xyz_SSA          (0:imax-1, 1:jmax, 1:kmax) )!single scahharing albedo =0.0
  allocate (xyz_AF           (0:imax-1, 1:jmax, 1:kmax) )!asynmetory factor = 0
  allocate (xyr_OptDep       (0:imax-1, 1:jmax, 0:kmax) )!opt depth from TOA
  allocate (xy_SurfAlbedo    (0:imax-1, 1:jmax)         )!albedo =0
  allocate (xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax) )!planck function integral kv
  allocate (xy_SurfPFInted   (0:imax-1, 1:jmax)         )!planck function in land surface
  allocate (xy_SurfDPFDTInted(0:imax-1, 1:jmax)         )!=0 planck function の温度微分
  allocate (xyr_RadUwFlux    (0:imax-1, 1:jmax, 0:kmax) )!FluxUp
  allocate (xyr_RadDwFlux    (0:imax-1, 1:jmax, 0:kmax) )!FluxDn
  allocate (xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )!
  allocate (xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )


  allocate( r_Press(0:kmax) )
  !allocate( z_Press(1:kmax) )
  allocate( r_Temp (0:kmax) )
  !allocate( z_Temp (0:kmax) )
  allocate( WaveNum(1:NumRangeWaveNum)            )
  !allocate( z_AbsorpCo(1:kmax, 1:NumRangeWaveNum) )
  !allocate( ScatCS(1:NumRangeWaveNum)             )
  !allocate( z_DelOptDep(1:kmax, 1:NumRangeWaveNum) )
  allocate( r_RadUwFlux(0:kmax, 1:NumRangeWaveNum) )
  allocate( r_RadDwFlux(0:kmax, 1:NumRangeWaveNum) )
  allocate( r_FluxUp(0:kmax) )
  allocate( r_FluxDn(0:kmax) )
  allocate( z_FluxConv(1:kmax) )

  allocate( r_OptDepTOA(0:kmax, 1:NumRangeWaveNum) )
  !== read vertical profile & optical depth
  !
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'p', r_Press)
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'wn', WaveNum, range='wn=1:300000')
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'temp', r_Temp)
  call HistoryGet('./output_OptDepTOA/OptDepTOA_ICRCCM_'//charData//'_'//RunCode//'.nc', 'OptDep', r_OptDepTOA, range='wn=1:300000')

 ! open(fot, file='test20140528')


  ! No Rayliegh Schattering
  !

  !! Calculate Flux
  !

  ! dcpam5
  xyz_SSA          (:, :, :) = 0.0_8 !single scahharing albedo =0.0
  xyz_AF           (:, :, :) = 0.0_8 !asynmetory factor = 0
  !xyr_OptDep       (:, :, :) !opt depth from TOA
  xy_SurfAlbedo    (:, :)    = 0.0_8 !albedo =0
  !xyr_PFInted      (:, :, :) !planck function integral kv
  !xy_SurfPFInted   (:, :) !planck function in land surface
  xy_SurfDPFDTInted(:, :)    = 0.0_8 !=0 planck function の温度微分


  ! dcpam5 radiation code
  do k = 1, NumRangeWaveNum

    !do kk = 1, kmax
    !  xyz_SSA(:,:,kk) = ScatCS(k) /(z_AbsorpCo(kk,k) + ScatCS(k))
    !  if(xyz_SSA(0,1,kk) == 1.0_8) then ! 2014/04/21
    !    xyz_SSA(:,:,kk) = 1.0_8 - 10.0**(-16) !2014/04/22 !0.999999_8
    !  end if
    !end do


    do kk = 0, kmax
      xyr_OptDep(:,:,kk) = r_OptDepTOA(kk, k)
    end do

!    xyr_OptDep(:,:,kmax) = 0.0_8
!    do kk = kmax-1, 0, -1
!      xyr_OptDep(:,:,kk) = xyr_OptDep(:,:,kk+1) + z_DelOptDep(kk+1,k)
!    end do


    if(Flag_LandF == 1) then  ! 2014/04/28
      xy_SurfPFInted(:, :) = 3.1415926535897932_8 * integral_planck(WaveNum(k) - 0.5_8, WaveNum(k) + 0.5_8, r_Temp(0))
    elseif(Flag_LandF == 0) then
      xy_SurfPFInted(:, :) = 0.0_8
    end if

    do kk = 0, kmax
      xyr_PFInted(:, :,kk) = 3.1415926535897932_8 * integral_planck(WaveNum(k) - 0.5_8, WaveNum(k) + 0.5_8, r_Temp(kk))
    end do

    ! dcpam5 radiation code
    call RadRTETwoStreamAppLW(                           &
      & xyz_SSA, xyz_AF,                                 & ! (in)
      & xyr_OptDep,                                      & ! (in)
      & xy_SurfAlbedo,                                   & ! (in)
      & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
      & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
      & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
      & )

    r_RadUwFlux(0:kmax,k) = xyr_RadUwFlux(0,1,0:kmax)
    r_RadDwFlux(0:kmax,k) = xyr_RadDwFlux(0,1,0:kmax)

    !write(fot, *) 'xyr_OptDep'
    !do kk = 0, kmax
    !  write(fot, *) xyr_OptDep(:,:,kk)
    !end do

    !write(fot, *) 'xyr_PFInted'
    !do kk = 0, kmax
    !  write(fot, *) xyr_PFInted(:,:,kk)
    !end do

    !write(fot, *) 'xy_SurfPFInted'
    !write(fot, *) xy_SurfPFInted

    !write(fot, *) 'xyr_RadUwFlux'
    !do kk = 0, kmax
    !  write(fot, *) xyr_RadUwFlux(0,1,kk)
    !end do

    !write(fot, *) 'xyr_RadDwFlux'
    !do kk = 0, kmax
    !  write(fot, *) xyr_RadDwFlux(0,1,kk)
    !end do

  end do

  !open(fo, file='Flux_ICRCCM_'//charData//'_'//RunCode//'')
  !write(fo, *) '# Flux; ICRCCM; '//DataBase//','//RunCode//''
  !write(fo, *) '# wavenumber[m-1], FluxUp(r_p=0, r_p=32), FluxDn(r_p=0, r_p=32)'

  !do k = NumMinWaveNum, NumRangeWaveNum
  !  write(fo, '(67E22.12e3)') WaveNum(k), r_RadUwFlux(0:32,k), r_RadDwFlux(0:32,k)
  !end do

!stop
  call HistoryCreate(                                    &
    & file='Flux_ICRCCM_'//charData//'_'//RunCode//'.nc', &
    & title='Flux; ICRCCM; '//DataBase//','//RunCode//'', &
    & source='main_Flux_ICRCCM_dcpam.f90;'//DataBase//','//RunCode//'', &
    & institution='DCRTM',       &
    & dims=(/'p ', 'wn'/), dimsizes=(/     33, 300000/),      &
    & longnames=(/'pressure  ', 'wavenumber'/),       &
    & units=(/'Pa ', 'm-1'/)                    )

  call HistoryPut('p ',r_Press)
  call HistoryPut('wn',WaveNum)

  call HistoryAddVariable( &
    & varname='temp', dims=(/'p'/), &
    & longname='temperature', units='K', xtype='double')
  call HistoryPut('temp',r_Temp)

  call HistoryAddVariable( &
    & varname='FluxUp', dims=(/'p ', 'wn'/), &
    & longname='upward flux', units='W m-2 m', xtype='double')
  call HistoryPut('FluxUp',r_RadUwFlux)

  call HistoryAddVariable( &
    & varname='FluxDn', dims=(/'p ', 'wn'/), &
    & longname='downward flux', units='W m-2 m', xtype='double')
  call HistoryPut('FluxDn',r_RadDwFlux)


  call HistoryClose



  deallocate (xyz_SSA           )!single scahharing albedo =0.0
  deallocate (xyz_AF            )!asynmetory factor = 0
  deallocate (xyr_OptDep        )!opt depth from TOA
  deallocate (xy_SurfAlbedo     )!albedo =0
  deallocate (xyr_PFInted       )!planck function integral kv
  deallocate (xy_SurfPFInted    )!planck function in land surface
  deallocate (xy_SurfDPFDTInted )!=0 planck function の温度微分
  deallocate (xyr_RadUwFlux     )!FluxUp
  deallocate (xyr_RadDwFlux     )!FluxDn
  deallocate (xyra_DelRadUwFlux )!
  deallocate (xyra_DelRadDwFlux )

  deallocate( r_Press     )
  !deallocate( z_Press     )
  deallocate( r_Temp      )
  !deallocate( z_Temp      )
  deallocate( WaveNum     )
  !deallocate( z_AbsorpCo  )
  !deallocate( ScatCS      )
  !deallocate( z_DelOptDep )
  deallocate( r_RadUwFlux )
  deallocate( r_RadDwFlux )
  deallocate( r_FluxUp    )
  deallocate( r_FluxDn    )
  deallocate( z_FluxConv  )

 
end program main_flux_icrccm_dcpam
