function [ncells,firstcell]=equare(dlat)

%======computes the number of cells in a lattitude band
%======and the first cellnumber of a lattitude band

%% EQUAL-AREA COMPUTATIONS
%% totcel   total number of E.A. boxes
%% ncells(720) number of E.A. boxes per lat zone

maxlat=180/dlat;
maxlon=360/dlat;

ncells = zeros(maxlat,1);
firstcell = zeros(maxlat,1);

%% common /equcom/tocell,ncells
rearth = 6371.2;
rcelat=(dlat*pi)/180;
totcel=0;
%% calculate height and area of equatorial zone
hezon=rearth*sin(rcelat);
aezon=2*pi*rearth*hezon;
%% calculate area of equatorial cell
aecell=(aezon*dlat)/360;
%% print*,'aecell',aecell
%% compute longitude zones for each latitude zone
maxlt2=maxlat/2;
for lat=1:maxlt2

  xlatb=(lat-1)*dlat;
  xlate=xlatb+dlat;
  rlatb=(2*pi*xlatb)/360;
  rlate=(2*pi*xlate)/360;
  %% calculate heights of latb,late,zone
  htb=rearth*sin(rlatb);
  hte=rearth*sin(rlate);
  htzone=hte-htb;
  azone=2*pi*rearth*htzone;
  %% calculate number of cells
  rcells=azone/aecell;
  icellr=fix(rcells+0.5);
  %% augment total # grid cells (both hemispheres)
  totcel=totcel+2*icellr;
  rcellr=icellr;
  dlongr=360/rcellr;
  acell=azone/rcellr;
  asq=azone/maxlon;
  %% create table of longitudes
  lat1=lat+maxlt2;
  lat2=maxlt2+1-lat;
  ncells(lat1)=icellr;
  ncells(lat2)=icellr;

end
    
%% search for the first cellnumber in each lat band
firstcell(1)=1;
for i=2:maxlat
  firstcell(i)=firstcell(i-1)+ncells(i-1);
end
    
