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 | |
---|
20 | ; needs to be a catted file ordered time, lev, lat, lon (or time, lat, lon for PS) |
---|
21 | filename="output_dcmip2016_regular.nc" |
---|
22 | ;======================================================================================= |
---|
23 | |
---|
24 | maxLat=90.0 |
---|
25 | minLat=0.0 |
---|
26 | |
---|
27 | ;======================================================================================= |
---|
28 | |
---|
29 | a = addfile(filename,"r") |
---|
30 | |
---|
31 | ; figure out what is surface level for wind speed... |
---|
32 | if (top_to_bottom) then |
---|
33 | sfcLev=dimsizes(a->lev)-1 ; dimsizes(a->lev)-1 subsets bottom level |
---|
34 | else |
---|
35 | sfcLev=0 |
---|
36 | end if |
---|
37 | |
---|
38 | do varid=1,2 |
---|
39 | |
---|
40 | if( varid .eq. 1) then |
---|
41 | whichVar = "PS" |
---|
42 | PS=a->PS(:,{minLat:maxLat},:) ; just take out NH since we know storm is there |
---|
43 | PS=PS/100. ; convert from Pa to hPa/mb |
---|
44 | PS@long_name="Surface pressure" |
---|
45 | PS@units="hPa" |
---|
46 | var_vtime = dim_min_n_Wrap(PS,(/1,2/)) ; for now, just find min over lat/lon domain at each time |
---|
47 | else if (varid .eq. 2) then |
---|
48 | whichVar = "WIND" |
---|
49 | UBOT=a->U(:,sfcLev,{minLat:maxLat},:) |
---|
50 | VBOT=a->V(:,sfcLev,{minLat:maxLat},:) |
---|
51 | WIND=UBOT |
---|
52 | WIND=sqrt(UBOT^2+VBOT^2) |
---|
53 | WIND@long_name="Lowest level wind" |
---|
54 | WIND@units="m/s" |
---|
55 | var_vtime = dim_max_n_Wrap(WIND,(/1,2/)) |
---|
56 | else |
---|
57 | print("unsupported var choice... exiting...") |
---|
58 | exit |
---|
59 | end if |
---|
60 | end if |
---|
61 | |
---|
62 | wks = gsn_open_wks (out_type,whichVar) ; send graphics to PNG file |
---|
63 | |
---|
64 | res = True ; plot mods desired |
---|
65 | res@gsnDraw=False |
---|
66 | res@gsnFrame=False |
---|
67 | res@tiMainString = "DCMIP2 "+model+" "+whichVar+" v. time" ; add title |
---|
68 | |
---|
69 | res@tiXAxisString = "Time (days)" |
---|
70 | res@tiYAxisString = var_vtime@long_name+" ("+var_vtime@units+")" |
---|
71 | |
---|
72 | res@xyLineThicknessF=5.0 |
---|
73 | res@xyLineColor="Blue" |
---|
74 | |
---|
75 | time = var_vtime&time_counter/86400 |
---|
76 | plot = gsn_csm_xy (wks,time,var_vtime,res) ; create plot |
---|
77 | draw(plot) |
---|
78 | frame(wks) |
---|
79 | |
---|
80 | end do |
---|
81 | |
---|
82 | end |
---|