Recipe: Fortran program as OPeNDAP client: Difference between revisions

From GEOS-5
Jump to navigation Jump to search
Pchakrab (talk | contribs)
Move from http to https
 
(26 intermediate revisions by one other user not shown)
Line 3: Line 3:
== Problem ==
== Problem ==


By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server, we want to read the surface temperature data inside the box bound by latitudes 25<sup>o</sup>N, 50<sup>o</sup>N and longitudes -130<sup>o</sup>W, -65<sup>o</sup>W for 2006/Sep/18, 9z and compute its min/max.
By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server
https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km
 
we want to read the surface temperature data inside the box bound by latitudes 25<sup>o</sup>N, 50<sup>o</sup>N and longitudes -130<sup>o</sup>W, -65<sup>o</sup>W for 2006/Sep/18, 9z and compute its min/max.


== Solution ==
== 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 <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
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
Line 19: Line 15:
return "yes", we have a compatible NetCDF-4 library. Here <code>nc-config</code> is a utility bundled with NetCDF-4 package.
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 ==
The metadata for the collection <code>inst01hr_3d_T_Cv</code> is available at
https://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv.info


We first download the collection inst01hr_3d_T_Cv for 2006-Sep-18/0900z as described [[G5NR Data Access Guide#ftp (global data)|here]]. The downloaded file is <code>c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4</code>. The code below does the following
==== Code ====
* reads global temperature data and computes its max/min
This code accesses the collection <code>inst01hr_3d_T_Cv</code> from the OPeNDAP server, reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min.
* reads a subset of the data and computes the subset's max/min


==== Code ====
NOTE:
# Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
# While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase. This is reflected in the call to <code>nf90_inq_varid</code>.
# Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly <code>inst</code> files are available starting at 2005/15/15, 2200z. Our desired time, 2006/09/18, 0900z is then the 11772<sup>th</sup> step.


<syntaxhighlight lang="fortran" line>
<syntaxhighlight lang="fortran" line>
program g5nr_reader
program g5nr_reader_dap


   use netcdf          ! for reading the NR files
   use netcdf          ! for reading the NR files
Line 40: Line 39:
   ! 4D array: (lon,lat,lev,time)
   ! 4D array: (lon,lat,lev,time)
   ! ------------------------------------
   ! ------------------------------------
   real, allocatable :: T(:,:,:,:)   ! global data
   real, allocatable :: T(:,:,:,:)
  real, allocatable :: Tsub(:,:,:,:) ! subset of global data


   ! Miscellaneous
   ! Miscellaneous
Line 47: Line 45:
   integer :: ierr
   integer :: ierr
   integer :: im, jm, lm
   integer :: im, jm, lm
  real :: minlat, minlon, maxlat, maxlon
  integer :: im_sub, jm_sub
  integer :: imin, imax, jmin, jmax
   integer :: ncid, varid
   integer :: ncid, varid
   integer :: start(4), count(4)
   integer :: start(4), count(4)
  integer :: imin, imax, jmin, jmax
  real :: minlat, minlon, maxlat, maxlon
  ! We replace filename by opendal url
  ! ----------------------------------
  T_file = "https://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"


   ! file name and dimensions
   ! file name and (subset) dimensions
   ! ------------------------
   ! bounding box:
   im = 720
  !  lons = -130:0.5:-65
   jm = 361
  !  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
   lm = 72
  T_file = "c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4"


   ! Open file and get var id
   ! Open url and get var id
   ! ------------------------
   ! -----------------------
   call check(nf90_open(T_file,NF90_NOWRITE,ncid), "opening T file")
   call check(nf90_open(T_file,NF90_NOWRITE,ncid), "opening T file")
   call check(nf90_inq_varid(ncid,"T",varid), "getting T varid")
   call check(nf90_inq_varid(ncid,"t",varid), "getting T varid")


   ! Read global temperature data
   ! Read temperature data
   ! ----------------------------
   ! ---------------------
   allocate(T(im,jm,lm,1))            ! global 4D array with 1 time level
   allocate(T(im,jm,lm,1))            ! global 4D array with 1 time level
   start = [1, 1, 1, 1]
   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
   count = [im, jm, lm, 1]            ! 1 time level, 3D (lon,lat,lev) array
   write(*,*) 'Reading global T...'
   write(*,*) 'Reading T...'
   call check(nf90_get_var(ncid,varid,T,start=start,count=count), "reading T")
   call check(nf90_get_var(ncid,varid,T,start=start,count=count), "reading T")
   write(*,*)'T: ', maxval(T),minval(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
   ! close file, release memory
   call check(nf90_close(ncid), "closing T file")
   call check(nf90_close(ncid), "closing T file")
   deallocate(T)
   deallocate(T)
  deallocate(Tsub)


   ! All done
   ! All done
Line 117: Line 105:
   end subroutine check
   end subroutine check


end program g5nr_reader
end program g5nr_reader_dap
</syntaxhighlight>
</syntaxhighlight>


Line 124: Line 112:
We use the utility <code>nf-config</code> (included in the NetCDF-4 installation) to identify the linking rules
We use the utility <code>nf-config</code> (included in the NetCDF-4 installation) to identify the linking rules


For a typical NetCDF-4 installation, the above code, <code>g5nr_reader.f90</code> can be compiled and linked to the NetCDF-4 library via
For a typical NetCDF-4 installation, the above code, <code>g5nr_reader_dap.f90</code> 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`
  gfortran -o g5nr_reader_dap.x `nf-config --fflags` g5nr_reader_dap.f90 `nf-config --flibs`


creating the executable <code>g5nr_reader.x</code>.
creating the executable <code>g5nr_reader_dap.x</code>.


NOTE:
NOTE:
Line 136: Line 124:
Running the executable,
Running the executable,


  ./g5nr_reader.x
  ./g5nr_reader_dap.x


produces the output
produces the output
  Reading global T...
  Reading T (subset)...
  T:    315.6512       180.3667   
  T (subset):    305.6048       191.6956
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
 
<syntaxhighlight lang="fortran" line>
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)
== Discussion ==
  ! ------------------------------------
  real, allocatable :: T(:,:,:,:)    ! global data
  real, allocatable :: Tsub(:,:,:,:) ! subset of global data


  ! Miscellaneous
== See Also ==
  ! -------------
  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
# File Spec: [[File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf]]
  ! ------------------------
# [[Recipe: Fortran program to read data from downloaded file]]
  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
== No Warranty ==
  ! ------------------------
  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
</syntaxhighlight>
 
===== Discussion =====
 
The four (4) differences from the case where an individual downloaded file was being read are
 
# <code>Line 31</code>: Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
# <code>Line 36</code>: 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 <code>nf90_inq_varid</code>.
# <code>Line 41</code>: Via the OPeNDAP URL, we now have access to all times for which data exists. The hourly <code>inst</code> files are available starting at 2005/15/15 2200z. Our desired time, 2006/09/18 0900z is then the 11772<sup>th</sup> file.
# <code>Line 63</code>: 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 <code>g5nr_reader_dap.x</code>.
 
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...
== Copyright ==
T:    315.6512      180.3667   
Reading T over US
Tsub:    302.8958      183.1974

Latest revision as of 10:24, 9 April 2019

Back to G5NR Data Access Guide.

Problem

By accessing the collection inst01hr_3d_T_Cv via the OPeNDAP server

https://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

https://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, reads a subset of the temperature data (all levels inside the bounding box specified above) and computes its max/min.

NOTE:

  1. Instead of reading a downloaded NetCDF-4 file, we read an OPeNDAP URL.
  2. While in the downloaded file, the temperature variable appears in the uppercase (T), on the OPeNDAP server, this variable is in lowercase. This is reflected in the call to nf90_inq_varid.
  3. 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 step.
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 = "https://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_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 T (subset)...
T (subset):    305.6048       191.6956

Discussion

See Also

  1. File Spec: File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf
  2. Recipe: Fortran program to read data from downloaded file

No Warranty

Copyright