program main
  use hdf5
  implicit none

  integer, parameter :: nMaxScan = 9000

  integer :: iscan, iray, iz, l

  integer :: ierror, file_id, group_id, dset_id

  character*200 :: filename, groupname, dsetname, productname

  !real, dimension (49, 8000) :: lat = -99999, lon = -99999
  !real, dimension (176, 49, 8000) :: precipRate = -99999

  real, allocatable, dimension (:,:) :: lon, lat, precipRateES
  real, allocatable, dimension (:,:,:) :: precipRate, zFactorCorrected
  real, allocatable, dimension (:,:,:,:) :: paramDSD

  integer, parameter :: nMaxRay = 49

  integer*8, dimension(2) :: dim2d = (/nMaxRay, nMaxScan/)
  integer*8, dimension(3) :: dim3d = (/176, nMaxRay, nMaxScan/)
  integer*8, dimension(4) :: dim4d = (/2, 176, nMaxRay, nMaxScan/)

  real :: Nw, Dm, D, ND, ND2, mu

  allocate( lat( nMaxRay, nMaxScan) ) ; lat = -99999
  allocate( lon( nMaxRay, nMaxScan) ) ; lon = -99999
  allocate( precipRateES( nMaxRay, nMaxScan) ) ; precipRateES = -99999
  allocate( precipRate( 176, nMaxRay, nMaxScan) ) ; precipRate = -99999
  allocate( zFactorCorrected( 176, nMaxRay, nMaxScan) ) ; zFactorCorrected = -99999
  allocate( paramDSD( 2, 176, nMaxRay, nMaxScan) ) ; paramDSD = -99999

  productname = "NS"

  call getarg(1, filename)
  if(iargc() /= 1) call exit(1)

  call h5open_f(ierror)
  !print *, ierror

  call h5fopen_f( filename, H5F_ACC_RDONLY_F, file_id, ierror)
  !print *, file_id, ierror

!-----------------------------------------------------------------------
  groupname = "/"//trim(productname)
  call h5gopen_f(file_id, groupname, group_id, ierror)
  !print *, group_id, ierror

  dsetname = "Latitude"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, lat,  dim2d, ierror)
  !print *, ierror

  call h5dclose_f(dset_id, ierror)
  !print *, ierror

  dsetname = "Longitude"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, lon, dim2d, ierror)
  !print *, ierror

  call h5dclose_f(dset_id, ierror)
  !print *, ierror

  call h5gclose_f(group_id, ierror)
  !print *, ierror
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
  groupname = "/"//trim(productname)//"/SLV"
  call h5gopen_f(file_id, groupname, group_id, ierror)
  !print *, group_id, ierror

  dsetname = "precipRate"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, precipRate,  dim3d, ierror)
  !print *, ierror

  dsetname = "zFactorCorrected"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, zFactorCorrected,  dim3d, ierror)

  call h5dclose_f(dset_id, ierror)
  !print *, ierror

  dsetname = "precipRateESurface"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, precipRateES,  dim2d, ierror)
  !print *, ierror

  call h5dclose_f(dset_id, ierror)
  !print *, ierror

  dsetname = "paramDSD"
  call h5dopen_f(group_id, dsetname, dset_id, ierror)
  !print *, dset_id, ierror

  call h5dread_f(dset_id, H5T_NATIVE_REAL, paramDSD,  dim4d, ierror)
  !print *, ierror

  call h5dclose_f(dset_id, ierror)
  !print *, ierror

  call h5gclose_f(group_id, ierror)
!-----------------------------------------------------------------------

  call h5fclose_f(file_id, ierror)
  !print *, ierror

  call h5close_f(ierror)
  !print *, ierror

  open(10,file="out_2d.txt", status="replace")

  open(20,file="out_3d.txt", status="replace")

  do iscan = 1, nMaxScan
    do iray = 1, 49
      if(lon(iray,iscan) == -99999) cycle

        write(10, "(i6,i4, 3f12.3)" ) &
        iscan, iray, lon(iray,iscan), lat(iray, iscan), &
        precipRateES(iray,iscan)

      !if(precipRateES(iray,iscan) > 0) then
      !  print  "(i6,i4, 3f12.3)", &
      !  iscan, iray, lon(iray,iscan), lat(iray, iscan), &
      !  precipRateES(iray,iscan)
      !endif

      do iz = 1, 176
        if(precipRate(iz,iray,iscan) <= 0.0) cycle

        Nw = 10**(paramDSD(1, iz, iray, iscan)/10.)
        Nw = paramDSD(1, iz, iray, iscan)
        Dm = paramDSD(2, iz, iray, iscan)
        mu = 3.0


!        if(precipRate(iz,iray,iscan) <= 50.0) cycle

        write(20, "(i6,i4,i4, 6f12.3)") &
        iscan, iray, iz, lon(iray,iscan), lat(iray, iscan), &
        precipRate(iz, iray, iscan), paramDSD(1:2, iz, iray, iscan)
cycle
        do l = 1, 1000
           D = l/100.

           ND = Nw*D**(mu)*exp(- (4 + mu)*D/Dm)

           ND2 = Nw*6*((mu + 4)**(mu + 4))/((4**4)*gamma(mu + 4))*((D/Dm)**mu)*exp(-(4 + mu)*D/Dm)
!           print *, real(D), ND, ND2, Nw, Dm !, precipRate(iz,iray,iscan)
!           print *, real(D), ND, ND2, precipRate(iz,iray,iscan), zFactorCorrected(iz,iray,iscan)
stop
        enddo
stop

      enddo
    enddo
  enddo

  close(10)
  close(20)


  deallocate(lon, lat, precipRateES, precipRate, paramDSD)

end program main