!= (火星大気向け) Non-LTE 放射モデル
!
!= Non-NLTE radiation model (for the Mars' atmosphere)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: rad_15m_NLTE.f90,v 1.1 2012/11/15 03:30:10 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module rad_15m_NLTE
  !
  != (火星大気向け) Non-LTE 放射モデル
  !
  != Non-NLTE radiation model (for the Mars' atmosphere)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! (火星大気向け) Non-LTE 放射モデル
  !
  ! Non-NLTE radiation model (for the Mars' atmosphere)
  !
  !== References
  !
  ! See Takahashi et al. (2003) for reference
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#rad_15m_NLTE_nml
  !

  ! モジュール引用 ; USE statements
  !

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
    &                 STRING     ! 文字列.       Strings.

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数.
                             ! Number of vertical level

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
!!$  public :: rad15mNLTE
  public :: rad15mNLTEMergeHR
  public :: rad15mNLTEInit

!!$  public :: rad15mNLTECalckMin
!!$  public :: rad15mNLTECalcWeight

  ! 公開変数
  ! Public variables
  !
  logical, save, public:: rad_15m_NLTE_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag

  ! 非公開変数
  ! Private variables
  !

  ! nl15fn               : maximum number of factors for 15 micron Non-LTE
  !                      : radiative cooling rate calculation
  ! nl15sn               : "reduced" optical depth for 15 micron Non-LTE
  !                      : radiative cooling rate calculation
  ! nl15fa               : parameter for 15 micron Non-LTE
  !                      : radiative cooling rate calculation

  integer , parameter :: nTable15mNLTE = 70
  real(DP)            :: a_Table15mNLTEsn( nTable15mNLTE )
                                        ! sigma N (non dimension)
  real(DP)            :: a_Table15mNLTEfa( nTable15mNLTE )
                                        ! E(tau)


  data a_Table15mNLTEsn &
    & / &
    & 0.1000e-6, 0.1000e-1, 0.1389e-1, 0.1931e-1, 0.2683e-1, &
    & 0.3728e-1, 0.5179e-1, 0.7197e-1, 0.1000e0 , 0.1389e0 , &
    & 0.1931e0 , 0.2683e0 , 0.3728e0 , 0.5179e0 , 0.7197e0 , &
    & 0.1000e1 , 0.1389e1 , 0.1931e1 , 0.2683e1 , 0.3728e1 , &
    & 0.5179e1 , 0.7197e1 , 0.1000e2 , 0.1389e2 , 0.1931e2 , &
    & 0.2683e2 , 0.3728e2 , 0.5179e2 , 0.7197e2 , 0.1000e3 , &
    & 0.1389e3 , 0.1931e3 , 0.2683e3 , 0.3728e3 , 0.5179e3 , &
    & 0.7197e3 , 0.1000e4 , 0.1389e4 , 0.1931e4 , 0.2683e4 , &
    & 0.3728e4 , 0.5179e4 , 0.7197e4 , 0.1000e5 , 0.1389e5 , &
    & 0.1931e5 , 0.2683e5 , 0.3728e5 , 0.5179e5 , 0.7197e5 , &
    & 0.1000e6 , 0.1389e6 , 0.1931e6 , 0.2683e6 , 0.3728e6 , &
    & 0.5179e6 , 0.7197e6 , 0.1000e7 , 0.1389e7 , 0.1931e7 , &
    & 0.2683e7 , 0.3728e7 , 0.5179e7 , 0.7197e7 , 0.1000e8 , &
    & 0.1389e8 , 0.1931e8 , 0.2683e8 , 0.3728e8 , 0.5179e8   &
    & /


  data a_Table15mNLTEfa &
    & / &
    & 0.5000e0 , 0.4993e0 , 0.4991e0 , 0.4987e0 , 0.4982e0 , &
    & 0.4975e0 , 0.4966e0 , 0.4953e0 , 0.4936e0 , 0.4913e0 , &
    & 0.4883e0 , 0.4844e0 , 0.4796e0 , 0.4734e0 , 0.4656e0 , &
    & 0.4557e0 , 0.4432e0 , 0.4279e0 , 0.4090e0 , 0.3864e0 , &
    & 0.3595e0 , 0.3284e0 , 0.2933e0 , 0.2550e0 , 0.2149e0 , &
    & 0.1750e0 , 0.1373e0 , 0.1041e0 , 0.7679e-1, 0.5566e-1, &
    & 0.4009e-1, 0.2892e-1, 0.2096e-1, 0.1530e-1, 0.1127e-1, &
    & 0.8387e-2, 0.6329e-2, 0.4851e-2, 0.3783e-2, 0.3005e-2, &
    & 0.2431e-2, 0.1999e-2, 0.1671e-2, 0.1415e-2, 0.1211e-2, &
    & 0.1044e-2, 0.9018e-3, 0.7798e-3, 0.6735e-3, 0.5801e-3, &
    & 0.4982e-3, 0.4267e-3, 0.3645e-3, 0.3108e-3, 0.2646e-3, &
    & 0.2249e-3, 0.1911e-3, 0.1622e-3, 0.1376e-3, 0.1167e-3, &
    & 0.9891e-4, 0.8382e-4, 0.7101e-4, 0.6013e-4, 0.5091e-4, &
    & 0.4307e-4, 0.3643e-4, 0.3079e-4, 0.2600e-4, 0.2194e-4  &
    & / 


  character(*), parameter:: module_name = 'rad_15m_NLTE'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_15m_NLTE.f90,v 1.1 2012/11/15 03:30:10 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  !**************************************************************************
  !     subroutine nlteradiation
  !     calculate Non-LTE radiative cooling rate
  !**************************************************************************
  ! pi
  ! amu      : atomic mass unit
  ! day      : one solar day (sec)
  ! grav     : gravitational acceleration
  ! im       : number of the vertical layers
  ! jm       : number of the meridional grids
  ! km       : number of the zonal grids
  ! ntemp    : neutral temperature (K)
  ! press    : pressure at midpoint of layer (Pa)
  ! cp       : specific heat at constant pressure
  ! q15m1    : infrared heating rate in 15 micron band below about 80 km
  ! rho      : mass density at midpoint of layers
  !**************************************************************************

  subroutine rad15mNLTE(                &
    & xyz_Press, xyz_Temp, xyz_VirTemp, &
    & xyz_DTempDt15mNLTE                &
    & )

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air

    real(DP), intent(in ) :: xyz_Press   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_VirTemp (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_DTempDt15mNLTE(0:imax-1, 1:jmax, 1:kmax)


    ! Local variables
    !
    real(DP) :: xyz_Rho(0:imax-1, 1:jmax, 1:kmax)
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: l

    real(DP) :: tau
    real(DP) :: e1, e2
    real(DP) :: ltau
    real(DP) :: ramda
    real(DP), parameter :: a10 = 2.44d0
    real(DP), parameter :: g10 = 2.0d0
    real(DP) :: kco2, ko
    real(DP) :: co2nd, ond
    real(DP) :: tmpfac
    real(DP) :: NLTECR
!!$    real(DP) :: xyz_Weight(0:imax-1, 1:jmax, 1:kmax)


    ! Variables for Reference Pressure
    !
    real(DP) :: p0, t0, ond0, f0

    ! kmin     : maximum index used for calculation of
    !          : atmospheric radiative cooling rate
!!$    integer  :: xy_kMin(0:imax-1, 1:jmax)

    real(DP), parameter :: AMU = 1.6605655d-27


    ! 初期化
    ! Initialization
    !
    if ( .not. rad_15m_NLTE_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

!!$    call calcimin(im,iup,jm,km,press,imin)
!!$    xy_kMin = kmax

!!$    call rad15mNLTECalcWeight( &
!!$      & xyz_Press,                   &
!!$!      &  xy_kMin,                    &
!!$      & xyz_Weight                   &
!!$      & )
!!$    xyz_Weight = 1.0_DP

    ! Set Variables for Reference Pressure
    !
    p0 = 1.2d-3 * 1.0d-6 * 1.0d5
    f0 = 0.5d0 * 1.0d-2
!            f0=1.0d0*1.0d-2


    xyz_Rho = xyz_Press / ( GasRDry * xyz_VirTemp )

    do j = 1, jmax
      do i = 0, imax-1

        ! Number Density is in CGS Unit
        !
        if( p0 <= xyz_Press(i,j,kmax) ) then
          t0 = xyz_Temp(i,j,kmax)
          ond0 = xyz_Rho(i,j,kmax) / ( 44.0d0 * AMU ) &
            * ( p0 / xyz_Press(i,j,kmax) ) &
            * f0 &
            * 1.0d-6
        else if( p0 <= xyz_Press(i,j,1) ) then
          search_p : do l = kmax-1, 2, -1
            if( p0 <= xyz_Press(i,j,l) ) exit search_p
          end do search_p
          t0 = ( xyz_Temp(i,j,l+1) - xyz_Temp(i,j,l) )        &
            &  / log( xyz_Press(i,j,l+1) / xyz_Press(i,j,l) ) &
            &  * log( p0 / xyz_Press(i,j,l) )                 &
            &  + xyz_Temp(i,j,l)
          ond0 =                                   &
            & xyz_Temp(i,j,l) / t0                 &
            &  * xyz_Rho(i,j,l) / ( 44.0d0 * AMU ) &
            &  * ( p0 / xyz_Press(i,j,l) )         &
            &  * f0                                &
            &  * 1.0d-6
        else
          write( 6, * ) 'Reference pressure or pressure is inappropriate.'
          write( 6, * ) 'Unexpected Error'
          write( 6, * ) 'Stop'
          stop
        endif

        do k = 1, kmax
          ! cgs unit
          co2nd = xyz_Rho(i,j,k) / ( 44.0d0 * AMU ) &
            & * 1.0d-6

!               kco2=1.0d-15
          ! from Lunt et al., 1985
          if( xyz_Temp(i,j,k) < 175.0_DP ) then
            kco2 = 4.2d-12 * exp( -2988.0d0 /175.0d0 + 303930.0d0 / ( 175.0d0 * 175.0d0 ) )
          else
            kco2 = 4.2d-12 * exp( -2988.0d0 / xyz_Temp(i,j,k) + 303930.0d0 / ( xyz_Temp(i,j,k) * xyz_Temp(i,j,k) ) )
          endif

          ! ond0 has already been set in cgs unit
          ond = t0 / xyz_Temp(i,j,k) &
            & * ond0                 &
            & * ( ( xyz_Press(i,j,k) /p0 )**(16.0d0/44.0d0) )

!               ko=1.5d-11*dexp(-800.0d0/ntemp(i,j,k))
          ! from Lopez-Valvelde and Lopez-Puertas, 1994, 
          ! Bougher et al., 1994
          ko = 3.0d-12
          tau = 6.4d-15 * xyz_Press(i,j,k) &
            & / ( 44.0d0 * AMU * Grav )    &
            & * 1.0d-4


!!$  integer , parameter :: nTable15mNLTE = 70
!!$  real(DP),           :: a_Table15mNLTEsn( nTable15mNLTE )
!!$                                        ! sigma N (non dimension)
!!$  real(DP),           :: a_Table15mNLTEfa( nTable15mNLTE )
!!$                                        ! E(tau)

!----------------------
          search_sn_1 : do l = 2, nTable15mNLTE-1
            if ( tau < a_Table15mNLTEsn(l) ) exit search_sn_1
          end do search_sn_1
          e1 =   ( a_Table15mNLTEfa(l) - a_Table15mNLTEfa(l-1) ) &
            &  / ( a_Table15mNLTEsn(l) - a_Table15mNLTEsn(l-1) ) &
            &  * ( tau - a_Table15mNLTEsn(l-1) )                 &
            &  + a_Table15mNLTEfa(l-1)
          if ( e1 > 0.5d0 ) e1 = 0.5d0
          if ( e1 < 0.0d0 ) e1 = 0.0d0
!----------------------
          search_sn_2 : do l = 2, nTable15mNLTE-1
            if ( ( tau / 2.0d0 ) <  a_Table15mNLTEsn(l) ) exit search_sn_2
          end do search_sn_2
          e2 =   ( a_Table15mNLTEfa(l) - a_Table15mNLTEfa(l-1) ) &
            &  / ( a_Table15mNLTEsn(l) - a_Table15mNLTEsn(l-1) ) &
            &  * ( ( tau / 2.0d0 ) - a_Table15mNLTEsn(l-1) )     &
            &  + a_Table15mNLTEfa(l-1)
          if ( e2 > 0.5d0 ) e2 = 0.5d0
          if ( e2 < 0.0d0 ) e2 = 0.0d0
!----------------------
          ltau = e1 + e2
          ramda = a10 / ( a10 + kco2 * co2nd + ko * ond )
          tmpfac = 0.5d0 * ramda * ltau / ( 1.0d0 - ramda + 0.5d0 * ramda * ltau )

          NLTECR = 1.33d-13 * g10 * exp( -960.0d0 / xyz_Temp(i,j,k) ) * co2nd &
            & * ( kco2 * co2nd + ko * ond ) * tmpfac &
            & * 1.0d-1 &
            & / ( xyz_Rho(i,j,k) * CpDry )

!!$          xyz_DTempDt15mNLTE(i,j,k) = &
!!$            & ( 1.0d0 - xyz_Weight(i,j,k) ) * ( -NLTECR )
          xyz_DTempDt15mNLTE(i,j,k) = - NLTECR
        end do

      end do
    end do

  end subroutine rad15mNLTE

  !----------------------------------------------------------------------------

  subroutine rad15mNLTEMergeHR(         &
    & xyz_Press, xyz_Temp, xyz_VirTemp, &
    & xyz_DTempDt                       &
    & )

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN                 ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.

    real(DP), intent(in   ) :: xyz_Press  (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in   ) :: xyz_VirTemp(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(inout) :: xyz_DTempDt(0:imax-1, 1:jmax, 1:kmax)


    ! Local variables
    !
    real(DP) :: xyz_Weight        (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DTempDt15mNLTE(0:imax-1, 1:jmax, 1:kmax)


    ! 初期化
    ! Initialization
    !
    if ( .not. rad_15m_NLTE_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    call rad15mNLTECalcWeight( &
      & xyz_Press,                   &
!      &  xy_kMin,                    &
      & xyz_Weight                   &
      & )
!!$    xyz_Weight = 1.0_DP

    call rad15mNLTE(                      &
      & xyz_Press, xyz_Temp, xyz_VirTemp, &
      & xyz_DTempDt15mNLTE                &
      & )

    xyz_DTempDt =                                              &
      &   xyz_Weight             * xyz_DTempDt                 &
      & + ( 1.0d0 - xyz_Weight ) * xyz_DTempDt15mNLTE


    call HistoryAutoPut( TimeN, 'DTempDt15mNLTE'   , xyz_DTempDt15mNLTE )
    call HistoryAutoPut( TimeN, 'DTempDtRadLMerged', xyz_DTempDt )


  end subroutine rad15mNLTEMergeHR

  !----------------------------------------------------------------------------

  subroutine rad15mNLTECalckMin( &
    & xyz_Press, &
    & xy_kMin    &
    & )

    real(DP), intent(in ) :: xyz_Press(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(out) :: xy_kMin  (0:imax-1, 1:jmax)

    integer :: i
    integer :: j
    integer :: k

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_15m_NLTE_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    ! set minimum index used for calculation of
    ! atmospheric radiative cooling rate
!    imin=im-25

    do j = 1, jmax
      do i = 0, imax-1
        find_kmin : do k = kmax, 1, -1
          if ( xyz_Press(i,j,k) > 1.0d-2 ) exit find_kmin
        end do find_kmin
        xy_kMin(i,j) = k
      end do
    end do

  end subroutine rad15mNLTECalckMin

  !----------------------------------------------------------------------------

  subroutine rad15mNLTECalcWeight( &
    & xyz_Press,                   &
!!$    & xy_kMin,          &
    & xyz_Weight                   &
    & )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    real(DP), intent(in ) :: xyz_Press (0:imax-1, 1:jmax, 1:kmax)
!!$    integer , intent(in ) :: xy_kMin   (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyz_Weight(0:imax-1, 1:jmax, 1:kmax)


    ! Local variables
    !
    integer :: i
    integer :: j
    integer :: k

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_15m_NLTE_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1

!               weight(i,j,k)=(atan(2.0d0 &
!                    *dlog(dsqrt(press(i,j,k)*press(i+1,j,k)) &
!                    /dsqrt(press(imin+4,j,k)*press(imin+4+1,j,k)))) &
!                    *1.2d0 &
!                    +pi/2.0d0)/pi

          xyz_Weight(i,j,k) = &
            & ( atan( 2.0d0 * log( xyz_Press(i,j,k) / ( 1.0d-2 * exp( 2.0d0 ) ) ) ) &
            &   * 1.2d0 + Pi / 2.0d0 ) &
            & / Pi
          xyz_Weight(i,j,k) = max( xyz_Weight(i,j,k), 0.0d0 )
          xyz_Weight(i,j,k) = min( xyz_Weight(i,j,k), 1.0d0 )
        end do
      end do
    end do


  end subroutine rad15mNLTECalcWeight

  !----------------------------------------------------------------------------

  subroutine Rad15mNLTEInit

    ! ファイル入出力補助
    ! File I/O support
    !
!!$    use dc_iounit, only: FileOpen

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
!!$    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & AxnameX, &
      & AxnameY, &
      & AxnameZ, &
      & AxnameR, &
      & AxnameT

    ! 宣言文 ; Declaration statements
    !

!!$    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
!!$                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
!!$                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
!!$    namelist /rad_15m_NLTE_nml/ &
!!$      & SolarConst
          !
          ! デフォルト値については初期化手続 "rad_15m_NLTE#Rad15mNLTEInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_15m_NLTE#Rad15mNLTEInit" for the default values.
          !

    if ( rad_15m_NLTE_inited ) return

    ! デフォルト値の設定
    ! Default values settings
    !
!!$    SolarConst      = 1380.0_DP / 1.52_DP**2

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
!!$    if ( trim(namelist_filename) /= '' ) then
!!$      call FileOpen( unit_nml, &          ! (out)
!!$        & namelist_filename, mode = 'r' ) ! (in)
!!$
!!$      rewind( unit_nml )
!!$      read( unit_nml,                     & ! (in)
!!$        & nml = rad_15m_NLTE_nml,         & ! (out)
!!$        & iostat = iostat_nml )             ! (out)
!!$      close( unit_nml )
!!$
!!$      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!!$    end if


    ! Initialization of modules used in this module
    !


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DTempDt15mNLTE',     &
      & (/ AxnameX, AxnameY, AxnameZ, AxnameT /),      &
      & 'radiative calculation by NLTE process at 15 micron meter', 'K s-1' )
    call HistoryAutoAddVariable( 'DTempDtRadLMerged',  &
      & (/ AxnameX, AxnameY, AxnameZ, AxnameT /),      &
      & 'radiative calculation in long wave merged with NLTE heating rate', 'K s-1' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
!!$    call MessageNotify( 'M', module_name, 'SolarConst = %f', d = (/ SolarConst /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_15m_NLTE_inited = .true.

  end subroutine Rad15mNLTEInit

!!$    !**************************************************************************
!!$    subroutine readnirfac(fn,nirfn,nirfp,nirfac)
!!$
!!$      use vtype
!!$
!!$      implicit none
!!$
!!$      character(len=128) :: fn
!!$      integer(i4b) :: nirfn
!!$      real(dp) :: nirfp(nirfn)
!!$      real(dp) :: nirfac(nirfn)
!!$
!!$
!!$      ! Local variables
!!$      !
!!$      integer(i4b) :: i
!!$      character(len=128) :: tmpl
!!$
!!$
!!$      open(70,file='./'//fn,status='unknown')
!!$      read(70,'(a128)') tmpl
!!$      do i=1,nirfn
!!$         read(70,*) nirfp(i),nirfac(i)
!!$         nirfp(i)=nirfp(i)*1.0d-6*1.0d5
!!$      enddo
!!$      close(70)
!!$
!!$      return
!!$    end subroutine readnirfac
!!$
!!$    !**************************************************************************
!!$    !     subroutine nirhrcorrect
!!$    !     correct near infrared heating rate
!!$    !**************************************************************************
!!$    ! im       : number of the vertical layers
!!$    ! jm       : number of the meridional grids
!!$    ! km       : number of the zonal grids
!!$    ! press    : pressure at midpoint of layer (Pa)
!!$    ! nirfn    : maximum number of near infrared heating rate
!!$    !          : factor
!!$    ! nirfp    : pressure for table of near infrared heating rate
!!$    !          : correct factor
!!$    ! nirfa    : near infrared heating rate correct factor
!!$    ! qnir     : near infrared heating rate
!!$    ! corsw    : artificial correction switch, 
!!$    !          : if sw is equal to 1 correction is down
!!$    !**************************************************************************
!!$
!!$    subroutine nirhrcorrect(im,jm,km,press,nirfn,nirfp,nirfac,qnir, &
!!$         corsw)
!!$
!!$      use vtype
!!$
!!$      implicit none
!!$
!!$      integer(i4b) :: im,jm,km
!!$      real(dp) :: press(im+1,jm,km)
!!$      integer(i4b) :: nirfn
!!$      real(dp) :: nirfp(nirfn)
!!$      real(dp) :: nirfac(nirfn)
!!$      real(dp) :: qnir(im,jm,km)
!!$      integer(i4b) :: corsw
!!$
!!$
!!$      ! Local variables
!!$      !
!!$      integer(i4b) :: i,j,k,l
!!$      real(dp) :: tmpp
!!$      real(dp) :: tmpfac
!!$
!!$
!!$      do k=1,km
!!$         do j=1,jm
!!$            do i=1,im
!!$               tmpp=dsqrt(press(i,j,k)*press(i+1,j,k))
!!$               if(tmpp.lt.nirfp(nirfn)) then
!!$
!!$!                  qnir(i,j,k)=qnir(i,j,k) &
!!$!                       /(nirfac(nirfn)*nirfp(nirfn)/tmpp)
!!$
!!$                  tmpfac=nirfac(nirfn) &
!!$                       *(nirfp(nirfn)/tmpp)*(nirfp(nirfn)/tmpp) &
!!$                       *(nirfp(nirfn)/tmpp)*(nirfp(nirfn)/tmpp)
!!$                  tmpfac=1.0d0/tmpfac
!!$               else if(tmpp.le.nirfp(1)) then
!!$                  do l=2,nirfn
!!$                     if(tmpp.gt.nirfp(l)) go to 100
!!$                  enddo
!!$100               continue
!!$                  tmpfac=(nirfac(l)-nirfac(l-1))/(nirfp(l)-nirfp(l-1)) &
!!$                       *(tmpp-nirfp(l-1))+nirfac(l-1)
!!$                  tmpfac=1.0d0/tmpfac
!!$                  if(tmpfac.gt.1.0d0) then
!!$                     write(*,*) 'Factor is greater than 1'
!!$                     write(*,*) 'Stop'
!!$                     stop
!!$                  endif
!!$               endif
!!$
!!$               ! This is artificial correction for 2.7 micron band
!!$               !
!!$               if((corsw.eq.1).and.(tmpp.lt.(1.2d-8*1.0d5))) then
!!$!                  qnir(i,j,k)=qnir(i,j,k)*dsqrt(tmpp/(1.2d-8*1.0d5))
!!$                  tmpfac=tmpfac*dsqrt(tmpp/(1.2d-8*1.0d5))
!!$               endif
!!$
!!$               qnir(i,j,k)=qnir(i,j,k)*tmpfac
!!$
!!$            enddo
!!$         enddo
!!$      enddo
!!$
!!$      return
!!$    end subroutine nirhrcorrect

end module rad_15m_NLTE