;*************************************************
; WRF: create time series plot
;************************************************
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/wrf/WRF_contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;************************************************
; open file and read in data
;************************************************
;fils = systemfunc ("ls wrfout_d03_*.nc") ; file paths
;f    = addfiles (fils, "r")   
;ListSetType (f, "cat")        ; concatenate (=default)
f    = addfile("wrfout_d03_2009_07_07_18_00_00.nc", "r")   
;************************************************
; Read all the times. Convert to numeric units for plottin.
;************************************************
;; Times = f[:]->Times                       ; Times is of type character
;; Time  = WRF_Times2Udunits_c(Times, 0)  ; convert to "hours since"
;; printVarSummary (Times)

;; inetcdf=018
;; jnetcdf=018
;; i=inetcdf-1
;; j=jnetcdf-1
;************************************************
; Import time series of perturbation geopotential (PH)
; at a (arbitrarily) specified grid point
;************************************************
;;  y     = f[:]->PH(:,8,118,118)                ; (Time, bottom_top, south_north, west_east)
;;  y     = f[:]->PBLH(:,i,j)                ; (Time, bottom_top, south_north, west_east)
;;  y     = f[:]->T2(:,i,j) - 273.15          ; (Time, bottom_top, south_north, west_east)
;;  y     = f[:]->TSLB(:,0,i,j) - 273.15          ; (Time, bottom_top, south_north, west_east)
;;   y@description = "Surface Temperature" 
;;y     = f[:]->HFX(:,i,j) - 273.15          ; (Time, bottom_top, south_north, west_east)
;;   y     = f[:]->PH(:,8,{43.69439478},{5.74614644}); (Time, bottom_top, south_north, west_east)
;Longitude en degrés décimaux :	5.74614644
;Latitude en degrés décimaux : 	43.69439478

;;printVarSummary (y)
;************************************************
; subtract the initial value to create anomalies
;************************************************
;; yAnom             = y-y(0)             ; anomalies from init time
;; yAnom@description = "Anomalies from initial Time" 
;; yAnom@units       = y@units
;************************************************
; For plot label, read the lat/lon location     
;************************************************
ip_locs = (/ "Pylone meteo Cadarache" /)

ip_lats = (/ 5.74614644/)

ip_lons = (/ 43.69439478/)


;; lat   = f[0]->XLAT(0,i,j)
;; lon   = f[0]->XLONG(0,i,j)
;;   lat   = f[0]->XLAT(0,{43.69439478},{5.74614644})
;;   lon   = f[0]->XLONG(0,{43.69439478},{5.74614644})

do ip = 1, 1    ; LOOP through above stations locations and


  ; Get ij point in model domain for location "ip"
  ; loc(0) is south-north (y) and loc(1) is west-east (x)
  ; > TH BUG :
  ;  loc = wrf_user_latlon_to_ij(a, ip_lats(ip), ip_lons(ip))
  ;  obsolete en NCL 5 cf NCL/NCL-5.1.1/ncl_ncarg-5.1.1/ni/src/examples/gsun/WRFUserARW.ncl
  opt = True
  opt@returnInt = True      ; Return real values (set to True for integer values - integer is default)
  opt@useTime = 0           ; times(it)

  loc = wrf_user_ll_to_ij(f, ip_lats(ip), ip_lons(ip), opt)

  y     = f[:]->HFX(:,loc(0),loc(1)) - 273.15          ; (Time, Lat, Lon)
  printVarSummary (y)

  yAnom             = y-y(0)             ; anomalies from init time
  yAnom@description = "Anomalies from initial Time" 
  yAnom@units       = y@units



  print("Attempting to plot: " + "Station - " + ip_locs(ip) )
  print("                    " + "at location: "+ ip_lats(ip) +" ; "+ ip_lons(ip) )
  print("                    " + "and (i,j)  = ("+ loc(0) +", "+ loc(1)+")" )

  if ( ismissing(loc(0)) .or. ismissing(loc(1)) )
    if ( FirstTime)
      print("                    " + "SKIP: Station outside model domain" )
    end if
  else
     
     					;************************************************
     					; create plots: Three slight variations.  
     					;************************************************
      wks1 = gsn_open_wks("x11" ,"WRF_xy1")          ; ps,pdf,x11,ncgm,eps
;;    wks2 = gsn_open_wks("x11" ,"WRF_xy2")          ; ps,pdf,x11,ncgm,eps
;;    wks3 = gsn_open_wks("x11" ,"WRF_xy3")          ; ps,pdf,x11,ncgm,eps
      
      res                       = True            ; plot mods desired
      res@gsnMaximize           = True            ; uncomment to maximize size
      res@xyLineThicknessF      = 3.0             ; make a bit thicker           
      res@tiMainString          = "("+inetcdf+","+jnetcdf+") "+lat+"N  "+fabs(lon)+"W"
      plot                      = gsn_csm_xy(wks1,Time,y,res)
      
      ;;   res@tiXAxisString         = Time@units      ; label bottom axis with units
      ;;   res@xyLineThicknessF      = 3.0             ; make a bit thicker           
      ;;   res@gsnYRefLine           = 0.0             ; draw a reference line
      ;;   plot                      = gsn_csm_xy(wks1,Time,yAnom,res)
      
      ;;   res@gsnAboveYRefLineColor = "red"           ; above ref line fill red
      ;;   res@gsnBelowYRefLineColor = "blue"          ; below ref line fill blue
      ;;   plot                      = gsn_csm_xy(wks1,Time,yAnom,res)
    end if
  end do     ; END OF LOCATIONS
end