Modified output so that data is appended after checkpoint, rather

than deleting any existing files.
* Starting a run from scratch will still wipe out any files that already
  exist in the output directory, as is usual for cactus 1d data.
* If you do a horizon find on the timestep when you are also checkpointing,
  then you'll end up with a duplicate point in the output, as the horizon
  is found again at restart.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@320 89daf98e-ef62-4674-b946-b8ff9de2216c

7 files changed, 243 insertions, 103 deletions

diff --git a/src/AHFinder_dat.F b/src/AHFinder_dat.F
index a7fc589..718d16e 100644
--- a/src/AHFinder_dat.F
+++ b/src/AHFinder_dat.F
@@ -59,7 +59,7 @@
       CCTK_REAL, allocatable, dimension (:,:) :: hflowc,cflowc,nflowc
       CCTK_REAL, allocatable, dimension (:,:) :: hflows,cflows,nflows
 
-      character*200 filestr
+      character(len=200) :: filestr
 
       data ahf_ncall / 0 /
       data status_old   / .false. /
diff --git a/src/AHFinder_dis.F b/src/AHFinder_dis.F
index 894271c..c18b0f8 100644
--- a/src/AHFinder_dis.F
+++ b/src/AHFinder_dis.F
@@ -40,7 +40,8 @@
       CCTK_REAL, allocatable, dimension(:) :: xa,ya,za
       CCTK_REAL, allocatable, dimension(:) :: txx,tyy,tzz,txy,txz,tyz
 
-      character*200 disf
+      character(len=200) :: disf
+      logical disf_exists
 
 
 !     **************************
@@ -571,22 +572,16 @@
 
  100  continue
 
-      disf = filestr(1:nfile) // "/ahf_d12" // ".tl"
-
       if (myproc.eq.0) then
-
-        if (ahf_ncall.eq.1) then
-           open(1,file=disf,form='formatted',status='replace')
-           write(1,"(A18)") '" Horizon distance'
+        inquire(file=disf, exist=disf_exists)
+        if (disf_exists) then
+          open(1,file=disf,form='formatted',status='old', position='append')
         else
-           open(1,file=disf,form='formatted',status='old',
-     .          position='append')
+          open(1,file=disf,form='formatted',status='replace')
+          write(1,"(A18)") '" Horizon distance'
         end if
-
         write(1,"(2ES14.6)") cctk_time,d12
-
         close(1)
-
       end if
 
 
diff --git a/src/AHFinder_flow.F b/src/AHFinder_flow.F
index 810019b..b5e4300 100644
--- a/src/AHFinder_flow.F
+++ b/src/AHFinder_flow.F
@@ -38,7 +38,7 @@
       CCTK_REAL zero,one
       CCTK_REAL dd,aux
 
-      character*200 logf
+      character(len=200) :: logf
 
 !     Description of variables:
 !
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 28ab392..d66c214 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -66,6 +66,7 @@
       CCTK_INT :: use_att, apower, nchars
       logical :: use_rot_att
       logical :: str_comp
+      logical :: gaussf_exists
       external str_comp
 
       character(len=32) :: tmpstr
@@ -1143,10 +1144,9 @@
             gaussf = filestr(1:nfile)//"/ahf.gauss"
          end if
 
-!        On first call replace file.
+         inquire(file=gaussf, exist=gaussf_exists)
 
          if (firstgau) then
-
             if (find3) then
                if (mfind.eq.3) then
                   firstgau = .false.
@@ -1154,49 +1154,20 @@
             else
                firstgau = .false.
             end if
+         end if
 
-            open(1,file=gaussf,form='formatted',
-     .      status='replace')
-
-            write(1,"(A20)") '# GAUSSIAN CURVATURE'
-            write(1,"(A1)")  '#'
-            write(1,"(A35)") '# The data is written in a loop as:'
-            write(1,"(A1)")  '#'
-            write(1,"(A17)") '# do i=1,ntheta'
-            write(1,"(A18)") '#    do j=1,nphi'
-            write(1,"(A27)") '#       write gaussian(i,j)'
-            write(1,"(A11)") '#    end do'
-            write(1,"(A8)")  '# end do'
-            write(1,"(A1)")  '#'
-            write(1,"(A40)") '# theta and phi are subdivided uniformly'
-            write(1,"(A26)") '# according to symmetries:'
-            write(1,"(A1)")  '#'
-            write(1,"(A36)") '# phi=[0,2 pi]   (refx=refy=.false.)'
-            write(1,"(A44)") '# phi=[0,pi]     (refx=.false., refy=.true.)'
-            write(1,"(A35)") '# phi=[0,pi/2]   (refx=refy=.true.)'
-            write(1,"(A1)")  '#'
-            write(1,"(A31)") '# theta=[0,pi]   (refz=.false.)'
-            write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)'
-            write(1,"(A1)")  '#'
-
-            write(1,"(A9,L1)")  '# refx = ',refx
-            write(1,"(A9,L1)")  '# refy = ',refy
-            write(1,"(A9,L1)")  '# refz = ',refz
-            write(1,"(A1)")     '#'
-            write(1,"(A11,I7)") '# ntheta = ',(ntheta+1)
-            write(1,"(A11,I7)") '# nphi   = ',(nphi+1)
+!        On first call replace file.
 
-            write(1,*)
-            write(1,"(A11,I4)")    '# Time step',cctk_iteration
-            write(1,"(A6,ES11.3)") '# Time',cctk_time
-            write(1,"(A6,I4)")     '# Call',ahf_ncall
+         if (.not.gaussf_exists) then
+
+           call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gaussf)
 
 !        On subsequent calls, append to file.
 
          else
 
-            open(1,file=gaussf,form='formatted',
-     .      status='old',position='append')
+            open(1,file=gaussf,form='formatted', status='old',
+     .          position='append')
 
             write(1,*)
             write(1,"(A11,I4)")    '# Time step',cctk_iteration
diff --git a/src/AHFinder_min.F b/src/AHFinder_min.F
index b458d85..1c4a044 100644
--- a/src/AHFinder_min.F
+++ b/src/AHFinder_min.F
@@ -32,7 +32,7 @@
       CCTK_REAL, dimension (NN)    :: P
       CCTK_REAL, dimension (NN,NN) :: XI
 
-      character*200 logf
+      character(len=200) :: logf
 
 !     Description of variables:
 !
diff --git a/src/AHFinder_output.F b/src/AHFinder_output.F
index 3f5b339..c9f4f8b 100644
--- a/src/AHFinder_output.F
+++ b/src/AHFinder_output.F
@@ -11,6 +11,167 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
+      subroutine AHFinder_InitOutput(CCTK_ARGUMENTS)
+      use AHFinder_dat
+      implicit none
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      
+      character(len=3) ::   ext
+      character(len=200) :: almf,areaf,massf,radf,circeqf,merip1f,merip2f
+      character(len=200) :: asymxf,asymyf,asymzf
+      character(len=200) :: disf
+      character(len=200) :: gaussf, gauss0f, gauss1f, gauss2f
+
+      if (myproc.ne.0) then
+        return
+      end if
+
+!     ****************************
+!     ***   BUILD FILE NAMES   ***
+!     ****************************
+
+      call CCTK_FortranString(nfile,out_dir,filestr)
+
+c     some compilers cannot trim an empty string, so we add at least one char
+      if (find3) then
+         ext = "_" // char(ichar('0') + mfind - 1) // "."
+      else
+         ext = "."
+      end if
+
+      almf = filestr(1:nfile) // "/ahf_coeff" // trim(ext) // "alm"
+      areaf = filestr(1:nfile) // "/ahf_area" // trim(ext) // "tl"
+      massf = filestr(1:nfile) // "/ahf_mass" // trim(ext) // "tl"
+      radf  = filestr(1:nfile) // "/ahf_rad" // trim(ext) // "tl"
+
+      circeqf = filestr(1:nfile) // "/ahf_circ_eq" // trim(ext) // "tl"
+      merip1f = filestr(1:nfile) // "/ahf_meri_p1" // trim(ext) // "tl"
+      merip2f = filestr(1:nfile) // "/ahf_meri_p2" // trim(ext) // "tl"
+
+      asymxf = filestr(1:nfile) // "/ahf_asymx" // trim(ext) // "tl"
+      asymyf = filestr(1:nfile) // "/ahf_asymy" // trim(ext) // "tl"
+      asymzf = filestr(1:nfile) // "/ahf_asymz" // trim(ext) // "tl"
+
+      disf = filestr(1:nfile) // "/ahf_d12" // trim(ext) // "tl"
+
+      gaussf = filestr(1:nfile) // "/ahf.gauss"
+      gauss0f = filestr(1:nfile) // "/ahf_0.gauss"
+      gauss1f = filestr(1:nfile) // "/ahf_1.gauss"
+      gauss2f = filestr(1:nfile) // "/ahf_2.gauss"
+
+
+      open(1,file=massf,form='formatted',status='replace')
+      write(1,"(A14)") '" Horizon mass'
+      close(1)
+
+      open(1,file=radf,form='formatted',status='replace')
+      write(1,"(A26)") '" Horizon radius'
+      close(1)
+
+      open(1,file=circeqf,form='formatted',status='replace')
+      write(1,"(A26)") '" Equatorial circumference'
+      close(1)
+
+      open(1,file=merip1f,form='formatted',status='replace')
+      write(1,"(A27)") '" Length of meridian, phi=0'
+      close(1)
+
+      open(1,file=merip1f,form='formatted',status='replace')
+      write(1,"(A27)") '" Length of meridian, phi=0'
+      close(1)
+
+      open(1,file=merip2f,form='formatted',status='replace')
+      write(1,"(A30)") '" Length of meridian, phi=pi/2'
+      close(1)
+
+      open(1,file=areaf,form='formatted',status='replace')
+      write(1,"(A14)") '" Horizon area'
+      close(1)
+
+      open(1,file=asymxf,form='formatted',status='replace')
+      write(1,"(A29)") '" Asymmetries on x reflection'
+      close(1)
+
+      open(1,file=asymyf,form='formatted',status='replace')
+      write(1,"(A29)") '" Asymmetries on y reflection'
+      close(1)
+
+      open(1,file=asymzf,form='formatted',status='replace')
+      write(1,"(A29)") '" Asymmetries on z reflection'
+      close(1)
+
+      open(1,file=almf,form='formatted',status='replace')
+      write(1,"(A21)")    '# Radial coefficients'
+      write(1,"(A1)")     '#'
+      close(1)
+
+      open(1,file=disf,form='formatted',status='replace')
+      write(1,"(A18)") '" Horizon distance'
+      close(1)
+
+      if (find3) then
+        call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss0f)
+        call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss1f)
+        call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss2f)
+      else
+        call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gaussf)
+      end if
+
+      end subroutine AHFinder_InitOutput
+
+
+
+      
+      subroutine AHFinder_WriteGaussHeader(CCTK_ARGUMENTS, filename)
+      use AHFinder_dat
+      implicit none
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      character(len=200) :: filename
+
+      open(1,file=filename,form='formatted', status='replace')
+      
+      write(1,"(A20)") '# GAUSSIAN CURVATURE'
+      write(1,"(A1)")  '#'
+      write(1,"(A35)") '# The data is written in a loop as:'
+      write(1,"(A1)")  '#'
+      write(1,"(A17)") '# do i=1,ntheta'
+      write(1,"(A18)") '#    do j=1,nphi'
+      write(1,"(A27)") '#       write gaussian(i,j)'
+      write(1,"(A11)") '#    end do'
+      write(1,"(A8)")  '# end do'
+      write(1,"(A1)")  '#'
+      write(1,"(A40)") '# theta and phi are subdivided uniformly'
+      write(1,"(A26)") '# according to symmetries:'
+      write(1,"(A1)")  '#'
+      write(1,"(A36)") '# phi=[0,2 pi]   (refx=refy=.false.)'
+      write(1,"(A44)") '# phi=[0,pi]     (refx=.false., refy=.true.)'
+      write(1,"(A35)") '# phi=[0,pi/2]   (refx=refy=.true.)'
+      write(1,"(A1)")  '#'
+      write(1,"(A31)") '# theta=[0,pi]   (refz=.false.)'
+      write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)'
+      write(1,"(A1)")  '#'
+
+      write(1,"(A9,L1)")  '# refx = ',refx
+      write(1,"(A9,L1)")  '# refy = ',refy
+      write(1,"(A9,L1)")  '# refz = ',refz
+      write(1,"(A1)")     '#'
+      write(1,"(A11,I7)") '# ntheta = ',(ntheta+1)
+      write(1,"(A11,I7)") '# nphi   = ',(nphi+1)
+
+      write(1,*)
+      write(1,"(A11,I4)")    '# Time step',cctk_iteration
+      write(1,"(A6,ES11.3)") '# Time',cctk_time
+      write(1,"(A6,I4)")     '# Call',ahf_ncall
+
+      close(1)
+
+      end subroutine AHFinder_WriteGaussHeader
+
+
+
 
       subroutine AHFinder_Output(CCTK_ARGUMENTS, report, status, horizon, mtype, intarea_h)
 
@@ -23,13 +184,16 @@
 
       integer mtype
       logical report, status, horizon
+      logical massf_exists, radf_exists, circeqf_exists
+      logical merip1f_exists, merip2f_exists, areaf_exists
+      logical asymxf_exists, asymyf_exists, asymzf_exists, almf_exists
       CCTK_REAL intarea_h
 
       integer l, m, ierror
-      character*3   ext
-      character*200 almf,areaf,massf,radf,circeqf,merip1f,merip2f
-      character*200 asymxf,asymyf,asymzf
-      character*200 options
+      character(len=3) ::   ext
+      character(len=200) :: almf,areaf,massf,radf,circeqf,merip1f,merip2f
+      character(len=200) :: asymxf,asymyf,asymzf
+      character(len=200) :: options
 
 
 !     ****************************
@@ -97,12 +261,13 @@ c     some compilers cannot trim an empty string, so we add at least one char
       if (myproc.eq.0) then
 
 !        Mass.
-
-        if (ahf_ncall.eq.1) then
-           open(1,file=massf,form='formatted',status='replace')
-           write(1,"(A14)") '" Horizon mass'
+        
+        inquire(file=massf, exist=massf_exists)
+        if (massf_exists) then
+          open(1,file=massf,form='formatted',status='old',position='append')
         else
-           open(1,file=massf,form='formatted',status='old',position='append')
+          open(1,file=massf,form='formatted',status='replace')
+          write(1,"(A14)") '" Horizon mass'
         end if
 
         if (status.and.report) then
@@ -116,11 +281,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 
 !        Radius.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=radf,form='formatted',status='replace')
-           write(1,"(A16)") '" Horizon radius'
+        inquire(file=radf, exist=radf_exists)
+        if (radf_exists) then
+          open(1,file=radf,form='formatted',status='old',position='append')
         else
-           open(1,file=radf,form='formatted',status='old',position='append')
+          open(1,file=radf,form='formatted',status='replace')
+          write(1,"(A26)") '" Horizon radius'
         end if
 
         if (status.and.report) then
@@ -134,11 +300,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 
 !        Equatorial circumference.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=circeqf,form='formatted',status='replace')
-           write(1,"(A26)") '" Equatorial circumference'
+        inquire(file=circeqf, exist=circeqf_exists)
+        if (circeqf_exists) then
+          open(1,file=circeqf,form='formatted',status='old',position='append')
         else
-           open(1,file=circeqf,form='formatted',status='old',position='append')
+          open(1,file=circeqf,form='formatted',status='replace')
+          write(1,"(A26)") '" Equatorial circumference'
         end if
 
         if (status.and.report) then
@@ -152,11 +319,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 
 !        Meridians.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=merip1f,form='formatted',status='replace')
-           write(1,"(A27)") '" Length of meridian, phi=0'
+        inquire(file=merip1f, exist=merip1f_exists)
+        if (merip1f_exists) then
+          open(1,file=merip1f,form='formatted',status='old',position='append')
         else
-           open(1,file=merip1f,form='formatted',status='old',position='append')
+          open(1,file=merip1f,form='formatted',status='replace')
+          write(1,"(A27)") '" Length of meridian, phi=0'
         end if
 
         if (status.and.report) then
@@ -168,11 +336,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 
         close(1)
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=merip2f,form='formatted',status='replace')
-           write(1,"(A30)") '" Length of meridian, phi=pi/2'
+        inquire(file=merip2f, exist=merip2f_exists)
+        if (merip2f_exists) then
+          open(1,file=merip2f,form='formatted',status='old', position='append')
         else
-           open(1,file=merip2f,form='formatted',status='old',position='append')
+          open(1,file=merip2f,form='formatted',status='replace')
+          write(1,"(A30)") '" Length of meridian, phi=pi/2'
         end if
 
         if (status.and.report) then
@@ -185,12 +354,13 @@ c     some compilers cannot trim an empty string, so we add at least one char
         close(1)
 
 !        Area.
-
-        if (ahf_ncall.eq.1) then
-           open(1,file=areaf,form='formatted',status='replace')
-           write(1,"(A14)") '" Horizon area'
+        
+        inquire(file=areaf, exist=areaf_exists)
+        if (areaf_exists) then
+          open(1,file=areaf,form='formatted',status='old',position='append')
         else
-           open(1,file=areaf,form='formatted',status='old',position='append')
+          open(1,file=areaf,form='formatted',status='replace')
+          write(1,"(A14)") '" Horizon area'
         end if
 
         if (status.and.report) then
@@ -206,11 +376,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 !        coefficients cs with even m indicate asymmetries
 !        on refelction x<->-x.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=asymxf,form='formatted',status='replace')
-           write(1,"(A29)") '" Asymmetries on x reflection'
+        inquire(file=asymxf, exist=asymxf_exists)
+        if (asymxf_exists) then
+          open(1,file=asymxf,form='formatted',status='old', position='append')
         else
-           open(1,file=asymxf,form='formatted',status='old',position='append')
+          open(1,file=asymxf,form='formatted',status='replace')
+          write(1,"(A29)") '" Asymmetries on x reflection'
         end if
 
         out_asymx = 0.0D0
@@ -233,11 +404,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 !        Asymmetries in y.  Any cs coefficient different from
 !        zero indicates asymmetries on reflection y<->-y.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=asymyf,form='formatted',status='replace')
-        write(1,"(A29)") '" Asymmetries on y reflection'
+        inquire(file=asymyf, exist=asymyf_exists)
+        if (asymyf_exists) then
+          open(1,file=asymyf,form='formatted',status='old', position='append')
         else
-           open(1,file=asymyf,form='formatted',status='old',position='append')
+          open(1,file=asymyf,form='formatted',status='replace')
+          write(1,"(A29)") '" Asymmetries on y reflection'
         end if
 
         out_asymy = 0.0D0
@@ -257,11 +429,12 @@ c     some compilers cannot trim an empty string, so we add at least one char
 !        or any coefficient cc and cs with odd (l-m) indicate
 !        asymmetries on reflection z<->-z.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=asymzf,form='formatted',status='replace')
-           write(1,"(A29)") '" Asymmetries on z reflection'
+        inquire(file=asymzf, exist=asymzf_exists)
+        if (asymzf_exists) then
+          open(1,file=asymzf,form='formatted',status='old',position='append')
         else
-           open(1,file=asymzf,form='formatted',status='old',position='append')
+          open(1,file=asymzf,form='formatted',status='replace')
+          write(1,"(A29)") '" Asymmetries on z reflection'
         end if
 
         out_asymz = 0.0D0
@@ -282,13 +455,14 @@ c     some compilers cannot trim an empty string, so we add at least one char
 
 !       Expansion coefficients.
 
-        if (ahf_ncall.eq.1) then
-           open(1,file=almf,form='formatted',status='replace')
-           write(1,"(A21)")    '# Radial coefficients'
-           write(1,"(A1)")     '#'
+        inquire(file=almf, exist=almf_exists)
+        if (almf_exists) then
+          open(1,file=almf,form='formatted',status='old',position='append')
+          write(1,*)
         else
-           open(1,file=almf,form='formatted',status='old',position='append')
-           write(1,*)
+          open(1,file=almf,form='formatted',status='replace')
+          write(1,"(A21)")    '# Radial coefficients'
+          write(1,"(A1)")     '#'
         end if
 
         if (status.and.report) then
diff --git a/src/AHFinder_pow.F b/src/AHFinder_pow.F
index fe53079..a928f76 100644
--- a/src/AHFinder_pow.F
+++ b/src/AHFinder_pow.F
@@ -77,7 +77,7 @@
       CCTK_REAL, dimension(1:N) :: P,PT,PTT,XIT
       CCTK_REAL, dimension(1:N,1:N) :: XI
 
-      character*200 logf
+      character(len=200) :: logf
 
       parameter (ZEPS=1.0D-10)
 
@@ -691,7 +691,7 @@
 
       CCTK_REAL, dimension(1:N) :: P
 
-      character*200 logf
+      character(len=200) :: logf
 
       save iter,logf