Recipe: Fortran program as OPeNDAP client
Back to G5NR Data Access Guide.
Problem
By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km
we want to read the surface temperature data inside the box bound by latitudes 25oN, 50oN and longitudes -130oW, -65oW for 2006/Sep/18, 9z and compute its min/max.
Solution
First, we ensure that our NetCDF-4 library has been built with Fortran and OPeNDAP support. If both the queries
> nc-config --has-f90 > nc-config --has-dap
return "yes", we have a compatible NetCDF-4 library. Here nc-config is a utility bundled with NetCDF-4 package.
The metadata for the collection inst01hr_3d_T_Cv is available at
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info
Code
This code accesses the collection inst01hr_3d_T_Cv from the OPeNDAP server and reads a subset of the temperature data.
NOTE: For hourly data, starting at 2005/May/15, 22z, the 11772th timestep corresponds to 2006/Sep/18, 9z.
program g5nr_reader_dap
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(:,:,:,:)
! Miscellaneous
! -------------
integer :: ierr
integer :: im, jm, lm
integer :: ncid, varid
integer :: start(4), count(4)
integer :: imin, imax, jmin, jmax
real :: minlat, minlon, maxlat, maxlon
! We replace filename by opendal url
! ----------------------------------
T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"
! file name and (subset) dimensions
! bounding box:
! lons = -130:0.5:-65
! lats = 25:0.5:50
! ---------------------------------
! array indices from the
! bounding lat/lon values
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
lm = 72
! Open url 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 temperature data
! ---------------------
allocate(T(im,jm,lm,1)) ! global 4D array with 1 time level
start = [imin, jmin, 1, 11772] ! time level 11772 corresponds to 2006/09/18, 9z
count = [im, jm, lm, 1] ! 1 time level, 3D (lon,lat,lev) array
write(*,*) 'Reading 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_dap
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.f90 can be compiled and linked to the NetCDF-4 library via
gfortran -o g5nr_reader.x `nf-config --fflags` g5nr_reader.f90 `nf-config --flibs`
creating the executable g5nr_reader.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
mpif90to link, even if your code does not use the MPI library.
Run
Running the executable,
./g5nr_reader.x
produces the output
Reading global T... T: 315.6512 180.3667 Reading T over US Tsub: 302.8958 183.1974
Access data via OPeNDAP server
Code
The code to read data from the OPeNDAP server is a minor modification (4 lines) of the above program and is reproduced below
program g5nr_reader_dap
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
real, allocatable :: Tsub(:,:,:,:) ! subset of global data
! Miscellaneous
! -------------
integer :: ierr
integer :: im, jm, lm
real :: minlat, minlon, maxlat, maxlon
integer :: im_sub, jm_sub
integer :: imin, imax, jmin, jmax
integer :: ncid, varid
integer :: start(4), count(4)
! file name and dimensions
! ------------------------
im = 720
jm = 361
lm = 72
T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv"
! 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, 11772] ! time level 11772 corresponds to 2006/09/18, 9z
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)
! Read temperature over US
! lons = -180:0.5:179.5
! lats = -90:0.5:90
! ------------------------
write(*,*) 'Reading T over US'
minlon = -130.
maxlon = -65.
minlat = 25.
maxlat = 50.
imin = nint((minlon + 180.)/0.5)
imax = nint((maxlon + 180.)/0.5)
jmin = nint((minlat + 90.)/0.5)
jmax = nint((maxlat + 90.)/0.5)
im_sub = imax-imin+1
jm_sub = jmax-jmin+1
allocate(Tsub(im_sub,jm_sub,lm,1)) ! global 4D array with 1 time level
start = [1, 1, 1, 11772] ! time level 11772 corresponds to 2006/09/18, 9z
count = [im_sub, jm_sub, lm, 1] ! 1 time level, 3D (lon,lat,lev) array
call check(nf90_get_var(ncid,varid,Tsub,start=start,count=count), "reading Tsub")
write(*,*)'Tsub: ', maxval(Tsub),minval(Tsub)
! close file, release memory
call check(nf90_close(ncid), "closing T file")
deallocate(T)
deallocate(Tsub)
! 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_dap
Discussion
The four (4) differences from the case where an individual downloaded file was being read are
Line 31: Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.Line 36: While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase (see http://opendap.nccs.nasa.gov:9090/dods/OSSE/GEOS-5.12/BETA9/0.5000_deg/inst/inst01hr_3d_T_Cv). This is reflected in the call tonf90_inq_varid.Line 41: Via the OPeNDAP URL, we now have access to all times for which data exists. The hourlyinstfiles are available starting at 2005/15/15 2200z. Our desired time, 2006/09/18 0900z is then the 11772th file.Line 63: Same as the difference above, but for retrieving a subset of data.
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_dap.f90 can be compiled and linked to the NetCDF-4 library via
gfortran -o g5nr_reader_dap.x `nf-config --fflags` g5nr_reader_dap.f90 `nf-config --flibs`
creating the executable g5nr_reader_dap.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_dap.x
produces the output
Reading global T... T: 315.6512 180.3667 Reading T over US Tsub: 302.8958 183.1974