program FP ********************************************************************* * * * This program, named as FP (fp.f), including subroutines * * subroutines init, ellipse (fp1.f), and hyperbolic (fp2.f) * * calcultes the spatial and temporal evolution of ionospheric * * F region perturbations generated by gravity waves. * * * ********************************************************************* c pf.f parameter (nz=126,nx=106) parameter (pi=3.141592654) real n,n0,deltan REAL ETIME, ELAPSED(2) common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /wx2/dt,t,tt 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 /wx8/epc,dperp common /wx9/n0(-5:nz),deltan(-5:nz,-5:nx) common /wx10/produ(-5:nz),recom(-5:nz) common /wx12/avern(1:121) call init t=0.0 tt=1201.0 do while (t.lt.tt) call hyperbolic call ellipse if (epc.gt.100.0) then go to 118 end if ********************************************************************* * * * Calculate average density and relative perturbations. * * * ********************************************************************* do 618 i=1,moz 618 avern(i)=0.0 do 620 i=1,moz do 622 j=1,mox-1 622 avern(i)=avern(i)+n(i,j) avern(i)=avern(i)/real(mox-1) do 620 j=1,mox deltan(i,j)=(n(i,j)-avern(i))/avern(i) 620 continue ********************************************************************* * * * Write the data into two files: 'inn1.dat' and 'inp1.dat'. * * 'inn1.dat' saves the density n(i,j), * * 'inp1.dat' saves the potential phi1(i,j). * * The data are written in the same form as that by which * * n(i,j) and phi1(i,j) are generated. * * * ********************************************************************* open(status='unknown',unit=9,file='inn1.dat') do 60 i=1,moz do 60 j=1,mox 60 write(9,*)n(i,j) close(9) open(status='unknown',unit=9,file='inp1.dat') do 65 i=1,moz do 65 j=1,mox 65 write(9,*)phi(i,j) close(9) ********************************************************************* * * * Write the data into three files prepared for plotting. * * 'outn1' saves the density n(i,j), * * 'outd1' saves the relative density perturbation * * deltan(i,j), * * 'outp1' saves the potential phi(i,j). * * The data are written in the form that the plotting * * software GRAFTOOL requires. * * Note that the number in FORMAT must be the range for 'i'. * * * ********************************************************************* open(status='unknown',unit=9,file='outn1') write(9,200)0.0,(i,i=11,111) 200 format(f3.1,101i12) do 7 j=1,mox 7 write(9,210)j,(n(i,j),i=11,111) 210 format(i3,101e12.4) close(9) open(status='unknown',unit=9,file='outd1') write(9,200)0.0,(i,i=11,111) do 76 j=1,mox 76 write(9,210)j,(deltan(i,j),i=11,111) close(9) open(status='unknown',unit=9,file='outa1') do 762 i=1,moz 762 write(9,220)i,avern(i) 220 format(i3,e12.4) close(9) open(status='unknown',unit=9,file='outp1') write(9,200)0.0,(i,i=11,111) do 8 j=1,mox 8 write(9,210)j,(phi(i,j),i=11,111) close(9) ******************************************************************** * * * Write the data into files for desired times. * * * ******************************************************************** if (t.gt.995.0.and.t.lt.1005.0) then open(status='unknown',unit=9,file='out1.dat') WRITE(9,*) 'TITLE="fp.f ~ t =',T,' s"' WRITE(9,*) 'VARIABLES="Latitude","Altitude","Density", $ "Perturbations","Potential"' WRITE(9,*) 'ZONE F=POINT I=',mox,',J=101' DO 711, J=1,MOX DO 711, I=11,111 XX=(DX*REAL(J-1)-LX)/1000. ZZ=160.0+DZ*REAL(I-1)/1000. 711 WRITE(9,*)XX,ZZ,N(I,J),DELTAN(I,J),PHI(I,J) CLOSE(9) end if if (t.gt.1995.0.and.t.lt.2005.0) then open(status='unknown',unit=9,file='out2.dat') write(9,*)0.0,(i,i=11,111) do 712 j=1,mox 712 write(9,*)j,(n(i,j),i=11,111) close(9) end if if (t.gt.2995.0.and.t.lt.3005.0) then open(status='unknown',unit=9,file='out3.dat') write(9,*)0.0,(i,i=11,111) do 713 j=1,mox 713 write(9,*)j,(n(i,j),i=11,111) close(9) end if if (t.gt.3995.0.and.t.lt.4005.0) then open(status='unknown',unit=9,file='out4.dat') write(9,*)0.0,(i,i=11,111) do 714 j=1,mox 714 write(9,*)j,(n(i,j),i=11,111) close(9) end if ********************************************************************* * * * Write down the time step and the total time. * * * ********************************************************************* t=t+dt write(8,*)'epc=',epc write(8,*)'dt=' write(8,*)dt write(8,*)'t=' write(8,*)t close(8) open (status='old',access='append',unit=8,file='s1') end do WRITE(6,*)'Elapsed Time is ',ETIME(ELAPSED),' seconds' 118 stop end ******************************************************************** ******************************************************************** subroutine init ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * ********************************************************************* parameter (nz=126,nx=106) parameter (pi=3.141592654) real n,n0,deltan common n(-5:nz,-5:nx),phi(-5:nz,-5:nx) common /wx1/dz,dx,z0,x0,moz,mox common /wx2/dt,t,tt 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 /wx8/epc,dperp common /wx9/n0(-5:nz),deltan(-5:nz,-5:nx) common /wx10/produ(-5:nz),recom(-5:nz) open (status='unknown',unit=8,file='s1') moz=121 mox=101 z0=480000.0 x0=600000.0 dz=z0/real(moz-1) dx=x0/real(mox-1) b=4.5e-5 gg=0.000001 sini=sin(45.0*pi/180.0) cosi=cos(45.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/300000.0 zk=-pi*2.0/100000.0 zki=5.2e-6 uki=2.75e-6 ********************************************************************* * * * Give height profiles of production rate 'produ(i)' and * * recombination coefficient 'recom(i)'. * * * ********************************************************************* do 155 i=1,moz c produ0=300000000.0 recom0=0.00015 produ0=3.0e-16 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) then recom(i)=recom0 end if if (produ(i).gt.produ0) then produ(i)=produ0 end if 155 continue ********************************************************************* * * * 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 2 i=-5,nz do 2 j=-5,nx 2 phi(i,j)=0.0 vin0=6.5 ven0=210.0 do 162 i=1,10 vin(i)=vin0*exp(0.14*real(11-i)) 162 ven(i)=ven0*exp(0.14*real(11-i)) do 164 i=11,moz vin(i)=vin0*exp(-0.05*real(i-11)) 164 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 10 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)*rvi(i)+cosi*cosi)/(1.0+rvi(i)*rvi(i)) riz1(i)=(rvi(i)*rvi(i)+sini*sini)/(1.0+rvi(i)*rvi(i)) ri2(i)=sini*cosi/(1.0+rvi(i)*rvi(i)) rex1(i)=(rve(i)*rve(i)+cosi*cosi)/(1.0+rve(i)*rve(i)) rez1(i)=(rve(i)*rve(i)+sini*sini)/(1.0+rve(i)*rve(i)) re2(i)=sini*cosi/(1.0+rve(i)*rve(i)) ti(i)=8254.81*tem/(mi*vin(i)) te(i)=1.51567e7*tem/(me*ven(i)) 10 continue ********************************************************************* * * * Read the background density profile. * * * ********************************************************************* open(status='unknown',unit=12,file='midn2a.dat') do 385 i=1,moz 385 read(12,*)n0(i) close(12) do 40 i=1,moz do 40 j=-5,nx n(i,j)=n0(i) 40 continue ********************************************************************* * * * When ki=1, this program restarts calculations with * * previously obtained data, saved in 'inn1.dat' and * * 'inp1.dat', as initial profiles. * * For this case, the time 't' must be changed to the * * value at which the program stopped previously. * * * ********************************************************************* ki=10 if (ki.eq.1) then open(status='unknown',unit=9,file='inn1.dat') do 60 i=1,moz do 60 j=1,mox 60 read(9,*)n(i,j) close(9) open(status='unknown',unit=9,file='inp1.dat') do 65 i=1,moz do 65 j=1,mox 65 read(9,*)phi(i,j) close(9) end if ********************************************************************* * * * Extend boundary condition in the vertical direction. * * * ********************************************************************* do 70 i=-5,0 vin(i)=vin(1) ven(i)=ven(1) produ(i)=produ(1) recom(i)=recom(1) 70 n0(i)=n0(1) do 80 i=moz+1,nz vin(i)=vin(moz) ven(i)=ven(moz) produ(i)=produ(moz) recom(i)=recom(moz) 80 n0(i)=n0(moz) do 90 i=-5,0 do 90 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 90 continue do 100 i=moz+1,nz do 100 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 100 continue ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* do 110 j=-5,0 do 110 i=-5,nz n(i,j)=n(i,mox+j-1) 110 phi(i,j)=phi(i,mox+j-1) do 120 j=mox+1,nx do 120 i=-5,nz n(i,j)=n(i,j-mox+1) 120 phi(i,j)=phi(i,j-mox+1) write(8,*)'end of init' return end ********************************************************************** ********************************************************************** include "fp1.for" include "fp2.for"