function yp=fopc_b(t,y,flag,tu,ts,name,ns,np)
% Subroutine for FOPC; bkwd. integ. of adjoint eqns 
% and Qg=int(HuPhi*HuPhi')dt;          8/97, 2/2/98
%
la=formm(y([1:np*ns]),ns);
u=interp1(tu(:,1),tu(:,2:end),t); 
s=interp1(ts(:,1),ts(:,2:end),t);
[fs,fu]=feval(name,u,s,t,3);
HuPhi=la'*fu; Qgd=-HuPhi*HuPhi'; 
lad=-fs'*la; yp=[formm(lad,'c'); forms(Qgd)];


