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
;; Temperature in celsuis  ;;
 temp = (/-50.,-40.,-30., -20.,-10.,  0., 10., 20., 30., 40., 50. /)
;==============================================
;==========Hyposometric Equation ==============
;==============================================
R=287.14
T0=15+273.
P0=1012 ;mb
g= 9.8 ;m/s

pres1 = (/1000.,900.,800.,700.,600.,500.,400.,300.,200.,100.,50/) 
temp1 = (/30   ,25  ,20  ,10  ,0   ,-10 ,-20 ,-30 ,-40 ,-50 ,-60/)

H =new((/dimsizes(pres1)/),"float")
z1=new((/dimsizes(pres1)/),"float")

do i=0,dimsizes(pres1)-1
  H(i)=R*(temp1(i)+273.15)/g   
  z1(i)=H(i)*log(P0/pres1(i))
  print("height is "+z1(i)+" at pressure "+pres1(i))
end do
;==============================================
wks                        = gsn_open_wks("png","with temp")        ; Open an ps workstation.
plot_res                   = True
plot_res@gsnMaximize       = True                           ; must include w/ Paper Orientation
plot_res@vpWidthF          = 0.9	                    ; set width of plot 
plot_res@vpHeightF         = 0.8  	                    ; set height of plot
plot_res@vpXF              = 0.1	                    ; set left hand side start point of plot
plot_res@vpYF              = 0.9	                    ; set top start point of plot
plot_res@gsnDraw           = False
plot_res@gsnFrame          = False
plot_res@tiXAxisString     = "Temperature"                  ; x-axis label
plot_res@tiYAxisString     = "log Pressure(Height) "        ; y-axis label
plot_res@tiXAxisFont        = "triplex_italic"
plot_res@tiXAxisFontHeightF = 0.02
plot_res@tiYAxisFont        = "triplex_italic"
plot_res@tiYAxisFontHeightF = 0.02
plot_res@tiMainFont        = "triplex_italic"
plot_res@tiMainFontHeightF = 0.02
plot_res@cnLineLabelsOn    = True                           ; no contour line labels
plot_res@cnLevelSpacingF   = 10
plot_res@cnInfoLabelFont   = "triplex_italic"
plot_res@cnInfoLabelOn     = False
plot_res@trXMinF           = min(temp)
plot_res@trXMaxF           = max(temp)
plot_res@trYMinF           = min(z1)
plot_res@trYMaxF           = max(z1)
plot_res@trGridType        = "TriangularMesh"
plot_res@trYAxisType       = "LinearAxis"  
;=============================================
;============ the skewed temperature =========
;=============================================
T=new((/dimsizes(z1),dimsizes(temp)/),"float")
T!0 = "z1"
T!1 = "temp"
T&z1 = z1
T&temp = temp 
;plot_res@sfXArray              = temp            ; X axes data points
;plot_res@sfYArray              = z1         ; Y axes data points
do i=0,dimsizes(z1)-1
 do j=0,dimsizes(temp)-1
   T(i,j) = temp(j)-10*i
 end do
end do

plot_res@cnLineColor     = "green"
plot_res@cnLineLabelFontColor = "green"
plot = gsn_contour(wks,T,plot_res)     
;==================================================
;=========== Drawing the pressure isobars =========
;=========== The Hyposometric equation. ===========
;==================================================
plot_text                 = True
plot_text@txFont          = 21
plot_text@txFontAspectF   = 1.9
plot_text@txFontHeightF   = 0.02
plot_text@txFontColor     = "red"
plot_line                 = True
plot_line@gsLineColor     = "red"

plot1=new(dimsizes(z1),"graphic")
plot2=new(dimsizes(z1),"graphic")
x1=(/temp(0),temp(dimsizes(temp)-1)/)
do i=0,dimsizes(z1)-1
  y1=(/z1(i),z1(i)/)
  print("y1 is "+y1)
  plot1(i)=gsn_add_polyline(wks,plot,x1,y1,plot_line)
  ;note that we increase the z1(i) by 250 as height is up to 20.000 , we use the 
  ; height  of axes itself .
  plot2(i)=gsn_add_text(wks,plot,pres1(i)+"mb",temp(0)+10,z1(i)+250,plot_text)
end do
draw(plot)   
frame(wks)
end