module cparam integer nx,ny,nz,nxy,nxyz parameter (nx=102,ny=102,nz=21,nxy=nx*ny,nxyz=nx*ny*nz) end module cparam c................................................................. module cfield use cparam real(kind=2) vrotx(nx,ny,nz),vroty(nx,ny,nz),vx(nx,ny,nz), 1 vy(nx,ny,nz),vz(nx,ny,nz),eta(nx,ny,nz), 2 alp0(nx,ny,nz),om(nx,ny,nz),rho(nx,ny,nz) end module cfield c.................................................................... module cdata use cparam real(kind=2) x(nx),y(ny),z(nz),dx,dy,dz,rdx,rdy,rdz,ralp,rgal, 1 romeg,rm,xmin,ymin,zmin,xmax,ymax,zmax,pi,alam,asp, 2 etarat,rsph(nx,ny,nz),rcp(nx,ny,nz),th(nx,ny,nz), 3 ph(nx,ny,nz),dtau,rbrandt,alpb,tyrs,psc,omega0 4 ,feta integer nquad,imag,irf,iom,idia,ip,k0,irho end module cdata c.................................................................... c.....version with periodic bcs in x- and y-directions c.....and fourier analysis c.....modified to have eta propto grad(v), 04/01/01 c.....and alpha-quenching prop to rho, 10/01/01 C...TAKE CARE WITH NONCIRCULAR VELS!!!! program main use cdata use cfield use cparam implicit none real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz),tau, 1 alp(nx,ny,nz),fm(nx,ny,nz),vm(nx,ny,nz), 2 bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz),scl 3 ,etaphys,tscale integer nstep,npr,irst,nsave common/mask/fm,vm open(2,file='OUT',status='unknown') c..... open(1,file='in',status='old') read(1,*)nstep,npr,irst,nsave,dtau read(1,*)ralp,romeg,rm,alam,asp,rgal read(1,*)alpb,nquad,ip,irf,idia,irho read(1,*)imag,iom,rbrandt,etarat,scl read(1,*)feta close(1) c..... print*,nstep,npr,irst,nsave,dtau print*,ralp,romeg,rm,alam,asp,rgal print*,alpb,nquad,ip,irf,idia,irho print*,imag,iom,rbrandt,etarat,scl print*,nx,ny,nz print*,feta c..... write(2,*)nstep,npr,irst,nsave,dtau write(2,*)ralp,romeg,rm,alam,asp,rgal write(2,*)alpb,nquad,ip,irf,idia,irho write(2,*)imag,iom,rbrandt,etarat,scl write(2,*)nx,ny,nz write(2,*)feta c..... if(iom.ne.1.and.imag.eq.10)then print*,'iom=',iom,' with imag=',imag stop end if c..... call setbinf(bx,by,bz) call setmask(fm,vm) call init(alp,vm) if(imag.eq.10)then etaphys=omega0*(alam*rgal*1d3*psc)**2/romeg rm=1d5*alam*rgal*1d3*psc/etaphys tscale=(alam*rgal*1d3*psc)**2/etaphys/tyrs write(6,10)rm,etaphys,tscale write(2,10)rm,etaphys,tscale end if call setaxyz(irst,tau,ax,ay,az,bx,by,bz,scl) c call max(ax,'ax ') c call max(ay,'ay ') c call max(az,'az ') c call max(alp,'alp ') c call max(eta,'eta ') c call bfield(ax,ay,az,bx,by,bz) c call max(bx,'bx ') c call max(by,'by ') c call max(bz,'bz ') c..... call integrk(nstep,npr,ax,ay,az,bx,by,bz,alp,tau,nsave) close(2) call diagn(ax,ay,az,bx,by,bz) 10 format('imag=10, rm recalculated as',1p,e10.2,' etaphys=', 1 e10.2,' tscale=',e10.2,' yrs') end program main c.......................................................... subroutine setbinf(bx,by,bz) use cparam implicit none real(kind=2) bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz) bx=1d20 by=1d20 bz=1d20 end c.............................................................. subroutine setaxyz(irst,tau,ax,ay,az,bx,by,bz,scl) use cdata use cfield use cparam implicit none integer irst,nxx,nyy,nzz,i,j,k real(kind=2) tau,aph(nx,ny,nz),ax(nx,ny,nz),ay(nx,ny,nz), 1 az(nx,ny,nz) 2 ,bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz),e,er,scl if(irst.eq.1)then open(1,file='O.DAT',form='unformatted',status='old') read(1)tau,nxx,nyy,nzz if((nx.ne.nxx).or.ny.ne.nyy.or.(nz.ne.nzz))then print*,'inconsistent dimensions on read',nx,ny,nz,nxx,nyy,nzz stop end if read(1)ax read(1)ay read(1)az close(1) return else tau=0d0 az=0d0 do k=1,nz do j=2,ny-1 do i=2,nx-1 c ax(i,j,k)=(z(k)**2*float(2-ip)+z(k)*float(ip-1))*(asp-abs(z(k))) c ay(i,j,k)=(z(k)**2*float(2-ip)+z(k)*float(ip-1))*(asp-abs(z(k))) c ax(i,j,k)=-y(j)*z(k)*float(ip-1)-y(j)*float(2-ip) c ay(i,j,k)=x(i)*z(k)*float(ip-1)+x(i)*float(2-ip) ax(i,j,k)=-sin(pi*y(j))*z(k)*float(ip-1)-sin(pi*y(j))*float(2-ip) ay(i,j,k)=sin(pi*x(i))*z(k)*float(ip-1)+sin(pi*x(i))*float(2-ip) end do end do do i=2,nx-1 ax(i,1,k)=ax(i,ny-1,k) ay(i,1,k)=ay(i,ny-1,k) ax(i,ny,k)=ax(i,2,k) ay(i,ny,k)=ay(i,2,k) end do do j=2,ny-1 ax(1,j,k)=ax(nx-1,j,k) ay(1,j,k)=ay(nx-1,j,k) ax(nx,j,k)=ax(2,j,k) ay(nx,j,k)=ay(2,j,k) end do end do ! k loop call bcs(ax,ay,az) ax=ax*scl ay=ay*scl az=az*scl call bfield(ax,ay,az,bx,by,bz) do k=1,nz write(19,20)k,z(k),ax(11,11,k),ay(11,11,k),az(11,11,k) 1 ,bx(11,11,k),bz(11,11,k) end do end if call energ(e,er,bx,by,bz) print*,'initial energy e=',e 20 format(i4,f8.4,1p,5e10.2) end c........................................................... subroutine init(alp,vm) use cdata use cfield use cparam implicit none integer i,j,k,kk real(kind=2) fzv,vx0(nx,ny),vy0(nx,ny),fomz(nz),eta0(nx,ny), 1 scl,alp(nx,ny,nz),brandt,ruom,vm(nx,ny,nz) 2 ,gvx(nx,ny),gvy(nx,ny) c.... psc=3.09d18 tyrs=3.16d7 omega0=1d5/(1d3*psc) ruom=1d0 c.....set grid quantities if((nquad-1)*(nquad-2).ne.0)then print*,'cannot have nquad=',nquad stop end if if(nquad.eq.1.and.((ip-1)*(ip-2).ne.0))then print*,'cannot have nquad=1 with ip=',ip stop end if c..... if(nquad.eq.2)then kk=(nz+1)/2+1 k0=kk-1 else kk=2 k0=1 end if pi=2d0*asin(1d0) scl=(asp-alam)/2d0 xmin=-1d0 ymin=-1d0 xmax=1d0 ymax=1d0 zmax=asp dx=(xmax-xmin)/float(nx-3) dy=(ymax-ymin)/float(ny-3) if(nquad.eq.1)then zmin=0d0 dz=(zmax-zmin)/float(nz-1) else zmin=-asp dz=(zmax-zmin)/float(nz-1) end if rdx=1d0/dx rdy=1d0/dy rdz=1d0/dz do i=1,nx x(i)=xmin+float(i-1)*dx-dx end do do j=1,ny y(j)=ymin+float(j-1)*dy-dy end do do k=1,nz z(k)=zmin+float(k-1)*dz end do do k=1,nz do j=1,ny do i=1,nx rcp(i,j,k)=sqrt(x(i)**2+y(j)**2) rsph(i,j,k)=sqrt(x(i)**2+y(j)**2+z(k)**2) th(i,j,k)=acos(z(k)/rsph(i,j,k)) ph(i,j,k)=atan2(y(j),x(i)) end do end do end do c..... c.....set velocities, diffusion, alpha vx0=0d0 vy0=0d0 do k=1,nz fomz(k)=1d0 end do if(iom.eq.0.or.imag.ne.10)then do k=1,nz do j=1,ny do i=1,nx om(i,j,k)=brandt(rcp(i,j,k))*fomz(k) if(i*j*k.eq.1)print*,'Brandt law with imag,iom=',imag,iom if(rcp(i,j,k).gt.1d0)om(i,j,k)=brandt(1d0)*fomz(k) if(irf.eq.1)om(i,j,k)=om(i,j,k)-fomz(k)*brandt(ruom) om(i,j,k)=om(i,j,k)*romeg end do end do end do else if(iom.eq.1)then if(i*j*k.eq.1)print*,'iom=1, omega set from vx0,vy0' else print*,'not set up for iom=',iom stop end if c.....end of Brandt section if(imag.eq.10)then call calvppe(vx0,vy0) write(19,*)'calculate max, min of raw vx0,vy0 in km/s' print*,'calculate max, min of raw vx0,vy0 in km/s' call max1(vx0,'vx0 ') call max1(vy0,'vy0 ') call gradv(gvx,gvy,vx0,vy0) call strip(vx0,vy0,om) do k=1,nz do j=1,ny do i=1,nx om(i,j,k)=om(i,j,1)*fomz(k) end do end do end do end if c.....calculate vrot in km/s, then nondimensionalize call calvrot call max(vrotx,'vrotx ') call max(vroty,'vroty ') vrotx=vrotx*rm*alam vroty=vroty*rm*alam c..... call seteta(eta0,gvx,gvy) if(irho.eq.1)call readden(rho) do k=1,nz do j=1,ny do i=1,nx if(irho.eq.0)rho(i,j,k)=1d0 if(abs(z(k)).gt.alam)then alp0(i,j,k)=0d0 eta(i,j,k)=1d0+(etarat-1d0)*(1d0-exp(-(abs(z(k))-alam)**2/scl**2)) else alp0(i,j,k)=sin(pi*z(k)/alam) eta(i,j,k)=1d0 end if eta(i,j,k)=eta0(i,j)*eta(i,j,k) if(rcp(i,j,k).gt.1d0)alp0(i,j,k)=0d0 alp0(i,j,k)=alp0(i,j,k)*ralp*alam alp(i,j,k)=alp0(i,j,k) end do !i loop end do !j loop end do !k loop c do k=1,nz if(abs(z(k)).lt.alam)then fzv=1d0 else fzv=exp(-((abs(z(k))-alam)/(asp-alam))**2) end if do j=1,ny do i=1,nx vx(i,j,k)=vx0(i,j)*fzv vy(i,j,k)=vy0(i,j)*fzv end do end do end do !k loop write(19,*)'max and min of raw vx,vy' print*,'max and min of raw vx,vy' call max(vx,'vx ') call max(vy,'vy ') vx=vx*rm*alam vy=vy*rm*alam vz=vz*rm*alam do k=1,nz do j=1,ny do i=1,nx if(i.ne.1.and.i.ne.nx)then vx(i,j,k)=vx(i,j,k)-alam**2*float(idia) 1 *.5d0*(eta(i+1,j,k)-eta(i-1,j,k))*rdx end if if(j.ne.1.and.j.ne.ny)then vy(i,j,k)=vy(i,j,k)-alam**2*float(idia) 1 *.5d0*(eta(i,j+1,k)-eta(i,j-1,k))*rdy end if if(k.ne.1.and.k.ne.nz)then vz(i,j,k)=vz(i,j,k)-alam**2*float(idia) 1 *.5d0*(eta(i,j,k+1)-eta(i,j,k-1))*rdz end if end do end do end do do j=1,ny do i=1,nx vz(i,j,1)=vz(i,j,1)-0.5d0*float(idia)*alam**2*rdz* 1 (-eta(i,j,1)*1.5d0+2d0*eta(i,j,2)-0.5d0*eta(i,j,3)) vz(i,j,nz)=vz(i,j,nz)-0.5d0*float(idia)*alam**2*rdz* 1 (1.5d0*eta(i,j,nz)-2d0*eta(i,j,nz-1) 2 +0.5d0*eta(i,j,nz-2)) end do end do c.....now apply radial masking vx=vx*vm vy=vy*vm vz=vz*vm vrotx=vrotx*vm vroty=vroty*vm c..... c.....write diagnostics write(19,10)x write(19,10)y write(19,10)z do k=1,nz write(19,20)k,z(k),rcp(nx/3,ny/3,k),rsph(nx/3,ny/3,k), 1 alp(nx/3,ny/3,k),eta(nx/3,ny/3,k), 2 om(nx/3,ny/3,k),vrotx(nx/3,ny/3,k),vroty(nx/3,ny/3,k) 3 ,vx(nx/3,ny/3,k) end do write(19,30)kk,z(kk) do i=1,nx write(19,20)i,x(i),rcp(i,ny/3,kk ),rsph(i,ny/3,kk ), 1 alp(i,ny/3,kk ),eta(i,ny/3,kk ), 2 om(i,ny/3,kk ),vrotx(i,ny/3,kk ),vroty(i,ny/3,kk ) 3 ,vx(i,ny/3,kk ) end do call max(alp,'alp ') call max(eta,'eta ') 10 format(1p,10e10.2) 20 format(i4,f8.4,1p,8e10.2) 30 format('x-wise slice at',i4,f8.4) end c........................................................... subroutine derivxyz(ax,ay,az,rx,ry,rz,bx,by,bz,alp,tau) use cdata use cfield use cparam implicit none real(kind=2) gx,gy,gz,fx,fy,fz,hx,hy,hz,ax,ay,az,rx,ry,rz, 1 alp,bx,by,bz,dv,tau integer i,j,k dimension gx(nx,ny,nz),gy(nx,ny,nz),gz(nx,ny,nz),fx(nx,ny,nz), 1 fy(nx,ny,nz),fz(nx,ny,nz),hx(nx,ny,nz),hy(nx,ny,nz), 2 hz(nx,ny,nz),alp(nx,ny,nz),bx(nx,ny,nz), 3 by(nx,ny,nz),bz(nx,ny,nz),dv(nx,ny,nz) dimension ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz),rx(nx,ny,nz), 1 ry(nx,ny,nz),rz(nx,ny,nz) call bfield(ax,ay,az,bx,by,bz) if(alpb.gt.1d-10)call nalp(alp,alp0,bx,by,bz,rho,alpb) call mult1(gx,gy,gz,bx,by,bz) call mult2(hx,hy,hz,bx,by,bz) call mult3(fx,fy,fz,alp,bx,by,bz) call div(ax,ay,az,dv) do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 rx(i,j,k)=eta(i,j,k)*((ax(i+1,j,k)+ax(i-1,j,k))*rdx**2+ 1 (ax(i,j+1,k)+ax(i,j-1,k))*rdy**2+ 2 (ax(i,j,k+1)+ax(i,j,k-1))*rdz**2- 3 2d0*ax(i,j,k)*(rdx**2+rdy**2+rdz**2)- x (dv(i+1,j,k)-dv(i-1,j,k))*rdx*0.5d0)*alam**2 4 +fx(i,j,k)+gx(i,j,k)+hx(i,j,k) ry(i,j,k)=eta(i,j,k)*((ay(i+1,j,k)+ay(i-1,j,k))*rdx**2+ 1 (ay(i,j+1,k)+ay(i,j-1,k))*rdy**2+ 2 (ay(i,j,k+1)+ay(i,j,k-1))*rdz**2- 3 2d0*ay(i,j,k)*(rdx**2+rdy**2+rdz**2)- x (dv(i,j+1,k)-dv(i,j-1,k))*rdy*0.5d0)*alam**2 4 +fy(i,j,k)+gy(i,j,k)+hy(i,j,k) rz(i,j,k)=eta(i,j,k)*((az(i+1,j,k)+az(i-1,j,k))*rdx**2+ 1 (az(i,j+1,k)+az(i,j-1,k))*rdy**2+ 2 (az(i,j,k+1)+az(i,j,k-1))*rdz**2- 3 2d0*az(i,j,k)*(rdx**2+rdy**2+rdz**2)- x (dv(i,j,k+1)-dv(i,j,k-1))*rdz*0.5d0)*alam**2 4 +fz(i,j,k)+gz(i,j,k)+hz(i,j,k) end do !i loop end do !j loop do i=2,nx-1 rx(i,1,k)=rx(i,ny-1,k) ry(i,1,k)=ry(i,ny-1,k) rz(i,1,k)=rz(i,ny-1,k) rx(i,ny,k)=rx(i,2,k) ry(i,ny,k)=ry(i,2,k) rz(i,ny,k)=rz(i,2,k) end do do j=2,ny-1 rx(1,j,k)=rx(nx-1,j,k) ry(1,j,k)=ry(nx-1,j,k) rz(1,j,k)=rz(nx-1,j,k) rx(nx,j,k)=rx(2,j,k) ry(nx,j,k)=ry(2,j,k) rz(nx,j,k)=rz(2,j,k) end do end do !k loop end c.................................................................... subroutine mult1(gx,gy,gz,bx,by,bz) use cfield use cparam implicit none real(kind=2) gx(nx,ny,nz),gy(nx,ny,nz),gz(nx,ny,nz), 1 bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz) c do l=1,nxyz c gx(l)=vroty(l)*bz(l) c gy(l)=-vrotx(l)*bz(l) c gz(l)=vrotx(l)*by(l)-vroty(l)*bx(l) c NEED TO REDIMENSION ARRAYS c INTEGER l c end do gx=vroty*bz gy=-vrotx*bz gz=vrotx*by-vroty*bx end c........................................................................ subroutine mult2(gx,gy,gz,bx,by,bz) use cfield use cparam implicit none real(kind=2) gx(nx,ny,nz),gy(nx,ny,nz),gz(nx,ny,nz) real(kind=2) bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz) gx=vy*bz-vz*by gy=-vx*bz+vz*bx gz=vx*by-vy*bx end c....................................................................... subroutine mult3(gx,gy,gz,alp,bx,by,bz) use cparam implicit none real(kind=2) gx(nx,ny,nz),gy(nx,ny,nz),gz(nx,ny,nz), 1 bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz), 1 alp(nx,ny,nz) gx=alp*bx gy=alp*by gz=alp*bz end c....................................................................... subroutine bfield(ax,ay,az,bx,by,bz) use cdata use cparam implicit none real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz) real(kind=2) bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz) integer i,j,k do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 bx(i,j,k)=-(ay(i,j,k+1)-ay(i,j,k-1))*0.5d0*rdz 1 +(az(i,j+1,k)-az(i,j-1,k))*0.5d0*rdy by(i,j,k)=+(ax(i,j,k+1)-ax(i,j,k-1))*0.5d0*rdz- 1 (az(i+1,j,k)-az(i-1,j,k))*rdx*0.5d0 bz(i,j,k)=(ay(i+1,j,k)-ay(i-1,j,k))*0.5d0*rdx- 1 (ax(i,j+1,k)-ax(i,j-1,k))*0.5d0*rdy end do end do end do c..... do k=2,nz-1 do j=2,ny-1 bx(1,j,k)=bx(nx-1,j,k) by(1,j,k)=by(nx-1,j,k) bz(1,j,k)=bz(nx-1,j,k) bx(nx,j,k)=bx(2,j,k) by(nx,j,k)=by(2,j,k) bz(nx,j,k)=bz(2,j,k) end do end do do k=2,nz-1 do i=2,nx-1 bx(i,1,k)=bx(i,ny-1,k) by(i,1,k)=by(i,ny-1,k) bz(i,1,k)=bz(i,ny-1,k) bx(i,ny,k)=bx(i,2,k) by(i,ny,k)=by(i,2,k) bz(i,ny,k)=bz(i,2,k) end do end do k=1 if(nquad.eq.1)then do j=2,ny-1 do i=2,nx-1 if(ip.eq.2)then bx(i,j,k)=(2d0*bx(i,j,k+1)-0.5d0*bx(i,j,k+2))/1.5d0 by(i,j,k)=(2d0*by(i,j,k+1)-0.5d0*by(i,j,k+2))/1.5d0 bz(i,j,k)=0d0 else print*,'revise bcs for ip=1 and nquad=1' stop end if bz(i,j,k)=(ay(i+1,j,k)-ay(i-1,j,k))*0.5d0*rdx- 1 (ax(i,j+1,k)-ax(i,j-1,k))*0.5d0*rdy end do end do else ! nquad print*,'b.c. on k=1 needs fixing for nquad=2' stop end if ! nquad k=nz do j=2,ny-1 do i=2,nx-1 bx(i,j,k)=-(0.5d0*ay(i,j,k-2)-2d0*ay(i,j,k-1)+ 2 1.5d0*ay(i,j,k))*rdz 1 +(az(i,j+1,k)-az(i,j-1,k))*0.5d0*rdy by(i,j,k)=+(0.5d0*ax(i,j,k-2)-2d0*ax(i,j,k-1)+ 2 1.5d0*ax(i,j,k))*rdz- 1 (az(i+1,j,k)-az(i-1,j,k))*rdx*0.5d0 bz(i,j,k)=(ay(i+1,j,k)-ay(i-1,j,k))*0.5d0*rdx- 1 (ax(i,j+1,k)-ax(i,j-1,k))*0.5d0*rdy end do end do c.....values at edges set as 1d20 initially do i=2,nx-1 bx(i,1,1)=bx(i,ny-1,1) by(i,1,1)=by(i,ny-1,1) bz(i,1,1)=bz(i,ny-1,1) bx(i,ny,1)=bx(i,2,1) by(i,ny,1)=by(i,2,1) bz(i,ny,1)=bz(i,2,1) end do do j=2,ny-1 bx(1,j,1)=bx(nx-1,j,1) by(1,j,1)=by(nx-1,j,1) bz(1,j,1)=bz(nx-1,j,1) bx(nx,j,1)=bx(2,j,1) by(nx,j,1)=by(2,j,1) bz(nx,j,1)=bz(2,j,1) end do end c............................................................. subroutine calvrot c.....calculates rotational velocity in km/s use cdata use cfield use cparam implicit none vrotx=-om*sin(ph)*rcp*rgal vroty=om*cos(ph)*rcp*rgal end c.......................................................... subroutine calvppe(vx0,vy0) use cdata use cfield use cparam implicit none integer nx1,ny1,i,j,one real(kind=2) vx0(nx,ny),vy0(nx,ny),vx00(nx-2,ny-2),vy00(nx-2,ny-2) 1 ,tau open(1,file='ppe_vxy_100.data',status='old') if(nx.le.320)then read(1,*)vx00 read(1,*)vy00 read(1,*)nx1,ny1 else read(1)vx00 read(1)vy00 read(1)nx1,ny1 end if if(nx1.ne.nx-2.or.ny1.ne.ny-2)then print*,'inconsistent data with imag = 10' print*,nx,ny,nx1,ny1 stop end if write(2,30)nx1,ny1 write(6,30)nx1,ny1 close(1) do j=2,ny-1 do i=2,nx-1 vx0(i,j)=vx00(i-1,j-1) vy0(i,j)=vy00(i-1,j-1) end do end do do j=2,ny vx0(1,j)=vx0(nx-1,j) vy0(1,j)=vy0(nx-1,j) vx0(nx,j)=vx0(2,j) vy0(nx,j)=vy0(2,j) end do do i=2,nx-1 vx0(i,1)=vx0(i,ny-1) vy0(i,1)=vy0(i,ny-1) vx0(i,nx)=vx0(i,2) vy0(i,nx)=vy0(i,2) end do vx0(1,1)=(2d0*vx0(2,1)-vx0(3,1)+2d0*vx0(1,2)-vx0(1,3))/2d0 vy0(1,1)=(2d0*vy0(2,1)-vy0(3,1)+2d0*vy0(1,2)-vy0(1,3))/2d0 vx0(nx,ny)=(2d0*vx0(nx-1,1)-vx0(nx-2,1)+2d0*vx0(1,ny-1)- 1 vx0(1,ny-2))/2d0 vy0(nx,ny)=(2d0*vy0(nx-1,1)-vy0(nx-2,1)+2d0*vy0(1,ny-1)- 1 vy0(1,ny-2))/2d0 vy0(1,ny)=(2d0*vy0(2,ny)-vy0(3,ny)+2d0*vy0(1,ny-1)- 1 vy0(1,ny-2))/2d0 vx0(1,ny)=(2d0*vx0(2,ny)-vx0(3,ny)+2d0*vx0(1,ny-1)- 1 vx0(1,ny-2))/2d0 vx0(nx,1)=(2d0*vx0(2,nx-1)-vx0(nx-2,1)+2d0*vx0(nx,2)- 1 vx0(nx,3))/2d0 vy0(nx,1)=(2d0*vy0(2,nx-1)-vy0(nx-2,1)+2d0*vy0(nx,2)- 1 vy0(nx,3))/2d0 one=1 tau=0d0 open(7,file='fort.OV',status='unknown') write(7,*)tau,nx,ny,one write(7,*)vx0 write(7,*)vy0 close(7) 30 format('read velocity data in z=0 with nx,ny=',2i4) end c................................................................ subroutine integrk(nstep,npr,ax,ay,az,bx,by,bz,alp,tau, 1 nsave) use cdata use cfield use cparam implicit none integer nstep,npr,n,kk,nsave,ienerg real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz), 1 bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz), 2 alp(nx,ny,nz),tau,p1,p2,e,erat 3 ,er,ef,et,e0,e1,e2,e3 if(nquad.eq.2)then kk=(nz+1)/2+1 else kk=2 end if do n=1,nstep call rk2(ax,ay,az,bx,by,bz,alp,dtau,tau) if(n/npr*npr.eq.n)then call bfield(ax,ay,az,bx,by,bz) call energ(e,erat,bx,by,bz) if(nquad.eq.2)then call parity(p1,p2,bx,by,bz) else p2=float(2*ip-3) p1=e end if write(6,10)n,tau,e,p1,p2,erat,bx(3*nx/4,3*ny/4,kk), 1 e0,e1 write(2,10)n,tau,e,p1,p2,erat,bx(3*nx/4,3*ny/4,kk), 1 e0,e1 call energxy(bx,by,bz,er,ef,et,e0,e1,e2,e3,ienerg) write(20,20)n,tau,e,erat,er,ef,e0,e1,e2,e3 open(29,file='snap',status='unknown') write(29,10)n,tau,e,p1,p2,erat,e0,e1,e2,e3 close(29) end if if((n/nsave)*nsave.eq.n.and.(n.ne.nstep))then open(1,file='O1.DAT',status='unknown',form='unformatted') write(1)tau,nx,ny,nz write(1)ax write(1)ay write(1)az close(1) print*,'intermediate data written at tau=',tau end if end do open(1,file='O.DAT',status='unknown',form='unformatted') write(1)tau,nx,ny,nz write(1)ax write(1)ay write(1)az close(1) open(1,file='OF.DAT',status='unknown',form='unformatted') write(1)tau,nx,ny,nz write(1)bx write(1)by write(1)bz close(1) ienerg=2 call energxy(bx,by,bz,er,ef,et,e0,e1,e2,e3,ienerg) 10 format(i8,f9.4,1p,2e11.4,0p,f15.12,1p,e10.3,4e10.2) 20 format(i8,f9.4,1p,8e10.2) end c............................................................... subroutine energ(e,erat,bx,by,bz) use cdata use cparam implicit none integer n c real(kind=2) e,bx(nxyz),by(nxyz),bz(nxyz),fm(nxyz),erat,em,eaz real(kind=2) e,bx(nx,ny,nz),by(nx,ny,nz),bz(nxyz), 1 fm(nxyz),erat,em,eaz,br(nxyz),bf(nxyz),vm(nxyz) common/mask/fm,vm eaz=0d0 em=0d0 call mult4(bx,by,ph,br,bf) do n=1,nxyz c e=e+(bx(n)**2+by(n)**2+bz(n)**2)*fm(n) em=em+(br(n)**2+bz(n)**2)*fm(n) eaz=eaz+(bf(n)**2)*fm(n) end do e=eaz+em if(eaz.gt.1d-12)then erat=em/eaz else erat=99d0 end if e=e*0.5d0*dx*dy*dz end c....................................................................... subroutine mult4(bx,by,ph,br,bf) use cparam implicit none integer k real(kind=2) bx(nxyz),by(nxyz),ph(nxyz),br(nxyz),bf(nxyz) do k=1,nxyz br(k)=bx(k)*cos(ph(k))+by(k)*sin(ph(k)) bf(k)=-bx(k)*sin(ph(k))+by(k)*cos(ph(k)) end do end c....................................................................... subroutine setmask(fm,vm) use cdata use cparam implicit none real(kind=2)fm(nx,ny,nz),vm(nx,ny,nz),xx,yy,cv,rr integer i,j,k cv=0.9d0 write(6,10)cv write(19,10)cv write(2,10)cv 10 format('vels go to zero for r gt',f8.4) fm=0d0 do k=2,nz-1 do j=3,ny-2 do i=3,nx-2 fm(i,j,k)=1d0 end do end do end do do k=1,nz do j=2,ny-1 do i=2,nx-1,nx-3 fm(i,j,k)=0.5d0 end do end do do j=2,ny-1,ny-3 do i=2,nx-1 fm(i,j,k)=0.5d0 end do end do end do !k loop if(nquad.eq.1)then do j=2,ny-1 do i=2,nx-1 fm(i,j,1)=0.5d0 end do end do end if c do k=1,nz do j=1,ny do i=1,nx if(x(i).gt.1d0)then xx=x(i)-2d0 else if(x(i).lt.-1d0)then xx=x(i)+2d0 else xx=x(i) end if if(y(j).gt.1d0)then yy=y(j)-2d0 else if(y(j).lt.-1d0)then yy=y(j)+2d0 else yy=y(j) end if rr=sqrt(xx**2+yy**2) if(rr.gt.cv)then vm(i,j,k)=exp(-2d0*((rr-cv)/(1d0-cv))**2) else vm(i,j,k)=1d0 end if end do end do end do end c....................................................................... subroutine nalp(alp,alp0,bx,by,bz,rho,alpb) use cparam implicit none real(kind=2) alp(nxyz),alp0(nxyz),bx(nxyz),by(nxyz),bz(nxyz), 1 alpb,rho(nxyz) integer n do n=1,nxyz alp(n)=alp0(n)/(1d0+alpb*(bx(n)**2+by(n)**2+bz(n)**2)/ 1 rho(n)) end do end c....................................................................... function brandt(r) use cdata use cparam implicit none real(kind=2) r,brandt brandt=1d0/sqrt(1d0+(r/rbrandt)**2) end c......................................................................... subroutine eveodd(a,aeve,aodd) use cparam implicit none real(kind=2) a(nx,ny,nz),aeve(nx,ny,nz),aodd(nx,ny,nz) integer i,j,k do k=1,nz do j=1,ny do i=1,nx aeve(i,j,k)=0.5d0*(a(i,j,k)+a(i,j,nz-k+1)) aodd(i,j,k)=0.5d0*(a(i,j,k)-a(i,j,nz-k+1)) end do end do end do end c.......................................................................... subroutine parity(p1,p2,bx,by,bz) use cdata use cparam implicit none real(kind=2) bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz), 1 bxeve(nx,ny,nz),byeve(nx,ny,nz),bzeve(nx,ny,nz), 2 bxodd(nx,ny,nz),byodd(nx,ny,nz),bzodd(nx,ny,nz), 3 p1,p2,eeve,eodd,er call eveodd(bx,bxeve,bxodd) call eveodd(by,byeve,byodd) call eveodd(bz,bzeve,bzodd) call energ(eeve,er,bxeve,byeve,bzodd) call energ(eodd,er,bxodd,byodd,bzeve) p1=eeve+eodd if(p1.lt.1d-14)then p2=99d0 else p2=(eeve-eodd)/p1 end if end c.......................................................................... subroutine bcs(ax,ay,az) use cdata use cparam implicit none real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz) integer i,j,k k=nz do j=1,ny do i=1,nx ax(i,j,k)=(2d0*ax(i,j,k-1)-0.5d0*ax(i,j,k-2))/1.5d0 ay(i,j,k)=(2d0*ay(i,j,k-1)-0.5d0*ay(i,j,k-2))/1.5d0 az(i,j,k)=0d0 end do end do k=1 if(nquad.eq.2)then do j=1,ny do i=1,nx ax(i,j,k)=(2d0*ax(i,j,2)-0.5d0*ax(i,j,3))/1.5d0 ay(i,j,k)=(2d0*ay(i,j,2)-0.5d0*ay(i,j,3))/1.5d0 az(i,j,k)=0d0 end do end do else if(ip.eq.1)then c odd parity do j=1,ny do i=1,nx az(i,j,k)=0d0 ax(i,j,k)=(2d0*ax(i,j,k+1)-0.5d0*ax(i,j,k+2))/1.5d0 ay(i,j,k)=(2d0*ay(i,j,k+1)-0.5d0*ay(i,j,k+2))/1.5d0 end do end do else c even parity do j=1,ny do i=1,nx az(i,j,k)=(2d0*az(i,j,k+1)-0.5d0*az(i,j,k+2))/1.5d0 ax(i,j,k)=0d0 ay(i,j,k)=0d0 end do end do end if end if end c.......................................................................... subroutine div(ax,ay,az,dv) use cdata use cparam implicit none integer i,j,k real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz), 1 dv(nx,ny,nz),ddx(nx,ny,nz),ddy(nx,ny,nz), 2 ddz(nx,ny,nz) do k=1,nz do j=1,ny do i=2,nx-1 ddx(i,j,k)=(ax(i+1,j,k)-ax(i-1,j,k))*rdx*0.5d0 end do ddx(1,j,k)=(ax(2,j,k)-ax(nx-2,j,k))*rdx*0.5d0 ddx(nx,j,k)=(ax(3,j,k)-ax(nx-1,j,k))*rdx*0.5d0 end do ! end j loop do j=2,ny-1 do i=1,nx ddy(i,j,k)=(ay(i,j+1,k)-ay(i,j-1,k))*rdy*0.5d0 end do end do ! end j loop do i=1,nx ddy(i,1,k)=(ay(i,2,k)-ay(i,ny-2,k))*rdy*0.5d0 ddy(i,ny,k)=(ay(i,3,k)-ay(i,ny-1,k))*rdy*0.5d0 end do ! end i loop end do ! end k loop do j=1,ny do i=1,nx ddz(i,j,1)=(-1.5d0*az(i,j,1)+2d0*az(i,j,2)-0.5d0*az(i,j,3))*rdz ddz(i,j,nz)=(1.5d0*az(i,j,nz)-2d0*az(i,j,nz-1)+0.5d0*az(i,j,nz-2)) 1 *rdz end do end do do k=2,nz-1 do j=1,ny do i=1,nx ddz(i,j,k)=(az(i,j,k+1)-az(i,j,k-1))*rdz*0.5d0 end do end do end do do k=1,nz do j=1,ny do i=1,nx dv(i,j,k)=ddx(i,j,k)+ddy(i,j,k)+ddz(i,j,k) end do end do end do end c............................................................ subroutine rk2(ax,ay,az,bx,by,bz,alp,dt,t) use cparam implicit none c.....on entry ax,ay,az contains current function values, on exit, the updated c.....values. derivxyz calculates derivatives using current ax,ay,az at time t real(kind=2) ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz), 1 rx(nx,ny,nz),ry(nx,ny,nz),rz(nx,ny,nz), 2 sx1(nx,ny,nz),sy1(nx,ny,nz),sz1(nx,ny,nz), 3 rx1(nx,ny,nz),ry1(nx,ny,nz),rz1(nx,ny,nz), x bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz), 4 alp(nx,ny,nz),dt,t integer i,j,k call derivxyz(ax,ay,az,rx,ry,rz,bx,by,bz,alp,t) do k=1,nz do j=1,ny do i=1,nx sx1(i,j,k)=ax(i,j,k)+dt*rx(i,j,k) sy1(i,j,k)=ay(i,j,k)+dt*ry(i,j,k) sz1(i,j,k)=az(i,j,k)+dt*rz(i,j,k) end do end do end do call bcs(sx1,sy1,sz1) t=t+dt call derivxyz(sx1,sy1,sz1,rx1,ry1,rz1,bx,by,bz,alp,t) do k=1,nz do j=1,ny do i=1,nx ax(i,j,k)=ax(i,j,k)+0.5d0*dt*(rx(i,j,k)+rx1(i,j,k)) ay(i,j,k)=ay(i,j,k)+0.5d0*dt*(ry(i,j,k)+ry1(i,j,k)) az(i,j,k)=az(i,j,k)+0.5d0*dt*(rz(i,j,k)+rz1(i,j,k)) end do end do end do call bcs(ax,ay,az) return end c................................................................ subroutine max(a,ch) use cparam implicit none integer i,j,k,mx,my,mz real(kind=2) a(nx,ny,nz),st character*4 ch st=0d0 mx=0 my=0 mz=0 do k=1,nz do j=1,ny do i=1,nx if(abs(a(i,j,k)).gt.st)then mx=i my=j mz=k st=abs(a(i,j,k)) end if end do end do end do write(19,*)ch,' max=',st,' at',mx,my,mz print*,ch,' max=',st,' at',mx,my,mz end c................................................................ subroutine max1(a,ch) use cparam implicit none integer i,j,mx,my real(kind=2) a(nx,ny),st character*4 ch st=0d0 mx=0 my=0 do j=1,ny do i=1,nx if(abs(a(i,j)).gt.st)then mx=i my=j st=abs(a(i,j)) end if end do end do write(19,*)ch,' max=',st,' at',mx,my print*,ch,' max=',st,' at',mx,my end c........................................................................ subroutine diagn(ax,ay,az,bx,by,bz) use cparam implicit none integer i,j,k,ix,jy,kz real(kind=2)ax(nx,ny,nz),ay(nx,ny,nz),az(nx,ny,nz),bx(nx,ny,nz), 1 by(nx,ny,nz),bz(nx,ny,nz) ix=(nx+1)/2 jy=(ny+1)/2 kz=nz/3 write(18,10)jy,kz do i=1,nx write(18,20)i,ax(i,jy,kz),ay(i,jy,kz),az(i,jy,kz), 1 bx(i,jy,kz),by(i,jy,kz),bz(i,jy,kz) end do write(18,30)ix,kz do j=1,ny write(18,20)j,ax(ix,j,kz),ay(ix,j,kz),az(ix,j,kz), 1 bx(ix,j,kz),by(ix,j,kz),bz(ix,j,kz) end do ix=ix+2 jy=jy+2 write(18,40)ix,jy do k=1,nz write(18,20)k,ax(ix,jy,k),ay(ix,jy,k),az(ix,jy,k), 1 bx(ix,jy,k),by(ix,jy,k),bz(ix,jy,k) end do 10 format('values of a,b at j,k=',2i4) 30 format('values of a,b at i,k=',2i4) 40 format('values of a,b at i,j=',2i4) 20 format(i4,1p,6e10.2) end c................................................................................ subroutine fouran(bf,a0,a1,a2,a3,b1,b2,b3,il,ih,ni,nj,dphi) use cdata use cparam implicit none integer il,ih,ni,nj,j,i real(kind=2) a0(ni),a1(ni),a2(ni),b1(ni),b2(ni),bf(ni,nj) 1 ,a3(ni),b3(ni),co1,co2,co3,si1,si2,si3,dphi common/four/co1(401),co2(401),co3(401),si1(401),si2(401),si3(401) do i=il,ih a0(i)=0.0 a1(i)=0.0 a2(i)=0.0 b1(i)=0.0 a3(i)=0.0 b3(i)=0.0 b2(i)=0.0 end do do i=il,ih do j=1,nj-1 a0(i)=a0(i)+bf(i,j)*dphi*0.5/pi a1(i)=a1(i)+bf(i,j)*co1(j)*dphi/pi a2(i)=a2(i)+bf(i,j)*co2(j)*dphi/pi b1(i)=b1(i)+bf(i,j)*si1(j)*dphi/pi a3(i)=a3(i)+bf(i,j)*co3(j)*dphi/pi b3(i)=b3(i)+bf(i,j)*si3(j)*dphi/pi b2(i)=b2(i)+bf(i,j)*si2(j)*dphi/pi end do end do return end c...................................................................... subroutine enfour(a0,a1,a2,a3,b1,b2,b3,e0,e1,e2,e3, 1 r,dr,ni) use cparam implicit none c.....calculates phi-energy in separate fourier modes integer i,ni real(kind=2) a0(ni),a1(ni),a2(ni),a3(ni),b1(ni),b2(ni),b3(ni) real(kind=2) r(ni),pi,e0,e1,e2,e3,dr pi=2d0*asin(1d0) e0=0.0d0 e1=0.0d0 e2=0.0d0 e3=0.0d0 do i=1,ni-1 e0=e0+r(i)*a0(i)**2+r(i+1)*a0(i+1)**2 e1=e1+r(i)*(b1(i)**2+a1(i)**2)+r(i+1)*(b1(i+1)**2+a1(i+1)**2) e2=e2+r(i)*(b2(i)**2+a2(i)**2)+r(i+1)*(b2(i+1)**2+a2(i+1)**2) e3=e3+r(i)*(b3(i)**2+a3(i)**2)+r(i+1)*(b3(i+1)**2+a3(i+1)**2) end do e0=e0*dr*0.5d0*2.d0*pi e1=e1*dr*0.5d0*pi e2=e2*dr*0.5d0*pi e3=e3*dr*0.5d0*pi c.....0.25 changed to 0.5 on 15/03/00 return end c...................................................................... subroutine energxy(bx,by,bz,er,ef,et,e0,e1,e2,e3,ienerg) use cdata use cparam implicit none integer ni,nj,k,ienerg,ix,iy,j,i parameter(ni=nx,nj=ny) real(kind=2) bx(nx,ny,nz),by(nx,ny,nz),bz(nx,ny,nz) real(kind=2) a0(ni),a1(ni),a2(ni),a3(ni),b1(ni),b2(ni),b3(ni) real(kind=2) br(ni,nj),bf(ni,nj),rr(ni) real(kind=2) dph,rmin,dr,phi,rrr,xx,yy,fx,fy,bbx1,bbx2,bby1,bby2, 1 bbx,bby,et,e0,e1,e2,e3,er0,er1,er2,er3,er,ef 2 ,pitch,pd if(nj.gt.401)then print*,'nj set gt 401. co1 etc incorrectly dimensioned',nj stop end if if(nquad.eq.1)then k=1 else k=(nz+1)/2 end if dph=2d0*pi/float(nj-1) rmin=0d0 dr=(xmax-rmin)/float(ni-1) if(ienerg.eq.0)then call coeffs(dph,nj) end if do j=1,nj phi=float(j-1)*dph do i=1,ni rrr=rmin+float(i-1)*dr rr(i)=rrr xx=rrr*cos(phi) yy=rrr*sin(phi) ix=(xx-xmin+1d-10)/dx+2 iy=(yy-ymin+1d-10)/dy+2 if(ix.lt.2.or.ix.gt.nx-1)then print*,'ix out of range',ix,xx,rr,phi,i,j end if if(iy.lt.2.or.iy.gt.ny-1)then print*,'iy out of range',iy,yy,rr,phi,i,j end if if(ix.eq.nx-1)ix=nx-2 if(iy.eq.ny-1)iy=ny-2 fx=(xx-x(ix))/dx fy=(yy-y(iy))/dx if((fx.gt.1d0+1d-12).or.(fx.lt.0d0-1d-12))then print*,'energxy: fx,ix,x(ix),i,j',fx,ix,x(ix),i,j,rrr,phi,xx,yy end if if((fy.gt.1d0+1d-12).or.(fy.lt.0d0-1d-12))then print*,'energxy: fy,iy,y(iy),i,j',fy,iy,y(iy),i,j,rrr,phi,xx,yy end if bbx1=bx(ix,iy,k)*(1d0-fx)+fx*bx(ix+1,iy,k) bbx2=bx(ix,iy+1,k)*(1d0-fx)+fx*bx(ix+1,iy+1,k) bby1=by(ix,iy,k)*(1d0-fx)+fx*by(ix+1,iy,k) bby2=by(ix,iy+1,k)*(1d0-fx)+fx*by(ix+1,iy+1,k) bbx=bbx1*(1d0-fy)+bbx2*fy bby=bby1*(1d0-fy)+bby2*fy br(i,j)=bbx*cos(phi)+bby*sin(phi) bf(i,j)=-bbx*sin(phi)+bby*cos(phi) end do end do if(ienerg.eq.0)ienerg=1 call fouran(bf,a0,a1,a2,a3,b1,b2,b3,1,ni,ni,nj,dph) call enfour(a0,a1,a2,a3,b1,b2,b3,e0,e1,e2,e3,rr,dr,ni) call fouran(br,a0,a1,a2,a3,b1,b2,b3,1,ni,ni,nj,dph) call enfour(a0,a1,a2,a3,b1,b2,b3,er0,er1,er2,er3,rr,dr,ni) e0=e0*alam e1=e1*alam e2=e2*alam e3=e3*alam er0=er0*alam er1=er1*alam er2=er2*alam er3=er3*alam ef=e0+e1+e2+e3 er=er0+er1+er2+er3 e0=e0+er0 e1=e1+er1 e2=e2+er2 e3=e3+er3 et=e0+e1+e2+e3 if(ienerg.eq.2)then do j=1,nj-5,(nj/4) phi=float(j-1)*dph write(21,10)j,dph do i=1,ni pitch=atan2(bf(i,j),br(i,j)) pd=pitch*180d0/pi write(21,20)i,rr(i),br(i,j),bf(i,j),pitch,pd end do end do end if 10 format('for j,phi=',i4,f8.4,' rad,Br,Bphi,pitch,pitch(deg)') 20 format(i4,f8.4,1p,4e11.3) return end c.................................................................... subroutine coeffs(dphi,nj) use cdata use cparam implicit none real(kind=2) co1,co2,co3,si1,si2,si3,phi,dphi integer j,nj common/four/co1(401),co2(401),co3(401),si1(401),si2(401),si3(401) do j=1,nj phi=float(j-1)*dphi co1(j)=cos(phi) co2(j)=cos(2d0*phi) co3(j)=cos(3d0*phi) si3(j)=sin(3d0*phi) si1(j)=sin(phi) si2(j)=sin(2d0*phi) end do return end c...................................................................... subroutine strip(vx0,vy0,om) c.....returns om(i,j,z=0) in km/s/kpc and vx0,vy0 in km/s. c.....after subtracting vcirc use cdata use cparam implicit none integer ni,nj,ix,iy,j,i,one parameter(ni=nx/2,nj=ny) real(kind=2) om(nx,ny,nz),tau real(kind=2) a0(ni),a1(ni),a2(ni),a3(ni),b1(ni),b2(ni),b3(ni) real(kind=2) vf0(ni,nj),rr(ni),vx0(nx,ny),vy0(nx,ny) real(kind=2) dph,rmin,dr,phi,rrr,xx,yy,fx,fy,bbx1,bbx2,bby1,bby2, 1 bbx,bby,vff,omega dph=2d0*pi/float(nj-1) rmin=0d0 dr=(xmax-rmin)/float(ni-1) call coeffs(dph,nj) do j=1,nj phi=float(j-1)*dph do i=1,ni rrr=rmin+float(i-1)*dr rr(i)=rrr xx=rrr*cos(phi) yy=rrr*sin(phi) ix=(xx-xmin+1d-10)/dx+2 iy=(yy-ymin+1d-10)/dy+2 if(ix.lt.2.or.ix.gt.nx-1)then print*,'ix out of range',ix,xx,rr,phi,i,j end if if(iy.lt.2.or.iy.gt.ny-1)then print*,'iy out of range',iy,yy,rr,phi,i,j end if if(ix.eq.nx-1)ix=nx-2 if(iy.eq.ny-1)iy=ny-2 fx=(xx-x(ix))/dx fy=(yy-y(iy))/dx if((fx.gt.1d0+1d-12).or.(fx.lt.0d0-1d-12))then print*,'strip: fx,ix,x(ix),i,j',fx,ix,x(ix),i,j,rrr,phi,xx,yy end if if((fy.gt.1d0+1d-12).or.(fy.lt.0d0-1d-12))then print*,'strip: fy,iy,y(iy),i,j',fy,iy,y(iy),i,j,rrr,phi,xx,yy end if bbx1=vx0(ix,iy)*(1d0-fx)+fx*vx0(ix+1,iy) bbx2=vx0(ix,iy+1)*(1d0-fx)+fx*vx0(ix+1,iy+1) bby1=vy0(ix,iy)*(1d0-fx)+fx*vy0(ix+1,iy) bby2=vy0(ix,iy+1)*(1d0-fx)+fx*vy0(ix+1,iy+1) bbx=bbx1*(1d0-fy)+bbx2*fy bby=bby1*(1d0-fy)+bby2*fy vf0(i,j)=-bbx*sin(phi)+bby*cos(phi) end do end do call fouran(vf0,a0,a1,a2,a3,b1,b2,b3,1,ni,ni,nj,dph) do j=1,nx do i=1,ny rrr=rcp(i,j,1) if(rrr.gt.1d0)then vff=a0(ni) else ix=(rrr-rmin+1d-10)/dr+1 fx=(rrr-rr(ix))/dr if(ix.eq.ni)ix=ni-1 if((fx.gt.1d0+1d-12).or.(fx.lt.0d0-1d-12))then print*,'strip1: fx,ix,x(ix),i,',fx,ix,x(ix),i,rrr,phi end if vff=(1d0-fx)*a0(ix)+fx*a0(ix+1) end if c.....remove rotation w.r.t. rotating frome from vx0,vy0: now applied c.....separately via om. vx0(i,j)=vx0(i,j)-(-vff*sin(ph(i,j,1))) vy0(i,j)=vy0(i,j)-vff*cos(ph(i,j,1)) if(rrr.gt.0d0)then om(i,j,1)=vff/(rcp(i,j,1)*rgal) c.....om(i,j,k) in units of km/s/kpc. Only used explicitly to c.....calculate vrotx,vroty, when rgal factor explicitly applied to give km/s else om(i,j,1)=1d8 end if end do end do call max(om,'om ') c one=1 tau=0d0 open(7,file='fort.OV1',status='unknown') write(7,*)tau,nx,ny,one write(7,*)vx0 write(7,*)vy0 close(7) ix=0 do j=1,ny do i=1,nx if(om(i,j,1).eq.1d8)then ix=ix+1 if(i.ne.1.and.(i.ne.nx).and.j.ne.1.and.(j.ne.ny))then om(i,j,1)=0.25d0*(om(i+1,j,1)+om(i-1,j,1)+om(i,j-1,1) 1 +om(i,j+1,1)) else print*,'problem with om(i,j,1):',i,j,om(i,j,1) end if end if end do end do print*,'r=0 and om(i,j,1) reset at',ix,' pts' open(23,file='omega.dat',status='unknown') write(23,10) do i=1,ni if(rr(i).gt.0d0)then omega=a0(i)/(rr(i)*rgal) else if (i.eq.1)then omega=(2d0*a0(2)/rr(2)-0.5d0*a0(3)/rr(3))/1.5d0/rgal else print*,'error in omega(r)',i,rr(i-1),a0(i-1),rr(i),a0(i), 1 rr(i+1),a0(i+1) stop end if write(23,20)i,rr(i),a0(i),omega end do close(23) 10 format('r and fourier analysed vrot(r) in km/s', 1 ' om(r) in km/s/kpc') 20 format(i4,f8.4,1p,2e11.3) end c................................................................... subroutine seteta(eta0,gvx,gvy) use cdata use cparam implicit none real(kind=2)eta0(nx,ny),gvx(nx,ny),gvy(nx,ny),gm integer i,j,nk2,nk3,nx1,ny1,l gm=0d0 do j=1,ny do i=1,nx if(gvx(i,j)+gvy(i,j).gt.gm)gm=gvx(i,j)+gvy(i,j) end do end do do j=1,ny do i=1,nx eta0(i,j)=1d0+feta*(gvx(i,j)+gvy(i,j))/gm end do end do ny1=ny-2 nx1=nx-2 open(11,file='E2DAT',status='unknown') write(11,220) write(11,230)ny1,ymin,ymax write(11,240)nx1,xmin,xmax write(11,250) write(11,260) write(11,270) nk2=nx/3 nk3=2*nx/3 do 3 j=2,ny-1 write(11,*)(eta0(l,j),l=2,nk2) write(11,*)(eta0(l,j),l=nk2+1,nk3) write(11,*)(eta0(l,j),l=nk3+1,nx-1) 3 continue write(11,30) close(11) print*,'write to E2DAT done',feta call max1(eta0,'eta0') return 30 format('$ END') 220 format('$ DATA=CONTOUR') 230 format('%ny=',i3,' ymin=',f9.5,' ymax=',f9.5) 240 format('%nx=',i3,' xmin=',f9.5,' xmax=',f9.5) 250 format('%xlabel="x"') 260 format('%ylabel="y"') 270 format('%toplabel= "eta0 contours','"') end c.............................................................. subroutine gradv(gvx,gvy,vx0,vy0) use cdata use cparam implicit none integer i,j real(kind=2) vx0(nx,ny),vy0(nx,ny),gvx(nx,ny),gvy(nx,ny) do j=2,ny-1 do i=2,nx-1 if(rcp(i,j,1).lt.1d0)then gvx(i,j)=abs((vy0(i+1,j)-vy0(i-1,j))*rdx) gvy(i,j)=abs((vx0(i,j+1)-vx0(i,j-1))*rdy) else gvx(i,j)=0d0 gvy(i,j)=0d0 end if end do end do do j=1,ny gvx(1,j)=0d0 gvy(1,j)=0d0 gvx(nx,j)=0d0 gvy(nx,j)=0d0 end do do i=1,nx gvx(i,1)=0d0 gvy(i,1)=0d0 gvx(i,ny)=0d0 gvy(i,ny)=0d0 end do end c.................................................................... subroutine readden(rho) use cdata use cparam implicit none integer nx1,ny1,i,j,k,l,nk2,nk3 real(kind=2)dmax,dmin,denn,den(nx,ny),den0(nx-2,ny-2), 1 rho(nx,ny,nz) open(1,file='../ppe_den_100.data',status='old') read(1,*)den0 read(1,*)nx1,ny1 close(1) if(nx-2.ne.nx1.or.ny-2.ne.ny1)then print*,'inconsistent dimensions for density',nx,ny,nx1,ny1 stop end if do j=2,ny-1 do i=2,nx-1 den(i,j)=den0(i-1,j-1) end do end do do j=2,ny den(1,j)=den(nx-1,j) den(nx,j)=den(2,j) end do do i=2,nx-1 den(i,1)=den(i,ny-1) den(i,nx)=den(i,2) end do den(1,1)=(2d0*den(2,1)-den(3,1)+2d0*den(1,2)-den(1,3))/2d0 den(nx,ny)=(2d0*den(nx-1,1)-den(nx-2,1)+2d0*den(1,ny-1)- 1 den(1,ny-2))/2d0 den(1,ny)=(2d0*den(2,ny)-den(3,ny)+2d0*den(1,ny-1)- 1 den(1,ny-2))/2d0 den(nx,1)=(2d0*den(2,nx-1)-den(nx-2,1)+2d0*den(nx,2)- 1 den(nx,3))/2d0 dmax=-1d10 dmin=1d10 do j=1,nx do i=1,ny if(den(i,j).gt.dmax)dmax=den(i,j) if(den(i,j).lt.dmin)dmin=den(i,j) end do end do write(2,30)dmin,dmax write(6,30)dmin,dmax do j=1,ny do i=1,nx denn=den(i,j) den(i,j)=den(i,j)-dmax end do end do open(11,file='PDEN',status='unknown') write(11,220) write(11,230)ny1,ymin,ymax write(11,240)nx1,xmin,xmax write(11,250) write(11,260) write(11,270) nk2=nx/3 nk3=2*nx/3 do 3 j=2,ny-1 write(11,*)(den(l,j),l=2,nk2) write(11,*)(den(l,j),l=nk2+1,nk3) write(11,*)(den(l,j),l=nk3+1,nx-1) 3 continue write(11,130) close(11) do k=1,nz do j=1,ny do i=1,nx rho(i,j,k)=1d1**den(i,j) end do end do end do 10 format(2i4) 20 format(1p,8e10.3) 30 format('log min, max densities',1p,2e11.3,' normalize with max') c40 format('% xmin=',f8.3,' xmax=',f8.3) c50 format('% ymin=',f8.3,' ymax=',f8.3) c60 format('% zmin=',f8.3,' zmax=',f8.3) 130 format('$ END') 220 format('$ DATA=CONTOUR') 230 format('%ny=',i3,' ymin=',f9.5,' ymax=',f9.5) 240 format('%nx=',i3,' xmin=',f9.5,' xmax=',f9.5) 250 format('%xlabel="x"') 260 format('%ylabel="y"') 270 format('%toplabel= " ','"') end