        SUBROUTINE OUTPUT

**********************************************************************
*                                                                    *
*       This subroutine writes the data for density and potential    *
*       for the initial calculation and every tenth iteration of     *
*       the do while (t.le.tt) loop.                                 *
*       The format for the files is Tecplot (www.tecplot.com). It    *
*       also generates a script file for each iteration.             *
*                                                                    *
**********************************************************************

        PARAMETER (NZ=126,NX=106)
        REAL N,N0,ZZ,XX
        REAL ETIME, ELAPSED(2)
        CHARACTER OUTFILE*9, ALPHA*10
        COMMON N(-5:NZ,-5:NX),PHI(-5:NZ,-5:NX)
        COMMON /INTS/II,JJ,KK
        COMMON /WX1/DZ,DX,Z0,X0,MOZ,MOX
        COMMON /WAVELENGTHS/LX,LZ
        COMMON /TIMES/DT,T,TT
        REAL AVERN(MOZ), DELTAN(MOZ,MOX)

********************************************************************
*                                                                  *
*       String Operation to number each file used                  *
*                                                                  *
********************************************************************

        ALPHA='0123456789'
        JJ=JJ+1
        II=1
        IF (JJ.EQ.11) THEN
           JJ=1
           KK=KK+1
           END IF
        OUTFILE(1:3)='out'
        OUTFILE(4:4)=ALPHA(KK:KK)
        OUTFILE(5:5)=ALPHA(JJ:JJ)
        OUTFILE(6:9)='.dat'

*********************************************************************
*                                                                   *
*       Calculates average densities and relative perturbations     *
*                                                                   *
*********************************************************************

        DO 404 I=1,MOZ
           DO 404 J=1,MOX
 404          DELTAN(I,J)=(N(I,J)-Ne(I))/Ne(I)*100.

*********************************************************************
*                                                                   *
*       Writes data to a Tecplot-friendly file (www.tecplot.com).   *
*       Data is printed "x-axis, altitude, density, perturbations,  *
*       potential".                                                 *
*                                                                   *
*********************************************************************

        OPEN(STATUS='UNKNOWN',UNIT=9,FILE=OUTFILE)
        WRITE(9,*) 'TITLE="fp.f ~ t =',T,' s"'
        WRITE(9,*) 'VARIABLES="Latitude","Altitude","Density (m<sup>-3</
     $sup>)","<greek>d</greek>N/N (%)","Potential"'
        WRITE(9,*) 'ZONE F=POINT I=',mox,',J=101'
        DO 420, J=1,MOX
           DO 420, I=11,111
              XX=(DX*REAL(J-1)-LX)/1000.
              ZZ=160.0+DZ*REAL(I-1)/1000.
 420          WRITE(9,*)XX,ZZ,N(I,J),DELTAN(I,J),PHI(I,J)
              CLOSE(9) 
        RETURN
        END


