Recipe: Fortran program as OPeNDAP client: Difference between revisions

From GEOS-5
Jump to navigation Jump to search
← Older editNewer edit →
Pchakrab (talk | contribs)
1,048 edits
Pchakrab (talk | contribs)
1,048 edits
Line 7: Line 7:
== Solution ==
== Solution ==


First, we ensure that out NetCDF-4 library has been built with Fortran, HDF5 and OPeNDAP support. All three queries
The OPeNDAP server for G5NR data is located at
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg


The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at
http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info
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-f90
> nc-config --has-hdf5
  > nc-config --has-dap
  > nc-config --has-dap


should return yes.
return "yes", we have a compatible NetCDF-4 library. Here <code>nc-config</code> is a utility bundled with NetCDF-4 package.


== Read downloaded nc4 file ==
== Read downloaded nc4 file ==

Revision as of 13:46, 29 October 2014

Back to G5NR Data Access Guide.

Problem

By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server, 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

The OPeNDAP server for G5NR data is located at 

http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg

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

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.

Read downloaded nc4 file

We first download the collection inst01hr_3d_T_Cv for 2006-Sep-18/0900z as described here. The downloaded file is c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4. The code below does the following

  • reads global temperature data and computes its max/min
  • reads a subset of the data and computes the subset's max/min

Code

program g5nr_reader

  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 = "c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4"

  ! 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)

  ! 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, 1]
  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

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:

  1. You can use your favorite Fortran compiler instead of gfortran.
  2. 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.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

  1. Line 31: Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
  2. 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 to nf90_inq_varid.
  3. Line 41: Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly inst files are available starting at 2005/15/15 2200z. Our desired time, 2006/09/18 0900z is then the 11772th file.
  4. 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:

  1. You can use your favorite Fortran compiler instead of gfortran.
  2. 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
Retrieved from "https://geos5.org/wiki/index.php?title=Recipe:_Fortran_program_as_OPeNDAP_client&oldid=3013"