;;9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
;;   Lat Cen    3.750   Lon Cen    6.250   Count    1
;;         JAN     FEB     MAR     APR     MAY     JUN     JUL     AUG     SEP     OCT     NOV     DEC
;;MEAN    50.2    75.5   146.0   204.6   244.3   271.9   324.0   250.0   321.7   265.9   118.2    62.7
;;STDV    16.2    26.0    36.3    37.3    27.3    33.2    78.7    76.0    61.8    40.3    28.0    19.4
;;1920    47.6    76.2   107.6   191.4   249.5   270.6   345.8   345.5   286.4   297.7   117.7    56.7
;;1921    55.9    61.1   156.0   262.1   237.2   265.9   247.3   352.0   351.2   271.8   108.4    75.8
;;....
;;2013    46.6    65.6   180.6   208.1   234.4   311.4   258.3   212.3   391.0   296.8   145.8    74.3
;;2014    35.7    65.4   119.1   204.4   232.4   279.5   281.1   123.4   427.7   236.0   147.8    84.7
;;9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
;;
;; This type of processing is all about book keeping!
;;
;;================================================================================
;---Define time period

      yrStrt = 1920
      yrLast = 2014
      yyyymm = yyyymm_time(yrStrt,yrLast,"integer")         ; yyyymm(time) 
      ntim   = dimsizes(yyyymm)                             ; all years & months
      printVarSummary(yyyymm)
      print("---")

;---Create grid

      grdMinLat =  6.250
      grdMaxLat = 43.750
      grdDlat   =  2.5
      nlat      = toint((grdMaxLat-grdMinLat)/grdDlat) + 1
      lat       = fspan(grdMinLat, grdMaxLat, nlat)
      print(lat)
      print("---")

      grdMinLon =-23.750
      grdMaxLon = 11.250
      grdDlon   =  2.5
      mlon      = toint((grdMaxLon-grdMinLon)/grdDlon) + 1
      lon       = fspan(grdMinLon, grdMaxLon, mlon)
      print(lon)
      print("---")

      print("ntim="+ntim+"  nlat="+nlat+"  mlon="+mlon)
      print("---")

      pmsg   = -999.0
      prc    = new( (/ntim,nlat,mlon/),"float",pmsg)     ; all entries initialized to 'pmsg' => _FillValue

;---Create work array

      nmos   = 12
      nyrs   = yrLast-yrStrt+1
      work   = new( (/nyrs,nmos/),"float",pmsg) 

;---Read data

      dir_rain   = "./" 
      file_name  = "Rain_estimate_2.5DegGrid.txt"
      path_name  = dir_rain + file_name
      dstr       = asciiread(path_name, -1, "string")
                                                            ; separator for new grid point
      grdPtStrt  = "9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999"
      indStrt    = ind(dstr.eq.grdPtStrt)                   ; indices 
      nGrdPt     = dimsizes(indStrt)                        ; number of grid point with data
      nGrdPt     = nGrdPt-1                                 ; do not count the last 'grdPtStrt'
      print("nGrdPt="+nGrdPt)
      print("---")
      print(indStrt)
      print("---")

      indSkip    = 5
      prtflg     = True
      
      do n=0,nGrdPt-1
         latPt  = tofloat( str_get_field(dstr(indStrt(n)+1), 3, " ") )
         lonPt  = tofloat( str_get_field(dstr(indStrt(n)+1), 6, " ") )
         if (latPt.ge.grdMinLat .and. latPt.le.grdMaxLat .and. \
             lonPt.ge.grdMinLon .and. lonPt.le.grdMaxLon) then
             nl     = ind(latPt.eq.lat)
             ml     = ind(lonPt.eq.lon)
             nStrt  = indStrt(n)+indSkip
             nLast  = nStrt + nyrs -1
             strBlk = dstr(nStrt:nLast)     ; data in string block     
             if (prtflg) then               ; debug
                 print("latPt="+latPt+"  lonPt="+lonPt)
                 print("nl="+nl+"  ml="+ml)
                 print("---")
                 print(strBlk)
                 print("---")
             end if

             do nmo=0,nmos-1
                work(:,nmo) = tofloat( str_get_field(strBlk, 2+nmo, " ") ) 
             end do

             if (prtflg) then               ; debug
                ;opt = True
                ;opt@fout = "RainEstimate."+sprinti("%0.3i", n)+".txt"
                 write_matrix (work, "12f7.1" , False)
                 print("---")
                 prtflg = False
             end if

             prc(:,nl,ml) = ndtooned(work)
         end if
      end do

; Add meta data

      lat!0     = "lat"
      lon!0     = "lon"
      lat@long_name = "latitude"
      lon@long_name = "longitude"
      lat@units = "degrees_north"
      lon@units = "degrees_east"

; Create time for netCDF

      yyyy      = yyyymm/100
      mm        = yyyymm - yyyy*100
      dd        = conform(yyyy, 1, -1)
      hh        = conform(yyyy, 0, -1)
      mn        = conform(yyyy, 0, -1)
      sc        = conform(yyyy, 0, -1)
      tunits    = "days since 1900-01-01 00:00:00"

      time      = cd_inv_calendar(yyyy,mm,dd,hh,mn,sc,tunits, 0)
      time!0    = "time"
      time&time =  time 
      time@units=  tunits

      yyyymm!0  = "time"
      yyyymm&time := time

      prc!0     = "time"
      prc!1     = "lat"
      prc!2     = "lon"
      prc&time  =  time
      prc&lat   =  lat
      prc&lon   =  lon

      prc@long_name = "precipitation"
      prc@units     = "mm"

      printVarSummary(prc)
      printMinMax(prc,0)
      

;---netCDF

     dirNc = "./"
     filNc = "PRECIP.nc"
     pthNc = dirNc + filNc

     system("/bin/rm -f "+pthNc)       ; remove any pre-existing file
     ncdf = addfile(pthNc, "c")        ; open output netCDF file

    ;===================================================================
    ; create global attributes of the file (optional)
    ;===================================================================
       fAtt               = True            ; assign file attributes
       fAtt@title         = "APPO PRECIP: Ascii to netCDF"
       fAtt@source_file   =  file_name         
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

    ;===================================================================
    ; make time an UNLIMITED dimension; recommended  for most applications
    ;===================================================================
       filedimdef(ncdf,"time",-1,True) 

    ;===================================================================
    ; output variables directly; NCL will call appropriate functions
    ; to write the meta data associated with each variable
    ;===================================================================
        ncdf->PRECIP = prc
        ncdf->YYYYMM = yyyymm