% function sample_nc_diagnostics



clear all; close all; clc % <-- clears memory, figures and command line



% specify files to open

cy6g_mif = '/project/p1382/svinogra/datasets/ncom-letkf-cycling-datasets/mif_0735120_analyzed_001_cy6g.nc';

cy6g_mof = '/project/p1382/svinogra/datasets/ncom-letkf-cycling-datasets/mof_0735120_forecast_001_cy6g.nc';



% open/read files from netcdf

cy6g_mif_nc = netcdf(cy6g_mif,'nowrite'); % change to netcdf.open!

cy6g_mof_nc = netcdf(cy6g_mof,'nowrite');



% reading grid and creating mesh for 2d plots

lon = cy6g_mif_nc{'lon',1}(:);  % deg East

lat = cy6g_mif_nc{'lat',1}(:);  % deg North

depth = cy6g_mif_nc{'depth',1}(:);  % m, down

[mesh.lon, mesh.lat] = meshgrid(lon, lat);

[IM,JM] = size(mesh.lon');



%specifying a vertical slice

xslice = [120,130]; % from 120E to 130E  

yslice = [26.3, 26.3270]; % getting just one grid line in latitude

[jdx, idx] = find ( ...

    xslice(1) <= mesh.lon & mesh.lon <= xslice(2) & ...

    yslice(1) <= mesh.lat & mesh.lat <= yslice(2));

jdx = jdx(1);

[slice.x, slice.z] = meshgrid(mesh.lon(1,idx), depth); % meshing the slice



% time

tim = cy6g_mif_nc{'time',1}(:); % hours since 2004.01.02 06:00 UTC

tau = cy6g_mif_nc{'tau',1}(:);  % hours since analysis



% fields cy6g

cy6g_mif_surf_el = cy6g_mif_nc{'surf_el',1}(:,:);  % extracting SSH

cy6g_mof_surf_el = cy6g_mof_nc{'surf_el',1}(:,:); 

cy6g_inc_el = cy6g_mif_surf_el - cy6g_mof_surf_el; % increments 



cy6g_mif_water_ptem = cy6g_mif_nc{'water_ptem',1}(:,:,:); % extracting temperature

cy6g_mof_water_ptem = cy6g_mof_nc{'water_ptem',1}(:,:,:); 

cy6g_inc_ptem = squeeze(cy6g_mif_water_ptem(1,:,:)) - squeeze(cy6g_mof_water_ptem(1,:,:)); % increments

slice_inc_ptem = squeeze(cy6g_mif_water_ptem(:,jdx,idx)) - squeeze(cy6g_mof_water_ptem(:,jdx,idx)); % slice of increments



cy6g_mif_water_u = cy6g_mif_nc{'water_u',1}(:,:,:); % extracting velocities

cy6g_mof_water_u = cy6g_mof_nc{'water_u',1}(:,:,:); 

cy6g_mif_water_v = cy6g_mif_nc{'water_v',1}(:,:,:); 

cy6g_mof_water_v = cy6g_mof_nc{'water_v',1}(:,:,:); 

cy6g_inc_vel = sqrt(squeeze(cy6g_mif_water_u(1,:,:)).^2 +squeeze(cy6g_mif_water_v(1,:,:)).^2) - ...

    sqrt(squeeze(cy6g_mof_water_u(1,:,:)).^2 +squeeze(cy6g_mof_water_v(1,:,:)).^2); % magnitudes of increments

slice_inc_uv = sqrt(squeeze(cy6g_mif_water_u(:,jdx,idx)).^2 +squeeze(cy6g_mif_water_v(:,jdx,idx)).^2) - ...

    sqrt(squeeze(cy6g_mof_water_u(:,jdx,idx)).^2 +squeeze(cy6g_mof_water_v(:,jdx,idx)).^2); % slics of magnitudes of increments





data = read_dof('./dof_0735120_analyzed_001_cy6g.txt'); % <-- reads dof 



% plotting the fields

latlim = [16 35]; 

clim = 'auto';



% % % SSH mif

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, cy6g_mif_surf_el); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mif surf el');

% 

% % % SSH mof

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, cy6g_mof_surf_el); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mof surf el');



% % SSH increments

cy6g_inc_el(find(cy6g_inc_el==0))=NaN; % set 0 values to NaN

figure('Position',[200 200 800 700]); hold on

pcolor(mesh.lon, mesh.lat, cy6g_inc_el); 

set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

title('cy6g mif SSH - cy6g mof SSH');

plot(data.lon,data.lat,'k.'); % <-- plots data locations

caxis([-2 2]); % <-- color scale, do caxis 'auto'; for automatic

text(118.2, 16.5, ['min=' num2str( min(min(cy6g_inc_el)) ) ', max= ' num2str( max(max(cy6g_inc_el)) )]);



% % % SST

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, squeeze(cy6g_mif_water_ptem(1,:,:))); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mif SST');

% 

% % % SST

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, squeeze(cy6g_mof_water_ptem(1,:,:))); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mof SST');



% % SST

figure('Position',[200 200 800 700]); hold on

pcolor(mesh.lon, mesh.lat, cy6g_inc_ptem); 

set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

title('cy6g mif SST - cy6g mof SST');

plot(data.lon,data.lat,'k.');

caxis([-2 2])

plot(mesh.lon(jdx,idx), mesh.lat(jdx,idx),'r.');

text(118.2, 16.5, ['min=' num2str( min(min(cy6g_inc_ptem)) ) ', max= ' num2str( max(max(cy6g_inc_ptem)) )]);



% 

figure('Position',[100 100 1500 300]); hold on;

pcolor(slice.x, -slice.z, slice_inc_ptem); shading flat; colorbar; set(gca,'YLim',[-500 0]);

[C,h] = contour(slice.x, -slice.z, slice_inc_ptem,[-10:1:10]); set(h,'LineColor','k');

xlabel('LONGITUDE, ^oE'); ylabel('DEPTH, M'); title(['cy6g mif - cy6g mof: ptem @ LAT = ' num2str(yslice(1)) '^oN']);

caxis([-8.5 6.1]);

caxis([-2 2]);

text(120.2, -450, ['CL=1, min=' num2str( min(min(slice_inc_ptem)) ) ', max= ' num2str( max(max(slice_inc_ptem)) )]);



% % % SSUV

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, sqrt(squeeze(cy6g_mif_water_u(1,:,:)).^2 +squeeze(cy6g_mif_water_v(1,:,:)).^2) ); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mif SSUV');

% 

% % % SSUV

% figure('Position',[200 200 1800/3 350]); hold on

% pcolor(mesh.lon, mesh.lat, sqrt(squeeze(cy6g_mof_water_u(1,:,:)).^2 +squeeze(cy6g_mof_water_v(1,:,:)).^2) ); 

% set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

% title('cy6g mof SSUV');



% % SSUV

figure('Position',[200 200 800 700]); hold on

pcolor(mesh.lon, mesh.lat, cy6g_inc_vel); 

set(gca,'YLim',latlim); shading flat; caxis(clim); colorbar;  

title('cy6g mif SSUV - cy6g mof SSUV');

plot(data.lon,data.lat,'k.');

plot(mesh.lon(jdx,idx), mesh.lat(jdx,idx),'r.');

caxis([-0.4 0.4]);

text(118.2, 16.5, ['min=' num2str( min(min(cy6g_inc_vel)) ) ', max= ' num2str( max(max(cy6g_inc_vel)) )]);



% % slice 

figure('Position',[100 100 1500 300]); hold on;

pcolor(slice.x, -slice.z, slice_inc_uv); shading flat; colorbar; set(gca,'YLim',[-500 0]);

[C,h] = contour(slice.x, -slice.z, slice_inc_uv,[-10:1:10]); set(h,'LineColor','k');

xlabel('LONGITUDE, ^oE'); ylabel('DEPTH, M'); title(['cy6g mif - cy6g mof: uv @ LAT = ' num2str(yslice(1)) '^oN']);

caxis([-0.55 1.2]);

caxis([-0.2 0.2]);

text(120.2, -450, ['CL=1, min=' num2str( min(min(slice_inc_uv)) ) ', max= ' num2str( max(max(slice_inc_uv)) )]);





