;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; Concepts illustrated:
;   - Drawing a scatter plot on a map
;   - Changing the marker color and size of stations on a map
;   - Plotting station locations using markers
;   - Manually creating a legend using markers and text
;
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;
; This example creates some station data, and then plots each value by coloring and sizing it
; depending on which range of values it falls in. It also plot the number of sea breeze days per year for
; all stations closer to the coast than 100 km
;
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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" 
load "landseabreeze_days04.ncl"                            
 
  begin 

	fils = systemfunc ("ls /media/amadou/COULIBALY12/DATA/CFSR/cfsr_data/wnd10m.gdas.20*.nc")
        f    = addfiles (fils, "r")              							; note the "s" of addfile
	time = f[:]->time
	u    = f[:]->U_GRD_L103
	v    = f[:]->V_GRD_L103

	lat = f[:]->lat
	lon = f[:]->lon
	
printVarSummary(lon)
printVarSummary(lat)
printVarSummary(u)
        
        time@units = "hours since 2001-01-01 00:00:00.0 +0:00"
 	ntim = dimsizes(time)
        iSize = dimsizes(time)
        date = cd_calendar(time, 0)
 	;printVarSummary(date)

;--> read the wind speed and direction from zonal and meridional winds
	wind_spd = sqrt(u^2 + v^2)
	wind_dir = wind_direction(u,v,0) 
;printVarSummary(wind_dir)

	wind_spd@coordinates = "time lat lon"
	wind_dir@coordinates = "time lat lon"

	wind_spd!0 = "time"
	wind_spd!1 = "lat"
	wind_spd!2 = "lon"
	wind_dir!0 = "time"
	wind_dir!1 = "lat"
	wind_dir!2 = "lon"

printVarSummary(wind_spd)

	wind_dir&time = time
	wind_dir&lat = lat
	wind_dir&lon = lon

	wind_spd&time = time
	wind_spd&lat = lat
	wind_spd&lon = lon

;printVarSummary(wind_dir)

	wind_spd1 = wind_spd(:, {6.16}, {1.25})
	wind_dir1 = wind_dir(:, {6.16}, {1.25})

        min_lon = min(lon)
        max_lon = max(lon)

        min_lat = min(lat)
        max_lat = max(lat)

;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Sea breeze dates %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         output1 = new((/12/), float)    ; Array which will contain the number of sea breeze days per month for all stations

          wdr = wind_dir1

	  out2 = LSBdaySelection(wdr, time) ; function selecting land sea breeze days
;%%%%%%%%%%%%%%%%%%%%% Number of Sea breeze days per month %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                do n = 0, 11

                     index = ind(out2(:,1) .ne. out2@_FillValue .and. out2(:,1) .eq. n+1)

                     if (.not.any(ismissing(index)) .and. dimsizes(index) .gt. 1) then

                             output1(n) = dimsizes(index)

                      end if 
                             delete(index)                        
                 end do  

             print(output1)
;%%%%%%%%%%%%%%%%%%%%%%%%%%  Monthly total occurrence of Land-Sea Breeze days %%%%%%%%%%%%%%%%%%%%%%%
          x = ispan(0,11,1)                     
          n = 12
;******************************************************************************************************
;                 create the plot
;******************************************************************************************************   
	  wks = gsn_open_wks("eps","Mon_Freq_Occur_CFSR_Lome")    ; open pdf file

         ; plot = new (5,graphic)

	  res                       = True                    ; plot mods desired         
	  res@gsnFrame              = False                   ; don't advance frame yet
          res@gsnDraw               = False 
	  res@gsnXYBarChart         = True                    ; turn on bar chart
	  res@gsnXYBarChartBarWidth = 0.5                     ; change bar widths
          res@gsnXYBarChartColors   = "blue"
	  res@gsnYRefLine           = 0
          
	  res@tmXBOn                = True                    ; turn off tickmarks at bot
	  res@trYMinF               = 0                       ; bring bars down to zero
	  res@trXMinF               = -1                      ; adds space on either end
	  res@trXMaxF               = 12                      ; of the 1st and last bars
	  res@tmXBMode              = "Explicit"              ; explicit labels
	  res@tmXBValues            = fspan(0,11,n)
	  res@tmXBLabels            = (/"J","F","M","A","M","J","J","A","S","O","N","D"/)
          res@tiYAxisString         = "Number of Land-Sea Breeze Days"
          res@tiXAxisString         = "Month"
          res@tiXAxisFontHeightF    = 0.015 
          res@tiYAxisFontHeightF    = 0.015 
        
	  plot = gsn_csm_xy(wks,x,output1(:),res)  
          ;plot(1) = gsn_csm_xy(wks,x,output1(7,:),res)
          ;plot(2) = gsn_csm_xy(wks,x,output1(12,:),res) 
          ;plot(3) = gsn_csm_xy(wks,x,output1(23,:),res) 
          ;plot(4) = gsn_csm_xy(wks,x,output1(25,:),res)   

          resP                        = True
	  ;resP@txString               = "Monthly Total Occurrence of Land-Sea Breeze Days"
	  resP@gsnPanelFigureStrings  = (/"Lome","Lome","Accra","Abidjan"/)  ; add strings to panel
	  resP@gsnPanelLabelBar       = True  
	  resP@gsnPaperOrientation    = "portrait"
	  resP@lbLabelFontHeightF     = 0.01
	  resP@gsnPanelBottom         = 0.05                   ; add space at bottom
	  resPgsnPanelTop             = 1.0
	  resP@gsnPanelScalePlotIndex = 0
	 ; resP@gsnPanelDebug          = True
	  resP@gsnMaximize            = True
	   
	  gsn_panel(wks,plot,(/1,1/),resP)          

  end