subroutine out(kappa,l,jdim,jmax,ndim,n,x,t,u,io,task)

c General information:

c     Template for writing subroutines out that perform output for the
c     unsteady heat conduction in a bar.  Follow this template to create
c     new subroutines out that produce different output.

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

c Usage information:

c     Subroutine out looks for its parameters in a file 'out.in' of
c     the form:

c        <Number of desired output times N>   <comment>
c        <Desired output time 1>              <comment>
c        <Desired output time 2>              <comment>
c        ...
c        <Desired output time N>              <comment>

c        For better readability, comment lines may be inserted anywhere
c     in file out.in.  All lines are considered comments, unless the
c     first nonblank character is a digit, sign, or a decimal point.

c     Subroutine out can be called in two different ways depending on
c     the chosen value of parameter task:

c     1. task = 0:
c        This is the standard call. In this case, subroutine out should
c     output the solution u(j,n) (j = 0, 1, ..., jmax) if output at, or
c     during the time step before, time t(n) was requested in the input
c     file out.in.
c        During such a call, subroutine out may want to interpolate to
c     the exact requested time using other computed time levels n2 in
c     addition to level n.  In doing this, subroutine out may assume
c     that the time levels u(j,n2) are stored cyclically, so that n2 =
c     n, n-1, n-2, ..., 1, 0, ndim, ndim-1, ..., n+2, n+1 are the last
c     ndim time levels computed, in order of decreasing time. Subroutine
c     out may assume that all these earlier time planes are correctly
c     defined as long as t(n2) is less than t(n).
c         This also means that if there are undefined time levels in
c     array u when calling subroutine out, the user must make sure that
c     the corresponding times are not less than t(n).

c     2. Any other value:
c        If task is positive, subroutine out should merely initialize
c     itself and exit, otherwise it should do nothing at all.

c Arguments:

c     Avoid typos:
      implicit none

c     Input: conduction constant (restrictions may be needed):
      double precision kappa

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: last time level computed:
      integer n

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

c     Input: array of the mesh times t:
      double precision t(0:ndim)

c     Input: temperature values:
      double precision u(0:jdim,0:ndim)
c        Subroutine out may assume that u(j,n2) is defined for time
c     level n2 = n, and for any other n2 for which t(n2) < t(n).

c     Input: I/O unit to write informational messages to:
      integer io
c        Subroutine out should only write to this unit if io is positive.

c     Input: task to perform:
      integer task
c        See the usage information above for more information on task.

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

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        fclose(module,filnam,io) closes the I/O unit "io" on which
c     file "filnam" is open.
      external fclose

c        fcreat(module,filnam,io) creates the new output file "filnam",
c     opening it on I/O unit "io". It selects "io" if initialized to zero.
      external fcreat

c        fopen(module,filnam,io) opens the existing file "filnam" on
c     I/O unit "io". It selects "io" if initialized to zero.
      external fopen

c        ioerr(module,io,errnum,text1,text2,text3) kills the program
c     after a fatal I/O error on I/O unit "io", printing the lines
c     "text1", "text2", and "text3", and possibly other info based on
c     Fortran IOSTAT error code "errnum".
      external ioerr
      integer errnum

c        iread(module,text,ioin,io,show) reads integer "iread" from I/O
c     unit "ioin". It writes it and its description "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      integer iread
      external iread

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        rread(module,text,ioin,io,show) reads number "rread" from I/O
c     unit "ioin". It writes it and its description "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      double precision rread
      external rread

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

c        warn(module,io,text1,text2,text3) writes the warning lines
c     "text1", "text2" and "text3" to the screen, and to I/O unit "io"
c     if nonzero.
      external warn

c Local variables:

c     ***************************************************************
c     ******* Declare any local variables that out needs here *******
c     ***************************************************************

c     An integer to keep track of whether out has been initialized:
      integer init

c     I/O unit for the input file:
      integer ioin
c     I/O unit for the temperature profiles plot file:
      integer iou

c     Mesh point index:
      integer j
c     The lowest value of j to output:
      integer j0
c     The maximum value of j to output:
      integer jmax2
c     Spacing of the j-values while printing some typical mesh points:
      integer jinc

c     Whether we warned for insufficient mesh points (we warn only once):
      integer nomesh

c     The minimum and maximum values of u:
      double precision minu,maxu
c     The j-position of the minimum and maximum:
      integer jminu,jmaxu

c     Whether we warned for pending overflow (we warn only once):
      integer isbig

c     Total number of requested times at which to do output:
      integer reqmax
      save reqmax
c     Dimensioned last output time:
      integer reqdim
      parameter (reqdim=99)
c     Index over the times at which output is requested:
      integer req
      save req
c     Index req as a character string:
      character*2 creq
c     Requested output times:
      double precision treq(reqdim)
      save treq

c     Step size while sorting the requested output times:
      integer step
c     Whether the sort changed anything:
      logical changd

c     Previously computed time level, used for interpolating the
c     requested time:
      integer n2
c     Relative weight to give the n2 time level while interpolating:
      double precision shift

c     Temporary:
      double precision tmp

c     Magnitude of u beyond which we warn for pending overflow:
      double precision big
      parameter (big=1.d15)

c     Constants, defined to simplify changes in precision:
      double precision zero,one
      parameter (zero=0.d0,one=1.d0)

c     Data statements:
      data init/0/
      data iou/0/
      data nomesh/0/
      data isbig/0/

c Executable statements:

c Ignore all negative tasks:

      if(task.lt.0)return

c Initialization during the first time that out is called:

      if(init.eq.0)then

c        Show which subroutine is being initialized:
         call cwrite('out',
     &   'Initialization of output subroutine out:',io,0)
c        ************************************************************
c        *********** Briefly describe subroutine out here ***********
c        ************************************************************
         call cwrite('out',
     &   '- template/out.f Copyright 1996 Leon van Dommelen.',io,0)
         call cwrite('out',
     &   '*** ERROR: You are using the template subroutine!',io,0)

c        Open the input file:
         ioin=0
         call fopen('out','out.in',ioin)

c        Read the number of times to output from the input file:
         reqmax=iread('out','Number of output times',ioin,io,1)
         if(reqmax.lt.0)call fatal('out',
     &   'The number of output times cannot be negative.',' ',' ')
         if(reqmax.gt.reqdim)then
            call iwrite('out',
     &      '    out:   Requested number of output times',reqmax,io,1)
            call iwrite('out',
     &      '         Dimensioned number of output times',reqdim,io,1)
            call fatal('out','Program limits exceeded.',' ',' ')
         endif

c        Read the times to output:
         if(reqmax.gt.0)then

c           Read the times from the input file using function rread:
            if(reqmax.gt.99)call fatal('out',
     &      'Code limit for the number of output times is 99.',' ',' ')
            do 120 req=1,reqmax
               creq=' '//char(48+req-req/10*10)
               if(req.gt.9)creq(1:1)=char(48+req/10-req/100*10)
               treq(req)=rread('out','Output time '//creq,ioin,0,0)
 120        continue

c           Sort the times in ascending order:
            step=reqmax
 140        step=(step+1)/2
 141        changd=.false.
            do 150 req=1,reqmax-step
               if(treq(req).le.treq(req+step))goto 150
               tmp=treq(req)
               treq(req)=treq(req+step)
               treq(req+step)=tmp
               changd=.true.
 150        continue
            if(changd)goto 141
            if(step.gt.1)goto 140

c           Write the sorted times to the output file:
            call cwrite('out','Requested output times, sorted:',io,0)
            if(io.gt.0)write(io,160,err=810,iostat=errnum)
     &       (treq(req),req=1,reqmax)
 160        format(5g15.7)

c           Create the temperature-profiles plot file:
            call fcreat('out','out_u.tmp',iou)

         endif

c        Close the input file:
         call fclose('out','out.in',ioin)

c        Initialize the first requested time to output:
         req=1

c        *************************************************************
c        ******* Read in any other parameters needed over here *******
c        *************************************************************

c        All done with the initialization:
         init=1

      endif

c Further ignore nonzero tasks:

      if(task.ne.0)return

c Check the sanity of some of the arguments:

      if(jdim.le.0 .or. ndim.lt.0 .or.
     & jmax.le.0 .or. n.lt.0 .or.
     & jmax.gt.jdim .or. n.gt.ndim)call fatal('out',
     & 'Program error: invalid arguments issued to subroutine out.',
     & 'Fix the program code.',' ')

      if(l.le.zero)then
         call rwrite('out','    out: Length of the bar l',l,io,1)
         call fatal('out','Length of the bar must be positive.',' ',' ')
      endif

c Find plotting limits for which all points are on the bar:

c     Find the lower limit j0 to plot:
      j0=-1
 210  j0=j0+1
      if(x(j0).lt.zero .and. j0.lt.jmax)goto 210

c     Find the upper limit jmax2 to plot:
      jmax2=jmax+1
 220  jmax2=jmax2-1
      if(x(jmax2).gt.l .and. jmax2.gt.0)goto 220

c     Check that there is an interval to plot:
      if(jmax2-j0.lt.1)then
         call iwrite('out',
     &   '    out: Number of mesh segments to plot',jmax2-j0,io,1)
         call fatal('out','No mesh segments to plot.',' ',' ')
      endif

c     Check whether there are interior points to plot:
      if(jmax2-j0.lt.2 .and. nomesh.eq.0)then
         call iwrite('out',
     &   '    out: Number of mesh intervals to plot',jmax2-j0,io,1)
         call warn('out',io,'No interior mesh points.',' ',' ')
         nomesh=1
      endif

c Check for pending overflow:

c     Find the minimum and maximum u values at level n:
      minu=u(j0,n)
      jminu=j0
      maxu=u(j0,n)
      jmaxu=j0
      do 250 j=j0+1,jmax2
         if(u(j,n).lt.minu)then
            minu=u(j,n)
            jminu=j
         endif
         if(u(j,n).gt.maxu)then
            maxu=u(j,n)
            jmaxu=j
         endif
 250  continue

c     Warn if too big and we did not warn before:
      if(max(-minu,maxu).gt.big .and. isbig.eq.0)then
         call rwrite('out',
     &   '     out: Minimum temperature',minu,io,1)
         call iwrite('out',
     &   '                   mesh point',jminu,io,1)
         call rwrite('out',
     &   '          Maximum temperature',maxu,io,1)
         call iwrite('out',
     &   '                   mesh point',jmaxu,io,1)
         call warn('out',io,'Pending overflow.',' ',' ')
         isbig=1
      endif

c Output all requested times:

c     See whether we want to output the current time

 500  if(req.gt.reqmax)goto 699
      if(treq(req).gt.t(n))goto 699

c     Output the time and limits of the solution:

      if(io.gt.0)write(io,550,err=820,iostat=errnum)
     & treq(req),minu,jminu,maxu,jmaxu
 550  format(/
     &'Time =',g13.5/
     &'Minimum temperature',g13.5,' at mesh point',i9/
     &'Maximum temperature',g13.5,' at mesh point',i9/
     &'Typical mesh point values:')

c     Output some typical mesh point values:

      jinc=max((jmax2-j0+2)/4,1)
      if(io.gt.0)write(io,551,err=821,iostat=errnum)
     & (x(j),j=j0,min(jmax2-1,j0+3*jinc),jinc),x(jmax2)
      if(io.gt.0)write(io,552,err=822,iostat=errnum)
     & (u(j,n),j=j0,min(jmax2-1,j0+3*jinc),jinc),u(jmax2,n)
 551  format('x: ',5g15.7)
 552  format('T: ',5g15.7)

c     Determine the interpolation to use to the requested time:

      n2=n-1
      if(n2.lt.0)n2=ndim
      if(t(n2).lt.t(n))then
         shift=(treq(req)-t(n))/(t(n2)-t(n))
      else
         n2=n
         shift=zero
      endif

c     Write the current temperature profile to the plot file (the plot
c     file is formatted for the 'xyplot' plotting program):

      write(iou,560,err=830,iostat=errnum)jmax2-j0+1,treq(req)
 560  format(i9,'  point temperature profile (x,T) at time t =',g15.7)
      write(iou,561,err=831,iostat=errnum)
     & (x(j),u(j,n)+shift*(u(j,n2)-u(j,n)),j=j0,jmax2)
 561  format(g12.4,',',g12.4)

c     ****************************************************************
c     **************** Perform additional output here ****************
c     ****************************************************************

c     Next requested time:

 690  req=req+1
      if(req.gt.reqmax)goto 699
      if(treq(req).eq.treq(req-1))goto 690
      goto 500
 699  if(req.gt.reqmax .and. iou.gt.0)then
         call fclose('out','out_u.tmp',iou)
         iou=0
      endif
      goto 900

c Exit:

c     Jump here for various I/O errors:
 810  call ioerr('out',io,errnum,
     &'Error writing to the output file,',
     &'using format 160 in the source code of out.',
     &'Check disk space and the format.')
 820  call ioerr('out',io,errnum,
     &'Error writing to the output file,',
     &'using format 550 in the source code of out.',
     &'Check disk space and the format.')
 821  call ioerr('out',io,errnum,
     &'Error writing to the output file,',
     &'using format 551 in the source code of out.',
     &'Check disk space and the format.')
 822  call ioerr('out',io,errnum,
     &'Error writing to the output file,',
     &'using format 552 in the source code of out.',
     &'Check disk space and the format.')
 830  call ioerr('out',iou,errnum,
     &'Error writing to the temperature-profiles plot file,',
     &'using format 560 in the source code of out.',
     &'Check disk space and the format.')
 831  call ioerr('out',iou,errnum,
     &'Error writing to the temperature profiles plot file,',
     &'using format 561 in the source code of out.',
     &'Check disk space and the format.')

c     Jump here when done:
 900  return
      end