; readqnblendwind.ncl
; NCL program to read binary QSCAT/NCEP blended wind data files.
; This routine can be adapted to read either the zonal and
; meridional wind data ("uv") or the windstress curl data ("curl").
;
; Contributed by: Tianyi Fan, University of Colorado, Boulder; 2007-07-12
;

begin
 nlat = 353
 nlon = 720

 readuv = True   ; to read either uv or curlfiles
;  readuv = False

; ===============================
;  open file
; ===============================
 if (readuv) then
   fils = "./QSCAT/uv.200012.bln"
 else
   fils = "./QSCAT/curl.200012.bln"
 end if   
 print("intput" + fils)

 
  numRec = 30*4  ; number of record to read in, mamybe 31*4 or 
 
 day = new( numRec, "float")
 u = new((/ numRec, nlat, nlon/), "float")
 v = new((/ numRec, nlat, nlon/), "float")
 curl = new((/ numRec, nlat, nlon/), "float")
 
 if (readuv) then
   recl = 7+2*nlon*nlat      ; recl length
 else
   recl = 5+nlon*nlat
 end if  
 dum_all = new(recl, "float")
 
 print("Assigning coordinate varibal information")
 
 u!0 = "time"
 u&time = ispan(1, numRec, 1)
 u!1 = "lat"
 u&lat = fspan(-88., 88., nlat)
 u&lat@units = "degrees_north"
 u!2 = "lon"
 u&lon = fspan(0.5,360.,nlon)
 u&lon@units = "degrees_east"
 
 v!0 = "time" 
 v&time = ispan(1, numRec, 1)
 v!1 = "lat"
 v&lat = fspan(-88., 88., nlat)
 v&lat@units = "degrees_north"
 v!2 = "lon"
 v&lon = fspan(0.5,360.,nlon)
 v&lon@units = "degrees_east"


 curl!0 = "time" 
 curl&time = ispan(1, numRec, 1)
 curl!1 = "lat"
 curl&lat = fspan(-88, 88, nlat)
 curl&lat@units = "degrees_north"
 curl!2 = "lon"
 curl&lon = fspan(0.5,360.,nlon)
 curl&lon@units = "degrees_east"
 
; ===============================
;  read data
; =============================== 
  if(readuv) then
    do i = 0, numRec -1   
      dum_all    = fbindirread(fils, i, recl, "float")       
;      idum1      = dum_all(0)
      day(i)     = dum_all(1)
;      idum2      = dum_all(2)
;      idum3      = dum_all(3)
      u(i, :,:)  = onedtond(dum_all(4:4+nlon*nlat-1),(/nlat, nlon/) )
;      idum4      = dum_all(4+nlon*nlat)
;      idum5      = dum_all(5+nlon*nlat)
      v(i, :,:)  = onedtond(dum_all(6+nlon*nlat:6+nlon*nlat*2-1),(/nlat, nlon/) )
;      idum6      = dum_all(6+nlon*nlat*2)
    end do
  else
    do i = 0, numRec -1
      dum_all    = fbindirread(fils, i, recl, "float")  
      idum1      = dum_all(0)
      day(i)     = dum_all(1)
      idum2      = dum_all(2)
      idum3      = dum_all(3)
      curl(i,:,:)= onedtond(dum_all(4:4+nlon*nlat-1),(/nlat, nlon/) ) 
      idum4      = dum_all(4+nlon*nlat)        
    end do
  end if
  

;============================
;  print
; ===========================
  if (readuv) then
    do i = 0, 5
      ip1 = 20
      ip2 = 24
      jp  = 19
      print("day = " + day(i) + "  print data at i = " + (ip1+1) + "  to  " + (ip2+1) + ", j = " + (jp+1)  )
      print(" u = " + u(i, jp, ip1:ip2))
      print(" v = " + v(i, jp, ip1:ip2))
      
      jp  = 299
      print("day = " + day(i) + "  print data at i = " + (ip1+1) + "  to  " + (ip2+1) + ", j = " + (jp+1)  )
      print(" u = " + u(i, jp, ip1:ip2))
      print(" v = " + v(i, jp, ip1:ip2))
    end do  
  else
    do i = 0, 5
      ip1 = 20
      ip2 = 24
      jp  = 19
      print("day = " + day(i) + "  print data at i = " + (ip1+1) + "  to  " + (ip2+1) + ", j = " + (jp+1)  )
      print(" c = " + curl(i, jp, ip1:ip2))
      
      jp  = 299
      print("day = " + day(i) + "  print data at i = " + (ip1+1) + "  to  " + (ip2+1) + ", j = " + (jp+1)  )
      print(" c = " + curl(i, jp, ip1:ip2))
    end do    
  end if
 
  
end