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" 
;script;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

begin
;set for the necessary varibles;;;;;;;;;;;;;;;;;;
  type     = "pdf"                     
  wksname  = "GPM"
  wks      = gsn_open_wks(type,wksname)
  gsn_define_colormap(wks,"precip3_16lev")
;;;;;;read in GPM files and get the cordinations;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 iFiles = systemfunc("ls /public/home/GlobeDisRadi/HRB2015/data/gpm20150624_regional/*.nc*")
; print(iFiles)
 nFiles = dimsizes(iFiles)  

 iFile    = addfile(iFiles(0),"r")
 gpmLat  = iFile->lat(:)
 gpmLon  = iFile->lon(:)
 nGpmLat = dimsizes(gpmLat)
 nGpmLon = dimsizes(gpmLon)
 minGpmLat = min(gpmLat)
 maxGpmLat = max(gpmLat)
 minGpmLon = min(gpmLon)
 maxGpmLon = max(gpmLon)

printVarSummary(gpmLat)
printVarSummary(gpmLon)

;;;;;;;;;get  gpm pre aray and ascii values;;;;;;;;;;;;;;;

gpmPre = new((/nFiles,nGpmLon,nGpmLat/),"float")
gpmTime = new(nFiles,string)
gpmAscii = new(nGpmLat,string)
do i =0,nFiles-1
      tmpFile = addfile(iFiles(i),"r")
      gpmPre(i,:,:) = tmpFile->precipitationCal(:,:)
     
     fn = stringtochar(iFiles(i))
  ;print(fn)
     asciiTxt = chartostring(fn(81:109))+".txt"

   do j =0,nGpmLat-1
      gpmAscii(i) = str_concat(sprintf("%10.3f",gpmPre(i,:,j)))
   end do 
   
      asciiwrite(asciiTxt,gpmAscii)
end do

if(any(ismissing(gpmPre))) then
    print("Your data contains some missing values. Beware.")
end if

;;;;;;;;;;set the resource plot;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  res = True
  res@gsnDraw = False
  res@gsnFrame= False     ;overlay, 必须设置
  res@gsnAddCyclic = False

  res@mpFillOn               = False
  res@mpOutlineOn            = True
  res@mpOutlineBoundarySets  = "National"   ; turn on country boundaries
  res@mpNationalLineColor    = "Black"

  res@gsnCenterString         ="GPM daily accumulated precipitation"
  res@gsnCenterStringFontHeightF      = 0.02
  res@gsnCenterStringOrthogonalPosF   = 0.14
  res@gsnLeftString           ="Precipitation (mm/day)"

  res@cnFillOn = True
  res@cnLinesOn = False
  res@cnLineLabelsOn = False
  res@cnInfoLabelOn = False; set false will not set an information label
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF = 1
  res@cnMaxLevelValF = 200
  res@cnLevelSpacingF = 5

  res@lbBoxLinesOn = False ; Label bar
  res@lbLabelAutoStride     = True            ; let NCL determine label spacing
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do i =0,1; nFiles-1
   fn = stringtochar(iFiles(i))
  ;print(fn)
   gpmTime = chartostring(fn(81:109))
    res@gsnRightString        =  gpmTime

 plot= gsn_csm_contour_map(wks,gpmPre(i,:,:),res)

;;;;add the basin polylines;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  polylinefile        = "Huaihe_Basin.shp"
  plres               = True
  plres@gsLineColor   ="red"  ;  RGB:30 144 255
  plres@gsLineThicknessF =1.5
 poly11 = gsn_add_shapefile_polylines(wks,plot,polylinefile,plres)

 draw (plot)

 frame(wks)
end do

end