* program MIDF ********************************************************************* * * * This program, named as MIDF (midf.f), including * * subroutines init, ellipse, and hyperbolic, calculates * * the spatial and temporal evolution of midlatitude * * F region irregularities generated by gravity waves. * * * ********************************************************************* parameter (nz=126,nx=106) 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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) common /wx8/epc,dperp common /wx9/n0(-5:nz),deltan(-5:nz,-5:nx) common /wx10/produ(-5:nz),recom(-5:nz) call init t=0.0 tt=100.0 do while (t.lt.tt) call hyperbolic call ellipse if (epc.gt.100.0) go to 118 ********************************************************************* * * * Write down the time and other data needed. * * * ********************************************************************* t=t+dt open (status='old',access='append',unit=8,file='s1') write(8,*)'epc=',epc write(8,*)'dt=',dt write(8,*)'t=',t write(8,*) close(8) end do 118 stop ********************************************************************* * * * Calculate the relative density perturbations. * * * ********************************************************************* do 620 i=1,moz do 620 j=1,mox 620 deltan(i,j)=(n(i,j)-n0(i))/n0(i) ********************************************************************* * * * 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. * * The data are written in the form that the plotting * * software Tecplot requires. * * * ********************************************************************* open(status='unknown',unit=9,file='midf.dat') write(9,*) 'TITLE="Ionospheric Simulation - midf.f"' write(9,*) 'VARIABLES="Altitude","Latitude","Density", *"Perturbations","Potential"' write(9,*) 'ZONE F=POINT I=',mox,',J=101' do 7 j=1,mox do 7 i=11,111 7 write(9,*)dz*real(i-1)/1000.,j,n(i,j),deltan(i,j),phi(i,j) close(9) end ******************************************************************** ******************************************************************** subroutine init ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * ********************************************************************* parameter (nz=126,nx=106) 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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) 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 dperp=1.0 eox=1.91e-8 eoy=1.5e-8 eoz=1.91e-8 uox=-0.00000001 uebx=0.00000001 u1x0=-18.0 u1z0=6.0 w=3.14159265*2.0/2400.0 xk=-3.14159265*2.0/300000.0 zk=-3.14159265*2.0/100000.0 ********************************************************************* * * * Give height profiles of production rate 'produ(i)' and * * recombination coefficient 'recom(i)'. * * * ********************************************************************* do 155 i=1,moz produ0=333000000.0 recom0=0.0002 hh=dz*real(i-1)/1000.0 produ(i)=produ0*exp(-abs(hh-300.0)/60.0) recom(i)=recom0*exp(-2.0*(hh-300.0)/60.0) 155 if (recom(i).gt.recom0) recom(i)=recom0 ********************************************************************* * * * 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 cosi=0.7071 sini=0.7071 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.058*real(i-11)) 164 ven(i)=ven0*exp(-0.058*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=0.15 mi=18.0 me=1.0 omegai=4310.6/mi omegae=7.91421e6 rvi(i)=vin(i)/omegai rve(i)=ven(i)/omegae ri1(i)=(rvi(i)*rvi(i)+cosi*cosi)/(1.0+rvi(i)*rvi(i)) ri2(i)=sini*cosi/(1.0+rvi(i)*rvi(i)) re1(i)=(rve(i)*rve(i)+cosi*cosi)/(1.0+rve(i)*rve(i)) re2(i)=sini*cosi/(1.0+rve(i)*rve(i)) ti(i)=8254.81*tem/(mi*vin(i)) 10 te(i)=1.51567e7*tem/(me*ven(i)) ********************************************************************* * * * Read the background density profile midn2a.dat * * * ********************************************************************* 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 40 n(i,j)=n0(i) ********************************************************************* * * * 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)-(n(1,j+1)-n(1,j-1))*zk*dz/(xk*dx) n(-1,j)=n(1,j)-(n(0,j+1)-n(0,j-1))*zk*dz/(xk*dx) n(-2,j)=n(0,j)-(n(-1,j+1)-n(-1,j-1))*zk*dz/(xk*dx) n(-3,j)=n(-1,j)-(n(-2,j+1)-n(-2,j-1))*zk*dz/(xk*dx) n(-4,j)=n(-2,j)-(n(-3,j+1)-n(-3,j-1))*zk*dz/(xk*dx) n(-5,j)=n(-3,j)-(n(-4,j+1)-n(-4,j-1))*zk*dz/(xk*dx) phi(0,j)=phi(2,j)-phi(1,j+1)+phi(1,j-1) phi(-1,j)=phi(1,j)-phi(0,j+1)+phi(0,j-1) phi(-2,j)=phi(0,j)-phi(-1,j+1)+phi(-1,j-1) phi(-3,j)=phi(-1,j)-phi(-2,j+1)+phi(-2,j-1) phi(-4,j)=phi(-2,j)-phi(-3,j+1)+phi(-3,j-1) 90 phi(-5,j)=phi(-3,j)-phi(-4,j+1)+phi(-4,j-1) do 100 i=moz+1,nz do 100 j=1,mox n(i,j)=n(i-2,j)+(n(i-1,j+1)-n(i-1,j-1))*zk*dz/(xk*dx) 100 phi(i,j)=phi(i-2,j)+phi(i-1,j+1)-phi(i-1,j-1) ********************************************************************* * * * 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' write(8,*) 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) real n 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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) common /wx8/epc,dperp common /wx10/produ(-5:nz),recom(-5:nz) common /wx21/cn(-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) dimension phi2(0:nx) open (status='old',access='append',unit=8,file='s1') write(8,*)'beginning of ellipse' write(8,*) close(8) ********************************************************************* * * * Give basic parameters. * * * ********************************************************************* re=2.0e-8 nit=0 dx1=1.0/(dx*dx) dz1=1.0/(dz*dz) 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 520 i=1,moz do 520 j=-5,nx if (phi(i,j).lt.phimin) phi(i,j)=phimin 520 if (phi(i,j).gt.phimax) phi(i,j)=phimax ********************************************************************* * * * Calculate the values of relavant quantities. * * * ********************************************************************* cosi=0.7071 sini=0.7071 do 802 i=-5,nz do 802 j=-5,nx z=dz*real(i-1)-z0/2.0 x=dx*real(j-1)-x0/2.0 u1z(i,j)=u1z0*cos(xk*x+zk*z-w*t) 802 u1x(i,j)=u1x0*cos(xk*x+zk*z-w*t) do 806 i=-5,nz cn(i)=-ri1(i)*ti(i)+re1(i)*te(i) 806 dn(i)=ri2(i)*ti(i)-re2(i)*te(i) do 808 i=-5,nz do 808 j=-5,nx fx(i,j)=ri1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * -ri2(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * -eoy*sini/((1.0+rvi(i)*rvi(i))*b) * -re1(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) * +re2(i)*(u1z(i,j)-eoz/(rve(i)*b)) * +eoy*sini/((1.0+rve(i)*rve(i))*b) fz(i,j)=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * +ri1(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * -eoy*cosi/((1.0+rvi(i)*rvi(i))*b) * +re2(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) * -re1(i)*(u1z(i,j)-eoz/(rve(i)*b)) * +eoy*cosi/((1.0+rve(i)*rve(i))*b) 808 continue ********************************************************************* * * * Determine the values of phi(i,j) at the lower boundary. * * * ********************************************************************* 5 do 10 j=1,mox 10 phi2(j)=phi(2,j) do 20 j=1,mox 20 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 40 i=3,moz-2 do 50 j=1,mox aphi=n(i,j)*(ri1(i)/(rvi(i)*b)+re1(i)/(rve(i)*b)) bphi=-n(i,j)*(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) ax=(ri1(i)/(rvi(i)*b)+re1(i)/(rve(i)*b)) * *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphi) * -(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*aphi) az=(n(i+1,j)*ri1(i+1)/(rvi(i+1)*b) * +n(i+1,j)*re1(i+1)/(rve(i+1)*b) * -n(i-1,j)*ri1(i-1)/(rvi(i-1)*b) * -n(i-1,j)*re1(i-1)/(rve(i-1)*b))/(2.0*dz*aphi) * -(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) * *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphi) abij=bphi/(2.0*dx*dz*aphi) a1ij=dx1+ax/(2.0*dx) c1ij=dx1-ax/(2.0*dx) d1ij=dz1+az/(2.0*dz) e1ij=dz1-az/(2.0*dz) s1=(cn(i)/aphi)*(dx1*(n(i,j+1)-2.0*n(i,j)+n(i,j-1)) * +dz1*(n(i+1,j)-2.0*n(i,j)+n(i-1,j))) s2=(dn(i+1)-dn(i-1))*(n(i,j+1)-n(i,j-1))/ * (4.0*dx*dz*aphi) s3=(cn(i+1)-cn(i-1))*(n(i+1,j)-n(i-1,j))/ * (4.0*dz*dz*aphi) s4=dn(i)*(n(i+1,j+1)-n(i+1,j-1)-n(i-1,j+1) * +n(i-1,j-1))/(2.0*dx*dz*aphi) s5=(n(i,j+1)*fx(i,j+1)-n(i,j-1)*fx(i,j-1))/(2.0*dx*aphi) s6=(n(i+1,j)*fz(i+1,j)-n(i-1,j)*fz(i-1,j))/(2.0*dz*aphi) sij=s1+s2+s3+s4+s5+s6 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 50 continue ********************************************************************* * * * Limit the amplitude of the potential phi2(j). * * * ********************************************************************* do 580 j=1,mox if (phi2(j).lt.phimin) phi2(j)=phimin 580 if (phi2(j).gt.phimax) phi2(j)=phimax ********************************************************************* * * * Put the values of phi2(j) to phi(i,j). * * * ********************************************************************* do 60 j=1,mox eps1=amax1(eps1,abs((phi2(j)-phi(i,j))*r)) 60 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) 40 if (eps1.gt.epc) epc=eps1 ********************************************************************* * * * 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 19 if (eps1.lt.epmax) then epmax=eps1 else go to 19 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 90 j=1,mox phi(2,j)=phi(4,j)-phi(3,j+1)+phi(3,j-1) phi(1,j)=phi(3,j)-phi(2,j+1)+phi(2,j-1) phi(moz-1,j)=phi(moz-3,j)+phi(moz-2,j+1)-phi(moz-2,j-1) 90 phi(moz,j)=phi(moz-2,j)+phi(moz-1,j+1)-phi(moz-1,j-1) 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 one in the present calculations. * * * ********************************************************************* if (it.eq.1) then a1=a2 a2=alog10(eps1/eps0) if (abs(a2-a1).lt.0.006) then it=2 if (eps1/eps0.ge.1.0) then r=r-0.1 else r=2.0/(1.0+sqrt(1.0-eps1/eps0)) end if end if end if if (r.gt.1.5) r=1.5 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 19 if (nit.gt.50) go to 19 c write(*,*)'eps1= ',eps1 c write(*,*)'r= ',r if (eps1.gt.re) then eps0=eps1 go to 5 end if ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* 19 do 120 j=-5,0 do 120 i=1,moz 120 phi(i,j)=phi(i,mox+j-1) do 130 j=mox+1,nx do 130 i=1,moz 130 phi(i,j)=phi(i,j-mox+1) open (status='old',access='append',unit=8,file='s1') write(8,*)'epc= ',epc write(8,*)'nit= ',nit write(8,*)'end of ellipse' write(8,*) close(8) return end ********************************************************************** ********************************************************************** subroutine hyperbolic ********************************************************************** * * * This subroutine, named as midf2.f, is to solve the * * hyperbolic equation of the density n(i,j). * * * ********************************************************************** parameter (nz=126,nx=106) real n,nt 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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) 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) dimension nt(-5:nz,-5:nx) re=0.05 eps0=1.0 epmax=1.0e18 nit=0 dxz=dx*dz open (status='old',access='append',unit=8,file='s1') write(8,*)'beginning of hyperbolic' write(8,*) close(8) ********************************************************************** * * * Determine the time step 'dt'. * * * ********************************************************************** do 416 i=-5,nz do 416 j=-5,nx z=dz*real(i-1)-z0/2.0 x=dx*real(j-1)-x0/2.0 u1z(i,j)=u1z0*cos(xk*x+zk*z-w*t) 416 u1x(i,j)=u1x0*cos(xk*x+zk*z-w*t) dt=0.0 552 eps1=0.0 do 501 i=2,moz-1 do 501 j=1,mox cosi=0.7071 sini=0.7071 vix=uebx+ri1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * -eoy*sini/((1.0+rvi(i)*rvi(i))*b) * -ri2(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * -ri1(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) * +ri2(i)*ti(i)*(n(i+1,j)-n(i-1,j))/(2.0*dz*n(i,j)) vix=abs(vix/dx) viz=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * -eoy*cosi/((1.0+rvi(i)*rvi(i))*b) * +ri1(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) * -ri1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) * +ri2(i)*ti(i)*(n(i,j+1)-n(i,j-1))/(2.0*dx*n(i,j)) viz=abs(viz/dz) 501 dt=amax1(dt,vix,viz) dt=0.35/dt if (dt.gt.10.0) dt=10.0 if (dt.lt.1.0) dt=1.0 ********************************************************************** * * * Calculate the time advanced solution for lower order flux. * * * ********************************************************************** do 400 i=-3,nz-2 do 400 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*ri1(i)*ti(i)*(n(i,j+1)-2.0*n(i,j)+n(i,j-1))/(dx*dx) * +dt*ri1(i)*ti(i)*(n(i+1,j)-2.0*n(i,j)+n(i-1,j))/(dz*dz) * +dt*dperp*(n(i,j+1)-2.0*n(i,j)+n(i,j-1))/(dx*dx) * +dt*dperp*(n(i+1,j)-2.0*n(i,j)+n(i-1,j))/(dz*dz) 400 continue ********************************************************************** * * * Calculate the time advanced solution for higher order flux * * in the inner region (i=3:moz-2, j=1:mox). * * * ********************************************************************** do 10 i=3,moz-2 do 10 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 10 eps1=amax1(eps1,abs(nt(i,j)-n(i,j))) ********************************************************************** * * * If the calculation becomes divergent, that is, 'eps1' * * begins to increase, stop calculation. * * * ********************************************************************** c write(*,*)'eps1=',eps1 c if (eps1.lt.epmax) then c epmax=eps1 c else c goto 192 c end if ********************************************************************** * * * Put the values of nt(i,j) to n(i,j). * * * ********************************************************************** do 20 i=3,moz-2 do 20 j=1,mox 20 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 202 j=1,mox n(2,j)=n(4,j)-(n(3,j+1)-n(3,j-1))*zk*dz/(xk*dx) n(1,j)=n(3,j)-(n(2,j+1)-n(2,j-1))*zk*dz/(xk*dx) n(moz-1,j)=n(moz-3,j)+(n(moz-2,j+1)-n(moz-2,j-1))*zk*dz/ * (xk*dx) 202 n(moz,j)=n(moz-2,j)+(n(moz-1,j+1)-n(moz-1,j-1))*zk*dz/(xk*dx) ********************************************************************** * * * Restrict sudden change of n(i,j) caused by numerical error. * * * ********************************************************************** do 206 i=1,moz do 206 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 206 continue ********************************************************************** * * * If the relative error 'abs(1.0-eps0/eps1)' is smaller than * * the value required, stop calculation. * * If iteration number is more than 4, stop calculation. * * Otherwise, return to the beginning and calculate again. * * Note that iteration in this subroutine causes terrorble * * errors, so that no iteration is used in calculations. * * * ********************************************************************** c nit=nit+1 c if (abs(1.0-eps0/eps1).lt.0.05) then c goto 192 c else c if (nit.gt.4) then c goto 192 c else c eps0=eps1 c goto 552 c end if c end if ********************************************************************** * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************** 192 do 120 j=-5,0 do 120 i=1,moz 120 n(i,j)=n(i,mox+j-1) do 130 j=mox+1,nx do 130 i=1,moz 130 n(i,j)=n(i,j-mox+1) open (status='old',access='append',unit=8,file='s1') write(8,*)'end of hyperbolic' write(8,*) close(8) return end ********************************************************************** * * * The functions used above. * * * ********************************************************************** function c1(i,j) if (af(i,j).ge.0.0) then c1=amin1(r1(i,j+1),r0(i,j)) else c1=amin1(r1(i,j),r0(i,j+1)) end if return end ********************************************************************** function c2(i,j) if (ag(i,j).ge.0.0) then c2=amin1(r1(i+1,j),r0(i,j)) else c2=amin1(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)=amax1(0.0,ag(i-1,j))-amin1(0.0,ag(i,j))+ * amax1(0.0,af(i,j-1))-amin1(0.0,af(i,j)) q(i,j)=(wmax(i,j)-dntd(i,j))*dxz if (p(i,j).gt.0.0) then r1=amin1(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)=amax1(0.0,ag(i,j))-amin1(0.0,ag(i-1,j))+ * amax1(0.0,af(i,j))-amin1(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=amin1(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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) cosi=0.7071 sini=0.7071 f=uebx+ri1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * -eoy*sini/((1.0+rvi(i)*rvi(i))*b) * -ri2(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * -ri1(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) * +ri2(i)*ti(i)*(n(i+1,j)-n(i-1,j))/(2.0*dz*n(i,j)) f=f*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,tem,eox,eoy,eoz,uox,uebx common /wx4/u1x0,u1z0,xk,zk,w common /wx5/vin(-5:nz),rvi(-5:nz),ri1(-5:nz),ri2(-5:nz) common /wx6/ven(-5:nz),rve(-5:nz),re1(-5:nz),re2(-5:nz) common /wx7/ti(-5:nz),te(-5:nz) common /wx34/u1x(-5:nz,-5:nx),u1z(-5:nz,-5:nx) cosi=0.7071 sini=0.7071 g=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) * -eoy*cosi/((1.0+rvi(i)*rvi(i))*b) * +ri1(i)*(u1z(i,j)+eoz/(rvi(i)*b)) * +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) * -ri1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) * +ri2(i)*ti(i)*(n(i,j+1)-n(i,j-1))/(2.0*dx*n(i,j)) g=g*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 /wx2/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 /wx2/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 /wx2/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 /wx2/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))) then ag=0.0 end if 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))) then af=0.0 end if 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)=amax1(n(i,j),dntd(i,j)) wmax=amax1(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)=amin1(n(i,j),dntd(i,j)) wmin=amin1(wb(i-1,j),wb(i,j),wb(i+1,j),wb(i,j-1), * wb(i,j+1)) return end ********************************************************************** **********************************************************************