module defary
    implicit none
    integer,parameter::NDmax=100
    character(len=200)::strcom ! Comment
    integer::nd;             ! Numbers of calculation
    integer::IE_inp(1:NDmax) ! 0:Inner pressure, 1:External pressure
    real(8)::PP_inp(1:NDmax) ! Water pressure
    real(8)::TT_inp(1:NDmax) ! Temperature change (+:increse of temperature)
    real(8)::aa_inp(1:NDmax) ! Internal radius of pressure tunnel
    real(8)::bb_inp(1:NDmax) ! Excavation radius
    real(8)::rr_inp(1:NDmax) ! Model edge radius
    real(8)::cc_inp(1:NDmax) ! Cover for reinforcement
    real(8)::ta_inp(1:NDmax) ! Equivalent thickness of inner reinforcement
    real(8)::tb_inp(1:NDmax) ! Equivalent thickness of outer reinforcement
    real(8)::Ec_inp(1:NDmax) ! Elastic modulus of concrete
    real(8)::nc_inp(1:NDmax) ! Poisson ratio of concrete
    real(8)::ac_inp(1:NDmax) ! Thermal expansion coefficient of concrete
    real(8)::Es_inp(1:NDmax) ! Elastic modulus of reinforcement
    real(8)::ns_inp(1:NDmax) ! Poisson ratio of reinforcement
    real(8)::as_inp(1:NDmax) ! Thermal expansion coefficient of reinforcement
    real(8)::Eg_inp(1:NDmax) ! Elastic modulus of bed rock
    real(8)::ng_inp(1:NDmax) ! Poisson ratio of bed rock
end module defary

!===========================
program f90_rcpt
!===========================
! Stress analysis for Circular Pressure Tunnel
!     using multiple cylindrical theorem
    use defary
    implicit none
    character(len=50)::fnameR,fnameW

    call getarg(1,fnameR)
    call getarg(2,fnameW)
    call main_part(fnameR,fnameW)

stop

contains
!------------------------------------------------------------------------------
subroutine main_part(fnameR,fnameW)
!------------------------------------------------------------------------------
    character(len=50),intent(in)::fnameR,fnameW
    character(len=200)::dat
    character(len=20)::dummy
    integer::i,j,k,kk,io
    integer::ncomma(1:16)
    integer::iwork
    real(8)::rwork

    ! Data input
    open(11,file=fnameR,status='old')
        read(11,*) strcom
        kk=0
        do
            read(11,'(a)',iostat=io) dat
            if(io<0)exit

            if(dat(1:1)/='#')then
                k=0
                i=1
                do
                    j=index(dat(i:),',')
                    if(j==0)exit
                    k=k+1
                    ncomma(k)=i+j-1
                    i=i+j
                    if(16<=k)exit
                end do
                kk=kk+1
                dummy=dat(1:ncomma(1)-1)              ;read(dummy,*) iwork; IE_inp(kk)=iwork
                dummy=dat(ncomma( 1)+1:ncomma( 2)-1)  ;read(dummy,*) rwork; PP_inp(kk)=rwork
                dummy=dat(ncomma( 2)+1:ncomma( 3)-1)  ;read(dummy,*) rwork; TT_inp(kk)=rwork
                dummy=dat(ncomma( 3)+1:ncomma( 4)-1)  ;read(dummy,*) rwork; aa_inp(kk)=rwork
                dummy=dat(ncomma( 4)+1:ncomma( 5)-1)  ;read(dummy,*) rwork; bb_inp(kk)=rwork
                dummy=dat(ncomma( 5)+1:ncomma( 6)-1)  ;read(dummy,*) rwork; rr_inp(kk)=rwork
                dummy=dat(ncomma( 6)+1:ncomma( 7)-1)  ;read(dummy,*) rwork; cc_inp(kk)=rwork
                dummy=dat(ncomma( 7)+1:ncomma( 8)-1)  ;read(dummy,*) rwork; ta_inp(kk)=rwork
                dummy=dat(ncomma( 8)+1:ncomma( 9)-1)  ;read(dummy,*) rwork; tb_inp(kk)=rwork
                dummy=dat(ncomma( 9)+1:ncomma(10)-1)  ;read(dummy,*) rwork; Ec_inp(kk)=rwork
                dummy=dat(ncomma(10)+1:ncomma(11)-1)  ;read(dummy,*) rwork; nc_inp(kk)=rwork
                dummy=dat(ncomma(11)+1:ncomma(12)-1)  ;read(dummy,*) rwork; ac_inp(kk)=rwork
                dummy=dat(ncomma(12)+1:ncomma(13)-1)  ;read(dummy,*) rwork; Es_inp(kk)=rwork
                dummy=dat(ncomma(13)+1:ncomma(14)-1)  ;read(dummy,*) rwork; ns_inp(kk)=rwork
                dummy=dat(ncomma(14)+1:ncomma(15)-1)  ;read(dummy,*) rwork; as_inp(kk)=rwork
                dummy=dat(ncomma(15)+1:ncomma(16)-1)  ;read(dummy,*) rwork; Eg_inp(kk)=rwork
                dummy=dat(ncomma(16)+1:len(trim(dat)));read(dummy,*) rwork; ng_inp(kk)=rwork
            end if
        end do
    close(11)
    nd=kk

    ! Making output file
    open(12, file=fnameW, status='replace')
        write(12,'(a)') trim(adjustl(strcom))
        write(12,'(a)') '*Input data'
        write(12,'(a)') 'k,IE,PP,TT,aa,bb,rr,cc,ta,tb,Ec,nc,ac,Es,ns,as,Eg,ng'
        do i=1,nd
            write(12,'(i5,",",i5,16(",",f10.3))') k,IE_inp(i),&
            PP_inp(i),TT_inp(i),aa_inp(i),bb_inp(i),rr_inp(i),cc_inp(i),ta_inp(i),tb_inp(i),&
            Ec_inp(i),nc_inp(i),ac_inp(i),Es_inp(i),ns_inp(i),as_inp(i),Eg_inp(i),ng_inp(i)
        end do
        write(12,'(a)') '*Output data'
        write(12,'(a)') 'k,IE,Eg,dT,sr_c,st_c,sr_si1,st_si1,sr_si2,st_si2,sr_so1,st_so1,sr_so2,st_so2,sr_g,st_g,ua,ub'
    close(12)

    ! k     : Sequence
    ! IE    : Inner pressure (=0) or Outer pressure (=1)
    ! PP    : Water pressure
    ! dT    : Temperature change
    ! sr_c  : Stress in the r-direction at inner surface concrete
    ! st_c  : Stress in the t-direction at inner surface concrete
    ! sr_si : Stress in the r-direction of inner reinforcement
    ! st_si : Stress in the t-direction of inner reinforcement
    ! sr_so : Stress in the r-direction of outer reinforcement
    ! st_so : Stress in the t-direction of outer reinforcement
    ! sr_g  : Stress in the r-direction at outer surface concrete or bed rock
    ! st_g  : Stress in the t-direction at outer surface concrete or bed rock
    ! ua    : Radial displacement of inner surface concrete
    ! ub    : Radial displacement of outer surface concrete

    ! Execution of calculation
    do k=1,nd
        select case(IE_inp(k))
            case(0)
                if(tb_inp(k)<=0)call main_sng_pin(fnameW,k)
                if(0<tb_inp(k)) call main_dbl_pin(fnameW,k)
            case(1)
                if(tb_inp(k)<=0)call main_sng_pex(fnameW,k)
                if(0<tb_inp(k)) call main_dbl_pex(fnameW,k)
        end select
    end do

end subroutine main_part
!------------------------------------------------------------------------------
subroutine main_sng_pin(fnameW,k)
!------------------------------------------------------------------------------
    character(len=50),intent(in)::fnameW
    integer,intent(in)::k
    integer,parameter::n=8
    integer::i,j
    real(8)::a(1:n,1:n+1)
    real(8)::ra,r1,r2,rb,rr
    real(8)::Eg,ng,Es,ns,alps,Ec,alpc
    real(8)::P0,DT
    real(8)::C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt
    real(8)::ug,sigrg,sigtg
    real(8)::us2,sigrs2,sigts2
    real(8)::us1,sigrs1,sigts1
    real(8)::uc,sigrc,sigtc=0.0D0

    ra  =aa_inp(k)
    r1  =aa_inp(k)+cc_inp(k)
    r2  =r1+ta_inp(k)
    rb  =bb_inp(k)
    rr  =rr_inp(k)
    Eg  =Eg_inp(k)
    ng  =ng_inp(k)
    Es  =Es_inp(k)
    ns  =ns_inp(k)
    alps=as_inp(k)
    Ec  =Ec_inp(k)
    alpc=ac_inp(k)
    P0  =PP_inp(k) ! Positive : internal pressure
    DT  =TT_inp(k)

    do i=1,n
        do j=1,n+1
            a(i,j)=0.0D0
        end do
    end do

    a(1,1)=rr
    a(1,2)= 1.0D0/rr
    a(2,1)= Eg/(1.0D0+ng)/(1.0D0-2.0D0*ng)
    a(2,2)=-Eg/(1.0D0+ng)/rb/rb
    a(2,3)= 0.0D0
    a(2,4)=-Ec/rb
    a(3,1)= rb
    a(3,2)= 1.0D0/rb
    a(3,3)=-1.0D0
    a(3,4)=-log(rb)
    a(4,3)= 0.0D0
    a(4,4)= Ec/r2
    a(4,5)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a(4,6)= Es/(1.0D0+ns)/r2/r2
    a(5,3)= 1.0D0
    a(5,4)= log(r2)
    a(5,5)=-r2
    a(5,6)=-1.0D0/r2
    a(6,5)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a(6,6)=-Es/(1.0D0+ns)/r1/r1
    a(6,7)= 0.0D0
    a(6,8)=-Ec/r1
    a(7,5)= r1
    a(7,6)= 1.0D0/r1
    a(7,7)=-1.0D0
    a(7,8)=-log(r1)
    a(8,7)= 0.0D0
    a(8,8)= Ec/ra;
    a(1,9)= 0.0D0
    a(2,9)= 0.0D0
    a(3,9)= alpc*DT*(rb-r2)
    a(4,9)=-Es*alps*DT/(1.0D0-ns)*(r2*r2-r1*r1)/2.0D0/r2/r2
    a(5,9)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r2*r2-r1*r1)/2.0D0/r2
    a(6,9)= 0.0D0
    a(7,9)= alpc*DT*(r1-ra)
    a(8,9)=-P0

    call STEGJ(n,a)

    C1=a(1,9);C2=a(2,9);xx=rb;xa=rb;E=Eg;nu=ng;alp=0.0D0;TT=0.0D0
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,ug,sigrg,sigtg)

    C1=a(5,9);C2=a(6,9);xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us2,sigrs2,sigts2)

    C1=a(5,9);C2=a(6,9);xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us1,sigrs1,sigts1)

    C1=a(7,9);C2=a(8,9);xx=ra;xa=ra;E=Ec;alp=alpc;TT=DT
    call CON(C1,C2,xx,xa,E,alp,DT,uc,sigrc)

    open(12, file=fnameW, position='append')
        write(12,'(i5,",",i5,2(",",f10.0),14(",",f15.3))') &
        k,IE_inp(k),Eg,DT,sigrc,sigtc,sigrs1,sigts1,sigrs2,sigts2,0.0D0,0.0D0,0.0D0,0.0D0,sigrg,sigtg,uc,ug
    close(12)

    write(6,'(a)') '===Single: External pressure==='
    write(6,'(3(e15.7))') uc,sigrc,sigtc
    write(6,'(3(e15.7))') us1,sigrs1,sigts1
    write(6,'(3(e15.7))') us2,sigrs2,sigts2
    write(6,'(3(e15.7))') ug,sigrg,sigtg

end subroutine main_sng_pin
!------------------------------------------------------------------------------
subroutine main_dbl_pin(fnameW,k)
!------------------------------------------------------------------------------
    character(len=50),intent(in)::fnameW
    integer,intent(in)::k
    integer,parameter::n=12
    integer i,j
    real(8)::a(1:n,1:n+1)
    real(8)::ra,r1,r2,r3,r4,rb,rr
    real(8)::Eg,ng,Es,ns,alps,Ec,alpc
    real(8)::P0,DT
    real(8)::C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt
    real(8)::ug,sigrg,sigtg
    real(8)::us4,sigrs4,sigts4
    real(8)::us3,sigrs3,sigts3
    real(8)::us2,sigrs2,sigts2
    real(8)::us1,sigrs1,sigts1
    real(8)::uc,sigrc,sigtc=0.0D0

    ra  =aa_inp(k)
    r1  =ra+cc_inp(k)
    r2  =r1+ta_inp(k)
    rb  =bb_inp(k)
    r4  =rb-cc_inp(k)
    r3  =r4-tb_inp(k)
    rr  =rr_inp(k)
    Eg  =Eg_inp(k)
    ng  =ng_inp(k)
    Es  =Es_inp(k)
    ns  =ns_inp(k)
    alps=as_inp(k)
    Ec  =Ec_inp(k)
    alpc=ac_inp(k)
    P0  =PP_inp(k) ! Positive : internal pressure
    DT  =TT_inp(k)

    do i=1,n
        do j=1,n+1
            a(i,j)=0.0D0
        end do
    end do
    a( 1, 1)= rr
    a( 1, 2)= 1.0D0/rr
    a( 2, 1)= Eg/(1.0D0+ng)/(1.0D0-2.0D0*ng)
    a( 2, 2)=-Eg/(1.0D0+ng)/rb/rb
    a( 2, 3)= 0.0D0
    a( 2, 4)=-Ec/rb
    a( 3, 1)= rb
    a( 3, 2)= 1.0D0/rb
    a( 3, 3)=-1.0D0
    a( 3, 4)=-log(rb)
    a( 4, 3)= 0.0D0
    a( 4, 4)= Ec/r4
    a( 4, 5)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 4, 6)= Es/(1.0D0+ns)/r4/r4
    a( 5, 3)= 1.0D0
    a( 5, 4)= log(r4)
    a( 5, 5)=-r4
    a( 5, 6)=-1.0D0/r4
    a( 6, 5)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 6, 6)=-Es/(1.0D0+ns)/r3/r3
    a( 6, 7)= 0.0D0
    a( 6, 8)=-Ec/r3
    a( 7, 5)= r3
    a( 7, 6)= 1.0D0/r3
    a( 7, 7)=-1.0D0
    a( 7, 8)=-log(r3)
    a( 8, 7)= 0.0D0
    a( 8, 8)= Ec/r2
    a( 8, 9)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 8,10)= Es/(1.0D0+ns)/r2/r2
    a( 9, 7)= 1.0D0
    a( 9, 8)= log(r2)
    a( 9, 9)=-r2
    a( 9,10)=-1.0D0/r2
    a(10, 9)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a(10,10)=-Es/(1.0D0+ns)/r1/r1
    a(10,11)= 0.0D0
    a(10,12)=-Ec/r1
    a(11, 9)= r1
    a(11,10)= 1.0D0/r1
    a(11,11)=-1.0D0
    a(11,12)=-log(r1)
    a(12,11)= 0.0D0
    a(12,12)= Ec/ra
    a( 1,13)= 0.0D0
    a( 2,13)= 0.0D0
    a( 3,13)= alpc*DT*(rb-r2)
    a( 4,13)=-Es*alps*DT/(1.0D0-ns)*(r4*r4-r3*r3)/2.0D0/r4/r4
    a( 5,13)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r4*r4-r3*r3)/2.0D0/r4
    a( 6,13)= 0.0D0
    a( 7,13)= alpc*DT*(r3-r2)
    a( 8,13)=-Es*alps*DT/(1.0D0-ns)*(r2*r2-r1*r1)/2.0D0/r2/r2
    a( 9,13)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r2*r2-r1*r1)/2.0D0/r2
    a(10,13)= 0.0D0
    a(11,13)= alpc*DT*(r1-ra)
    a(12,13)=-P0

    call STEGJ(n,a)

    C1=a(1,13);C2=a(2,13);xx=rb;xa=rb;E=Eg;nu=ng;alp=0.0D0;TT=0.0D0
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,ug,sigrg,sigtg)

    C1=a(5,13);C2=a(6,13);xx=r4;xa=r3;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us4,sigrs4,sigts4)

    C1=a(5,13);C2=a(6,13);xx=r3;xa=r3;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us3,sigrs3,sigts3)

    C1=a(9,13);C2=a(10,13);xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us2,sigrs2,sigts2)

    C1=a(9,13);C2=a(10,13);xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us1,sigrs1,sigts1)

    C1=a(11,13);C2=a(12,13);xx=ra;xa=ra;E=Ec;alp=alpc;TT=DT
    call CON(C1,C2,xx,xa,E,alp,DT,uc,sigrc)

    open(12, file=fnameW, position='append')
        write(12,'(i5,",",i5,2(",",f10.0),14(",",f15.3))') &
        k,IE_inp(k),Eg,DT,sigrc,sigtc,sigrs1,sigts1,sigrs2,sigts2,sigrs3,sigts3,sigrs4,sigts4,&
        sigrg,sigtg,uc,ug
    close(12)

    write(6,'(a)') '===Double: External pressure==='
    write(6,'(3(e15.7))') uc,sigrc,sigtc
    write(6,'(3(e15.7))') us1,sigrs1,sigts1
    write(6,'(3(e15.7))') us2,sigrs2,sigts2
    write(6,'(3(e15.7))') us3,sigrs3,sigts3
    write(6,'(3(e15.7))') us4,sigrs4,sigts4
    write(6,'(3(e15.7))') ug,sigrg,sigtg

end subroutine main_dbl_pin
!------------------------------------------------------------------------------
subroutine main_sng_pex(fnameW,k)
!------------------------------------------------------------------------------
    character(len=50),intent(in)::fnameW
    integer,intent(in)::k
    integer,parameter::n=6
    integer::i,j
    real(8)::a(1:n,1:n+1)
    real(8)::ra,r1,r2,rb,r
    real(8)::Eg=0.0,Es,ns,alps,Ec,nc,alpc
    real(8)::Pe,DT
    real(8)::C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt
    real(8)::ug,sigrg,sigtg
    real(8)::us2,sigrs2,sigts2
    real(8)::us1,sigrs1,sigts1
    real(8)::uc,sigrc,sigtc

    ra  =aa_inp(k)
    r1  =aa_inp(k)+cc_inp(k)
    r2  =r1+ta_inp(k)
    rb  =bb_inp(k)
    Es  =Es_inp(k)
    ns  =ns_inp(k)
    alps=as_inp(k)
    Ec  =Ec_inp(k)
    nc  =nc_inp(k)
    alpc=ac_inp(k)
    Pe  =PP_inp(k) ! Positive : outer pressure
    DT  =TT_inp(k)

    do i=1,n
        do j=1,n+1
            a(i,j)=0.0D0
        end do
    end do
    a(1,1)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a(1,2)=-Ec/(1.0D0+nc)/rb/rb
    a(2,1)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a(2,2)=-Ec/(1.0D0+nc)/r2/r2
    a(2,3)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a(2,4)= Es/(1.0D0+ns)/r2/r2;
    a(3,1)= r2
    a(3,2)= 1.0D0/r2
    a(3,3)=-r2
    a(3,4)=-1.0D0/r2;
    a(4,3)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a(4,4)=-Es/(1.0D0+ns)/r1/r1
    a(4,5)=-Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a(4,6)= Ec/(1.0D0+nc)/r1/r1
    a(5,3)= r1
    a(5,4)= 1.0D0/r1
    a(5,5)=-r1
    a(5,6)=-1.0D0/r1
    a(6,5)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a(6,6)=-Ec/(1.0D0+nc)/ra/ra;
    a(1,7)=-Pe
    a(2,7)=-Es*alps*DT/(1.0D0-ns)*(r2*r2-r1*r1)/2.0D0/r2/r2
    a(3,7)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r2*r2-r1*r1)/2.0D0/r2
    a(4,7)=-Ec*alpc*DT/(1.0D0-nc)*(r1*r1-ra*ra)/2.0D0/r1/r1
    a(5,7)= (1.0D0+nc)/(1.0D0-nc)*alpc*DT*(r1*r1-ra*ra)/2.0D0/r1
    a(6,7)= 0.0D0

    call STEGJ(n,a)

    C1=a(1,7);C2=a(2,7);xx=rb;xa=r2;E=Ec;nu=nc;alp=alpc;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,ug,sigrg,sigtg)

    C1=a(3,7);C2=a(4,7);xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us2,sigrs2,sigts2)

    C1=a(3,7);C2=a(4,7);xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us1,sigrs1,sigts1)

    C1=a(5,7);C2=a(6,7);xx=ra;xa=ra;E=Ec;nu=nc;alp=alpc;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,uc,sigrc,sigtc)

    open(12, file=fnameW, position='append')
        write(12,'(i5,",",i5,2(",",f10.0),14(",",f15.3))') &
        k,IE_inp(k),Eg,DT,sigrc,sigtc,sigrs1,sigts1,sigrs2,sigts2,0.0D0,0.0D0,0.0D0,0.0D0,sigrg,sigtg,uc,ug
    close(12)

    write(6,'(a)') '===Single: External pressure==='
    write(6,'(3(e15.7))') uc,sigrc,sigtc
    write(6,'(3(e15.7))') us1,sigrs1,sigts1
    write(6,'(3(e15.7))') us2,sigrs2,sigts2
    write(6,'(3(e15.7))') ug,sigrg,sigtg

end subroutine main_sng_pex
!------------------------------------------------------------------------------
subroutine main_dbl_pex(fnameW,k)
!------------------------------------------------------------------------------
    character(len=50),intent(in)::fnameW
    integer,intent(in)::k
    integer,parameter::n=10
    integer::i,j
    real(8)::a(1:n,1:n+1)
    real(8)::ra,r1,r2,r3,r4,rb,rr
    real(8)::Eg=0.0,Es,ns,alps,Ec,nc,alpc
    real(8)::Pe,DT
    real(8)::C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt
    real(8)::ug,sigrg,sigtg
    real(8)::us4,sigrs4,sigts4
    real(8)::us3,sigrs3,sigts3
    real(8)::us2,sigrs2,sigts2
    real(8)::us1,sigrs1,sigts1
    real(8)::uc,sigrc,sigtc

    ra  =aa_inp(k)
    r1  =ra+cc_inp(k)
    r2  =r1+ta_inp(k)
    rb  =bb_inp(k)
    r4  =rb-cc_inp(k)
    r3  =r4-tb_inp(k)
    Es  =Es_inp(k)
    ns  =ns_inp(k)
    alps=as_inp(k)
    Ec  =Ec_inp(k)
    nc  =nc_inp(k)
    alpc=ac_inp(k)
    Pe  =PP_inp(k) ! Positive : outer pressure
    DT  =TT_inp(k)

    do i=1,n
        do j=1,n+1
            a(i,j)=0.0D0
        end do
    end do
    a( 1, 1)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a( 1, 2)=-Ec/(1.0D0+nc)/rb/rb
    a( 2, 1)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a( 2, 2)=-Ec/(1.0D0+nc)/r4/r4
    a( 2, 3)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 2, 4)= Es/(1.0D0+ns)/r4/r4
    a( 3, 1)= r4
    a( 3, 2)= 1.0D0/r4
    a( 3, 3)=-r4
    a( 3, 4)=-1.0D0/r4;
    a( 4, 3)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 4, 4)=-Es/(1.0D0+ns)/r3/r3
    a( 4, 5)=-Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a( 4, 6)= Ec/(1.0D0+nc)/r3/r3;
    a( 5, 3)= r3
    a( 5, 4)= 1.0D0/r3
    a( 5, 5)=-r3
    a( 5, 6)=-1.0D0/r3;
    a( 6, 5)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a( 6, 6)=-Ec/(1.0D0+nc)/r2/r2
    a( 6, 7)=-Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 6, 8)= Es/(1.0D0+ns)/r2/r2
    a( 7, 5)= r2
    a( 7, 6)= 1.0D0/r2
    a( 7, 7)=-r2
    a( 7, 8)=-1.0D0/r2;
    a( 8, 7)= Es/(1.0D0+ns)/(1.0D0-2.0D0*ns)
    a( 8, 8)=-Es/(1.0D0+ns)/r1/r1
    a( 8, 9)=-Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a( 8,10)= Ec/(1.0D0+nc)/r1/r1
    a( 9, 7)= r1
    a( 9, 8)= 1.0D0/r1
    a( 9, 9)=-r1
    a( 9,10)=-1.0D0/r1;
    a(10, 9)= Ec/(1.0D0+nc)/(1.0D0-2.0D0*nc)
    a(10,10)=-Ec/(1.0D0+nc)/ra/ra
    a( 1,11)=-Pe
    a( 2,11)=-Es*alps*DT/(1.0D0-ns)*(r4*r4-r3*r3)/2.0D0/r4/r4
    a( 3,11)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r4*r4-r3*r3)/2.0D0/r4
    a( 4,11)=-Ec*alpc*DT/(1.0D0-nc)*(r3*r3-r2*r2)/2.0D0/r3/r3
    a( 5,11)= (1.0D0+nc)/(1.0D0-nc)*alpc*DT*(r3*r3-r2*r2)/2.0D0/r3
    a( 6,11)=-Es*alps*DT/(1.0D0-ns)*(r2*r2-r1*r1)/2.0D0/r2/r2
    a( 7,11)= (1.0D0+ns)/(1.0D0-ns)*alps*DT*(r2*r2-r1*r1)/2.0D0/r2
    a( 8,11)=-Ec*alpc*DT/(1.0D0-nc)*(r1*r1-ra*ra)/2.0D0/r1/r1
    a( 9,11)= (1.0D0+nc)/(1.0D0-nc)*alpc*DT*(r1*r1-ra*ra)/2.0D0/r1
    a(10,11)= 0.0D0

    call STEGJ(n,a)

    C1=a(1,11);C2=a(2,11);xx=rb;xa=r4;E=Ec;nu=nc;alp=alpc;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,ug,sigrg,sigtg)

    C1=a(3,11);C2=a(4,11);xx=r4;xa=r3;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us4,sigrs4,sigts4)

    C1=a(3,11);C2=a(4,11);xx=r3;xa=r3;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us3,sigrs3,sigts3)

    C1=a(7,11);C2=a(8,11);xx=r2;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us2,sigrs2,sigts2)

    C1=a(7,11);C2=a(8,11);xx=r1;xa=r1;E=Es;nu=ns;alp=alps;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,us1,sigrs1,sigts1)

    C1=a(9,11);C2=a(10,11);xx=ra;xa=ra;E=Ec;nu=nc;alp=alpc;TT=DT
    call ELAS(C1,C2,xx,xa,E,nu,alp,TT,uc,sigrc,sigtc)

    open(12, file=fnameW, position='append')
        write(12,'(i5,",",i5,2(",",f10.0),14(",",f15.3))') &
        k,IE_inp(k),Eg,DT,sigrc,sigtc,sigrs1,sigts1,sigrs2,sigts2,sigrs3,sigts3,sigrs4,sigts4,&
        sigrg,sigtg,uc,ug
    close(12)

    write(6,'(a)') '===Double: External pressure==='
    write(6,'(3(e15.7))') uc,sigrc,sigtc
    write(6,'(3(e15.7))') us1,sigrs1,sigts1
    write(6,'(3(e15.7))') us2,sigrs2,sigts2
    write(6,'(3(e15.7))') us3,sigrs3,sigts3
    write(6,'(3(e15.7))') us4,sigrs4,sigts4
    write(6,'(3(e15.7))') ug,sigrg,sigtg

end subroutine main_dbl_pex
!------------------------------------------------------------------------------
subroutine ELAS(C1,C2,xx,xa,E,nu,alp,TT,u,sigr,sigt)
!------------------------------------------------------------------------------
    real(8),intent(in)::C1,C2,xx,xa,E,nu,alp,TT
    real(8),intent(out)::u,sigr,sigt
    u=(1.0D0+nu)/(1.0D0-nu)*alp*TT*(xx*xx-xa*xa)/2.0D0/xx+C1*xx+C2/xx
    sigr=-E*alp*TT/(1.0D0-nu)*(xx*xx-xa*xa)/2.0D0/xx/xx+E/(1.0D0+nu)/(1.0D0-2.0D0*nu)*C1-E/(1.0D0+nu)*C2/xx/xx
    sigt=-E*alp*TT/(1.0D0-nu)*(xx*xx+xa*xa)/2.0D0/xx/xx+E/(1.0D0+nu)/(1.0D0-2.0D0*nu)*C1+E/(1.0D0+nu)*C2/xx/xx
end subroutine ELAS
!------------------------------------------------------------------------------
subroutine CON(C1,C2,xx,xa,E,alp,TT,u,sigr)
!------------------------------------------------------------------------------
    real(8),intent(in)::C1,C2,xx,xa,E,alp,TT
    real(8),intent(out)::u,sigr
    u=C1+C2*log(xx)+alp*TT*(xx-xa)
    sigr=E*C2/xx
end subroutine CON
!------------------------------------------------------------------------------
    subroutine STEGJ(n,ary)
!------------------------------------------------------------------------------
        !Gauss-Jordan法による連立一次方程式の解法
        !河西朝雄著「C 言語によるはじめてのアルゴリズム入門」
        integer,intent(in)::n
        real(8),intent(inout)::ary(1:n,1:n+1)

        integer::i,j,k,s
        real(8)::p,d,max,dumy

        do k=1,n
            max=0.0D0
            s=k
            do j=k,n
                if(abs(ary(j,k))>max)then
                    max=abs(ary(j,k))
                    s=j
                end if
            end do
            do j=1,n+1
                dumy=ary(k,j)
                ary(k,j)=ary(s,j)
                ary(s,j)=dumy
            end do
            p=ary(k,k)
            do j=k,n+1
                ary(k,j)=ary(k,j)/p
            end do
            do i=1,n
                if(i/=k)then
                    d=ary(i,k)
                    do j=k,n+1
                        ary(i,j)=ary(i,j)-d*ary(k,j)
                    end do
                end if
            end do
        end do
    end subroutine STEGJ
!------------------------------------------------------------------------------
end program f90_rcpt
!------------------------------------------------------------------------------
