Add SOILHGT to Rapid Refresh NCEP hybrid Vtable. metgrid complains when SOILHGT is...
[WPS.git] / util / gfs.ncl
blob6ad9c7af6a2fdaa922a089e7d5e6a8e7052a5d32
1 ;***************************************************************************************************************
2 ; This is a sample NCL script to process GRIB data, plot surface variable (SKINTEMP and MSLP) and upper-air 
3 ; variables (U, V, T, RH) in pressure level.  
5 ;**************************************************************************************************************
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
9 begin
10   grb_file = addfile( "gfs_150416_12_00.grb","r")
11   names = getfilevarnames(grb_file)  ; Get the variable names in the
12   print(names)                       ; GRIB file and print them out.
14 ;  atts = getfilevaratts(grb_file,names(0)) ; Get the variable attributes and
15 ;  dims = getfilevardims(grb_file,names(0)) ; dimension names from the GRIB
16   pres= grb_file->lv_ISBL0
17   npres=dimsizes(pres)
18   presRH= grb_file->lv_ISBL6
19   npresRH=dimsizes(presRH)
20   soil= grb_file->lv_DBLL10_l0
21   nsoil=dimsizes(soil)
22   wks = gsn_open_wks("x11","gfs") ; Open an X11 workstation.
23   gsn_define_colormap(wks,"wh-bl-gr-ye-re")
25 ;----------- MCHEN Begin plot -----------------------------------------
26   resources             = True
27   resources@cnFillOn                    = True                 ; turn on color
28   resources@cnLinesOn                   = False                ; no contour lines
29   resources@gsnSpreadColors             = True                 ; use full color map
30   resources@lbLabelAutoStride           = True                    ; every other label
31   resources@cnLevelSelectionMode        = "ManualLevels"       ; manual levels
33   do i = 0,dimsizes(names)-1
34     names_char = stringtochar(names(i))
35 ; PMSL
36     if(names(i).eq."PRMSL_P0_L101_GLL0") then   
37       MSLP = grb_file->PRMSL_P0_L101_GLL0
38         if(isatt(MSLP,"units").eq.True)  
39           resources@tiMainString = "MSLP ( " + MSLP@units  + ")" 
40         else
41           resources@tiMainString = "PRMSL_P0_L101_GLL0"
42         end if               
43     plot = gsn_csm_contour_map_ce(wks,MSLP,resources)
44     end if
45 ; SKINTEMP
46     if(names(i).eq."TMP_P0_L1_GLL0") then   
47       SKINTEMP = grb_file->TMP_P0_L1_GLL0
48         if(isatt(SKINTEMP,"units").eq.True)  
49           resources@tiMainString = "SKINTEMP ( " + SKINTEMP@units + ")"
50         else
51           resources@tiMainString = "TMP_P0_L1_GLL0"
52         end if               
53     plot = gsn_csm_contour_map_ce(wks,SKINTEMP,resources)
54     end if
55 ; T2
56     if(names(i).eq."TMP_P0_L103_GLL0") then   
57       T2 = grb_file->TMP_P0_L103_GLL0
58         if(isatt(T2,"units").eq.True)  
59           resources@tiMainString = "T2    (" + T2@units + ")"
60         else
61           resources@tiMainString = "TMP_P0_L103_GLL0"
62         end if               
63     plot = gsn_csm_contour_map_ce(wks,T2(0,:,:),resources)
64     end if
65 ; U10
66     if(names(i).eq."UGRD_P0_L103_GLL0") then   
67       U10 = grb_file->UGRD_P0_L103_GLL0
68         if(isatt(T2,"units").eq.True)  
69           resources@tiMainString = "U10   (" + U10@units + ")"
70         else
71           resources@tiMainString = "UGRD_P0_L103_GLL0"
72         end if               
73     plot = gsn_csm_contour_map_ce(wks,U10(0,:,:),resources)
74     end if
75 ; V10
76     if(names(i).eq."VGRD_P0_L103_GLL0") then   
77       V10 = grb_file->VGRD_P0_L103_GLL0
78         if(isatt(T2,"units").eq.True)  
79           resources@tiMainString = "V10,UNIT: " + V10@units + ")"
80         else
81           resources@tiMainString = "VGRD_P0_L103_GLL0"
82         end if               
83     plot = gsn_csm_contour_map_ce(wks,V10(0,:,:),resources)
84     end if
85 ;  SOILM in 4 levels
86     if(names(i).eq."SOILW_P0_2L106_GLL0") then   
87       SOILM = grb_file->SOILW_P0_2L106_GLL0
88         if(isatt(SOILM,"units").eq.True)  
89           MainString = "SOILM, UNIT:" + SOILM@units + ")"
90         end if               
91      do nlev=0,nsoil-1 
92           resources@tiMainString = MainString + " At " + soil(nlev)  + " m"
93        plot = gsn_csm_contour_map_ce(wks,SOILM(nlev,:,:),resources)
94      end do
95     end if
96 ;  SOILT in 4 levels
97     if(names(i).eq."TMP_P0_2L106_GLL0") then   
98       SOILT = grb_file->TMP_P0_2L106_GLL0
99         if(isatt(SOILT,"units").eq.True)  
100           MainString = "SOILT  (" + SOILT@units + ")"
101         end if               
102      do nlev=0,nsoil-1 
103           resources@tiMainString = MainString + " At " + soil(nlev) + " m" 
104        plot = gsn_csm_contour_map_ce(wks,SOILT(nlev,:,:),resources)
105      end do
106     end if
107 ;  T in Pres Level
108  ;   if(names_char(0:15).eq."TMP_P0_L100_GLL0") then   
109     if(names(i).eq."TMP_P0_L100_GLL0") then   
110       T = grb_file->TMP_P0_L100_GLL0
111         if(isatt(T,"units").eq.True)  
112           MainString = "TEMPERATURE ( " + T@units + ")"
113         end if               
114      do nlev=0,npres-1 
115           resources@tiMainString = MainString + " At " + pres(nlev)/100 +"hpa"
116        plot = gsn_csm_contour_map_ce(wks,T(nlev,:,:),resources)
117      end do
118     end if
119 ;  RH in Pres Level
120   ;  if(names_char(0:14).eq."RH_P0_L100_GLL0") then   
121     if(names(i).eq."RH_P0_L100_GLL0") then   
122       RH = grb_file->RH_P0_L100_GLL0
123         if(isatt(RH,"units").eq.True)  
124           MainString = "RH  (" + RH@units  + ")"
125         end if               
126      do nlev=0,npresRH-1 
127           resources@tiMainString = MainString + " At " + presRH(nlev)/100 +"hpa"
128        plot = gsn_csm_contour_map_ce(wks,RH(nlev,:,:),resources)
129      end do
130     end if
131 ; UGRD in pressure
132    ; if(names_char(0:16).eq."UGRD_P0_L100_GLL0") then   
133     if(names(i).eq."UGRD_P0_L100_GLL0") then   
134       U = grb_file->UGRD_P0_L100_GLL0
135         if(isatt(U,"units").eq.True)  
136           MainString = "U  ( " + U@units + ")"
137         end if               
138      do nlev=0,npres-1 
139           resources@tiMainString = MainString + " At " + pres(nlev)/100 +"hpa"
140        plot = gsn_csm_contour_map_ce(wks,U(nlev,:,:),resources)
141      end do
142     end if
143 ; VGRD in pressure
144   ;  if(names_char(0:11).eq."VGRD_P0_L100") then   
145     if(names(i).eq."VGRD_P0_L100_GLL0") then   
146       V = grb_file->VGRD_P0_L100_GLL0
147         if(isatt(V,"units").eq.True)  
148           MainString = "V ( " + V@units  + ")"
149         end if               
150      do nlev=0,npres-1 
151           resources@tiMainString = MainString + " At " + pres(nlev)/100 +"hpa"
152        plot = gsn_csm_contour_map_ce(wks,V(nlev,:,:),resources)
153      end do
154     end if
155     delete(names_char)
156   end do