Index: /issm/trunk-jpl/src/m/coordsystems/ll2utm.m
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/ll2utm.m	(revision 27626)
+++ /issm/trunk-jpl/src/m/coordsystems/ll2utm.m	(revision 27627)
@@ -18,6 +18,6 @@
 %	meters)	and F is flattening of the user-defined ellipsoid.
 %
-%	LL2UTM(...,ZONE) forces the UTM ZONE (scalar integer) instead of
-%	automatic set.
+%	LL2UTM(...,ZONE) forces the UTM ZONE (scalar integer or same size as
+%   LAT and LON) instead of automatic set.
 %
 %	[X,Y,ZONE]=LL2UTM(...) returns also the computed UTM ZONE (negative
@@ -37,13 +37,13 @@
 %		   Notes Techniques NT/G 76, janvier 1995.
 %
-%	Acknowledgments: Mathieu.
+%	Acknowledgments: Mathieu, Frederic Christen.
 %
 %
 %	Author: Francois Beauducel, <beauducel@ipgp.fr>
 %	Created: 2003-12-02
-%	Updated: 2014-08-24
-
-
-%	Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%	Updated: 2019-05-29
+
+
+%	Copyright (c) 2001-2019, François Beauducel, covered by BSD License.
 %	All rights reserved.
 %
@@ -93,12 +93,14 @@
 end
 
-if isnumeric(varargin{1}) & size(varargin{1},2) == 2
+if nargin > 1 && isnumeric(varargin{1}) && isnumeric(varargin{2}) ...
+		&& (all(size(varargin{1})==size(varargin{2})) ...
+		|| isscalar(varargin(1)) || isscalar(varargin{2}))
+	lat = varargin{1};
+	lon = varargin{2};
+	v = 2;
+elseif isnumeric(varargin{1}) && size(varargin{1},2) == 2
 	lat = varargin{1}(:,1);
 	lon = varargin{1}(:,2);
 	v = 1;
-elseif nargin > 1 & isnumeric(varargin{2})
-	lat = varargin{1};
-	lon = varargin{2};
-	v = 2;
 else
 	error('Single input argument must be a 2-column matrix [LAT,LON].')
@@ -109,11 +111,23 @@
 end
 
+if any(abs(lat)>90)
+	error('LAT absolute values must be lower than 90.')
+end
+
+% checks for DATUM and/or ZONE syntax
+% NOTE: the following strategy works in any case except if ZONE argument
+% has a size of 1x2 (in that case it will be interpreted as a DATUM). To
+% force the ZONE syntax with 2 elements, just use ZONE(:) to make a colum
+% vector of 2x1.
 for n = (v+1):nargin
 	% LL2UTM(...,DATUM)
-	if ischar(varargin{n}) | (isnumeric(varargin{n}) & numel(varargin{n})==2)
+	if ischar(varargin{n}) || (isnumeric(varargin{n}) ...
+			&& all(size(varargin{n})==[1,2]))
 		datum = varargin{n};
 	% LL2UTM(...,ZONE)
-	elseif isnumeric(varargin{n}) & isscalar(varargin{n})
-			zone = round(varargin{n});
+	elseif isnumeric(varargin{n}) && (isscalar(varargin{n}) ...
+			|| (isscalar(lat) || all(size(varargin{n})==size(lat))) ...
+			&& (isscalar(lon) || all(size(varargin{n})==size(lon))))
+		zone = round(varargin{n});
 	else
 		error('Unknown argument #%d. See documentation.',n)
@@ -142,5 +156,5 @@
 	F0 = round((l1*D0 + 183)/6);
 else
-	F0 = zone;
+	F0 = abs(zone);
 end
 
@@ -148,5 +162,5 @@
 E1 = sqrt((A1*A1 - B1*B1)/(A1*A1));
 P0 = 0/D0;
-L0 = (6*F0 - 183)/D0;		% UTM origin longitude (rad)
+L0 = (6*F0 - 183)/D0;	% UTM origin longitude (rad)
 Y0 = 1e7*(p1 < 0);		% UTM false northern (m)
 N = K0*A1;
@@ -166,9 +180,8 @@
 % same size as x/y in case of crossed zones
 if nargout > 2
-	fu = unique(F0.*sign(lat));
+   	f = F0.*sign(lat);
+	fu = unique(f);
 	if isscalar(fu)
 		f = fu;
-	else
-		f = F0;
 	end
 end
Index: /issm/trunk-jpl/src/m/coordsystems/utm2ll.m
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/utm2ll.m	(revision 27626)
+++ /issm/trunk-jpl/src/m/coordsystems/utm2ll.m	(revision 27627)
@@ -5,5 +5,5 @@
 %	(in degrees). Default datum is WGS84.
 %
-%	X and Y can be scalars, vectors or matrix. Outputs LAT and LON will
+%	X, Y and F can be scalars, vectors or matrix. Outputs LAT and LON will
 %	have the same size as inputs.
 %
@@ -32,8 +32,12 @@
 %	Author: Francois Beauducel, <beauducel@ipgp.fr>
 %	Created: 2001-08-23
-%	Updated: 2014-04-20
-
-
-%	Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%	Updated: 2019-05-29
+
+%	Revision history:
+%
+%	[2019-05-29]
+%	  - fix an issue when X or Y are matrices.
+
+%	Copyright (c) 2001-2019, François Beauducel, covered by BSD License.
 %	All rights reserved.
 %
@@ -74,10 +78,13 @@
 end
 
-if all([numel(x),numel(y)] > 1) && any(size(x) ~= size(y))
-	error('X and Y must be the same size or scalars.')
-end
-
-if ~isnumeric(f) || ~isscalar(f) || f ~= round(f)
-	error('ZONE must be a scalar integer.')
+% checks if input arguments have compatible sizes using unique(complex)
+sz = [size(x);size(y);size(f)];
+sz = complex(sz(:,1),sz(:,2));
+if length(unique(sz(sz~=complex(1,1)))) > 1
+	error('X, Y and ZONE must be scalar or vector/matrix of the same size.')
+end
+
+if ~isnumeric(f) || any(f ~= round(f))
+	error('ZONE must be integer value.')
 end
 
@@ -101,4 +108,9 @@
 end
 
+% calculations are made on column vectors
+x = x(:);
+y = y(:);
+f = f(:);
+
 % constants
 D0 = 180/pi;	% conversion rad to deg
@@ -133,5 +145,5 @@
 p0 = NaN;
 n = 0;
-while any(isnan(p0) | abs(p - p0) > eps) & n < maxiter
+while any(isnan(p0) | abs(p - p0) > eps) && n < maxiter
 	p0 = p;
 	es = E1*sin(p0);
@@ -141,8 +153,10 @@
 
 if nargout < 2
-	lat = D0*[p(:),l(:)];
+	lat = D0*[p,l];
 else
-	lat = p*D0;
-	lon = l*D0;
+	% reshapes vectors to x/y/f original size
+	sz = max([size(x);size(y);size(f)]);
+	lat = reshape(p*D0,sz);
+	lon = reshape(l*D0,sz);
 end