%
%  If you uncomment the "C!" lines in LTE_linear.f, LTE_nonlinear.f,
%  LTEice_linear.f or LTEice_nonlinear.f, it will generate files pde2d.m 
%  and pde2d.rdm.
%  To see MATLAB plots, replace pde2d.m with this program and run it, 
%  it will read the solution from pde2d.rdm and produce plots of U,V,Eta 
%  (and Q, if applicable) during the last period and the geographical 
%  distribution of the time-averaged energy dissipation rate.
%
 fid = fopen('pde2d.rdm');
 % Read solution from pde2d.rdm
 NSAVE = fscanf(fid,'%g',1);
 NEQN = fscanf(fid,'%g',1);
 NP1 = fscanf(fid,'%g',1);
 NP2 = fscanf(fid,'%g',1);
 if (NP1<1 || NP2<1)
    error('NP1 or NP2 < 1')
 end
 L0 = 0;
 T = zeros(NSAVE+1,1);
 P1 = zeros(NP2+1,NP1+1);
 P2 = zeros(NP2+1,NP1+1);
 X = zeros(NP2+1,NP1+1);
 Y = zeros(NP2+1,NP1+1);
 U = zeros(NP2+1,NP1+1,NSAVE+1,3,NEQN);
 for i=0:NP1
 for j=0:NP2
    P1(j+1,i+1) = fscanf(fid,'%g',1);
    P2(j+1,i+1) = fscanf(fid,'%g',1);
    X(j+1,i+1) = fscanf(fid,'%g',1);
    Y(j+1,i+1) = fscanf(fid,'%g',1);
 end
 end
 for l=0:NSAVE
 T(l+1) = fscanf(fid,'%g',1);
 for ieq=1:NEQN
 for ider=1:3
 for i=0:NP1
 for j=0:NP2
    U(j+1,i+1,l+1,ider,ieq) = fscanf(fid,'%g',1);
 end
 end
 end
 end
 end
 xmin = min(min(X(:,:)));
 xmax = max(max(X(:,:)));
 ymin = min(min(Y(:,:)));
 ymax = max(max(Y(:,:)));
 hx = 0.1*(xmax-xmin);
 hy = 0.1*(ymax-ymin);
 % Surface plots of each variable
             for IEQ=1:NEQN
 % Plot U_IEQ,  if IDER=1
 %      Ux_IEQ, if IDER=2
 %      Uy_IEQ, if IDER=3
 IDER = 1;
 umin = min(min(min(U(:,:,181:201,IDER,IEQ))));
 umax = max(max(max(U(:,:,181:201,IDER,IEQ))));
 if (umax==umin); umax = umin+1; end
 for L=180:5:200
    figure
    surf(X,Y,U(:,:,L+1,IDER,IEQ))
    axis([xmin xmax ymin ymax umin umax])
    xlabel('Longitude (rad.)','FontSize',15)
    ylabel('Latitude (rad.)','FontSize',15)
    if (IEQ==1)
       zlabel('U (m/sec)','FontSize',15)
    end   
    if (IEQ==2)
       zlabel('V (m/sec)','FontSize',15)
    end
    if (IEQ==3)
       zlabel('\eta (m)','FontSize',15)
    end
    if (IEQ==4)
       zlabel('Q (Pa)','FontSize',15)
    end
    title([' T = ',num2str(T(L+1))],'FontSize',15)
    view(-37.5,30.0)
 end
             end
 % Surface and contour plots of dissipation rate
 L = NSAVE;
 IEQ = 2;
 IDER = 3;
 umin = min(min(U(:,:,L+1,IDER,IEQ)));
 umax = max(max(U(:,:,L+1,IDER,IEQ)));
 if (umax==umin); umax = umin+1; end
 figure
 surf(X,Y,U(:,:,L+1,IDER,IEQ))
 axis([xmin xmax ymin ymax umin umax])
 xlabel('Longitude (rad.)','FontSize',15)
 ylabel('Latitude (rad.)','FontSize',15)
 zlabel('Dissipation rate (W/m^2)','FontSize',15)
 view(-37.5,30.0)

 figure
 [C h] = contour(X,Y,U(:,:,L+1,IDER,IEQ));          
 clabel(C,h)               
 xlabel('Longitude (rad.)','FontSize',15)
 ylabel('Latitude (rad.)','FontSize',15)
 title('Dissipation rate (W/m^2)','FontSize',15)
