dimension v(1000,1000),d(4,4),dt(4,4),c(1000,1000)
    *,x(1000,1000),z(1000,1000),zt(1000,1000)
    *,ztdt(1000,1000),dx(1000,1000)
    *,ztdtc(1000,1000),ztdtcdx(1,1)
       real v,d,dt,c,x,z,zt,vv,dx,ddx,ddz,h,l
    *,xxx,zzz,xxxx,zzzz
    h=1
           L=1
    read*,xxxx,zzzz
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc节点的速度
          do i=1,5
       do j=1,5
           v(i,j)=1600
       enddo
    enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc插值：求任意一点的速度
if(xxxx.ge.1.and.xxxx.le.2.and.zzzz.ge.1.and.zzzz.le.2)then
                 i=int(xxxx/h)+1
          j=int(zzzz/l)+1
          xxx=(i-1)*h
          zzz=(j-1)*l
           d(1,1)=-1
            d(1,2)=2
           d(1,3)=-1
           d(1,4)=0
           d(2,1)=3
           d(2,2)=-5
           d(2,3)=0
           d(2,4)=2
           d(3,1)=-3
           d(3,2)=4
           d(3,3)=1
           d(3,4)=0
           d(4,1)=1
           d(4,2)=-1
           d(4,3)=0
           d(4,4)=0
           dt(1,1)=-1
            dt(1,2)=3
           dt(1,3)=-3
           dt(1,4)=1
           dt(2,1)=2
           dt(2,2)=-5
           dt(2,3)=4
           dt(2,4)=-1
           dt(3,1)=-1
           dt(3,2)=0
           dt(3,3)=1
            dt(3,4)=0
           dt(4,1)=0
           dt(4,2)=2
           dt(4,3)=0
           dt(4,4)=0                              
          c(1,1)=v(i-1,j-1)           
          c(1,2)=v(i,j-1)
          c(1,3)=v(i+1,j-1)
           c(1,4)=v(i+2,j-1)
          c(2,1)=v(i-1,j)
          c(2,2)=v(i,j)
          c(2,3)=v(i+1,j)
            c(2,4)=v(i+2,j)
          c(3,1)=v(i-1,j+1)
          c(3,2)=v(i,j+1)
          c(3,3)=v(i+1,j+1)
          c(3,4)=v(i+2,j+1)
          c(4,1)=v(i-1,j+2)
          c(4,2)=v(i,j+2)
          c(4,3)=v(i+1,j+2)  
          c(4,4)=v(i+2,j+2)                
          ddx=(xxxx-xxx)/h
                 ddz=(zzzz-zzz)/l
                 x(1,1)=ddx**3
          x(2,1)=ddx**2
            x(3,1)=ddx**1
                 x(4,1)=1
          z(1,1)=ddz**3
          z(2,1)=ddz**2
                 z(3,1)=ddz**1
                 z(4,1)=1
          zt(1,1)=z(1,1)
          zt(1,2)=z(2,1)
          zt(1,3)=z(3,1)
          zt(1,4)=z(4,1)    
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 两个矩阵相乘公式，ii为行，jj为列。       
c    ztdt=zt*dt
      do ii=1,1
        do jj=1,4
         ztdt(ii,jj)=0
          do k=1,4
            ztdt(ii,jj)=ztdt(ii,jj)+zt(ii,k)*dt(k,jj)
        enddo
      enddo
    enddo
c      dx=d*x    
      do ii=1,4
        do jj=1,1
         dx(ii,jj)=0
        do k=1,4
            dx(ii,jj)=dx(ii,jj)+d(ii,k)*x(k,jj)
        enddo
      enddo
    enddo
c    ztdtc=zt*dt*c=ztdt*c
      do ii=1,1
        do jj=1,4
         ztdtc(ii,jj)=0
        do k=1,4
            ztdtc(ii,jj)=ztdtc(ii,jj)+ztdt(ii,k)*c(k,jj)
        enddo
      enddo
    enddo
c     ztdtcdx=zt*dt*c*dx=ztdtc*dx    
      do ii=1,1
        do jj=1,1
         ztdtcdx(ii,jj)=0
        do k=1,4
            ztdtcdx(ii,jj)=ztdtcdx(ii,jj)+ztdtc(ii,k)*dx(k,jj)
           vv=1/4*ztdtcdx(ii,jj)           
        enddo
      enddo
    enddo
else
     vv=8888    
endif    
print*,vv
end