source: configs/dcmip2016/PLOTS/testcase/CYCLONE/include/radialAvg.ncl @ 438

Last change on this file since 438 was 437, checked in by dubos, 8 years ago

Plots for DCMIP2016
File size: 12.8 KB
Line 
1; This file contains various functions relating to TC plotting in radial space
2; Colin Zarzycki (zarzycki@ucar.edu)
3; History:
4; 6/5/2016 - Released for DCMIP-2016
5;
6; contains: radialAvg3D, radialAvg2D, calcWindComponents, radialInt3D, calc_vertmax_from_rad
7;
8;
9; radialAvg3D
10; Usage: radialAvg3D(data,lat,lon,lev,psminlat,psminlon,outerRad,mergeInnerBins)
11; data (numeric): 3-D (lev,lat,lon) of variable to transform to radial space
12; lat (numeric): 1-D array of lats (deg)
13; lon (numeric): 1-D array of lons (deg)
14; lev (numeric): 1-D array of levs (Pa or m)
15; psminlat (numeric): lat of TC center (deg)
16; psminlon (numeric): lon of TC center (deg)
17; outerRad (numeric) : outermost radius to calculate (km) (typically 500-800km)
18; mergeInnerBins (logical) :merge bins at storm center for sampling purpose? (generally True)
19;
20; ex: radialAvg3D(v_theta,lat,lon,v_theta&lev,psminlat,psminlon,500.,True)
21;
22;
23; radialAvg2D
24; Usage: radialAvg2D(data,lat,lon,psminlat,psminlon,outerRad,mergeInnerBins)
25; data (numeric): 2-D (lat,lon) of variable to transform to radial space
26; lat (numeric): 1-D array of lats (deg)
27; lon (numeric): 1-D array of lons (deg)
28; psminlat (numeric): lat of TC center (deg)
29; psminlon (numeric): lon of TC center (deg)
30; outerRad (numeric) : outermost radius to calculate (km) (typically 500-800km)
31; mergeInnerBins (logical) :merge bins at storm center for sampling purpose? (generally True)
32;
33; ex:
34;
35;
36; calcWindComponents
37; Usage: radialAvg3D(U,V,lat,lon,psminlat,psminlon)
38; U (numeric): 2-D (lat,lon) or 3-D (lev, lat, lon) or 4-D (time, lev, lat, lon) zonal wind
39; V (numeric): 2-D (lat,lon) or 3-D (lev, lat, lon) or 4-D (time, lev, lat, lon) zonal wind
40; lat (numeric): 1-D array of lats (deg)
41; lon (numeric): 1-D array of lons (deg)
42; psminlat (numeric): lat of TC center (deg)
43; psminlon (numeric): lon of TC center (deg)
44;
45; ex: calcWindComponents(U_Z,V_Z,lat,lon,psminlat,psminlon)
46;
47;
48; radialInt3D
49; Usage: radialInt3D(data,lat,lon,psminlat,psminlon,outerRad)
50; data (numeric): 3-D (lev,lat,lon)
51; lat (numeric): 1-D array of lats (deg)
52; lon (numeric): 1-D array of lons (deg)
53; lev (numeric): 1-D array of levs (Pa or m)
54; psminlat (numeric): lat of TC center (deg)
55; psminlon (numeric): lon of TC center (deg)
56; outerRad (numeric) : outermost radius to integrate (km) (typically 500-800km)
57;
58; ex: radialInt3D(tofloat(DTCOND),tofloat(lat),tofloat(lon),tofloat(lev),psminlat,psminlon,radius_3D)
59
60undef ( "radialAvg3D" )
61function radialAvg3D ( \
62  data[*][*][*]    : numeric, \
63  lat[*]           : numeric, \
64  lon[*]           : numeric, \
65  lev[*]           : numeric, \
66  psminlat         : numeric, \
67  psminlon         : numeric, \
68  outerRad         : numeric, \
69  mergeInnerBins   : logical \
70)
71local None
72
73; lat (deg)
74; lon (deg)
75; center_lat(deg)
76; center_lon(deg)
77; outerRad (km)
78
79begin
80
81deltaLat = lat(1)-lat(0)
82deltaLon = lon(1)-lon(0)
83deltaMax = max((/deltaLon, deltaLat/))
84deltaMin = min((/deltaLon, deltaLat/))
85kmInDeg = 111.32
86kmGrid = kmInDeg * deltaMax
87print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
88
89nlat = dimsizes(lat)
90nlon = dimsizes(lon)
91nlev = dimsizes(lev)
92
93pi = 3.14159
94d2r = pi/180.
95r2d = 180./pi
96lonr = lon*d2r
97latr = lat*d2r
98copy_VarCoords(lon,lonr)
99copy_VarCoords(lat,latr)
100
101psminLocY = ind_nearest_coord(psminlat,lat,0)
102psminLocX = ind_nearest_coord(psminlon,lon,0)
103
104;;;;; DO RADIAL AVERAGES
105
106loopWidth = toint(ceil((outerRad*1.10) / (kmInDeg * deltaMin)))
107; for standard lat/lon grid
108;1.25 gives n ~= 5 in inner grid if numMerge = 2
109;1.4 gives n ~= 9 in inner grid if numMerge = 2
110timesGrid = 1.1 ; we want each radius bin to be timeGrid times kmGrid
111
112nx = toint(outerRad/(timesGrid*kmGrid))
113print("Number of bins is equal to "+nx)
114 
115if (mergeInnerBins) then
116;  origRadiusArr = fspan(0,outerRad,nx)
117;  radiusArr = new(dimsizes(origRadiusArr)-1,double)
118;  radiusArr(0) = origRadiusArr(0)
119;  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(2:dimsizes(origRadiusArr)-1)
120  numMerge = 2 ; numMerge is the number of innermost radial bins we want to merge
121  print("Merging innermost "+numMerge+" bins because radial average is so small.")
122  numMergeMinusOne = numMerge-1
123  origRadiusArr = fspan(0,outerRad,nx+numMergeMinusOne)
124  radiusArr = new((dimsizes(origRadiusArr)-numMergeMinusOne),double)
125  radiusArr(0) = origRadiusArr(0)
126  print("Going from "+dimsizes(origRadiusArr)+" to "+dimsizes(radiusArr)+" bins...")
127  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(numMerge:dimsizes(origRadiusArr)-1)
128else
129  print("Not merging any innermost bins -- be careful that your inner bins have > 1 pt.")
130  radiusArr = fspan(0,outerRad,nx)
131end if
132
133numRadBins = dimsizes(radiusArr)
134
135rad_thevar_hit  = new((/nlev,numRadBins/),integer)
136rad_thevar_cum  = new((/nlev,numRadBins/),double)
137
138; NEEDS to be 1 or 0
139rad_thevar_hit!1 = "radius"
140rad_thevar_cum!1 = "radius"
141rad_thevar_hit&radius = radiusArr
142rad_thevar_cum&radius = radiusArr
143
144rad_thevar_hit!0 = "lev"
145rad_thevar_cum!0 = "lev"
146rad_thevar_hit&lev = lev
147rad_thevar_cum&lev = lev
148
149rad_thevar_hit = 0
150rad_thevar_cum = 0
151
152; NEED TO CORRECT NLEV
153print("Starting loop")
154do i = psminLocY-loopWidth,psminLocY+loopWidth
155  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
156  do j = psminLocX-loopWidth,psminLocX+loopWidth
157    do k = 0,dimsizes(lev)-1
158      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
159      if (gcdist .le. outerRad)
160        rad_thevar_hit(k,{gcdist}) = rad_thevar_hit(k,{gcdist}) + 1
161        rad_thevar_cum(k,{gcdist}) = rad_thevar_cum(k,{gcdist}) + data(k,i,j)
162      end if
163    end do
164  end do
165end do
166
167print("Minimum number of hits per gridbox: "+min(rad_thevar_hit))
168print("Maximum number of hits per gridbox: "+max(rad_thevar_hit))
169rad_thevar = rad_thevar_cum/rad_thevar_hit
170
171copy_VarCoords(rad_thevar_hit,rad_thevar)
172
173rad_thevar@units = data@units
174rad_thevar@long_name=data@long_name
175
176return(rad_thevar)
177
178end
179
180; ====================================================================================================
181
182undef ( "radialAvg2D" )
183function radialAvg2D ( \
184  data[*][*]       : numeric, \
185  lat[*]           : numeric, \
186  lon[*]           : numeric, \
187  psminlat         : numeric, \
188  psminlon         : numeric, \
189  outerRad         : numeric, \
190  mergeInnerBins   : logical \
191)
192local None
193
194; lat (deg)
195; lon (deg)
196; center_lat(deg)
197; center_lon(deg)
198; outerRad (km)
199
200begin
201
202deltaLat = lat(1)-lat(0)
203deltaLon = lon(1)-lon(0)
204deltaMax = max((/deltaLon, deltaLat/))
205deltaMin = min((/deltaLon, deltaLat/))
206kmInDeg = 111.32
207kmGrid = kmInDeg * deltaMax
208print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
209
210nlat = dimsizes(lat)
211nlon = dimsizes(lon)
212
213pi = 3.14159
214d2r = pi/180.
215r2d = 180./pi
216lonr = lon*d2r
217latr = lat*d2r
218copy_VarCoords(lon,lonr)
219copy_VarCoords(lat,latr)
220
221psminLocY = ind_nearest_coord(psminlat,lat,0)
222psminLocX = ind_nearest_coord(psminlon,lon,0)
223
224;;;;; DO RADIAL AVERAGES
225
226loopWidth = toint(ceil((outerRad*1.10) / (kmInDeg * deltaMin)))
227; for standard lat/lon grid
228;1.25 gives n ~= 5 in inner grid if numMerge = 2
229;1.4 gives n ~= 9 in inner grid if numMerge = 2
230timesGrid = 1.1 ; we want each radius bin to be timeGrid times kmGrid
231
232nx = toint(outerRad/(timesGrid*kmGrid))
233print("Number of bins is equal to "+nx)
234 
235if (mergeInnerBins) then
236;  origRadiusArr = fspan(0,outerRad,nx)
237;  radiusArr = new(dimsizes(origRadiusArr)-1,double)
238;  radiusArr(0) = origRadiusArr(0)
239;  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(2:dimsizes(origRadiusArr)-1)
240  numMerge = 2 ; numMerge is the number of innermost radial bins we want to merge
241  print("Merging innermost "+numMerge+" bins because radial average is so small.")
242  numMergeMinusOne = numMerge-1
243  origRadiusArr = fspan(0,outerRad,nx+numMergeMinusOne)
244  radiusArr = new((dimsizes(origRadiusArr)-numMergeMinusOne),double)
245  radiusArr(0) = origRadiusArr(0)
246  print("Going from "+dimsizes(origRadiusArr)+" to "+dimsizes(radiusArr)+" bins...")
247  radiusArr(1:dimsizes(radiusArr)-1) = origRadiusArr(numMerge:dimsizes(origRadiusArr)-1)
248else
249  print("Not merging any innermost bins -- be careful that your inner bins have > 1 pt.")
250  radiusArr = fspan(0,outerRad,nx)
251end if
252
253numRadBins = dimsizes(radiusArr)
254
255rad_thevar_hit  = new((/numRadBins/),integer)
256rad_thevar_cum  = new((/numRadBins/),double)
257
258; NEEDS to be 1 or 0
259rad_thevar_hit!0 = "radius"
260rad_thevar_cum!0 = "radius"
261rad_thevar_hit&radius = radiusArr
262rad_thevar_cum&radius = radiusArr
263
264rad_thevar_hit = 0
265rad_thevar_cum = 0
266
267; NEED TO CORRECT NLEV
268print("Starting loop")
269do i = psminLocY-loopWidth,psminLocY+loopWidth
270  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
271  do j = psminLocX-loopWidth,psminLocX+loopWidth
272      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
273      if (gcdist .le. outerRad)
274        rad_thevar_hit({gcdist}) = rad_thevar_hit({gcdist}) + 1
275        rad_thevar_cum({gcdist}) = rad_thevar_cum({gcdist}) + data(i,j)
276      end if
277  end do
278end do
279
280print("Minimum number of hits per gridbox: "+min(rad_thevar_hit))
281print("Maximum number of hits per gridbox: "+max(rad_thevar_hit))
282rad_thevar = rad_thevar_cum/rad_thevar_hit
283
284copy_VarCoords(rad_thevar_hit,rad_thevar)
285
286rad_thevar@units = data@units
287rad_thevar@long_name=data@long_name
288
289return(rad_thevar)
290
291end
292
293; ====================================================================================================
294
295undef ( "calcWindComponents" )
296function calcWindComponents ( \
297  U[*][*][*]       : numeric, \
298  V[*][*][*]       : numeric, \
299  lat[*]           : numeric, \
300  lon[*]           : numeric, \
301  psminlat         : numeric, \
302  psminlon         : numeric \
303)
304local None
305begin
306
307pi = 3.14159
308d2r = pi/180.
309r2d = 180./pi
310lonr = lon*d2r
311latr = lat*d2r
312copy_VarCoords(lon,lonr)
313copy_VarCoords(lat,latr)
314
315deltalong = lonr(:)-(psminlon*d2r)
316
317arr1 = conform(U,sin(deltalong),2) * conform(U,cos(latr),1)
318arr2 = cos(psminlat*d2r)*conform(U,sin(latr),1) - sin(psminlat*d2r)*conform(U,cos(latr),1)*conform(U,cos(deltalong),2)
319dir_angle_r = atan2(arr1,arr2)
320; need to conform here
321dir_angle_d = dir_angle_r*r2d
322
323copy_VarCoords(U,dir_angle_r)
324copy_VarCoords(U,dir_angle_d)
325
326eps = 10.e-6
327phi_r = atan2(U,V)
328phi_d = phi_r*r2d
329
330copy_VarCoords(U,phi_r)
331copy_VarCoords(U,phi_d)
332
333WIND = U
334WIND = sqrt(U^2+V^2)
335
336v_theta = WIND*sin(dir_angle_r - phi_r)
337v_rad = -WIND*cos(dir_angle_r - phi_r)
338
339copy_VarCoords(U,v_theta)
340copy_VarCoords(U,v_rad)
341
342v_theta@units = U@units
343v_theta@long_name="Tangential wind"
344v_rad@units = U@units
345v_rad@long_name="Radial wind"
346
347return( [/ v_rad, v_theta /] )
348
349end
350
351
352undef ( "radialInt3D" )
353function radialInt3D ( \
354  data[*][*][*]    : numeric, \
355  lat[*]           : numeric, \
356  lon[*]           : numeric, \
357  lev[*]           : numeric, \
358  psminlat         : numeric, \
359  psminlon         : numeric, \
360  outerRad         : numeric \
361)
362local None
363
364; lat (deg)
365; lon (deg)
366; center_lat(deg)
367; center_lon(deg)
368; outerRad (km)
369
370begin
371
372deltaLat = lat(1)-lat(0)
373deltaLon = lon(1)-lon(0)
374deltaMax = max((/deltaLon, deltaLat/))
375deltaMin = min((/deltaLon, deltaLat/))
376kmInDeg = 111.32
377kmGrid = kmInDeg * deltaMax
378print("The max lat/lon km grid spacing at equator is "+kmGrid+" km")
379
380nlat = dimsizes(lat)
381nlon = dimsizes(lon)
382nlev = dimsizes(lev)
383
384pi = 3.14159
385d2r = pi/180.
386r2d = 180./pi
387lonr = lon*d2r
388latr = lat*d2r
389copy_VarCoords(lon,lonr)
390copy_VarCoords(lat,latr)
391
392psminLocY = ind_nearest_coord(psminlat,lat,0)
393psminLocX = ind_nearest_coord(psminlon,lon,0)
394
395; Get cosine weights
396rad    = 4.0*atan(1.0)/180.0
397clat   = cos(lat*rad)
398
399loopWidth = toint(ceil((outerRad*1.25) / (kmInDeg * deltaMin)))
400; make sure we are more than 1 to do a few rows of buffer cells
401
402rad_thevar_hit  = new((/nlev/),float)
403rad_thevar_cum  = new((/nlev/),float)
404
405rad_thevar_hit!0 = "lev"
406rad_thevar_cum!0 = "lev"
407rad_thevar_hit&lev = lev
408rad_thevar_cum&lev = lev
409
410rad_thevar_hit = 0
411rad_thevar_cum = 0
412
413; NEED TO CORRECT NLEV
414print("Starting loop")
415do i = psminLocY-loopWidth,psminLocY+loopWidth
416  print("Doing outer loop #: "+(i-psminLocY+loopWidth)+" of "+(loopWidth*2)+"   ("+((tofloat(i-psminLocY+loopWidth)/tofloat(loopWidth*2))*100.)+"%)")
417  do j = psminLocX-loopWidth,psminLocX+loopWidth
418    do k = 0,dimsizes(lev)-1
419      gcdist = tofloat(gc_latlon(psminlat,psminlon,lat(i),lon(j),2,4))
420      if (gcdist .le. outerRad)
421        rad_thevar_hit(k) = rad_thevar_hit(k) + data(k,i,j)*clat(i)
422        rad_thevar_cum(k) = rad_thevar_cum(k) + clat(i)
423      end if
424    end do
425  end do
426end do
427
428rad_thevar = rad_thevar_hit/rad_thevar_cum
429
430copy_VarCoords(rad_thevar_hit,rad_thevar)
431
432rad_thevar@units = data@units
433rad_thevar@long_name=data@long_name
434
435return(rad_thevar)
436
437end
438
439
440
441undef ( "calc_vertmax_from_rad" )
442function calc_vertmax_from_rad ( \
443  vrad[*][*]       : numeric \
444)
445local None
446
447begin
448
449lev = vrad&lev
450nlev = dimsizes(lev)
451rad = vrad&radius
452nrad = dimsizes(rad)
453
454rmw = new((/nlev/),"float")
455
456do k = 0,dimsizes(lev)-1
457  maxix = maxind(vrad(k,:))
458  rmw(k) = tofloat(rad(maxix))
459end do
460
461rmw!0 = "lev"
462rmw&lev = lev
463
464return(rmw)
465
466end
Note: See TracBrowser for help on using the repository browser.