% program tpeval         
% tpeval.m  
%*******************************************************************
%                                                                  *
%                 PROGRAM TPEVAL                                   *
%                                                                  *
%*******************************************************************
% basic program: READS INPIT FILES FROM QGSBL MODEL OUTPUT
%& subplot 1: plots nrms difference between simulated and referenced data:
% output from QGSBL model;
% &subplot 2: plots predicted (from model) and observed (CTD) SLA
% at BATS site;
% &subplot 3: plots scatter diagram between predicted (from model)
% and observed (ERS 2 data inside "diamond") SLA;
% &subplot 4: plots scatter diagram between predicted (from model)
% and observed (CTD) SLA time series at BATS site;
% predicted SLA are taken at the same day, as CTD data.
% QGSBL model was run beginning Oct 27, 1992 (day 300) until
% Sep 17, 1998 (day 2450), so the total length of the run is 2151 days.
% 
% gap occured in ERS1/2 mission betwen December 23,1993 and April 29,1995
% ( days 722 to 1214 ).
% START EVALUATION 
% load source files: T/P and ERS data
% ADD PATH   ADD PATH    ADD PATH     ADD PATH
     addpath /usr/people/valery/TPD_1.0/data/
% *********************************************
     load tp36ers2a13_1;            
     va1=tp36ers2a13_1;
     clear tp36ers2a13_1;
     load tp36ers2a13_2;            
     va2=tp36ers2a13_2;
     clear tp36ers2a13_2;

% set up subplots position
    h1=subplot(411);
    get(h1,'Position');
    h2=subplot(412);
    get(h2,'Position');
    h3=subplot(413);
    get(h3,'Position');
    h4=subplot(414);
      set(h1,'Position',[0.1300    0.7012    0.7750    0.2238]);
      set(h2,'Fontweigh','bold','Fontsize',10,'Position',[0.1300    0.4056    0.7750    0.2238]);
      set(h3,'Position',[0.1300    0.1100    0.3375    0.2238]);
      set(h4,'Fontweigh','bold','Fontsize',10,'Position',[0.5675    0.1100    0.3375    0.2238]);
    drawnow
    set(gcf,'DefaultAxesFontSize',6)
%  load and transform from PSI to SLA; calculate RMS; 
%  (SLA=PSI/2.3162 cm)
      load dh3a15b5f.vbr
      b=dh3a15b5f;
      b=b/2.3162; 

      load dh3a15b5f.vbc
      a=dh3a15b5f;
      clear dh3a15b5f
      a=a/2.3162;
% RMS
    rsl=sqrt(mean((a(1:2151,:)-b(1:2151,:)).^2,2));
% dimensional mean
%    rms=mean(rsl);
% referenced SLA (streamfunction) rms
    sf=sqrt(mean((b(1:2151,:)).^2,2));
% normalized RMS difference
    v=rsl./mean(sf)*100;
% normaliaed mean
    nmrs=mean(v)
% draw plot
    xx=300:2450;
    set(gcf,'CurrentAxes',h1)
    plot(xx,v,'r-');
    hold on
    h5=area(xx,v);
% draw gap in ERS data
    xz=722:1:1214;
    yz(1:493)=85;
    plot(xz,yz,'-','Linewidth',1.5)
    hold off
    tt=text(820,90,['NO ERS DATA'],'Fontsize',6)
    text(370,90,['NRMS=',num2str(nmrs,'%5.1f'),'%'],'HorizontalAlignment','left','fontsize',6);
    axis([300 2500 0 100]);
    grid on
    set(gca,'XTick',[366 731 1096 1461 1827 2192])
    set(gca,'XTickLabel',[1993 1994 1995 1996 1997 1998])
    ylabel('NRMS difference (%)','fontweigh','bold','fontsize',10);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% TIME SERIES AT BATS SITE
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Find and extract streamfunction value at BATS site (inside "diamond");
% compare simulated and referenced SLA values to the observed (dynamic heigh
% from CTD data). 
       t = 0;                   % Loop counter
%   set coordinate of BATS site and coordinate system origin
       f=pi*64./180;
% BATS coordinates (degrees)
       x0=-64.17;                                                                                                                                    
       y0=31.74;
% small domain coordinate system origin (relative)
       yn=10.2728;
       xn=23.4591;
% create "small" domain
       x=ones(25,31);
       y=ones(25,31);

       for k=1:1:25
       for l=1:1:31
       x(k,l)=(xn)+(-(k-1)*sin(f)+(l-1)*cos(f))*0.366;
       y(k,l)=(yn)+((k-1)*cos(f)+(l-1)*sin(f))*0.366;
       end
       end
% calculate geographical coordinates
       xi=-70+(x-.85)*10/34;     
       yi=28.25+(y-1)*8/32;
% "diamond" perimeter
       xv=[xi(1,1);xi(1,31);xi(25,31);xi(25,1);xi(1,1)];
       yv=[yi(1,1);yi(1,31);yi(25,31);yi(25,1);yi(1,1)];
% double check BATS site choice
       sa=find(abs(xi-x0)<=.05 & abs(yi-y0)<=.05)
       xi(sa)
       yi(sa)
       [r,q]=find(abs(xi-x0)<=.05 & abs(yi-y0)<=.05)
       xi(r,q)
       yi(r,q)
% load dynheig
       load dynh20.dat;
       c=dynh20;
% remove mean (dyn height in cm!)
       c(:,4)=(c(:,4)-mean(c(:,4),1));
       c(1:1:42,:)=[];
% extract data
     sbats=ones(2151,3);
     xl=300:2450;
% LOOP: create PSI time series 
     for i = [1:2151]   
     a1=a(i,:);
     a1=reshape(a1,31,25);
     a1=a1';

     b1=b(i,:);
     b1=reshape(b1,31,25);
     b1=b1';

     sbats(i,1)=i+299;
     sbats(i,3)=b1(r,q);
     sbats(i,2)=a1(r,q);
     end
% remove mean
    sbats(:,2)=sbats(:,2)-mean(sbats(:,2));
% 
     dyn=c;
     ta=ones(1,3);
     tb=ones(1,3);
% start extraction
     for k=[1:size(dyn(:,1))]
     ch=dyn(k,3);
     rr=find(sbats(:,1)==ch);
     tb=[tb;sbats(rr,:)];
     end
%
     tb(1,:)=[];
     dyn(72:76,:)=[];
%    cr1=corrcoef(sbats(:,2),sbats(:,3))
     cr2=corrcoef(dyn(:,4),tb(:,2))
% regression (X and y in cm)
     X=[ones(length(dyn),1),dyn(:,4)];
     [bu,bint,r,rint,stat]=regress(tb(:,2),X);
% DRAW SUBPLOT 2

     set(gcf,'CurrentAxes',h2) 

     plot(xl,sbats(:,2),'r-','Linewidth',0.5);
     hold on
     plot(c(:,3),c(:,4),'xk-','Linewidth',1);
     hold off
     axis([300 2500 -20 30]);
     grid on
     set(gca,'XTick',[366 731 1096 1461 1827 2192])
     set(gca,'XTickLabel',[1993 1994 1995 1996 1997 1998])
     
     set(get(h2,'Ylabel'),'String','Predicted \DeltaD (cm)')
     xlabel('Time (years)','fontweigh','bold','fontsize',10);
     legend('simulated (QG model)','observed (CTD)');
% DRAW SUBPLOT 4 (cm)

     rbm=sqrt(mean((dyn(:,4)-tb(:,2)).^2));
     set(gcf,'CurrentAxes',h4);
     plot(dyn(:,4),tb(:,2),'ok'); 
     axis square
     axis([-40 50 -40 50]);
     g1=refline;
     g2=refline(1,0);
     set(g2,'Linestyle','--')
  text(-36,40,['RMS=',num2str(rbm,'%5.1f'),' cm'],'HorizontalAlignment','left','fontsize',7);
  text(29,-25,['r=',num2str(cr2(1,2),'%0.2g')],'HorizontalAlignment','left','fontsize',6);
  text(29,-30,['a=',num2str(bu(1),'%5.2f')],'HorizontalAlignment','left','fontsize',6);
  text(29,-35,['b=',num2str(bu(2),'%5.2f')],'HorizontalAlignment','left','fontsize',6); 
  set(get(h4,'Ylabel'),'String','Predicted \DeltaD (cm)')
  set(get(h4,'Xlabel'),'String','Observed \DeltaD (cm)')
  set(gcf,'PaperPosition',[0.75 1 7 9]);
%
%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$44
% t is a given JD since 01/01/92; i=1 corresponds to t=i+299=300
% find out ERS data from total sets
     n=0; ib=901; ie=2151;
     sla=ones(1,3);
     slar=ones(1,3);
     sers =ones(1,5);
% LOOP !!!!!!
     for i=[ib:ie]
     t=i+299;
     tmin=t-n;
     tmax=t+n+1;
% find ERS points at  a given day
     r1=find(va1(:,1)>=tmin & va1(:,1)<=tmax & va1(:,5)==15.0);
     r2=find(va2(:,1)>=tmin & va2(:,1)<=tmax & va2(:,5)==15.0);
     v1=va1(r1,:);
     v2=va2(r2,:);
% find ERS pointS INSIDE "diamond" box;
     f1=find(v1(:,3)>=-66.0 & v1(:,3)<=-62.0 & v1(:,2)>=30.0 & v1(:,2)<=34.0);
     vv1=v1(f1,:);
     f2=find(v2(:,3)>=-66.0 & v2(:,3)<=-62.0 & v2(:,2)>=30.0 & v2(:,2)<=34.0);
     vv2=v2(f2,:);
     u1=inpolygon(vv1(:,3),vv1(:,2),xv,yv);
     u2=inpolygon(vv2(:,3),vv2(:,2),xv,yv);
     er1=vv1(u1,:);
     er2=vv2(u2,:);
     er2=[er1;er2];
    clear v1;
    clear v2;
% find ERS sla at selected alongtrack points
     if isempty(er2)==0
% find corresponding points from a  model grid:
% predicted sla
      a1=a(i,:);
      a1=a1';
      as=reshape(a1,31,25);
      as=as';
% referenced sla
      b1=b(i,:);
      b1=b1';
      bs=reshape(b1,31,25);
      bs=bs';
% interpolate sla and slaref values from model grid onto ERS points
      xe=er2(:,3);
      ye=er2(:,2);
      am2=griddata(xi,yi,as,xe,ye);
      bm2=griddata(xi,yi,bs,xe,ye);
      sla=[sla; xe ye am2];
      slar=[slar; xe ye bm2];
      sers=[sers; er2]; 
      clear er1;
      clear er2;
      clear am2;
      clear bm2;
      clear xe;
      clear ye;
      clear f1;
      clear f2;
% end if isempty
      end 
% end of loop
      end 
%
      sla(1,:)=[];
      slar(1,:)=[];
      sers(1,:)=[];
      sers(:,4)=0.1*sers(:,4);
      c1=corrcoef(sla(:,3),sers(:,4))
%     
      X=[ones(length(sers),1),sers(:,4)];
      [bu,bint,r,rint,stat]=regress(sla(:,3),X);
%
      set(gcf,'CurrentAxes',h3);
      plot(sers(:,4),sla(:,3),'.b','Markersize',2.5)
      g1=refline;
      set(g1,'color','k')
      g2=refline(1,0);
      set(g2,'color','k','Linestyle','--')
      axis square
      axis([-40 50 -40 50]);

      rm=sqrt(mean((sla(:,3)-sers(:,4)).^2));
 
   text(-36,40,['RMS=',num2str(rm,'%5.1f'),' cm'],'HorizontalAlignment','left','fontsize',7);
   text(29,-25,['r=',num2str(c1(1,2),'%0.2g')],'HorizontalAlignment','left','fontsize',6);
   text(29,-30,['a=',num2str(bu(1),'%5.2f')],'HorizontalAlignment','left','fontsize',6);
   text(29,-35,['b=',num2str(bu(2),'%5.2f')],'HorizontalAlignment','left','fontsize',6); 
   ylabel('Predicted SLA (cm)','fontweigh','bold','fontsize',10)
   xlabel('Observed SLA (cm)','fontweigh','bold','fontsize',10)
%44444444444444444444444444444444444444444444444444444444444444444
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now, calculate ERS vs OA
      figure
      c3=corrcoef(slar(:,3),sers(:,4))
%     
      X=[ones(length(sers),1),sers(:,4)];
      [bu,bint,r,rint,stat]=regress(slar(:,3),X);
%
      plot(sers(:,4),slar(:,3),'.b','Markersize',2.5);
      g1=refline;
%
      set(g1,'color','k')
      g2=refline(1,0);
      set(g2,'color','k','Linestyle','--')
      axis square
      axis([-40 50 -40 50]);

      rmr=sqrt(mean((slar(:,3)-sers(:,4)).^2));
 
   text(-36,40,['RMS=',num2str(rmr,'%5.1f'),' cm'],'HorizontalAlignment','left','fontsize',7);
   text(29,-25,['r=',num2str(c3(1,2),'%0.2g')],'HorizontalAlignment','left','fontsize',6);
   text(29,-30,['a=',num2str(bu(1),'%5.2f')],'HorizontalAlignment','left','fontsize',6);
   text(29,-35,['b=',num2str(bu(2),'%5.2f')],'HorizontalAlignment','left','fontsize',6); 
   xlabel('Observed SLA (cm)','fontweigh','bold');
   ylabel('OA  SLA (cm)','fontweigh','bold');
 

 
