mxbl(1:3)=round((.5/.58)*mx);
mxbl(2)=mx+2;
mybl(1:3)=round((.23/.58)*my);
mybl(2)=round((.12/.58)*my);
ixstart=round((.02/.58)*mx)+3;
 if mod(mybl(1),2)==1
     mybl(3)=mybl(3)+1;
     mybl(1)=mybl(1)+1;
 end
 
 ixend=mxbl(2)-mxbl(1)-ixstart;

       for ix=1:mxbl(1)+1
       for iy=1:mybl(1)+1
	      vx=(2.0*epack/mxbl(1))*(ix-(mxbl(1)/2.0)-1);
	      vy=(2.0*epack/mybl(1))*(iy-(mybl(1)/2.0)-1);
	      eta1=log(vx+sqrt(1+vx^2))/log(epack+sqrt(1+epack^2));
              eta2=log(vy+sqrt(1+vy^2))/log(epack+sqrt(1+epack^2));
	  N1=(1/4.0)*(1.0-eta1)*(1.0-eta2);
	  N2=(1/4.0)*(1.0+eta1)*(1.0-eta2);
	  N3=(1/4.0)*(1.0+eta1)*(1.0+eta2);
	  N4=(1/4.0)*(1.0-eta1)*(1.0+eta2);
	  
      xy(3).rz(2,ix,iy)= .06*N1+.03*N2+.29*N3+.29*N4;
	  xy(3).rz(1,ix,iy)=  .02*N1+.52*N2+.4*N3+.25*N4;
	  
% Inverse Jacobian terms
          
	  dzdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.06*(1-eta2)+.03*(1-eta2)+.29*(1+eta2)-.29*(1+eta2));
	  drdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.02*(1-eta2)+.52*(1-eta2)+.4*(1+eta2)-.25*(1+eta2));
	  dzdiy=((2*epack/mybl(3))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.06*(1-eta1)-.03*(1+eta1)+.29*(1+eta1)+.29*(1-eta1));
	  drdiy=((2*epack/mybl(3))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.02*(1-eta1)-.52*(1+eta1)+.4*(1+eta1)+.25*(1-eta1));
        
          
		
          inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
	  xy(3).dixdr(ix,iy)=Jac(1,1);
	  xy(3).dixdz(ix,iy)=Jac(1,2);
	  xy(3).diydr(ix,iy)=Jac(2,1);
	  xy(3).diydz(ix,iy)=Jac(2,2); 
      
      
      
	 		       	       
%end jacobian terms		       
		       	  	  
	  
	  xy(1).rz(2,ix,iy)= -.29*N1-.29*N2-.03*N3-.06*N4;
	  xy(1).rz(1,ix,iy)= .25*N1+.4*N2+.52*N3+.02*N4;
	  
% Inverse Jacobian terms
          
	  dzdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.29*(1-eta2)-.29*(1-eta2)-.03*(1+eta2)+.06*(1+eta2));
	  drdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.25*(1-eta2)+.4*(1-eta2)+.52*(1+eta2)-.02*(1+eta2));
	  dzdiy=((2*epack/mybl(3))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.29*(1-eta1)+.29*(1+eta1)-.03*(1+eta1)-.06*(1-eta1));
	  drdiy=((2*epack/mybl(3))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.25*(1-eta1)-.4*(1+eta1)+.52*(1+eta1)+.02*(1-eta1));
		
      inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
	  xy(1).dixdr(ix,iy)=Jac(1,1);
	  xy(1).dixdz(ix,iy)=Jac(1,2);
	  xy(1).diydr(ix,iy)=Jac(2,1);
	  xy(1).diydz(ix,iy)=Jac(2,2);	   					
	 		       	       
%end jacobian terms	
	  
        end
      end       
      
      for ix=1:mxbl(2)+1        
	for iy=1:mybl(2)+1
	  
	  if ((ix<=ixstart+mxbl(1)) & (ix>=ixstart)) ;
	      vx=(2.0*epack/mxbl(1))*((ix-3)-(mxbl(1)/2.0)- ...
                                                     (.02/.58)*mx);
	     vy=(2.0*epack/mybl(2))*((iy-1)-(mybl(2)/2.0));

	     eta1=log(vx+sqrt(1+vx^2))/log(epack+sqrt(1+epack^2));
	     eta2=log(vy+sqrt(1+vy^2))/log(epack+sqrt(1+epack^2));	     
	  N1=(1/4.0)*(1.0-eta1)*(1.0-eta2);
	  N2=(1/4.0)*(1.0+eta1)*(1.0-eta2);
	  N3=(1/4.0)*(1.0+eta1)*(1.0+eta2);
	  N4=(1/4.0)*(1.0-eta1)*(1.0+eta2);
	  xy(2).rz(1,ix,iy)=.02*N1+.52*N2+.52*N3+.02*N4;
	  xy(2).rz(2,ix,iy)= -.06*N1-.03*N2+.03*N3+.06*N4;
	  
	  % Inverse Jacobian terms
          
	  dzdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.06*(1-eta2)-.03*(1-eta2)+.03*(1+eta2)-.06*(1+eta2));
	  drdix=((2*epack/mxbl(3))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.02*(1-eta2)+.52*(1-eta2)+.52*(1+eta2)-.02*(1+eta2));
	  dzdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.06*(1-eta1)+.03*(1+eta1)+.03*(1+eta1)+.06*(1-eta1));
	  drdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.02*(1-eta1)-.52*(1+eta1)+.52*(1+eta1)+.02*(1-eta1));
		
          inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
	  xy(2).dixdr(ix,iy)=Jac(1,1);
	  xy(2).dixdz(ix,iy)=Jac(1,2);
	  xy(2).diydr(ix,iy)=Jac(2,1);
	  xy(2).diydz(ix,iy)=Jac(2,2);	   					
	 		       	       
         %end jacobian terms	
	  
	  
	  end
	  
	  if (ix<=ixstart)

             vx=(epack/(.01*mx/.58+1))*((ix-1)-(.01*mx/.58+1));
	     vy=(2.0*epack/mybl(2))*((iy-1)-(mybl(2)/2.0));

             eta1=log(vx+sqrt(1+vx^2))/log(epack+sqrt(1+epack^2));  
	     eta2=log(vy+sqrt(1+vy^2))/log(epack+sqrt(1+epack^2));
	     N1=(1/4.0)*(1.0-eta1)*(1.0-eta2);
	     N2=(1/4.0)*(1.0+eta1)*(1.0-eta2);
	     N3=(1/4.0)*(1.0+eta1)*(1.0+eta2);
	     N4=(1/4.0)*(1.0-eta1)*(1.0+eta2);
             xy(2).rz(1,ix,iy)=.02*N2+.02*N3;
             xy(2).rz(2,ix,iy)=-.06*N1-.06*N2+.06*N3+.06*N4;
	     
	  % Inverse Jacobian terms
          
	  dzdix=((epack/(.01*mx/.58))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.06*(1-eta2)-.06*(1-eta2)+.06*(1+eta2)-.06*(1+eta2));
	  drdix=((epack/(.01*mx/.58))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.02*(1-eta2)+.02*(1+eta2));
	  dzdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.06*(1-eta1)+.06*(1+eta1)+.06*(1+eta1)+.06*(1-eta1));
	  drdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.02*(1+eta1)+.02*(1+eta1));
		
          
          
        inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
      if (ix~=ixstart)
	  xy(2).dixdr(ix,iy)=Jac(1,1);
	  xy(2).dixdz(ix,iy)=Jac(1,2);
	  xy(2).diydr(ix,iy)=Jac(2,1);
	  xy(2).diydz(ix,iy)=Jac(2,2);
      else
      xy(2).dixdr(ix,iy)=(Jac(1,1)+xy(2).dixdr(ix,iy))/2;
	  xy(2).dixdz(ix,iy)=(Jac(1,2)+xy(2).dixdz(ix,iy))/2;
	  xy(2).diydr(ix,iy)=(Jac(2,1)+xy(2).diydr(ix,iy))/2;
	  xy(2).diydz(ix,iy)=(Jac(2,2)+xy(2).diydz(ix,iy))/2;
      end
	 		       	       
         %end jacobian terms	

	  elseif (ix>=ixstart+mxbl(1)) 
             vx=(epack/(.03*mx/.58))*((ix-3)-(.55*mx/.58));
             vy=(2.0*epack/mybl(2))*((iy-1)-(mybl(2)/2.0));

             eta1=log(vx+sqrt(1+vx^2))/log(epack+sqrt(1+epack^2)); 
	     eta2=log(vy+sqrt(1+vy^2))/log(epack+sqrt(1+epack^2));
	     N1=(1/4.0)*(1.0-eta1)*(1.0-eta2);
	     N2=(1/4.0)*(1.0+eta1)*(1.0-eta2);
	     N3=(1/4.0)*(1.0+eta1)*(1.0+eta2);
	     N4=(1/4.0)*(1.0-eta1)*(1.0+eta2);
	     xy(2).rz(1,ix,iy)=.52*N1+.58*N2+.58*N3+.52*N4;
	     xy(2).rz(2,ix,iy)=-.03*N1-.03*N2+.03*N3+.03*N4; 
	     
	  % Inverse Jacobian terms
          
	  dzdix=((epack/(.03*mx/.58))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.03*(1-eta2)-.03*(1-eta2)+.03*(1+eta2)-.03*(1+eta2));
	  drdix=((epack/(.03*mx/.58))/(sqrt(1+vx^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.52*(1-eta2)+.58*(1-eta2)+.58*(1+eta2)-.52*(1+eta2));
	  dzdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(.03*(1-eta1)+.03*(1+eta1)+.03*(1+eta1)+.03*(1-eta1));
	  drdiy=((2*epack/mybl(2))/(sqrt(1+vy^2)*log(epack+sqrt(1+epack^2))))* ...
	        (1/4)*(-.52*(1-eta1)-.58*(1+eta1)+.58*(1+eta1)+.52*(1-eta1));
		
        
        
        inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
      if (ix~=ixstart+mxbl(1))
	  xy(2).dixdr(ix,iy)=Jac(1,1);
	  xy(2).dixdz(ix,iy)=Jac(1,2);
	  xy(2).diydr(ix,iy)=Jac(2,1);
	  xy(2).diydz(ix,iy)=Jac(2,2);
  else
     xy(2).dixdr(ix,iy)=(Jac(1,1)+xy(2).dixdr(ix,iy))/2;
	  xy(2).dixdz(ix,iy)=(Jac(1,2)+xy(2).dixdz(ix,iy))/2;
	  xy(2).diydr(ix,iy)=(Jac(2,1)+xy(2).diydr(ix,iy))/2;
	  xy(2).diydz(ix,iy)=(Jac(2,2)+xy(2).diydz(ix,iy))/2;
  end	 		       	       
         %end jacobian terms	
	       

     end
       end
      end	
      
      
% grid seam points

   for ix=2:mxbl(1)
       dzdix=(xy(1).rz(2,ix+1,mybl(1)+1)-xy(1).rz(2,ix-1,mybl(1)+1))/2;
       dzdiy=(-xy(2).rz(2,ix+ixstart-1,3)+8*xy(2).rz(2,ix+ixstart-1,2)- ...
                  8*xy(1).rz(2,ix,mybl(1))+xy(1).rz(2,ix,mybl(1)-1))/12;
       drdix=(xy(1).rz(1,ix+1,mybl(1)+1)-xy(1).rz(1,ix-1,mybl(1)+1))/2;
       drdiy=(-xy(2).rz(1,ix+ixstart-1,3)+8*xy(2).rz(1,ix+ixstart-1,2)- ...
                  8*xy(1).rz(1,ix,mybl(1))+xy(1).rz(1,ix,mybl(1)-1))/12;
  

       inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
      
      xy(1).dixdr(ix,mybl(1)+1)=Jac(1,1);
	  xy(1).dixdz(ix,mybl(1)+1)=Jac(1,2);
	  xy(1).diydr(ix,mybl(1)+1)=Jac(2,1);
	  xy(1).diydz(ix,mybl(1)+1)=Jac(2,2);
      
       dzdix=(xy(3).rz(2,ix+1,1)-xy(3).rz(2,ix-1,1))/2;
       dzdiy=(-xy(3).rz(2,ix,3)+8*xy(3).rz(2,ix,2)- ...
                8*xy(2).rz(2,ix+ixstart-1,mybl(2))+xy(2).rz(2,ix+ixstart-1,mybl(2)-1))/12;
       drdix=(xy(3).rz(1,ix+1,1)-xy(3).rz(1,ix-1,1))/2;
       drdiy=(-xy(3).rz(1,ix,3)+8*xy(3).rz(1,ix,2)- ...
                8*xy(2).rz(1,ix+ixstart-1,mybl(2))+xy(2).rz(1,ix+ixstart-1,mybl(2)-1))/12;
  

       inv_Jac=[drdix drdiy
	           dzdix dzdiy];
	  Jac=inv(inv_Jac);
      
      xy(3).dixdr(ix,1)=Jac(1,1);
	  xy(3).dixdz(ix,1)=Jac(1,2);
	  xy(3).diydr(ix,1)=Jac(2,1);
	  xy(3).diydz(ix,1)=Jac(2,2);
  end

      
      
% second order jacobian terms

	
    for ibl=[1 3 2]
	
	for ix=2:mxbl(ibl)
	  for iy=1:mybl(ibl)+1
                        

    	if (ibl==1 & iy==mybl(ibl)+1)
%	xy(ibl).dixdz(ix,iy)=(xy(ibl).dixdz(ix,iy)+...
%	                     xy(2).dixdz(ix+ixstart-1,1))/2;
%	xy(ibl).diydz(ix,iy)=(xy(ibl).diydz(ix,iy)+...
%	                     xy(2).diydz(ix+ixstart-1,1))/2;
%	xy(ibl).dixdr(ix,iy)=(xy(ibl).dixdr(ix,iy)+...
%	                     xy(2).dixdr(ix+ixstart-1,1))/2;
%	xy(ibl).diydr(ix,iy)=(xy(ibl).diydr(ix,iy)+...
%	                     xy(2).diydr(ix+ixstart-1,1))/2;	     			       
	   
	xy(ibl).dixdrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).dixdr(ix+1,iy)-xy(ibl).dixdr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(2).dixdr(ix+ixstart-1,2)-xy(ibl).dixdr(ix,iy-1))/2);
	xy(ibl).dixdzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).dixdz(ix+1,iy)-xy(ibl).dixdz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(2).dixdz(ix+ixstart-1,2)-xy(ibl).dixdz(ix,iy-1))/2);	    	        
	xy(ibl).diydrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).diydr(ix+1,iy)-xy(ibl).diydr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(2).diydr(ix+ixstart-1,2)-xy(ibl).diydr(ix,iy-1))/2);
	xy(ibl).diydzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).diydz(ix+1,iy)-xy(ibl).diydz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(2).diydz(ix+ixstart-1,2)-xy(ibl).diydz(ix,iy-1))/2);   	     
	end     
	if (ibl==3 & iy==1)
%	xy(ibl).dixdz(ix,iy)=(xy(ibl).dixdz(ix,iy)+...
%	                     xy(2).dixdz(ix+ixstart-1,mybl(2)+1))/2;
%	xy(ibl).diydz(ix,iy)=(xy(ibl).diydz(ix,iy)+...
%	                     xy(2).diydz(ix+ixstart-1,mybl(2)+1))/2;
%	xy(ibl).dixdr(ix,iy)=(xy(ibl).dixdr(ix,iy)+...
%	                     xy(2).dixdr(ix+ixstart-1,mybl(2)+1))/2;
%	xy(ibl).diydr(ix,iy)=(xy(ibl).diydr(ix,iy)+...
%	                     xy(2).diydr(ix+ixstart-1,mybl(2)+1))/2;	     			       
	   
	xy(ibl).dixdrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).dixdr(ix+1,iy)-xy(ibl).dixdr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(ibl).dixdr(ix,iy+1)-xy(2).dixdr(ix+ixstart-1,mybl(2)))/2);
	xy(ibl).dixdzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).dixdz(ix+1,iy)-xy(ibl).dixdz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(ibl).dixdz(ix,iy+1)-xy(2).dixdz(ix+ixstart-1,mybl(2)))/2);	    	        
	xy(ibl).diydrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).diydr(ix+1,iy)-xy(ibl).diydr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(ibl).diydr(ix,iy+1)-xy(2).diydr(ix+ixstart-1,mybl(2)))/2);
	xy(ibl).diydzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).diydz(ix+1,iy)-xy(ibl).diydz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(ibl).diydz(ix,iy+1)-xy(2).diydz(ix+ixstart-1,mybl(2)))/2);   	     
	end     	 
    
    if (ibl==2 & iy==1 & ix>ixstart & ix<ixstart+mxbl(1))
    xy(ibl).dixdz(ix,iy)=xy(1).dixdz(ix-ixstart+1,mybl(1)+1);
	xy(ibl).diydz(ix,iy)=xy(1).diydz(ix-ixstart+1,mybl(1)+1);
	xy(ibl).dixdr(ix,iy)=xy(1).dixdr(ix-ixstart+1,mybl(1)+1);
	xy(ibl).diydr(ix,iy)=xy(1).diydr(ix-ixstart+1,mybl(1)+1);
    
       if (ix~=ixstart+mxbl(1))
    xy(ibl).dixdzz(ix,iy)=xy(1).dixdzz(ix-ixstart+1,mybl(1)+1);
	xy(ibl).diydzz(ix,iy)=xy(1).diydzz(ix-ixstart+1,mybl(1)+1);
	xy(ibl).dixdrr(ix,iy)=xy(1).dixdrr(ix-ixstart+1,mybl(1)+1);
	xy(ibl).diydrr(ix,iy)=xy(1).diydrr(ix-ixstart+1,mybl(1)+1);
      end
    end
    if (ibl==2 & iy==mybl(2)+1 & ix>ixstart & ix<ixstart+mxbl(1))
    xy(ibl).dixdz(ix,iy)=xy(3).dixdz(ix-ixstart+1,1);
	xy(ibl).diydz(ix,iy)=xy(3).diydz(ix-ixstart+1,1);
	xy(ibl).dixdr(ix,iy)=xy(3).dixdr(ix-ixstart+1,1);
	xy(ibl).diydr(ix,iy)=xy(3).diydr(ix-ixstart+1,1);
    
      if (ix~=ixstart+mxbl(1))
    xy(ibl).dixdzz(ix,iy)=xy(3).dixdzz(ix-ixstart+1,1);
	xy(ibl).diydzz(ix,iy)=xy(3).diydzz(ix-ixstart+1,1);
	xy(ibl).dixdrr(ix,iy)=xy(3).dixdrr(ix-ixstart+1,1);
	xy(ibl).diydrr(ix,iy)=xy(3).diydrr(ix-ixstart+1,1);
      end
    end
          
    if (iy~=1 & iy~=mybl(ibl)+1)
	
	xy(ibl).dixdrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).dixdr(ix+1,iy)-xy(ibl).dixdr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(ibl).dixdr(ix,iy+1)-xy(ibl).dixdr(ix,iy-1))/2);
	xy(ibl).dixdzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).dixdz(ix+1,iy)-xy(ibl).dixdz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(ibl).dixdz(ix,iy+1)-xy(ibl).dixdz(ix,iy-1))/2);	    	        
	xy(ibl).diydrr(ix,iy)=xy(ibl).dixdr(ix,iy)*...
	     ((xy(ibl).diydr(ix+1,iy)-xy(ibl).diydr(ix-1,iy))/2)+...
	                      xy(ibl).diydr(ix,iy)*...
	     ((xy(ibl).diydr(ix,iy+1)-xy(ibl).diydr(ix,iy-1))/2);
	xy(ibl).diydzz(ix,iy)=xy(ibl).dixdz(ix,iy)*...
	     ((xy(ibl).diydz(ix+1,iy)-xy(ibl).diydz(ix-1,iy))/2)+...
	                      xy(ibl).diydz(ix,iy)*...
	     ((xy(ibl).diydz(ix,iy+1)-xy(ibl).diydz(ix,iy-1))/2);
	end  

    
   
	
    end
   end
    
 end