PROGRAM FP ********************************************************************* * * * This program, named as FP including subroutines * * init, ellipse, hyperbolic, output, and final calculates the * * spatial and temporal evolution of ionospheric F region * * perturbations generated by gravity waves. * * * * Version 0.3 - Last Modified 3/23/04 * * * ********************************************************************* PARAMETER (NZ=126,NX=106) PARAMETER (PI=3.141592654) REAL N,N0,ZZ,XX REAL ETIME, ELAPSED(2) COMMON /INTS/II,JJ,KK common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) 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/n0(-5:nz) common /wx10/produ(-5:nz),recom(-5:nz) CALL INIT T=0.0 II=0 JJ=0 KK=1 DO WHILE (T.LE.TT) CALL HYPERBOLIC CALL ELLIPSE IF (EPC.GT.100.0) THEN WRITE(6,*) 'Divergent Calculation - Program Aborted!' STOP END IF CALL OUTPUT END DO CALL FINAL END ******************************************************************** ******************************************************************** 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,n0 real junk,o,o2,n2,rho,temp common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) 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/n0(-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,*) WRITE(6,*)'Enter vertical wavelength in km' READ(5,*)LZ LZ=LZ*1000.0 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=pi*2.0/2400.0 xk=-pi*2.0/lx zk=-pi*2.0/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='msis.dat') READ(12,*) CJUNK DO 150, I=1,MOZ READ(12,*)JUNK,O,N2,O2,RHO,TEMP,JUNK N0(I)=O+N2+O2 150 continue CLOSE(12) do 155, i=1,moz do 155, j=-5,nx 155 n(i,j)=n0(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. * * * ********************************************************************* 170 write(6,*)'Running Simulation' 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) c open (status='old',access='append',unit=8,file='s1') c write(8,*)'end of init' c close(8) return end ********************************************************************** ********************************************************************** subroutine ellipse ********************************************************************* * * * This subroutine is to solve the ellipse * * equation of the potential phi(i,j). * * * ********************************************************************* parameter (nz=126,nx=106) parameter (pi=3.141592654) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) 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 /wx10/produ(-5:nz),recom(-5:nz) common /wx21/cnx(-5:nz),cnz(-5:nz),dn(-5:nz) common /wx22/fx(-5:nz,-5:nx),fz(-5:nz,-5:nx) common /wx23/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) common /wx35/u1x0h(-5:nz),u1z0h(-5:nz),zkr(-5:nz) dimension phi2(0:nx) ********************************************************************* * * * Give basic parameters. * * * ********************************************************************* re=2.0e-8 nit=0 dx1=1.0/(dx*dx) dz1=1.0/(dz*dz) c b1=-0.5/(dx1+dz1) it=1 a2=100.0 r=1.0 eps0=1.0 epc=0.0 epmax=10000.0 ********************************************************************* * * * Limit the amplitude of the potential phi(i,j). * * * ********************************************************************* phimax=15.0 phimin=-15.0 do 200, i=1,moz do 200, j=-5,nx if (phi(i,j).lt.phimin) phi(i,j)=phimin 200 if (phi(i,j).gt.phimax) phi(i,j)=phimax ********************************************************************* * * * Calculate the amplitude of the gravity wave. * * The vertical wavelength and amplitude of the gravity wave * * vary with height. * * * ********************************************************************* do 210, i=-5,nz do 210, j=-5,nx c u1x0h=-u1z0h*zkr/xk z=dz*real(i-1)-z0/2.0 x=dx*real(j-1)-x0/2.0 u1z(i,j)=u1z0h(i)*cos(xk*x+zkr(i)*z-w*t) 210 u1x(i,j)=u1x0h(i)*cos(xk*x+zkr(i)*z-w*t) ********************************************************************* * * * Calculate the values of relavant quantities. * * * ********************************************************************* do 220, i=-5,nz cnx(i)=-2.0*(cosi**2)*ti(i)/(1.0+rvi(i)**2) cnz(i)=-2.0*(sini**2)*ti(i)/(1.0+rvi(i)**2) 220 dn(i)=2.0*ri2(i)*ti(i) do 222, i=-5,nz do 222, j=-5,nx fx(i,j)=rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -rex1(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) $ +re2(i)*(uoz+u1z(i,j)-eoz/(rve(i)*b)) $ -eoy*sini/((1.0+rvi(i)**2)*b) $ +eoy*sini/((1.0+rve(i)**2)*b) $ +sini*cosi*gg/((1.0+rvi(i)**2)*vin(i)) fz(i,j)=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ +re2(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) $ -rez1(i)*(uoz+u1z(i,j)-eoz/(rve(i)*b)) $ -eoy*cosi/((1.0+rvi(i)**2)*b) $ +eoy*cosi/((1.0+rve(i)**2)*b) $ -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i)) 222 continue ********************************************************************* * * * Determine the values of phi(i,j) at the lower boundary. * * * ********************************************************************* 230 do 232, j=1,mox 232 phi2(j)=phi(2,j) do 234, j=1,mox 234 phi(1,j)=phi2(j)*r+(1.0-r)*phi(1,j) phi(1,mox+1)=phi(1,2) phi(1,0)=phi(1,mox-1) ********************************************************************* * * * Calculate phi(i,j) in the inner region (i=3:moz-2, j=1:mox) * * * ********************************************************************* eps1=0.0 do 265, i=3,moz-2 aven=0.0 do 240, j=1,mox-1 240 aven=aven+n(i,j) aven=aven/real(mox-1) do 245, j=1,mox aphix=n(i,j)*(rix1(i)/(rvi(i)*b)+rex1(i)/(rve(i)*b)) aphiz=n(i,j)*(riz1(i)/(rvi(i)*b)+rez1(i)/(rve(i)*b)) bphi=-n(i,j)*(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) ax=(rix1(i)/(rvi(i)*b)+rex1(i)/(rve(i)*b)) $ *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphix) $ -(n(i+1,j)*ri2(i+1)/(rvi(i+1)*b) $ +n(i+1,j)*re2(i+1)/(rve(i+1)*b) $ -n(i-1,j)*ri2(i-1)/(rvi(i-1)*b) $ -n(i-1,j)*re2(i-1)/(rve(i-1)*b))/(2.0*dz*aphix) az=(n(i+1,j)*riz1(i+1)/(rvi(i+1)*b) $ +n(i+1,j)*rez1(i+1)/(rve(i+1)*b) $ -n(i-1,j)*riz1(i-1)/(rvi(i-1)*b) $ -n(i-1,j)*rez1(i-1)/(rve(i-1)*b))/(2.0*dz*aphix) $ -(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) $ *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphix) s14=2.0*ti(i)*(cosi*xk-sini*zkr(i))*(cosi*xk-sini*zkr(i)) s14=s14/((1.0+rvi(i)*rvi(i))*aphix) s14=s14*(n(i,j)-aven)/aven s5=(n(i,j+1)*fx(i,j+1)-n(i,j-1)*fx(i,j-1))/(2.0*dx*aphix) s6=(n(i+1,j)*fz(i+1,j)-n(i-1,j)*fz(i-1,j))/(2.0*dz*aphix) sij=s14+s5+s6 dzz1=dz1*aphiz/aphix b1=-0.5/(dx1+dzz1) abij=bphi/(2.0*dx*dz*aphix) a1ij=dx1+ax/(2.0*dx) c1ij=dx1-ax/(2.0*dx) d1ij=dzz1+az/(2.0*dz) e1ij=dzz1-az/(2.0*dz) phi2(j)=b1*(sij-abij*phi(i+1,j+1)+abij*phi(i+1,j-1) $ +abij*phi(i-1,j+1)-abij*phi(i-1,j-1) $ -a1ij*phi(i,j+1)-c1ij*phi(i,j-1) $ -d1ij*phi(i+1,j)-e1ij*phi(i-1,j)) ********************************************************************* * * * Restrict the sudden change of phi2(i,j). * * * ********************************************************************* phiap=5.0*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) phian=-5.0*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) if (t.gt.20.0) then if (phi2(j).gt.phiap) then phi2(j)=0.25*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) end if if (phi2(j).lt.phian) then phi2(j)=0.25*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) end if end if 245 continue ********************************************************************* * * * Limit the amplitude of the potential phi2(j). * * * ********************************************************************* do 250, j=1,mox if (phi2(j).lt.phimin) then phi2(j)=phimin end if if (phi2(j).gt.phimax) then phi2(j)=phimax end if 250 continue ********************************************************************* * * * Put the values of phi2(j) to phi(i,j). * * * ********************************************************************* do 260, j=1,mox eps1=max(eps1,abs((phi2(j)-phi(i,j))*r)) 260 phi(i,j)=phi2(j)*r+(1.0-r)*phi(i,j) phi(i,mox+1)=phi(i,2) phi(i,0)=phi(i,mox-1) if (eps1.gt.epc) then epc=eps1 end if 265 continue ********************************************************************* * * * If 'epc' is bigger than 100, stop the calculation. * * If the calculation becomes divergent, that is, 'eps1' * * begins to increase, stop the calculation. * * * ********************************************************************* nit=nit+1 if (epc.gt.100.0) go to 280 if (eps1.lt.epmax) then epmax=eps1 else go to 280 end if ********************************************************************* * * * Determine the values of phi(2,j), phi(1,j), phi(moz-1,j), * * and phi(moz,j) for j=1,mox under suitable boundary * * conditions. These values are not calculated in the * * calculations for the inner region. * * * ********************************************************************* do 270, j=1,mox phi(2,j)=phi(4,j)-(cosi/sini)*(phi(3,j+1)-phi(3,j-1))*dz/dx phi(1,j)=phi(3,j)-(cosi/sini)*(phi(2,j+1)-phi(2,j-1))*dz/dx phi(moz-1,j)=phi(moz-3,j) $ +(cosi/sini)*(phi(moz-2,j+1)-phi(moz-2,j-1))*dz/dx phi(moz,j)=phi(moz-2,j) $ +(cosi/sini)*(phi(moz-1,j+1)-phi(moz-1,j-1))*dz/dx 270 continue phi(moz,mox+1)=phi(moz,2) phi(moz,0)=phi(moz,mox-1) ********************************************************************* * * * Calculate the relaxation factor 'r'. * * In fact, 'r' is taken to be one in the present calculations * * * ********************************************************************* c if (it.eq.1) then c a1=a2 c a2=alog10(eps1/eps0) c if (abs(a2-a1).lt.0.006) then c it=2 c if (eps1/eps0.ge.1.0) then c r=r-0.1 c else c r=2.0/(1.0+sqrt(1.0-eps1/eps0)) c end if c end if c end if c if (r.gt.1.5) r=1.5 c if (r.lt.1.0) r=1.0 r=1.0 ********************************************************************* * * * Determine the accuracy of the calculated data. * * If the absolute error 'abs(esp1-eps0)' or the relative * * error 'abs(1.0-eps0/eps1)' is smaller than the value * * required, stop calculation. * * If iteration number is more than 50, stop calculation. * * Otherwise, return to the beginning and calculate again. * * * ********************************************************************* c if (abs(eps1-eps0).lt.re) then if (abs(1.0-eps0/eps1).lt.re) go to 280 if (nit.gt.50) go to 280 c write(8,*)'eps1= ',eps1 c write(8,*)'r= ',r if (eps1.gt.re) then eps0=eps1 go to 230 end if ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* 280 do 282, j=-5,0 do 282, i=1,moz 282 phi(i,j)=phi(i,mox+j-1) do 284, j=mox+1,nx do 284, i=1,moz 284 phi(i,j)=phi(i,j-mox+1) c open (status='old',access='append',unit=8,file='s1') c write(8,*)'epc= ',epc c write(8,*)'nit= ',nit c write(8,*)'end of ellipse' c close(8) return end ********************************************************************** ********************************************************************** subroutine hyperbolic ********************************************************************** * * * This subroutine is to solve the * * hyperbolic equation of the density n(i,j). * * * ********************************************************************** parameter (nz=126,nx=106) parameter (pi=3.141592654) real n,nt common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) 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 /wx10/produ(-5:nz),recom(-5:nz) common /wx32/dntd(-5:nz,-5:nx) common /wx33/dxz common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) common /wx35/u1x0h(-5:nz),u1z0h(-5:nz),zkr(-5:nz) dimension nt(-5:nz,-5:nx) re=0.05 eps0=1.0 epmax=1.0e18 dxz=dx*dz ********************************************************************* * * * Calculate the amplitude of the gravity wave. * * The vertical wavelength and amplitude of the gravity wave * * vary with height. * * * ********************************************************************* do 300, i=-5,nz do 300, j=-5,nx c u1x0h=-u1z0h*zkr/xk z=dz*real(i-1)-z0/2.0 x=dx*real(j-1)-x0/2.0 u1z(i,j)=u1z0h(i)*cos(xk*x+zkr(i)*z-w*t) 300 u1x(i,j)=u1x0h(i)*cos(xk*x+zkr(i)*z-w*t) ********************************************************************** * * * Determine the time step 'dt'. * * * ********************************************************************** dt=0.0 310 eps1=0.0 do 312, i=2,moz-1 do 312, j=1,mox vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*sini/((1.0+rvi(i)**2)*b) $ +sini*cosi*gg/((1.0+rvi(i)**2)*vin(i)) $ -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) vix1=abs(vix1/dx) viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*cosi/((1.0+rvi(i)**2)*b) $ -(rvi(i)*rvi(i)+sini**2)*gg/((1.0+rvi(i)**2)*vin(i)) $ +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) viz1=abs(viz1/dz) 312 dt=max(dt,vix1,viz1) dt=0.35/dt if (dt.gt.10.0) then dt=10.0 end if if (dt.lt.1.0) then dt=1.0 end if ********************************************************************** * * * Calculate the time advanced solution for lower order flux. * * * ********************************************************************** do 322, i=-3,nz-2 aven=0.0 do 320, j=1,mox-1 320 aven=aven+n(i,j) aven=aven/real(mox-1) do 322, j=-3,nx-2 dntd(i,j)=n(i,j)+dt*(produ(i)-recom(i)*n(i,j)) $ -(fl(i,j)-fl(i,j-1)+gl(i,j)-gl(i-1,j))/dxz $ -dt*((cosi*xk-sini*zkr(i))**2)*(n(i,j)-aven) $ *ti(i)/(1.0+rvi(i)**2) 322 continue ********************************************************************** * * * Calculate the time advanced solution for higher order flux * * in the inner region (i=3:moz-2, j=1:mox). * * * ********************************************************************** do 330, i=3,moz-2 do 330, j=1,mox nt(i,j)=dntd(i,j)-(c1(i,j)*af(i,j)-c1(i,j-1)*af(i,j-1) $ +c2(i,j)*ag(i,j)-c2(i-1,j)*ag(i-1,j))/dxz 330 eps1=max(eps1,abs(nt(i,j)-n(i,j))) ********************************************************************** * * * Put the values of nt(i,j) to n(i,j). * * * ********************************************************************** do 340, i=3,moz-2 do 340, j=1,mox 340 n(i,j)=nt(i,j) ********************************************************************** * * * Calculate the values of n(2,j), n(1,j), n(moz-1,j), and * * n(moz,j) for j=1,mox under suitable boundary conditions. * * These values are not calculated in the calculations for * * the inner region. * * * ********************************************************************** do 350, j=1,mox n(2,j)=n(4,j)-(cosi/sini)*(n(3,j+1)-n(3,j-1))*dz/dx n(1,j)=n(3,j)-(cosi/sini)*(n(2,j+1)-n(2,j-1))*dz/dx n(moz-1,j)=n(moz-3,j) $ +(cosi/sini)*(n(moz-2,j+1)-n(moz-2,j-1))*dz/dx n(moz,j)=n(moz-2,j) $ +(cosi/sini)*(n(moz-1,j+1)-n(moz-1,j-1))*dz/dx 350 continue ********************************************************************** * * * Restrict sudden change of n(i,j) caused by numerical error. * * * ********************************************************************** do 360, i=1,moz do 360, j=1,mox if (n(i,j).lt.0.0) then n(i,j)=0.25*(n(i+1,j)+n(i-1,j)+n(i,j+1)+n(i,j-1)) end if 360 continue ********************************************************************** * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************** 370 do 372, j=-5,0 do 372, i=1,moz 372 n(i,j)=n(i,mox+j-1) do 374, j=mox+1,nx do 374, i=1,moz 374 n(i,j)=n(i,j-mox+1) c open (status='old',access='append',unit=8,file='s1') c write(8,*)'end of hyperbolic' c close(8) return end ********************************************************************** * * * The functions used above. * * * ********************************************************************** function c1(i,j) if (af(i,j).ge.0.0) then c1=min(r1(i,j+1),r0(i,j)) else c1=min(r1(i,j),r0(i,j+1)) end if return end ********************************************************************** function c2(i,j) if (ag(i,j).ge.0.0) then c2=min(r1(i+1,j),r0(i,j)) else c2=min(r1(i,j),r0(i+1,j)) end if return end ********************************************************************** function r1(i,j) parameter (nz=126,nx=106) common /wx32/dntd(-5:nz,-5:nx) common /wx33/dxz p(i,j)=max(0.0,ag(i-1,j))-min(0.0,ag(i,j))+ $ max(0.0,af(i,j-1))-min(0.0,af(i,j)) q(i,j)=(wmax(i,j)-dntd(i,j))*dxz if (p(i,j).gt.0.0) then r1=min(1.0,q(i,j)/p(i,j)) else r1=0.0 end if return end ********************************************************************** function r0(i,j) parameter (nz=126,nx=106) common /wx32/dntd(-5:nz,-5:nx) common /wx33/dxz p(i,j)=max(0.0,ag(i,j))-min(0.0,ag(i-1,j))+ $ max(0.0,af(i,j))-min(0.0,af(i,j-1)) q(i,j)=(dntd(i,j)-wmin(i,j))*dxz if (p(i,j).gt.0.0) then r0=min(1.0,q(i,j)/p(i,j)) else r0=0.0 end if return end ********************************************************************** function f(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx common /wx4/u1x0,u1z0,xk,zk,w,zki,uki 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 /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*sini/((1.0+rvi(i)*rvi(i))*b) $ +sini*cosi*gg/((1.0+rvi(i)*rvi(i))*vin(i)) $ -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) f=vix1*n(i,j) return end ********************************************************************** function g(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /wx3/b,gg,tem,eox,eoy,eoz,uox,uoz,uebx common /wx4/u1x0,u1z0,xk,zk,w,zki,uki 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 /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*cosi/((1.0+rvi(i)*rvi(i))*b) $ -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i)) $ +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) g=viz1*n(i,j) return end ********************************************************************** function fl(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /times/dt,t,tt fl=(0.5*dt*(f(i,j+1)+f(i,j))-0.25*dx*(n(i,j+1)-n(i,j)))*dz return end ********************************************************************** function gl(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /times/dt,t,tt gl=(0.5*dt*(g(i+1,j)+g(i,j))-0.25*dz*(n(i+1,j)-n(i,j)))*dx return end ********************************************************************** function fh(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /times/dt,t,tt fh=(7.0*(f(i,j+1)+f(i,j))-(f(i,j+2)+f(i,j-1)) $ )*dt*dz/12.0 return end ********************************************************************** function gh(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /times/dt,t,tt gh=(7.0*(g(i+1,j)+g(i,j))-(g(i+2,j)+g(i-1,j)) $ )*dt*dx/12.0 return end ********************************************************************** function ag(i,j) parameter (nz=126,nx=106) common /wx32/dntd(-5:nz,-5:nx) ag=gh(i,j)-gl(i,j) t1=ag*(dntd(i+1,j)-dntd(i,j)) t2=ag*(dntd(i+2,j)-dntd(i+1,j)) t3=ag*(dntd(i,j)-dntd(i-1,j)) if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) ag=0.0 return end ********************************************************************** function af(i,j) parameter (nz=126,nx=106) common /wx32/dntd(-5:nz,-5:nx) af=fh(i,j)-fl(i,j) t1=af*(dntd(i,j+1)-dntd(i,j)) t2=af*(dntd(i,j+2)-dntd(i,j+1)) t3=af*(dntd(i,j)-dntd(i,j-1)) if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) af=0.0 return end ********************************************************************** function wmax(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx32/dntd(-5:nz,-5:nx) wa(i,j)=max(n(i,j),dntd(i,j)) wmax=max(wa(i-1,j),wa(i,j),wa(i,j-1),wa(i,j+1), $ wa(i+1,j)) return end ********************************************************************** function wmin(i,j) parameter (nz=126,nx=106) real n common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx32/dntd(-5:nz,-5:nx) wb(i,j)=min(n(i,j),dntd(i,j)) wmin=min(wb(i-1,j),wb(i,j),wb(i+1,j),wb(i,j-1), $ wb(i,j+1)) return end ********************************************************************** ********************************************************************** SUBROUTINE OUTPUT ********************************************************************** * * * This subroutine writes the data for density and potential * * for 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) PARAMETER (PI=3.141592654) REAL N,N0,ZZ,XX REAL ETIME, ELAPSED(2) CHARACTER OUTFILE*9, ALPHA*10 REAL AVERN(MOZ), DELTAN(MOZ,MOX) 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/n0(-5:nz) common /wx10/produ(-5:nz),recom(-5:nz) ******************************************************************** * * * 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 400 I=1,MOZ 400 AVERN(I)=0.0 DO 404 I=1,MOZ DO 402 J=1,MOX-1 402 AVERN(I)=AVERN(I)+N(I,J) AVERN(I)=AVERN(I)/REAL(MOX-1) DO 404 J=1,MOX 404 DELTAN(I,J)=(N(I,J)-AVERN(I))/AVERN(I) ********************************************************************* * * * 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-3)","Perturbations","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) T=T+DT RETURN END *************************************************************************** *************************************************************************** SUBROUTINE FINAL ********************************************************************* * * * Write the data into 'in1.dat' * * The data are written in the same form as that by which * * n(i,j) and phi1(i,j) are generated. * * * ********************************************************************* c INSERT('var.inc') PARAMETER (NZ=126,NX=106) PARAMETER (PI=3.141592654) 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 /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/n0(-5:nz) common /wx10/produ(-5:nz),recom(-5:nz) OPEN(STATUS='UNKNOWN',UNIT=9,FILE='in1.dat') WRITE(9,*)T DO 500, I=1,MOZ DO 500, J=1,MOX 500 WRITE(9,*)N(I,J),PHI(I,J) CLOSE(9) WRITE(6,*) WRITE(6,*)'Program is finished' WRITE(6,*) WRITE(6,*)'Elapsed Time is ',ETIME(ELAPSED),' seconds' WRITE(6,*) RETURN END