% {{{ Initialize 
% Input grid
nx=960;
ny=1280;

pn1='/nobackupp2/dmenemen/llc_4320/regions/CalSWOT/';
pn2='~dmenemen/llc_4320/regions/CalSWOT/run_template1/MITgcm/run/';

% Steric height computed by Hector
Dir='/nobackup/htorresg/hires/stericH/global/';

% ECMWF surface pressure



cd ~dmenemen/llc_4320/regions/CalSWOT/run_template1/animations
% }}}

% {{{ Compare Eta, PhiBot, and steric

% {{{ Extract and plot time series 

% Index to extract
IX=400; IY=600;

% {{{ Extract time series 
STR=zeros(167,1);
AP=zeros(167,1);
AT=zeros(167,1);
ETA=zeros(167,1);
PHB=zeros(167,1);
for ts=144:144:24168, mydisp(ts)
    
    st=hdf5read([Dir,'steric_height_llc_' myint2str(ts/144,4) '.h5'],'stericH');
    STR(ts/144)=st(IX,IY);
    
    skip=round((datenum(2012,6,1)-datenum(2012,1,1))*24+ts/144);
    tmp=readbin(Pres,[length(lon) length(lat)],1,'real*4',skip);
    ap=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;
    AP(ts/144)=ap(IX,IY);
    
    skip1=floor(skip/6);
    skip2=ceil(skip/6);
    if skip2==skip1, skip2=skip2+1; end
    tmp1=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip1);
    tmp2=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip2);
    tmp=(skip2-skip/6)*tmp1+(skip/6-skip1)*tmp2;
    at=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;
    AT(ts/144)=at(IX,IY);

    dy=ts2dte(ts,25,2012,6,1,30);
    eg=readbin([pn1 'Eta/Eta_960x1280.' dy],[nx ny]);
    ETA(ts/144)=eg(IX,IY);
    
    pg=readbin([pn1 'PhiBot/PhiBot_960x1280.' dy],[nx ny])/9.81;
    PHB(ts/144)=pg(IX,IY);
end
% }}}

% {{{ Plot time series 

Tme=2:167;
figure(2)
clf
plot(Tme,STR(Tme)-STR(1), ...
     Tme,ETA(Tme)-ETA(1), ...
     Tme,PHB(Tme)-PHB(1), ...
     Tme,AP(Tme)-AP(1), ...
     Tme,AT(Tme)-AT(1),'linewidth',3)
legend('steric','Eta','PhiBot','APtidal','Atmos')
figure(3)
clf
plot(Tme,STR(Tme)-STR(1), ...
     Tme,ETA(Tme)-ETA(1)-PHB(Tme)+PHB(1)-STR(Tme)+STR(1)+AP(Tme)-AP(1), ...
     Tme,ETA(Tme)-ETA(1)-PHB(Tme)+PHB(1)+AP(Tme)-AP(1),'linewidth',3)
legend('steric','Eta-PhiBot-steric-APtidal','Eta-PhiBot-APtidal')

% }}}

% }}}

% {{{ Extract fields at hour 1 (time step 144).
ts=144;
% {{{ Steric height computed by Hector .
fnm=[Dir,'steric_height_llc_' myint2str(ts/144,4) '.h5'];
st1=hdf5read(fnm,'stericH');
% }}}
% {{{ ECMWF surface pressure 
XC=readbin([pn1 'grid/XC_960x1280'],[nx ny]);
XC(find(XC<0))=XC(find(XC<0))+360;
YC=readbin([pn1 'grid/YC_960x1280'],[nx ny]);
Pres='/nobackupp8/dmenemen/forcing/ECMWF_operational/EOG_pres_tide_2012';
Atm='/nobackupp8/dmenemen/forcing/ECMWF_operational/EOG_pres_2012';
lon=0:.140625:359.9;
dy= [0 .1394 .14018 .14039 .1404695 .140496 .1405145 .1405275 .1405375 ...
     .1405455 .140552 .1405575 .140562 .1405655 .140568 .1405695 ...
     ones(1,1249)*.14057 ...
     .1405695 .140568 .1405655 .140562 .1405575 .140552 .1405455 ...
     .1405375 .1405275 .1405145 .140496 .1404695 .14039 .14018 .1394];
lat=-89.8924+cumsum(dy);
skip=round((datenum(2012,6,1)-datenum(2012,1,1))*24+ts/144);
tmp=readbin(Pres,[length(lon) length(lat)],1,'real*4',skip);
ap1=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;
skip1=floor(skip/6);
skip2=ceil(skip/6);
tmp1=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip1);
tmp2=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip2);
tmp=(skip2-skip/6)*tmp1+(skip/6-skip1)*tmp2;
at1=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;
% }}}
% {{{ Sea surface height and bottom pressure
dy=ts2dte(ts,25,2012,6,1,30);
e1g=readbin([pn1 'Eta/Eta_960x1280.' dy],[nx ny]);
p1g=readbin([pn1 'PhiBot/PhiBot_960x1280.' dy],[nx ny])/9.81;
% }}}
% {{{ set land to nan
msk=readbin([pn1 'grid/hFacC_960x1280x87'],[nx ny]);
ix=1:959; iy=1:1279;
st1(find(~msk(ix,iy)))=nan;
ap1(find(~msk))=nan;
at1(find(~msk))=nan;
e1g(find(~msk))=nan;
p1g(find(~msk))=nan;
% }}}
% }}}

% {{{ Plot fields
for ts=288:144:24168, mydisp(ts)
    
    st=hdf5read([Dir,'steric_height_llc_' myint2str(ts/144,4) '.h5'],'stericH');
    STR(ts/144)=st(IX,IY);
    
    skip=round((datenum(2012,6,1)-datenum(2012,1,1))*24+ts/144);
    tmp=readbin(Pres,[length(lon) length(lat)],1,'real*4',skip);
    ap=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;
    
    skip1=floor(skip/6);
    skip2=ceil(skip/6);
    if skip2==skip1, skip2=skip2+1; end
    tmp1=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip1);
    tmp2=readbin(Atm,[length(lon) length(lat)],1,'real*4',skip2);
    tmp=(skip2-skip/6)*tmp1+(skip/6-skip1)*tmp2;
    at=interp2(lat,lon,tmp,YC,XC)/9.81/1027.5;

    dy=ts2dte(ts,25,2012,6,1,30);
    eg=readbin([pn1 'Eta/Eta_960x1280.' dy],[nx ny]);
    
    pg=readbin([pn1 'PhiBot/PhiBot_960x1280.' dy],[nx ny])/9.81;

    clf
    subplot(331), mypcolor(st'-st1'); thincolorbar
    title(['steric height (hr' int2str(ts/144) ' - hr1)'])
    subplot(332), mypcolor(eg(ix,iy)'-e1g(ix,iy)'); thincolorbar
    title(['Eta (hr' int2str(ts/144) ' - hr1)'])
    subplot(333), mypcolor(pg(ix,iy)'-p1g(ix,iy)'); thincolorbar
    title(['PhiBot (hr' int2str(ts/144) ' - hr1)'])
    subplot(334), mypcolor(eg(ix,iy)'-e1g(ix,iy)'-pg(ix,iy)'+p1g(ix,iy)'); thincolorbar
    title(['Eta - PhiBot (hr' int2str(ts/144) ' - hr1)'])
    subplot(335), mypcolor(eg(ix,iy)'-e1g(ix,iy)'-pg(ix,iy)'+p1g(ix,iy)'-st'+st1'); thincolorbar
    title(['Eta - PhiBot - steric (hr' int2str(ts/144) ' - hr1)'])
    subplot(336), mypcolor(ap1'-ap'); thincolorbar
    title(['atmospheric pressure (hr1 - hr' int2str(ts/144) ')'])
    subplot(337), mypcolor(eg(ix,iy)'-e1g(ix,iy)'-pg(ix,iy)'+p1g(ix,iy)'-st(ix,iy)'+st1(ix,iy)'-ap1(ix,iy)'+ap(ix,iy)'); thincolorbar
    title(['Eta - PhiBot - steric - atm (hr' int2str(ts/144) ' - hr1)'])
    subplot(338), mypcolor(eg(ix,iy)'-e1g(ix,iy)'-pg(ix,iy)'+p1g(ix,iy)'-ap1(ix,iy)'+ap(ix,iy)'); thincolorbar
    title(['Eta - PhiBot - atm (hr' int2str(ts/144) ' - hr1)'])
    subplot(339), mypcolor(at1'-at'); thincolorbar
    title(['atmospheric pressure no tides (hr1 - hr' int2str(ts/144) ')'])
    pause(.1)
end
% }}}

% }}}

% {{{ Compare Eta - PhiBot 

cd ~dmenemen/llc_4320/regions/CalSWOT/run_template1/animations
ts=144;
dy=ts2dte(ts,25,2012,6,1,30);
e1g=readbin([pn1 'Eta/Eta_960x1280.' dy],[nx ny]);
p1g=readbin([pn1 'PhiBot/PhiBot_960x1280.' dy],[nx ny]);
e1r=readbin([pn2 'Eta.' myint2str(ts,10) '.data'],[nx ny]);
p1r=readbin([pn2 'PhiBot.' myint2str(ts,10) '.data'],[nx ny]);
msk=readbin([pn1 'grid/hFacC_960x1280x87'],[nx ny]);
clf reset
orient tall
wysiwyg
cx=[-.5 1];
for ts=144:144:24168
    dy=ts2dte(ts,25,2012,6,1,30);
    eg=readbin([pn1 'Eta/Eta_960x1280.' dy],[nx ny]);
    er=readbin([pn2 'Eta.' myint2str(ts,10) '.data'],[nx ny]);
    pg=readbin([pn1 'PhiBot/PhiBot_960x1280.' dy],[nx ny]);
    pr=readbin([pn2 'PhiBot.' myint2str(ts,10) '.data'],[nx ny]);
 
    eg=(eg-(pg-p1g)/9.81).*msk-mean(e1g(find(msk)));
    er=(er-e1r+e1g-(pr-p1r)/9.81).*msk-mean(e1g(find(msk)));
    
    clf
    subplot(311)
    mypcolor(eg');
    caxis(cx)
    thincolorbar
    title(['Eta minus PhiBot from global 1/48 on hour ' int2str(ts/144)])
 
    subplot(312)
    mypcolor(er');
    caxis(cx)
    thincolorbar
    title(['Eta minus PhiBot from regional 1/48 on hour ' int2str(ts/144)])

    subplot(313)
    mypcolor(er'-eg');
    caxis([-1 1]/10)
    thincolorbar
    title('Regional minus global Eta')
        
    pause(1)
end

% }}}

% {{{ Compare Eta 
cd ~dmenemen/llc_4320/regions/CalSWOT/run_template1/animations
mkdir frames
cd frames
fn1=dir([pn1 'Eta/Eta*']);
fn2=dir([pn2 'Eta*data']);
frm=0;
cx=[-2 2];
for h=1:min(length(fn1),length(fn2))
    e1=readbin([pn1 'Eta/' fn1(h).name],[nx ny]);
    for t=1:6;
        clf
        subplot(121)
        mypcolor(e1');
        caxis(cx)
        thincolorbar
        title(['1/48, hour ' int2str(h-1)])
        e2=readbin([pn2 fn2((h-1)*6+t).name],[nx ny]);
        subplot(122)
        mypcolor(e2');
        caxis(cx)
        thincolorbar
        title(['regional, hour ' int2str(h-1) ':' myint2str((t-1)*10)])
        eval(['print -djpeg frame' myint2str(frm,4)])
        frm=frm+1;
    end
end
!module load imagemagick; convert -delay 20 *.jpg Eta.gif
!mv *.gif ..
cd ..
!rm -rf frames
% }}}

% {{{ Compare U 

fn1=dir([pn1 'U/U*']);
fn2=dir([pn2 'U*data']);
cx=[-1 1]/2;
eval(['cd ' pn1])
mkdir frame2
cd frame2
!rm *
frm=0;
for h=1:length(fn1)
    e1=readbin([pn1 'U/' fn1(h).name],[nx ny]);
    for t=1:6;
        clf
        subplot(121)
        mypcolor(e1');
        caxis(cx)
        thincolorbar
        title(['1/48, hour ' int2str(h-1)])
        e2=readbin([pn2 fn2((h-1)*6+t).name],[nx ny]);
        subplot(122)
        mypcolor(e2');
        caxis(cx)
        thincolorbar
        title(['regional, hour ' int2str(h-1) ':' myint2str((t-1)*10)])
        eval(['print -djpeg frame' myint2str(frm,4)])
        frm=frm+1;
    end
end
% cd ~/llc_4320/regions/CalSWOT/frame2
% module load imagemagick
% convert -delay 20 *.jpg U.gif
% mv *.gif ..

% compare V
fn1=dir([pn1 'V/V*']);
fn2=dir([pn2 'V*data']);
cx=[-1 1]/2;
eval(['cd ' pn1])
mkdir frame3
cd frame3
!rm *
frm=0;
for h=1:length(fn1)
    e1=readbin([pn1 'V/' fn1(h).name],[nx ny]);
    for t=1:6;
        clf
        subplot(121)
        mypcolor(e1');
        caxis(cx)
        thincolorbar
        title(['1/48, hour ' int2str(h-1)])
        e2=readbin([pn2 fn2((h-1)*6+t).name],[nx ny]);
        subplot(122)
        mypcolor(e2');
        caxis(cx)
        thincolorbar
        title(['regional, hour ' int2str(h-1) ':' myint2str((t-1)*10)])
        eval(['print -djpeg frame' myint2str(frm,4)])
        frm=frm+1;
    end
end
% cd ~/llc_4320/regions/CalSWOT/frame3
% module load imagemagick
% convert -delay 20 *.jpg V.gif
% mv *.gif ..

% compare Theta
fn1=dir([pn1 'Theta/Theta*']);
fn2=dir([pn2 'T*data']);
cx=[7.5 20.5];
eval(['cd ' pn1])
mkdir frames
cd frames
!rm *
frm=0;
for h=1:length(fn1)
    e1=readbin([pn1 'Theta/' fn1(h).name],[nx ny]);
    clf
    subplot(121)
    mypcolor(e1');
    caxis(cx)
    thincolorbar
    title(['1/48, hour ' int2str(h-1)])
    e2=readbin([pn2 fn2((h-1)*6+1).name],[nx ny]);
    subplot(122)
    mypcolor(e2');
    caxis(cx)
    thincolorbar
    title(['regional, hour ' int2str(h-1)])
    eval(['print -djpeg frame' myint2str(frm,4)])
    frm=frm+1;
end
% cd ~/llc_4320/regions/CalSWOT/frames
% module load imagemagick
% convert -delay 20 *.jpg T.gif
% mv *.gif ..

% compare Salt
fn1=dir([pn1 'Salt/Salt*']);
fn2=dir([pn2 'S*data']);
cx=[32 36];
eval(['cd ' pn1])
mkdir frame4
cd frame4
!rm *
frm=0;
for h=1:length(fn1)
    e1=readbin([pn1 'Salt/' fn1(h).name],[nx ny]);
    clf
    subplot(121)
    mypcolor(e1');
    caxis(cx)
    thincolorbar
    title(['1/48, hour ' int2str(h-1)])
    e2=readbin([pn2 fn2((h-1)*6+1).name],[nx ny]);
    subplot(122)
    mypcolor(e2');
    caxis(cx)
    thincolorbar
    title(['regional, hour ' int2str(h-1)])
    eval(['print -djpeg frame' myint2str(frm,3)])
    frm=frm+1;
end
% cd ~/llc_4320/regions/CalSWOT/frame4
% module load imagemagick
% convert -delay 20 *.jpg S.gif
% mv *.gif ..

% }}}

% {{{ ECMWF forcing 
fnm=dir('/nobackup/dmenemen/forcing/ECMWF_operational/EOG_pres_tide_2012');
lon=0:.140625:359.9;
lat=-89.8924;
lat_inc=[.1394,.14018,.14039,.1404695,.140496,.1405145,.1405275,.1405375, ...
         .1405455,.140552,.1405575,.140562,.1405655,.140568,.1405695, ...
         ones(1,1249)*.14057, ...
         .1405695,.140568,.1405655,.140562,.1405575,.140552,.1405455, ...
         .1405375,.1405275,.1405145,.140496,.1404695,.14039,.14018,.1394];
 precip_nlon       = 2560,
 precip_nlat       = 1280,
 % }}}