function [lat,lon] = grid2ll(xgrid,ygrid)
%   grid2ll converts a point defined by the  Universal Transverse Mercator
%  projection on the GLOBEC grid to a latitude and longitude.
%   
%   j. klinck, Feb 2001 (klinck@ccpo.odu.edu)

global SECS COSROT SINROT ONCE
global UTMXOFF UTMYOFF XOFF YOFF ROTANGLE IZONE ISPHERE
global AXIS BXIS AK0 RADSEC

%      set up transformation variables
if isempty(ONCE)
%         need to do only once
  ONCE=1;
  setuputm;
end
[m,n]=size(xgrid);
if n==1
   xgrid=xgrid';
   ygrid=ygrid';
end

x=1000*(xgrid-XOFF);
y=1000*(ygrid-YOFF);

xutm=  COSROT*x - SINROT*y + UTMXOFF;
yutm=  SINROT*x + COSROT*y + UTMYOFF;

[lat,lon,conv]=utm2ll(xutm, yutm);
lat=lat/SECS;
lon=-lon/SECS;
i=find(lon < -180);
lon(i)=lon(i)+360;
conv=conv/SECS;

function [lat,lon,conv]=utm2ll(x,y);
%   convert from x,y in a universal transverse mercator
%   to lat, lon based on the parameters given in SETUPUTM.
%     universal transverse mercator conversion
%
%     rewritten 6-16-83 by j.f. waananen using formulation
%     by j.p. snyder in bulletin 1532, pages 68-69
%
%     convert utm.f to convert from utm to lat lon (idir > 0)
%      2/5/2001  j. klinck

global SECS COSROT SINROT ONCE
global UTMXOFF UTMYOFF XOFF YOFF ROTANGLE IZONE ISPHERE
global AXIS BXIS AK0 RADSEC

a = AXIS(ISPHERE);
b = BXIS(ISPHERE);
es= (a^2-b^2)/(a^2);

%     compute utm zone(izone) and central meridian in seconds for
%     geodetic to utm conversion where zone is not input.

if (IZONE == 0 )
  disp(' *************   error exit from utm. ');
  disp('  zone is not given for conversion from ');
  disp('  x,y to lat,lon. ');
  exit
end

if(abs(IZONE) > 30)
   iutz = abs(IZONE)-30;
   cm=((iutz*6.0d0)-3.0d0)*(-3600.0d0);
else
   iutz = 30-abs(IZONE);
   cm=((iutz*6.0d0)+3.0d0)*3600.0d0;
end

if(IZONE < 0)
  yy=y-1.0d7;
else
  yy = y;
end
xx = x - 5.0d5;
em = yy/AK0;
um = em./(a.*(1.d0-(es/4.d0)-(3.d0*es.*es/6.4d1) - ...
       (5.d0*es.*es.*es/2.56d2)));
e1 = (1.d0-sqrt(1.d0-es))./(1.d0+sqrt(1.d0-es));

phi1 = um+((3.d0*e1/2.d0)-(2.7d1*e1.^3/3.2d1)).*sin(2.d0*um) + ...
        ((2.1d1*e1.*e1/1.6d1)-(5.5d1*e1.^4/3.2d1)).*sin(4.d0*um) + ...
        (1.51d2*e1.^3/9.6d1).*sin(6.d0*um);

en = a./sqrt(1.0d0-es.*sin(phi1).^2);
t  = tan(phi1).^2;
epri = es./(1.d0-es);
c  = epri.*cos(phi1).^2;
r  = (a.*(1.d0-es))./((1.d0-es.*sin(phi1).^2).^1.5d0);
d  = xx./(en*AK0);

phi = phi1 - (en.*tan(phi1)./r) .* ((d.*d/2.d0) - ...
      (5.d0+3.d0*t+10.d0*c-4.d0*c.*c-9.d0*epri).*d.^4/2.4d1 ...
      + (6.1d1+9.d1*t+2.98d2*c+4.5d1*t.*t ...
       -2.52d2*epri-3.d0*c.*c).*d.^6/7.2d2);

alam = (cm/RADSEC)-(d-(1.d0+2.d0*t+c).*d.^3/6.d0 + ...
       (5.d0-2.d0*c+2.8d1*t -3.d0*c.*c+8.d0*epri+2.4d1*t.*t)...
       .*d.^5/1.2d2)./cos(phi1);

lat = phi*RADSEC;
lon = alam*RADSEC;
dlam = -(lon-cm)./RADSEC;

%    test to see if within definition of utm projection
%  ---- latitude within 84 degrees
i=find(abs(lat) > 302400);
if (~isempty(i))
  disp(' *************   error exit from utm. ');
  disp('   latitude value is poleward of 84 degrees');
  disp('   utm is not valid.');
  disp('   calculation has continued but values ');
  disp('   may not be valid.');
    lat(i)/3600.d0
    lon(i)/3600.d0
end
%  ---- delta longitude within 0.16 radians of central meridian
%    print disabled to avoid extra output. Some points on the south west
%    side of the GLOBEC grid are just outside of the permissible area.
%i= (abs(dlam) > 1.6d-1);
%if (~isempty(i))
%  disp(' *************   error exit from utm. ');
%  disp('  d lon not within 0.16 radians of central merid.');
%  disp('   calculation has continued but values ');
%  disp('   may not be valid.');
%    lat(i)/3600.d0
%    lon(i)/3600.d0
%end
%     compute convergence

conv = dlam.*(sin(phi)+1.9587d-12*(dlam.^2).*sin(phi).*cos(phi.*phi));
conv = conv*RADSEC;

%   clean up dummy variables
clear xx yy em um e1 phi1 en t c r d epri phi alam dlam iutz cm

%   end of grid2ll.m
