SUBROUTINE INIT ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * * Version 0.5 - Last Modified 7/14/04 * * * ********************************************************************* INCLUDE "fpvar.inc" INTEGER INPUT CHARACTER CJUNK,INFILE*12 REAL LX,LZ REAL TIME,LAT,LONG,YEAR,MONTH,DAY,DATE REAL LAT0,LONG0,MLAT,POLE,TILT REAL JUNK,O,N2,O2,RHO,TEMP(ST:NZ) REAL N0,NE,DIP,HH,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 1200 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=1200. 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 = 2.*LX/REAL(MOX-1) LAT0=79.3 LONG0=71.4 B = 4.5E-5 GG = 0.000001 GRAV=9.8 BOLTZ=1.381E-23 DPERP = 1. EOX = 1.91E-8 EOY = 1.8E-8 EOZ = 1.91E-8 UOX = 0.00000001 UOZ = 0.00000001 UEBX = 0.00000002 U1X0 = -12. U1Z0 = 4. W = 2.*PI/2400. XK = -2.*PI/LX ZK = -2.*PI/LZ ZKI = 5.2E-6 UKI = 2.75E-6 ********************************************************************** * * * Allows amplitude and vertical wavenumber to vary with height.* * * ********************************************************************** DO 101, I=ST,NZ ZZ=DZ*REAL(I-11) ZKR(I)=ZK*EXP(-ZKI*ZZ) U1Z0H(I)=U1Z0*EXP(UKI*ZZ) 101 U1X0H(I)=U1X0*EXP(UKI*ZZ) ********************************************************************* * * * Reads the background profile from data generated by the * * MSIS-90 model. Density is converted to per cubic meters. * * HH=kT/mg * * * ********************************************************************* c WRITE(6,*)'Enter the input file name:' c READ(5,*)INFILE infile='input.dat' OPEN(12,FILE=INFILE) READ(12,*) CJUNK,YEAR,CJUNK,MONTH,CJUNK,DAY,CJUNK,TIME READ(12,*) CJUNK,LAT,CJUNK,LONG READ(12,*) CJUNK READ(12,*)CJUNK DO 112, I=1,MOZ READ(12,*)JUNK,O,N2,O2,RHO(I),TEMP(I) N0(I)=O+N2+O2 112 HH(I)=BOLTZ*N0(I)*TEMP(I)/(RHO(I)*GRAV) READ(12,*) CJUNK DO 113, I=1,MOZ READ(12,*)JUNK,NE(I) 113 NE(I)=NE(I)*1.E6 CLOSE(12) TEMP(0)=TEMP(1) TEMP(MOZ+1)=TEMP(MOZ) DO 115, I=1,MOZ HH(I)=1./(1./HH(I)+(TEMP(I+1)-TEMP(I-1))/(2.*DZ*TEMP(I))) DO 115, J=ST,NX 115 N(I,J)=NE(I) TLOCAL=TIME*60.*60. LAT=TORAD(LAT) LAT0=TORAD(LAT0) LONG=TORAD(LONG) LONG0=TORAD(LONG0) MLAT=ASIN(SIN(LAT)*SIN(LAT0)+COS(LAT)*COS(LAT0)*COS(LONG-LONG0)) DIP=TODEG(ATAN(2.*TAN(MLAT))) LAT=TODEG(LAT) LONG=TODEG(LONG) write(6,*) year,month,day,TIME write(6,*) lat,long write(6,*) DIP SINI = SIN(TORAD(DIP)) COSI = COS(TORAD(DIP)) IF(MONTH.EQ.1.)DATE=DAY IF(MONTH.EQ.2.)DATE=DAY+31. IF(MONTH.EQ.3.)DATE=DAY+59. IF(MONTH.EQ.4.)DATE=DAY+90. IF(MONTH.EQ.5.)DATE=DAY+120. IF(MONTH.EQ.6.)DATE=DAY+151. IF(MONTH.EQ.7.)DATE=DAY+181. IF(MONTH.EQ.8.)DATE=DAY+212. IF(MONTH.EQ.9.)DATE=DAY+243. IF(MONTH.EQ.10.)DATE=DAY+273. IF(MONTH.EQ.11.)DATE=DAY+304. IF(MONTH.EQ.12.)DATE=DAY+334. SOLSTICE=172. TILT=23.45 POLE=TILT*COS((DATE-SOLSTICE)*2.*PI/365.) SEASONAL=COS(TORAD(LAT-POLE)) ********************************************************************* * * * 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) NE(I)=NE(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) NE(I)=NE(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*TIM/(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