program FP ********************************************************************* * * * This program, named as FP including subroutines * * subroutines init, ellipse, and hyperbolic, calculates * * the spatial and temporal evolution of ionospheric * * F region perturbations generated by gravity waves. * * * ********************************************************************* parameter (nz=126,nx=106) real n,n0,deltan,zz,xx character*9 outfile 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),deltan(-5:nz,-5:nx) common /wx10/produ(-5:nz),recom(-5:nz) common /wx12/avern(1:121) call init t=0.0 ii=0 jj=0 do while (t.le.tt) call hyperbolic call ellipse if (epc.gt.100.0) go to 999 ********************************************************************* * * * Calculate average density and relative perturbations. * * * ********************************************************************* c do 10 i=1,moz c 10 avern(i)=0.0 c do 12 i=1,moz c do 12 j=1,mox-1 c 12 avern(i)=avern(i)+n(i,j) c avern(i)=avern(i)/real(mox-1) c do 14 j=1,mox c 14 deltan(i,j)=(n(i,j)-avern(i))/avern(i) ******************************************************************** * * * Write the data into files for desired times. * * * ******************************************************************** if (ii.eq.10) then jj=jj+1 ii=0 else ii=ii+1 go to 30 end if open(status='unknown',unit=99,file='junk') write(99,*)jj,'out.dat' backspace(99) read(99,*) outfile close(99) open(status='unknown',unit=9,file=outfile) write(9,*) 'TITLE="fp.f ~ t =',t,'"' write(9,*) 'VARIABLES="Latitude","Altitude","Density"' write(9,*) 'ZONE F=POINT I=',mox,',J=101' do 20 j=1,mox do 20 i=11,111 zz=160.0+dz*real(i-1)/1000. xx=(dx*real(j-1)-lx)/1000. 20 write(9,*)xx,zz,n(i,j) close(9) ********************************************************************* * * * Write down the time step and the total time. * * * ********************************************************************* 30 open (status='old',access='append',unit=8,file='s1') t=t+dt write(8,*)'epc=',epc write(8,*)'dt=',dt write(8,*)'t=',t close(8) end do ********************************************************************* * * * 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='in1.dat') write(9,*)t do 60 i=1,moz do 60 j=1,mox 60 write(9,*)n(i,j),phi(i,j) close(9) write(6,*)'Program is finished' 999 stop end ******************************************************************** ******************************************************************** subroutine init ********************************************************************* * * * This subroutine gives the initial values and background * * parameters. * * * ********************************************************************* parameter (nz=126,nx=106) integer input real n,n0,deltan 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),deltan(-5:nz,-5:nx) common /wx10/produ(-5:nz),recom(-5:nz) *********************************************************************** * * * Initiates software, and asks several user input options * * * *********************************************************************** open (status='unknown',unit=8,file='s1') close(8) 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(45.0*3.14159/180.0) cosi=cos(45.0*3.14159/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=3.14159265*2.0/2400.0 xk=-3.14159265*2.0/lx zk=-3.14159265*2.0/lz zki=5.2e-6 uki=2.75e-6 100 write(6,*)'Input Method:' write(6,*)' (1) - New Calculations' write(6,*)' (2) - Old Calculations' read(5,*) input if (input.eq.1) then go to 110 else if (input.eq.2) then go to 160 else write(6,*)'Not a valid option' write(6,*) go to 100 end if end if ********************************************************************* * * * Give height profiles of production rate 'produ(i)' and * * recombination coefficient 'recom(i)'. * * * ********************************************************************* 110 write(6,*)'Day or Night?' write(6,*)' (1) - Day' write(6,*)' (2) - Night' read(5,*) input if (input.eq.1) then produ0=300000000.0 else if (input.eq.2) then produ0=3.0e-16 else write(6,*)'Not a valid option' write(6,*) go to 110 end if end if 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 if (produ(i).gt.produ0) produ(i)=produ0 120 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 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)*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)) 140 te(i)=1.51567e7*tem/(me*ven(i)) ********************************************************************* * * * Read the background density profile. * * * ********************************************************************* open(status='unknown',unit=12,file='midn2a.dat') do 150 i=1,moz 150 read(12,*)n0(i) close(12) do 155 i=1,moz do 155 j=-5,nx 155 n(i,j)=n0(i) go to 170 ********************************************************************* * * * 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 open(status='unknown',unit=9,file='in1.dat') read(9,*)t do 162 i=1,moz do 162 j=1,mox 162 read(9,*)n(i,j),phi(i,j) close(9) ********************************************************************* * * * 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) open (status='old',access='append',unit=8,file='s1') write(8,*)'end of init' 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 /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) 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 zz=dz*real(i-11) zkr=zk*exp(-zki*zz) u1z0h=u1z0*exp(uki*zz) u1x0h=u1x0*exp(uki*zz) 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*cos(xk*x+zkr*z-w*t) 210 u1x(i,j)=u1x0h*cos(xk*x+zkr*z-w*t) ********************************************************************* * * * Calculate the values of relavant quantities. * * * ********************************************************************* do 220 i=-5,nz cnx(i)=-2.0*cosi*cosi*ti(i)/(1.0+rvi(i)*rvi(i)) cnz(i)=-2.0*sini*sini*ti(i)/(1.0+rvi(i)*rvi(i)) 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)*rvi(i))*b) * +eoy*sini/((1.0+rve(i)*rve(i))*b) * +sini*cosi*gg/((1.0+rvi(i)*rvi(i))*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)*rvi(i))*b) * +eoy*cosi/((1.0+rve(i)*rve(i))*b) * -(rvi(i)*rvi(i)+sini*sini)*gg/((1.0+rvi(i)*rvi(i))*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 zz=dz*real(i-11) zkr=zk*exp(-zki*zz) 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)*(cosi*xk-sini*zkr) 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=amax1(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 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 c 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 open (status='old',access='append',unit=8,file='s1') write(8,*)'eps1= ',eps1 write(8,*)'r= ',r close(8) 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) open (status='old',access='append',unit=8,file='s1') write(8,*)'epc= ',epc write(8,*)'nit= ',nit write(8,*)'end of ellipse' 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) 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) 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 wavelength and amplitude of the gravity wave * * vary with height. * * * ********************************************************************** do 300 i=-5,nz do 300 j=-5,nx zz=dz*real(i-11) zkr=zk*exp(-zki*zz) u1z0h=u1z0*exp(uki*zz) u1x0h=u1x0*exp(uki*zz) 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*cos(xk*x+zkr*z-w*t) 300 u1x(i,j)=u1x0h*cos(xk*x+zkr*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)*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) 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)*rvi(i))*b) * -(rvi(i)*rvi(i)+sini*sini)*gg/((1.0+rvi(i)*rvi(i))*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=amax1(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 zz=dz*real(i-11) zkr=zk*exp(-zki*zz) 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)*(cosi*xk-sini*zkr)*(n(i,j)-aven) * *ti(i)/(1.0+rvi(i)*rvi(i)) 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=amax1(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) open (status='old',access='append',unit=8,file='s1') write(8,*)'end of hyperbolic' 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,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)*rvi(i)+sini*sini)*gg/((1.0+rvi(i)*rvi(i))*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)=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 ********************************************************************** **********************************************************************