load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  netCDF = True
  PLOT   = True

; Read data into a big 1D string array
  diri  = "./"              ; input directory
  fili  = systemfunc("cd "+diri+" ; ls 2007*.dat" )
  nfili = dimsizes(fili)
  print(fili)

; plot resources that do no change
  res = True
  res@gsnAddCyclic = False 
  
  ;;;;;;;;;;;;;;set for the map
  ;res@mpDataBaseVersion = "HighRes"
  ;res@mpDataSetName     = "Earth..4"
  res@mpGridAndLimbOn = False ; Turn off lat/lon lines
  res@mpOutlineOn     = True   ; Turn on map outlines
  res@mpGeophysicalLineColor = "Black"; Change the outline line color
  res@mpGeophysicalLineThicknessF= 2.; double the thickness of geophysical boundaries
  res@mpOutlineSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/) 
  
;;;;;;;;;;;;;;;ZOOM
  res@mpMinLatF = 30 
  res@mpMaxLatF = 34 
  res@mpMinLonF = 114 
  res@mpMaxLonF = 122 
  res@pmTickMarkDisplayMode = "Always"    
  
;;;;;;;;;;;;;; contours
  res@cnFillOn = True 
  res@cnLinesOn = False     ; Turn off contour lines 
  res@cnLevelSelectionMode = "ExplicitLevels" 
  res@cnLevels = (/1, 2, 6, 10, 15, 20, 30, 40, 50, 70, 90, 110, 130, 150, 200/)
;;;;;;;;;;;;;; label bar
  res@lbBoxLinesOn = False ; Turn off label bar box lines 

  do nf=0,nfili-1
    data := asciiread(diri+fili(nf),-1,"string")
  	
  ; Count the number of fields, just to show it can be done.
    nfields = str_fields_count(data(1)," ")
    print("number of fields = " + nfields)
    
  ;
  ; Skip first row of "data" because it's just a header line. (optional)
  ;
  ; Use a space (" ") as a delimiter in str_get_field. The first
  ; field is field=1 (unlike str_get_cols, in which the first column
  ; is column=0).
  ;
   ;lat  := stringtofloat(str_get_field(data(1:), 2," ")) (optional)
   ;lon  := stringtofloat(str_get_field(data(1:), 3," ")) (optional)
   ;pwv  := stringtofloat(str_get_field(data(1:), 4," ")) (optional)
    
  ; No skipping first coloumn
    lat  := stringtofloat(str_get_field(data(1:), 2," ")) 
    lon  := stringtofloat(str_get_field(data(1:), 1," ")) 
    rain := stringtofloat(str_get_field(data(1:), 3," ")) 
  
    lon@long_name   = "longitude"
    lon@units       = "degrees_east"
    lat@long_name   = "latitude"
    lat@units       = "degrees_north"
  
    rain@_FillValue = -255.0
    rain@long_name  = "Precipitation"
    rain@units      = "mm/???" 
  
    rain!0 = "npts"       ; give all the same dimension name
    lat!0  = "npts"       ; not required but better
    lon!0  = "npts"
  
    printVarSummary(rain)
  
    fname_root = str_get_field(fili(nf),1,".")
  
    if (netCDF) then
  ;***************************************************
  ; Open output netcdf file. Remove first just in case.
  ;***************************************************
       ncdf_out_fname = fname_root + ".nc"
       system("rm -f " + ncdf_out_fname) ; remove any pre-existing file 
       ncdf_out = addfile(ncdf_out_fname,"c")   
  
       fAtt = True
       fAtt@title  = "Precipitation"  
       fAtt@source_file   =  fili(nf)
       fAtt@Conventions   = "None"   
       fAtt@creation_date = systemfunc ("date")        
       fileattdef( ncdf_out, fAtt )   
   
       ncdf_out->lat  =  lat
       ncdf_out->lon  =  lon
       ncdf_out->rain =  rain
    end if
  
    if (PLOT) then
        wks = gsn_open_wks("x11", fname_root)
       ;gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
       ;gsn_define_colormap(wks,"WhViBlGrYeOrRe")
       ;gsn_define_colormap(wks,"WhBlGrYeRe")
        gsn_define_colormap(wks,"precip2_17lev")
        ;;;;;;;;;;;;;; locations
        res@sfXArray = lon              ; rain is one-dimensional
        res@sfYArray = lat
        
        ;;;;;;;;;;;;;;set for the map
        ;res@mpDataBaseVersion       = "HighRes"
        ;res@mpDataSetName           = "Earth..4"
        res@mpGridAndLimbOn = False                         ; Turn off lat/lon lines
        res@mpOutlineOn           = True         ; Turn on map outlines
        res@mpGeophysicalLineColor = "Black"      ; Change the outline line color
        res@mpGeophysicalLineThicknessF= 2.      ; double the thickness of geophysical boundaries
        res@mpOutlineSpecifiers    = (/"China","Taiwan","Disputed area between India and China","India:Arunachal Pradesh"/) 
        
        ;;;;;;;;;;;;;;;ZOOM
        res@mpMinLatF = 30 
        res@mpMaxLatF = 34 
        res@mpMinLonF = 114 
        res@mpMaxLonF = 122 
        res@pmTickMarkDisplayMode = "Always"    
        
        ;;;;;;;;;;;;;; locations
        res@sfXArray = lon              ; rain is one-dimensional
        res@sfYArray = lat
  
        ;;;;;;;;;;;;;; contours
        res@cnFillOn = True 
        res@cnLinesOn = False           ; Turn off contour lines 
        res@cnLevelSelectionMode = "ExplicitLevels" 
        res@cnLevels = (/1, 2, 6, 10, 15, 20, 30, 40, 50, 70, 90, 110, 130, 150, 200/)
        ;;;;;;;;;;;;;; label bar
        res@lbBoxLinesOn = False ; Turn off label bar box lines 
  
        res@tiMainString    = fili(nf)
       ;res@gsnCenterString = fili(nf)
        plot = gsn_csm_contour_map(wks, rain, res)
    end if
   end do     ; nf loop
end