; Reads the best storm tracks of tropical cyclones from G5NR
; and from the GFS CTL and RO runs
; and computes an average track and intensity error
; as well as individual storm tracks for each case
;
;
; Andrew Kren
; June 16, 2017

begin

; begin user's settings
;---------------------------------------------------------

path_g5nr = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/G5NR/"
path_ctl = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/CTL/"
path_exp = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/C2REF/"

; list files in directory to determine how many dropsondes to read
fils_g5nr = systemfunc("csh -c 'cd " + path_g5nr + " ; ls *.tc'")
fils_ctl = systemfunc("csh -c 'cd " + path_ctl + " ; ls *.txt'")
fils_exp = systemfunc("csh -c 'cd " + path_exp + " ; ls *.txt'")

output_file = "storm_tracks_G5NR"  ; output figure file
output_type = "png" ; "x11", "ps", etc

input_file_g5nr = path_g5nr + fils_g5nr
input_file_ctl =  path_ctl + fils_ctl
input_file_exp =  path_exp + fils_exp

scale_factor = 0.1
scale_factor_west = -0.1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; end of user's settings
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;
; This procedure adds a line to a plot, making sure that it is returned
; to a unique variable name, and that this variable is retained even
; outside this procedure call.
;
procedure add_line(wks,plot,xval,yval,line_color)
local plres, str, y
begin
  plres = True
  plres@gsLineColor = line_color   ; Set the line color.
  plres@gsLineThicknessF = 4

  
  str = unique_string("dum")   ; "unique_string" will return a unique
                                    ; string every time it is called from
                                    ;  within a single NCL session.
;
; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.
;
  plot@$str$ = gsn_add_polyline(wks, plot, xval, yval, plres)
end

procedure add_mark(wks,plot,xval,yval,mark_color)
local res_mark, str, y
begin
   res_mark = True

   res_mark@gsMarkerIndex    = 1       ; polymarker style
   res_mark@gsMarkerSizeF    = 0.012   ; polymarker size

  str = unique_string("duma")   ; "unique_string" will return a unique
                                    ; string every time it is called from
                                    ;  within a single NCL session.
;
; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.
;
  plot@$str$ = gsn_add_polymarker(wks, plot, xval, yval, res_mark)
end

; create arrays to store storm track data
g5nr_pres = new((/25,100/),"float")
g5nr_name = new((/25/),"string")
g5nr_lat = g5nr_pres
g5nr_lon = g5nr_pres
g5nr_ymd = new((/25,100/),"string")

ctl_pres = new((/25,3000/),"float")
ctl_name = new((/25/),"string")
ctl_lat = ctl_pres
ctl_lon = ctl_pres
ctl_ymd = new((/25,3000/),"string")
ctl_hr = ctl_ymd

exp_pres = new((/25,3000/),"float")
exp_name = new((/25/),"string")
exp_lat = ctl_pres
exp_lon = ctl_pres
exp_ymd = new((/25,3000/),"string")
exp_hr = exp_ymd

track_error_ctl = new((/25,20,65/),"float") ; saved track errors and intensity errors, 65 corresponds to 0 to 384 hours of 6hr output
intensity_error_ctl = track_error_ctl
track_error_exp = track_error_ctl
intensity_error_exp = track_error_ctl
fcst_hr_array_ctl = new((/65/),"float") ; number of forecasts in each 6hr forecast segment for statistics
fcst_hr_array_exp = fcst_hr_array_ctl
ctl_lat_array = track_error_ctl
ctl_lon_array = track_error_ctl
exp_lat_array = track_error_ctl
exp_lon_array = track_error_ctl

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; loop over each of the 25 storms, plotting the storm track at each iteration
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

print("Reading G5NR storm tracks...")
do i=0,dimsizes(fils_g5nr)-1 

; read in data for this storm
 data1 = readAsciiTable(input_file_g5nr(i),19,"string",0) ; 19 column array
 
 data_1d = ndtooned(data1)
 
 ; loop over data_1d and read data using str_get_cols
 count = num(.not.ismissing(data_1d))
 
 lat = tofloat(str_get_cols(data_1d,116,118)) * scale_factor ; latitudes
 lon = tofloat(str_get_cols(data_1d,121,124)) ; longitudes
 deg_we = str_get_cols(data_1d,125,125)
 
 do j=0,count-1
 
 if (deg_we(j).eq."E") then
   lon(j) = lon(j) * scale_factor      ; deg W
 end if
 if (deg_we(j).eq."W") then
   lon(j) = lon(j) * scale_factor_west ; deg E
 end if
 
 end do
 
 pres = tofloat(str_get_cols(data_1d,135,138)) ; pressure in hPa at time and location
 storm_name = str_get_cols(data_1d,92,95) ; Storm name in TC vitals file
 storm_cat = tofloat(str_get_cols(data_1d,176,177)) ; Storm category, e.g. cat 1 hurricane, 2, 3, 0 is TD/TC
 ymd = str_get_cols(data_1d,102,109) ; YMD
 hr = str_get_cols(data_1d,111,112)  ; HH
 date = ymd ; initialize
 month = ymd
 day = ymd
 
 do j=0,count-1
 
 temp_concat = (/ymd(j),hr(j)/)
 date(j) = str_concat(temp_concat)
 month(j) = str_get_cols(date(j),4,5)
 day(j) = str_get_cols(date(j),6,7)
 
 end do
 
 ; store G5NR data into arrays for later use
 size_good = num(.not.ismissing(pres)) ; number of good elements in array since when reading some are missing due to reading multiple lines in file
 
 g5nr_pres(i,0:size_good-1) = pres(0:size_good-1)
 g5nr_name(i) = storm_name(0)
 g5nr_lat(i,0:size_good-1) = lat(0:size_good-1)
 g5nr_lon(i,0:size_good-1) = lon(0:size_good-1)
 g5nr_ymd(i,0:size_good-1) = date(0:size_good-1)
 
 
delete([/data1,data_1d,lat,lon,deg_we,pres,storm_name,storm_cat,ymd,hr,date,month,day,temp_concat/])
 
end do ; end loop over 25 G5NR storms

; read the control and experiment storm track files

; CTL
print("Reading CTL and EXP storm tracks...")
do i=0,dimsizes(fils_ctl)-1

; open file and read
data1 = asciiread(input_file_ctl(i),-1,"string") ; 20 column array of values

delim = ","
ctl_ymd(i,0:dimsizes(data1)-1) = str_get_field(data1,3,delim)
ctl_hr(i,0:dimsizes(data1)-1) = str_get_field(data1,6,delim)
ctl_lat(i,0:dimsizes(data1)-1) = tofloat(str_get_cols(str_get_field(data1,7,delim),0,3)) * scale_factor
lon = str_get_cols(str_get_field(data1,8,delim),0,4)
ctl_pres(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,10,delim))
ctl_name(i) = str_get_cols(input_file_ctl(i),48,51)

deg_we = str_get_cols(str_get_field(data1,8,delim),5,5)

do j=0,dimsizes(data1)-1 ; assign correct scale factor to longitudes based on degW or degE
 
 if (deg_we(j).eq."E") then
   lon(j) = tofloat(lon(j)) * scale_factor      ; deg W
 end if
 if (deg_we(j).eq."W") then
   lon(j) = tofloat(lon(j)) * scale_factor_west ; deg E
 end if
 
 end do
 
ctl_lon(i,0:dimsizes(data1)-1) = tofloat(lon)

delete([/data1,lon,deg_we/])

; EXP 
; open file and read
data1 = asciiread(input_file_exp(i),-1,"string") ; 20 column array of values

delim = ","
exp_ymd(i,0:dimsizes(data1)-1) = str_get_field(data1,3,delim)
exp_hr(i,0:dimsizes(data1)-1) = str_get_field(data1,6,delim)
exp_lat(i,0:dimsizes(data1)-1) = tofloat(str_get_cols(str_get_field(data1,7,delim),0,3)) * scale_factor
lon = str_get_cols(str_get_field(data1,8,delim),0,4)
exp_pres(i,0:dimsizes(data1)-1) = tofloat(str_get_field(data1,10,delim))
exp_name(i) = str_get_cols(input_file_exp(i),50,53)

deg_we = str_get_cols(str_get_field(data1,8,delim),5,5)

do j=0,dimsizes(data1)-1 ; assign correct scale factor to longitudes based on degW or degE
 
 if (deg_we(j).eq."E") then
   lon(j) = tofloat(lon(j)) * scale_factor      ; deg W
 end if
 if (deg_we(j).eq."W") then
   lon(j) = tofloat(lon(j)) * scale_factor_west ; deg E
 end if
 
 end do
 
exp_lon(i,0:dimsizes(data1)-1) = tofloat(lon)

delete([/data1,lon,deg_we/])

end do ; end loop over control and experiment forecasts



; compute storm track and intensity error over each storm individually (maps for individual ones) and save for computing statistics over all storms, including t-test
print("Computing track and intensity errors for CTL...")
; first for CTL only since the dates between CTL and EXP are not always the same so treat separately
do i=0,dimsizes(g5nr_name)-1 

; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm

storm_index = ind(ctl_name.eq.g5nr_name(i))

if (.not.ismissing(storm_index)) then ; then storm exists in both NR and CTL
print("Analyzing storm: " + g5nr_name(i) + "...")

; find first instance of matching dates between starting NR track and starting GFS run date, if exists
fcst_dates = get_unique_values(ctl_ymd(storm_index,:)) ; fcst init dates from GFS
fcst_runs = count_unique_values(ctl_ymd(storm_index,:)) ; number of GFS runs for this particular TC storm

fcst_hours = get_unique_values(ctl_hr(storm_index,:))
print(fcst_hours)
exit

do j=0,fcst_runs-1 ; loop over forecast runs

check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j)))

if (.not.ismissing(check)) then
start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR 
start_index = check ; index start of nature run
print("Analyzing GFS run " + start_date + "...")

; loop over time periods of NR and GFS, computing and saving storm track and intensity error
fcst_indices = ind(str_squeeze(ctl_ymd(storm_index,:)).eq.start_date) ; find all indices where the GFS run is initialized at the NR starting date 

if (.not.all(ismissing(fcst_indices))) then ; loop over fcsts of GFS, computing storm track and intensity error from NR

; conform length of NR tracks or GFS run to each other to account for runs to persist longer or shorter than the best track (if GFS forecast dies/goes beyond best track)
length_gfs = dimsizes(fcst_indices) ; number of forecast hours in this GFS run
length_nr  = num(.not.ismissing(g5nr_ymd(i,start_index::))) ; number of dates in NR from start of GFS run to the end of best track dataset for this storm


if (length_gfs.gt.length_nr) then ; means there are more gfs runs than best tracks
start_ind_nr = start_index
start_ind_gfs = fcst_indices(0)
end_ind_nr = num(.not.ismissing(g5nr_ymd(i,::)))-1 ; last date index of G5NR best track file
end_ind_gfs = (end_ind_nr - start_ind_nr) + fcst_indices(0) ; last date index of G5NR best track file
end if


if (length_nr.ge.length_gfs) then ; means there are more nr tracks than gfs runs
start_ind_nr = start_index
start_ind_gfs = fcst_indices(0)
end_ind_gfs = fcst_indices(num(.not.ismissing(fcst_indices))-1) ; last date index of GFS track file
end_ind_nr = (end_ind_gfs - start_ind_gfs) + start_index ; last date index of G5NR best track file
end if


; track and intensity error
end_2 = end_ind_gfs - start_ind_gfs

track_error_ctl(i,j,0:end_2) = abs(gc_latlon(ctl_lat(storm_index,start_ind_gfs:end_ind_gfs),ctl_lon(storm_index,start_ind_gfs:end_ind_gfs),g5nr_lat(i,start_ind_nr:end_ind_nr),g5nr_lon(i,start_ind_nr:end_ind_nr),2,4)) ; returns km distance
intensity_error_ctl(i,j,0:end_2) = abs( ctl_pres(storm_index,start_ind_gfs:end_ind_gfs) - g5nr_pres(i,start_ind_nr:end_ind_nr) )

ctl_lat_array(i,j,0:end_2) = ctl_lat(storm_index,start_ind_gfs:end_ind_gfs)
ctl_lon_array(i,j,0:end_2) = ctl_lon(storm_index,start_ind_gfs:end_ind_gfs)

delete(fcst_indices)

end if ; end if for checking missing values of fcst_indices


else ; then no matching date with NR and GFS for this init. cycle

print("No matching NR date with GFS Init Run of " + fcst_dates(j) + "...")

end if ; end if for checking matching dates

end do ; end loop over GFS forecast runs


delete(fcst_dates)
delete(fcst_runs)

else ; else for checking if both storms exist in NR and CTL

print("Skipping storm: " + g5nr_name(i) + "...")



end if ; end if for checking if both storms exist in NR and CTL



end do ; end loop over computing stats over each storm of NR

count_storms = 0
;--------------------------------------------------------------------------------------------------------
print("Computing track and intensity errors for EXP...")
; now for EXP only since the dates between CTL and EXP are not always the same so treat separately
do i=0,dimsizes(g5nr_name)-1 

; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm

storm_index = ind(exp_name.eq.g5nr_name(i))

if (.not.ismissing(storm_index)) then ; then storm exists in both NR and EXP
count_storms = count_storms + 1
print("Analyzing storm: " + g5nr_name(i) + "...")

; find first instance of matching dates between starting NR track and starting GFS run date, if exists
fcst_dates = get_unique_values(exp_ymd(storm_index,:)) ; fcst init dates from GFS
fcst_runs = count_unique_values(exp_ymd(storm_index,:)) ; number of GFS runs for this particular TC storm


do j=0,fcst_runs-1 ; loop over forecast runs

check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j)))

if (.not.ismissing(check)) then
start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR 
start_index = check ; index start of nature run
print("Analyzing GFS run " + start_date + "...")

; loop over time periods of NR and GFS, computing and saving storm track and intensity error
fcst_indices = ind(str_squeeze(exp_ymd(storm_index,:)).eq.start_date) ; find all indices where the GFS run is initialized at the NR starting date 

if (.not.all(ismissing(fcst_indices))) then ; loop over fcsts of GFS, computing storm track and intensity error from NR

; conform length of NR tracks or GFS run to each other to account for runs to persist longer or shorter than the best track (if GFS forecast dies/goes beyond best track)
length_gfs = dimsizes(fcst_indices) ; number of forecast hours in this GFS run
length_nr  = num(.not.ismissing(g5nr_ymd(i,start_index::))) ; number of dates in NR from start of GFS run to the end of best track dataset for this storm


if (length_gfs.gt.length_nr) then ; means there are more gfs runs than best tracks
start_ind_nr = start_index
start_ind_gfs = fcst_indices(0)
end_ind_nr = num(.not.ismissing(g5nr_ymd(i,::)))-1 ; last date index of G5NR best track file
end_ind_gfs = (end_ind_nr - start_ind_nr) + fcst_indices(0) ; last date index of G5NR best track file
end if


if (length_nr.ge.length_gfs) then ; means there are more nr tracks than gfs runs
start_ind_nr = start_index
start_ind_gfs = fcst_indices(0)
end_ind_gfs = fcst_indices(num(.not.ismissing(fcst_indices))-1) ; last date index of GFS track file
end_ind_nr = (end_ind_gfs - start_ind_gfs) + start_index ; last date index of G5NR best track file
end if


; track and intensity error
end_2 = end_ind_gfs - start_ind_gfs

track_error_exp(i,j,0:end_2) = abs(gc_latlon(exp_lat(storm_index,start_ind_gfs:end_ind_gfs),exp_lon(storm_index,start_ind_gfs:end_ind_gfs),g5nr_lat(i,start_ind_nr:end_ind_nr),g5nr_lon(i,start_ind_nr:end_ind_nr),2,4)) ; returns km distance
intensity_error_exp(i,j,0:end_2) = abs( exp_pres(storm_index,start_ind_gfs:end_ind_gfs) - g5nr_pres(i,start_ind_nr:end_ind_nr) )

exp_lat_array(i,j,0:end_2) = exp_lat(storm_index,start_ind_gfs:end_ind_gfs)
exp_lon_array(i,j,0:end_2) = exp_lon(storm_index,start_ind_gfs:end_ind_gfs)

delete(fcst_indices)

end if ; end if for checking missing values of fcst_indices


else ; then no matching date with NR and GFS for this init. cycle

print("No matching NR date with GFS Init Run of " + fcst_dates(j) + "...")

end if ; end if for checking matching dates

end do ; end loop over GFS forecast runs


delete(fcst_dates)
delete(fcst_runs)

else ; else for checking if both storms exist in NR and EXP

print("Skipping storm: " + g5nr_name(i) + "...")



end if ; end if for checking if both storms exist in NR and EXP



end do ; end loop over computing stats over each storm of NR


fcst_hr = ispan(0,384,6)

;;;;;;;;;;;;;;;;;;;;;;;;
;; plots results for each storm and GFS forecast separately
;;;;;;;;;;;;;;;;;;;;;;;;


print("Plotting results as function of Storm and GFS run separately...")
; now for EXP only since the dates between CTL and EXP are not always the same so treat separately
do i=0,dimsizes(g5nr_name)-1 

; first make sure the G5NR name is also found in the CTL experiment, if not, skip the storm

storm_index = ind(exp_name.eq.g5nr_name(i))

if (.not.ismissing(storm_index)) then ; then storm exists in both NR and EXP

print("plotting storm: " + g5nr_name(i) + "...")

; find first instance of matching dates between starting NR track and starting GFS run date, if exists
fcst_dates = get_unique_values(exp_ymd(storm_index,:)) ; fcst init dates from GFS
fcst_runs = count_unique_values(exp_ymd(storm_index,:)) ; number of GFS runs for this particular TC storm
 

do j=0,fcst_runs-1 ; loop over forecast runs

check = ind(str_squeeze(g5nr_ymd(i,:)).eq.str_squeeze(fcst_dates(j)))

if (.not.ismissing(check)) then
start_date = str_squeeze(g5nr_ymd(i,check)) ; starting date for NR 
start_index = check ; index start of nature run
print("Plotting GFS run " + start_date + "...")

if (.not.all(ismissing(track_error_ctl(storm_index,j,:))) .and. .not.all(ismissing(track_error_ctl(storm_index,j,:))) ) then ; then make track and intensity error plots

if (num(.not.ismissing(track_error_ctl(storm_index,j,:))).gt.1) then ; need 2 elements to plot 


; count number of non-missing values, conform each to the smaller array
count_ctl = num(.not.ismissing(track_error_ctl(storm_index,j,:)))
count_exp = num(.not.ismissing(track_error_exp(storm_index,j,:)))

if (count_ctl.le.count_exp) then

x_axis_max = count_ctl -1 ; subtract one based on array indexing

end if

if (count_exp.lt.count_ctl) then

x_axis_max = count_exp -1 ; subtract one based on array indexing

end if

max1 = max(track_error_ctl(storm_index,j,:))
max2 = max(track_error_exp(storm_index,j,:))
maxes = (/max1,max2/)
max_track = max(maxes)

max1 = max(intensity_error_ctl(storm_index,j,:))
max2 = max(intensity_error_exp(storm_index,j,:))
maxes = (/max1,max2/)
max_intensity = max(maxes)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;     
; plot errors as a function of forecast lead time
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


wks   = gsn_open_wks ("png","track_error_fcst_hour-" + g5nr_name(i) + "_" + start_date)

 res                   = True                      ; plot mods desired
 res@tiMainString      = "Track Error (" + g5nr_name + "): " + start_date
 res@trYMaxF = max_track + 100.0
 res@trYMinF = 0               
 res@gsnFrame          = False                      ; don't advance frame yet
 colors_plot           = (/"Black","Red"/)           ; change line color
 res@gsnMaximize = True
 res@gsnDraw = False
 
  res@vpWidthF         = 0.8                      ; set width and height
  res@vpHeightF        = 0.3
  res@vpXF             = 0.15
  res@vpYF             = 0.9

 res@trXMinF = 0
 res@trXMaxF = x_axis_max

 x_axis = ispan(0,x_axis_max,1)
 res@tmXBMode        = "Explicit"                     ; Force labels where we want them, later as function of forecast hour
 res@tmXBValues = ispan(0,x_axis_max,2)
 res@tmXBLabels      = fcst_hr(0:x_axis_max:2)
 res@tiYAxisString   = "Track Error (km)"
 res@tiXAxisString   = "Forecast lead time (hr)"
 res@tiXAxisFontHeightF = 0.015
 res@tiYAxisFontHeightF = 0.015
 res@tmXBLabelFontHeightF = 0.012
 res@tmYLLabelFontHeightF = 0.015
 res@xyLineThicknessF = 4.0
 res@tiMainFontHeightF = 0.021
 res@xyLineColor = colors_plot(0)

 plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl(storm_index,j,0:x_axis_max),res)

 res@xyLineColor = colors_plot(1)
 res@xyLineThicknessF = 4.0
 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp(storm_index,j,0:x_axis_max),res)


overlay(plot0,plot1)

xtxt2 = (/0.19,0.19/)
ytxt2 = (/0.64,0.62/)
labels2 = (/"CTL","C2REF"/)

 ; Attach a legend
 txres               = True
 txres@txFontHeightF = 0.016

  do k=0,dimsizes(labels2)-1
  txres@txFontColor = colors_plot(k)
  gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres)
  end do
  
  draw(plot0)
  frame(wks)
  
  delete(res)
  
wks   = gsn_open_wks ("png","intensity_error_fcst_hour-" + g5nr_name(i) + "_" + start_date)


 res                   = True                      ; plot mods desired
 res@tiMainString      = "Intensity Error (" + g5nr_name + "): " + start_date
 res@trYMaxF = max_intensity + 5
 res@trYMinF = 0               
 res@gsnFrame          = False                      ; don't advance frame yet
 colors_plot           = (/"Black","Red"/)           ; change line color
 res@gsnMaximize = True
 res@gsnDraw = False
 
  res@vpWidthF         = 0.8                      ; set width and height
  res@vpHeightF        = 0.3
  res@vpXF             = 0.15
  res@vpYF             = 0.9

 res@trXMinF = 0
 res@trXMaxF = x_axis_max

 x_axis = ispan(0,x_axis_max,1)
 res@tmXBMode        = "Explicit"                     ; Force labels where we want them, later as function of forecast hour
 res@tmXBValues = ispan(0,x_axis_max,2)
 res@tmXBLabels      = fcst_hr(0:x_axis_max:2)
 res@tiYAxisString   = "Intensity Error (hPa)"
 res@tiXAxisString   = "Forecast lead time (hr)"
 res@tiXAxisFontHeightF = 0.015
 res@tiYAxisFontHeightF = 0.015
 res@tmXBLabelFontHeightF = 0.012
 res@tmYLLabelFontHeightF = 0.015
 res@xyLineThicknessF = 4.0
 res@tiMainFontHeightF = 0.021
 res@xyLineColor = colors_plot(0)

 plot0 = gsn_csm_xy(wks,x_axis,intensity_error_ctl(storm_index,j,0:x_axis_max),res)

 res@xyLineColor = colors_plot(1)
 res@xyLineThicknessF = 4.0
 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp(storm_index,j,0:x_axis_max),res)


overlay(plot0,plot1)

xtxt2 = (/0.19,0.19/)
ytxt2 = (/0.64,0.62/)
labels2 = (/"CTL","C2REF"/)

 ; Attach a legend
 txres               = True
 txres@txFontHeightF = 0.016

  do k=0,dimsizes(labels2)-1
  txres@txFontColor = colors_plot(k)
  gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres)
  end do
  
  draw(plot0)
  frame(wks)
  
  delete(res)
  delete(x_axis)
  
  
  ;;;; storm track map result
  
  ; determine bounds of track map
  maxes_lon = max(g5nr_lon(i,:))
  maxes_lat = max(g5nr_lat(i,:))
  
  mins_lon = min(g5nr_lon(i,:))
  mins_lat = min(g5nr_lat(i,:))
  
  wks=gsn_open_wks("png","track_map_" + g5nr_name(i) + "_" + start_date)

 res = True
 res@gsnMaximize = True   ; maximize plot in frame
 res@gsnFrame = False     ; don't advance frame
 res@gsnDraw = False      ; don't draw plot right away, wait till later
 
 res@mpFillOn = False
 res@mpPerimOn = True
 res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
 res@pmTickMarkDisplayMode = "Always" 
 res@mpLimitMode = "LatLon"        ; select subregion
 res@mpMinLatF   = mins_lat - 10.0
 res@mpMaxLatF   = maxes_lat + 10.0           
 res@mpMinLonF   = mins_lon - 10.0
 res@mpMaxLonF   = maxes_lon + 10.0
 

 res@tmYROn      = False     ; turn off right and top tickmarks
 res@tmXTOn      = False

 res@tiMainString      = " " + g5nr_name + ": " + start_date
 res@tiMainFontHeightF = 0.02

 map = gsn_csm_map_ce(wks,res)                  ; create map

 ; Set up some legend resources.
   lgres                    = True
   lgres@lgLineColors       = (/"black","red","green"/)
   lgres@lgLineThicknessF   = 2. 
   lgres@lgLabelFontHeightF = .08            ; set the legend label font thickness
   lgres@vpWidthF           = 0.15           ; width of legend (NDC)
   lgres@vpHeightF          = 0.1            ; height of legend (NDC)
   lgres@lgMonoDashIndex    = True   
   lgres@lgPerimColor       = "orange"       ; draw the box perimeter in orange
   lgres@lgPerimThicknessF = 2.0            ; thicken the box perimeter
   labels = (/"G5NR","CTL","C2REF"/)

 ; Create the legend.
   lbid   = gsn_create_legend(wks,3,labels,lgres)         ; create legend

; Set up resources to attach legend to map.
   amres = True
   amres@amParallelPosF   =  -0.35 ;0.27 	         ; positive move legend to the right
   amres@amOrthogonalPosF = -0.35                 ; positive move the legend down
   annoid1 = gsn_add_annotation(map,lbid,amres)   ; attach legend to plot

 ; Add trajectory lines.
   pres                  = True               ; polyline resources
   pres@gsLineThicknessF = 3.0                ; line thickness
   pres@gsLineColor      = "black"
   line1 = gsn_add_polyline(wks,map,g5nr_lon(i,start_index::),g5nr_lat(i,start_index::),pres)      ; draw the traj

   pres                  = True               ; polyline resources
   pres@gsLineColor      = "red"
   line2 = gsn_add_polyline(wks,map,ctl_lon_array(storm_index,j,:),ctl_lat_array(storm_index,j,:),pres)      ; draw the traj

   pres                  = True               ; polyline resources
   pres@gsLineColor      = "green"
   line3 = gsn_add_polyline(wks,map,exp_lon_array(storm_index,j,:),exp_lat_array(storm_index,j,:),pres)      ; draw the traj


; Add markers to the trajectories.
   mres                = True         ; marker resources for best track
   mres@gsMarkerIndex  = 16           ; marker style (filled circle)
   mres@gsMarkerSizeF  = 8.0          ; marker size
   mres@gsMarkerColor  = "black"      ; maker color
   markers1 = gsn_add_polymarker(wks,map,g5nr_lon(i,start_index::),g5nr_lat(i,start_index::),mres)
   mres@gsMarkerColor  = "red"      ; maker color
   markers2 = gsn_add_polymarker(wks,map,ctl_lon_array(storm_index,j,:),ctl_lat_array(storm_index,j,:),mres)
   mres@gsMarkerColor  = "green"      ; maker color
   markers3 = gsn_add_polymarker(wks,map,exp_lon_array(storm_index,j,:),exp_lat_array(storm_index,j,:),mres)


   draw(map)                                           
   frame(wks)    
   delete(wks)
  delete(res)
  delete(mres)
  delete(pres)
  delete(lgres)
  
  delete(maxes_lon)
  delete(maxes_lat)
  delete(mins_lon)
  delete(mins_lat)
  
end if

end if ; end if for checking missing values

end if

end do ; end loop over GFS forecast runs


delete(fcst_dates)
delete(fcst_runs)

else ; else for checking if both storms exist in NR and EXP

print("Skipping storm: " + g5nr_name(i) + "...")



end if ; end if for checking if both storms exist in NR and EXP



end do ; end loop over computing stats over each storm of NR




  
  
  
  
  

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


; average the track and intensity errors for both experiments
print("Computing average track and intensity error statistics...")

track_error_ctl_st_fh = dim_avg_n(track_error_ctl,1) ; average as function of storm and forecast hour
track_error_exp_st_fh = dim_avg_n(track_error_exp,1)
intensity_error_ctl_st_fh = dim_avg_n(intensity_error_ctl,1)
intensity_error_exp_st_fh = dim_avg_n(intensity_error_exp,1)

track_error_ctl_fh = dim_avg_n(track_error_ctl,(/0,1/)) ; average as function of forecast hour only
track_error_exp_fh = dim_avg_n(track_error_exp,(/0,1/))

track_error_ctl_fh_var = dim_variance_n(track_error_ctl,(/0,1/)) ; average as function of forecast hour only
track_error_exp_fh_var = dim_variance_n(track_error_exp,(/0,1/))

intensity_error_ctl_fh = dim_avg_n(intensity_error_ctl,(/0,1/))
intensity_error_exp_fh = dim_avg_n(intensity_error_exp,(/0,1/))

; total average over all storms and hours, one value here

track_error_ctl_avg = avg(track_error_ctl) ; average as function of forecast hour only
track_error_exp_avg = avg(track_error_exp)
intensity_error_ctl_avg = avg(intensity_error_ctl)
intensity_error_exp_avg = avg(intensity_error_exp)

print("Total storms analyzed: " + count_storms)

; determine number of non-missing values at each forecast hour in each experiment


do i=0,64 ; loop over forecast hours (output every 6 hours)
    fcst_hr_array_ctl(i) = num(.not.ismissing(track_error_ctl(:,:,i)))
end do

do i=0,64 ; loop over forecast hours (output every 6 hours)
    fcst_hr_array_exp(i) = num(.not.ismissing(track_error_exp(:,:,i)))
end do

; overlay the number of samples in each forecast hour
fcst_hr_array_avg = round(((fcst_hr_array_ctl + fcst_hr_array_exp) / 2.0),3) ; average number between the two experiments

; compute significance over full avg and over forecast hour
siglvl = 0.05
ave_ctl = track_error_ctl_avg
ave_exp = track_error_exp_avg

var_ctl = variance(track_error_ctl)
var_exp = variance(track_error_exp)

s_ctl = num(.not.ismissing(track_error_ctl))
s_exp = num(.not.ismissing(track_error_exp))

iflag = False
tval_opt = False
prob_avg = ttest(ave_ctl,var_ctl,s_ctl,ave_exp,var_exp,s_exp,iflag,tval_opt)

delete([/s_ctl,s_exp,ave_ctl,ave_exp,var_ctl,var_exp/])

ave_ctl = track_error_ctl_fh(0:50)
ave_exp = track_error_exp_fh(0:50)

var_ctl = track_error_ctl_fh_var(0:50)
var_exp = track_error_exp_fh_var(0:50)

s_ctl = fcst_hr_array_ctl(0:50)
s_exp = fcst_hr_array_exp(0:50)

iflag = False
tval_opt = False
prob_fh = ttest(ave_ctl,var_ctl,s_ctl,ave_exp,var_exp,s_exp,iflag,tval_opt)

; plot results

print("Plotting results...")


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Storm track errors as histogram, along with wind and intensity errors
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

colors_plot = (/"red","blue"/)
names_short_1 = (/"A","B","","A","B"/) ; names for legend
names_short_2 = (/"A = CTL","B = C2REF"/) ; names for legend
xtxt = (/0.25,0.25/)
ytxt = (/0.5,0.47/)

print("Making track, wind, and intensity error plots as histograms...")

x_axis_1 = ispan(2,4,2)
x_axis_2 = ispan(8,10,2)
y_axis_values_1 = (/track_error_ctl_avg,track_error_exp_avg/) ; y values for the bar charts 
y_axis_values_2 = (/intensity_error_ctl_avg,intensity_error_exp_avg/)

max_y1=max(y_axis_values_1)
max_y2=max(y_axis_values_2)

wks = gsn_open_wks(output_type,"track-intensity-error-all-storms") ; create workstation
 res = True
 res@gsnFrame = False
 res@gsnDraw = False
 res@tiMainString  = "Track and Intensity Error (" + count_storms + " TCs)"
 res@gsnXYBarChart = True            ; Create bar plot
 ;res@tiMainOffsetYF = 0.006
 res@tiMainFontHeightF = 0.02
 res@tmXTMode       = "Explicit" ; Define own tick mark labels.
 res@tmXTValues     = x_axis_1
 res@tmXTOn         = False             ; turn on the top tickmarks
 res@tmXBOn         = True
 res@tmXBLabelsOn = True
 res@tmXBLabelJust = "BottomCenter"
 res@gsnScale = True
 res@tmYROn         = False
 res@tmYRMinorOn    = True
 res@tmYLOn         = True
 res@tmYLMinorOn    = True
 res@gsnXYBarChartOutlineThicknessF = 1.5
 res@tmXTMinorOn    = False
 ;res@tiYAxisOffsetXF = 0.035
 res@tiYAxisString = "Track Error (km)"
 res@tiYAxisFontHeightF = 0.018
 res@trYMaxF         = 400   ; max value on y-axis
 res@trYMinF         = 300
 res@gsnXYBarChartBarWidth = 1.5                 ; change bar widths 0.75 originally
 res@gsnXYBarChartColors = colors_plot(0)
 res@tmXBMode                = "Explicit"
 res@trXMinF               = 0                    ; adds space on either end
 res@trXMaxF               = 12                    ; of the 1st and last bars
 
 res@tmXBValues               = (/2,4/)
 res@tmXBLabels               = names_short_1(0:1)
 res@tmXBLabelJust            = "BottomCenter"

plot = gsn_csm_xy(wks,x_axis_1,y_axis_values_1,res)

bres = True
bres@tmXBMode                = "Explicit"
bres@gsnFrame = False
bres@gsnDraw = False
bres@gsnXYBarChart = True
bres@gsnXYBarChartColors = colors_plot(1)
bres@tmYROn = True
bres@tmYRMinorOn = True
bres@tmYLOn = False
bres@trYMinF = 20
bres@tmXBValues = (/8,10/)
bres@tiMainString = ""
bres@trYMaxF = 33.0 ; hPa
bres@tmYRLabelsOn = True
bres@tiYAxisFontHeightF = 0.018
bres@gsnXYBarChartOutlineThicknessF = 1.5
bres@tmYROn = True
bres@gsnXYBarChartBarWidth = 1.5                 ; change bar widths 0.75 originally
bres@tmXBOn         = True
bres@trXMinF               = 0                    ; adds space on either end
bres@trXMaxF               = 12                    ; of the 1st and last bars
bres@tiYAxisString = "Intensity Error (hPa)"
bres@tiYAxisSide         = "Right"
bres@tiYAxisAngleF       = 270
bres@tmXBLabels = names_short_1(0:1)
bres@tmXBLabelJust            = "BottomCenter"
bres@tmXTOn = False
draw(plot)

;**********************************************************
; add text labels
;**********************************************************
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.016                ; default size is HUGE!
  txres@txAngleF            = 0                  ; text angle
  txres@txJust              = "CenterCenter"         ; puts text on top of bars;
gsn_text(wks,plot,names_short_2(0),2.4,390,txres) ; add labels
gsn_text(wks,plot,names_short_2(1),2.54,387,txres) ; add labels

 plot2 = gsn_csm_xy(wks,x_axis_2,y_axis_values_2,bres)

;---Create a blank plot in order to get additional X tickmark labels
  vres = True
  vres@tmXBMode             = "Explicit"
  vres@tmXBValues           = (/3.1,9.1/)
  vres@tmXBLabels           = (/"Position Error","Magnitude Error"/)
  vres@tmXBLabelDeltaF      = 0.5
  vres@tmXBLabelFontHeightF = 0.018
  vres@tmYLOn = False
  vres@tmYROn = False
; vres@tmXBOn = False
  vres@tmXTOn = False
  blank_plot = gsn_csm_blank_plot(wks,vres)
  overlay(plot2,blank_plot)

draw(plot2)
frame(wks)

delete(bres)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;     
; plot errors as a function of forecast lead time
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


delete(res)
delete(txres)
delete(colors_plot)

wks   = gsn_open_wks ("png","track_error_fcst_hour-all-storms")

 res                   = True                      ; plot mods desired
 res@gsnXYBarChart = False
 res@tiMainString      = "Track Error (" + count_storms + " TCs)"
 res@trYMaxF = 1800
 res@trYMinF = 0               
 res@gsnFrame          = False                      ; don't advance frame yet
 colors_plot           = (/"Black","Red"/)           ; change line color
 res@gsnDraw = False
 res@gsnFrame = False

 res@gsnScale = True
 res@trXMinF = 0
 res@trXMaxF = 50

 x_axis = ispan(0,50,1)
 res@tmXBMode        = "Explicit"                     ; Force labels where we want them, later as function of forecast hour
 res@tmXBValues = ispan(0,50,4)
 res@tmXBLabels      = ispan(0,300,24)
 res@tiYAxisString   = "Track Error (km)"
 res@tiXAxisString   = "Forecast lead time (hr)"
 res@tiXAxisFontHeightF = 0.015
 res@tmYROn = False
 res@tmYROn = False
 res@tmYRMinorOn = False
 res@tiYAxisFontHeightF = 0.015
 res@tmXBLabelFontHeightF = 0.015
 res@tmYLLabelFontHeightF = 0.015
 res@xyLineThicknessF = 4.0
 res@tiMainFontHeightF = 0.021
 res@xyLineColor = colors_plot(0)

 plot0 = gsn_csm_xy(wks,x_axis,track_error_ctl_fh(0:50),res)

 res@xyLineColor = colors_plot(1)
 res@xyLineThicknessF = 4.0
 plot1 = gsn_csm_xy(wks,x_axis,track_error_exp_fh(0:50),res)


bres = True
bres@tmXBMode                = "Explicit"
bres@gsnFrame = False
bres@gsnDraw = False
bres@gsnXYBarChart = True
bres@gsnXYBarChartColors = "grey78"
bres@tmYROn = True
bres@tmYRMinorOn = True
bres@tmYLOn = False
bres@trYMinF = 0
bres@tiMainString = ""
bres@trYMaxF = 150.0 ; number of samples in each forecast hour
bres@tmYRLabelsOn = True
bres@tiYAxisFontHeightF = 0.015
bres@gsnXYBarChartOutlineThicknessF = 1.5
bres@tmYROn = True
bres@gsnXYBarChartBarWidth = 0.7                 ; change bar widths 0.75 originally
bres@tmXBOn         = False
bres@trXMinF               = 0                    ; adds space on either end
bres@trXMaxF               = 50                    ; of the 1st and last bars
bres@tiYAxisString = "# of Samples"
bres@tiYAxisSide         = "Right"
bres@tiYAxisAngleF       = 270
bres@tmXTOn = False


plot2 = gsn_csm_xy(wks,x_axis,fcst_hr_array_avg(0:50),bres)


xtxt2 = (/0.3,0.3/)
ytxt2 = (/0.75,0.72/)
labels2 = (/"CTL","C2REF"/)


overlay(plot0,plot1)

draw(plot2)
  
  ; Attach a legend
 txres               = True
 txres@txFontHeightF = 0.016

  do k=0,dimsizes(labels2)-1
  txres@txFontColor = colors_plot(k)
  gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres)
  end do
  
  draw(plot0)
  
  frame(wks)
  
  delete(res)
  delete(bres)

  
wks   = gsn_open_wks ("png","intensity_error_fcst_hour-all-storms")

 res                   = True                      ; plot mods desired
 res@gsnXYBarChart = False
 res@tiMainString      = "Intensity Error (" + count_storms + " TCs)"
 res@trYMaxF = 40
 res@trYMinF = 0               
 res@gsnFrame          = False                      ; don't advance frame yet
 colors_plot           = (/"Black","Red"/)           ; change line color
 res@gsnDraw = False
 res@gsnFrame = False

 res@gsnScale = True
 res@trXMinF = 0
 res@trXMaxF = 50

 x_axis = ispan(0,50,1)
 res@tmXBMode        = "Explicit"                     ; Force labels where we want them, later as function of forecast hour
 res@tmXBValues = ispan(0,50,4)
 res@tmXBLabels      = ispan(0,300,24)
 res@tiYAxisString   = "Intensity Error (hPa)"
 res@tiXAxisString   = "Forecast lead time (hr)"
 res@tiXAxisFontHeightF = 0.015
 res@tmYROn = False
 res@tmYROn = False
 res@tmYRMinorOn = False
 res@tiYAxisFontHeightF = 0.015
 res@tmXBLabelFontHeightF = 0.015
 res@tmYLLabelFontHeightF = 0.015
 res@xyLineThicknessF = 4.0
 res@tiMainFontHeightF = 0.021
 res@xyLineColor = colors_plot(0)

 plot0 = gsn_csm_xy(wks,x_axis,intensity_error_ctl_fh(0:50),res)

 res@xyLineColor = colors_plot(1)
 res@xyLineThicknessF = 4.0
 plot1 = gsn_csm_xy(wks,x_axis,intensity_error_exp_fh(0:50),res)


bres = True
bres@tmXBMode                = "Explicit"
bres@gsnFrame = False
bres@gsnDraw = False
bres@gsnXYBarChart = True
bres@gsnXYBarChartColors = "grey78"
bres@tmYROn = True
bres@tmYRMinorOn = True
bres@tmYLOn = False
bres@trYMinF = 0
bres@tiMainString = ""
bres@trYMaxF = 150.0 ; number of samples in each forecast hour
bres@tmYRLabelsOn = True
bres@tiYAxisFontHeightF = 0.015
bres@gsnXYBarChartOutlineThicknessF = 1.5
bres@tmYROn = True
bres@gsnXYBarChartBarWidth = 0.7                 ; change bar widths 0.75 originally
bres@tmXBOn         = False
bres@trXMinF               = 0                    ; adds space on either end
bres@trXMaxF               = 50                    ; of the 1st and last bars
bres@tiYAxisString = "# of Samples"
bres@tiYAxisSide         = "Right"
bres@tiYAxisAngleF       = 270
bres@tmXTOn = False


plot2 = gsn_csm_xy(wks,x_axis,fcst_hr_array_avg(0:50),bres)


xtxt2 = (/0.3,0.3/)
ytxt2 = (/0.75,0.72/)
labels2 = (/"CTL","C2REF"/)


overlay(plot0,plot1)

draw(plot2)
  
  ; Attach a legend
 txres               = True
 txres@txFontHeightF = 0.016

  do k=0,dimsizes(labels2)-1
  txres@txFontColor = colors_plot(k)
  gsn_text_ndc(wks,labels2(k),xtxt2(k),ytxt2(k),txres)
  end do
  
  draw(plot0)
  
  frame(wks)
  
  delete(res)
  delete(bres)





print("-------Track Program completed--------")
end