subroutine ellipse ********************************************************************* * * * This subroutine, named as fp1.f, is to solve the ellips * * equation of the potential phi(i,j). * * * ********************************************************************* c fp1.f 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 /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 /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 520 i=1,moz do 520 j=-5,nx if (phi(i,j).lt.phimin) then phi(i,j)=phimin end if if (phi(i,j).gt.phimax) then phi(i,j)=phimax end if 520 continue ********************************************************************* * * * Calculate the amplitude of the gravity wave. * * The vertical wavelength and amplitude of the gravity wave * * vary with height. * * * ********************************************************************* do 802 i=-5,nz do 802 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) u1x(i,j)=u1x0h*cos(xk*x+zkr*z-w*t) 802 continue ********************************************************************* * * * Calculate the values of relavant quantities. * * * ********************************************************************* do 806 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)) dn(i)=2.0*ri2(i)*ti(i) 806 continue do 808 i=-5,nz do 808 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)) 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 zz=dz*real(i-11) zkr=zk*exp(-zki*zz) aven=0.0 do 508 j=1,mox-1 508 aven=aven+n(i,j) aven=aven/real(mox-1) do 50 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 50 continue ********************************************************************* * * * Limit the amplitude of the potential phi2(j). * * * ********************************************************************* do 580 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 580 continue ********************************************************************* * * * 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) if (eps1.gt.epc) then epc=eps1 end if 40 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) then go to 19 end if 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)-(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 90 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) then r=1.5 end if if (r.lt.1.0) then r=1.0 end if 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) then go to 19 end if if (nit.gt.50) then go to 19 end if 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) c write(*,*)'epc= ',epc c write(*,*)'nit= ',nit c write(*,*)'end of ellipse' return end ********************************************************************** **********************************************************************