gdal_transform can be used to transform the grid X,Y indices into meters using on the projection info from A Guide to NSIDC's Polar Stereographic Projection | National Snow and Ice Data Center
the end result is similar to the example 28, except the linear scale annotation on the map axes make no sense when using stereographic projection.
#!/bin/sh
rm -f gmt.*
gdal_translate -of NETCDF -a_ullr -3850000 -5350000 3750000 5850000 -a_srs EPSG:3413 NETCDF:"RDEFT4_20101021.nc":sea_ice_thickness grid_translated.nc
# set ice thickness below 0 to NaN
gmt grdclip grid_translated.nc -Sb0/NaN -Ggrid_translated.nc
# define map scale
SCALE=1:75000000
gmt begin RDEFT4_20101021_sea_ice_thickness1 png
# plot Cartesian grid data, in meters, scale as defined above
gmt grdimage grid_translated.nc -Jx$SCALE
# coastline, stereographic projection, origin at -45 degrees lon/ 90 degrees lat.
# Scale (defined above) valid at a standard latitude 70 degrees (no distortion at 70 degreed north)
# -Js-45/90/70/1:75000000
gmt coast -Rgrid_translated.nc -Js-45/90/70/$SCALE -Bxa90-45g15 -Byg15 -Di -W0.25p #-Gwhite
gmt colorbar -DJRM+w5c/0.5c+o1.5c/0+mc -F+p+i -Bxa -By+lm
gmt end
