;; -*- mode: emacs-lisp; coding: emacs-mule; -*-;  -*- mode: 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/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"


ana=addfile("sfc_toa_fine.nc","r")
Absp= ana->Absp_fine_no_ds1
sfc= ana->sfc_fine_no_ds1
delete(ana)

ana=addfile("vertical_profile_fine_ocean.nc","r")
lev3=ana->lev
Absp3D_ocean=ana->Absp3D_fine_no_ds_tot
lat=ana->lat
lon=ana->lon
delete(ana)

ana=addfile("vertical_profile_fine_land.nc","r")
Absp3D_land=ana->Absp3D_fine_no_ds_tot
delete(ana)


ana=addfile("vertical_profile_fine_tot.nc","r")
Absp3D_tot=ana->Absp3D_fine_no_ds_tot
delete(ana)



jlat  = dimsizes( lat )
rad    = 4.0*atan(1.0)/180.0
re     = 6371220.0
rr     = re*rad
dlon   = abs(lon(2)-lon(1))*rr
dx     = dlon*cos(lat*rad)                                     ;lat can have variable spacing
dy     = new ( jlat, typeof(dx))                                                          ; close enough
dy(0)  = abs(lat(2)-lat(1))*rr
dy(1:jlat-2)  = abs(lat(2:jlat-1)-lat(0:jlat-3))*rr*0.5   
dy(jlat-1)    = abs(lat(jlat-1)-lat(jlat-2))*rr
area   = dx*dy                                ; cell area function of latitude only
clat   = cos(lat*rad)


tot= wgt_areaave(Absp3D_tot,clat,1.0,0)*24*3600
ocean= wgt_areaave(Absp3D_ocean,clat,1.0,0)*24*3600
land= wgt_areaave(Absp3D_land,clat,1.0,0)*24*3600

data2      = new((/3,25/),double)

 data2(0,:) = tot
 data2(1,:) = ocean
 data2(2,:) =land



wks = gsn_open_wks("eps","atm_sfc_fine")
;gsn_define_colormap(wks,"BlAqGrYeOrRe")
;gsn_define_colormap(wks,"BlWhRe")
 plot = new(3,graphic)   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  res                     = True               ; plot mods desired
  res@cnFillOn            = True              ; turn on color   
  res@gsnSpreadColors     = True               ; use full color map
  res@cnLinesOn           = False              ; no contour lines
  res@cnLineLabelsOn      = False              ; no line labels


 ; res@gsnAddCyclic        = False 
  res@gsnDraw             = False           ; don't draw
  res@gsnFrame            = False           ; don't advance frame
  res@cnInfoLabelOn       = False           ; turn off cn info label
  res@cnFillOn            = True            ; turn on color
  res@gsnSpreadColors     = True            ; spread out color table
  res@gsnSpreadColorStart = 2               ; color to startt

  res@lbLabelBarOn        = True          ; turn off individual cb's

;  res@cnLevelSelectionMode =  "AutomaticLevels"   ; AutomaticLevels




res@cnLevelSelectionMode =  "ManualLevels"   ; AutomaticLevels
res@cnMinLevelValF       = -20.              ; min level
res@cnMaxLevelValF       =  20.               ; max level
res@cnLevelSpacingF      =  2               ; interval


   
res@lbTitleString    = "W/m^2"                ; title string
res@lbTitlePosition  = "Right"              ; title position
res@lbTitleFontHeightF= .025                ; make title smaller
res@lbTitleDirection = "Across"             ; title direction


res@gsnStringFontHeightF = 0.025

res@gsnLeftString       = "Solar radiation absorption in the atmosphere"         ; change left string
plot(0) = gsn_csm_contour_map_ce(wks,Absp,res)
res@gsnLeftString       = "Surface net flux"         ; change left string
plot(1) = gsn_csm_contour_map_ce(wks,sfc,res)

;*********************************************************************************
res2                   = True                       ; plot mods desired
res2@lbLabelFontHeightF =.022
res2@gsnStringFontHeightF = 0.025

res2@gsnDraw = False
res2@gsnFrame = False
;res2@tiMainString      = "Vertical profile of atmospheric heating rate"             ; add title
res2@tiYAxisString = "Pressure (hPa)"  ; xaxis string
res2@tiXAxisString = "Atmospheric heating rate (K/day)";Pressure levels (hPa)"

res2@gsnXRefLine       = 0.    
res2@gsnXRefLineDashPatterns = 4

res2@trYMinF     = 100		; set minimum Y-axis value
res2@trYMaxF     =  1000		; set maximum Y-axis value
res2@trXMinF     = -0.02      ; set minimum X-axis value
;res2@trXMaxF     = 0.09      ; set maximum X-axis value
 res2@trYReverse = True                     ; reverse Y-axis
res2@xyDashPattern  = 0                         ; Make curves all solid



res2@xyLineThicknesses = (/2.0, 2.0, 2.0/)
res2@xyLineColors   = (/"blue","red","green"/) ; 3 different colors
; add a legend
 res2@pmLegendDisplayMode    = "Always"              ; turn on legend
 
 res2@pmLegendSide           = "Top"                 ; Change location of 
 res2@pmLegendParallelPosF   = .90                   ; move units right
 res2@pmLegendOrthogonalPosF = -0.6                  ; more neg = down
 
 res2@pmLegendWidthF         = 0.06                  ; Change width and
 res2@pmLegendHeightF        = 0.25                  ; height of legend.
 res2@lgLabelFontHeightF     = .015                   ; change font height
 res2@lgPerimOn              = False                 ; no box around
; labels for the legend
res2@xyExplicitLegendLabels = (/"Total","Oceans","Land"/)
;plot  = gsn_csm_xy (wks,Absp3D_fine_no_ds5,lev3,res) ; create plot
plot(2)  = gsn_csm_xy (wks,data2,lev3,res) ; create plot
;************************************************
; create panel
;************************************************
resP                  = True                    ; modify the panel plot
resP@gsnPanelLabelBar = False                  ; add common colorbar
resP@gsnPanelFigureStrings = (/"a)","b)","c)"/)      ; add strings to panel
gsn_panel(wks,plot,(/2,2/),resP)