Skip to content

Instantly share code, notes, and snippets.

@wjcapehart

wjcapehart/fnl_netcdf.f90

Created October 17, 2018 21:54
Show Gist options
  • Save wjcapehart/f777b5ae06c4d65ad4a561257aecd43c to your computer and use it in GitHub Desktop.
Save wjcapehart/f777b5ae06c4d65ad4a561257aecd43c to your computer and use it in GitHub Desktop.
Download ZIP
Fortran Program for Reading ASCII Model Data, and Writing [then Reading] NetCDF
Raw
fnl_netcdf.f90
program fnl_netcdf
! Takes two ascii fields of global data, Temp and Hgt @ 500 hPa
! And creates geostrophic u and v winds @ 500 hPa
!
! These fields are written as a netcdf file
! Build Instructions
! compile with:
! gfortran -I${NETCDFINC} -L{NETCDFLIB} -lnetcdff make_netcdf.f90
!
! executable will be named a.out unless you request a different name
! Geospatial information for the FNL data
! Data will read as (LON, LAT) with the LON (x) read in fastest
! NX = 360
! NY = 181
! DLON = +1 degree_east
! DLAT = -1 degree_north (it starts at 90 degrees_north!)
! LAT(001,001) = 90 N (!)
! LON(001,001) = 0 E
! LAT(360,180) = -90 N
! LON(360,180) = 355 E (i.e., 1 W)
! NZ(NP) = 1 (yes, use it, please! You'll see why!)
! PRESSURE = 500hPa (just one level)
! NT = 1
! Time = 0 seconds from 2007-03-01t00:00:00z
! Earth's Rad = 6,371,000 m (estimated by earth's volume)
!
! Fields and Units
! TMP = Air Temperature (K)
! UGRD = Grid-Relative Real Wind U-Component (m s-1)
! VGRD = Grid-Relative Real Wind V-Component (m s-1)
! HGT = Geopotential Height above Sea Level (m)
! Input files
! FNL_2007-03-01_00_TMP-500mb.txt
! FNL_2007-03-01_00_HGT-500mb.txt
! FNL_2007-03-01_00_RH-500mb.txt
! FNL_2007-03-01_00_UGRD-500mb.txt
! FNL_2007-03-01_00_VGRD-500mb.txt
! All of these files have one data value on each line.
! Things to be aware of
! We are working in NX,NY,NZ (or NP), NT as the array order
! but the netcdf file will reverse this order. This will be
! confusing.
!
!
! Also explore the GFS Dataset
! Geospatial information for the GFS data
! Data will read as (LON, LAT) with the LON (x) read in fastest
! NX = 720
! NY = 361
! DLON = +0.5 degree_east
! DLAT = -0.5 degree_north (it starts at 90 degrees_north!)
! LAT(001,001) = 90.0 N (!)
! LON(001,001) = 0.0 E
! LAT(360,180) = -90.0 N
! LON(360,180) = 355.5 E (i.e., 0.5 W)
! NZ(NP) = 1 (yes, use it, please! You'll see why!)
! PRESSURE = 500hPa (just one level)
! NT = 1
! Time = 0 seconds from 2009-03-30t12:00:00z
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
USE NETCDF ! the netcdf module is at ${NETCDFINC}/NETCDF.mod
IMPLICIT NONE
INTEGER, PARAMETER :: NLON = 360
INTEGER, PARAMETER :: NLAT = 181
INTEGER, PARAMETER :: NP = 1
INTEGER, PARAMETER :: NT = 1
CHARACTER (LEN=*), PARAMETER :: ORIG_FILEDIR = &
"/academic/atm_519/data/FNL_ASCII/"
!
! Coordinates for Writring
!
REAL, DIMENSION(NLON) :: longitude ! lambda : degree_east
REAL, DIMENSION(NLAT) :: latitude ! phi : degree_north
REAL, DIMENSION(NP) :: pressure ! p : Pa
REAL, DIMENSION(NT) :: time ! t : s'
!
! Variables for Writring
!
REAL, DIMENSION(NLON,NLAT,NP,NT) :: tmp ! T : K
REAL, DIMENSION(NLON,NLAT,NP,NT) :: hgt ! z : gpm
REAL, DIMENSION(NLON,NLAT,NP,NT) :: ugrd ! u : m s-2
REAL, DIMENSION(NLON,NLAT,NP,NT) :: vgrd ! v : m s-1
!
! Variables for Reading
!
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: tmp_2 ! T : K
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: hgt_2 ! z : gpm
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: ugrd_2 ! u : m s-2
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: vgrd_2 ! v : m s-1
!
! Geospatial Data
!
REAL, PARAMETER :: PI = 3.14159265359 ! π
REAL, PARAMETER :: RE = 6.371e6 ! Spherical Radius of the Earth (m)
REAL, PARAMETER :: LON_0 = 0. ! Longitude of first grid cell : degree_east
REAL, PARAMETER :: LAT_0 = 90. ! Latitude of first grid cell : degree_north
REAL, PARAMETER :: DLAT = -1 ! degree_north (starts at the north pole
REAL, PARAMETER :: DLON = 1 ! degree_east
REAL :: dx, dy ! dx and dy (calculated on the go)
INTEGER :: i, j, k, m ! we don't use "L" for obvious reasons!
REAL, PARAMETER :: FILL_VALUE = 9.96921e+36 ! Fill (or Missing) Value
INTEGER :: NLON_2, NLAT_2, NP_2, NT_2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! NetCDF Identifiers
integer :: ncstat ! generic netcdf status return variable
INTEGER :: netcdf_id_file ! netcdf file id
INTEGER :: netcdf_id_file_2 ! netcdf file id backend read
INTEGER :: netcdf_dim_lon ! netcdf lon (NX) dimension ID
INTEGER :: netcdf_dim_lat ! netcdf lat (NY) dimension ID
INTEGER :: netcdf_dim_pres ! netcdf pres (NZ) dimension ID
INTEGER :: netcdf_dim_time ! netcdf time (NT) dimension ID
INTEGER :: netcdf_dim_lon_2 ! netcdf lon (NX) dimension ID backend read
INTEGER :: netcdf_dim_lat_2 ! netcdf lat (NY) dimension ID backend read
INTEGER :: netcdf_dim_pres_2 ! netcdf pres (NZ) dimension ID backend read
INTEGER :: netcdf_dim_time_2 ! netcdf time (NT) dimension ID backend read
INTEGER :: netcdf_id_lon ! netcdf lon variable ID
INTEGER :: netcdf_id_lat ! netcdf lat variable ID
INTEGER :: netcdf_id_pres ! netcdf pres variable ID
INTEGER :: netcdf_id_time ! netcdf time variable ID
INTEGER :: netcdf_id_temp ! netcdf temperature variable ID
INTEGER :: netcdf_id_hgt ! netcdf geopot. height variable ID
INTEGER :: netcdf_id_ugrd ! netcdf grid-x wind variable ID
INTEGER :: netcdf_id_vgrd ! netcdf grid-y wind variable ID
INTEGER :: netcdf_id_temp_2 ! netcdf temperature variable ID backend read
INTEGER :: netcdf_id_hgt_2 ! netcdf geopot height variable ID backend rd
INTEGER :: netcdf_id_ugrd_2 ! netcdf grid-x wind variable ID backend read
INTEGER :: netcdf_id_vgrd_2 ! netcdf grid-y wind variable ID backend read
INTEGER, DIMENSION(4) :: netcdf_dims_4d ! NX,NY,NZ,NT netcdf dim bundle
INTEGER, DIMENSION(4) :: netcdf_dims_4d_start ! 1,1,1,1 array
INTEGER, DIMENSION(4) :: netcdf_dims_4d_end ! NX,NY,NZ,NT array
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! spatiotemporal information
!
DO i = 1, NLON
longitude(i) = LON_0 + (i-1) * DLON
END DO
DO j = 1, NLAT
latitude(j) = LAT_0 + (j-1) * DLAT
END DO
pressure(1) = 500.0e2 ! Pa
time(1) = 0 !
dy = 2. * RE * PI / 360.0 * DLAT ! 1 degree lat in meters
! dx = 2 * (RE*cos(lat*PI/180.)) * PI / 360.0 * DLON ! 1 degree lon in meters
!
! input of ascii grid data
!
! (Note that we have little-endian fortran binary files too in this dir!)
!
! These files have a .fbin as the suffix
!
! To read the fortran binary data change the "form" to "unformatted"
! and change the READ statemet from "READ(1,*) ..." to just READ(1)..."
!
OPEN(1, file=trim(ORIG_FILEDIR)//"FNL_2007-03-01_00_HGT-500mb.txt", form="formatted")
READ(1, *) hgt
CLOSE(1)
OPEN(1, file=trim(ORIG_FILEDIR)//"FNL_2007-03-01_00_TMP-500mb.txt", form="formatted")
READ(1, *) tmp
CLOSE(1)
OPEN(1, file=trim(ORIG_FILEDIR)//"FNL_2007-03-01_00_UGRD-500mb.txt", form="formatted")
READ(1, *) ugrd
CLOSE(1)
OPEN(1, file=trim(ORIG_FILEDIR)//"FNL_2007-03-01_00_VGRD-500mb.txt", form="formatted")
READ(1, *) vgrd
CLOSE(1)
print*, maxval(hgt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Encode NetCDF
!
! 1) open netcdf file
! 2) define the contests of the file
! 2a) create dimensions
! 2b) create variables and their attributes
! 2c) create global attributes if any
! 3) end definition mode
! 4) drop the data into the file
! 5) close netcdf file
!
! STOOPID PROGRAMER WARNING!!!
!
! For speed you will be tempted to do a lot of cutting and pasting.
! While there is nothing wrong with this in-and-of-itself, it can
! leave you vulnerable to error when you cut and paste large blocks of
! code and forget to modify everything.
!
!
!
! 1 -- Open NetCDF File and "clobber" any file that is currently in its spot
!
ncstat = NF90_CREATE(path = "./FNL_2007-03-01_00-500mb.nc", &
cmode = NF90_CLOBBER, &
ncid = netcdf_id_file)
! 2 -- Define Contents of File
!
! 2-a Create Dimensions
!
ncstat = NF90_DEF_DIM(ncid = netcdf_id_file, &
name = "time", &
len = NT, &
dimid = netcdf_dim_time)
ncstat = NF90_DEF_DIM(ncid = netcdf_id_file, &
name = "pressure", &
len = NP, &
dimid = netcdf_dim_pres)
ncstat = NF90_DEF_DIM(ncid = netcdf_id_file, &
name = "lon", &
len = NLON, &
dimid = netcdf_dim_lon)
ncstat = NF90_DEF_DIM(ncid = netcdf_id_file, &
name = "lat", &
len = NLAT, &
dimid = netcdf_dim_lat)
!
! 2-b (continued) For Multi-dimensional data we should create arrays
! for the data. Since we are working in fortran we do this
! nx,ny,nz,nt even though the netCDF libraries will write it
! nt,nz,ny,nx
!
! Here it is a good idea for you to create arrays for each multi-D
! combination (nx,ny), (nx,ny,nz), (nx,ny,nt), (nx,ny,nz,nt) that
! you will need. In this example we only have one: (nx,ny,nz,nt)
!
! For each combination you should have arrays
! - one with the netCSF *dimension* id "handles"
! - one with the starting point for your dimensions (all 1's)
! - one with the ending point for your dimensions (NX,NY,NZ,NT)
netcdf_dims_4d = (/ netcdf_dim_lon, &
netcdf_dim_lat, &
netcdf_dim_pres, &
netcdf_dim_time /)
netcdf_dims_4d_start = (/ 1, 1, 1, 1 /)
netcdf_dims_4d_end = (/ NLON, NLAT, NP, NT /)
!
! 2-b Create Variables and their Attributes
! (We do the data at a later step)
!
! I now like to start with my coordinates
!
! time
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "time", &
xtype = NF90_REAL, &
dimids = netcdf_dim_time, &
varid = netcdf_id_time)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_time, &
name = "long_name", &
values = "time")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_time, &
name = "standard_name", &
values = "time")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_time, &
name = "units", &
values = "seconds from 2007-03-01 00:00:00")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_time, &
name = "_CoordinateAxisType", &
values = "time")
! isobaric height
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "isobaric", &
xtype = NF90_REAL, &
dimids = netcdf_dim_pres, &
varid = netcdf_id_pres)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "long_name", &
values = "isobaric height")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "standard_name", &
values = "pressure")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "units", &
values = "Pa")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "_CoordinateAxisType", &
values = "Pressure")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "_CoordinateZisPositive", &
values = "down")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
name = "positive", &
values = "down")
! latitude
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "lat", &
xtype = NF90_REAL, &
dimids = netcdf_dim_lat, &
varid = netcdf_id_lat)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lat, &
name = "long_name", &
values = "latitude")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lat, &
name = "standard_name", &
values = "latitude")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lat, &
name = "units", &
values = "degree_north")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lat, &
name = "_CoordinateAxisType", &
values = "Lat")
! longitude
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "lon", &
xtype = NF90_REAL, &
dimids = netcdf_dim_lon, &
varid = netcdf_id_lon)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lon, &
name = "long_name", &
values = "longitude")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lon, &
name = "standard_name", &
values = "longitude")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lon, &
name = "units", &
values = "degree_east")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_lon, &
name = "_CoordinateAxisType", &
values = "Lon")
!
! 2-b (continued) Now lets do the actual data
!
! air temperature
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "temperature", &
xtype = NF90_REAL, &
dimids = netcdf_dims_4d, &
varid = netcdf_id_temp)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
name = "long_name", &
values = "temperature @ isobaric")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
name = "standard_name", &
values = "air_temperature")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
name = "units", &
values = "K")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
name = "coordinates", &
values = "time isobaric lat lon")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
name = "_FillValue", &
values = FILL_VALUE)
! geopotential_height
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "height", &
xtype = NF90_REAL, &
dimids = netcdf_dims_4d, &
varid = netcdf_id_hgt)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
name = "long_name", &
values = "geopotential height @ isobaric")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
name = "standard_name", &
values = "geopotential_height")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
name = "units", &
values = "m")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
name = "coordinates", &
values = "time isobaric lat lon")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
name = "_FillValue", &
values = FILL_VALUE)
! grid_eastward_wind
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "u", &
xtype = NF90_REAL, &
dimids = netcdf_dims_4d, &
varid = netcdf_id_ugrd)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
name = "long_name", &
values = "x-wind component @ isobaric")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
name = "standard_name", &
values = "grid_eastward_wind")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
name = "units", &
values = "m s-1")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
name = "coordinates", &
values = "time isobaric lat lon")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
name = "_FillValue", &
values = FILL_VALUE)
! grid_northward_wind
ncstat = NF90_DEF_VAR(ncid = netcdf_id_file, &
name = "v", &
xtype = NF90_REAL, &
dimids = netcdf_dims_4d, &
varid = netcdf_id_vgrd)
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
name = "long_name", &
values = "y-wind component @ isobaric")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
name = "standard_name", &
values = "grid_northward_wind")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
name = "units", &
values = "m s-1")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
name = "coordinates", &
values = "time isobaric lat lon")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
name = "_FillValue", &
values = FILL_VALUE)
!
! 2-c And now comes the global definitions
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = NF90_GLOBAL, &
name = "comment1", &
values = "My First NetCDF File")
ncstat = NF90_PUT_ATT(ncid = netcdf_id_file, &
varid = NF90_GLOBAL, &
name = "comment2", &
values = "Data Taken From NCAR Global FNL Analyses")
!
! 3 -- We now leave "definition mode." This closes off the part of the file
! that includes the header information. Once closed no more metatdata
! can be added.
!
ncstat = NF90_ENDDEF(ncid = netcdf_id_file)
!
! 4 -- And now we add the data. This includes the coordinates and the data,
! itself. The space allocated for the file has already been made
! and the data can be inserted in any order (but I like to keep it in
! the same order as I defined it to keep the code tidy - as you can
! see, this is messy enough!)
!
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_time, &
values = time, &
start = (/ 1 /) , &
count = (/ NT /) )
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_pres, &
values = pressure, &
start = (/ 1 /) , &
count = (/ NP /) )
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_lat, &
values = latitude, &
start = (/ 1 /) , &
count = (/ NLAT /) )
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_lon, &
values = longitude, &
start = (/ 1 /) , &
count = (/ NLON /) )
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_hgt, &
values = hgt, &
start = netcdf_dims_4d_start, &
count = netcdf_dims_4d_end)
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_temp, &
values = tmp, &
start = netcdf_dims_4d_start, &
count = netcdf_dims_4d_end)
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_ugrd, &
values = ugrd, &
start = netcdf_dims_4d_start, &
count = netcdf_dims_4d_end)
ncstat = NF90_PUT_VAR(ncid = netcdf_id_file, &
varid = netcdf_id_vgrd, &
values = vgrd, &
start = netcdf_dims_4d_start, &
count = netcdf_dims_4d_end)
!
! 5 -- Finally we close the file
!
ncstat = NF90_CLOSE(ncid = netcdf_id_file)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Now let's read that NetCDF
!
! 1) open netcdf file
! 2) read the variable
! 3) close netcdf file
!
!
! 1 -- Open the file for read-only status
!
ncstat = NF90_OPEN(path = "./FNL_2007-03-01_00-500mb.nc", &
mode = NF90_NOWRITE, &
ncid = netcdf_id_file_2)
!
! 2 -- Get Dimension Information (assuming we don't know our dimensions)
!
ncstat = NF90_INQ_DIMID(ncid = netcdf_id_file_2, &
name = "time", &
dimid = netcdf_dim_time_2)
ncstat = NF90_INQ_DIMID(ncid = netcdf_id_file_2, &
name = "pressure", &
dimid = netcdf_dim_pres_2)
ncstat = NF90_INQ_DIMID(ncid = netcdf_id_file_2, &
name = "lat", &
dimid = netcdf_dim_lat_2)
ncstat = NF90_INQ_DIMID(ncid = netcdf_id_file_2, &
name = "lon", &
dimid = netcdf_dim_lon_2)
ncstat = NF90_INQUIRE_DIMENSION(ncid = netcdf_id_file_2, &
dimid = netcdf_dim_time_2, &
len = NT_2)
ncstat = NF90_INQUIRE_DIMENSION(ncid = netcdf_id_file_2, &
dimid = netcdf_dim_pres_2, &
len = NP_2)
ncstat = NF90_INQUIRE_DIMENSION(ncid = netcdf_id_file_2, &
dimid = netcdf_dim_lat_2, &
len = NLAT_2)
ncstat = NF90_INQUIRE_DIMENSION(ncid = netcdf_id_file_2, &
dimid = netcdf_dim_lon_2, &
len = NLON_2)
!
! 3 -- allocate variables
!
WRITE(*,*) "ALLOCATING HEIGHT FOR DIMENSIONS", NLON_2, NLAT_2, NP_2, NT_2
ALLOCATE( hgt_2(NLON_2, NLAT_2, NP_2, NT_2) )
!
! 4 -- get the requested variable
!
ncstat = NF90_INQ_VARID(ncid = netcdf_id_file_2, &
name = "height", &
varid = netcdf_id_hgt_2)
ncstat = NF90_GET_VAR(ncid = netcdf_id_file_2, &
varid = netcdf_id_hgt_2, &
values = hgt_2)
!
! 5 -- close netcdf file
!
ncstat = NF90_CLOSE(ncid = netcdf_id_file_2)
!
! Now let's look at that variable
!
WRITE(*,*) "Statistics for Geopotential Height"
write(*,*) " Maximum: ", MINVAL(hgt_2)
write(*,*) " Average: ", (SUM(hgt_2) / ( NT_2 * NP_2 * NLON_2 * NLAT_2 ) )
write(*,*) " Minimum: ", MAXVAL(hgt_2)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! close up and clean up
!
DEALLOCATE( hgt_2 )
WRITE(*,*) "Task Complete!"
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end program fnl_netcdf
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment