% um_read_nc
% reads a netCDF file
% inputs: file location, UM pole lat and lon 
% These will be used in the call to um_xy_nc.m which
% computes the actual lat/lon (x_p, y_p); given rotated lat/lon (x,y)
% 
% Input variable names are STASH dependent, so will need checked!
% Output variable names are defined here.
% Variables generally q for 4D variable, or q_sfc at a specified level
% Field values are stored as 4D arrays e.g. q(time, levels, y, x)
%
% written by Ian Renfrew Aug 2006
% based on read_nc by Jeff Chagnon, 2006

% ***** set file name -----------------------------------------
dira = 'C:\data\um\';
jobname = 'xcbmh';
filename = [jobname '.nc0'];
infile = [dira jobname '\' filename];
disp(['reading ' infile]);
cd([dira jobname]);

% ***** set UM grid details ------------------------------------------------
polelat = 37.5; polelon = 177.5 ;      % use for 12 km standard UK grid
nc = netcdf(infile);                   % read file into stucture nc

% get dimensions  --- Note: these naming conventions (i.e., the arguments)
% may not be appropriate for every converted pp file
% (x,y) are equatorial lat,lon; x = [no_columns,1], y=[no_rows, 1]
% (x_p,y_p) are rotated lat,lon; xp =[no_columns, no_rows]
x = nc{'x'};                %-360-.5*(180-polelon); degrees east
y = nc{'y'};                %+90-polelat; degrees north
x_1 = nc{'x_1'};            % second horizontal grid
y_1 = nc{'y_1'};            % second horizontal grid
um_xy_nc;                   % need polelat and polelon to convert from eq coord to rotpol coord

% check that all data have been interpolated onto the same (A) grid
if (max(abs(x(:,:)-x_1(:,:))) > 1e-4) disp('UM grids require interpolation - use um2nc.tcl'); end;
if (max(abs(y(:,:)-y_1(:,:))) > 1e-4) disp('UM grids require interpolation - use um2nc.tcl'); end;

% dump metadata to the screen to find out what variables are available
nc_dump(infile);

% height dimensions
hheight_rho = nc{'hybrid_ht'};        % Model hybrid heights levels (Mrho)
hheight_theta = nc{'hybrid_ht_1'};    % Model hybrid height levels (Mtheta)
theta_levels = nc{'theta_1'};         % Constant theta levels (K)
p_levels = nc{'p'};                   % Constant pressure levels (mbar)
z_rho = nc{'ht_1'};                   % Geometric height of Mrho levels (3D)
z_theta = nc{'ht_2'};                 % Geometric height of Mtheta levels (3D)

% time dimensions: units are days since start time
t = nc{'t'};             
t_1 = nc{'t_1'};   
t_2 = nc{'t_2'}; 
t_3 = nc{'t_3'}; 

% atmospheric variables: model levels, all times
u = nc{'x-wind'};           % u on Mrho z levels m/s
v = nc{'y-wind'};           % v on Mrho z levels m/s
theta = nc{'theta'};        % pot temp on Mtheta z_1 levels K
q = nc{'q'};                % specific humidity on Mtheta z_1 levels, kg/kg
w = nc{'dz_dt'};            % w on Mtheta levels m/s
pv = nc{'pv_1'};            % potential vorticity on Mrho levels K m^2 s^-1 kg^-1

% atmospheric variables: model levels, not available at t=0
qcl = nc{'QCL'};            % cloud liquid water, z_1 levels, kg/kg
qcf = nc{'QCF'};            % cloud ice content, z_1 levels, kg/kg

% atmospheric variables: theta levels
pv_theta = nc{'pv'};        % potential vorticity K m^2 s^-1 kg^-1

% atmospheric variables: pressure levels
gph = nc{'ht_3'};           % geopotential height m

% single level variables: constant with time
lsm = nc{'lsm'};            % land surface mask [1=land] at surface
orography = nc{'ht'};       % orography, surface

% single level variables: all times
temp_sfc = nc{'temp'};      % temperature at surface, K 
mslp = nc{'p_1'};           % mean sea level pressure, mbar

% single level variables: not available at t=0
u_10m = nc{'x-wind_1'};     % u at 10 m, m/s
v_10m = nc{'y-wind_1'};     % v at 10 m, m/s
temp_1pt5m = nc{'temp_1'};  % temperature at 1.5 m
q_1pt5m = nc{'q_1'};        % specific humidity at 1.5 m
low_cloud = nc{'field33'};  % low cloud amount (0 to 1)
med_cloud = nc{'field32'};  % medium cloud amount (0 to 1)
high_cloud = nc{'field31'}; % high cloud amount (0 to 1)

% single level variables: means and accumulations
lsrain = nc{'lsrain'};      % large scale rain, surface, kg m^-2
lssnow = nc{'lssnow'};      % large scale snow, surface, kg m^-2
cvrain = nc{'cvrain'};      % convective rain, surface, kg m^-2
cvsnow = nc{'cvsnow'};      % convective snow, surface, kg m^-2
shf_sfc = nc{'sh'};         % surface sensible heat flux
lhf_sfc = nc{'lh'};         % surface latent heat flux
taux_sfc = nc{'taux'};      % x-comp of mean sea surface stress, N/m^2
tauy_sfc = nc{'tauy'};      % y-comp of mean sea surface stress, N/m^2

% derived diagnostics
% wind_speed = sqrt(u(:,:,:,:).*u(:,:,:,:) + v.*v);
tau_sfc = sqrt(taux_sfc(:,:,:,:).*taux_sfc(:,:,:,:) + tauy_sfc.*tauy_sfc);
wind_speed_10m = sqrt(u_10m(:,:,:,:).*u_10m(:,:,:,:) + v_10m.*v_10m);