diff --git a/CMakeLists.txt b/CMakeLists.txt index ec9721ba6..718ba11b4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,7 +43,6 @@ add_library(fv3atm cpl/module_cap_cpl.F90 io/FV3GFS_io.F90 io/module_write_netcdf.F90 - io/module_write_netcdf_parallel.F90 io/module_fv3_io_def.F90 io/module_write_internal_state.F90 io/module_wrt_grid_comp.F90 diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 30959e625..f7f518f3a 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -1,7 +1,8 @@ -#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 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) write(0,*) "file: ", __FILE__, " line: ", __LINE__, trim(nf90_strerror(status)); \ if (status /= nf90_noerr) call ESMF_Finalize(endflag=ESMF_END_ABORT) module module_write_netcdf @@ -10,46 +11,63 @@ module module_write_netcdf use netcdf use module_fv3_io_def,only : ideflate, nbits, & output_grid,dx,dy,lon1,lat1,lon2,lat2 + use mpi implicit none private public write_netcdf + logical :: par + + interface quantize_array + module procedure quantize_array_3d + module procedure quantize_array_4d + end interface + contains !---------------------------------------------------------------------------------------- - subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d, rc) + subroutine write_netcdf(wrtfb, filename, & + use_parallel_netcdf, mpi_comm, mype, & + ichunk2d, jchunk2d, ichunk3d, jchunk3d, kchunk3d, rc) ! - type(ESMF_FieldBundle), intent(in) :: fieldbundle type(ESMF_FieldBundle), intent(in) :: wrtfb character(*), intent(in) :: filename + logical, intent(in) :: use_parallel_netcdf 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 :: i,j,t, istart,iend,jstart,jend + integer :: im, jm, lm + + integer, dimension(:), allocatable :: fldlev + + real(ESMF_KIND_R4), dimension(:,:), pointer :: array_r4 + real(ESMF_KIND_R4), dimension(:,:,:), pointer :: array_r4_cube + real(ESMF_KIND_R4), dimension(:,:,:), pointer :: array_r4_3d + real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: array_r4_3d_cube - 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(ESMF_KIND_R8), dimension(:,:), pointer :: array_r8 + real(ESMF_KIND_R8), dimension(:,:,:), pointer :: array_r8_cube + real(ESMF_KIND_R8), dimension(:,:,:), pointer :: array_r8_3d + real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: array_r8_3d_cube - real(8) x(im),y(jm) + real(8), dimension(:), allocatable :: x,y integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound + integer, dimension(:), allocatable :: start_idx type(ESMF_Field), allocatable :: fcstField(:) type(ESMF_TypeKind_Flag) :: typekind type(ESMF_TypeKind_Flag) :: attTypeKind type(ESMF_Grid) :: wrtgrid type(ESMF_Array) :: array + type(ESMF_DistGrid) :: distgrid - integer :: attcount + integer :: attCount character(len=ESMF_MAXSTR) :: attName, fldName integer :: varival @@ -58,35 +76,92 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ic real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval - character(128) :: time_units - - integer :: ncerr + integer :: ncerr,ierr 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 :: im_dimid, jm_dimid, tile_dimid, pfull_dimid, phalf_dimid, time_dimid + integer :: im_varid, jm_varid, tile_varid, lon_varid, lat_varid + integer, dimension(:), allocatable :: dimids_2d, dimids_3d integer, dimension(:), allocatable :: varids logical shuffle - call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) + logical :: is_cubed_sphere + integer :: rank, deCount, localDeCount, dimCount, tileCount + integer :: my_tile, start_i, start_j + integer, dimension(:,:), allocatable :: minIndexPDe, maxIndexPDe + integer, dimension(:,:), allocatable :: minIndexPTile, maxIndexPTile + integer, dimension(:), allocatable :: deToTileMap, localDeToDeMap + logical :: do_io + integer :: par_access +! + is_cubed_sphere = .false. + tileCount = 0 + my_tile = 0 + start_i = -10000000 + start_j = -10000000 + + par = use_parallel_netcdf + do_io = par .or. (mype==0) + + call ESMF_FieldBundleGet(wrtfb, 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, & + call ESMF_FieldBundleGet(wrtfb, 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) + call ESMF_FieldGet(fcstField(i), dimCount=fieldDimCount, array=array, rc=rc); ESMF_ERR_RETURN(rc) + if (fieldDimCount > 3) then write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" - stop + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! use first field to determine tile number, grid size, start index etc. + if (i == 1) then + call ESMF_ArrayGet(array, & + distgrid=distgrid, & + dimCount=dimCount, & + deCount=deCount, & + localDeCount=localDeCount, & + tileCount=tileCount, & + rc=rc); ESMF_ERR_RETURN(rc) + + allocate(minIndexPDe(dimCount,deCount)) + allocate(maxIndexPDe(dimCount,deCount)) + allocate(minIndexPTile(dimCount, tileCount)) + allocate(maxIndexPTile(dimCount, tileCount)) + call ESMF_DistGridGet(distgrid, & + minIndexPDe=minIndexPDe, maxIndexPDe=maxIndexPDe, & + minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, & + rc=rc); ESMF_ERR_RETURN(rc) + + allocate(deToTileMap(deCount)) + allocate(localDeToDeMap(localDeCount)) + call ESMF_ArrayGet(array, & + deToTileMap=deToTileMap, & + localDeToDeMap=localDeToDeMap, & + rc=rc); ESMF_ERR_RETURN(rc) + + is_cubed_sphere = (tileCount == 6) + my_tile = deToTileMap(localDeToDeMap(1)+1) + im = maxIndexPTile(1,1) + jm = maxIndexPTile(2,1) + start_i = minIndexPDe(1,localDeToDeMap(1)+1) + start_j = minIndexPDe(2,localDeToDeMap(1)+1) + if (.not. par) then + start_i = 1 + start_j = 1 + end if end if + if (fieldDimCount > gridDimCount) then allocate(ungriddedLBound(fieldDimCount-gridDimCount)) allocate(ungriddedUBound(fieldDimCount-gridDimCount)) @@ -104,301 +179,559 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ic 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) + ! for serial output allocate 'global' arrays + if (.not. par) then + allocate(array_r4(im,jm)) + allocate(array_r8(im,jm)) + allocate(array_r4_3d(im,jm,lm)) + allocate(array_r8_3d(im,jm,lm)) + if (is_cubed_sphere) then + allocate(array_r4_cube(im,jm,tileCount)) + allocate(array_r8_cube(im,jm,tileCount)) + allocate(array_r4_3d_cube(im,jm,lm,tileCount)) + allocate(array_r8_3d_cube(im,jm,lm,tileCount)) + end if end if - call add_dim(ncid, "time", time_dimid, wrtgrid, rc) + ! create netcdf file and enter define mode + if (do_io) then + + if (par) then + ncerr = nf90_create(trim(filename),& + cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),& + comm=mpi_comm, info = MPI_INFO_NULL, ncid=ncid); 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) + end if - call get_global_attr(wrtfb, ncid, rc) + ! disable auto filling. + ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) - 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. + ! define dimensions [grid_xt, grid_yta ,(pfull/phalf), (tile), time] + 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) + if (lm > 1) then + call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, rc) + call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, rc) + end if + if (is_cubed_sphere) then + ncerr = nf90_def_dim(ncid, "tile", tileCount, tile_dimid); NC_ERR_STOP(ncerr) + end if + call add_dim(ncid, "time", time_dimid, wrtgrid, rc) + + ! define coordinate variables + ncerr = nf90_def_var(ncid, "grid_xt", NF90_DOUBLE, im_dimid, im_varid); 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_put_att(ncid, jm_varid, "cartesian_axis", "Y"); NC_ERR_STOP(ncerr) + if (is_cubed_sphere) then + ncerr = nf90_def_var(ncid, "tile", NF90_INT, tile_dimid, tile_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, tile_varid, "long_name", "cubed-spehere face"); NC_ERR_STOP(ncerr) + end if + + ! coordinate variable attributes based on output_grid type + if (trim(output_grid) == 'gaussian_grid' .or. & + trim(output_grid) == 'global_latlon' .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) + end if + + ! define longitude variable + if (is_cubed_sphere) then + ncerr = nf90_def_var(ncid, "lon", NF90_DOUBLE, [im_dimid,jm_dimid,tile_dimid], lon_varid); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, "lon", NF90_DOUBLE, [im_dimid,jm_dimid ], lon_varid); NC_ERR_STOP(ncerr) + end if + 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) + + ! define latitude variable + if (is_cubed_sphere) then + ncerr = nf90_def_var(ncid, "lat", NF90_DOUBLE, [im_dimid,jm_dimid,tile_dimid], lat_varid); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, "lat", NF90_DOUBLE, [im_dimid,jm_dimid ], lat_varid); NC_ERR_STOP(ncerr) + end if + 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) + + if (par) then + ncerr = nf90_var_par_access(ncid, im_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + ncerr = nf90_var_par_access(ncid, lon_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + ncerr = nf90_var_par_access(ncid, jm_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + ncerr = nf90_var_par_access(ncid, lat_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + if (is_cubed_sphere) then + ncerr = nf90_var_par_access(ncid, tile_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + end if + end if + + + call get_global_attr(wrtfb, ncid, rc) + + + ! define variables (fields) + if (is_cubed_sphere) then + allocate(dimids_2d(4)) + allocate(dimids_3d(5)) + dimids_2d = [im_dimid,jm_dimid, tile_dimid,time_dimid] + if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid,tile_dimid,time_dimid] + else + allocate(dimids_2d(3)) + allocate(dimids_3d(4)) + dimids_2d = [im_dimid,jm_dimid, time_dimid] + if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid, time_dimid] + end if + + do i=1, fieldCount + call ESMF_FieldGet(fcstField(i), name=fldName, rank=rank, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + + par_access = NF90_INDEPENDENT + ! define variables + if (rank == 2) 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, & + dimids_2d, varids(i), & + shuffle=.true.,deflate_level=ideflate); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_2d, varids(i), & + shuffle=.true.,deflate_level=ideflate,& + chunksizes=[ichunk2d,jchunk2d,tileCount,1]); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_2d, varids(i), & + shuffle=.true.,deflate_level=ideflate,& + chunksizes=[ichunk2d,jchunk2d, 1]); NC_ERR_STOP(ncerr) + end if + end if + ! compression filters require collective access. + par_access = NF90_COLLECTIVE 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) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_2d, varids(i)); NC_ERR_STOP(ncerr) + end if + else if (typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & + dimids_2d, varids(i)); NC_ERR_STOP(ncerr) + else + write(0,*)'Unsupported typekind ', typekind + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + else if (rank == 3) 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. + end if + if (ichunk3d < 0 .or. jchunk3d < 0 .or. kchunk3d < 0) then + ! let netcdf lib choose chunksize + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_3d, varids(i), & + shuffle=shuffle,deflate_level=ideflate); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_3d, varids(i), & + shuffle=shuffle,deflate_level=ideflate,& + chunksizes=[ichunk3d,jchunk3d,kchunk3d,tileCount,1]); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_3d, varids(i), & + shuffle=shuffle,deflate_level=ideflate,& + chunksizes=[ichunk3d,jchunk3d,kchunk3d, 1]); NC_ERR_STOP(ncerr) + end if + end if + ! compression filters require collective access. + par_access = NF90_COLLECTIVE 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 + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + dimids_3d, varids(i)); NC_ERR_STOP(ncerr) + end if + else if (typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & + dimids_3d, varids(i)); NC_ERR_STOP(ncerr) 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) + write(0,*)'Unsupported typekind ', typekind + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + else + write(0,*)'Unsupported rank ', rank + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + if (par) then + ncerr = nf90_var_par_access(ncid, varids(i), par_access); NC_ERR_STOP(ncerr) + end if - if ( index(trim(attName),"ESMF") /= 0 ) then - cycle - endif + ! define variable attributes + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + attnestflag=ESMF_ATTNEST_OFF, count=attCount, & + rc=rc); ESMF_ERR_RETURN(rc) - if (attTypeKind==ESMF_TYPEKIND_I4) then + do j=1,attCount call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & - name=trim(attName), value=varival, & + attnestflag=ESMF_ATTNEST_OFF, attributeIndex=j, & + name=attName, typekind=attTypeKind, & 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) + if (index(trim(attName),"ESMF") /= 0) then + cycle + end if - 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 + 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 when using NF90_NETCDF4 + ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr) + end if + + 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) - 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) == 'global_latlon' .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 do ! j=1,attCount + + if (is_cubed_sphere) then + ncerr = nf90_put_att(ncid, varids(i), 'coordinates', 'lon lat'); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'grid_mapping', 'cubed_sphere'); NC_ERR_STOP(ncerr) + end if + + end do ! i=1,fieldCount + + ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr) + end if + ! end of define mode + + ! + ! write dimension variables and lon,lat variables + ! + if (allocated(start_idx)) deallocate(start_idx) + if (is_cubed_sphere) then + allocate(start_idx(3)) + start_idx = [start_i, start_j, my_tile] + else + allocate(start_idx(2)) + start_idx = [start_i, start_j] end if -! end of define mode + ! write lon (lon_varid) + if (par) then + call ESMF_GridGetCoord(wrtGrid, coordDim=1, farrayPtr=array_r8, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_var(ncid, lon_varid, values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + else + call ESMF_GridGetCoord(wrtGrid, coordDim=1, array=array, rc=rc); ESMF_ERR_RETURN(rc) + if (is_cubed_sphere) then + do t=1,tileCount + call ESMF_ArrayGather(array, array_r8_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (do_io) then + ncerr = nf90_put_var(ncid, lon_varid, values=array_r8_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_ArrayGather(array, array_r8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, lon_varid, values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + end if + endif + end if - ! 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 + istart = lbound(array_r8,1); iend = ubound(array_r8,1) + jstart = lbound(array_r8,2); jend = ubound(array_r8,2) + + ! write grid_xt (im_varid) + if (do_io) then + allocate (x(im)) if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'global_latlon' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, im_varid, values=array_r8(:,jstart), start=[istart], count=[iend-istart+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) + end do + 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) + end do + ncerr = nf90_put_var(ncid, im_varid, values=x); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'cubed_sphere_grid') then + do i=1,im + x(i) = i + end do + ncerr = nf90_put_var(ncid, im_varid, values=x); NC_ERR_STOP(ncerr) + else + write(0,*)'unknown output_grid ', trim(output_grid) + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + + ! write lat (lat_varid) + if (par) then + call ESMF_GridGetCoord(wrtGrid, coordDim=2, farrayPtr=array_r8, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_var(ncid, lat_varid, values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + else + call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc) + if (is_cubed_sphere) then + do t=1,tileCount + call ESMF_ArrayGather(array, array_r8_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (do_io) then + ncerr = nf90_put_var(ncid, lat_varid, values=array_r8_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_ArrayGather(array, array_r8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, lat_varid, values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + end if endif - ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) - endif + end if - 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 + ! write grid_yt (jm_varid) + if (do_io) then + allocate (y(jm)) if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'global_latlon' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, jm_varid, values=array_r8(istart,:), start=[jstart], count=[jend-jstart+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) + end do + 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 + end do + ncerr = nf90_put_var(ncid, jm_varid, values=y); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'cubed_sphere_grid') then + do j=1,jm + y(j) = j + end do + ncerr = nf90_put_var(ncid, jm_varid, values=y); NC_ERR_STOP(ncerr) + else + write(0,*)'unknown output_grid ', trim(output_grid) + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + end if + + ! write tile (tile_varid) + if (do_io .and. is_cubed_sphere) then + ncerr = nf90_put_var(ncid, tile_varid, values=[1,2,3,4,5,6]); NC_ERR_STOP(ncerr) + end if + ! write variables (fields) do i=1, fieldCount - call ESMF_FieldGet(fcstField(i),name=fldName,typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + call ESMF_FieldGet(fcstField(i),name=fldName,rank=rank,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 + if (rank == 2) then + + if (allocated(start_idx)) deallocate(start_idx) + if (is_cubed_sphere) then + allocate(start_idx(4)) + start_idx = [start_i,start_j,my_tile,1] + else + allocate(start_idx(3)) + start_idx = [start_i,start_j, 1] 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 + if (par) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r4, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_var(ncid, varids(i), values=array_r4, start=start_idx); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + do t=1,tileCount + call ESMF_ArrayGather(array, array_r4_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_FieldGather(fcstField(i), array_r4, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r4, start=start_idx); NC_ERR_STOP(ncerr) + end if + end if + 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 + if (par) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r8, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_var(ncid, varids(i), values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + do t=1,tileCount + call ESMF_ArrayGather(array, array_r8_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_FieldGather(fcstField(i), array_r8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (do_io) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8, start=start_idx); NC_ERR_STOP(ncerr) + end if + end if + end if end if - end if + else if (rank == 3) then - end do + if (allocated(start_idx)) deallocate(start_idx) + if (is_cubed_sphere) then + allocate(start_idx(5)) + start_idx = [start_i,start_j,1,my_tile,1] + else + allocate(start_idx(4)) + start_idx = [start_i,start_j,1, 1] + end if - if (ideflate > 0 .and. nbits > 0 .and. mype == 0) then + if (typekind == ESMF_TYPEKIND_R4) then + if (par) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r4_3d, rc=rc); ESMF_ERR_RETURN(rc) + if (ideflate > 0 .and. nbits > 0) then + dataMax = maxval(array_r4_3d) + dataMin = minval(array_r4_3d) + call mpi_allreduce(mpi_in_place,dataMax,1,mpi_real4,mpi_max,mpi_comm,ierr) + call mpi_allreduce(mpi_in_place,dataMin,1,mpi_real4,mpi_min,mpi_comm,ierr) + call quantize_array(array_r4_3d, dataMin, dataMax, nbits, compress_err(i)) + call mpi_allreduce(mpi_in_place,compress_err(i),1,mpi_real4,mpi_max,mpi_comm,ierr) + end if + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d, start=start_idx); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + do t=1,tileCount + call ESMF_ArrayGather(array, array_r4_3d_cube(:,:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (mype==0) then + if (ideflate > 0 .and. nbits > 0) then + call quantize_array(array_r4_3d_cube, minval(array_r4_3d_cube), maxval(array_r4_3d_cube), nbits, compress_err(i)) + end if + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_FieldGather(fcstField(i), array_r4_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (mype==0) then + if (ideflate > 0 .and. nbits > 0) then + call quantize_array(array_r4_3d, minval(array_r4_3d), maxval(array_r4_3d), nbits, compress_err(i)) + end if + ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d, start=start_idx); NC_ERR_STOP(ncerr) + end if + end if + end if + else if (typekind == ESMF_TYPEKIND_R8) then + if (par) then + call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=array_r8_3d, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_3d, start=start_idx); NC_ERR_STOP(ncerr) + else + if (is_cubed_sphere) then + call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc) + do t=1,tileCount + call ESMF_ArrayGather(array, array_r8_3d_cube(:,:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc) + end do + if (mype==0) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_3d_cube, start=start_idx); NC_ERR_STOP(ncerr) + end if + else + call ESMF_FieldGather(fcstField(i), array_r8_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (mype==0) then + ncerr = nf90_put_var(ncid, varids(i), values=array_r8_3d, start=start_idx); NC_ERR_STOP(ncerr) + end if + end if + end if + end if ! end typekind + + else + + write(0,*)'Unsupported rank ', rank + call ESMF_Finalize(endflag=ESMF_END_ABORT) + + end if ! end rank + + end do ! end fieldCount + + if (ideflate > 0 .and. nbits > 0 .and. do_io) 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 + end if + end do ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) - endif + end if - deallocate(arrayr4) - deallocate(arrayr8) - deallocate(arrayr4_3d,arrayr4_3d_save) - deallocate(arrayr8_3d) + if (.not. par) then + deallocate(array_r4) + deallocate(array_r8) + deallocate(array_r4_3d) + deallocate(array_r8_3d) + if (is_cubed_sphere) then + deallocate(array_r4_cube) + deallocate(array_r8_cube) + deallocate(array_r4_3d_cube) + deallocate(array_r8_3d_cube) + end if + end if + + if (do_io) then + deallocate(dimids_2d) + deallocate(dimids_3d) + end if deallocate(fcstField) deallocate(varids) deallocate(compress_err) - if (mype==0) then - ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) + if (do_io) 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 @@ -406,21 +739,19 @@ subroutine get_global_attr(fldbundle, ncid, rc) integer, intent(out) :: rc ! local variable - integer :: i, attcount + 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, & + attnestflag=ESMF_ATTNEST_OFF, count=attCount, & rc=rc); ESMF_ERR_RETURN(rc) do i=1,attCount @@ -458,7 +789,7 @@ subroutine get_global_attr(fldbundle, ncid, rc) end do end subroutine get_global_attr -! + !---------------------------------------------------------------------------------------- subroutine get_grid_attr(grid, prefix, ncid, varid, rc) type(ESMF_Grid), intent(in) :: grid @@ -468,7 +799,7 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) integer, intent(out) :: rc ! local variable - integer :: i, attcount, n, ind + integer :: i, attCount, n, ind integer :: ncerr character(len=ESMF_MAXSTR) :: attName type(ESMF_TypeKind_Flag) :: typekind @@ -479,16 +810,14 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) character(len=ESMF_MAXSTR) :: varcval ! call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & - attnestflag=ESMF_ATTNEST_OFF, Count=attcount, & + 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)//":") @@ -507,10 +836,10 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) 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 + ! FIXME: _FillValue must be cast to var type when using + ! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue. ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr8val); NC_ERR_STOP(ncerr) - endif + end if else if (typekind==ESMF_TYPEKIND_CHARACTER) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & @@ -525,6 +854,7 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) end subroutine get_grid_attr +!---------------------------------------------------------------------------------------- subroutine add_dim(ncid, dim_name, dimid, grid, rc) integer, intent(in) :: ncid character(len=*), intent(in) :: dim_name @@ -533,75 +863,127 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) integer, intent(out) :: rc ! local variable - integer :: i, attcount, n, dim_varid + integer :: 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) + if (trim(dim_name) == "time") then + ! using an unlimited dim requires collective mode (NF90_COLLECTIVE) + ! for parallel writes, which seems to slow things down on hera. + !ncerr = nf90_def_dim(ncid, trim(dim_name), NF90_UNLIMITED, dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_dim(ncid, trim(dim_name), 1, 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) + 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_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) + 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_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 + if (par) then + ncerr = nf90_var_par_access(ncid, dim_varid, NF90_INDEPENDENT); NC_ERR_STOP(ncerr) + 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" + subroutine quantize_array_3d(array, dataMin, dataMax, nbits, compress_err) + + real(4), dimension(:,:,:), intent(inout) :: array + real(4), intent(in) :: dataMin, dataMax + integer, intent(in) :: nbits + real(4), intent(out) :: compress_err + + real(4) :: scale_fact, offset + real(4), dimension(:,:,:), allocatable :: array_save + ! 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). + scale_fact = (dataMax - dataMin) / (2**nbits-1) + offset = dataMin + if (scale_fact > 0.) then + allocate(array_save, source=array) + array = scale_fact*(nint((array_save - offset) / scale_fact)) + offset + ! compute max abs compression error + compress_err = maxval(abs(array_save-array)) + deallocate(array_save) + else + ! field is constant + compress_err = 0. 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 subroutine quantize_array_3d + + subroutine quantize_array_4d(array, dataMin, dataMax, nbits, compress_err) + + real(4), dimension(:,:,:,:), intent(inout) :: array + real(4), intent(in) :: dataMin, dataMax + integer, intent(in) :: nbits + real(4), intent(out) :: compress_err + + real(4) :: scale_fact, offset + real(4), dimension(:,:,:,:), allocatable :: array_save + + ! 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). + scale_fact = (dataMax - dataMin) / (2**nbits-1) + offset = dataMin + if (scale_fact > 0.) then + allocate(array_save, source=array) + array = scale_fact*(nint((array_save - offset) / scale_fact)) + offset + ! compute max abs compression error + compress_err = maxval(abs(array_save-array)) + deallocate(array_save) + else + ! field is constant + compress_err = 0. + end if + end subroutine quantize_array_4d +!---------------------------------------------------------------------------------------- end module module_write_netcdf diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 deleted file mode 100644 index 0506d794a..000000000 --- a/io/module_write_netcdf_parallel.F90 +++ /dev/null @@ -1,627 +0,0 @@ -#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_parallel - - use esmf - use netcdf - use module_fv3_io_def,only : ideflate, nbits, & - output_grid,dx,dy,lon1,lat1,lon2,lat2 - use mpi - - implicit none - private - public write_netcdf_parallel - - contains - -#ifdef NO_PARALLEL_NETCDF -!---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(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, ichunk2d, jchunk2d, & - ichunk3d, jchunk3d, kchunk3d - integer, optional,intent(out) :: rc - print *,'in stub write_netcdf_parallel - model not built with parallel netcdf support, return' - end subroutine write_netcdf_parallel -#else -!---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(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, ichunk2d, jchunk2d, & - ichunk3d, jchunk3d, kchunk3d - integer, optional,intent(out) :: rc -! -!** local vars - integer :: i,j,m,n,k,istart,iend,jstart,jend,i1,i2,j1,j2,k1,k2 - integer :: lm - - integer, dimension(:), allocatable :: fldlev - real(ESMF_KIND_R4), dimension(:,:), pointer :: arrayr4 - real(ESMF_KIND_R8), dimension(:,:), pointer :: arrayr8 - real(ESMF_KIND_R4), dimension(:,:,:), pointer :: arrayr4_3d,arrayr4_3d_save - real(ESMF_KIND_R8), dimension(:,:,:), pointer :: 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 :: totalLBound2d(2),totalUBound2d(2),totalLBound3d(3),totalUBound3d(3) - - 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,ierr - 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(:)) - -! create netcdf file for parallel access - - ncerr = nf90_create(trim(filename),& - cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),& - comm=mpi_comm, info = MPI_INFO_NULL, ncid=ncid); NC_ERR_STOP(ncerr) -! disable auto filling. - 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_var_par_access(ncid, im_varid, NF90_INDEPENDENT) - ncerr = nf90_def_var(ncid, "lon", NF90_DOUBLE, (/im_dimid,jm_dimid/), lon_varid); NC_ERR_STOP(ncerr) - !ncerr = nf90_var_par_access(ncid, lon_varid, NF90_INDEPENDENT) - 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_var_par_access(ncid, jm_varid, NF90_INDEPENDENT) - ncerr = nf90_def_var(ncid, "lat", NF90_DOUBLE, (/im_dimid,jm_dimid/), lat_varid); NC_ERR_STOP(ncerr) - ncerr = nf90_var_par_access(ncid, lat_varid, NF90_INDEPENDENT) - 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/)); NC_ERR_STOP(ncerr) - endif - ! compression filters require collective access. - ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) - else - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) - ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) - 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) - ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) - 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 - ! compression filters require collective access. - ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) - else - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) - ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) - 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) - ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) - 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 when using NF90_NETCDF4 - 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) == 'global_latlon' .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 of define mode - - ! write grid_xt, grid_yt values - call ESMF_GridGetCoord(wrtGrid, coordDim=1, farrayPtr=arrayr8, rc=rc); ESMF_ERR_RETURN(rc) - istart = lbound(arrayr8,1); iend = ubound(arrayr8,1) - jstart = lbound(arrayr8,2); jend = ubound(arrayr8,2) - !print *,'in write netcdf mpi dim 1',istart,iend,jstart,jend,shape(arrayr8),minval(arrayr8(:,jstart)),maxval(arrayr8(:,jstart)) - - if (trim(output_grid) == 'gaussian_grid' .or. & - trim(output_grid) == 'global_latlon' .or. & - trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,jstart),start=(/istart/), count=(/iend-istart+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, start=(/istart,jstart/)); NC_ERR_STOP(ncerr) - - call ESMF_GridGetCoord(wrtGrid, coordDim=2, farrayPtr=arrayr8, rc=rc); ESMF_ERR_RETURN(rc) - !print *,'in write netcdf mpi dim 2',istart,iend,jstart,jend,shape(arrayr8),minval(arrayr8(istart,:)),maxval(arrayr8(istart,:)) - if (trim(output_grid) == 'gaussian_grid' .or. & - trim(output_grid) == 'global_latlon' .or. & - trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(istart,:),start=(/jstart/),count=(/jend-jstart+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, start=(/istart,jstart/)); NC_ERR_STOP(ncerr) - - 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_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr4, totalLBound=totalLBound2d, totalUBound=totalUBound2d,rc=rc); ESMF_ERR_RETURN(rc) - !print *,'field name=',trim(fldName),'bound=',totalLBound2d,'ubound=',totalUBound2d - ncerr = nf90_put_var(ncid, varids(i), values=arrayr4, start=(/totalLBound2d(1),totalLBound2d(2),1/)); NC_ERR_STOP(ncerr) - else if (typekind == ESMF_TYPEKIND_R8) then - call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr8, totalLBound=totalLBound2d, totalUBound=totalUBound2d,rc=rc); ESMF_ERR_RETURN(rc) - ncerr = nf90_put_var(ncid, varids(i), values=arrayr8, start=(/totalLBound2d(1),totalLBound2d(2),1/)); NC_ERR_STOP(ncerr) - end if - else if (fldlev(i) > 1) then - if (typekind == ESMF_TYPEKIND_R4) then - call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr4_3d, totalLBound=totalLBound3d, totalUBound=totalUBound3d,rc=rc); ESMF_ERR_RETURN(rc) - if (ideflate > 0 .and. nbits > 0) then - i1=totalLBound3d(1);i2=totalUBound3d(1) - j1=totalLBound3d(2);j2=totalUBound3d(2) - k1=totalLBound3d(3);k2=totalUBound3d(3) - dataMax = maxval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) - dataMin = minval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) - call mpi_allreduce(mpi_in_place,dataMax,1,mpi_real4,mpi_max,mpi_comm,ierr) - call mpi_allreduce(mpi_in_place,dataMin,1,mpi_real4,mpi_min,mpi_comm,ierr) - ! 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) - scale_fact = (dataMax - dataMin) / (2**nbits-1); offset = dataMin - if (scale_fact > 0.) then - allocate(arrayr4_3d_save(i1:i2,j1:j2,k1:k2)) - arrayr4_3d_save(i1:i2,j1:j2,k1:k2)=arrayr4_3d(i1:i2,j1:j2,k1:k2) - arrayr4_3d = scale_fact*(nint((arrayr4_3d_save - offset) / scale_fact)) + offset - ! compute max abs compression error. - compress_err(i) = & - maxval(abs(arrayr4_3d_save(i1:i2,j1:j2,k1:k2)-arrayr4_3d(i1:i2,j1:j2,k1:k2))) - deallocate(arrayr4_3d_save) - call mpi_allreduce(mpi_in_place,compress_err(i),1,mpi_real4,mpi_max,mpi_comm,ierr) - !print *,'field name=',trim(fldName),dataMin,dataMax,compress_err(i) - else - ! field is constant - compress_err(i) = 0. - endif - endif - ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3),1/)); NC_ERR_STOP(ncerr) - else if (typekind == ESMF_TYPEKIND_R8) then - call ESMF_FieldGet(fcstField(i), localDe=0, farrayPtr=arrayr8_3d, totalLBound=totalLBound3d, totalUBound=totalUBound3d,rc=rc); ESMF_ERR_RETURN(rc) - !print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d - ncerr = nf90_put_var(ncid, varids(i), values=arrayr8_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3),1/)); NC_ERR_STOP(ncerr) - end if - - end if !end fldlev(i) - - end do ! end fieldCount - - if (ideflate > 0 .and. nbits > 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(fcstField) - deallocate(varids) - deallocate(compress_err) - - ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) - !call mpi_barrier(mpi_comm,ierr) - !print *,'netcdf parallel close, finished write_netcdf_parallel' - - end subroutine write_netcdf_parallel -#endif - -!---------------------------------------------------------------------------------------- - 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 when using - ! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue. - 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 - ! using an unlimited dim requires collective mode (NF90_COLLECTIVE) - ! for parallel writes, which seems to slow things down on hera. - !ncerr = nf90_def_dim(ncid, trim(dim_name), NF90_UNLIMITED, dimid); NC_ERR_STOP(ncerr) - ncerr = nf90_def_dim(ncid, trim(dim_name), 1, 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) - ncerr = nf90_var_par_access(ncid, dim_varid, NF90_INDEPENDENT) - 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) - ncerr = nf90_var_par_access(ncid, dim_varid, NF90_INDEPENDENT) - 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 - -end module module_write_netcdf_parallel diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 70257b8d6..97c9c8885 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -43,7 +43,6 @@ module module_wrt_grid_comp use module_write_netcdf, only : write_netcdf use physcons, only : pi => con_pi use inline_post, only : inline_post_run, inline_post_getattr - use module_write_netcdf_parallel, only : write_netcdf_parallel ! !----------------------------------------------------------------------- ! @@ -278,6 +277,9 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) if (ideflate < 0) ideflate=0 call ESMF_ConfigGetAttribute(config=CF,value=nbits,default=0,label ='nbits:',rc=rc) + if (lprnt) then + print *,'ideflate=',ideflate,' nbits=',nbits + end if ! nbits quantization level for lossy compression (must be between 1 and 31) ! 1 is most compression, 31 is least. If outside this range, set to zero ! which means use lossless compression. @@ -1576,6 +1578,11 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) if (trim(output_grid) == 'cubed_sphere_grid') then wbeg = MPI_Wtime() + if (trim(output_file(nbdl)) == 'netcdf_parallel') then + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .true., wrt_mpi_comm,wrt_int_state%mype, & + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) + else call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & convention="NetCDF", purpose="FV3", & status=ESMF_FILESTATUS_REPLACE, & @@ -1591,6 +1598,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + end if wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' actual netcdf Write Time is ',wend-wbeg & @@ -1602,8 +1610,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) if (trim(output_file(nbdl)) == 'netcdf') then wbeg = MPI_Wtime() - call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .false., wrt_mpi_comm,wrt_int_state%mype, & ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then @@ -1620,9 +1628,9 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) line=__LINE__, file=__FILE__)) return #endif wbeg = MPI_Wtime() - call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & - trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& - ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .true., wrt_mpi_comm,wrt_int_state%mype, & + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & @@ -1658,8 +1666,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) if (trim(output_file(nbdl)) == 'netcdf') then wbeg = MPI_Wtime() - call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .false., wrt_mpi_comm,wrt_int_state%mype, & ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then @@ -1676,9 +1684,9 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) line=__LINE__, file=__FILE__)) return #endif wbeg = MPI_Wtime() - call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & - trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& - ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .true., wrt_mpi_comm,wrt_int_state%mype, & + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & @@ -1714,8 +1722,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) if (trim(output_file(nbdl)) == 'netcdf' .and. nbits==0) then wbeg = MPI_Wtime() - call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .false., wrt_mpi_comm,wrt_int_state%mype, & ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (mype == lead_write_task) then @@ -1732,9 +1740,9 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) line=__LINE__, file=__FILE__)) return #endif wbeg = MPI_Wtime() - call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & - trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& - ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) + call write_netcdf(wrt_int_state%wrtFB(nbdl),trim(filename), & + .true., wrt_mpi_comm,wrt_int_state%mype, & + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg &