diff -uNr GMT4.0/share/pslib/PSL_prologue.ps GMT4.0_barb/share/pslib/PSL_prologue.ps
--- GMT4.0/share/pslib/PSL_prologue.ps	2004-01-13 04:18:13.000000000 +0900
+++ GMT4.0_barb/share/pslib/PSL_prologue.ps	2005-03-13 08:17:16.000000000 +0900
@@ -39,6 +39,8 @@
 /A4 {0 exch M 0 D D D D D 0 D P S} bind def
 /A5 {0 exch M 0 D D D D D 0 D P V K F U N} bind def
 /A6 {0 exch M 0 D D D D D 0 D P V K F U S} bind def
+/Bt {0 M D D P F} bind def
+/Bl {0 M D S} bind def
 /C0 {0 360 arc V A F U N} bind def
 /C1 {0 360 arc V A F U S} bind def
 /C2 {0 360 arc V C F U N} bind def
diff -uNr GMT4.0/src/GMT.in GMT4.0_barb/src/GMT.in
--- GMT4.0/src/GMT.in	2004-01-03 12:30:02.000000000 +0900
+++ GMT4.0_barb/src/GMT.in	2005-03-13 08:17:16.000000000 +0900
@@ -73,6 +73,7 @@
 	grd2cpt		Make a color palette table from a grdfile
 	grdblend	Blend several partially over-lapping grdfiles onto one grid
 	grd2xyz		Conversion from 2-D gridded file to table data
+	grdbarb	        Plot wind barbs from grdfiles
 	grdclip		Limit the z-range in gridded data
 	grdcontour	Contouring of 2-D gridded data
 	grdcut		Cut a sub-region from a grd file
diff -uNr GMT4.0/src/Makefile GMT4.0_barb/src/Makefile
--- GMT4.0/src/Makefile	2004-08-25 09:19:15.000000000 +0900
+++ GMT4.0_barb/src/Makefile	2005-03-13 08:17:16.000000000 +0900
@@ -108,7 +108,7 @@
 	  psbasemap.c psclip.c pscoast.c pshistogram.c psimage.c pslegend.c psmask.c \
 	  psrose.c psscale.c pstext.c pscontour.c pswiggle.c psxy.c psxyz.c \
 	  sample1d.c spectrum1d.c splitxyz.c surface.c trend1d.c trend2d.c \
-	  triangulate.c xyz2grd.c
+	  triangulate.c xyz2grd.c grdbarb.c
 
 first:	all
 
@@ -149,7 +149,7 @@
 all:		init libs $(PROGS)
 
 init:		gmtmacros gmt_notposix.h gmt_nan.h
-		
+
 gmtmacros:
 		if [ `cat makegmt.macros | wc -c` = 0 ]; then \
 			echo "src/makegmt.macros is empty - you must rerun configure in the main GMT directory"; \
diff -uNr GMT4.0/src/grdbarb.c GMT4.0_barb/src/grdbarb.c
--- GMT4.0/src/grdbarb.c	1970-01-01 09:00:00.000000000 +0900
+++ GMT4.0_barb/src/grdbarb.c	2005-03-13 08:16:50.000000000 +0900
@@ -0,0 +1,409 @@
+/*--------------------------------------------------------------------
+ *	$Id: grdvector.c,v 1.14 2004/08/17 01:48:06 pwessel Exp $
+ *
+ *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
+ *	See COPYING file for copying and redistribution conditions.
+ *
+ *	This program is free software; you can redistribute it and/or modify
+ *	it under the terms of the GNU General Public License as published by
+ *	the Free Software Foundation; version 2 of the License.
+ *
+ *	This program is distributed in the hope that it will be useful,
+ *	but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *	GNU General Public License for more details.
+ *
+ *	Contact info: gmt.soest.hawaii.edu
+ *--------------------------------------------------------------------*/
+/*
+   grdvector reads 2 grdfiles that contains the 2 components of a vector
+   field (cartesian or polar) and plots vectors at the grid positions.
+   This is basically a short-hand for using grd2xyz | psxy -SV and is
+   more convenient for such plots on a grid.
+
+   Author:	Paul Wessel
+   Date:	12-JUN-1995
+   Revised:	15-FEB-2000
+   Version:	4
+
+   grdbarb ( based source code is grdvector.c ) reads 2 grd files that
+   contains the 2 components of a vector field ( cartesian or polar )
+   and plots wind barbs at the grid positions. This is basically a
+   short-hand for using grd2xyz | psxy -SU and is more convenient for
+   such plots on a grid.
+
+   Author:      KATO Masay <kato@ep.sci.hokudai.ac.jp>
+   Date:        13-MAR-2005
+   Version:     0.6 ( for GMT4.0 )
+ 
+ */
+
+
+#include "gmt.h"
+
+float *r, *theta;
+
+main (int argc, char **argv)
+{
+
+	int	i, j, n = 0, nm, nx, ny, ij, i0, j0, di, dj, off;
+	int     dummy[4], knot_or_ms = 0;
+	
+	BOOLEAN convert_angles = FALSE, get_rgb = FALSE, cartesian = TRUE, shrink = FALSE, set_fill = FALSE;
+	BOOLEAN error = FALSE, center = FALSE, outline = FALSE, azimuth = FALSE, inc_set = FALSE;
+	BOOLEAN clip = TRUE, got_fix_length = FALSE;
+	
+	char *file[2], *cpt, txt_a[32], txt_b[32], txt_c[32], unit = 0;
+	
+	double dx2, dy2, b_width = 0.1, b_length = 0.2, b_angle = 120., b_scale = 5.;
+	double v_w, h_l, h_w, v_shrink, v_norm = 0.0, tmp, x, y, plot_x, plot_y, x_off, y_off;
+	double west, east, south, north, x2, y2, scale = 1.0, fix_length = 0.0;
+	double data_west, data_east, data_south, data_north, value, c, s;
+	
+	struct GRD_HEADER h[2];
+	struct GMT_FILL fill;
+	struct GMT_PEN pen;
+
+	GMT_init_pen (&pen, GMT_PENWIDTH);
+	GMT_init_fill (&fill, -1, -1, -1);
+	west = east = south = north = 0.0;
+	dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
+	di = dj = 1;
+	i0 = j0 = 0;
+	dx2 = dy2 = 0.0;
+	
+	argc = GMT_begin (argc, argv);
+
+	for (i = 1; i < argc; i++) {
+		if (argv[i][0] == '-') {
+			switch (argv[i][1]) {
+				/* Common parameters */
+			
+				case 'B':
+				case 'J':
+				case 'K':
+				case 'O':
+				case 'P':
+				case 'R':
+				case 'U':
+				case 'V':
+				case 'X':
+				case 'x':
+				case 'Y':
+				case 'y':
+				case 'c':
+				case '\0':
+					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
+					break;
+				
+				/* Supplemental parameters */
+			
+				case 'A':
+					cartesian = FALSE;
+					break;
+				case 'E':
+					center = TRUE;
+					break;
+				case 'I':	/* Only use gridnodes dx2,dy2 apart */
+					GMT_getinc (&argv[i][2], &dx2, &dy2);
+					inc_set = TRUE;
+					break;
+				case 'N':	/* Do not clip at border */
+					clip = FALSE;
+					break;
+				case 'Q':
+					if (argv[i][2] && argv[i][3] != 'n') {
+						if (sscanf (&argv[i][2], "%[^/]/%[^/]/%lf/%lf", txt_a, txt_b, &b_angle, &b_scale) != 4) {
+						        fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Could not decode barbwidth/barblength/barbangle/barbscale\n", GMT_program);
+							error++;
+						}
+						else {
+							b_width  = GMT_convert_units (txt_a, GMT_INCH);
+							b_length = GMT_convert_units (txt_b, GMT_INCH);
+						}
+					}
+					break;
+				case 'T':
+					convert_angles = TRUE;
+					break;
+				case 'W':		/* Set line attributes */
+					if (argv[i][2] && GMT_getpen (&argv[i][2], &pen)) {
+						GMT_pen_syntax ('W');
+						error++;
+					}
+					outline = TRUE;
+					break;
+				case 'Z':
+					azimuth = TRUE;
+					break;
+				case 'C':
+				case 'G':
+				case 'S':
+					break;
+				default:
+					error = TRUE;
+					GMT_default_error (argv[i][1]);
+					break;
+			}
+		}
+		else if (n < 2)
+			file[n++] = argv[i];
+		else
+			n++;
+	}
+	
+	if (argc == 1 || GMT_quick) {
+		fprintf (stderr, "grdbarb %s - Plot wind barb fields from grdfiles\n\n", GMT_VERSION);
+		fprintf (stderr, "usage: grdbarb compx.grd compy.grd -J<params> -R<west/east/south/north> [-A]\n");
+		fprintf (stderr, "\t[-B<tickinfo>] [-E] [-I<dx/dy>] [-K] [-O] [-P] [-Q<params>] [-N] [-T]\n");
+		fprintf (stderr, "\t[-U[<label>]] [-V] [-W<pen>] [-X<x_shift>] [-Y<y_shift>] [-Z] [-c<ncopies>] [-:]\n\n");
+		
+		if (GMT_quick) exit (EXIT_FAILURE);
+		
+		fprintf (stderr, "\tcompx & compy are grdfiles with the 2 vector components.\n");
+		GMT_explain_option ('j');
+		fprintf (stderr, "\n\tOPTIONS:\n");
+		fprintf (stderr, "\t-A means grdfiles have polar (r, theta) components [Default is Cartesian]\n");
+		GMT_explain_option ('b');
+		fprintf (stderr, "\t-E cEnter vectors on grid nodes [Default draws from grid node]\n");
+		fprintf (stderr, "\t-I plots only those nodes that are <dx/dy> apart [Default is all nodes]\n");
+		GMT_explain_option ('K');
+		fprintf (stderr, "\t-N Do Not clip vectors that exceed the map boundaries [Default will clip]\n");
+		GMT_explain_option ('O');
+		GMT_explain_option ('P');
+		fprintf (stderr, "\t-Q Select vector plot.\n");
+		fprintf (stderr, "\t   Optionally, specify vector parameters\n");
+		fprintf (stderr, "\t   <params> are barbwidth/barblength/barbangle/barbscale [Default is 0.4i/1.0i/120/5]\n");
+		fprintf (stderr, "\t   Append n<size> which will cause vectors shorter than <size> to be\n");
+		fprintf (stderr, "\t     scaled down\n");
+		GMT_explain_option ('R');
+		fprintf (stderr, "\t-T means azimuth should be converted to angles based on map projection\n");
+		GMT_explain_option ('U');
+		GMT_explain_option ('V');
+		fprintf (stderr, "\t-W sets pen attributes [width = %gp, color = (%d/%d/%d), texture = solid line].\n", 
+			pen.width, pen.rgb[0], pen.rgb[1], pen.rgb[2]);
+		GMT_explain_option ('X');
+		fprintf (stderr, "\t-Z means the angles provided are azimuths rather than direction\n");
+		GMT_explain_option ('c');
+		GMT_explain_option ('.');
+		exit (EXIT_FAILURE);
+	}
+
+	if (inc_set && (dx2 <= 0.0 || dy2 <= 0.0)) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify positive increments\n", GMT_program);
+		error++;
+	}
+	/*
+	if (scale == 0.0 && !got_fix_length) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Scale must be nonzero\n", GMT_program);
+		error++;
+	}
+	if (fix_length <= 0.0 && got_fix_length) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -Sl option:  Length must be positive\n", GMT_program);
+		error++;
+	}
+	if (got_fix_length && shrink) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -Sl, -Q options:  Cannot use -Q..n<size> with -Sl\n", GMT_program);
+		error++;
+	}
+	*/
+	if (azimuth && cartesian) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  Azimuths not valid input for Cartesian data\n", GMT_program);
+		error++;
+	}
+	if (get_rgb && !cpt) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify a color palette table\n", GMT_program);
+		error++;
+	}
+	if (!(set_fill || outline || get_rgb)) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -W\n", GMT_program);
+		error++;
+	}
+	if (n != 2) {
+		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify two input grdfiles\n", GMT_program);
+		error++;
+	}
+	
+	if (error) exit (EXIT_FAILURE);
+
+	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */
+
+	if (get_rgb) GMT_read_cpt (cpt);
+
+	if (!(strcmp (file[0], "=") || strcmp (file[1], "="))) {
+		fprintf (stderr, "%s: Piping of grdfiles not supported!\n", GMT_program);
+		exit (EXIT_FAILURE);
+	}
+
+	if (GMT_read_grd_info (file[0], &h[0])) {
+		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file[0]);
+		exit (EXIT_FAILURE);
+	}
+	
+	if (GMT_read_grd_info (file[1], &h[1])) {
+		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file[1]);
+		exit (EXIT_FAILURE);
+	}
+	
+	if (!(h[0].nx == h[1].nx && h[0].ny == h[1].ny && h[0].x_min == h[1].x_min && h[0].y_min == h[1].y_min 
+		&& h[0].x_inc == h[1].x_inc && h[0].y_inc == h[1].y_inc)) {
+		fprintf (stderr, "%s: files %s and %s does not match!\n", GMT_program, file[0], file[1]);
+		exit (EXIT_FAILURE);
+	}
+	off = (h[0].node_offset) ? 0 : 1;
+		
+	/* Determine what wesn to pass to map_setup */
+
+	if (!project_info.region_supplied) {
+		west = h[0].x_min;
+		east = h[0].x_max;
+		south = h[0].y_min;
+		north = h[0].y_max;
+	}
+
+	GMT_map_setup (west, east, south, north);
+
+	/* Determine the wesn to be used to read the grdfile */
+
+	if (GMT_grd_setregion (&h[0], &data_west, &data_east, &data_south, &data_north)) {	/* No grid to plot; just do empty map and exit */
+		ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
+			gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
+			gmtdefs.dpi, GMT_INCH , gmtdefs.paper_width, gmtdefs.page_rgb, gmtdefs.encoding.name, GMT_epsinfo (argv[0]));
+		GMT_echo_command (argc, argv);
+		if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
+		GMT_map_basemap ();
+		ps_plotend (gmtdefs.last_page);
+		GMT_end (argc, argv);
+	}
+
+	/* Read data */
+
+	nx = irint ( (data_east - data_west) / h[0].x_inc) + off;
+	ny = irint ( (data_north - data_south) / h[0].y_inc) + off;
+	nm = nx * ny;
+	r = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
+	theta = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
+	if (GMT_read_grd (file[0], &h[0], r, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
+		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[0]);
+		exit (EXIT_FAILURE);
+	}
+	if (GMT_read_grd (file[1], &h[1], theta, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
+		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, file[1]);
+		exit (EXIT_FAILURE);
+	}
+	/* if (!got_fix_length) scale = 1.0 / scale; */
+
+	switch (unit) {	/* Adjust for possible unit selection */
+		case 'C':
+		case 'c':
+			scale *= GMT_u2u[GMT_CM][GMT_INCH];
+			fix_length *= GMT_u2u[GMT_CM][GMT_INCH];
+			break;
+		case 'I':
+		case 'i':
+			scale *= GMT_u2u[GMT_INCH][GMT_INCH];
+			fix_length *= GMT_u2u[GMT_INCH][GMT_INCH];
+			break;
+		case 'M':
+		case 'm':
+			scale *= GMT_u2u[GMT_M][GMT_INCH];
+			fix_length *= GMT_u2u[GMT_M][GMT_INCH];
+			break;
+		case 'P':
+		case 'p':
+			scale *= GMT_u2u[GMT_PT][GMT_INCH];
+			fix_length *= GMT_u2u[GMT_PT][GMT_INCH];
+			break;
+		default:
+			scale *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH];
+			fix_length *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH];
+			break;
+	}
+	
+	ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
+		gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
+		gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, gmtdefs.encoding.name, GMT_epsinfo (argv[0]));
+	GMT_echo_command (argc, argv);
+	if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
+
+	GMT_setpen (&pen);
+        if (clip) GMT_map_clip_on (GMT_no_rgb, 3);
+
+	if (dx2 != 0.0 && dy2 != 0.0) {
+		dj = irint (dy2 / h[0].y_inc);
+		di = irint (dx2 / h[0].x_inc);
+		tmp = ceil (h[0].y_max / dy2) * dy2;
+		if (tmp > h[0].y_max) tmp -= dy2;
+		j0 = irint ((h[0].y_max - tmp) / h[0].y_inc);
+		tmp = floor (h[0].x_min / dx2) * dx2;
+		if (tmp < h[0].x_min) tmp += dx2;
+		i0 = irint ((tmp - h[0].x_min) / h[0].x_inc);
+	}
+	
+	for (j = j0; j < h[1].ny; j += dj) {
+		y = h[0].y_max - j *h[0].y_inc;
+		for (i = i0; i < h[1].nx; i += di) {
+		
+			ij = j * h[0].nx + i;
+			if (GMT_is_fnan (r[ij]) || GMT_is_fnan (theta[ij])) continue;	/* Cannot plot NaN-vectors */
+			
+			value = r[ij];
+			
+			if (cartesian) {
+				value = hypot (theta[ij], r[ij]);
+				if (value == 0.0) continue;
+				theta[ij] = (float)(R2D * atan2 (theta[ij], r[ij]));
+				r[ij] = (float)value;
+			}
+			else if (r[ij] < 0.0) {
+				r[ij] = -r[ij];
+				theta[ij] += 180.0;
+			}
+			else if (r[ij] == 0.0) continue;
+			
+			if (get_rgb) GMT_get_rgb24 (value, fill.rgb);
+			
+			x = h[0].x_min + i * h[0].x_inc;
+			GMT_geo_to_xy (x, y, &plot_x, &plot_y);
+			
+			if (convert_angles) {
+				if (!azimuth) theta[ij] = (float)90.0 - theta[ij];
+				GMT_azim_to_angle (x, y, 0.1, (double)theta[ij], &tmp);
+				theta[ij] = (float)(tmp * D2R);
+			}
+			else
+				theta[ij] *= (float)D2R;
+			
+			if (got_fix_length)
+				r[ij] = (float)fix_length;
+			else
+				r[ij] *= (float)scale;
+			
+			sincos (theta[ij], &s, &c);
+			x2 = plot_x + r[ij] * c;
+			y2 = plot_y + r[ij] * s;
+			
+			if (center) {
+				x_off = 0.5 * (x2 - plot_x);
+				y_off = 0.5 * (y2 - plot_y);
+				plot_x -= x_off;
+				plot_y -= y_off;
+				x2 -= x_off;
+				y2 -= y_off;
+			}
+
+			ps_barb (plot_x, plot_y, x2, y2, b_width, b_length, b_angle, b_scale, knot_or_ms, fill.rgb, outline);
+		}
+	}
+	
+        if (clip) GMT_map_clip_off ();
+
+	GMT_map_basemap ();
+
+	ps_plotend (gmtdefs.last_page);
+
+	GMT_free ((void *)r);
+	GMT_free ((void *)theta);
+	
+	GMT_end (argc, argv);
+}
diff -uNr GMT4.0/src/pslib.c GMT4.0_barb/src/pslib.c
--- GMT4.0/src/pslib.c	2004-09-04 05:26:45.000000000 +0900
+++ GMT4.0_barb/src/pslib.c	2005-03-13 08:17:16.000000000 +0900
@@ -109,6 +109,7 @@
  *	ps_transrotate		: Translates and rotates the coordinate system
  *	ps_triangle		: Plots a triangle and [optionally] fills it
  *	ps_vector		: Draws an vector as specified
+ *      ps_barb                 : Draws a barb as specified
  *	ps_words		: Plots a text paragraph
  *
  *
@@ -124,6 +125,10 @@
  *
  * The environmental variable GMTHOME must be set to the directory that holds the subdirectory
  *   share/pattern where all the pattern Sun raster files are stored
+ *
+ * Date:        13-MAR-2005
+ * Modified:    KATO Masaya <kato@ep.sci.hokudai.ac.jp>
+ *
  */
 
 #include "pslib_inc.h"
@@ -2419,6 +2424,76 @@
 	 ps_vector (*xtail, *ytail, *xtip, *ytip, *tailwidth, *headlength, *headwidth, *headshape, rgb, *outline);
 }
 
+void ps_barb (double xtail, double ytail, double xtip, double ytip, double barbwidth, double barblength, double barbangle, double barbscale, int knot_or_ms, int rgb[], int outline)
+{
+	/* Will make sure that arrow has a finite width in PS coordinates */
+	
+	double angle, barb_tmp, length;
+	int w0, w1, dx, dy, nloop, number_of_barb, loop, itri;
+	int num_of_tri, num_of_ful, num_of_hlf;
+	int num_of_TRI, num_of_FUL, num_of_BLB;
+	
+	length = hypot((xtail-xtip), (ytail-ytip));
+	if (length == 0.) return;
+
+	angle = atan2 ((ytip-ytail),(xtip-xtail)) * R2D + 180.;
+	fprintf (ps.fp, "V %d %d T", irint (xtail * ps.scale), irint (ytail * ps.scale));
+	if (angle != 0.0) fprintf (ps.fp, " %lg R", angle);
+
+	w0 = irint (barblength * ps.scale / 2.);
+	w1 = irint (0.25 * (barblength * ps.scale - w0));
+	dx = irint(- barbwidth * ps.scale * cos(barbangle / R2D));
+	dy = irint(- barbwidth * ps.scale * sin(barbangle / R2D));
+
+	nloop = 1;
+	num_of_tri = (int) (length/barbscale/5. + 0.0001);
+	barb_tmp   = length - barbscale * (double) ( 5 * num_of_tri );
+	num_of_ful = (int)(barb_tmp/barbscale + 0.0001);
+	barb_tmp   = barb_tmp - barbscale * (double)num_of_ful;
+	num_of_hlf = (int)(barb_tmp * 2. / barbscale + 0.0001);
+	number_of_barb = num_of_tri + num_of_ful + num_of_hlf;
+
+	num_of_TRI = 0;
+	num_of_FUL = num_of_ful;
+	num_of_BLB = 0;
+	itri = (num_of_tri == 0) ? 0 : 1;
+	itri *= w1;
+
+	fprintf (ps.fp, " %d 0 0 Bl", irint(barblength * ps.scale));
+
+	for(loop=1;loop <= number_of_barb; loop++) {
+	  if(num_of_tri != 0){
+	    if(num_of_FUL == 0){
+	      fprintf(ps.fp, " %d %d %d %d %d Bt",dx,dy,-2*w1,0,irint(barblength*ps.scale)+w1+2*num_of_TRI*w1);
+	    }else{
+	      fprintf(ps.fp, " %d %d %d %d %d Bt",dx,dy,-2*w1,0,irint(barblength*ps.scale)+2*w1+2*num_of_TRI*w1);
+	    }
+	    num_of_tri--;
+	    num_of_TRI++;
+	  }else if(num_of_ful != 0){
+	    fprintf(ps.fp, " %d %d %d Bl",dx,dy,irint(barblength*ps.scale)-w1*num_of_BLB-itri);
+	    num_of_ful--;
+	    num_of_BLB++;
+	  }else{
+	    if(num_of_BLB == 0){
+	      fprintf(ps.fp, " %d %d %d Bl",irint(dx/2),irint(dy/2),irint(barblength*ps.scale)-w1-itri);
+	    }else{
+	      fprintf(ps.fp, " %d %d %d Bl",irint(dx/2),irint(dy/2),irint(barblength*ps.scale)-w1*num_of_BLB-itri);
+	      num_of_hlf = 0;
+	    }
+	  }
+	}	  
+	fprintf(ps.fp, " U \n");
+}
+
+/* fortran interface */
+/*
+void ps_barb_ (double *xtail, double *ytail, double *xtip, double *ytip, double *barbwidth, double *barblength, double *barbwidth, double *barbscale, int *knot_or_ms, int *rgb, int *outline)
+{
+	 ps_barb (*xtail, *ytail, *xtip, *ytip, *barbwidth, *barblength, *barbwidth, *barbscale, *knot_or_ms, rgb, *outline);
+}
+*/
+
 /* Support functions used in ps_* functions.  No Fortran bindings needed */
 
 void get_uppercase (char *new, char *old)
diff -uNr GMT4.0/src/pslib.h GMT4.0_barb/src/pslib.h
--- GMT4.0/src/pslib.h	2004-06-04 13:24:27.000000000 +0900
+++ GMT4.0_barb/src/pslib.h	2005-03-13 08:17:16.000000000 +0900
@@ -124,6 +124,7 @@
 EXTERN_MSC void ps_transrotate(double x, double y, double angle);
 EXTERN_MSC void ps_triangle(double x, double y, double side, int rgb[], int outline);
 EXTERN_MSC void ps_vector(double xtail, double ytail, double xtip, double ytip, double tailwidth, double headlength, double headwidth, double headshape, int rgb[], int outline);
+EXTERN_MSC void ps_barb(double xtail, double ytail, double xtip, double ytip, double barbwidth, double barblength, double barbangle, double barbscale, int knot_or_ms, int rgb[], int outline);
 EXTERN_MSC unsigned char *ps_loadraster (char *file, struct rasterfile *header, BOOLEAN invert, BOOLEAN monochrome, BOOLEAN colorize, int f_rgb[], int b_rgb[]);
 EXTERN_MSC void ps_words (double x, double y, char **text, int n_words, double line_space, double par_width, int par_just, int font, double font_size, double angle, int rgb[3], int justify, int draw_box, double x_off, double y_off, double x_gap, double y_gap, int boxpen_width, char *boxpen_texture, int boxpen_offset, int boxpen_rgb[], int vecpen_width, char *vecpen_texture, int vecpen_offset, int vecpen_rgb[], int boxfill_rgb[3]);
 EXTERN_MSC void ps_setline(int linewidth);
diff -uNr GMT4.0/src/psxy.c GMT4.0_barb/src/psxy.c
--- GMT4.0/src/psxy.c	2004-09-03 03:15:23.000000000 +0900
+++ GMT4.0_barb/src/psxy.c	2005-03-13 08:17:16.000000000 +0900
@@ -26,6 +26,9 @@
  * Date:	15-SEP-2001
  * Version:	4
  *
+ * Date:        13-MAR-2005
+ * Modified:    KATO Masaya <kato@ep.sci.hokudai.ac.jp>
+ *
  */
 
 #include "gmt.h"
@@ -54,6 +57,7 @@
 #define PENTAGON	19
 #define OCTAGON		20
 #define CUSTOM		21
+#define BARB            22
 #define FRONT		-100
 #define QUOTED_LINE	-200
 
@@ -99,12 +103,17 @@
 
 	int convert_angles;	/* If 2, convert azimuth to angle on map, 1 special case for -JX, 0 plain case */
 	BOOLEAN read_vector;	/* if TRUE must read vector attributes */
+        BOOLEAN read_barb;      /* if TRUE must read wind barb attributes */
 	BOOLEAN shrink;		/* If TRUE, shrink vector attributes for small lengths */
 	double v_norm;		/* shrink when lengths are smaller than this */
 	double v_shrink;	/* Required scale factor */
 	double v_width;		/* Width of vector stem in inches */
 	double h_length;	/* Length of vector head in inches */
 	double h_width;		/* Width of vector head in inches */
+	double b_width;		/* Width of wind barb in inches */
+	double b_length;	/* Length of wind barb in inches */
+	double b_angle;		/* Angle of wind barb in degree */
+	double b_scale;		/* Scale of each wind barbs */
 	int v_just;		/* How to justify vector: head point given (3), head (2), center(1), tail (0 - Default) */
 	BOOLEAN v_double_heads;		/* If TRUE, Add 8 (|= 8) to outline to specify double-headed vector (FALSE is single-headed) */
 
@@ -167,15 +176,23 @@
 	S.symbol = S.n_required = S.convert_angles = 0;
 	S.size = S.base = S.v_norm = S.v_shrink = 0.0;
 	S.read_size = S.user_unit  = S.shrink = S.read_vector = FALSE;
+	S.read_barb = FALSE;
 	S.v_width = 0.03;
 	memset ((void *)S.string, 0, (size_t)(16 * sizeof (char)));
 	S.h_length = 0.12;
 	S.h_width = 0.1;
 	S.font_no = gmtdefs.annot_font[0];
 
+	S.b_width  = 0.1;
+	S.b_length = 0.2;
+	S.b_angle  = 120.;
+	S.b_scale  = 5.;
+
 	if ((S.u = gmtdefs.measure_unit) == GMT_CM) {
 		S.v_width = 0.075 / 2.54;	S.h_length = 0.3 / 2.54;
 		S.h_width = 0.25 / 2.54; error_width = 0.25 / 2.54;
+		S.b_width = 0.25 / 2.54;
+		S.b_length = 0.50 / 2.54;
 	}
 	
 	/* Check and interpret the command line arguments */
@@ -346,7 +363,7 @@
 		GMT_explain_option ('P');
 		fprintf (stderr, "\t-S to select symbol type and symbol size (in %s).  Choose between\n", GMT_unit_names[gmtdefs.measure_unit]);
 		fprintf (stderr, "\t   -(xdash), st(a)r, (b)ar, (c)ircle, (d)iamond, (e)llipse, (f)ront, octa(g)on, (h)exagon, (i)nvtriangle,\n");
-		fprintf (stderr, "\t   (k)ustom, (l)etter, pe(n)tagon, (p)oint, (q)uoted line, (r)ect, (s)quare, (t)riangle, (v)ector, (w)edge, (x)cross, (y)dash\n");
+		fprintf (stderr, "\t   (k)ustom, (l)etter, pe(n)tagon, (p)oint, (q)uoted line, (r)ect, (s)quare, (t)riangle, (u)barb, (v)ector, (w)edge, (x)cross, (y)dash\n");
 		fprintf (stderr, "\t   If no size is specified, psxy expects the 3rd column to have sizes.\n");
 		fprintf (stderr, "\t   If no symbol is specified, psxy expects the last column to have symbol code.\n");
 		fprintf (stderr, "\t   [Note: if -C is selected then 3rd means 4th column, etc.]\n");
@@ -380,6 +397,10 @@
 		fprintf (stderr, "\t     Furthermore, <size> means arrowwidth/headlength/headwith [Default is %g/%g/%g]\n", S.v_width, S.h_length, S.h_width);
 		fprintf (stderr, "\t     If -SV rather than -Sv is selected, psxy will expect azimuth and length\n");
 		fprintf (stderr, "\t     and convert azimuths based on the chosen map projection\n");
+		fprintf (stderr, "\t   Barbs: the direction and length must be in input columns 3 and 4.\n");
+		fprintf (stderr, "\t     Furthermore, <size> means barbwidth/barblength/barbangle/barbscale [Default is %lg/%lg/%lg/%lg]\n", S.b_width, S.b_length, S.b_angle, S.b_scale);
+		fprintf (stderr, "\t     If -SU rather than -Su is selected, psxy will expect azimuth and length\n");
+		fprintf (stderr, "\t     and convert azimuths based on the chosen map projection\n");
 		fprintf (stderr, "\t     Insert h(head), b(balance point), or t(ail) after -Sv|V to justify vector w.r.t. input (x,y).\n");
 		fprintf (stderr, "\t     Insert s(egment) if (x,y) is tail and columns 3 and 4 holds head (x,y).\n");
 		fprintf (stderr, "\t     Upper case H, B, T, S will give double-headed vector [Default is t]\n");
@@ -407,7 +428,7 @@
 	
 	if (S.symbol > 0 && cpt_given) get_rgb = TRUE;	/* Need to read z-vales from input data file */
 	
-	if (error_bars && (S.read_vector || S.symbol == ELLIPSE || S.symbol == FRONT || S.symbol == QUOTED_LINE)) {
+	if (error_bars && (S.read_vector || S.read_barb || S.symbol == ELLIPSE || S.symbol == FRONT || S.symbol == QUOTED_LINE)) {
 		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Incompatible with -Se, -Sf, -Sq, -Sv\n", GMT_program);
 		error++;
 	}
@@ -652,6 +673,10 @@
 					x_len = 0.5 * in[two+S.read_size];
 					y_len = 0.5 * in[three+S.read_size];
 				}
+				else if (S.read_barb){
+					direction = in[two];
+					length = in[three];
+				}
 				
 				/* Skip zero-size symbols */
 				
@@ -664,6 +689,10 @@
 					if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Vector length <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read);
 					continue;
 				}
+				if (S.read_barb && length <= 0.0) {
+					if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Vector length <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read);
+					continue;
+				}
 				if (S.symbol == ELLIPSE && (major <= 0.0 || minor <= 0.0)) {
 					if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Ellipse axes <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read);
 					continue;
@@ -838,6 +867,16 @@
 					case RECT:
 						ps_rect (plot_x - x_len, plot_y - y_len, plot_x + x_len, plot_y + y_len, fill.rgb, outline);
 						break;
+					case BARB:
+						if (S.convert_angles) {
+							GMT_azim_to_angle (in[0], in[1], 0.1, direction, &tmp);
+							direction = tmp;
+						}
+						sincos (direction * D2R, &s, &c);
+						x2 = plot_x + length * c;
+						y2 = plot_y + length * s;
+						ps_barb (plot_x, plot_y, x2, y2, S.b_width, S.b_length, S.b_angle, S.b_scale, 0, fill.rgb, outline);
+						break;
 					case CUSTOM:
 						GMT_draw_custom_symbol (plot_x, plot_y, 0.0, S.size, S.custom, &pen, &fill, outline);
 						break;
@@ -1272,6 +1311,24 @@
 		case 't':
 			p->symbol = TRIANGLE;
 			break;
+		case 'U':
+			p->convert_angles = TRUE;
+		case 'u':
+			p->symbol = BARB;
+			len = strlen(text) - 1;
+			if (text[1]) {
+				if (isalpha ((int)text[len]) && isalpha ((int)text[len-1])) {
+					p->u = GMT_get_unit (text[len]);
+					text[len] = 0;
+				}
+				sscanf (&text[1], "%[^/]/%[^/]/%lf/%lf", txt_a, txt_b, &p->b_angle, &p->b_scale);
+				p->b_width  = GMT_convert_units (txt_a, GMT_INCH);
+				p->b_length = GMT_convert_units (txt_b, GMT_INCH);
+			}
+			p->read_barb = TRUE;
+			p->n_required = 2;
+			check = FALSE;
+			break;
 		case 'V':
 			p->convert_angles = 1;
 		case 'v':