load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
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/shea_util.ncl"
;************************************************
;----------------------------------------------------------------------
; This procedure attaches lat/lon labels to a masked lambert plot
;
; You will likely need to change lat_values and/or lon_values to
; contain the locations where you want lat/lon labels.
;----------------------------------------------------------------------
procedure add_lc_labels(wks,map,minlat,maxlat,minlon,maxlon)
local lat_values, nlat, lat1_ndc, lat2_ndc, lon1_ndc, lon2_ndc,slope,txres, \
lon_values, PI, RAD_TO_DEG, dum_lft, dum_rgt, dum_bot
begin

  PI         = 3.14159
  RAD_TO_DEG = 180./PI

;---Pick some "nice" values for the latitude labels.
  lat_values = ispan(toint(minlat),toint(maxlat),10) * 1.
  nlat       = dimsizes(lat_values)
;
; We need to get the slope of the left and right min/max longitude lines.
; Use NDC coordinates to do this.
;
  lat1_ndc = new(1,float)
  lon1_ndc = new(1,float)
  lat2_ndc = new(1,float)
  lon2_ndc = new(1,float)
  datatondc(map,minlon,lat_values(0),lon1_ndc,lat1_ndc)
  datatondc(map,minlon,lat_values(nlat-1),lon2_ndc,lat2_ndc)
  slope_lft = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc)

  datatondc(map,maxlon,lat_values(0),lon1_ndc,lat1_ndc)
  datatondc(map,maxlon,lat_values(nlat-1),lon2_ndc,lat2_ndc)
  slope_rgt = (lat2_ndc-lat1_ndc)/(lon2_ndc-lon1_ndc)

;---Set some text resources
  txres               = True
  txres@txFontHeightF = 0.03
  txres@txPosXF       = 0.8
;
; Loop through lat values, and attach labels to the left and
; right edges of the masked LC plot.  The labels will be 
; rotated to fit the line better.
;
  dum_lft = new(nlat,graphic)    ; Dummy array to hold attached strings.
  dum_rgt = new(nlat,graphic)    ; Dummy array to hold attached strings.
  do n=0,nlat-1
; Add extra white space to labels.
    lat_label_rgt = "                  "  +lat_values(n) 

;---Check if North, South, or Zero
    if(lat_values(n).lt.0) then
     lat_label_lft = lat_values(n) + "S                "
      lat_label_rgt = lat_label_rgt + "S"
    end if
    if(lat_values(n).gt.0) then
      lat_label_lft = lat_values(n) + "N                "
      lat_label_rgt = lat_label_rgt + "N"
    end if
    if(lat_values(n).eq.0) then
      lat_label_lft = lat_values(n) + "~S~o~N~              "
    end if

;---Left label
    txres@txAngleF = RAD_TO_DEG * atan(slope_lft) - 90
    dum_lft(n) = gsn_add_text(wks,map,lat_label_lft,minlon,lat_values(n),txres)

;---Right label
    txres@txAngleF = RAD_TO_DEG * atan(slope_rgt) + 90
    dum_rgt(n) = gsn_add_text(wks,map,lat_label_rgt,maxlon,lat_values(n),txres)
  end do

;----------------------------------------------------------------------
; Now do longitude labels. These are harder because we're not
; adding them to a straight line.
;
; Loop through lon values, and attach labels to the bottom edge of the
; masked LC plot.
;
  delete(txres@txPosXF)
  txres@txPosYF = -5.0

;---Pick some "nice" values for the longitude labels.
  lon_values = ispan(toint(minlon+10),toint(maxlon-10),10) * 1.
  nlon       = dimsizes(lon_values)

  dum_bot = new(nlon,graphic)    ; Dummy array to hold attached strings.

  do n=0,nlon-1
;
; For each longitude label, we need to figure out how much to rotate
; it, so get the approximate slope at that point.
;
    datatondc(map,lon_values(n)-0.25,minlat,lon1_ndc,lat1_ndc)
    datatondc(map,lon_values(n)+0.25,minlat,lon2_ndc,lat2_ndc)
    slope_bot = (lat1_ndc-lat2_ndc)/(lon1_ndc-lon2_ndc)
    txres@txAngleF  = atan(slope_bot) * RAD_TO_DEG
;
; Create longitude label. Add extra carriage returns to
; move label away from plot.
;
;---Check if East, West, or Zero
    lon_label_bot = "                 " + abs(lon_values(n))
    if(lon_values(n).lt.0) then
      lon_label_bot = lon_label_bot + "W"
    end if
    if(lon_values(n).gt.0) then
      lon_label_bot = lon_label_bot + "E"
    end if

;---Attach to map.
    dum_bot(n) = gsn_add_text(wks,map,lon_label_bot,lon_values(n),minlat,txres)
  end do
end


;----------------------------------------------------------------------
; Main code.
;----------------------------------------------------------------------
begin
;************************************************
; read in data1 file for binary
;************************************************
  nlat = 61
  mlon = 61
  NMOD = 8      ; principal modes
  model = "mfpspa"
  MODEL = "MFPSPA"
  season = "MAM"
  season1 = "mam"
  pModes = "08"
  pModes1 = 08

  nf   = new((/61,61/),float)
  nf01 = new((/5,61,61/),float)
  nf4d = new((/NMOD,5,61,61/),float)
  mf   = new((/61,61/),float)
  mf01 = new((/5,61,61/),float)
  mf4d = new((/NMOD,5,61,61/),float)
  kf   = new((/61,61/),float)
  kf01 = new((/5,61,61/),float)
  kf4d = new((/NMOD,5,61,61/),float)
  var1 = "mslp"
  var2 = "z500"
  var3 = "evs"

  fname = "output/"+model+"_"+ season + "_" + pModes
  wks = gsn_open_wks("pdf",fname)                   ; open a pdf file

  path1 = model+"_"+ var1 +"_"+ season1 +".score.5x61x61"
  path2 = model+"_"+ var2 +"_"+ season1 +".score.5x61x61"
  path3 = model+"_"+ var3 +"_"+ season1 +".score.5x61x61"
  N=0
  do nmod = 0,NMOD-1,1
  do ntim = 0,4,1
    nf = fbindirread(path1, N, (/nlat,mlon/), "float")
    nf01(ntim,:,:) = nf(:,:)
    nf4d(nmod,ntim,:,:)=nf01(ntim,:,:)
    mf = fbindirread(path2, N, (/nlat,mlon/), "float")
    mf01(ntim,:,:) = mf(:,:)
    mf4d(nmod,ntim,:,:)=mf01(ntim,:,:)
    kf = fbindirread(path3, N, (/nlat,mlon/), "float")
    kf01(ntim,:,:) = kf(:,:)
    kf4d(nmod,ntim,:,:)=kf01(ntim,:,:)
    N=N+1
  end do
  end do
;        print(nf01(1,:,:))                                 ; show in the screen

;************************************************ 
; Create lon and lat coordinates variables 
;************************************************
  minlat = 10.
  maxlat = 40.
  minlon = 100.
  maxlon = 140.

  lat = new(61, float) 
  lat = fspan( 10, 40, 61 ) 
  lat!0 = "lat" 
  lat@long_name = "latitude" 
  lat@units = "degrees-north" 
  lat&lat = lat 

  lon = new(61,float) 
  lon = fspan( 100, 140, 61 ) 
  lon!0 = "lon" 
  lon@long_name = "longitude" 
  lon@units = "degrees-east" 
  lon&lon = lon 

  nlat@units = "degrees_north"
  mlon@units = "degrees_east"

  nf!0 = "lat"                ; name dimensions
  nf!1 = "lon"
  nf&lat = fspan( 10, 40, 61 )
  nf&lon = fspan( 100, 140, 61 )
  mf!0 = "lat"                ; name dimensions
  mf!1 = "lon"
  mf&lat = fspan( 10, 40, 61 )
  mf&lon = fspan( 100, 140, 61 )
  kf!0 = "lat"                ; name dimensions
  kf!1 = "lon"
  kf&lat = fspan( 10, 40, 61 )
  kf&lon = fspan( 100, 140, 61 )

;  print(nf&lon)

;************************************************
; create plot
;************************************************
;  gsn_define_colormap(wks,"BlAqGrYeOrRe")        ; choose color map
;  gsn_define_colormap(wks,"gsdtol") 
;************************************************
; resource list for first data array
;************************************************
  res          =    True                          ; plot mods desired
  res@mpProjection          = "LambertConformal"  ; choose projection
  res@mpFillOn              = True               ; turn off gray map
  res@mpMinLatF             = minlat
  res@mpMaxLatF             = maxlat
  res@mpMinLonF             = minlon
  res@mpMaxLonF             = maxlon
;  res@mpGeophysicalLineColor         = "black"       ; color of cont. outlines
;  res@mpGeophysicalLineThicknessF    = 1             ; thickness of outlines

  res@cnInfoLabelOn         = False               ; turn off info label
  res@cnFillOn              = False                ; turn on color
;  res@cnFillMode            = "RasterFill"
  res@cnLinesOn             = True               ; turn off contour lines 
  res@cnLineLabelsOn        = False
  res@cnLineThicknessF      = 1.7                  ; thickness of contour
;  res@cnMissingValPerimOn   = True                ; turn on the missing value perimeter

;  res@cnLevelSelectionMode  = "ManualLevels"      ; set manual contour levels
;  res@cnMinLevelValF        =  -0.2                ; set min contour level
;  res@cnMaxLevelValF        =  0.2                ; set max contour level
;  res@cnLevelSpacingF       =  0.05                ; set contour spacing

  res@gsnContourNegLineDashPattern = 1       ; sets negative contours to dash pattern 1

  tiMain = (/"Day-2","Day-1","Day+0","Day+1","Day+2"/)
  res@tiMainOn              = True
  res@tiMainFont            = "helvetica-bold"   ; set Font form
  res@tiMainFontHeightF     = 0.04               ; title font size

  res@tiYAxisFontHeightF    = 0.07            ; Increase size of fonts
  res@tiYAxisAngleF         = 60
  res@tiYAxisFont           = "helvetica-bold"

;  res@gsnSpreadColors       = True                ; use full colormap 
;  res@gsnSpreadColorStart   = 0                   ; start at color 
;  res@gsnSpreadColorEnd     = 95                  ; end at color 

  res@gsnDraw               = False               ; don't draw
  res@gsnFrame              = False               ; don't advance frame
;  res@gsnMaximize           = True                ; maximize image
;  res@gsnSpreadColors       = True                ; use full range of colors

  res@sfXCStartV            = minlon
  res@sfXCEndV              = maxlon
  res@sfYCStartV            = minlat
  res@sfYCEndV              = maxlat
  
  res@gsnMaskLambertConformal = True            ; turn on lc masking
  res@mpGridAndLimbOn         = True  
  res@mpGridLatSpacingF       = 10
  res@mpGridLonSpacingF       = 5
  res@gsnAddCyclic            = False

  plot = new(15,graphic)                          ; create graphics array

  ntim=0
  do m = 0,4
     res@tiMainString          =  tiMain(m)
     res@cnLevelSelectionMode  = "ManualLevels"      ; set manual contour levels
     res@cnMinLevelValF        =  -3.0                ; set min contour level
     res@cnMaxLevelValF        =  3.0                ; set max contour level
     res@cnLevelSpacingF       =  0.6                ; set contour spacing
;  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
;  res@cnLevels            = (/-0.4,-0.2,-0.12,-0.08,-0.04,-0.0,0.04,0.08,0.12,0.2,0.25,0.3,0.35,0.4/)

     if m.eq.0
        res@tiYAxisString         = "EVS"
     else
        res@tiYAxisString         = ""
     end if
     plot(m) =gsn_csm_contour_map(wks,kf4d(pModes1-1,ntim,:,:) ,res)        ; create plot
;     plot(m) = ShadeGtContour(plot(m),0,17) ;shade contours >0
;---Attach latitude labels
;     add_lc_labels(wks,plot(m),minlat,maxlat,minlon,maxlon)
     ntim=ntim+1
  end do
  ntim=0
  res@tiMainOn              = False
  do m = 5,9
     res@tiMainString          =  tiMain(m-5)
     res@cnLevelSelectionMode  = "ManualLevels"      ; set manual contour levels
     res@cnMinLevelValF        =  -3.6                ; set min contour level
     res@cnMaxLevelValF        =  3.6                ; set max contour level
     res@cnLevelSpacingF       =  0.4                ; set contour spacing
;  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
;  res@cnLevels            = (/-0.4,-0.2,-0.12,-0.08,-0.04,-0.0,0.04,0.08,0.12,0.2,0.28,0.3,0.38,0.4/)
     if m.eq.5
        res@tiYAxisString         = "Z500"
     else
        res@tiYAxisString         = ""
     end if
     plot(m) =gsn_csm_contour_map(wks,mf4d(pModes1-1,ntim,:,:) ,res)        ; create plot
;     plot(m) = ShadeGtContour(plot(m),0,17) ;shade contours >0
;---Attach latitude labels
;     add_lc_labels(wks,plot(m),minlat,maxlat,minlon,maxlon)
     ntim=ntim+1
  end do
  ntim=0
  do m = 10,14
     res@tiMainString          =  tiMain(m-10)
     res@cnLevelSelectionMode  = "ManualLevels"      ; set manual contour levels
     res@cnMinLevelValF        =  -3.6                ; set min contour level
     res@cnMaxLevelValF        =  3.6                ; set max contour level
     res@cnLevelSpacingF       =  0.4                ; set contour spacing
     if m.eq.10
        res@tiYAxisString         = "MSLP"
     else
        res@tiYAxisString         = ""
     end if
     plot(m) =gsn_csm_contour_map(wks,nf4d(pModes1-1,ntim,:,:) ,res)        ; create plot
;     plot(m) = ShadeGtContour(plot(m),0,17) ;shade contours >0
;---Attach latitude labels
;     add_lc_labels(wks,plot(m),minlat,maxlat,minlon,maxlon)
     ntim=ntim+1
  end do

;*******************************************
; draw panel plot with title
;*******************************************
  pres                   = True                       ; mod panel plot
  pres@gsnFrame          = False
;  pres@gsnPanelRight     = 0.5
  pres@gsnMaximize    = True                ; maximize plots
 
  pres@txString          = MODEL+"_" + season + "_" + pModes           ; add common title
  pres@txFont            = "helvetica-bold" 
  pres@txFontHeightF     = 0.025                       ; font size
  pres@txPosXF           = 0.5
;  pres@pmLabelBarParallelPosF = 0.015
;  pres@gsnPanelXF        = 10.0

  gsn_panel(wks,plot(0:14),(/3,5/),pres)
  frame(wks)

  end