c.................................................................... c program von3daniken c c 3D cartesian code c c includes: hard edge boundary conditions (reflecting x,y,z) c polytropic equation of state c.................................................................... include 'von_3defs.h' save c......initialization call param_in call data_in call timestep call set_gravity3d c......evolution cycle do istep=nstep,nsteps c........Perturbation...................... c if(istep.eq.1) then c call perturbation c end if call gravity3d call source call transport call interrupt call timestep end do call interrupt end subroutine param_in c......Read the input deck include 'von_3defs.h' save open (unit=1,file='parameters',status='old') c......istart: 0 = from scratch, or i = binary model number read(1,*) istrt c......Trun: expected run time read(1,*) trun c......outint: time interval for output read(1,*) outint c......starting step number (for positioning timestep file) read(1,*) nstart_step c......nsteps: Maximum ending step number read(1,*) nsteps c......cf: courant factor read(1,*) cf c......x1min: read(1,*) x1min c......x1max: read(1,*) x1max c......x2min: read(1,*) x2min c......x2max: read(1,*) x2max c......x3min: read(1,*) x3min c......x3max: read(1,*) x3max c......cv: coefficient for artvisc (negative if inactive) read(1,*) cv c......gamma: polytropic index read(1,*) gamma c......pK: polytropic constant read(1,*) pK c......ndim: 2=thin disk geometry, 3=thick disk geometry read(1,*) ndim c......eps: grav softening read(1,*) eps close(1) c......position 'time.record' file to the proper line open(unit=11,file='time.record',status='unknown') do i=1,nstart_step-1 read(11,*) end do return end subroutine source c......source terms in the momentum equation. include 'von_3defs.h' c......local viscosity arrays. dimension q1(-1:n1+2,-1:n2+2,-1:n3+2) dimension q2(-1:n1+2,-1:n2+2,-1:n3+2) dimension q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......Pressure from a polytropic equation of state do i=1,n1 do j=1,n2 do k=1,n3 P(i,j,k)=e(i,j,k)*(gamma-1) end do end do end do call x1bi(p,0,1,n1,n2,n3) call x1bo(p,0,1,n1,n2,n3) call x2bi(p,0,1,n1,n2,n3) call x2bo(p,0,1,n1,n2,n3) call x3bi(p,0,1,n1,n2,n3) call x3bo(p,0,1,n1,n2,n3) c......1 (x) direction force terms do i=1,n1 do j=1,n2 do k=1,n3 v1(i,j,k)=v1(i,j,k)+ + dt*((-(P(i,j,k)-P(i-1,j,k))/dx1b(i))/ + ((d(i,j,k)+d(i-1,j,k))/2.) + -(phi(i,j,k)-phi(i-1,j,k))/dx1b(i)) end do end do end do call x1bo(v1,1,-1,n1,n2,n3) call x1bi(v1,1,-1,n1,n2,n3) call x1oz(v1,n1,n2,n3) call x1iz(v1,n1,n2,n3) call x2bo(v1,0,1,n1,n2,n3) call x2bi(v1,0,1,n1,n2,n3) call x3bo(v1,0,1,n1,n2,n3) call x3bi(v1,0,1,n1,n2,n3) c...... (y) direction force terms do i=1,n1 do j=1,n2 do k=1,n3 v2(i,j,k)=v2(i,j,k)+ + dt*((-(P(i,j,k)-P(i,j-1,k))/dx2b(j))/ + ((d(i,j,k)+d(i,j-1,k))/2.) + -(phi(i,j,k)-phi(i,j-1,k))/dx2b(j)) end do end do end do call x2bo(v2,1,-1,n1,n2,n3) call x2bi(v2,1,-1,n1,n2,n3) call x2oz(v2,n1,n2,n3) call x2iz(v2,n1,n2,n3) call x1bo(v2,0,1,n1,n2,n3) call x1bi(v2,0,1,n1,n2,n3) call x3bo(v2,0,1,n1,n2,n3) call x3bi(v2,0,1,n1,n2,n3) c...... (z) direction force terms do i=1,n1 do j=1,n2 do k=1,n3 v3(i,j,k)=v3(i,j,k)+ + dt*((-(P(i,j,k)-P(i,j,k-1))/dx3b(k))/ + ((d(i,j,k)+d(i,j,k-1))/2.) + -(phi(i,j,k)-phi(i,j,k-1))/dx3b(k)) end do end do end do call x3bo(v3,1,-1,n1,n2,n3) call x3bi(v3,1,-1,n1,n2,n3) call x3oz(v3,n1,n2,n3) call x3iz(v3,n1,n2,n3) call x2bo(v3,0,1,n1,n2,n3) call x2bi(v3,0,1,n1,n2,n3) call x1bo(v3,0,1,n1,n2,n3) call x1bi(v3,0,1,n1,n2,n3) c......determine the artificial viscosity c......note that artificial viscosity corrections to the velocity c......are not computed along the grid boundaries c*.....in the 1D case (n2=n3=1) 'j' and 'k' should go from 1 to n2, c*.....and 1 to 'n3' respectively, otherwise the artificial viscosity c*.....is zero everywhere! if (cv.gt.0.) then do i=1,n1 do j=1,n2 do k=1,n3 if((v1(i+1,j,k)-v1(i,j,k)).lt.0.) then q1(i,j,k)=cv*d(i,j,k)*(v1(i+1,j,k)-v1(i,j,k))**2 else q1(i,j,k)=0. end if if((v2(i,j+1,k)-v2(i,j,k)).lt.0.) then q2(i,j,k)=cv*d(i,j,k)*(v2(i,j+1,k)-v2(i,j,k))**2 else q2(i,j,k)=0. end if if((v3(i,j,k+1)-v3(i,j,k)).lt.0.) then q3(i,j,k)=cv*d(i,j,k)*(v3(i,j,k+1)-v3(i,j,k))**2 else q3(i,j,k)=0. end if end do end do end do cc........apply viscous pressure to update the velocities c cc*.......here 'j' and 'k' has to go from 1,n2 and 1,n3 c*.......(or else don't have update) do i=1,n1 do j=1,n2 do k=1,n3 v1(i,j,k)=v1(i,j,k)-dt*(q1(i,j,k)-q1(i-1,j,k))/ + (dx1b(i)*((d(i,j,k)+d(i-1,j,k))/2.)) v2(i,j,k)=v2(i,j,k)-dt*(q2(i,j,k)-q2(i,j-1,k))/ + (dx2b(j)*((d(i,j,k)+d(i,j-1,k))/2.)) v3(i,j,k)=v3(i,j,k)-dt*(q3(i,j,k)-q3(i,j,k-1))/ + (dx3b(k)*((d(i,j,k)+d(i,j,k-1))/2.)) cc*.....update of energy e(i,j,k)=e(i,j,k)-dt*(q1(i,j,k)*(v1(i+1,j,k)- + v1(i,j,k))/dx1a(i)+q2(i,j,k)* + (v2(i,j+1,k)-v2(i,j,k))/dx2a(j)+ + q3(i,j,k)*(v3(i,j,k+1)-v3(i,j,k))/dx3a(k)) end do end do end do end if cc*....compressional heating (gamma low!!) do i=1,n1 do j=1,n2 do k=1,n3 e(i,j,k)=((1.-(dt/2.)*(gamma-1)*((v1(i+1,j,k)- + v1(i,j,k))/dx1a(i)+(v2(i,j+1,k)-v2(i,j,k))/dx2a(j)+ + (v3(i,j,k+1)-v3(i,j,k))/dx3a(k)))/ + (1.+(dt/2.)*(gamma-1)*((v1(i+1,j,k)-v1(i,j,k))/dx1a(i)+ + (v2(i,j+1,k)-v2(i,j,k))/dx2a(j)+(v3(i,j,k+1)- + v3(i,j,k))/dx3a(k))))* + e(i,j,k) end do end do end do return end subroutine transport c......1. density advection, 2. momentum advection, c......3. energy advection include 'von_3defs.h' c......transport related common variables: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......Transform velocities to momenta do i=1,n1 do j=1,n2 do k=1,n3 q1(i,j,k)=((d(i-1,j,k)+d(i,j,k))/2.)*v1(i,j,k) q2(i,j,k)=((d(i,j-1,k)+d(i,j,k))/2.)*v2(i,j,k) q3(i,j,k)=((d(i,j,k-1)+d(i,j,k))/2.)*v3(i,j,k) end do end do end do c......change ordering of advection calls every timestep.... if ((istep-(istep/3)*3).eq.1) then call rho_x_advect call energy_x_advect call xmoment_x_advect call ymoment_x_advect call zmoment_x_advect call rho_y_advect call energy_y_advect call xmoment_y_advect call ymoment_y_advect call zmoment_y_advect call rho_z_advect call energy_z_advect call xmoment_z_advect call ymoment_z_advect call zmoment_z_advect endif if ((istep-(istep/3)*3).eq.2) then call rho_y_advect call energy_y_advect call xmoment_y_advect call ymoment_y_advect call zmoment_y_advect call rho_z_advect call energy_z_advect call xmoment_z_advect call ymoment_z_advect call zmoment_z_advect call rho_x_advect call energy_x_advect call xmoment_x_advect call ymoment_x_advect call zmoment_x_advect endif if ((istep-(istep/3)*3).eq.0) then call rho_z_advect call energy_z_advect call xmoment_z_advect call ymoment_z_advect call zmoment_z_advect call rho_x_advect call energy_x_advect call xmoment_x_advect call ymoment_x_advect call zmoment_x_advect call rho_y_advect call energy_y_advect call xmoment_y_advect call ymoment_y_advect call zmoment_y_advect endif c......transform momenta back to velocities c*.....avoid division by zero do i=1,n1 do j=1,n2 do k=1,n3 if(((d(i-1,j,k)+d(i,j,k))/2.).eq.0.) then v1(i,j,k)=0. else v1(i,j,k)=q1(i,j,k)/((d(i-1,j,k)+d(i,j,k))/2.) endif if(((d(i,j-1,k)+d(i,j,k))/2.).eq.0.) then v2(i,j,k)=0. else v2(i,j,k)=q2(i,j,k)/((d(i,j-1,k)+d(i,j,k))/2.) endif if(((d(i,j,k-1)+d(i,j,k))/2.).eq.0.) then v3(i,j,k)=0. else v3(i,j,k)=q3(i,j,k)/((d(i,j,k-1)+d(i,j,k))/2.) endif end do end do end do call x1bo(v1,1,-1,n1,n2,n3) call x1bi(v1,1,-1,n1,n2,n3) call x1oz(v1,n1,n2,n3) call x1iz(v1,n1,n2,n3) call x2bo(v1,0,1,n1,n2,n3) call x2bi(v1,0,1,n1,n2,n3) call x3bo(v1,0,1,n1,n2,n3) call x3bi(v1,0,1,n1,n2,n3) call x2bo(v2,1,-1,n1,n2,n3) call x2bi(v2,1,-1,n1,n2,n3) call x2oz(v2,n1,n2,n3) call x2iz(v2,n1,n2,n3) call x1bo(v2,0,1,n1,n2,n3) call x1bi(v2,0,1,n1,n2,n3) call x3bo(v2,0,1,n1,n2,n3) call x3bi(v2,0,1,n1,n2,n3) call x3bo(v3,1,-1,n1,n2,n3) call x3bi(v3,1,-1,n1,n2,n3) call x3oz(v3,n1,n2,n3) call x3iz(v3,n1,n2,n3) call x2bo(v3,0,1,n1,n2,n3) call x2bi(v3,0,1,n1,n2,n3) call x1bo(v3,0,1,n1,n2,n3) call x1bi(v3,0,1,n1,n2,n3) return end c................................................................. c subroutine grid_set c c................................................................ include 'von_3defs.h' save c......define grid radii c......set up x zones (evenly spaced) dx=(x1max-x1min)/real(n1) do i=-1,n1+2 x1a(i)=x1min+dx*real(i-1) x1b(i)=x1a(i)+dx/2. dx1a(i)=dx dx1b(i)=dx end do c......set up y zones (evenly spaced) dy=(x2max-x2min)/real(n2) do i=-1,n2+2 x2a(i)=x2min+dy*real(i-1) x2b(i)=x2a(i)+dy/2. dx2a(i)=dy dx2b(i)=dy end do c......set up z zones (evenly spaced) dz=(x3max-x3min)/real(n3) do i=-1,n3+2 x3a(i)=x3min+dz*real(i-1) x3b(i)=x3a(i)+dz/2. dx3a(i)=dz dx3b(i)=dz end do c......compute volume element array, distance from center array c c......With n2 ** 3 equal cells, dv should not be an array !! c do i=1,n1 do j=1,n2 do k=1,n3 dv(i,j,k)=dx1a(i)*dx2a(j)*dx3a(k) end do end do end do return end c.................................................................. c subroutine boundary c c.................................................................. include 'von_3defs.h' save call x1bi(d,0,1,n1,n2,n3) call x1bo(d,0,1,n1,n2,n3) call x2bi(d,0,1,n1,n2,n3) call x2bo(d,0,1,n1,n2,n3) c*.....for the 1D test, the inner boundary condition is called c*.... twice. (Have to find a better solution for this..) call x2bi(d,0,1,n1,n2,n3) call x3bo(d,0,1,n1,n2,n3) call x1bi(p,0,1,n1,n2,n3) call x1bo(p,0,1,n1,n2,n3) call x2bi(p,0,1,n1,n2,n3) call x2bo(p,0,1,n1,n2,n3) call x3bi(p,0,1,n1,n2,n3) call x3bo(p,0,1,n1,n2,n3) call x1bi(e,0,1,n1,n2,n3) call x1bo(e,0,1,n1,n2,n3) call x2bi(e,0,1,n1,n2,n3) call x2bo(e,0,1,n1,n2,n3) call x3bi(e,0,1,n1,n2,n3) call x3bo(e,0,1,n1,n2,n3) call x1bo(v1,1,-1,n1,n2,n3) call x1bi(v1,1,-1,n1,n2,n3) call x1oz(v1,n1,n2,n3) call x1iz(v1,n1,n2,n3) call x2bo(v1,0,1,n1,n2,n3) call x2bi(v1,0,1,n1,n2,n3) call x3bo(v1,0,1,n1,n2,n3) call x3bi(v1,0,1,n1,n2,n3) call x2bo(v2,1,-1,n1,n2,n3) call x2bi(v2,1,-1,n1,n2,n3) call x2oz(v2,n1,n2,n3) call x2iz(v2,n1,n2,n3) call x1bo(v2,0,1,n1,n2,n3) call x1bi(v2,0,1,n1,n2,n3) call x3bo(v2,0,1,n1,n2,n3) call x3bi(v2,0,1,n1,n2,n3) call x3bo(v3,1,-1,n1,n2,n3) call x3bi(v3,1,-1,n1,n2,n3) call x3oz(v3,n1,n2,n3) call x3iz(v3,n1,n2,n3) call x2bo(v3,0,1,n1,n2,n3) call x2bi(v3,0,1,n1,n2,n3) call x1bo(v3,0,1,n1,n2,n3) call x1bi(v3,0,1,n1,n2,n3) return end c...................................................................... c subroutine timestep c c..................................................................... include 'von_3defs.h' save c......Establish timestep criteria dt=0 dt1=1.e+20 c*.....in the original version c*.....if (cv.gt.0.) then if (cv.lt.0.) then c.. no artificial viscosity in timestep limiter do i=1,n1 do j=1,n2 do k=1,n3 csound=sqrt(gamma*p(i,j,k)/d(i,j,k)) dt1=dx1a(i)/csound sample=dx2a(j)/csound if (sample.lt.dt1) dt1=sample sample=dx3a(k)/csound if(sample.lt.dt1) dt1=sample dt2=dx1a(i)/(max(abs(v1(i,j,k)),1.d-16)) dt3=dx2a(j)/(max(abs(v2(i,j,k)),1.d-16)) dt4=dx3a(k)/(max(abs(v3(i,j,k)),1.d-16)) sample=(1./dt1**2+1./dt2**2+1./dt3**2+1./dt4**2) c............. Find which zone limits timestep: if(sample.gt.dt) then dt=sample iflag=i jflag=j kflag=k end if c.............................................. end do end do end do dt=cf/sqrt(dt) else c.. include artificial viscosity in the timestep limiter do i=1,n1 do j=1,n2 do k=1,n3 if (d(i,j,k).lt.0.) then write(11,*) 'Negative density zone:',i,j,k stop end if csound=sqrt(gamma*p(i,j,k)/d(i,j,k)) dt1=dx1a(i)/csound sample=dx2a(j)/csound if(sample.lt.dt1) dt1=sample sample=dx3a(k)/csound if(sample.lt.dt1) dt1=sample dt2=dx1a(i)/(max(abs(v1(i,j,k)),1.d-16)) dt3=dx2a(j)/(max(abs(v2(i,j,k)),1.d-16)) dt4=dx3a(k)/(max(abs(v3(i,j,k)),1.d-16)) dt5=dx1a(i)/(4*cv*(max(abs(v1(i+1,j,k) + -v1(i,j,k)),1.d-16))) sample=dx2a(j)/(4*cv*(max(abs(v2(i,j+1,k) + -v2(i,j,k)),1.d-16))) if(sample.lt.dt5) dt5=sample sample=dx3a(k)/(4*cv*(max(abs(v3(i,j,k+1) + -v3(i,j,k)),1.d-16))) if(sample.lt.dt5) dt5=sample sample=(1./dt1**2+1./dt2**2+1./dt3**2 + +1./dt4**2+1./dt5**2) if(sample.gt.dt) then dt=sample iflag=i jflag=j kflag=k end if end do end do end do dt=cf/sqrt(dt) end if if(istep.ge.1) then time=time+dt end if c......timestep file update on each timestep.......... write(11,2110) istep,time,dt,iflag,jflag,kflag 2110 format(i6,1p2e12.4,3i4) c.............................................. return end c............................................................... c subroutine data_in c c.............................................................. include 'von_3defs.h' save open(unit=2,file='image',status='unknown',form='unformatted') if(istrt.eq.0) then call grid_set call set_ic else do i=1,istrt read(2) x1a,x2a,x3a, + x1b,x2b,x3b, + dx1a,dx2a,dx3a, + dx1b,dx2b,dx3b, + dv, + d,v1,v2,v3, + time,dt,istart,nstep end do time_restart=time istart_restart=istart end if return end c................................................................... c subroutine set_ic c c................................................................... include 'von_3defs.h' save c......routine to set up variables with their initial values c*.....set initial values for the energy and pressure nstep=1 time=0. do i=-1,n1+2 do j=-1,n2+2 do k=-1,n3+2 v1(i,j,k)=0. v2(i,j,k)=0. v3(i,j,k)=0. d(i,j,k)= 0. e(i,j,k)=0. p(i,j,k)=0. phi(i,j,k)=0. Plm(i,j,k)=0. end do end do end do c*....advection test initial conditions c do j=1,n2 c do i=1,n1 c do k=1,n3 c if (i.ge.5.and.i.le.54) then c v1(i,j,k)=1. c v2(i,j,k)=0. c v3(i,j,k)=0. c d(i,j,k)=1. c e(i,j,k)=1. c p(i,j,k)=e(i,j,k)*(gamma-1) c else c v1(i,j,k)=0. c v2(i,j,k)=0. c v3(i,j,k)=1. c d(i,j,k)=0. c e(i,j,k)=0. c p(i,j,k)=e(i,j,k)*(gamma-1) c end if c end do c end do c end do c c*.....Sod shock tube test initial conditions c do i=1,n1 c do j=1,n2 c do k=1,n3 c if (k.ge.1.and.k.le.50) then c v1(i,j,k)=0. c v2(i,j,k)=0. c v3(i,j,k)=0. c d(i,j,k)=1. c e(i,j,k)=1./(gamma-1) c p(i,j,k)=e(i,j,k)*(gamma -1) c else c v1(i,j,k)=0. c v2(i,j,k)=0. c v3(i,j,k)=0. c d(i,j,k)=0.125 c e(i,j,k)=0.1/(gamma-1) c p(i,j,k)=e(i,j,k)*(gamma -1) c end if c end do c end do c end do c....potential solver test initial cond-point charge c d(n1/2,n2/2,n3/2)=1. c d(10,10,10)=0.7 c d(2,3,3)=0.5 c d(20,12,11)=1.2 c d(30,32,5)=0.4 c d(2,3,5)=0.3 c d(3,3,4)=0.6 c d(1,10,15)=1.1 c d(15,15,17)=0.1 c d(7,7,1)=1. c write(*,*) d(n1/2,n2/2,n3/2) c.......write(*,*) (n1/2),(n2/2),(n3/2),d(n1/2,n2/2,n3/2) c do i=1,n1 c do j=1,n2 c do k=1,n3 c Plm(i,j,k)=(-1.*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(n1/2)))**2+ c + (dx*(real(j)-real(n2/2)))**2+(dx*(real(k)-real(n3/2)))**2)+ c + (-0.7*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(10)))**2+ c + (dx*(real(j)-real(10)))**2+(dx*(real(k)-real(10)))**2)+ c + (-0.5*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(2)))**2+ c + (dx*(real(j)-real(3)))**2+(dx*(real(k)-real(3)))**2)+ c + (-1.2*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(20)))**2+ c + (dx*(real(j)-real(12)))**2+(dx*(real(k)-real(11)))**2)+ c + (-0.4*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(30)))**2+ c + (dx*(real(j)-real(32)))**2+(dx*(real(k)-real(5)))**2)+ c + (-0.3*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(2)))**2+ c + (dx*(real(j)-real(3)))**2+(dx*(real(k)-real(5)))**2)+ c + (-0.6*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(3)))**2+ c + (dx*(real(j)-real(3)))**2+(dx*(real(k)-real(4)))**2)+ c + (-1.1*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(1)))**2+ c + (dx*(real(j)-real(10)))**2+(dx*(real(k)-real(15)))**2)+ c + (-0.1*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(15)))**2+ c + (dx*(real(j)-real(15)))**2+(dx*(real(k)-real(17)))**2)+ c + (-1*dx*dx*dx)/sqrt c + (eps**2+(dx*(real(i)-real(7)))**2+ c + (dx*(real(j)-real(7)))**2+(dx*(real(k)-real(1)))**2) c c end do c end do c end do c cc......Full call to the boundary conditions call boundary call interrupt return end c......................................................................... c subroutine interrupt c c......................................................................... include 'von_3defs.h' c......local character variables for file manipulation character*3 sstring character*8 filestem character*10 nstring character*12 target_out DATA nstring/'0123456789'/ save c......start counting if (istep.eq.1.or.istep.eq.2) then istart=1 endif nsnap=istart c..... Binary dump if (time.gt.trun) then write(2) x1a,x2a,x3a, + x1b,x2b,x3b, + dx1a,dx2a,dx3a, + dx1b,dx2b,dx3b, + dv, + d,v1,v2,v3, + time,dt,istart,nstep end if c.......Ascii file dump. Stored in File den.xxx if(time.eq.0. .or. (time-time_restart).gt. + (istart-istart_restart)*outint) then c if ((time-time_restart).gt.(istart-istart_restart)*outint) then sstring(1:1)=nstring(1+nsnap/100:1+nsnap/100) istring=1+MOD(nsnap,100)/10 sstring(2:2)=nstring(istring:istring) istring=1+MOD(nsnap,10) sstring(3:3)=nstring(istring:istring) filestem='pse.' target_out=filestem(1:4)//sstring(1:3) istart=istart+1 open (unit=3,file=target_out,status='unknown') write(3,*) time do i=1,n1 write(3,2109) x1a(i) end do do j=1,n2 write(3,2109) x2a(j) end do do k=1,n3 write(3,2109) x3a(k) end do c..write out a 'slice' (n1/2,j,n3/2) of the errors c do i=1,n1 c do j=1,n2 c do k=1,n3 c write(3,2109) abs(phi(n1/2,j,n3/2)- Plm(n1/2,j,n3/2)), c + abs((phi(n1/2,j,n3/2)- Plm(n1/2,j,n3/2))/ c + Plm(n1/2,j,n3/2)) c end do c end do c write(3,*) '' c end do c*........write out energies do i=1,n1 do j=1,n2 do k=1,n3 write(3,2109) e(i,j,k) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 write(3,2109) v1(i,j,k) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 write(3,2109) v2(i,j,k) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 write(3,2109) v3(i,j,k) end do end do end do close(3) end if c2109 format(1p1e12.4) c*......new format to handle very low densities c2109 format(1p1e12.4e3) 2109 format(1p1e26.13e3,1p1e26.13e3) 2110 format(9(1p1e14.6)) c.......Mass sum checker c sum=0. c do i=1,n1 c do j=1,n2 c do k=1,n3 c sum=sum+d(i,j,k)*dv(i,j,k) c end do c end do c end do c write(*,*) sum c....................................... if (time.gt.trun) then write(*,*) 'End of Simulation' stop end if return end c....................................................................... double precision FUNCTION RAN2(IDUM) c....................................................................... c.....random number generator (from numerical recipies) implicit real*8(a-h,o-z), integer*4(i-n) PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6) DIMENSION IR(97) DATA IFF /0/ save IF(IDUM.LT.0.OR.IFF.EQ.0)THEN IFF=1 IDUM=MOD(IC-IDUM,M) DO 11 J=1,97 IDUM=MOD(IA*IDUM+IC,M) IR(J)=IDUM 11 CONTINUE IDUM=MOD(IA*IDUM+IC,M) IY=IDUM ENDIF J=1+(97*IY)/M IF(J.GT.97.OR.J.LT.1)PAUSE IY=IR(J) RAN2=IY*RM IDUM=MOD(IA*IDUM+IC,M) IR(J)=IDUM RETURN end c................................................................ c subroutine perturbation c c................................................................ c......Supplies a perturbation to the density include 'von_3defs.h' save do i=1,1 do j=1,n2 do k=1,n3 c........... random perturbation............................. c d(i,j,k)=d(i,j,k)*(1.+.001*(2*ran2(idum)-1)) c............................................................... end do end do end do c.........Full call to boundary condition. call boundary return end c............................................................... c subroutine rho_x_advect c c............................................................... c......density advection in x direction include 'von_3defs.h' c......Transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......set up interpolated densities at x zone interfaces do i=1,n1+1 do j=1,n2 do k=1,n3 if(v1(i,j,k).gt.0.) then dqpi12=(d(i,j,k)-d(i-1,j,k))/dx1b(i) dqmi12=(d(i-1,j,k)-d(i-2,j,k))/dx1b(i-1) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds1(i,j,k)=d(i-1,j,k)+(dx1a(i-1)-v1(i,j,k)*dt)*dqi else dqpi12=(d(i+1,j,k)-d(i,j,k))/dx1b(i+1) dqmi12=(d(i,j,k)-d(i-1,j,k))/dx1b(i) if((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds1(i,j,k)=d(i,j,k)-(dx1a(i)+v1(i,j,k)*dt)*dqi end if end do end do end do c......Face centered fluxes in the x direction do i=1,n1+1 do j=1,n2 do k=1,n3 ds1(i,j,k)=ds1(i,j,k)*v1(i,j,k) Fm1(i,j,k)=ds1(i,j,k)*dx3a(k)*dx2a(j) end do end do end do call x2bi(ds1,0,1,n1,n2,n3) call x2bi(fm1,0,1,n1,n2,n3) call x2bo(ds1,0,1,n1,n2,n3) call x2bo(fm1,0,1,n1,n2,n3) call x1bi(ds1,1,1,n1,n2,n3) call x1bi(fm1,1,-1,n1,n2,n3) call x1iz(fm1,n1,n2,n3) call x1bo(ds1,1,1,n1,n2,n3) call x1bo(fm1,1,-1,n1,n2,n3) call x1oz(fm1,n1,n2,n3) call x3bi(ds1,0,1,n1,n2,n3) call x3bi(fm1,0,1,n1,n2,n3) call x3bo(ds1,0,1,n1,n2,n3) call x3bo(fm1,0,1,n1,n2,n3) c......update the density do i=1,n1 do j=1,n2 do k=1,n3 d(i,j,k)=d(i,j,k)-(dt/dv(i,j,k))*(Fm1(i+1,j,k)-Fm1(i,j,k)) end do end do end do call x1bi(d,0,1,n1,n2,n3) call x1bo(d,0,1,n1,n2,n3) call x2bi(d,0,1,n1,n2,n3) call x2bo(d,0,1,n1,n2,n3) call x3bi(d,0,1,n1,n2,n3) call x3bo(d,0,1,n1,n2,n3) return end c..................................................................... c subroutine rho_y_advect c c..................................................................... c......density advection in r direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......density advection (2 -- (r) direction) c......set up interpolated densities at r zone interfaces do i=1,n1 do j=1,n2+1 do k=1,n3 if (v2(i,j,k).gt.0.) then dqpi12=(d(i,j,k)-d(i,j-1,k))/dx2b(j) dqmi12=(d(i,j-1,k)-d(i,j-2,k))/dx2b(j-1) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds2(i,j,k)=d(i,j-1,k)+(dx2a(j-1)-v2(i,j,k)*dt)*dqi else dqpi12=(d(i,j+1,k)-d(i,j,k))/dx2b(j+1) dqmi12=(d(i,j,k)-d(i,j-1,k))/dx2b(j) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds2(i,j,k)=d(i,j,k)-(dx2a(j)+v2(i,j,k)*dt)*dqi end if end do end do end do c......Face centered fluxes in the R direction do i=1,n1 do j=1,n2+1 do k=1,n3 ds2(i,j,k)=ds2(i,j,k)*v2(i,j,k) Fm2(i,j,k)=ds2(i,j,k)*dx1a(i)*dx3a(k) end do end do end do call x1bi(ds2,0,1,n1,n2,n3) call x1bi(fm2,0,1,n1,n2,n3) call x1bo(ds2,0,1,n1,n2,n3) call x1bo(fm2,0,1,n1,n2,n3) call x2bi(ds2,1,1,n1,n2,n3) call x2bi(fm2,1,-1,n1,n2,n3) call x2iz(fm2,n1,n2,n3) call x2bo(ds2,1,1,n1,n2,n3) call x2bo(fm2,1,-1,n1,n2,n3) call x2oz(fm2,n1,n2,n3) call x3bi(ds2,0,1,n1,n2,n3) call x3bi(fm2,0,1,n1,n2,n3) call x3bo(ds2,0,1,n1,n2,n3) call x3bo(fm2,0,1,n1,n2,n3) c......Update the density do i=1,n1 do j=1,n2 do k=1,n3 d(i,j,k)=d(i,j,k)-(dt/dv(i,j,k))*(Fm2(i,j+1,k)-Fm2(i,j,k)) end do end do end do call x1bi(d,0,1,n1,n2,n3) call x1bo(d,0,1,n1,n2,n3) call x2bi(d,0,1,n1,n2,n3) call x2bo(d,0,1,n1,n2,n3) call x3bi(d,0,1,n1,n2,n3) call x3bo(d,0,1,n1,n2,n3) return end c..................................................................... c subroutine rho_z_advect c c..................................................................... c......density advection in z direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c.. set up interpolated densities on the z zone faces do i=1,n1 do j=1,n2 do k=1,n3+1 if(v3(i,j,k).gt.0.) then dqpi12=(d(i,j,k)-d(i,j,k-1))/(dx3b(k)) dqmi12=(d(i,j,k-1)-d(i,j,k-2))/(dx3b(k-1)) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds3(i,j,k)=d(i,j,k-1)+(dx3a(k-1)-v3(i,j,k)*dt) + *dqi else dqpi12=(d(i,j,k+1)-d(i,j,k))/(dx3b(k+1)) dqmi12=(d(i,j,k)-d(i,j,k-1))/(dx3b(k)) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. ds3(i,j,k)=d(i,j,k)-(dx3a(k)+v3(i,j,k)*dt) + *dqi end if end do end do end do c......Compute mass flux densities (by re-labeling ds1(i,j,k) etc.) c......Compute face centered fluxes in the phi direction do i=1,n1 do j=1,n2 do k=1,n3+1 ds3(i,j,k)=ds3(i,j,k)*v3(i,j,k) Fm3(i,j,k)=ds3(i,j,k)*dx1a(i)*dx2a(j) end do end do end do call x2bi(ds3,0,1,n1,n2,n3) call x2bi(fm3,0,1,n1,n2,n3) call x2bo(ds3,0,1,n1,n2,n3) call x2bo(fm3,0,1,n1,n2,n3) call x3bi(ds3,1,1,n1,n2,n3) call x3bi(fm3,1,-1,n1,n2,n3) call x3iz(fm3,n1,n2,n3) call x3bo(ds3,1,1,n1,n2,n3) call x3bo(fm3,1,-1,n1,n2,n3) call x3oz(fm3,n1,n2,n3) call x1bi(ds3,0,1,n1,n2,n3) call x1bi(fm3,0,1,n1,n2,n3) call x1bo(ds3,0,1,n1,n2,n3) call x1bo(fm3,0,1,n1,n2,n3) c......Update the density do i=1,n1 do j=1,n2 do k=1,n3 d(i,j,k)=d(i,j,k)-(dt/dv(i,j,k))*(Fm3(i,j,k+1)-Fm3(i,j,k)) end do end do end do call x1bi(d,0,1,n1,n2,n3) call x1bo(d,0,1,n1,n2,n3) call x2bi(d,0,1,n1,n2,n3) call x2bo(d,0,1,n1,n2,n3) call x3bi(d,0,1,n1,n2,n3) call x3bo(d,0,1,n1,n2,n3) return end c........................................................................ c subroutine xmoment_x_advect c c........................................................................ c......x momentum advection in the x direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......Extrapolate velocity to control volume interfaces do i=0,n1 do j=1,n2 do k=1,n3 if (((Fm1(i,j,k)+Fm1(i+1,j,k))/2.).gt.0.) then dqmi12=(v1(i,j,k)-v1(i-1,j,k))/dx1a(i-1) dqpi12=(v1(i+1,j,k)-v1(i,j,k))/dx1a(i) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c........donor cell c dqi=0. c.................. vs(i,j,k)=v1(i,j,k)+(dx1b(i)-((v1(i,j,k)+v1(i+1,j,k))/2.) + *dt)*dqi else dqmi12=(v1(i+1,j,k)-v1(i,j,k))/dx1a(i) dqpi12=(v1(i+2,j,k)-v1(i+1,j,k))/dx1a(i+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k)=v1(i+1,j,k)-(dx1b(i+1)+((v1(i,j,k)+v1(i+1,j,k))/2.) + *dt)*dqi end if end do end do end do call x1bi(vs,1,-1,n1,n2,n3) call x1iz(vs,n1,n2,n3) call x1bo(vs,1,-1,n1,n2,n3) call x1oz(vs,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) c......construct flux array for x momentum advection in x direction do i=1,n1 do j=1,n2 do k=1,n3 Area=dx3a(k)*dx2a(j) F(i,j,k)=vs(i,j,k)*(.5*(ds1(i,j,k)+ds1(i+1,j,k)))*Area if(i.eq.1) F(0,j,k)=F(1,j,k) dvg(i,j,k)=dx1b(i)*Area end do end do end do c......adjust x momentum due to flux in x direction do i=1,n1 do j=1,n2 do k=1,n3 q1(i,j,k)=q1(i,j,k)-(dt/dvg(i,j,k))*(F(i,j,k)-F(i-1,j,k)) end do end do end do return end c....................................................................... c subroutine ymoment_x_advect c c....................................................................... c......y momentum advection in the x direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......extrapolate r velocity to the control volume interface do i=0,n1 do j=1,n2 do k=1,n3 if ((Fm1(i+1,j,k)+Fm1(i+1,j-1,k)).gt.0.) then dqmi12=(v2(i,j,k)-v2(i-1,j,k))/dx1b(i) dqpi12=(v2(i+1,j,k)-v2(i,j,k))/dx1b(i+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i+1,j,k)=v2(i,j,k)+(dx1a(i)-((v1(i+1,j,k)+ + v1(i+1,j-1,k))/2.)*dt)*dqi else dqmi12=(v2(i+1,j,k)-v2(i,j,k))/dx1b(i+1) dqpi12=(v2(i+2,j,k)-v2(i+1,j,k))/dx1b(i+2) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i+1,j,k)=v2(i+1,j,k)-(dx1a(i+1)+((v1(i+1,j,k)+ + v1(i+1,j-1,k))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,1,1,n1,n2,n3) call x1iz(vs,n1,n2,n3) call x1bo(vs,1,1,n1,n2,n3) call x1oz(vs,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) c......construct flux array for y momentum advection in the x direction do i=1,n1+1 do j=1,n2 do k=1,n3 Area=dx3a(k)*dx2b(j) F(i,j,k)=vs(i,j,k)*.5*(ds1(i,j,k)+ds1(i,j-1,k))*Area dvg(i,j,k)=Area*dx1a(i) end do end do end do c......Update specific y momentum due to advection in the x direction do i=1,n1 do j=1,n2 do k=1,n3 q2(i,j,k)=q2(i,j,k)-(dt/dvg(i,j,k))*(F(i+1,j,k)-F(i,j,k)) end do end do end do return end c.................................................................. c subroutine zmoment_x_advect c c.................................................................. c......z momentum advection in the x direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......z momentum advection in the x direction c......extrapolate z velocity to the control volume interface: do i=0,n1 do j=1,n2 do k=1,n3 if((Fm1(i+1,j,k)+Fm1(i+1,j,k-1)).gt.0.) then dqmi12=(v3(i,j,k)-v3(i-1,j,k))/dx1b(i) dqpi12=(v3(i+1,j,k)-v3(i,j,k))/dx1b(i+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i+1,j,k)=v3(i,j,k)+(dx1a(i)-((v1(i+1,j,k)+ + v1(i+1,j,k-1))/2.)*dt)*dqi else dqmi12=(v3(i+1,j,k)-v3(i,j,k))/dx1b(i+1) dqpi12=(v3(i+2,j,k)-v3(i+1,j,k))/dx1b(i+2) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i+1,j,k)=v3(i+1,j,k)-(dx1a(i+1)+((v1(i+1,j,k) + +v1(i+1,j,k-1))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,1,1,n1,n2,n3) call x1iz(vs,n1,n2,n3) call x1bo(vs,1,1,n1,n2,n3) call x1oz(vs,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) do i=1,n1+1 do j=1,n2 do k=1,n3 Area=dx3b(k)*dx2a(j) F(i,j,k)=vs(i,j,k)*.5*(ds1(i,j,k)+ds1(i,j,k-1))*area dvg(i,j,k)=Area*dx1a(i) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q3(i,j,k)=q3(i,j,k)-(dt/dvg(i,j,k))* + (F(i+1,j,k)-F(i,j,k)) end do end do end do return end c..................................................................... c subroutine xmoment_y_advect c c.................................................................... c......x momentum advection in the y direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......extrapolate x velocity to the control volume interface do i=1,n1 do j=0,n2 do k=1,n3 if((Fm2(i,j+1,k)+Fm2(i-1,j+1,k)).gt.0.) then dqmi12=(v1(i,j,k)-v1(i,j-1,k))/dx2b(j) dqpi12=(v1(i,j+1,k)-v1(i,j,k))/dx2b(j+1) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j+1,k)=v1(i,j,k)+(dx2a(j)-((v2(i,j+1,k)+ + v2(i-1,j+1,k))/2.)*dt)*dqi else dqmi12=(v1(i,j+1,k)-v1(i,j,k))/dx2b(j+1) dqpi12=(v1(i,j+2,k)-v1(i,j+1,k))/dx2b(j+2) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j+1,k)=v1(i,j+1,k)-(dx2a(j+1)+((v2(i,j+1,k)+ + v2(i-1,j+1,k))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,1,1,n1,n2,n3) call x2iz(vs,n1,n2,n3) call x2bo(vs,1,1,n1,n2,n3) call x2oz(vs,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) c......Construct flux array for x momentum advection in y direction do i=1,n1 do j=1,n2+1 do k=1,n3 Area=dx3a(k)*dx1b(i) F(i,j,k)=vs(i,j,k)*.5*(ds2(i,j,k)+ds2(i-1,j,k))*Area dvg(i,j,k)=dx3a(k)*dx1b(i)*dx2a(j) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q1(i,j,k)=q1(i,j,k)-(dt/dvg(i,j,k))*(F(i,j+1,k)-F(i,j,k)) end do end do end do return end c................................................................. c subroutine ymoment_y_advect c c................................................................. c......y momentum advection in the y direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c... Extrapolate y velocity to y control volume surface do i=1,n1 do j=0,n2 do k=1,n3 if((Fm2(i,j,k)+Fm2(i,j+1,k)).gt.0.) then dqmi12=(v2(i,j,k)-v2(i,j-1,k))/dx2a(j-1) dqpi12=(v2(i,j+1,k)-v2(i,j,k))/dx2a(j) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k)=v2(i,j,k)+(dx2b(j)-((v2(i,j,k)+v2(i,j+1,k)) + /2.)*dt)*dqi else dqmi12=(v2(i,j+1,k)-v2(i,j,k))/dx2a(j) dqpi12=(v2(i,j+2,k)-v2(i,j+1,k))/dx2a(j+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. endif c.........donor cell c dqi=0. c.................. vs(i,j,k)=v2(i,j+1,k)-(dx2b(j+1)+((v2(i,j,k) + +v2(i,j+1,k))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,1,-1,n1,n2,n3) call x2iz(vs,n1,n2,n3) call x2bo(vs,1,-1,n1,n2,n3) call x2oz(vs,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) c.......construct flux array for y momentum advection in y direction do i=1,n1 do j=1,n2 do k=1,n3 Area=dx3a(k)*dx1a(i) F(i,j,k)=vs(i,j,k)*(.5*(ds2(i,j,k)+ds2(i,j+1,k)))*Area if(j.eq.1) F(i,0,k)=F(i,1,k) dvg(i,j,k)=dx2b(j)*dx3a(k)*dx1a(i) end do end do end do c.......Update specific y momentum due to y momentum advection do i=1,n1 do j=1,n2 do k=1,n3 q2(i,j,k)=q2(i,j,k)-(dt/dvg(i,j,k))*(F(i,j,k)-F(i,j-1,k)) end do end do end do return end c....................................................................... c subroutine zmoment_y_advect c c....................................................................... c......z momentum advection in the y direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c........z momentum in the y direction c........extrapolate z velocity to control volume interface do i=1,n1 do j=0,n2 do k=1,n3 if((Fm2(i,j+1,k)+Fm2(i,j+1,k-1)).gt.0.) then dqmi12=(v3(i,j,k)-v3(i,j-1,k))/dx2b(j) dqpi12=(v3(i,j+1,k)-v3(i,j,k))/dx2b(j+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j+1,k)=v3(i,j,k)+(dx2a(j)-((v2(i,j+1,k)+ + v2(i,j+1,k-1))/2.)*dt)*dqi else dqmi12=(v3(i,j+1,k)-v3(i,j,k))/dx2b(j+1) dqpi12=(v3(i,j+2,k)-v3(i,j+1,k))/dx2b(j+2) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j+1,k)=v3(i,j+1,k)-(dx2a(j+1)+ + ((v2(i,j+1,k)+v2(i,j+1,k-1))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,1,1,n1,n2,n3) call x2iz(vs,n1,n2,n3) call x2bo(vs,1,1,n1,n2,n3) call x2oz(vs,n1,n2,n3) call x3bi(vs,0,1,n1,n2,n3) call x3bo(vs,0,1,n1,n2,n3) c........construct flux array for z momentum advection in y direction do i=1,n1 do j=1,n2+1 do k=1,n3 Area=dx1a(i)*dx3b(k) F(i,j,k)=vs(i,j,k)*.5*(ds2(i,j,k)+ds2(i,j,k-1))*Area dvg(i,j,k)=dx1a(i)*dx3b(k)*dx2a(j) end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q3(i,j,k)=q3(i,j,k)-(dt/dvg(i,j,k))*(F(i,j+1,k)-F(i,j,k)) end do end do end do return end c........................................................................ c subroutine xmoment_z_advect c c........................................................................ c......x momentum advection in the z direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save do i=1,n1 do j=1,n2 do k=0,n3 if((Fm3(i,j,k+1)+Fm3(i-1,j,k+1)).gt.0.) then dqmi12=(v1(i,j,k)-v1(i,j,k-1))/dx3b(k) dqpi12=(v1(i,j,k+1)-v1(i,j,k))/dx3b(k+1) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k+1)=v1(i,j,k)+(dx3a(k)-((v3(i,j,k+1) + +v3(i-1,j,k+1))/2.)*dt)*dqi else dqmi12=(v1(i,j,k+1)-v1(i,j,k))/dx3b(k+1) dqpi12=(v1(i,j,k+2)-v1(i,j,k+1))/dx3b(k+2) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k+1)=v1(i,j,k+1)-(dx3a(k+1)+((v3(i,j,k+1) + +v3(i-1,j,k+1))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,1,1,n1,n2,n3) call x3iz(vs,n1,n2,n3) call x3bo(vs,1,1,n1,n2,n3) call x3oz(vs,n1,n2,n3) c........construct flux array for x momentum advection in the z direction do i=1,n1 do j=1,n2 do k=1,n3+1 Area=dx1b(i)*dx2a(j) dvg(i,j,k)=dx1b(i)*dx3a(k)*dx2a(j) F(i,j,k)=vs(i,j,k)*.5*(ds3(i,j,k)+ds3(i-1,j,k))*area end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q1(i,j,k)=q1(i,j,k)-(dt/dvg(i,j,k))*(F(i,j,k+1)-F(i,j,k)) end do end do end do return end c......................................................................... c subroutine ymoment_z_advect c c........................................................................ include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c........y momentum advection in the z direction do i=1,n1 do j=1,n2 do k=0,n3 if((Fm3(i,j,k+1)+Fm3(i,j-1,k+1)).gt.0.) then dqmi12=(v2(i,j,k)-v2(i,j,k-1))/dx3b(k) dqpi12=(v2(i,j,k+1)-v2(i,j,k))/dx3b(k+1) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k+1)=v2(i,j,k)+(dx3a(k)-((v3(i,j-1,k+1) + +v3(i,j,k+1))/2.)*dt)*dqi else dqmi12=(v2(i,j,k+1)-v2(i,j,k))/dx3b(k+1) dqpi12=(v2(i,j,k+2)-v2(i,j,k+1))/dx3b(k+2) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k+1)=v2(i,j,k+1)-(dx3a(k+1)+ + ((v3(i,j-1,k+1)+v3(i,j,k+1))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,1,1,n1,n2,n3) call x3iz(vs,n1,n2,n3) call x3bo(vs,1,1,n1,n2,n3) call x3oz(vs,n1,n2,n3) c........construct flux array for y momentum advection in the z direction do i=1,n1 do j=1,n2 do k=1,n3+1 Area=dx1a(i)*dx2b(j) dvg(i,j,k)=dx2b(j)*dx3a(k)*dx1a(i) F(i,j,k)=vs(i,j,k)*.5*(ds3(i,j,k)+ds3(i,j-1,k))*Area end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q2(i,j,k)=q2(i,j,k)-(dt/dvg(i,j,k))*(F(i,j,k+1)-F(i,j,k)) end do end do end do return end c....................................................................... c subroutine zmoment_z_advect c c....................................................................... c......z momentum advection in the z direction include 'von_3defs.h' c......transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save do i=1,n1 do j=1,n2 do k=0,n3 if((Fm3(i,j,k)+Fm3(i,j,k+1)).gt.0.) then dqmi12=(v3(i,j,k)-v3(i,j,k-1))/dx3a(k-1) dqpi12=(v3(i,j,k+1)-v3(i,j,k))/dx3a(k) if((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c..... vs(i,j,k)=v3(i,j,k)+(dx3b(k)-((v3(i,j,k)+ + v3(i,j,k+1))/2.)*dt)*dqi else dqmi12=(v3(i,j,k+1)-v3(i,j,k))/dx3a(k) dqpi12=(v3(i,j,k+2)-v3(i,j,k+1))/dx3a(k+1) if ((dqmi12*dqpi12).gt.0.) then dqi=(dqmi12*dqpi12)/(dqmi12+dqpi12) else dqi=0. end if c.........donor cell c dqi=0. c.................. vs(i,j,k)=v3(i,j,k+1)-(dx3b(k+1)+((v3(i,j,k)+ + v3(i,j,k+1))/2.)*dt)*dqi end if end do end do end do call x1bi(vs,0,1,n1,n2,n3) call x1bo(vs,0,1,n1,n2,n3) call x2bi(vs,0,1,n1,n2,n3) call x2bo(vs,0,1,n1,n2,n3) call x3bi(vs,1,-1,n1,n2,n3) call x3iz(vs,n1,n2,n3) call x3bo(vs,1,-1,n1,n2,n3) call x3oz(vs,n1,n2,n3) c.. construct flux array for z momentum in the z direction do i=1,n1 do j=1,n2 do k=0,n3 Area=dx1a(i)*dx2a(j) dvg(i,j,k)=dx2a(j)*dx3b(k)*dx1a(i) if(k.eq.1) F(i,j,0)=F(i,j,1) F(i,j,k)=vs(i,j,k)*.5*(ds3(i,j,k) + +ds3(i,j,k+1))*Area end do end do end do do i=1,n1 do j=1,n2 do k=1,n3 q3(i,j,k)=q3(i,j,k)-(dt/dvg(i,j,k))*(F(i,j,k)-F(i,j,k-1)) end do end do end do return end c*............................................................... c subroutine energy_x_advect c c*............................................................... c*......energy advection in x direction include 'von_3defs.h' c......Transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......set up interpolated energies at x zone interfaces do i=1,n1+1 do j=1,n2 do k=1,n3 if(v1(i,j,k).gt.0.) then dqpi12=(e(i,j,k)-e(i-1,j,k))/dx1b(i) dqmi12=(e(i-1,j,k)-e(i-2,j,k))/dx1b(i-1) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es1(i,j,k)=e(i-1,j,k)+(dx1a(i-1)-v1(i,j,k)*dt)*dqi else dqpi12=(e(i+1,j,k)-e(i,j,k))/dx1b(i+1) dqmi12=(e(i,j,k)-e(i-1,j,k))/dx1b(i) if((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es1(i,j,k)=e(i,j,k)-(dx1a(i)+v1(i,j,k)*dt)*dqi end if end do end do end do c......Face centered fluxes in the x direction do i=1,n1+1 do j=1,n2 do k=1,n3 es1(i,j,k)=es1(i,j,k)*v1(i,j,k) Fe1(i,j,k)=es1(i,j,k)*dx3a(k)*dx2a(j) end do end do end do call x2bi(es1,0,1,n1,n2,n3) call x2bi(fe1,0,1,n1,n2,n3) call x2bo(es1,0,1,n1,n2,n3) call x2bo(fe1,0,1,n1,n2,n3) call x1bi(es1,1,1,n1,n2,n3) call x1bi(fe1,1,-1,n1,n2,n3) call x1iz(fe1,n1,n2,n3) call x1bo(es1,1,1,n1,n2,n3) call x1bo(fe1,1,-1,n1,n2,n3) call x1oz(fe1,n1,n2,n3) call x3bi(es1,0,1,n1,n2,n3) call x3bi(fe1,0,1,n1,n2,n3) call x3bo(es1,0,1,n1,n2,n3) call x3bo(fe1,0,1,n1,n2,n3) c......update the energy do i=1,n1 do j=1,n2 do k=1,n3 e(i,j,k)=e(i,j,k)-(dt/dv(i,j,k))*(Fe1(i+1,j,k)-Fe1(i,j,k)) end do end do end do call x1bi(e,0,1,n1,n2,n3) call x1bo(e,0,1,n1,n2,n3) call x2bi(e,0,1,n1,n2,n3) call x2bo(e,0,1,n1,n2,n3) call x3bi(e,0,1,n1,n2,n3) call x3bo(e,0,1,n1,n2,n3) return end c*............................................................... c subroutine energy_y_advect c c*............................................................... c*......energy advection in y direction include 'von_3defs.h' c......Transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......set up interpolated energies at y zone interfaces do i=1,n1 do j=1,n2+1 do k=1,n3 if(v2(i,j,k).gt.0.) then dqpi12=(e(i,j,k)-e(i,j-1,k))/dx2b(j) dqmi12=(e(i,j-1,k)-e(i,j-2,k))/dx2b(j-1) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es2(i,j,k)=e(i,j-1,k)+(dx2a(j-1)-v2(i,j,k)*dt)*dqi else dqpi12=(e(i,j+1,k)-e(i,j,k))/dx2b(j+1) dqmi12=(e(i,j,k)-e(i,j-1,k))/dx2b(j) if((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es2(i,j,k)=e(i,j,k)-(dx2a(j)+v2(i,j,k)*dt)*dqi end if end do end do end do c......Face centered fluxes in the y direction do i=1,n1 do j=1,n2+1 do k=1,n3 es2(i,j,k)=es2(i,j,k)*v2(i,j,k) Fe2(i,j,k)=es2(i,j,k)*dx3a(k)*dx1a(i) end do end do end do call x1bi(es2,0,1,n1,n2,n3) call x1bi(fe2,0,1,n1,n2,n3) call x1bo(es2,0,1,n1,n2,n3) call x1bo(fe2,0,1,n1,n2,n3) call x2bi(es2,1,1,n1,n2,n3) call x2bi(fe2,1,-1,n1,n2,n3) call x2iz(fe2,n1,n2,n3) call x2bo(es2,1,1,n1,n2,n3) call x2bo(fe2,1,-1,n1,n2,n3) call x2oz(fe2,n1,n2,n3) call x3bi(es2,0,1,n1,n2,n3) call x3bi(fe2,0,1,n1,n2,n3) call x3bo(es2,0,1,n1,n2,n3) call x3bo(fe2,0,1,n1,n2,n3) c......update the energy do i=1,n1 do j=1,n2 do k=1,n3 e(i,j,k)=e(i,j,k)-(dt/dv(i,j,k))*(Fe2(i,j+1,k)-Fe2(i,j,k)) end do end do end do call x1bi(e,0,1,n1,n2,n3) call x1bo(e,0,1,n1,n2,n3) call x2bi(e,0,1,n1,n2,n3) call x2bo(e,0,1,n1,n2,n3) call x3bi(e,0,1,n1,n2,n3) call x3bo(e,0,1,n1,n2,n3) return end c*............................................................... c subroutine energy_z_advect c c*............................................................... c*......energy advection in z direction include 'von_3defs.h' c......Transport related common variable: common /trans_array/ F(-1:n1+2,-1:n2+2,-1:n3+2), + dvg(-1:n1+2,-1:n2+2,-1:n3+2), + q1(-1:n1+2,-1:n2+2,-1:n3+2), + q2(-1:n1+2,-1:n2+2,-1:n3+2), + q3(-1:n1+2,-1:n2+2,-1:n3+2) save c......set up interpolated energies at z zone interfaces do i=1,n1 do j=1,n2 do k=1,n3+1 if(v3(i,j,k).gt.0.) then dqpi12=(e(i,j,k)-e(i,j,k-1))/dx3b(k) dqmi12=(e(i,j,k-1)-e(i,j,k-2))/dx3b(k-1) if ((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es3(i,j,k)=e(i,j,k-1)+(dx3a(k-1)-v3(i,j,k)*dt)*dqi else dqpi12=(e(i,j,k+1)-e(i,j,k))/dx3b(k+1) dqmi12=(e(i,j,k)-e(i,j,k-1))/dx3b(k) if((dqpi12*dqmi12).gt.0.) then dqi=(dqpi12*dqmi12)/(dqpi12+dqmi12) else dqi=0. end if es3(i,j,k)=e(i,j,k)-(dx3a(k)+v3(i,j,k)*dt)*dqi end if end do end do end do c......Face centered fluxes in the z direction do i=1,n1 do j=1,n2 do k=1,n3+1 es3(i,j,k)=es3(i,j,k)*v3(i,j,k) Fe3(i,j,k)=es3(i,j,k)*dx1a(i)*dx2a(j) end do end do end do call x1bi(es3,0,1,n1,n2,n3) call x1bi(fe3,0,1,n1,n2,n3) call x1bo(es3,0,1,n1,n2,n3) call x1bo(fe3,0,1,n1,n2,n3) call x3bi(es3,1,1,n1,n2,n3) call x3bi(fe3,1,-1,n1,n2,n3) call x3iz(fe3,n1,n2,n3) call x3bo(es3,1,1,n1,n2,n3) call x3bo(fe3,1,-1,n1,n2,n3) call x3oz(fe3,n1,n2,n3) call x2bi(es3,0,1,n1,n2,n3) call x2bi(fe3,0,1,n1,n2,n3) call x2bo(es3,0,1,n1,n2,n3) call x2bo(fe3,0,1,n1,n2,n3) c......update the energy do i=1,n1 do j=1,n2 do k=1,n3 e(i,j,k)=e(i,j,k)-(dt/dv(i,j,k))*(Fe3(i,j,k+1)-Fe3(i,j,k)) end do end do end do call x1bi(e,0,1,n1,n2,n3) call x1bo(e,0,1,n1,n2,n3) call x2bi(e,0,1,n1,n2,n3) call x2bo(e,0,1,n1,n2,n3) call x3bi(e,0,1,n1,n2,n3) call x3bo(e,0,1,n1,n2,n3) return end c........................................................................ c subroutine x1bi(var,i_add,i_sign,n1,n2,n3) c c........................................................................ c*.......this line was missing in the original version c*.......in each boundary routines implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' ic=-1 do i=-1,0 ic=ic+1 do j=-1,n2+2 do k=-1,n3+2 var(i,j,k)=i_sign*var(2-ic+i_add,j,k) end do end do end do return end c........................................................................ c subroutine x2bi(var,i_add,i_sign,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' jc=-1 do j=-1,0 jc=jc+1 do i=-1,n1+2 do k=-1,n3+2 var(i,j,k)=i_sign*var(i,2-jc+i_add,k) end do end do end do return end c........................................................................ c subroutine x3bi(var,i_add,i_sign,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' kc=-1 do k=-1,0 kc=kc+1 do i=-1,n1+2 do j=-1,n2+2 var(i,j,k)=i_sign*var(i,j,2-kc+i_add) end do end do end do return end c........................................................................ c subroutine x1iz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do j=-1,n2+2 do k=-1,n3+2 var(1,j,k)=0. end do end do return end c........................................................................ c subroutine x2iz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do i=-1,n1+2 do k=-1,n3+2 var(i,1,k)=0. end do end do return end c........................................................................ c subroutine x3iz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do i=-1,n1+2 do j=-1,n2+2 var(i,j,1)=0. end do end do return end c........................................................................ c subroutine x1bo(var,i_add,i_sign,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' ic=-1 do i=n1+1,n1+2 ic=ic+1 do j=-1,n2+2 do k=-1,n3+2 var(i,j,k)=i_sign*var(n1-ic+i_add,j,k) end do end do end do return end c........................................................................ c subroutine x2bo(var,i_add,i_sign,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' jc=-1 do j=n2+1,n2+2 jc=jc+1 do i=-1,n1+2 do k=-1,n3+2 var(i,j,k)=i_sign*var(i,n2-jc+i_add,k) end do end do end do return end c........................................................................ c subroutine x3bo(var,i_add,i_sign,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine implements reflecting inner boundary condition c........on standard array variable `var' kc=-1 do k=n3+1,n3+2 kc=kc+1 do i=-1,n1+2 do j=-1,n2+2 var(i,j,k)=i_sign*var(i,j,n3-kc+i_add) end do end do end do return end c........................................................................ c subroutine x1oz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do j=1,n2 do k=1,n3 var(n1+1,j,k)=0. end do end do return end c........................................................................ c subroutine x2oz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do i=1,n1 do k=1,n3 var(i,n2+1,k)=0. end do end do return end c........................................................................ c subroutine x3oz(var,n1,n2,n3) c c........................................................................ implicit real*8(a-h,o-z), integer*4 (i-n) dimension var(-1:n1+2,-1:n2+2,-1:n3+2) c........routine zeros face-centered variables along inner boundary do i=1,n1 do j=1,n2 var(i,j,n3+1)=0. end do end do return end c..................End Routine........................................