; Calculate the climalogical mean tropopause height
; using the monthly temperatue 
; from NCEP/NCAR reanalysis1


begin

dir  = "/home/wangn/data/NCEP1/"
tf   = addfile(dir + "air.mon.mean.nc", "r")


temp = tf->air  
temp = temp + 273.15                               ; units conversion (from degC to K) 
ta   = dim_avg_n_Wrap(temp, 0)                     ; calculate the climatological mean                  
pre  = ta&level                                    ; get the pressure levels (hPa)

trop = trop_wmo(pre(::-1), ta(lat|:,lon|:,level|::-1), 0, False)
copy_VarCoords(dim_avg_n_Wrap(ta,0), trop)
trop@long_name  = "Tropopause pressure"
trop@_FillValue = ta@_FillValue
trop@units = "hPa"


wks    = gsn_open_wks("x11","ncep_clim_trop")

cmap   = read_colormap_file("GMT_no_green")

res  = True

res@vpWidthF          = 0.7
res@vpHeightF         = 0.4 

res@gsnFrame          = False              ; don't draw yet
res@gsnDraw           = False              ; don't advance frame

res@cnLinesOn         = False
res@cnFillOn          = True
res@cnFillPalette     = cmap 
res@cnMissingValFillColor = "white"

res@cnLevelSelectionMode = "ManualLevels"       ; manually select levels
res@cnMinLevelValF       = 80                   ; min level
res@cnMaxLevelValF       = 300
res@cnLevelSpacingF      = 20
;------------------------------------- 

res2  = True

res2@gsnFrame  = True
res2@gsnDraw   = True 
res2@gsnZonalMeanXMinF    = 0      ; min X-axis value for zonal mean plot  
res2@gsnZonalMeanXMaxF    = 300      ; max X-axis value for zonal mean plot  


plot1 = gsn_csm_contour_map(wks, trop, res)
plot2 = gsn_csm_attach_zonal_means(wks, plot1, trop, res2)

end