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

c General information:

c     Subroutine out performs the output for computations of the
c     unsteady heat conduction in a bar.

c     Subroutine out requires a separate exact solution;
c           function uexa(x,t,kappa,l,io,task).
c     The numerical data are compared to this exact solution.

c     The Fourier coefficients are determined if possible.
c     These are found by simple integration of the mesh point values;
c     as a result, both the numerical and the 'exact' coefficients
c     will have aliasing errors.

c     Most output is performed to the file open on I/O unit io. But
c     subroutine out will also create a plot file of the temperature
c     profiles, cexa_u.tmp, and plot files of the Fourier coefficients
c     at the requested times called cexa_c01.tmp, cexa_c02.tmp, ...  The
c     currently implemented plot file format is xyplot-compatible.

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 'cexa.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 cexa.in.  All lines are considered comments, unless
c     the 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 outputs
c     the solution u(j,n) (j = 0, 1, ..., jmax) if output at, or during
c     the time step before, time t(n) was requested in the input file.
c        During this call, out assumes that all input parameters have
c     meaningful values.  In particular, subroutine out assumes that it
c     can find the solution for the previous time step at n2 = n-1, or
c     at n2 = ndim in case n = 0.  If the solution at n2 is undefined,
c     force out to ignore it by making t(n2) greater or equal to t(n).

c     2. Any other value:
c        If task is positive, subroutine out merely initializes itself
c     and exits, otherwise it does nothing at all.

c Arguments:

c     Avoid typos:
      implicit none

c     Input: conduction constant (any restrictions are those of uexa):
      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, at least 2 (3 if periodic):
      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        If any time level n2 not equal to n is undefined, make sure
c     that the corresponding t(n2) is not less than t(n).

c     Input: I/O unit to write informational messages to:
      integer io
c        Subroutine out only writes 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     Subroutine out requires a separate exact solution;
c           function uexa(x,t,kappa,l,io,task).
c     to generate the exact temperature values.
c     For an example exact solution see:
c           ../../u_exa/triangle/u_exa.f
c     For a template for writing new exact solutions see:
c           ../../u_exa/template/u_exa.f
      double precision uexa
      external uexa

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     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     I/O unit for the Fourier coefficients plot files:
      integer ioc

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     The number of points to plot for the exact solution:
      integer jexa

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     Identifiers for the possible properties of the exact solution uexa:
c     Symmetry around the left boundary:
      integer syml
      parameter (syml=1)
c     Symmetry around the right boundary:
      integer symr
      parameter (symr=2)
c     Anti-symmetry around the left boundary:
      integer asyml
      parameter (asyml=4)
c     Anti-symmetry around the right boundary:
      integer asymr
      parameter (asymr=8)
c     Periodicity:
      integer perlr
      parameter (perlr=16)
c     Actual sum of the properties:
      integer symmet
c     Actual individual properties
      integer psyml,psymr,pasyml,pasymr,pperlr

c     First wave number to produce a Fourier coefficient for:
      integer k0
c     Increment in wave number from one mode to the next:
      integer kinc
c     Final wave number:
      integer kmax
c     The number of wave numbers:
      integer kcount
c     Wave number index:
      integer k
c     Common factor, (= 2 pi/period), in the arguments of the modes:
      double precision argfac
c     Point count in the integration of the coefficients:
      double precision count
c     Exact and numerical u-values in integrating the coefficients:
      double precision un,ue
c     Integrated numerical Fourier cosine and sine coefficients:
      double precision cnum,snum
c     Integrated analytical Fourier cosine and sine coefficients:
      double precision cexa,sexa

c     Current error, with its maximum and RMS values:
      double precision err,errmax,errrms

c     The name of the Fourier coefficients plot file:
      character*12 filnam

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,two,four,o4
      parameter (zero=0.d0,one=1.d0,two=2.d0,four=4.d0,o4=.25d0)
      double precision pi
      save pi

c     Data statements:
      data init/0/
      data isbig/0/
      data iou/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)
         call cwrite('out',
     &   '- evaluates errors using an exact solution uexa;',io,0)
         call cwrite('out',
     &   '- integrates the Fourier coefficients;',io,0)
         call cwrite('out',
     &   '- cexa/out.f Copyright 1996 Leon van Dommelen.',io,0)

c        Value of pi:
         pi=four*atan(one)

c        Open the input file:
         ioin=0
         call fopen('out','cexa.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','cexa_u.tmp',iou)

         endif

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

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

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)then
         call iwrite('out',
     &   '    out: Number of mesh intervals to plot',jmax2-j0,io,1)
         call fatal('out','No interior mesh points.',' ',' ')
      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):

c     Write the numerical temperature profile (the reversed sign on the
c     number of points causes xyplot to draw the points as symbols):
      write(iou,560,err=830,iostat=errnum)-(jmax2-j0+1),treq(req)
 560  format(i10,'  *(-1) point temperature profile (x,T) at 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     Select the number of points to plot for the exact solution:
      jexa=jmax2-j0
 562  if(jexa.ge.100)goto 563
      jexa=jexa*2
      goto 562
 563  continue

c     Draw the exact solution as a continuous line, using jexa+1 points:
      write(iou,564,err=832,iostat=errnum)jexa+1,treq(req)
 564  format(i9,'  point exact temperature profile (x,T) at t =',g15.7)
      write(iou,565,err=833,iostat=errnum)
     & (l*j/jexa,uexa(l*j/jexa,treq(req),kappa,l,io,0),j=0,jexa)
 565  format(g12.4,',',g12.4)

c     Find the errors in the mesh point values:

      errmax=zero
      errrms=zero

      do 580 j=j0,jmax2
c        Pointwise error:
         err=abs(u(j,n)-uexa(x(j),t(n),kappa,l,io,0))
c        Update the maximum error:
         if(err.gt.errmax)errmax=err
c        Give the error its weighting factor in the RMS sum:
         err=err/sqrt((jmax2-j0+one))
c        Update the RMS error while avoiding unnecessary overflows:
         tmp=max(err,errrms)
         if(tmp.gt.zero)errrms=tmp*sqrt((errrms/tmp)**2+(err/tmp)**2)
 580  continue

      call rwrite('out',
     &'Maximum error in the mesh point values',errmax,io,1)
      call rwrite('out',
     &'    RMS error in the mesh point values',errrms,io,0)

c     Find the wave numbers to compute Fourier coefficients for:

c     First inquire the exact solution about its properties:
      symmet=-1
      tmp=uexa(x(j0),t(n),kappa,l,io,symmet)
c     Divide away any unrecognized parameters in the returned symmet:
      symmet=symmet/syml*syml
      symmet=symmet-symmet/(2*perlr)*(2*perlr)
c     Extract the properties from large to small:
      pperlr=symmet/perlr*perlr
      symmet=symmet-pperlr
      pasymr=symmet/asymr*asymr
      symmet=symmet-pasymr
      pasyml=symmet/asyml*asyml
      symmet=symmet-pasyml
      psymr=symmet/symr*symr
      symmet=symmet-psymr
      psyml=symmet/syml*syml
      symmet=symmet-psyml

c     Wave numbers for the periodic case:
      if(pperlr.eq.perlr)then
         k0=0
         kinc=1
         kmax=(jmax2-j0-1)/2
         kcount=2*kmax+1
         if(psyml.eq.syml)kcount=kcount-kmax
         if(pasyml.eq.asyml)kcount=kcount-kmax-1
         argfac=two*pi/l
         goto 610
      endif

c     Wave numbers for two symmetry boundary conditions:
      if(psyml+psymr.eq.syml+symr)then
         k0=0
         kinc=1
         kmax=(2*(jmax2-j0)-1)/2
         kcount=kmax+1
         if(pasyml.eq.asyml)kcount=0
         argfac=pi/l
         goto 610
      endif
        
c     Wave numbers for two anti-symmetry boundary conditions:
      if(pasyml+pasymr.eq.asyml+asymr)then
         k0=1
         kinc=1
         kmax=(2*(jmax2-j0)-1)/2
         kcount=kmax
         if(psyml.eq.syml)kcount=0
         argfac=pi/l
         goto 610
      endif
        
c     Wave numbers for symmetry and anti-symmetry boundary conditions:
      if(psyml+pasymr.eq.syml+asymr .or. pasyml+psymr.eq.asyml+symr)then
         k0=1
         kinc=2
         kmax=(4*(jmax2-j0)-1)/2
         kcount=1+(kmax-1)/2
         if(pasyml+psyml.eq.asyml+syml)kcount=0
         argfac=pi/(two*l)
         goto 610
      endif
        
c     The exact solution has an unimplemented symmetry. Skip the coefficients:
      goto 650

c     Jump here after the wave numbers have been found:
 610  continue

c     Check whether we have something to do:
      if(kcount.le.0)then
         call warn('out',io,'All Fourier modes are trivial.',' ',' ')
         goto 650
      endif

c     Find the Fourier coefficients:
c        Fast Fourier transforms would be more efficient to do this,
c     but more complicated, less easy to follow.  We just crunch.

c     Initialize the RMS and maximum errors in the coefficients:

      errmax=zero
      errrms=zero

c     Create the Fourier coefficients plot file:

c     Generate a unique file name:
      if(req.gt.99)call fatal('out',
     &'Code limit for the number of output times is 99.',' ',' ')
      filnam='cexa_c00.tmp'
      filnam(7:7)=char(48+req/10-req/100*10)
      filnam(8:8)=char(48+req-req/10*10)

c     Create the file:
      ioc=0
      call fcreat('out',filnam,ioc)

c     Write the number of numerical coefficients (the reversed sign on
c     the number of points causes xyplot to draw the points as symbols):
      write(ioc,620,err=840,iostat=errnum)-kcount,treq(req)
 620  format(i10,'  *(-1) numerical Fourier coefficients at t =',g15.7)

c     Compare the numerical and exact Fourier coefficients:

      do 630 k=k0,kmax,kinc

c        Integrate the Fourier coefficients:

c        Initialize the count of the number of integration points:
         count=zero
c        Initialize the numerical Fourier coefficients:
         cnum=zero
         snum=zero
c        Initialize the exact Fourier coefficients:
         cexa=zero
         sexa=zero

c        Integrate the coefficients:
         do 625 j=j0,jmax2
c           Increment the count of the integration points:
            count=count+one
c           Exact temperature value:
            ue=uexa(x(j),treq(req),kappa,l,io,0)
c           Numerical temperature value interpolated to time treq(req):
            un=u(j,n)+shift*(u(j,n2)-u(j,n))
c           Perform a correction if the point is an end point:
            if(x(j).eq.zero .or. x(j).eq.l)then
               count=count-(one/two)
               ue=ue/two
               un=un/two
            endif
c           Increment the numerical coefficients:
            cnum=cnum+un*cos(k*argfac*x(j))
            snum=snum+un*sin(k*argfac*x(j))
c           Increment the exact coefficients:
            cexa=cexa+ue*cos(k*argfac*x(j))
            sexa=sexa+ue*sin(k*argfac*x(j))
 625     continue

c        Normalize the coefficients
         cnum=cnum/count
         snum=snum/count
         cexa=cexa/count
         sexa=sexa/count

c        Process the nontrivial numerical coefficients:

c        Process the cosine coefficient:
         if(pasyml.ne.asyml)then
c           Write it to the plot file (shift it slightly):
            write(ioc,626,err=850,iostat=errnum)k+o4,cnum
 626        format(g12.4,',',g12.4)
c           Compute the error:
            err=abs(cnum-cexa)
c           Update the maximum error:
            if(err.gt.errmax)errmax=err
c           Give the error its weighting factor in the RMS sum:
            err=err/sqrt(kcount+zero)
c           Update the RMS error while avoiding unnecessary overflows:
            tmp=max(errrms,err)
            if(tmp.gt.zero)errrms=tmp*sqrt((errrms/tmp)**2+(err/tmp)**2)
         endif

c        Process the sine coefficient:
         if(k.gt.0 .and. psyml.ne.syml)then
c           Write it to the plot file:
            write(ioc,627,err=851,iostat=errnum)k+zero,snum
 627        format(g12.4,',',g12.4)
c           Compute the error:
            err=abs(snum-sexa)
c           Update the maximum error:
            if(err.gt.errmax)errmax=err
c           Give the error its weighting factor in the RMS sum:
            err=err/sqrt(kcount+zero)
c           Update the RMS error while avoiding unnecessary overflows:
            tmp=max(errrms,err)
            if(tmp.gt.zero)errrms=tmp*sqrt((errrms/tmp)**2+(err/tmp)**2)
         endif

 630  continue

      call rwrite('out',
     &'Maximum error in the Fourier coefficients',errmax,io,0)
      call rwrite('out',
     &'    RMS error in the Fourier coefficients',errrms,io,0)

c     Compute the exact Fourier coefficients once more for plotting:

      do 640 k=k0,kmax,kinc

c        Initialize the count of the number of integration points:
         count=zero
c        Initialize the exact Fourier coefficients:
         cexa=zero
         sexa=zero

c        Integrate the coefficients:
         do 635 j=j0,jmax2
c           Increment the count of the integration points:
            count=count+one
c           Exact temperature value:
            ue=uexa(x(j),treq(req),kappa,l,io,0)
c           Perform a correction if the point is an end point:
            if(x(j).eq.zero .or. x(j).eq.l)then
               count=count-(one/two)
               ue=ue/two
            endif
c           Increment the exact coefficients:
            cexa=cexa+ue*cos(k*argfac*x(j))
            sexa=sexa+ue*sin(k*argfac*x(j))
 635     continue

c        Normalize the coefficients
         cexa=cexa/count
         sexa=sexa/count

c        Write the nontrivial exact coefficients to the plot file in
c        a format that causes xyplot to draw a spike

c        Write the cosine coefficient:
         if(pasyml.ne.asyml)
     &    write(ioc,636,err=860,iostat=errnum)k+o4,k+o4,cexa
 636     format('2'/g12.4,',0.'/g12.4,',',g12.4)

c        Write the sine coefficient:
         if(k.gt.0 .and. psyml.ne.syml)
     &    write(ioc,637,err=861,iostat=errnum)k+zero,k+zero,sexa
 637     format('2'/g12.4,',0.'/g12.4,',',g12.4)

 640  continue

c     Close the plot file:

      call fclose('out',filnam,ioc)
      call cwrite('out',
     &'Created Fourier coefficients plot file '//filnam,io,1)

c     Done with Fourier coefficients:

 650  continue

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','cexa_u.tmp',iou)
         call cwrite('out',
     &   'Created temperature-profiles plot file cexa_u.tmp',io,1)
         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.')
 832  call ioerr('out',iou,errnum,
     &'Error writing to the temperature profiles plot file,',
     &'using format 564 in the source code of out.',
     &'Check disk space and the format.')
 833  call ioerr('out',iou,errnum,
     &'Error writing to the temperature profiles plot file,',
     &'using format 565 in the source code of out.',
     &'Check disk space and the format.')
 840  call ioerr('out',ioc,errnum,
     &'Error writing to the Fourier coefficients plot file,',
     &'using format 620 in the source code of out.',
     &'Check disk space and the format.')
 850  call ioerr('out',ioc,errnum,
     &'Error writing to the Fourier coefficients plot file,',
     &'using format 626 in the source code of out.',
     &'Check disk space and the format.')
 851  call ioerr('out',ioc,errnum,
     &'Error writing to the Fourier coefficients plot file,',
     &'using format 627 in the source code of out.',
     &'Check disk space and the format.')
 860  call ioerr('out',ioc,errnum,
     &'Error writing to the Fourier coefficients plot file,',
     &'using format 636 in the source code of out.',
     &'Check disk space and the format.')
 861  call ioerr('out',ioc,errnum,
     &'Error writing to the Fourier coefficients plot file,',
     &'using format 637 in the source code of out.',
     &'Check disk space and the format.')

c     Jump here when done:
 900  return
      end