subroutine hyperbolic ********************************************************************** * * * This subroutine, named as fp2.f, is to solve the * * hyperbolic equation of the density n(i,j). * * * ********************************************************************** c fp2.f 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 /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 /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 416 i=-5,nz do 416 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) 416 u1x(i,j)=u1x0h*cos(xk*x+zkr*z-w*t) ********************************************************************** * * * Determine the time step 'dt'. * * * ********************************************************************** dt=0.0 552 eps1=0.0 do 501 i=2,moz-1 do 501 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) 501 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 400 i=-3,nz-2 zz=dz*real(i-11) zkr=zk*exp(-zki*zz) aven=0.0 do 408 j=1,mox-1 408 aven=aven+n(i,j) aven=aven/real(mox-1) 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*(cosi*xk-sini*zkr)*(cosi*xk-sini*zkr)*(n(i,j)-aven) * *ti(i)/(1.0+rvi(i)*rvi(i)) 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 eps1=amax1(eps1,abs(nt(i,j)-n(i,j))) 10 continue ********************************************************************** * * * Put the values of nt(i,j) to n(i,j). * * * ********************************************************************** do 20 i=3,moz-2 do 20 j=1,mox n(i,j)=nt(i,j) 20 continue ********************************************************************** * * * 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)-(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 202 continue ********************************************************************** * * * 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 ********************************************************************** * * * 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) c write(*,*)'end of hyperbolic' 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 /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 ********************************************************************** **********************************************************************