function [D,vg] = calc_dynH_fixlev(lon,lat,z,TEMP,SAL,ref_dep)
%This uses Ilson'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)
%dbstop('43')
addpath('/automount1/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)
 %       if ii==27,keyboard,end
        sal = squeeze(SAL(ii,:,:)); temp = squeeze(TEMP(ii,:,:));
       d = calcula_DH(sal,temp,pressura,dx,nn);
       D(ii,:,:) = d;%.*mask;
    end
else
   % mask = SAL ~= 0 & ~isnan(SAL);
    D = calcula_DH(SAL,TEMP,pressura,dx,nn);%.*mask;
end
D = D * 10;
if nargout > 1
[vg] = calcula_vel_from_DH(lon,lat,D);
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;
 %Pn(Pn==0) = nan;
 
% 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
 %ky = find(~kn);
 
% 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
  %    disp(num2str(iy))
   if (ke>1) %~=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

[fi,fj]    = find(gpan==0.0);
     f     = find(gpan==0.0);
gpan(f)    = NaN*ones(size(f));
gpan(kp,:) = zeros(1,F);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% comeca loop para completar e extrapolar  an. geopotencial
%
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

lf = length(fi);
%fi = flipud(fi);  % inverte para comecar pelas estacoes mais fundas 
%fj = flipud(fj);
[as,sa]=sort(abs(fj-F/2));
fi = fi(sa); fj = fj(sa);
%fj = [sort(fj(fj>(f/2))) fj(fj<(f/2))]

      %for l=lf:-1:1 %this was before, but started from outside since fi,fj
      %are inverted
      for l=1:lf
%      lj=fj(l); lj2=min(F,fj(l)+2); lj1=min(F-1,fj(l)+1);
%                ljm2=max(1,fj(l)-2); ljm1=max(2,fj(l)-1);
%REPLACED BY: ESPETACULO
      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 isnan(gpan(fi(l),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
       
       %gpan(fi(l):kp,lj) = gpan(fi(l):kp,lj1)-dgpan*dx(lj);
   %     gpan(fi(l):kp,lj) = gpan(fi(l):kp,lj1)+signo*dgpan*dx(lj);     %commented Jan19
      
       % completa integracao somando o utlimo valor ao resto 
       % da coluna que possui dados
     %  if ~isnan(gpan(fi(l),lj))
       gpan(1:fi(l)-1,lj) = gpan(1:fi(l)-1,lj) + gpan(fi(l),lj);
       %else
           %keyboard
  %         gpan(1:kp,lj) = interp1([1:fi(l)-1 kp],gpan([1:fi(l)-1 kp],lj),1:kp);
  %         gpan(kp:-1:1,lj) = cumsum(gpan(kp:-1:1,lj));
    %   end
      end
      
      %IT was good for 800m, but gets weird afterwards. 
      %Need to reference all at the surface, calculate the dgpan and then
      %extrapolate
      %uu = (~isnan(gpan(1,:)));
%uu_int = find(~uu&~isnan(mask(1,:)));
% gpan(1,uu_int) = interp1(find(uu),gpan(1,uu),uu_int,'nearest','extrap');



return

function [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));
%dx = dx(:,ones(1,nz));
%[dy]=sw_dist(lon,lat,'km')*1e3;dy(end+1) = dy(end);%dy = dy';
latmat=lat;%*ones(1,length(lon));

% O=7.2921e-05;
% f=2*O*sin(2*pi*latmat/360);
% g=9.81;

%MAKE IT LIKE load_data.m
g  = sw_g(lat);g=g(:,ones(1,nz));
pi = 3.14159;
O  = 2.*pi/86400.;
f  = 2*O*sin(latmat*pi/180.);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*dyE./f;




return
