subroutine initx(l,jdim,jmax,ndim,io,x,t,u)

c General information:

c     Subroutine initx creates a mesh of x-values for the computation
c     of the unsteady heat conduction in a bar of length l.

c     Setting equally spaced x-values is not such a big deal, but initx
c     will automatically add mesh segments if the boundary conditions
c     require it.

c     Copyright 1996 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 12/16/96

c Arguments:

c     Avoid typos:
      implicit none

c     Input: length of the bar (0 < l):
      double precision l

c     Input: declared maximum mesh point index in arrays x and u:
      integer jdim

c     Input: actual maximum mesh point index:
      integer jmax

c     Input: declared maximum time level index in arrays t and u:
      integer ndim

c     Input: I/O unit of an already open output file or zero:
      integer io
c        Subroutine initx only writes to this unit if io is positive.

c     Output: array of the x-positions of the mesh points:
      double precision x(0:jdim)

c     Unused but needed as argument: array of the mesh times t:
      double precision t(0:ndim)

c     Unused but needed as argument: temperature values:
      double precision u(0:jdim,0:ndim)

c External variables and info for compiling or changing subroutine initx:

c     Subroutine initx requires a separate subroutine:
c           bc1(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u)
c     to tell it the left hand boundary condition.
c     For an example subroutine bc1 see:
c           ../../bc_1/diri_exa/bc_1.f
c     For a template for writing new subroutines bc1 see:
c           ../../bc_1/template/bc_1.f
      external bc1

c     Subroutine initx requires a separate subroutine:
c           bc2(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u)
c     to tell it the right hand boundary condition.
c     For an example subroutine bc2 see:
c           ../../bc_2/diri_exa/bc_2.f
c     For a template for writing new subroutines bc2 see:
c           ../../bc_2/template/bc_2.f
      external bc2

c     The following utility routines from ../../lib/util.f were used:

c        cwrite(module,text,io,show) writes line "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      external cwrite

c        fatal(module,text1,text2,text3) kills the program after a
c     fatal error, printing the lines "text1", "text2" and "text3".
      external fatal

c        iwrite(module,text,val,io,show) writes integer "val" and its
c     description "text" to I/O unit "io" if nonzero, and to the screen
c     if "show" is nonzero.
      external iwrite

c Local variables:

c     Numerical codes for the properties of the boundary-condition
c     subroutines bc1 for the left end of the bar and bc2 for the right
c     end.  Subroutines bc1 and bc2 must return the sum of the
c     applicable codes when queried by subroutine initx:

c     Integer code for a nonlinear boundary condition:
      integer nonlin
      parameter (nonlin=1)
c        Ignored by init_x, but of interest for implicit schemes.  If
c     bc1 or bc2 returns this code, it means that the boundary value
c     returned by bc1, respectively bc2, varies nonlinearly with changes
c     in the values of the interior and/or previous mesh points.

c     Integer code for a periodic boundary conditions:
      integer period
      parameter (period=2)
c        The boundary conditions are periodic, with the mesh point j = 0
c     physically identical to the mesh point j = jmax-1 and the mesh
c     point j = jmax physically identical to the mesh point j = 1.
c     Since periodic boundary conditions must apply at both sides, if
c     bc1 returns this code, bc2 must also return it.  Otherwise initx
c     will notify the user to select a different pair of subroutines.

c     Integer code for a mirror-point boundary treatment:
      integer mirror
      parameter (mirror=4)
c        The end mesh point is a mirror point one mesh cell outside the
c     actual bar.  In other words, the end of the bar is at j = 1,
c     respectively at j = jmax - 1.

c     Integer code for a staggered boundary treatment:
      integer staggr
      parameter (staggr=8)
c        The boundary is halfway in between the end mesh point and the
c     next interior point.  In other words, the end of the bar is at
c     j = 0.5, respectively at j = jmax - 0.5.

c     Total number of currently defined codes:
      integer codmax
      parameter (codmax=4)
c        Change this value when you add more codes!

c     Sums of the actual codes as returned by bc1 and bc2:
      integer bc1cod,bc2cod

c     Integers that contain the individual codes of bc1 and bc2
      integer bc1non,bc2non
      integer bc1per,bc2per
      integer bc1mir,bc2mir
      integer bc1stg,bc2stg

c     Mesh point index:
      integer j

c     We define a second mesh point value k = 2 j so that we can
c     refer to the mid-points of the mesh.  This allows us to deal
c     more easily with the bar ends in staggered meshes:

c     The k-value of the left end of the bar x = 0 (default value 0):
      integer kleft

c     The k-value of the right end of the bar x = l (default value 2 jmax):
      integer kright

c     The difference between the above two k-values:
      integer kdiff

c     An integer we will use to indicate an undefined value:
      integer none
      parameter (none=-123456789)

c     Constants defined to make changing precision easier:
      double precision zero
      parameter (zero=0.d0)

c Data statements:

c     Mark the ends of the bar as unknown:
      data kleft/none/,kright/none/,kdiff/none/
c        We will keep all three variables at all times current.

c Executable statements:

c     See what the properties of the boundary conditions are:

c     Get the properties of the left boundary condition:
      bc1cod=-1
      call bc1(zero,l,jdim,jmax,ndim,0,x,t,io,bc1cod,u)

c     Check for consistency:
      if(bc1cod.ge.2**codmax)call fatal('initx',
     &'Unimplemented left boundary condition:',
     &'initx does not know where to put the mesh points.',
     &'You may need to change the source code of initx.')

c     Take the properties apart from large to small:
      bc1stg=bc1cod/staggr*staggr
      bc1cod=bc1cod-bc1stg
      bc1mir=bc1cod/mirror*mirror
      bc1cod=bc1cod-bc1mir
      bc1per=bc1cod/period*period
      bc1cod=bc1cod-bc1per
      bc1non=bc1cod

c     Get the properties of the right boundary condition:
      bc2cod=-1
      call bc2(zero,l,jdim,jmax,ndim,0,x,t,io,bc2cod,u)

c     Check for consistency:
      if(bc2cod.ge.2**codmax)call fatal('initx',
     &'Unimplemented right boundary condition:',
     &'initx does not know where to put the mesh points.',
     &'You may need to change the source code of initx.')

c     Take the properties apart from large to small:
      bc2stg=bc2cod/staggr*staggr
      bc2cod=bc2cod-bc2stg
      bc2mir=bc2cod/mirror*mirror
      bc2cod=bc2cod-bc2mir
      bc2per=bc2cod/period*period
      bc2cod=bc2cod-bc2per
      bc2non=bc2cod

c Handle periodicity:

c     Cannot have a single periodic boundary condition at one end only:
      if(bc1per+bc2per.eq.period)call fatal('initx',
     &'Both boundary conditions must be periodic or none.',
     &'Select a different pair of boundary conditions.',' ')

c     Periodicity requires a mesh that is one mesh interval longer than
c     the physical bar length.  This translates to two half intervals;
c     so the difference kdiff between the left and right end positions
c     must be 2 jmax - 2 instead of 2 jmax.

      if(bc1per+bc2per.eq.2*period)then

c        Check for conflicts:
         if(kdiff.ne.none .and. kdiff.ne.2*jmax-2)call fatal('initx',
     &   'Boundary conditions bc1 and bc2 have conflicting',
     &   'requirements for the length of the mesh.',
     &   'Please select a different pair of boundary conditions.')

c        Set the length of the bar in mesh half-point positions:
         kdiff=2*jmax-2

c        Update kleft and kright for the fact that we now know kdiff:
         if(kleft.ne.none)kright=kleft+kdiff
         if(kright.ne.none)kleft=kright-kdiff

      endif

c Handle mirror and staggered positions:

c     A mesh point cannot be both a mirror point and staggered:
      if(bc1mir+bc1stg.eq.mirror+staggr)call fatal('initx',
     &'Subroutine bc1 returns impossible properties.',
     &'Correct the source code of subroutine bc1.',' ')
      if(bc2mir+bc2stg.eq.mirror+staggr)call fatal('initx',
     &'Subroutine bc2 returns impossible properties.',
     &'Correct the source code of subroutine bc2.',' ')

c     Handle a mirror point position request of bc1:

      if(bc1mir.eq.mirror)then

c        Check for conflicts:
         if(kleft.ne.none .and. kleft.ne.2)call fatal('initx',
     &   'Subroutines bc1 and bc2 have incompatible demands for the',
     &   'mesh-position of the left end of the bar.',
     &   'Please select a different pair of boundary conditions.')

c        Set the desired end location:
         kleft=2

c        Update kdiff for the fact that we now know kleft:
         if(kright.ne.none)kdiff=kright-kleft

c        Update kright for the fact that we now know kleft:
         if(kdiff.ne.none)kright=kleft+kdiff

      endif

c     Handle a mirror point position request of bc2:

      if(bc2mir.eq.mirror)then

c        Check for conflicts:
         if(kright.ne.none .and. kright.ne.2*jmax-2)call fatal('initx',
     &   'Subroutines bc1 and bc2 have incompatible demands for the',
     &   'mesh-position of the right end of the bar.',
     &   'Please select a different pair of boundary conditions.')

c        Set the desired end location:
         kright=2*jmax-2

c        Update kdiff for the fact that we now know kright:
         if(kleft.ne.none)kdiff=kright-kleft

c        Update kleft for the fact that we now know kright:
         if(kdiff.ne.none)kleft=kright-kdiff

      endif

c     Handle a stagger position request of bc1:

      if(bc1stg.eq.staggr)then

c        Check for conflicts:
         if(kleft.ne.none .and. kleft.ne.1)call fatal('initx',
     &   'Subroutines bc1 and bc2 have incompatible demands for the',
     &   'mesh-position of the left end of the bar.',
     &   'Please select a different pair of boundary conditions.')

c        Set the desired end location:
         kleft=1

c        Update kdiff for the fact that we now know kleft:
         if(kright.ne.none)kdiff=kright-kleft

c        Update kright for the fact that we now know kleft:
         if(kdiff.ne.none)kright=kleft+kdiff

      endif

c     Handle a stagger position request of bc2:

      if(bc2stg.eq.staggr)then

c        Check for conflicts:
         if(kright.ne.none .and. kright.ne.2*jmax-1)call fatal('initx',
     &   'Subroutines bc1 and bc2 have incompatible demands for the',
     &   'mesh-position of the right end of the bar.',
     &   'Please select a different pair of boundary conditions.')

c        Set the desired end location:
         kright=2*jmax-1

c        Update kdiff for the fact that we now know kright:
         if(kleft.ne.none)kdiff=kright-kleft

c        Update kleft for the fact that we now know kright:
         if(kdiff.ne.none)kleft=kright-kdiff

      endif

c Choose the still undefined quantities:

c     Define the left end:
      if(kleft.eq.none)then

c        Set the default left end location:
         kleft=0

c        Update kright for the fact that we now know kleft:
         if(kdiff.ne.none)kright=kleft+kdiff

      endif

c     Define the right end:
      if(kright.eq.none)kright=2*jmax

c Create the mesh:

c     Check the mesh for feasibility:
      if(kright-kleft.lt.2)then
         call iwrite('initx',
     &   '    initx: Number of half mesh-intervals along the bar',
     &    kright-kleft,io,1)
         call fatal('initx',
     &   'At least one full mesh interval is required.',
     &   'Increase the value of jmax.',' ')
      endif

c     Set the mesh points, being careful for round-off errors
      do 850 j=0,jmax
         x(j)=(2*j-kleft+zero)/(kright-kleft+zero)*l
 850  continue

      call cwrite('initx',
     &'A mesh has been initialized by initx:',io,0)
      call cwrite('initx',
     &'- init_x.f Copyright 1996 Leon van Dommelen.',io,0)
      call iwrite('initx',
     &'- index of the last mesh point',jmax,io,0)
      if(kleft/2*2.eq.kleft)then
         call iwrite('initx',
     &   '- left end of the bar at mesh point',kleft/2,io,0)
      else
         call iwrite('initx',
     &   '- left end of the bar a half unit right of mesh point',
     &    kleft/2,io,0)
      endif
      if(kright/2*2.eq.kright)then
         call iwrite('initx',
     &   '- right end of the bar at mesh point',kright/2,io,0)
      else
         call iwrite('initx',
     &   '- right end of the bar a half unit left of mesh point',
     &    (kright+1)/2,io,0)
      endif

c     All done:
      goto 900

c Exit:

c     Jump here when done:
 900  return
      end