
module set_1d_profile

  use vtype_module

  implicit none

  private



  public :: Set1DProfileAtm
  public :: Set1DProfileTemp


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

contains

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

  subroutine Set1DProfileTemp(      &
    & InFileName, PressName, TempName, &
    & kmax, z_Sigma,       &
    & Ps, z_Press, SurfTemp, z_Temp  &
    & )

    use ni3_module
    use interpolate, only: Interpolate1D

    character(*), intent(in) :: InFileName
    character(*), intent(in) :: PressName
    character(*), intent(in) :: TempName

    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: z_Sigma  (1:kmax)
    real(DP), intent(out) :: Ps
    real(DP), intent(out) :: z_Press  (1:kmax)
    real(DP), intent(out) :: SurfTemp
    real(DP), intent(out) :: z_Temp   (1:kmax)


    !
    ! local variables
    !
    integer               :: NcID

    integer               :: Inkmax
    real(DP), allocatable :: a_InPress(:)
    real(DP), allocatable :: a_InTemp (:)

    real(DP), allocatable :: a_Tmp(:)

    integer               :: k


    call ni3_open( InFileName, "read", NcID )

    call ni3_inq_dimlen( NcID, PressName, Inkmax )

    allocate( a_InPress   (Inkmax) )
    allocate( a_InTemp    (Inkmax) )
    allocate( a_Tmp       (Inkmax) )

    call ni3_get_var( NcID, PressName , a_InPress )
    call ni3_get_var( NcID, TempName  , a_InTemp  )

    call ni3_close( NcID )


!!$    Ps       = a_InPress(1)
!!$    SurfTemp = a_InTemp (1)
    if ( a_InPress(1) < a_InPress(Inkmax) ) then
      a_Tmp = a_InPress
      do k = 1, Inkmax
        a_InPress(k) = a_Tmp(Inkmax-k+1)
      end do
      a_Tmp = a_InTemp
      do k = 1, Inkmax
        a_InTemp(k) = a_Tmp(Inkmax-k+1)
      end do
    end if
    Ps       = a_InPress(1)
    SurfTemp = a_InTemp (1)

    z_Press = Ps * z_Sigma


    call Interpolate1D(       &
      & Inkmax, a_InPress, a_InTemp,    &
      & kmax, z_Press,                  &
      & z_Temp                          &
      & )


    deallocate( a_InPress    )
    deallocate( a_InTemp     )
    deallocate( a_Tmp        )


  end subroutine Set1DProfileTemp

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

  subroutine Set1DProfileAtm(      &
    & InFileName, PressName, TempName, H2OVapName, O3Name, &
    & kmax, z_Sigma,       &
    & Ps, z_Press, SurfTemp, z_Temp, z_QH2OVap, z_QO3     &
    & )

    use ni3_module
    use interpolate, only: Interpolate1D

    character(*), intent(in) :: InFileName
    character(*), intent(in) :: PressName
    character(*), intent(in) :: TempName
    character(*), intent(in) :: H2OVapName
    character(*), intent(in) :: O3Name

    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: z_Sigma  (1:kmax)
    real(DP), intent(out) :: Ps
    real(DP), intent(out) :: z_Press  (1:kmax)
    real(DP), intent(out) :: SurfTemp
    real(DP), intent(out) :: z_Temp   (1:kmax)
    real(DP), intent(out) :: z_QH2OVap(1:kmax)
    real(DP), intent(out) :: z_QO3    (1:kmax)


    !
    ! local variables
    !
    integer               :: NcID

    integer               :: Inkmax
    real(DP), allocatable :: a_InPress(:)
    real(DP), allocatable :: a_InTemp (:)
    real(DP), allocatable :: a_InQH2O (:)
    real(DP), allocatable :: a_InQO3  (:)

    real(DP), allocatable :: a_InLogValue(:)



    call ni3_open( InFileName, "read", NcID )

    call ni3_inq_dimlen( NcID, PressName, Inkmax )

    allocate( a_InPress   (Inkmax) )
    allocate( a_InTemp    (Inkmax) )
    allocate( a_InQH2O    (Inkmax) )
    allocate( a_InQO3     (Inkmax) )
    allocate( a_InLogValue(Inkmax) )

    call ni3_get_var( NcID, PressName , a_InPress )
    call ni3_get_var( NcID, TempName  , a_InTemp  )
    call ni3_get_var( NcID, H2OVapName, a_InQH2O  )
    call ni3_get_var( NcID, O3Name    , a_InQO3   )

    call ni3_close( NcID )


    Ps       = a_InPress(1)
    SurfTemp = a_InTemp (1)

    z_Press = Ps * z_Sigma


    call Interpolate1D(       &
      & Inkmax, a_InPress, a_InTemp,    &
      & kmax, z_Press,                  &
      & z_Temp                          &
      & )


    if ( any( a_InQH2O <= 0.0_DP ) ) then
      stop 'QH2O contains values <= 0.'
    end if
    a_InLogValue = log( a_InQH2O )
    !
    call Interpolate1D(        &
      & Inkmax, a_InPress, a_InLogValue, &
      & kmax, z_Press,                   &
      & z_QH2OVap                        &
      & )
    !
    z_QH2OVap = exp( z_QH2OVap )

    if ( any( a_InQO3 <= 0.0_DP ) ) then
      stop 'QO3 contains values <= 0.'
    end if
    a_InLogValue = log( a_InQO3 )
    !
    call Interpolate1D(        &
      & Inkmax, a_InPress, a_InLogValue, &
      & kmax, z_Press,                   &
      & z_QO3                            &
      & )
    !
    z_QO3 = exp( z_QO3 )


    deallocate( a_InPress    )
    deallocate( a_InTemp     )
    deallocate( a_InQH2O     )
    deallocate( a_InQO3      )
    deallocate( a_InLogValue )


  end subroutine Set1DProfileAtm

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

  subroutine Set1DProfileInterpolate( &
    & NLev, a_Press, a_Array,         &
    & kmax, z_Press,                  &
    & z_Array                         &
    & )

    integer , intent(in ) :: NLev
    real(DP), intent(in ) :: a_Press  (1:NLev)
    real(DP), intent(in ) :: a_Array  (1:NLev)
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: z_Press(1:kmax)
    real(DP), intent(out) :: z_Array(1:kmax)


    !
    ! local variables
    !
    integer :: k
    integer :: kk


    do k = 1, kmax

      if( z_Press(k) <= a_Press(NLev) ) then
        z_Array(k) = a_Array(NLev)
      else
        search_loop : do kk = 2, NLev
          if( a_Press( kk ) < z_Press(k) ) exit search_loop
        end do search_loop
        if( kk > NLev ) &
          stop 'Unexpected error in setting temperature profile'
        z_Array(k) =                                     &
          &   ( a_Array( kk ) - a_Array( kk-1 ) )        &
          & / ( log( a_Press(kk  ) / a_Press( kk-1 ) ) ) &
          & * ( log( z_Press(k   ) / a_Press( kk-1 ) ) ) &
          & + a_Array( kk-1 )
      end if

    end do


  end subroutine Set1DProfileInterpolate

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

end module set_1d_profile
