;********************************************************************
; This script reads netcdf files containing AI mask and inflates
; the mask to include surrounding non-cloudy areas and writes output
; to separate netcdf files for 00 time step and 30 timestep
;********************************************************************


files = systemfunc("ls /discover/nobackup/aakbari1/TC/DiffusionMaskV1.1/mask_202009*.nc4")
dim1 = dimsizes(files)

print(typeof(files))

filname = stringtochar(files(0))
;print(filname)
;print(dimsizes(filname))

yr = chartostring(filname(54:57))
mnth = chartostring(filname(58:59))
day = chartostring(filname(60:61))
hr = chartostring(filname(62:63))

yr_float = tofloat(yr)
mnth_float = tofloat(mnth)
day_float = tofloat(day)
hr_float = tofloat(hr)

print(yr)
print(mnth)
print(day)
print(hr)

f = addfile(files(0), "r")


;g = f=>GROUPNAME
;print(g)
aimask = f->mask(0,:,:)
lat  = f->latitude
lon = f->longitude
time = f->time

print(dimsizes(aimask))
print(dimsizes(lat))
print(dimsizes(lon))
print(dimsizes(time))

printMinMax(lat,0)
printMinMax(lon,0)

lat@units = "degrees_north"
lon@units = "degrees_east"

aimask!0 = "lat"
aimask!1 = "lon"
aimask&lat = lat
aimask&lon = lon

dummyvar = new(dimsizes(aimask),"float")
dummyvar = 1

plotvar = mask(dummyvar,aimask,1)   ; Return the values of ts where oro=1

plotvar!0 = "lat"
plotvar!1 = "lon"
plotvar&lat = lat
plotvar&lon = lon

;*****************************************************
; Read TC location info
;*****************************************************
tcfile = readAsciiTable("/discover/nobackup/mganesha/AImask/2020/2020_HurricaneList.txt",8,"float",0)
print(dimsizes(tcfile))
;print(tcfile(:,0))
;print(tcfile(:,7))

year = tcfile(:,0)
mm = tcfile(:,1)
dd = tcfile(:,2)
hh = tcfile(:,3)
lattc = tcfile(:,4)
lontc = tcfile(:,5)

;print(year)
;print(lattc)
;print(lontc)

dim2 = dimsizes(lattc)
print(dim2)
tcname = asciiread("/discover/nobackup/mganesha/AImask/2020/2020_HurricaneList.txt",dim2,"string")

tclist = new(dim2,"string")
do i = 0,dim2-1
  tempname = stringtochar(tcname(i))
  dim3 = dimsizes(tempname)
  tclist(i) = chartostring(tempname(46:dim3-1))

  delete(tempname)
  delete(dim3)
end do

;print(tclist)

if (any(mnth_float.eq.mm.and.day_float.eq.dd)) then
  ind1 = ind(mnth_float.eq.mm.and.day_float.eq.dd)
  print(ind1)
  print("These are the TCs observed on "+mnth+"-"+day+"-"+yr)
;  print(tclist(ind1))

  ind2 = ind(hh(ind1).eq.hr_float) 

  print(tclist(ind1(ind2)))

  lat_tc = lattc(ind1(ind2))
  lon_tc = lontc(ind1(ind2))
end if
  
;*****************************************************
;plot the AI mask
;*****************************************************

wks = gsn_open_wks("eps","Map_with_TCs_"+mnth+"-"+day+"-"+yr+"_"+hr+"Z")    

  res             = True                              ; plot mods desired
  res@gsnFrame    = False                             ; don't advance yet
  res@gsnMaximize = True
  res@cnFillOn             = True     ; turn on color fill
  res@cnLinesOn            = False    ; turn off contour lines
  res@cnConstFEnableFill = True
  res@cnConstFLabelOn  = False

  res@mpMinLatF = -60
  res@mpMaxLatF = 60
  res@mpMinLonF = -180
  res@mpMaxLonF = 180
  res@tiMainString = "AI masks on "+mnth+"-"+day+"-"+yr+" at "+hr+"Z"

  res@lbLabelBarOn = False

  res@gsnDraw     = False
  res@gsnFrame    = False

;  map = gsn_csm_map(wks,res)                ; Create the map plot, don't draw
  plot = gsn_csm_contour_map(wks,plotvar,res)                ; Create the map plot, don't draw

;Add TC locations
  
  dum = new(dimsizes(ind2),graphic)
  pmres = True
  pmres@gsMarkerColor = "blue"
  pmres@gsMarkerIndex = 16
  pmres@gsMarkerSizeF = 0.01
  pmres@gsMarkerOpacityF = 0.3

  do i=0,dimsizes(ind2)-1

    dum(i) = gsn_add_polymarker(wks, plot, lon_tc(i), lat_tc(i), pmres)
  end do
 
  draw(plot)
  frame(wks)