program main

  use vtype_module
  use ni3_module
  use fi_module

  implicit none

  character(extstr)     :: CtlFN = "lsf.conf"
  integer               :: FileUnit
  integer               :: iostat

  integer               :: NLev
  real(dp), allocatable :: p_SWUpFluxObs(:)
  real(dp), allocatable :: p_SWDnFluxObs(:)
  real(dp)              :: MissingValue

  integer , parameter   :: NParam = 4
  real(dp), allocatable :: pa_SWUpFluxModelPro(:,:)
  real(dp), allocatable :: pa_SWDnFluxModelPro(:,:)
  real(dp)              :: a_SWUpFluxModelTOA( 0:NParam )
  real(dp)              :: a_SWDnFluxModelTOA( 0:NParam )
  real(dp)              :: a_LWUpFluxModelTOA( 0:NParam )
  real(dp)              :: a_LWDnFluxModelTOA( 0:NParam )

  real(dp), allocatable :: pa_DSWUpFluxModelDParam(:,:)
  real(dp), allocatable :: pa_DSWDnFluxModelDParam(:,:)
  real(dp)              :: a_DSWUpFluxModelTOADParam( NParam )
  real(dp)              :: a_DSWDnFluxModelTOADParam( NParam )
  real(dp)              :: a_DLWUpFluxModelTOADParam( NParam )
  real(dp)              :: a_DLWDnFluxModelTOADParam( NParam )

  character(2*extstr)   :: NcFNObs
  character(2*extstr)   :: a_NcFNModelPro( 0:NParam )
  character(2*extstr)   :: a_NcFNModelTOA( 0:NParam )

  real(dp)              :: aa_Param( 0:1, NParam )

  character(2*extstr)   :: NcFNModelPro( 0:NParam )
  character(2*extstr)   :: NcFNModelTOA( 0:NParam )
  real(dp)              :: Param( 0:1, NParam )

  real(dp)              :: Gamma
  real(dp)              :: aa_Mat( NParam, NParam )
  real(dp)              :: a_Vec ( NParam )

  real(dp)              :: SWUpFlux
  real(dp)              :: SWDnFlux
  real(dp)              :: LWUpFlux
  real(dp)              :: LWDnFlux

  integer               :: NcID
  character(extstr)     :: Mode
  character(extstr)     :: Name
  character(extstr)     :: AttName

  integer               :: j
  integer               :: ks
  integer               :: ke
  integer               :: k
  integer               :: n

  namelist /file_nml/ &
    & NcFNModelPro, NcFNModelTOA
  namelist /param_nml/ &
    & Param


  NcFNObs        = "data/Venus_SWFluxProfile_Tomasko-etal1980.nc"

  NcFNModelPro = ""
  NcFNModelTOA = ""

  Param = 0.0d0


  ! Read a namelist file

  call fi_open( CtlFN, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = file_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in file_nml is something wrong.'
  end if
  write( 6, nml = file_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = param_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in param_nml is something wrong.'
  end if
  write( 6, nml = param_nml )
  !
  close( FileUnit )


  a_NcFNModelPro = NcFNModelPro
  a_NcFNModelTOA = NcFNModelTOA
  aa_Param       = Param


  Mode = "read"
  call ni3_open( NcFNObs, Mode, NcID )
  Name = "plev"
  call ni3_inq_dimlen( NcID, Name, NLev )
  allocate( p_SWUpFluxObs( NLev ) )
  allocate( p_SWDnFluxObs( NLev ) )
  Mode = "read"
  call ni3_open( NcFNObs, Mode, NcID )
  Name = "upflux"
  call ni3_get_var( NcID, Name, p_SWUpFluxObs )
  Name = "dnflux"
  call ni3_get_var( NcID, Name, p_SWDnFluxObs )
  AttName = "missing_value"
  call ni3_get_att( NcID, Name, AttName, MissingValue )
  call ni3_close( NcID )

  do k = 1, NLev
    if ( p_SWUpFluxObs(k) /= MissingValue ) exit
  end do
  ks = k
  do k = ks, NLev
    if ( p_SWUpFluxObs(k) == MissingValue ) exit
  end do
  ke = k-1


  allocate( pa_SWUpFluxModelPro( NLev, 0:NParam ) )
  allocate( pa_SWDnFluxModelPro( NLev, 0:NParam ) )


  do n = 0, NParam
    call readncfile( &
      & a_NcFNModelPro(n), a_NcFNModelTOA(n), &
      & NLev, &
      & pa_SWUpFluxModelPro(:,n), pa_SWDnFluxModelPro(:,n), &
      & a_SWUpFluxModelTOA(n), a_SWDnFluxModelTOA(n), &
      & a_LWUpFluxModelTOA(n), a_LWDnFluxModelTOA(n) &
      & )
  end do


  allocate( pa_DSWUpFluxModelDParam( NLev, NParam ) )
  allocate( pa_DSWDnFluxModelDParam( NLev, NParam ) )

  do n = 1, NParam
    do k = 1, ks-1
      pa_DSWUpFluxModelDParam( k, n ) = MissingValue
      pa_DSWDnFluxModelDParam( k, n ) = MissingValue
    end do
    do k = ks, ke
      pa_DSWUpFluxModelDParam( k, n ) = &
        &   ( pa_SWUpFluxModelPro( k, n ) - pa_SWUpFluxModelPro( k, 0 ) ) &
        & / ( aa_Param           ( 1, n ) - aa_Param           ( 0, n ) )
      pa_DSWDnFluxModelDParam( k, n ) = &
        &   ( pa_SWDnFluxModelPro( k, n ) - pa_SWDnFluxModelPro( k, 0 ) ) &
        & / ( aa_Param           ( 1, n ) - aa_Param           ( 0, n ) )
    end do
    do k = ke+1, NLev
      pa_DSWUpFluxModelDParam( k, n ) = MissingValue
      pa_DSWDnFluxModelDParam( k, n ) = MissingValue
    end do
    !
    a_DSWUpFluxModelTOADParam( n ) = &
        &   ( a_SWUpFluxModelTOA( n )    - a_SWUpFluxModelTOA( 0 )    ) &
        & / ( aa_Param          ( 1, n ) - aa_Param          ( 0, n ) )
    a_DSWDnFluxModelTOADParam( n ) = &
        &   ( a_SWDnFluxModelTOA( n )    - a_SWDnFluxModelTOA( 0 )    ) &
        & / ( aa_Param          ( 1, n ) - aa_Param          ( 0, n ) )
    a_DLWUpFluxModelTOADParam( n ) = &
        &   ( a_LWUpFluxModelTOA( n )    - a_LWUpFluxModelTOA( 0 )    ) &
        & / ( aa_Param          ( 1, n ) - aa_Param          ( 0, n ) )
    a_DLWDnFluxModelTOADParam( n ) = &
        &   ( a_LWDnFluxModelTOA( n )    - a_LWDnFluxModelTOA( 0 )    ) &
        & / ( aa_Param          ( 1, n ) - aa_Param          ( 0, n ) )
  end do



!!$  do n = NParam, NParam
!!$    do k = 1, ks-1
!!$      pa_DSWUpFluxModelDParam( k, n ) = MissingValue
!!$      pa_DSWDnFluxModelDParam( k, n ) = MissingValue
!!$    end do
!!$    do k = ks, ke
!!$      pa_DSWUpFluxModelDParam( k, n ) = 0.0d0
!!$      pa_DSWDnFluxModelDParam( k, n ) = 0.0d0
!!$    end do
!!$    do k = ke+1, NLev
!!$      pa_DSWUpFluxModelDParam( k, n ) = MissingValue
!!$      pa_DSWDnFluxModelDParam( k, n ) = MissingValue
!!$    end do
!!$    !
!!$    a_DSWUpFluxModelTOADParam( n ) = 0.0d0
!!$    a_DSWDnFluxModelTOADParam( n ) = 0.0d0
!!$    a_DLWUpFluxModelTOADParam( n ) = 0.0d0
!!$    a_DLWDnFluxModelTOADParam( n ) = 0.0d0
!!$  end do



  Gamma = ke - ks + 1
!!$  Gamma = Gamma / 4.0
!!$  Gamma = 20.0
!!$  Gamma = 10.0
  Gamma =  5.0
!!$  Gamma =  1.0
!!$  Gamma = 0
  write( 6, * ) 'Gamma = ', Gamma
  ! row
  do n = 1, NParam

    ! column
    do j = 1, NParam
      aa_Mat(n,j) = 0.0d0
      do k = ks, ke
        aa_Mat(n,j) = aa_Mat(n,j) &
          & + pa_DSWUpFluxModelDParam(k,n) * pa_DSWUpFluxModelDParam(k,j) &
          & + pa_DSWDnFluxModelDParam(k,n) * pa_DSWDnFluxModelDParam(k,j)
      end do
      aa_Mat(n,j) = aa_Mat(n,j) &
        & + Gamma &
        &   * ( &
        &         ( a_DSWUpFluxModelTOADParam(j) - a_DSWDnFluxModelTOADParam(j) ) &
        &       + ( a_DLWUpFluxModelTOADParam(j) - a_DLWDnFluxModelTOADParam(j) ) &
        &     ) &
        &   * ( &
        &         ( a_DSWUpFluxModelTOADParam(n) - a_DSWDnFluxModelTOADParam(n) ) &
        &       + ( a_DLWUpFluxModelTOADParam(n) - a_DLWDnFluxModelTOADParam(n) ) &
        &     )
    end do

    a_Vec(n) = 0.0d0
    do k = ks, ke
      a_Vec(n) = a_Vec(n) &
        & + ( pa_SWUpFluxModelPro(k,0) - p_SWUpFluxObs(k) ) &
        &     * pa_DSWUpFluxModelDParam(k,n)                &
        & + ( pa_SWDnFluxModelPro(k,0) - p_SWDnFluxObs(k) ) &
        &     * pa_DSWDnFluxModelDParam(k,n)
    end do
    a_Vec(n) = a_Vec(n) &
      & + Gamma &
      &   * (                                                     &
      &       + ( a_SWUpFluxModelTOA(0) - a_SWDnFluxModelTOA(0) ) &
      &       + ( a_LWUpFluxModelTOA(0) - a_LWDnFluxModelTOA(0) ) &
      &     )                                                     &
      &   * ( &
      &         ( a_DSWUpFluxModelTOADParam(n) - a_DSWDnFluxModelTOADParam(n) ) &
      &       + ( a_DLWUpFluxModelTOADParam(n) - a_DLWDnFluxModelTOADParam(n) ) &
      &     )
    a_Vec(n) = - a_Vec(n)

  end do


  ! No adjustment for k_{abs}
!!$  n = NParam
!!$  do j = 1, NParam
!!$    aa_Mat(n,j) = 0.0d0
!!$  end do
!!$  aa_Mat(n,n) = 1.0d0
!!$  a_Vec(n) = aa_Param(0,n) - aa_Param(0,n)

  ! No adjustment for cloud
!!$  do n = 1, 3
!!$    do j = 1, NParam
!!$      aa_Mat(n,j) = 0.0d0
!!$    end do
!!$    aa_Mat(n,n) = 1.0d0
!!$    a_Vec(n) = aa_Param(0,n) - aa_Param(0,n)
!!$  end do


  call ludecomp_prep_simple( NParam, aa_Mat )
  call ludecomp_solve_simple( NParam, aa_Mat, a_Vec )

  write( 6, * ) 'Best fit parameters:'
  write( 6, * ) '  cloud mode 1 & UV absorber: alpha1 = ', a_Vec(1) + aa_Param(0,1)
  write( 6, * ) '  cloud mode 2 & 2 prime    : alpha2 = ', a_Vec(2) + aa_Param(0,2)
  write( 6, * ) '  cloud mode 3              : alpha3 = ', a_Vec(3) + aa_Param(0,3)
  write( 6, * ) '  k_abs                     : alpha4 = ', a_Vec(4) + aa_Param(0,4)

  do k = ks, ke

    SWUpFlux = pa_SWUpFluxModelPro(k,0)
    SWDnFlux = pa_SWDnFluxModelPro(k,0)
    do n = 1, NParam
      SWUpFlux = SWUpFlux + pa_DSWUpFluxModelDParam(k,n) * a_Vec(n)
      SWDnFlux = SWDnFlux + pa_DSWDnFluxModelDParam(k,n) * a_Vec(n)
    end do

    write( 60, * ) k, SWUpFlux, p_SWUpFluxObs(k), SWDnFlux, p_SWDnFluxObs(k)
  end do

  SWUpFlux = a_SWUpFluxModelTOA(0)
  SWDnFlux = a_SWDnFluxModelTOA(0)
  LWUpFlux = a_LWUpFluxModelTOA(0)
  LWDnFlux = a_LWDnFluxModelTOA(0)
  do n = 1, NParam
    SWUpFlux = SWUpFlux + a_DSWUpFluxModelTOADParam(n) * a_Vec(n)
    SWDnFlux = SWDnFlux + a_DSWDnFluxModelTOADParam(n) * a_Vec(n)
    LWUpFlux = LWUpFlux + a_DLWUpFluxModelTOADParam(n) * a_Vec(n)
    LWDnFlux = LWDnFlux + a_DLWDnFluxModelTOADParam(n) * a_Vec(n)
  end do

  write( 6, * ) 0, SWUpFlux, SWDnFlux, LWUpFlux, LWDnFlux, (SWUpFlux-SWDnFlux)+(LWUpFlux-LWDnFlux)


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

contains

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

  subroutine readncfile( &
    & NcFNModelPro, NcFNModelTOA, &
    & NLev, &
    & p_SWUpFluxModelPro, p_SWDnFluxModelPro, &
    & SWUpFluxModelTOA, SWDnFluxModelTOA, LWUpFluxModelTOA, LWDnFluxModelTOA &
    & )

    use ni3_module

    character(*), intent(in ) :: NcFNModelPro
    character(*), intent(in ) :: NcFNModelTOA
    integer     , intent(in ) :: NLev
    real(dp)    , intent(out) :: p_SWUpFluxModelPro( NLev )
    real(dp)    , intent(out) :: p_SWDnFluxModelPro( NLev )
    real(dp)    , intent(out) :: SWUpFluxModelTOA
    real(dp)    , intent(out) :: SWDnFluxModelTOA
    real(dp)    , intent(out) :: LWUpFluxModelTOA
    real(dp)    , intent(out) :: LWDnFluxModelTOA


    ! Local variables
    !
    integer               :: NcID
    character(extstr)     :: Mode
    character(extstr)     :: Name

    real(dp)              :: p_Flux( NLev )



    Mode = "read"
    call ni3_open( NcFNModelPro, Mode, NcID )
    Name = "RadUwFluxSW"
    call ni3_get_var( NcID, Name, p_Flux )
    p_SWUpFluxModelPro = p_Flux
    Name = "RadDwFluxSW"
    call ni3_get_var( NcID, Name, p_Flux )
    p_SWDnFluxModelPro = p_Flux
    call ni3_close( NcID )

    Mode = "read"
    call ni3_open( NcFNModelTOA, Mode, NcID )
    Name = "RadUwFluxSW"
    call ni3_get_var( NcID, Name, p_Flux )
    SWUpFluxModelTOA = p_Flux(NLev)
    Name = "RadDwFluxSW"
    call ni3_get_var( NcID, Name, p_Flux )
    SWDnFluxModelTOA = p_Flux(NLev)
    Name = "RadUwFluxLW"
    call ni3_get_var( NcID, Name, p_Flux )
    LWUpFluxModelTOA = p_Flux(NLev)
    Name = "RadDwFluxLW"
    call ni3_get_var( NcID, Name, p_Flux )
    LWDnFluxModelTOA = p_Flux(NLev)
    call ni3_close( NcID )


  end subroutine readncfile

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

end program main
