Running

      program main
      implicit real *8 (a-h,o-z)
      parameter (nnmax=10000,nemax=10000,ncase_max=10)
    parameter (nn = 10000)
      parameter (iabaqus =1, zero=0.0d0) ! set it to 0 for pc software
      character (len=1) :: inter
      character (len=10) :: fname(ncase_max)
      dimension :: inode(nnmax),coord(2,nnmax)
      dimension :: ielem(4,nemax),inelem(4,nemax)
      dimension :: kp(11) !! previous it is kp(11), 08/09
      dimension :: num_nd(nnmax),coord_nd(2,nnmax)
      dimension :: larray_ele(8,nemax)
      dimension :: key_point(11)
    dimension :: nodeflag(nnmax) !!this array is for node's material identification
      dimension :: num_nd_T(nnmax),coord_nd_T(2,nnmax)
      dimension :: larray_ele_T(8,nemax)
      dimension :: key_point_T(11)

    REAL*8 VR(NN),PEAK(NN)
    REAL*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
    REAL*8 XFE1,XPE1
    REAL*8 FC1,PC1,BC1

    REAL*8 C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    REAL*8 XFE2,XPE2,FC2,PC2,BC2

    REAL*8 DELTAT,XTERM,DIFF
    real*8 AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1&    ,Tliquid1,Tsolid1
    real*8 AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2&    ,Tliquid2,Tsolid2
    REAL*8 XA(NN),XF(NN),XP(NN),XB(NN),XM(NN),GS(NN),G(NN),HV(NN)&  ,XL(NN)

    COMMON/COMP/C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
    COMMON/COMP2/C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    COMMON/TEMPS/AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1&  ,Tliquid1,Tsolid1
    COMMON/TEMPS2/AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2& ,Tliquid2,Tsolid2
    COMMON/STRUCT/XA,XF,XP,XB,XM,GS,G,HV,XL
    COMMON/KCONST/fc1,pc1,bc1
    COMMON/KCONST2/fc2,pc2,bc2
    COMMON/BASELINE/XFE1,XPE1
    COMMON/BASELINE2/XFE2,XPE2
    common/nodetype/nodeflag

      include 'common.f'
      fname(1) = 'case1.dat'
      fname(2) = 'case2.dat'
      fname(3) = 'case3.dat'
      fname(4) = 'case4.dat'
      fname(5) = 'case5.dat'
      fname(6) = 'case6.dat'
      fname(7) = 'case7.dat'
      fname(8) = 'case8.dat'
      fname(9) = 'case9.dat'
      fname(10) = 'case10.dat'
    OLDMESHER = 1 ! SET THIS TO 1 IF YOU WANTTO USE THE NEW MESH

      open ( 9,file='data.txt')
      open (12,file='file.log')
      open (97,file='HT_STAT.DAT')
      open (98,file='HT_ERR.DAT')

      read (9,*) ncase
      if (ncase .gt. ncase_max) then
         print *, 'Too many cases, precess stoped'
         stop
      end if

 1    format(1x,A,I2,A,I2,A)
 3    format(A12,A12,5x,A12,A12)
      do 900 icase=1,ncase

         open (10,file=fname(icase))

         print 1, '***** CASE #',icase,&        ' OUT OF ',ncase, ' CASE(S) *****'
         write(12,1) '***** CASE #',icase,&        ' OUT OF ',ncase, ' CASE(S) *****'
         n_node = 0
         n_elem = 0

         do i=1,nnmax
            num_nd(i) = 0
            coord_nd(1,i) = zero
            coord_nd(2,i) = zero
            inode(i) = 0
            coord(1,i) = zero
            coord(2,i) = zero
         enddo

         do i=1,nemax
            do j=1,8
               larray_ele(j,i) = 0
            enddo
            do j=1,4
               ielem(j,i) = 0
               inelem(j,i) = 0
            enddo
         enddo

         do i=1,11
            key_point(i) = 0
         enddo

         print 1, 'Read the parameters/data into the program.'
         call setup

         if (id_conf .eq. 0) then

            call mesh00(thickness_p,inode,coord,nnode, &           ielem,inelem,nelem,key_point)

         elseif (id_conf .eq. 1) then

            d_pipe = d_in + 2*thickness_p
            vl_pipe = 40.0*max(thickness_p,thickness_s)

            if (vl_pipe .gt. 500.0d0) vl_pipe=500.0d0

            IF ( OLDMESHER .LE. 0.5d0 ) THEN

             call mesh01(d_pipe,thickness_p,thickness_s,gap, &  area_w,angle_w,vl_pipe,nnode,nelem,& num_nd,coord_nd,larray_ele,key_point,r_cc)
            ELSE
             call mesh_FINE01(d_pipe,thickness_p,thickness_s,gap,&  area_w,angle_w,vl_pipe,nnode,nelem,& num_nd,coord_nd,larray_ele,key_point,r_cc,iabaqus)
            END IF
    print *, "Mesh completed!"
         elseif (id_conf .eq. 2) then

            d_branch = d_branch_in
            vl_runpipe = 20.0d0*thickness_s
            vl_branch = 20.0d0*thickness_s

            IF ( OLDMESHER .LE. 0.5d0 ) THEN
                call mesh02(d_branch,thickness_p,thickness_s,&  gap,gap2,area_w,angle_wrtb,vl_runpipe,& vl_branch,size_w,nnode,nelem,num_nd,coord_nd,&  larray_ele,key_point,r_cc)
              ELSE

                call mesh02_FINE(d_branch,thickness_p,thickness_s, &  gap,gap2,area_w,angle_wrtb,vl_runpipe,&  vl_branch,size_w,nnode,nelem,num_nd,coord_nd,& larray_ele,key_point,r_cc,iabaqus)
            END IF

         elseif (id_conf .eq. 3) then

            d_pipe = d_in + 2.*thickness_p ! outside radius of the pipe
            vl_pipe = 40.0d0*thickness_p
            call mesh03(d_pipe,thickness_p,area_w,vl_pipe,& size_w,nnode,nelem,num_nd,&    coord_nd,larray_ele,key_point,r_cc)
         else

            print 1, "Invalid configuration, analysis stoped"
            stop

         endif

         if (key_point(1) .eq. 0) then

            print 1, 'Analysis of case #',icase,&           ' has failed!'
            print *, "An error has occurred.  The weld geometry ", &           "can not be setup properly."
            print *, "Please review the error message for details."
            write(12,1) 'Analysis of case #',icase,&        ' has failed!'
            write(10,3) 'NA','NA','NA','NA'
            if (icase .lt. ncase) then
               print *, "Do you wish to continue the analysis ", &              "for the rest cases?"
               print *, "Please enter 'y' or 'Y' to continue!"
               print *, "Or enter any other key to terminate!"
               read(*,'(A)') inter
               if ((inter .eq. 'y') .or. (inter .eq. 'Y')) then
                  goto 800
               else
                  stop
               endif
            else
               goto 800
            endif

         endif

         print *, "There are",nnode," nodes in the FEA model."
         print *, "There are",nelem," elements in the FEA model."

         if (id_conf .ne. 0) then

            print *, "Resort node numbers and coordinates."
            do i=1,nnode
               inode(i) = num_nd(i)
               coord(1,i) = coord_nd(1,i)
               coord(2,i) = coord_nd(2,i)
            enddo

            print *, "Resort element numbers and connectivity."
            do i=1,nelem
               ielem(1,i) = larray_ele(1,i)
               nn1 = larray_ele(2,i)
               nn2 = larray_ele(3,i)
               nn3 = larray_ele(4,i)
               nn4 = larray_ele(5,i)

               ielem(2,i) = larray_ele(6,i)
               ielem(3,i) = larray_ele(7,i)
               ielem(4,i) = larray_ele(8,i)
               do ii=1,nnode
                  if (inode(ii) .eq. nn1) inelem(1,i)=ii
                  if (inode(ii) .eq. nn2) inelem(2,i)=ii
                  if (inode(ii) .eq. nn3) inelem(3,i)=ii
                  if (inode(ii) .eq. nn4) inelem(4,i)=ii
               enddo
            enddo

          do i=1,nelem
            idKey = ielem(2,i)
            do ii = 1,4
                innn = inelem(ii,i)
                nodeflag(innn) = idKey
            end do
          end do
         endif

         do i=1,11
            if (key_point(i) .ne. 0) then
               kpp = key_point(i)
               kpn = i
               do ii=1,nnode
                  if (inode(ii) .eq. kpp) then
                     kp(i) = ii
                     exit
                  endif
               enddo
            endif
         enddo
    CONTINUE

         if ((icase .eq. ncase) .and. (iabaqus .eq. 1)) then
            call modelout(nnode,inode,coord,nelem, &           ielem,inelem,kpn,key_point)
         endif

    if ( oldmesher .gt. 0.5d0 ) then
        call PrepareForMicrostructureCalculation(nnode)
    end if

         call solver (1,nnode,inode,coord, &     nelem,ielem,inelem,kpn,kp,istat)

         if (istat .eq. 0) then
            write(*,1) 'SOLUTIONS FOR CASE #',ICASE, &           ' HAVE FAILED TO CONVERGE!'
            write(*,1) 'NO RESULTS IS AVAILABLE FOR THIS CASE!'
            write(98,1) 'SOLUTION FOR CASE #',ICASE,&           ' HAS FAILED TO CONVERGE!'
            write(98,1) 'NO RESULTS IS AVAILABLE FOR THIS CASE!'
            write(10,3) 'NA','NA','NA','NA'
            if (icase .lt. ncase) then
               print *, "Do you wish to continue the analysis ",&              "for the rest cases?"
               print *, "please enter 'y' or 'Y to continue!"
               read(*,'(A)') inter
               if ((inter .eq. 'y') .or. (inter .eq. 'Y')) then
               else
                  stop
               endif
            endif
         endif

 800     close(10)
         write(*,*)
         write(12,*)
 900  continue

      write(97,*) 'OK'

      stop
      end

      subroutine modelout(nnode,inode,coord,nelem, &     ielem,inelem,kpn,key_point)

      implicit double precision (a-h,o-z)
      character (len=8) :: elset
      dimension :: inode(*),coord(2,*)
      dimension :: ielem(4,*),inelem(4,*)
      dimension :: key_point(kpn)
      dimension :: nelweld(nelem),nelapp(nelem),nelbase(nelem)

      include 'common.f'
      open(15,file='node.inp')
      open(16,file='elem.inp')
      open(17,file='convects')
      open(18,file='parm.f')
      open(19,file='sets.inp')

      write(17,'(A)') '*FILM,OP=NEW'

      do i=1,nnode
         if (coord(1,i) .le. 1.0d-6) coord(1,i)=0.0d0
         write (15,5) i,(coord(ii,i),ii=1,2)
      enddo
 5    format(i5,2(', ',ES15.7))
      close(15)

      iw=0
      iapp=0
      ib=0
      do i=1,nelem
         write(16,15) ielem(1,i),(inelem(j,i),j=1,4)
         if (ielem(2,i) .eq. 3) then
            iw=iw+1
            nelweld(iw)=ielem(1,i)
         elseif (ielem(2,i) .eq. 2) then
            iapp=iapp+1
            nelapp(iapp)=ielem(1,i)
         elseif (ielem(2,i) .eq. 1) then
            ib=ib+1
            nelbase(ib)=ielem(1,i)
         endif

         if (ielem(3,i) .eq. 2) then ! forced convection
            write(17,25) ielem(1,i),ielem(4,i),&           t0_f,2.5E-4
         elseif (ielem(3,i) .eq. 1) then ! free convection
            write(17,25) ielem(1,i),ielem(4,i),&           t_env,h_env
         endif
      enddo
 15   format(i5,9(:',',i5))
 20   format(i5,', F',I1,'NU')
 25   format(I8,', F',I1,',',ES12.5,', ',ES12.5)

      if (iw .gt. 0) then
         elset='EWELD'
         call elsetwrite(elset,iw,nelweld)
      endif
      if(iapp .gt. 0) then
         if (id_conf .eq. 1) elset='SLEEVE'
         if (id_conf .eq. 2) elset='BRANCH'
         call elsetwrite(elset,iapp,nelapp)
      endif
      if(ib .gt. 0) then
         elset='PIPE'
         call elsetwrite(elset,ib,nelbase)
      endif

      write(18,35) 'block data kparm'
      write(18,35) 'implicit real*8 (a-h,o-z)'
      write(18,35) 'parameter (nprecd=2)'
      write(18,35) 'common / block1 / heat, speed, ',&     'rxx, ryy, dzz'
      write(18,35) 'common / block2 / xcent, ycent'
      write(18,35) 'common / block3 / insd'
      write(18,35) 'common / block4 / d_in, id_f, t0_f, frate'

      write(18,45) 'data  heat / ', heat_in
      write(18,45) 'data speed / ', speed
      write(18,45) 'data   dzz / ', r_cc
      write(18,55) 'data  insd / ', key_point(4)
      write(18,45) 'data  d_in / ', d_in
      write(18,55) 'data  id_f / ', id_f
      write(18,45) 'data  t0_f / ', t0_f
      write(18,45) 'data frate / ', velocity
      write(18,35) 'end'
 35   format(6x,A,A)
 45   format(6x,A,ES12.5,' /')
 55   format(6x,A,I5,' /')

      close(15)
      close(16)
      close(17)
      close(18)

      return
      end

      subroutine elsetwrite(name,iel,nelem)
      implicit real*8 (a-h,o-z)
      character (len=8) name
      dimension nelem(iel)

      write(19,'(A,A,A)') '*ELSET, ELSET=',name

      il=iel/10
      ir=mod(iel,10)

      if (il .gt. 0) then
         do i=1,il
            write(19,3) (nelem(j),j=10*(i-1)+1,10*i)
         enddo
      endif
      if (ir .gt. 0) then
         write(19,3) (nelem(j),j=10*il+1,iel)       !9/10/2009
      endif
 3    format(1x,i5,9(:',',1x,i4))

      return
      end

      subroutine setup
      implicit double precision (a-h,o-z)

      include 'common.f'

      pi = acos(-1.0d0)
      emissivity = 0.1
      h_env = 8.4d-5

      read(9,*) id_case, id_conf
      read(9,*) id_m(1),d_in,thickness_p,carbon_eq,t0_p
      read(9,*) id_f,velocity,pressure,t0_f,id_vl
      if (id_conf .eq. 0) then  ! heat sink capacity
         read(9,*) heat_in,igas,fl,oriface,dheight,&        efficiency,angle
         read(9,*)
      else

         read(9,*) id_m(3),id_ele,current,voltage,speed,&        efficiency,angle_w
         read(9,*) id_m(2),d_branch_in,thickness_s,gap,gap2,&        t0_s,angle_wrtb
         if (gap2 .lt. 10.d-3) gap2 = gap
      endif
      read(9,*) t_env

      if (carbon_eq .gt. 1.0d-3) then
         read(9,*) w_b,w_c,w_cr,w_cu,w_mn,w_mo
         read(9,*) w_ni,w_nb,w_n,w_si,w_v
         read(9,*) w_b2,w_c2,w_cr2,w_cu2,w_mn2,w_mo2
         read(9,*) w_ni2,w_nb2,w_n2,w_si2,w_v2
      endif

    read(9,*)   DMESHSIZE,oldmesher, gs01,gs02

    if ( carbon_eq .le. 1.0d-3) then
        oldmesher = 0.0
    end if

      if (id_vl .eq. 1) then
         t0 = 288.7               ! temperature in standard cond.
         p0 = 1.03250d5           ! pressure in standard cond. (1 atm)
         call zcompcalc(id_f,t0,p0,zcomp0)
         pressure = pressure+p0
         call zcompcalc(id_f,t0_f,pressure,zcomp)
         zfact = p0*t0_f/(t0*pressure)*zcomp/zcomp0
         velocity = velocity*zfact
      endif

      ts(1) = 1700.0d0            ! carbon steel
      tl(1) = 1780.0d0            ! carbon steel
      hf(1) = 2.72d2              ! carbon steel
      ts(2) = 1720.0d0            ! austenitic stainless steel
      tl(2) = 1780.0d0            ! austenitic stainless steel
      hf(2) = 2.47d2              ! austenitic stainless steel
      ts(3) = 1720.0d0            ! duplex stainless steel
      tl(3) = 1780.0d0            ! duplex stainless steel
      hf(3) = 2.47d2              ! duplex stainless steel

      if (id_conf .eq. 0) then

         tcr1 = 5.231d2         ! 250 deg C
         tcr2 = 3.731d2         ! 100 deg C
         tcr3 = 4.481d2         ! 175 deg C

         if (heat_in .ne. 0.0) then
            efficiency=0.40
            heat_in=1.05506d3*heat_in*efficiency
         else

         endif

         area_w = pi*25.4**2
      else

         tcr1 = 1.073d3         ! 800 deg C
         tcr2 = 7.731d2         ! 500 deg C
         tcr3 = 8.108d2         ! 1000 deg F
         heat_in = current*voltage
         speed_ips = speed/25.4

         do i=6,12,6
            write(i,5) "Welding Parameters:"
            write(i,5) "Current = ",current, " A"
            write(i,5) "Voltage = ",voltage, " V"
            write(i,5) "Speed = ",speed, " mm/sec"
            write(i,5) "Heat Input =", heat_in/speed_ips/1000.,&           " kJ/in"
            write(i,5) "Arc Efficiency =", efficiency
         enddo
 5       format(1x,A,F6.2,1x,A)

         heat_in = heat_in * efficiency/0.75d0

         if (id_ele .eq. 1) then ! electrode type E7018
            area_w = -3.556d-3 + 8.754d-7*heat_in/speed_ips
         elseif (id_ele .eq. 2) then ! electrode type E7010
            area_w = 1.109d-3 + 8.283d-7*heat_in/speed_ips
         elseif (id_ele .eq. 3) then !  electrode type stainless
            area_w = -1.65d-3 + 8.495d-7*heat_in/speed_ips
         else
            write(12,*) 'Invalid electrode type, analysis stopped!'
            stop
         endif

         heat_in = current*voltage*efficiency
         area_w = area_w*25.4*25.4
         write(12,*) 'weld area = ',area_w
         if (area_w .le. 0.0d0) then
            write(98,*) 'Estimated weld cross section area is ',area_w, &           ' mm'
            write(98,*) 'Information about welding heat-input is '
            write(98,*) 'incorrect.  Go back and change them!'
            stop
         endif
         endif
      IF(ID_CONF.NE.1) GO TO 2111
      XWELDLEG = 2.0D0*AREA_W**0.5D0
      WELDXXX=THICKNESS_S+GAP
      IF(XWELDLEG.LE.WELDXXX) GO TO 2111
      write(12,*) 'WELD LEG LENGTH ',XWELDLEG,&           ' mm'
      write(12,*) 'GREATER THAN SLEEVE PLUS GAP ', WELDXXX, 'mm'
      write(12,*) 'incorrect.  CHANGE INPUTS!'
      write(97,*) 'WELD LEG LENGTH ',XWELDLEG, &           ' mm'
      write(97,*) 'GREATER THAN SLEEVE PLUS GAP ', WELDXXX, 'mm'
      write(97,*) 'incorrect.  CHANGE INPUTS!'
      STOP

 2111 write(*,*) 'Parameter reading and processing completed.'

      return
      end

      subroutine init(nelem,ielem,inelem,&     nnode,t,told,f,fold)
      implicit double precision (a-h,o-z)

      dimension :: ielem(4,*),inelem(4,*)
      dimension :: t(nnode),told(nnode)
      dimension :: f(nnode),fold(nnode)

      include 'common.f'

      do i=1,nelem

         id=ielem(2,i)

         if (id .eq. 1) then
            do ii=1,4
               told(inelem(ii,i))=t0_p
               t(inelem(ii,i))=t0_p
            enddo
         elseif (id .eq. 2) then
            do ii=1,4
               told(inelem(ii,i))=t0_s
               t(inelem(ii,i))=t0_s
            end do
         elseif (id .eq. 3) then
            do ii=1,4
               told(inelem(ii,i))=t_env
               t(inelem(ii,i))=t_env
            end do
         endif
      enddo

      do i=1,nnode
         f(i)=0.0d0
         fold(i)=0.0d0
      enddo

      return
      end

      subroutine jacob(x,vj,r_a,i_gaus,dphi)
      implicit real *8 (a-h,o-z)
      dimension x(2,4), a(2,2),r_a(2,2)
      dimension dphi(4,4,2)

      do ii=1,2
         do jj=1,2
            a(ii,jj)=0.0d0
            r_a(ii,jj)=0.0d0
         enddo
      enddo

      do i=1,4                  ! Gaussian points
         do ii=1,2              ! local coordinates
            do jj=1,2           ! global coordinates
               a(ii,jj)=x(jj,i)*dphi(i,i_gaus,ii)+a(ii,jj)
            enddo
         enddo
      enddo

      vj=a(1,1)*a(2,2)-a(1,2)*a(2,1)
      if (vj .lt. 0.0d0) then
         print *, 'Nagative Jacobian.  Analysis Stopped!'
         stop
      end if

      r_a(1,1)=a(2,2)/vj
      r_a(2,1)=-a(2,1)/vj
      r_a(1,2)=-a(1,2)/vj
      r_a(2,2)=a(1,1)/vj

      return
      end

      subroutine phi(iflag,value,i,j,k)
      implicit real*8 (a-h,o-z)
      parameter (r_gau=0.577350269189626)
      dimension matrix(4,2)

      data ((matrix(ii,jj),ii=1,4),jj=1,2)&     /-1,1,1,-1,-1,-1,1,1/

      if (iflag .eq. 0) then
         value=(1.0d0+matrix(j,1)*matrix(i,1)*r_gau)&        *(1.0d0+matrix(j,2)*matrix(i,2)*r_gau)
      elseif (iflag .eq. 1) then
         if (k .eq. 1) then
            k1=2
         endif
         if (k .eq. 2) then
            k1=1
         endif
         value=matrix(i,k)&        *(1.0d0+matrix(j,k1)*matrix(i,k1)*r_gau)
      else
         value=0.0d0
      endif

      value=value*0.25d0

      return
      end

      subroutine dflux(flux,time,anc,pmax,weldtime)
      implicit real*8 (a-h,o-z)

      include 'common.f'

      if (time .gt. weldtime) then
         flux = 0.0d0
         return
      endif

      dist = time*speed
      if (dist .le. r_cc) then
         dz = 1.0d0 - dist/r_cc
         flux=pmax*exp(-3*dz*dz)
      else
         dz = (dist-r_cc)/(anc*r_cc)
         flux=pmax*exp(-3*dz*dz)
      endif

      return
      end

      subroutine sflux(flux,x)
      implicit real*8 (a-h,o-z)

      include 'common.f'

      flux = -heat_in/area_w

      return
      end


      subroutine fflux(flux,tsurface)
      implicit real*8 (a-h,o-z)

      parameter (simga = 5.67d-14)

      include 'common.f'

      flux = h_env*(tsurface-t_env) + sigma*emissivity*&     (tsurface**4-t_env**4)

      return
      end

      function cp(temp,ID)
      implicit real*8 (a-h,o-z)

      if (id .eq. 1)  then

         cp = cp_steel1(temp)

      elseif (id .eq. 2) then

         cp = cp_steel2(temp)

      elseif (id .eq. 3) then

         cp = cp_steel3(temp)

      endif

      end

      function cp_steel1(temp)
      implicit real *8 (a-h,o-z)
      dimension t(11),value(11)

      data (value(i),t(i),i=1,11) /&  486.0d-3, 100.0d0,&  498.0d-3, 200.0d0,&  515.0d-3, 300.0d0,&  536.0d-3, 400.0d0,&  557.0d-3, 500.0d0,&     586.0d-3, 600.0d0,&  619.0d-3, 700.0d0,& 691.0d-3, 800.0d0,&  695.0d-3, 900.0d0,& 691.0d-3, 100.0d1,&  700.0d-3, 300.0d1/

      t_c=temp-273.15d0
      if (t_c .le. t(1)) then
         cp_steel1=value(1)
      elseif (t_c .ge. t(11)) then
         cp_steel1=value(11)
      else
         do i=1,10
            if ((t_c .gt. t(i)) .and. (t_c .le. t(i+1))) exit
         enddo
         cp_steel1=value(i)+(value(i+1)-value(i))&        *(t_c-t(i))/(t(i+1)-t(i))
      endif

      return
      end

      function cp_steel2(temp)
      implicit real *8 (a-h,o-z)

      parameter (np=26)
      dimension t(np),value(np)

      data (value(i),t(i),i=1,np)/
     &  5.104E-01,  373., &   5.272E-01, 423.,&  5.314E-01,   473.,&  5.397E-01,   523.,&   5.481E-01,   573.,&    5.607E-01,   623.,&    5.690E-01,   673.,&  5.816E-01,   723.,&   5.941E-01,   773.,
     &  6.276E-01,   823.,&   6.485E-01,   873.,&      6.318E-01,   923.,&  6.234E-01,   973.,&  6.276E-01,  1023.,&  6.402E-01,  1073.,&  6.485E-01,  1123.,&   6.402E-01,  1173.,&      6.485E-01,  1223.,&   6.485E-01,  1273.,&   6.527E-01,  1323.,&  6.611E-01,  1373.,&   6.653E-01,  1423.,&   6.736E-01,  1473.,&   6.778E-01,  1523.,&  6.820E-01,  1573.,&   7.000E-01,  3000./

      if (temp .le. t(1)) then
         cp_steel2 = value(1)
      elseif (temp .ge. t(np)) then
         cp_steel2 = value(np)
      else
         do i=1,np-1
            if ((temp .gt. t(i)) .and. (temp .le. t(i+1))) exit
         enddo
         fac = (temp-t(i))/(t(i+1)-t(i))
         cp_steel2 = value(i) + (value(i+1)-value(i))*fac
      endif

      return
      end

      function cp_steel3(temp)
      implicit real *8 (a-h,o-z)
      dimension t(11),value(11)

      data (t(i),value(i),i=1,11) /
     &     100.0d0, 486.0d-3,
     &     200.0d0, 498.0d-3,
     &     300.0d0, 515.0d-3,
     &     400.0d0, 536.0d-3,
     &     500.0d0, 557.0d-3,
     &     600.0d0, 586.0d-3,
     &     700.0d0, 619.0d-3,
     &     800.0d0, 691.0d-3,
     &     900.0d0, 695.0d-3,
     &     100.0d1, 691.0d-3,
     &     300.0d1, 700.0d-3/

      t_c=temp-273.15d0
      if (t_c .lt. t(1)) then
         cp_steel3=value(1)
      elseif (t_c .ge. t(8)) then
         cp_steel3=value(8)
      else
         do i=1,7
            if ((t_c .gt. t(i)) .and. (t_c .le. t(i+1))) exit
         enddo
         cp_steel3=value(i)+(value(i+1)-value(i)) &        *(t_c-t(i))/(t(i+1)-t(i))
      endif

      return
      end

      function dens(ID)
      implicit real *8 (a-h,o-z)

      if (id .eq. 1)  then      ! carbon steel
         dens = 7.820E-03
      elseif (id .eq. 2) then   ! austenitic stainless steel
         dens = 7.820E-03
      elseif (id .eq. 3) then   ! duplex stainless steel
         dens = 7.820E-03
      endif

      end
      subroutine h(xx,temp,h_con,id)
      implicit real*8 (a-h,o-z)

      if (id .eq. 1) then       ! air
         call h_air(xx,temp,h_con)
      elseif (id .eq. 2) then   ! gas methane
         call h_methane1(xx,temp,h_con)
      elseif (id .eq. 3) then   ! hydrogen
         call h_hydrogen(xx,temp,h_con)
      elseif (id .eq. 4) then   ! steam
         call h_steam(xx,temp,h_con)
      elseif (id .eq. 5) then   ! ethylene
         call h_ethyl1(xx,temp,h_con)

      elseif (id .eq. 11) then   ! ethyl alcohol
         call h_ethyl2(xx,temp,h_con)
      elseif (id .eq. 12) then   ! water
         call h_water(xx,temp,h_con)
      elseif (id .eq. 13) then   ! ammonia
         call h_ammonia(xx,temp,h_con)
      elseif (id .eq. 14) then   ! propane
         call h_propane(xx,temp,h_con)
      elseif (id .eq. 15) then  ! fuel oil
         call h_foil(xx,temp,h_con)
      elseif (id .eq. 16) then  ! gasoline
         call h_gasoline(xx,temp,h_con)
      elseif (id .eq. 17) then  ! liquid methane
         call h_methane2(xx,temp,h_con)
      elseif (id .eq. 18) then  ! crude oil
         call h_coil(xx,temp,h_con)
      elseif (id .eq. 19) then  ! motor oil
         call h_moil(xx,temp,h_con)
      else
         print *, 'FLOW MATERIAL PROPERTIES NOT AVAILABLE.'
         print *, 'ANALYSIS STOPED!'
      endif

      return
      end

      subroutine h_air(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'
      integer iflag
      save iflag,pr,R_e,c_p,density,viscosity,tk_air,beta

      d_in_m=d_in/1.0d3
      if (iflag .eq. 0) then
         iflag=1

         beta=1.0d0/t0_f
         pr=0.7d0
         c_p=1.0d3
         viscosity=0.01716d-3*(t0_f/273.16d0)**1.5d0&        *383.72d0/(t0_f+110.56d0)
         density=pressure*29/8134.0d0/t0_f
         R_e=density*d_in_m*velocity/viscosity
         tk_air=viscosity*c_p/pr
      endif

      if (R_e .gt. 5000) then

         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.318d0*(d_in/xx)**0.75d0
         h_f=factor*vnu*tk_air/d_in_m
      else if (r_e.le.2300) then

         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_air/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.318d0*(d_in/xx)**0.75d0
         h1=factor*vnu*tk_air/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_air/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      endif

      delt=temp-t0_f
      if (delt .le. 0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr .ge. 1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_air/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_air/d_in_m
         endif
      endif
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_ammonia(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      pr=2.0d0
      beta=2.45d-3
      density=703.69d0+(t0_f-273.15d0-50.0d0) &     *(564.33d0-703.69d0)/100.0d0
      c_p=4.635d3
      tk_ammonia=0.547d0+(t0_f-273.15d0-50.0d0) &     *(0.476d0-0.547d0)/100.0d0
      viscosity=pr*tk_ammonia/c_p
      r_e=density*d_in_m*velocity/viscosity
      if (R_e .gt. 5000) then

         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_ammonia/d_in_m
      else if (r_e.le.2300) then

         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_ammonia/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_ammonia/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_ammonia/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_ammonia/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_ammonia/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_coil(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      beta=0.70d-3
      density=899.00d0+(t0_f-273.15d0)&     *(805.89d0-889.00d0)/160.0d0
      c_p=2.000d3
      tk_coil=0.140d0
      viscosity=0.00428d0+(t0_f-273.15d0)&     *(0.056d-4-0.00428d0)/160.0d0
      r_e=density*d_in_m*velocity/viscosity
      pr=viscosity*c_p/tk_coil
      if (R_e.gt.5000) then     ! turbulent flow
                                ! vnu=0.023*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_coil/d_in_m
      else if (r_e.le.2300) then ! laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_coil/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_coil/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_coil/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_coil/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_coil/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_ethyl1(xx,temp,h)
      implicit real*8 (a-h,o-z)
      include 'common.f'

      integer iflag
      save iflag, density, c_p, tk_C2H4, viscosity, R_e, Pr, beta

      d_in_m = d_in / 1.0D+3
      if (iflag .eq. 0) then
        iflag = 1
        beta = 1.0D0 / t0_f
        density = 28.0D0*pressure / (8.314D+3*t0_f)
        c_p = 1.535D+3
        Pr = 0.75D0
        viscosity = 14.3D-6
        tk_C2H4 = c_p * viscosity / Pr
        R_e = density * d_in_m * velocity / viscosity
      end if

      if (R_e > 5000.0D0) then
        vnu = 0.023D0 * (R_e**0.8D0) * (Pr**0.33333D0)
        h_f = vnu * tk_C2H4 / d_in_m
      elseif (r_e <= 2300.0D0) then
        vnu = 3.65D0 + (0.0668D0*d_in*R_e*Pr/xx)
     1      /(1.0D0+0.04D0*(d_in*R_e*Pr/xx)**0.66667D0)
        h_f = vnu * tk_C2H4 / d_in_m
      else
        vnu = 0.023D0 * (R_e**0.8D0) * (Pr**0.33333D0)
        h1 = vnu * tk_C2H4 / d_in_m
        vnu = 3.65D0 + (0.0668D0*d_in*R_e*Pr/xx)
     1      /(1.0D0+0.04D0*(d_in*R_e*Pr/xx)**0.66667D0)
        h2 = vnu * tk_C2H4 / d_in_m
        h_f = h2 + (h1-h2)*(R_e-2300.0D0)/(5000.0D0-2300.0D0)
      endif

      delt = temp - t0_f
      if (delt <= 0.0D0) then
        gr = 0.0D0
        h_n = 0.0D0
      else
        gr = d_in_m**3.0D0 * beta * delt * density *
     1      9.8D0 * density / viscosity / viscosity
        if (gr >= 1.0D+9) then
            vnu = 0.16D0 * (gr*Pr)**0.33333d0
            h_n = vnu * tk_C2H4 / d_in_m
         else
            vnu = 0.59D0 * (gr*Pr)**0.25D0
            h_n = vnu * tk_C2H4 / d_in_m
         end if
      end if

      h = (h_f*h_f*h_f + h_n*h_n*h_n)**0.33333D0
      h = h*1.0D-6
      return
      end

      subroutine h_ethyl2(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      beta=1.10d-3
      density=806.00d0+(t0_f-273.15d0) &     *(383.00d0-806.00d0)/240.0d0
      c_p=1.884d3+(4664.0d0-1884.0d0) &     *(t0_f-373.15d0)/260.0d0
      tk_ethyl2=0.180d0
      viscosity=470.0d-4+(t0_f-273.15d0)&     *(5.03d-4-470.0d-4)/170.0d0
      pr=viscosity*c_p/tk_ethyl2
      r_e=density*d_in_m*velocity/viscosity
      if (R_e.gt.5000) then     !turbulent flow
                                ! vnu=0.023*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_ethyl2/d_in_m
      else if (r_e.le.2300) then !laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_ethyl2/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_ethyl2/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_ethyl2/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_ethyl2/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_ethyl2/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_foil(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'
      integer iflag
      save iflag

      if (iflag .eq. 0) then
         iflag=1
         density=850.0d0
         c_p=200.0d0
         tk=0.14d0
      end if

      h=2300.0d-6

      return
      end

      subroutine h_gasoline(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      beta=1.14d-3
      density=718.00d0+(t0_f-273.15d0)&     *(337.00d0-718.00d0)/290.0d0
      c_p=2.177d3
      tk_gasoline=0.140d0
      viscousity=7.00d-4+(t0_f-273.15d0)&     *(2.08d-4-7.0d-4)/120.0d0
      pr=viscousity*c_p/tk_gasoline
      r_e=density*d_in_m*velocity/viscousity
      if (R_e.gt.5000) then     ! turbulent flow
                                ! vnu=0.023*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_gasoline/d_in_m
      else if (r_e.le.2300) then ! laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_gasoline/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_gasoline/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_gasoline/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density* &        9.8d0*density/viscousity/viscousity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_gasoline/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_gasoline/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_hydrogen(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      beta=1.0d0/t0_f
      pr=0.7d0
      c_p=15.5d3
      viscosity=0.008411d-3*(t0_f/273.16d0)**1.5d0 &     *369.81d0/(t0_f+96.667d0)
      density=pressure*2/8134.0d0/t0_f
      R_e=density*d_in_m*velocity/viscosity
      tk_hyd=viscosity*c_p/pr
      if (R_e.gt.5000) then     ! turbulent flow
                                ! vnu=0.023*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.318d0*(d_in/xx)**0.75d0
         h_f=factor*vnu*tk_hyd/d_in_m
      else if (r_e.le.2300) then ! laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_hyd/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.318d0*(d_in/xx)**0.75d0
         h1=factor*vnu*tk_hyd/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_hyd/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_hyd/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_hyd/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_methane1(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'
      integer, save :: iflag=0
      save pr,R_e,c_p,density,viscosity,tk_methane,beta

      d_in_m=d_in/1.0d3
      if (iflag .eq. 0) then
         iflag=1

         beta=1.0d0/temp
         c_p=2165.4d0+(t0_f-273.0d0)*3202.0d0/1200.0d0
         viscosity=7.4d-6+(t0_f-193.0d0)*15.24d-6/579.0d0
         density=pressure*16./8134.0d0/t0_f
         tk_methane=0.01291d0+(t0_f-123.15d0)*0.02431d0/2.0d2
         R_e=density*d_in_m*velocity/viscosity
         pr=viscosity*c_p/tk_methane
      end if

      if (R_e.gt.5000) then     ! turbulent flow

         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         h_f=vnu*tk_methane/d_in_m
         h_f=h_f*(1.0d0+0.318d0*(d_in/xx)**0.75d0)
      else if (r_e.le.2300) then ! laminar flow

         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_methane/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         h1=vnu*tk_methane/d_in_m
         h1=h1*(1.0d0+0.318d0*(d_in/xx)**0.75d0)
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_methane/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_methane/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_methane/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      viscosity_w=7.4d-6+(temp-193.0d0)*15.24d-6/579.0d0
      h=h*1.0d-6*(viscosity/viscosity_w)**0.14d0

      return
      end

      subroutine h_methane2(x,temp,h)
      implicit real *8 (a-h,o-z)
      PARAMETER (Ta=536.7,Tb=534.9,g=32.17,D=0.8645,V=16.9,
     $     Pf=600.0,a11=32.654,a21=-2.43437E-1,a31=7.01127E-4,
     $     a41=-8.96041E-7,a51=4.28695E-10,     a12=-9.24685E-2,
     $     a22=7.04526E-4,a32=-2.02503E-6,a42=2.5938E-9,
     $     a52=-1.24701E-12,  a13=3.88139E-5,a23=-3.01735E-7,
     $     a33=8.8031E-10,a43=-1.141E-12,a53=5.54013E-16)

      t=temp*1.8+32.0+497.5
      f1=a11+a21*Tb+a31*Tb**2.0+a41*Tb**3.0+a51*Tb**4.0
      f2=a12+a22*Tb+a32*Tb**2.0+a42*Tb**3.0+a52*Tb**4.0
      f3=a13+a23*Tb+a33*Tb**2.0+a43*Tb**3.0+a53*Tb**4.0
      Z=f1+f2*(Pf+14.7)+f3*(Pf+14.7)**2.0
      Cp=0.4366E-3*Tb+0.3004
      cod=0.1255E-7*Tb**2.0+0.3214E-4*Tb-0.6064E-3
      u=4.966E-7*Tb**1.5/(Tb+295.2)
      belta=2.0/(T+Tb)
      rou=144.0*(Pf+14.7)/(96.4*Tb*Z)
      uw=4.966E-7*T**1.5/(T+295.2)
      Re=rou*V*D/u
      Pr=3600.0*u*Cp/cod

      IF (Re .GE. 5000.0) THEN
         C1=3600.0*Cp*V*rou*0.023*(Re**-0.2)*
     $        (Pr**(-2.0/3.0))*(u/uw)**0.14
         Hf=(1+0.318*(D/X)**0.75)*C1
      ELSE IF (Re .LE. 2100.0) THEN
         C2=1.615*(Re*Pr/(X/D))**(1.0/3.0)
         Hf=(cod/D)*((C2**5.0+4.364**5)**0.2)*(u/uw)**0.14
      ELSE
         C1=3600.0*Cp*V*rou*0.023*(5000.0**-0.2)*
     $        (Pr**(-2.0/3.0))*(u/uw)**0.14
         HT=(1+0.318*(D/X)**0.75)*C1
         C2=1.615*(2100.0*Pr/(X/D))**(1.0/3.0)
         HL=(cod/D)*((C2**5.0+4.364**5)**0.2)*(u/uw)**0.14
         Hf=HL+(Re-2100.0)*(HT-HL)/2900.0
      END IF

      Tc=(T+Tb)/2.0
      f1n=a11+a21*Tc+a31*Tc**2.0+a41*Tc**3.0+a51*Tc**4.0
      f2n=a12+a22*Tc+a32*Tc**2.0+a42*Tc**3.0+a52*Tc**4.0
      f3n=a13+a23*Tc+a33*Tc**2.0+a43*Tc**3.0+a53*Tc**4.0
      Zn=f1n+f2n*(Pf+14.7)+f3n*(Pf+14.7)**2.0
      Cpn=0.4366E-3*Tc+0.3004
      con=0.1255E-7*Tc**2.0+0.3214E-4*Tc-0.6064E-3
      un=4.966E-7*Tc**1.5/(Tc+295.2)
      ron=144.0*(Pf+14.7)/(96.4*Tc*Zn)
      Gr=(D**3.0)*g*belta*(T-Tb)*(ron/un)**2.0
      Prn=3600.0*un*Cpn/con
      GP=Gr*Prn
      IF (GP .GT. 3.0E10) THEN
         an=0.5
         C=0.0362
      ELSE
         an=0.25
         C=0.27
      END IF

      Hn=(con/D)*C*GP**an

      h=5.6788E-6*(Hf**3.0+Hn**3.0)**(1.0/3.0)

      h=10000.0d-6
      return
      end

      subroutine h_moil(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'

      d_in_m=d_in/1.0d3
      beta=0.70d-3
      density=899.00d0+(t0_f-273.15d0)&     *(805.89d0-889.00d0)/160.0d0
      c_p=2.000d3
      tk_moil=0.140d0
      viscosity=7962.2d-4+(t0_f-293.15d0)&     *(65.02d-4-7962.2d-4)/120.0d0
      r_e=density*d_in_m*velocity/viscosity
      pr=viscosity*c_p/tk_moil
      if (R_e.gt.5000) then

         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_moil/d_in_m
      else if (r_e.le.2300) then

         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_moil/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_moil/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_moil/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_moil/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_moil/d_in_m
         end if
      end if

      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_propane(xx,temp,h)
      implicit real *8 (a-h,o-z)

      h=230.0d-6

      return
      end

      subroutine h_steam(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'
      integer iflag
      save iflag,pr,R_e,c_p,density,viscosity,tk_steam,beta

      d_in_m=d_in/1.0d3
      if (iflag.eq.0) then
         iflag=1
         beta=1.0d0/t0_f
         c_p=2.0d3
         pr=0.85d0
         viscosity=134.0d-7+(t0_f-400.0d0)*144.2d-7/4.0d2
         density=pressure*18/8134.0d0/t0_f
         R_e=density*d_in_m*velocity/viscosity
         tk_steam=viscosity*c_p/pr
      end if
      if (R_e.gt.5000) then     ! turbulent flow
                                ! vnu=0.026*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         h_f=vnu*tk_steam/d_in_m
      else if (r_e.le.2300) then ! laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_steam/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         h1=vnu*tk_steam/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_steam/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density* &        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_steam/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_steam/d_in_m
         end if
      end if
      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine h_water(xx,temp,h)
      implicit real *8 (a-h,o-z)
      include 'common.f'
      integer iflag
      save iflag,pr,R_e,c_p,density,viscosity,tk_water,beta

      d_in_m=d_in/1.0d3
      if (iflag.eq.0) then
         iflag=1

         beta=1.8d-4
         c_p=4.2d3
         t_c=t0_f-273.16d0-8.435d0
         vvv=0.021482d0*(t_c+(t_c*t_c+8078.4)**0.5d0)-1.2d0
         viscosity=1.0d-3/vvv
         density=1.0d3
         tk_water=0.65d0
         pr=c_p*viscosity/tk_water
         R_e=density*d_in_m*velocity/viscosity
      end if
      if (R_e.gt.5000) then     ! turbulent flow
                                ! vnu=0.026*Re^0.8*Pr^0.333
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h_f=1.0d0*vnu*tk_water/d_in_m
      else if (r_e.le.2300) then ! laminar flow
                                ! vnu=3.65+...........
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx) &        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h_f=vnu*tk_water/d_in_m
      else
         vnu=0.023d0*r_e**0.8d0*pr**0.33333d0
         factor=1.0d0+0.125d0*(d_in/xx)
         h1=factor*vnu*tk_water/d_in_m
         vnu=3.65d0+(0.0668d0*d_in*r_e*pr/xx)&        /(1.0d0+0.04d0*(d_in*r_e*pr/xx)**0.66667d0)
         h2=vnu*tk_water/d_in_m
         h_f=h2+(h1-h2)*(r_e-2300.0d0)/(5000.0d0-2300.0d0)
      end if

      delt=temp-t0_f
      if (delt.le.0.0d0) then
         gr=0.0d0
         h_n=0.0d0
      else
         gr=d_in_m**3*beta*delt*density*&        9.8d0*density/viscosity/viscosity
         if (gr.ge.1.0d9) then
            vnu=0.16d0*(gr*pr)**0.333333d0
            h_n=vnu*tk_water/d_in_m
         else
            vnu=0.59d0*(gr*pr)**0.25d0
            h_n=vnu*tk_water/d_in_m
         end if
      end if

      h=(h_f*h_f*h_f+h_n*h_n*h_n)**0.33333d0
      h=h*1.0d-6

      return
      end

      subroutine mesh00(thickness,inode,coord,nnode,&     ielem,inelem,nelem,kp)

      implicit double precision (a-h,o-z)

      dimension :: inode(*),coord(2,*)
      dimension :: ielem(4,*),inelem(4,*)
      dimension :: kp(*)

      pi=acos(-1.0)
      ndiv=16
      dxmin=25.4/ndiv           ! minimum element size in radial direction
      dymin=1.                  ! minimum element size in thickness direction
      plength=40.*thickness
      ndy=nint(thickness/dymin)
      dy=thickness/ndy
      dxmax=20.*dxmin
      ndiv1=ndiv+4

      dist1=ndiv1*dxmin
      q=1.25
      ndiv2=int(log(dxmax/dxmin)/log(q))
      dist2=dxmin*q*(q**ndiv2-1.)/(q-1.)

      if (dist1+dist2 .gt. plength) then
         dist3=0.
         dist2=plength-dist1
         ndiv2=int(log(1.+dist2*(q-1.)/dxmin/q)/log(q))
      else
         dist3=plength-dist1-dist2
         ndiv3=1+int(dist3/dxmax)
      endif

      kp(1)=1+ndy
      kp(2)=(1+ndiv)*(1+ndy)

      do j=1,1+ndiv1+ndiv2+ndiv3
         if (j .le. 1+ndiv1) then
            xx=(j-1)*dxmin
         elseif (j .le. 1+ndiv1+ndiv2) then
            jj=j-1-ndiv1
            xx=ndiv1*dxmin+dxmin*q*(q**jj-1.)/(q-1.)
         else
            jj=j-1-ndiv1-ndiv2
            xx=dist1+dist2+jj*dist3/ndiv3
         endif
         do i=1,1+ndy
            nd=(j-1)*(1+ndy)+i
            inode(nd)=nd
            coord(1,nd)=xx
            coord(2,nd)=(i-1)*thickness/ndy
c            write(1,3) nd,coord(1,nd),coord(2,nd)
         enddo
      enddo
 3    format(i8,2(:', 'ES12.5))

      do j=1,ndiv1+ndiv2+ndiv3
         do i=1,ndy
            iel=(j-1)*ndy+i
            ielem(1,iel)=iel
            ielem(2,iel)=1
            if (i .eq. 1) then
               ielem(3,iel)=2
               ielem(4,iel)=1
            elseif (i .eq. ndy) then
               ielem(3,iel)=1
               ielem(4,iel)=3
            else
               ielem(3,iel)=0
               ielem(4,iel)=0
            endif

            inelem(1,iel)=(j-1)*(1+ndy)+i
            inelem(2,iel)=inelem(1,iel)+1+ndy
            inelem(3,iel)=inelem(2,iel)+1
            inelem(4,iel)=inelem(1,iel)+1
         enddo
      enddo
 5    format(1x,i4,4(:', 'i4))

      nnode=nd
      nelem=iel

      return
      end
      subroutine mesh01(d_pipe,t_pipe,t_sleeve,gap,weldarea,
     1     angle_w,vl_pipe,nnode,nelem,inode,
     2     coord,larray_ele,key_point,r_cc)

      implicit double precision (a-h,o-z)

      dimension :: key_point(*),larray_ele(8,*)
      dimension :: inode(*),coord(2,*)
      dimension :: x(1000),y(1000),n_repeat(10,2)
      dimension :: co_f(2),co_l(2),dist(4),larray2D(1000,5)
      common /coord/x,y

      meshtest = 0

      if(oldmesher.eq.0) meshtest=1

      bf_t_pipe=1.3
      bf_sleeve_top=1.3
      bf_pipe_left=1.5
      bf_pipe_right=1.5
      ne_pipe=9
      ne_pipe_left=19
      ne_pipe_right=19
      ne_sleeve_top=5

      pi=acos(-1.0d0)
      to_rad = pi/180.d0
      ang_set = 30.d0*to_rad
      angle_w = angle_w*to_rad

      weldxx=sqrt(2.0*weldarea)
      if (weldxx .gt. (gap+7.*t_sleeve/8.)) then
         weldxx=t_sleeve+gap
      endif
      weldyy=2.0*weldarea/weldxx

      r_cc=0.5*sqrt(weldxx*weldxx+weldyy*weldyy)
      write(*,3) 'Estimated cross-sectional area of weld deposit =', &     weldarea,'mm^2'
      write(*,3) 'Meshing in progress!'
 3    format(1x,A,F6.2,1x,A)

      key_point(1)=1
      key_point(2)=5
      key_point(3)=9
      key_point(4)=81

      x(1) = 0.d0
      y(1) = 0.d0
      x(9) = -weldyy
      y(9) = -weldxx
      x(81) = x(1)
      y(81) = y(9)

      x(5) = 0.5d0*(x(1)+x(9))
      y(5) = 0.5d0*(y(1)+y(9))
      x(41) = 0.5d0*(x(1)+x(81))
      y(41) = 0.5d0*(y(1)+y(81))

      if (ang_set .lt. angle_w) ang_set=angle_w
      x(45) = (dtan(ang_set)*x(41)+x(5)-y(41)+y(5)) &       /(dtan(ang_set)+1.d0)
      y(45) = y(5)-x(45)+x(5)
      x(49) = x(45)
      y(49) = y(9)

      nb = 1
      ne = 81
      nd = 10
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ng_line(nb,ne,nd,co_f,co_l,bf)

      do 20 i=1,4
         dist(i) = dabs(dabs(y(71-10*(i-1))-y(81))-gap)
 20   continue
      dist_mini = t_sleeve
      do 40 i=1,4
         if (dist(i) .lt. dist_mini) then
            dist_mini = dist(i)
            id_mini = i
         endif
 40   continue
      key_point(5) = 71-10*(id_mini-1)
      key_point(6) = 6
      key_point(7) = 7
      key_point(8) = 8

      if (gap .ge. 0.5d0*dabs(y(71)-y(81))) then
         y(71-10*(id_mini-1)) = -weldxx+gap
      else
         y(71) = -weldxx+0.5d0*dabs(y(71)-y(81))
      end if
      nb = 1
      ne = 9
      nd = 1
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ng_line(nb,ne,nd,co_f,co_l,bf)
      nb = 9
      ne = 49
      nd = 10
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ng_line(nb,ne,nd,co_f,co_l,bf)
      nb = 5
      ne = 45
      nd = 10
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ng_line(nb,ne,nd,co_f,co_l,bf)

      do 60 i=1,4
         nb = 11+(i-1)*10
         ne = 15+(i-1)*10
         nd = 1
         bf = 1.0
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 60   continue

      do 80 i=1,4
         nb = 15+(i-1)*10
         ne = 19+(i-1)*10
         nd = 1
         bf = 1.0
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 80   continue

      x(55) = x(46)
      y(55) = y(46)
      x(65) = x(47)
      y(65) = y(47)
      x(75) = x(48)
      y(75) = y(48)
      x(85) = x(49)
      y(85) = y(49)

      do 100 i=1,4
         nb = 51+(i-1)*10
         ne = 55+(i-1)*10
         nd = 1
         bf = 1.0
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 100  continue


      do k=1,5
         nd_delta = 10*(k-1)
         nd_count_delta = 9*(k-1)
         do i=1+nd_delta,9+nd_delta
            nd_count = i-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*d_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count
      do k=1,4
         nd_delta = 10*(k-1)
         nd_count_delta = 4*(k-1)
         do i=51+nd_delta,54+nd_delta
            nd_count = nnode+i-50-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*d_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count

      do 120 k=1,4
         n_repeat(k,1) = 55+10*(k-1)
         n_repeat(k,2) = 46+(k-1)
 120  continue

      n_ele = 1
      nd1 = 1
      nd2 = nd1+1
      nd3 = nd2+10
      nd4 = nd1+10
      nele1row = 8
      nddt1row = 1
      nedt1row = 1
      nelerows = 4
      nddtrows = 10
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,&     nedt1row,nelerows,nddtrows,nedtrows,larray2D)
      id_mat = 3
      do 140 j=1,32
         if (j .le. 8) then
            id_location = 1
            id_face = 1
         else
            id_location = 0
            id_face = 0
         end if
         do m=1,5
            larray_ele(m,j) = larray2D(j,m)
         enddo
         larray_ele(6,j) = id_mat
         larray_ele(7,j) = id_location
         larray_ele(8,j) = id_face
 140  continue

      n_ele = nele1row*nelerows+n_ele
      nd1 = 41
      nd2 = nd1+1
      nd3 = nd2+10
      nd4 = nd1+10
      nele1row = 4
      nddt1row = 1
      nedt1row = 1
      nelerows = 4
      nddtrows = 10
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row, &     nedt1row,nelerows,nddtrows,nedtrows,larray2D)
      do 160 j=1,16
      do 160 i=2,5
      do 160 k=1,4
         if (larray2D(j,i) .eq. n_repeat(k,1)) then
            larray2D(j,i)=n_repeat(k,2)
         endif
160   continue

      do 180 j=1,16
         if (id_mini .eq. 1 .and. j .eq. 13) then
            id_location = 3
            id_face = 4
         elseif (id_mini .eq. 2 .and. (j .eq. 9 .or. j .eq. 13)) then
            id_location = 3
            id_face = 4
         elseif (id_mini .eq. 3 .and. &           (j .eq. 5 .or. j .eq. 9 .or. j .eq. 13)) then
            id_location = 3
            id_face = 4
         elseif ((j .eq. 1 .or. j .eq. 5 .or. j .eq. 9 .or. j .eq. 13)&           .and. (id_mini .eq. 4)) then
            id_location = 3
            id_face = 4
         else
            id_location = 0
            id_face = 0
         endif
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 180  continue

200   dr=1.d0
      do k=1,ne_pipe-1,1
         dr=dr+bf_t_pipe**k
      enddo
      w0=t_pipe/dr
      if (w0 .lt. abs(y(48)-y(49))) then
         ne_pipe = ne_pipe-1
         goto 200
      end if

      nd_grade = 10*ne_pipe
      key_point(9) = 101+nd_grade
      key_point(10) = 105+nd_grade
      key_point(11) = 109+nd_grade

      do i=1,5
         x(100+i) = x(9+(i-1)*10)
         y(100+i) = y(9+(i-1)*10)
      enddo
      do i=1,4
         x(105+i) = x(85-i)
         y(105+i) = y(85-i)
      enddo

      do i=1,9
         x(100+i+ne_pipe*10) = x(100+i)
         y(100+i+ne_pipe*10) = y(100+i)-t_pipe
      enddo

      do 220 i=1,9
         nb = 100+i
         ne = nb+ne_pipe*10
         nd = 10
         bf = bf_t_pipe
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 220  continue

      do k=1,5
         n_repeat(k,1) = 101+(k-1)
         n_repeat(k,2) = 9+10*(k-1)
      enddo
      do k=6,9
         n_repeat(k,1) = 101+(k-1)
         n_repeat(k,2) = 84-(k-6)
      enddo

      do k=1,ne_pipe
         nd_delta = 10*(k-1)
         nd_count_delta = 9*(k-1)
         do i=111+nd_delta,119+nd_delta
            nd_count = nnode+i-110-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*d_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count
      n_ele = nele1row*nelerows+n_ele
      nd1 = 101
      nd2 = nd1+10
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = ne_pipe
      nddt1row = 10
      nedt1row = 1
      nelerows = 8
      nddtrows = 1
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&           nelerows,nddtrows,nedtrows,larray2D)
      do 240 j=1,nele1row*nelerows
      do 240 i=2,5
      do 240 k=1,9
         if (larray2D(j,i) .eq. n_repeat(k,1)) then
            larray2D(j,i)=n_repeat(k,2)
         endif
 240  continue
c
      id_mat = 1
      do 260 j=1,nele1row*nelerows
         if (mod(j,ne_pipe) .eq. 0) then
            id_location = 2
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 260  continue

 300  dr=1.d0
      do k=1,ne_pipe_left-1,1
         dr=dr+bf_pipe_left**k
      enddo
      w0=(0.5d0*vl_pipe-weldyy)/dr
      if (w0 .lt. abs(x(101)-x(102))) then
         ne_pipe_left = ne_pipe_left-1
         goto 300
      end if

      do i=1,ne_pipe+1
         x(200+1+20*(i-1)) = x(100+1+10*(i-1))
         y(200+1+20*(i-1)) = y(100+1+10*(i-1))
      enddo

      do i=1,ne_pipe+1
         x(200+1+ne_pipe_left+20*(i-1)) = -0.5d0*vl_pipe
         y(200+1+ne_pipe_left+20*(i-1)) = y(100+1+10*(i-1))
      enddo

      do 320 i=1,ne_pipe+1
         nb = 200+1+20*(i-1)
         ne = nb+ne_pipe_left
         nd = 1
         bf = bf_pipe_left
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 320  continue

      n_repeat(1,1) = 201
      n_repeat(1,2) = 9
      do k=2,ne_pipe+1
         n_repeat(k,1) = 201+20*(k-1)
         n_repeat(k,2) = 101+10*(k-1)
      enddo

      do k=1,ne_pipe+1
         nd_delta = 20*(k-1)
         nd_count_delta = ne_pipe_left*(k-1)
         do i=202+nd_delta,202+(ne_pipe_left-1)+nd_delta
            nd_count = nnode+i-201-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*D_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count

      n_ele = nele1row*nelerows+n_ele
      nd1 = 201
      nd2 = nd1+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_pipe_left
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_pipe
      nddtrows = 20
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&           nelerows,nddtrows,nedtrows,larray2D)
      do 340 j=1,nele1row*nelerows
      do 340 i=2,5
      do 340 k=1,ne_pipe+1
         if (larray2D(j,i) .eq. n_repeat(k,1)) then
            larray2D(j,i)=n_repeat(k,2)
         endif
 340  continue

      id_mat = 1
      do 360 j=1,nele1row*nelerows
         if (j .le. ne_pipe_left) then
            id_location = 1
            id_face = 1
         else if (j .gt. ne_pipe_left*(ne_pipe-1)) then
            id_location = 2
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 360  continue

 400  dr=1.d0
      do k=1,ne_pipe_right-1,1
         dr=dr+bf_pipe_right**k
      enddo
      w0=(0.5d0*vl_pipe-weldyy*dtan(angle_w))/dr
      if (w0 .lt. abs(x(108)-x(109))) then
         ne_pipe_right = ne_pipe_right-1
         goto 400
      end if

      do i=1,ne_pipe+1
         x(400+1+20*(i-1)) = x(100+9+10*(i-1))
         y(400+1+20*(i-1)) = y(100+9+10*(i-1))
      enddo

      do i=1,ne_pipe+1
         x(400+1+ne_pipe_right+20*(i-1)) = 0.5d0*vl_pipe
         y(400+1+ne_pipe_right+20*(i-1)) = y(100+9+10*(i-1))
      enddo

      do 420 i=1,ne_pipe+1
         nb = 400+1+20*(i-1)
         ne = nb+ne_pipe_right
         nd = 1
         bf = bf_pipe_right
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 420  continue

      n_repeat(1,1) = 401
      n_repeat(1,2) = 81
      do k=2,ne_pipe+1
         n_repeat(k,1) = 401+20*(k-1)
         n_repeat(k,2) = 109+10*(k-1)
      enddo

      do k=1,ne_pipe+1
         nd_delta = 20*(k-1)
         nd_count_delta = ne_pipe_right*(k-1)
         do i=402+nd_delta,402+(ne_pipe_right-1)+nd_delta
            nd_count = nnode+i-401-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*D_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count

      n_ele = nele1row*nelerows+n_ele
      nd2 = 401
      nd1 = nd2+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_pipe_right
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_pipe
      nddtrows = 20
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row, &           nelerows,nddtrows,nedtrows,larray2D)
      do 440 j=1,nele1row*nelerows
      do 440 i=2,5
      do 440 k=1,ne_pipe+1
         if (larray2D(j,i) .eq. n_repeat(k,1)) then
            larray2D(j,i)=n_repeat(k,2)
         endif
 440  continue

      id_mat = 1
      do 460 j=1,nele1row*nelerows
         if (j .le. ne_pipe_right) then
            id_location = 3
            id_face = 1
         else if (j .gt. ne_pipe_right*(ne_pipe-1)) then
            id_location = 2
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 460  continue
      do i=1,8-(id_mini-1)
         x(601+20*(i-1)) = x(71-10*(id_mini-1)-10*(i-1))
         y(601+20*(i-1)) = y(71-10*(id_mini-1)-10*(i-1))
      enddo

      do i=1,8-(id_mini-1)
         x(601+ne_pipe_right+20*(i-1)) = 0.5d0*vl_pipe
         y(601+ne_pipe_right+20*(i-1)) = y(71-10*(id_mini-1)-10*(i-1))
      enddo

      do 520 i=1,8-(id_mini-1)
         nb = 601+20*(i-1)
         ne = nb+ne_pipe_right
         nd = 1
         bf = bf_pipe_right
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 520  continue

      do k=1,8-(id_mini-1)
         n_repeat(k,1) = 601+20*(k-1)
         n_repeat(k,2) = 71-10*(id_mini-1)-10*(k-1)
      enddo

      do k=1,9-id_mini
         nd_delta = 20*(k-1)
         nd_count_delta = ne_pipe_right*(k-1)
         do i=602+nd_delta,602+(ne_pipe_right-1)+nd_delta
            nd_count = nnode+i-601-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*D_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count

      n_ele = nele1row*nelerows+n_ele
      nd1 = 601
      nd2 = nd1+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_pipe_right
      nddt1row = 1
      nedt1row = 1
      nelerows = 7-(id_mini-1)
      nddtrows = 20
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&           nelerows,nddtrows,nedtrows,larray2D)
      do 540 j=1,nele1row*nelerows
      do 540 i=2,5
      do 540 k=1,8-(id_mini-1)
         if (larray2D(j,i) .eq. n_repeat(k,1)) then
            larray2D(j,i)=n_repeat(k,2)
         endif
 540  continue

      id_mat = 2
      do 560 j=1,nele1row*nelerows
         if (j .le. ne_pipe_right) then
            id_location = 3
            id_face = 1
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 560  continue

 600  dr=1.d0
      do k=1,ne_sleeve_top-1,1
         dr=dr+bf_sleeve_top**k
      enddo
      w0=(t_sleeve+gap-weldxx)/dr
      if (w0 .lt. dabs(y(1)-y(11)) .and. ne_sleeve_top .gt. 1) then
         ne_sleeve_top = ne_sleeve_top-1
         goto 600
      end if
      if (ne_sleeve_top .lt. 1) ne_sleeve_top = 1
      do 620 i=1,ne_pipe_right+1
         id_old = 601+20*(7-(id_mini-1))+(i-1)
         id_new = id_old+20*ne_sleeve_top
         x(id_new) = x(id_old)
         y(id_new) = y(601)+t_sleeve
 620  continue

      do 640 i=1,ne_pipe_right+1
         nb = 600+1+20*(7-(id_mini-1))+(i-1)
         ne = nb+20*ne_sleeve_top
         nd = 20
         bf = bf_sleeve_top
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ng_line(nb,ne,nd,co_f,co_l,bf)
 640  continue

      n_repeat(1,1) = 601+20*(7-(id_mini-1))
      n_repeat(1,2) = 1

      do k=1,ne_sleeve_top
         nd_delta = 20*(k-1)+20*(9-id_mini)
         nd_count_delta = (ne_pipe_right+1)*(k-1)
         do i=601+nd_delta,601+ne_pipe_right+nd_delta
            nd_count = nnode+i-600-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = y(i)+0.5d0*D_pipe+weldxx
            coord(2,nd_count) = -x(i)+0.5d0*vl_pipe
         enddo
      enddo
      nnode = nd_count

      n_ele = nele1row*nelerows+n_ele
      nd1 = 601+20*(7-(id_mini-1))
      nd2 = nd1+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_pipe_right
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_sleeve_top
      nddtrows = 20
      nedtrows = nele1row
      call elg4(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&           nelerows,nddtrows,nedtrows,larray2D)
      do 660 j=1,nele1row*nelerows
      do 660 i=2,5
         if (larray2D(j,i) .eq. n_repeat(1,1)) then
            larray2D(j,i)=n_repeat(1,2)
         endif
 660  continue

      id_mat = 2
      do 680 j=1,nele1row*nelerows
         if (mod(j,ne_pipe_right) .eq. 1) then
            id_location = 1
            id_face = 4
         else if (j .gt. (ne_sleeve_top-1)*ne_pipe_right) then
            id_location = 1
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if
         if ((mod(j,ne_pipe_right) .eq. 1) .and. &        (j .gt. (ne_sleeve_top-1)*ne_pipe_right)) then
            do i=1,5
               larray_ele(i,j+n_ele-1) = larray2D(j,i)
            enddo

            larray_ele(6,j+n_ele-1) = id_mat
            larray_ele(7,j+n_ele-1) = id_location
            larray_ele(8,j+n_ele-1) = id_face
         else
            do i=1,5
               larray_ele(i,j+n_ele-1) = larray2D(j,i)
            enddo
            larray_ele(6,j+n_ele-1) = id_mat
            larray_ele(7,j+n_ele-1) = id_location
            larray_ele(8,j+n_ele-1) = id_face
         end if
 680  continue
      nelem = n_ele-1+nele1row*nelerows

      if (meshtest .eq. 0) goto 999

      open(15,file='node.inp')
      open(16,file='elem.inp')
      do i=1,nnode
         write(15,4) inode(i),coord(1,i),coord(2,i)
      enddo
      do i=1,nelem
         write(16,8) (larray_ele(ii,i),ii=1,5)
      enddo
      close(15)
      close(16)
      stop
 4    format(i4,3(:', ',ES12.5))
 8    format(i4,4(:', ',i4))

 999  write(*,*) 'Meshing completed!'

      return
      end


      subroutine mesh02(d_branch,t_runpipe,t_branch,&  rootgap,vland,weldarea,electro_ang,&  vl_runpipe,vl_branch,size_w,nnode,nelem,inode, &     coord,larray_ele,key_point,r_cc)

      implicit real*8 (a-h,o-z)

      dimension :: key_point(*),larray_ele(8,*)
      dimension :: inode(*),coord(2,*)
      dimension :: x(1000),y(1000),n_repeat(10,2)
      dimension :: larray2D(1000,5),rfw(5)

      data (rfw(i),i=1,5)/0.27, 0.41, 0.6, 0.8, 1.0/

      common /coord/x,y

      if_weldsize_reset = 0
      meshtest = 0
      t_branch_cri=0.1d0
      bf_t_runpipe=1.5d0
      bf_runpipe_left=1.5d0
      bf_runpipe_right=1.5d0
      bf_branch_length=1.5d0
      bf_branch_thick=1.2d0
      ne_runpipe=9
      ne_runpipe_left=19
      ne_runpipe_right=19
      ne_branch_length=15
      ne_branch_thick=12
      zero=0.0d0
      pi=acos(-1.0d0)
      to_rad = pi/180.d0
      electro_ang = electro_ang*to_rad
      wprep_ang = 0.5 * electro_ang

      write(*,3) "Estimated cross-sectional area ", &     "of the weld deposit =", weldarea,"mm^2."
 3    format(1x,A,A,ES9.2,1x,A,A)

      angle1=electro_ang
      angle2=wprep_ang
      if (rootgap .lt. 0.5*vland) rootgap=0.5*vland

      vxmin=0.3*vland
      vymin=vxmin*tan(angle2)
      weldmin=rootgap*(vland+vxmin) &     +0.5d0*vxmin*vymin &     +0.5d0*(vymin+rootgap)**2/tan(angle1)

      if (weldarea .le. weldmin) then
         vx=zero
         vy=zero
         rootgap=sqrt(weldarea)
         vland=rootgap
         size_w=zero
         r_cc=0.5*rootgap
         key_point(1)=83
         key_point(2)=43
         key_point(3)=3
      else
         a=0.5*tan(angle2)*(1.+tan(angle2)/tan(angle1))
         b=rootgap*(1.+tan(angle2)/tan(angle1))
         c=rootgap*vland+0.5*rootgap*rootgap/tan(angle1)-weldarea
         vx=0.5*(-b+sqrt(b*b-4.*a*c))/a

         if (vland+vx .gt. t_branch-t_branch_cri) then
            vx=t_branch-vland
            vy=vx*tan(angle2)
            area3=weldarea-rootgap*vland-0.5*vx*(2.*rootgap+vy)
            weldsize=vx+2.*area3/(vy+rootgap)
            electro_ang=atan((vy+rootgap)/(weldsize-vx))
         else
            vy=vx*tan(angle2)
            weldsize=vx+(vy+rootgap)/tan(angle1)
         endif
         size_w=(vy+rootgap)/sin(electro_ang)
         r_cc=0.5*size_w
         key_point(1)=89
         key_point(2)=49
         key_point(3)=9
      endif

      x0=0.5*d_branch-t_branch
      write(*,3) 'Meshing in progress!'

      id_mat = 3
      nd_count=0
      do 10 j=1,3
         do i=1,3
            n=i+(j-1)*40
            x(n)=0.5*(i-1)*vland
            y(n)=0.5*(j-1)*rootgap
            nd_count=nd_count+1
            inode(nd_count)=n
            coord(1,nd_count) = x(n)+x0
            coord(2,nd_count) = y(n)
         enddo
 10   continue

      nelem = 1
      nd1 = 1
      nd2 = nd1+1
      nd3 = nd2+40
      nd4 = nd1+40
      nele1row = 2
      nddt1row = 1
      nedt1row = 1
      nelerows = 2
      nddtrows = 40
      nedtrows = 2
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,&     nedt1row,nelerows,nddtrows,nedtrows,larray2D)
      do 20 j=1,4
         if ((j .eq. 1) .or. (j .eq. 3)) then
            id_location = 1
            id_face = 4
         else
            id_location = 0
            id_face = 0
         endif
         do i=1,5
            larray_ele(i,j) = larray2D(j,i)
         enddo
         larray_ele(6,j) = id_mat
         larray_ele(7,j) = id_location
         larray_ele(8,j) = id_face
 20   continue
      nnode = nd_count
      nelem=5

      if (vx .lt. 1.0d-3) then
         vx=t_branch-vland
         vy=vx*tan(angle2)
         weldsize=t_branch-vland
      endif

      do i=1,5
         x(4+i) = vland+rfw(i)*weldsize
         y(4+i) = y(1)
         x(84+i) = vland+rfw(i)*vx
         y(84+i) = rootgap+rfw(i)*vy
      enddo
      do n=5,9
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo
      do n=85,89
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo
      nnode = nd_count

      if (weldarea .le. weldmin) goto 90
      do i=1,5
         x(4+i) = vland+rfw(i)*weldsize
         y(4+i) = y(1)
         x(84+i) = vland+rfw(i)*vx
         y(84+i) = rootgap+rfw(i)*vy
      enddo

      x(45) = 0.5*(x(5)+x(85))
      y(45) = 0.5*(y(5)+y(85))
      x(25) = 0.5*(x(5)+x(45))
      y(25) = 0.5*(y(5)+y(45))
      x(65) = 0.5*(x(85)+x(45))
      y(65) = 0.5*(y(85)+y(45))

      x(44) = 0.5*(x(43)+x(45))
      y(44) = 0.5*(y(43)+y(45))
      x(24) = 0.25*(x(3)+x(5)+x(43)+x(45))
      y(24) = 0.25*(y(3)+y(5)+y(43)+y(45))
      x(64) = 0.25*(x(83)+x(85)+x(43)+x(45))
      y(64) = 0.25*(y(83)+y(85)+y(43)+y(45))

      nnum=7
      ninc=10
      bias=1.0
      do i=0,2
         node1 = 7+i
         node2 = 87+i
         nstart=node1+ninc
         call nfill(node1,node2,nnum,nstart,ninc,bias)
      enddo

      do k=0,3
         x(16+k*20)=0.25*(x(5+k*20)+x(7+k*20)+x(25+k*20)+x(27+k*20))
         y(16+k*20)=0.25*(y(5+k*20)+y(7+k*20)+y(25+k*20)+y(27+k*20))
      enddo
      x(46) = 0.5*(x(45)+x(47))
      y(46) = 0.5*(y(45)+y(47))
      do j=1,3
         n=24+(j-1)*20
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo

      do j=1,5
         n = 5+(j-1)*20
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo

      do j=1,2
         n = 6+(j-1)*10
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo

      do j=1,3
         n = 36+(j-1)*10
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo

      do j=1,2
         n = 76+(j-1)*10
         nd_count = nd_count+1
         inode(nd_count) = n
         coord(1,nd_count) = x(n)+x0
         coord(2,nd_count) = y(n)
      enddo

      nnode=nd_count
      do k=1,9
         nd_delta = (k-1)*10
         nd_count_delta = 3*(k-1)
         do n=7+nd_delta,9+nd_delta
            nd_count = nnode+n-6-nd_delta+nd_count_delta
            inode(nd_count) = n
            coord(1,nd_count) = x(n)+x0
            coord(2,nd_count) = y(n)
         enddo
      enddo
      nnode = nd_count

      call elgen(5,3,5,25,24,larray2D)
      call elgen(6,3,24,44,43,larray2D)
      call elgen(7,24,25,45,44,larray2D)
      call elgen(8,43,44,64,83,larray2D)
      call elgen(9,44,45,65,64,larray2D)
      call elgen(10,83,64,65,85,larray2D)
      call elgen(11,5,6,16,25,larray2D)
      call elgen(12,6,7,17,16,larray2D)
      call elgen(13,25,16,17,27,larray2D)
      call elgen(14,25,27,37,36,larray2D)
      call elgen(15,25,36,46,45,larray2D)
      call elgen(16,36,37,47,46,larray2D)
      call elgen(17,45,46,56,65,larray2D)
      call elgen(18,46,47,57,56,larray2D)
      call elgen(19,65,56,57,67,larray2D)
      call elgen(20,65,67,77,76,larray2D)
      call elgen(21,65,76,86,85,larray2D)
      call elgen(22,76,77,87,86,larray2D)

      do 40 j=5,22
         do i=1,5
            larray_ele(i,j)=larray2D(j,i)
         enddo
         larray_ele(6,j) = id_mat
         larray_ele(7,j) = 0
         larray_ele(8,j) = 0
 40   continue

      nelem = 23
      nd1 = 7
      nd2 = nd1+1
      nd3 = nd2+10
      nd4 = nd1+10
      nele1row = 2
      nddt1row = 1
      nedt1row = 1
      nelerows = 8
      nddtrows = 10
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,&     nedt1row,nelerows,nddtrows,nedtrows,larray2D)
      do 80 j=1,16
         if (mod(j,2) .eq. 0) then
            id_location = 1
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
 80   continue
      nelem = nelem+nele1row*nelerows

 90   id_mat = 1
 100  dr=1.0
      do k=1,ne_runpipe-1,1
         dr=dr+bf_t_runpipe**k
      enddo
      w0=t_runpipe/dr
      dist_comp = bf_t_runpipe*min(0.5*vland,abs(x(5)-x(6)))
      if ((w0 .lt. dist_comp) .and. (ne_runpipe .gt. 1)) then
         ne_runpipe = ne_runpipe-1
         goto 100
      end if

      nd_grade = 10*ne_runpipe
      key_point(9) = 101+nd_grade
      key_point(10) = 105+nd_grade
      key_point(11) = 108+nd_grade

      do i=1,3
         x(100+i) = x(i)
         y(100+i) = y(i)
      enddo
      do i=1,5
         x(103+i) = x(4+i)
         y(103+i) = y(4+i)
      enddo

      do i=1,8
         x(100+i+nd_grade) = x(100+i)
         y(100+i+nd_grade) = y(100+i)-t_runpipe
      enddo

      nnum = ne_runpipe-1
      ninc = 10
      bias = bf_t_runpipe
      do i=1,8
         node1 = 100+i
         node2 = node1+nd_grade
         nstart = node1+ninc
         call nfill(node1,node2,nnum,nstart,ninc,bias)
      enddo

      do k=1,8
         n_repeat(k,1) = 100+k
         if (k .le. 3) then
            n_repeat(k,2) = k
         else
            n_repeat(k,2) = k+1
         endif
      enddo

      do k=1,ne_runpipe
         nd_delta = 10*(k-1)
         nd_count_delta = 8*(k-1)
         do n=111+nd_delta,118+nd_delta
            nd_count = nnode+n-110-nd_delta+nd_count_delta
            inode(nd_count) = n
            coord(1,nd_count) = x(n)+x0
            coord(2,nd_count) = y(n)
         enddo
      enddo
      nnode = nd_count
      nd1 = 101
      nd2 = nd1+10
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = ne_runpipe
      nddt1row = 10
      nedt1row = 1
      nelerows = 7
      nddtrows = 1
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,&     nedt1row,nelerows,nddtrows,nedtrows,larray2D)
      do j=1,nele1row*nelerows
         do i=2,5
            do k=1,8
               if (larray2D(j,i) .eq. n_repeat(k,1))&              larray2D(j,i)=n_repeat(k,2)
            enddo
         enddo
      enddo

      do j=1,nele1row*nelerows
         if (mod(j,ne_runpipe) .eq. 0) then
            id_location = 2
            id_face = 2
         else
            id_location = 0
            id_face = 0
         endif

         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
      enddo
      nelem=nelem+nele1row*nelerows

 200  dr=1.d0
      do k=1,ne_runpipe_left-1,1
         dr=dr+bf_runpipe_left**k
      enddo
      w0=(0.5*d_branch-t_branch)/dr
      w0c=(1.+0.5*(bf_runpipe_left-1.))*abs(x(2)-x(1))
      if ((w0 .lt. w0c) .and. (ne_runpipe_left .gt. 1)) then
         ne_runpipe_left = ne_runpipe_left-1
         goto 200
      end if

      do i=1,ne_runpipe+1
         x(201+20*(i-1)) = x(101+10*(i-1))
         y(201+20*(i-1)) = y(101+10*(i-1))
      enddo

      do i=1,ne_runpipe+1
         x(201+ne_runpipe_left+20*(i-1))=-0.5*d_branch+t_branch
         y(201+ne_runpipe_left+20*(i-1))=y(101+10*(i-1))
      enddo

      nnum = ne_runpipe_left-1
      ninc = 1
      bias = bf_runpipe_left
      do i=1,ne_runpipe+1
         node1 = 201+20*(i-1)
         node2 = node1+ne_runpipe_left
         nstart = node1+ninc
         call nfill(node1,node2,nnum,nstart,ninc,bias)
      enddo

      n_repeat(1,1) = 201
      n_repeat(1,2) = 1
      do k=2,ne_runpipe+1
         n_repeat(k,1) = 201+20*(k-1)
         n_repeat(k,2) = 101+10*(k-1)
      enddo

      do k=1,ne_runpipe+1
         nd_delta = 20*(k-1)
         nd_count_delta = ne_runpipe_left*(k-1)
         do n=202+nd_delta,202+(ne_runpipe_left-1)+nd_delta
            nd_count = nnode+n-201-nd_delta+nd_count_delta
            inode(nd_count) = n
            coord(1,nd_count) = x(n)+x0
            coord(2,nd_count) = y(n)
         enddo
      enddo
      nnode = nd_count

      nd1 = 201
      nd2 = nd1+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_runpipe_left
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_runpipe
      nddtrows = 20
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&     nelerows,nddtrows,nedtrows,larray2D)
      do j=1,nele1row*nelerows
         do i=2,5
            do k=1,ne_runpipe+1
               if (larray2D(j,i) .eq. n_repeat(k,1)) &              larray2D(j,i)=n_repeat(k,2)
            enddo
         enddo
      enddo

      do j=1,nele1row*nelerows
         if (j .le. ne_runpipe_left) then
            id_location = 1
            id_face = 1
         else if (j .gt. ne_runpipe_left*(ne_runpipe-1)) then
            id_location = 2
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if

         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
      enddo
      nelem = nelem+nele1row*nelerows

 300  dr=1.d0
      do k=1,ne_runpipe_right-1,1
         dr=dr+bf_runpipe_right**k
      enddo
      w0=(vl_runpipe+t_branch-weldsize-vland)/dr
      w0c=(1.d0+0.5d0*(bf_runpipe_right-1.d0))*abs(x(8)-x(9))
      if ((w0 .lt. w0c) .and. (ne_runpipe_right .gt. 1)) then
         ne_runpipe_right = ne_runpipe_right-1
         goto 300
      end if

      do i=1,ne_runpipe+1
         x(401+20*(i-1)) = x(108+10*(i-1))
         y(401+20*(i-1)) = y(108+10*(i-1))
      enddo

      do i=1,ne_runpipe+1
         x(401+ne_runpipe_right+20*(i-1)) = vl_runpipe+t_branch
         y(401+ne_runpipe_right+20*(i-1)) = y(108+10*(i-1))
      enddo

      nnum = ne_runpipe_right-1
      ninc = 1
      bias = bf_runpipe_right
      do i=1,ne_runpipe+1
         node1 = 401+20*(i-1)
         node2 = node1+ne_runpipe_right
         nstart = node1+ninc
         call nfill(node1,node2,nnum,nstart,ninc,bias)
      enddo

      n_repeat(1,1) = 401
      n_repeat(1,2) = 9
      do k=2,ne_runpipe+1
         n_repeat(k,1) = 401+20*(k-1)
         n_repeat(k,2) = 108+10*(k-1)
      enddo

      do k=1,ne_runpipe+1
         nd_delta = 20*(k-1)
         nd_count_delta = ne_runpipe_right*(k-1)
         do n=402+nd_delta,402+(ne_runpipe_right-1)+nd_delta
            nd_count = nnode+n-401-nd_delta+nd_count_delta
            inode(nd_count) = n
            coord(1,nd_count) = x(n)+x0
            coord(2,nd_count) = y(n)
         enddo
      enddo
      nnode = nd_count

      nd2 = 401
      nd1 = nd2+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_runpipe_right
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_runpipe
      nddtrows = 20
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row, &     nelerows,nddtrows,nedtrows,larray2D)
      do j=1,nele1row*nelerows
         do i=2,5
            do k=1,ne_runpipe+1
               if (larray2D(j,i) .eq. n_repeat(k,1)) &              larray2D(j,i)=n_repeat(k,2)
            enddo
         enddo
      enddo

      do j=1,nele1row*nelerows
         if (j .le. ne_runpipe_right) then
            id_location = 1
            id_face = 1
         else if (j .gt. ne_runpipe_right*(ne_runpipe-1)) then
            id_location = 2
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if

         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
      enddo
      nelem = nelem+nele1row*nelerows

      key_point(4) = 761
      key_point(5) = 781
      key_point(6) = 402
      key_point(7) = 403
      key_point(8) = 404

      do 320 i=1,3
         n_repeat(i,1) = 601+20*(i-1)
         n_repeat(i,2) = 81+(i-1)
         x(n_repeat(i,1)) = x(n_repeat(i,2))
         y(n_repeat(i,1)) = y(n_repeat(i,2))
 320  continue

      do 340 i=4,8
         n_repeat(i,1) = 601+20*(i-1)
         n_repeat(i,2) = 82+(i-1)
         x(n_repeat(i,1)) = x(n_repeat(i,2))
         y(n_repeat(i,1)) = y(n_repeat(i,2))
 340  continue

      fl_branch_local = 1.d0+vx/(t_branch-vland)
      vl_branch_local = rootgap+fl_branch_local* &     dtan(wprep_ang)*(t_branch-vland)
      do 360 i=1,8
         x(601+4+20*(i-1)) = x(601+20*(i-1))
         y(601+4+20*(i-1)) = vl_branch_local
 360  continue

      nnum=3
      ninc = 1
      bias = bf_branch_length
      do 380 i=1,8
         node1 = 601+20*(i-1)
         node2 = node1+4
         nstart = node1+1
         call nfill(node1,node2,nnum,nstart,ninc,bias)
 380  continue

 400  dr=1.d0
      do k=1,ne_branch_length-1,1
         dr=dr+bf_branch_length**k
      enddo
      w0=(vl_branch-vl_branch_local)/dr
      w0c=(1.d0+0.5d0*(bf_branch_length-1.d0))*dabs(y(605)-y(604))
      if ((w0 .lt. w0c) .and. (ne_branch_length .gt. 1)) then
         ne_branch_length = ne_branch_length-1
         goto 400
      end if

      do i=1,8
         ii=601+4+ne_branch_length+20*(i-1)
         x(ii) = x(605+20*(i-1))
         y(ii) = vl_branch
      enddo

      nnum = ne_branch_length-1
      ninc = 1
      bias = bf_branch_length
      do 420 i=1,8
         node1 = 605+20*(i-1)
         node2 = node1+ne_branch_length
         nstart = node1+1
         call nfill(node1,node2,nnum,nstart,ninc,bias)
 420  continue

      do 440 k=1,8
         nd_delta = 20*(k-1)
         nd_count_delta = (ne_branch_length+4)*(k-1)
         do i=602+nd_delta,602+(ne_branch_length+3)+nd_delta
            nd_count = nnode+i-601-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = x(i)+0.5d0*d_branch-t_branch
            coord(2,nd_count) = y(i)
         enddo
 440  continue
      nnode = nd_count

      nd1 = 601
      nd2 = nd1+20
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = 4
      nddt1row = 1
      nedt1row = 1
      nelerows = 7
      nddtrows = 20
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row, &     nelerows,nddtrows,nedtrows,larray2D)
      do j=1,nele1row*nelerows
         do i=2,5
            do k=1,8
               if (larray2D(j,i) .eq. n_repeat(k,1))&              larray2D(j,i)=n_repeat(k,2)
            enddo
         enddo
      enddo

      do 460 j=1,nele1row*nelerows
         if (j .le. 4) then
            id_location = 1
            id_face = 4
         else if (j .gt. 24 .and. if_weldsize_reset .eq. 1) then
            id_location = 1
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if

         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
 460  continue
      nelem = nelem+nele1row*nelerows

      nd1 = 605
      nd2 = nd1+20
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = ne_branch_length
      nddt1row = 1
      nedt1row = 1
      nelerows = 7
      nddtrows = 20
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row, &     nelerows,nddtrows,nedtrows,larray2D)

      do j=1,nele1row*nelerows
         if (j .le. ne_branch_length) then
            id_location = 1
            id_face = 4
         else if (j .gt. ne_branch_length*6 .and.  &           if_weldsize_reset .eq. 1) then
            id_location = 1
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if

         do i=1,5
            larray_ele(i,j+nelem-1) = larray2D(j,i)
         enddo
         larray_ele(6,j+nelem-1) = id_mat
         larray_ele(7,j+nelem-1) = id_location
         larray_ele(8,j+nelem-1) = id_face
      enddo
      nnode=nd_count
      nelem = nele1row*nelerows+nelem
      nel0=nelem
      w0=(t_branch-vland-vx)
      if (w0 .lt. 1.0d-3) then
         nelem=nelem-1
         print *, 'Meshing completed!'
         return
      endif
 500  dr=1.d0
      do k=1,ne_branch_thick-1,1
         dr=dr+bf_branch_thick**k
      enddo
      w0=(t_branch-vland-vx)/dr
      w0c=(1.d0+0.5d0*(bf_branch_thick-1.d0))*dabs(x(741)-x(721))
      if ((w0 .lt. w0c) .and. (ne_branch_thick .gt. 1)) then
         ne_branch_thick = ne_branch_thick-1
         goto 500
      endif

      x(741+20*ne_branch_thick) = t_branch
      y(741+20*ne_branch_thick) = rootgap+ &     (t_branch-vland)*dtan(wprep_ang)

      n_repeat(1,1) = 741
      n_repeat(1,2) = 89
      node1 = 741
      node2 = node1+20*ne_branch_thick
      nnum=ne_branch_thick-1
      ninc = 20
      nstart = node1+ninc
      bias = bf_branch_thick
      call nfill(node1,node2,nnum,nstart,ninc,bias)
      x(745+20*ne_branch_thick) = t_branch
      y(745+20*ne_branch_thick) = y(745)
      node1 = 745
      node2 = node1+20*ne_branch_thick
      ninc=20
      nstart = node1+ninc
      call nfill(node1,node2,nnum,nstart,ninc,bias)
      x(745+20*ne_branch_thick+ne_branch_length) = t_branch
      y(745+20*ne_branch_thick+ne_branch_length) = vl_branch

      node1 = 745+ne_branch_length
      node2 = node1+20*ne_branch_thick
      nstart = node1+ninc
      call nfill(node1,node2,nnum,nstart,ninc,bias)

      bias = bf_branch_length
      do 520 i=1,ne_branch_thick
         node1 = 741+20*i
         node2 = 745+20*i
         ninc = 1
         nstart=node1+1
         bias = bf_branch_length
         nnum=3
         call nfill(node1,node2,nnum,nstart,ninc,bias)
 520  continue

      do 540 i=1,ne_branch_thick
         node1 = 745+20*i
         node2 = node1+ne_branch_length
         ninc = 1
         nstart=node1+1
         nnum=ne_branch_length-1
         bias = bf_branch_length
         call nfill(node1,node2,nnum,nstart,ninc,bias)
 540  continue

      do 560 k=1,ne_branch_thick
         nd_delta = 20*(k-1)
         nd_count_delta = (ne_branch_length+5)*(k-1)
         do i=761+nd_delta,761+(ne_branch_length+4)+nd_delta
            nd_count = nnode+i-760-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = x(i)+0.5d0*d_branch-t_branch
            coord(2,nd_count) = y(i)
         enddo
 560  continue
      nnode = nd_count

      nd1 = 741
      nd2 = nd1+20
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = 4+ne_branch_length
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_branch_thick
      nddtrows = 20
      nedtrows = nele1row
      call elg4(nelem,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row, &     nelerows,nddtrows,nedtrows,larray2D)

      do j=1,nele1row*nelerows
         do i=2,5
            if (larray2D(j,i) .eq. n_repeat(1,1))&           larray2D(j,i)=n_repeat(1,2)
         enddo
      enddo

      do j=1,nele1row*nelerows
         if (mod(j,(ne_branch_length+4)) .eq. 1) then
            id_location = 1
            id_face = 1
         else if (j .gt. (ne_branch_length+4)*(ne_branch_thick-1)) then
            id_location = 1
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if

         if ((mod(j,(ne_branch_length+4)) .eq. 1) .and. &        (j .gt. (ne_branch_length+4)*(ne_branch_thick-1))) then
            do i=1,5
               larray_ele(i,j+nelem-1) = larray2D(j,i)
            enddo
            larray_ele(6,j+nelem-1) = id_mat
            larray_ele(7,j+nelem-1) = id_location
            larray_ele(8,j+nelem-1) = id_face
         else
            do i=1,5
               larray_ele(i,j+nelem-1) = larray2D(j,i)
            enddo
            larray_ele(6,j+nelem-1) = id_mat
            larray_ele(7,j+nelem-1) = id_location
            larray_ele(8,j+nelem-1) = id_face
         endif
      enddo
      nelem = nelem-1+nele1row*nelerows

      if (meshtest .eq. 0) goto 999
      open(15,file='node.inp')
      open(16,file='elem.inp')
      do i=1,nnode
         write(15,4) inode(i),coord(1,i),coord(2,i)
      enddo
      do i=1,nelem
         write(16,8) (larray_ele(ii,i),ii=1,5)
      enddo
      close(15)
      close(16)
      stop
 4    format(i4,3(:', ',ES12.5))
 8    format(i4,4(:', ',i4))

 999  write(*,*) 'Meshing completed!'
      return
      end

      subroutine mesh03(d_pipe,t_pipe,weldarea,vl_pipe,size_w, &     nnode,nelem,inode,coord,larray_ele,key_point,r_cc)
      implicit double precision (a-h,o-z)

      dimension :: key_point(*),larray_ele(8,*)
      dimension :: inode(*),coord(2,*)
      dimension :: node_symm(30),x(1000),y(1000),n_repeat(10,2)
      dimension :: co_f(2),co_l(2),larray2d(1000,5)
      common /coord/x,y
      pi=acos(-1.0)
      to_rad = pi/180.d0
      r_pipe=0.5*d_pipe
      bf_t_pipe=1.3
      bf_pipe_left=1.5
      ne_pipe=9
      ne_pipe_left=19
      ne_pipe_right=10
      ang_set=30.*to_rad
      ang_midnode=30.*to_rad
      f_node41=0.4
      f_node49=0.6
      whratio=3.2
      alpha=atan(0.5*whratio)
      ang_cir=pi-2.*alpha
      theta=2.*ang_cir
      ra_cir=sqrt(2.*weldarea/(theta-sin(theta)))
      weldwidth=2.*ra_cir*sin(ang_cir)
      weldheight=weldwidth/whratio
      weldwidth_h = 0.5*weldwidth
      ang_sector = 0.125*ang_cir
      write(*,3) 'Estimated cross-sectional area of weld deposit =', &     weldarea,'mm^2'
      write(*,3) 'Estimated weld bead width =',weldwidth,'mm'
      write(*,3) 'Estimated weld bead height =',weldheight,'mm'
      write(*,3) 'Meshing in progress!'
 3    format(1x,A,F6.2,1x,A)

      r_cc=sqrt(weldheight**2+weldwidth_h**2)
      size_w=ra_cir*theta
      key_point(1)=1
      key_point(2)=81

      x(1) = r_pipe+weldheight
      y(1) = 0.0d0
      x(9) = r_pipe
      y(9) = weldwidth_h
      x(81) = r_pipe
      y(81) = 0.0d0
      do 20 i=2,8,1
         x(i) = r_pipe + weldheight &        - ra_cir*(1.d0-dcos(ang_sector*dfloat(i-1)))
         y(i) = ra_cir*dsin(ang_sector*dfloat(i-1))
 20   continue

      x(41) = r_pipe + weldheight - f_node41*weldheight
      y(41) = 0.0d0
      x(49) = r_pipe
      y(49) = f_node49*weldwidth_h
      x(45) = 0.5d0*(x(49)+x(5))
      y(45) = y(49)-dabs(x(49)-x(45))*dtan(ang_midnode)

      nb = 45
      ne = 49
      nd = 1
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ngline(nb,ne,nd,co_f,co_l,bf)

      x(55) = x(46)
      y(55) = y(46)
      x(65) = x(47)
      y(65) = y(47)
      x(75) = x(48)
      y(75) = y(48)
      x(85) = x(49)
      y(85) = y(49)
      nb = 41
      ne = 45
      nd = 1
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ngline(nb,ne,nd,co_f,co_l,bf)

      do 40 i=1,5
         nb = 1+(i-1)
         ne = 41+(i-1)
         nd = 10
         bf = 1.1
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ngline(nb,ne,nd,co_f,co_l,bf)
 40   continue

      nb = 81
      ne = 85
      nd = 1
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ngline(nb,ne,nd,co_f,co_l,bf)

      do 60 i=1,4
         nb = 41+(i-1)
         ne = 81+(i-1)
         nd = 10
         bf = 1.0
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ngline(nb,ne,nd,co_f,co_l,bf)
 60   continue

      nb = 45
      ne = 49
      nd = 1
      bf = 1.0
      co_f(1) = x(nb)
      co_f(2) = y(nb)
      co_l(1) = x(ne)
      co_l(2) = y(ne)
      call ngline(nb,ne,nd,co_f,co_l,bf)

      do 80 i=1,4
         nb = 6+(i-1)
         ne = 46+(i-1)
         nd = 10
         bf = 1.1
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ngline(nb,ne,nd,co_f,co_l,bf)
 80   continue

      do 100 i=1,9,1
         node_symm(i) = 1+10*(i-1)
 100  continue

      do 120 k=1,4
         n_repeat(k,1) = 55+10*(k-1)
         n_repeat(k,2) = 46+(k-1)
 120  continue

      n_ele = 1
      nd1 = 1
      nd2 = nd1+1
      nd3 = nd2+10
      nd4 = nd1+10
      nele1row = 8
      nddt1row = 1
      nedt1row = 1
      nelerows = 4
      nddtrows = 10
      nedtrows = nele1row
      call elg(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,nedt1row,&     nelerows,nddtrows,nedtrows,larray2d)
      id_mat = 3
      do 140 j=1,32
         if (j .le. 8) then
            id_location = 1
            id_face = 1
         else
            id_location = 0
            id_face = 0
         end if
         do m=1,5
            larray_ele(m,j) = larray2d(j,m)
         enddo
         larray_ele(6,j) = id_mat
         larray_ele(7,j) = id_location
         larray_ele(8,j) = id_face
 140  continue

      n_ele = nele1row*nelerows+n_ele
      nd1 = 41
      nd2 = nd1+1
      nd3 = nd2+10
      nd4 = nd1+10
      nele1row = 4
      nddt1row = 1
      nedt1row = 1
      nelerows = 4
      nddtrows = 10
      nedtrows = nele1row
      call elg(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row, &     nedt1row,nelerows,nddtrows,nedtrows,larray2d)
      do 160 j=1,16
      do 160 i=2,5
      do 160 k=1,4
         if (larray2d(j,i) .eq. n_repeat(k,1))  &        larray2d(j,i)=n_repeat(k,2)
 160  continue

      id_location = 0
      id_face = 0
      do 180 j=1,16
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2d(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 180  continue

 200  dr=1.d0
      do k=1,ne_pipe-1,1
         dr=dr+bf_t_pipe**k
      enddo
      w0=t_pipe/dr
      if (w0 .lt. abs(x(48)-x(49))) then
         ne_pipe = ne_pipe-1
         goto 200
      endif

      nd_grade = 10*ne_pipe
      key_point(3) = 101+nd_grade
      key_point(4) = 109+nd_grade
      do i=1,ne_pipe
         node_symm(i+9) = 109+10*i
      enddo

      do 220 i=1,5
         x(100+i) = x(9+(i-1)*10)
         y(100+i) = y(9+(i-1)*10)
 220  continue

      do 240 i=1,4
         x(105+i) = x(85-i)
         y(105+i) = y(85-i)
 240  continue

      do 260 i=1,9
         x(100+i+ne_pipe*10) = x(100+i)-t_pipe
         y(100+i+ne_pipe*10) = y(100+i)
 260  continue

      do 280 i=1,9
         nb = 100+i
         ne = nb+ne_pipe*10
         nd = 10
         bf = bf_t_pipe
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ngline(nb,ne,nd,co_f,co_l,bf)
 280  continue

      do k=1,5
         n_repeat(k,1) = 101+(k-1)
         n_repeat(k,2) = 9+10*(k-1)
      enddo
      do k=6,9
         n_repeat(k,1) = 101+(k-1)
         n_repeat(k,2) = 84-(k-6)
      enddo

      n_ele = nele1row*nelerows+n_ele
      nd1 = 101
      nd2 = nd1+10
      nd3 = nd2+1
      nd4 = nd1+1
      nele1row = ne_pipe
      nddt1row = 10
      nedt1row = 1
      nelerows = 8
      nddtrows = 1
      nedtrows = nele1row
      call elg(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,
     1     nedt1row,nelerows,nddtrows,nedtrows,larray2d)
      do 300 j=1,nele1row*nelerows
      do 300 i=2,5
      do 300 k=1,9
      if (larray2d(j,i) .eq. n_repeat(k,1)) &        larray2d(j,i)=n_repeat(k,2)
 300  continue

      id_mat = 1
      do 320 j=1,nele1row*nelerows
         if (mod(j,ne_pipe) .eq. 0) then
            id_location = 2
            id_face = 2
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2d(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 320  continue

 340  dr=1.d0
      do k=1,ne_pipe_left-1,1
         dr=dr+bf_pipe_left**k
      enddo
      w0=(0.5d0*vl_pipe-weldwidth_h)/dr
      if (w0 .lt. abs(y(101)-y(102))) then
         ne_pipe_left = ne_pipe_left-1
         goto 340
      endif

      do i=1,ne_pipe+1
         x(201+20*(i-1)) = x(101+10*(i-1))
         y(201+20*(i-1)) = y(101+10*(i-1))
      enddo

      do i=1,ne_pipe+1
         x(201+ne_pipe_left+20*(i-1)) = x(101+10*(i-1))
         y(201+ne_pipe_left+20*(i-1)) = 0.5d0*vl_pipe
      enddo

      do 360 i=1,ne_pipe+1
         nb = 201+20*(i-1)
         ne = nb+ne_pipe_left
         nd = 1
         bf = bf_pipe_left
         co_f(1) = x(nb)
         co_f(2) = y(nb)
         co_l(1) = x(ne)
         co_l(2) = y(ne)
         call ngline(nb,ne,nd,co_f,co_l,bf)
 360  continue

      n_repeat(1,1) = 201
      n_repeat(1,2) = 9
      do k=2,ne_pipe+1
         n_repeat(k,1) = 201+20*(k-1)
         n_repeat(k,2) = 101+10*(k-1)
      enddo

      n_ele = nele1row*nelerows+n_ele
      nd1 = 201
      nd2 = nd1+1
      nd3 = nd2+20
      nd4 = nd1+20
      nele1row = ne_pipe_left
      nddt1row = 1
      nedt1row = 1
      nelerows = ne_pipe
      nddtrows = 20
      nedtrows = nele1row
      call elg(n_ele,nd1,nd2,nd3,nd4,nele1row,nddt1row,&     nedt1row,nelerows,nddtrows,nedtrows ,larray2d)
      do 380 j=1,nele1row*nelerows
      do 380 i=2,5
      do 380 k=1,ne_pipe+1
      if (larray2d(j,i) .eq. n_repeat(k,1)) &        larray2d(j,i)=n_repeat(k,2)
 380  continue

      id_mat = 1
      do 400 j=1,nele1row*nelerows
         if (j .le. ne_pipe_left) then
            id_location = 1
            id_face = 1
         elseif (j .gt. ne_pipe_left*(ne_pipe-1)) then
            id_location = 2
            id_face = 3
         else
            id_location = 0
            id_face = 0
         end if
         do i=1,5
            larray_ele(i,j+n_ele-1) = larray2d(j,i)
         enddo
         larray_ele(6,j+n_ele-1) = id_mat
         larray_ele(7,j+n_ele-1) = id_location
         larray_ele(8,j+n_ele-1) = id_face
 400  continue

      nelem = n_ele-1+nele1row*nelerows

      do k=1,5
         nd_delta = 10*(k-1)
         nd_count_delta = 9*(k-1)
         do i=1+nd_delta,9+nd_delta
            nd_count = i-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = x(i)
            coord(2,nd_count) = y(i)
         enddo
      enddo
      nnode = nd_count

      do k=1,4
         nd_delta = 10*(k-1)
         nd_count_delta = 4*(k-1)
         do i=51+nd_delta,54+nd_delta
            nd_count = nnode+i-50-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = x(i)
            coord(2,nd_count) = y(i)
         enddo
      enddo
      nnode = nd_count

      do k=1,ne_pipe
         nd_delta = 10*(k-1)
         nd_count_delta = 9*(k-1)
         do i=111+nd_delta,119+nd_delta
            nd_count = nnode+i-110-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) = x(i)
            coord(2,nd_count) = y(i)
         enddo
      enddo
      nnode = nd_count

      do k=1,ne_pipe+1
         nd_delta = 20*(k-1)
         nd_count_delta = ne_pipe_left*(k-1)
         do i=202+nd_delta,202+(ne_pipe_left-1)+nd_delta
            nd_count = nnode+i-201-nd_delta+nd_count_delta
            inode(nd_count) = i
            coord(1,nd_count) =  x(i)
            coord(2,nd_count) =  y(i)
         enddo
      enddo
      nnode = nd_count

      write(*,3) "Meshing completed!"
      return
      end

      subroutine ngline(node_f,node_l,node_d,cord_f,cord_l,b_factor)
      implicit real*8 (a-h,o-z)
      parameter (one=1.0d0, none=1, ntwo=2)
      dimension cord_f(2),cord_l(2),x(1000),y(1000)

      common /coord/x,y

      nnode = (node_l-node_f)/node_d
      dr=one
      do 110 k=none,nnode-none,none
110   dr=dr+b_factor**k
      x0 = (cord_l(1)-cord_f(1))/dr
      y0 = (cord_l(2)-cord_f(2))/dr

      x(node_f+node_d) = cord_f(1)+x0
      y(node_f+node_d) = cord_f(2)+y0
      do 210 k=ntwo,nnode-none,none
      nd = node_f+node_d*k
      x(nd) = x(nd-node_d)+x0*b_factor**(k-none)
      y(nd) = y(nd-node_d)+y0*b_factor**(k-none)
210   enddo


      return
      end

      subroutine elg(n0,n1,n2,n3,n4, &     ne1rw,nd1rw,ned1rw,nrow,ndrw,nedrw,iarray)
      implicit real*8 (a-h,o-z)
      dimension iarray (1000,5)

      iarray(1,1) = n0
      iarray(1,2) = n1
      iarray(1,3) = n2
      iarray(1,4) = n3
      iarray(1,5) = n4
      do 200 i=2,ne1rw
      iarray(i,1) = iarray((i-1),1)+ned1rw
      do 200 l=2,5
200   iarray(i,l) = iarray((i-1),l)+nd1rw

      do 400 j=2,nrow
      do 400 i=1,ne1rw
      iarray(ne1rw*(j-1)+i,1) = iarray(ne1rw*(j-2)+i,1)+nedrw
      do 400 l=2,5,1
400   iarray(ne1rw*(j-1)+i,l) = iarray(ne1rw*(j-2)+i,l)+ndrw

      return
      end

      subroutine nfill (node1,node2,nnum,nstart,ninc,bias)
      implicit real*8 (a-h,o-z)
      parameter (ntest=0)
      dimension x(1000),y(1000)

      common /coord/x,y

      sum=1.0
      do k=1,nnum,1
         sum=sum+bias**k
      enddo
      dx = (x(node2)-x(node1))/sum
      dy = (y(node2)-y(node1))/sum

      x(nstart) = x(node1) + dx
      y(nstart) = y(node1) + dy
      if (node1 .eq. ntest) then
         print *, node1, x(node1),y(node1)
         print *, nstart, x(nstart),y(nstart)
      endif
      do k=2,nnum,1
         nd = nstart+ninc*(k-1)
         x(nd) = x(nd-ninc)+dx*bias**(k-1)
         y(nd) = y(nd-ninc)+dy*bias**(k-1)
         if (node1 .eq. ntest) print *, nd, x(nd), y(nd)
      enddo
      if (node1 .eq. ntest) print *, node2, x(node2),y(node2)

      return
      end

      subroutine elgen(nel,nd1,nd2,nd3,nd4,iarray)
      implicit real*8 (a-h,o-z)
      dimension iarray (1000,5)

      iarray(nel,1) = nel
      iarray(nel,2) = nd1
      iarray(nel,3) = nd2
      iarray(nel,4) = nd3
      iarray(nel,5) = nd4

      return
      end

      function tff(temp,id)
      implicit real*8 (a-h,o-z)
      include 'common.f'

      if (temp .lt. ts(id)) then
         tff=0.0d0
      elseif (temp .lt. tl(id)) then
         tff=(temp-ts(id))/(tl(id)-ts(id))
      else
         tff=1.0d0
      endif

      end

      function tk(temp,id)
      implicit real*8 (a-h,o-z)

      if (id .eq. 1)  then
         tk = tk_steel1(temp)
      elseif (id .eq. 2) then
         tk = tk_steel2(temp)
      elseif (id .eq. 3) then
     tk = tk_steel3(temp)
      endif

      end
      function tk_steel1(temp)
      implicit real *8 (a-h,o-z)
      dimension t(33),value(33)

      data (value(i),t(i),i=1,33) /
     &     5.119d-2,   0.0d0,
     &     5.146d-2,  50.0d0,
     &     5.104d-2, 100.0d0,
     &     4.979d-2, 150.0d0,
     &     4.853d-2, 200.0d0,
     &     4.977d-2, 250.0d0,
     &     4.435d-2, 300.0d0,
     &     4.348d-2, 350.0d0,
     &     4.268d-2, 400.0d0,
     &     4.100d-2, 450.0d0,
     &     3.933d-2, 500.0d0,
     &     3.766d-2, 550.0d0,
     &     3.556d-2, 600.0d0,
     &     3.389d-2, 650.0d0,
     &     3.180d-2, 700.0d0,
     &     2.866d-2, 750.0d0,
     &     2.552d-2, 800.0d0,
     &     2.552d-2, 850.0d0,
     &     2.636d-2, 900.0d0,
     &     2.678d-2, 950.0d0,
     &     2.720d-2,1000.0d0,
     &     2.803d-2,1050.0d0,
     &     2.845d-2,1100.0d0,
     &     2.929d-2,1150.0d0,
     &     2.971d-2,1200.0d0,
     &     3.051d-2,1250.0d0,
     &     3.111d-2,1300.0d0,
     &     3.171d-2,1350.0d0,
     &     3.231d-2,1400.0d0,
     &     3.291d-2,1450.0d0,
     &     3.500d-2,1500.0d0,
     &     6.500d-2,1550.0d0,
     &     1.250d-1,1600.0d0/

      t_c=temp-273.15d0
      if (t_c .le. t(1)) then
         tk_steel1=value(1)
      elseif (t_c .ge. t(33)) then
         tk_steel1=value(33)
      else
         i=int(t_c/5.0d1)+1
         tk_steel1=value(i)+(value(i+1)-value(i))&        *(t_c-t(i))/(t(i+1)-t(i))
      endif

      return
      end

      function tk_steel2(temp)
      implicit real *8 (a-h,o-z)

      integer, parameter :: np= 33

      double precision, dimension (np) :: t,value

      data (value(k),t(k),k=1,np) /
     1      1.590E-02,     0.,
     2      1.590E-02,    50.,
     3      1.631E-02,   100.,
     4      1.673E-02,   150.,
     5      1.715E-02,   200.,
     6      1.757E-02,   250.,
     7      1.841E-02,   300.,
     8      1.924E-02,   350.,
     9      2.008E-02,   400.,
     &      2.092E-02,   450.,
     1      2.175E-02,   500.,
     2      2.301E-02,   550.,
     3      2.384E-02,   600.,
     4      2.468E-02,   650.,
     5      2.552E-02,   700.,
     6      2.635E-02,   750.,
     7      2.677E-02,   800.,
     8      2.635E-02,   850.,
     9      2.677E-02,   900.,
     &      2.761E-02,   950.,
     1      2.803E-02,  1000.,
     2      2.844E-02,  1050.,
     3      2.886E-02,  1100.,
     4      2.928E-02,  1150.,
     5      2.970E-02,  1200.,
     6      3.051E-02,  1250.,
     7      3.111E-02,  1300.,
     8      3.171E-02,  1350.,
     9      3.231E-02,  1400.,
     &      3.291E-02,  1450.,
     1      3.500E-02,  1500.,
     2      6.500E-02,  1550.,
     3      1.250E-01,  1600./

      t_c = temp-273.15d0

      if (t_c .le. t(1))  then

         tk_steel2 = value(1)

      elseif (t_c .ge. t(np)) then

         tk_steel2 = value(6)

      else

         do i=1,np-1
            if ((t_c .gt. t(i)) .and. (t_c .le. t(i+1))) exit
         enddo
         tk_steel2 = value(i) + (value(i+1)-value(i)) &        *(t_c-t(i))/(t(i+1)-t(i))
      endif

      return
      end

      function tk_steel3(temp)
      implicit real *8 (a-h,o-z)

      dimension t(30),value(30)

      data (t(k),value(k),k=1,30) /
     1     100.0d0, 51.1d-3,
     2     150.0d0, 49.8d-3,
     3     200.0d0, 48.6d-3,
     4     250.0d0, 46.5d-3,
     5     300.0d0, 44.4d-3,
     6     350.0d0, 43.5d-3,
     7     400.0d0, 42.7d-3,
     8     450.0d0, 41.0d-3,
     9     500.0d0, 39.4d-3,
     &     550.0d0, 37.7d-3,
     1     600.0d0, 35.6d-3,
     2     650.0d0, 33.9d-3,
     3     700.0d0, 31.8d-3,
     4     750.0d0, 28.5d-3,
     5     800.0d0, 26.0d-3,
     6     850.0d0, 26.0d-3,
     7     900.0d0, 26.4d-3,
     8     950.0d0, 26.8d-3,
     9     100.0d1, 27.2d-3,
     &     105.0d1, 28.1d-3,
     1     110.0d1, 28.5d-3,
     2     115.0d1, 29.3d-3,
     3     120.0d1, 29.7d-3,
     4     125.0d1, 34.0d-3,
     5     130.0d1, 43.0d-3,
     6     135.0d1, 52.0d-3,
     7     140.0d1, 61.0d-3,
     8     145.0d1, 70.0d-3,
     9     150.0d1, 80.0d-3,
     &     300.0d1, 12.5d-2/

      t_c=temp-273.15d0
      if (t_c .le. t(1))  then
         tk_steel3=value(1)
      elseif (t_c .ge. t(30)) then
         tk_steel3=value(30)
      else
         do i=1,29
            if ((t_c .gt. t(i)) .and. (t_c .le. t(i+1))) exit
         enddo
         tk_steel3=value(i)+(value(i+1)-value(i)) &        *(t_c-t(i))/(t(i+1)-t(i))
      endif
      return
      end
      subroutine yurioka (t85,hv)
      implicit real*8 (a-h,o-z)

      include 'common.f'

      fn=50d0*(0.02d0-w_n)
      if (w_b .le. 1.0) then
         deltah=0.0d0
      elseif (w_b .le. 2.5) then
         deltah=0.03d0*fn
      elseif (w_b .le. 4.0) then
         deltah=0.06d0*fn
      else
         deltah=0.09d0*fn
      endif

      if (w_c .le. 0.3) then
         cp=w_c
      else
         cp=w_c/6.+0.25d0
      endif

      ce_i=cp+w_si/24.+w_mn/6.+w_cu/15.+w_ni/12.&     +w_cr*(1.-0.16*sqrt(w_cr))/8.+w_mo/4.+deltah
      tau_mart=exp(10.60*ce_i-4.80)
      h_mart=884.*w_c*(1.-0.3*w_c*w_c)+294.

      ce_ii=w_c+w_si/24.+w_mn/5.+w_cu/10.+w_ni/18.+w_cr/5. &     +w_mo/2.5+w_v/5.+w_nb/3.
      if (ce_ii .le. 0.75) then
         h_bain=197.0*ce_ii+117.
      else
         h_bain=145.0+130.d0*tanh(2.65*ce_ii-0.69)
      endif

      ce_iii=cp+w_mn/3.6+w_cu/20.+w_ni/9.+w_cr/5.+w_mo/4.
      tau_bain=exp(6.20*ce_iii+0.74)

      if (t85 .le. tau_mart) then
         hv=h_mart
      elseif (t85 .gt. tau_bain) then
         hv=h_bain
      else
         x=4.*log(t85/tau_mart)/log(tau_bain/tau_mart)-2.0d0
         hv=0.5d0*(h_mart+h_bain) &        -0.5d0*(h_mart-h_bain)*atan(x)/atan(2.0d0)
      endif

      return
      end
      subroutine yurioka2(t85,hv2)
      implicit real*8 (a-h,o-z)

      include 'common.f'

      fn=50d0*(0.02d0-w_n2)
      if (w_b2 .le. 1.0) then
         deltah=0.0d0
      elseif (w_b2 .le. 2.5) then
         deltah=0.03d0*fn
      elseif (w_b2 .le. 4.0) then
         deltah=0.06d0*fn
      else
         deltah=0.09d0*fn
      endif

      if (w_c2 .le. 0.3) then
         cp=w_c2
      else
         cp=w_c2/6.+0.25d0
      endif

      ce_i=cp+w_si2/24.+w_mn2/6.+w_cu2/15.+w_ni2/12. &     +w_cr2*(1.-0.16*sqrt(w_cr2))/8.+w_mo2/4.+deltah
      tau_mart=exp(10.60*ce_i-4.80)
      h_mart=884.*w_c2*(1.-0.3*w_c2*w_c2)+294.

      ce_ii=w_c2+w_si2/24.+w_mn2/5.+w_cu2/10.+w_ni2/18.+w_cr2/5.&     +w_mo2/2.5+w_v2/5.+w_nb2/3.
      if (ce_ii .le. 0.75) then
         h_bain=197.0*ce_ii+117.
      else
         h_bain=145.0+130.d0*tanh(2.65*ce_ii-0.69)
      endif

      ce_iii=cp+w_mn2/3.6+w_cu2/20.+w_ni2/9.+w_cr2/5.+w_mo2/4.
      tau_bain=exp(6.20*ce_iii+0.74)

      if (t85 .le. tau_mart) then
         hv2=h_mart
      elseif (t85 .gt. tau_bain) then
         hv2=h_bain
      else
         x=4.*log(t85/tau_mart)/log(tau_bain/tau_mart)-2.0d0
         hv2=0.5d0*(h_mart+h_bain)&        -0.5d0*(h_mart-h_bain)*atan(x)/atan(2.0d0)
      endif

      return
      end

      subroutine zcompcalc(id,tbulk,press,zcomp)
      implicit real*8 (a-h,o-z)
      dimension tc(5),pc(5)

      data (tc(i),i=1,5)/126.2, 191.1, 33.3, 647.4, 282.4/
      data (pc(i),i=1,5)/33.5, 45.8, 12.8, 218.3, 50.5/
      tr=tbulk/tc(id)
      pr=9.869233d-6*press/pc(id)

      a=1.0
      b=-(1.0+0.125*pr/tr)/3.
      c=9.*pr/64./tr/tr
      d=-27.*pr*pr/512./tr/tr/tr

      q=a*c-b*b
      r=0.5*(3.*a*b*c-a*a*d)-b*b*b
      qr=q*q*q+r*r

      s1=(r+sqrt(qr))**(1./3.)
      s2=(r-sqrt(qr))**(1./3.)
      zcomp=(s1+s2-b)/a

      return
      end

      subroutine elg4(n0,n1,n2,n3,n4,ne1rw,nd1rw,ned1rw,&     nrow,ndrw,nedrw,iarray)
      dimension iarray (1000,5)

      iarray(1,1) = n0
      iarray(1,2) = n1
      iarray(1,3) = n2
      iarray(1,4) = n3
      iarray(1,5) = n4

      do 200 i=2,ne1rw
      iarray(i,1) = iarray((i-1),1)+ned1rw
      do 200 l=2,5
200   iarray(i,l) = iarray((i-1),l)+nd1rw

      do 400 j=2,nrow
      do 400 i=1,ne1rw
      iarray(ne1rw*(j-1)+i,1) = iarray(ne1rw*(j-2)+i,1)+nedrw
      do 400 l=2,5,1
400   iarray(ne1rw*(j-1)+i,l) = iarray(ne1rw*(j-2)+i,l)+ndrw

      return
      end

      subroutine ng_line(node_f,node_l,node_d,cord_f,cord_l,b_factor)

      implicit real*8 (a-h,o-z)
      dimension cord_f(2),cord_l(2),x(1000),y(1000)

      common /coord/x,y
      data none/1/,ntwo/2/
      data one/1.d0/

      nnode = (node_l-node_f)/node_d
      dr=one
      do 110 k=none,nnode-none,none
110   dr=dr+b_factor**k
      x0 = (cord_l(1)-cord_f(1))/dr
      y0 = (cord_l(2)-cord_f(2))/dr

      x(node_f+node_d) = cord_f(1)+x0
      y(node_f+node_d) = cord_f(2)+y0
      do 210 k=ntwo,nnode-none,none
      nd = node_f+node_d*k
      x(nd) = x(nd-node_d)+x0*b_factor**(k-none)
      y(nd) = y(nd-node_d)+y0*b_factor**(k-none)
210   continue

      return
      end

      subroutine assemble(nnode,inode,coord,nelem,ielem,&     inelem,ael,x_g,stiff)

      implicit real*8 (a-h,o-z)
      dimension :: inode(nnode),coord(2,nnode)
      dimension :: ielem(4,nelem),inelem(4,nelem),x_g(2,4,nelem)
      dimension :: ael(4,nelem),stiff(4,4,nelem)
      dimension :: x(2,4)
      dimension :: phi0(4,4),dphi(4,4,2),n_e(4)
      dimension :: r_a(2,2),r_a_e(2,2,4),vj_e(4,nelem)
    real*8 value

    value = 0.0d0
      do i=1,4
         do j=1,4
            call phi(0,value,i,j,0)
            phi0(i,j)=value
            do k=1,2
               call phi(1,value,i,j,k)
               dphi(i,j,k)=value
            enddo
         enddo
      enddo

      do i=1,nelem
         do ii=1,4
            n_e(ii)=inelem(ii,i)
         enddo

         do jj=1,4
            do ii=1,2
               x(ii,jj)=coord(ii,n_e(jj))
            enddo
         enddo

         do ii=1,4
            do jj=1,2
               sum=0.0d0
               do kk=1,4
                  sum=sum+x(jj,kk)*phi0(kk,ii)
               enddo
               x_g(jj,ii,i)=sum
            enddo
         enddo

         do ii=1,4
            call jacob(x,vj,r_a,ii,dphi)
            vj_e(ii,i)=vj
            do jj=1,2
               do kk=1,2
                  r_a_e(jj,kk,ii)=r_a(jj,kk)
               enddo
            enddo
         enddo

         do ii=1,4
            ael(ii,i)=x_g(1,1,i)*phi0(ii,1)*vj_e(1,i)
     &           +x_g(1,2,i)*phi0(ii,2)*vj_e(2,i)
     &           +x_g(1,3,i)*phi0(ii,3)*vj_e(3,i)
     &           +x_g(1,4,i)*phi0(ii,4)*vj_e(4,i)
         enddo

         do ii=1,4
            do jj=ii,4
               sum=0.0d0
               do kk=1,4
                  dphidx_i=r_a_e(1,1,kk)*dphi(ii,kk,1)
     &                 +r_a_e(1,2,kk)*dphi(ii,kk,2)
                  dphidx_j=r_a_e(1,1,kk)*dphi(jj,kk,1)
     &                 +r_a_e(1,2,kk)*dphi(jj,kk,2)
                  dphidy_i=r_a_e(2,1,kk)*dphi(ii,kk,1)
     &                 +r_a_e(2,2,kk)*dphi(ii,kk,2)
                  dphidy_j=r_a_e(2,1,kk)*dphi(jj,kk,1)
     &                 +r_a_e(2,2,kk)*dphi(jj,kk,2)
                  dphidx=dphidx_i*dphidx_j
                  dphidy=dphidy_i*dphidy_j
                  sum=x_g(1,kk,i)*(dphidx+dphidy)*vj_e(kk,i)+sum
               enddo
               stiff(ii,jj,i)=sum
               stiff(jj,ii,i)=stiff(ii,jj,i)
            enddo
         enddo
      enddo

      return
      end
      subroutine post1(iaba,kpn,cctime,crate,t_kp_max,kp)
      implicit real*8 (a-h,o-z)
      dimension cctime(kpn),crate(kpn),t_kp_max(kpn),kp(kpn)
      include 'common.f'

      i1=0
      i2=0
      i3=0
      if (cctime(1).gt.1.0d-3) i1=1
      if (cctime(2).gt.1.0d-3) i2=1
      if (cctime(3).gt.1.0d-3) i3=1
      if ((i1+i2+i3) .eq. 0) then
         write(10,*)  'Model fails to give a cooling time'
      else
         coolt=(cctime(1)*i2*2.0d0
     &        +cctime(2)*i2*4.0d0
     &        +cctime(3)*i3*2.0d0)
     &        /(i1*2.0d0+i2*4.0d0+i3*2.0d0)
Comment Check
         write(12,*) 'cctime(1)=',cctime(1)
         write(12,*) 'cctime(2)=',cctime(2)
         write(12,*) 'cctime(3)=',cctime(3)
         write(12,*) 'cooltime=',coolt

      end if

      crate_max=0.0d0
      do i=1,3
         if (dabs(crate(i)) .gt. 1.0d-3) then
            crate_max=dmax1(crate_max,dabs(crate(i)))
         endif
      end do
      write(12,*) 'crate(1)=',crate(1)
      write(12,*) 'crate(2)=',crate(2)
      write(12,*) 'crate(3)=',crate(3)
      write(12,*) 'max. cool rate=',crate_max
      if (crate_max .lt. 1.0d-3) then
         write(10,*) 'model fails to give a cooling rate'
      end if

      if (carbon_eq .gt. 1.0d-3) then
         call yurioka(cctime(3),hv)
       call yurioka2(cctime(1),hv2)

       if ( id_conf .eq. 3 ) then

        write(10,5) cctime(1),cctime(1),crate_max,t_kp_max(kpn),hv,0.0  !12/16/2009 !02/24/2010
       else

         write(10,5) cctime(3),cctime(1),crate_max,t_kp_max(kpn),hv,hv2
       end if
      else
         write(10,7) coolt,crate_max,t_kp_max(kpn),'Not Requested'
      endif

      do i=6,12,6
         write(i,3) 'ANALYSIS COMPLETED SUCCESSFULLY!'
         write(i,3), 'T85 TIME (sec)','COOL. RATE (K/sec)', &        'MAX. TEMP (C)','MAX. HAZ HARDNESS (Hv)'
         if (carbon_eq .gt. 1.0d-3) then
            write(i,5) coolt,crate_max,t_kp_max(kpn)-273.,hv
         else
            write(i,7) coolt,crate_max,t_kp_max(kpn)-273.,&           'Not Requested'
         endif
      enddo

    if ( id_conf .eq. 3 ) then
    write(10,8) CCTIME(1),CCTIME(1),CCTIME(1)           !02/24/2010
    else
    write(10,8) CCTIME(1),CCTIME(2),CCTIME(3)
    end if

 3    format(1X,A,2x,A,2x,A,2x,A)
 5    format(1x,G12.5,5x,G12.5,5x,G12.5,5x,G12.5,6x,G12.5,6x,G12.5)         !08/26/2009
 7    format(1x,G12.5,5x,G12.5,5x,G12.5,5x,A)
 8    FORMAT(1X,3g12.5)
      return
      end

      subroutine post3(iaba,kpn,cctime,crate,t_kp_max,kp,hdmax1,hdmax2)
      implicit real*8 (a-h,o-z)
      dimension cctime(kpn),crate(kpn),t_kp_max(kpn),kp(kpn)
      include 'common.f'

      i1=0
      i2=0
      i3=0
      if (cctime(1).gt.1.0d-3) i1=1
      if (cctime(2).gt.1.0d-3) i2=1
      if (cctime(3).gt.1.0d-3) i3=1
      if ((i1+i2+i3) .eq. 0) then
         write(10,*)  'Model fails to give a cooling time'
      else
         coolt=(cctime(1)*i2*2.0d0
     &        +cctime(2)*i2*4.0d0
     &        +cctime(3)*i3*2.0d0)
     &        /(i1*2.0d0+i2*4.0d0+i3*2.0d0)
         write(12,*) 'cctime(1)=',cctime(1)
         write(12,*) 'cctime(2)=',cctime(2)
         write(12,*) 'cctime(3)=',cctime(3)
         write(12,*) 'cooltime=',coolt
      end if

      crate_max=0.0d0
      do i=1,3
         if (dabs(crate(i)) .gt. 1.0d-3) then
            crate_max=dmax1(crate_max,dabs(crate(i)))
         endif
      end do
      write(12,*) 'crate(1)=',crate(1)
      write(12,*) 'crate(2)=',crate(2)
      write(12,*) 'crate(3)=',crate(3)
      write(12,*) 'max. cool rate=',crate_max
      if (crate_max .lt. 1.0d-3) then
         write(10,*) 'model fails to give a cooling rate'
      end if

      if (carbon_eq .gt. 1.0d-3) then
         call yurioka(cctime(3),hv)
       call yurioka2(cctime(1),hv2)
       if ( id_conf .eq. 3 ) then
          write(10,5) cctime(1),cctime(1),crate_max,t_kp_max(kpn),      !12/16/2009 !02/24/2010&            hdmax1,0.0
       else
      if(id_conf.eq.1) kpn=10
          write(10,5) cctime(3),cctime(1),crate_max,t_kp_max(kpn),&        hdmax1,hdmax2            !08/26/2009
       end if
      else
         write(10,7) coolt,crate_max,t_kp_max(kpn),'Not Requested'
      endif

      do i=6,12,6
         write(i,3) 'ANALYSIS COMPLETED SUCCESSFULLY!'
         write(i,3), 'T85 TIME (sec)','COOL. RATE (K/sec)',&        'MAX. TEMP (C)','MAX. HAZ HARDNESS (Hv)'
         if (carbon_eq .gt. 1.0d-3) then
            write(i,5) coolt,crate_max,t_kp_max(kpn)-273.,hv
         else
            write(i,7) coolt,crate_max,t_kp_max(kpn)-273.,&           'Not Requested'
         endif
      enddo

      if ( id_conf .eq. 3 ) then
       WRITE(10,*) "Yurioka: ", hv,0.0
         testt=1.0
       WRITE(10,*) "Kirkaldy: ", hdmax1,0.0
    else
       WRITE(10,*) "Yurioka: ", hv, hv2
       WRITE(10,*) "Kirkaldy: ", hdmax1, hdmax2
    end if

    if ( id_conf .eq. 3 ) then
    write(10,8) CCTIME(1),CCTIME(1),CCTIME(1)           !02/24/2010
    else
    write(10,8) CCTIME(1),CCTIME(2),CCTIME(3)
    end if

 3    format(1X,A,2x,A,2x,A,2x,A)
 5    format(1x,G12.5,5x,G12.5,5x,G12.5,5x,G12.5,6x,2G12.5)
 7    format(1x,G12.5,5x,G12.5,5x,G12.5,5x,A)
 8    FORMAT(1X,3g12.5)
    return
    end

      subroutine post2(iaba,kpn,cctime,crate,t_kp_max,kp)
      implicit real*8 (a-h,o-z)
      dimension cctime(kpn),crate(kpn),t_kp_max(kpn),kp(kpn)
      include 'common.f'

      write(10,7) cctime(1),abs(crate(1)),t_kp_max(2)

      do i=6,12,6
         write(i,3) 'ANALYSIS COMPLETED SUCCESSFULLY!'
         write(i,3), 'COOL. TIME (sec)','COOL. RATE (K/sec)', &        'MAX. TEMP (C)'
         write(i,7) cctime(1),abs(crate(1)),t_kp_max(2)-273.
      enddo

 3    format(1X,A,2x,A,2x,A,2x,A)
 5    format(1x,G12.5,5x,G12.5,5x,G12.5,6x,G12.5)
 7    format(3x,G12.5,5x,G12.5,5x,G12.5,5x,A)

      return

      end

      subroutine PrepareForMicrostructureCalculation(NNODE)
      implicit real*8 (a-h,o-z)
    parameter ( nn = 10000 )
    include "common.f"

    real*8 C, Mn, Si, Ni, Cr, Mo, Cu, W, V, Nb, Ti, Al, N
    REAL*8 C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    real*8 AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1
     &  ,Tliquid1,Tsolid1
    real*8 AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2
     &  ,Tliquid2,Tsolid2
    real*8 GS0,G0,hv0,gsmax
    real*8 XFE1,XPE1, xfe2,xpe2
    REAL*8 XA(NN),XF(NN),XP(NN),XB(NN),XM(NN),GS(NN),G(NN),HV(NN),XL(NN)
    real*8 fc1,pc1,bc1,FC2,PC2,BC2
    real*8 xcompos(13), krate(3)
    dimension nodeflag(nn)

    COMMON/COMP/C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
    COMMON/COMP2/C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    COMMON/TEMPS/AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1
     &  ,Tliquid1,Tsolid1
    COMMON/TEMPS2/AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2
     &  ,Tliquid2,Tsolid2
    COMMON/STRUCT/XA,XF,XP,XB,XM,GS,G,HV,XL
    COMMON/KCONST/fc1,pc1,bc1
    COMMON/KCONST2/fc2,pc2,bc2
    COMMON/BASELINE/XFE1,XPE1
    COMMON/BASELINE2/XFE2,XPE2
    COMMON/KINI/GS0,G0,hv0,gsmax
    common/nodetype/nodeflag

    XCOMPOS(1) = w_c
    XCOMPOS(2) = w_mn
    XCOMPOS(3) = w_si
    XCOMPOS(4) = w_ni
    XCOMPOS(5) = w_cr
    XCOMPOS(6) = w_mo
    XCOMPOS(7) = w_cu
    XCOMPOS(8) = 0.0d0                   !!!!w_w, no w_w, hard-code it
    XCOMPOS(9) = w_v
    XCOMPOS(10) = w_nb
    XCOMPOS(11) = 0.0d0                     !!!!w_ti, no w_ti, hard-code it
    XCOMPOS(12) = 0.0d0                  !!!!w_al, no w_al, hard-code it
    XCOMPOS(13) = w_n
    C = w_c
    Mn = w_mn
    Si = w_si
    Ni = w_ni
    Cr = w_cr
    Mo = w_mo
    Cu = w_cu
    W = 0.0d0                   !!!!w_w, no w_w, hard-code it
    V = w_v
    Nb = w_nb
    Ti = 0.0d0                  !!!!w_ti, no w_ti, hard-code it
    Al = 0.0d0                  !!!!w_al, no w_al, hard-code it
    N = w_n

    CALL THERMODYNAMICS(XCOMPOS,XFE1,XPE1,TMS1,TBS1,AE1_1,
     &  AE3_1,TDISS1,AE4_1,Tsolid1,Tliquid1,as1_1,krate)
    fc1 = krate(1)
    pc1 = krate(2)
    bc1 = krate(3)
    AC1_1=AE1_1+0.d0
    AC3_1=AE3_1+0.d0
    IF (TDISS1 .LT. AC3_1) THEN
       TDISS1=AC3_1+0.d0
    ELSE
       TDISS1=AC3_1+100.d0
    ENDIF

    XCOMPOS(1) = w_c2
    XCOMPOS(2) = w_mn2
    XCOMPOS(3) = w_si2
    XCOMPOS(4) = w_ni2
    XCOMPOS(5) = w_cr2
    XCOMPOS(6) = w_mo2
    XCOMPOS(7) = w_cu2
    XCOMPOS(8) = 0.0d0                   !!!!w_w, no w_w, hard-code it
    XCOMPOS(9) = w_v2
    XCOMPOS(10) = w_nb2
    XCOMPOS(11) = 0.0d0                     !!!!w_ti, no w_ti, hard-code it
    XCOMPOS(12) = 0.0d0                  !!!!w_al, no w_al, hard-code it
    XCOMPOS(13) = w_n2
    C2 = w_c2
    Mn2 = w_mn2
    Si2 = w_si2
    Ni2 = w_ni2
    Cr2 = w_cr2
    Mo2 = w_mo2
    Cu2 = w_cu2
    W2 = 0.0d0                   !!!!w_w, no w_w, hard-code it
    V2 = w_v2
    Nb2 = w_nb2
    Ti2 = 0.0d0                     !!!!w_ti, no w_ti, hard-code it
    Al2 = 0.0d0                  !!!!w_al, no w_al, hard-code it
    N2 = w_n2

    CALL THERMODYNAMICS(XCOMPOS,XFE2,XPE2,TMS2,TBS2,AE1_2,AE3_2,TDISS2
     &  ,AE4_2,Tsolid2,Tliquid2,as1_2,krate)
    fc2 = krate(1)
    pc2 = krate(2)
    bc2 = krate(3)
    AC1_2=AE1_2+0.d0
    AC3_2=AE3_2+0.d0
    IF (TDISS2 .LT. AC3_2) THEN
       TDISS2=AC3_2+0.d0
    ELSE
       TDISS2=AC3_2+100.d0
    ENDIF


    G01=1.0D0+2.0D0*DLOG(0.254D0/GS01)/DLOG(2.0D0)
    G02=1.0D0+2.0D0*DLOG(0.254D0/GS02)/DLOG(2.0D0)
    hv0=165.D0
    gsmax=0.3D0

    DO I=1,NNODE,1
       HV(I)=hv0
       if (nodeflag(i) .eq. 1 ) then
        GS(I)=GS02*1.0d-3
        G(I)=G02
       else
        GS(I) = GS01*1.0d-3
        G(I) = G01
       end if
       XF(I)=XFE
       XP(I)=1.0D0-XPE
       XB(I)=0.0D0
       XM(I)=0.0D0
       XA(I)=0.0D0
       XL(I)=0.0D0
      END DO
    return
    end

      subroutine MicrostructureTransition(time1, time2, temp1,&     temp2, peak, vr,nnode)
      implicit double precision (a-h,o-z)
    parameter (nn=10000)
    REAL*8 temp1(*), temp2(*), peak(nn), vr(nn)
    REAL*8 tempC1(NN), tempC2(NN)

    real*8 time1, time2, deltat
    real*8 C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
    REAL*8 C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    real*8 AC1,AC3,TDISS,AE4,AE3,AE1,TBS,TMS,as1,Tliquid,Tsolid
    real*8 AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2&    ,Tliquid2,Tsolid2
    real*8 XFE1,XPE1,XFE2,XPE2
    REAL*8 XA(NN),XF(NN),XP(NN),XB(NN),XM(NN),GS(NN),G(NN),HV(NN)&  ,XL(NN)
    real*8 fc1,pc1,bc1,FC2,PC2,BC2
    real*8 GS0,G0,hv0,gsmax

    COMMON/COMP/C,Mn,Si,Ni,Cr,Mo,Cu,W,V,Nb,Ti,Al,N
    COMMON/COMP2/C2,Mn2,Si2,Ni2,Cr2,Mo2,Cu2,W2,V2,Nb2,Ti2,Al2,N2
    COMMON/TEMPS/AC1_1,AC3_1,TDISS1,AE4_1,AE3_1,AE1_1,TBS1,TMS1,as1_1&  ,Tliquid1,Tsolid1
    COMMON/TEMPS2/AC1_2,AC3_2,TDISS2,AE4_2,AE3_2,AE1_2,TBS2,TMS2,as1_2& ,Tliquid2,Tsolid2
    COMMON/STRUCT/XA,XF,XP,XB,XM,GS,G,HV,XL
    COMMON/KCONST/fc1,pc1,bc1
    COMMON/KCONST2/fc2,pc2,bc2
    COMMON/BASELINE/XFE1,XPE1
    COMMON/BASELINE2/XFE2,XPE2
      COMMON/KINI/GS0,G0,hv0,gsmax


    deltat = (time2-time1)/1.0d2
    DO I = 1, NNODE
        TEMPC1(I) = TEMP1(I) - 273.15D0
        TEMPC2(I) = TEMP2(I) - 273.15D0
    END DO
    CALL TRANSFORMATION(TIME1,TEMPC1,TIME2,TEMPC2,DELTAT,PEAK,VR,NNODE)

    return
    end