% This computes the OAX grid file using the "old" method of
% scale definition

scxy=loadgrid('../MESH/sab_clim2');
scll=loadgrid('../MESH/sab_clim2ll');
scxy=scll;

trx0=-7.672028358888956e+06;
try0=3.128247751710321e+06;
lo0=-80.75;
la0=31.25;
[scxy.x,scxy.y]=convll2m(scll.x,scll.y,lo0,la0);
scxy.x=scxy.x-trx0;
scxy.y=scxy.y-try0;

data=load('super_smoothed_z_sab_clim2.dat');
smoothz=data(:,1);


XX=scxy.x;
YY=scxy.y;

gh=grad(scxy.e,scxy.x,scxy.y,smoothz);
hx=gh(:,1);hy=gh(:,2);

ZZ=scxy.z;
[sm,a]=vecmagdir(hx,hy);

%Make the horiz. correlation scales 
% bottom-depth dependent

i70=find(ZZ<=70);
i70to600=find(ZZ>70&ZZ<=600);
i600plus=find(ZZ>600);

% Make the oax .grid file (correlation scales)
Cx=NaN*(ones(size(ZZ)));
Cy=Cx;
S0=Cx;

Cx(i70)=3;
Cy(i70)=10;
S0(i70)=100;

Cx(i70to600)=1;
Cy(i70to600)=10;
S0(i70to600)=100;

Cx(i600plus)=10;
Cy(i600plus)=10;
S0(i600plus)=100;


Sxp=S0.*(1+Cx.*sqrt(sm))*1000;
Syp=S0.*(1+Cy.*sqrt(sm))*1000;

% Load neighbor file for sab_clim mesh
[x,y,bc,z,nbs]=read_nei('/afs/isis/depts/marine/opnml/ARCHIVE/GRIDS/sab_clim2/sab_clim2ll.nei');
for i=1:5
 Sxp=avg_over_nei(scll,nbs,Sxp);
 Syp=avg_over_nei(scll,nbs,Syp);
end

%z800=load('sab_clim2_smooth_800max.bat.BOB');
%write_bat(z800,'sab_clim2_smooth_800max.bat')
z800=scxy.z;

nnv=21;
ZZ=flipud(sinegrid(z800,zeros(size(scxy.z)),nnv,1.)');
XX=repmat(scxy.x',[nnv 1]);
YY=repmat(scxy.y',[nnv 1]);

SZ=NaN*ones(size(YY));
SM=repmat(sm',[nnv 1]);
A=repmat(a',[nnv 1]);
SXP=repmat(Sxp',[nnv 1]);
SYP=repmat(Syp',[nnv 1]);

ZBOT=repmat(z800',[nnv 1]);

%i100=find(ZBOT<=100);
%i500=find(ZBOT>100&ZBOT<=500);
%i500plus=find(ZBOT>500);
%SZ(i100)=25;
%SZ(i500)=50;
%SZ(i500plus)=100;
%find(isnan(SZ))

% This is the "old" method for vert scales, 
% based on the bvertical node depth, not the total
% water column depth
i1=find(ZZ>-100);
SZ(i1)=10;
i2=find( ZZ<=-100 & ZZ>-500);
SZ(i2)=25;
i3=find( ZZ<=-500);
SZ(i3)=50;

find(isnan(SZ))


%ST=30*ones(size(YY));       % Time correlation scale (days)
%T=45*ones(size(YY));        % Time for estimation (Year-day)


%  oaxgrid output for sigma grid
out=[XX(:) YY(:) abs(ZZ(:)) A(:) SXP(:) SYP(:) SZ(:)];
save clim_fem_sg_notime.oaxgrid out -ascii



% Make fem oax file on level surfaces.
zlevs=[0:3:12 20:20:100 150:50:500 600:100:800];
ZLEVS=repmat(zlevs,[1 length(XX)]);
ZLEVS=ZLEVS(:);
nzlevs=length(zlevs);
XXX=repmat(scxy.x,[1 nzlevs])';
XXX=XXX(:);
YYY=repmat(scxy.y,[1 nzlevs])';
SZ=NaN*ones(size(YYY));
YYY=YYY(:);
SM=repmat(sm,[1 nzlevs])';
SM=SM(:);
A=repmat(a,[1 nzlevs])';
A=A(:);
SXP=repmat(Sxp,[1 nzlevs])';
SXP=SXP(:);
SYP=repmat(Syp,[1 nzlevs])';
SYP=SYP(:);

for i=1:length(zlevs)
   if zlevs(i)<100
      SZ(i,:)=10;
   elseif  zlevs(i)>=100 & zlevs(i)<500
      SZ(i,:)=25;
   else
      SZ(i,:)=50;
   end
end
SZ=SZ(:);

%ZBOT=repmat(scxy.z,[1 nzlevs])';
%i100=find(ZBOT<=100);
%i500=find(ZBOT>100&ZBOT<=500);
%i500plus=find(ZBOT>500);
%SZ(i100)=25;
%SZ(i500)=50;
%SZ(i500plus)=100;
%find(isnan(SZ))

out=[XXX(:) YYY(:) abs(ZLEVS(:)) A(:) SXP(:) SYP(:) SZ(:)];
save clim_fem_zlevs_notime.oaxgrid out -ascii