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

program f90_toroidal
    use defpi
    implicit none
    integer,parameter::nd=360
    integer::NODT,NELT,MATEL,KOZ,KOR,NF
    real(8)::dt,aa,bb,z0,pp,tt
    real(8)::theta,cs,sn,zz,rr,dl,r0,pz1,pr1,pz2,pr2
    real(8)::z(1:nd*2),r(1:nd*2),fz(1:nd),fr(1:nd)
    integer::i,j,ne

    NODT=2*nd
    NELT=nd
    MATEL=1
    KOZ=0
    KOR=0
    NF=nd
    dt=360.0D0/dble(nd)
    aa=2000.0D0
    bb=4000.0D0
    z0=5000.0D0
    pp=1.0D0
    tt=10.0D0

    write(6,*) 'Toroidal'
    write(6,*) NODT,NELT,MATEL,KOZ,KOR,NF,1
    write(6,*) 2.0D5,0.3,2.5,0,1.0D-5,100000
    
    do ne=1,NELT-1
        write(6,*) ne,ne+1,ne+1+nd,ne+nd,1
    end do
    ne=NELT
    write(6,*) ne,1,1+nd,NODT,1

    do i=1,nd
        theta=dt*dble(i-1)
        z(i)=z0+aa*sin(theta/180.0D0*pi)
        r(i)=bb+aa*cos(theta/180.0D0*pi)
        z(i+nd)=z0+(aa+tt)*sin(theta/180.0D0*pi)
        r(i+nd)=bb+(aa+tt)*cos(theta/180.0D0*pi)
    end do
    do i=1,NODT
        write(6,*) z(i),r(i),0
    end do

    ne=1
        i=ne
        j=ne+1
        if(ne==nd)j=1
        zz=0.5D0*(z(i)+z(j))
        rr=0.5D0*(r(i)+r(j))
        cs=((rr-bb)/sqrt((zz-z0)**2+(rr-bb)**2))
        sn=((zz-z0)/sqrt((zz-z0)**2+(rr-bb)**2))
        dl=sqrt((z(i)-z(j))**2+(r(i)-r(j))**2)
        r0=0.5D0*(r(i)+r(j))
        pz1=pp*dl*sn*r0
        pr1=pp*dl*cs*r0
        fz(1)=0.5D0*pz1
        fr(1)=0.5D0*pr1
    do ne=2,nd
        i=ne
        j=ne+1
        if(ne==nd)j=1
        zz=0.5D0*(z(i)+z(j))
        rr=0.5D0*(r(i)+r(j))
        cs=((rr-bb)/sqrt((zz-z0)**2+(rr-bb)**2))
        sn=((zz-z0)/sqrt((zz-z0)**2+(rr-bb)**2))
        dl=sqrt((z(i)-z(j))**2+(r(i)-r(j))**2)
        r0=0.5D0*(r(i)+r(j))
        pz2=pp*dl*sn*r0
        pr2=pp*dl*cs*r0
        fz(i)=0.5D0*(pz1+pz2)
        fr(i)=0.5D0*(pr1+pr2)
        pz1=pz2
        pr1=pr2
    end do
    fz(1)=fz(1)+0.5D0*pz2
    fr(1)=fr(1)+0.5D0*pr2
    do i=1,NF
        write(6,*) i,fz(i),fr(i)
    end do

end program f90_toroidal