% eval_grid
%
% calculate the grid of 3x3 degree points that the Argo array
% will sample.  Also calculate a higher resolution grid that
% covers the same region, in order to require a float within 3 deg.
% of all open-ocean regions.

for index=1:2;
    
    if (1)  % TPOS 
        if (index==1) d=2.1; else d=2.1/2;end;
        long=(125:d:300)'; %(d/2:d:360-d/2)'; 
        latg=-10:d:10;
    else  % global
        if (index==1) d=3; else d=3/2; end;
        long=(d/2:d:360-d/2)'; 
        latg=-60:d:70;
         if (index==2) long=[long;360]; end;
    end;
    long=long*ones(1,length(latg));
    latg=ones(size(long,1),1)*latg;
    
    % remove points in:
    % Hudson Bay
    i=[258   257   294   294]; j=[79    49    51    79];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    i=[293   299   304   294   288];j=[46    46    52    54    50];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Med, Red and Black seas
    i=[-4    11    56    56    -1]; j=[40    49    48    20    21];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    i=[42    46    36    31];j=[11    18    29    25];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    i=[354 354 360 360];j=[30 43 43 30];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Indonesian Seas
    i=[133   106    99   100   122   126   142];
    j=[-11    -7     8    15    14     0    -5];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    q=find(long>135 & long<143 & latg>-18 & latg<-8);
    long(q)=[]; latg(q)=[];
    % East Sea
    i=[ 104   121   120   101]; j=[13    14    30    29];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    i=[116   128   132   122   108];
    j=[27    26    34    45    36];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Japan Sea
    i=[127   127   141   143   142   139   131];
    j=[36    44    52    52    43    36    33];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Baltic
    i=[ 14    14    24    33    28]; j=[52    64    67    61    52];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Okhotsk
    i=[ 133   143   159   166   148];j=[58    43    55    64    66];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Arctic
    i=[25    27   316   328   328   328];
    j=[90    65    65    71    72    90];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
    % Indonesian Seas
    i=[165   157   150   137   130   128   122   121   130   143   157   164];
    j=[-11    -4    -1     0     6    11    10     2    -6   -10   -11   -13];
    q=find(isinpoly(long,latg,i,j)); long(q)=[]; latg(q)=[];
  
    
    load /phodnet/drifter/lumpkin/bathystuff/GINA/bathyworld.mat
    lon2=[lon2,lon2+360];
    z2=[z2,z2];
    
    h=interp2(lon2,lat2,z2,long,latg);
    q=find(h<-2000);
    long=long(q)';
    latg=latg(q)'; 
    
    % remove points close to land 
    lon2m=(lon2'*ones(size(lat2)))';
	lat2m=(ones(length(lon2),1)*lat2)';
	q=find(z2>=0);
    del=10*ones(size(long));
    for j=1:length(long);
        del(j)=sqrt(min( (lon2m(q)-long(j)).^2 + (lat2m(q)-latg(j)).^2));
    end;  
    q=find(del>1); 
    long=long(q); latg=latg(q);
    
    plot(long,latg,'.');
    overlay; axis([0 360 -90 90]);
    title(num2str(length(long)),'fontsize',12);
    
    pause(1)
    
    if (index==1) long2=long; latg2=latg; end;

end;
 
plot(long,latg,'.',long2,latg2,'ok'); 
 
save eval_grid long* latg*  


