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 "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"

begin

; Open model level output file
export_text     = False
plot_map        = False
plot_timeseries = False

 File = addfile( "./out.nc","r")
 time= File->time
 gust=File->ws_10m
 mslp=File->slp
 topo=File->Z_sfc
 lat=File->lat
 lon=File->lon
 hour=File->hour

; mask over land points
; not really needed if mslp is used 
 mslp=mask(mslp, topo.ge.5 , False)

; create plot

ResC = True
ResC@gsnDraw=False
ResC@gsnFrame = False
ResC@mpLimitMode = "LatLon"
ResC@mpMaxLonF    = 17                      ; specify the plot domain
ResC@mpMinLonF    = 11                      ;                         
ResC@mpMinLatF    = 34                     ;
ResC@mpMaxLatF    = 39                     ;
ResC@mpDataBaseVersion  = "HighRes"
;ResC@gsnAddCyclic = False

 ResC@gsnMaximize   = True         ; Maximize plot in frame. 
 
   wks = gsn_open_wks("x11","track") 
   cmap = read_colormap_file("circular_2")
   res_lb                    = True
   res_lb@vpWidthF           = 0.50
   res_lb@vpHeightF          = 0.10
   res_lb@lbPerimOn          = False            ; Turn off perimeter.
   res_lb@lbOrientation      = "Horizontal"     ; Default is vertical.
   ;res_lb@lbLabelStride      = 10    
   res_lb@lbLabelAlignment   = "InteriorEdges"  ; Default is "BoxCenters".
   res_lb@lbFillColors       = cmap(2:,:)      ; Colors for boxes.
   res_lb@lbMonoFillPattern  = True             ; Fill them all solid.
   res_lb@lbLabelFontHeightF = 0.015 
 
 lres= True
 lres@gsLineThicknessF = 3.0

 mres=True
 mres@gsMarkerIndex = 16
 mres@gsMarkerSizeF = 0.005
 
 txres=True
 txres@txFontHeightF = 0.01
  
 ilat=new(dimsizes(time), integer) 
 ilon=new(dimsizes(time), integer) 
 lat_min=new(dimsizes(time), double)
 lon_min=new(dimsizes(time), double)
 mslp_min=new(dimsizes(time), double)
 gust_max=new(dimsizes(time), double)
 copy_VarAtts(gust,gust_max)

start=0

;radius of lon-lat to search max. gust
;depends only on the resolution
lon_radius=0.05
lat_radius=0.05

i=start
 dims = dimsizes(mslp(i,:,:))
 x1d = ndtooned(mslp(i,:,:))      ; convert 2D array to 1D for use in minind
 inds = ind_resolve(minind (x1d), dims)    ; convert 1D array back to 2D 
 ilat(i) = inds(0,0)        ; select the latitude index where the X array is at its' minimum  
 ilon(i) = inds(0,1)        ; select the longitude index where the X array is at its' minimum  
 lat_min(i) = lat(ilat(i),ilon(i))   ; insert the latitude index into the lat coordinate variable
 lon_min(i) = lon(ilat(i),ilon(i))   ; insert the longitude index into the lon coordinate variable
 mslp_min(i)= mslp(i,ilat(i),ilon(i))
; gust_max(i)=max(gust(i,0,{lat_min(i):lat_min(i)+lat_radius},{lon_min(i):lon_min(i)+lon_radius}))
 delete(dims)
 delete(x1d)
 delete(inds)

i=start+1
 dims = dimsizes(mslp(i,:,:))
 x1d = ndtooned(mslp(i,:,:))      ; convert 2D array to 1D for use in minind
 inds = ind_resolve(minind (x1d), dims)    ; convert 1D array back to 2D 
 ilat(i) = inds(0,0)        ; select the latitude index where the X array is at its' minimum  
 ilon(i) = inds(0,1)        ; select the longitude index where the X array is at its' minimum  
 lat_min(i) = lat(ilat(i),ilon(i))   ; insert the latitude index into the lat coordinate variable
 lon_min(i) = lon(ilat(i),ilon(i))   ; insert the longitude index into the lon coordinate variable
 mslp_min(i)= mslp(i,ilat(i),ilon(i))
 ;gust_max(i)=max(gust(i,0,{lat_min(i):lat_min(i)+lat_radius},{lon_min(i):lon_min(i)+lon_radius}))
 delete(dims)
 delete(x1d)
 delete(inds)

; box for the minimum search around the first guess
 box=4
;

do i=start+2,dimsizes(time)-1
 ilat_c=2*ilat(i-1)-ilat(i-2)
 ilon_c=2*ilon(i-1)-ilon(i-2)
 
 dims = dimsizes(mslp(i, ilat_c-box:ilat_c+box , ilon_c-box:ilon_c+box ))
 x1d = ndtooned(mslp(i,  ilat_c-box:ilat_c+box , ilon_c-box:ilon_c+box  ))      ; convert 2D array to 1D for use in minind
 inds = ind_resolve(minind (x1d), dims)    ; convert 1D array back to 2D 
 ilat(i) =  inds(0,0)-box+ilat_c       ; select the latitude index where the X array is at its' minimum  
 ilon(i) =  inds(0,1)-box+ilon_c        ; select the longitude index where the X array is at its' minimum
 lat_min(i) = lat(ilat(i),ilon(i))   ; insert the latitude index into the lat coordinate variable
 lon_min(i) = lon(ilat(i),ilon(i))   ; insert the longitude index into the lon coordinate variable
 mslp_min(i)= mslp(i,ilat(i),ilon(i))
; gust_max(i)=max(gust(i,0,{lat_min(i):lat_min(i)+lat_radius},{lon_min(i):lon_min(i)+lon_radius}))
 delete(dims)
 delete(x1d)
 delete(inds)
 delete(ilat_c)
 delete(ilon_c)
end do 

if plot_map then
wks = gsn_open_wks("x11","track") 
print("Plotting map")
map=gsn_csm_map(wks,ResC)

cnLevels=fspan(min(mslp_min),max(mslp_min),10)
do j=0,dimsizes(mslp_min)-2
   lres@gsLineColor=GetFillColor(cnLevels,cmap,avg( (/mslp_min(j),mslp_min(j+1)/)))
   gsn_polyline(wks,map,(/lon_min(j),lon_min(j+1)/),(/lat_min(j),lat_min(j+1)/),lres) 
end do
marker1 = gsn_add_polymarker(wks,map,lon_min,lat_min,mres)
text1 = gsn_add_text(wks, map,hour,lon_min-0.01, lat_min+0.1,txres)
gsn_labelbar_ndc(wks,dimsizes(cnLevels)+1,floor(cnLevels),0.30,0.30,res_lb)

draw(map)
frame(wks)
end if

if export_text then 
; write coordinates and mslp in a text file
alist=[/hour_acc,lon_min,lat_min,mslp_min,gust_max/]
write_table(diri+"_track.txt","w",alist,"%2i%10.8f%10.8f%8.4f%4.2f")
end if

draw(plot)
frame(wks2)

if plot_timeseries then
wks3 = gsn_open_wks("x11","mslp_wind") 
resxy=True
resxy@vpWidthF       = 0.9
resxy@vpHeightF      = 0.5
resxy@gsnDraw = False
resxy@gsnFrame = False
resxy@trXMinF = start
resxy@trXMaxF = limit
resxy@xyLineColors="black"
resxy@xyLineThicknessF = 4
resxy@tiYAxisFontHeightF = 0.02
resxy@gsnMaximize   = True
resxy2=resxy
resxy2@xyLineColors="red"

timeseries=gsn_csm_xy2(wks3, hour_acc, mslp_min, gust_max,resxy,resxy2)

draw(timeseries)
frame(wks3)
end if 
 

end