1 | ;======================================================================================= |
---|
2 | ; This NCL |
---|
3 | ; |
---|
4 | ; Written by Colin Zarzycki (zarzycki@ucar.edu) |
---|
5 | ; Version 0.1 (6/5/2016) - DCMIP-2016 release |
---|
6 | ;======================================================================================= |
---|
7 | |
---|
8 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
---|
9 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
---|
10 | begin |
---|
11 | |
---|
12 | ;======================================================================================= |
---|
13 | |
---|
14 | |
---|
15 | ;whichVar="PS" ; PS or WIND |
---|
16 | top_to_bottom=True ; is data top to bottom? (false for bottom to top) |
---|
17 | model="DYNAMICO" |
---|
18 | out_type="png" |
---|
19 | hemisphere = (/"SOUTH","NORTH"/) |
---|
20 | |
---|
21 | physdir = (/ "variant/dry/", "variant/moist_kessler/", "variant/moist_full_physiq/", "variant/moist_full_physiq_0.5deg/", "variant/moist_full_physiq_30min/" /) |
---|
22 | ; needs to be a catted file ordered time, lev, lat, lon (or time, lat, lon for PS) |
---|
23 | filename="output_dcmip2016_regular.nc" |
---|
24 | ;======================================================================================= |
---|
25 | |
---|
26 | maxLat=90.0 |
---|
27 | eqLat=0.0 |
---|
28 | minLat=-90.0 |
---|
29 | |
---|
30 | ;======================================================================================= |
---|
31 | ;; SET UP THE PANEL PLOT WITH LEFT AND RIGHT COLUMNS |
---|
32 | wks_type = "png" |
---|
33 | ; Open the workspace on the current file |
---|
34 | wks = gsn_open_wks(wks_type,"BAROCLINIC_WAVE_BC_physCompare.png") |
---|
35 | |
---|
36 | lgres = True |
---|
37 | lgres@lgLineColors = (/"Red","Orange","Blue","Green","Black"/) |
---|
38 | lgres@lgLineThicknessF = 5.0 ; legend line thickness |
---|
39 | lgres@lgItemType = "Lines" ; show lines only (default) |
---|
40 | lgres@lgLineDashSegLenF = 0.5 |
---|
41 | lgres@lgMonoDashIndex = True |
---|
42 | lgres@lgLabelFontHeightF = .04 ; set the legend label font thickness |
---|
43 | lgres@vpWidthF = 0.25 ; width of legend (NDC) |
---|
44 | lgres@vpHeightF = 0.15 ; height of legend (NDC) |
---|
45 | lgres@lgPerimColor = "Black" ; draw the box perimeter in orange |
---|
46 | lgres@lgPerimThicknessF = 2.0 ; thicken the box perimeter |
---|
47 | lgres@lgBoxBackground = "White" |
---|
48 | lgres@lgPerimFill = 0 |
---|
49 | lgres@lgPerimFillColor = "White" |
---|
50 | lgres@lgBoxMinorExtentF = 0.25 |
---|
51 | |
---|
52 | plotBCSF = new(2,graphic) |
---|
53 | do hemi = 0,1 |
---|
54 | ; Loop over dry, Kessler, and full physics cases |
---|
55 | do physCase = 0,4 |
---|
56 | infile = (/ physdir(physCase), filename /) |
---|
57 | rfile = str_concat(infile) |
---|
58 | a = addfile(rfile,"r") |
---|
59 | |
---|
60 | ; Get the minimum pressure |
---|
61 | whichVar = "PS" |
---|
62 | if (hemi .eq. 0) then |
---|
63 | PS=a->PS(:,{minLat:eqLat},:) |
---|
64 | end if |
---|
65 | if (hemi .eq. 1) then |
---|
66 | PS=a->PS(:,{eqLat:maxLat},:) |
---|
67 | end if |
---|
68 | PS=PS/100. ; convert from Pa to hPa/mb |
---|
69 | PS@long_name="Surface pressure" |
---|
70 | PS@units="hPa" |
---|
71 | var_vtime = dim_min_n_Wrap(PS,(/1,2/)) ; for now, just find min over lat/lon domain at each time |
---|
72 | |
---|
73 | time = var_vtime&time_counter / 86400 |
---|
74 | |
---|
75 | ; CREATE THE PLOTS AND SEND THEM OUT TO .PDF FILES |
---|
76 | if (physCase .eq. 0) then |
---|
77 | ;; SET THE RESOURCES ATTRIBUTES FOR THE XY PLOTS |
---|
78 | rescom = True |
---|
79 | ; Graphics resources |
---|
80 | rescom@gsnDraw = False |
---|
81 | rescom@gsnFrame = False |
---|
82 | rescom@gsnMaximize = True |
---|
83 | rescom@gsnPaperOrientation = "portrait" |
---|
84 | ; Plot area resources |
---|
85 | rescom@tmXMajorGrid = True |
---|
86 | rescom@tmYMajorGrid = True |
---|
87 | rescom@tmXMajorGridLineDashPattern = 2 |
---|
88 | rescom@tmYMajorGridLineDashPattern = 2 |
---|
89 | rescom@tiMainString = "DCMIP1 "+model+" "+whichVar+" v. time "+hemisphere(hemi) ; add title |
---|
90 | rescom@tiXAxisString = "Time (days)" |
---|
91 | rescom@tiYAxisString = var_vtime@long_name+" ("+var_vtime@units+")" |
---|
92 | rescom@xyLineColor = "Red" |
---|
93 | rescom@xyDashPattern = 0 |
---|
94 | rescom@xyLineThicknessF = 5.0 |
---|
95 | plotBCSF(hemi) = gsn_csm_xy(wks,time,var_vtime,rescom) |
---|
96 | end if |
---|
97 | |
---|
98 | if (physCase .gt. 0) then |
---|
99 | dataid = create "xyData" coordArraysClass defaultapp |
---|
100 | "caYArray": var_vtime |
---|
101 | "caXArray": time |
---|
102 | end create |
---|
103 | |
---|
104 | dataspec = NhlAddData(plotBCSF(hemi),"xyCoordData",dataid) |
---|
105 | |
---|
106 | if (physCase .eq. 1) then |
---|
107 | lcolor = "Orange" |
---|
108 | end if |
---|
109 | if (physCase .eq. 2) then |
---|
110 | lcolor = "Blue" |
---|
111 | end if |
---|
112 | if (physCase .eq. 3) then |
---|
113 | lcolor = "Green" |
---|
114 | end if |
---|
115 | if (physCase .eq. 4) then |
---|
116 | lcolor = "Black" |
---|
117 | end if |
---|
118 | setvalues dataspec |
---|
119 | "xyLineColor" : lcolor |
---|
120 | "xyDashPattern" : 0 |
---|
121 | "xyLineThicknessF" : 5.0 |
---|
122 | end setvalues |
---|
123 | end if |
---|
124 | delete(PS) |
---|
125 | end do |
---|
126 | |
---|
127 | ; Add a legend and draw the plot |
---|
128 | lbid = gsn_create_legend(wks,5,(/ "Dry","Kessler","Full Physics", "Full Physics 1/2 deg", "Full physics 30min" /),lgres) |
---|
129 | amres = True |
---|
130 | amres@amParallelPosF = -0.3 |
---|
131 | amres@amOrthogonalPosF = 0.3 |
---|
132 | annoid1 = gsn_add_annotation(plotBCSF(hemi),lbid,amres) ; add legend to plot |
---|
133 | |
---|
134 | ; Delete some stuff for the next plot in the panel |
---|
135 | delete(lbid) |
---|
136 | delete(amres) |
---|
137 | delete(dataid) |
---|
138 | delete(dataspec) |
---|
139 | delete(annoid1) |
---|
140 | delete(rescom) |
---|
141 | end do |
---|
142 | ; SET PANEL RESOURCES |
---|
143 | resP = True ; modify the panel plot |
---|
144 | resP@gsnMaximize = True ; maximize plots |
---|
145 | resP@gsnPaperOrientation = "landscape" |
---|
146 | ;resP@gsnFrame = False ; don't advance panel plot |
---|
147 | ;resP@amJust = "TopLeft" ; put in top left of plot |
---|
148 | ;resP@gsnPaperMargin = 0.5 |
---|
149 | ;resP@gsnAttachBorderOn = False |
---|
150 | ;resP@gsnPanelLabelBar = True |
---|
151 | ;resP@lbLabelFontHeightF = 0.02 |
---|
152 | ;resP@lbLabelStride = 2 |
---|
153 | ;resP@lbTitleOn = True |
---|
154 | ;resP@lbTitleString = "Potential Temperature (K)" |
---|
155 | ;resP@lbTitleFontHeightF= 0.02 |
---|
156 | ;resP@pmLabelBarOrthogonalPosF = -0.03 |
---|
157 | ;resP@pmLabelBarWidthF = 0.8 |
---|
158 | ;resP@pmLabelBarHeightF = 0.075 |
---|
159 | |
---|
160 | resP@gsnPanelLeft = 0.05 |
---|
161 | resP@gsnPanelRight = 0.95 |
---|
162 | gsn_panel(wks,plotBCSF,(/1,2/),resP) |
---|
163 | end |
---|