; read GEFS data from multiple directories
; requires (or cancel) a file that has a list of times
; contact barry.h.lynn@gmail.com if you have comments, improvements to suggest
; thank you to ncl-talk participants and the NCL on-line help
; thank you to a couple of folk who posted shell scripts for reading NCEP files
; for example: "www.meteo-blog.net"   


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/contrib/cd_string.ncl"
external GET_IJ "./get_ij.so"

undef("get_ij_ncl(lat2d,lon2d,lat_beg,lon_be")
function get_ij_ncl(lat2d,lon2d,lat_beg,lon_beg)
local nj, ni
begin
  i_loc = 1
  j_loc = 1
  nj = dimsizes(lat2d(:,0))
  ni = dimsizes(lat2d(0,:))
  GET_IJ::get_ij(lat2d,lon2d,lat_beg,lon_beg,i_loc,j_loc,nj,ni)
  return((/i_loc,j_loc/))
end

undef("average_subset(x,lat,lon,xy_start,xy_en")
function average_subset(x,lat,lon,xy_start,xy_end)
begin
  x_ave_new = x(xy_start(1):xy_end(1),xy_start(0):xy_end(0))
  x_ave_new!0 = "lat"
  x_ave_new!1 = "lon"
  x_ave_new&lat= lat(xy_start(1):xy_end(1))
  x_ave_new&lon= lon(xy_start(0):xy_end(0))
  return(x_ave_new)
end

undef("set_res1_list")
function set_res1_list()
begin
  res                     = True
  res@mpDataBaseVersion = "Ncarg4_1"           ; choose more recent database
  res@mpDataSetName = "Earth..4"               ; high resolution
  res@mpOutlineBoundarySets = "National"       ; National borders
  res@mpNationalLineColor = "Black"            ; National borders color
  res@mpNationalLineThicknessF    = 3.0
  res@mpGeophysicalLineThicknessF = 2.0
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnFillOn            = True            ; turn on color
  res@cnLinesOn           = False           ; no contour lines
  res@gsnLeftString       = "700 mb Wind Vectors"         ; change left string
  res@gsnRightString      = "km/h"         ; assign right string
  res@mpFillOn            = False           ; no map fill
  res@cnLevelSelectionMode = "ExplicitLevels"    ; set explicit contour levels
  res@cnLevels    = (/ -34.,-32.,-30.,-28.,-26.,\
                       -24.,-22.,-20.,-18.,-16.,-14.,\
                       -12.,-10.,-8.,-6.,-4.,-2.,0.,\
                       2.,4.,6.,8.,10.,12.,14.,16.,18.,\
                       20.,22.,24.,26.,28.,30.,32./)
  res@cnFillColors = (/ 10,12,14,15,17,19,22,25,31,35,39,43,48,53,\
                        61,64,67,71,89,95,99,104,108,111,115,123,\
                        130,135,143,148,150,155,159,163,165/)
                         ; set the colors to be used
; res1 label bar
  res@cnInfoLabelOn = False                  ; turn off contour info label
  res@lbAutoManage          = False          ; control label bar
  res@pmLabelBarDisplayMode = "Always"       ; turns on label bar
  res@lbOrientation         = "Vertical"
  res@pmLabelBarSide        = "Right"
  res@lbLabelAutoStride =   True
  res@pmLabelBarWidthF      = 0.20
  res@pmLabelBarHeightF     = 0.7
  res@lbLabelFontHeightF    = .025
  res@lbPerimOn             = False
  res@lbLabelFont      = "Helvetica"         ; label font
  res@lbTitleOn        =  True               ; turn on title
  res@lbTitleString    = "Temp [~S~o~N~C]"        ; title string
  res@lbTitlePosition  = "Top"               ; title position
  res@lbTitleFontHeightF= .020               ; make title smaller
  res@lbTitleDirection = "Across"            ; title direction
  res@lbTitleOffsetF  = 0.002 ; label bar title position
  res@lbTopMarginF = 0.001
  res@gsnAddCyclic          = True           ; data already has cyclic point
  return(res)
end

undef("set_res2_list")
function set_res2_list()
begin
  res                     = True
  res@mpDataBaseVersion = "Ncarg4_1"           ; choose more recent database
  res@mpDataSetName = "Earth..4"               ; high resolution
  res@mpOutlineBoundarySets = "National"       ; National borders
  res@mpNationalLineColor = "Black"            ; National borders color
  res@mpNationalLineThicknessF    = 3.0
  res@mpGeophysicalLineThicknessF = 2.0
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnFillOn            = True            ; turn on color
  res@cnFillPalette       = "gui_default"   ; set color map
  res@cnLinesOn           = False           ; no contour lines
  res@gsnLeftString       = ""         ; change left string
  res@gsnRightString      = ""        ; assign right string
  res@mpFillOn            = False           ; no map fill
  res@lbTitleString    = "RH [%]"        ; title string
  res@cnInfoLabelOn = False                  ; turn off contour info label
  res@lbAutoManage          = False          ; control label bar
  res@pmLabelBarDisplayMode = "Always"       ; turns on label bar
  res@lbOrientation         = "Vertical"
  res@pmLabelBarSide        = "Right"
  res@lbLabelAutoStride =   True
  res@pmLabelBarWidthF      = 0.20
  res@pmLabelBarHeightF     = 0.7
  res@lbLabelFontHeightF    = .025
  res@lbPerimOn             = False
  res@lbLabelFont      = "Helvetica"         ; label font
  res@lbTitleOn        =  True               ; turn on title
  res@lbTitleString    = "RH [%]"        ; title string
  res@lbTitlePosition  = "Top"               ; title position
  res@lbTitleFontHeightF= .020               ; make title smaller
  res@lbTitleDirection = "Across"            ; title direction
  res@lbTitleOffsetF  = 0.002
  res@lbTopMarginF = 0.001
  return(res)
end

undef("set_res3_list")
function set_res3_list()
begin
  res   = True
  res@gsnAddCyclic          = False          ; we will manually zoom in.
  res@gsnDraw              = False             ; do not draw the plot
  res@gsnFrame             = False             ; do not advance the frame
  res@gsnLeftString       = "700 mb Geopotential Height"         ; change left string
  res@gsnRightString      = "dam"        ; assign right string
  res@cnConstFLabelPerimOn           = False
  res@cnLevelSelectionMode           = "ManualLevels"
  res@cnMinLevelValF                 = 250
  res@cnMaxLevelValF                 = 350
  res@cnLevelSpacingF                = 5
  res@cnLineLabelAngleF              = 0.0
  res@cnLineLabelFontHeightF         = 0.01
  res@cnLineLabelBackgroundColor     = -1
  res@gsnContourLineThicknessesScale = 1.5
  res@tfDoNDCOverlay = True
  return(res)
end

undef("set_vecres_list")
function set_vecres_list()
begin
  res                  = True            ; vector only resources
  res@gsnDraw          = False           ; don't draw
  res@gsnFrame         = False           ; don't advance frame
  res@vcGlyphStyle     = "CurlyVector"   ; curly vectors
  res@vcRefMagnitudeF  = 20              ; define vector ref mag
  res@vcRefLengthF     = 0.045           ; define length of vec ref
  res@gsnRightString   = " "             ; turn off right string
  res@gsnLeftString    = " "             ; turn off left string
  res@tiXAxisString    = " "             ; turn off axis label 
  res@vcRefAnnoOrthogonalPosF = 0.01    ; move ref vector into plot
  res@gsnAddCyclic          = False          ; we will manually zoom in.
  return(res)
end

;***************************************

begin
  plot = new(2,graphic)                           ; create a plot array
; define number of time variables
  n_t_vars=5
  n_l_times = 0
; ascii file
  a_asc = (/"gefs_time.file"/)
;---Read the data in as a small1D array
  values = asciiread(a_asc,-1,"float")
; set how many time steps
  n_l_times = 39
  nvars  = 5
  nxy    = 1
  nblk   = nxy * nvars
; temporary
  ny=1
  nx=1
; create variables to hold the ascii file variables.
;
  n_t_vars=5
  time_array     = new((/n_t_vars/),integer)
  year     = new((/n_l_times,ny,nx/),integer)
  month    = new((/n_l_times,ny,nx/),integer)
  day      = new((/n_l_times,ny,nx/),integer)
  hour     = new((/n_l_times,ny,nx/),integer)
  minute     = new((/n_l_times,ny,nx/),integer)
  do it=0,n_l_times-1
      nstart = it*(nblk+1)+1
      nend   = nstart+nblk-1
      year(it,:,:)    = toint(reshape(values(nstart:nend:nvars),(/ny,nx/)))
      month(it,:,:)    = toint(reshape(values(nstart+1:nend:nvars),(/ny,nx/)))
      day(it,:,:)    = toint(reshape(values(nstart+2:nend:nvars),(/ny,nx/)))
      hour(it,:,:)    = toint(reshape(values(nstart+3:nend:nvars),(/ny,nx/)))
      minute(it,:,:)    = toint(reshape(values(nstart+4:nend:nvars),(/ny,nx/)))
  end do
;      print("year = " + year(:,0,0))
;      print("day = " + day(:,0,0))
  DATADir = "/home/cust100021_vol1/barry/cust100021_vol2/GEFS"
; find directories with this date, etc
  cmd = "find " + DATADir + "/GEFS* -type d -print"
  dirWRF = systemfunc (cmd)
; print(dirWRF)
; this is the total number of directories within the listed directory header
; it can be just the number of, for example, "GEFS_XX" directories or also all subdirectories within
  n_dirWRF = dimsizes(dirWRF)
; print("n_dirWRF" + n_dirWRF)
; set name of output file
  outfile="gefs_700mb_maps"
  wks = gsn_open_wks("pdf",outfile)                 ; open wk station
  gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200")
; http://www.ncl.ucar.edu/Document/Functions/Built-in/systemfunc.shtml
; to retrieve a list of files from another directory
  diri = dirWRF(0) + "/"
  dir_files = systemfunc("ls " + diri + "ge*pgrb2b*grb")
  n_files = dimsizes(dir_files)
; print("n_files = " + n_files)

; Create various resource lists
  res1 = set_res1_list()
  res2 = set_res2_list()
  res3 = set_res3_list()
  vecres = set_vecres_list()

;; Map Projection (min and max Lon and Lat)
  lat_beg = 20
  lat_end = 50
  lon_beg = 20
  lon_end = 50

;   res1@mpProjection          = "LambertConformal"
;   The GEFS is not!
  res1@mpMinLatF            = lat_beg
  res1@mpMaxLatF            = lat_end
  res1@mpMinLonF            = lon_beg
  res1@mpMaxLonF            = lon_end
  res1@mpLimitMode="LatLon"
  res2@mpMinLatF            = lat_beg
  res2@mpMaxLatF            = lat_end
  res2@mpMinLonF            = lon_beg
  res2@mpMaxLonF            = lon_end
  res2@mpLimitMode="LatLon"

  resP                = True                ; panel only resources
  resP@gsnMaximize    = True                ; maximize plots

;;; Now that we have identified the number of files and number of directories
  all_files     = new((/n_files,n_dirWRF/),string)
  do i_dir = 0,n_dirWRF-1
     diri = dirWRF(i_dir) + "/"
     all_files(:,i_dir) = systemfunc("ls " + diri + "ge*pgrb2b*grb")
  end do
; print("all_files = " + all_files)
; loop over number of files within each directory
; do n=0,n_files -1,1
  do n=0,38,2
    print("n times steps = " + n)
; loop over directories
    n_dir = 0
    do i_dir = 0,n_dirWRF-1
     check = ismissing(all_files(n,i_dir))
;    print(check)
     if (check.eq.True)then
        print("missing!!!")
      end if
     if (check.eq.False)then
      n_dir = n_dir + 1
;     print("all_files(n,i_dir) = " + all_files(n,i_dir))
      fin = addfile(all_files(n,i_dir),"r")

;     print("lv0 = " + fin->lv_ISBL0)
;     print("lv8 = " + fin->lv_ISBL8)
;     print("lv10 = " + fin->lv_ISBL10)
      t = fin->TMP_P1_L100_GLL0(17,:,:)
      rh = fin->RH_P1_L100_GLL0(17,:,:)
      z_up = fin->HGT_P1_L100_GLL0(16,:,:)
      z_dn = fin->HGT_P1_L100_GLL0(17,:,:)
      z =0.5*(z_up+z_dn)
;     print("z ave = "+ z(20,20))
      copy_VarCoords(z_up,z)                         ; copy coord vars to speed
      u_up = fin->UGRD_P1_L100_GLL0(13,:,:)
      u_dn = fin->UGRD_P1_L100_GLL0(14,:,:)
      v_up = fin->VGRD_P1_L100_GLL0(13,:,:)
      v_dn = fin->VGRD_P1_L100_GLL0(14,:,:)
      u =0.5*(u_up+u_dn)
      v =0.5*(v_up+v_dn)
      copy_VarCoords(u_up,u)                         ; copy coord vars to speed
      copy_VarCoords(v_up,v)                         ; copy coord vars to speed
      delete([/u_up,v_up,u_dn,v_dn/])                ; remove vars we no longer need
      lon = fin->lon_0(:)
      lat = fin->lat_0(:)
      if (i_dir.eq.0)then
        t_ave = t
        rh_ave = rh
        z_ave = z
        u_ave = u
        v_ave = v
      end if
      if (i_dir.ne.0)then
        t_ave = t + t_ave
        rh_ave = rh + rh_ave
        z_ave = z + z_ave
        u_ave = u + u_ave
        v_ave = v + v_ave
      end if
      end if
;  make plot
      if (i_dir.eq.n_dirWRF-1)then
        print("n_dir = " + n_dir)
        t_ave = t_ave/n_dir
        rh_ave = rh_ave/n_dir
        z_ave = z_ave/n_dir
        u_ave = u_ave/n_dir
        v_ave = v_ave/n_dir
        t_ave=t_ave-273.16  ; convert to C
        z_ave=z_ave/10.  ; convert to dam
;       printMinMax(z_ave,False)
        u = u * 3.6
        v = v * 3.6
        copy_VarCoords(u,u_ave)                         ; copy coord vars to speed
        copy_VarCoords(v,v_ave)                         ; copy coord vars to speed
        copy_VarCoords(z,z_ave)                         ; copy coord vars to speed
        copy_VarCoords(t,t_ave)                         ; copy coord vars to speed
        copy_VarCoords(rh,rh_ave)                         ; copy coord vars to speed
; map perameters
        new_time_units = "hours since 1800-01-01 00:00"
        new_time_array = cd_inv_calendar( year(n,0,0),  month(n,0,0), \
               day(n,0,0), hour(n,0,0), minute(n,0,0), 0, new_time_units,0)
;       print("new_time_array = " + new_time_array)
        format = "" ;; defaults to "%H%M IST %d %c %Y"
        format = "%H%M IST %d %c %Y"
        new_time_array = new_time_array + 2.0
        time_var=cd_string(new_time_array,format)
        res1@tiMainString         = time_var
        res2@tiMainString         = time_var
                                                ; this must also be set for any zoom
; must manuall excerpt variable
        dims2d = dimsizes(z)
        lat2d = conform_dims(dims2d,lat,0)  
        lon2d = conform_dims(dims2d,lon,1)

        xy_start = get_ij_ncl(lat2d,lon2d,lat_beg,lon_beg)
        xy_end   = get_ij_ncl(lat2d,lon2d,lat_end,lon_end)
;       print("xy_start(0) = " + xy_start(0))
;       print("xy_start(1) = " + xy_start(1))
;       print("xy_end(0)   = " + xy_end(0))
;       print("xy_end(1)   = " + xy_end(1))

        u_ave_new = average_subset(u_ave,lat,lon,xy_start,xy_end)
        v_ave_new = average_subset(v_ave,lat,lon,xy_start,xy_end)
        z_ave_new = average_subset(z_ave,lat,lon,xy_start,xy_end)
        delete([/u_ave,v_ave,z_ave/])

        plot(0) = gsn_csm_contour_map(wks,t_ave,res1)
        plot(1) = gsn_csm_contour_map(wks,rh_ave,res2)
        plotB   = gsn_csm_vector(wks,u_ave_new(::4,::4),v_ave_new(::4,::4),vecres)
        plotD   = gsn_csm_contour(wks,z_ave_new,res3)
        overlay(plot(0),plotB)
        overlay(plot(1),plotD)
        gsn_panel(wks,plot,(/2,1/),resP)          ; now draw as one plot
      end if
    end do
  end do
end