function [nx,ny,nz] = grad_locally_referenced(sr,ct,pr,lon,lat,zlev,ra)
%
% FUNCTION [nx, ny, nz] = grad_neutral_density(sr,ct,lon,lat,zlev,ra)
%
% Compute the gradient of the locally referenced potential density
%
% INPUTS:
%   sr(npi,npj,npk) : 3-dimensional reference composition salinity (g/kg)
%   ct(npi,npj,npk) : 3-dimensional Conservative Temperature (deg C)
%   pr(npi,npj,npk) : 3-dimensional pressure field (dbar) 
%   lon(npi)        : 1-dimensional longitude array (deg)
%   lat(npj)        : 1-dimensional latitude array (deg)
%   zlev(npk)       : 1-dimensional array of positive depth (meters)
%   ra              : scalar Earth Radius (meters)
%
% OUTPUTS:
%   nx(npi,npj,npk) : zonal gradient of locally referenced potential dens
%   ny(npi,npj,npk) : meridional gradient of locally referenced pot dens
%   nz(npi,npj,npk) : vertical gradient of locally referenced pot dens
%
% 
%
% Author: Remi Tailleux, University of Reading, 14 October 2018
%                                      Modified: 15 September 2020
% -------------------------------------------------------------------------
[npi,npj,npk] = size(sr); 
nx = zeros(npi,npj,npk); ny = zeros(npi,npj,npk); nz = zeros(npi,npj,npk);
dx = zeros(npj,npk); 
valid = sr > 0. & sr < 45.; 
sr(~valid) = NaN; ct(~valid) = NaN; pr(~valid) = NaN; 

% Compute mid density
rho_mid = gsw_rho_CT_exact(sr,ct,pr) -1000. ; 

% Compute grid size elements assuming uniform grid size 
dlon = (lon(2)-lon(1))*pi/180.; 
for j=1:npj
    dx(j,:) = ra*cos(lat(j)*pi/180.)*dlon; 
end  
dy = ra*(lat(2)-lat(1))*pi/180.;

% Compute the vertical gradients 
% ------------------------------
% Interior points 
for k=2:npk-1
    rho_upper = ...
        gsw_rho_CT_exact(sr(:,:,k-1),ct(:,:,k-1),pr(:,:,k)) - 1000.;
    rho_lower = ...
        gsw_rho_CT_exact(sr(:,:,k+1),ct(:,:,k+1),pr(:,:,k)) - 1000.;
    dz1 = zlev(k)-zlev(k-1);
    dz2 = zlev(k+1) - zlev(k);
    dz = dz1+dz2;
    nz(:,:,k) = (dz2/dz)*(rho_upper - rho_mid(:,:,k))./dz1 + ...
                (dz1/dz)*(rho_mid(:,:,k) - rho_lower)./dz2 ; 
end
% Top level
% pr_mid = 0.5*(pr(:,:,1)+pr(:,:,2));
pr_mid = pr(:,:,1); 
rho_upper = ...
    gsw_rho_CT_exact(sr(:,:,1),ct(:,:,1),pr_mid) - 1000.;
rho_lower = ...
    gsw_rho_CT_exact(sr(:,:,2),ct(:,:,2),pr_mid) - 1000. ;
nz(:,:,1) = (rho_upper - rho_lower)/(zlev(2)-zlev(1));
% Lowest level
% pr_mid = 0.5*(pr(:,:,npk-1)+pr(:,:,npk));
pr_mid = pr(:,:,npk); 
rho_upper = ...
    gsw_rho_CT_exact(sr(:,:,npk-1),ct(:,:,npk-1),pr_mid) - 1000.;
rho_lower = ...
    gsw_rho_CT_exact(sr(:,:,npk),ct(:,:,npk),pr_mid) - 1000.;
nz(:,:,npk) = (rho_upper - rho_lower)/(zlev(npk)-zlev(npk-1)); 

% Compute the latitudinal gradients 
% ---------------------------------
% Interior points 
for j=2:npj-1
    rho_north = ...
        gsw_rho_CT_exact(sr(:,j+1,:),ct(:,j+1,:),pr(:,j,:)) - 1000. ;
    rho_south = ...
        gsw_rho_CT_exact(sr(:,j-1,:),ct(:,j-1,:),pr(:,j,:)) - 1000. ; 
    ny(:,j,:) = (rho_north - rho_south)/(2.*dy);
end
% Northermost point
pr_mid = 0.5*(pr(:,npj-1,:)+pr(:,npj,:));
rho_north = ...
    gsw_rho_CT_exact(sr(:,npj,:),ct(:,npj,:),pr_mid) - 1000.;
rho_south = ...
    gsw_rho_CT_exact(sr(:,npj-1,:),ct(:,npj-1,:),pr_mid) - 1000.;
ny(:,npj,:) = (rho_north - rho_south)/dy;
% Southernmost point
pr_mid = 0.5*(pr(:,1,:) + pr(:,2,:));
rho_north = ...
    gsw_rho_CT_exact(sr(:,2,:),ct(:,2,:),pr_mid) - 1000.;
rho_south = ...
    gsw_rho_CT_exact(sr(:,1,:),ct(:,2,:),pr_mid) - 1000. ;
ny(:,1,:) = (rho_north - rho_south)/dy;

% Compute the longitudinal gradients
% ----------------------------------
% Interior points
for i=2:npi-1
    rho_east = gsw_rho_CT_exact(squeeze(sr(i+1,:,:)), ...
        squeeze(ct(i+1,:,:)),squeeze(pr(i,:,:))) - 1000. ;
    rho_west = gsw_rho_CT_exact(squeeze(sr(i-1,:,:)), ...
        squeeze(ct(i-1,:,:)), squeeze(pr(i,:,:))) - 1000. ;
    nx(i,:,:) = (rho_east - rho_west)./(2.*dx);
end
disp('size rho_east and west')
disp(size(rho_east))
disp(size(dx)) 
disp(max(dx(:)))
% Easternmost point
rho_east = gsw_rho_CT_exact(squeeze(sr(1,:,:)), ...
    squeeze(ct(1,:,:)), squeeze(pr(npi,:,:))) - 1000. ;
rho_west = gsw_rho_CT_exact(squeeze(sr(npi-1,:,:)), ...
    squeeze(ct(npi-1,:,:)), squeeze(pr(npi,:,:))) - 1000. ;
nx(npi,:,:) = (rho_east-rho_west)./(2.*dx);
% Westernmost point
rho_east = gsw_rho_CT_exact(squeeze(sr(2,:,:)), ...
    squeeze(ct(2,:,:)), squeeze(pr(1,:,:))) - 1000. ;
rho_west = gsw_rho_CT_exact(squeeze(sr(npi,:,:)), ...
    squeeze(ct(npi,:,:)), squeeze(pr(1,:,:))) - 1000. ;
nx(1,:,:) = (rho_east - rho_west)./(2.*dx); 

% Remove the North Pole point to NaN
%----------------------------
nx(:,npj,:) = NaN; ny(:,npj,:)=NaN; nz(:,npj,:)=NaN; 




