From 45d0ce2771a022ebbbfcb08005f8f8c67950a6c4 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Wed, 8 Jan 2020 02:38:52 +0000 Subject: [PATCH 01/32] netcdf parallel writing and lon/lat in netcdf file --- io/makefile | 5 +- io/module_write_nemsio.F90 | 2 +- io/module_write_netcdf.F90 | 11 +- io/module_write_netcdf_parallel.F90 | 610 ++++++++++++++++++++++++++++ io/module_wrt_grid_comp.F90 | 71 +++- io/post_gfs.F90 | 2 +- 6 files changed, 674 insertions(+), 27 deletions(-) create mode 100644 io/module_write_netcdf_parallel.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_write_nemsio.F90 b/io/module_write_nemsio.F90 index 3afd66789..e51f64a52 100644 --- a/io/module_write_nemsio.F90 +++ b/io/module_write_nemsio.F90 @@ -51,7 +51,7 @@ subroutine nemsio_first_call(fieldbundle, imo, jmo, & integer, intent(in) :: wrt_mype, wrt_ntasks, wrt_mpi_comm integer, intent(in) :: wrt_nbdl, mybdl integer, intent(in) :: inidate(7) - real, intent(in) :: lat(:), lon(:) + real(8), intent(in) :: lat(:), lon(:) integer, optional,intent(out) :: rc !** local vars diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 1fce3d8b9..52cf740e6 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -42,7 +42,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save real(8), dimension(:,:,:), allocatable :: arrayr8_3d - real(8) rad2dg,x(im),y(jm) + real(8) x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound @@ -71,7 +71,6 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! !! ! - rad2dg = 45./atan(1.0) call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) @@ -247,7 +246,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(:,1) ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) @@ -271,7 +270,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc) @@ -279,7 +278,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then - ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:) ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) @@ -303,7 +302,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif do i=1, fieldCount diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 new file mode 100644 index 000000000..c2ea784fd --- /dev/null +++ b/io/module_write_netcdf_parallel.F90 @@ -0,0 +1,610 @@ +!#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) \ + 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 + +!---------------------------------------------------------------------------------------- + subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) +! + type(ESMF_FieldBundle), intent(in) :: fieldbundle + type(ESMF_FieldBundle), intent(in) :: wrtfb + character(*), intent(in) :: filename + integer, intent(in) :: mpi_comm + integer, intent(in) :: mype + integer, intent(in) :: im, jm + integer, optional,intent(out) :: rc +! +!** local vars + integer :: i,j,m,n,k,istart,iend,jstart,jend + 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) rad2dg,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 :: start1d(1), count1d(1), start(2), count(2), start3d(3), count3d(3) + integer :: totalLBound2d(2), totalUBound2d(2),totalLBound3d(3),totalUBound3d(3) + + integer :: varival + real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax + real(8) :: varr8val + character(len=ESMF_MAXSTR) :: varcval + + character(128) :: time_units + + integer :: ncerr + integer :: ncid + integer :: oldMode + integer :: im_dimid, jm_dimid, pfull_dimid, phalf_dimid, time_dimid + integer :: im_varid, jm_varid, lm_varid, time_varid, lon_varid, lat_varid + integer, dimension(:), allocatable :: varids +! +!! +! + rad2dg = 45./atan(1.0) + + call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) + + allocate(fldlev(fieldCount)) ; fldlev = 0 + allocate(fcstField(fieldCount)) + allocate(varids(fieldCount)) + + call ESMF_FieldBundleGet(fieldbundle, fieldList=fcstField, grid=wrtGrid, & +! itemorderflag=ESMF_ITEMORDER_ADDORDER, & + rc=rc); ESMF_ERR_RETURN(rc) + + call ESMF_GridGet(wrtgrid, dimCount=gridDimCount, rc=rc); ESMF_ERR_RETURN(rc) + + do i=1,fieldCount + call ESMF_FieldGet(fcstField(i), dimCount=fieldDimCount, rc=rc); ESMF_ERR_RETURN(rc) + if (fieldDimCount > 3) then + write(0,*)"write_netcdf: Only 2D and 3D fields are supported!" + stop + end if + if (fieldDimCount > gridDimCount) then + allocate(ungriddedLBound(fieldDimCount-gridDimCount)) + allocate(ungriddedUBound(fieldDimCount-gridDimCount)) + call ESMF_FieldGet(fcstField(i), & + ungriddedLBound=ungriddedLBound, & + ungriddedUBound=ungriddedUBound, rc=rc); ESMF_ERR_RETURN(rc) + fldlev(i) = ungriddedUBound(fieldDimCount-gridDimCount) - & + ungriddedLBound(fieldDimCount-gridDimCount) + 1 + deallocate(ungriddedLBound) + deallocate(ungriddedUBound) + else if (fieldDimCount == 2) then + fldlev(i) = 1 + end if + end do + + lm = maxval(fldlev(:)) + +! allocate(arrayr4(im,jm)) +! allocate(arrayr8(im,jm)) +! allocate(arrayr4_3d(im,jm,lm),arrayr4_3d_save(im,jm,lm)) +! allocate(arrayr8_3d(im,jm,lm)) + +! create netcdf file and enter define mode +!not build with pnetcdf +! if (ideflate == 0) then +! ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), & +! comm=mpi_comm, info = MPI_INFO_NULL,ncid=ncid); NC_ERR_STOP(ncerr) +! ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) +! else +! ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), & +! ncid=ncid); NC_ERR_STOP(ncerr) +! ! if compression on use HDF5 format with default _FillValue +! endif +! + ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_NETCDF4), & + comm=mpi_comm, info = MPI_INFO_NULL, ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) + + ! define dimensions + ncerr = nf90_def_dim(ncid, "grid_xt", im, im_dimid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_dim(ncid, "grid_yt", jm, jm_dimid); NC_ERR_STOP(ncerr) + ! define coordinate variables + ncerr = nf90_def_var(ncid, "grid_xt", NF90_DOUBLE, im_dimid, im_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "lon", NF90_DOUBLE, (/im_dimid,jm_dimid/), lon_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, lon_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, lon_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "cartesian_axis", "X"); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "grid_yt", NF90_DOUBLE, jm_dimid, jm_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, "lat", NF90_DOUBLE, (/im_dimid,jm_dimid/), lat_varid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, lat_varid, "long_name", "T-cell latitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, lat_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "cartesian_axis", "Y"); NC_ERR_STOP(ncerr) + + if (lm > 1) then + call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, rc) + call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, rc) + end if + + call add_dim(ncid, "time", time_dimid, wrtgrid, rc) + + call get_global_attr(wrtfb, ncid, rc) + + do i=1, fieldCount + call ESMF_FieldGet(fcstField(i), name=fldName, typekind=typekind, rc=rc); ESMF_ERR_RETURN(rc) + + ! define variables + if (fldlev(i) == 1) then + if (typekind == ESMF_TYPEKIND_R4) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) + if (ideflate > 0) then + ! shuffle filter on for lossless compression + ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate) + NC_ERR_STOP(ncerr) + 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 + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) + if (ideflate > 0) then + ! shuffle filter off since lossy compression used + ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate) + NC_ERR_STOP(ncerr) + endif + else if (typekind == ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) + else + write(0,*)'Unsupported typekind ', typekind + stop + end if + end if + + ! define variable attributes + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + attnestflag=ESMF_ATTNEST_OFF, Count=attcount, & + rc=rc); ESMF_ERR_RETURN(rc) + + do j=1,attCount + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + attnestflag=ESMF_ATTNEST_OFF, attributeIndex=j, & + name=attName, typekind=attTypeKind, itemCount=n, & + rc=rc); ESMF_ERR_RETURN(rc) + + if ( index(trim(attName),"ESMF") /= 0 ) then + cycle + endif + + if (attTypeKind==ESMF_TYPEKIND_I4) then + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + name=trim(attName), value=varival, & + rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_att(ncid, varids(i), trim(attName), varival); NC_ERR_STOP(ncerr) + + else if (attTypeKind==ESMF_TYPEKIND_R4) then + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + name=trim(attName), value=varr4val, & + rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr4val); NC_ERR_STOP(ncerr) + + else if (attTypeKind==ESMF_TYPEKIND_R8) then + call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & + name=trim(attName), value=varr8val, & + rc=rc); ESMF_ERR_RETURN(rc) + if (trim(attName) /= '_FillValue' .or. ideflate == 0) 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 + + 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) + +! call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) +! if (mype==0) then + if (trim(output_grid) == 'gaussian_grid' .or. trim(output_grid) == 'regional_latlon') then + istart = lbound(arrayr8,1) + iend = ubound(arrayr8,1) + jstart = lbound(arrayr8,2) + jend = ubound(arrayr8,2) + print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend + start1d(1) = istart + count1d(1) = iend - istart + 1 + ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(istart:iend,1),start=start1d, count=count1d ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'rotated_latlon') then + do i=1,im + x(i) = lon1 + (lon2-lon1)/(im-1) * (i-1) + enddo + ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'lambert_conformal') then + do i=1,im + x(i) = dx * (i-1) + enddo + ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + endif + start = (/istart,jstart/) + count = (/iend-istart+1,jend-jstart+1/) + ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8, start=start, count=count ); NC_ERR_STOP(ncerr) + + call ESMF_GridGetCoord(wrtGrid, coordDim=2, farrayPtr=arrayr8, rc=rc); ESMF_ERR_RETURN(rc) +! call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) + if (trim(output_grid) == 'gaussian_grid' .or. trim(output_grid) == 'regional_latlon') then + istart = lbound(arrayr8,1) + iend = ubound(arrayr8,1) + jstart = lbound(arrayr8,2) + jend = ubound(arrayr8,2) + print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend + start1d(1) = jstart + count1d(1) = jend - jstart + 1 + ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:),start=start1d,count=count1d ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'rotated_latlon') then + do j=1,jm + y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1) + enddo + ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'lambert_conformal') then + do j=1,jm + y(j) = dy * (j-1) + enddo + ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + endif + start = (/istart,jstart/) + count = (/iend-istart+1, jend-jstart+1/) + ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8, start=start, count=count ); 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 + count3d=(/totalUBound2d(1)-totalLBound2d(1)+1, & + totalUBound2d(2)-totalLBound2d(2)+1, 1/) + ncerr = nf90_put_var(ncid, varids(i), values=arrayr4, start=(/totalLBound2d(1),totalLBound2d(2),1/),count=count3d); 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) + count3d=(/totalUBound2d(1)-totalLBound2d(1)+1, & + totalUBound2d(2)-totalLBound2d(2)+1, 1/) + ncerr = nf90_put_var(ncid, varids(i), values=arrayr8, start=(/totalLBound2d(1),totalLBound2d(2),1/),count=count3d); 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) + print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d + 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 = maxval(abs(arrayr4_3d_save-arrayr4_3d)) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + endif + count3d = (/totalUBound3d(1)-totalLBound3d(1)+1, & + totalUBound3d(2)-totalLBound3d(2)+1, & + totalUBound3d(3)-totalLBound3d(3)+1/) + ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3)/),count=count3d ); 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 + count3d = (/totalUBound3d(1)-totalLBound3d(1)+1, & + totalUBound3d(2)-totalLBound3d(2)+1, & + totalUBound3d(3)-totalLBound3d(3)+1/) + ncerr = nf90_put_var(ncid, varids(i), values=arrayr8_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3)/),count=count3d ); NC_ERR_STOP(ncerr) + end if + + end if !end fldlev(i) + + end do ! end fieldCount + +! if(allocated(arrayr4))deallocate(arrayr4) +! if(allocated(arrayr8))deallocate(arrayr8) +! if(allocated(arrayr4_3d))deallocate(arrayr4_3d) +! if(allocated(arrayr4_3d_save))deallocate(arrayr4_3d_save) +! if(allocated(arrayr8_3d))deallocate(arrayr8_3d) + + deallocate(fcstField) + deallocate(varids) + + if (mype==0) then + ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) + end if + + end subroutine write_netcdf_parallel + +!---------------------------------------------------------------------------------------- + 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' .or. ideflate == 0) 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 + ncerr = nf90_def_dim(ncid, trim(dim_name), NF90_UNLIMITED, dimid); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_dim(ncid, trim(dim_name), n, dimid); NC_ERR_STOP(ncerr) + end if + + if (typekind==ESMF_TYPEKIND_R8) then + ncerr = nf90_def_var(ncid, dim_name, NF90_REAL8, dimids=(/dimid/), varid=dim_varid); NC_ERR_STOP(ncerr) + allocate(valueListR8(n)) + call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & + name=trim(dim_name), valueList=valueListR8, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, dim_varid, values=valueListR8 ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + deallocate(valueListR8) + else if (typekind==ESMF_TYPEKIND_R4) then + ncerr = nf90_def_var(ncid, dim_name, NF90_REAL4, dimids=(/dimid/), varid=dim_varid); NC_ERR_STOP(ncerr) + allocate(valueListR4(n)) + call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & + name=trim(dim_name), valueList=valueListR4, rc=rc); ESMF_ERR_RETURN(rc) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4 ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + deallocate(valueListR4) + else + write(0,*)'Error in module_write_netcdf.F90(add_dim) unknown typekind for ',trim(dim_name) + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + call get_grid_attr(grid, dim_name, ncid, dim_varid, rc) + + end subroutine add_dim +! +!---------------------------------------------------------------------------------------- + subroutine nccheck(status) + use netcdf + implicit none + integer, intent (in) :: status + + if (status /= nf90_noerr) then + write(0,*) status, trim(nf90_strerror(status)) + stop "stopped" + end if + end subroutine nccheck + + elemental real function quantized(dataIn, nbits, dataMin, dataMax) + integer, intent(in) :: nbits + real(4), intent(in) :: dataIn, dataMin, dataMax + real(4) offset, scale_fact + ! convert data to 32 bit integers in range 0 to 2**nbits-1, then cast + ! cast back to 32 bit floats (data is then quantized in steps + ! proportional to 2**nbits so last 32-nbits in floating + ! point representation should be zero for efficient zlib compression). + scale_fact = (dataMax - dataMin) / (2**nbits-1); offset = dataMin + quantized = scale_fact*(nint((dataIn - offset) / scale_fact)) + offset + end function quantized + +end module module_write_netcdf_parallel diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index aad3991d0..6e6bc3da8 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -43,6 +43,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 ! !----------------------------------------------------------------------- ! @@ -183,7 +184,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) character(128) :: FBlist_outfilename(100), outfile_name character(128),dimension(:,:), allocatable :: outfilename real(8), dimension(:), allocatable :: slat - real, dimension(:), allocatable :: lat, lon, axesdata + real(8), dimension(:), allocatable :: lat, lon real(ESMF_KIND_R8), dimension(:,:), pointer :: lonPtr, latPtr real(ESMF_KIND_R8) :: rot_lon, rot_lat real(ESMF_KIND_R8) :: geo_lon, geo_lat @@ -358,13 +359,13 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc) wrt_int_state%latstart = lat(1) wrt_int_state%latlast = lat(jmo) do j=1,imo - lon(j) = 360./real(imo) *real(j-1) + lon(j) = 360.d0/real(imo,8) *real(j-1,8) enddo wrt_int_state%lonstart = lon(1) wrt_int_state%lonlast = lon(imo) do j=lbound(latPtr,2),ubound(latPtr,2) do i=lbound(lonPtr,1),ubound(lonPtr,1) - lonPtr(i,j) = 360./real(imo) * (i-1) + lonPtr(i,j) = 360.d0/real(imo,8) * real(i-1,8) latPtr(i,j) = lat(j) enddo enddo @@ -1434,6 +1435,17 @@ 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_parallel') then + + 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,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) == 'netcdf_esmf') then wbeg = MPI_Wtime() @@ -1622,13 +1634,14 @@ subroutine recover_fields(file_bundle,rc) character(100) fieldName,uwindname,vwindname type(ESMF_Field), allocatable :: fcstField(:) real(ESMF_KIND_R8), dimension(:,:), pointer :: lon, lat + real(ESMF_KIND_R8), dimension(:,:), pointer :: lonloc, latloc real(ESMF_KIND_R4), dimension(:,:), pointer :: pressfc real(ESMF_KIND_R4), dimension(:,:), pointer :: uwind2dr4,vwind2dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: uwind3dr4,vwind3dr4 real(ESMF_KIND_R4), dimension(:,:,:), pointer :: cart3dPtr2dr4 real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: cart3dPtr3dr4 real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: cart3dPtr3dr8 - save lon, lat + save lonloc, latloc real(ESMF_KIND_R8) :: coslon, sinlon, sinlat ! ! get filed count @@ -1648,9 +1661,18 @@ subroutine recover_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return - lon = lon * pi/180. -! print *,'in 3DCartesian2wind, lon dim=',lbound(lon,1),ubound(lon,1),lbound(lon,2),ubound(lon,2), & -! 'lon=',lon(lbound(lon,1),lbound(lon,2)), lon(ubound(lon,1),ubound(lon,2)) + allocate(lonloc(lbound(lon,1):ubound(lon,1),lbound(lon,2):ubound(lon,2))) + istart = lbound(lon,1) + iend = ubound(lon,1) + jstart = lbound(lon,2) + jend = ubound(lon,2) +!$omp parallel do default(none) shared(lon,lonloc,jstart,jend,istart,iend) & +!$omp private(i,j) + do j=jstart,jend + do i=istart,iend + lonloc(i,j) = lon(i,j) * pi/180. + enddo + enddo CALL ESMF_LogWrite("call recover field get coord 2",ESMF_LOGMSG_INFO,rc=RC) @@ -1658,6 +1680,19 @@ subroutine recover_fields(file_bundle,rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + allocate(latloc(lbound(lat,1):ubound(lat,1),lbound(lat,2):ubound(lat,2))) + istart = lbound(lat,1) + iend = ubound(lat,1) + jstart = lbound(lat,2) + jend = ubound(lat,2) +!$omp parallel do default(none) shared(lat,latloc,jstart,jend,istart,iend) & +!$omp private(i,j) + do j=jstart,jend + do i=istart,iend + latloc(i,j) = lat(i,j) * pi/180.d0 + enddo + enddo + lat = lat * pi/180. ! print *,'in 3DCartesian2wind, lat dim=',lbound(lat,1),ubound(lat,1),lbound(lat,2),ubound(lat,2), & ! 'lat=',lat(lbound(lon,1),lbound(lon,2)), lat(ubound(lon,1),ubound(lon,2)) @@ -1718,18 +1753,18 @@ subroutine recover_fields(file_bundle,rc) ! update u , v wind !$omp parallel do default(shared) private(i,j,k,coslon,sinlon,sinlat) do k=kstart,kend -!!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lon,lat,cart3dPtr3dr4,jstart,jend,istart,iend,k) & -!!$omp private(i,j,coslon,sinlon,sinlat) +!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lonloc,latloc,cart3dPtr3dr4,jstart,jend,istart,iend,k) & +!$omp private(i,j,coslon,sinlon,sinlat) do j=jstart, jend do i=istart, iend - coslon = cos(lon(i,j)) - sinlon = sin(lon(i,j)) - sinlat = sin(lat(i,j)) + coslon = cos(lonloc(i,j)) + sinlon = sin(lonloc(i,j)) + sinlat = sin(latloc(i,j)) uwind3dr4(i,j,k) = cart3dPtr3dr4(1,i,j,k) * coslon & + cart3dPtr3dr4(2,i,j,k) * sinlon vwind3dr4(i,j,k) =-cart3dPtr3dr4(1,i,j,k) * sinlat*sinlon & + cart3dPtr3dr4(2,i,j,k) * sinlat*coslon & - + cart3dPtr3dr4(3,i,j,k) * cos(lat(i,j)) + + cart3dPtr3dr4(3,i,j,k) * cos(latloc(i,j)) enddo enddo enddo @@ -1749,18 +1784,18 @@ subroutine recover_fields(file_bundle,rc) call ESMF_FieldGet(ufield, localDe=0, farrayPtr=uwind2dr4,rc=rc) call ESMF_FieldGet(vfield, localDe=0, farrayPtr=vwind2dr4,rc=rc) ! update u , v wind -!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lon,lat,cart3dPtr2dr4,jstart,jend,istart,iend) & +!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lonloc,latloc,cart3dPtr2dr4,jstart,jend,istart,iend) & !$omp private(i,j,k,coslon,sinlon,sinlat) do j=jstart, jend do i=istart, iend - coslon = cos(lon(i,j)) - sinlon = sin(lon(i,j)) - sinlat = sin(lat(i,j)) + coslon = cos(lonloc(i,j)) + sinlon = sin(lonloc(i,j)) + sinlat = sin(latloc(i,j)) uwind2dr4(i,j) = cart3dPtr2dr4(1,i,j) * coslon & + cart3dPtr2dr4(2,i,j) * sinlon vwind2dr4(i,j) =-cart3dPtr2dr4(1,i,j) * sinlat*sinlon & + cart3dPtr2dr4(2,i,j) * sinlat*coslon & - + cart3dPtr2dr4(3,i,j) * cos(lat(i,j)) + + cart3dPtr2dr4(3,i,j) * cos(latloc(i,j)) enddo enddo endif diff --git a/io/post_gfs.F90 b/io/post_gfs.F90 index 98d4fef50..7559f52ba 100644 --- a/io/post_gfs.F90 +++ b/io/post_gfs.F90 @@ -1146,7 +1146,7 @@ subroutine set_postvars_gfs(wrt_int_state,mpicomp,setvar_atmfile, & do i=ista, iend if (arrayr42d(i,j) /= spval) then !set range within (0,1) - sr(i,j) = min(1.,max(0.,sr(i,j))) + sr(i,j) = min(1.,max(0.,arrayr42d(i,j))) else sr(i,j) = spval endif From 6b9b8651c5cca43f5951a18e787b3f036b5855d9 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 8 Jan 2020 22:00:06 +0000 Subject: [PATCH 02/32] some changes to get it to (almost) run on hera --- io/module_write_netcdf_parallel.F90 | 63 +++++++++++++---------------- 1 file changed, 29 insertions(+), 34 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index c2ea784fd..1c0d2ab7e 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -116,20 +116,12 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i ! 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 -!not build with pnetcdf -! if (ideflate == 0) then -! ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), & -! comm=mpi_comm, info = MPI_INFO_NULL,ncid=ncid); NC_ERR_STOP(ncerr) -! ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) -! else -! ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), & -! ncid=ncid); NC_ERR_STOP(ncerr) -! ! if compression on use HDF5 format with default _FillValue -! endif -! - ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_NETCDF4), & +! 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 @@ -193,6 +185,8 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i stop end if end if + ! time dimension is unlimited, so must use collective access. + !ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) ! define variable attributes call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & @@ -225,7 +219,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, & rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then + 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 @@ -249,8 +243,6 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i ! write grid_xt, grid_yt values call ESMF_GridGetCoord(wrtGrid, coordDim=1, farrayPtr=arrayr8, rc=rc); ESMF_ERR_RETURN(rc) -! call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) -! if (mype==0) then if (trim(output_grid) == 'gaussian_grid' .or. trim(output_grid) == 'regional_latlon') then istart = lbound(arrayr8,1) iend = ubound(arrayr8,1) @@ -286,22 +278,22 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i start = (/istart,jstart/) count = (/iend-istart+1,jend-jstart+1/) ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8, start=start, count=count ); NC_ERR_STOP(ncerr) + print *,'done write lon' call ESMF_GridGetCoord(wrtGrid, coordDim=2, farrayPtr=arrayr8, rc=rc); ESMF_ERR_RETURN(rc) -! call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc) if (trim(output_grid) == 'gaussian_grid' .or. trim(output_grid) == 'regional_latlon') then - istart = lbound(arrayr8,1) - iend = ubound(arrayr8,1) - jstart = lbound(arrayr8,2) - jend = ubound(arrayr8,2) - print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend - start1d(1) = jstart - count1d(1) = jend - jstart + 1 - ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:),start=start1d,count=count1d ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + istart = lbound(arrayr8,1) + iend = ubound(arrayr8,1) + jstart = lbound(arrayr8,2) + jend = ubound(arrayr8,2) + print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend + start1d(1) = jstart + count1d(1) = jend - jstart + 1 + ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:),start=start1d,count=count1d ); NC_ERR_STOP(ncerr) + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'rotated_latlon') then do j=1,jm y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1) @@ -324,6 +316,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i start = (/istart,jstart/) count = (/iend-istart+1, jend-jstart+1/) ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8, start=start, count=count ); NC_ERR_STOP(ncerr) + print *,'done write lat' do i=1, fieldCount @@ -345,7 +338,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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) - print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d + print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d if (ideflate > 0 .and. nbits > 0) then ! Lossy compression if nbits>0. ! The floating point data is quantized to improve compression @@ -396,9 +389,8 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i deallocate(fcstField) deallocate(varids) - if (mype==0) then ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) - end if + print *,'netcdf parallel close, finished write_netcdf_parallel' end subroutine write_netcdf_parallel @@ -509,7 +501,7 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) else if (typekind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then + 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) @@ -551,7 +543,8 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) 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) + !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 @@ -562,6 +555,7 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) 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_var_par_access(ncid, dim_varid, NF90_COLLECTIVE) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR8 ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR8) @@ -571,6 +565,7 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) 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_var_par_access(ncid, dim_varid, NF90_COLLECTIVE) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4 ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR4) From ccdbf2743613ac7c16c5b5b3533e6aab54a7ff7d Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 9 Jan 2020 17:46:19 +0000 Subject: [PATCH 03/32] lat should be in degrees - remove conversion to radians. --- io/module_wrt_grid_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 6e6bc3da8..2390c0b32 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1693,7 +1693,7 @@ subroutine recover_fields(file_bundle,rc) enddo enddo - lat = lat * pi/180. + !lat = lat * pi/180. ! print *,'in 3DCartesian2wind, lat dim=',lbound(lat,1),ubound(lat,1),lbound(lat,2),ubound(lat,2), & ! 'lat=',lat(lbound(lon,1),lbound(lon,2)), lat(ubound(lon,1),ubound(lon,2)) first_getlatlon = .false. From b3e52a0079b0a4b7f0f7481d51566eb4969c25bd Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 9 Jan 2020 17:47:07 +0000 Subject: [PATCH 04/32] import updates from jswhit/fv3atm --- io/module_write_netcdf.F90 | 116 +++++++++++++++++++------------------ 1 file changed, 60 insertions(+), 56 deletions(-) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 52cf740e6..2a05bf337 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -42,7 +42,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save real(8), dimension(:,:,:), allocatable :: arrayr8_3d - real(8) x(im),y(jm) + real(8) rad2dg,x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound @@ -56,7 +56,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc character(len=ESMF_MAXSTR) :: attName, fldName integer :: varival - real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax + real(4) :: varr4val, scale_fact, offset, dataMin, dataMax + real(4), allocatable, dimension(:) :: compress_err real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval @@ -71,9 +72,11 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! !! ! + rad2dg = 45./atan(1.0) 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)) @@ -115,15 +118,9 @@ 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(NF90_CLOBBER,NF90_64BIT_OFFSET), & - ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr) - else - ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), & - ncid=ncid); NC_ERR_STOP(ncerr) - ! if compression on use HDF5 format with default _FillValue - endif + ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_NETCDF4), & + 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) @@ -155,12 +152,14 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! define variables if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter on for lossless compression - ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate) - NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i), & + shuffle=.true.,deflate_level=ideflate,& + chunksizes=(/im,jm,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) endif else if (typekind == ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & @@ -171,13 +170,15 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter off since lossy compression used - ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate) - NC_ERR_STOP(ncerr) - endif + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & + shuffle=.false.,deflate_level=ideflate,& + chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + 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) @@ -218,8 +219,8 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, & rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then - ! FIXME: _FillValue must be cast to var type when using NF90_NETCDF4 + if (trim(attName) /= '_FillValue') then + ! FIXME: _FillValue must be cast to var type for recent versions of netcdf ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr) endif @@ -235,6 +236,25 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end do ! i=1,fieldCount + ! write grid_xt, grid_yt attributes + if (trim(output_grid) == 'gaussian_grid' .or. & + trim(output_grid) == 'regional_latlon') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'rotated_latlon') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) + else if (trim(output_grid) == 'lambert_conformal') then + ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) + endif + ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr) end if @@ -247,28 +267,16 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'rotated_latlon') then do i=1,im x(i) = lon1 + (lon2-lon1)/(im-1) * (i-1) enddo ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do i=1,im x(i) = dx * (i-1) enddo ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif @@ -279,28 +287,16 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc if (trim(output_grid) == 'gaussian_grid' .or. & trim(output_grid) == 'regional_latlon') then ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'rotated_latlon') then do j=1,jm y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1) enddo ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do j=1,jm y(j) = dy * (j-1) enddo ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr) endif @@ -343,11 +339,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax) ! compute max abs compression error, save as a variable ! attribute. - compress_err = maxval(abs(arrayr4_3d_save-arrayr4_3d)) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d)) endif ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/1,1,1/),count=(/im,jm,lm,1/) ); NC_ERR_STOP(ncerr) end if @@ -362,6 +354,17 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc end do + if (ideflate > 0 .and. nbits > 0 .and. mype == 0) then + ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) + do i=1, fieldCount + if (compress_err(i) > 0) then + ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err(i)); NC_ERR_STOP(ncerr) + ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) + endif + enddo + ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + endif + deallocate(arrayr4) deallocate(arrayr8) deallocate(arrayr4_3d,arrayr4_3d_save) @@ -369,6 +372,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) @@ -483,9 +487,9 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc) else if (typekind==ESMF_TYPEKIND_R8) then call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", & name=trim(attName), value=varr8val, rc=rc); ESMF_ERR_RETURN(rc) - if (trim(attName) /= '_FillValue' .or. ideflate == 0) then - ! FIXME: _FillValue must be cast to var type when using - ! NF90_NETCDF4. Until this is fixed, using netCDF default _FillValue. + if (trim(attName) /= '_FillValue') then + ! FIXME: _FillValue must be cast to var type for recent versions + ! of netcdf ncerr = nf90_put_att(ncid, varid, trim(attName(ind+1:len(attName))), varr8val); NC_ERR_STOP(ncerr) endif From a201d6c11fe8df9cb87ba7f46dacb61108a24117 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 9 Jan 2020 17:47:28 +0000 Subject: [PATCH 05/32] more bug fixes - now works on hera --- io/module_write_netcdf_parallel.F90 | 183 ++++++++++++---------------- 1 file changed, 78 insertions(+), 105 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 1c0d2ab7e..1f488caee 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -43,7 +43,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i real(ESMF_KIND_R4), dimension(:,:,:), pointer :: arrayr4_3d,arrayr4_3d_save real(ESMF_KIND_R8), dimension(:,:,:), pointer :: arrayr8_3d - real(8) rad2dg,x(im),y(jm) + real(8) x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound @@ -55,29 +55,26 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i integer :: attcount character(len=ESMF_MAXSTR) :: attName, fldName - integer :: start1d(1), count1d(1), start(2), count(2), start3d(3), count3d(3) - integer :: totalLBound2d(2), totalUBound2d(2),totalLBound3d(3),totalUBound3d(3) + integer :: totalLBound2d(2),totalUBound2d(2),totalLBound3d(3),totalUBound3d(3) integer :: varival - real(4) :: varr4val, scale_fact, compress_err, offset, dataMin, dataMax + real(4) :: varr4val, scale_fact, offset, dataMin, dataMax + real(4), allocatable, dimension(:) :: compress_err real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval 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, dimension(:), allocatable :: varids ! -!! -! - rad2dg = 45./atan(1.0) - 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)) @@ -111,11 +108,6 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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 for parallel access ncerr = nf90_create(trim(filename),& @@ -129,12 +121,16 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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) @@ -154,12 +150,14 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i ! define variables if (fldlev(i) == 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter on for lossless compression - ncerr = nf90_def_var_deflate(ncid, varids(i), 1, 1, ideflate) - NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i), & + shuffle=.true.,deflate_level=ideflate,& + chunksizes=(/im,jm,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + else + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) endif else if (typekind == ESMF_TYPEKIND_R8) then ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, & @@ -170,12 +168,14 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i end if else if (fldlev(i) > 1) then if (typekind == ESMF_TYPEKIND_R4) then - ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & - (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr) if (ideflate > 0) then - ! shuffle filter off since lossy compression used - ncerr = nf90_def_var_deflate(ncid, varids(i), 0, 1, ideflate) - NC_ERR_STOP(ncerr) + ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & + (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & + shuffle=.true.,deflate_level=ideflate,& + chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); 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, & @@ -185,8 +185,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i stop end if end if - ! time dimension is unlimited, so must use collective access. - !ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) + ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) ! define variable attributes call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & @@ -236,87 +235,66 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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 - istart = lbound(arrayr8,1) - iend = ubound(arrayr8,1) - jstart = lbound(arrayr8,2) - jend = ubound(arrayr8,2) - print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend - start1d(1) = istart - count1d(1) = iend - istart + 1 - ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(istart:iend,1),start=start1d, count=count1d ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + 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) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do i=1,im x(i) = dx * (i-1) enddo ncerr = nf90_put_var(ncid, im_varid, values=x ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - start = (/istart,jstart/) - count = (/iend-istart+1,jend-jstart+1/) - ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8, start=start, count=count ); NC_ERR_STOP(ncerr) - print *,'done write lon' + 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 - istart = lbound(arrayr8,1) - iend = ubound(arrayr8,1) - jstart = lbound(arrayr8,2) - jend = ubound(arrayr8,2) - print *,'in write netcdf mpi, istart=',istart,iend,'jstart=',jstart,jend - start1d(1) = jstart - count1d(1) = jend - jstart + 1 - ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:),start=start1d,count=count1d ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) + 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) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) else if (trim(output_grid) == 'lambert_conformal') then do j=1,jm y(j) = dy * (j-1) enddo ncerr = nf90_put_var(ncid, jm_varid, values=y ); NC_ERR_STOP(ncerr) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) endif - start = (/istart,jstart/) - count = (/iend-istart+1, jend-jstart+1/) - ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8, start=start, count=count ); NC_ERR_STOP(ncerr) - print *,'done write lat' + ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8, start=(/istart,jstart/)); NC_ERR_STOP(ncerr) do i=1, fieldCount @@ -325,20 +303,16 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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 - count3d=(/totalUBound2d(1)-totalLBound2d(1)+1, & - totalUBound2d(2)-totalLBound2d(2)+1, 1/) - ncerr = nf90_put_var(ncid, varids(i), values=arrayr4, start=(/totalLBound2d(1),totalLBound2d(2),1/),count=count3d); NC_ERR_STOP(ncerr) + !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) - count3d=(/totalUBound2d(1)-totalLBound2d(1)+1, & - totalUBound2d(2)-totalLBound2d(2)+1, 1/) - ncerr = nf90_put_var(ncid, varids(i), values=arrayr8, start=(/totalLBound2d(1),totalLBound2d(2),1/),count=count3d); NC_ERR_STOP(ncerr) + 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) - print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d + !print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d if (ideflate > 0 .and. nbits > 0) then ! Lossy compression if nbits>0. ! The floating point data is quantized to improve compression @@ -357,40 +331,37 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i arrayr4_3d = quantized(arrayr4_3d_save, nbits, dataMin, dataMax) ! compute max abs compression error, save as a variable ! attribute. - compress_err = maxval(abs(arrayr4_3d_save-arrayr4_3d)) - ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'max_abs_compression_error', compress_err); NC_ERR_STOP(ncerr) - ncerr = nf90_put_att(ncid, varids(i), 'nbits', nbits); NC_ERR_STOP(ncerr) - ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr) - endif - count3d = (/totalUBound3d(1)-totalLBound3d(1)+1, & - totalUBound3d(2)-totalLBound3d(2)+1, & - totalUBound3d(3)-totalLBound3d(3)+1/) - ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3)/),count=count3d ); NC_ERR_STOP(ncerr) + compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d)) + endif + ncerr = nf90_put_var(ncid, varids(i), values=arrayr4_3d, start=(/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 - count3d = (/totalUBound3d(1)-totalLBound3d(1)+1, & - totalUBound3d(2)-totalLBound3d(2)+1, & - totalUBound3d(3)-totalLBound3d(3)+1/) - ncerr = nf90_put_var(ncid, varids(i), values=arrayr8_3d, start=(/totalLBound3d(1),totalLBound3d(2),totalLBound3d(3)/),count=count3d ); NC_ERR_STOP(ncerr) + !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(allocated(arrayr4))deallocate(arrayr4) -! if(allocated(arrayr8))deallocate(arrayr8) -! if(allocated(arrayr4_3d))deallocate(arrayr4_3d) -! if(allocated(arrayr4_3d_save))deallocate(arrayr4_3d_save) -! if(allocated(arrayr8_3d))deallocate(arrayr8_3d) + 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) - print *,'netcdf parallel close, finished write_netcdf_parallel' + !call mpi_barrier(mpi_comm,ierr) + !print *,'netcdf parallel close, finished write_netcdf_parallel' end subroutine write_netcdf_parallel @@ -543,6 +514,8 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) 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 @@ -551,21 +524,21 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc) 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_var_par_access(ncid, dim_varid, NF90_COLLECTIVE) 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_var_par_access(ncid, dim_varid, NF90_COLLECTIVE) ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4 ); NC_ERR_STOP(ncerr) ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr) deallocate(valueListR4) From 656c5466d852976767159ad50b110955cacdecf6 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 9 Jan 2020 18:03:31 +0000 Subject: [PATCH 06/32] specify collective access if compression turned on. --- io/module_write_netcdf_parallel.F90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 1f488caee..dc47bce9b 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -155,13 +155,17 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i (/im_dimid,jm_dimid,time_dimid/), varids(i), & shuffle=.true.,deflate_level=ideflate,& chunksizes=(/im,jm,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + ! 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 @@ -173,19 +177,22 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & shuffle=.true.,deflate_level=ideflate,& chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + ! 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 - ncerr = nf90_var_par_access(ncid, varids(i), NF90_INDEPENDENT) ! define variable attributes call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", & From ade86ed515d3676ceeab7f4c859f636550980f22 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 9 Jan 2020 22:26:54 +0000 Subject: [PATCH 07/32] use classic model --- io/module_write_netcdf.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 2a05bf337..5b948f092 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -118,8 +118,9 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! create netcdf file and enter define mode if (mype==0) then - ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_NETCDF4), & - ncid=ncid); NC_ERR_STOP(ncerr) + 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 From 49ffb5e4f7f5444e24cc5d5993194e9011cc31ea Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 10 Jan 2020 14:12:50 +0000 Subject: [PATCH 08/32] bug fixes for parallel IO with compression --- io/module_write_netcdf_parallel.F90 | 71 +++++++++++++++-------------- 1 file changed, 36 insertions(+), 35 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index dc47bce9b..23bda3391 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -34,7 +34,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i integer, optional,intent(out) :: rc ! !** local vars - integer :: i,j,m,n,k,istart,iend,jstart,jend + integer :: i,j,m,n,k,istart,iend,jstart,jend,i1,i2,j1,j2,k1,k2 integer :: lm integer, dimension(:), allocatable :: fldlev @@ -58,8 +58,9 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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(4) :: varr4val, scale_fact, offset, dataMin, dataMax, dataMin1, & + dataMax1 + real(4), allocatable, dimension(:) :: compress_err,compress_err1 real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval @@ -75,6 +76,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) allocate(compress_err(fieldCount)); compress_err=-999. + allocate(compress_err1(fieldCount)); compress_err1=0. allocate(fldlev(fieldCount)) ; fldlev = 0 allocate(fcstField(fieldCount)) allocate(varids(fieldCount)) @@ -319,26 +321,37 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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) - !print *,'field name=',trim(fldName),'bound=',totalLBound3d,'ubound=',totalUBound3d 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)) + i1=totalLBound3d(1);i2=totalUBound3d(1) + j1=totalLBound3d(2);j2=totalUBound3d(2) + k1=totalLBound3d(3);k2=totalUBound3d(3) + dataMax1 = maxval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) + dataMin1 = minval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) + call mpi_allreduce(dataMax1,dataMax,1,mpi_real4,mpi_max,mpi_comm,ierr) + call mpi_allreduce(dataMin1,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,mold=arrayr4_3d) + arrayr4_3d_save=arrayr4_3d + arrayr4_3d=scale_fact*(nint((arrayr4_3d-offset)/scale_fact))+offset + ! compute max abs compression error. + compress_err1(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(compress_err1(i),compress_err(i),1,mpi_real4,mpi_max,mpi_comm,ierr) + print *,'field name=',trim(fldName),dataMin,dataMax,compress_err(i),& + minval(arrayr4_3d(i1:i2,j1:j2,k1:k2)),maxval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) 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 @@ -364,7 +377,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i deallocate(fcstField) deallocate(varids) - deallocate(compress_err) + deallocate(compress_err,compress_err1) ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) !call mpi_barrier(mpi_comm,ierr) @@ -570,16 +583,4 @@ subroutine nccheck(status) end if end subroutine nccheck - elemental real function quantized(dataIn, nbits, dataMin, dataMax) - integer, intent(in) :: nbits - real(4), intent(in) :: dataIn, dataMin, dataMax - real(4) offset, scale_fact - ! convert data to 32 bit integers in range 0 to 2**nbits-1, then cast - ! cast back to 32 bit floats (data is then quantized in steps - ! proportional to 2**nbits so last 32-nbits in floating - ! point representation should be zero for efficient zlib compression). - scale_fact = (dataMax - dataMin) / (2**nbits-1); offset = dataMin - quantized = scale_fact*(nint((dataIn - offset) / scale_fact)) + offset - end function quantized - end module module_write_netcdf_parallel From 134ebc0dc51570edacfda3f1e91549f4b91dda81 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 10 Jan 2020 14:47:08 +0000 Subject: [PATCH 09/32] fix calculation of max compression error --- io/module_write_netcdf_parallel.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 23bda3391..33d2b4e87 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -342,16 +342,15 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i ! rescaling (plus it allows the ability to adjust the packing ! range) scale_fact = (dataMax - dataMin) / (2**nbits-1); offset = dataMin - allocate(arrayr4_3d_save,mold=arrayr4_3d) - arrayr4_3d_save=arrayr4_3d - arrayr4_3d=scale_fact*(nint((arrayr4_3d-offset)/scale_fact))+offset + 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(i1:i2,j1:j2,k1:k2)=scale_fact*(nint((arrayr4_3d(i1:i2,j1:j2,k1:k2)-offset)/scale_fact))+offset ! compute max abs compression error. compress_err1(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(compress_err1(i),compress_err(i),1,mpi_real4,mpi_max,mpi_comm,ierr) - print *,'field name=',trim(fldName),dataMin,dataMax,compress_err(i),& - minval(arrayr4_3d(i1:i2,j1:j2,k1:k2)),maxval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) + !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 From 5fce869943c97ba1f45ab5649f6258e53ec32da7 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 10 Jan 2020 15:56:19 +0000 Subject: [PATCH 10/32] turn off shuffle filter --- io/module_write_netcdf_parallel.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 33d2b4e87..eced14c67 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -177,7 +177,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i if (ideflate > 0) then ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & - shuffle=.true.,deflate_level=ideflate,& + shuffle=.false.,deflate_level=ideflate,& chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) From 30516fc6c45458f19f453ab76a1e7425e0f8bbc4 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 10 Jan 2020 17:15:55 +0000 Subject: [PATCH 11/32] code simplification --- io/module_write_netcdf_parallel.F90 | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index eced14c67..66b126f8f 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -58,9 +58,8 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i integer :: totalLBound2d(2),totalUBound2d(2),totalLBound3d(3),totalUBound3d(3) integer :: varival - real(4) :: varr4val, scale_fact, offset, dataMin, dataMax, dataMin1, & - dataMax1 - real(4), allocatable, dimension(:) :: compress_err,compress_err1 + real(4) :: varr4val, scale_fact, offset, dataMin, dataMax + real(4), allocatable, dimension(:) :: compress_err real(8) :: varr8val character(len=ESMF_MAXSTR) :: varcval @@ -76,7 +75,6 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) allocate(compress_err(fieldCount)); compress_err=-999. - allocate(compress_err1(fieldCount)); compress_err1=0. allocate(fldlev(fieldCount)) ; fldlev = 0 allocate(fcstField(fieldCount)) allocate(varids(fieldCount)) @@ -325,10 +323,10 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i i1=totalLBound3d(1);i2=totalUBound3d(1) j1=totalLBound3d(2);j2=totalUBound3d(2) k1=totalLBound3d(3);k2=totalUBound3d(3) - dataMax1 = maxval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) - dataMin1 = minval(arrayr4_3d(i1:i2,j1:j2,k1:k2)) - call mpi_allreduce(dataMax1,dataMax,1,mpi_real4,mpi_max,mpi_comm,ierr) - call mpi_allreduce(dataMin1,dataMin,1,mpi_real4,mpi_min,mpi_comm,ierr) + 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 @@ -344,12 +342,12 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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(i1:i2,j1:j2,k1:k2)=scale_fact*(nint((arrayr4_3d(i1:i2,j1:j2,k1:k2)-offset)/scale_fact))+offset + arrayr4_3d = scale_fact*(nint((arrayr4_3d_save - offset) / scale_fact)) + offset ! compute max abs compression error. - compress_err1(i) = & + 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(compress_err1(i),compress_err(i),1,mpi_real4,mpi_max,mpi_comm,ierr) + 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) @@ -376,7 +374,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i deallocate(fcstField) deallocate(varids) - deallocate(compress_err,compress_err1) + deallocate(compress_err) ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr) !call mpi_barrier(mpi_comm,ierr) From 930d622bde526a96b2810177df623ca765f86e75 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 17 Jan 2020 15:45:39 +0000 Subject: [PATCH 12/32] remove debug print --- io/module_write_netcdf.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 5b948f092..3be91df60 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -341,6 +341,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc ! compute max abs compression error, save as a variable ! attribute. compress_err(i) = maxval(abs(arrayr4_3d_save-arrayr4_3d)) + !print *,'field name=',trim(fldName),dataMin,dataMax,compress_err(i) 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 From 6cb687d4558a12baed86987a631cadc3de82e01f Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Wed, 22 Jan 2020 20:50:11 +0000 Subject: [PATCH 13/32] don't use parallel IO for 2d file (since it seems to increase run time) --- io/module_wrt_grid_comp.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 2390c0b32..ab075db5b 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1438,8 +1438,13 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) else if (trim(output_file) == 'netcdf_parallel') then 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,rc) + if (nbdl == 1) then ! only use parallel write for 3D file + call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & + trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) + else + call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & + wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) + endif wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & From c883717a384c0494e1c55397874cdccbf69ee7b4 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 23 Jan 2020 22:29:19 +0000 Subject: [PATCH 14/32] allow multiple values of output_file, as long as they all start with 'netcdf' --- fv3_cap.F90 | 38 ++++++++++++++++++++++++++++++++----- io/module_fv3_io_def.F90 | 2 +- io/module_wrt_grid_comp.F90 | 37 ++++++++++++++++-------------------- 3 files changed, 50 insertions(+), 27 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index a099ae32e..70594031f 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -235,7 +235,7 @@ 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 logical :: isPetLocal logical :: OPENED @@ -346,8 +346,38 @@ 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)) + CALL ESMF_ConfigFindLabel(CF,'output_file:',rc=RC) + ierr = 0 + do i=1,num_files + CALL ESMF_ConfigGetAttribute(config=CF,value=output_file(i), rc=rc) + ! if only 1 output_file specified, copy to all elements + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + ierr = 1 + if (i .eq. 1) then + return + else + output_file(i) = output_file(i-1) + endif + endif + ! specifying multiple values of output_file will only work if + ! they all start with 'netcdf' - need to add a check for this + ! Mixing netcdf and nemsio will not work. + if (i > 1 .and. ierr == 0) then + if (output_file(i)(1:6) .ne. 'netcdf' .or. & + output_file(i-1)(1:6) .ne. 'netcdf') then + write(0,*)'fv3_cap.F90:muliple values of output_file must all begin with netcdf' + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + endif + enddo + 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 +389,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/module_fv3_io_def.F90 b/io/module_fv3_io_def.F90 index 17e9f308d..55249b14c 100644 --- a/io/module_fv3_io_def.F90 +++ b/io/module_fv3_io_def.F90 @@ -17,13 +17,13 @@ module module_fv3_io_def integer :: num_files character(255) :: app_domain character(255) :: output_grid - character(255) :: output_file integer :: imo,jmo 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_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index ab075db5b..9309d7707 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1122,14 +1122,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) ! !----------------------------------------------------------------------- @@ -1363,7 +1363,7 @@ 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 + 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 +1413,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,7 +1424,7 @@ 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), & @@ -1435,23 +1435,18 @@ 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_parallel') then + else if (trim(output_file(nbdl)) == 'netcdf_parallel') then wbeg = MPI_Wtime() - if (nbdl == 1) then ! only use parallel write for 3D file - call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & - trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) - else - call write_netcdf(file_bundle,wrt_int_state%wrtFB(nbdl),trim(filename), & - wrt_mpi_comm,wrt_int_state%mype,imo,jmo,rc) - endif + call write_netcdf_parallel(file_bundle,wrt_int_state%wrtFB(nbdl), & + trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,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) == 'netcdf_esmf') then + else if (trim(output_file(nbdl)) == 'netcdf_esmf') then wbeg = MPI_Wtime() call ESMFproto_FieldBundleWrite(gridFB, filename=trim(filename), & @@ -1498,7 +1493,7 @@ 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), & From ffb33e008ca8cbf141dedc8cdbf4b41dd47dca23 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 24 Jan 2020 02:34:27 +0000 Subject: [PATCH 15/32] use default chunksize for 2d vars --- io/module_write_netcdf_parallel.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 66b126f8f..110c6b3a5 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -153,8 +153,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + shuffle=.true.,deflate_level=ideflate); NC_ERR_STOP(ncerr) ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else From 89184bddee1f5157b66244bcd743c146b08b7d6c Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 24 Jan 2020 17:18:57 +0000 Subject: [PATCH 16/32] delete commented out macro ESMF_ERR_ABORT --- io/module_write_netcdf.F90 | 4 ---- io/module_write_netcdf_parallel.F90 | 4 ---- 2 files changed, 8 deletions(-) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 3be91df60..aca6350f5 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) \ diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 110c6b3a5..5fe52b0bf 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.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) \ From 5af9ccf158b31879b01dfe0b9e0895ba452e109f Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 24 Jan 2020 17:20:27 +0000 Subject: [PATCH 17/32] delete rad2dg --- io/module_write_netcdf.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index aca6350f5..c8fc04ecf 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -38,7 +38,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save real(8), dimension(:,:,:), allocatable :: arrayr8_3d - real(8) rad2dg,x(im),y(jm) + real(8) x(im),y(jm) integer :: fieldCount, fieldDimCount, gridDimCount integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound @@ -65,10 +65,6 @@ 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 -! -!! -! - rad2dg = 45./atan(1.0) call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc) From b49ef58a94c41c58c18bb0cbf670cf7b9b484cdb Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sat, 25 Jan 2020 16:07:44 +0000 Subject: [PATCH 18/32] add module_write_netcdf_parallel.F90 --- io/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index 5eecdd9ea..c5ba4c240 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -12,6 +12,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 From 2e63f49f59de840577cd6225dc39e78b3dac8d7d Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sat, 25 Jan 2020 17:17:48 +0000 Subject: [PATCH 19/32] add option to build without parallel netcdf (-DNO_PARALLEL_NETCDF) --- io/CMakeLists.txt | 9 ++++++++- io/makefile | 12 +++++++++++- io/module_wrt_grid_comp.F90 | 7 +++++++ 3 files changed, 26 insertions(+), 2 deletions(-) diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index c5ba4c240..bf3a292b4 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -5,6 +5,13 @@ else() add_definitions(-DNO_INLINE_POST) endif() +if(PARALLEL_NETCDF) + set(NETCDF_PARALLEL_SRC model_write_netcdf_parallel.F90) +else() + set(NETCDF_PARALLEL_SRC model_write_netcdf_parallel_stub.F90) + add_definitions(-DNO_PARALLEL_NETCDF) +endif() + add_library( io @@ -12,7 +19,7 @@ add_library( FV3GFS_io.F90 module_write_nemsio.F90 module_write_netcdf.F90 - module_write_netcdf_parallel.F90 + ${NETCDF_PARALLEL_SRC} module_fv3_io_def.F90 module_write_internal_state.F90 module_wrt_grid_comp.F90 diff --git a/io/makefile b/io/makefile index b6a3946b0..adcb6cc38 100644 --- a/io/makefile +++ b/io/makefile @@ -29,6 +29,14 @@ POST_SRC = \ ./post_nems_routines.F90 endif +ifneq (,$(findstring NO_PARALLEL_NETCDF,$(CPPDEFS))) +NETCDF_PARALLEL_SRC = \ + ./module_write_netcdf_parallel_stub.F90 +else +NETCDF_PARALLEL_SRC = \ + ./module_write_netcdf_parallel.F90 +endif + SRCS_f = SRCS_f90 = @@ -41,7 +49,7 @@ SRCS_F90 = \ $(POST_SRC) \ ./module_write_nemsio.F90 \ ./module_write_netcdf.F90 \ - ./module_write_netcdf_parallel.F90 \ + $(NETCDF_PARALLEL_SRC) \ ./module_fv3_io_def.F90 \ ./module_write_internal_state.F90 \ ./module_wrt_grid_comp.F90 @@ -73,6 +81,8 @@ module_write_netcdf.o: 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_netcdf_parallel_stub.o: module_write_netcdf_parallel_stub.F90 + $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) -c module_write_netcdf_parallel_stub.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_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 9309d7707..7c65306a4 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1437,6 +1437,13 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) else if (trim(output_file(nbdl)) == 'netcdf_parallel') then +#ifdef NO_PARALLEL_NETCDF + rc = ESMF_RC_NOT_IMPL + print *,'netcdf_parallel post 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,rc) From 4a45855378122872a29be204ddfc6a84cf607576 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sat, 25 Jan 2020 17:37:44 +0000 Subject: [PATCH 20/32] fix typo --- io/module_wrt_grid_comp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 7c65306a4..33f4fa58f 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1439,7 +1439,7 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) #ifdef NO_PARALLEL_NETCDF rc = ESMF_RC_NOT_IMPL - print *,'netcdf_parallel post not available on this machine' + print *,'netcdf_parallel not available on this machine' if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, file=__FILE__)) return #endif From 14dfcab11c5a849e754a472d7e58d75c3ced137a Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sun, 26 Jan 2020 02:21:22 +0000 Subject: [PATCH 21/32] stub file for building without parallel netcdf lib --- io/module_write_netcdf_parallel_stub.F90 | 26 ++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 io/module_write_netcdf_parallel_stub.F90 diff --git a/io/module_write_netcdf_parallel_stub.F90 b/io/module_write_netcdf_parallel_stub.F90 new file mode 100644 index 000000000..3fe32e063 --- /dev/null +++ b/io/module_write_netcdf_parallel_stub.F90 @@ -0,0 +1,26 @@ +module module_write_netcdf_parallel + + use esmf + implicit none + private + public write_netcdf_parallel + + contains + +!---------------------------------------------------------------------------------------- + subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) +! + type(ESMF_FieldBundle), intent(in) :: fieldbundle + type(ESMF_FieldBundle), intent(in) :: wrtfb + character(*), intent(in) :: filename + integer, intent(in) :: mpi_comm + integer, intent(in) :: mype + integer, intent(in) :: im, jm + integer, optional,intent(out) :: rc + + print *,'in stub write_netcdf_parallel - model not built with parallel netcdf support, return' + + end subroutine write_netcdf_parallel + + +end module module_write_netcdf_parallel From dab1dd89212bf8bf6719043073f4c1884848860e Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sun, 26 Jan 2020 20:38:43 +0000 Subject: [PATCH 22/32] allow chunksizes for 2d arrays to be set in model_configure (ichunk2d,jchunk2d) Default is size of array on each write task. --- fv3_cap.F90 | 7 ++++++- io/module_fv3_io_def.F90 | 1 + io/module_write_netcdf_parallel.F90 | 9 +++++---- io/module_wrt_grid_comp.F90 | 17 ++++++++++++++--- 4 files changed, 26 insertions(+), 8 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index 70594031f..c58b898ac 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -37,7 +37,7 @@ 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,& write_fsyncflag, nsout_io, & cen_lon, cen_lat, ideflate, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -307,6 +307,11 @@ 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) + + ! zlib compression flag call ESMF_ConfigGetAttribute(config=CF,value=ideflate,default=0,label ='ideflate:',rc=rc) if (ideflate < 0) ideflate=0 diff --git a/io/module_fv3_io_def.F90 b/io/module_fv3_io_def.F90 index 55249b14c..2b326f8b8 100644 --- a/io/module_fv3_io_def.F90 +++ b/io/module_fv3_io_def.F90 @@ -18,6 +18,7 @@ module module_fv3_io_def character(255) :: app_domain character(255) :: output_grid integer :: imo,jmo + integer :: ichunk2d,jchunk2d integer :: nbdlphys integer :: nsout_io, iau_offset, ideflate, nbits real :: cen_lon, cen_lat, lon1, lat1, lon2, lat2, dlon, dlat diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 5fe52b0bf..584521e52 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -19,14 +19,14 @@ module module_write_netcdf_parallel contains !---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) + subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d, jchunk2d, rc) ! type(ESMF_FieldBundle), intent(in) :: fieldbundle type(ESMF_FieldBundle), intent(in) :: wrtfb character(*), intent(in) :: filename integer, intent(in) :: mpi_comm integer, intent(in) :: mype - integer, intent(in) :: im, jm + integer, intent(in) :: im, jm, ichunk2d, jchunk2d integer, optional,intent(out) :: rc ! !** local vars @@ -149,7 +149,8 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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); NC_ERR_STOP(ncerr) + shuffle=.false.,deflate_level=ideflate,& + chunksizes=(/ichunk2d,jchunk2d,1/)); NC_ERR_STOP(ncerr) ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else @@ -171,7 +172,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, & (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), & shuffle=.false.,deflate_level=ideflate,& - chunksizes=(/im,jm,1,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr) + chunksizes=(/im,jm,1,1/)); NC_ERR_STOP(ncerr) ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index c70b55d28..e55b70bb4 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -34,7 +34,7 @@ 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,& nsout => nsout_io, & cen_lon, cen_lat, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -1444,10 +1444,21 @@ 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 #endif - + ! set default chunksizes for 2d arrays + if (ichunk2d <= 0) then + 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 + 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 (wrt_int_state%mype == 0) then + print *,'ichunk2d,jchunk2d',ichunk2d,jchunk2d + 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,rc) + trim(filename), wrt_mpi_comm,wrt_int_state%mype,imo,jmo,ichunk2d,jchunk2d,rc) wend = MPI_Wtime() if (lprnt) then write(*,'(A,F10.5,A,I4.2,A,I2.2)')' parallel netcdf Write Time is ',wend-wbeg & From a6c36fbee098ab2ce6e479f55a2e0ba126d7f374 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 02:41:14 +0000 Subject: [PATCH 23/32] add ichunk3d,jchunk3d,kchunk3d to specify 3d chunksizes. Default is now ichunk3d,jchunk3d same as array size on each PE, kchunk3d=nlevs This results in the fastest writes on hera. --- fv3_cap.F90 | 4 ++++ io/module_fv3_io_def.F90 | 2 +- io/module_write_netcdf_parallel.F90 | 7 ++++--- io/module_wrt_grid_comp.F90 | 24 +++++++++++++++++++++--- 4 files changed, 30 insertions(+), 7 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index c58b898ac..7d2ba2c35 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -38,6 +38,7 @@ module fv3gfs_cap_mod lead_wrttask, last_wrttask, & output_grid, output_file, & imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,& + ichunk3d,jchunk3d,kchunk3d, & write_fsyncflag, nsout_io, & cen_lon, cen_lat, ideflate, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -310,6 +311,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! 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) diff --git a/io/module_fv3_io_def.F90 b/io/module_fv3_io_def.F90 index 2b326f8b8..b31ab666b 100644 --- a/io/module_fv3_io_def.F90 +++ b/io/module_fv3_io_def.F90 @@ -18,7 +18,7 @@ module module_fv3_io_def character(255) :: app_domain character(255) :: output_grid integer :: imo,jmo - integer :: ichunk2d,jchunk2d + 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 diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 584521e52..8d4dc9b7e 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -19,14 +19,15 @@ module module_write_netcdf_parallel contains !---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d, jchunk2d, rc) + 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 + integer, intent(in) :: im, jm, ichunk2d, jchunk2d, & + ichunk3d, jchunk3d, kchunk3d integer, optional,intent(out) :: rc ! !** local vars @@ -172,7 +173,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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) + chunksizes=(/ichunk3d,jchunk3d,kchunk3d,1/)); NC_ERR_STOP(ncerr) ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index e55b70bb4..9068a2ac0 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -35,6 +35,7 @@ module module_wrt_grid_comp n_group, num_files, app_domain, & filename_base, output_grid, output_file, & imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,& + ichunk3d,jchunk3d,kchunk3d, & nsout => nsout_io, & cen_lon, cen_lat, & lon1, lat1, lon2, lat2, dlon, dlat, & @@ -1171,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) @@ -1444,7 +1445,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 #endif - ! set default chunksizes for 2d arrays + ! set default chunksizes if (ichunk2d <= 0) then ichunk2d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 call mpi_bcast(ichunk2d,1,mpi_integer,0,wrt_mpi_comm,rc) @@ -1453,12 +1454,29 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) 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 + 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 + 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 + 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) + 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 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,rc) + 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 & From 1cf5a6c9c5253853d801eab5bf7775354aa186e2 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 17:25:51 +0000 Subject: [PATCH 24/32] fix typo --- io/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index bf3a292b4..f25645128 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -6,9 +6,9 @@ else() endif() if(PARALLEL_NETCDF) - set(NETCDF_PARALLEL_SRC model_write_netcdf_parallel.F90) + set(NETCDF_PARALLEL_SRC module_write_netcdf_parallel.F90) else() - set(NETCDF_PARALLEL_SRC model_write_netcdf_parallel_stub.F90) + set(NETCDF_PARALLEL_SRC module_write_netcdf_parallel_stub.F90) add_definitions(-DNO_PARALLEL_NETCDF) endif() From a3f03061e4a4a3544e653b42a64b61242acdea08 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 17:33:24 +0000 Subject: [PATCH 25/32] put ifdefs in module_write_netcdf_parallel.F90 so no stub file needed --- io/CMakeLists.txt | 7 ++----- io/makefile | 10 +--------- io/module_write_netcdf_parallel.F90 | 15 +++++++++++++++ 3 files changed, 18 insertions(+), 14 deletions(-) diff --git a/io/CMakeLists.txt b/io/CMakeLists.txt index f25645128..e2cedc43d 100644 --- a/io/CMakeLists.txt +++ b/io/CMakeLists.txt @@ -5,10 +5,7 @@ else() add_definitions(-DNO_INLINE_POST) endif() -if(PARALLEL_NETCDF) - set(NETCDF_PARALLEL_SRC module_write_netcdf_parallel.F90) -else() - set(NETCDF_PARALLEL_SRC module_write_netcdf_parallel_stub.F90) +if(NOT PARALLEL_NETCDF) add_definitions(-DNO_PARALLEL_NETCDF) endif() @@ -19,7 +16,7 @@ add_library( FV3GFS_io.F90 module_write_nemsio.F90 module_write_netcdf.F90 - ${NETCDF_PARALLEL_SRC} + 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 adcb6cc38..e204f1de8 100644 --- a/io/makefile +++ b/io/makefile @@ -29,14 +29,6 @@ POST_SRC = \ ./post_nems_routines.F90 endif -ifneq (,$(findstring NO_PARALLEL_NETCDF,$(CPPDEFS))) -NETCDF_PARALLEL_SRC = \ - ./module_write_netcdf_parallel_stub.F90 -else -NETCDF_PARALLEL_SRC = \ - ./module_write_netcdf_parallel.F90 -endif - SRCS_f = SRCS_f90 = @@ -49,7 +41,7 @@ SRCS_F90 = \ $(POST_SRC) \ ./module_write_nemsio.F90 \ ./module_write_netcdf.F90 \ - $(NETCDF_PARALLEL_SRC) \ + ./module_write_netcdf_parallel.F90 \ ./module_fv3_io_def.F90 \ ./module_write_internal_state.F90 \ ./module_wrt_grid_comp.F90 diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 8d4dc9b7e..335d0339d 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -18,6 +18,20 @@ module module_write_netcdf_parallel contains +#ifdef NO_PARALLEL_NETCDF +!---------------------------------------------------------------------------------------- + subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) + type(ESMF_FieldBundle), intent(in) :: fieldbundle + type(ESMF_FieldBundle), intent(in) :: wrtfb + character(*), intent(in) :: filename + integer, intent(in) :: mpi_comm + integer, intent(in) :: mype + integer, intent(in) :: im, jm + integer, optional,intent(out) :: rc + print *,'in stub write_netcdf_parallel - model not built with parallel netcdf support, return' + end subroutine write_netcdf_parallel +end module module_write_netcdf_parallel +#else !---------------------------------------------------------------------------------------- subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d, jchunk2d, ichunk3d, jchunk3d, kchunk3d, rc) ! @@ -378,6 +392,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i !print *,'netcdf parallel close, finished write_netcdf_parallel' end subroutine write_netcdf_parallel +#endif !---------------------------------------------------------------------------------------- subroutine get_global_attr(fldbundle, ncid, rc) From 1525cc05e6ca414513337f3706c4c47c71118218 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 17:34:06 +0000 Subject: [PATCH 26/32] don't need this file anymore --- io/module_write_netcdf_parallel_stub.F90 | 26 ------------------------ 1 file changed, 26 deletions(-) delete mode 100644 io/module_write_netcdf_parallel_stub.F90 diff --git a/io/module_write_netcdf_parallel_stub.F90 b/io/module_write_netcdf_parallel_stub.F90 deleted file mode 100644 index 3fe32e063..000000000 --- a/io/module_write_netcdf_parallel_stub.F90 +++ /dev/null @@ -1,26 +0,0 @@ -module module_write_netcdf_parallel - - use esmf - implicit none - private - public write_netcdf_parallel - - contains - -!---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) -! - type(ESMF_FieldBundle), intent(in) :: fieldbundle - type(ESMF_FieldBundle), intent(in) :: wrtfb - character(*), intent(in) :: filename - integer, intent(in) :: mpi_comm - integer, intent(in) :: mype - integer, intent(in) :: im, jm - integer, optional,intent(out) :: rc - - print *,'in stub write_netcdf_parallel - model not built with parallel netcdf support, return' - - end subroutine write_netcdf_parallel - - -end module module_write_netcdf_parallel From 51d3dc1d7dc5e0814fa4ad454f1c48fe218c0ace Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 17:42:44 +0000 Subject: [PATCH 27/32] remove module_write_netcdf_parallel_stub.o target --- io/makefile | 2 -- 1 file changed, 2 deletions(-) diff --git a/io/makefile b/io/makefile index e204f1de8..b6a3946b0 100644 --- a/io/makefile +++ b/io/makefile @@ -73,8 +73,6 @@ module_write_netcdf.o: 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_netcdf_parallel_stub.o: module_write_netcdf_parallel_stub.F90 - $(FC) $(CPPDEFS) $(CPPFLAGS) $(FPPFLAGS) $(FFLAGS) $(OTHERFLAGS) $(OTHER_FFLAGS) $(ESMF_INC) -c module_write_netcdf_parallel_stub.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 From 8ddde8975396e1025d1e621fad2457ef5cd78576 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 18:31:16 +0000 Subject: [PATCH 28/32] use specified chunksizes for serial IO. If chunksize parameter negative, let netcdf library choose defaults. --- io/module_write_netcdf.F90 | 53 +++++++++++++++------- io/module_write_netcdf_parallel.F90 | 38 ++++++++++++---- io/module_wrt_grid_comp.F90 | 68 ++++++++++++++++------------- 3 files changed, 106 insertions(+), 53 deletions(-) diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 33bf17a7a..270e36ec6 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -18,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 @@ -26,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 @@ -65,6 +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) @@ -146,28 +148,49 @@ 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/),cache_size=40*im*jm); 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/),cache_size=40*im*jm); 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) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 335d0339d..55dd3c19b 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -82,6 +82,7 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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) @@ -162,10 +163,18 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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=.false.,deflate_level=ideflate,& - chunksizes=(/ichunk2d,jchunk2d,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/)); NC_ERR_STOP(ncerr) + endif ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else @@ -184,10 +193,23 @@ subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, i 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=(/ichunk3d,jchunk3d,kchunk3d,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 ! compression filters require collective access. ncerr = nf90_var_par_access(ncid, varids(i), NF90_COLLECTIVE) else diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index 9068a2ac0..fb5dc3d7b 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1365,6 +1365,40 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) file_bundle = wrt_int_state%wrtFB(nbdl) endif + ! set default chunksizes for netcdf output + ! (use 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 + 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 + 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 + 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 + 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 + 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) + 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 @@ -1430,7 +1464,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) 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 & @@ -1445,34 +1480,6 @@ 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 #endif - ! set default chunksizes - if (ichunk2d <= 0) then - 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 - 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 - 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 - 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 - 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) - 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 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,& @@ -1534,7 +1541,8 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) 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 & From 98948bddf6571abf61e4966c803d7b9f8371a4af Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Mon, 27 Jan 2020 18:43:09 +0000 Subject: [PATCH 29/32] update comments --- fv3_cap.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index 7d2ba2c35..72997d32d 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -360,7 +360,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ierr = 0 do i=1,num_files CALL ESMF_ConfigGetAttribute(config=CF,value=output_file(i), rc=rc) - ! if only 1 output_file specified, copy to all elements + ! if only 1 output_file specified, copy to all elements, set ierr=1 if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then ierr = 1 if (i .eq. 1) then @@ -370,9 +370,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif endif ! specifying multiple values of output_file will only work if - ! they all start with 'netcdf' - need to add a check for this + ! they all start with 'netcdf'. ! Mixing netcdf and nemsio will not work. if (i > 1 .and. ierr == 0) then + ! if multiple values of output_file specified, check that + ! they all start with 'netcdf'. if (output_file(i)(1:6) .ne. 'netcdf' .or. & output_file(i-1)(1:6) .ne. 'netcdf') then write(0,*)'fv3_cap.F90:muliple values of output_file must all begin with netcdf' From 45d35d5836cad0c29809e22aebabfd5d8c9d9a83 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Thu, 30 Jan 2020 16:51:14 +0000 Subject: [PATCH 30/32] get output_file without esmf call error --- fv3_cap.F90 | 46 ++++++++++++++++++++-------------------------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index 72997d32d..df8779e54 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -238,6 +238,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) integer :: mpi_comm_atm integer :: i, j, k, io_unit, urc, ierr integer :: petcount, mype + integer :: num_output_file logical :: isPetLocal logical :: OPENED character(ESMF_MAXSTR) :: name @@ -355,33 +356,26 @@ 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 + allocate(output_file(num_files)) - CALL ESMF_ConfigFindLabel(CF,'output_file:',rc=RC) - ierr = 0 - do i=1,num_files - CALL ESMF_ConfigGetAttribute(config=CF,value=output_file(i), rc=rc) - ! if only 1 output_file specified, copy to all elements, set ierr=1 - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then - ierr = 1 - if (i .eq. 1) then - return - else - output_file(i) = output_file(i-1) - endif - endif - ! specifying multiple values of output_file will only work if - ! they all start with 'netcdf'. - ! Mixing netcdf and nemsio will not work. - if (i > 1 .and. ierr == 0) then - ! if multiple values of output_file specified, check that - ! they all start with 'netcdf'. - if (output_file(i)(1:6) .ne. 'netcdf' .or. & - output_file(i-1)(1:6) .ne. 'netcdf') then - write(0,*)'fv3_cap.F90:muliple values of output_file must all begin with netcdf' - call ESMF_Finalize(endflag=ESMF_END_ABORT) - endif - endif - enddo + 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' + ndif if(mype == 0) then print *,'af nems config,num_files=',num_files do i=1,num_files From a0f9ffa7063893dfaa01fda59649aabc4c326c25 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Thu, 30 Jan 2020 20:46:12 +0000 Subject: [PATCH 31/32] syntax fix --- fv3_cap.F90 | 2 +- io/module_wrt_grid_comp.F90 | 20 +++++++++++++------- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/fv3_cap.F90 b/fv3_cap.F90 index df8779e54..c211a9044 100644 --- a/fv3_cap.F90 +++ b/fv3_cap.F90 @@ -375,7 +375,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) output_file(1:num_files) = output_file(1) else output_file(1:num_files) = 'netcdf' - ndif + endif if(mype == 0) then print *,'af nems config,num_files=',num_files do i=1,num_files diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index fb5dc3d7b..a7177fe22 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -1366,31 +1366,37 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc) endif ! set default chunksizes for netcdf output - ! (use use MPI decomposition size). + ! (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 - ichunk2d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 + 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 - jchunk2d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 + 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 - ichunk3d = wrt_int_state%lon_end-wrt_int_state%lon_start+1 + 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 - jchunk3d = wrt_int_state%lat_end-wrt_int_state%lat_start+1 + 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 - call ESMF_FieldBundleGet(wrt_int_state%wrtFB(nbdl), grid=wrtgrid) - call ESMF_AttributeGet(wrtgrid, convention="NetCDF", purpose="FV3", & + 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 From 715f1aeb9de2a8ecce53135138c37c8f58cc4e73 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Sat, 1 Feb 2020 16:49:08 +0000 Subject: [PATCH 32/32] fix stub interface in module_write_netcdf_parallel.F90 for cmake build --- io/module_write_netcdf_parallel.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/io/module_write_netcdf_parallel.F90 b/io/module_write_netcdf_parallel.F90 index 55dd3c19b..7d7dbb0f7 100644 --- a/io/module_write_netcdf_parallel.F90 +++ b/io/module_write_netcdf_parallel.F90 @@ -20,17 +20,17 @@ module module_write_netcdf_parallel #ifdef NO_PARALLEL_NETCDF !---------------------------------------------------------------------------------------- - subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc) + 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 + 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 -end module module_write_netcdf_parallel #else !---------------------------------------------------------------------------------------- subroutine write_netcdf_parallel(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, ichunk2d, jchunk2d, ichunk3d, jchunk3d, kchunk3d, rc)