% Replace line below with location of your root directory
cd ~dmenemen/llc_4320/regions/Amundsen2

% Size of domain
nx=1031;
ny=916;
nz=83;

% Read horizontal grid information.
% The grid orientation is not fully aligned with latitude and
% longitude lines so you need 2D arrays to define the grid.
% (XC,YC) is longitude/latitude coordinates of centers of grid cells
% (XG,YG) is longitude/latitude coordinates of southwest corners
XC=readbin('grid/XC_1031x916',[nx ny]);
YC=readbin('grid/YC_1031x916',[nx ny]);
XG=readbin('grid/XG_1031x916',[nx ny]);
YG=readbin('grid/YG_1031x916',[nx ny]);

% Read verticla grid information.
load grid/thk90
dpt=dpt90(1:nz);   % depth of vertical level in m
thk=thk90(1:nz);   % thickness of vertical level in m
clear *90

% Read potential temperature.
% The model output is hourly, one file per hour.
% The first 10 digits of filename denote the model
% time step (ts), and I give example below how to
% convert time step to a date (dte).
ts=1491840;
dte=ts2dte(ts,25,2011,9,10);
fnm=['Theta/' myint2str(ts,10) '_Theta_13325.1965.1_1031.916.83'];
fld=readbin(fnm,[nx ny nz]);

% Example plotting of level k.
k=30;
clf
mypcolor(fld(:,:,k)');
colorbar
title(['Theta (deg C) at ' int2str(dpt(10)) ' m depth on ' dte])