program main

  use vtype_module
  use ni3_module
  use ca_module
  use netcdf
  use saturate

  implicit none

  interface
    subroutine findfu( fn, ios, fu, mode )
      use vtype_module
      implicit none
      character(len=*), intent(in )           :: fn
      integer(i4b)    , intent(out)           :: ios, fu
      character(len=*), intent(in ), optional :: mode
    end subroutine findfu
  end interface

  character(len=extstr)            :: mode

  character(len=extstr)            :: ncfn_q
  character(len=extstr)            :: ncfn_t
  character(len=extstr)            :: ncfn_out
  integer                          :: ncid_q
  integer                          :: ncid_t
  integer                          :: ncid_out
  character(len=extstr)            :: varname_q
  character(len=extstr)            :: varname_t
  character(len=extstr)            :: varname_out

  character(len=extstr)            :: name
!!$  character(len=extstr)            :: stdname
!!$  character(len=extstr)            :: units

  integer(i4b)                       :: NDims
  character(len=extstr), allocatable :: a_DimNames(:)

  integer                          :: imax
  integer                          :: jmax
  integer                          :: kmax
  integer                          :: tmax

  real(dp)                         :: FillValue
  real(dp)                         :: Temp1
  real(dp)                         :: Temp2
  character(len=256)               :: MixPhaseType

  real(dp)           , allocatable :: x_Lon(:)
  real(dp)           , allocatable :: y_Lat(:)
  real(dp)           , allocatable :: z_VLev(:)
  real(dp)           , allocatable :: xyz_Q (:,:,:)
  real(dp)           , allocatable :: xyz_T (:,:,:)
  real(dp)           , allocatable :: xyz_P (:,:,:)
  real(dp)           , allocatable :: xyz_RH(:,:,:)
  real(dp)           , allocatable :: a_Time(:)


  logical                          :: ex
  integer                          :: ios

!!$  integer                          :: DayStart
!!$  integer                          :: DayEnd

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


  character(extstr)                :: fn

  character(extstr)                :: ctlfn = "calcrh.cntl"
  integer                          :: ctlfu


  namelist /file/  ncfn_q, varname_q, ncfn_t, varname_t, ncfn_out
  namelist /qs/ Temp1, Temp2, MixPhaseType


  ncfn_q    = "q"
  varname_q = "q"
  ncfn_t    = "t"
  varname_t = "t"
  ncfn_out  = "out.nc"

  Temp1 = 0.0
  Temp2 = 0.0
  MixPhaseType = 'Lin'

  fn = ctlfn
  inquire( file = fn, exist = ex )
  if( .not. ex ) then
    write( 6, * ) "Control file ", trim( ctlfn ), " does not exit."
    write( 6, * ) "Control file should contain namelists listed below."
    write( 6, * ) "    &input"
    write( 6, * ) "        dir       = '.'"
    write( 6, * ) "    &end"
    write( 6, * ) "    &lander"
    write( 6, * ) "        LanderSym = 'VL1'"
    write( 6, * ) "    &end"
    write( 6, * ) "    &output"
    write( 6, * ) "        outfn     = 'vl1_dmps.dat'"
    write( 6, * ) "    &end"
    write( 6, * ) "    &dayinfo"
    write( 6, * ) "        DayStart  =     1"
    write( 6, * ) "        DayEnd    = 10000"
    write( 6, * ) "    &end"
    stop
  end if
  call findfu( fn, ios, ctlfu )
  if( ios /= 0 ) then
    write( 6, * ) 'STOP: ', ios
    stop
  end if
  open( ctlfu, file = fn, status = 'unknown' )

  rewind( ctlfu )

  loop_namelist : do

    read( ctlfu, nml = file, iostat = ios )
    write( 6, nml = file )
    if ( ios /= 0 ) exit
    rewind( ctlfu )
    read( ctlfu, nml = qs  , iostat = ios )
    write( 6, nml = qs )
    if ( ios /= 0 ) exit
    write( 6, * ) 'Read   ', trim( ncfn_q   )
    write( 6, * ) 'Read   ', trim( ncfn_t   )
    write( 6, * ) 'Output ', trim( ncfn_out )


    !----------------------------------------------------------
    !
    ! Preparing for ps file
    !
    mode = 'read'
    call ni3_open( ncfn_q, mode, ncid_q )
    mode = 'read'
    call ni3_open( ncfn_t, mode, ncid_t )
    !
    !----------------------------------------------------------

    !
    ! Number of dimensions are checked.
    !
    call ni3_inq_var( ncid_q, varname_q, ndims = NDims )
    allocate( a_DimNames( NDims ) )
    call ni3_inq_vardimnames( ncid_q, varname_q, NDims, a_DimNames )

    call ni3_inq_dimlen( ncid_q, a_DimNames(1)    , imax )
    call ni3_inq_dimlen( ncid_q, a_DimNames(2)    , jmax )
    call ni3_inq_dimlen( ncid_q, a_DimNames(NDims), tmax )
    allocate( x_Lon (0:imax-1) )
    allocate( y_Lat (1:jmax  ) )
    allocate( a_Time(1:tmax  ) )
    call ni3_get_var( ncid_q, a_DimNames(1)    , x_Lon  )
    call ni3_get_var( ncid_q, a_DimNames(2)    , y_Lat  )
    call ni3_get_var( ncid_q, a_DimNames(NDims), a_Time )
    call ni3_inq_dimlen( ncid_q, a_DimNames(3), kmax )
    allocate( z_VLev(1:kmax) )
    call ni3_get_var( ncid_q, a_DimNames(3), z_VLev )

    allocate( xyz_Q (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_T (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_P (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_RH(0:imax-1, 1:jmax, 1:kmax) )


    !----------------------------------------------------------
    !
    ! Preparing for output file
    !
    mode = 'new'
    call ni3_open( ncfn_out, mode, ncid_out )

    name = a_DimNames(1)
    call ni3_set_dim( ncid_out, name, NF90_REAL, x_Lon )
    call ni3_cp_atts( ncid_q, ncid_out, name )

    name = a_DimNames(2)
    call ni3_set_dim( ncid_out, name, NF90_REAL, y_Lat )
    call ni3_cp_atts( ncid_q, ncid_out, name )

    name = a_DimNames(3)
    call ni3_set_dim( ncid_out, name, NF90_REAL, z_VLev )
    call ni3_cp_atts( ncid_q, ncid_out, name )

    name = a_DimNames(NDims)
    call ni3_def_dim( ncid_out, name, NF90_REAL, NF90_UNLIMITED )
    call ni3_cp_atts( ncid_q, ncid_out, name )

!!$  name    = "time"
!!$  stdname = "time"
!!$  units   = "days since 0000-01-01"
!!$  call ni3_def_dim( ncid_out, name, NF90_REAL, NF90_UNLIMITED, stdname = stdname, units  = units )
!!$  call ni3_cp_atts( ncid_ice, ncid_out, name )


    varname_out = 'RH'
    name = varname_out
    call ni3_def_var( ncid_out, name, NF90_REAL, NDims, a_DimNames, &
      & longname = 'relative humidity', units = '1' )
!!$    call ni3_cp_atts( ncid_q, ncid_out, name )

    if ( ni3_chk_att( ncid_q, varname_q, "missing_value" ) ) then
      call ni3_get_att( ncid_q  , varname_q  , "missing_value", FillValue )
      call ni3_put_att( ncid_out, varname_out, "missing_value", real(FillValue) )
      call ni3_put_att( ncid_out, varname_out, "_FillValue"   , real(FillValue) )
    else if ( ni3_chk_att( ncid_q, varname_q, "_FillValue" ) ) then
      call ni3_get_att( ncid_q  , varname_q  , "_FillValue"   , FillValue )
      call ni3_put_att( ncid_out, varname_out, "missing_value", real(FillValue) )
      call ni3_put_att( ncid_out, varname_out, "_FillValue"   , real(FillValue) )
    else
      FillValue = -999.0
      call ni3_put_att( ncid_out, varname_out, "missing_value", real(FillValue) )
      call ni3_put_att( ncid_out, varname_out, "_FillValue"   , real(FillValue) )
    end if
    !
    !----------------------------------------------------------


    do t = 1, tmax

      print *, t, tmax

      call ni3_get_varss( ncid_q, varname_q, t, xyz_Q )
      call ni3_get_varss( ncid_t, varname_t, t, xyz_T )

      do k = 1, kmax
        xyz_P(:,:,k) = z_VLev(k)
      end do
!!$      call xyz_CalcQVapSat( 0, imax-1, 1, jmax, 1, kmax, xyz_T, xyz_P, FillValue, xyz_RH )
      call xyz_CalcQVapSat( Temp1, Temp2, MixPhaseType, xyz_T, xyz_P, FillValue, xyz_RH )

      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_Q(i,j,k) /= FillValue ) then
              xyz_RH(i,j,k) = xyz_Q(i,j,k) / xyz_RH(i,j,k)
            end if
          end do
        end do
      end do

      call ni3_put_varss( ncid_out, "time"     , t, a_Time(t) )
      call ni3_put_varss( ncid_out, varname_out, t, xyz_RH    )

    end do


    deallocate( a_DimNames )
    deallocate( x_Lon )
    deallocate( y_Lat )
    deallocate( a_Time )
    deallocate( z_VLev )
    deallocate( xyz_Q )
    deallocate( xyz_T )
    deallocate( xyz_P )
    deallocate( xyz_RH )


    call ni3_close( ncid_q   )
    call ni3_close( ncid_t   )
    call ni3_close( ncid_out )

  end do loop_namelist

  close( ctlfu )

end program main
