Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

netcdf write optimizations #19

Closed
wants to merge 4 commits into from
Closed

netcdf write optimizations #19

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1265f04
netcdf write optimizations
jswhit2 Dec 5, 2019
735c54c
set chunksize to be 2d horizontal grid slice
jswhit2 Dec 5, 2019
8f0cd89
don't write _FillValue (use default to avoid _FillValue/datatype mism…
jswhit2 Dec 6, 2019
499950c
Merge remote-tracking branch 'upstream/develop' into develop
jswhit2 Dec 20, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 58 additions & 48 deletions io/module_write_netcdf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
character(len=ESMF_MAXSTR) :: attName, fldName

integer :: varival
real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax
real(4) :: varr4val, scale_fact, offset, dataMin, dataMax
real(4), allocatable, dimension(:) :: compress_err
real(8) :: varr8val
character(len=ESMF_MAXSTR) :: varcval

Expand All @@ -75,6 +76,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc

call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc)

allocate(compress_err(fieldCount)); compress_err=-999.
allocate(fldlev(fieldCount)) ; fldlev = 0
allocate(fcstField(fieldCount))
allocate(varids(fieldCount))
Expand Down Expand Up @@ -117,13 +119,14 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
if (mype==0) then

if (ideflate == 0) then
ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
ncerr = nf90_create(trim(filename), &
cmode=IOR(IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),NF90_SHARE), &
ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
else
ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), &
ncid=ncid); NC_ERR_STOP(ncerr)
! if compression on use HDF5 format with default _FillValue
ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
endif

! define dimensions
Expand Down Expand Up @@ -156,12 +159,14 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
! define variables
if (fldlev(i) == 1) then
if (typekind == ESMF_TYPEKIND_R4) then
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
if (ideflate > 0) then
! shuffle filter on for lossless compression
ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate)
NC_ERR_STOP(ncerr)
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,time_dimid/), varids(i), &
shuffle=.true.,deflate_level=ideflate,&
chunksizes=(/im,jm,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr)
else
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
endif
else if (typekind == ESMF_TYPEKIND_R8) then
ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
Expand All @@ -172,13 +177,15 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
end if
else if (fldlev(i) > 1) then
if (typekind == ESMF_TYPEKIND_R4) then
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
if (ideflate > 0) then
! shuffle filter off since lossy compression used
ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate)
NC_ERR_STOP(ncerr)
endif
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), &
shuffle=.false.,deflate_level=ideflate,&
chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Must the chunksizes be /im,jm,1,1/?

Having chunksizes of 1 performs badly. Having chunksizes that are the entire extent of the dimension can be fast for smaller sizes, but larger sizes should be broken up.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, I just set it to that because it seemed like the most common access pattern (reading a 2-d horizontal grid).

else
ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
(/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
endif
else if (typekind == ESMF_TYPEKIND_R8) then
ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
(/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
Expand Down Expand Up @@ -219,8 +226,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", &
name=trim(attName), value=varr8val, &
rc=rc); ESMF_ERR_RETURN(rc)
if (trim(attName) /= '_FillValue' .or. ideflate == 0) then
! FIXME: _FillValue must be cast to var type when using NF90_NETCDF4
if (trim(attName) /= '_FillValue') then
! FIXME: _FillValue must be cast to var type for recent versions of netcdf
ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr)
endif

Expand All @@ -236,6 +243,25 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc

end do ! i=1,fieldCount

! write grid_xt, grid_yt attributes
if (trim(output_grid) == 'gaussian_grid' .or. &
trim(output_grid) == 'regional_latlon') then
ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'rotated_latlon') then
ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'lambert_conformal') then
ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
endif

ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr)
end if

Expand All @@ -248,28 +274,16 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
if (trim(output_grid) == 'gaussian_grid' .or. &
trim(output_grid) == 'regional_latlon') then
ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(:,1) ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'rotated_latlon') then
do i=1,im
x(i) = lon1 + (lon2-lon1)/(im-1) * (i-1)
enddo
ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'lambert_conformal') then
do i=1,im
x(i) = dx * (i-1)
enddo
ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
endif
ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
endif
Expand All @@ -280,28 +294,16 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
if (trim(output_grid) == 'gaussian_grid' .or. &
trim(output_grid) == 'regional_latlon') then
ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:) ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'rotated_latlon') then
do j=1,jm
y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1)
enddo
ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
else if (trim(output_grid) == 'lambert_conformal') then
do j=1,jm
y(j) = dy * (j-1)
enddo
ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
endif
ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
endif
Expand Down Expand Up @@ -344,11 +346,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax)
! compute max abs compression error, save as a variable
! attribute.
compress_err = maxval(abs(arrayr4_3d_save-arrayr4_3d))
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d))
endif
ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/1,1,1/),count=(/im,jm,lm,1/) ); NC_ERR_STOP(ncerr)
end if
Expand All @@ -363,13 +361,25 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc

end do

if (ideflate > 0 .and. nbits > 0 .and. mype == 0) then
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
do i=1, fieldCount
if (compress_err(i) > 0) then
ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err(i)); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr)
endif
enddo
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
endif

deallocate(arrayr4)
deallocate(arrayr8)
deallocate(arrayr4_3d,arrayr4_3d_save)
deallocate(arrayr8_3d)

deallocate(fcstField)
deallocate(varids)
deallocate(compress_err)

if (mype==0) then
ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr)
Expand Down Expand Up @@ -484,9 +494,9 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc)
else if (typekind==ESMF_TYPEKIND_R8) then
call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
name=trim(attName), value=varr8val, rc=rc); ESMF_ERR_RETURN(rc)
if (trim(attName) /= '_FillValue' .or. ideflate == 0) then
! FIXME: _FillValue must be cast to var type when using
! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue.
if (trim(attName) /= '_FillValue') then
! FIXME: _FillValue must be cast to var type for recent versions
! of netcdf
ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr8val); NC_ERR_STOP(ncerr)
endif

Expand Down