% This script examines some of the terms balanced
% in geostrophic and thermal wind momentum balances
% at 26.5N in the ECCOv4r2 latlon_avg
% it is assumed that this file is placed in the same location
% as the data files, which can be downloaded from
% http://fox-kemper.com/data/ecco_for_las/version_4/release2/
% The specific files that are needed are
% oceCURLTAUbyF.nc  (wind stress/f curl)
% oceCURLTAU.nc  (wind stress curl)
% oceTAUE.nc (zonal wind stress)
% llgrid.mat   (grid)
% WVELMASS.nc  (vertical velocity)
% PHIHYD.nc (pressure anomaly)
% NVELMASS.nc (northward velocity)
% ETAN.nc  (sea surface height)
% RHOAnoma.nc (density anomaly)

load llgrid.mat

fname='oceCURLTAUbyF.nc'
readncfile

fname='oceCURLTAU.nc'
readncfile

fname='WVELMASS.nc'
readncfile

fname='PHIHYD.nc'
readncfile

fname='ETAN.nc'
readncfile

fname='RHOAnoma.nc'
readncfile

lonc=lon;
latc=lat;

fname='oceTAUE.nc'
readncfile

lone=lon;

% Note that all of these variables are located at cell center, so lon, lat
% can be overwritten without worry.

fname='NVELMASS.nc'
readncfile

latn=lat;

rho0=1029;
fcor=2*2*pi/(24*3600)*sind(26.5);
beta=2*2*pi/(24*3600)*cosd(26.5)./6370./1000;
dx=mean(mean(diff(lon*114*1000*cosd(26.5))));

ndxc=find(abs(latc(1,:)-26.5780)<0.1);
ndxn=find(abs(latn(1,:)-26.5780)<0.7);

% Focus on only our section at 26.5N

% And we make a 1d longitude for ease of use later.
lon=lonc(:,1);

% Now we choose our longitudes

% Examining this list, we see that indices 100 to 166 
% are uninterrupted Atlantic values.  The westernmost is located
% at 26.5N, 80W, which is right offshore of Boynton Beach, FL
[(ETAN(:,ndxc)) lon (1:length(lon))']


% Note that the 3d variables need 3 indices, with only lat & lon set
% For time-dependent fields, e.g., the climatology, youll have to
% edit these to get the time index, too.

nd=100:166;

oceCURLTAU=squeeze(oceCURLTAU(nd,ndxc));
oceTAUE=squeeze(oceTAUE(99:166,ndxc));
oceTAUE=0.5*(oceTAUE(2:end,:)+oceTAUE(1:end-1,:)); % average to midpoint
oceCURLTAUbyF=squeeze(oceCURLTAUbyF(nd,ndxc));
ETAN=squeeze(ETAN(nd,ndxc));
PHIHYD=squeeze(PHIHYD(nd,ndxc,:));
WVELMASS=squeeze(WVELMASS(nd,ndxc,:));
RHOAnoma=squeeze(RHOAnoma(nd,ndxc,:));

% This one we average the neighboring values onto our latitude.
NVELMASS=squeeze(mean(NVELMASS(nd,ndxn,:),2));
lon=lon(nd);

figure(1)
clf
plot(lon,ETAN)
xlabel('longitude (deg)')
ylabel('Sea Surface Height (m)')


figure(2)
subplot(211)
pcolor(lon,llgrid.RC(1:29),NVELMASS(:,1:29)');
shading('flat');
caxis([-.3 .3])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward Velocity (m/s)')

subplot(212)
pcolor(lon,llgrid.RC(30:50),NVELMASS(:,30:50)');
shading('flat');
caxis([-.04 .04])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward Velocity (m/s)')

PRESS=9.81*ETAN*ones(1,50)-PHIHYD;

NVELGEO=diff(PRESS)./mean(diff(lon*114*1000*cosd(26.5)))./fcor;
lonh=0.5*(lon(2:end)+lon(1:end-1));

subplot(212)
pcolor(lonh,llgrid.RC(30:50),NVELGEO(:,30:50)');
shading('flat');
caxis([-.04 .04])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward Geostrophic Velocity (m/s)')


H=400;  %Assumed depth
NVELGEOSURF=9.81*diff(ETAN)./dx./fcor*H;

NVELINT=nansum(NVELMASS(:,:).*(ones(67,1)*llgrid.DRF'),2);

NVELEKMAN=oceTAUE./fcor./rho0;

NVELSVER=oceCURLTAU./rho0./beta;

NVELSVERPLUS=fcor.*oceCURLTAUbyF./rho0./beta;

figure(3)
clf
plot(lon,NVELINT,lonh,NVELGEOSURF,lon,NVELEKMAN,lon,NVELSVER,lon,NVELSVERPLUS);

xlabel('?')
ylabel('?')
legend('?','??','???','????','?????')


dvdztw=-9.81*diff(RHOAnoma)./dx./fcor./rho0;

figure(4)
subplot(211)
pcolor(lonh,llgrid.RC(1:29),dvdztw(:,1:29)');
shading('flat');
caxis([-5e-5 5e-5])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward Thermal Wind Velocity Shear (1/s)')

subplot(212)
pcolor(lonh,llgrid.RC(30:50),dvdztw(:,30:50)');
shading('flat');
caxis([-5e-5 5e-5])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward Thermal Wind Velocity Shear (1/s)')

figure(5)
zh=0.5*(llgrid.RC(2:end)+llgrid.RC(1:end-1));
dvdz=diff(NVELMASS,[],2)./(ones(67,1)*diff(llgrid.RC)');

subplot(211)
pcolor(lon,zh(1:29),dvdz(:,1:29)');
shading('flat');
caxis([-5e-5 5e-5])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward ECCO Velocity Shear (1/s)')

subplot(212)
pcolor(lon,zh(30:49),dvdz(:,30:49)');
shading('flat');
caxis([-5e-5 5e-5])
colorbar
xlabel('longitude')
ylabel('depth (m)')
title('Northward ECCO Wind Velocity Shear (1/s)')
