aboutsummaryrefslogtreecommitdiff
path: root/src/Extract.F
diff options
context:
space:
mode:
Diffstat (limited to 'src/Extract.F')
-rw-r--r--src/Extract.F814
1 files changed, 814 insertions, 0 deletions
diff --git a/src/Extract.F b/src/Extract.F
new file mode 100644
index 0000000..1aa8534
--- /dev/null
+++ b/src/Extract.F
@@ -0,0 +1,814 @@
+c/*@@
+c @file Extract.F
+c @date 28th October 1997
+c @author Gab Allen
+c @desc Waveform extraction
+c
+c @enddesc
+c@@*/
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c/*@@
+c @routine Extract
+c @date 28th October 1997
+c @author Gab Allen
+c @desc Entry point for waveform extraction routines
+c
+c @enddesc
+c @calls D3_extract
+c @calledby No idea
+c @history
+c
+c @endhistory
+c@@*/
+
+
+ SUBROUTINE Extract(CCTK_FARGUMENTS)
+
+ USE D3_extract_int
+
+ IMPLICIT NONE
+
+ DECLARE_CCTK_FARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
+
+c Non-Cactus input variables for D3_extract
+
+ INTEGER ::
+ & igrid,lmode,mmode,Psi_power=1
+ INTEGER ::
+ & do_ADMmass(2)
+ CCTK_REAL ::
+ & orig(3),radius
+ CCTK_REAL ::
+ & x_1d(cctk_lsh(1)),y_1d(cctk_lsh(2)),z_1d(cctk_lsh(3))
+
+
+c Output variables from D3_extract
+
+ CCTK_REAL ::
+ & dtaudt,ADMmass(2),mass,rsch,momentum(3),spin(3)
+ CCTK_REAL,ALLOCATABLE ::
+ & Qodd(:,:,:), Qeven(:,:,:)
+
+
+c Local variables
+
+ INTEGER ::
+ & ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep,
+ & il,im,it,ndet,idet,nx_parfile,nchar,iconv,ioutput
+ INTEGER,SAVE :: openfiles = 1
+ INTEGER,PARAMETER ::
+ & max_detectors = 9,number_timecoords = 2
+ CCTK_REAL,PARAMETER ::
+ & two = 2.0D0
+ CCTK_REAL ::
+ & time,r1,r2,dr,r_max,xmin,xmax,ymin,ymax,zmin,zmax,gf_min,gf_max
+ CCTK_REAL ::
+ & detector(max_detectors)
+ CHARACTER*100 ::
+ & out1,out2,massfile,rschfile,ADMmassfile1,ADMmassfile2,
+ & momentumfile1,momentumfile2,momentumfile3,spinfile1,spinfile2,spinfile3
+ CHARACTER*3 ::
+ & timestring
+
+ LOGICAL :: use_proper,use_coordinate,do_output
+ LOGICAL,SAVE :: firsttime
+ CCTK_REAL,DIMENSION(10),SAVE :: proper_time
+ CCTK_REAL,SAVE :: last_time
+ INTEGER, save :: open_file_level(10) = 0
+
+ INTEGER myproc,ierr
+ CCTK_REAL dx,dy,dz
+ character*200 filestr
+ character*80 infoline
+
+c ------------------------------------------------------------------
+
+c Get the output directory sorted out
+ call CCTK_FortranString(nchar,outdir,filestr)
+ if (firsttime) then
+ call CCTK_mkdir(ierr,filestr);
+ firsttime = .FALSE.
+ end if
+
+c ------------------------------------------------------------------
+
+ myproc = CCTK_MyProc(cctkGH)
+ dx = cctk_delta_space(1)
+ dy = cctk_delta_space(2)
+ dz = cctk_delta_space(3)
+
+c ------------------------------------------------------------------
+c
+c 0. Check to see if Extract should have been called
+c
+c ------------------------------------------------------------------
+
+ IF (MODULO(int(cctk_iteration),int(itout)) .NE. 0) RETURN
+
+ if (verbose == 1) then
+ write(infoline,'(A24,G12.7)')
+ & 'Calling Extract at time ',cctk_time
+ call CCTK_INFO(infoline)
+ end if
+
+c ------------------------------------------------------------------
+c
+c 1. Initial Stuff
+c
+c ------------------------------------------------------------------
+
+
+c Get the number of polar divisions
+
+ IF (MOD(int(Nt),2) == 1 .OR. MOD(int(Np),2) == 1 .OR.
+ & Nt < 0 .OR. Np < 0) THEN
+ call CCTK_WARN(0,"Error in Nt or Np in Extract")
+ END IF
+
+
+c Get the origin of spherical symmetry
+
+ orig(1) = origin_x
+ orig(2) = origin_y
+ orig(3) = origin_z
+
+
+c Set the value of igrid
+
+ IF (CCTK_EQUALS(domain,"octant")) THEN
+ igrid = 1
+ ELSEIF (CCTK_EQUALS(domain,"bitant")) THEN
+ igrid = 3
+ ELSE
+ igrid = 0
+ END IF
+
+#ifdef THORN_CARTOON_2D
+ IF (contains("cartoon_active","yes") .NE. 0) THEN
+ igrid = 2
+ Np = 1
+ END IF
+#endif
+
+c Create 1D coordinate arrays
+
+ do ix = 1, cctk_lsh(1)
+ x_1d(ix) = x(ix,1,1)
+ enddo
+ do iy = 1, cctk_lsh(2)
+ y_1d(iy) = y(1,iy,1)
+ enddo
+ do iz = 1, cctk_lsh(3)
+ z_1d(iz) = z(1,1,iz)
+ enddo
+
+c See if the ADM mass should be calculated and output
+ IF (doADMmass == 1) THEN
+ do_ADMmass(1) = 1
+ IF (use_conformal == 1) THEN
+ do_ADMmass(2) = 1
+ ELSE
+ do_ADMmass(2) = 0
+ ENDIF
+ ELSE
+ do_ADMmass(:) = 0
+ END IF
+
+c Set array containing gtt ... for the moment I am lazy and
+c just make it the lapse and do not include the shift!
+ g00 = alp
+
+c What kind of timecoordinate do I want to use
+ if (CCTK_EQUALS(timecoord,"both")) then
+ use_proper = .true.
+ use_coordinate = .true.
+ if (cctk_iteration == 0) then
+ proper_time = 0
+ last_time = cctk_time
+ end if
+ elseif (CCTK_EQUALS(timecoord,"proper")) then
+ use_proper = .true.
+ use_coordinate = .false.
+ if (cctk_iteration == 0) then
+ proper_time = 0
+ last_time = cctk_time
+ end if
+ elseif (CCTK_EQUALS(timecoord,"coord")) then
+ use_proper = .false.
+ use_coordinate = .true.
+ else
+ use_proper = .false.
+ use_coordinate = .true.
+ endif
+
+c ------------------------------------------------------------------
+c
+c 2. Sort out which modes to extract
+c
+c ------------------------------------------------------------------
+
+ lmode = l_mode
+
+ IF (all_modes == 1) THEN
+ lmin = 2
+ lmax = lmode
+c Get m_mode in this case too, or the T3E has a junk value (PW)
+ mmode = 0
+ ALLOCATE(Qodd(2,2:lmode,0:lmode),Qeven(2,2:lmode,0:lmode))
+ ELSE
+ lmin = lmode
+ lmax = lmode
+ mmode = m_mode
+ ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode))
+ ENDIF
+
+
+c Check do not have bad values for the modes
+
+ IF (lmode < 2 .OR. mmode > lmode .OR. mmode < 0) THEN
+ call CCTK_WARN(0,"Error in lmode,mmode in Extract")
+ END IF
+ IF (lmode > 9) THEN
+ call CCTK_WARN(4,"Filenames will probably be crazy in Extract")
+ END IF
+
+c If using an octant, do not do odd modes
+
+ IF (igrid == 1 .OR. igrid == 2) THEN
+ lstep = 2 ; mstep = 2
+ ELSE
+ lstep = 1 ; mstep = 1
+ END IF
+
+c If we have cartoon, then jump past all the m-modes in loop
+ IF (igrid == 2) mstep = lmode+1
+
+
+c ------------------------------------------------------------------
+c
+c 3. Find maximum radius of extraction on grid
+c
+c ------------------------------------------------------------------
+
+ call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,"x")
+ call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,"y")
+ call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z")
+
+ IF (igrid == 2) THEN
+ r_max = MIN(xmax-two*Dx-orig(1),
+ & zmax-two*Dz-orig(3))
+ ELSE IF (igrid == 1) THEN
+ r_max = MIN(zmax-two*Dx-orig(1),
+ & ymax-two*Dy-orig(2),
+ & zmax-two*Dz-orig(3))
+ ELSE IF (igrid == 3) THEN
+ r_max = MIN(xmax-two*Dx-orig(1),
+ & ymax-two*Dy-orig(2),
+ & zmax-two*Dz-orig(3),
+ & DABS(xmin+two*Dx-orig(1)),
+ & DABS(ymin+two*Dy-orig(2)))
+ ELSE
+ r_max = MIN(xmax-two*Dx-orig(1),
+ & ymax-two*Dy-orig(2),
+ & zmax-two*Dz-orig(3),
+ & DABS(xmin+two*Dx-orig(1)),
+ & DABS(ymin+two*Dy-orig(2)),
+ & DABS(zmin+two*Dz-orig(3)))
+ END IF
+
+
+
+
+c ------------------------------------------------------------------
+c
+c 4. Extract Cauchy initial data for a linear wave equation
+c
+c ------------------------------------------------------------------
+
+ test_Cauchy: IF (Cauchy == 1) THEN
+
+ it = Cauchy_timestep
+
+
+c check_timestep: IF (cctk_iteration == it) THEN
+
+ check_timestep: IF (open_file_level(1) .eq. 0) THEN
+
+c Open output files
+ write (*,*) "Open output file on level ", 1
+ open_file_level(1)=1;
+
+ test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN
+
+ out1 = filestr(1:nchar)//"/rsch_ini.rl"
+ OPEN(UNIT= 77,file = out1,STATUS="unknown")
+ out1 = filestr(1:nchar)//"/mass_ini.rl"
+ OPEN(UNIT= 78,file = out1,STATUS="unknown")
+ IF (do_ADMmass(1) == 1) THEN
+ out1 = filestr(1:nchar)//"/ADMmass_ini.rl"
+ OPEN(UNIT= 79,file = out1,STATUS="unknown")
+ ENDIF
+ IF (do_ADMmass(2) == 1) THEN
+ out1 = filestr(1:nchar)//"/ADMmassc_ini.rl"
+ OPEN(UNIT= 80,file = out1,STATUS="unknown")
+ ENDIF
+
+ IF (do_momentum == 1) THEN
+ out1 = filestr(1:nchar)//"/momentum_x_ini.rl"
+ OPEN(UNIT=781,file = out1,STATUS="unknown")
+ out1 = filestr(1:nchar)//"/momentum_y_ini.rl"
+ OPEN(UNIT=782,file = out1,STATUS="unknown")
+ out1 = filestr(1:nchar)//"/momentum_z_ini.rl"
+ OPEN(UNIT=783,file = out1,STATUS="unknown")
+ ENDIF
+ IF (do_spin == 1) THEN
+ out1 = filestr(1:nchar)//"/spin_x_ini.rl"
+ OPEN(UNIT=784,file = out1,STATUS="unknown")
+ out1 = filestr(1:nchar)//"/spin_y_ini.rl"
+ OPEN(UNIT=785,file = out1,STATUS="unknown")
+ out1 = filestr(1:nchar)//"/spin_z_ini.rl"
+ OPEN(UNIT=786,file = out1,STATUS="unknown")
+ ENDIF
+
+ loop_l1: DO il = lmin,lmax,lstep
+
+ IF (all_modes == 0) THEN
+ mmin = mmode ; mmax = mmode
+ ELSE
+ mmin = 0 ; mmax = il
+ END IF
+
+ loop_m1: DO im = mmin,mmax,lstep
+
+
+ fn1 = il*10+im
+ out1 = filestr(1:nchar)//"/Qeven_ini_"//
+ & CHAR(il+48)//CHAR(im+48)//".rl"
+ OPEN(UNIT=fn1,FILE=out1,STATUS="unknown")
+
+c Only print odd-parity if full grid
+ IF (igrid == 0) THEN
+ fn2 = 100+il*10+im
+ out2 = filestr(1:nchar)//"/Qodd_ini_"//
+ & CHAR(il+48)//CHAR(im+48)//".rl"
+ OPEN(UNIT=fn2,FILE=out2,STATUS="unknown")
+ END IF
+
+ END DO loop_m1
+
+ END DO loop_l1
+
+ END IF test_myproc1
+
+
+c Find range of extraction radii
+
+ r1 = Cauchy_r1
+ r2 = r_max
+ dr = Cauchy_dr
+
+
+
+c Do extraction at each radius
+
+ IF (verbose == 1) THEN
+ WRITE(*,*)
+ WRITE(*,*) "Extracting Cauchy initial data"
+ WRITE(*,*) "------------------------------"
+ WRITE(*,*) "From r =",r1
+ WRITE(*,*) "To r =",r2
+ WRITE(*,*) "In steps of ",dr
+ WRITE(*,*)
+ END IF
+
+ radius = r1
+
+ extract_at_each_radius: DO WHILE (radius < r2)
+
+ CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,
+ & igrid,orig,myproc,Nt,Np,all_modes,lmode,
+ & mmode,x_1d,y_1d,z_1d,Dx,Dy,Dz,Psi_power,Psi,g00,
+ & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz,
+ & radius,ADMmass,momentum,spin,mass,rsch,
+ & Qodd,Qeven,temp3d,dtaudt)
+
+ IF (verbose == 1) THEN
+ WRITE(*,*) "Extracted at r =",radius
+ WRITE(*,*) "Schwarzschild radius =",rsch
+ WRITE(*,*) "Schwarzschild mass =",mass
+ if (do_ADMmass(1) == 1) then
+ WRITE(*,*) "ADM mass =",ADMmass(1)
+ end if
+ if (do_ADMmass(2) == 1) then
+ WRITE(*,*) "ADM mass (conformal) =",ADMmass(2)
+ end if
+ if (do_momentum == 1) then
+ WRITE(*,*) "momentum = (",momentum(1),",",
+ & momentum(2),",",momentum(3),")"
+ end if
+ if (do_spin == 1) then
+ WRITE(*,*) "spin = (",spin(1),",",
+ & spin(2),",",spin(3),")"
+ end if
+ ENDIF
+
+
+c Write to file at each radius
+
+ test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN
+
+ WRITE( 77,*) radius,rsch
+ WRITE( 78,*) radius,mass
+ IF (do_ADMmass(1) == 1) WRITE( 79,*) radius,ADMmass(1)
+ IF (do_ADMmass(2) == 1) WRITE( 80,*) radius,ADMmass(2)
+ IF (do_momentum == 1) then
+ WRITE(781,*) radius,momentum(1)
+ WRITE(782,*) radius,momentum(2)
+ WRITE(783,*) radius,momentum(3)
+ end if
+
+ IF (do_spin == 1) THEN
+ WRITE(784,*) radius,spin(1)
+ WRITE(785,*) radius,spin(2)
+ WRITE(786,*) radius,spin(3)
+ END IF
+
+ loop_l2: DO il = lmin,lmax,lstep
+
+ IF (all_modes == 0) THEN
+ mmin = mmode ; mmax = mmode
+ ELSE
+ mmin = 0 ; mmax = il
+ END IF
+
+ loop_m2: DO im = mmin,mmax,mstep
+
+ fn1 = il*10+im
+ WRITE(fn1,*) rsch,Qeven(:,il,im)
+
+c Only print odd-parity if full grid
+ IF (igrid == 0) THEN
+ fn2 = 100+il*10+im
+ WRITE(fn2,*) rsch,Qeven(:,il,im)
+ END IF
+
+ END DO loop_m2
+
+ END DO loop_l2
+
+ END IF test_myproc2
+
+ radius = radius+dr
+
+ END DO extract_at_each_radius
+
+
+c Now close all the files
+
+ test_myproc3: IF (myproc == 0) THEN
+
+ CLOSE( 77) ! Schwarzschild radius
+ CLOSE( 78) ! Schwarzschild mass
+ IF (do_ADMmass(1) == 1) CLOSE( 79) ! ADM mass
+ IF (do_ADMmass(2) == 1) CLOSE( 80) ! ADM mass
+
+ IF (do_momentum == 1) THEN
+ CLOSE(781) ! momentum
+ CLOSE(782) ! momentum
+ CLOSE(783) ! momentum
+ end if
+
+ IF (do_spin == 1) THEN
+ CLOSE(784) ! spin
+ CLOSE(785) ! spin
+ CLOSE(786) ! spin
+ END IF
+
+ loop_l3: DO il = lmin,lmax,lstep
+
+ IF (all_modes == 0) THEN
+ mmin = mmode ; mmax = mmode
+ ELSE
+ mmin = 0 ; mmax = il
+ END IF
+
+ loop_m3: DO im = mmin,mmax,mstep
+
+ fn1 = il*10+im
+ CLOSE(fn1) ! Qeven_ini_lm.dat
+
+ IF (igrid == 0) THEN
+ fn2 = 100+il*10+im
+ CLOSE(fn2) ! Qodd_ini_lm.dat
+ END IF
+
+ END DO loop_m3
+
+ END DO loop_l3
+
+ END IF test_myproc3
+
+ END IF check_timestep
+
+ END IF test_Cauchy
+
+
+c ------------------------------------------------------------------
+c
+c 5. Extract waveforms at detectors
+c
+c ------------------------------------------------------------------
+
+c Cannot use the conformal equation for ADM mass now
+ do_ADMmass(2) = 0
+
+ ndet = num_detectors
+
+ test_detectors: IF (ndet > 0) THEN
+
+ IF (ndet > 9) THEN
+ call CCTK_WARN(0,"Too many detectors in Extract")
+ END IF
+
+ detector(1) = detector1
+ detector(2) = detector2
+ detector(3) = detector3
+ detector(4) = detector4
+ detector(5) = detector5
+ detector(6) = detector6
+ detector(7) = detector7
+ detector(8) = detector8
+ detector(9) = detector9
+
+ DO idet = ndet+1, max_detectors
+ detector(idet) = 0.0D0
+ END DO
+
+ DO idet = 1, ndet
+ IF (detector(idet) > r_max) THEN
+ IF (openfiles == 1 .OR. cctk_iteration == it) THEN
+ call CCTK_WARN(8,"Removing detectors outside coordinate range")
+c IF (iconv == 0) STOP
+ END IF
+ ndet =idet-1
+ GOTO 404
+ ENDIF
+ END DO
+ 404 CONTINUE
+
+ IF (verbose == 1) THEN
+ IF (openfiles == 1) THEN
+ WRITE(*,*)
+ WRITE(*,*) "Extracting at detectors"
+ WRITE(*,*) "-----------------------"
+ WRITE(*,*) "Using ",ndet," detectors at"
+ DO idet=1,ndet
+ WRITE(*,*) "r = ",detector(idet)
+ ENDDO
+ END IF
+ END IF
+
+ detector_loop: DO idet = 1,ndet
+
+ radius = detector(idet)
+
+ IF (verbose == 1) THEN
+ WRITE(*,*) "Calling extract at radius ... ",radius
+ END IF
+
+ CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,
+ & igrid,orig,myproc,Nt,
+ & Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d,
+ & Dx,Dy,Dz,Psi_power,Psi,g00,
+ & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz,
+ & radius,ADMmass,momentum,spin,mass,rsch,
+ & Qodd,Qeven,temp3d,dtaudt)
+
+ test_verbose:
+ & IF (verbose == 1) THEN
+
+ WRITE(*,*)
+ WRITE(*,*) "Detector ",idet," ..."
+ WRITE(*,*) "Schwarzschild radius =",rsch
+ WRITE(*,*) "Schwarzschild mass =",mass
+
+ loop_l4: DO il = lmin,lmax,lstep
+
+ IF (all_modes == 0) THEN
+ mmin = mmode ; mmax = mmode
+ ELSE
+ mmin = 0 ; mmax = il
+ END IF
+
+ loop_m4: DO im = mmin,mmax,mstep
+
+ WRITE(*,*) filestr(1:nchar)//"/Qeven_"//
+ & CHAR(il+48)//CHAR(im+48),Qeven(:,il,im)
+ WRITE(*,*) filestr(1:nchar)//"/Qodd_"//
+ & CHAR(il+48)//CHAR(im+48),Qodd(:,il,im)
+
+ END DO loop_m4
+
+ END DO loop_l4
+
+ END IF test_verbose
+
+c Output to files
+
+ test_myproc4: IF (myproc == 0) THEN
+
+ DO ioutput=1,number_timecoords
+
+ do_output = .false.
+ timestring = "ERR"
+ time = 0D0
+ IF (ioutput == 1 .AND. use_coordinate) THEN
+ timestring = ".tl"
+ do_output = .true.
+ time = cctk_time
+ ELSEIF (ioutput == 2 .AND. use_proper) THEN
+ timestring = ".ul"
+ do_output = .true.
+ time = proper_time(idet) + (cctk_time-last_time)*dtaudt
+ proper_time(idet) = time
+ if (idet == ndet) last_time = cctk_time
+ ENDIF
+
+ IF (do_output) THEN
+
+c Output extracted radius and mass
+ rschfile = filestr(1:nchar)//"/rsch_"//"R"//
+ & CHAR(idet+48)//timestring
+ massfile = filestr(1:nchar)//"/mass_"//"R"//
+ & CHAR(idet+48)//timestring
+ ADMmassfile1 = filestr(1:nchar)//"/ADMmass_"//"R"//
+ & CHAR(idet+48)//timestring
+ ADMmassfile2 = filestr(1:nchar)//"/ADMmassc_"//"R"//
+ & CHAR(idet+48)//timestring
+ momentumfile1 = filestr(1:nchar)//"/momentum_x_"//"R"//
+ & CHAR(idet+48)//timestring
+ momentumfile2 = filestr(1:nchar)//"/momentum_y_"//"R"//
+ & CHAR(idet+48)//timestring
+ momentumfile3 = filestr(1:nchar)//"/momentum_z_"//"R"//
+ & CHAR(idet+48)//timestring
+ spinfile1 = filestr(1:nchar)//"/spin_x_"//"R"//
+ & CHAR(idet+48)//timestring
+ spinfile2 = filestr(1:nchar)//"/spin_y_"//"R"//
+ & CHAR(idet+48)//timestring
+ spinfile3 = filestr(1:nchar)//"/spin_z_"//"R"//
+ & CHAR(idet+48)//timestring
+
+ IF (openfiles == 1) THEN
+ OPEN(UNIT= 77,FILE=rschfile,STATUS="unknown")
+ OPEN(UNIT= 78,FILE=massfile,STATUS="unknown")
+ IF (do_ADMmass(1) == 1)
+ & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="unknown")
+ IF (do_ADMmass(2) == 1)
+ & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="unknown")
+ IF (do_momentum == 1) THEN
+ OPEN(UNIT=781,FILE=momentumfile1,STATUS="unknown")
+ OPEN(UNIT=782,FILE=momentumfile2,STATUS="unknown")
+ OPEN(UNIT=783,FILE=momentumfile3,STATUS="unknown")
+ END IF
+
+ IF (do_spin == 1) THEN
+ OPEN(UNIT=784,FILE=spinfile1,STATUS="unknown")
+ OPEN(UNIT=785,FILE=spinfile2,STATUS="unknown")
+ OPEN(UNIT=786,FILE=spinfile3,STATUS="unknown")
+ END IF
+ ELSE
+ OPEN(UNIT= 77,FILE=rschfile,STATUS="old",
+ & POSITION="append")
+ OPEN(UNIT= 78,FILE=massfile,STATUS="old",
+ & POSITION="append")
+ IF (do_ADMmass(1) == 1)
+ & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="old",
+ & POSITION="append")
+ IF (do_ADMmass(2) == 1)
+ & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="old",
+ & POSITION="append")
+ IF (do_momentum == 1) THEN
+ OPEN(UNIT=781,FILE=momentumfile1,STATUS="old",
+ & POSITION="append")
+ OPEN(UNIT=782,FILE=momentumfile2,STATUS="old",
+ & POSITION="append")
+ OPEN(UNIT=783,FILE=momentumfile3,STATUS="old",
+ & POSITION="append")
+ END IF
+
+ IF (do_spin == 1)THEN
+ OPEN(UNIT=784,FILE=spinfile1,STATUS="old",
+ & POSITION="append")
+ OPEN(UNIT=785,FILE=spinfile2,STATUS="old",
+ & POSITION="append")
+ OPEN(UNIT=786,FILE=spinfile3,STATUS="old",
+ & POSITION="append")
+ END IF
+
+ ENDIF
+ WRITE( 77,21) time,rsch
+ WRITE( 78,21) time,mass
+ IF (do_ADMmass(1) == 1) WRITE( 79,21) time,ADMmass(1)
+ IF (do_ADMmass(2) == 1) WRITE( 80,21) time,ADMmass(2)
+ IF (do_momentum == 1) THEN
+ WRITE(781,*) cctk_time,momentum(1)
+ WRITE(782,*) cctk_time,momentum(2)
+ WRITE(783,*) cctk_time,momentum(3)
+ END IF
+ IF (do_spin == 1) THEN
+ WRITE(784,*) cctk_time,spin(1)
+ WRITE(785,*) cctk_time,spin(2)
+ WRITE(786,*) cctk_time,spin(3)
+ END IF
+
+ CLOSE( 77)
+ CLOSE( 78)
+ IF (do_ADMmass(1) == 1) CLOSE( 79)
+ IF (do_ADMmass(2) == 1) CLOSE( 80)
+
+ IF (do_momentum == 1) THEN
+ CLOSE(781)
+ CLOSE(782)
+ CLOSE(783)
+ END IF
+ IF (do_spin == 1) THEN
+ CLOSE(782)
+ CLOSE(783)
+ CLOSE(784)
+ END IF
+
+
+c Output gauge invariant variables
+
+ loop_l5: DO il = lmin,lmax,lstep
+
+ IF (all_modes == 0) THEN
+ mmin = mmode ; mmax = mmode
+ ELSE
+ mmin = 0 ; mmax = il
+ END IF
+
+ loop_m5: DO im = mmin,mmax,mstep
+
+ fn1 = 11
+ fn2 = 12
+
+ out1 = filestr(1:nchar)//"/Qeven_"//"R"
+ & //CHAR(idet+48)//"_"//CHAR(il+48)//
+ & CHAR(im+48)//timestring
+ out2 = filestr(1:nchar)//"/Qodd_"//"R"
+ & //CHAR(idet+48)//"_"//CHAR(il+48)//
+ & CHAR(im+48)//timestring
+
+c Write even parity waveforms
+ IF (openfiles == 1) THEN
+ OPEN(UNIT=fn1,FILE=out1,STATUS="unknown")
+ ELSE
+ OPEN(UNIT=fn1,FILE=out1,STATUS="old",
+ & POSITION="append")
+ ENDIF
+
+ WRITE(fn1,20) time,Qeven(:,il,im)
+
+ CLOSE(fn1)
+
+c Only write odd parity waveforms if full grid
+ IF (igrid == 0) THEN
+ IF (openfiles == 1) THEN
+ OPEN(UNIT=fn2,FILE=out2,STATUS="unknown")
+ ELSE
+ OPEN(UNIT=fn2,FILE=out2,STATUS="old",
+ & POSITION="append")
+ ENDIF
+
+ WRITE(fn2,20) time,Qodd(:,il,im)
+
+ CLOSE(fn2)
+ END IF
+
+ END DO loop_m5
+
+ END DO loop_l5
+
+ ENDIF
+ END DO
+ END IF test_myproc4
+
+ END DO detector_loop
+
+ END IF test_detectors
+
+ DEALLOCATE(Qodd,Qeven)
+
+ openfiles = 0
+
+ 20 format(e24.15,e24.15,e24.15)
+ 21 format(e24.15,e24.15)
+
+ END SUBROUTINE Extract
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 17:14:11 +0000