GFSv16 netcdf post ficein cpl (ufs-community#48)

* fv3atm issue ufs-community#37: fix the real(8) lat/lon in netcdf file
* fv3atm ufs-community#35: Reducing background vertical diffusivities in the inversion layers
* fv3atm ufs-community#24: bug in gfsphysics/physics/moninedmf_hafs.f
* fv3atm ufs-community#18: Optimize netcdf write component and bugfix for post and samfdeepcnv.f
* set (0-1) bounds for ficein_cpl
* remove cache_size due to lower netcdf verion 4.5.1 on mars
* Change ice falling to 0.9 in gfsphysics/physics/gfdl_cloud_microphys.F90

Author: DusanJovic-NOAA

.gitmodules

@@ -5,6 +5,8 @@
 [submodule "ccpp/framework"]
 	path = ccpp/framework
 	url = https://github.com/NCAR/ccpp-framework
+	branch = master
 [submodule "ccpp/physics"]
 	path = ccpp/physics
 	url = https://github.com/NCAR/ccpp-physics
+	branch = master

atmos_model.F90

@@ -1720,7 +1720,7 @@ subroutine assign_importdata(rc)
                   IPD_Data(nb)%Coupling%ficein_cpl(ix)   = zero
                   if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then
                     if (datar8(i,j) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then
-                      IPD_Data(nb)%Coupling%ficein_cpl(ix) = datar8(i,j)
+                      IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(datar8(i,j),one))
 !                     if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) == one) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points
                       IPD_Data(nb)%Sfcprop%slmsk(ix)         = 2.                                        !slmsk=2 crashes in gcycle on partial land points
                       IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4.

ccpp/framework

@@ -3266,7 +3266,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
             else
                 tc (k) = tk (k) - tice
                 vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee
-                vti (k) = vi0 * exp (log_10 * vti (k)) * 0.8 
+                vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9
                 vti (k) = min (vi_max, max (vf_min, vti (k)))
             endif
         enddo

ccpp/physics

@@ -1360,7 +1360,11 @@ subroutine moninedmf_hafs(ix,im,km,ntrac,ntcw,dv,du,tau,rtg,      &
           tem = 0.5 * (diss(i,k-1)+diss(i,k))
           tem  = max(tem, 0.)
           ttend = tem / cp
-          tau(i,k) = tau(i,k) + 0.5*ttend
+          if (alpha .gt. 0.0) then
+            tau(i,k) = tau(i,k) + 0.5*ttend
+          else
+            tau(i,k) = tau(i,k) + 0.7*ttend  ! in HWRF/HMON, use 0.7
+          endif
         enddo
       enddo
 !

gfsphysics/physics/gfdl_cloud_microphys.F90

@@ -1547,22 +1547,22 @@ subroutine samfdeepcnv(im,ix,km,delt,itc,ntc,ntk,ntr,delp,
       enddo
       enddo
       do i = 1, im
-        betamn = betas
-        if(islimsk(i) == 1) betamn = betal
-        if(ntk > 0) then
-          betamx = betamn + dbeta
-          if(tkemean(i) > tkemx) then
-            beta = betamn
-          else if(tkemean(i) < tkemn) then
-            beta = betamx
+        if(cnvflg(i)) then
+          betamn = betas
+          if(islimsk(i) == 1) betamn = betal
+          if(ntk > 0) then
+            betamx = betamn + dbeta
+            if(tkemean(i) > tkemx) then
+              beta = betamn
+            else if(tkemean(i) < tkemn) then
+              beta = betamx
+            else
+              tem = (betamx - betamn) * (tkemean(i) - tkemn)
+              beta = betamx - tem  / dtke
+            endif
           else
-            tem = (betamx - betamn) * (tkemean(i) - tkemn)
-            beta = betamx - tem  / dtke
+            beta = betamn
           endif
-        else
-          beta = betamn
-        endif
-        if(cnvflg(i)) then
           dz  = (sumx(i)+zi(i,1))/float(kbcon(i))
           tem = 1./float(kbcon(i))
           xlamd(i) = (1.-beta**tem)/dz

gfsphysics/physics/module_sf_noahmp_glacier.f90

@@ -148,7 +148,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
      &                     epsi,    beta,   chx,    cqx,
      &                     rdt,     rdz,    qmin,   qlmin,
      &                     rimin,   rbcr,   rbint,  tdzmin,
-     &                     rlmn,    rlmn1,  rlmx,   elmx,
+     &                     rlmn,    rlmn1,  rlmn2,  
+     &                     rlmx,    elmx,
      &                     ttend,   utend,  vtend,  qtend,
      &                     zfac,    zfmin,  vk,     spdk2,
      &                     tkmin,   tkminx, xkzinv, xkgdx,
@@ -172,7 +173,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
       parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
       parameter(vk=0.4,rimin=-100.)
       parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
-      parameter(rlmn=30.,rlmn1=5.,rlmx=300.,elmx=300.)
+      parameter(rlmn=30.,rlmn1=5.,rlmn2=10.)
+      parameter(rlmx=300.,elmx=300.)
       parameter(prmin=0.25,prmax=4.0)
       parameter(pr0=1.0,prtke=1.0,prscu=0.67)
       parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
@@ -698,8 +700,9 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
 !         if(tem1 > 1.e-5) then
           tem1 = tvx(i,k+1)-tvx(i,k)
           if(tem1 > 0.) then
-             xkzo(i,k)  = min(xkzo(i,k),xkzinv)
-             xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
+             xkzo(i,k)  = min(xkzo(i,k), xkzinv)
+             xkzmo(i,k) = min(xkzmo(i,k), xkzinv)
+             rlmnz(i,k) = min(rlmnz(i,k), rlmn2)
           endif
         enddo
       enddo

gfsphysics/physics/module_sf_noahmplsm.f90

@@ -285,7 +285,7 @@ subroutine FV3GFS_IPD_checksum (Model, IPD_Data, Atm_block)
        temp2d(i,j,55) = IPD_Data(nb)%Coupling%visbmui(ix)
        temp2d(i,j,56) = IPD_Data(nb)%Coupling%visdfui(ix)
        temp2d(i,j,57) = IPD_Data(nb)%Coupling%sfcdsw(ix)
-       temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcnsw(ix)
+       temp2d(i,j,58) = IPD_Data(nb)%Coupling%sfcnsw(ix)
        temp2d(i,j,59) = IPD_Data(nb)%Coupling%sfcdlw(ix)
        temp2d(i,j,60) = IPD_Data(nb)%Grid%xlon(ix)
        temp2d(i,j,61) = IPD_Data(nb)%Grid%xlat(ix)

gfsphysics/physics/module_wrf_utl.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

gfsphysics/physics/moninedmf_hafs.f

@@ -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
 
@@ -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,10 +72,10 @@ 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))
@@ -117,13 +118,13 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
     if (mype==0) then
 
     if (ideflate == 0) then
-        ncerr = nf90_create(trim(filename), cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
+        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)
-        ! if compression on use HDF5 format with default _FillValue
+        ncerr = nf90_set_fill(ncid, NF90_NOFILL, oldMode); NC_ERR_STOP(ncerr)
     endif
 
     ! define dimensions
@@ -156,28 +157,32 @@ 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/)); 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, &
-                               (/im_dimid,jm_dimid,time_dimid/), varids(i)); NC_ERR_STOP(ncerr)
+                              (/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)
+             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)
+           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, &
@@ -219,8 +224,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
 
@@ -236,6 +241,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,63 +271,39 @@ 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_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(:,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
-       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)
     call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc)
     if (mype==0) then
        if (trim(output_grid) == 'gaussian_grid' .or. &
            trim(output_grid) == 'regional_latlon') then
-          ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*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)
+          ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:)  ); NC_ERR_STOP(ncerr)
        else if (trim(output_grid) == 'rotated_latlon') then
           do j=1,jm
              y(j) = lat1 + (lat2-lat1)/(jm-1) * (j-1)
           enddo
           ncerr = nf90_put_var(ncid, jm_varid, values=y  ); NC_ERR_STOP(ncerr)
-          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=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr)
     endif
 
     do i=1, fieldCount
@@ -344,11 +344,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
@@ -363,6 +359,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)
@@ -484,9 +491,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