************************************************************************ * This code evolves binary black hole spacetimes using the ADM 3+1 * formalism. A initial spacelike hypersurface is formed by using * an exact form for a known spacetime, or by solving the time-independent * Einstein equations, also known as the constraint equations. This * initial hypersurface is evolved using the time-dependent Einstein * equations. In this manner, a foliation of the full spacetime can * be computed. * The hypersurfaces are freely evolved, i.e.--the constraint * equations are not used in the solution of hypersurfaces subsequent * to the initial one. The evolution is a three-level update scheme, * using time-centered and space-centered derivatives. At the boundaries * of the computational domain, the outer edges of the grid and the * boundary of the excised singularity, we use one-sided differences * to compute spatial derivatives. * Two important subroutines used are the apparent * horizon solver and a generic interpolator/extrapolator. The program * is designed in a modular fasion to allow a variety of finite differencing * schemes and conforms to a DAGH compatible data structure. * * Randall R. Correll * Mijan F. Huq * Scott A. Klasky * Richard A. Matzner * * COPYRIGHT 1995-1996 * * Modification to evolve linearized Teukolsky waves (MFH 1-96) ************************************************************************ program t2 implicit none c Problem parameters integer nx, ny, nz, nxm, nym, nzm, nl, nlmax, nplot integer i, j, k parameter(nxm=65, nym=65, nzm=65) parameter(nlmax=6*(nxm*nxm+nym*nym+nzm*nzm)) real*8 xmin, xmax, $ ymin, ymax, $ zmin, zmax real*8 tmin, tmax integer hdfflag, asciiflag c Black hole parameters real*8 mass real*8 bha(2) real*8 bhr(2,3), bhp(2,3), bhs(2,3) c Grid specifications real*8 x(nxm), y(nym), z(nzm) real*8 dx,dy,dz c Main loop parameters integer it,nt,itsave real*8 dt, dtdx, time real*8 dtsave, timesave real*8 tol c Dynamic field variables real*8 g11(nxm,nym,nzm), g12(nxm,nym,nzm), g13(nxm,nym,nzm), * g22(nxm,nym,nzm), g23(nxm,nym,nzm), g33(nxm,nym,nzm), * g11m1(nxm,nym,nzm), g12m1(nxm,nym,nzm), g13m1(nxm,nym,nzm), * g22m1(nxm,nym,nzm), g23m1(nxm,nym,nzm), g33m1(nxm,nym,nzm), * gup11(nxm,nym,nzm), gup12(nxm,nym,nzm), gup13(nxm,nym,nzm), * gup22(nxm,nym,nzm), gup23(nxm,nym,nzm), gup33(nxm,nym,nzm) real*8 k11(nxm,nym,nzm), k12(nxm,nym,nzm), k13(nxm,nym,nzm), * k21(nxm,nym,nzm), k22(nxm,nym,nzm), k23(nxm,nym,nzm), * k31(nxm,nym,nzm), k32(nxm,nym,nzm), k33(nxm,nym,nzm), * k11m1(nxm,nym,nzm), k12m1(nxm,nym,nzm), k13m1(nxm,nym,nzm), * k21m1(nxm,nym,nzm), k22m1(nxm,nym,nzm), k23m1(nxm,nym,nzm), * k31m1(nxm,nym,nzm), k32m1(nxm,nym,nzm), k33m1(nxm,nym,nzm) c Gauge variables real*8 beta1(nxm,nym,nzm), beta2(nxm,nym,nzm), * beta3(nxm,nym,nzm), * alpha(nxm,nym,nzm) c Spatial tensors real*8 r11(nxm,nym,nzm), r12(nxm,nym,nzm), r13(nxm,nym,nzm), * r22(nxm,nym,nzm), r23(nxm,nym,nzm), r33(nxm,nym,nzm) real*8 c111(nxm,nym,nzm), c112(nxm,nym,nzm), c113(nxm,nym,nzm), * c122(nxm,nym,nzm), c123(nxm,nym,nzm), c133(nxm,nym,nzm), * c211(nxm,nym,nzm), c212(nxm,nym,nzm), c213(nxm,nym,nzm), * c222(nxm,nym,nzm), c223(nxm,nym,nzm), c233(nxm,nym,nzm), * c311(nxm,nym,nzm), c312(nxm,nym,nzm), c313(nxm,nym,nzm), * c322(nxm,nym,nzm), c323(nxm,nym,nzm), c333(nxm,nym,nzm) c Temporary arrays for calculations within subroutines real*8 t11(nxm,nym,nzm), t12(nxm,nym,nzm), t13(nxm,nym,nzm), $ t21(nxm,nym,nzm), t22(nxm,nym,nzm), t23(nxm,nym,nzm), $ t31(nxm,nym,nzm), t32(nxm,nym,nzm), t33(nxm,nym,nzm), $ temp1(nxm,nym,nzm),temp2(nxm,nym,nzm), $ temp3(nxm,nym,nzm) c Variables to output constraint monitors real*8 hamnorm, mom1norm, mom2norm, mom3norm c Characteristic functions real*8 mask(nxm,nym,nzm),maskold(nxm,nym,nzm),d(3*nlmax) integer m(3*nlmax) integer xi(nlmax),yj(nlmax),zk(nlmax) real*8 fx(nlmax),fy(nlmax),fz(nlmax) c Miscellaneous real*8 pi, coeff integer lb(3), smth c Linearized wave parameters real*8 amp, width amp = 1.0e-5 width=1.0 lb(1) = 0 lb(2) = 0 lb(3) = 0 smth = 8 coeff = 0.5 ************************************************************************ write(6,*) 'Begin TCODE run.' write(6,*) c Initialize variables pi = 4.0e0*atan(1.0e0) time = 0.0e0 it = 0 c To write hdf files make hdfflag = .true. hdfflag = 0 c To write ascii files make asciiflag = .true. asciiflag = 0 call input(nx,ny,nz,nt,nplot, $ xmin,xmax,ymin,ymax,zmin,zmax, $ hdfflag,asciiflag, $ bha,bhr,bhp,bhs,tol) dx = (xmax-xmin)/(nx-1) dy = (ymax-ymin)/(ny-1) dz = (zmax-zmin)/(nz-1) write(6,*) '...Make grid' call makegrid(nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax, $ dx,dy,dz,x,y,z,lb) call makegrid3D(nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax, $ dx,dy,dz,t21,t22,t23,t13,lb) call horizon(nx,ny,nz,g11,g12,g13,g22,g23,g33, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ x,y,z,bha,bhr,mask) call copy(nx,ny,nz,mask,maskold) c write(6,*) '...Initialize linearized wave data' c x,y,z,r are 3d arrays for this routine so pass t21, t22, t23, t13 c call linearwave(nx, ny, nz, amp, width, c * t21, t22, t23, t13, g11, g12, g13, g22, g23, g33, c * k11, k12, k13, k21, k22, k23, k31, k32, k33, c * temp1, temp2, temp3) c call initdata(nx,ny,nz,g11,g12,g13,g22,g23,g33, c & gup11,gup12,gup13,gup22,gup23,gup33, c & k11,k12,k13,k21,k22,k23,k31,k32,k33, c & alpha, beta1,beta2,beta3, c & t11,temp1,temp2,temp3,x,y,z,bha,bhr) call init_ascii(nx,ny,nz,g11,g12,g13,g22,g23,g33, & gup11,gup12,gup13,gup22,gup23,gup33, & k11,k12,k13,k21,k22,k23,k31,k32,k33, & alpha, beta1, beta2, beta3, x, y, z, bha,bhr, * t11,temp1, temp2, temp3) call remask(nx,ny,nz,g11,g12,g13,g22,g23,g33, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ mask,maskold) write(6,*) '...Compute excise mask and boundary points' call flags(nx,ny,nz,nlmax,xi,yj,zk,fx,fy,fz,mask,m,d,nl) write(6,*) '...Compute gauge' call setgauge(nx,ny,nz,mask, $ x,y,z,bha,bhr, $ alpha,beta1,beta2,beta3) write(6,*) '...Compute connections and Ricci' write(6,*) '.....inverse metric' call inverse(g11,g12,g13,g22,g23,g33, * gup11,gup12,gup13,gup22, gup23, gup33,temp1, * mask,nx, ny, nz) write(6,*) '.....connection coefficinets' call connect(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, dx, nx, ny, nz,nl,m,d) write(6,*) '.....Ricci Tensor' call riccidd(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, nx, ny, nz) call timestep(dx,dy,dz,dt) dtsave = dt timesave = time itsave = it dt = -0.5e0*dt write(6,*) '...Take Euler step backward to generate n+1 level' dt = 0.5e0*dt call copy(nx,ny,nz,g11,g11m1) call copy(nx,ny,nz,g12,g12m1) call copy(nx,ny,nz,g13,g13m1) call copy(nx,ny,nz,g22,g22m1) call copy(nx,ny,nz,g23,g23m1) call copy(nx,ny,nz,g33,g33m1) call evolve_g(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1,temp2, temp3, r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) call copy(nx,ny,nz,k11,k11m1) call copy(nx,ny,nz,k12,k12m1) call copy(nx,ny,nz,k13,k13m1) call copy(nx,ny,nz,k21,k21m1) call copy(nx,ny,nz,k22,k22m1) call copy(nx,ny,nz,k23,k23m1) call copy(nx,ny,nz,k31,k31m1) call copy(nx,ny,nz,k32,k32m1) call copy(nx,ny,nz,k33,k33m1) call evolve_k(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1, temp2, temp3, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) call swap(nx,ny,nz,g11,g11m1,temp1) call swap(nx,ny,nz,g12,g12m1,temp1) call swap(nx,ny,nz,g13,g13m1,temp1) call swap(nx,ny,nz,g22,g22m1,temp1) call swap(nx,ny,nz,g23,g23m1,temp1) call swap(nx,ny,nz,g33,g33m1,temp1) call swap(nx,ny,nz,k11,k11m1,temp1) call swap(nx,ny,nz,k12,k12m1,temp1) call swap(nx,ny,nz,k13,k13m1,temp1) call swap(nx,ny,nz,k21,k21m1,temp1) call swap(nx,ny,nz,k22,k22m1,temp1) call swap(nx,ny,nz,k23,k23m1,temp1) call swap(nx,ny,nz,k31,k31m1,temp1) call swap(nx,ny,nz,k32,k32m1,temp1) call swap(nx,ny,nz,k33,k33m1,temp1) dt = 2.0e0*dt c Then continue with standard three-level scheme. write(6,*) '...Complete initial backward timestep' time = time + dt call inverse(g11,g12,g13,g22,g23,g33, * gup11,gup12,gup13,gup22, gup23, gup33,temp1, * mask,nx, ny, nz) call connect(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, dx, nx, ny, nz,nl,m,d) call riccidd(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, nx, ny, nz) call evolve_g(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1,temp2, temp3, r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) call evolve_k(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1, temp2, temp3, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) c Restore times and time levels dt = dtsave time = timesave it = itsave write(6,*) '...Re-initialize t=0 data' call makegrid3D(nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax, $ dx,dy,dz,t21,t22,t23,t13,lb) c x,y,z,r are 3d arrays for this routine so pass t21, t22, t23, t13 c call linearwave(nx, ny, nz, amp, width, c * t21, t22, t23, t13, g11, g12, g13, g22, g23, g33, c * k11, k12, k13, k21, k22, k23, k31, k32, k33, c * temp1, temp2, temp3) c call initdata(nx,ny,nz,g11,g12,g13,g22,g23,g33, c & gup11,gup12,gup13,gup22,gup23,gup33, c & k11,k12,k13,k21,k22,k23,k31,k32,k33, c & alpha, beta1,beta2,beta3, c & t11,temp1,temp2,temp3,x,y,z,bha,bhr) call init_ascii(nx,ny,nz,g11,g12,g13,g22,g23,g33, & gup11,gup12,gup13,gup22,gup23,gup33, & k11,k12,k13,k21,k22,k23,k31,k32,k33, & alpha, beta1, beta2, beta3, x, y, z, bha,bhr, * t11,temp1, temp2, temp3) call remask(nx,ny,nz,g11,g12,g13,g22,g23,g33, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ mask,maskold) c call initdata(nx,ny,nz,g11,g12,g13,g22,g23,g33, c $ k11,k12,k13,k21,k22,k23,k31,k32,k33, c $ x,y,z,bha,bhr,temp1) call inverse(g11,g12,g13,g22,g23,g33, * gup11,gup12,gup13,gup22, gup23, gup33,temp1, * mask,nx, ny, nz) call connect(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, dx, nx, ny, nz,nl,m,d) call riccidd(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, nx, ny, nz) c Begin main loop, time-centered, space-centered. write(6,*) '...Begin main loop, forward in time' c Output constraint monitors write(6,*) 'Iter Time HamErr Mom1Err Mom2Err Mom3Err' call constraint(nx,ny,nz,dx,dy,dz,mask, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ hamnorm,mom1norm,mom2norm,mom3norm, $ c111,c112,c113,c122,c123,c133, $ c211,c212,c213,c222,c223,c233, $ c311,c312,c313,c322,c323,c333, $ r11,r12,r13,r22,r23,r33, $ gup11,gup12,gup13,gup22,gup23,gup33, $ t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1, $ dx,m,d,nl) write(6,'(i5,f9.5,1P,4e12.3)') it,time, $ hamnorm,mom1norm,mom2norm,mom3norm write(2,'(i5,f9.5,1P,4e12.3)') it,time, $ hamnorm,mom1norm,mom2norm,mom3norm do it = 1, nt call setgauge(nx,ny,nz,mask, $ x,y,z,bha,bhr, $ alpha,beta1,beta2,beta3) call timestep(dx,dy,dz,dt) time = time + dt call evolve_g(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1,temp2, temp3, r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) call evolve_k(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11m1,k12m1,k13m1, $ k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1, $ k11, k12, k13, $ k21, k22, k23, $ k31, k32, k33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ beta1, beta2, beta3, alpha, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ temp1, temp2, temp3, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, dt,nx, ny, nz) call swap(nx,ny,nz,g11,g11m1,temp1) call swap(nx,ny,nz,g12,g12m1,temp1) call swap(nx,ny,nz,g13,g13m1,temp1) call swap(nx,ny,nz,g22,g22m1,temp1) call swap(nx,ny,nz,g23,g23m1,temp1) call swap(nx,ny,nz,g33,g33m1,temp1) call swap(nx,ny,nz,k11,k11m1,temp1) call swap(nx,ny,nz,k12,k12m1,temp1) call swap(nx,ny,nz,k13,k13m1,temp1) call swap(nx,ny,nz,k21,k21m1,temp1) call swap(nx,ny,nz,k22,k22m1,temp1) call swap(nx,ny,nz,k23,k23m1,temp1) call swap(nx,ny,nz,k31,k31m1,temp1) call swap(nx,ny,nz,k32,k32m1,temp1) call swap(nx,ny,nz,k33,k33m1,temp1) call copy(nx,ny,nz,mask,maskold) call horizon(nx,ny,nz,g11,g12,g13,g22,g23,g33, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ x,y,z,bha,bhr,mask) call remask(nx,ny,nz,g11,g12,g13,g22,g23,g33, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ mask,maskold) call remask(nx,ny,nz,g11m1,g12m1,g13m1,g22m1,g23m1,g33m1, $ k11m1,k12m1,k13m1,k21m1,k22m1,k23m1, $ k31m1,k32m1,k33m1,mask,maskold) call flags(nx,ny,nz,nlmax,xi,yj,zk,fx,fy,fz,mask,m,d,nl) call inverse(g11,g12,g13,g22,g23,g33, * gup11,gup12,gup13,gup22, gup23, gup33,temp1, * mask,nx, ny, nz) call connect(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, dx, nx, ny, nz,nl,m,d) call riccidd(g11, g12, g13, g22, g23, g33, $ gup11, gup12, gup13, gup22, gup23, gup33, $ c111, c112, c113, c122, c123, c133, $ c211, c212, c213, c222, c223, c233, $ c311, c312, c313, c322, c323, c333, $ t11, t12, t13, t21, t22, t23, t31, t32, t33, $ r11, r12, r13, r22, r23, r33, $ m,d,nl, dx, nx, ny, nz) c output constraint monitors call constraint(nx,ny,nz,dx,dy,dz,mask, $ k11,k12,k13,k21,k22,k23,k31,k32,k33, $ hamnorm,mom1norm,mom2norm,mom3norm, $ c111,c112,c113,c122,c123,c133, $ c211,c212,c213,c222,c223,c233, $ c311,c312,c313,c322,c323,c333, $ r11,r12,r13,r22,r23,r33, $ gup11,gup12,gup13,gup22,gup23,gup33, $ t11,t12,t13,t21,t22,t23,t31,t32,t33,temp1, $ dx,m,d,nl) write(6,'(i5,f9.5,1P,4e12.3)') it,time, $ hamnorm,mom1norm,mom2norm,mom3norm write(2,'(i5,f9.5,1P,4e12.3)') it,time, $ hamnorm,mom1norm,mom2norm,mom3norm call dumpznamed(nx, ny, nz, (nz-1)/2 - 1, $ 'g11',it, x, y, z, g11) c call dumpznamed(nx, ny, nz, (nz-1)/2 - 1, c $ 'g22',it, x, y, z, g22) c call dumpznamed(nx, ny, nz, (nz-1)/2 - 1, c $ 'g33',it, x, y, z, g33) if(hamnorm .gt. tol)goto 1919 end do 1919 continue c Terminate run stop end