diff --git a/fv3_cap.F90 b/fv3_cap.F90 index a099ae32e..c211a9044 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -37,7 +37,8 @@ module fv3gfs_cap_mod wrttasks_per_group, n_group, & lead_wrttask, last_wrttask, & output_grid, output_file, & - imo, jmo, write_nemsioflip, & + imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,& + ichunk3d,jchunk3d,kchunk3d, & write_fsyncflag, nsout_io, & cen_lon, cen_lat, ideflate, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -235,8 +236,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer,dimension(6) :: date, date_init integer :: mpi_comm_atm - integer :: i, j, k, io_unit, urc + integer :: i, j, k, io_unit, urc, ierr integer :: petcount, mype + integer :: num_output_file logical :: isPetLocal logical :: OPENED character(ESMF_MAXSTR) :: name @@ -307,6 +309,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call ESMF_ConfigGetAttribute(config=CF,value=iau_offset,default=0,label ='iau_offset:',rc=rc) if (iau_offset < 0) iau_offset=0 + ! chunksizes for netcdf_parallel + call ESMF_ConfigGetAttribute(config=CF,value=ichunk2d,default=0,label ='ichunk2d:',rc=rc) + call ESMF_ConfigGetAttribute(config=CF,value=jchunk2d,default=0,label ='jchunk2d:',rc=rc) + call ESMF_ConfigGetAttribute(config=CF,value=ichunk3d,default=0,label ='ichunk3d:',rc=rc) + call ESMF_ConfigGetAttribute(config=CF,value=jchunk3d,default=0,label ='jchunk3d:',rc=rc) + call ESMF_ConfigGetAttribute(config=CF,value=kchunk3d,default=0,label ='kchunk3d:',rc=rc) + + ! zlib compression flag call ESMF_ConfigGetAttribute(config=CF,value=ideflate,default=0,label ='ideflate:',rc=rc) if (ideflate < 0) ideflate=0 @@ -346,8 +356,33 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) CALL ESMF_ConfigGetAttribute(config=CF,value=filename_base(i), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return enddo - if(mype == 0) print *,'af nems config,num_files=',num_files, & - 'filename_base=',filename_base + + allocate(output_file(num_files)) + num_output_file = ESMF_ConfigGetLen(config=CF, label ='output_file:',rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (num_files == num_output_file) then + CALL ESMF_ConfigGetAttribute(CF,valueList=output_file,label='output_file:', & + count=num_files, rc=RC) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + do i = 1, num_files + if(output_file(i) /= "netcdf" .and. output_file(i) /= "netcdf_parallel") then + write(0,*)"fv3_cap.F90: only netcdf and netcdf_parallel are allowed for multiple values of output_file" + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + enddo + else if ( num_output_file == 1) then + CALL ESMF_ConfigGetAttribute(CF,valuelist=output_file,label='output_file:', count=1, rc=RC) + output_file(1:num_files) = output_file(1) + else + output_file(1:num_files) = 'netcdf' + endif + if(mype == 0) then + print *,'af nems config,num_files=',num_files + do i=1,num_files + print *,'num_file=',i,'filename_base= ',trim(filename_base(i)),& + ' output_file= ',trim(output_file(i)) + enddo + endif ! ! variables for alarms call ESMF_ConfigGetAttribute(config=CF, value=nfhout, label ='nfhout:', rc=rc) @@ -359,10 +394,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! variables for I/O options call ESMF_ConfigGetAttribute(config=CF, value=output_grid, label ='output_grid:',rc=rc) - call ESMF_ConfigGetAttribute(config=CF, value=output_file, label ='output_file:',rc=rc) if (mype == 0) then print *,'output_grid=',trim(output_grid) - print *,'output_file=',trim(output_file) end if write_nemsioflip =.false. write_fsyncflag =.false. diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index 5eecdd9ea..e2cedc43d 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -5,6 +5,10 @@ else() add_definitions(-DNO_INLINE_POST) endif() +if(NOT PARALLEL_NETCDF) + add_definitions(-DNO_PARALLEL_NETCDF) +endif() + add_library( io @@ -12,6 +16,7 @@ add_library( FV3GFS_io.F90 module_write_nemsio.F90 module_write_netcdf.F90 + module_write_netcdf_parallel.F90 module_fv3_io_def.F90 module_write_internal_state.F90 module_wrt_grid_comp.F90 diff --git a/io/makefile b/io/makefile index 82e61f68f..b6a3946b0 100644 --- a/io/makefile +++ b/io/makefile @@ -41,6 +41,7 @@ SRCS_F90 = \ $(POST_SRC) \ ./module_write_nemsio.F90 \ ./module_write_netcdf.F90 \ + ./module_write_netcdf_parallel.F90 \ ./module_fv3_io_def.F90 \ ./module_write_internal_state.F90 \ ./module_wrt_grid_comp.F90 @@ -69,7 +70,9 @@ post_nems_routines.o: post_nems_routines.F90 module_write_nemsio.o: module_write_nemsio.F90 $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) $(NEMSIOINC) -c module_write_nemsio.F90 module_write_netcdf.o: module_write_netcdf.F90 - $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) $(NEMSIOINC) -c module_write_netcdf.F90 + $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) -c module_write_netcdf.F90 +module_write_netcdf_parallel.o: module_write_netcdf_parallel.F90 + $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) -c module_write_netcdf_parallel.F90 module_write_internal_state.o: module_write_internal_state.F90 $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) -c module_write_internal_state.F90 module_wrt_grid_comp.o: module_wrt_grid_comp.F90 diff --git a/io/module_fv3_io_def.F90 b/io/module_fv3_io_def.F90 index 17e9f308d..b31ab666b 100644 --- a/io/module_fv3_io_def.F90 +++ b/io/module_fv3_io_def.F90 @@ -17,13 +17,14 @@ module module_fv3_io_def integer :: num_files character(255) :: app_domain character(255) :: output_grid - character(255) :: output_file integer :: imo,jmo + integer :: ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d integer :: nbdlphys integer :: nsout_io, iau_offset, ideflate, nbits real :: cen_lon, cen_lat, lon1, lat1, lon2, lat2, dlon, dlat real :: stdlat1, stdlat2, dx, dy character(255),dimension(:),allocatable :: filename_base + character(255),dimension(:),allocatable :: output_file ! integer,dimension(:),allocatable :: lead_wrttask, last_wrttask ! diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index fc59c75c9..270e36ec6 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -1,7 +1,3 @@ -!#define ESMF_ERR_ABORT(rc) \ -!if (rc /= ESMF_SUCCESS) write(0,*) 'rc=',rc,__FILE__,__LINE__; \ -! if (ESMF_LogFoundError(rc, msg=ESMF_LOGERR_PASSTHRU, 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) \ @@ -22,7 +18,7 @@ module module_write_netcdf contains !---------------------------------------------------------------------------------------- - subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) + 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 @@ -30,6 +26,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc 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 @@ -69,9 +66,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc 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) @@ -117,15 +112,10 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! create netcdf file and enter define mode if (mype==0) then - if (ideflate == 0) then - ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),NF90_SHARE), & - ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) - else - ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), & - ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) - endif + 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) @@ -158,32 +148,53 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then if (ideflate > 0) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,time_dimid/), varids(i), & - shuffle=.true.,deflate_level=ideflate, & - chunksizes=(/im,jm,1/)); NC_ERR_STOP(ncerr) + 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) + 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) + 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 + write(0,*)'Unsupported typekind ', typekind + stop end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then if (ideflate > 0) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & - shuffle=.false.,deflate_level=ideflate, & - chunksizes=(/im,jm,1,1/)); NC_ERR_STOP(ncerr) + ! 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 + (/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) @@ -224,7 +235,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, & rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' ) then + 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 @@ -241,23 +252,23 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end do ! i=1,fieldCount -! write grid_xt, grid_yt attributes + ! 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) + 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) + 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) + 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) @@ -377,6 +388,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc deallocate(fcstField) deallocate(varids) + deallocate(compress_err) if (mype==0) then ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 new file mode 100644 index 000000000..7d7dbb0f7 --- /dev/null +++ b/io/module_write_netcdf_parallel.F90 @@ -0,0 +1,617 @@ +#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) == '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) == '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) == '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 + 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) + 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 84769eaea..a7177fe22 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -34,7 +34,8 @@ module module_wrt_grid_comp use module_fv3_io_def, only : num_pes_fcst,lead_wrttask, last_wrttask, & n_group, num_files, app_domain, & filename_base, output_grid, output_file, & - imo, jmo, write_nemsioflip, & + imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,& + ichunk3d,jchunk3d,kchunk3d, & nsout => nsout_io, & cen_lon, cen_lat, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -43,6 +44,7 @@ module module_wrt_grid_comp use module_write_netcdf, only : write_netcdf use physcons, only : pi => con_pi use post_gfs, only : post_run_gfs, post_getattr_gfs + use module_write_netcdf_parallel, only : write_netcdf_parallel ! !----------------------------------------------------------------------- ! @@ -1122,14 +1124,14 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) !----------------------------------------------------------------------- ! call ESMF_LogWrite("before initialize for nemsio file", ESMF_LOGMSG_INFO, rc=rc) - if (trim(output_grid) == 'gaussian_grid' .and. trim(output_file) == 'nemsio') then -! if (lprnt) write(0,*) 'in wrt initial, befnemsio_first_call wrt_int_state%FBcount=',wrt_int_state%FBcount - do i= 1, wrt_int_state%FBcount - call nemsio_first_call(wrt_int_state%wrtFB(i), imo, jmo, & - wrt_int_state%mype, ntasks, wrt_mpi_comm, & - wrt_int_state%FBcount, i, idate, lat, lon, rc) - enddo - endif + do i= 1, wrt_int_state%FBcount + if (trim(output_grid) == 'gaussian_grid' .and. trim(output_file(i)) == 'nemsio') then +! if (lprnt) write(0,*) 'in wrt initial, befnemsio_first_call wrt_int_state%FBcount=',wrt_int_state%FBcount + call nemsio_first_call(wrt_int_state%wrtFB(i), imo, jmo, & + wrt_int_state%mype, ntasks, wrt_mpi_comm, & + wrt_int_state%FBcount, i, idate, lat, lon, rc) + endif + enddo call ESMF_LogWrite("after initialize for nemsio file", ESMF_LOGMSG_INFO, rc=rc) ! !----------------------------------------------------------------------- @@ -1170,7 +1172,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) type(ESMF_Time) :: currtime type(ESMF_TypeKind_Flag) :: datatype type(ESMF_Field) :: field_work - type(ESMF_Grid) :: grid_work, fbgrid + type(ESMF_Grid) :: grid_work, fbgrid, wrtgrid type(ESMF_Array) :: array_work type(ESMF_State),save :: stateGridFB type(optimizeT), save :: optimize(4) @@ -1363,7 +1365,47 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) file_bundle = wrt_int_state%wrtFB(nbdl) endif - if ( trim(output_file) == 'nemsio' ) then + ! set default chunksizes for netcdf output + ! (use MPI decomposition size). + ! if chunksize parameter set to negative value, + ! netcdf library default is used. + if (output_file(nbdl)(1:6) == 'netcdf') then + if (ichunk2d == 0) then + if( wrt_int_state%mype == 0 ) & + ichunk2d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 + call mpi_bcast(ichunk2d,1,mpi_integer,0,wrt_mpi_comm,rc) + endif + if (jchunk2d == 0) then + if( wrt_int_state%mype == 0 ) & + jchunk2d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 + call mpi_bcast(jchunk2d,1,mpi_integer,0,wrt_mpi_comm,rc) + endif + if (ichunk3d == 0) then + if( wrt_int_state%mype == 0 ) & + ichunk3d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 + call mpi_bcast(ichunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) + endif + if (jchunk3d == 0) then + if( wrt_int_state%mype == 0 ) & + jchunk3d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 + call mpi_bcast(jchunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) + endif + if (kchunk3d == 0 .and. nbdl == 1) then + if( wrt_int_state%mype == 0 ) then + call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), grid=wrtgrid) + call ESMF_AttributeGet(wrtgrid, convention="NetCDF", purpose="FV3", & + attnestflag=ESMF_ATTNEST_OFF, name='pfull', & + itemCount=kchunk3d, rc=rc) + endif + call mpi_bcast(kchunk3d,1,mpi_integer,0,wrt_mpi_comm,rc) + endif + if (wrt_int_state%mype == 0) then + print *,'ichunk2d,jchunk2d',ichunk2d,jchunk2d + print *,'ichunk3d,jchunk3d,kchunk3d',ichunk3d,jchunk3d,kchunk3d + endif + endif + + if ( trim(output_file(nbdl)) == 'nemsio' ) then filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nemsio' else filename = trim(wrt_int_state%wrtFB_names(nbdl))//'f'//trim(cfhour)//'.nc' @@ -1413,7 +1455,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) else if (trim(output_grid) == 'gaussian_grid') then - if (trim(output_file) == 'nemsio') then + if (trim(output_file(nbdl)) == 'nemsio') then wbeg = MPI_Wtime() call write_nemsio(file_bundle,trim(filename),nf_hours, nf_minutes, & @@ -1424,18 +1466,37 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif - else if (trim(output_file) == 'netcdf') then + else 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,rc) + wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg & ,' at Fcst ',NF_HOURS,':',NF_MINUTES endif - else if (trim(output_file) == 'netcdf_esmf') then + else if (trim(output_file(nbdl)) == 'netcdf_parallel') then + +#ifdef NO_PARALLEL_NETCDF + rc = ESMF_RC_NOT_IMPL + print *,'netcdf_parallel not available on this machine' + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + 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) + wend = MPI_Wtime() + if (lprnt) then + write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & + ,' at Fcst ',NF_HOURS,':',NF_MINUTES + endif + + else if (trim(output_file(nbdl)) == 'netcdf_esmf') then wbeg = MPI_Wtime() call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & @@ -1482,11 +1543,12 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) write(*,'(A,F10.5,A,I4.2,A,I2.2)')' mask_fields time is ',wend-wbeg endif - if (trim(output_file) == 'netcdf') then + 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,rc) + wrt_mpi_comm,wrt_int_state%mype,imo,jmo,& + ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d,rc) wend = MPI_Wtime() if (mype == lead_write_task) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' netcdf Write Time is ',wend-wbeg &