begin 
    ;---- change the time attribute after applying ncrcat operator and subsetting with ncea operator
    ;---- need to do this to get the correct time sequence

      fi   = addfile ("./pr_1x1.nc", "r") 

      file_out = "pr_1x1_test.nc"
      system ("rm " + file_out)  ; remove any previously exist file
    ;---- read variables
      lat = fi->lat		
      lon = fi->lon
      time = fi->time
    ;---- get dimension sizes
      ntim = dimsizes(time)	; 24000 values; 1901-01 to 2100-12
      nlat = dimsizes(lat)	; 180 values
      nlon = dimsizes(lon)	; 360

    ;---- open output file and explicitly declare file definition mode.
      fout = addfile(file_out, "c")
      setfileoption(fout,"DefineMode",True)
     
    ;---- predefine the coordinate variables and their dimensionality
      dimNames = (/"time", "lat", "lon"/)  
      dimSizes = (/ -1   , nlat,   nlon /) 
      dimUnlim = (/ True , False, False /)   
      filedimdef(fout,dimNames,dimSizes,dimUnlim)

    ;---- create global attributes of the file
      filevardef( fout, "lat", "double", "lat" )
      latAtt=0
      latAtt@units     = "degrees_north"
      latAtt@axis      = "Y"
      latAtt@bounds    = "lat_bnds"
      latAtt@long_name = "latitude"
      filevarattdef( fout, "lat", latAtt )

      filevardef( fout, "lon", "double", "lon" )
      lonAtt=0
      lonAtt@units     = "degrees_east"
      lonAtt@axis      = "X"
      lonAtt@bounds    = "lon_bnds"
      lonAtt@long_name = "longitude"
      filevarattdef( fout, "lon", lonAtt )

      filevardef( fout, "time", "double", "time" )
      timeAtt=0
      timeAtt@units     = "days since 1901-01-01 00:00:00"
      timeAtt@calendar  = "noleap"
      timeAtt@axis      = "T"
      timeAtt@bounds    = "time_bnds"
      timeAtt@long_name = "time"
      filevarattdef( fout, "time", timeAtt )

      filevardef( fout, "pr", "float", (/"time","lat","lon"/) )
      prAtt=0
      prAtt@units     = "kg m-2 s-1"
      prAtt@missing_value   = 1e+20
      prAtt@FillValue      = 1e+20
      prAtt@long_name = "Precipitation"
      filevarattdef( fout, "pr", prAtt )
   
   ;---- exit file definition mode.
      setfileoption(fout,"DefineMode",False)
   ;---- output 
       fout->time   = (/time/)     
       fout->lat    = (/lat/)
       fout->lon    = (/lon/) 
       fout->pr      = (/pr/)

    end