%

cp=computer;

% a Matlab routine to calculate transport
%
whsec=menu('Which Section?','GSC','SFlank','Both');
%
disp('Note: common time for 1985 deployment is 241.135 - 310.677');
disp('     while sept-oct runs from 244-305.');
disp('     while Apr-May runs from 115-129.5');
disp(' ');
disp('Note: common time for 1978 deployment is 272.25 to 352.125')
disp('       (in yrday1_gmt convention)')
disp(' ');
disp('Note: common time for 1982 deployment is 223.0417 to  265.4583')
disp('Note: common time for 1995 deployment is 34 to 136');
%
start=menu('enter initial time in yeardays: ','241.135','244',...
      '272.25','115','223.0417','34');
if start==1;
  start=241.135;
elseif start==2;
  start=244;
elseif start==3;
  start=272.25;
elseif start==4;
  start=115;
elseif start==5;
  start=223.0417;
elseif start==6;
  start=34;
end;

tend=menu('enter end time in yeardays: ','310.677','305',...
  '129.5','352.125','265.4583','136');
if tend==1;
  tend=310.677;
elseif tend==2;
  tend=305;
elseif tend==3;
  tend=129.5;
elseif tend==4;
  tend=352.125;
elseif tend==5;
  tend=265.4583;
elseif tend==6
  tend=136;
end;
%
%
if whsec==1;
%
% first for the GSC and then for the southern flank
% Note that some lines are commented out depending
% on what season the transport is for (AM or SO).
%
  whd=menu('which deployment? ','whoi/usgs 1978','nmfs 1985');
  if whd==1;
    load /data4/jmanning/moor/gsc/mm1.mtl
    load /data4/jmanning/moor/gsc/mm2.mtl
    load /data4/jmanning/moor/gsc/mn1.mtl
    load /data4/jmanning/moor/gsc/mn2.mtl
    load /data4/jmanning/moor/gsc/mo2.mtl
    load /data4/jmanning/moor/gsc/mo3.mtl
%
% Find the records that fit the common time period
%
%
    im1=find(mm1(:,2)>=start&mm1(:,2)<=tend);
    im2=find(mm2(:,2)>=start&mm2(:,2)<=tend);
    in1=find(mn1(:,2)>=start&mn1(:,2)<=tend);
    in2=find(mn2(:,2)>=start&mn2(:,2)<=tend);
    io2=find(mo2(:,2)>=start&mo2(:,2)<=tend);
    io3=find(mo3(:,2)>=start&mo3(:,2)<=tend);
    %
    % multiply northward velocity by weighting = f(dx,dz)
    %
    t1=mm1(im1,4)*3720;
    t2=mm2(im2,4)*3720;
    t3=mn1(in1,4)*1572;
    t4=mn2(in2,4)/10*3930;
    t5=mo2(io2,4)*1110;
    t6=mo3(io3,4)*2220;
    t7=mo2(io2,4)*1443;
%
% define yearday  given the latest deployment
%
    yd=mm2(im2,2);


  else;
    load /data4/jmanning/moor/gsc/m21a.mtl
    load /data4/jmanning/moor/gsc/m23.mtl
    load /data4/jmanning/moor/gsc/m22.mtl
    load /data4/jmanning/moor/gsc/m31cfix.dat
    load /data4/jmanning/moor/gsc/m32.mtl
    load /data4/jmanning/moor/gsc/m411.mtl
    load /data4/jmanning/moor/gsc/m412.mtl
    load /data4/jmanning/moor/gsc/m51cfix.dat
    load /data4/jmanning/moor/gsc/m52cfix.dat
    load /data4/jmanning/moor/gsc/m53cfix.dat
    load /data4/jmanning/moor/gsc/m62cfix.dat
    load /data4/jmanning/moor/gsc/m64cfix.dat
%
% Find the records that fit the common time period
%
%
    i21a=find(m21a(:,2)>=start&m21a(:,2)<=tend);
    i23=find(m23(:,2)>=start&m23(:,2)<=tend);
    i22=find(m22(:,2)>=start&m22(:,2)<=tend);
    i31c=find(m31cfix(:,2)>=start&m31cfix(:,2)<=tend);
    i32=find(m32(:,2)>=start&m32(:,2)<=tend);
    i411=find(m411(:,2)>=start&m411(:,2)<=tend);
    i412=find(m412(:,2)>=start&m412(:,2)<=tend);
    i51c=find(m51cfix(:,2)>=start&m51cfix(:,2)<=tend);
    i52c=find(m52cfix(:,2)>=start&m52cfix(:,2)<=tend);
    i53c=find(m53cfix(:,2)>=start&m53cfix(:,2)<=tend);
    i62c=find(m62cfix(:,2)>=start&m62cfix(:,2)<=tend);
    i64c=find(m64cfix(:,2)>=start&m64cfix(:,2)<=tend);
  %
%
% define yearday 
%
    yd=m31cfix(i31c,2);

% assign weights for different sites so that they add up to one
     w2=[0;0;0;.1;.5;.7;.8];
     w3=[.2;.3;.5;.8;.5;.3;.2];
     w4=[.8;.7;.5;.1;0;0;0];
     
% multiply northward velocity by weighting = f(dx,dz)
%
    inc=0;
    for is=1:6;
      for iw=1:6;

% note that each site has a weight associated with it
% that is dependent on the "is" and "iw" value

 
        tc2=m23(i23,4)*10*w2(iw)*40*10;
        tc3=m22(i22,4)*10*w2(iw)*40*10;
        tc4=m31cfix(i31c,4)*10*w3(iw)*20*10;
        tc5=tc4*0.79;
        tc6=m412(i412,4)*10*w4(iw)*40*10;
        tf1=m412(i412,3)*10*w4(is)*40*10;
        tf2=m51cfix(i51c,3)*10*w3(is)*30*10;
        tf3=m52cfix(i52c,3)*10*w3(is)*10*10;
        tf4=m62cfix(i62c,3)*10*w2(is)*40*10;
 

  %
  % sum up transport for each hour at this location
  %
%        for lag=24:24:372;
          lag=4*24;
          clear tc tf;
          for k=lag+1:length(tc3);

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%  edit the following line according to run parameters

            tc(k)=tc3(k)+tc4(k)+tc5(k)+tc6(k);
            tf(k)=tf1(k-lag)+tf2(k-lag)+tf3(k-lag)+tf4(k-lag);
 
          end;

% calculate cumulative transport
          if is==1;
            tfm(is,iw)=mean(tf);
          else;
            tfm(is,iw)=mean(tf)+sum(tfm(1:is-1,iw));
          end;
          if iw==1;
            tcm(is,iw)=mean(tc);
          else;
            tcm(is,iw)=mean(tc)+sum(tcm(is,1:iw-1));
          end;
          pct(is,iw)=tcm(is,iw)/tfm(is,iw)*100.;

%          ik=lag/24;
%          keep(ik)=mean(mean(pct));
%        end;
%         inc=inc+1;
%         idk(inc)=find(keep==min(keep))
      end;
    end;
    disp('best lag = ');
    mean(idk)

    disp('contour pct');
    pause;
    
  end;


  subplot(111);
  %
  subplot(221);
  title('raw inchannel transport')
  ylabel('m**3/s');
  plot(yd,t(1:length(yd)))
  title(num2str(mean(t)));
  pause(5);
  %
  % filter out tides
  %
  nyqf=12;
  cof=(24/33)/nyqf;
  [b,a]=butter(3,cof);
  f=filtfilt(b,a,t);
  ix=find(yd>274&yd<349);
  subplot(222);
  plot(yd(ix),f(ix));
  title('low pass inchannel transport');
  %
  subplot(223)
  f1=filtfilt(b,a,t1);
  f2=filtfilt(b,a,t2);
  f3=filtfilt(b,a,t3);
  f4=filtfilt(b,a,t4);
  f5=filtfilt(b,a,t5);
  f6=filtfilt(b,a,t6);
  text(min(yd(ix))+1,-1,'         contributions');
  xlabel('yearday');
  plot(yd(ix),f1(ix),yd(ix),f2(ix),yd(ix),f3(ix),yd(ix),f4(ix),yd(ix),f5(ix),yd(ix),f6(ix));
  keyboard
  %
  % save filtered inchannel as "tich"
  %
  tich=f;
  %
end;
%
if whsec==2;
    %
    % now start on the southern flank transport
    %
  wy=menu('Which year?','1982','1985','1995');
  if wy==1;
     cp=computer;

     if cp(1:3)=='SGI';
        load /data5/jmanning/moor/moor82/wcrhr

     else;

        load f:\moor\moor82\wcrhr;

     end;   

%   calculate transport associated with each instrument

    wib=menu('What Lateral Extent or delimiting factor?',...
                                    '120m isobath only',...
                                    'entire transect',...
                                    'f(temperature)',...
                                    'f(salt)',...
                                    'shear site 5',...
                                    'f(salt)&shear site 5')

    %   for each case define the transport components;
    %   # of instruments    %   (noi)
    %   instrument codes for legend
    if wib==1;
      noi=4;
      leg=['Total';'   51';'   61';'   62';'   63'];
      t51=y(:,11)/100*60*42000;
      t61=y(:,13)/100*20*32500;
      t62=y(:,15)/100*40*32500;
      t63=y(:,17)/100*40*32500;
      yt=-1*[t51 t61 t62 t63];
      tr=-1*(t51+t61+t62+t63);
    elseif wib==2;
      noi=8;
      leg=['Total';'   51';'   61';'   62';'   63';'   32';'   33';'   34';'   71'];
      t51=y(:,11)/100*60*42000;
      t61=y(:,13)/100*20*32500;
      t62=y(:,15)/100*40*32500;
      t63=y(:,17)/100*40*32500;
      t31=y(:,5)/100*60*14000;
      t32=y(:,7)/100*140*14000;
      t33=y(:,9)/100*120*14000;
      t71=y(:,19)/100*320*5000;
      yt=-1*[t51 t61 t62 t63 t31 t32 t33 t71];
      tr=-1*(t51+t61+t62+t63+t31+t32+t33+t71);
    elseif wib==3;
      yc=y;
      load /data5/jmanning/moor/moor82/wcrthr
      noi=8;
      leg=['Total';'   51';'   61';'   62';'   63';'   31';'   32';'   33';'   71'];
      % the following are the indices for instruments wanted (first temp &
      % then current)
      tid=[6 7 8 9 3 4 5 10];
      cid=[11 13 15 17 5 7 9 19];
      %      limt=input('What max temp?');
      % we set max tmeps exceptable according to AFAP offshore climatolgy
      % adjusted for depth
%      limt=[22 22 18 13 22 18 13 22]-5;
      limt=[30 30 30 30 30 30 30 30];
      depw=[60 20 40 60 140 120 320];
      xdis=[42000 32500 32500 14000 14000 14000 5000];
      for j=1:8;
        yc(find(y(:,tid(j))>limt(j)),cid(j))=zeros(length(find(y(:,tid(j))>limt(j))),1);
      end;
      t51=yc(:,11)/100*60*42000;
      t61=yc(:,13)/100*20*32500;
      t62=yc(:,15)/100*40*32500;
      t63=yc(:,17)/100*40*32500;
      t31=yc(:,5)/100*60*14000;
      t32=yc(:,7)/100*140*14000;
      t33=yc(:,9)/100*120*14000;
      t71=yc(:,19)/100*320*5000;
      yt=-1*[t51 t61 t62 t63 t31 t32 t33 t71];
      tr=-1*(t51+t61+t62+t63+t31+t32+t33+t71);
      dtdy=y(:,10)-y(:,7); %cross bank temp gradient
    elseif wib==4;
      noi=6;
      leg=['Total';'   51';'   61';'   62';'   63';'   32';'   71'];
      t51=y(:,11)/100*60*42000;
      t61=y(:,13)/100*20*32500;
      t62=y(:,15)/100*((40*32500)-(5000*24)/2-9000*23);
      t63=y(:,17)/100*(16*18000)/2;
      t31=y(:,5)/100*((60*14000)-(14000*11)/2-14000*10);
      t71=y(:,19)/100*50*5000;
      yt=-1*[t51 t61 t62 t63 t31 t71];
      tr=-1*(t51+t61+t62+t63+t31+t71);
    elseif wib==5;
      noi=8;
      leg=['Total';'   51';'   61';'   62';'   63';'   32';'   33';'   34';'   71'];
      nlev=52;

      if cp(1:3)=='SGI';
         load /home3/ocn/jmanning/wcr/m5shear

      else;

         load e:\wcr\m5shear

      end;   
      dv=depthvel';
      wgt=(-1*dv(:,2))/18.1; % velocity normalized by maximum
      for i=1:nlev;
        t51(:,i)=y(:,11)/100*1*42000*wgt(i);
      end;
      t51old=y(:,11)/100*60*42000;
      t61=y(:,13)/100*20*32500;
      t62=y(:,15)/100*40*32500;
      t63=y(:,17)/100*40*32500;
      t31=y(:,5)/100*60*14000;
      t32=y(:,7)/100*140*14000;
      t33=y(:,9)/100*120*14000;
      t71=y(:,19)/100*320*5000;
      yt=-1*[sum(t51')' t61 t62 t63 t31 t32 t33 t71];
      tr=-1*(sum(t51')'+t61+t62+t63+t31+t32+t33+t71);
    elseif wib==6;
      noi=6;
      leg=['Total';'   51';'   61';'   62';'   63';'   32';'   71'];
      nlev=52;

      if cp(1:3)=='SGI';
         load /home3/ocn/jmanning/wcr/m5shear

      else;

         load e:\wcr\m5shear;

      end;   
      dv=depthvel';
      wmeth=menu('Which shear method?','relative to 18.1',...
             'shifted eachtime');
      if wmeth==1;
        wgt=-1*dv(:,2)/18.2; % velocity normalized by maximum
        for i=1:nlev;
          t51(:,i)=y(:,11)/100*1*42000*wgt(i);
        end;
      elseif wmeth==2;
        t51=ones(length(y),nlev+2);
        % shift each velocity
        for i=1:length(y);
          wgt=dv(:,2)+(y(i,11)+18.0955);
          t51(i,1:nlev)=wgt'/100*1*42000;
        % add top & bottom levels
          t51(i,nlev+1)=wgt(nlev)/100*4.5*42000;
          t51(i,nlev+2)=wgt(1)/100*3.5*42000;
        end;
        
      end;
      t51old=y(:,11)/100*60*42000;
      t61=y(:,13)/100*20*32500;
      t62=y(:,15)/100*((40*32500)-(5000*24)/2-9000*23);
      t63=y(:,17)/100*(16*18000)/2;
      t31=y(:,5)/100*((60*14000)-(14000*11)/2-14000*10);
      t71=y(:,19)/100*50*5000;
      yt=-1*[sum(t51')' t61 t62 t63 t31 t71];
      tr=-1*(sum(t51')'+t61+t62+t63+t31+t71);
    end;
    disp(['mean sflank transport = ',num2str(meannnan(tr))]);
    disp(['std sflank transport = ',num2str(std(tr))]);
    %
    yd=t;
    t=tr;
    %
    % filter out tides
    %
    f=plfilt(t);
    if wib==4;fshf=f;end;
    plot(yd,f,'linewidth',3,'color','g');
    hold on;
    if whsec==4;
      plot(yd,fshf,'--','linewidth',3,'color','g');
    end;
    title('Alongbank Transport');
    text(min(yd),maxnnan(f),['Mean Transport = ',num2str(mean(t),5),'m3/s'],...
      'verticalalignment','bottom');
    text(min(yd),maxnnan(f)-100000,['      Std Dev  = ',num2str(std(t),5),'m3/s']);
    %
    hold on;

% plot the individual detided contributions

    load ccodemat
    

    for i=1:noi;
      yf(:,i)=plfilt(yt(:,i));
      if i==8;

         plot(yd,yf(:,i),'--','color','k')

      else;

         cc=ccode(i,:);
        plot(yd,yf(:,i),'color',cc)
      end;
    end;

    cd c:\matlab\toolbox\matlab\graph2d
    legend(leg,-1) 
    ylabel('m3/s')
    yd2calab(82)
    %
    % save filtered transport to a .mat file
    %
%    eval(['save /data5/jmanning/moor/moor82/trnsprt82',num2str(wib),'  yd yt yf t f yf']);


  elseif wy==2;

    load /data4/jmanning/moor/gsc/m411.mtl
    load /data4/jmanning/moor/gsc/m412.mtl
    load /data4/jmanning/moor/gsc/m51cfix.dat
    load /data4/jmanning/moor/gsc/m52cfix.dat
    load /data4/jmanning/moor/gsc/m53cfix.dat
    load /data4/jmanning/moor/gsc/m62cfix.dat
    load /data4/jmanning/moor/gsc/m64cfix.dat
    %
    % Find the records that fit the common time period
    %
    i411=find(m411(:,2)>=start&m411(:,2)<=tend);
    i412=find(m412(:,2)>=start&m412(:,2)<=tend);
    i51c=find(m51cfix(:,2)>=start&m51cfix(:,2)<=tend);
    i52c=find(m52cfix(:,2)>=start&m52cfix(:,2)<=tend);
    i53c=find(m53cfix(:,2)>=start&m53cfix(:,2)<=tend);
    i62c=find(m62cfix(:,2)>=start&m62cfix(:,2)<=tend);
    i64c=find(m64cfix(:,2)>=start&m64cfix(:,2)<=tend);
    %

    %clear t1 t2 t3 t4 t5 t6  t f1 f2 f3 f4 f5 f6 f

    %
    % the following is the original weights used
    %
    %  t1=m411(i411,3)*5200;
    %  t1=m412(i412,3)*5200;
    %  t2=m51cfix(i51c,3)*4218;
    %  t3=m52cfix(i52c,3)*3774;
    %  t4=m53cfix(i53c,3)*4884;
    %  t5=m62cfix(i62c,3)*9640;
    %  t6=m64cfix(i64c,3)*14460;

    % the following is the "equivalent weights"
    % area is same as above for the three records
          area=7350.7

    %  t1=m411(i411,3)*5200;
       t1=m412(i412,3)*area;
       t2=m51cfix(i51c,3)*area*2/3;
       t3=m52cfix(i52c,3)*area*1/3;
    %  t4=m53cfix(i53c,3)*4884;
       t5=m62cfix(i62c,3)*area;
    %  t6=m64cfix(i64c,3)*14460;
    %
    for k=1:length(t2);
     t(k)=t1(k)+t2(k)+t3(k)+t5(k);
    end;
    %
    disp(['mean sflank transport = ',num2str(mean(t))]);
    disp(['std sflank transport = ',num2str(std(t))]);
    %
    yd=m51cfix(i51c,2);
    subplot(111);
    %
    subplot(221);
    title('raw along-bank transport')
    plot(yd,t)
    ylabel(' transport m**3/s)')
    title(num2str(mean(t)));
    %
    % filter out tides
    %
    nyqf=12;
    cof=(24/33)/nyqf;
    [b,a]=butter(3,cof);
    f=filtfilt(b,a,t);
    ix=find(yd>246&yd<308);
    subplot(222);
    plot(yd(ix),f(ix));
    title('low pass along-bank transport');
    %
    subplot(223)
    f1=filtfilt(b,a,t1);
    f2=filtfilt(b,a,t2);
    f3=filtfilt(b,a,t3);
    f4=filtfilt(b,a,t4);
    f5=filtfilt(b,a,t5);
    f6=filtfilt(b,a,t6);
    text(250,-1500000,'         contributions');
    plot(yd(ix),f1(ix),yd(ix),f2(ix),yd(ix),f3(ix),yd(ix),f4(ix),yd(ix),f5(ix),yd(ix),f6(ix));
    %plot(yd(ix),f1(ix),yd(ix),f2(ix),yd(ix),f3(ix),yd(ix),f4(ix));
    %
    keyboard;
    %
    %
    % save filtered transport to a .mat file
    %
    save fallequiv  yd(ix) tich(ix) f(ix)
    %
  elseif wy==3;
    if cp(1:3)=='SGI';
       load /data5/jmanning/moor/moor95/st2f

    else;

       load f:\moor\moor95\st2f

    end;

%   calculate transport associated with each instrument

    wib=menu('What Lateral Extent or delimiting factor?',...
                                    'entire transect',...
                                    'shear site 5');

    %   for each case define the transport components;
    %   # of instruments    %   (noi)
    %   instrument codes for legend
    if wib==1;
      noi=3;
      leg=[' Total';'   15m';'   35m';'   62m'];
      %get the portion wanted
      id1=find(st21f(:,1)>=start&st21f(:,1)<=tend);
      id2=find(st22f(:,1)>=start&st22f(:,1)<=tend);
      id3=find(st23f(:,1)>=start&st23f(:,1)<=tend);
      t51=st21f(id1,2)/100*20*42000;
      t52=st22f(id2,2)/100*20*42000;
      t53=st23f(id3,2)/100*20*42000;
      yt=-1*[t51 t52 t53];
      tr=-1*(t51+t52+t53);
    elseif wib==2;
      yt=-1*[t51 t61 t62 t63 t31 t32 t33 t71];
      tr=-1*(t51+t61+t62+t63+t31+t32+t33+t71);
    end;
    disp(['mean sflank transport = ',num2str(meannnan(tr))]);
    disp(['std sflank transport = ',num2str(stdnnan(tr))]);
    %
    yd=st21f(:,1);
    plot(yd,tr,'linewidth',2,'color','g');
    hold on;
    title('Alongbank Transport at ST2');
    text(min(yd),maxnnan(tr),['Mean Transport = ',num2str(meannnan(tr),5),'m3/s'],...
      'verticalalignment','bottom');
    text(min(yd),maxnnan(tr)-50000,['      Std Dev  = ',num2str(stdnnan(tr),5),'m3/s']);
    %
    hold on;

% plot the individual detided contributions

    load e:\atsee\ccodemat


    for i=1:noi;
      if i==8;
        plot(yd,yt(:,i),'--','color','k')
     else;

        cc=ccode(i,:);

        plot(yd,yt(:,i),'color',cc)
      end;
    end;
    ylabel('m3/s')
    axis([34 136 -400000 600000])
%    load /data5/jmanning/wind/stress/b1195ws
%    for i=1:length(taux);
%      [nx(i),ny(i)]=rot(-45,taux(i),tauy(i));
%    end;
%    plot(t,-1*nx*1000000,'g','linewidth',2)
%    legend('Westward Transport @ST2','Southeastward Wind Stress @44011')

    cd c:\matlab\toolbox\matlab\graph2d
    legend(leg) 
    yd2calab(95)
    %
    % save filtered transport to a .mat file
    %
%    eval(['save /data5/jmanning/moor/moor95/trnsprt95',num2str(wib),'  yd yt tr']);


  end;

end;
%
if whsec==3;
  load fallequiv;

  % Generate Wind &  percent in-channel transport figure
  %
  hold on;
  %
  subplot(3,1,1);
  load /data4/jmanning/wind/b1185ws
  iw=find(t>start&t<tend&isnan(taux)~=1);
    stickplotjm(t(iw),taux(iw),tauy(iw),[244 310],'Pa')
  ylabel('Wind Stress at Buoy #11');
  set(gca,'xticklabels',[ ]);
  set(gca,'units','normalized','position',[.13 .70 .775 .28]);
  hold on;
  %
  subplot(3,1,2);
  plot(yd,tich/1000000,'g');
  hold on;
  h0=plot(yd,f*(-1)/1000000,'c');
  set(h0,'linewidth',2);
  h1=text(250,-.1,'  In-Channel (northward)');
  set(h1,'color','g');
  h2=text(250,1.0,' Along-Bank (westward)');
  set(h2,'color','c');
  %text(250,.6,'Not Including Moorings #2 and #6');
  set(gca,'xticklabels',[ ]);
  ylabel('Transport (Sv)');
  axis([244 310 -.25 1.25]);
  set(gca,'units','normalized','position',[.13 .38 .775 .28]);

  subplot(3,1,3);
  pc=tich./(f.*(-1));
  plot(yd,pc*100,'r');
  text(250,-30,'(In-Channel/Along-Bank Transport) * 100');
  ylabel('Percent Recirculation');
  axis([244 310 -50 150]);
  set(gca,'units','normalized','position',[.13 .06 .775 .28]);

  set(gca,'xticklabels',[]);
  text(244,-60,'Sept 7');
  text(259,-60,'17');
  text(269,-60,'27');
  text(275,-60,'Oct  7');
  text(289,-60,'17');
  text(299,-60,'27');
  text(307,-60,'Nov 6');
  text(244,-75,'1985');

  set(gcf,'paperposition',[.25,1,8.25,9.5]);
  %
  %
  %fileout='c:\gbr\trnsprt\sepoct.dat';
  %
  %disp(' ');
  %disp([' Please wait while writing output as ',fileout])
  %disp(' ');
  %
  %   for numrow=1:length(ff);
  %    fprintf(fileout,'%10.4f %10.4f %10.4f\n',yd(numrow),ff(numrow),fff(numrow));
%   end;
%
%plot(yd,ff,yd,fff);
end;