% Example matlab program using gcmfaces to display state estimate.

% This example 
% 0) matlab set up 
% 1) Reads in the v4 grid, 
% 2) Reads in the two horizontal components of velocity over the 
%    entire model domain at a particular month, and vertical level 
% 3) Computes curl of UV velocities
% 4) Plots the curl computed in 3)

%pre-requisites : 
% 1) download and expand 
%      ancillary_data/gcmfaces_MITprof_r1.tar
%      ancillary_data/gcmfaces_MITprof_sample_input_r1.tar
%      ancillary_data/GRID_r1.tar
%    in directory where present routine was put
% 2) have Matlab 2012a or newer that has built-in support of OPeNDAP

%note : has not yet been tested with MS windows

% -------------------------------------------------------------------------
% Step 0a) 
% Include gcmfaces package in path. 
% This example assumes the package is in directory ./gcmfaces/
% Change according to where the package actually is. 

if ~isdir('gcmfaces'); error('missing gcmfaces/ directory'); end;
path(genpath([pwd '/gcmfaces/']),path);

% -------------------------------------------------------------------------
% Step 0b)
% It's often a good idea to clear and close all before proceeding.

clear all
close all

% -------------------------------------------------------------------------
% Step 0c)
% Initialize the gcmfaces package.

gcmfaces_global;

% -------------------------------------------------------------------------
% Step 1) 
% Set grid parameters and read in grid information for v4.
% This example assumes the grid files are in directory ./GRID/
% Change according to where the grid files actually are.

dirGrid='GRID/';   % Directory expanded from GRID_r1.tar
if ~isdir(dirGrid); error('missing GRID/ directory'); end;
nF=5; fileFormat='compact';  %format specification
grid_load(dirGrid,nF,fileFormat);

% -------------------------------------------------------------------------
% Step 2)
% Read in U and V velocities. The velocities are in model native 
% coordinate and are NOT necessarily zonal or meridional components.

%This expample reads in U and V in netcdf format from MIT's OPeNDAP server 
%(need Matlab 2012a or newer that has built-in support of OPeNDAP)
uname = 'http://mit.ecco-group.org/opendap/ecco_for_las/version_4/release1/nctiles/UVELMASS/UVELMASS';
vname = 'http://mit.ecco-group.org/opendap/ecco_for_las/version_4/release1/nctiles/VVELMASS/VVELMASS';

%Beyond this example, downloading the netcdf files to a local directory is highly 
%recommended, since it will then be much faster to read in. Then one would use :
%uname = 'UVELMASS/UVELMASS';
%vname = 'VVELMASS/VVELMASS';

fprintf('reading UVELMASS\n');
fld_u = read_nctiles(uname,'UVELMASS',1,1);
fprintf('reading VVELMASS\n');
fld_v = read_nctiles(vname,'VVELMASS',1,1);

% -------------------------------------------------------------------------
% Step 3)
% Compute curl of UV velocities.

uvcurl=calc_UV_curl(fld_u,fld_v,1,1);


% -------------------------------------------------------------------------
% Step 4)
% Plot the curl 

[x,y,z]=convert2pcol(mygrid.XC,mygrid.YC,uvcurl);   
figure;
pcolor(x,y,z);shading flat;
axis([-180 180 -90 90])
colorbar; caxis([-6e-6 6e-6]);

