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


program f90_tunnel
    use defpi
    implicit none
    character(len=50)::fnameW
    character(len=50)::fnameF
    character(len=50)::dummy
    integer::ksp
    real(8)::DIA,Tc,Tx,W,ts,tp,theta
    real(8)::B,H,Ae,Ac,Ls,DAc,DLs

    call getarg(1,fnameW)                  !Script for GMT execution
    call getarg(2,fnameF)                  !Output image filename
    call getarg(3,dummy);read(dummy,*) ksp !0: no-steel lining, 1: steel lining
    call getarg(4,dummy);read(dummy,*) DIA !Internal diameter of waterway
    call getarg(5,dummy);read(dummy,*) Tc  !Thickness of concrete lining
    call getarg(6,dummy);read(dummy,*) Tx  !Additional thickness of invert concrete
    call getarg(7,dummy);read(dummy,*) W   !Width of invert concrete
    call getarg(8,dummy);read(dummy,*) ts  !Thickness of shotcrete
    call getarg(9,dummy);read(dummy,*) tp  !Thickness from design line to pay line

    theta=CAL_THETA(DIA,Tc,Tx,W)
    ! Excavated area (Pay line)
        B=DIA+2.0D0*Tc+2.0D0*ts+2.0D0*tp
        H=DIA+2.0D0*Tc+Tx+ts+2.0D0*tp
        Ae=CAL_AREA(B,H,theta)
        Ae=Ae*1.0D-6
    ! Area of concrete lining (Concrete)
        B=DIA+2.0D0*Tc
        H=DIA+2.0D0*Tc+Tx
        Ac=CAL_AREA(B,H,theta)-0.25D0*pi*DIA*DIA
        B=DIA+2.0D0*Tc+2.0D0*tp
        H=DIA+2.0D0*Tc+Tx+tp
        DAc=(CAL_AREA(B,H,theta)-0.25D0*pi*DIA*DIA)-Ac
        Ac=Ac*1.0D-6
        DAc=DAc*1.0D-6
    ! Length of shot crete (Shotcrete)
        B=DIA+2.0D0*Tc+2.0D0*ts
        H=DIA+2.0D0*Tc+Tx+ts
        Ls=CAL_LENG(B,H,theta)
        B=DIA+2.0D0*Tc+2.0D0*ts+2.0D0*tp
        H=DIA+2.0D0*Tc+Tx+ts+tp
        DLs=CAL_LENG(B,H,theta)-Ls
        Ls=Ls*1.0D-3
        DLs=DLs*1.0D-3
    write(6,'("Ae=",f6.2," Ac=",f6.2," DAc=",f6.2," Ls=",f6.2," DLs=",f6.2)') Ae,Ac,DAc,Ls,DLs

    call MAKEBAT(fnameW,fnameF,ksp,DIA,Tc,Tx,W,ts,tp,theta,Ae,Ac,Ls,DAc,DLs)

    call system('bash '//fnameW)

stop

contains

    real(8) function CAL_AREA(B,H,theta)
        real(8),intent(in)::B,H,theta
        real(8)::rr,hh,ang,A0,A1,A2

        ang=theta/180.0D0*pi
        rr=0.5D0*B
        hh=H-rr-rr*sin(ang)
        A0=(2.0D0*rr*cos(ang)-hh*tan(ang))*hh ! Area of trapezoid with height of hh
        A1=rr*rr*(ang+0.5D0*pi)               ! Area of sector
        A2=rr*rr*sin(ang)*cos(ang)            ! Area of triangle
        CAL_AREA=A0+A1+A2
    end function CAL_AREA


    real(8) function CAL_LENG(B,H,theta)
        real(8),intent(in)::B,H,theta
        real(8)::rr,hh,ang,L0,L1

        ang=theta/180.0D0*pi
        rr=0.5D0*B
        hh=H-rr-rr*sin(ang)
        L0=rr*(pi+2.0D0*ang) ! Length of sector
        L1=hh/cos(ang)     ! Length of straight line
        CAL_LENG=L0+2.0D0*L1
    end function CAL_LENG


    real(8) function CAL_THETA(DIA,Tc,Tx,W)
        real(8),intent(in)::DIA,Tc,Tx,W
        ! Set angle (degree)
        real(8)::r,f,x1,x2,x3,f1,f2,f3

        r=0.5D0*DIA+Tc
        x1=0.0D0
        x2=90.0D0
        x3=0.5D0*(x1+x2)
        do while(abs(x1-x2)>1.0D-6)
            f1=CAL_F(x1,W,r,Tx)
            f2=CAL_F(x2,W,r,Tx)
            f3=CAL_F(x3,W,r,Tx)
            if(f1*f3<0.0D0)then
                x2=x3
                x3=0.5D0*(x1+x2)
            end if
            if(f3*f2<0.0D0)then
                x1=x3
                x3=0.5D0*(x1+x2)
            end if
        end do
        CAL_THETA=0.5D0*(x1+x2)
    end function CAL_THETA


    real(8) function CAL_F(ang,W,r,Tx)
        real(8),intent(in)::ang,W,r,Tx

        CAL_F=W-(2.0D0*r*cos(ang/180.0D0*pi)-2.0D0*(r+Tx-r*sin(ang/180.0D0*pi))*tan(ang/180.0D0*pi))
    end function CAL_F


    subroutine MAKEBAT(fnameW,fnameF,ksp,DIA,Tc,Tx,W,ts,tp,theta,Ae,Ac,Ls,DAc,DLs)
        character(len=50),intent(in)::fnameW,fnameF
        integer,intent(in)::ksp
        real(8),intent(in)::DIA,Tc,Tx,W,ts,tp,theta
        real(8),intent(in)::Ae,Ac,Ls,DAc,DLs
        real(8)::x1,y1,x2,y2,x3,y3,x4,y4,h,r,dt
        real(8)::ds,dss,dl,dsh,phi
        integer::fsize
        real(8)::sang
        character(len=50)::dummy
        character(len=100)::linebuf
        integer::i,nd

        fsize=8
        ds=1000.0D0
        dss=200.0D0
        dl=0.5D0*ds
        dsh=dss*2.0D0
        nd=360
        h=0.5D0*DIA+Tc+Tx-DIA*sin(theta*pi/180.0D0)

        !Outer line of shotcrete
        open(12, file='_temp1.txt',status='replace')
        r=0.5D0*DIA+Tc+ts
        x1=0.5D0*W+ts/cos(theta*pi/180.0D0)
        y1=-(0.5D0*DIA+Tc+Tx)
        x2= r*cos(theta/180.0D0*pi)
        y2=-r*sin(theta/180.0D0*pi)
        x3=-x2
        y3= y2
        x4=-x1
        y4= y1
        write(12,'(2f10.3)') x1,y1
        write(12,'(2f10.3)') x2,y2
        dt=(180.0D0+2.0D0*theta)/180.0D0*pi/dble(nd)
        do i=1,nd-1
            phi=-theta/180.0D0*pi+dt*dble(i)
            write(12,'(2f10.3)') r*cos(phi),r*sin(phi)
        end do
        write(12,'(2f10.3)') x3,y3
        write(12,'(2f10.3)') x4,y4
        write(12,'(2f10.3)') x1,y1
        close(12)

        !Outer line of concrete lining
        open(12, file='_temp2.txt',status='replace')
        r=0.5D0*DIA+Tc
        x1=0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx)
        x2= r*cos(theta/180.0D0*pi)
        y2=-r*sin(theta/180.0D0*pi)
        x3=-x2
        y3= y2
        x4=-x1
        y4= y1
        write(12,'(2f10.3)') x1,y1
        write(12,'(2f10.3)') x2,y2
        dt=(180.0D0+2.0D0*theta)/180.0D0*pi/dble(nd)
        do i=1,nd-1
            phi=-theta/180.0D0*pi+dt*dble(i)
            write(12,'(2f10.3)') r*cos(phi),r*sin(phi)
        end do
        write(12,'(2f10.3)') x3,y3
        write(12,'(2f10.3)') x4,y4
        write(12,'(2f10.3)') x1,y1
        close(12)

        !Internal surface line
        if(ksp==1)then
            open(12, file='_tempS.txt',status='replace')
            r=0.5D0*(DIA+150.0D0)
            dt=360.0D0/180.0D0*pi/dble(nd)
            do i=0,nd
                phi=0.0D0+dt*dble(i)
                write(12,'(2f10.3)') r*cos(phi),r*sin(phi)
            end do
            close(12)
        end if

        !Internal surface line
        open(12, file='_temp3.txt',status='replace')
        r=0.5D0*(DIA)
        dt=360.0D0/180.0D0*pi/dble(nd)
        do i=0,nd
            phi=0.0D0+dt*dble(i)
            write(12,'(2f10.3)') r*cos(phi),r*sin(phi)
        end do
        close(12)

        !Center line, spring line
        open(12, file='_temp4.txt',status='replace')
        write(12,'(">")')
        x1=0.0D0
        y1=-(0.5D0*DIA+Tc+Tx+0.5D0*ds)
        write(12,'(2f10.3)') x1,y1
        x1=0.0D0
        y1=0.5D0*DIA+Tc+0.5D0*ds
        write(12,'(2f10.3)') x1,y1
        write(12,'(">")')
        x1=-(0.5D0*DIA+Tc+ts+0.5D0*ds)
        y1=0.0D0
        write(12,'(2f10.3)') x1,y1
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds+dss
        y1=0.0D0
        write(12,'(2f10.3)') x1,y1
        close(12)

        !Additional line
        open(12, file='_temp5.txt',status='replace')
        !Horizontal line
        x1=0.0D0
        y1=0.5D0*DIA+Tc
        x2=0.5D0*DIA+Tc+ts+2.0D0*ds+dss
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.0D0
        y1=0.5D0*DIA
        x2=0.5D0*DIA+Tc+ts+1.0D0*ds+dss
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.0D0
        y1=0.0D0
        x2= (0.5D0*DIA+Tc)*cos(theta/180.0D0*pi)
        y2=-(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1= (0.5D0*DIA+Tc)*cos(theta/180.0D0*pi)
        y1=-(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        x2=0.5D0*DIA+Tc+ts+1.0D0*ds+dss
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.0D0
        y1=-0.5D0*DIA
        x2=0.5D0*DIA+Tc+ts+1.0D0*ds+dss
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx)
        x2=0.5D0*DIA+Tc+ts+2.0D0*ds+dss
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        !Vertical
        x1=-(0.5D0*DIA+Tc)
        y1=0.0D0
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx+2.0D0*ds+dss)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=-0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx)
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx+1.0D0*ds+dss)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx)
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx+1.0D0*ds+dss)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*DIA+Tc
        y1=0.0D0
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx+2.0D0*ds+dss)
        write(12,'(4f10.3)') x1,y1,x2,y2
        close(12)

        !Arrow line (both)
        open(12, file='_temp6.txt',status='replace')
        !Vertical
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds
        y1=0.5D0*DIA
        x2=x1
        y2=0.0D0
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds
        y1=0.0D0
        x2=x1
        y2=-(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds
        y1=-(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        x2=x1
        y2=-0.5D0*DIA
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*DIA+Tc+ts+2.0D0*ds
        y1=0.5D0*DIA+Tc
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx)
        write(12,'(4f10.3)') x1,y1,x2,y2
        !Horizontal
        x1=-(0.5D0*DIA+Tc)
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds)
        x2=-0.5D0*W
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=-0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds)
        x2=0.5D0*W
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*W
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds)
        x2=0.5D0*DIA+Tc
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=-(0.5D0*DIA+Tc)
        y1=-(0.5D0*DIA+Tc+Tx+2.0D0*ds)
        x2=0.5D0*DIA+Tc
        y2=y1
        write(12,'(4f10.3)') x1,y1,x2,y2
        close(12)

        !Arrow line (single)
        open(12, file='_temp7.txt',status='replace')
        x1=0.0D0
        y1=0.0D0
        x2=0.5D0*DIA*cos(0.25*pi)
        y2=0.5D0*DIA*sin(0.25*pi)
        write(12,'(4f10.3)') x1,y1,x2,y2
        !Vertical
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds
        y1=0.5D0*DIA+Tc+dl
        x2=x1
        y2=0.5D0*DIA+Tc
        write(12,'(4f10.3)') x1,y1,x2,y2
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds
        y1=-(0.5D0*DIA+Tc+Tx+dl)
        x2=x1
        y2=-(0.5D0*DIA+Tc+Tx)
        write(12,'(4f10.3)') x1,y1,x2,y2
        x2=-(0.5D0*DIA+Tc+ts)*cos(0.25D0*pi)
        y2= (0.5D0*DIA+Tc+ts)*sin(0.25D0*pi)
        x1=x2-0.5D0*ds
        y1=y2+0.5D0*ds
        write(12,'(4f10.3)') x1,y1,x2,y2
        if(ksp==1)then
            x2=-(0.5D0*DIA+75.0D0)*sin(0.1D0*pi)
            y2= (0.5D0*DIA+75.0D0)*cos(0.1D0*pi)
            x1=x2-(Tc+ts+1.0D0*ds)*sin(0.1D0*pi)
            y1=y2+(Tc+ts+1.0D0*ds)*cos(0.1D0*pi)
            write(12,'(4f10.3)') x1,y1,x2,y2
        end if
        close(12)

        !Character
        open(12, file='_temp8.txt',status='replace')
        !Vertical
        x2=Tc
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds+dsh
        y1=0.5D0*DIA+0.5D0*x2
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=0.5D0*DIA
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds-dsh
        y1=0.5D0*x2
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds+dsh
        y1=-0.5D0*x2
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=0.5D0*DIA-(0.5D0*DIA+Tc)*sin(theta/180.0D0*pi)
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds-dsh
        y1=-0.5D0*DIA+0.5D0*x2
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=Tc+Tx
        x1=0.5D0*DIA+Tc+ts+1.0D0*ds+dsh
        y1=-(0.5D0*DIA+0.5*x2)
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=Tc+DIA+Tc+Tx
        x1=0.5D0*DIA+Tc+ts+2.0D0*ds+dsh
        y1=0.5D0*DIA+Tc-0.5D0*x2
        sang=90.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        !Horizontal
        x2=0.5D0*DIA+Tc-0.5D0*W
        x1=-0.5D0*W-0.5D0*x2
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds+dsh)
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=W
        x1=0.0D0
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds+dsh)
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=0.5D0*DIA+Tc-0.5D0*W
        x1=0.5D0*W+0.5D0*x2
        y1=-(0.5D0*DIA+Tc+Tx+1.0D0*ds+dsh)
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=DIA+2.0D0*Tc
        x1=0.0D0
        y1=-(0.5D0*DIA+Tc+Tx+2.0D0*ds+dsh)
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        !Thickness of shotcrete
        x1=-(0.5D0*DIA+Tc+ts)*cos(0.25D0*pi)-0.5D0*ds-dsh
        y1= (0.5D0*DIA+Tc+ts)*sin(0.25D0*pi)+0.5D0*ds+dsh
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'("t="i3)') nint(ts)
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x1=x1
        y1=y1+1.2D0*dsh
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            dummy='Shotcrete'
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        if(ksp==1)then
            x2=-(0.5D0*DIA+75.0D0)*sin(0.1D0*pi)
            y2= (0.5D0*DIA+75.0D0)*cos(0.1D0*pi)
            x1=x2-(Tc+ts+1.0D0*ds)*sin(0.1D0*pi)
            y1=y2+(Tc+ts+1.0D0*ds)*cos(0.1D0*pi)+dsh
            sang=0.0D0
                write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
                dummy='Steel pipe'
                linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))
                write(12,'(a)') trim(adjustl(linebuf))
        end if
        !angle
        x2=0.5D0*DIA
        x1=0.5D0*DIA/2.0D0*cos(0.25D0*pi)-dsh*sin(0.25D0*pi)
        y1=0.5D0*DIA/2.0D0*sin(0.25D0*pi)+dsh*cos(0.25D0*pi)
        sang=45.0D0
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(i6)') nint(x2)
            linebuf=trim(adjustl(linebuf))//' r='//trim(adjustl(dummy))
            write(12,'(a)') trim(adjustl(linebuf))
        x2=theta
        x1= 0.5D0*DIA/2.0D0*cos(theta/180.0D0*pi)-dsh*sin(theta/180.0D0*pi)
        y1=-0.5D0*DIA/2.0D0*sin(theta/180.0D0*pi)-dsh*cos(theta/180.0D0*pi)
        sang=-x2
            write(linebuf,'(2f10.3,i3,f10.3," MC")') x1,y1,fsize,sang
            write(dummy,'(f6.3)') x2
            linebuf=trim(adjustl(linebuf))//' '//trim(adjustl(dummy))//'@%12%\260@%%'
            write(12,'(a)') trim(adjustl(linebuf))
        !Quantity
        x1=0.0D0
        y1=0.5D0*DIA+Tc+ts+tp+2.0D0*ds
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," ML")') x1,y1,fsize,sang
            write(dummy,'(f6.2)') Ls
            linebuf=trim(adjustl(linebuf))//' Ls='//trim(adjustl(dummy))
            write(dummy,'(f6.2)') DLs
            linebuf=trim(adjustl(linebuf))//'m@+3@+ (@~D@~Ls='//trim(adjustl(dummy))//'m@+3@+)'
            write(12,'(a)') trim(adjustl(linebuf))
        y1=y1+2.0D0*dsh
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," ML")') x1,y1,fsize,sang
            write(dummy,'(f6.2)') Ac
            linebuf=trim(adjustl(linebuf))//' Ac='//trim(adjustl(dummy))
            write(dummy,'(f6.2)') DAc
            linebuf=trim(adjustl(linebuf))//'m@+3@+ (@~D@~Ac='//trim(adjustl(dummy))//'m@+3@+)'
            write(12,'(a)') trim(adjustl(linebuf))
        y1=y1+2.0D0*dsh
        sang=0.0D0
            write(linebuf,'(2f10.3,i3,f10.3," 0 ML")') x1,y1,fsize,sang
            write(dummy,'(f6.2)') Ae
            linebuf=trim(adjustl(linebuf))//' Ae='//trim(adjustl(dummy))//'m@+3@+'
            write(12,'(a)') trim(adjustl(linebuf))
        close(12)

        !Make script file
        open(12, file=fnameW,status='replace')
        linebuf='fig='//fnameF
        write(12,'(a)') trim(adjustl(linebuf))
        write(12,'(a)') 'range=-15000/15000/-15000/15000'
        write(12,'(a)') 'scale=15/15'
        write(12,'(a)') 'gmt set MAP_VECTOR_SHAPE 1'
        write(12,'(a)') 'gmt psxy _temp1.txt -R$range -JX$scale -W1 -P -K > $fig'
!        write(12,'(a)') 'gmt psxy _temp2.txt -R -J -Gp400/62 -W1 -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp2.txt -R -J -G#eeeeee -W1 -O -K >> $fig'
        if(ksp==1)write(12,'(a)') 'gmt psxy _tempS.txt -R -J -W1 -Gblack -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp3.txt -R -J -W1 -Gwhite -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp4.txt -R -J -W0.5,10_3_3_3:0 -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp5.txt -R -J -Sv0.2+s -W0.5 -G0 -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp6.txt -R -J -Sv0.2+s+b+e+a45 -W0.5 -G0 -O -K >> $fig'
        write(12,'(a)') 'gmt psxy _temp7.txt -R -J -Sv0.2+s+e+a45 -W0.5 -G0 -O -K >> $fig'
        write(12,'(a)') 'gmt pstext _temp8.txt -R -J -F+f+a+j -N -O -K >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        close(12)

    end subroutine MAKEBAT

end program f90_tunnel
