module defmat
    implicit none
    ! Characteristics of Concrete
    real(8),parameter::epscu=0.0035D0
    real(8),parameter::gmc=1.5D0
    real(8)::fcu
    ! Characteristics of Reinforcement bar
    real(8),parameter::Es=2.0D5
    real(8),parameter::fy=460.0D0
    real(8),parameter::gms=1.05D0
end module defmat

program f90_RCLS
    ! Ultimate strength of RC member
    use defmat
    implicit none
    integer,parameter::nn=100
    integer::mb,mh,ms
    real(8)::FN,MX,MY,bb,hh
    real(8),allocatable::xc(:),yc(:),rc(:)
    real(8),allocatable::xs(:),ys(:),area(:),rs(:)
    real(8)::ecx,ecy,theta
    real(8)::da,xx,yy,dx,dy
    real(8)::rcmax,rsmax,eps,e0,de,eu,ey
    real(8)::sum_a,sum_m
    real(8)::wa,wm,wwa,wwm
    character(len=50)::fnameR
    integer::i,j,k,iflag

    ! data input
    call getarg(1,fnameR)
    open(11,file=fnameR,status='old')
    read(11,*) fcu,FN,MX,MY,bb,hh,mb,mh,ms
    allocate(xc(1:mb*mh),yc(1:mb*mh),rc(1:mb*mh))
    allocate(xs(1:ms),ys(1:ms),area(1:ms),rs(1:ms))
    do k=1,ms
        read(11,*) xs(k),ys(k),area(k)
    end do
    close(11)

    ! Eccentricity and angle of rotation
    ecy=MX/FN
    ecx=MY/FN
    theta=acos(ecy/sqrt(ecx*ecx+ecy*ecy))

    ! Coordinates of concrete elements and re-bars
    dx=bb/dble(mb)
    dy=hh/dble(mh)
    da=dx*dy
    k=0
    do i=1,mh
        yy=0.5D0*hh+0.5D0*dy-dy*dble(i)
    do j=1,mb
        xx=-0.5D0*bb-0.5D0*dx+dx*dble(j)
        k=k+1
        xc(k)=xx
        yc(k)=yy
    end do
    end do
    ! Arm for moment calculation
    rcmax=0.0D0
    do k=1,mh*mb
        rc(k)=xc(k)*sin(theta)+yc(k)*cos(theta)
        if(rcmax<rc(k))rcmax=rc(k)
    end do
    rsmax=0.0D0
    do k=1,ms
        rs(k)=xs(k)*sin(theta)+ys(k)*cos(theta)
        if(rs(k)<rsmax)rsmax=rs(k)
    end do

    ! Calcculation of axial force and moment
    eu=epscu      ! Ultimate strain of concrete
    ey=-fy/gms/Es ! Yield strain of reinforcement
    de=(eu-ey)/dble(nn)
    e0=eu+de
    iflag=0
    do i=0,nn
        e0=e0-de
        sum_a=0.0D0
        sum_m=0.0D0
        do k=1,mh*mb
            if(iflag==0)then
                eps=e0+(eu-e0)/rcmax*rc(k)
            else
                eps=e0+(ey-e0)/rsmax*rc(k)
            end if
            sum_a=sum_a+SIGC(eps)*da
            sum_m=sum_m+SIGC(eps)*da*rc(k)
        end do
        do k=1,ms
            if(iflag==0)then
                eps=e0+(eu-e0)/rcmax*rs(k)
            else
                eps=e0+(ey-e0)/rsmax*rs(k)
            end if
            sum_a=sum_a+(SIGS(eps)-SIGC(eps))*area(k)
            sum_m=sum_m+(SIGS(eps)-SIGC(eps))*area(k)*rs(k)
        end do
        if(iflag==0.and.e0+(eu-e0)/rcmax*rsmax<=ey)iflag=1
        if(0.0D0<=wa.and.sum_a<=0.0D0)then
            wwa=0.0D0
            wwm=wm-(wm-sum_m)/(wa-sum_a)*wa
            write(6,'(a10,2(e15.7))') 'M0',wwa/bb/hh,wwm/bb/hh/hh
        else
            write(6,'(i10,2(e15.7))') i,sum_a/bb/hh,sum_m/bb/hh/hh
        end if
        wa=sum_a
        wm=sum_m
    end do

stop

contains

    real(8) function SIGC(eps_c)
        ! Stress-strain curve of concrete
        real(8),intent(in)::eps_c
        real(8)::sig,sigy,epsy

        sigy=0.67D0*fcu/gmc
        epsy=2.4D-4*sqrt(fcu/gmc)
        sig=-sigy/epsy/epsy*(eps_c-epsy)*(eps_c-epsy)+sigy
        if(epsy<eps_c)sig=sigy
        if(epscu<eps_c)sig=0.0D0
        if(eps_c<=0.0D0)sig=0.0D0
        SIGC=sig
    end function SIGC

    real(8) function SIGS(eps_s)
        ! Stress-strain cyrve of reinforcement
        real(8),intent(in)::eps_s
        real(8)::sig

        sig=Es*eps_s
        if(fy/gms<sig)sig=fy/gms
        if(sig<-fy/gms)sig=-fy/gms
        SIGS=sig
    end function SIGS

end program f90_RCLS
