Recipe: Fortran program as OPeNDAP client: Difference between revisions

From GEOS-5
Jump to navigation Jump to search
Pchakrab (talk | contribs)
Pchakrab (talk | contribs)
Line 191: Line 191:
end program g5nr_reader_dap
end program g5nr_reader_dap
</syntaxhighlight>
</syntaxhighlight>
The three (3) differences are
# <code>T_file</code>: instead of reading a downloaded NetCDF-4 file, we read a URL.
# <code>start</code>:
# <code>nf90_inq_varid</code>


==== Compile and run ====
==== Compile and run ====

Revision as of 09:46, 27 October 2014

Back to G5NR Data Access Guide.

First, we ensure that out NetCDF-4 library has been built with Fortran, HDF5 and OPeNDAP support. All three queries

> nc-config --has-f90
> nc-config --has-hdf5
> nc-config --has-dap

should return yes.

Read downloaded nc4 file

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 program to read this file is

Code

program g5nr_reader

  use netcdf           ! for reading the NR files  

  implicit none

  ! File name                                                                                                                 
  ! ---------                                                                                                                 
  character(len=256) :: T_file

  ! Global, 4D array: (lon,lat,lev,time)                                                                                      
  ! ------------------------------------ 
  real, allocatable :: T(:,:,:,:)

  ! Miscellaneous                                                                                                              
  ! -------------                                                                                                              
  integer :: ierr
  integer :: im, jm, lm
  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"

  ! Allocate the Global 4-D array with only 1 time level                                                                      
  ! ----------------------------------------------------                                                                      
  allocate(T(im,jm,lm,1))

  ! Hypercube for reading one 1 array
  ! ---------------------------------
  start = (/  1,  1,  1, 1 /)
  count = (/ im, jm, lm, 1  /)   ! 1 time level, 3D (lon,lat,lev) array                                                         

  ! Read the data file                                                                                                         
  ! ------------------                                                                                                         
  write(*,*)'Reading T'
  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_get_var(ncid,varid,T,start=start,count=count), "reading T")
  call check(nf90_close(ncid), "closing T file")

  ! Print min/max of arrays                                                                                                    
  ! -----------------------                                                                                                    
  write(*,*)'T: ', maxval(T),minval(T)

  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

Compile and link

We use the utility nf-config (included in the NetCDF-4 installation) to identify the linker 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 T
T:    315.6512       180.3667

Access data via OPeNDAP server

Code

The code to read data from the OPeNDAP server is a minor modification (3 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

  ! Global, 4D array: (lon,lat,lev,time)
  ! ------------------------------------
  real, allocatable :: T(:,:,:,:)

  ! Miscellaneous
  ! -------------
  integer :: ierr
  integer :: im, jm, lm
  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"

  ! Allocate the Global 4-D array with only 1 time level
  ! ----------------------------------------------------
  allocate(T(im,jm,lm,1))

  ! Hypercube for reading one 1 array
  ! ---------------------------------
  start = (/  1,  1,  1, 11772 /) ! time level 11748 corresponds to 2006/09/18, 9z
  count = (/ im, jm, lm, 1  /)    ! 1 time level, 3D (lon,lat,lev) array

  ! Read the data file
  ! ------------------
  write(*,*)'Reading T'
  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_get_var(ncid,varid,T,start=start,count=count), "reading T")
  call check(nf90_close(ncid), "closing T file")

  ! Print min/max of arrays
  ! -----------------------
  write(*,*)'T: ', maxval(T),minval(T)

  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

The three (3) differences are

  1. T_file: instead of reading a downloaded NetCDF-4 file, we read a URL.
  2. start:
  3. nf90_inq_varid

Compile and run