module defpi
implicit none
real(8),parameter::pi=3.14159265358979323846D0
end module defpi

program f90_navier
    use defpi
    implicit none
    integer,parameter::num=1001
    integer::mmax,nmax
    real(8),parameter::E=200000.0D0
    real(8),parameter::po=0.3D0
    real(8)::D,h,q,a,b
    real(8)::x,y,coef,qmn,work
    real(8)::w,Mxx,Myy,Mxy,Qx,Qy,Vx,Vy,R
    integer::i,m,n
    real(8)::array(1:num,1:num*10)

    mmax=num
    h=0.1
    D=E*h*h*h/12.0D0/(1.0D0-po*po)
    a=1.0D0
    q=1.0D0

    write(6,'(a8,a8,2(a8),5(a8))') 'b/a','w','Mx','My','Qx','Qy','Vx','Vy','R'
    do i=0,9
        b=a+a*dble(i)*0.1D0
        nmax=int(b/a*dble(mmax)/2.0D0)*2+1
        call LOAD(num,mmax,nmax,D,q,a,b,array)
        x=0.5D0*a
        y=0.5D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,w,Mxx,Myy,work,work,work,work,work)
        x=0.0D0*a
        y=0.5D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,work,Qx,work,Vx,work)
        x=0.5D0*a
        y=0.0D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,work,work,Qy,work,Vy)
        x=1.0D0*a
        y=1.0D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,Mxy,work,work,work,work)
        R=2.0D0*Mxy
        write(6,'(f8.3,f8.5,2(f8.4),5(f8.3))') b/a,w*D/q/a**4.0D0,Mxx/q/a/a,Myy/q/a/a,Qx/q/a,Qy/q/a,Vx/q/a,Vy/q/a,R/q/a/a
    end do

    do i=2,10
        b=a*dble(i)
        nmax=int(b/a*dble(mmax))
        call LOAD(num,mmax,nmax,D,q,a,b,array)
        x=0.5D0*a
        y=0.5D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,w,Mxx,Myy,work,work,work,work,work)
        x=0.0D0*a
        y=0.5D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,work,Qx,work,Vx,work)
        x=0.5D0*a
        y=0.0D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,work,work,Qy,work,Vy)
        x=1.0D0*a
        y=1.0D0*b
        call CALC(num,mmax,nmax,D,po,array,a,b,x,y,work,work,work,Mxy,work,work,work,work)
        R=2.0D0*Mxy
        write(6,'(f8.3,f8.5,2(f8.4),5(f8.3))') b/a,w*D/q/a**4.0D0,Mxx/q/a/a,Myy/q/a/a,Qx/q/a,Qy/q/a,Vx/q/a,Vy/q/a,R/q/a/a
    end do

stop

contains

    subroutine LOAD(num,mmax,nmax,D,q,a,b,array)
        real(8),intent(in)::D,q,a,b
        integer,intent(in)::num,mmax,nmax
        real(8),intent(out)::array(1:num,1:num*10)
        integer::m,n
        real(8)::qmn,coef

        do m=1,mmax
        do n=1,nmax
            array(m,n)=0.0D0
        end do
        end do
        do m=1,mmax,2
        do n=1,nmax,2
            qmn=16.0D0*q/pi/pi/dble(m)/dble(n)
!            qmn=8.0D0*q/pi/pi/dble(m)/dble(n)*dble(-1)**(dble(m-1))
            coef=D*pi**4.0D0*(dble(m)*dble(m)/a/a+dble(n)*dble(n)/b/b)**2.0D0
            array(m,n)=qmn/coef
        end do
        end do
    end subroutine LOAD


    subroutine CALC(num,mmax,nmax,D,po,array,a,b,x,y,w,Mxx,Myy,Mxy,Qx,Qy,Vx,Vy)
        real(8),intent(in)::D,po,a,b,x,y
        integer,intent(in)::num,mmax,nmax
        real(8),intent(in)::array(1:num,1:num*10)
        real(8),intent(out)::w,Mxx,Myy,Mxy,Qx,Qy,Vx,Vy
        integer::m,n
        real(8)::dMxydx,dMxydy
        real(8)::snx,sny,csx,csy,dmm,dnn,dmn

        w=0.0D0
        Mxx=0.0D0
        Myy=0.0D0
        Mxy=0.0D0
        Qx=0.0D0
        Qy=0.0D0
        dMxydx=0.0D0
        dMxydy=0.0D0
        do m=1,mmax,2
        do n=1,nmax,2
            snx=sin(dble(m)*pi*x/a)
            sny=sin(dble(n)*pi*y/b)
            csx=cos(dble(m)*pi*x/a)
            csy=cos(dble(n)*pi*y/b)
            dmm=dble(m)*dble(m)
            dnn=dble(n)*dble(n)
            dmn=dble(m)*dble(n)
            w=w+array(m,n)*snx*sny
            Mxx=Mxx+array(m,n)*(dmm/a/a+po*dnn/b/b)*snx*sny
            Myy=Myy+array(m,n)*(po*dmm/a/a+dnn/b/b)*snx*sny
            Mxy=Mxy+array(m,n)*dmn/a/b*csx*csy
            Qx=Qx+array(m,n)*(dmm/a/a+dnn/b/b)*dble(m)/a*csx*sny
            Qy=Qy+array(m,n)*(dmm/a/a+dnn/b/b)*dble(n)/b*snx*csy
            dMxydx=dMxydx-array(m,n)*dmn/a/b*dble(m)*pi/a*snx*csy
            dMxydy=dMxydy-array(m,n)*dmn/a/b*dble(n)*pi/b*csx*sny
        end do
        end do
        Mxx=D*pi*pi*Mxx
        Myy=D*pi*pi*Myy
        Mxy=(1.0D0-po)*D*pi*pi*Mxy
        Qx=D*pi*pi*pi*Qx
        Qy=D*pi*pi*pi*Qy
        Vx=Qx-(1.0D0-po)*D*pi*pi*dMxydy
        Vy=Qy-(1.0D0-po)*D*pi*pi*dMxydx
    end subroutine CALC


end program f90_navier
