SUBROUTINE INIT ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * ********************************************************************* PARAMETER (NZ=126,NX=106) PARAMETER (PI=3.141592654) INTEGER INPUT CHARACTER CJUNK REAL N,Ne,junk,o,o2,n2 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 /times/dt,t,tt common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx common /wx4/u1x0,u1z0,xk,zk,w,zki,uki COMMON /WAVELENGTHS/LX,LZ common /wx5/vin(-5:nz),ven(-5:nz) common /wx61/rvi(-5:nz),rix1(-5:nz),riz1(-5:nz),ri2(-5:nz) common /wx62/rve(-5:nz),rex1(-5:nz),rez1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz),sini,cosi common /wx8/epc,dperp common /wx9/ne(-5:nz) common /wx10/produ(-5:nz),recom(-5:nz) common /wx23/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) common /wx35/u1x0h(-5:nz),u1z0h(-5:nz),zkr(-5:nz) *********************************************************************** * * * Initiates software, and asks several user input options * * * *********************************************************************** WRITE(6,*) WRITE(6,*) WRITE(6,*)'Gravity Wave Simulator' WRITE(6,*) WRITE(6,*) WRITE(6,*)'Enter time of simulation in seconds:' READ(5,*)TT WRITE(6,*) WRITE(6,*)'Enter horizontal wavelength in km' READ(5,*)LX LX=LX*1000.0 WRITE(6,*)'Enter vertical wavelength in km' READ(5,*)LZ LZ=LZ*1000.0 *********************************************************************** * * * Initializes common variables used throughout the simulation. * * * *********************************************************************** T = 0.0 II = 10 JJ = 0 KK = 1 MOZ = 121 MOX = 101 Z0 = 480000.0 X0 = 2.0*LX DZ = Z0/REAL(MOZ-1) DX = X0/REAL(MOX-1) B = 4.5E-5 GG = 0.000001 SINI = SIN(68.0*PI/180.0) COSI = COS(68.0*PI/180.0) 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 * * * ********************************************************************** 100 WRITE(6,*) WRITE(6,*)'Enter the Case Number:' READ(5,*) INPUT IF (INPUT.EQ.1) THEN PRODU0=300000000.0 DO 101, I=-5,NZ ZKR(I)=ZK U1Z0H(I)=U1Z0 101 U1X0H(I)=U1X0 ELSE IF (INPUT.EQ.2) THEN PRODU0=300000000.0 DO 102, I=-5,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=300000000.0 DO 103, I=-5,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=-5,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=300000000.0 DO 105, I=-5,NZ ZKR(I)=ZK U1Z0H(I)=0.0 105 U1X0H(I)=0.0 END IF END IF END IF END IF ********************************************************************* * * * Give height profiles of production rate 'produ(i)' and * * recombination coefficient 'recom(i)'. * * * ********************************************************************* 110 recom0=0.00015 do 120, i=1,moz hh=160.0+dz*real(i-1)/1000.0 produ(i)=produ0*exp(-(hh-300.0)/60.0) recom(i)=recom0*exp(-1.75*(hh-300.0)/60.0) if (recom(i).gt.recom0) recom(i)=recom0 120 if (produ(i).gt.produ0) produ(i)=produ0 ********************************************************************* * * * 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=-5,nz do 130, j=-5,nx 130 phi(i,j)=0.0 vin0=6.5 ven0=210.0 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)) ********************************************************************* * * * 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=1,moz 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)) ********************************************************************* * * * Read the background density profile. * * * ********************************************************************* OPEN(STATUS='UNKNOWN',UNIT=12,FILE='iri.dat') READ(12,*) cjunk DO 150, I=1,MOZ READ(12,*)JUNK,NE(I) 150 NE(I)=NE(I)*1.E6 CLOSE(12) do 155, i=1,moz do 155, j=-5,nx 155 n(i,j)=ne(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 in the vertical direction. * * * ********************************************************************* do 172, i=-5,0 vin(i)=vin(1) ven(i)=ven(1) produ(i)=produ(1) recom(i)=recom(1) 172 n0(i)=n0(1) do 174, i=moz+1,nz vin(i)=vin(moz) ven(i)=ven(moz) produ(i)=produ(moz) recom(i)=recom(moz) 174 n0(i)=n0(moz) do 176, i=-5,0 do 176, j=1,mox n(0,j)=n(2,j)-(cosi/sini)*(n(1,j+1)-n(1,j-1))*dz/dx n(-1,j)=n(1,j)-(cosi/sini)*(n(0,j+1)-n(0,j-1))*dz/dx n(-2,j)=n(0,j)-(cosi/sini)*(n(-1,j+1)-n(-1,j-1))*dz/dx n(-3,j)=n(-1,j)-(cosi/sini)*(n(-2,j+1)-n(-2,j-1))*dz/dx n(-4,j)=n(-2,j)-(cosi/sini)*(n(-3,j+1)-n(-3,j-1))*dz/dx n(-5,j)=n(-3,j)-(cosi/sini)*(n(-4,j+1)-n(-4,j-1))*dz/dx phi(0,j)=phi(2,j) $ -(cosi/sini)*(phi(1,j+1)-phi(1,j-1))*dz/dx phi(-1,j)=phi(1,j) $ -(cosi/sini)*(phi(0,j+1)-phi(0,j-1))*dz/dx phi(-2,j)=phi(0,j) $ -(cosi/sini)*(phi(-1,j+1)-phi(-1,j-1))*dz/dx phi(-3,j)=phi(-1,j) $ -(cosi/sini)*(phi(-2,j+1)-phi(-2,j-1))*dz/dx phi(-4,j)=phi(-2,j) $ -(cosi/sini)*(phi(-3,j+1)-phi(-3,j-1))*dz/dx phi(-5,j)=phi(-3,j) $ -(cosi/sini)*(phi(-4,j+1)-phi(-4,j-1))*dz/dx 176 continue do 178, i=moz+1,nz do 178, 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 178 continue ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* do 180, j=-5,0 do 180, i=-5,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=-5,nz n(i,j)=n(i,j-mox+1) 182 phi(i,j)=phi(i,j-mox+1) ********************************************************************* * * * Creates troubleshootin file named "trouble.dat" * * * ********************************************************************* OPEN(STATUS='UNKNOWN',UNIT=8,FILE='trouble.dat') CLOSE(8) WRITE(6,*)'Running Simulation' RETURN END