;=======================================================================================
; This NCL 
;
; Written by Colin Zarzycki (zarzycki@ucar.edu)
; Version 0.1 (6/5/2016) - DCMIP-2016 release
;=======================================================================================

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
begin

;=======================================================================================


;whichVar="PS"          ; PS or WIND
top_to_bottom=True     ; is data top to bottom? (false for bottom to top)
model="DYNAMICO"
out_type="png"
hemisphere = (/"SOUTH","NORTH"/)

physdir = (/ "variant/dry/", "variant/moist_kessler/", "variant/moist_full_physiq/", "variant/moist_full_physiq_0.5deg/",  "variant/moist_full_physiq_30min/" /)
; needs to be a catted file ordered time, lev, lat, lon (or time, lat, lon for PS)
filename="output_dcmip2016_regular.nc"
;=======================================================================================

maxLat=90.0
eqLat=0.0
minLat=-90.0

;=======================================================================================
;; SET UP THE PANEL PLOT WITH LEFT AND RIGHT COLUMNS
  wks_type = "png"
; Open the workspace on the current file
  wks = gsn_open_wks(wks_type,"BAROCLINIC_WAVE_BC_physCompare.png")

  lgres = True
  lgres@lgLineColors     = (/"Red","Orange","Blue","Green","Black"/)
  lgres@lgLineThicknessF   = 5.0            ; legend line thickness
  lgres@lgItemType         = "Lines"        ; show lines only (default)
  lgres@lgLineDashSegLenF  = 0.5
  lgres@lgMonoDashIndex = True
  lgres@lgLabelFontHeightF = .04            ; set the legend label font thickness
  lgres@vpWidthF           = 0.25            ; width of legend (NDC)
  lgres@vpHeightF          = 0.15            ; height of legend (NDC)
  lgres@lgPerimColor       = "Black"        ; draw the box perimeter in orange
  lgres@lgPerimThicknessF  = 2.0            ; thicken the box perimeter
  lgres@lgBoxBackground = "White"
  lgres@lgPerimFill = 0
  lgres@lgPerimFillColor = "White"
  lgres@lgBoxMinorExtentF = 0.25

plotBCSF = new(2,graphic)
do hemi = 0,1
	; Loop over dry, Kessler, and full physics cases
	do physCase = 0,4
		infile = (/ physdir(physCase), filename /)
		rfile = str_concat(infile)
		a = addfile(rfile,"r")

	; Get the minimum pressure
		whichVar = "PS"
		if (hemi .eq. 0) then
			PS=a->PS(:,{minLat:eqLat},:)
		end if
		if (hemi .eq. 1) then
			PS=a->PS(:,{eqLat:maxLat},:)
		end if
		PS=PS/100.                                 ; convert from Pa to hPa/mb
		PS@long_name="Surface pressure"
		PS@units="hPa"
		var_vtime = dim_min_n_Wrap(PS,(/1,2/))     ; for now, just find min over lat/lon domain at each time

		time = var_vtime&time_counter / 86400

	  ; CREATE THE PLOTS AND SEND THEM OUT TO .PDF FILES
		if (physCase .eq. 0) then
			;; SET THE RESOURCES ATTRIBUTES FOR THE XY PLOTS
			rescom = True
			; Graphics resources
			rescom@gsnDraw = False
			rescom@gsnFrame = False
			rescom@gsnMaximize = True
			rescom@gsnPaperOrientation = "portrait"
			; Plot area resources
			rescom@tmXMajorGrid = True
			rescom@tmYMajorGrid = True
			rescom@tmXMajorGridLineDashPattern = 2
			rescom@tmYMajorGridLineDashPattern = 2
			rescom@tiMainString = "DCMIP1 "+model+" "+whichVar+" v. time "+hemisphere(hemi) ; add title
			rescom@tiXAxisString = "Time (days)"  
			rescom@tiYAxisString = var_vtime@long_name+" ("+var_vtime@units+")"
			rescom@xyLineColor = "Red"
			rescom@xyDashPattern = 0
			rescom@xyLineThicknessF = 5.0
			plotBCSF(hemi) = gsn_csm_xy(wks,time,var_vtime,rescom)
		end if

		if (physCase .gt. 0) then
			dataid = create "xyData" coordArraysClass defaultapp
				"caYArray": var_vtime
				"caXArray": time
			end create

			dataspec = NhlAddData(plotBCSF(hemi),"xyCoordData",dataid)

			if (physCase .eq. 1) then
				lcolor = "Orange"
			end if 
			if (physCase .eq. 2) then
				lcolor = "Blue"
			end if
			if (physCase .eq. 3) then
				lcolor = "Green"
			end if
			if (physCase .eq. 4) then
				lcolor = "Black"
			end if
			setvalues dataspec
				"xyLineColor"      : lcolor
				"xyDashPattern"    : 0
				"xyLineThicknessF" : 5.0
			end setvalues
		end if
		delete(PS)
	end do

	; Add a legend and draw the plot
	lbid = gsn_create_legend(wks,5,(/ "Dry","Kessler","Full Physics", "Full Physics 1/2 deg", "Full physics 30min" /),lgres)
	amres = True
	amres@amParallelPosF   = -0.3
	amres@amOrthogonalPosF = 0.3
	annoid1 = gsn_add_annotation(plotBCSF(hemi),lbid,amres)   ; add legend to plot

	; Delete some stuff for the next plot in the panel
	delete(lbid)
	delete(amres)
	delete(dataid)
	delete(dataspec)
	delete(annoid1)
	delete(rescom)
end do
  ; SET PANEL RESOURCES
  resP                        = True                ; modify the panel plot
  resP@gsnMaximize            = True                ; maximize plots
  resP@gsnPaperOrientation    = "landscape"
  ;resP@gsnFrame               = False               ; don't advance panel plot
  ;resP@amJust                 = "TopLeft"           ; put in top left of plot
  ;resP@gsnPaperMargin         = 0.5
  ;resP@gsnAttachBorderOn = False
  ;resP@gsnPanelLabelBar = True
  ;resP@lbLabelFontHeightF = 0.02
  ;resP@lbLabelStride = 2
  ;resP@lbTitleOn = True 
  ;resP@lbTitleString = "Potential Temperature (K)"
  ;resP@lbTitleFontHeightF= 0.02
  ;resP@pmLabelBarOrthogonalPosF = -0.03
  ;resP@pmLabelBarWidthF      = 0.8
  ;resP@pmLabelBarHeightF     = 0.075
  
  resP@gsnPanelLeft = 0.05
  resP@gsnPanelRight = 0.95
  gsn_panel(wks,plotBCSF,(/1,2/),resP)
end