% Compute mean and standard deviation of December 2011 surface Theta and KPPhbl.
% Extract surface Theta and KPPhbl at Argo locations during December 2011.

% initialize
clear all, close all
cd ~dmenemen/llc_4320/regions/Argo

% load Argo and keep only December 2011 profiles
load mldinfo_all
it=find(profiledate>datenum(2011,12,1)&profiledate<datenum(2012,1,1));
da_mld       =da_mld       (it);
da_mlpd      =da_mlpd      (it);
da_mls       =da_mls       (it);
da_mlt       =da_mlt       (it);
dt_mld       =dt_mld       (it);
dt_mlpd      =dt_mlpd      (it);
dt_mls       =dt_mls       (it);
dt_mlt       =dt_mlt       (it);
floatnumber  =floatnumber  (it);
names        =names        (it);
profiledate  =profiledate  (it);
profilelat   =profilelat   (it);
profilelon   =profilelon   (it);
profilenumber=profilenumber(it);
status       =status       (it);
ta_mld       =ta_mld       (it);
tt_mld       =tt_mld       (it);

% get llc4320 grid
nx=4320;
prec='real*4';
gdir='~dmenemen/llc_4320/grid/';
YC=quikread_llc([gdir 'YC.data'],nx,1,prec); % latitude from -90 to 90
XC=quikread_llc([gdir 'XC.data'],nx,1,prec); % longitude from -180 to 180

% iterate over all model time steps to extract fields
N=0;
SST=0*it;
HBL=0*it;
meanSST=zeros(nx,nx*13,'single');
mnsqSST=zeros(nx,nx*13,'single');
meanHBL=zeros(nx,nx*13,'single');
mnsqHBL=zeros(nx,nx*13,'single');
pin='~dmenemen/llc_4320/MITgcm/run/';
for ts=283536:144:390528
    N=N+1;

    % extract Argo profiles indices
    tme=datenum(ts2dte(ts,25,2011,9,10));
    it=find(profiledate>(tme-1/24)&profiledate<=tme);
    ix=it;
    for i=1:length(it)
        [Y I]=min((XC(:)-profilelon(it(i))).^2+(YC(:)-profilelat(it(i))).^2);
        ix(i)=I;
    end

    % extract mean, mean squared, and Argo profile SST
    fin=[pin 'Theta.' myint2str(ts,10) '.data'];
    disp(fin)
    fld=quikread_llc(fin,nx,1,prec);
    meanSST=meanSST+fld.^2;
    mnsqSST=mnsqSST+fld.^2;
    if ~isempty(ix)
        SST(it)=fld(ix);
    end

    % extract mean, mean squared, and Argo profile SST
    fin=[pin 'KPPhbl.' myint2str(ts,10) '.data'];
    disp(fin)
    fld=quikread_llc(fin,nx,1,prec);
    meanHBL=meanHBL+fld.^2;
    mnsqHBL=mnsqHBL+fld.^2;
    if ~isempty(ix)
        HBL(it)=fld(ix);
    end
end

meanSST=meanSST/N;
mnsqSST=mnsqSST/N;
meanHBL=meanHBL/N;
mnsqHBL=mnsqHBL/N;

clear I N Y a* fi* g* i* nx pi* pre* tm* ts
    
save December2011