function [G]=check_metrics(gname);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Copyright (c) 2003 ROMS/TOMS Adjoint Group                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [grid]=check_metrics(gname)                                      %
%                                                                           %
% This routine checks grid metrics.                                         %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    gname       GRID NetCDF file name (character string).                  %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    G           Grid structure array with several quantities.              %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Read in grid quatities.

h=nc_read(gname,'h');
[Lp,Mp]=size(h);
L=Lp-1; Lm=L-1;
M=Mp-1; Mm=M-1;

pm=nc_read(gname,'pm');
pn=nc_read(gname,'pn');

dx=1./pm;
dy=1./pn;

scale=0.000001;


% Derive quantities at u-, v-, psi-points using:  pm, pn 

G.dxr=dx;
G.dyr=dy;
G.areaR=scale*sum(sum(G.dxr(2:L,2:M).*G.dyr(2:L,2:M)));


G.dxu(1:L,1:Mp)=2./(pm(2:Lp,1:Mp)+pm(1:L,1:Mp));
G.dyu(1:L,1:Mp)=2./(pn(2:Lp,1:Mp)+pn(1:L,1:Mp));
G.areaU=scale*(sum(sum(G.dxu(2:Lm,2:M).*G.dyu(2:Lm,2:M)))+ ...
               sum(0.5.*G.dxu(1,2:M).*G.dyu(1,2:M))+ ...
	       sum(0.5.*G.dxu(L,2:M).*G.dyu(L,2:M)));

G.dxv(1:Lp,1:M)=2./(pm(1:Lp,2:Mp)+pm(1:Lp,1:M));
G.dyv(1:Lp,1:M)=2./(pn(1:Lp,2:Mp)+pn(1:Lp,1:M));
G.areaV=scale*(sum(sum(G.dxv(2:L,2:Mm).*G.dyv(2:L,2:Mm)))+ ...
               sum(0.5.*G.dxv(2:L,1).*G.dyv(2:L,1))+ ...
	       sum(0.5.*G.dxv(2:L,M).*G.dyv(2:L,M)));

G.dxp(1:L,1:M)=4./(pm(1:L,1:M )+pm(2:Lp,1:M )+ ...
                   pm(1:L,2:Mp)+pm(2:Lp,2:Mp));
G.dyp(1:L,1:M)=4./(pn(1:L,1:M )+pn(2:Lp,1:M )+ ...
                   pn(1:L,2:Mp)+pn(2:Lp,2:Mp));
G.areaP=scale*(sum(sum(G.dxp(2:Lm,2:Mm).*G.dyp(2:Lm,2:Mm)))+ ...
               sum(0.5.*G.dxp(2:Lm,1).*G.dyp(2:Lm,1))+ ...
	       sum(0.5.*G.dxp(2:Lm,M).*G.dyp(2:Lm,M))+ ...
	       sum(0.5.*G.dxp(1,2:Mm).*G.dyp(1,2:Mm))+ ...
	       sum(0.5.*G.dxp(L,2:Mm).*G.dyp(L,2:Mm))+ ...
	       0.25*G.dxp(1,1)*G.dyp(1,1)+ ...
	       0.25*G.dxp(L,1)*G.dyp(L,1)+ ...
	       0.25*G.dxp(1,M)*G.dyp(1,M)+ ...
	       0.25*G.dxp(L,M)*G.dyp(L,M));
	       
disp(' ');
disp(' Domain area (km^2) computed using inverse metrics pm, pn (old way): ');
disp(' ');
disp(['    Interior rho-points = ', num2str(G.areaR)]);
disp(['    Interior psi-points = ', num2str(G.areaP)]);
disp(['    Interior   u-points = ', num2str(G.areaU)]);
disp(['    Interior   v-points = ', num2str(G.areaV)]);


% Derive quantities at u-, v-, psi-points using: dx, dy

G.DXu(1:L,1:Mp)=0.5.*(dx(2:Lp,1:Mp)+dx(1:L,1:Mp));
G.DYu(1:L,1:Mp)=0.5.*(dy(2:Lp,1:Mp)+dy(1:L,1:Mp));
G.AreaU=scale*(sum(sum(G.DXu(2:Lm,2:M).*G.DYu(2:Lm,2:M)))+ ...
               sum(0.5.*G.DXu(1,2:M).*G.DYu(1,2:M))+ ...
	       sum(0.5.*G.DXu(L,2:M).*G.DYu(L,2:M)));

G.DXv(1:Lp,1:M)=0.5.*(dx(1:Lp,2:Mp)+dx(1:Lp,1:M));
G.DYv(1:Lp,1:M)=0.5.*(dy(1:Lp,2:Mp)+dy(1:Lp,1:M));
G.AreaV=scale*(sum(sum(G.DXv(2:L,2:Mm).*G.DYv(2:L,2:Mm)))+ ...
               sum(0.5.*G.DXv(2:L,1).*G.DYv(2:L,1))+ ...
	       sum(0.5.*G.DXv(2:L,M).*G.DYv(2:L,M)));

G.DXp(1:L,1:M)=0.25.*(dx(1:L,1:M )+dx(2:Lp,1:M )+ ...
                      dx(1:L,2:Mp)+dx(2:Lp,2:Mp));
G.DYp(1:L,1:M)=0.25.*(dy(1:L,1:M )+dy(2:Lp,1:M )+ ...
                      dy(1:L,2:Mp)+dy(2:Lp,2:Mp));
G.AreaP=scale*(sum(sum(G.DXp(2:Lm,2:Mm).*G.DYp(2:Lm,2:Mm)))+ ...
               sum(0.5.*G.DXp(2:Lm,1).*G.DYp(2:Lm,1))+ ...
	       sum(0.5.*G.DXp(2:Lm,M).*G.DYp(2:Lm,M))+ ...
	       sum(0.5.*G.DXp(1,2:Mm).*G.DYp(1,2:Mm))+ ...
	       sum(0.5.*G.DXp(L,2:Mm).*G.DYp(L,2:Mm))+ ...
	       0.25*G.DXp(1,1)*G.DYp(1,1)+ ...
	       0.25*G.DXp(L,1)*G.DYp(L,1)+ ...
	       0.25*G.DXp(1,M)*G.DYp(1,M)+ ...
	       0.25*G.DXp(L,M)*G.DYp(L,M));
		  
disp(' ');
disp(' Domain area (km^2) computed using grid spacing dx, dy (new way): ');
disp(' ');
disp(['    Interior rho-points = ', num2str(G.areaR)]);
disp(['    Interior psi-points = ', num2str(G.AreaP)]);
disp(['    Interior   u-points = ', num2str(G.AreaU)]);
disp(['    Interior   v-points = ', num2str(G.AreaV)]);

return
