Use netcdf
Integer :: ncid, id_lat_dim, id_lon_dim, id_t_dim, varid, Nlat, Nlon, Ntime, ifile
Integer :: i,j, i0, i1, j0, j1, ic, Nx, Ny
Real, Allocatable :: Lat(:), Lon(:), Time(:), X(:,:,:)
Real, Allocatable :: Lat2(:), Lon2(:), Time2(:), Y(:,:,:)
Integer :: Nlat2, Nlon2, Ntime2
Character (Len=255) :: FileName, VarName, LatName, LonName, TName
Real :: Lat0, Lat1, Lon0, Lon1
Real, Allocatable :: XLat(:), XLon(:), XDat(:,:)
!------------------------------------------------------------------------------------

FileName = "/home/tanpv/SEACLID/Data/APH_025_1981_Prec_Mon.nc"
status = nf90_open(trim(FileName), nf90_NoWrite, ncid)
Call get_data_info("precip","latitude","longitude","time") !,X) !, Lat, Lon, Nlat, Nlon)
status =  nf90_close(ncid)
Where (X==-99.9) X=-1

print*, "Nlat, Nlon, Ntime=", Nlat, Nlon, Ntime
print*, "Lat="
print*, (Lat(j),j=1, Nlat)
print*, "Lon="
print*, (Lon(i),i=1, Nlon)
print*, "Time="
print*, (Time(i),i=1, Ntime)
!-------------------------------------------------------------------
Lon0=102; Lon1=110; Lat0=8; Lat1=23

Call get_domain_info (Lon0, Lon1, Lat0, Lat1, i0, i1, j0, j1, Nx, Ny)
print*, "i0, i1, j0, j1, Nx, Ny=",i0, i1, j0, j1, Nx, Ny

pause

Do itime=1,Ntime
!print*,X(:,:,itime)
Call get_data_domain(itime,i0, i1, j0, j1, Nx,Ny, XDat, XLon, XLat)

print*, "XDat, Nx, Ny, Mon=", Nx, Ny, itime+1
Do j=Ny, 1, -1 ! Ny
   print "(40F5.1)", (XDat(i,j),i=1, Nx)
Enddo
pause
Enddo

Call rd_ncdf_3d (trim(FileName),lon2,lat2,Time2, Nlon2, Nlat2, Ntime2,Y,"precip")
print*,">>>>>>>>>>>>>>>>"
print*, X-Y

!------------------------------------------------------------------------------------
CONTAINS

Subroutine get_domain_info (Lon0, Lon1, Lat0, Lat1, i0, i1, j0, j1, Nx, Ny)
Real :: Lon0, Lon1, Lat0, Lat1
Integer :: i,i0, i1, j0, j1, Nx, Ny
i=0
Do while (i<Nlon)
  i=i+1;  if (Lon(i)>=Lon0) then; i0=i; exit; endif
Enddo
i=Nlon+1
Do while (i>0)
  i=i-1;  if (Lon(i)<=Lon1) then; i1=i; exit; endif
Enddo
!---------
i=0
Do while (i<Nlat)
  i=i+1;  if (Lat(i)>=Lat0) then; j0=i; exit; endif
Enddo
i=Nlat+1
Do while (i>0)
  i=i-1;  if (Lat(i)<=Lat1) then; j1=i; exit; endif
Enddo
Nx = i1-i0+1; Ny = j1-j0+1
End Subroutine get_domain_info
!------------------------------------------------------------------------------------

Subroutine get_data_domain (itime,i0, i1, j0, j1, Nx,Ny, XDat, XLon, XLat) ! (Lat0, Lat1, Lon0, Lon1, XLat(:), XLon(:), XDat(:,:), Ny, Nx)
Integer :: i,j, i0, i1, j0, j1, itime, Nx, Ny
Real, Allocatable :: XDat(:,:),XLon(:), XLat(:)
if (Allocated(XDat)) Deallocate(XDat)
if (Allocated(XLon)) Deallocate(XLon)
if (Allocated(XLat)) Deallocate(XLat)
Allocate(XLat(Ny),XLon(Nx),XDat(Nx,Ny))
Do i=i0, i1
   XLon(i-i0+1) = Lon(i)
   Do j=j0, j1
      XLat(j-j0+1) = Lat(j)
      XDat(i-i0+1,j-j0+1) = X(i,j,itime)
   Enddo
Enddo
End Subroutine get_data_domain
!------------------------------------------------------------------------------------

Subroutine get_data_info (VarName, LatName, LonName, TName)!, X) !, Lat, Lon, Nlat, Nlon)
Character (Len=*) :: VarName, LatName, LonName, TName
  ! Get Lat Dim
  status = nf90_inq_dimid(ncid,LatName,id_lat_dim)
  status = nf90_inquire_dimension(ncid, id_lat_dim,len=Nlat)
  ! Get Lon Dim
  status = nf90_inq_dimid(ncid,LonName, id_lon_dim)
  status = nf90_inquire_dimension(ncid, id_lon_dim,len=Nlon)
  ! Get Time Dim
  status = nf90_inq_dimid(ncid,TName, id_t_dim)
  status = nf90_inquire_dimension(ncid, id_t_dim,len=Ntime)

  Allocate (X(Nlon,Nlat,Ntime),Lon(Nlon),Lat(Nlat), Time(Ntime))

  status = nf90_inq_varid(ncid,LatName, varid) ! Get Lat
  status = nf90_get_var(ncid, varid, Lat)
  status = nf90_inq_varid(ncid,LonName, varid) ! Get Lon
  status = nf90_get_var(ncid, varid, Lon)

  status = nf90_inq_varid(ncid,TName, varid)   ! Get Time
  status = nf90_get_var(ncid, varid, Time)
  status = nf90_inq_varid(ncid,VarName, varid) ! Get Var
  status = nf90_get_var(ncid, varid, X(:,:,:))
EndSubroutine get_data_info
!------------------------------------------------------------------------------------
Include "sub_rd_ncdf.F90"

END

!Include "sub_rd_ncdf.F90"