1 | ;======================================================================================= |
---|
2 | ; This NCL code calculates radially-averaged tangential and radial wind components |
---|
3 | ; as well as T anomaly for DCMIP test case #2 (cyclone) |
---|
4 | ; this code requires the accompanying function set "radialAvg.ncl" |
---|
5 | ; |
---|
6 | ; Usage: User should modify "user options" for their particular data set. Currently, |
---|
7 | ; U, V, T, PS are required as variables. |
---|
8 | ; If variables are on constant Z surfaces, life is easy. |
---|
9 | ; |
---|
10 | ; Grepping for "MODELSPEC" will point to possible areas of the code in need of modification |
---|
11 | ; for model specific output |
---|
12 | ; |
---|
13 | ; Written by Colin Zarzycki (zarzycki@ucar.edu) |
---|
14 | ; Version 0.1 (6/5/2016) - DCMIP-2016 release |
---|
15 | ;======================================================================================= |
---|
16 | |
---|
17 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
---|
18 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
---|
19 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
---|
20 | load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" |
---|
21 | |
---|
22 | begin |
---|
23 | |
---|
24 | ;======================================================================================= |
---|
25 | ; User options |
---|
26 | ;======================================================================================= |
---|
27 | ;filename = "/glade/scratch/zarzycki/cam5_work/dcmip2/kessler_bryan/dcmip2.cam.h0.2000-01-05-00000.nc" |
---|
28 | |
---|
29 | filename="output_dcmip2016_regular.nc" |
---|
30 | data_on_constant_Z=False ; is data already on CONSTANT Z surfaces? |
---|
31 | convert_hybridP_to_Z=False ; is data on hybrid pressure levels? |
---|
32 | hasTimeIx=True ; does file have time index? |
---|
33 | timeStep=79 ; If yes, what time index do you want to plot? |
---|
34 | model="DYNAMICO" ; used for mainStr, but also for model specific if statements |
---|
35 | mainStr=model+" Day 8" ; main string for plot titles |
---|
36 | out_type="png" ; output format, popular options are x11, png, pdf, eps |
---|
37 | |
---|
38 | Uname="U" ; Variable name of zonal wind (typically "U" or "ua") |
---|
39 | Vname="V" ; Variable name of meridional wind (typically "V" or "ua") |
---|
40 | PSname="PS" ; Variable name of surface pressure (typically "PS","PSL","ps",or "slp") |
---|
41 | Tname="T" ; Variable name of air temperature (typically "T" or "ta") |
---|
42 | Zname="Z" ; Variable name of height (typically "Z" or "za") |
---|
43 | PRECTname="PRECL" ; Variable name of height (typically "Z" or "za") |
---|
44 | |
---|
45 | ;======================================================================================= |
---|
46 | |
---|
47 | ;======================================================================================= |
---|
48 | ; Other settings, required to remain constant for intercomparison |
---|
49 | ; Generally, don't touch unless experimenting/debugging |
---|
50 | ;======================================================================================= |
---|
51 | stride=1 ; subset in horizontal, leave @ 1 unless debugging |
---|
52 | minLat=0. ; max lat to read in (deg) |
---|
53 | maxLat=80. ; min lat to read in (deg) |
---|
54 | thetaWindMax=60.0 ; max theta wind to plot (m/s) |
---|
55 | radMinMax=18.0 ; min/max for radial plot (m/s) |
---|
56 | tAnomMax=15.0 ; max for tAnom plots (K) |
---|
57 | offAnomDist=5.0 ; offset distance for anom calcs (deg) |
---|
58 | zintmin=20.0 ; min height for z interpolation (m) |
---|
59 | zintmax=18000.0 ; max height for z interpolation (m) |
---|
60 | nzint=100 ; num points for z interpolation |
---|
61 | Zmidpt=2500.0 ; half of Z points BELOW this height, half ABOVE |
---|
62 | ;======================================================================================= |
---|
63 | |
---|
64 | ;======================================================================================= |
---|
65 | ; Input data checks |
---|
66 | ;======================================================================================= |
---|
67 | if (data_on_constant_Z .and. convert_hybridP_to_Z) then |
---|
68 | print("Both data_on_constant_Z and convert_hybridP_to_Z cannot be true, exiting...") |
---|
69 | exit |
---|
70 | end if |
---|
71 | |
---|
72 | ;======================================================================================= |
---|
73 | ; Get file, coordinate variables, and U, V, PSL |
---|
74 | ;======================================================================================= |
---|
75 | print("Loading data from file...") |
---|
76 | |
---|
77 | f = addfile(filename,"r") |
---|
78 | |
---|
79 | lat = f->lat({minLat:maxLat:stride}) |
---|
80 | lon = f->lon(::stride) |
---|
81 | lev = f->lev(:) |
---|
82 | nlat = dimsizes(lat) |
---|
83 | nlon = dimsizes(lon) |
---|
84 | nlev = dimsizes(lev) |
---|
85 | |
---|
86 | if (hasTimeIx) then |
---|
87 | U = f->$Uname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
88 | V = f->$Vname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
89 | T = f->$Tname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
90 | PS = f->$PSname$(timeStep,{minLat:maxLat:stride},::stride) |
---|
91 | PRECT = f->$PRECTname$(timeStep,{minLat:maxLat:stride},::stride) |
---|
92 | else |
---|
93 | U = f->$Uname$(:,{minLat:maxLat:stride},::stride) |
---|
94 | V = f->$Vname$(:,{minLat:maxLat:stride},::stride) |
---|
95 | T = f->$Tname$(:,{minLat:maxLat:stride},::stride) |
---|
96 | PS = f->$PSname$({minLat:maxLat:stride},::stride) |
---|
97 | PRECT = f->$PRECTname$({minLat:maxLat:stride},::stride) |
---|
98 | end if |
---|
99 | ; If U and V are not m/s, please convert here |
---|
100 | U@units="m/s" |
---|
101 | V@units="m/s" |
---|
102 | T@units="K" |
---|
103 | T@long_name="Temperature" |
---|
104 | ;======================================================================================= |
---|
105 | |
---|
106 | ;======================================================================================= |
---|
107 | ; Find center of storm by minimizing PS |
---|
108 | ;======================================================================================= |
---|
109 | print("Finding minimum PS location...") |
---|
110 | a = new((/nlat,nlon/),typeof(PS)) |
---|
111 | a(:,:) = PS(:,:) |
---|
112 | a1D = ndtooned(a) |
---|
113 | dsizes_a = dimsizes(a) |
---|
114 | a_indices = ind_resolve(minind(a1D),dsizes_a) ; Resolve 1D indices to original array |
---|
115 | psminlat = lat(a_indices(0,0)) |
---|
116 | psminlon = lon(a_indices(0,1)) |
---|
117 | delete([/a,a1D,dsizes_a,a_indices/]) |
---|
118 | print("... PS minimum found at lat: "+psminlat+" lon: "+psminlon) |
---|
119 | ;======================================================================================= |
---|
120 | |
---|
121 | |
---|
122 | ;======================================================================================= |
---|
123 | ; Calculate temperature anomaly |
---|
124 | ;======================================================================================= |
---|
125 | Tanom = T |
---|
126 | Tref = T(:,{psminlat+offAnomDist},{psminlon+offAnomDist}) |
---|
127 | Tanom = T - conform(T,Tref,0) |
---|
128 | Tanom@long_name="Temperature anomaly" |
---|
129 | ;======================================================================================= |
---|
130 | |
---|
131 | |
---|
132 | ;======================================================================================= |
---|
133 | ; Do Z calculations/interpolations if necessary |
---|
134 | ;======================================================================================= |
---|
135 | if (.not. data_on_constant_Z) then |
---|
136 | ;======================================================================================= |
---|
137 | ; Convert from hybrid levels to Z (CAM) |
---|
138 | ;======================================================================================= |
---|
139 | if (convert_hybridP_to_Z) |
---|
140 | print("Converting from hybrid P to Z...") |
---|
141 | ; If need_to_convert_P_to_Z is true, variables T,Q,hyam,hybm,hyai,hybm,P0 needed |
---|
142 | ; all hybrid coefs and TV need to be TOP TO BOTTOM! |
---|
143 | |
---|
144 | hyam=f->hyam |
---|
145 | hybm=f->hybm |
---|
146 | hyai=f->hyai |
---|
147 | hybi=f->hybi |
---|
148 | P0=f->P0 |
---|
149 | Tconv = f->T(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
150 | Qconv = f->Q(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
151 | |
---|
152 | ; Calculate TV from T (K) and Q (kg/kg) |
---|
153 | TV=Tconv*(1.+0.61*Qconv) |
---|
154 | |
---|
155 | ; PHIS is nlatxnlon array = 0.0 for DCMIP Test 2 |
---|
156 | PHIS=PS |
---|
157 | PHIS=0.0 |
---|
158 | |
---|
159 | Z = cz2ccm(PS,PHIS,TV,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1)) |
---|
160 | Z@units="m" |
---|
161 | |
---|
162 | delete([/Tconv,Qconv,TV,PHIS/]) |
---|
163 | end if |
---|
164 | |
---|
165 | ;======================================================================================= |
---|
166 | ; Example of calculating Z at each level |
---|
167 | ; MODELSPEC: (if you do not put output directly on either Z levels OR hybrid levels, you will |
---|
168 | ; need to modify the code here (Tempest shown as example). |
---|
169 | ;======================================================================================= |
---|
170 | if (model .eq. "tempest") then |
---|
171 | Zs=f->Zs |
---|
172 | Ztop=30000.0 |
---|
173 | Z = U |
---|
174 | do kk = 0,nlev-1 |
---|
175 | do ii = 0,nlat-1 |
---|
176 | do jj = 0,nlon-1 |
---|
177 | Z(kk,ii,jj) = Zs(ii,jj) + lev(kk)*(Ztop - Zs(ii,jj)) |
---|
178 | end do |
---|
179 | end do |
---|
180 | end do |
---|
181 | Z@units="m" |
---|
182 | end if |
---|
183 | |
---|
184 | if (model .eq. "DYNAMICO") then |
---|
185 | print("DYNAMICO : computing Z from geopotential ...") |
---|
186 | ; NB : geopotential is defined at level interfaces, average from interfaces to full levels |
---|
187 | Phi=f->PHI |
---|
188 | Z=U ; create Z with correct shape |
---|
189 | do kk = 0,nlev-1 |
---|
190 | do ii = 0,nlat-1 |
---|
191 | do jj = 0,nlon-1 |
---|
192 | Z(kk,ii,jj) = (.5/9.81)*(Phi(timeStep,kk,ii,jj)+Phi(timeStep,kk+1,ii,jj)) |
---|
193 | end do |
---|
194 | end do |
---|
195 | end do |
---|
196 | Z@units="m" |
---|
197 | print("...done") |
---|
198 | end if |
---|
199 | |
---|
200 | |
---|
201 | ;======================================================================================= |
---|
202 | ; If all else fails, try to load Z directly from file |
---|
203 | ;======================================================================================= |
---|
204 | if(.not. isdefined("Z")) then |
---|
205 | if (isfilevar(f, Zname)) then |
---|
206 | print("Found Z on file...") |
---|
207 | Z = f->Z |
---|
208 | else |
---|
209 | print("FATAL: Z needs to be either loaded from a file or otherwise defined before continuing...") |
---|
210 | exit |
---|
211 | end if |
---|
212 | end if |
---|
213 | ;======================================================================================= |
---|
214 | |
---|
215 | end if |
---|
216 | |
---|
217 | ;======================================================================================= |
---|
218 | ; Generate Z levels of interest |
---|
219 | ;======================================================================================= |
---|
220 | print("Generating constant Z levels to interpolate to") |
---|
221 | ZlevPBL = fspan(zintmin,Zmidpt,nzint/2) |
---|
222 | ZlevABL = fspan(Zmidpt,zintmax,nzint/2) |
---|
223 | |
---|
224 | Zlev=(/100.,1000.,1500.,2500.,5000.,10000.,15000./) |
---|
225 | Zlev@units = "m" |
---|
226 | Zlev!0 = "lev" |
---|
227 | Zlev&lev = Zlev |
---|
228 | |
---|
229 | ;======================================================================================= |
---|
230 | ; Interpolate lat/lon variables to constant Z levels |
---|
231 | ;======================================================================================= |
---|
232 | print("Interpolating to Z surfaces") |
---|
233 | U_Z = int2p_n_Wrap(Z,U,Zlev,2,0) |
---|
234 | V_Z = int2p_n_Wrap(Z,V,Zlev,2,0) |
---|
235 | Tanom_Z = int2p_n_Wrap(Z,Tanom,Zlev,2,0) |
---|
236 | ;======================================================================================= |
---|
237 | |
---|
238 | WIND_Z = U_Z |
---|
239 | WIND_Z = sqrt(U_Z^2+V_Z^2) |
---|
240 | WIND_Z@long_name="Horizontal wind" |
---|
241 | |
---|
242 | PRECT = PRECT*8.64e7 ; convert m/s to mm/day |
---|
243 | PRECT@long_name="Precipitation rate" |
---|
244 | PRECT@units="mm/d" |
---|
245 | |
---|
246 | print("Plotting...") |
---|
247 | |
---|
248 | wks = gsn_open_wks (out_type,"horiz_crossx") ; send graphics to PNG file |
---|
249 | contour = new(4,"graphic") |
---|
250 | |
---|
251 | gsn_define_colormap(wks,"BlAqGrYeOrReVi200") |
---|
252 | |
---|
253 | res = True |
---|
254 | res@gsnDraw = False |
---|
255 | res@gsnFrame = False |
---|
256 | res@gsnSpreadColors = True ; Span full color map |
---|
257 | res@cnFillOn = True ; Turn on contour fill |
---|
258 | res@cnLinesOn = False |
---|
259 | res@cnLineLabelsOn = False |
---|
260 | res@cnInfoLabelOn = False |
---|
261 | res@gsnAddCyclic = True |
---|
262 | |
---|
263 | latWidth = 10.0 |
---|
264 | lonWidth = 10.0 |
---|
265 | res@mpOutlineOn = False |
---|
266 | res@mpMaxLatF = psminlat + latWidth |
---|
267 | res@mpMinLatF = psminlat - latWidth |
---|
268 | res@mpMinLonF = psminlon - lonWidth |
---|
269 | res@mpMaxLonF = psminlon + lonWidth |
---|
270 | |
---|
271 | res@cnLevelSelectionMode = "ManualLevels" |
---|
272 | res@cnLevelSpacingF = 5.0 |
---|
273 | res@cnMinLevelValF = 10.0 |
---|
274 | res@cnMaxLevelValF = 80.0 |
---|
275 | res@tiMainString="1500m Horiz. Wind" |
---|
276 | contour(0) = gsn_csm_contour_map(wks,WIND_Z({1500.0},:,:),res) ; create the plot |
---|
277 | |
---|
278 | res@cnLevelSelectionMode = "ManualLevels" |
---|
279 | res@cnLevelSpacingF = 10.0 |
---|
280 | res@cnMinLevelValF = -60.0 |
---|
281 | res@cnMaxLevelValF = 60.0 |
---|
282 | res@tiMainString="100m U-Wind" |
---|
283 | contour(1) = gsn_csm_contour_map(wks,U_Z({100.0},:,:),res) ; create the plot |
---|
284 | |
---|
285 | res@cnLevelSelectionMode = "ManualLevels" |
---|
286 | res@cnLevelSpacingF = 1.0 |
---|
287 | res@cnMinLevelValF = 0.0 |
---|
288 | res@cnMaxLevelValF = 10.0 |
---|
289 | res@tiMainString="5000m T anom." |
---|
290 | contour(2) = gsn_csm_contour_map(wks,Tanom_Z({5000.0},:,:),res) ; create the plot |
---|
291 | |
---|
292 | res@cnLevelSelectionMode = "ManualLevels" |
---|
293 | res@cnLevelSpacingF = 200.0 |
---|
294 | res@cnMinLevelValF = 200.0 |
---|
295 | res@cnMaxLevelValF = 4000.0 |
---|
296 | res@tiMainString="Precip. rate" |
---|
297 | contour(3) = gsn_csm_contour_map(wks,PRECT(:,:),res) ; create the plot |
---|
298 | |
---|
299 | resP = True ; modify the panel plot |
---|
300 | |
---|
301 | gsn_panel(wks,contour,(/2,2/),resP) ; now draw as one plot |
---|
302 | |
---|
303 | end |
---|