SUBROUTINE INIT ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * * Version 0.4 - Last Modified 5/25/04 * * * ********************************************************************* INCLUDE "fpvar.inc" INTEGER INPUT REAL JUNK,O(MOZ),N2(MOZ),O2(MOZ),RHO,TEMP,OPT REAL N0,ANGLE,HH(ST:NZ),BOLTZ,GRAV *********************************************************************** * * * Initiates software, and asks several user input options * * * *********************************************************************** WRITE(6,*) WRITE(6,*) WRITE(6,*) 'Gravity Wave Simulator' WRITE(6,*) WRITE(6,*) WRITE(6,*) 'Menu:' WRITE(6,*) ' (1) Test simulation with zero seconds.' WRITE(6,*) ' (2) Test simulation with 1000 seconds.' WRITE(6,*) ' (3) More Options.' READ(5,*) CHOICE 100 IF (CHOICE.EQ.1) THEN TMAX=0. LX=300000. LZ=100000. ELSE IF (CHOICE.EQ.2) THEN TMAX=1000. LX=300000. LZ=100000. ELSE IF (CHOICE.EQ.3) THEN WRITE(6,*) 'Enter time of simulation in seconds:' READ(5,*) TMAX WRITE(6,*) WRITE(6,*) 'Enter horizontal wavelength in km:' READ(5,*) LX WRITE(6,*) LX = LX*1000. WRITE(6,*) 'Enter vertical wavelength in km:' READ(5,*) LZ WRITE(6,*) LZ = LZ*1000. ELSE WRITE(6,*) 'You have chosen poorly, choose again.' READ(5,*) CHOICE GOTO 100 END IF END IF END IF *********************************************************************** * * * Initializes common variables used throughout the simulation. * * * *********************************************************************** T = 0. II = 10 JJ = 0 KK = 1 MOZ = 121 MOX = 101 Z0 = 480000. X0 = 2.*LX DZ = Z0/REAL(MOZ-1) DX = X0/REAL(MOX-1) B = 4.5E-5 GG = 0.000001 GRAV=9.8 BOLTZ=1.381E-23 DPERP = 1.0 EOX = 1.91E-8 EOY = 1.8E-8 EOZ = 1.91E-8 UOX = 0.00000001 UOZ = 0.00000001 UEBX = 0.00000002 U1X0 = -12.0 U1Z0 = 4.0 W = 2.0*PI/2400.0 XK = -2.0*PI/LX ZK = -2.0*PI/LZ ZKI = 5.2E-6 UKI = 2.75E-6 ********************************************************************** * * * Allows user to select the case from the Huang paper. * * 1 -> Uniform Gravity Wave, Day * * 2 -> Amplitude varies w/ height, Day * * 3 -> Amplitude, Wavelength varies, Day * * 4 -> Amplitude, Wavelength varies, Night * * all else -> No Gravity wave applied to the atmosphere * * * ********************************************************************** WRITE(6,*) WRITE(6,*)'Enter the Case Number:' READ(5,*) INPUT IF (INPUT.EQ.1) THEN PRODU0=3.E8 DO 101, I=ST,NZ ZKR(I)=ZK U1Z0H(I)=U1Z0 101 U1X0H(I)=U1X0 ELSE IF (INPUT.EQ.2) THEN PRODU0=3.E8 DO 102, I=ST,NZ ZZ=DZ*REAL(I-11) ZKR(I)=ZK U1Z0H(I)=U1Z0*EXP(UKI*ZZ) 102 U1X0H(I)=U1X0*EXP(UKI*ZZ) ELSE IF (INPUT.EQ.3) THEN PRODU0=3.E8 DO 103, I=ST,NZ ZZ=DZ*REAL(I-11) ZKR(I)=ZK*EXP(-ZKI*ZZ) U1Z0H(I)=U1Z0*EXP(UKI*ZZ) 103 U1X0H(I)=U1X0*EXP(UKI*ZZ) ELSE IF (INPUT.EQ.4) THEN PRODU0=3E-16 DO 104, I=ST,NZ ZZ=DZ*REAL(I-11) ZKR(I)=ZK*EXP(-ZKI*ZZ) U1Z0H(I)=U1Z0*EXP(UKI*ZZ) 104 U1X0H(I)=U1X0*EXP(UKI*ZZ) ELSE PRODU0=3.E8 DO 105, I=ST,NZ ZKR(I)=ZK U1Z0H(I)=0.0 105 U1X0H(I)=0.0 END IF END IF END IF END IF ********************************************************************* * * * Reads the background profile from data generated by * * the MSIS-90 model. * * * ********************************************************************* OPEN(STATUS='UNKNOWN',UNIT=12,FILE='msis.dat') READ(12,*) ANGLE ANGLE=90-ANGLE DO 110, I=1,MOZ READ(12,*)ALT(I),O(I),N2(I),O2(I),RHO,TEMP,JUNK N0(I)=O(I)+N2(I)+O2(I) 110 HH(I)=BOLTZ*N0(I)*TEMP/(RHO*GRAV) CLOSE(12) DO 115, I=1,MOZ DO 115, J=ST,NX 115 N(I,J)=N0(I) SINI = SIN(ANGLE*PI/180.0) COSI = COS(ANGLE*PI/180.0) ********************************************************************* * * * Give height profiles of production rate 'produ(i)' and * * recombination coefficient 'recom(i)'. * * * ********************************************************************* OPEN(12,FILE='rate.dat', STATUS='unknown') WRITE(12,*)'TITLE="Rates"' WRITE(12,*)'VARIABLES="Altitude (km)","Prod","Recom","Opt"' WRITE(12,*)'ZONE I=',MOZ RECOM0=0.00015 DO 120, I=1,MOZ OPT=2.-EXP(-(ALT(I)-325.)/HH(I)) IF(OPT.LT.0.) OPT=0. PRODU(I)=PRODU0*EXP(-(ALT(I)-325.0)/HH(I))*OPT RECOM(I)=RECOM0*EXP(-1.75*(ALT(I)-325.0)/HH(I)) 120 WRITE(12,*) ALT(I),PRODU(I),RECOM(I),OPT CLOSE(12) ********************************************************************* * * * The initial potential profile is taken to be zero. * * The ion-neutral collision frequency vin(i) and electron- * * neutral collision frequency ven(i) are given. * * * ********************************************************************* DO 130, I=ST,NZ DO 130, J=ST,NX 130 PHI(I,J)=0.0 VIN0=6.5 VEN0=210. DO 132, I=1,10 VIN(I)=VIN0*EXP(0.14*REAL(11-I)) 132 VEN(I)=VEN0*EXP(0.14*REAL(11-I)) DO 134, I=11,MOZ VIN(I)=VIN0*EXP(-0.05*REAL(I-11)) 134 VEN(I)=VEN0*EXP(-0.05*REAL(I-11)) ********************************************************************* * * * Extends vertical boundaries for collision frequencies * * and production/recombination rates by copying edge values * * * ********************************************************************* DO 136, I=ST,0 VIN(I)=VIN(1) VEN(I)=VEN(1) N0(I)=N0(1) PRODU(I)=PRODU(1) 136 RECOM(I)=PRODU(1) DO 138, I=MOZ+1,NZ VIN(I)=VIN(MOZ) VEN(I)=VEN(MOZ) N0(I)=N0(MOZ) PRODU(I)=PRODU(MOZ) 138 RECOM(I)=RECOM(MOZ) ********************************************************************* * * * Give the height profiles of rvi(i), rve(i), ri1(i), * * ri2(i), re1(i), re2(i), ti(i), and te(i). * * * ********************************************************************* do 140, i=ST,nz tem=1500.0 mi=18.0 me=1.0 omegai=4310.6/mi omegae=7.91421e6 rvi(i)=vin(i)/omegai rve(i)=ven(i)/omegae rix1(i)=(rvi(i)**2+cosi**2)/(1.0+rvi(i)**2) riz1(i)=(rvi(i)**2+sini**2)/(1.0+rvi(i)**2) ri2(i)=sini*cosi/(1.0+rvi(i)**2) rex1(i)=(rve(i)**2+cosi**2)/(1.0+rve(i)**2) rez1(i)=(rve(i)**2+sini**2)/(1.0+rve(i)**2) re2(i)=sini*cosi/(1.0+rve(i)**2) ti(i)=8254.81*tem/(mi*vin(i)) 140 te(i)=1.51567e7*tem/(me*ven(i)) ********************************************************************* * * * Optional rotuine to restart calculations with * * previously obtained data, saved in 'in1.dat' * * as initial profiles. * * For this case, the time 't' must be changed to the * * value at which the program stopped previously. * * * ********************************************************************* 160 WRITE(6,*) WRITE(6,*)'Input Method:' WRITE(6,*)' (1) - New Calculations' WRITE(6,*)' (2) - Old Calculations' READ(5,*) INPUT IF (INPUT.EQ.2) THEN OPEN(STATUS='unknown',UNIT=9,FILE='in1.dat') READ(9,*)T DO 165, I=1,MOZ DO 165, J=1,MOX 165 READ(9,*)N(I,J),PHI(I,J) CLOSE(9) END IF ********************************************************************* * * * Extend boundary condition for n and phi in the vertical * * direction. * * * * Calculations based on neighbouring cells. * * * ********************************************************************* DO 170, I=0,ST,-1 DO 170, J=1,MOX N(I,J)=N(I+2,J)-(COSI/SINI)*(N(I+1,J+1)-N(I+1,J-1))*DZ/DX PHI(I,J)=PHI(I+2,J) $ -(COSI/SINI)*(PHI(I+1,J+1)-PHI(I+1,J-1))*DZ/DX 170 CONTINUE DO 172, I=MOZ+1,NZ DO 172, J=1,MOX N(I,J)=N(I-2,J)+(COSI/SINI)*(N(I-1,J+1)-n(I-1,J-1))*DZ/DX PHI(I,J)=PHI(I-2,J) $ +(COSI/SINI)*(PHI(I-1,J+1)-PHI(I-1,J-1))*DZ/DX 172 CONTINUE ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* DO 180, J=ST,0 DO 180, I=ST,NZ N(I,J)=N(I,MOX+J-1) 180 PHI(I,J)=PHI(I,MOX+J-1) DO 182, J=MOX+1,NX DO 182, I=ST,NZ N(I,J)=N(I,J-MOX+1) 182 PHI(I,J)=PHI(I,J-MOX+1) ********************************************************************* * * * Creates troubleshooting file named "trouble.dat" * * * ********************************************************************* OPEN(UNIT=8,FILE='trouble.dat') CLOSE(8) ********************************************************************* * * * End of Subroutine * * * ********************************************************************* WRITE(6,*) WRITE(6,*)'Running Simulation' RETURN END