function [D,vg,ug] = calc_dynH_fixlev(lon,lat,z,TEMP,SAL,ref_dep)
%By Marlos Goes's method over topography
%Interpolates to the desired section [lon,lat] and 
%calculates the Dynamic height for a section from the 
%temperature e salinity data, referenced in 800 m;
%INPUT: longitude,latitude, DEPTH, temperature, salinity
%D = calc_dynH(lon,lat,z,temp,sal,ref_dep,fig)
%OUTPUT: Dynamic height (in cm)
addpath('/data/mgoes/matlab/Seawater/')
disp('calc_dynH_fixlev')

%Decibars to Pascals
if (nargin<=5 || isempty(ref_dep)),ref_dep = 400;end
[dummy,nn] = min(abs(z-ref_dep));
disp(['Reference level is:', num2str(z(nn))])
dx = sw_dist(lat',lon','km')*1e3;
dx = [dx' dx(end)];

         [Y,Z] = meshgrid(lat,z);
      pressura = Z;%sw_pres(Z,Y);  %Z
      
if ndims(TEMP) ==3,
  %  mask = squeeze(SAL(1,:,:)) ~= 0 & ~isnan(squeeze(SAL(1,:,:)));
    D = zeros(size(TEMP));
    for ii = 1:size(TEMP,1)
        sal = squeeze(SAL(ii,:,:)); temp = squeeze(TEMP(ii,:,:));
       d = calcula_DH(sal,temp,pressura,dx,nn);
       D(ii,:,:) = d;%.*mask;
    end
else
    mask = double(SAL ~= 0 & ~isnan(SAL));
    mask(mask==0)=nan;
    D = calcula_DH(SAL,TEMP,pressura,dx,nn);%.*mask;
end
D = D/10;% * 10;
if nargout > 1
   [ug,vg] = calcula_vel_from_DH(lon,lat,D);
   vg=vg'.*mask;ug=ug'.*mask; %ADDED MAR19
end

function gpan = calcula_DH(Sn,Tn,Pn,dx,kp)
 db2pa = 1e4;
 mask = Sn ~= 0 & ~isnan(Sn);
  mask = double(mask);mask(mask==0)=nan;
  [nn,F] = size(Sn);
  gpans = NaN*ones(nn,F);
  gpan = gpans;

 Pn = Pn.*mask;
 
% encontrando as estacoes com profundidades ate nivel adotado (kp)

kn = (not(isnan(Pn(kp,:))));
ky = find(~kn & ~isnan(mask(1,:)));  %NEW CONSTRAIN TO SURFACE MASK
 
% Calculo da anomalia do geopotencial  referenciado na superficie
% para estacoes com profundidade ate o nivel adotado.

gpans(:,kn) = sw_gpan(Sn(:,kn),Tn(:,kn),Pn(:,kn));

% anomalia do geopotencial referenciado na profundidade de referencia
gpan(:,kn) = ones(nn,1)*gpans(kp,kn)- gpans(:,kn);
         f = find(gpan==0.0);
gpan(f)    = NaN*ones(size(f));

% Calculo da anomalia do geopotencial  referenciado na superficie
% para estacoes com profundidade que nao atingem o nivel adotado.

for iy = ky  %1:kn(1)-1
    ke = find(isnan(Pn(:,iy))); ke = ke(1)-1; %Escolhe ultimo indice valido
   if (ke~=0),
     gpans(1:ke,iy) = sw_gpan(Sn(1:ke,iy),Tn(1:ke,iy),Pn(1:ke,iy));
      gpan(1:ke,iy) = ones(ke,1)*gpans(ke,iy) - gpans(1:ke,iy);
   else
   %      gpan(:,iy) = nan;   %COMMENTED JAN 19
   end
 %  gpan(1:ke,iy) = gpan(1:ke,iy)+ interp1(find(kn),gpans(kp,kn),iy,'nearest','extrap');
end
  
  
  
 % return
% fj -- indice das estacoes que nao alcancam a prof. referencia
% fi -- indice das profundidades maximas que as estacoes fj alcan@am
fj2   = find([kn(1) kn(2:end-1)&(kn(1:end-2)&kn(3:end)) kn(end)]); %find(~isnan(gpan(kp-1,:)));
[fi,fj]    = find(gpan==0.0);
     f     = find(gpan==0.0);
gpan(f)    = NaN*ones(size(f));
gpan(kp,:) = zeros(1,F);
f2 = abs(fj2-fj)';
fk = min(f2);
[~,fki]=sort(fk);
fj=fj(fki);
fi=fi(fki);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% comeca loop para completar e extrapolar  an. geopotencial
%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
lf = length(fi);
      %startp loop on shallow stations
      for l=1:lf
      lj=fj(l); lj2=min(F,fj(l)+2); lj1=min(F,fj(l)+1);
                ljm2=max(1,fj(l)-2); ljm1=max(1,fj(l)-1);
       % cria vetor com gradientes de an. geopotencial 
       %   do par imediatamente anterior de estacoes
       % entre a prof max da estacao em extrapolacao e
       % a profundidade de referencia

       dgpan = gpan(fi(l):kp,lj2)-gpan(fi(l):kp,lj1);  % diferencas avancadas
       dgpan = 1/2/dx(lj1)*dgpan;
        signo=-1;
        gpan(fi(l):kp,lj) = gpan(fi(l):kp,lj1)+signo*dgpan*dx(lj);
       
%       dgpan = gpan(fi(l):kp,lj2)-gpan(fi(l):kp,lj1);  % diferencas centradas
%       dgpan = 1/(dx(lj1)+dx(lj2))*dgpan;
          
       if any(isnan(gpan(fi(l):kp,lj))), signo=1; 
          dgpan = gpan(fi(l):kp,ljm1)-gpan(fi(l):kp,ljm2);  % diferencas atrasadas
          dgpan = 1/2/dx(ljm1)*dgpan;
          gpan(fi(l):kp,lj) = gpan(fi(l):kp,ljm1)+signo*dgpan*dx(lj);
       end
       % cria vetor de valores extrapolados para a estacao l
       % o nivel fi(lj) e kp
       % completa integracao somando o utlimo valor ao resto 
       % da coluna que possui dados
       gpan(1:fi(l)-1,lj) = gpan(1:fi(l)-1,lj) + gpan(fi(l),lj);
      end
return

function [u,v]=calcula_vel_from_DH(lon,lat,D)

if (length(lon)==size(lon,2))
    lon=lon';
end;
if (length(lat)==size(lat,2))
    lat=lat';
end;
nz = size(D,1);
[area,dx,dy]=sw_area(lon,lat);
dy=diag(dy);dx = diag(dx);
widthT = sqrt(dx.^2+dy.^2);
dx = widthT(:,ones(1,nz));
latmat=lat;%*ones(1,length(lon));


%MAKE IT LIKE load_data.m
g  = sw_g(lat);g=g(:,ones(1,nz));  % g=9.81;
pi = 3.14159;
O  = 2.*pi/86400.; % O=7.2921e-05;
f  = 2*O*sin(latmat*pi/180.);   % f=2*O*sin(2*pi*latmat/360);
f=f(:,ones(1,nz));

q=find(abs(latmat)<.001);f(q)=NaN;


[dxE,dyE]=gradient(D);
dxE=dxE'./dx;
%dyE=dyE'./dy;

v = g.*dxE./f;

u=-g.*dxE./f;

return
