; use NCL to calculate and plot atmosheric meridional overturning stream function
; Eli: modified from Wanying's
; to run from terminal and produce a pdf:
; $ module load ncl
; $ ncl Hadley_stream_function.ncl
;
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"		; functions required to
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"           	; plot.  include before
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"         	; begin
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
begin

; plot control
showbottom=1000;
showtop=100;

;H=1
a = 6.37122e06 ;  radius of the earth
PI=3.1415926
g=9.8
bottom=""

; read in data
filename_V="/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/atm/proc/tseries/monthly/V/"+"b.e11.B1850C5CN.f09_g16.005.cam.h0.V.040001-049912.nc"
filename_PS="/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/atm/proc/tseries/monthly/PS/"+"b.e11.B1850C5CN.f09_g16.005.cam.h0.PS.040001-049912.nc"
f=addfile(filename_V,"r");
f_PS=addfile(filename_PS,"r");
lat=f->lat
pres_=f->lev
indpres=ind(pres_.gt.showtop .and. pres_.lt.showbottom)
pres=pres_(indpres)
phi = lat*PI/180.0     ; latitude in radians
nlat=dimsizes(lat)
V=f->V(0,indpres,:,:)
PS=f_PS->PS(0,:,:)

;;; Calculate meridional circulation
psi=zonal_mpsi(V,lat,pres*100,PS) ; multiply 100 because pres is in mb
;unit kg/s

psi!0="lev"
psi&lev=pres
psi&lev@units="mb" ; added by Eli
psi!1="lat"
psi&lat=lat
psi@units="kg/s"
psi@long_name="Meridional Streamfunction"
psi=(/psi/1e09/)
psi@units="svp"

;***********************************
;plot
;***********************************

  wks   = gsn_open_wks ("pdf", "Output/streamfunction_NCL.pdf")        ; open workstation
  gsn_define_colormap(wks,"NCV_blu_red")
  res                      = True
  res@tiMainString      = "Meridional Circulation:"; title
  res@cnFillOn             = True            ; turn on color fill
  res@cnLinesOn            = False            ; turn lines on/off         ; True is default
  res@cnLineLabelsOn       = False            ; turn line labels on/off   ; True is default
  res@gsnMajorLatSpacing = 30
  res@cnLevelSelectionMode = "ManualLevels"       ; manually select levels
  res@trYMaxF               = showbottom
  res@trYMinF               = showtop
  res@gsnFrame=True
  res@tmYLMode = "Explicit"              ; explicitly set YR tick marks
  res@tmYLValues = (/850.,500.,300.,150.,100.,70.,50.,30.,20.,10.,5./)
  res@tmYLLabels = (/850.,500.,300.,150.,100.,70.,50.,30.,20.,10.,5./)
  res@lbLabelFontHeightF=0.03
  res@tmXBLabelFontHeightF=0.03
  res@tmYLLabelFontHeightF=0.03
  res@tmYRLabelFontHeightF=0.03
  res@tiYAxisFontHeightF=0.03
  res@pmLabelBarHeightF    = 0.1
  dlim       = -150.                 ; min level
  ulim       = 150.                  ; max leve
  space      = 10.                  ; contour spacing
  res@cnMinLevelValF       = dlim                 ; min level
  res@cnMaxLevelValF       = ulim                 ; max leve
  res@cnLevelSpacingF      = space                 ; contour spacing

  plot  = gsn_csm_pres_hgt(wks,psi,res)

end