%%========= Sample code to read the MIZMAS model results =============
%%========== Written by Wenli Zhong, OUC-UW, May 29, 2018 =============
%%=======================================================================
%% ================= Load the Model grids ====================
nx=168; ny=180; nz=40;  nt=365;  % longitude: nx=168;  latitude: ny=180;  level: nz=40;  Year days: nt=365;

gridname='grid.dat.rot_168_180';
depthname='dz.dta40';
maskname='levels_40_t_168_180';

grid=load(gridname);
lontt=zeros([nx*ny 1]);  lattt=zeros([nx*ny 1]);
fg=0;
for i=1:size(grid,1)/2
    for j=1:size(grid,2)
        fg=fg+1;
        lontt(fg,1)=grid(i,j);
        lattt(fg,1)=grid(i+size(grid,1)/2,j);
    end
end
mlon=reshape(lontt,[nx,ny]);  mlat=reshape(lattt,[nx,ny]); % convert the griddata to model longitude and latitude

dz=load(depthname);
dep=cumsum(dz/100);   % The real depth of the each model layer 
%% ===============================================================

for year=1992:2017
    
    %% Ice concentration data
    icefilename=strcat('aiday_168_180.H',num2str(year));
    fid1 = fopen(icefilename,'r');
    preidata = fread(fid1,'float32');       % Read the ice concentration data
    idata = reshape(preidata,[nx ny nt]);   % fraction of ice concerntration
    fclose(fid1);
    clear preidata
    
    %% Surface wind field
    windfilename=strcat('sur_wind_168_180.',num2str(year));
    fid2 = fopen(windfilename,'r');
    prewind = fread(fid2,[nx*ny*2 nt],'float32');   % Read the surface wind field data
    preUw=prewind(1:size(prewind,1)/2,:);                 % decomposed the velocity into u and v
    preVw=prewind(size(prewind,1)/2+1:size(prewind,1),:);
    windUU=reshape(preUw,[nx,ny,nt]);        % don't write it this way: UU=reshape(uv(1:size(uv,1)/2,:),[nx,ny,nz,nt]); because it's more time consuming!!!
    windVV=reshape(preVw,[nx,ny,nt]);
    fclose(fid2);
    clear prewind preUw preVw
    
    %% Surface geostrophic current from SSH field
    sshugfilename=strcat('ugday_168_180.H',num2str(year));
    fid3 = fopen(sshugfilename,'r');
    preUVg = fread(fid3,[nx*ny*2 nt],'float32');   % Read the Surface geostrophic current
    preUg=preUVg(1:size(preUVg,1)/2,:);                 % decomposed the velocity into u and v
    preVg=preUVg(size(preUVg,1)/2+1:size(preUVg,1),:);
    UUg_ssh = reshape(preUg,[nx ny nt]);
    VVg_ssh = reshape(preVg,[nx ny nt]);
    fclose(fid3);
    clear preUVg preUg preVg
    
    %% Sea ice motion vectors
    icevecfilename=strcat('uiday_168_180.H',num2str(year));
    fid4 = fopen(icevecfilename,'r');
    preicevec = fread(fid4,[nx*ny*2 nt],'float32');   % Read the Sea ice motion vectors
    preUi=preicevec(1:size(preicevec,1)/2,:);                 % decomposed the velocity into u and v
    preVi=preicevec(size(preicevec,1)/2+1:size(preicevec,1),:);
    iceUU = reshape(preUi,[nx ny nt]);
    iceVV = reshape(preVi,[nx ny nt]);
    fclose(fid4);
    clear preicevec preUi preVi
    
    %% Ocean temperature
    tofilename=strcat('to_168_180.H',num2str(year));    
    fid5 = fopen(tofilename,'r');     % Read the ocean temperature data
    to = fread(fid5,'float32');
    Temp = reshape(to,[nx,ny,nz,nt]);
    fclose(fid5);
    
    %% Ocean Salinity
    sofilename=strcat('so_168_180.H',num2str(year));
    fid6 = fopen(sofilename,'r');     % Read the ocean salinity data
    so = fread(fid6,'float32');
    Salt = reshape(so,[nx,ny,nz,nt]);
    fclose(fid6);
    
    %% Ocean current (U,V)
    uvfilename=strcat('uo_168_180.H',num2str(year));  
    fid = fopen(uvfilename,'r');   % Read the ocean current data
    uv = fread(fid,[nx*ny*nz*2 nt],'float32');
    fclose(fid);
    
    preU=uv(1:size(uv,1)/2,:);             % decomposed the velocity into u and v
    preV=uv(size(uv,1)/2+1:size(uv,1),:);
    UU=reshape(preU,[nx ny nz nt]);        % don't write it this way: UU=reshape(uv(1:size(uv,1)/2,:),[nx,ny,nz,nt]); because it's more time consuming!!!
    VV=reshape(preV,[nx ny nz nt]);
    clear preU preV uv
    
end





