Matlab plot -- vertical temperature distribution

clc;clear;close;

yr=2019;
start_date='2019-01-01 00:00:00';

%% obs vertical temperature distribution
% load observation data
ncload('/Users/chuyanzhao/Desktop/gls_observation_data/obs_glerl_mooring/glerl_lake_superior_temperature_mooring_2018_2020.nc');

% adjust time to target time period
tmp_time_obs=datetime(datestr(time/86400+datenum(1970,1,1)));
time_obs_start=datenum(start_date-tmp_time_obs(1))*24;
time_obs_end=time_obs_start+365*24;
tmp_year=time(time_obs_start:time_obs_end);
time_obs=(tmp_year-tmp_year(1))/3600/24;

% read observation depth
depth_obs=-z;

% read observation temperature
temp_obs=sea_water_temperature(time_obs_start:time_obs_end,:);
tmp_ind=temp_obs<-1000;
temp_obs(tmp_ind)=NaN;
lat_obs=[lat;lat];
lon_obs=[lon;lon];

% set colorbar
cmin=0;cinc=1;cmax=20;
ncc=round([cmax-cmin]/cinc);
ccmap=flipud(spectral(ncc));

% plot obs vertical temperature
figure(101);clf;hold on;box on;
    set(gca,'tickdir','out');
    contourf(time_obs,depth_obs,temp_obs',[cmin:cinc:cmax],'LineColor','none');
    caxis([cmin cmax]);
    colormap(ccmap);
    colorbar('location','eastoutside');
    axis([0 366 -190 -5]);
    lcb=colorbar('location','eastoutside');
    ylabel(lcb,'LST (\circC)');
    text(10,-180,['Vertical temp profile by mooring'],'background','w');
    xlabel(['Julian day ' num2str(yr)]);
    ylabel(['Depth (m)'])

    saveas(gcf,['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/figure/HighRes/fig_vertical_temp_obs_' num2str(yr) '.png'],'png'); 

%% find the nearest node in fvcom to the mooring location

% read fvcom grids
ncload(['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/raw_data/HighRes/gls_lst_2019_HighRes.nc'],'lon','lat');
lat_mod=lat;
lon_mod=lon-360;

% find the nearest node of the mooring location
xy_node=[lat_mod lon_mod];
xy_obs=[lat_obs lon_obs];

loc=find_nearest_group(xy_node,xy_obs);

pos=loc(1,1)

%% fvcom vertical temperature distribution

% generate the raw data by linux command and load
ncload(['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/raw_data/HighRes/gls_ver_temp_2019_HighRes.nc']);
temp_mod=temp;

[nt,nn]=size(temp_mod);

tmp_time_mod=[1:nt];
time_mod=tmp_time_mod/24;
depth_mod=h.*siglay(:);

cmin=0;cinc=1;cmax=20.0;
ncc=round([cmax-cmin]/cinc);
ccmap=flipud(spectral(ncc));

figure(201);clf;hold on;box on;
    set(gca,'tickdir','out');
    contourf(time_mod,depth_mod,temp_mod',[cmin:cinc:cmax],'LineColor','none');
    caxis([cmin cmax]);
    colormap(ccmap);
    axis([0 366 -190 -5]);
    lcb=colorbar('location','eastoutside');
    ylabel(lcb,'LST (\circC)');
    text(10,-180,['Vertical temp profile by fvcom'],'background','w');
    xlabel(['Julian day ' num2str(yr)]);
    ylabel(['Depth (m)'])

    saveas(gcf,['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/figure/HighRes/fig_vertical_temp_mod_' num2str(yr) '.png'],'png'); 

%%  Difference between fvcom and observation -- interpolate fvcom temp to mooring location
temp_mod2obs= interp1(depth_mod,temp_mod',depth_obs);

% calculate the difference
temp_diff=temp_mod2obs-temp_obs';

% set colorbar
color_scale1=[
178,24,43
214,96,77
244,165,130
253,219,199
247,247,247
247,247,247
209,229,240
146,197,222
67,147,195
33,102,172]/255;

cmin=-5;cinc=2;cmax=5.0;
ccmap=flipud(color_scale1);

% plot
figure(301);clf;hold on;box on;
    set(gca,'tickdir','out');
    contourf(time_obs,depth_obs,temp_diff,[cmin:cinc:cmax],'LineColor','none');
    caxis([cmin cmax]);
    colormap(ccmap);
    colorbar('location','eastoutside');
    axis([0 366 -180 -5]);
    lcb=colorbar('location','eastoutside');
    ylabel(lcb,'LST (\circC)');
    text(10,-170,['Vertical temp difference profile (fvcom2obs)'],'background','w');
    xlabel(['Julian day ' num2str(yr)]);
    ylabel(['Depth (m)'])

saveas(gcf,['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/figure/fig_vertical_temp_diff_fvcom2obs' num2str(yr) '.png'],'png'); 

%%  Difference between fvcom and observation -- interpolate mooring temp to fvcom layer
temp_obs2mod= interp1(depth_obs,temp_obs',depth_mod);

% calculate the difference
temp_diff=temp_mod'-temp_obs2mod;

% set colorbar
cmin=-5;cinc=2;cmax=5.0;
ccmap=flipud(color_scale1);

% plot
figure(401);clf;hold on;box on;
    set(gca,'tickdir','out');
    contourf(time_mod,depth_mod,temp_diff,[cmin:cinc:cmax],'LineColor','none');
    caxis([cmin cmax]);
    colormap(ccmap);
    colorbar('location','eastoutside');
    axis([0 366 -180 -5]);
    lcb=colorbar('location','eastoutside');
    ylabel(lcb,'LST (\circC)');
    text(10,-170,['Vertical temp difference profile (obs2fvcom)'],'background','w');
    xlabel(['Julian day ' num2str(yr)]);
    ylabel(['Depth (m)'])

saveas(gcf,['/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/2019/figure/fig_vertical_temp_mod_obs2fvcom' num2str(yr) '.png'],'png');