!= continumm absorption calculation module
!
!= continuum absorption calculation module
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20160713 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]


module mod_contabs

  ! 物理定数設定1
  ! physical constants setteing 1
  !
  use constants0, only: GasRUniv

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

  !
  ! continuum absorption calculation
  !
!  use mod_contabs_mtckd, only: calc_contabs_mtckd

  implicit none

  private

  public :: cont_absorpcs
  public :: cont_absorpk

contains

  subroutine cont_absorpk(          & !2016/06/21
    & Flag_CONT, mmax, kmax,        & !(in )
    & r_Press, r_Temp, rm_MixRatio, & !(in )
    & m_MolWt, MinWN, MaxWN, ResWN, & !(in )
    & a_alpha,                      & !(in )
    & rmw_CAbsorpk                  ) !(out)
    !
    ! 
    ! 
    !
    integer, intent(in ) :: Flag_CONT
    integer, intent(in ) :: mmax
    integer, intent(in ) :: kmax
    real(8), intent(in ) :: r_Press    (0:kmax) ![Pa]
    real(8), intent(in ) :: r_Temp     (0:kmax) ![K]
    real(8), intent(in ) :: rm_MixRatio(0:kmax, 0:  50)
    real(8), intent(in ) :: m_MolWt    (        0:mmax)
    real(8), intent(in ) :: MinWN ![m-1]
    real(8), intent(in ) :: MaxWN ![m-1]
    real(8), intent(in ) :: ResWN ![m-1]
    real(8), intent(in ) :: a_alpha(1:50) ![kg2 m-4 molecule-1]
!    real(8), intent(out) :: w_WaveNum (0:int((MaxWN-MinWN)/ResWN))
    real(8), intent(out) :: rmw_CAbsorpk(0:kmax,0:mmax,0:int((MaxWN-MinWN)/ResWN))
    integer :: k_w, wmax
    integer :: k_m, k_r
    integer :: k_mol
    real(8) :: contVMR(1:5)
    real(8) ::  w_WaveNum (0:int((MaxWN-MinWN)/ResWN))

    wmax = int((MaxWN - MinWN)/ResWN)

    call cont_absorpcs(               &
      & Flag_CONT, mmax, kmax,        & !(in )
      & r_Press, r_Temp, rm_MixRatio, & !(in )
      & m_MolWt, MinWN, MaxWN, ResWN, & !(in )
      & a_alpha,                      & !(in )
      & rmw_CAbsorpk                  ) !(out)

    do k_m = 0, mmax
      do k_w = 0, wmax
      rmw_CAbsorpk(0:kmax,k_m,k_w) = rmw_CAbsorpk(0:kmax,k_m,k_w) *  AvogadroNum/m_MolWt(k_m)
      end do
    end do

  end subroutine cont_absorpk

  subroutine cont_absorpcs(         &
    & Flag_CONT, mmax, kmax,        & !(in )
    & r_Press, r_Temp, rm_MixRatio, & !(in )
    & m_MolWt, MinWN, MaxWN, ResWN, & !(in )
    & a_alpha,                      & !(in )
    & rmw_CAbsorpCS                 ) !(out)
    !
    ! 
    ! 
    !
    integer, intent(in ) :: Flag_CONT
    integer, intent(in ) :: mmax
    integer, intent(in ) :: kmax
    real(8), intent(in ) :: r_Press    (0:kmax) ![Pa]
    real(8), intent(in ) :: r_Temp     (0:kmax) ![K]
    real(8), intent(in ) :: rm_MixRatio(0:kmax, 0:  50)
    real(8), intent(in ) :: m_MolWt    (        0:mmax)
    real(8), intent(in ) :: MinWN ![m-1]
    real(8), intent(in ) :: MaxWN ![m-1]
    real(8), intent(in ) :: ResWN ![m-1]
    real(8), intent(in ) :: a_alpha(1:50) ![kg2 m-4 molecule-1]
!    real(8) :: a_alpha(0:50) ![kg2 m-4 molecule-1]
!    real(8), intent(out) :: w_WaveNum (0:int((MaxWN-MinWN)/ResWN))
    real(8), intent(out) :: rmw_CAbsorpCS(0:kmax,0:mmax,0:int((MaxWN-MinWN)/ResWN))
    integer :: k_w, wmax
    integer :: k_m, k_r
    integer :: k_mol
    real(8) :: contVMR(1:5)
    real(8) ::  w_WaveNum (0:int((MaxWN-MinWN)/ResWN))
    integer :: fi = 21
    character(99) :: CharDMMY


    wmax = int((MaxWN - MinWN)/ResWN)

    if( Flag_CONT == 0 ) then

      rmw_CAbsorpCS = 0.0_8

    elseif( Flag_CONT == 1) then
      !
      ! MT_CKD continuum model
      !
      print *, 'This DCRTM version DOES NOT contain MT_CKD continuum model'
      print *, 'You do NOT use Flag_CONT == 1: stop'
      stop

    elseif( Flag_CONT == 3) then
      !
      ! H2O のみ連続吸収を持つ
      ! only H2O has continuum
      !
      print *, 'Now CONSTRACTING: stop'
      stop

    elseif( Flag_CONT == 2) then


!      open(fi, file='./src/optdep/contabs_RhoSquareCS_alpha')
!      read(fi, * ) CharDMMY

!      do k_m = 1, mmax
!        read(fi, * ) CharDMMY
!        read(fi, * ) a_alpha(k_m)
!        print *, 'alpha=', a_alpha(k_m)
!      end do

!      close(fi)

      do k_m = 1, mmax
        do k_r = 0, kmax
          do k_w = 0, wmax
            call contabs_RhoSquareCS( &
              & r_Press(k_r), r_Temp(k_r), m_MolWt(k_m), &
              & a_alpha(k_m), rmw_CAbsorpCS(k_r,k_m,k_w))
          end do
        end do
      end do

    else

      print *, '"This program does not support Flag_CONT:"', Flag_CONT
      print *, 'program stop'
      stop

    end if

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

  end subroutine cont_absorpcs

  subroutine contabs_RhoSquareCS( &
    & Press, Temp, MolWt, alpha,  &
    & AbsorpCS)
    !
    ! continuum absorption is proportional to square of gas density
    ! Pollack et al. (1993)
    !
    real(8), intent(in ) :: Press
    real(8), intent(in ) :: Temp
    real(8), intent(in ) :: MolWt
    real(8), intent(in ) :: alpha ![kg2 m-4 molecule-1]
    real(8), intent(out) :: AbsorpCS

    AbsorpCS = alpha * Press * Press * MolWt * MolWt/ &
      & (GasRUniv * GasRUniv * Temp * Temp)

  end subroutine contabs_RhoSquareCS

end module mod_contabs
