-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathAnomalies_Figure_advection_T.ncl
182 lines (138 loc) · 6.58 KB
/
Anomalies_Figure_advection_T.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
begin
;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"
; uq_his=hus_his*ua_his
; uq_rcp=hus_rcp*ua_rcp
; vq_his=hus_his*va_his
; vq_rcp=hus_rcp*va_rcp
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;; Nat-Hist;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
months = (/"Jan.","Feb.","Mar.","Apr.","May","Jun.","Jul.","Aug.","Sept.","Oct.","Nov.","Dec."/)
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++
;++++++++++++ advect varible of temperature ++++++++++++
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++
file_precip_0 = addfile("GPH_T_U_V_500hpa.nc","r")
t = file_precip_0->t (time|:, expver|0, latitude|160:360,longitude|240:720)
t := t(:,::-1,:)
u = file_precip_0->u (time|:, expver|0, latitude|160:360,longitude|240:720)
u := u(:,::-1,:)
v = file_precip_0->v (time|:, expver|0, latitude|160:360,longitude|240:720)
v := v(:,::-1,:)
lat = file_precip_0->latitude(160:360)
lon = file_precip_0->longitude(240:720)
gridType = 0 ; global gaussian grid
opt_adv = 0 ; return only the advected variable
long_name = "advection of temperature"
units = "K/s"
Tadv = advect_variable(short2flt(u), short2flt(v), short2flt(t),gridType,long_name,units,opt_adv)
temp := Tadv(latitude|:,longitude|:,time|0:503)
temp!0 = "lat"
temp!1 = "lon"
temp!2 = "time"
tp_clm := clmMonLLT (short2flt(temp)) ;Calculate LTMM's
Tadv_anomalies = calcMonAnomLLT(short2flt(temp), tp_clm) ;Calculate Anomalies from Means
temp := u(latitude|:,longitude|:,time|0:503)
temp!0 = "lat"
temp!1 = "lon"
temp!2 = "time"
tp_clm := clmMonLLT (short2flt(temp)) ;Calculate LTMM's
u5_anomalies = calcMonAnomLLT(short2flt(temp), tp_clm) ;Calculate Anomalies from Means
temp := v(latitude|:,longitude|:,time|0:503)
temp!0 = "lat"
temp!1 = "lon"
temp!2 = "time"
tp_clm := clmMonLLT (short2flt(temp)) ;Calculate LTMM's
v5_anomalies = calcMonAnomLLT(short2flt(temp), tp_clm) ;Calculate Anomalies from Means
; region_min := dim_min( omega_anomalies)
; region_mean := dim_avg( omega_anomalies)
; region_max := dim_max( omega_anomalies)
region_mean := dim_avg(Tadv)
print (region_mean(100,200))
;*******************************************************************
; create plot for total preciptation and wind_850hpa anomalies
;*******************************************************************
res = True ; plot mods desired
res@gsnDraw = False
res@gsnFrame = False
; res@mpProjection = "CylindricalEqualArea"; choose projection
res@mpFillOn = False ; turn off map fill
res@cnFillOn = True ; turn on color
res@cnLinesOn = False ; turn off contour lines
res@cnFillPalette = "GMT_polar" ; set color map
res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF = -300. ; set min contour level
res@cnMaxLevelValF = 300. ; set max contour level
res@cnLevelSpacingF = 50. ; set contour spacing
res@gsnAddCyclic = False ; regional plot
; res@mpMinLatF = 20 ; min lat to mask
; res@mpMaxLatF = 80 ; max lat to mask
; res@mpMinLonF = -90 ; min lon to mask
; res@mpMaxLonF = 40 ; max lon to mask
res@mpMaxLatF = max(lat) ; choose subregion
res@mpMinLatF = min(lat)
res@mpMaxLonF = 150
res@mpMinLonF = min(lon)
vcres = True
vcres@gsnDraw = False
vcres@gsnFrame = False
vcres@vcGlyphStyle = "CurlyVector" ; curly vectors
;---Set up some vector resources.
vcres@vcLevelSelectionMode = "ManualLevels"
vcres@vcMinLevelValF = -1.5
vcres@vcMaxLevelValF = 1.5
vcres@vcLevelSpacingF = 0.5 ;
; vcres@vcLevelPalette = "amwg_blueyellowred" ; assign color map to vectors
; vcres@lbLabelBarOn = False
;---Vector lengths and color
; vcres@vcFillArrowsOn = True
; vcres@vcLineArrowThicknessF = 2.0
vcres@vcMinFracLengthF = 0.33
vcres@vcMinMagnitudeF = 10
vcres@vcMinDistanceF = 0.03
; vcres@vcMonoFillArrowFillColor = False
; vcres@vcMonoLineArrowColor = False
vcres@vcRefLengthF = 0.045
vcres@vcRefMagnitudeF = 3.0
vcres@vcRefAnnoOrthogonalPosF = -0.12
vcres@vcRefAnnoParallelPosF = 0.997
vcres@vcRefAnnoFontHeightF = 0.015
vcres@lbTitleString = "wind flux"
vcres@lbTitleOffsetF = -0.25
vcres@lbTitleFontHeightF = 0.02
vcres@lbLabelFontHeightF = 0.015
vcres@lbLabelAutoStride = True
;---Make sure vectors are drawn in "Postdraw" phase.
vcres@vcVectorDrawOrder = "Postdraw"
;*********************************************************************
; create plot for tempature advction and wind_500hpa anomalies
;*********************************************************************
res@gsnDraw = False
res@cnMinLevelValF = -2.e-5 ; set min contour level
res@cnMaxLevelValF = 2.e-5 ; set max contour level
res@cnLevelSpacingF = 5.e-6 ; set contour spacing
vcres@vcRefLengthF = 0.045
vcres@vcRefMagnitudeF = 4.0
vcres@vcRefAnnoOrthogonalPosF = -0.12
vcres@vcRefAnnoParallelPosF = 0.997
vcres@vcRefAnnoFontHeightF = 0.015
wks3 = gsn_open_wks("pdf","Temperature_advection_wind_500hpa_anomalies") ; send graphics to PNG file
gsn_define_colormap(wks3, "GMT_polar")
; gsn_reverse_colormap(wks3)
plot = new(12, graphic)
do month = 0, 11
res@tiMainString = months(month) + "(mm/day)"
plotA= gsn_csm_contour_map(wks3, Tadv_anomalies (:,:,month+492), res ) ; draw second plot
plotB = gsn_csm_vector(wks3,u5_anomalies(:,:,month+492), v5_anomalies(:,:,month+492),vcres) ; create a default plot
overlay(plotA, plotB)
plot(month) = plotA
; maximize_output(wks3,True)
end do
resP = True
resP@gsnMaximize = True ; maximize plots
resP@gsnPanelLabelBar = True
resP@lbAutoManage = False
resP@lbLabelFontHeightF = 0.0075
gsn_panel(wks3,plot,(/4,3/), resP)
end;