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 | load "radialAvg.ncl" |
---|
22 | |
---|
23 | begin |
---|
24 | |
---|
25 | ;======================================================================================= |
---|
26 | ; User options |
---|
27 | ;======================================================================================= |
---|
28 | filename = "output_dcmip2016_regular.nc" |
---|
29 | data_on_constant_Z=False ; is data already on CONSTANT Z surfaces? |
---|
30 | convert_hybridP_to_Z=False ; is data on hybrid pressure levels? |
---|
31 | hasTimeIx=True ; does file have time index? |
---|
32 | timeStep=64 ; If yes, what time index do you want to plot? |
---|
33 | model="DYNAMICO" ; used for mainStr, but also for model specific if statements |
---|
34 | mainStr=model+" Day 8" ; main string for plot titles |
---|
35 | outType="png" ; output format, popular options are x11, png, pdf, eps |
---|
36 | |
---|
37 | Uname="U" ; Variable name of zonal wind (typically "U" or "ua") |
---|
38 | Vname="V" ; Variable name of meridional wind (typically "V" or "ua") |
---|
39 | Rhoname="Rho" ; Variable name of Rho |
---|
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 | ;======================================================================================= |
---|
44 | |
---|
45 | ;======================================================================================= |
---|
46 | ; Other settings, required to remain constant for intercomparison |
---|
47 | ; Generally, don't touch unless experimenting/debugging |
---|
48 | ;======================================================================================= |
---|
49 | stride=1 ; subset in horizontal, leave @ 1 unless debugging |
---|
50 | minLat=0. ; max lat to read in (deg) |
---|
51 | maxLat=80. ; min lat to read in (deg) |
---|
52 | thetaWindMax=60.0 ; max theta wind to plot (m/s) |
---|
53 | radMinMax=18.0 ; min/max for radial plot (m/s) |
---|
54 | tAnomMax=15.0 ; max for tAnom plots (K) |
---|
55 | offAnomDist=5.0 ; offset distance for anom calcs (deg) |
---|
56 | zintmin=20.0 ; min height for z interpolation (m) |
---|
57 | zintmax=18000.0 ; max height for z interpolation (m) |
---|
58 | nzint=100 ; num points for z interpolation |
---|
59 | Zmidpt=2500.0 ; half of Z points BELOW this height, half ABOVE |
---|
60 | ;======================================================================================= |
---|
61 | |
---|
62 | ;======================================================================================= |
---|
63 | ; Input data checks |
---|
64 | ;======================================================================================= |
---|
65 | if (data_on_constant_Z .and. convert_hybridP_to_Z) then |
---|
66 | print("Both data_on_constant_Z and convert_hybridP_to_Z cannot be true, exiting...") |
---|
67 | exit |
---|
68 | end if |
---|
69 | |
---|
70 | ;======================================================================================= |
---|
71 | ; Get file, coordinate variables, and U, V, PSL |
---|
72 | ;======================================================================================= |
---|
73 | print("Loading data from file...") |
---|
74 | |
---|
75 | f = addfile(filename,"r") |
---|
76 | |
---|
77 | lat = f->lat({minLat:maxLat:stride}) |
---|
78 | lon = f->lon(::stride) |
---|
79 | lev = f->lev(:) |
---|
80 | nlat = dimsizes(lat) |
---|
81 | nlon = dimsizes(lon) |
---|
82 | nlev = dimsizes(lev) |
---|
83 | |
---|
84 | if (hasTimeIx) then |
---|
85 | U = f->$Uname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
86 | V = f->$Vname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
87 | T = f->$Tname$(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
88 | if (model .eq. "DYNAMICO") then |
---|
89 | P0 = f->P(timeStep,0,{minLat:maxLat:stride},::stride) |
---|
90 | T0 = f->$Tname$(timeStep,0,{minLat:maxLat:stride},::stride) |
---|
91 | Rho = P0/(287*T0) ; not quite exact (FIXME) |
---|
92 | Rho@units="kg/m3" |
---|
93 | else |
---|
94 | Rho = f->$Rhoname$(timeStep,0,{minLat:maxLat:stride},::stride) |
---|
95 | end if |
---|
96 | PS = f->$PSname$(timeStep,{minLat:maxLat:stride},::stride) |
---|
97 | else |
---|
98 | U = f->$Uname$(:,{minLat:maxLat:stride},::stride) |
---|
99 | V = f->$Vname$(:,{minLat:maxLat:stride},::stride) |
---|
100 | T = f->$Tname$(:,{minLat:maxLat:stride},::stride) |
---|
101 | Rho = f->$Rhoname$(0,{minLat:maxLat:stride},::stride) |
---|
102 | ;PS = f->$PSname$({minLat:maxLat:stride},::stride) |
---|
103 | end if |
---|
104 | ; If U and V are not m/s, please convert here |
---|
105 | U@units="m/s" |
---|
106 | V@units="m/s" |
---|
107 | T@units="K" |
---|
108 | T@long_name="Temperature" |
---|
109 | ;======================================================================================= |
---|
110 | |
---|
111 | ;======================================================================================= |
---|
112 | ; Find center of storm by minimizing PS |
---|
113 | ;======================================================================================= |
---|
114 | print("Finding minimum PS location...") |
---|
115 | a = new((/nlat,nlon/),typeof(PS)) |
---|
116 | a(:,:) = PS(:,:) |
---|
117 | a1D = ndtooned(a) |
---|
118 | dsizes_a = dimsizes(a) |
---|
119 | a_indices = ind_resolve(minind(a1D),dsizes_a) ; Resolve 1D indices to original array |
---|
120 | psminlat = lat(a_indices(0,0)) |
---|
121 | psminlon = lon(a_indices(0,1)) |
---|
122 | delete([/a,a1D,dsizes_a,a_indices/]) |
---|
123 | print("... minimum found at lat: "+psminlat+" lon: "+psminlon) |
---|
124 | ;======================================================================================= |
---|
125 | |
---|
126 | |
---|
127 | ;======================================================================================= |
---|
128 | ; Calculate temperature anomaly |
---|
129 | ;======================================================================================= |
---|
130 | Tanom = T |
---|
131 | Tref = T(:,{psminlat+offAnomDist},{psminlon+offAnomDist}) |
---|
132 | Tanom = T - conform(T,Tref,0) |
---|
133 | Tanom@long_name="Temperature anomaly" |
---|
134 | ;======================================================================================= |
---|
135 | |
---|
136 | |
---|
137 | ;======================================================================================= |
---|
138 | ; Do Z calculations/interpolations if necessary |
---|
139 | ;======================================================================================= |
---|
140 | if (.not. data_on_constant_Z) then |
---|
141 | ;======================================================================================= |
---|
142 | ; Convert from hybrid levels to Z (CAM) |
---|
143 | ;======================================================================================= |
---|
144 | if (convert_hybridP_to_Z) |
---|
145 | print("Converting from hybrid P to Z...") |
---|
146 | ; If need_to_convert_P_to_Z is true, variables T,Q,hyam,hybm,hyai,hybm,P0 needed |
---|
147 | ; all hybrid coefs and TV need to be TOP TO BOTTOM! |
---|
148 | |
---|
149 | hyam=f->hyam |
---|
150 | hybm=f->hybm |
---|
151 | hyai=f->hyai |
---|
152 | hybi=f->hybi |
---|
153 | P0=f->P0 |
---|
154 | Tconv = f->T(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
155 | Qconv = f->Q(timeStep,:,{minLat:maxLat:stride},::stride) |
---|
156 | |
---|
157 | ; Calculate TV from T (K) and Q (kg/kg) |
---|
158 | TV=Tconv*(1.+0.61*Qconv) |
---|
159 | |
---|
160 | ; PHIS is nlatxnlon array = 0.0 for DCMIP Test 2 |
---|
161 | PHIS=Rho(0,:,:) |
---|
162 | PHIS=0.0 |
---|
163 | |
---|
164 | Z = cz2ccm(Rho,PHIS,TV,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1)) |
---|
165 | Z@units="m" |
---|
166 | |
---|
167 | delete([/Tconv,Qconv,TV,PHIS/]) |
---|
168 | end if |
---|
169 | |
---|
170 | ;======================================================================================= |
---|
171 | ; Example of calculating Z at each level |
---|
172 | ; MODELSPEC: (if you do not put output directly on either Z levels OR hybrid levels, you will |
---|
173 | ; need to modify the code here (Tempest shown as example). |
---|
174 | ;======================================================================================= |
---|
175 | if (model .eq. "DYNAMICO") then |
---|
176 | print("DYNAMICO : computing Z from geopotential ...") |
---|
177 | ; NB : geopotential is defined at level interfaces, average from interfaces to full levels |
---|
178 | Phi=f->PHI |
---|
179 | Z=U ; create Z with correct shape |
---|
180 | do kk = 0,nlev-1 |
---|
181 | do ii = 0,nlat-1 |
---|
182 | do jj = 0,nlon-1 |
---|
183 | Z(kk,ii,jj) = (.5/9.81)*(Phi(timeStep,kk,ii,jj)+Phi(timeStep,kk+1,ii,jj)) |
---|
184 | end do |
---|
185 | end do |
---|
186 | end do |
---|
187 | Z@units="m" |
---|
188 | print("...done") |
---|
189 | end if |
---|
190 | |
---|
191 | ;======================================================================================= |
---|
192 | ; If all else fails, try to load Z directly from file |
---|
193 | ;======================================================================================= |
---|
194 | if(.not. isdefined("Z")) then |
---|
195 | if (isfilevar(f, Zname)) then |
---|
196 | print("Found Z on file...") |
---|
197 | Z = f->Z |
---|
198 | else |
---|
199 | print("FATAL: Z needs to be either loaded from a file or otherwise defined before continuing...") |
---|
200 | exit |
---|
201 | end if |
---|
202 | end if |
---|
203 | |
---|
204 | ;======================================================================================= |
---|
205 | ; Generate uniform Z levs to interpolate to |
---|
206 | ; Extra density (half of pts) is specified in PBL, controlled by Zmidpt |
---|
207 | ;======================================================================================= |
---|
208 | print("Generating constant Z levels to interpolate to") |
---|
209 | ZlevPBL = fspan(zintmin,Zmidpt,nzint/2) |
---|
210 | ZlevABL = fspan(Zmidpt,zintmax,nzint/2) |
---|
211 | |
---|
212 | Zlev=new(nzint-1,"float") |
---|
213 | Zlev(0:(nzint/2)-1)=ZlevPBL |
---|
214 | Zlev((nzint/2)-1:nzint-2)=ZlevABL |
---|
215 | Zlev@units = "m" |
---|
216 | Zlev!0 = "lev" |
---|
217 | Zlev&lev = Zlev |
---|
218 | delete([/ZlevPBL,ZlevABL/]) |
---|
219 | ;======================================================================================= |
---|
220 | |
---|
221 | ;======================================================================================= |
---|
222 | ; Interpolate lat/lon variables to constant Z levels |
---|
223 | ;======================================================================================= |
---|
224 | print("Interpolating from Z model levels to constant Z surfaces") |
---|
225 | U_Z = int2p_n_Wrap(Z,U,Zlev,2,0) |
---|
226 | V_Z = int2p_n_Wrap(Z,V,Zlev,2,0) |
---|
227 | Tanom_Z = int2p_n_Wrap(Z,Tanom,Zlev,2,0) |
---|
228 | ;======================================================================================= |
---|
229 | else |
---|
230 | print("Data already on Z surfaces, copying to new vars") |
---|
231 | U_Z = U |
---|
232 | V_Z = V |
---|
233 | Tanom_Z = Tanom |
---|
234 | end if |
---|
235 | |
---|
236 | ;======================================================================================= |
---|
237 | ; Calculate radial and tangential wind components on native grid |
---|
238 | ;======================================================================================= |
---|
239 | print("Separating winds into radial and tangential components given Rho min loc...") |
---|
240 | vComps = calcWindComponents(U_Z,V_Z,lat,lon,psminlat,psminlon) |
---|
241 | v_rad = vComps[0] |
---|
242 | v_theta = vComps[1] |
---|
243 | delete(vComps) |
---|
244 | ;======================================================================================= |
---|
245 | |
---|
246 | |
---|
247 | ;======================================================================================= |
---|
248 | ; Perform radial integrations |
---|
249 | ;======================================================================================= |
---|
250 | print("Do radial integrations...") |
---|
251 | rad_v_theta = radialAvg3D(v_theta,lat,lon,v_theta&lev,psminlat,psminlon,800.,False) |
---|
252 | rad_v_rad = radialAvg3D(v_rad ,lat,lon,v_rad&lev ,psminlat,psminlon,800.,False) |
---|
253 | rad_t_anom = radialAvg3D(Tanom_Z,lat,lon,Tanom_Z&lev ,psminlat,psminlon,800.,False) |
---|
254 | ;======================================================================================= |
---|
255 | |
---|
256 | |
---|
257 | print(v_theta&lev) |
---|
258 | print(rad_v_theta&radius) |
---|
259 | |
---|
260 | print("Plotting...") |
---|
261 | |
---|
262 | res = True |
---|
263 | res@gsnDraw = False |
---|
264 | res@gsnFrame = False |
---|
265 | res@gsnSpreadColors = True ; Span full color map |
---|
266 | res@cnFillOn = True ; Turn on contour fill |
---|
267 | res@cnLinesOn = False |
---|
268 | res@cnLineLabelsOn = False |
---|
269 | res@cnInfoLabelOn = False |
---|
270 | res@tiYAxisString = "Height (m)" |
---|
271 | res@gsnYAxisIrregular2Linear = True |
---|
272 | res@gsnXAxisIrregular2Linear = True |
---|
273 | |
---|
274 | wks = gsn_open_wks (outType,"t_anom_"+model) |
---|
275 | gsn_define_colormap(wks,"matlab_jet") |
---|
276 | res_rad = res |
---|
277 | res_rad@cnLevelSelectionMode = "ExplicitLevels" |
---|
278 | res_rad@cnLevels = ispan(0,toint(tAnomMax),1) |
---|
279 | plot = gsn_csm_contour(wks,rad_t_anom,res_rad) |
---|
280 | draw(plot) |
---|
281 | frame(wks) |
---|
282 | |
---|
283 | delete(plot) |
---|
284 | delete(wks) |
---|
285 | delete(res_rad) |
---|
286 | |
---|
287 | wks = gsn_open_wks (outType,"v_theta_"+model) |
---|
288 | gsn_define_colormap(wks,"BlAqGrYeOrReVi200") |
---|
289 | res_rad = res |
---|
290 | res_rad@tiMainString = mainStr+" tangential wind" |
---|
291 | res_rad@cnLevelSelectionMode = "ExplicitLevels" |
---|
292 | res_rad@cnLevels = fspan(0.0,thetaWindMax,21) |
---|
293 | plot = gsn_csm_contour(wks,rad_v_theta,res_rad) |
---|
294 | draw(plot) |
---|
295 | frame(wks) |
---|
296 | |
---|
297 | delete(plot) |
---|
298 | delete(wks) |
---|
299 | delete(res_rad) |
---|
300 | |
---|
301 | wks = gsn_open_wks (outType,"v_rad_"+model) |
---|
302 | gsn_define_colormap(wks,"hotcolr_19lev") |
---|
303 | res_rad = res |
---|
304 | res_rad@tiMainString = mainStr+" radial wind" |
---|
305 | res_rad@cnLevelSelectionMode = "ExplicitLevels" |
---|
306 | res_rad@cnLevels = ispan(toint(-radMinMax),toint(radMinMax),1) |
---|
307 | plot = gsn_csm_contour(wks,rad_v_rad,res_rad) |
---|
308 | draw(plot) |
---|
309 | frame(wks) |
---|
310 | |
---|
311 | delete(plot) |
---|
312 | delete(wks) |
---|
313 | delete(res_rad) |
---|
314 | |
---|
315 | |
---|
316 | |
---|
317 | end |
---|