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/csm/shea_util.ncl"
;************************************
begin
;Reading the control file

        fname = "mom-ctrl-slc.dat"
        lines = asciiread(fname,-1,"string")
        para = asciiread(fname,-1,"double")

        ip = (str_get_field(lines(1:),1,","))   ;Variable from first coloumn
        minlat = para(0)
        maxlat = para(1)
        minlon = para(2)
        maxlon = para(3)
        time = para(4)
        plot_width = para(5)
        plot_height = para(6)
print(maxlon)
        plot_time = (str_get_field(lines(1:),6,","))
print(ip)
print(plot_time)

aa = systemfunc("ls Forecast*.nc")
print(aa)
init_time = str_get_cols(aa, 8, 18)
print(init_time)
        yy1 = str_get_cols(init_time, 0, 3)
        mm1 = str_get_cols(init_time, 4,5)
        dd1 = str_get_cols(init_time, 6, 7)
        hr1 = str_get_cols(init_time, 8,9)
print(yy1)
print(mm1)
print(dd1)
print(hr1)

        yy2 = str_get_cols(plot_time, 0, 3)
        mm2 = str_get_cols(plot_time, 4,5)
        dd2 = str_get_cols(plot_time, 6, 7)
        hr2 = str_get_cols(plot_time, 8,9)
print(yy2)
print(mm2)
print(dd2)
print(hr2)


total_days = stringtointeger(dd2)-stringtointeger(dd1)
total_hrs = stringtointeger(hr2)-stringtointeger(hr1)

print(total_days)
print(total_hrs)
time_cell = ((24*total_days + total_hrs)/6)-1 ;subtract 1,since time starts from 6th hrs in mom data
print(time_cell)


;************************************
 in = addfile(aa,"r")
 var = in->$ip$

printVarSummary(var)

var1 = lonPivot(var,181)
printVarSummary(var1)

xlat=var1&YT_J
xlon=var1&XT_I

print(xlon)

minlocY = ind_nearest_coord (minlat, xlat, 0)        ;minlat cell number in the data file
maxlocY = ind_nearest_coord (maxlat, xlat, 0)        ;maxlat cell number in the data file
minlocX = ind_nearest_coord (minlon, xlon, 0)        ;minlong cell number in the data file
maxlocX = ind_nearest_coord (maxlon, xlon, 0)        ;maxlong cell number in the data file

print(minlocY)
print(maxlocY)
print(minlocX)
print(maxlocX)


;************************************
type = "png"
type@wkWidth = plot_width
type@wkHeight = plot_height
plt_name_1 = "test_1"
plt_name = "mom_slc_"+ip+"_"+time+"UTC"
wks = gsn_open_wks(type,plt_name_1)

;************************************

 res=True
 res@cnFillOn             = True     ; turn on color fill
 res@gsnAddCyclic = False
 res@mpFillOn             = True    ; gray continents
 res@mpOutlineOn          = True    ; continental outline
 res@mpGeophysicalLineThicknessF= 3      ; Increase the thickness of coastal boundaries
 res@gsnDraw              =False    ; do not draw picture
 res@gsnFrame             = False    ; do not advance frame
 res@gsnSpreadColors      = True     ; use full colormap
 res@lbLabelBarOn = True
 res@lbTitleOn = False                   ;Turn off the color bar text
 res@lbOrientation = "Vertical"          ;color bar orientation

 res@pmTickMarkDisplayMode= "Always"
 res@gsnLeftString = ""
 res@gsnRightString = ""
 res@gsnCenterString = ""+var@long_name+"("+var@units+")     Time :"+time+"UTC"
 plot = gsn_csm_contour_map(wks,var(time_cell,0,minlocY:maxlocY,minlocX:maxlocX),res)     ; create plot
 draw(plot)
 frame(wks)
;************************************

trim_plot= systemfunc("convert -trim "+plt_name_1+"."+type+" "+plt_name+"."+type)
rem_plot = systemfunc("rm " +plt_name_1+"."+type )


end