[437] | 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 |
---|