function [A,lon,lat,lev]=lect_regcm_sri25(flon,flat,flev,ivar,month,day)
% [A,lon,lat,lev]=lect_regcm_sri25(flon,flat,flev,ivar,month,day);
% function to read 25km RegCM daily outputs from Joshua Qian
% flon,flat,flev: pairs of coordinate values defining desired window
% ivar: =1 for u; =2 for v; =4 for q (see grads .ctl file) 
% month,day: pairs of values values defining desired time window

nlon=59; nlat=59; nlev=18;
x=72.5926:0.2248:85.6310;
y=[2.12   2.35   2.57   2.80   3.02   3.25   3.47   3.70   3.92   4.14 ...
   4.37   4.59   4.82   5.04   5.26   5.49   5.71   5.94   6.16   6.38 ...
   6.61   6.83   7.05   7.28   7.50   7.72   7.94   8.17   8.39   8.61 ...
   8.83   9.06   9.28   9.50   9.72   9.94  10.16  10.39  10.61  10.83 ...
  11.05  11.27  11.49  11.71  11.93  12.15  12.37  12.59  12.81  13.03 ...
  13.25  13.46  13.68  13.90  14.12  14.34  14.55  14.77  14.99];
z=[0.995 0.985 0.97 0.945 0.91 0.865 0.81 0.745 0.67 0.59 0.51 0.43 ...
  0.35 0.27 0.195 0.135 0.075 0.025];
[X,Y,Z]=meshgrid(x,y,z);
%X=X'; Y=Y'; Z=Z';
window=find(X >= flon(1) & X <= flon(2) & Y >= flat(1) & Y <= flat(2) ...
          & Z >= flev(2) & Z <= flev(1));
yb=1980:4:2006; % leap years
YR0=1978;   % first year minus 1
YEAR=1:28;
%YEAR=1:2;

cd /mstore/iri/modelling/jqian/RegCM3/tino/sri/postday
addpath /beluga/data/vincent/iri_html/matlab_function_version1

A=[];
for i=1:length(YEAR);
  disp(['year #',num2str(i),' ...']);
  fid=fopen(['ATM',num2str(YEAR(i)+YR0),'.day.nr.AVG.DAT'],'r');
  fseek(fid,(ivar-1)*nlev*nlon*nlat*4,'bof'); % starting byte posn in file
  a(:,1)=fread(fid,nlon*nlat*nlev,'float32','b');
  if isempty(find((i+YR0)==yb)); % non-leap yrs
    for j=2:365;
      fseek(fid,(7*nlev+5)*nlon*nlat*4,'cof'); % skip 7*18 lev + 5*1 lev
      a(:,j)=fread(fid,nlon*nlat*nlev,'float32','b');
    end
    [indice,l]=choose_date_daily([YEAR(i)+YR0,1,1,0],[YEAR(i)+YR0 YEAR(i)+YR0],month,day,0);
    a=a(window,indice);
    fclose(fid);
    [ra,ca]=size(a);
  else
    for j=2:366;
      fseek(fid,(7*nlev+5)*nlon*nlat*4,'cof'); % skip 7*18 lev + 5*1 lev
      a(:,j)=fread(fid,nlon*nlat*nlev,'float32','b');
    end
    fclose(fid);
    [indice,l]=choose_date_daily([YEAR(i)+YR0,1,1,0],[YEAR(i)+YR0 YEAR(i)+YR0],month,day,0);
    a=a(window,indice);
    [ra,ca]=size(a);
  end
  a2=a';
  [R,C]=size(A);
  A(R+1:R+ca,:)=a2;
  clear a a2
end   
lon=X(window);
lat=Y(window);
lev=Z(window);



