
!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

! XSUM_KERNELS_ASCII
!
! USAGE
!   NPROC bin/xsum_kernels_ascii INPUT_FILE OUTPUT_DIR
!
!
! COMMAND LINE ARGUMENTS
!   INPUT_FILE             - text file containing list of kernel directories
!   OUTPUT_PATH            - directory to which summed kernels are written
!
!
! DESCRIPTION
!   Sums kernels from directories specified in INPUT_FILE.
!   Writes the resulting sums to OUTPUT_DIR.
!
!   INPUT_FILE is a text file containing a list of absolute or relative paths to
!   kernel directories, one directory per line. Each kernel directory must
!   contain the following ASCII files, which are generated by SPECFEM2D when
!   adjoint simulation and ascii kernel options are specified:
!
!     proc000000_rhop_alpha_beta_kernel.dat
!     proc000000_rho_kappa_mu.dat
!
!   This routine is modeled after xsum_kernels in SPECFEM3D, in so far as
!   similarities between 2D and 3D packages allow.


program sum_kernels_ascii

  use postprocess_par, only: MAX_STRING_LEN, MAX_KERNEL_PATHS, MAX_KERNEL_NAMES, &
    IIN

  implicit none

  integer, parameter :: NARGS = 2

  integer :: ipath, npath, iker, nlines, i,j
  double precision, allocatable, dimension(:) :: kernel1, kernel2, kernel3
  double precision, allocatable, dimension(:) :: kernel_sum1, kernel_sum2, kernel_sum3
  double precision, allocatable, dimension(:,:) :: coord
  character(len=MAX_STRING_LEN) :: input_file, output_dir, kernel_name
  character(len=MAX_STRING_LEN) :: kernel_paths(MAX_KERNEL_PATHS)
  character(len=MAX_STRING_LEN) :: filename, arg(2), line
  integer :: ier
  double precision :: dummy1, dummy2


  write(*,*) 'Running XSUM_KERNELS_ASCII'
  write(*,*)

  if (command_argument_count() /= NARGS) then
    print *, 'USAGE: bin/xsum_kernels_ascii INPUT_FILE OUTPUT_DIR'
    print *, ''
    stop 'Please check command line arguments'
  endif

  ! parse command line arguments
  do i = 1, NARGS
    call get_command_argument(i,arg(i), status=ier)
  enddo

  read(arg(1),'(a)') input_file
  read(arg(2),'(a)') output_dir

  ! parse paths from INPUT_FILE
  npath=0
  open(unit = IIN, file = trim(input_file), status = 'old',iostat = ier)
  if (ier /= 0) then
     print *,'Error opening ',trim(input_file)
     stop 'Please check command line argument: INPUT_FILE'
  endif
  do while (1 == 1)
     read(IIN,'(a)',iostat=ier) line
     if (ier /= 0) exit
     npath = npath+1
     if (npath > MAX_KERNEL_PATHS) stop 'Error number of paths exceeds MAX_KERNEL_PATHS'
     kernel_paths(npath) = line
  enddo
  close(IIN)

  ! print status update
  write(*,*) 'Running XSUM_KERNELS_ASCII'
  write(*,*)

  write(*,*) 'Combining arrays from ',npath,' paths'
  write(*,*)

  ! allocate arrays
  call get_number_gll_points((kernel_paths(1)), nlines)
  allocate( coord(2,nlines) )
  allocate( kernel1(nlines) )
  allocate( kernel2(nlines) )
  allocate( kernel3(nlines) )
  allocate( kernel_sum1(nlines) )
  allocate( kernel_sum2(nlines) )
  allocate( kernel_sum3(nlines) )
  kernel1(:) = 0.0d0
  kernel2(:) = 0.0d0
  kernel3(:) = 0.0d0
  kernel_sum1(:) = 0.0d0
  kernel_sum2(:) = 0.0d0
  kernel_sum3(:) = 0.0d0

  ! sum kernels
  do iker = 1, 2

    if (iker == 1) write(kernel_name,'(a)') 'proc000000_rhop_alpha_beta_kernel.dat'
    if (iker == 2) write(kernel_name,'(a)') 'proc000000_rho_kappa_mu_kernel.dat'

    ! sum kernels
    do ipath = 1, npath
       filename = trim(kernel_paths(ipath)) //'/'// trim(kernel_name)
       open(unit=3,file=filename, status='old', action='read')
       do j = 1,nlines
          read(3,*) dummy1, dummy2, kernel1(j), kernel2(j), kernel3(j)
       enddo
       close(3)

       if (ipath == 1) then
          kernel_sum1 = kernel1
          kernel_sum2 = kernel2
          kernel_sum3 = kernel3
       else
         kernel_sum1(:)=kernel_sum1(:)+kernel1(:)
         kernel_sum2(:)=kernel_sum2(:)+kernel2(:)
         kernel_sum3(:)=kernel_sum3(:)+kernel3(:)
       endif
    enddo

    ! save result
    filename = trim(output_dir) //'/'// trim(kernel_name)
    open(unit=4,file=filename,status='unknown',action='write')
    do j = 1, nlines
          write(4,'(5e11.3)') coord(1,j),coord(2,j),kernel_sum1(j),kernel_sum2(j),kernel_sum3(j)
    enddo
    close(4)

  enddo

  write(*,*) 'Finished writing kernels '
  write(*,*)

end program sum_kernels_ascii


! ------------------------------------------------------------------------------

subroutine get_number_gll_points(kernel_path, nlines)

  use postprocess_par, only: MAX_STRING_LEN, MAX_LINES

  implicit none

  double precision :: dummy1, dummy2, dummy3, dummy4, dummy5

  character(len=MAX_STRING_LEN) :: kernel_path, filename

  integer :: j, ios, nlines

  filename = trim(kernel_path)//'/'//'proc000000_rhop_alpha_beta_kernel.dat'
  open(unit=3,file=filename,status='old',action='read')
  nlines = 0
  do j=1,MAX_LINES
     read(3,*,iostat=ios) dummy1, dummy2, dummy3, dummy4, dummy5
     if (ios /= 0) exit
     nlines=nlines+1
  enddo
  close(3)

end subroutine get_number_gll_points

