      program plate
      implicit real*8(a-h,o-z)
      parameter(ist=4500000)
c*****************************************************************
c      Mindlin plate analysis under a movin mass;
c*****************************************************************
      dimension ma(ist)
      data nbd,ir,lwe,lss,lsswsp/190,2,4,3,2/
      
      nbd1=nbd+1
      open(8,file='results.dat',status='unknown')
c
c     nbd   - bandwith of the stiffness matrix
c     ir=2  - length of Real/length of integer:  R*8/I*4
c
      do 1 i=1,ist
1       ma(i)=0

c     reading initial data .....................................
      open(1,file='plate.cmm',status='old')
      read(1,*)nx,ny,xlx,yly,dt,alfa,vjazdy,qtotal,gamma,xmasa
      read(1,*)e,xni,ro,t,xkz
      close(1)

      np=(nx+1)*(ny+1)
      ne=nx*ny
      nwb=(nx+ny)*2
      nszf=np*lss

c     nx, ny   - No. of spatial elements in x and y direction
c     nwb      - No. of fixed nodes
c     dt,h     - time step
c     e        - Young modulus
c     xni      - Poissone coefficient
c     ro       - mass density
c     alfa     - parametr (0,1] procedury calkowania rown. ruchu
c     vjazdy   - velocity of the moving mass
c     gamma    - coefficient of numerical damping
c     kz       - stiffness of foundation


c     memory partition............................................
c
c     1    - wpw
c     n1   - npe
c     n2   - ityp
c     n3   - iwb
c     n4   - r                (right-had side vector)
c     n5   - 
c     n6   - ef
c     n7   - q                (velocity vector)
c     n8   - u                (displacement vector)
c     n13  - skd | ind        (stiffness matrix) | temp. matrix
c     n14  - xmater
c     n15  - ind

      n1= 1+np*lsswsp*ir
      n2=n1+lwe*ne
      n3=n2+ne

      call dane(ma(1),ma(n1),ma(n3),ma(n2),np,ne,nx,ny,
     x     xlx,yly,dx,dy,nwbakt,lwe,lss,nbd00)

      nwb=nwbakt

      if(nbd00.gt.nbd)then
         write(*,*)'increase nbd ',nbd,' -->',nbd00
         stop
      else
         write(*,*)'decrease nbd ',nbd,' -->',nbd00
      endif

      n4=n3+2*nwb
      n6=n4+nszf*ir
      n7=n6+nszf*ir
      n8=n7+2*nszf*ir
      n13=n8+nszf*ir
      n14=n13+nszf*nbd1*ir
      n15=n14+20*ir
      n16=n15+ne
      n20=n16

      write(*,710)n16,ist
 710  format(' mem/ist=',i10,'/',i10)
      if(n16.gt.ist)stop 'lack of memory'
c------------------------------------------------------------------
      iq1=1
      iq2=2
      it=0

100   it=it+1
      czas=it*dt
      write(*,*)'step, time=',it,czas      
      call load(ma(1),ma(n4),np,lss,
     x     czas,vjazdy,qtotal,nx,ny,xlx,yly,dt,alfa,xmasa)
      
      call loadnodal(ma(n6),ma(n4),np,lss)

      call formk(ma(n13),ma(1),ma(n1),np,ne,dx,dy,ma(n4),lwe,lss,
     1     nszf,nbd,nbd1,dt,alfa,ma(n6),ma(n7),iq1,iq2,ma(n2),
     2     xmasa,ximasa,vjazdy,ma(n8),xlx,yly,
     3      e,xni,ro,t,xkz)
     
      call modyfwb(ma(n13),nbd,nbd1,np,lss,ma(n4),nwb,ma(n3))

      call dpbsv('U',nszf,nbd,1,ma(n13),nbd1,ma(n4),nszf,info)
      if(info.ne.0)stop 'singular matrix'

      call przemieszczenia(ma(n4),ma(n7),ma(n8),iq1,iq2,alfa,dt,
     x     gamma,np,lss)

      call silywezlowe(ma(1),ma(n1),np,ne,dx,dy,ma(n4),lwe,lss,
     1     nszf,nbd,nbd1,dt,alfa,ma(n6),ma(n2),ma(n8),ma(n14),
     2     xlx,yly,xmasa,vjazdy,e,xni,ro,t,xkz)

      call wydruki(ma(1),ma(n1),ma(n8),ma(n7),iq1,iq2,np,ne,lss,
     x     nx,ny,czas,vjazdy,xlx,wspx,e,xi)

      iq1=3-iq1
      iq2=3-iq2
      if(wspx.lt.1.2)go to 100
      end


      subroutine dane(wpw,npe,iwb,ityp,np,ne,nx,ny,xlx,yly,
     x     dx,dy,nwb,lwe,lss,nbd00)
      implicit real*8(a-h,o-z)
      dimension wpw(2,1),npe(lwe,1),iwb(2,1),ityp(1)

      dx=xlx/nx
      dy=yly/ny

      l=0
      do 9 i=0,nx-1
         do 9 j=0,ny-1
            l=l+1
            npe(1,l)=i*(ny+1)+j+1
            npe(2,l)=i*(ny+1)+j+2
            npe(3,l)=(i+1)*(ny+1)+j+1
            npe(4,l)=(i+1)*(ny+1)+j+2
            ityp(l)=1
 9    continue

      nwb=0
      do 10 i=1,ny+1
          nwb=nwb+1
          iwb(1,nwb)=i
          iwb(2,nwb)=100
          nwb=nwb+1
          iwb(1,nwb)=(nx)*(ny+1)+i
10        iwb(2,nwb)=100
      
      do 11 i=2,nx
          nwb=nwb+1
          iwb(1,nwb)=(ny+1)*(i-1)+1
          iwb(2,nwb)=100
          nwb=nwb+1
          iwb(1,nwb)=(ny+1)*(i)
11        iwb(2,nwb)=100
      
      ll=0
      do 12 i=0,nx
         do 12 j=0,ny
            ll=ll+1
            wpw(1,ll)=i*dx
            wpw(2,ll)=j*dy
 12   continue

      nbd00=0
      do 91 i=1,ne
         n1=iabs(npe(1,i)-npe(2,i))
         nbd00=max0(nbd00,n1)
         if(npe(3,i).gt.0)then
            n1=iabs(npe(2,i)-npe(3,i))
            nbd00=max0(nbd00,n1)
            n1=iabs(npe(3,i)-npe(4,i))
            nbd00=max0(nbd00,n1)
            n1=iabs(npe(4,i)-npe(1,i))
            nbd00=max0(nbd00,n1)
         endif
 91   continue

      nbd00=(nbd00+1)*lss
      write(*,*)'bandwidth =',nbd00

      do 71 i=1,nwb
         if(iwb(1,i).gt.np)stop 7
 71   continue

      end
      
      

      subroutine FORMK(skd,WPW,NPE,np,ne,dx,dy,r1,lwe,lss,
     1     nszf,nb,nbd1,h,alfa,ef,q,iq1,iq2,ityp,
     2     xmasa,ximasa,vm,u,xlx,yly,e,xni,ro,t,xkz)
      implicit real*8(a-h,o-z)

      dimension skd(nbd1,1),wpw(2,1),r1(6),dd(12,24),npe(lwe,1),
     1     q(2,1),wsp(2,4),nel(4),ef(1),ityp(1),u(1)
      common /wezly/x0,ktory_wez1,ktory_wez2,xkappa,sk1,sk2,ine

      lsh=lwe*lss

      do 300 j=1,nszf
         do 300 i=1,nbd1
 300        skd(i,j)=0d0


      DO 400 N=1,NE
         do 81 ii=1,12
            do 81 jj=1,24
 81            dd(ii,jj)=0d0

      nel(1)=npe(1,n)
      nel(2)=npe(2,n)
      nel(3)=npe(3,n)
      nel(4)=npe(4,n)
c----------------------------------------------
      call plate_elem(dd,dx,dy,h,alfa,1,e,xni,ro,t,xkz)
      lwelok=4

      do 7 i=1,lwelok
         nw1=nel(i)
         do 7 ii=1,lss
            il=nw1*lss-lss+ii
            ix=i*lss-lss+ii
            
            do 7 j=1,lwelok
               nw2=nel(j)
               do 7 jj=1,lss
                  ik=nw2*lss-lss+jj
                  iy=j*lss-lss+jj
                  r1(il)=r1(il)-dd(ix,iy)*q(iq1,nw2*lss-lss+jj)
                  if(il.le.ik)skd(nbd1+il-ik,ik)=
     x                 skd(nbd1+il-ik,ik)+dd(ix,iy+lsh)
 7             continue
 400   CONTINUE

      nel(1)=ktory_wez1
      nel(2)=ktory_wez2
      lsh=6

      call mass(dd,dx,dy,h,alfa,xmasa,vm)

      lwelok=2
      do 71 i=1,lwelok
         nw1=nel(i)
         do 71 ii=1,lss
            il=nw1*lss-lss+ii
            ix=i*lss-lss+ii
            
            do 71 j=1,lwelok
               nw2=nel(j)
               do 71 jj=1,lss
                  
                  ik=nw2*lss-lss+jj
                  iy=j*lss-lss+jj
                  r1(il)=r1(il)-dd(ix,iy)*q(iq1,nw2*lss-lss+jj)
                  if(il.le.ik)skd(nbd1+il-ik,ik)=
     x                 skd(nbd1+il-ik,ik)+dd(ix,iy+lsh)
 71            continue
       end
       

      subroutine modyfwb(b,nbd,nbd1,np,lss,r,nwb,iwb)

      implicit real*8(a-h,o-z)
      dimension b(nbd1,1),r(6),iwb(2,1)

      do 1 iii=1,nwb
         iwezel=iwb(1,iii)
         kod=iwb(2,iii)
         do 2 jj=1,lss
            icyfra=int(kod/10**(lss-jj))
            kod=kod-icyfra*10**(lss-jj)

            if(icyfra.eq.1)then
               nrl=iwezel*lss-lss+jj
               al=0d0

               do 400 ii=1,nbd1
                  lin=nrl-nbd1+ii
                  if(lin.gt.np*lss)stop 2
                  if(ii.gt.np*lss)stop 3
                  r(lin)=r(lin)-al*b(ii,nrl)
                  b(ii,nrl)=0d0
 400           continue
                     
               do 401 ii=1,min0(nbd,np*lss-nrl)
                  r(nrl+ii)=r(nrl+ii)-al*b(nbd1-ii,nrl+ii)
                  b(nbd1-ii,nrl+ii)=0d0                  
 401           continue

               b(nbd1,nrl)=1.d0
               r(nrl)=al
            endif

 2       continue
 1    continue

      end



      subroutine load(wpw,fzewn,np,lss,
     x     czas,vjazdy,qtotal,nx,ny,xlx,yly,h,alfa,xm0)
      implicit real*8(a-h,o-z)

      dimension wpw(2,1),fzewn(1)
      common /wezly/x0,ktory_wez1,ktory_wez2,xkappa,sk1,sk2,ine
      common sw,ktory_wez1m,ktory_wez2m

c     zeroing of external load vector

      do 2 i=1,np*lss
 2       fzewn(i)=0d0

      dx=xlx/nx
      dy=yly/ny

c     nodes which are subjected by a load: ktory_wez1, ktory_wez2
      x=vjazdy*czas
      if(x.le.xlx)then
         ine=int(czas*vjazdy/dx-1d-10)+1
         x0=(czas-h)*vjazdy-ine*dx+dx
         xkappa=(x0+vjazdy*alfa*h)/dx

         ktory_wez1=(ny+1)*(ine-1)+ny/2+1
         ktory_wez2=(ny+1)*(ine  )+ny/2+1

         sk1=1d0-xkappa
         sk2=xkappa

         fzewn(ktory_wez1*3-2)=-qtotal*sk1           
         fzewn(ktory_wez2*3-2)=-qtotal*sk2
      nnn=(ny+1)*(nx/2+1)-ny/2

      else
         w1=0.
         w2=0.
         ktory_wez1=1
         ktory_wez2=1
      endif
      
      fzewn(3*ktory_wez1-2)=fzewn(3*ktory_wez1-2)-sw*(1-xkappa)
      fzewn(3*ktory_wez2-2)=fzewn(3*ktory_wez2-2)-sw*xkappa
      end
      
      
      
      subroutine loadnodal(ef,r,np,lss)
      implicit real*8(a-h,o-z)
      dimension ef(1),r(1)

c     computation of right-hand side vector

      do 1 i=1,np*lss
 1       r(i)=r(i)-ef(i)
      end
      

      
      subroutine przemieszczenia(r,q,u,iq1,iq2,alfa,dt,gamma,np,lss)
      implicit real*8(a-h,o-z)
      dimension r(1),q(2,1),u(1)

      beta=1d0-alfa/(1d0+gamma)

      do 9 i=1,np*lss
         q(iq2,i)=r(i)
         u(i)=u(i)+(q(iq1,i)*(1d0-beta)+q(iq2,i)*beta)*dt
 9    continue
      end
      
      
      subroutine silywezlowe(WPW,NPE,np,ne,dx,dy,r1,lwe,lss,
     1     nszf,nbd,nbd1,h,alfa,ef,ityp,u,xmater,xlx,yly,
     2     xxm0,vm,e,xni,ro,t,xkz)
      implicit real*8(a-h,o-z)

      dimension wpw(2,1),r1(6),dd(12,24),npe(lwe,1),
     2  wsp(2,4),nel(4),ef(1),ityp(1),u(1),xmater(20)
      common /wezly/x0,ktory_wez1,ktory_wez2,xkappa,sk1,sk2,ine
      common sw,ktory_wez1m,ktory_wez2m

      lsh=lwe*lss
      itryb=2

      do 17 i=1,nszf
 17      ef(i)=0d0

      DO 400 N=1,NE
         nel(1)=npe(1,n)
         nel(2)=npe(2,n)
         nel(3)=npe(3,n)
         nel(4)=npe(4,n)
         
c     element matrix..........................................
      call plate_elem(dd,dx,dy,h,alfa,itryb,e,xni,ro,t,xkz)
         lwelok=4
         do 7 i=1,lwelok
            nw1=nel(i)
            do 7 ii=1,lss
               il=nw1*lss-lss+ii
               ix=i*lss-lss+ii

               do 7 j=1,lwelok
                  nw2=nel(j)
                  do 7 jj=1,lss

                     ik=nw2*lss-lss+jj
                     iy=j*lss-lss+jj
                     ef(il)=ef(il)+dd(ix,iy)*u(nw2*lss-lss+jj)
7        continue
400   CONTINUE

c.............. nodal forces of the mass ..................
      u1=u(3*ktory_wez1m-2)
      u2=u(3*ktory_wez2m-2)

      u3=u(3*ktory_wez1-2)
      u4=u(3*ktory_wez2-2)
      ktory_wez1m=ktory_wez1
      ktory_wez2m=ktory_wez2

      sw=xxm0*vm/dx/h*(u4-u3-u2+u1)

      end
       
       
       
      subroutine wydruki(wpw,npe,u,q,iq1,iq2,np,ne,lss,nx,ny,czas,
     x     vjazdy,xlx,wspx,e,xi)
      implicit real*8(a-h,o-z)
      character*24 dysk
      dimension wpw(2,1),npe(4,1),u(1),q(2,1)
      common /wezly/x0,ktory_wez1,ktory_wez2,xkappa,sk1,sk2,ine

      wspx=czas*vjazdy/xlx
      dx=xlx/nx

      nnn=(ny+1)*(nx/2+1)-ny/2
      w_sle=u(ktory_wez1*3-2)*sk1+u(ktory_wez2*3-2)*sk2
      write(8,991)czas*vjazdy/xlx,w_sle,u(nnn*3-2)
991      format(1pe15.6,1pe15.6,1pe15.6)

      end
      
      
      subroutine plate_elem(dd,dxbok,dybok,h,alfa,itryb,e,xni,ro,t,xkz)
      implicit real*8(a-h,o-z)
      dimension iksi(4),ieta(4)
      dimension dd(12,24)
      data iksi,ieta/1,-1,1,-1,-1,-1,1,1/
c...........................................
      xkappa=5./6.
      g=e/2./(1+xni)
      a=dxbok/2.
      b=dybok/2.
c...........................................
      wsp1=e*t**3/48./a/b/(1.-xni**2)
      wsp2=g*xkappa*t

      do 711 i=1,4
         do 711 j=1,4
            dd(3*i-2,3*j-2)=0d0
            dd(3*i-2,3*j-1)=0d0
            dd(3*i-2,3*j)=0d0
            dd(3*i-1,3*j-2)=0d0
            dd(3*i,3*j-2)=0d0
      
            dd(3*i-1,3*j-1)=((2*a**2*ieta(i)*ieta(j)*
     x      (iksi(i)*iksi(j)+3)+b**2*iksi(i)*iksi(j)*
     x      (1-xni)*(ieta(i)*ieta(j)+3))/6)*wsp1
            dd(3*i-1,3*j)=(-a*b*(2*ieta(i)*xni*iksi(j)+ieta(j)*
     x      iksi(i)*(1-xni))/2)*wsp1
            dd(3*i,3*j-1)=(a*b*(ieta(i)*iksi(j)*(xni-1)-2*
     x      ieta(j)*xni*iksi(i))/2)*wsp1
            dd(3*i,3*j)=((2*b**2*iksi(i)*iksi(j)*
     x      (ieta(i)*ieta(j)+3)-a**2*ieta(i)*ieta(j)*(xni-1)*
     x      (iksi(i)*iksi(j)+3))/6)*wsp1
c----------------------
            dd(3*i-2,3*j-2)=dd(3*i-2,3*j-2)+((a**2*ieta(i)*ieta(j)*
     x      (iksi(i)*iksi(j)+3)+b**2*iksi(i)*iksi(j)*
     x      (ieta(i)*ieta(j)+3))/(12*a*b))*wsp2
            dd(3*i-2,3*j-1)=dd(3*i-2,3*j-1)+(-a*ieta(i)*
     x      (iksi(i)*iksi(j)+3)/12)*wsp2
            dd(3*i-2,3*j)=dd(3*i-2,3*j)+(b*iksi(i)*
     x      (ieta(i)*ieta(j)+3)/12)*wsp2
            dd(3*i-1,3*j-2)=dd(3*i-1,3*j-2)+(-a*ieta(j)*(iksi(i)*
     x      iksi(j)+3)/12)*wsp2
            dd(3*i-1,3*j-1)=dd(3*i-1,3*j-1)+
     x      (a*b*(ieta(i)*ieta(j)+3)*(iksi(i)*iksi(j)+3)/36)*wsp2
            dd(3*i-1,3*j)=dd(3*i-1,3*j)+0d0
            dd(3*i,3*j-2)=dd(3*i,3*j-2)+(b*iksi(j)*
     x      (ieta(i)*ieta(j)+3)/12)*wsp2
            dd(3*i,3*j-1)=dd(3*i,3*j-1)+0d0
            dd(3*i,3*j)=dd(3*i,3*j)+(a*b*(ieta(i)*ieta(j)+3)*
     x      (iksi(i)*iksi(j)+3)/36)*wsp2
711      continue

c........................... elastic foundation
      do 73 i=1,12,3
73      dd(i,i)=dd(i,i)+xkz*dxbok*dybok/4.

c........................... further dynamic stiffness matrix
      if(itryb.eq.2)return

c     formulation of entire matrix [A|B] ...................

      xmnoz=h*alfa**2/2.d0
      do 16 i=1,12
         do 16 j=13,24
 16         dd(i,j)=dd(i,j-12)*xmnoz

      xmnoz=(1.d0-alfa/2.d0)*alfa*h
      do 17 i=1,12
         do 17 j=1,12
 17         dd(i,j)=dd(i,j)*xmnoz

c     addition of the mass matrix ............ [-1/h*...|+1/h*...]

      xmasa=dxbok*dybok*t*ro/4./h/108.

      do 71 ii=1,4
         do 71 jj=1,4
            if(ii.eq.jj)mnoz=48
            if(iabs(ii-jj).eq.1)mnoz=24
            if(iabs(ii-jj).eq.2)mnoz=12
            if(iabs(ii-jj).eq.3)mnoz=24
            dd(3*ii-2,   3*jj-2)=dd(3*ii-2,   3*jj-2)-xmasa*mnoz
            dd(3*ii-2,12+3*jj-2)=dd(3*ii-2,12+3*jj-2)+xmasa*mnoz
 71   continue

      xmasa=dxbok*dybok*t**3*ro/4./h/108.

      do 72 ii=1,4
         do 72 jj=1,4
            if(ii.eq.jj)mnoz=4
            if(iabs(ii-jj).eq.1)mnoz=2
            if(iabs(ii-jj).eq.2)mnoz=1
            if(iabs(ii-jj).eq.3)mnoz=2

            dd(3*ii-1,   3*jj-1)=dd(3*ii-1,   3*jj-1)-xmasa*mnoz
            dd(3*ii  ,   3*jj  )=dd(3*ii  ,   3*jj  )-xmasa*mnoz
            dd(3*ii-1,12+3*jj-1)=dd(3*ii-1,12+3*jj-1)+xmasa*mnoz
            dd(3*ii  ,12+3*jj  )=dd(3*ii  ,12+3*jj  )+xmasa*mnoz
 72   continue
      end


      subroutine mass(dd,dxbok,dybok,h,alfa,xxm0,vm)
      implicit real*8(a-h,o-z)
      dimension dd(12,24),xm(6,12),xc(6,12),xk(6,12)
      common /wezly/x0,ktory_wez1,ktory_wez2,xkappa,sk1,sk2,ine
      
c     formulation of matrices of mass element

      data xm,xc,xk/72*0d0,72*0d0,72*0d0/

      do 11 i=1,12
         do 11 j=1,24
 11         dd(i,j)=0d0

      alfa2=0.
c...........................formulation of matrix M
      xm(1,1)=-xxm0/h*(1.-xkappa)**2
      xm(1,4)=-xxm0/h*xkappa*(1.-xkappa)
      xm(4,1)=-xxm0/h*xkappa*(1.-xkappa)
      xm(4,4)=-xxm0/h*xkappa**2

      xm(1,7)= xxm0/h*(1.-xkappa)**2
      xm(1,10)=xxm0/h*xkappa*(1.-xkappa)
      xm(4,7)= xxm0/h*xkappa*(1.-xkappa)
      xm(4,10)=xxm0/h*xkappa**2

c...........................formulation of matrix C
      w=2d0*xxm0*vm/dxbok

      xc(1,1)=-w*(1-xkappa)*(1-alfa2)
      xc(1,4)= w*(1-xkappa)*(1-alfa2)
      xc(4,1)=-w*xkappa*(1-alfa2)
      xc(4,4)= w*xkappa*(1-alfa2)

      xc(1,7)= -w*(1-xkappa)*alfa2
      xc(1,10)= w*(1-xkappa)*alfa2
      xc(4,7)= -w*xkappa*alfa2
      xc(4,10)= w*xkappa*alfa2

      do 1 i=1,6
      do 1 j=1,12
1        dd(i,j)=xm(i,j)+xc(i,j)

      end

