!= 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_mtckd

!  implicit none

!  private

!  public :: calc_contabs 
!  public :: calc_contabs_mtckd

!contains

  subroutine calc_contabs_mtckd(   &
    & InPress, InTemp, ArrVMR_all, & !(in )
    & InWNmin, InWNmax, INWNres,   & !(in )
    & AbsorpCS                     ) !(out)
    !
    ! 
    ! 
    !
    real(8), intent(in ) :: InPress ![Pa]
    real(8), intent(in ) :: InTemp  ![K]
    real(8), intent(in ) :: ArrVMR_all(1:5)
    real(8), intent(in ) :: InWNmin ![m-1]
    real(8), intent(in ) :: InWNmax ![m-1]
    real(8), intent(in ) :: InWNres ![m-1]
    !real(8), intent(out) :: WaveNum (0:int((InWNmax-InWNmin)/InWNres))
    real(8), intent(out) :: AbsorpCS(0:int((InWNmax-InWNmin)/InWNres))
    real(8) :: Press
    real(8) :: WNmin
    real(8) :: WNmax
    real(8) :: WNres
    real(8) :: WaveNum (0:int((InWNmax-InWNmin)/InWNres))
    real(8) :: AbsorpCS_m(1:2, 1:5001)
    real(8) :: ArrVMR_dry(1:5)
    real(8) :: SumVMR
    real(8) :: N_wn
    integer :: Lmax
    integer :: k, k_w

    ! change VMR: VMR_all to VMR_dry
    !
!    SumVMR = 0.0_8
!    do k = 1, 5
!      SumVMR = SumVMR + ArrVMR_all(k)
!    end do

!    if( SumVMR > 1.0_8 ) then
!      print *, 'mixing ratio is over the range: SumVMR:', SumVMR
!      print *, 'stop (call_mtckd)'
!      stop
!    end if

    ! mt_ckd trick??
    !
    ! we need to chane the VMR of N2 for "correct" continuum calculation.
    ! for example, vmr(H2O) = 0.01, any other molecules are no absorption
    ! i.e., ArrVMR_all(2:5) = 0, MT_CKD does not work (output is bad value).
    ! this case, MT_CKD needs ArrVMR(5) = 1. - vmr(H2O)
    !

    if( ArrVMR_all(1) == 1.0_8 ) then ! All modeclules are H2O
      ArrVMR_dry(1) = 1.0_8
      ArrVMR_dry(2:5) = 0.0_8
    elseif( ArrVMR_all(1) == 0.5_8 ) then !
      !ArrVMR_all(:) = ArrVMR_all(:) * 10.0_8
      ArrVMR_dry(:) = ArrVMR_all(:)*10.0_8/(1.0_8 - ArrVMR_all(1)*10.0_8)
    else
      ArrVMR_dry(:) = ArrVMR_all(:)/(1.0_8 - ArrVMR_all(1))
    end if

    !print *, ArrVMR_dry

    ! change scale: SI to MT_CKD(mbr, cm-1)
    !
    Press = InPress * 1e-2_8
    WNmin = InWNmin * 1e-2_8
    WNmax = InWNmax * 1e-2_8
    WNres = InWNres * 1e-2_8

!    N_wn = 1.0_8 + (WNmax - WNmin)/WNres
    N_wn = 1.0_8 + int( (WNmax - WNmin)/WNres )
    Lmax = int( (N_wn - 1.0_8)/5000.0_8 ) + 1
!    print *, (WNmax - WNmin)/WNres
!    print *, 'N_wn:', N_wn, 'Lmax:', Lmax
!    print *, 'int(N_wn)', int(N_wn)
!    print *, int(N_wn)-5000*(Lmax-1)-1

    do k = 1, Lmax

      if( k == Lmax ) then

        call MT_CKD(Press, InTemp, ArrVMR_dry, &
          & WNmin + (k-1)*5000*WNres, WNmax, &
          & WNres, AbsorpCS_m)

        WaveNum(5000*(Lmax-1):int((WNmax-WNmin)/WNres)) = &
          & AbsorpCS_m(1,1:int(N_wn) - 5000*(Lmax-1)) * 100.0_8 ![m-1]

        !print *, '==check'
        !print *, 5000*(Lmax-1), int((WNmax-WNmin)/WNres)
        !print *, 1, int(N_wn) - 5000*(Lmax-1)

        !do k_w = 1, int(N_wn) - 5000*(Lmax-1)!-1
        !  print *, AbsorpCS_m(1,k_w)
        !end do

        if( ArrVMR_all(1) == 0.5_8 ) then
          AbsorpCS(5000*(Lmax-1):int((WNmax-WNmin)/WNres)) = &
            & AbsorpCS_m(2,1:int(N_wn) - 5000*(Lmax-1)) * 1e-4_8 *0.10_8![m2/molecule]
        else
          AbsorpCS(5000*(Lmax-1):int((WNmax-WNmin)/WNres)) = &
            & AbsorpCS_m(2,1:int(N_wn) - 5000*(Lmax-1)) * 1e-4_8 ![m2/molecule]
        end if

      else

        call MT_CKD(Press, InTemp, ArrVMR_dry, &
          & WNmin + (k-1)*5000*WNres, WNmin + (5000*k-1)*WNres, &
          & WNres, AbsorpCS_m)

        WaveNum ((k-1)*5000:5000*k-1) = AbsorpCS_m(1,1:5000) * 100.0_8 ![m-1]

        if( ArrVMR_all(1) == 0.5_8 ) then
          AbsorpCS((k-1)*5000:5000*k-1) = AbsorpCS_m(2,1:5000) * 1e-4_8 *0.10_8![m2/molecule]
        else
          AbsorpCS((k-1)*5000:5000*k-1) = AbsorpCS_m(2,1:5000) * 1e-4_8 ![m2/molecule]
        end if

      end if

    end do

  end subroutine calc_contabs_mtckd

!end module mod_contabs_mtckd

  Include 'sub_cntnm_progr_MO.f90' ! 2014/06/03
