      PROGRAM FP

********************************************************************* 
*                                                                   *
*       This program, named as FP including subroutines             *
*       init(fp0.5a.f), ellipse(fp0.5b.f), hyperbolic(fp0.5c.f),    *
*       and final(fp0.5d.f), and a list of functions (fp0.5e.f      *
*       calculates the spatial and temporal evolution of            *
*       ionospheric F-region perturbations generated by gravity     *
*       waves. This program also requires a user-specified profile  *
*       generated from MSIS-90 and IRI-2001 (see sample in          *
*       input.dat) and a common variable list (var.inc).            *
*                                                                   *
*       Version 0.5 - Last Modified 7/14/04                         *
*                                                                   *
*********************************************************************  

        INCLUDE "fpvar.inc"

        CALL INIT
        DO WHILE (T.LE.TMAX)

*********************************************************************
*                                                                   *
*       Initializes the gravity waves for each time iteration.      *
*                                                                   *
*********************************************************************

           DO 10, I=ST,NZ
              DO 10, J=ST,NX
                 Z=DZ*REAL(I-1)-Z0/2.
                 X=DX*REAL(J-1)-X0/2.
                 U1Z(I,J)=U1Z0H(I)*COS(XK*X+ZKR(I)*Z-W*T)
 10              U1X(I,J)=U1X0H(I)*COS(XK*X+ZKR(I)*Z-W*T)

*********************************************************************
*                                                                   *
*       Initializes production and recombination curves for         *
*       local time slot.                                            *
*                                                                   *
*********************************************************************

           RECOM0=0.00015
           PRODU0=3.E8
           PRODU0=PRODU0*SIN(PI*TLOCAL/(24.*60.*60.))*SEASONAL
	   WRITE(6,*)PRODU0
           DO 20, I=1,MOZ
              ZZ=160.+DZ*REAL(I-1)/1000.
              OPT=2.-EXP(-(ZZ-325.)/HH(I))
              IF(OPT.LT.0.) OPT=0.
              PRODU(I)=PRODU0*EXP(-(ZZ-325.)/HH(I))*OPT
 20           RECOM(I)=RECOM0*EXP(-1.75*(ZZ-325.)/HH(I))

*********************************************************************
*                                                                   *
*       Calculates new values for N and PHI                         *
*                                                                   *
*********************************************************************

           CALL HYPERBOLIC
           CALL ELLIPSE
           IF (EPC.GT.100.0) THEN
              WRITE(6,*) 'Divergent Calculation - Program Aborted!'
              STOP
              END IF

*********************************************************************
*                                                                   *
*         Writes loop values to a troubleshooting file              *
*         named "trouble.dat".                                      *
*                                                                   *
*********************************************************************

           OPEN (STATUS='OLD',ACCESS='APPEND',UNIT=8,FILE='trouble.dat')
           T=T+DT
	   TLOCAL=TLOCAL+DT
    	   WRITE(8,*)'epc=',EPC
           WRITE(8,*)'dt=',DT
           WRITE(8,*)'t=',T
           CLOSE(8)
           END DO
	CALL FINAL
        STOP
        END

**********************************************************************
*                                                                    *
*       These are the required files for the 4 subroutines           *
*                                                                    *
*       a - INIT                                                     *
*       b - ELLIPSE                                                  *
*       c - HYPERBOLIC                                               *
*       d - FINAL                                                    *
*       e - functions                                                *
*                                                                    *
**********************************************************************

        INCLUDE "fp0.5a.f"
        INCLUDE "fp0.5b.f"
        INCLUDE "fp0.5c.f"
        INCLUDE "fp0.5d.f"
        INCLUDE "fp0.5e.f"