!= continumm absorption calculation module
!
!= continuum absorption calculation module
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20161125 $
! 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)

    rmw_CAbsorpCS = 0.0_8

    if( Flag_CONT == 0 ) then

      rmw_CAbsorpCS = 0.0_8

    elseif( Flag_CONT == 1) then
      !
      ! MT_CKD continuum model
      !

      contVMR = 0.0_8
  if(.false.) then !2016/10/21
      if(mmax >= 22) then !N2
        do k_r = 0, kmax
          contVMR(5) = contVMR(5) + rm_MixRatio(k_r,22)
        end do
      else
        contVMR(5) = 0.0_8
      end if 

      if(mmax >=  7) then !O2
        do k_r = 0, kmax
          contVMR(4) = contVMR(4) + rm_MixRatio(k_r,7)
        end do
      else
        contVMR(4) = 0.0_8
      end if 

      if(mmax >=  3) then !O3
        do k_r = 0, kmax
          contVMR(3) = contVMR(3) + rm_MixRatio(k_r,3)
        end do
      else
        contVMR(3) = 0.0_8
      end if 

      if(mmax >=  2) then !CO2
        do k_r = 0, kmax
          contVMR(2) = contVMR(2) + rm_MixRatio(k_r,2)
        end do
      else
        contVMR(2) = 0.0_8
      end if 

      if(mmax >=  1) then !H2O
        do k_r = 0, kmax
          contVMR(1) = contVMR(1) + rm_MixRatio(k_r,1)
        end do
      else
        contVMR(1) = 0.0_8
      end if 
  end if


      do k_r = 0, kmax  ! 2016/10/21

        contVMR(1) = rm_MixRatio(k_r,1) !2016/10/21

      if( contVMR(1) .ne. 0.0_8 ) then
        !
        ! H2O かつ, CO2, O3, O2, N2 のいずれかに連続吸収あり
        ! H2O and (CO2 or O3 or O2 or N2) have continuum absorption
        !

        ! k_m = 1 に格納

        !do k_r = 0, kmax

          contVMR(1) = rm_MixRatio(k_r, 1)
          contVMR(2) = rm_MixRatio(k_r, 2)
          contVMR(3) = rm_MixRatio(k_r, 3)
          contVMR(4) = rm_MixRatio(k_r, 7)
          contVMR(5) = rm_MixRatio(k_r,22)

          if( contVMR(1) > 0.0_8 ) then
            contVMR(5) = 1.0_8 - contVMR(1) - contVMR(2) - contVMR(3) - contVMR(4)
          end if

          call calc_contabs_mtckd(                &
            & r_Press(k_r), r_Temp(k_r), contVMR, & !(in )
            & MinWN, MaxWN, ResWN,                & !(in )
            & rmw_CAbsorpCS(k_r,1,:)              ) !(out)

          rmw_CAbsorpCS(k_r,1,:) = rmw_CAbsorpCS(k_r,1,:)/rm_MixRatio(k_r,1)
  
        !end do

      else
        !
        ! H2O 連続吸収なし, CO2, O3, O2, N2 のいずれかに連続吸収あり
        ! H2O has no continuum, but (CO2 or O3 or O2 or N2) has continuum
        !

        ! k_m = 2 に格納
                
        !do k_r = 0, kmax
  
          contVMR(1) = 0.0_8
          contVMR(2) = rm_MixRatio(k_r, 2)
          contVMR(3) = rm_MixRatio(k_r, 3)
          contVMR(4) = rm_MixRatio(k_r, 7)
          contVMR(5) = rm_MixRatio(k_r,22)

          call calc_contabs_mtckd(                &
            & r_Press(k_r), r_Temp(k_r), contVMR, & !(in )
            & MinWN, MaxWN, ResWN,                & !(in )
            & rmw_CAbsorpCS(k_r,2,:)              ) !(out)

          rmw_CAbsorpCS(k_r,2,:) = rmw_CAbsorpCS(k_r,2,:)/rm_MixRatio(k_r,2)
  
        !end do

      end if

      end do

      !print *, '==continuum'
      !do k_w = 0, wmax
      !  print *, k_w, w_WaveNum(k_w)
      !end do

    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

  Include 'mod_contabs_mtckd.f90' ! 
