function [u,v,adjointu,adjointv,direction]= ...
    grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,...
      index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);
%GRAD_J - compute the gradient of the misfit with respect to the inversed parameter
%
%   this routine is calculates the difference between the calculated and the given velocity
%   field (cost function), then calculates the adjoint state and uses it to determine the gradient of the
%   cost function to update the spatial distribution of viscosity.
%
%   Usage:
%      [u,v,adjointu,adjointv,direction]=grad_J(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);

      
disp(['      Recalculating model velocity with last iteration from velfinder']);
Rhs_parsed=P*(Rhs - F*specified_velocity);
Fprime=P*F*P';

solution=Fprime\Rhs_parsed;
solution=P'*solution + specified_velocity;
u=solution(1:nods);
v=solution(nods+1:2*nods);
disp('      Done');


disp(['      Calculating adjoints.']);
Du=[(u-vx_obs).*weighting.*weighting
   (v-vy_obs).*weighting.*weighting];

[nu2,nu3]=visc_grad(index,nel,alpha,beta,u,v,B_bar);
%G_calc;
ARhs_aprime=-P*(S*Du);
%ADprime=P0*(F-G_bar)*P0';
ADprime=P*F*P';

templambda=P'*(ADprime\ARhs_aprime);
adjointu=templambda(1:nods);
adjointv=templambda((nods+1):(nods*2));
disp('      Done');

disp('      Building the gradJ');
value=zeros(nel,27);
row=zeros(nel,27);
count=0;
for m=1:3
   for k=1:3

	   for l=1:3
         
		   count=count+1;
		  	row(:,count)=index(:,m);

     		
			value(:,count)=-area/3.*z_thick_bar.*nu3.*(...
              (  2.*alpha(:,k).*u(index(:,k))  +beta(:,k).*v(index(:,k))  )  .*2.*alpha(:,l).*adjointu(index(:,l))+...
           (  beta(:,k).*u(index(:,k)) + alpha(:,k).*v(index(:,k))  )  .*( beta(:,l).*adjointu(index(:,l)) +alpha(:,l).*adjointv(index(:,l)) )+...
           (  2*beta(:,k).*v(index(:,k))  +  alpha(:,k).*u(index(:,k))  ).*2.*  beta(:,l).*adjointv(index(:,l))...
           );
      	

      end
   end
end


direction=sparse(row,ones(size(row)),value,nods,1);
direction=full(direction);
disp('      Done');
