%atn 9.feb.2021
%example to look at Argo profiles in the Nordic Seas in 2015-2016.
%First, the directory profiles/ must be downloaded, and assumed to store in dir_profile

clear all;
warning off

load coast;	%coast line [long,lat], default of matlab, if does not exist, can comment out this line

addpath('./ASTE_mlibrary/');	%<-- change absolute path as needed
addpath('./ASTE_mlibrary/MITprof_toolbox/');

dir_profile='./';		%<-- change absolute path as needed

%example: pick two years of Argo to show how to concat files, 2014,2015
fileIn={'argo_feb2016_2014_llc270_18x18_ASTE_model.nc',...
        'argo_feb2017_2015_llc270_18x18_ASTE_model.nc'};

a0=MITprof_load([dir_profile fileIn{1}]);			%2014
a =MITprof_concat(a0,MITprof_load([dir_profile fileIn{2}]));	%combined 2014 and 2015
%a0 = 						a = 
%  struct with fields:                             struct with fields:
%                                                
%          prof_YYYYMMDD: [34389×1 double]                 prof_YYYYMMDD: [72781×1 double]
%            prof_HHMMSS: [34389×1 double]                   prof_HHMMSS: [72781×1 double]
%               prof_lon: [34389×1 double]                      prof_lon: [72781×1 double]
%               prof_lat: [34389×1 double]                      prof_lat: [72781×1 double]
%             prof_depth: [55×1 double]                       prof_depth: [55×1 double]
%             prof_descr: {34389×1 cell}                      prof_descr: {72781×1 cell}
%              prof_date: [34389×1 double]                     prof_date: [72781×1 double]
%             prof_basin: [34389×1 double]                    prof_basin: [72781×1 double]
%             prof_point: [34389×1 double]                    prof_point: [72781×1 double]
%                 prof_T: [34389×55 double]                       prof_T: [72781×55 double]
%           prof_Tweight: [34389×55 double]                 prof_Tweight: [72781×55 double]
%            prof_Testim: [34389×55 double]                  prof_Testim: [72781×55 double]
%              prof_Terr: [34389×55 double]                    prof_Terr: [72781×55 double]
%             prof_Tflag: [34389×55 double]                   prof_Tflag: [72781×55 double]
%                 prof_S: [34389×55 double]                       prof_S: [72781×55 double]
%           prof_Sweight: [34389×55 double]                 prof_Sweight: [72781×55 double]
%            prof_Sestim: [34389×55 double]                  prof_Sestim: [72781×55 double]
%              prof_Serr: [34389×55 double]                    prof_Serr: [72781×55 double]
%             prof_Sflag: [34389×55 double]                   prof_Sflag: [72781×55 double]
%       prof_interp_XC11: [34389×1 double]              prof_interp_XC11: [72781×1 double]
%       prof_interp_YC11: [34389×1 double]              prof_interp_YC11: [72781×1 double]
%     prof_interp_XCNINJ: [34389×1 double]            prof_interp_XCNINJ: [72781×1 double]
%     prof_interp_YCNINJ: [34389×1 double]            prof_interp_YCNINJ: [72781×1 double]
%          prof_interp_i: [34389×1 double]                 prof_interp_i: [72781×1 double]
%          prof_interp_j: [34389×1 double]                 prof_interp_j: [72781×1 double]
%        prof_interp_lon: [34389×1 double]               prof_interp_lon: [72781×1 double]
%        prof_interp_lat: [34389×1 double]               prof_interp_lat: [72781×1 double]
%    prof_interp_weights: [34389×1 double]           prof_interp_weights: [72781×1 double]
%             list_descr: {1189×1 cell}                       list_descr: {1489×1 cell}
%                     nd: 1189                                        nd: 1489
%                     np: 34389                                       np: 72781
%                     nr: 55                                          nr: 55                

%now compute (Tmodel minus Tobs) / sigma_T = (Tmodel minus Tobs) * sqrt(weightT)
%model fields are a.prof_[S,T]estim
%obs   fields are a.prof_[S,T]
%weights = 1/(sigma^2) are stored in a.prof_weight[T,S]

%izero=find(a.prof_weightT(:)==0);
%sigmaT=1./sqrt(a.prof_weightT);sigmaT(izero)=nan;	%make nan to avoid div by zero

misfitT=(a.prof_Testim-a.prof_T).*sqrt(a.prof_Tweight);
misfitS=(a.prof_Sestim-a.prof_S).*sqrt(a.prof_Sweight);

figure(1);clf;
plot(a.prof_lon,a.prof_lat,'.');grid;
if(exist('long','var'));hold on;plot(long,lat,'k-');hold off;end;

%limit to only within Nordic Seas, here APPROXIMATE by convenient lat/lon limits:
ij=find(a.prof_lat>64 & a.prof_lon>-20);

%pick a depth, based on a.prof_depth' =
%[5 15 25 35 45 55 65 75 85 95 105 115 125 135 145 155 165 175 185 200 220 240 260 280 300 320 340 360 380 400 
%420 440 460 480 500 550 600 650 700 750 800 850 900 950 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000]
iz=find(a.prof_depth==105); %pick 105m

figure(2);clf;colormap(jet(21));
subplot(121);
hh=scatter3(a.prof_lon(ij),a.prof_lat(ij),misfitT(ij,iz),10.*ones(length(ij),1),misfitT(ij,iz),'filled');
  grid;view(0,90);colorbar('southoutside');xlabel('Lon [degE]');ylabel('Lat [degN]');
  if(exist('long','var'));hold on;plot(long,lat,'k-');hold off;end;axis([-20 20 64 80]);
  title(['(T_{ASTEr1}-T_{Argo})/\sigma{T} (unitless), z=' num2str(a.prof_depth(iz)) 'm, 2014--2015'],'interpreter','tex');

subplot(122);
hh=scatter3(a.prof_lon(ij),a.prof_lat(ij),misfitS(ij,iz),10.*ones(length(ij),1),misfitS(ij,iz),'filled');
  grid;view(0,90);colorbar('southoutside');xlabel('Lon [degE]');ylabel('Lat [degN]');
  if(exist('long','var'));hold on;plot(long,lat,'k-');hold off;end;axis([-20 20 64 80]);
  title(['(S_{ASTEr1}-S_{Argo})/\sigma{S} (unitless), z=' num2str(a.prof_depth(iz)) 'm, 2014--2015'],'interpreter','tex');

figure(2);set(gcf,'paperunits','inches','paperposition',[0 0 12 7]);
fpr=['./ASTER1_Argo_NordicSeas_20142015_misfitTS.png'];print(fpr,'-dpng');
