%%0) setup:

%add matlab path and load ECCO v4 grid:
p = genpath([pwd filesep 'gcmfaces/']); addpath(p);
p = genpath([pwd filesep 'MITprof/']); addpath(p);
p = genpath([pwd filesep 'm_map/']); addpath(p);
grid_load; gcmfaces_global;
m=mygrid.mskC(:,:,1);%standard land mask

%input directories:
dirNctiles='release1/nctiles_climatology/';
dirBinary='release1/binary_inputs/';

%constant parameters:
parms.rhoconst =1029; %sea water density [in kg/m^3]
parms.rhoi     = 910; %sea ice density [in kg/m^3]
parms.rhosn    = 330; %snow density [in kg/m^3]


%%1) sea surface temperature and related variables:

%THETA:    potential temperature [in degree C]
%SALT:     salinity [in psu]
%SIarea:   fractional ice-covered area [0 to 1]
%MXLDEPTH: Mixed-Layer Depth [in m]
%
%note:     MXLDEPTH is based on the Kara, Rochford, and Hurlburt 2000 formula

if isdir(dirNctiles);
sst_ecco_clim=read_nctiles([dirNctiles 'THETA/THETA'],'THETA',[],1);
sss_ecco_clim=read_nctiles([dirNctiles 'SALT/SALT'],'SALT',[],1);
mld_ecco_clim=read_nctiles([dirNctiles 'MXLDEPTH/MXLDEPTH'],'MXLDEPTH');
ice_ecco_clim=read_nctiles([dirNctiles 'SIarea/SIarea'],'SIarea');
end;

%sst_reynolds:    Reynolds monthly mean sea surface temperature maps
%iicecover_nsidc: NSIDC monthly mean fractional ice-covered area [0 to 1]
%
%note:            references are provided in Forget et al 2015 (GMD)

if isdir(dirBinary);
sst_reynolds_1992=read_bin([dirBinary 'sst_reynolds/reynolds_oiv2_r1_1992']);
sst_reynolds_1992(sst_reynolds_1992==0)=NaN;%0's denote missing values in these files
ice_nsidc_1992=read_bin([dirBinary 'icecover_nsidc/nsidc79_monthly_1992']);
ice_nsidc_1992(ice_nsidc_1992==-9999)=NaN;%-9999's denote missing values in these files
end;


%%2) heat and freshwater fluxes into the ocean (and ice pack)

%oceFWflx: net surface Fresh-Water flux into the ocean [in kg/m^2/s]
%oceQnet:  net surface heat flux into the ocean [in W/m^2]
%SIatmFW:  net freshwater flux from atmosphere & land [in kg/m^2/s]
%SIatmQnt: net heat flux to atmosphere & land [in W/m^2]
%
%note:     the key difference between oceQnet and SIatmQnt (or oceFWflx and SIatmFW) is
%          that the former is computed at the bottom of the ice pack when there is sea-ice
%note:     the MITgcm sign conventions are opposite for SIatmQnt and oceQnet so below
%          I flip the sign of SIatmQnt upon reading it from file.

if isdir(dirNctiles);
qnet_to_oce=read_nctiles([dirNctiles 'oceQnet/oceQnet'],'oceQnet');
qnet_to_surf=-read_nctiles([dirNctiles 'SIatmQnt/SIatmQnt'],'SIatmQnt');
qnet_to_ice=qnet_to_surf-qnet_to_oce;

fwf_to_oce=read_nctiles([dirNctiles 'oceFWflx/oceFWflx'],'oceFWflx');
fwf_to_surf=read_nctiles([dirNctiles 'SIatmFW/SIatmFW'],'SIatmFW');
fwf_to_ice=fwf_to_surf-fwf_to_oce;
end;


%%3) sea surface height and bottom pressure:

%ETAN:     Surface Height Anomaly [in m]
%sIceLoad: sea-ice loading [in km/m^2]
%PHIBOT:   Bottom Pressure Anomaly [in m^2/s^2]
%
%note:      ETAN corresponds to the ocean-ice interface when there is seaice. The 
%           height contribution of seaice+snow (sIceLoad/parms.rhoconst) therefore
%           needs to be accounted for to allow comparison with altimetry.
%note:      bp is converted to equivalent sea level meters by dividing by g.

if isdir(dirNctiles);
ETAN=read_nctiles([dirNctiles 'ETAN/ETAN'],'ETAN');
sIceLoad=read_nctiles([dirNctiles 'sIceLoad/sIceLoad'],'sIceLoad');
PHIBOT=read_nctiles([dirNctiles 'PHIBOT/PHIBOT'],'PHIBOT');

ssh=ETAN+sIceLoad/parms.rhoconst;
bp=PHIBOT/9.81;
end;

%%4) near surface velocity and related vector variables:

%UVELMASS: first component of horizontal velocity [in m/s]
%VVELMASS: second component of horizontal velocity [in m/s]
%
%note: the loaded vector components are along grid line directions. To
%      allow for comparison with e.g. current meters they need to be 
%      converted to zonal and meridional components using calc_UEVNfromUXVY.   

UVELMASS=read_nctiles([dirNctiles 'UVELMASS/UVELMASS'],'UVELMASS',[],1);
VVELMASS=read_nctiles([dirNctiles 'VVELMASS/VVELMASS'],'VVELMASS',[],1);
oceTAUX=read_nctiles([dirNctiles 'oceTAUX/oceTAUX'],'oceTAUX');
oceTAUY=read_nctiles([dirNctiles 'oceTAUY/oceTAUY'],'oceTAUY');

[Ue,Vn]=calc_UEVNfromUXVY(UVELMASS,VVELMASS);
[Te,Tn]=calc_UEVNfromUXVY(oceTAUX,oceTAUY);