#define ESMF_ERR_RETURN(rc) if (ESMF_LogFoundError(rc, msg="Breaking out of subroutine", line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT) #define NC_ERR_STOP(status) \ if (status /= nf90_noerr) write(0,*) "line ", __LINE__, trim(nf90_strerror(status)); \ if (status /= nf90_noerr) call ESMF_Finalize(endflag=ESMF_END_ABORT) module module_write_netcdf use esmf use netcdf use module_fv3_io_def,only : ideflate, nbits, & output_grid,dx,dy,lon1,lat1,lon2,lat2 implicit none private public write_netcdf contains !---------------------------------------------------------------------------------------- subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d, rc) ! type(ESMF_FieldBundle), intent(in) :: fieldbundle type(ESMF_FieldBundle), intent(in) :: wrtfb character(*), intent(in) :: filename integer, intent(in) :: mpi_comm integer, intent(in) :: mype integer, intent(in) :: im, jm integer, intent(in) :: ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d integer, optional,intent(out) :: rc ! !** local vars integer :: i,j,m,n,k integer :: lm integer, dimension(:), allocatable :: fldlev real(4), dimension(:,:), allocatable :: arrayr4 real(8), dimension(:,:), allocatable :: arrayr8 real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save real(8), dimension(:,:,:), allocatable :: arrayr8_3d real(8) x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound type(ESMF_Field), allocatable :: fcstField(:) type(ESMF_TypeKind_Flag) :: typekind type(ESMF_TypeKind_Flag) :: attTypeKind type(ESMF_Grid) :: wrtgrid type(ESMF_Array) :: array integer :: attcount character(len=ESMF_MAXSTR) :: attName, fldName integer :: varival real(4) :: varr4val, scale_fact, offset, dataMin, dataMax real(4), allocatable, dimension(:) :: compress_err real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval character(128) :: time_units integer :: ncerr integer :: ncid integer :: oldMode integer :: im_dimid, jm_dimid, pfull_dimid, phalf_dimid, time_dimid integer :: im_varid, jm_varid, lm_varid, time_varid, lon_varid, lat_varid integer, dimension(:), allocatable :: varids logical shuffle 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)) call ESMF_FieldBundleGet(fieldbundle, fieldList=fcstField, grid=wrtGrid, & ! itemorderflag=ESMF_ITEMORDER_ADDORDER, & rc=rc); ESMF_ERR_RETURN(rc) call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc); ESMF_ERR_RETURN(rc) do i=1,fieldCount call ESMF_FieldGet(fcstField(i), dimCount=fieldDimCount, rc=rc); ESMF_ERR_RETURN(rc) if (fieldDimCount > 3) then write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" stop end if if (fieldDimCount > gridDimCount) then allocate(ungriddedLBound(fieldDimCount-gridDimCount)) allocate(ungriddedUBound(fieldDimCount-gridDimCount)) call ESMF_FieldGet(fcstField(i), & ungriddedLBound=ungriddedLBound, & ungriddedUBound=ungriddedUBound, rc=rc); ESMF_ERR_RETURN(rc) fldlev(i) = ungriddedUBound(fieldDimCount-gridDimCount) - & ungriddedLBound(fieldDimCount-gridDimCount) + 1 deallocate(ungriddedLBound) deallocate(ungriddedUBound) else if (fieldDimCount == 2) then fldlev(i) = 1 end if end do lm = maxval(fldlev(:)) allocate(arrayr4(im,jm)) allocate(arrayr8(im,jm)) allocate(arrayr4_3d(im,jm,lm),arrayr4_3d_save(im,jm,lm)) allocate(arrayr8_3d(im,jm,lm)) ! create netcdf file and enter define mode if (mype==0) then ncerr = nf90_create(trim(filename),& cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),& ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) ! define dimensions ncerr = nf90_def_dim(ncid, "grid_xt", im, im_dimid); NC_ERR_STOP(ncerr) ncerr = nf90_def_dim(ncid, "grid_yt", jm, jm_dimid); NC_ERR_STOP(ncerr) ! define coordinate variables ncerr = nf90_def_var(ncid, "grid_xt", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr) ncerr = nf90_def_var(ncid, "lon", NF90_DOUBLE, (/im_dimid,jm_dimid/), lon_varid); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, lon_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, lon_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, im_varid, "cartesian_axis", "X"); NC_ERR_STOP(ncerr) ncerr = nf90_def_var(ncid, "grid_yt", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr) ncerr = nf90_def_var(ncid, "lat", NF90_DOUBLE, (/im_dimid,jm_dimid/), lat_varid); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, lat_varid, "long_name", "T-cell latitude"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, lat_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, jm_varid, "cartesian_axis", "Y"); NC_ERR_STOP(ncerr) if (lm > 1) then call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, rc) call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, rc) end if call add_dim(ncid, "time", time_dimid, wrtgrid, rc) call get_global_attr(wrtfb, ncid, rc) do i=1, fieldCount call ESMF_FieldGet(fcstField(i), name=fldName, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) ! define variables if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then if (ideflate > 0) then if (ichunk2d < 0 .or. jchunk2d < 0) then ! let netcdf lib choose chunksize ! shuffle filter on for 2d fields (lossless compression) ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,time_dimid/), varids(i), & shuffle=.true.,deflate_level=ideflate); NC_ERR_STOP(ncerr) else ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,time_dimid/), varids(i), & shuffle=.true.,deflate_level=ideflate,& chunksizes=(/ichunk2d,jchunk2d,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) endif 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, & (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) else write(0,*)'Unsupported typekind ', typekind stop end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then if (ideflate > 0) then ! shuffle filter off for 3d fields using lossy compression if (nbits > 0) then shuffle=.false. else shuffle=.true. endif if (ichunk3d < 0 .or. jchunk3d < 0 .or. kchunk3d < 0) then ! let netcdf lib choose chunksize ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & shuffle=shuffle,deflate_level=ideflate); NC_ERR_STOP(ncerr) else ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & shuffle=shuffle,deflate_level=ideflate,& chunksizes=(/ichunk3d,jchunk3d,kchunk3d,1/)); NC_ERR_STOP(ncerr) endif 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) else write(0,*)'Unsupported typekind ', typekind stop end if end if ! define variable attributes call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, Count=attcount, & rc=rc); ESMF_ERR_RETURN(rc) do j=1,attCount call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, attributeIndex=j, & name=attName, typekind=attTypeKind, itemCount=n, & rc=rc); ESMF_ERR_RETURN(rc) if ( index(trim(attName),"ESMF") /= 0 ) then cycle endif if (attTypeKind==ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varival, & rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varids(i), trim(attName), varival); NC_ERR_STOP(ncerr) else if (attTypeKind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr4val, & rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr4val); NC_ERR_STOP(ncerr) else if (attTypeKind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, & rc=rc); ESMF_ERR_RETURN(rc) 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 else if (attTypeKind==ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varcval, & rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varids(i), trim(attName), trim(varcval)); NC_ERR_STOP(ncerr) end if end do ! j=1,attCount 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 ! end of define mode ! write grid_xt, grid_yt values call ESMF_GridGetCoord(wrtGrid, coordDim=1, array=array, rc=rc); ESMF_ERR_RETURN(rc) call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); 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) 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) endif ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc) call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); 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) 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) endif ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif do i=1, fieldCount call ESMF_FieldGet(fcstField(i),name=fldName,typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then call ESMF_FieldGather(fcstField(i), arrayr4, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then ncerr = nf90_put_var(ncid, varids(i), values=arrayr4, start=(/1,1,1/),count=(/im,jm,1/) ); NC_ERR_STOP(ncerr) end if else if (typekind == ESMF_TYPEKIND_R8) then call ESMF_FieldGather(fcstField(i), arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then ncerr = nf90_put_var(ncid, varids(i), values=arrayr8, start=(/1,1,1/),count=(/im,jm,1/) ); NC_ERR_STOP(ncerr) end if end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then call ESMF_FieldGather(fcstField(i), arrayr4_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then if (ideflate > 0 .and. nbits > 0) then ! Lossy compression if nbits>0. ! The floating point data is quantized to improve compression ! See doi:10.5194/gmd-10-413-2017. The method employed ! here is identical to the 'scaled linear packing' method in ! that paper, except that the data are scaling into an arbitrary ! range (2**nbits-1 not just 2**16-1) and are stored as ! re-scaled floats instead of short integers. ! The zlib algorithm does almost as ! well packing the re-scaled floats as it does the scaled ! integers, and this avoids the need for the client to apply the ! rescaling (plus it allows the ability to adjust the packing ! range). arrayr4_3d_save = arrayr4_3d dataMax = maxval(arrayr4_3d); dataMin = minval(arrayr4_3d) arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax) ! compute max abs compression error, save as a variable ! attribute. 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 else if (typekind == ESMF_TYPEKIND_R8) then call ESMF_FieldGather(fcstField(i), arrayr8_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (mype==0) then ncerr = nf90_put_var(ncid, varids(i), values=arrayr8_3d, start=(/1,1,1/),count=(/im,jm,lm,1/) ); NC_ERR_STOP(ncerr) end if end if end if 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) end if end subroutine write_netcdf !---------------------------------------------------------------------------------------- subroutine get_global_attr(fldbundle, ncid, rc) type(ESMF_FieldBundle), intent(in) :: fldbundle integer, intent(in) :: ncid integer, intent(out) :: rc ! local variable integer :: i, attcount integer :: ncerr character(len=ESMF_MAXSTR) :: attName type(ESMF_TypeKind_Flag) :: typekind integer :: varival real(ESMF_KIND_R4) :: varr4val real(ESMF_KIND_R4), dimension(:), allocatable :: varr4list real(ESMF_KIND_R8) :: varr8val real(ESMF_KIND_R8), dimension(:), allocatable :: varr8list integer :: itemCount character(len=ESMF_MAXSTR) :: varcval ! call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, Count=attcount, & rc=rc); ESMF_ERR_RETURN(rc) do i=1,attCount call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, attributeIndex=i, name=attName, & typekind=typekind, itemCount=itemCount, rc=rc); ESMF_ERR_RETURN(rc) if (typekind==ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varival, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), varival); NC_ERR_STOP(ncerr) else if (typekind==ESMF_TYPEKIND_R4) then allocate (varr4list(itemCount)) call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & name=trim(attName), valueList=varr4list, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), varr4list); NC_ERR_STOP(ncerr) deallocate(varr4list) else if (typekind==ESMF_TYPEKIND_R8) then allocate (varr8list(itemCount)) call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & name=trim(attName), valueList=varr8list, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), varr8list); NC_ERR_STOP(ncerr) deallocate(varr8list) else if (typekind==ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(fldbundle, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varcval, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, NF90_GLOBAL, trim(attName), trim(varcval)); NC_ERR_STOP(ncerr) end if end do end subroutine get_global_attr ! !---------------------------------------------------------------------------------------- subroutine get_grid_attr(grid, prefix, ncid, varid, rc) type(ESMF_Grid), intent(in) :: grid character(len=*), intent(in) :: prefix integer, intent(in) :: ncid integer, intent(in) :: varid integer, intent(out) :: rc ! local variable integer :: i, attcount, n, ind integer :: ncerr character(len=ESMF_MAXSTR) :: attName type(ESMF_TypeKind_Flag) :: typekind integer :: varival real(ESMF_KIND_R4) :: varr4val real(ESMF_KIND_R8) :: varr8val character(len=ESMF_MAXSTR) :: varcval ! call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, Count=attcount, & rc=rc); ESMF_ERR_RETURN(rc) !write(0,*)'grid attcount = ', attcount do i=1,attCount call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, attributeIndex=i, name=attName, & typekind=typekind, itemCount=n, rc=rc); ESMF_ERR_RETURN(rc) !write(0,*)'grid att = ',i,trim(attName), ' itemCount = ' , n if (index(trim(attName), trim(prefix)//":")==1) then ind = len(trim(prefix)//":") if (typekind==ESMF_TYPEKIND_I4) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varival, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varival); NC_ERR_STOP(ncerr) else if (typekind==ESMF_TYPEKIND_R4) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr4val, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr4val); NC_ERR_STOP(ncerr) 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') 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 else if (typekind==ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varcval, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), trim(varcval)); NC_ERR_STOP(ncerr) end if end if end do end subroutine get_grid_attr subroutine add_dim(ncid, dim_name, dimid, grid, rc) integer, intent(in) :: ncid character(len=*), intent(in) :: dim_name integer, intent(inout) :: dimid type(ESMF_Grid), intent(in) :: grid integer, intent(out) :: rc ! local variable integer :: i, attcount, n, dim_varid integer :: ncerr character(len=ESMF_MAXSTR) :: attName type(ESMF_TypeKind_Flag) :: typekind integer, allocatable :: valueListI(:) real(ESMF_KIND_R4), allocatable :: valueListR4(:) real(ESMF_KIND_R8), allocatable :: valueListR8(:) character(len=ESMF_MAXSTR), allocatable :: valueListC(:) ! call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & attnestflag=ESMF_ATTNEST_OFF, name=dim_name, & typekind=typekind, itemCount=n, rc=rc); ESMF_ERR_RETURN(rc) if ( trim(dim_name) == "time" ) then ncerr = nf90_def_dim(ncid, trim(dim_name), NF90_UNLIMITED, dimid); NC_ERR_STOP(ncerr) else ncerr = nf90_def_dim(ncid, trim(dim_name), n, dimid); NC_ERR_STOP(ncerr) end if if (typekind==ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, dim_name, NF90_REAL8, dimids=(/dimid/), varid=dim_varid); NC_ERR_STOP(ncerr) allocate(valueListR8(n)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dim_name), valueList=valueListR8, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR8 ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR8) else if (typekind==ESMF_TYPEKIND_R4) then ncerr = nf90_def_var(ncid, dim_name, NF90_REAL4, dimids=(/dimid/), varid=dim_varid); NC_ERR_STOP(ncerr) allocate(valueListR4(n)) call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(dim_name), valueList=valueListR4, rc=rc); ESMF_ERR_RETURN(rc) ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4 ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR4) else write(0,*)'Error in module_write_netcdf.F90(add_dim) unknown typekind for ',trim(dim_name) call ESMF_Finalize(endflag=ESMF_END_ABORT) end if call get_grid_attr(grid, dim_name, ncid, dim_varid, rc) end subroutine add_dim ! !---------------------------------------------------------------------------------------- subroutine nccheck(status) use netcdf implicit none integer, intent (in) :: status if (status /= nf90_noerr) then write(0,*) status, trim(nf90_strerror(status)) stop "stopped" end if end subroutine nccheck elemental real function quantized(dataIn, nbits, dataMin, dataMax) integer, intent(in) :: nbits real(4), intent(in) :: dataIn, dataMin, dataMax real(4) offset, scale_fact ! convert data to 32 bit integers in range 0 to 2**nbits-1, then cast ! cast back to 32 bit floats (data is then quantized in steps ! proportional to 2**nbits so last 32-nbits in floating ! point representation should be zero for efficient zlib compression). scale_fact = (dataMax - dataMin) / (2**nbits-1); offset = dataMin quantized = scale_fact*(nint((dataIn - offset) / scale_fact)) + offset end function quantized end module module_write_netcdf