GEOS-5

Recipe: Fortran program as OPeNDAP client: Difference between revisions

← Older edit
Pchakrab (talk | contribs)
1,048 edits
Move from http to https
 
(19 intermediate revisions by one other user not shown)
Line 4: Line 4:


By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server
By accessing the collection <code>inst01hr_3d_T_Cv</code> via the OPeNDAP server
  http://opendap.nccs.nasa.gov/dods/OSSE/G5NR/Ganymed/7km
  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.
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.
Line 16: Line 16:


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


==== Code ====
==== Code ====
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.
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_dap
program g5nr_reader_dap
Line 32: Line 39:
   ! 4D array: (lon,lat,lev,time)
   ! 4D array: (lon,lat,lev,time)
   ! ------------------------------------
   ! ------------------------------------
   real, allocatable :: T(:,:,:,:) ! subset of global data
   real, allocatable :: T(:,:,:,:)


   ! Miscellaneous
   ! Miscellaneous
Line 42: Line 49:
   integer :: imin, imax, jmin, jmax
   integer :: imin, imax, jmin, jmax
   real :: minlat, minlon, maxlat, maxlon
   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
   ! file name and (subset) dimensions
   ! bounding box:
   ! bounding box:
   !  lons = -180:0.5:179.5
   !  lons = -130:0.5:-65
   !  lats = -90:0.5:90
   !  lats =   25:0.5:50
   ! ---------------------------------
   ! ---------------------------------
  T_file = "http://opendap.nccs.nasa.gov:9090/dods/OSSE/G5NR/Ganymed/7km/0.5000_deg/inst/inst01hr_3d_T_Cv"
   ! array indices from the
  minlon = -130.
  maxlon = -65.
  minlat = 25.
  maxlat = 50.
   ! compute array indices from the
   ! bounding lat/lon values
   ! bounding lat/lon values
   imin = nint((minlon + 180.)/0.5)
   imin = nint((-130. + 180.)/0.5)
   imax = nint((maxlon + 180.)/0.5)
   imax = nint(( -65. + 180.)/0.5)
   jmin = nint((minlat +  90.)/0.5)
   jmin = nint(( 25 +  90.)/0.5)
   jmax = nint((maxlat +  90.)/0.5)
   jmax = nint(( 50 +  90.)/0.5)
   ! compute array sizes
   ! array sizes
   im = imax-imin+1
   im = imax-imin+1
   jm = jmax-jmin+1
   jm = jmax-jmin+1
   lm = 72
   lm = 72


   ! 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 = [imin, jmin, 1, 11772]           ! time level 11772 corresponds to 2006/09/18, 9z
   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 T...'
   write(*,*) 'Reading T...'
Line 106: 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 118: 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 ==
== Discussion ==


==== Code ====
== See Also ==


The code to read data from the OPeNDAP server is a minor modification (4 lines) of the above program and is reproduced below
# File Spec: [[File:G5NR-Ganymed-7km FileSpec-ON6-V1.0.pdf]]
# [[Recipe: Fortran program to read data from downloaded file]]


<syntaxhighlight lang="fortran" line>
== No Warranty ==
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
</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
Retrieved from "https://geos5.org/wiki/index.php?title=Recipe:_Fortran_program_as_OPeNDAP_client"