!2014/05/20
!2014/05/20

program Main_LinebyLine

  use gtool_history
  use mod_LinebyLine
  implicit none

  real :: cputime1, cputime2

  real(8) :: AvogadroNum_ = 6.02214129e+23_8 !CODATA 2010
  real(8) :: Grav_ = 9.8_8

  real(8) :: GasRUniv = 8.3144621_8 !CODATA [J mol-1 K-1] 
  character(10) :: DataBase
  integer :: kmax
  real(8) :: Press(1:2)
  real(8) :: Temp(1:2)
  real(8) :: MixRatioH(1:2, 1:50)
  real(8) :: MinWaveNum = 0.0_8
  real(8) :: MaxWaveNum = 100.0_8
  real(8) :: ResWaveNum = 1.0_8
  real(8), allocatable :: OptDepth(:)
  real(8), allocatable :: OptDepthM(:)
  real(8), allocatable :: r_Press(:)
  real(8), allocatable :: r_Temp(:)
  real(8), allocatable :: r_OptDepTOA(:,:)
  real(8), allocatable :: WaveNum(:)
  real(8), allocatable :: Absorp(:)
  integer :: k_mol
  real(8) :: MolAve
  real(8) :: VolMixRatio

  integer :: mfi = 10
  integer :: mfo = 11
  integer :: mfc = 12
  integer :: mk, mkk

  integer :: LayNum
  integer :: LayNum_s
  integer :: LayNum_e

  character( 4) :: RunCode
  character(30) :: OutFileName
  character(15) :: InFileName
  character( 4) :: Lnum4
  character(30) :: TAPE5_head
  real(8) :: TAPE5_ZPT(1:2, 1:3)
  real(8), allocatable :: TAPE5_x(:, :)
  real(8), allocatable :: TAPE6_v(:, :)
  real(8), allocatable :: MolecAmount(:, :)
  real(8), allocatable :: z_Press(:)
  real(8), allocatable :: z_Temp(:)
  real(8), allocatable :: z_MixRatioH(:,:)
  character(16) charScale
  integer :: NumMol
  real(8) :: MolWeight
  real(8) :: MolWeight_m01 = 18.010565e-3
  real(8) :: MolWeight_m03 = 47.9847453e-3
  integer :: mfot = 15
  integer :: mfit = 16
  

  !==read condition file
  namelist /optdeplay_icrccm_nml/ &
    & RunCode, &
    & InFileName,  &
    & DataBase,    &
    & MinWaveNum,  &
    & MaxWaveNum,  &
    & ResWaveNum,  &
    & LayNum_s,    &
    & LayNum_e

  call cpu_time(cputime1)

  open(mfc, file='OptDepLay_ICRCCM.nml')
  read(mfc, nml=optdeplay_icrccm_nml)
  close(mfc)
  print *, RunCode, InFileName, DataBase, MinWaveNum, MaxWaveNum, ResWaveNum, LayNum_s, LayNum_e

  kmax = Int( ( MaxWaveNum - MinWaveNum ) /ResWaveNum )

  allocate( OptDepth (0:kmax) )
  allocate( OptDepthM(0:kmax) )
  allocate( Absorp   (0:kmax) )
  allocate( WaveNum  (0:kmax) )
  allocate( r_OptDepTOA(LayNum_s-1:LayNum_e,0:kmax) )
  allocate( r_Press    (LayNum_s-1:LayNum_e) )
  allocate( r_Temp     (LayNum_s-1:LayNum_e) )

  do mk = 0, kmax
    WaveNum(mk) = MinWaveNum + ResWaveNum * mk
  end do

  open (mfi, file='./data/'//InFileName//'')

  do mk = 1, 3 ! calculate mixing ratio from molecular density data
   read (mfi, *)
  end do
  
  read (mfi, '(A30)') TAPE5_head
  !print *, TAPE5_head
  !print *, TAPE5_head(30:30)

  read(TAPE5_head(30:30), '(I1)') NumMol
  !print *, NumMol
  close(mfi)

  allocate( TAPE5_x(1:2, 1:NumMol) )
  allocate( TAPE6_v(LayNum_s:LayNum_e, 1:9+NumMol) )
  allocate( MolecAmount(LayNum_s:LayNum_e, 0:NumMol) )
  allocate( z_Press(LayNum_s:LayNum_e) )
  allocate( z_Temp(LayNum_s:LayNum_e) )
  allocate( z_MixRatioH(LayNum_s:LayNum_e, 1:NumMol) )

  open(mfit, file='./data/TAPE6_ZPTMA_'//RunCode//'')
  read(mfit, *)
  read(mfit, *)
  read(mfit, *)
  read(mfit, *)

  do mk = LayNum_s, LayNum_e
    read(mfit, *) TAPE6_v(mk,1:9+NumMol)
    z_Press(mk) = TAPE6_v(mk, 6) * 100.0_8
    z_Temp(mk) = TAPE6_v(mk, 7)
    MolecAmount(mk, 0:NumMol) = TAPE6_v(mk, 8:8+NumMol) * 1.0e+4_8
    z_MixRatioH(mk, 1:NumMol) = MolecAmount(mk, 1:NumMol)/MolecAmount(mk,0)
    !print *, z_Press(mk), z_Temp(mk)
    !print *, MolecAmount(mk,:)
    !print *, z_MixRatioH(mk,:)
  end do

  ! calculate molecular weight
  ! almost all atmosphere is N2:78%, O2:21%, Ar:1%
  MolAve = 0.78 * 28.006147e-3 + 0.21 * 31.989830e-3 + 0.01 * 39.948e-3 ![kg/mol]

  r_OptDepTOA(LayNum_e, :) = 0.0_8

  do LayNum = LayNum_e,  LayNum_s, -1
    print *, 'LayNum:', LayNum
    if(LayNum < 10) then
      write(Lnum4, '(A3, I1)') '000', LayNum
    elseif(LayNum < 100) then
      write(Lnum4, '(A2, I2)') '00', LayNum
    elseif(LayNum < 1000) then
      write(Lnum4, '(A1, I3)') '0', LayNum
    elseif(LayNum < 10000) then
      write(Lnum4, '(I4)') LayNum
    else
      print *, 'layer number is over, layer number:', LayNum
      print *, 'layer number: stop'
      stop
    end if

    OutFileName = 'OptDepLay_ICRCCM_'//RunCode//'_Lay'//Lnum4//''

    open (mfi, file='./data/'//InFileName//'')

    do mk = 1, 11 ! calculate mixing ratio from molecular density data
     read (mfi, *)
    end do
  
    if(LayNum > 1) then
      do mk = 1, (LayNum - 1) * 2
        read (mfi, *)
      end do
    end if

    read(mfi, '(3E14.8e2, A16)') TAPE5_ZPT(1, 1:3), charScale
    read(mfi, *) TAPE5_x(1, 1:NumMol)
    !print *, TAPE5_ZPT(1, 1:3)
    !print *, charScale
    !print *, charScale(11:11)
    !print *, charScale(12:12)
    !print *, TAPE5_x(1, 1:NumMol)

    read(mfi, '(3E14.8e2, A16)') TAPE5_ZPT(2, 1:3), charScale
    read(mfi, *) TAPE5_x(2, 1:NumMol)
    !print *, TAPE5_ZPT(2, 1:3)
    !print *, charScale
    !print *, charScale(11:11)
    !print *, charScale(12:12)
    !print *, TAPE5_x(2, 1:NumMol)

    close(mfi)

    if(charScale(6:6)=='A') then
      !print *, 'charScale(6:6) is A'
      Press(1:2) = TAPE5_ZPT(1:2, 2) * 100.0_8 ! [Pa]
      !print *, Press(1:2)
    else
      print *, 'scale of p is ??. charScale(6:6):', charScale(6:6)
      print *, 'stop'
      stop
    end if

    r_Press(LayNum-1) = Press(1)
    r_Press(LayNum  ) = Press(2)

    if(charScale(7:7)=='A') then
      !print *, 'charScale(7:7) is A'
      Temp(1:2) = TAPE5_ZPT(1:2, 3) ! [K]
      !print *, Temp(1:2)
    else
      print *, 'scale of p is ??. charScale(7:7):', charScale(6:6)
      print *, 'stop'
      stop
    end if

    r_Temp(LayNum-1) = Temp(1)
    r_Temp(LayNum  ) = Temp(2)

    MixRatioH(:, :) = 0.0_8

    do mkk = 1, NumMol

      if(charScale(10+mkk:10+mkk) == 'A') then
        ! ppmv
        !print *, 'charScale is A'
        MixRatioH(1:2, mkk) = TAPE5_x(1:2, mkk) * 1.0e-6_8
        !print *, mkk, MixRatioH(1:2, mkk)
      elseif(charScale(10+mkk:10+mkk) == 'D') then
        ! [gm m-3]
        !print *, 'charScale is D'
        if(mkk==1) then
          MolWeight = MolWeight_m01
        elseif(mkk==3) then
          MolWeight = MolWeight_m03
        else
          print *, 'mol weight is ??(stop)'
          stop
        end if

        MixRatioH(1:2, mkk) = TAPE5_x(1:2, mkk) * Temp(1:2) * GasRUniv/(Press(1:2)*MolWeight) *0.0010_8
        !print *, mkk, MixRatioH(1:2, mkk)
      else
        print *, 'scale of molecular amount is ??. mkk:', mkk, 'charScale:', charScale(10+mkk:10+mkk)
        print *, 'stop'
        stop
      end if

    end do

!    call MolecularWeight(Press, Temp, MixRatioH, MolAve)

    ! print *, "mol ave", MolAve

    OptDepth(:) = 0.0_8
    !print *, 'start calculate optical depth'
    do k_mol = 1, NumMol

      !OptDepthM(:) = 0.0_8
      !print *, 'if: MixRatioH .ne. 0'
      if( z_MixRatioH(LayNum, k_mol) .ne. 0.0_8 ) then
        !call mol_opticaldepth(DataBase, kmax, k_mol, Press, Temp, &
        !  & MixRatioH(:, k_mol), MinWaveNum, ResWaveNum, MolAve, &
        !  & OptDepthM)

        !print *, 'call mol_absorption, k_mol:', k_mol
        !print *, 'DataBase:', DataBase
        !print *, 'kmax:', kmax
        !print *, 'Press:', Press
        !print *, 'Temp:', Temp
        !print *, 'MixRatioH(:,k_mol):', MixRatioH(:, k_mol)
        !print *, 'MinWaveNum:', MinWaveNum
        !print *, 'ResWaveNum:', ResWaveNum
        !call mol_absorption(DataBase, kmax, k_mol, Press, Temp, &
        !  & MixRatioH(:, k_mol), MinWaveNum, ResWaveNum, &
        !  & Absorp(0:kmax))

        !VolMixRatio = (MixRatioH(1,k_mol) * Press(1)/Temp(1) + MixRatioH(2,k_mol) * Press(2)/Temp(2))/&
        !  & ( Press(1)/Temp(1) + Press(2)/Temp(2) )

        !OptDepthM(0:kmax) = VolMixRatio * abs(Press(1) - Press(2)) * AvogadroNum_/&
        !  & (Grav_ * MolAve) * Absorp(0:kmax)

        !print *, 'k_mol=', k_mol
        !print *, 'VolMixRatio=', VolMixRatio
        !print *, 'Molecular Amount=', MolecAmount(LayNum, k_mol) !VolMixRatio * abs(Press(1) - Press(2)) * AvogadroNum_/(Grav_ * MolAve)

        call linebyline(DataBase, kmax, k_mol, z_Press(LayNum), &
          & z_Temp(LayNum), z_MixRatioH(LayNum, k_mol), MinWaveNum, ResWaveNum, &
          & Absorp)

        OptDepthM(0:kmax) = MolecAmount(LayNum, k_mol) * Absorp(0:kmax)


      end if

!    open (mfot, file='test_OptDepLay03')
!    write(mfot, *) 'LayNum:', LayNum, 'k_mol', k_mol, 'MixRatioH:', MixRatioH(:,k_mol)
!    do mk = 0, kmax
!      write(mfot, *) int(MinWaveNum + ResWaveNum * mk), OptDepthM(mk)
!    end do

      OptDepth(:) = OptDepth(:) + OptDepthM(:)  
    end do
!  print *, "P=", Press, "T=", Temp
if(.false.) then
    ! 出力ファイルの設定==SET
    !open (mfo, file=OutFileName)
    open (mfo, file=OutFileName)
    write(mfo, '(A1,A4,A8,I4,A11,E10.4e2,A1,E10.4e2,A9,E9.3e2,A1,E9.3e2)') &
      & '#',RunCode,',LayNum:', LayNum, ',Press[Pa]:',Press(1),';',Press(2),',Temp[K]:',Temp(1),';', Temp(2)
    write(mfo, *) "# Wavenumber[m-1], OpticalDepth[1]"

    do mk = 0, kmax
      write(mfo, *) int(MinWaveNum + ResWaveNum * mk), OptDepth(mk)
    end do
  
    close(mfo)
end if
    r_OptDepTOA(LayNum-1, :) = r_OptDepTOA(LayNum, :) + OptDepth(:)

  end do ! LayNum

  call HistoryCreate(                                    &
    & file='OptDepTOA_ICRCCM_'//RunCode//'_LBLRTM-PT.nc', &
    & title='Optical Depth from TOA; ICRCCM; '//RunCode//'', &
    & source='main_linebyline_OptDepLay_ICRCCM_LBLRTM-MA.f90:'//RunCode//'', &
    & institution='DCRTM',       &
    & dims=(/'p ', 'wn'/), dimsizes=(/     33, 300001/),      &
    & 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='OptDep', dims=(/'p ', 'wn'/), &
    & longname='optical depth from TOA', units='1', xtype='double')
  call HistoryPut('OptDep',r_OptDepTOA)

  call HistoryClose


  call cpu_time(cputime2)
  print *, "k_mol;", k_mol, "cputime;", cputime2 - cputime1

end program Main_LinebyLine
