parallel netcdf IO (#41)

* netcdf parallel writing and lon/lat in netcdf file

* some changes to get it to (almost) run on hera

* lat should be in degrees - remove conversion to radians.

* import updates from jswhit/fv3atm

* more bug fixes - now works on hera

* specify collective access if compression turned on.

* use classic model

* bug fixes for parallel IO with compression

* fix calculation of max compression error

* turn off shuffle filter

* code simplification

* remove debug print

* don't use parallel IO for 2d file (since it seems to increase run time)

* allow multiple values of output_file, as long as they all start with 'netcdf'

* use default chunksize for 2d vars

* delete commented out macro ESMF_ERR_ABORT

* delete rad2dg

* add module_write_netcdf_parallel.F90

* add option to build without parallel netcdf (-DNO_PARALLEL_NETCDF)

* fix typo

* stub file for building without parallel netcdf lib

* allow chunksizes for 2d arrays to be set in model_configure (ichunk2d,jchunk2d)
Default is size of array on each write task.

* 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.

* fix typo

* put ifdefs in module_write_netcdf_parallel.F90 so no stub file needed

* don't need this file anymore

* remove module_write_netcdf_parallel_stub.o target

* use specified chunksizes for serial IO. If chunksize parameter negative, let netcdf library choose defaults.

* update comments

* get output_file without esmf call error

* syntax fix

* fix stub interface in module_write_netcdf_parallel.F90 for cmake build

Co-authored-by: junwang-noaa <37633869+junwang-noaa@users.noreply.github.com>

fv3_cap.F90

@@ -37,7 +37,8 @@ module fv3gfs_cap_mod
                                     wrttasks_per_group, n_group,             &
                                     lead_wrttask, last_wrttask,              &
                                     output_grid, output_file,                &
-                                    imo, jmo, write_nemsioflip,              &
+                                    imo,jmo,ichunk2d,jchunk2d,write_nemsioflip,&
+                                    ichunk3d,jchunk3d,kchunk3d,              &
                                     write_fsyncflag, nsout_io,               &
                                     cen_lon, cen_lat, ideflate,              &
                                     lon1, lat1, lon2, lat2, dlon, dlat,      &
@@ -235,8 +236,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
 
     integer,dimension(6)                   :: date, date_init
     integer                                :: mpi_comm_atm
-    integer                                :: i, j, k, io_unit, urc
+    integer                                :: i, j, k, io_unit, urc, ierr
     integer                                :: petcount, mype
+    integer                                :: num_output_file
     logical                                :: isPetLocal
     logical                                :: OPENED
     character(ESMF_MAXSTR)                 :: name
@@ -307,6 +309,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     call ESMF_ConfigGetAttribute(config=CF,value=iau_offset,default=0,label ='iau_offset:',rc=rc)
     if (iau_offset < 0) iau_offset=0
 
+    ! chunksizes for netcdf_parallel
+    call ESMF_ConfigGetAttribute(config=CF,value=ichunk2d,default=0,label ='ichunk2d:',rc=rc)
+    call ESMF_ConfigGetAttribute(config=CF,value=jchunk2d,default=0,label ='jchunk2d:',rc=rc)
+    call ESMF_ConfigGetAttribute(config=CF,value=ichunk3d,default=0,label ='ichunk3d:',rc=rc)
+    call ESMF_ConfigGetAttribute(config=CF,value=jchunk3d,default=0,label ='jchunk3d:',rc=rc)
+    call ESMF_ConfigGetAttribute(config=CF,value=kchunk3d,default=0,label ='kchunk3d:',rc=rc)
+
+    ! zlib compression flag
     call ESMF_ConfigGetAttribute(config=CF,value=ideflate,default=0,label ='ideflate:',rc=rc)
     if (ideflate < 0) ideflate=0
 
@@ -346,8 +356,33 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
         CALL ESMF_ConfigGetAttribute(config=CF,value=filename_base(i), rc=rc)
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
       enddo
-      if(mype == 0) print *,'af nems config,num_files=',num_files, &
-                            'filename_base=',filename_base
+
+      allocate(output_file(num_files))
+      num_output_file = ESMF_ConfigGetLen(config=CF, label ='output_file:',rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+      if (num_files == num_output_file) then
+        CALL ESMF_ConfigGetAttribute(CF,valueList=output_file,label='output_file:', &
+             count=num_files, rc=RC)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+        do i = 1, num_files
+          if(output_file(i) /= "netcdf" .and. output_file(i) /= "netcdf_parallel") then
+            write(0,*)"fv3_cap.F90: only netcdf and netcdf_parallel are allowed for multiple values of output_file"
+            call ESMF_Finalize(endflag=ESMF_END_ABORT)
+          endif
+        enddo
+      else if ( num_output_file == 1) then
+        CALL ESMF_ConfigGetAttribute(CF,valuelist=output_file,label='output_file:', count=1, rc=RC)
+        output_file(1:num_files) = output_file(1)
+      else
+        output_file(1:num_files) = 'netcdf'
+      endif
+      if(mype == 0) then
+        print *,'af nems config,num_files=',num_files
+        do i=1,num_files
+           print *,'num_file=',i,'filename_base= ',trim(filename_base(i)),&
+           ' output_file= ',trim(output_file(i))
+        enddo
+      endif
 !
 ! variables for alarms
       call ESMF_ConfigGetAttribute(config=CF, value=nfhout,   label ='nfhout:',   rc=rc)
@@ -359,10 +394,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
 
 ! variables for I/O options
       call ESMF_ConfigGetAttribute(config=CF, value=output_grid, label ='output_grid:',rc=rc)
-      call ESMF_ConfigGetAttribute(config=CF, value=output_file, label ='output_file:',rc=rc)
       if (mype == 0) then
         print *,'output_grid=',trim(output_grid)
-        print *,'output_file=',trim(output_file)
       end if
       write_nemsioflip =.false.
       write_fsyncflag  =.false.

io/CMakeLists.txt

@@ -5,13 +5,18 @@ else()
   add_definitions(-DNO_INLINE_POST)
 endif()
 
+if(NOT PARALLEL_NETCDF)
+  add_definitions(-DNO_PARALLEL_NETCDF)
+endif()
+
 add_library(
     io
 
     ffsync.F90
     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

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

io/module_fv3_io_def.F90

@@ -17,13 +17,14 @@ module module_fv3_io_def
   integer           :: num_files
   character(255)    :: app_domain
   character(255)    :: output_grid
-  character(255)    :: output_file
   integer           :: imo,jmo
+  integer           :: ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d
   integer           :: nbdlphys
   integer           :: nsout_io, iau_offset, ideflate, nbits
   real              :: cen_lon, cen_lat, lon1, lat1, lon2, lat2, dlon, dlat
   real              :: stdlat1, stdlat2, dx, dy
   character(255),dimension(:),allocatable :: filename_base
+  character(255),dimension(:),allocatable :: output_file
 !
   integer,dimension(:),allocatable     :: lead_wrttask, last_wrttask
 !

io/module_write_netcdf.F90

@@ -1,7 +1,3 @@
-!#define ESMF_ERR_ABORT(rc) \ 
-!if (rc /= ESMF_SUCCESS) write(0,*) 'rc=',rc,__FILE__,__LINE__; \
-!   if (ESMF_LogFoundError(rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
-
 #define ESMF_ERR_RETURN(rc) if (ESMF_LogFoundError(rc, msg="Breaking out of subroutine", line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
 
 #define NC_ERR_STOP(status) \
@@ -22,14 +18,15 @@ 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
     character(*), intent(in)           :: filename
     integer, intent(in)                :: mpi_comm
     integer, intent(in)                :: mype
     integer, intent(in)                :: im, jm
+    integer, intent(in)                :: ichunk2d,jchunk2d,ichunk3d,jchunk3d,kchunk3d
     integer, optional,intent(out)      :: rc
 !
 !** local vars
@@ -69,9 +66,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     integer :: im_dimid, jm_dimid, pfull_dimid, phalf_dimid, time_dimid
     integer :: im_varid, jm_varid, lm_varid, time_varid, lon_varid, lat_varid
     integer, dimension(:), allocatable :: varids
-!
-!!
-!
+    logical shuffle
 
     call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc)
 
@@ -117,15 +112,10 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 ! create netcdf file and enter define mode
     if (mype==0) then
 
-    if (ideflate == 0) then
-        ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),NF90_SHARE), &
-        ncid=ncid); NC_ERR_STOP(ncerr)
-        ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
-    else
-        ncerr = nf90_create(trim(filename), cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL), &
-        ncid=ncid); NC_ERR_STOP(ncerr)
-        ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
-    endif
+    ncerr = nf90_create(trim(filename),&
+            cmode=IOR(IOR(NF90_CLOBBER,NF90_NETCDF4),NF90_CLASSIC_MODEL),&
+            ncid=ncid); NC_ERR_STOP(ncerr)
+    ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
 
     ! define dimensions
     ncerr = nf90_def_dim(ncid, "grid_xt", im, im_dimid); NC_ERR_STOP(ncerr)
@@ -158,32 +148,53 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
       if (fldlev(i) == 1) then
         if (typekind == ESMF_TYPEKIND_R4) then
           if (ideflate > 0) then
-            ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                                (/im_dimid,jm_dimid,time_dimid/), varids(i), &
-                                 shuffle=.true.,deflate_level=ideflate, &
-                                 chunksizes=(/im,jm,1/)); NC_ERR_STOP(ncerr)
+             if (ichunk2d < 0 .or. jchunk2d < 0) then
+                ! let netcdf lib choose chunksize
+                ! shuffle filter on for 2d fields (lossless compression)
+                ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                        (/im_dimid,jm_dimid,time_dimid/), varids(i), &
+                        shuffle=.true.,deflate_level=ideflate); NC_ERR_STOP(ncerr)
+             else
+                ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                        (/im_dimid,jm_dimid,time_dimid/), varids(i), &
+                        shuffle=.true.,deflate_level=ideflate,&
+                        chunksizes=(/ichunk2d,jchunk2d,1/),cache_size=40*im*jm); NC_ERR_STOP(ncerr)
+             endif
           else
-            ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                                (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
+             ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+             (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
           endif
         else if (typekind == ESMF_TYPEKIND_R8) then
-          ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
-                              (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
+           ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
+                               (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
         else
-          write(0,*)'Unsupported typekind ', typekind
-          stop
+           write(0,*)'Unsupported typekind ', typekind
+           stop
         end if
       else if (fldlev(i) > 1) then
          if (typekind == ESMF_TYPEKIND_R4) then
            if (ideflate > 0) then
-             ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                                 (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), &
-                                 shuffle=.false.,deflate_level=ideflate, &
-                                 chunksizes=(/im,jm,1,1/)); NC_ERR_STOP(ncerr)
+             ! shuffle filter off for 3d fields using lossy compression
+             if (nbits > 0) then
+                shuffle=.false.
+             else
+                shuffle=.true.
+             endif
+             if (ichunk3d < 0 .or. jchunk3d < 0 .or. kchunk3d < 0) then
+                ! let netcdf lib choose chunksize
+                ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                        (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), &
+                        shuffle=shuffle,deflate_level=ideflate); NC_ERR_STOP(ncerr)
+             else
+                ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
+                       (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i), &
+                       shuffle=shuffle,deflate_level=ideflate,&
+                       chunksizes=(/ichunk3d,jchunk3d,kchunk3d,1/)); NC_ERR_STOP(ncerr)
+             endif
            else
              ncerr = nf90_def_var(ncid, trim(fldName), NF90_FLOAT, &
-                                 (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
-           endif
+                     (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
+           endif 
          else if (typekind == ESMF_TYPEKIND_R8) then
            ncerr = nf90_def_var(ncid, trim(fldName), NF90_DOUBLE, &
                                 (/im_dimid,jm_dimid,pfull_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
@@ -224,7 +235,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
            call ESMF_AttributeGet(fcstField(i), convention="NetCDF", purpose="FV3", &
                                   name=trim(attName), value=varr8val, &
                                   rc=rc); ESMF_ERR_RETURN(rc)
-           if (trim(attName) /= '_FillValue' ) then
+           if (trim(attName) /= '_FillValue') then
               ! FIXME:  _FillValue must be cast to var type for recent versions of netcdf
               ncerr = nf90_put_att(ncid, varids(i), trim(attName), varr8val); NC_ERR_STOP(ncerr)
            endif
@@ -241,23 +252,23 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 
     end do   ! i=1,fieldCount
 
-! write grid_xt, grid_yt attributes
+    ! write grid_xt, grid_yt attributes
     if (trim(output_grid) == 'gaussian_grid' .or. &
         trim(output_grid) == 'regional_latlon') then
-      ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
     else if (trim(output_grid) == 'rotated_latlon') then
-      ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "long_name", "rotated T-cell longiitude"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "long_name", "rotated T-cell latiitude"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees"); NC_ERR_STOP(ncerr)
     else if (trim(output_grid) == 'lambert_conformal') then
-      ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
-      ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "long_name", "x-coordinate of projection"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "long_name", "y-coordinate of projection"); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
     endif
 
     ncerr = nf90_enddef(ncid); NC_ERR_STOP(ncerr)
@@ -377,6 +388,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
 
     deallocate(fcstField)
     deallocate(varids)
+    deallocate(compress_err)
 
     if (mype==0) then
     ncerr = nf90_close(ncid=ncid); NC_ERR_STOP(ncerr)