%% read data
fid = 'JRA55/anl_p25.007_hgt.197901_197912';
data = ncdataset(fid);
lon0  = squeeze(data.data('lon'));
lat  = squeeze(data.data('lat'));
level = squeeze(data.data('isobaric'));
timemid = squeeze(data.data('time'));
hgtmid0 = double(squeeze(data.data('Geopotential_height_isobaric_Mixed_intervals_Average')));
fid = 'JRA55/anl_p25.011_tmp.197901_197912';
data = ncdataset(fid);
airmid0 = squeeze(data.data('Temperature_isobaric_Mixed_intervals_Average'));
air = zeros(144,73,37,528);
hgt = zeros(144,73,37,528);
time = zeros(528,1);
t = 1;
airmid = permute(airmid0,[4 3 2 1]);
hgtmid = permute(hgtmid0,[4 3 2 1]);
air(:,:,:,1+(t-1)*12:t*12) = airmid;
hgt(:,:,:,1+(t-1)*12:t*12) = hgtmid;
time(1+(t-1)*12:t*12) = timemid;
for t = 2:44
    fid = ['JRA55/anl_p25.007_hgt.',num2str(1978+t),'01_',num2str(1978+t),'12'];
    data = ncdataset(fid);
    hgtmid0 = double(squeeze(data.data('Geopotential_height_isobaric_Mixed_intervals_Average')));
    fid = ['JRA55/anl_p25.011_tmp.',num2str(1978+t),'01_',num2str(1978+t),'12'];
    data = ncdataset(fid);
    airmid0 = squeeze(data.data('Temperature_isobaric_Mixed_intervals_Average'));
    timemid = squeeze(data.data('time'));
    airmid = permute(airmid0,[4 3 2 1]);
    hgtmid = permute(hgtmid0,[4 3 2 1]);
    air(:,:,:,1+(t-1)*12:t*12) = airmid;
    hgt(:,:,:,1+(t-1)*12:t*12) = hgtmid;
    time(1+(t-1)*12:t*12) = timemid;
end
save jra55data air hgt

%% ==== n2
clear all 
clc;
fid = 'NCEPdata/land.nc';
ncdisp(fid)
land = ncread(fid,'land');% 1 is land; 0 is ocean
lat = ncread(fid,'lat');
lon0 = ncread(fid,'lon');
fid = 'NCEPdata/hgt.mon.mean.nc';
ncdisp(fid)
hgt = ncread(fid,'hgt');% 
level = ncread(fid,'level');
time = ncread(fid,'time');
[xm,ym,zm,tm] = size(hgt);
lon0(xm/2:xm) =  lon0(xm/2:xm)-360;
% ==== calculate the N2
%% SST\SH\LH 1900-1994 DAYLY
  load('masknor.mat')
% fid = 'NCEPdata/daily/interpolated_dlwrf.sfc.gauss.1979.nc';
% ncdisp(fid)
% data1 = ncread(fid,'dlwrf');
% fid = 'NCEPdata/daily/dlwrf.sfc.gauss.1979.nc';
% data2 = ncread(fid,'dlwrf');
% figure;contourf(data1(:,:,3))
% figure;contourf(data2(:,:,3))
for yr = 1989:1995
    fid = ['NCEPdata/daily/air.',num2str(yr),'.nc'];
    % ncdisp(fid)
    airYR = ncread(fid,'air');
    level = ncread(fid,'level');
    lon0NCEP = double(ncread(fid,'lon'));
    latNCEP = double(ncread(fid,'lat'));
    % fid= ['NCEPdata/daily/air.1988.nc';
    % air1988 = ncread(fid,'air');
    % === read vertical velocity
    fid = ['NCEPdata/daily/omega.',num2str(yr),'.nc'];
    % ncdisp(fid)
    omegaYR = ncread(fid,'omega');
    fid = ['NCEPdata/daily/hgt.',num2str(yr),'.nc'];
    % ncdisp(fid)
    hgtYR = ncread(fid,'hgt');
    [xm,ym,zm,tmYR] = size(hgtYR);
    % ---- POTENTIAL TMEPERATURE
    Rd = 287.058;
    gg0 = 9.8 ;% gravity accaleration;
    gg0R = 1.0/9.8;
    thetaYR = zeros(xm,ym,zm,tmYR);% potential temperature;
    N20YR = zeros(xm,ym,zm,tmYR);% the stability paramenter;
    Thro0YR = zeros(xm,ym,zm,tmYR);
    for k =1:zm
        if k==1
            thetaYR(:,:,k,:) = airYR(:,:,k,:).*(1000./level(k))^0.286;
        end
        if k < zm
            thetaYR(:,:,k+1,:) = airYR(:,:,k+1,:).*(1000./level(k+1))^0.286;
            %             N20(:,:,k,:) = (theta(:,:,k+1,:)-theta(:,:,k,:))./((hgt(:,:,k+1,:)-hgt(:,:,k,:)));
            N20YR(:,:,k,:) = gg0.*(thetaYR(:,:,k+1,:)-thetaYR(:,:,k,:))./...
                (thetaYR(:,:,k,:).*(hgtYR(:,:,k+1,:)-hgtYR(:,:,k,:)));
            Thro0YR(:,:,k,:) = level(k)*100/Rd/airYR(:,:,k,:);
        end
    end
    fnm = ['DAY50/resultN2',num2str(yr),'.mat'];
    save  (fnm, 'N20YR','thetaYR','hgtYR' )
    clear thetaYR hgtYR
    LTN20YR = squeeze(mean(N20YR(:,:,2:3,:),3));%700-850,850-925

    % === read guass data
    % READ SH AND LH DATA 192*94
    fid = ['NCEPdata/daily/interpolated_lhtfl.sfc.gauss.',num2str(yr),'.nc'];
    LHlrNCEPYR = double(ncread(fid,'lhtfl'));%
    %
    fid = ['NCEPdata/daily/interpolated_shtfl.sfc.gauss.',num2str(yr),'.nc'];
    SHlrNCEPYR = double(ncread(fid,'shtfl'));
    %
    fid= ['NCEPdata/daily/interpolated_uswrf.sfc.gauss.',num2str(yr),'.nc'];
    USlrNCEPYR = double(ncread(fid,'uswrf'));%
    %
    fid= ['NCEPdata/daily/interpolated_dswrf.sfc.gauss.',num2str(yr),'.nc'];
    DSlrNCEPYR = double(ncread(fid,'dswrf'));%
    %
    fid= ['NCEPdata/daily/interpolated_ulwrf.sfc.gauss.',num2str(yr),'.nc'];
    ULlrNCEPYR = double(ncread(fid,'ulwrf'));%
    %
    fid= ['NCEPdata/daily/interpolated_dlwrf.sfc.gauss.',num2str(yr),'.nc'];
    DLlrNCEPYR = double(ncread(fid,'dlwrf'));
    % ==== SST
    fid = ['NCEPdata/daily/interpolated_sst.day.mean.',num2str(yr),'.nc'];
    % Timesst(1069,:) is 197901
    sstLrYR = double(ncread(fid,'sst'));
    %% pick up the points and regional average original series
    [xmNCEP,ymNCEP,zmNCEP,tmNCEP] = size(omegaYR);
    LTN20NCEP_AAYR = ones(xmNCEP*7,tmYR).*nan;
    LTN20NCEP_BBYR = ones(xmNCEP*7,tmYR).*nan;
    SST_AAYR = ones(xmNCEP*7, tmYR).*nan;
    SST_BBYR = ones(xmNCEP*7, tmYR).*nan;
    air925_AAYR = ones(xmNCEP*7, tmYR).*nan;
    air925_BBYR = ones(xmNCEP*7, tmYR).*nan;
    air850_AAYR = ones(xmNCEP*7, tmYR).*nan;
    air850_BBYR = ones(xmNCEP*7, tmYR).*nan;
    omega_AAYR = ones(xmNCEP*7,tmYR).*nan;
    omega_BBYR = ones(xmNCEP*7,tmYR).*nan;
    SH_AAYR = ones(xmNCEP*7, tmYR).*nan;
    SH_BBYR = ones(xmNCEP*7, tmYR).*nan;
    LH_AAYR = ones(xmNCEP*7, tmYR).*nan;
    LH_BBYR = ones(xmNCEP*7, tmYR).*nan;
    UL_AAYR = ones(xmNCEP*7, tmYR).*nan;
    DL_AAYR = ones(xmNCEP*7, tmYR).*nan;
    US_AAYR = ones(xmNCEP*7, tmYR).*nan;
    DS_AAYR = ones(xmNCEP*7, tmYR).*nan;
    UL_BBYR = ones(xmNCEP*7, tmYR).*nan;
    DL_BBYR = ones(xmNCEP*7, tmYR).*nan;
    US_BBYR = ones(xmNCEP*7, tmYR).*nan;
    DS_BBYR = ones(xmNCEP*7, tmYR).*nan;
    kk1 = 1;
    kk2 = 1;
    for i = 1:xmNCEP
        for j = 1:25
            % === pick up the points with the trend
            ln1mNCEP = (squeeze(LTN20YR(i,j,1:tmYR)));
            ln2mNCEP = (squeeze(sstLrYR(i,j,1:tmYR)));
            ln3mNCEP = (squeeze(airYR(i,j,2,1:tmYR)));
            ln4mNCEP = (squeeze(airYR(i,j,3,1:tmYR)));
            ln5mNCEP = (squeeze(omegaYR(i,j,6,1:tmYR)));%500hPa
            ln6mNCEP = (squeeze(LHlrNCEPYR(i,j,1:tmYR)));
            ln7mNCEP = (squeeze(SHlrNCEPYR(i,j,1:tmYR)));
            ln8mNCEP = (squeeze(ULlrNCEPYR(i,j,1:tmYR)));
            ln9mNCEP = (squeeze(DLlrNCEPYR(i,j,1:tmYR)));
            ln10mNCEP = (squeeze(USlrNCEPYR(i,j,1:tmYR)));
            ln11mNCEP = (squeeze(DSlrNCEPYR(i,j,1:tmYR)));
            if (maskNor(i,j)==3)%norwegian
                LTN20NCEP_AAYR(kk1,:) = ln1mNCEP;
                SST_AAYR(kk1,:) = ln2mNCEP;
                air925_AAYR(kk1,:)  = ln3mNCEP;
                air850_AAYR(kk1,:)  = ln4mNCEP;
                omega_AAYR(kk1,:)  = ln5mNCEP;
                LH_AAYR(kk1,:)  = ln6mNCEP;
                SH_AAYR(kk1,:)  = ln7mNCEP;
                UL_AAYR(kk1,:) = ln8mNCEP;
                DL_AAYR(kk1,:) = ln9mNCEP;
                US_AAYR(kk1,:) = ln10mNCEP;
                DS_AAYR(kk1,:) = ln11mNCEP;
                kk1 = kk1+1;
            elseif (maskNor(i,j)==1)%iceland
                LTN20NCEP_BBYR(kk2,:) = ln1mNCEP;
                SST_BBYR(kk2,:)  = ln2mNCEP;
                air925_BBYR(kk2,:) = ln3mNCEP;
                air850_BBYR(kk2,:) = ln4mNCEP;
                omega_BBYR(kk2,:) = ln5mNCEP;
                LH_BBYR(kk2,:) = ln6mNCEP;
                SH_BBYR(kk2,:) = ln7mNCEP;
                UL_BBYR(kk2,:) = ln8mNCEP;
                DL_BBYR(kk2,:) = ln9mNCEP;
                US_BBYR(kk2,:) = ln10mNCEP;
                DS_BBYR(kk2,:) = ln11mNCEP;
                kk2 = kk2+1;
            end
        end
    end
    fnm = ['DAY50/resultAABB',num2str(yr),'.mat'];
    save  (fnm, 'LTN20NCEP_AAYR','LTN20NCEP_BBYR' , ...
        'SST_AAYR' ,'SST_BBYR', ...
        'air925_AAYR','air925_BBYR', ...
        'air850_AAYR','air850_BBYR', ...
        'omega_AAYR','omega_BBYR', ...
        'SH_AAYR','SH_BBYR','LH_AAYR'...
        ,'LH_BBYR','UL_AAYR','DL_AAYR', ...
        'US_AAYR','DS_AAYR','UL_BBYR', ...
        'DL_BBYR','US_BBYR','DS_BBYR')
end
%% READ sst and HEAT FLUX
tmERA5 = 516;
% ==== read omega
fid = 'NCEPdata/omega.mon.mean.nc';
ncdisp(fid)
time = ncread(fid,'time');
level = ncread(fid,'level');
omega =  double(ncread(fid,'omega'));
lon0NCEP = double(ncread(fid,'lon'));
latNCEP = double(ncread(fid,'lat'));
[xm,ym,zm,tm] = size(omega);
fid = 'NCEPdata/sst.mon.mean.nc';
ncdisp(fid)
timest = ncread(fid,'time');
Timesst = datestr((datenum('18900101','yyyymmdd')+timest),'yyyymm');
% Timesst(1069,:) is 197901
startLoc = [1 1  1069]; 
count  = [360 180 1590-1068]; 
stride = [1 1 1]; 
sstNCEP = double(ncread(fid,'sst',startLoc,count,stride));
latsst = double(ncread(fid,'lat'));
lon0sst = double(ncread(fid,'lon'));
%  === change grid 
[x,y] = meshgrid(lon0sst,latsst);
[xq,yq] = meshgrid(lon0NCEP,latNCEP);
sstLr = zeros(xm,ym,tmERA5);
for t = 1: tmERA5
    sst2mid = griddata(x,y,squeeze(sstNCEP(:,:,t))',xq,yq,'linear');
    sstLr(:,:,t) = sst2mid';
end
sstNCEP(sstNCEP>100) = nan;
sstLr(sstLr>100) = nan;
figure;contourf(squeeze(sstNCEP(:,:,7)));
figure;contourf(sstLr(:,:,7))
clear sstNCEP
% READ SH AND LH DATA 192*94
fid = 'NCEPdata/lhtfl.sfc.mon.mean.nc';
ncdisp(fid)
LH = double(ncread(fid,'lhtfl'));% 
time = ncread(fid,'time');
lonHt = double(ncread(fid,'lon'));
latHt = double(ncread(fid,'lat'));
[xmHT,ymHT,tmHT] = size(LH);
%
fid = 'NCEPdata/uswrf.sfc.mon.meanNEW.nc';
ncdisp(fid)
US = double(ncread(fid,'uswrf'));% 
%
fid = 'NCEPdata/dswrf.sfc.mon.meanNEW.nc';
ncdisp(fid)
DS = double(ncread(fid,'dswrf'));% 
%
fid = 'NCEPdata/ulwrf.sfc.mon.meanNEW.nc';
ncdisp(fid)
UL = double(ncread(fid,'ulwrf'));% 
%
fid = 'NCEPdata/dlwrf.sfc.mon.meanNEW.nc';
ncdisp(fid)
DL = double(ncread(fid,'dlwrf'));
%
fid = 'NCEPdata/shtfl.sfc.mon.mean.nc';
ncdisp(fid)
SH = double(ncread(fid,'shtfl'));
% === first, change the spatial resolution from 192*94 into 144*73
LHlrNCEP = zeros(xm,ym,tm);
SHlrNCEP = zeros(xm,ym,tm);
ULlrNCEP = zeros(xm,ym,tm);
DLlrNCEP = zeros(xm,ym,tm);
USlrNCEP = zeros(xm,ym,tm);
DSlrNCEP = zeros(xm,ym,tm);
%  === change grid 
[x,y] = meshgrid(lonHt,latHt);
[xq,yq] = meshgrid(lon0NCEP,latNCEP);
for t = 1: tm
    mid = griddata(x,y,squeeze(LH(:,:,t))',xq,yq,'linear');
    LHlrNCEP(:,:,t) = mid';
    mid = griddata(x,y,squeeze(SH(:,:,t))',xq,yq,'linear');
    SHlrNCEP(:,:,t) = mid';
    mid = griddata(x,y,squeeze(UL(:,:,t))',xq,yq,'linear');
    ULlrNCEP(:,:,t) = mid';
    mid = griddata(x,y,squeeze(DL(:,:,t))',xq,yq,'linear');
    DLlrNCEP(:,:,t) = mid';
    mid = griddata(x,y,squeeze(US(:,:,t))',xq,yq,'linear');
    USlrNCEP(:,:,t) = mid';
    mid = griddata(x,y,squeeze(DS(:,:,t))',xq,yq,'linear');
    DSlrNCEP(:,:,t) = mid';
end
clear SH LH UL DL US DS
save heatmonthly  LHlrNCEP SHlrNCEP ULlrNCEP DLlrNCEP USlrNCEP DSlrNCEP sstLr