Recipe: Fortran program to read data from downloaded file
Back to G5NR Data Access Guide.
Problem
We want to read a downloaded data file using Fortran.
Solution
For the purpose of this example, we assume that we have already downloaded the file c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4
from the ftp server. For more information about file naming conventions, and how to download a file from the ftp server, please follow the links in the #See Also section.
Here we use the NetCDF-4 library to read the downloaded file (since NetCDF-4 uses HDF-5 as the underlying data format, HDF-5 library can also be used directly). First, we ensure that out NetCDF-4 library has been built with Fortran and HDF5 support. If both the queries
> nc-config --has-f90 > nc-config --has-hdf5
return "yes", we have a compatible NetCDF-4 library. Here nc-config
is a utility bundled with NetCDF-4 package.
Code
This code reads the global temperature data from file and computes the maximum and minimum temperatures.
program g5nr_reader_global
use netcdf ! for reading the NR files
implicit none
! File name
! ---------
character(len=256) :: T_file
! 4D array: (lon,lat,lev,time)
! ----------------------------
real, allocatable :: T(:,:,:,:) ! global data
! Miscellaneous
! -------------
integer :: ierr
integer :: im, jm, lm
integer :: ncid, varid
integer :: start(4), count(4)
! file name and dimensions
! ------------------------
T_file = "c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4"
im = 720
jm = 361
lm = 72
! Open file and get var id
! ------------------------
call check(nf90_open(T_file,NF90_NOWRITE,ncid), "opening T file")
call check(nf90_inq_varid(ncid,"T",varid), "getting T varid")
! Read global temperature data
! ----------------------------
allocate(T(im,jm,lm,1)) ! global 4D array with 1 time level
start = [1, 1, 1, 1]
count = [im, jm, lm, 1] ! 1 time level, 3D (lon,lat,lev) array
write(*,*) 'Reading global T...'
call check(nf90_get_var(ncid,varid,T,start=start,count=count), "reading T")
write(*,*)'T: ', maxval(T),minval(T)
! close file, release memory
call check(nf90_close(ncid), "closing T file")
deallocate(T)
! All done
! --------
contains
subroutine check(status, loc)
integer, intent(in) :: status
character(len=*), intent(in) :: loc
if(status /= NF90_NOERR) then
write (*,*) "Error at ", loc
write (*,*) NF90_STRERROR(status)
end if
end subroutine check
end program g5nr_reader_global
Compile and link
We use the utility nf-config
(included in the NetCDF-4 installation) to identify the linking rules.
For a typical NetCDF-4 installation, the above code, g5nr_reader_global.f90
can be compiled and linked to the NetCDF-4 library via
gfortran -o g5nr_reader_global.x `nf-config --fflags` g5nr_reader_global.f90 `nf-config --flibs`
creating the executable g5nr_reader_global.x
.
NOTE:
- You can use your favorite Fortran compiler instead of gfortran.
- If the NetCDF-4 library was built with parallel I/O support, you will need to use
mpif90
to link, even if your code does not use the MPI library.
Run
Running the executable
> ./g5nr_reader_global.x
produces the output
Reading global T... T: 315.6512 180.3667
Discussions
Modifications to read a subset of the data
The above Fortran code can be easily modified to read a subset of the temperature data instead. If we want to read temperature inside the box bounded by latitudes 25oN, 50oN and longitudes -130oN, -65oN, we need to modify the values of the dimensions im
, jm
and the argument start
of the call to nf90_get_var
.
Our new dimensions are
! array indices from boundary lat/lon
imin = nint((-130. + 180.)/0.5)
imax = nint((- 65. + 180.)/0.5)
jmin = nint(( 25. + 90.)/0.5)
jmax = nint(( 50. + 90.)/0.5)
! array sizes
im = imax-imin+1
jm = jmax-jmin+1
The new start
argument to nf90_get_var
defining the first indices of the slice of array to be read, is
start = [imin, jmin, 1, 1]
The modified version of the code is g5nr_reader_subset.txt.
Running this modified version
./g5nr_reader_subset.x
produces the output
Reading Tsub Tsub: 305.6048 191.6956
See Also
- File Spec: File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf
- Recipe: Retrieve (global) data from FTP server
- Recipe: Fortran program as OPeNDAP client