program read_PMR
  use netcdf
  use,intrinsic :: iso_fortran_env
  use,intrinsic :: iso_c_binding

  integer, parameter :: len_char = 1024

  integer :: iarg
  character(len=len_char) :: carg, ofile

  type type_input
    character(len=len_char) :: filename
    character(len=len_char) :: freqname
  end type

  type(type_input) :: input
!
! PMR header + a
!
  type type_headerPMR
    character(len=1) :: OrbitDirection
    integer :: MeanAnomaly, MeanMotion, Eccentricity, PerigeeArgument, AscendingNodeLongitude, OrbitalInclination
    integer :: EpochTime
    integer :: OrbitNumber
    character(len=10) :: startDate
    character(len=12) :: startTime
    character(len=10) :: endDate
    character(len=12) :: endTime
    character(len=24) :: startDateTime
    character(len=24) :: endDateTime
    real :: scLatPoint(4)
    real :: scLonPoint(4)
  end type type_headerPMR

  type(type_headerPMR) :: header
!
! PMR
!
  real, allocatable :: LatitudePMR(:,:,:), LongitudePMR(:,:,:)
  integer(int8), allocatable :: flagPrecipPMR(:,:), dataQualityPMR(:,:)
  integer(int16), allocatable :: LandSurfaceTypePMR(:,:), qualityDataPMR(:,:), dayCountPMR(:)
  integer, allocatable :: msCountPMR(:)
!
! Common
!
  real, allocatable :: Latitude(:,:), Longitude(:,:), ellipsoidBinOffset(:,:), localZenithAngle(:,:)

  integer(int8), allocatable :: snowIceCover(:,:)
  integer(int16), allocatable :: binRealSurface(:,:), binFirstLatlon(:,:), binClutterFreeBottom(:,:), binStormTop(:,:)
  integer, allocatable :: LandSurfaceType(:,:), flagPrecip(:,:), qualityData(:,:)

  integer(int16), allocatable :: Year(:), DayOfYear(:), MilliSecond(:)
  integer(int8), allocatable :: Month(:), DayOfMonth(:), Hour(:), Minute(:), Second(:)
  real(real64), allocatable :: SecondOfDay(:), FractionalGranuleNumber(:)

  real, allocatable :: height(:,:,:), zFactorMeasured(:,:,:), sigmaZeroMeasured(:,:), elevation(:,:), snRatioAtRealSurface(:,:), heightStormTop(:,:)
!
  integer :: nray, nscan, nlev
  integer :: iray, iscan, ilev

  integer :: values(9)
  character(len=24) :: ctimefmt

  interface

    integer(c_int64_t) function mktime_c( values ) bind(c, name="mktime_c")
      import c_int, c_int64_t
      integer(c_int), dimension(9), intent(in) :: values
    end function

    integer(c_int64_t) function mkutctime_c( values ) bind(c, name="mkutctime_c")
      import c_int, c_int64_t
      integer(c_int), dimension(9), intent(in) :: values
    end function

    subroutine sub_gmtime_c( t, values ) bind(c, name="sub_gmtime_c")
      import c_int, c_int64_t
      integer(c_int64_t), value, intent(in) :: t
      integer(c_int), dimension(9), intent(out) :: values
    end subroutine

  end interface

  if(iargc() /= 2) then
    print *, "read_PMR.exe (filename) (freqname)"
    call exit(0)
  endif

  call getarg(1, input%filename) ! filename
  call getarg(2, input%freqname) ! Ku or Ka

  write(0,*) trim(input%filename), " ", trim(input%freqname)

  call init_PMR_L1()

  call read_PMR_L1(trim(input%filename), trim(input%freqname))

  call calc_PMR_L1()

  call free_PMR_L1()

contains

subroutine init_PMR_L1
  implicit none
end subroutine init_PMR_L1

subroutine read_PMR_L1(filename, freqname)
  use hdf5
  implicit none

  integer, parameter :: len_char_nc = 1024

  character(len=*), intent(in) :: filename, freqname
  character(len=len_char_nc) :: varname, grpname, subgrpname, dimname

  integer :: ncid, grpid, dimid, varid, kindid, subgrpid
  integer :: dimids(4), dimsize, iparent
  integer :: dimnum
  integer :: ijk

  integer :: istat

  integer :: ierror
  integer(HID_T) :: file_id, attr_id, data_dims
  INTEGER(HSIZE_T) :: dims(1)

  real(real32)  , parameter :: MISSING_FLOAT  = -9999.9_real32
  real(real64)  , parameter :: MISSING_DOUBLE = -9999.9_real64
  integer(int8) , parameter :: MISSING_BYTE   = -99_int8
  integer(int16), parameter :: MISSING_SHORT  = -9999_int16
  integer(int32), parameter :: MISSING_INT    = -9999_int32

! time
  integer*8 :: imktime
  integer :: idate, isod

! HDF5
  call h5open_f(ierror)

  call h5fopen_f(filename, H5F_ACC_RDONLY_F, file_id, ierror)

  call h5aopen_name_f(file_id, "Orbit Direction", attr_id, ierror)

  dims(:) = 1
  call h5aread_f(attr_id, H5T_NATIVE_CHARACTER, header%OrbitDirection, dims, ierror)

  call h5aclose_f(attr_id, ierror)

  call h5fclose_f(file_id, ierror)

  call h5close_f(ierror)


! NetCDF4
  call nc_check( nf90_open(trim(filename), nf90_nowrite, ncid) )

  call nc_check( nf90_get_att(ncid, nf90_global, "MeanAnomaly", header%MeanAnomaly) )
  call nc_check( nf90_get_att(ncid, nf90_global, "MeanMotion", header%MeanMotion) )
  call nc_check( nf90_get_att(ncid, nf90_global, "OrbitalInclination", header%OrbitalInclination) )
  call nc_check( nf90_get_att(ncid, nf90_global, "EpochTime", header%EpochTime) )

  call nc_check( nf90_get_att(ncid, nf90_global, "Orbit Number", header%OrbitNumber) )

  !call nc_check( nf90_get_att(ncid, nf90_global, "Orbit Direction", header%OrbitDirection) )

  call nc_check( nf90_get_att(ncid, nf90_global, "Observing Beginning Date", header%startDate) )
  call nc_check( nf90_get_att(ncid, nf90_global, "Observing Beginning Time", header%startTime) )
  call nc_check( nf90_get_att(ncid, nf90_global, "Observing Ending Date"   , header%endDate) )
  call nc_check( nf90_get_att(ncid, nf90_global, "Observing Ending Time"   , header%endTime) )

  call nc_check( nf90_get_att(ncid, nf90_global, "Orbit Point Latitude", header%scLatPoint) )
  call nc_check( nf90_get_att(ncid, nf90_global, "Orbit Point Longitude", header%scLonPoint) )

  header%startDateTime = header%startDate//"T"//header%startTime//"Z"
  header%endDateTime = header%endDate//"T"//header%endTime//"Z"

  grpname = "Geolocation"

  call nc_check( nf90_inq_grp_ncid(ncid, trim(grpname), grpid) )

  subgrpname = trim(freqname)

  call nc_check( nf90_inq_grp_ncid(grpid, trim(subgrpname), subgrpid) )

  call nc_check( nf90_inq_dimids(subgrpid, dimnum, dimids, iparent) )

  do ijk = 1, dimnum
    call nc_check( nf90_inquire_dimension(subgrpid, dimids(ijk), dimname, dimsize) )
    if(ijk == 1) nscan = dimsize
    if(ijk == 2) nray  = dimsize
    if(ijk == 4) nlev  = dimsize
  enddo

  allocate( LatitudePMR( 2, nray, nscan ), source = MISSING_FLOAT )
  allocate( LongitudePMR( 2, nray, nscan ), source = MISSING_FLOAT )

  allocate( ellipsoidBinOffset( nray, nscan ), source = MISSING_FLOAT )

  allocate( localZenithAngle( nray, nscan ), source = MISSING_FLOAT )

  allocate( elevation( nray, nscan ), source = MISSING_FLOAT )

  allocate( height( nlev, nray, nscan ), source = MISSING_FLOAT )

  allocate( LandSurfaceTypePMR( nray, nscan ), source = MISSING_SHORT )
  allocate( msCountPMR( nscan ), source = MISSING_INT )
  allocate( dayCountPMR( nscan ), source = MISSING_SHORT )

  allocate( LandSurfaceType( nray, nscan ), source = MISSING_INT )
  allocate( Latitude( nray, nscan ), source = MISSING_FLOAT )
  allocate( Longitude( nray, nscan ), source = MISSING_FLOAT )

  varname = "Latitude"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, LatitudePMR) )

  varname = "Longitude"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, LongitudePMR) )

  varname = "ellipsoidBinOffset"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, ellipsoidBinOffset) )

  varname = "localZenithAngle"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, localZenithAngle) )

  varname = "LandSurfaceType"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, LandSurfaceTypePMR) )

  varname = "elevation"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, elevation) )

  varname = "height"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, height) )

  varname = "dayCount"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, dayCountPMR) )

  varname = "msCount"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, msCountPMR) )

  grpname = "FLG"

  call nc_check( nf90_inq_grp_ncid(ncid, trim(grpname), grpid) )

  subgrpname = trim(freqname)

  call nc_check( nf90_inq_grp_ncid(grpid, trim(subgrpname), subgrpid) )

  allocate( qualityDataPMR( nray, nscan ), source = MISSING_SHORT )
  allocate( qualityData( nray, nscan ), source = MISSING_INT )

  varname = "qualityData"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, qualityDataPMR) )

  grpname = "PRE"

  call nc_check( nf90_inq_grp_ncid(ncid, trim(grpname), grpid) )

  subgrpname = trim(freqname)

  call nc_check( nf90_inq_grp_ncid(grpid, trim(subgrpname), subgrpid) )

  allocate( zFactorMeasured( nlev, nray, nscan ), source = MISSING_FLOAT )
  allocate( sigmaZeroMeasured( nray, nscan ), source = MISSING_FLOAT )
  allocate( snRatioAtRealSurface( nray, nscan ), source = MISSING_FLOAT )
  allocate( binRealSurface( nray, nscan ), source = MISSING_SHORT )
  allocate( binFirstLatlon( nray, nscan ), source = MISSING_SHORT )
  allocate( binClutterFreeBottom( nray, nscan ), source = MISSING_SHORT )
  allocate( flagPrecipPMR( nray, nscan ), source = MISSING_BYTE )
  allocate( binStormTop( nray, nscan ), source = MISSING_SHORT )
  allocate( heightStormTop( nray, nscan ), source = MISSING_FLOAT )
  allocate( snowIceCover( nray, nscan ), source = MISSING_BYTE )

  allocate( flagPrecip( nray, nscan ), source = MISSING_INT )

  varname = "binRealSurface"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, binRealSurface) )

  varname = "zFactorMeasured"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, zFactorMeasured) )

  varname = "sigmaZeroMeasured"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, sigmaZeroMeasured) )

  varname = "snRatioAtRealSurface"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, snRatioAtRealSurface) )

  varname = "binRealSurface"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, binRealSurface) )

  varname = "binFirstLatlon"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, binFirstLatlon) )

  varname = "binClutterFreeBottom"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, binClutterFreeBottom) )

  varname = "flagPrecip"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, flagPrecipPMR) )

  varname = "binStormTop"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, binStormTop) )

  varname = "heightStormTop"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, heightStormTop) )

  varname = "snowIceCover"
  call nc_check( nf90_inq_varid(subgrpid, trim(varname), varid) )
  call nc_check( nf90_get_var(subgrpid, varid, snowIceCover) )

  call nc_check( nf90_close(ncid) )

  ! adjust time

  allocate( Year( nscan ), source = MISSING_SHORT )
  allocate( DayOfYear( nscan ), source = MISSING_SHORT )

  allocate( Month( nscan ), source = MISSING_BYTE )
  allocate( DayOfMonth( nscan ), source = MISSING_BYTE )

  allocate( Hour( nscan ), source = MISSING_BYTE )
  allocate( Minute( nscan ), source = MISSING_BYTE )
  allocate( Second( nscan ), source = MISSING_BYTE )

  allocate( SecondOfDay( nscan ), source = MISSING_DOUBLE )
  allocate( MilliSecond( nscan ), source = MISSING_SHORT )

  allocate( FractionalGranuleNumber( nscan ), source = MISSING_DOUBLE )

block
  integer :: ia
  ia = 0
  if(header%OrbitDirection == "A") ia = 1
  if(header%OrbitDirection == "D") ia = 2

  FractionalGranuleNumber(:) = header%OrbitNumber + 0.5*(ia-1)
end block

  call setmktime(20000101, 12*3600, values(:)) ! 2000/01/01 12:00:00 UTC is base datetime ???

  imktime = mkutctime_c(values(:))

  do iscan = 1, nscan
    MilliSecond(iscan) = mod(int(msCountPMR(iscan)*1e-1), 1000)

    call sub_gmtime_c( int(imktime + dayCountPMR(iscan)*86400 + int(msCountPMR(iscan)*1e-4), int64), values(:))

    !call timeformat(values(:), ctimefmt) ; write(ctimefmt(21:23), "(i3.3)") MilliSecond(iscan)

    call mktime2datesod(values, idate, isod)

    Year(iscan) = idate/10000
    Month(iscan) = mod(idate/100,100)
    DayOfMonth(iscan) = mod(idate,100)

    Hour(iscan) = values(3)
    Minute(iscan) = values(2)
    Second(iscan) = values(1)

    SecondOfDay(iscan) = isod + MilliSecond(iscan)*1d-3

    !print *, int(iscan,int16), MilliSecond(iscan), " "//ctimefmt
  enddo

  Latitude(:,:) = LatitudePMR(2,:,:)
  Longitude(:,:) = LongitudePMR(2,:,:)

  landSurfaceType(:,:) = int( LandSurfaceTypePMR(:,:), int32 )
  flagPrecip(:,:) = int(flagPrecipPMR(:,:), int32 )

  qualityData(:,:) = int( qualityData(:,:), int32 )

  deallocate( LatitudePMR, LongitudePMR )
  deallocate( LandSurfaceTypePMR )
  deallocate( flagPrecipPMR )
  deallocate( qualityDataPMR )

  deallocate( dayCountPMR, msCountPMR )

end subroutine read_PMR_L1

subroutine calc_PMR_L1
  implicit none

  real :: dBzmSTH, dBzmCFB

  do iscan = 1, nscan
    call timeformat(values(:), ctimefmt) ; write(ctimefmt(21:23), "(i3.3)") MilliSecond(iscan)
    do iray = 1, nray

      if(binStormTop(iray,iscan) > 0) then
        dBzmSTH = zFactorMeasured(binStormTop(iray,iscan),iray,iscan)
      else
        dBzmSTH = -9999.9e0
      endif

      if(binClutterFreeBottom(iray,iscan) > 0) then
        dBzmCFB = zFactorMeasured(binClutterFreeBottom(iray,iscan),iray,iscan)
      else
        dBzmCFB = -9999.9e0
      endif

      print *, int(header%OrbitNumber,int16), int(iscan,int16), int(iray,int8), ctimefmt, Longitude(iray,iscan), Latitude(iray,iscan), sigmaZeroMeasured(iray,iscan), &
      flagPrecip(iray,iscan), heightStormTop(iray,iscan), dBzmSTH, dBzmCFB        
    enddo
  enddo

end subroutine calc_PMR_L1

subroutine free_PMR_L1
  implicit none

! PRE
  if( allocated( zFactorMeasured ) )  deallocate( zFactorMeasured )
  deallocate( sigmaZeroMeasured )
  deallocate( snRatioAtRealSurface )
  deallocate( binRealSurface )
  deallocate( binFirstLatlon )
  deallocate( binClutterFreeBottom )
  deallocate( flagPrecip )
  deallocate( binStormTop )
  deallocate( heightStormTop )
  deallocate( snowIceCover )

! FLG
  deallocate( qualityData )

! Geolocation
  if( allocated( height ) ) deallocate( height )
  deallocate( ellipsoidBinOffset, localZenithAngle, elevation, landSurfaceType)

  deallocate( Year, Month, DayOfMonth, DayOfYear, SecondOfDay, MilliSecond )
  deallocate( Hour, Minute, Second )

  deallocate( Longitude, Latitude )
  deallocate( FractionalGranuleNumber )

end subroutine free_PMR_L1

subroutine nc_check( istat )
  implicit none

  integer, intent(in) :: istat

  if(istat /= nf90_noerr) then
    print *, trim(nf90_strerror(istat))
    stop !"Stopped"
  endif

end subroutine nc_check

subroutine setmktime(idate, isod, values)
  implicit none
  integer, intent(in) :: idate, isod
  integer, intent(out) :: values(9)

  values(1) = mod(isod,60)
  values(2) = mod(int(isod/60),60)
  values(3) = mod(int(isod/3600),60)
  values(4) = mod(idate,100)
  values(5) = mod(int(idate/100),100) -1
  values(6) = int(idate/10000) - 1900
  values(7:9) = 0

end subroutine setmktime

subroutine mktime2datesod(values, idate, isod)
  implicit none
  integer, intent(in) :: values(9)
  integer, intent(out) :: idate, isod

  idate = ( values(6) + 1900 )*10000 + ( values(5) + 1)*100 + ( values(4) )
  isod = (values(3)*3600 + values(2)*60 + values(1) )

end subroutine mktime2datesod

subroutine timeformat(values, ctimefmt)
  implicit none
  integer, intent(in) :: values(9)
  character(len=*), intent(out) :: ctimefmt

  ctimefmt = ""

  write( ctimefmt( 1: 4), "(i4.4)") values(6) + 1900
  ctimefmt( 5: 5) = "-"
  write( ctimefmt( 6: 7), "(i2.2)") values(5) + 1
  ctimefmt( 8: 8) = "-"
  write( ctimefmt( 9:10), "(i2.2)") values(4)
  ctimefmt(11:11) = "T"
  write( ctimefmt(12:13), "(i2.2)") values(3)
  ctimefmt(14:14) = ":"
  write( ctimefmt(15:16), "(i2.2)") values(2)
  ctimefmt(17:17) = ":"
  write( ctimefmt(18:19), "(i2.2)") values(1)

  select case(len(ctimefmt))
    case(20)
      ctimefmt(20:20) = "Z"
    case(23)
      ctimefmt(11:11) = " "
      ctimefmt(20:23) = " UTC"
    case(24)
      ctimefmt(20:23) = ".000"
      ctimefmt(24:24) = "Z"
  end select

end subroutine timeformat

end program read_PMR