Moved the output code from the main function into a separate routine

located in AHFinder_output.F.
This routine does the AHFinder ASCII and IEEEIO output. It also invokes
an I/O method "IOAHFinderHDF5" to output AHFinder data in HDF5 file format.
This is implemented by thorn AlphaThorns/IOAHFinderHDF5.


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

1 files changed, 10 insertions, 416 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index 779d5a0..c3b4adf 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -48,7 +48,6 @@
       CCTK_REAL c0_temp(0:18)
       CCTK_REAL maxval(3),minval(3)
       CCTK_REAL rhor
-      CCTK_REAL asymx,asymy,asymz
 
       CCTK_REAL, parameter :: zero = 0.0D0
       CCTK_REAL, parameter :: half = 0.5D0
@@ -56,8 +55,7 @@
 
       CCTK_REAL, dimension (1:3) :: ahfmass
 
-      character*200 almf,logf,areaf,massf,radf,circeqf,merip1f,merip2f
-      character*200 asymxf,asymyf,asymzf
+      character*200 logf
 
       save ahfafter,ahfso,atime
 
@@ -314,93 +312,17 @@
       firstint = .true.
 
 
-!     **********************
-!     ***   FILE NAMES   ***
-!     **********************
-
-      call CCTK_FortranString(nfile,outdir,filestr)
-
-      if (find3) then
-            
-         if (mfind.eq.1) then
-
-            logf = filestr(1:nfile)//"/ahf_logfile_0"
-            almf = filestr(1:nfile)//"/ahf_coeff_0.alm"
-
-            areaf = filestr(1:nfile)//"/ahf_area_0.tl"
-            massf = filestr(1:nfile)//"/ahf_mass_0.tl"
-            radf  = filestr(1:nfile)//"/ahf_rad_0.tl"
-
-            circeqf = filestr(1:nfile)//"/ahf_circ_eq_0.tl"
-            merip1f = filestr(1:nfile)//"/ahf_meri_p1_0.tl"
-            merip2f = filestr(1:nfile)//"/ahf_meri_p2_0.tl"
-
-            asymxf = filestr(1:nfile)//"/ahf_asymx_0.tl"
-            asymyf = filestr(1:nfile)//"/ahf_asymy_0.tl"
-            asymzf = filestr(1:nfile)//"/ahf_asymz_0.tl"
-
-         else if (mfind.eq.2) then
-
-            logf = filestr(1:nfile)//"/ahf_logfile_1"
-            almf = filestr(1:nfile)//"/ahf_coeff_1.alm"
-
-            areaf = filestr(1:nfile)//"/ahf_area_1.tl"
-            massf = filestr(1:nfile)//"/ahf_mass_1.tl"
-            radf  = filestr(1:nfile)//"/ahf_rad_1.tl"
-
-            circeqf = filestr(1:nfile)//"/ahf_circ_eq_1.tl"
-            merip1f = filestr(1:nfile)//"/ahf_meri_p1_1.tl"
-            merip2f = filestr(1:nfile)//"/ahf_meri_p2_1.tl"
-
-            asymxf = filestr(1:nfile)//"/ahf_asymx_1.tl"
-            asymyf = filestr(1:nfile)//"/ahf_asymy_1.tl"
-            asymzf = filestr(1:nfile)//"/ahf_asymz_1.tl"
-
-         else
-
-            logf = filestr(1:nfile)//"/ahf_logfile_2"
-            almf = filestr(1:nfile)//"/ahf_coeff_2.alm"
-
-            areaf = filestr(1:nfile)//"/ahf_area_2.tl"
-            massf = filestr(1:nfile)//"/ahf_mass_2.tl"
-            radf  = filestr(1:nfile)//"/ahf_rad_2.tl"
-
-            circeqf = filestr(1:nfile)//"/ahf_circ_eq_2.tl"
-            merip1f = filestr(1:nfile)//"/ahf_meri_p1_2.tl"
-            merip2f = filestr(1:nfile)//"/ahf_meri_p2_2.tl"
-
-            asymxf = filestr(1:nfile)//"/ahf_asymx_2.tl"
-            asymyf = filestr(1:nfile)//"/ahf_asymy_2.tl"
-            asymzf = filestr(1:nfile)//"/ahf_asymz_2.tl"
-
-         end if
-
-      else
-
-         logf = filestr(1:nfile)//"/ahf_logfile"
-         almf = filestr(1:nfile)//"/ahf_coeff.alm"
-            
-         areaf = filestr(1:nfile)//"/ahf_area.tl"
-         massf = filestr(1:nfile)//"/ahf_mass.tl"
-         radf  = filestr(1:nfile)//"/ahf_rad.tl"
-
-         circeqf = filestr(1:nfile)//"/ahf_circ_eq.tl"
-         merip1f = filestr(1:nfile)//"/ahf_meri_p1.tl"
-         merip2f = filestr(1:nfile)//"/ahf_meri_p2.tl"
-
-         asymxf = filestr(1:nfile)//"/ahf_asymx.tl"
-         asymyf = filestr(1:nfile)//"/ahf_asymy.tl"
-         asymzf = filestr(1:nfile)//"/ahf_asymz.tl"
-
-      end if
-
-
 !     ***************************************
 !     ***   FIND OUT THE TYPE OF OUTPUT   ***
 !     ***************************************
 
+      call CCTK_FortranString(nfile,outdir,filestr)
       if (ahf_logfile.ne.0) then
          logfile = .true.
+         logf = filestr(1:nfile)//"/ahf_logfile"
+         if (find3) then
+            logf = logf // "_" // char('0' + mfind - 1)
+         end if
       else
          logfile = .false.
       end if
@@ -1473,300 +1395,6 @@
       end if
 
 
-!     ****************************
-!     ***   WRITE DATA FILES   ***
-!     ****************************
-
-      if (myproc.eq.0) then
-
-!        Expansion coefficients.
-         
-         if (ahf_ncall.eq.1) then
-            open(1,file=almf,form='formatted',
-     .      status='replace')
-            write(1,"(A21)")    '# Radial coefficients'
-            write(1,"(A1)")     '#'
-            write(1,"(A12,I4)") '# Time step ',cctk_iteration
-            write(1,"(A7,ES14.6)") '# Time ',cctk_time
-            write(1,"(A7,I4)")  '# Call ',ahf_ncall
-         else
-            open(1,file=almf,form='formatted',
-     .      status='old',position='append')
-            write(1,*)
-            write(1,"(A12,I4)") '# Time step ',cctk_iteration
-            write(1,"(A7,ES14.6)") '# Time ',cctk_time
-            write(1,"(A7,I4)")  '# Call ',ahf_ncall
-         end if
-
-         if (status.and.report) then
-            if (mtype.eq.1) then
-               if (horizon) then
-                  write(1,"(A30,I4)") '# Surface found: Outer horizon'
-               else
-                  write(1,"(A31,I4)") '# Surface found: Outer horizon?'
-               end if
-            else if (mtype.eq.2) then
-               if (horizon) then
-                  write(1,"(A30,I4)") '# Surface found: Inner horizon'
-               else
-                  write(1,"(A31,I4)") '# Surface found: Inner horizon?'
-               end if
-            else if (mtype.eq.3) then
-               if (horizon) then
-                  write(1,"(A34,I4)") '# Surface found: Marginal horizon?'
-               else
-                  write(1,"(A32,I4)") '# Surface found: Trapped surface'
-               end if
-            else if (mtype.eq.4) then
-               if (horizon) then
-                  write(1,"(A34,I4)") '# Surface found: Marginal horizon?'
-               else
-                  write(1,"(A30,I4)") '# Surface found: Not a horizon'
-               end if
-            end if
-         else
-            write(1,"(A18,I4)") '# No surface found'
-         end if
-
-         if (status.and.report) then
-
-            write(1,"(A15,3ES14.6)") '# centered on: ',xc,yc,zc
-
-            write(1,"(A1)")  '#'
-            write(1,"(A22)") '#     a_lm       l   m'
-            write(1,"(A1)")  '#'
-
-            write(1,"(ES14.6,2I4)") c0(0),0,0
-
-            do l=1,lmax
-               do m=l,1,-1
-                  write(1,"(ES14.6,2I4)") cs(l,m),l,-m
-               end do
-               write(1,"(ES14.6,2I4)") c0(l),l,0
-               do m=1,l
-                  write(1,"(ES14.6,2I4)") cc(l,m),l,m
-               end do
-            end do
-
-         end if
-
-         close(1)
-
-!        Area.
-
-         if (ahf_ncall.eq.1) then
-            open(1,file=areaf,form='formatted',
-     .      status='replace')
-            write(1,"(A14)") '" Horizon area'
-         else
-            open(1,file=areaf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,intarea_h
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Mass.
-
-         if (ahf_ncall.eq.1) then
-            open(1,file=massf,form='formatted',
-     .      status='replace')
-            write(1,"(A14)") '" Horizon mass'
-         else
-            open(1,file=massf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,
-     .      0.141047396D0*dsqrt(intarea_h)
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Radius.
-
-         if (ahf_ncall.eq.1) then
-            open(1,file=radf,form='formatted',
-     .      status='replace')
-            write(1,"(A16)") '" Horizon radius'
-         else
-            open(1,file=radf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,c0(0)
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Equatorial circumference.
-
-         if (ahf_ncall.eq.1) then
-            open(1,file=circeqf,form='formatted',
-     .      status='replace')
-            write(1,"(A26)") '" Equatorial circumference'
-         else
-            open(1,file=circeqf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,circ_eq
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Meridians.
-
-         if (ahf_ncall.eq.1) then
-            open(1,file=merip1f,form='formatted',
-     .      status='replace')
-            write(1,"(A27)") '" Length of meridian, phi=0'
-         else
-            open(1,file=merip1f,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,meri_p1
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         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'
-         else
-            open(1,file=merip2f,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,meri_p2
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Asymmetries in x.  Coefficients cc with odd m, and
-!        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'
-         else
-            open(1,file=asymxf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         asymx = 0.0D0
-
-         if (.not.refx) then
-            do l=1,lmax
-               do m=1,l
-                  if (mod(m,2).ne.0) then
-                     asymx = asymx + dabs(cc(l,m))
-                  else
-                     asymx = asymx + dabs(cs(l,m))
-                  end if
-               end do
-            end do
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,asymx/c0(0)
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        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'
-         else
-            open(1,file=asymyf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         asymy = 0.0D0
-
-         if (.not.refy) then
-            do l=1,lmax
-               do m=1,l
-                  asymy = asymy + dabs(cs(l,m))
-               end do
-            end do
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,asymy/c0(0)
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-!        Asymmetries in z. Any coefficient c0 with odd l,
-!        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'
-         else
-            open(1,file=asymzf,form='formatted',
-     .      status='old',position='append')
-         end if
-
-         asymz = 0.0D0
-
-         if (.not.refz) then
-            do l=1,lmax
-               if (mod(l,2).ne.0) asymz = asymz + dabs(c0(l))
-               do m=1,l
-                  if (mod(l-m,2).ne.0) then
-                     asymz = asymz + dabs(cc(l,m)) + dabs(cs(l,m))
-                  end if
-               end do
-            end do
-         end if
-
-         if (status.and.report) then
-            write(1,"(2ES14.6)") cctk_time,asymz/c0(0)
-         else
-            write(1,"(2ES14.6)") cctk_time,0.0D0
-         end if
-
-         close(1)
-
-      end if
-
-
 !     ******************************************************
 !     ***   PREPARING  {ahfgrid3,ahf_exp3}  FOR OUTPUT   ***
 !     ******************************************************
@@ -1775,45 +1403,11 @@
          call AHFinder_find3(CCTK_ARGUMENTS,mtype)
       end if
 
+!     ****************************
+!     ***   WRITE DATA FILES   ***
+!     ****************************
 
-!     *********************************
-!     ***   OUTPUT GRID FUNCTIONS   ***
-!     *********************************
-
-      if (ahf_2Doutput.ne.0) then
-         if (find3) then
-            if (mfind.eq.3) then
-               call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .              "ahfinder::ahfgrid3","IOFlexIO_2D","ahfgrid3")
-               call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .              "ahfinder::ahf_exp3","IOFlexIO_2D","ahf_exp3")
-            end if
-         else
-            call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .           "ahfinder::ahfgrid","IOFlexIO_2D","ahfgrid")
-            call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .           "ahfinder::ahf_exp","IOFlexIO_2D","ahf_exp")
-         end if
-      end if
-
-      if (ahf_3Doutput.ne.0) then
-
-         if (find3) then
-            if (mfind.eq.3) then
-               call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .              "ahfinder::ahfgrid","IOFlexIO_3D","ahfgrid3")
-               call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .              "ahfinder::ahf_exp","IOFlexIO_3D","ahf_exp3")
-            end if
-         else
-            call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .           "ahfinder::ahfgrid","IOFlexIO_3D","ahfgrid")
-            call CCTK_OutputVarAsByMethod(ierror,cctkGH,
-     .           "ahfinder::ahf_exp","IOFlexIO_3D","ahf_exp")
-         end if
-
-
-      end if
+      call AHFinder_Output (CCTK_ARGUMENTS, report, status, horizon, mtype, intarea_h)
 
 
 !     ****************