source: configs/dcmip2016/PLOTS/dcmip1-SFP0-trace.ncl @ 439

;=======================================================================================
; 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