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

program f90_natm
    !****************************************************************
    !Ground Response Curve for NATM
    !****************************************************************
    use defpi
    implicit none
    integer,parameter::nnn=11   !Number of time steps
    integer,parameter::mm=100   !Number of devided support pressure
    integer::i,j,k
    real(8)::a        !Radius of the tunnel
    real(8)::P0       !Initial ground pressure
    real(8)::E        !Elastic modulus of thr rock
    real(8)::po       !Poisson's ratio of thr rock
    real(8)::c        !Pure Shearing strength of the rock
    real(8)::phi      !Internal friction angle of the rock
    real(8)::zeta,qu  !Coefficient of passive earth pressure
    real(8)::Pin      !Internal pressure by support
    real(8)::R        !Radius of plastic zone
    real(8)::PR       !Radial stress at radius=R
    real(8)::ua       !Total inner wall displacement
    real(8)::ub       !Total displacement at end point of the rockbolt
    real(8)::Pp       !Ground pressure (working variable)
    real(8)::day3D    !Required time for 3D progress
    real(8)::t_sc     !Thickness of shotcrete
    real(8)::s_t,s_l  !installation pitch of rockbolt
    real(8)::L_rb     !length of rockbolt
    real(8)::A_sr     !Section area of Steel rib
    real(8)::s_sr     !Pitch of steel rib
    real(8)::P1(1:mm+1),u1(1:mm+1),R1(1:mm+1)
    real(8)::P2(1:mm+1),u2(1:mm+1),R2(1:mm+1)
    real(8)::day(1:nnn),rate(1:nnn),utu(1:nnn),utp(1:nnn)
    real(8)::psupport(1:nnn),usupport(1:nnn)
    real(8)::psu1,usu1,psu2,usu2,xx,yy
    real(8)::work,E_sc,weps,deps,wsig
    real(8)::tt,dt    !Time variables
    real(8)::Pcr,wqu
    character(len=20)::dummy,fnameW1,fnameW2

    !Data input
    call getarg( 1,dummy); read(dummy,*) a
    call getarg( 2,dummy); read(dummy,*) P0
    call getarg( 3,dummy); read(dummy,*) E
    call getarg( 4,dummy); read(dummy,*) po
    call getarg( 5,dummy); read(dummy,*) c
    call getarg( 6,dummy); read(dummy,*) phi
    call getarg( 7,dummy); read(dummy,*) day3D
    call getarg( 8,dummy); read(dummy,*) t_sc
    call getarg( 9,dummy); read(dummy,*) s_t
    call getarg(10,dummy); read(dummy,*) s_l
    call getarg(11,dummy); read(dummy,*) L_rb
    call getarg(12,dummy); read(dummy,*) A_sr
    call getarg(13,dummy); read(dummy,*) s_sr
    call getarg(14,fnameW1)
    call getarg(15,fnameW2)

    zeta=(1.0D0+sin(phi/180.0D0*pi))/(1.0D0-sin(phi/180.0D0*pi))
    qu=2.0D0*c*cos(phi/180.0D0*pi)/(1.0D0-sin(phi/180.0D0*pi))

    !Definition of time steps and ground pressure release rate
    do i=1,nnn
        select case(i)
        case(1); xx=0.00D0
        case(2); xx=0.25D0
        case(3); xx=0.50D0
        case(4); xx=0.75D0
        case(5); xx=1.00D0
        case(6); xx=1.25D0
        case(7); xx=1.50D0
        case(8); xx=2.00D0
        case(9); xx=2.50D0
        case(10);xx=3.00D0
        case(11);xx=3.50D0
        end select
        day(i)=xx/3.5D0*day3D
    end do
    rate(1)=0.500D0   !xx=0.00D0
    rate(2)=0.738D0   !xx=0.25D0
    rate(3)=0.860D0   !xx=0.50D0
    rate(4)=0.915D0   !xx=0.75D0
    rate(5)=0.945D0   !xx=1.00D0
    rate(6)=0.962D0   !xx=1.25D0
    rate(7)=0.973D0   !xx=1.50D0
    rate(8)=0.986D0   !xx=2.00D0
    rate(9)=0.993D0   !xx=2.50D0
    rate(10)=0.997D0   !xx=3.00D0
    rate(11)=1.000D0   !xx=3.50D0

    !Limit of the Support pressure and displacement
    tt=0.0D0
    dt=0.01D0
    weps=0.0D0
    usu1=0.0D0
    xx=PRB(0.09D0*L_rb,t_sc,s_t,s_l,L_rb)
    yy=PSR(0.09D0,A_sr,a,s_sr)
    psu1=xx+yy
    k=int(40.0D0/dt)
    open(12,file=fnameW2,status='replace')
    write(12,'(3e15.7)') tt,usu1,psu1
    do i=1,k
        work=SCS(tt+0.5D0*dt)
        E_sc=SCE(work)
        deps=(SCS(tt+dt)-SCS(tt))/E_sc
        weps=weps+deps
        wsig=SCS(tt+dt)
        psu1=PSC(a,wsig,t_sc)+xx+yy
        usu1=(weps)*a
        write(12,'(6e15.7)') tt+dt,usu1,psu1,PSC(a,wsig,t_sc),xx,yy
        do j=1,nnn
            if(day(j)-dt<=tt.and.tt<=day(j)+dt)then
                utu(j)=usu1
                utp(j)=psu1
            end if
        end do
        tt=tt+dt
    end do
    close(12)

    !Ground response curve and displacement estimation
    Pp=rate(1)*P0
!    Pcr=Pp-0.5D0*qu
!    wqu=2.0D0*Pp-Pcr*(1.0D0+zeta)
!    if(wqu<=0.0D0)wqu=0.0D0
    wqu=qu
    do i=1,mm+1
        Pin=Pp-Pp/dble(mm)*dble(i-1)
        if(wqu==0.0D0.and.Pin==0.0D0)Pin=1.0D-6
        call CALRU(Pp,Pin,wqu,zeta,a,E,po,R,ua)
        P1(i)=Pin
        R1(i)=R
        u1(i)=ua
    end do
    usupport(1)=u1(mm+1)
    psupport(1)=0.0D0
    usu1=usupport(1)
    psu1=psupport(1)
    do k=2,nnn
        Pp=rate(k)*P0
!        Pcr=Pp-0.5D0*qu
!        wqu=2.0D0*Pp-Pcr*(1.0D0+zeta)
!        if(wqu<=0.0D0)wqu=0.0D0
        do i=1,mm+1
            Pin=Pp-Pp/dble(mm)*dble(i-1)
            if(wqu==0.0D0.and.Pin==0.0D0)Pin=1.0D-6
            call CALRU(Pp,Pin,wqu,zeta,a,E,po,R,ua)
            P2(i)=Pin
            R2(i)=R
            u2(i)=ua
        end do
        if(t_sc+s_t*s_l+s_sr==0.0D0)then
            usupport(k)=u2(mm+1)  !Total displacement
            psupport(k)=0.0D0     !Support pressure
        else
            usu2=usu1+utu(k)
            psu2=psu1+utp(k)
            do i=1,mm
                call POINT(u2(i),P2(i),u2(i+1),P2(i+1),usupport(k-1),psupport(k-1),usu2,psu2,xx,yy)
                if((u2(i)<=xx.and.xx<=u2(i+1)).and.(P2(i+1)<=yy.and.yy<=P2(i)))exit
            end do
            usupport(k)=xx
            psupport(k)=yy
        end if
    end do

    open(12,file=fnameW1,status='replace')
    write(12,'(i5,6f10.3)') mm,u1(mm+1),usupport(nnn),psupport(nnn),utu(nnn),utp(nnn),day3D
    !Initial displacement and Final displacement of the ground
    do i=1,mm+1
        write(12,'(6e15.7)') u1(i),P1(i),R1(i)-a,u2(i),P2(i),R2(i)-a
    end do
    !Supoort displacement and support pressure
    do k=1,nnn
        write(12,'(6e15.7)') usupport(k),psupport(k),day(k),usupport(k)-usupport(1)
    end do
    close(12)

stop

contains

    subroutine CALRU(Pp,Pin,wqu,zeta,a,E,po,R,ua)
        real(8),intent(in)::Pp,Pin,wqu,zeta,a,E,po
        real(8),intent(out)::R,ua

        R=a
        ua=a*(1.0D0+po)/E*(Pp-Pin)
        if(Pin<(2.0D0*Pp-wqu)/(zeta+1.0D0))then
            R=a*(2.0D0/(zeta+1.0D0)*(Pp+wqu/(zeta-1.0D0))/(Pin+wqu/(zeta-1.0D0)))**(1.0D0/(zeta-1.0D0))
            ua=R*R/a*(1.0D0+po)/E*(zeta-1.0D0)/(zeta+1.0D0)*(Pp+wqu/(zeta-1.0D0))
        end if
    end subroutine CALRU


    real(8) function SCS(dd)
        !Strength of Shotcrete
        real(8),intent(in)::dd
        SCS=7.56D0*sqrt(dd)
    end function SCS


    real(8) function SCE(sig)
        !Elastic modulus of Shotcrete
        real(8),intent(in)::sig
        SCE=1826.0D0*sqrt(sig)
    end function SCE


    real(8) function PSC(a,sig_sc,t_sc)
        !Support pressure by Shotcrete
        real(8),intent(in)::a,sig_sc,t_sc
        PSC=0.0D0
        if(0.0D0<t_sc)PSC=sig_sc*t_sc/a
    end function PSC


    real(8) function PRB(dlenrb,t_sc,s_t,s_l,L_rb)
        !Support pressure by Rockbolt
        real(8),intent(in)::dlenrb,t_sc,s_t,s_l,L_rb
        real(8)::sig_rb,eps_rb
        real(8)::Es=200000.0D0
        real(8)::sigu=226.0D0
        real(8)::As=490.874D0 !D25
        eps_rb=dlenrb/L_rb
        sig_rb=Es*eps_rb
        if(sigu<sig_rb)sig_rb=sigu
        PRB=0.0D0
        if(0.0D0<s_t*s_l)PRB=sig_rb*As/s_t/s_l
        if(0.1D0<eps_rb)PRB=0.0D0
    end function PRB


    real(8) function PSR(eps,a,A_sr,s_sr)
        !Support pressure by Steel Rib
        real(8),intent(in)::eps,a,A_sr,s_sr
        real(8)::sig_sr
        real(8)::Es=200000.0D0
        real(8)::sigu=275.0D0
        sig_sr=Es*eps
        if(sigu<sig_sr)sig_sr=sigu
        PSR=0.0D0
        if(0.0D0<s_sr)PSR=sig_sr*A_sr/s_sr/a
        if(0.1D0<eps)PSR=0.0D0
    end function PSR


    subroutine POINT(x1,y1,x2,y2,x3,y3,x4,y4,xx,yy)
        real(8),intent(in)::x1,y1,x2,y2,x3,y3,x4,y4
        real(8),intent(out)::xx,yy
        real(8)::a2,b2,a1,b1
        a1=(y2-y1)/(x2-x1)
        b1=(x2*y1-x1*y2)/(x2-x1)
        a2=(y4-y3)/(x4-x3)
        b2=(x4*y3-x3*y4)/(x4-x3)
        xx=(b2-b1)/(a1-a2)
        yy=a1*xx+b1
    end subroutine POINT

end program f90_natm
