
HisFile='/home/arango/ocean/toms/adjoint/Test/SCB/Forward/scb20_his_01.nc';
GrdFile='/home/arango/ocean/toms/adjoint/Test/SCB/Data/scb20_grd.nc';

Rec=20;

% Read in Coriolis parameter.

f=nc_read(GrdFile,'f');
rmask=nc_read(GrdFile,'mask_rho');

% Read in velocity components, temperature and salinity.

u=nc_read(HisFile,'u',Rec);
v=nc_read(HisFile,'v',Rec);
T=nc_read(HisFile,'temp',Rec);
S=nc_read(HisFile,'salt',Rec);

% Compute depths at RHO-points.

[Zr]=depths(HisFile,GrdFile,1,0,Rec);
[Zw]=depths(HisFile,GrdFile,5,0,Rec);

[Lp,Mp,N]=size(Zr);
L=Lp-1; Lm=L-1;
M=Mp-1; Mm=M-1;

% Compute Brunt-Vaisala frequency.

[N2]=bvf(T,S,Zr,Zw);

% Compute velocity shear at horizontal RHO-points and vertical W-points.

Ur=0.5.*(u(1:Lm,1:Mp,1:N)+u(2:L ,1:Mp,1:N));
Vr=0.5.*(v(1:Lp,1:Mm,1:N)+v(1:Lp,2:M ,1:N));

for k=1:N-1,
  dudz(1:Lm,1:Mp,k+1)=(Ur(1:Lm,1:Mp,k+1)-Ur(1:Lm,1:Mp,k  ))./ ...
                      (Zr(2:L ,1:Mp,k+1)-Zr(2:L ,1:Mp,k  ));

  dvdz(1:Lp,1:Mm,k+1)=(Vr(1:Lp,1:Mm,k+1)-Vr(1:Lp,1:Mm,k  ))./ ...
                      (Zr(1:Lp,2:M ,k+1)-Zr(1:Lp,2:M ,k  ));
end,

%  Bottom and surface boundary conditions.

dudz(:,:,1)=0;   dudz(:,:,N+1)=0;
dvdz(:,:,1)=0;   dvdz(:,:,N+1)=0;

dudz2=dudz.^2;
dvdz2=dvdz.^2;

% Convert Coriolis parameter to a 3D array.

f=repmat(f,[1 1 N+1]);

% Compute desired quantity at interior RHO-points.

Q=f(2:Lm,2:Mm,:).*sqrt(dudz2(1:Lm-1,2:Mm,:) + ...
                       dvdz2(2:Lm,1:Mm-1,:));

Q(:,:,2:N)=Q(:,:,2:N)./N2(2:Lm,2:Mm,2:N);

% Plot

pcolor(Q(:,:,20)'); shading interp; colorbar;
