;================================================;
; Zooming Satellite Projection:
; -----------------------------
; Temperature contour over Maitri Station
; at South pole, Antarctica
;================================================;

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   

begin

;;;;;;;; Defining Station name, Station lat & long ;;;;;;;;;;;;;;

  name = "MAITRI"
  slat = -70.76583
  slon =  11.73583

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Reading the filename table                             ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  filenames = asciiread("/expr/saha/DATA/INPUTS/filename_grib.tbl",-1,"string")
  nhours = dimsizes(filenames)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Defining variable arrays                               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ;shour_arr  = new((/nhours/),string)
  ;syear_arr  = new((/nhours/),string)
  ;smonth_arr = new((/nhours/),string)
  ;sdate_arr  = new((/nhours/),string)

;************************************************
; read in data 
;************************************************

 do ihr = 0, nhours-1, 2

  in = addfile(filenames(ihr),"r")
  ;in = addfile(filenames,"r")

  ;Times = in->PRES_4_MWSL_10@initial_time
  ;tk2 = in->TMP_4_SFC_10(:,:)
  initime = in->TMP_P0_L1_GLL0@initial_time
  fcsttime = in->TMP_P0_L1_GLL0@forecast_time
  ;timeunits = in->TMP_P0_L1_GLL0@forecast_time_units 
  tk2 = in->TMP_P0_L1_GLL0(:,:)
  tc2 = tk2-273.15
  tc2@units = "DegC" 
  tc2@description = "Temp 2mtr"
 
  str = str_split(initime," ")
  date = str(0)
  ;str1 = str(1) 
  ;str2 = str_split(str1,"(")
  ;time = stringtoint(str_split(str2,")"))
  ;time = time+fcsttime

  ;print(date)
  ;print(time)

  ;print(initime)
  ;print(fcsttime)
  ;print(timeunits)
   
  ;xlat = in->lat_4
  xlat = in->lat_0
  ;xlat@units = "degrees_south"
  xlat@units = "degrees_north"

  ;xlon = in->lon_4
  xlon = in->lon_0
  ;xlon@units = "degrees_west"
  xlon@units = "degrees_east"

  tc2!0 = "xlat"
  tc2!1 = "xlon"

  tc2&xlat = xlat
  tc2&xlon = xlon

  ;printVarSummary(xlat)
  ;print(xlat)
  ;print(Times)

;************************************************
; create plot
;************************************************

  ;wks  = gsn_open_wks("ps","temp_maitri"+fcsttime)  ; open a ps file
  wks  = gsn_open_wks("x11","temp_maitri")  ; open a ps file
  gsn_define_colormap(wks,"testcmap")            ; choose colormap
  ;gsn_define_colormap(wks,"sunshine_9lev")            ; choose colormap

  res                       = True         ; plot mods desired
  res@gsnMaximize           = True
  res@gsnDraw               = False
  res@gsnFrame              = False
  res@gsnAddCyclic          = False

  res@cnFillOn              = True         ; turn on color
  res@cnLinesOn             = False        ; turn off contour lines
  ;res@cnLevelSpacingF       = 0.5          ; contour spacing
  ;res@cnFillDrawOrder       = "PreDraw"    ; fill first

  ;res@lbLabelStride         = 6            ; skip every other label
  res@lbBoxLinesOn          = False        ; turn off box between colors

  res@gsnSpreadColors       = True         ; use full colormap
  ;res@gsnSpreadColorStart   = 10           ; start at color 10
  ;res@gsnSpreadColorEnd     = 185          ; end at color 185
  ;res@gsnPolar              = "SH"         ; southern hemisphere

  ;res@mpProjection          = "Satellite"  ; choose map projection
  ;res@mpCenterLatF          = slat         ; choose center lat
  ;res@mpCenterLonF          = slon         ; choose center lon
  ;res@mpLimitMode           = "Angles"     ; required
  ;res@mpLeftAngleF          = -50.
  ;res@mpRightAngleF         = 60
  ;res@mpBottomAngleF        = -90 
  ;res@mpTopAngleF           = -60

  res@mpProjection          = "CylindricalEquidistant"  ; choose map projection
  res@mpCenterLatF          = slat         ; choose center lat
  res@mpCenterLonF          = slon         ; choose center lon
  res@mpLimitMode           = "LatLon"     ; required
  res@mpMinLatF             = -90.         ; min lat
  res@mpMaxLatF             = -60.         ; max lat
  res@mpMinLonF             = -50.         ; min lon
  res@mpMaxLonF             =  60.         ; max lon

  ;res@mpGeophysicalLineColor       = "white"
  ;res@mpGeophysicalLineThicknessF = 1.5

  if fcsttime .ge. 0 .and. fcsttime .le. 9 then
  res@tiMainString          = "MAITRI:INI TM:"+initime+" FCST TM:0"+fcsttime+"Hrs" 
  else
  res@tiMainString          = "MAITRI:INI TM:"+initime+" FCST TM:"+fcsttime+"Hrs" 
  end if

  res@pmTickMarkDisplayMode = "Always"     ; turn on automatic tickmarks
  res@mpGridAndLimbOn       = True         ; turn on lat/lon lines
  ;res@mpGridMaskMode        = "MaskLand"   ; Mask grid over land.

;;;;;;;;;;;;;;;; Attach a polymark ;;;;;;;;;;;;;;

  gsres                    = True
  ;gsres@gsMarkerColor      = "white"
  gsres@gsMarkerIndex      = 16
  ;gsres@gsMarkerSizeF      = 0.014
  ;gsres@gsMarkerThicknessF = 2.

;;;;;;;;;;;;;;;; Attach some text strings ;;;;;;;;;;

  txres               = True
  ;txres@txFontColor   = "white"
  txres@txFontHeightF = 0.015
  ;txres@txAngleF      = -45.

;;;;;;;;;;;;;;; Drawing of plot ;;;;;;;;;;;;;;;;;;;;

  plot = gsn_csm_contour_map_ce(wks,tc2,res)       
  ;plot = gsn_csm_contour_map_polar(wks,tc2,res)       

  mark = gsn_add_polymarker(wks,plot,slon,slat,gsres)

  text = new(1,graphic)
  text = gsn_add_text(wks,plot,name,slon,slat+2,txres)

 draw(plot)
 frame(wks)
 delete(wks)

 end do

end