text:
module test_function_object
  ! Object oriented approach for function
  ! evaluation.

  use utils, only : wp
  use num_integration

  implicit none

  abstract interface
    function f_interface(x)
      import :: wp
      real(wp), intent(in) :: x
      real(wp) :: f_interface
    end function f_interface
  end interface

  ! function_object is an object that can point
  ! to any function that follow the f_interface.
  ! This allows for the object to change its behavior
  ! in runtime.
  type, extends(fun) :: function_object
    procedure(f_interface), pointer, nopass :: f
  contains
    procedure :: eval => specialized_function
  end type function_object

  contains

  function specialized_function(self, x)
    class(function_object), intent(in) :: self
    real(wp), intent(in) :: x
    real(wp) :: specialized_function
    specialized_function = self%f(x)
  end function specialized_function

end module test_function_object

program romgerg_tests
use utils
use num_integration
use test_function_object

! Note that not every test is meant to pass,
! depending on the chosen number of iterations.
! Estimations by Romberg's using a low number
! of iterations should not give a decent
! approximation for higher order polynomials
! or complex functions.

implicit none

!---------------
! Variables
!---------------

integer :: i,j,l,m,n,total_tests,passed_tests
character (len=:), allocatable :: expression
real(wp) :: a(3),b(3)
integer :: k(3)
type(function_object) :: fobj

!---------------
! Logic
!---------------

total_tests = 1
passed_tests = 0
k = [5, 10, 15]
a = [-2.0, 0.0, 1.0]
b = [1.0, 1.0, 5.0]

write(*,*)
write(*,*) "!-------------------------------"
write(*,*) "! Testing integration by Romberg"
write(*,*) "!-------------------------------"
write(*,*)

do i = 1, size(k)
  write(*,*)
  write(*,*) "!----------------------------------"
  write(*,*) "! With k = ",k(i)
  write(*,*) "!----------------------------------"
  do j = 1, size(a)
    expression = trim("f(x)=1")
    fobj%f => poly0
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    expression = trim("f(x)=x+1")
    fobj%f => poly1
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    expression = trim("f(x)=3*x**2 + 2")
    fobj%f => poly2
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    expression = trim("f(x)=33.33*x**5 + 8.8*x**4 + x**3 + 4.5*x*2 + 2.123")
    fobj%f => poly5
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    expression = trim("f(x)=8*x**7/77 + (4*sin(2*x*pi) + 9*cos(11*x*pi))/5")
    fobj%f => ftrig
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    expression = trim("f(x)=e^x")
    fobj%f => fe
    call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)

    if ( a(j) > 0 ) then
      expression = trim("f(x)=ln(x)")
      fobj%f => fln
      call test_integration(fobj,a(j),b(j),k(i),expression,passed_tests,total_tests,romberg,close)
    end if
  end do
  write(*,*)
end do

write(*,*) "!---------------------------"
write(*,'(a,g14.6)') " !Total tests ", total_tests-1
write(*,'(a,g14.6)') " !Tests passed", passed_tests
write(*,*) "!---------------------------"

contains

  function poly0(x)
    real(wp), intent(in) :: x
    real(wp) :: poly0

    poly0 = 1.0_wp
  end function poly0

  function poly1(x)
    real(wp), intent(in) :: x
    real(wp) :: poly1

    poly1 = x + 1
  end function poly1

  function poly2(x)
    real(wp), intent(in) :: x
    real(wp) :: poly2

    poly2 = 3*x**2 + 2
  end function poly2

  function poly5(x)
    real(wp), intent(in) :: x
    real(wp) :: poly5

    poly5 = 33.33*x**5 + 8.8*x**4 + x**3 + 4.5*x*2 + 2.123
  end function poly5

  function ftrig(x)
    real(wp), intent(in) :: x
    real(wp) :: ftrig

    ftrig = 8*x**7/77 + (4*sin(2*x*pi) + 9*cos(11*x*pi))/5
  end function ftrig

  function fe(x)
    real(wp), intent(in) :: x
    real(wp) :: fe

    fe = exp(x)
  end function fe

  function fln(x)
    real(wp), intent(in) :: x
    real(wp) :: fln

    fln = log(x)
  end function fln

  subroutine test_integration(fobj,a,b,k,expression,passed_tests,total_tests,romberg,close)
    real(wp), intent(in) :: a,b
    real(wp) :: true,result,abserr
    real(wp), parameter :: epsabs = 0.0D+00
    real(wp), parameter :: epsrel = 0.001D+00
    integer, intent(inout) :: passed_tests, total_tests
    integer, intent(in) :: k
    integer, parameter :: key = 6
    integer, parameter :: limit = 500
    integer, parameter :: lenw = 4 * limit
    integer :: ier
    integer :: iwork(limit)
    integer :: last
    integer :: neval
    real(wp) :: work(lenw)

!    real(wp), external :: f
!    real(wp), external :: fcomp
    class(function_object), intent(in) :: fobj
    real(wp), external :: romberg
    logical, external :: close

    character (len=:), allocatable, intent(in) :: expression

    write(*,*)
    write(*,*) "Test number ", total_tests
    total_tests = total_tests + 1
    result = romberg(fobj, a, b, k)
    call dqag ( fobj%f , a, b, epsabs, epsrel, key, true, abserr, neval, ier,&
    limit, lenw, last, iwork, work )

    write (*, '(A,A)') '  Integrand is ', expression
    write ( *, '(A,F15.5)' ) '  Integral left endpoint A     ', a
    write ( *, '(A,F15.5)' ) '  Integral right endpoint B    ', b
    write (*, '(A,I15)') '  #Romberg iterations(k) is', k
    write ( *, '(A,F15.5)' ) '  Integral is                  ', true
    write ( *, '(A,F15.5)' ) '  Estimated integral is        ', result
    if ( close(result,true,1.0_wp*1e-5) ) then
      passed_tests = passed_tests + 1
      write(*,*) "!Test passed!"
    else
      write(*,*) "!-----------!"
      write(*,*) "!Test failed!"
      write(*,*) "!-----------!"
    end if
    write(*,*)

  end subroutine test_integration
end program romgerg_tests

max_stars_repo_name: AndreTGMello/numerical-analysis-course

max_issues_repo_name: AndreTGMello/numerical-analysis-course

max_forks_repo_name: AndreTGMello/numerical-analysis-course


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!> @brief parallelisation wrapper for "omp_old_save"
!> @param[in] level level number
!> @param[in] ismpl local sample index
      subroutine old_save(level,ismpl)

        use mod_OMP_TOOLS

        implicit none

        integer :: ismpl

        include 'OMP_TOOLS.inc'

        integer, intent (in) :: level

#ifdef fOMP
!$OMP parallel num_threads(Tlevel_1)
!$      call omp_binding(abs(ismpl))
#endif

        call omp_old_save(level,ismpl)

#ifdef fOMP
!$OMP end parallel
#endif
!
        return

      end subroutine old_save

!> @brief save current state as an old version
!> @param[in] level level number (which old version)
!> @param[in] ismpl local sample index
      subroutine omp_old_save(level,ismpl)

        use arrays
        use mod_genrl
        use mod_conc

        implicit none

        ! local sample index
        integer :: ismpl

        ! cgen level index
        integer, intent (in) :: level

        ! directional cell-indices
        integer :: i, j, k

        ! counter for concentration
        integer :: l

        ! Start position in array for process
        integer :: tpos

        ! Number of array elements for process
        integer :: tanz

        ! OpenMP partition, get tpos, tanz
        call omp_part(i0*j0*k0,tpos,tanz)

        ! Get i, j, k indices for tpos
        call ijk_m(tpos,i,j,k)

        ! Copy all variable arrays to old: var -> varold
        call dcopy(tanz,head(i,j,k,abs(ismpl)),1, headold(tpos,level,max(1,ismpl)), 1)
        call dcopy(tanz,temp(i,j,k,abs(ismpl)),1, tempold(tpos,level,max(1,ismpl)), 1)
        call dcopy(tanz,pres(i,j,k,abs(ismpl)),1, presold(tpos,level,max(1,ismpl)), 1)
        do l = 1, ntrans
          call dcopy(tanz,conc(i,j,k,l,abs(ismpl)),1, concold(tpos,l,level,max(1,ismpl)),1)
        end do

        return

      end subroutine omp_old_save


text:
UC Davis The university has commencement (graduation) ceremonies twice a year. Once in June after spring quarter and once in December after fall quarter. The commencement ceremonies are split up into multiple separate events for different colleges and even majors.

  Bars In Downtown Davis really like the June graduation.

 How to Graduate
After a student has http://registrar.ucdavis.edu/graduation/ filed for graduation, they may choose to participate in the commencement ceremony, for which they have to http://commencement.ucdavis.edu/ register also.

 Tickets
Students are mailed guest tickets in midMay.  The actual number received depends on the college and how many people are graduating. If available, extra ticket(s) can obtained with more details in the letter sent home to students.  Only guests need tickets, graduates do not need a ticket to get in.

Historically, the College of Letters and Science Spring Commencement has given five to six guest tickets per student, whereas the College of Ag gives closer to ten tickets.

Although university policy forbids selling tickets, giving them away is fine, so if you have extras, or need more, head over to UCD Commencement Ticket Swap ticket swap and help someone out!

 Doing it in style
While some people show up in jeans and sneakers, they realize that theyre a little underdressed. Old sneakers stand out when youre walking across the stage so they should be avoided if you care about things like that. A decent percentage of people decorate their hats. This is totally allowable because you get to keep it after you graduate. People do everything from their name in glitter to lights, toy Unitrans Buses and cows. Whether you choose to decorate your hat or wear a colorful flower lei, the idea is that youll stand out so that your family in the stands can more easily see you.
 General Commencement Information
Nearly all of the information that you need to know is covered in the http://commencement.ucdavis.edu/ Commencement Handbook.

The large commencement ceremonies take place indoors at The Pavilion at ARC. Parking can be found in the lots around the event center for $6 or for free if you park further away.

 Past Commencements
 Fall 2005
The Fall 2005 Letters and Science Commencement ceremony was held at 2pm on Sunday, December 18th. It was a very rainy day and The Pavilion at ARC The Pavilion was at nearly full capacity. The student speaker was Sapana Shende. Little of interest occurred. The only notable thing was the graduate who took a face plant on stage but he probably did it on purpose.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! Truong Dang
! 01-11-2016
program Laplace_2D
    implicit none
    interface
       function functionf2D(x,y,k)
          real(8)::functionf2D,y,x
          integer(4)::k
       end function functionf2D

       function exact_solution2D(x,y,k)
          real(8)::exact_solution2D,y,x
          integer(4)::k
       end function exact_solution2D
    end interface

    integer(4)::N,k1,j,i,t
    real(8)::aa,bb
    real(8)::h,k
    real(8),dimension(4)::M
    real(8),dimension(:),allocatable::x,xcp,y,ycp,F,u
    real(8),dimension(:,:),allocatable::B,A,C,D,B_temp,B_temp_inv
    real(8),dimension(:,:),allocatable::u_dis,u_exact

    real(8),dimension(4)::norm_l2,norm_l1

    real(8)::ai,bi,cj,dj,sij
    integer(4)::i1,i2

    N=2
    k1=3

    do j=1,4
       aa=0.0d0
       bb=1.0d0
       N=N*2
       allocate(x(1:N))
       allocate(xcp(1:N+1))
       allocate(y(1:N))
       allocate(ycp(1:N+1))
       allocate(B(1:(N-1)*(N-1),1:(N-1)*(N-1)))
       allocate(B_temp(1:(N-1)*(N-1),1:(N-1)*(N-1)))
       allocate(B_temp_inv(1:(N-1)*(N-1),1:(N-1)*(N-1)))
       allocate(A(1:N-1,1:N-1))
       allocate(C(1:N-1,1:N-1))
       allocate(D(1:N-1,1:N-1))
       allocate(F(1:(N-1)*(N-1)))
       allocate(u(1:(N-1)*(N-1)))
       allocate(u_dis(1:N+1,1:N+1))
       allocate(u_exact(1:N+1,1:N+1))

       M(j)=N

    !%%Make mesh points x(i+1/2)
    h=1.0/(N-1)
    x(1)=aa
    do i=2,N-1
        x(i)=(i-1)*h + x(1)
    enddo
    x(N)=bb

    !%%Make control points xcp(i)
    xcp=0.0d0
    xcp(1)=x(1)
    do i=2,N+1
        if (i==N+1) then
            xcp(i)=x(i-1)
        else
            xcp(i)=(x(i)+x(i-1))/2.0d0;
        endif
    enddo

    !%%Make mesh points y(i+1/2) i=0,1,...,N
    k=1.0d0/(N-1)
    y(1)=aa
    do i=2,N-1
        y(i)=(i-1)*k + y(1)
    enddo
    y(N)=bb

    !%%Make control points ycp(i) i=0,1,...,N+1
    ycp=0.0d0
    ycp(1)=y(1)
    do i=2,N+1
        if (i==N+1) then
            ycp(i)=y(i-1)
        else
            ycp(i)=(y(i)+y(i-1))/2.0d0
        endif
    enddo
    !%%

    B=0.0d0

    !%%
    ai=-1.0d0/h**2
    bi=-1.0d0/h**2
    cj=-1.0d0/k**2
    dj=-1.0d0/k**2
    sij=(ai+bi+cj+dj)
    !%%
    A=0.0d0
    do i=1,N-1
        if(i==1) then
            A(i,i)=sij
            A(i,i+1)=-bi
            elseif(i==N-1)then
                A(i,i)=sij
                A(i,i-1)=-ai
            else
                A(i,i)=sij
                A(i,i-1)=-ai
                A(i,i+1)=-bi
            endif

    enddo
    !%%

        C=0.0d0
        do i=1,N-1
            C(i,i)=-cj
        enddo

    !%%

        D=0.0d0
        do i=1,N-1
            D(i,i)=-dj
        enddo

    !%%

    do i=1,N-1
        if(i==1)then
            B((i-1)*(N-1)+1:i*(N-1),(i-1)*(N-1)+1:i*(N-1))=-A
            B((i-1)*(N-1)+1:i*(N-1),i*(N-1)+1:(i+1)*(N-1))=-D
        elseif(i==N-1)then
                B((i-1)*(N-1)+1:i*(N-1),(i-1)*(N-1)+1:i*(N-1))=-A
                B((i-1)*(N-1)+1:i*(N-1),(i-2)*(N-1)+1:(i-1)*(N-1))=-C
            else
                B((i-1)*(N-1)+1:i*(N-1),(i-1)*(N-1)+1:i*(N-1))=-A
                B((i-1)*(N-1)+1:i*(N-1),i*(N-1)+1:(i+1)*(N-1))=-D
                B((i-1)*(N-1)+1:i*(N-1),(i-2)*(N-1)+1:(i-1)*(N-1))=-C

        endif
    enddo
    !%A=(1/h^2)*A;
    F=0.0d0
    do i1=1,N-1
        do i2=1,N-1
            F((i1-1)*(N-1)+i2)=functionf2D(xcp(i2),ycp(i1),k1)
        enddo
    enddo
    !u=B\F;
    B_temp = B
    call gauss_2(B,F,u,(N-1)*(N-1))
    call matinv(B_temp,(N-1)*(N-1),B_temp_inv)
    !print*,B_temp_inv

    !do i1=1,(N-1)*(N-1)
    !        u(i1) = (1.0D0/B(i1,i1))*(F(i1) &
    !        -dot_product(B(i1,1:(N-1)*(N-1)),u)+B(i1,i1) &
    !        *u(i1))
    !enddo

    u_dis=0.0d0
    do i1=1,N-1
        do i2=1,N-1
            u_dis(i1+1,i2+1)=u((i1-1)*(N-1)+i2)
        enddo
    enddo

    open(unit=1,file='u_dis')
    do i1=1,N+1
        do i2=1,N+1
            write(1,*),xcp(i1),ycp(i2),u_dis(i1,i2)
        enddo
        write(1,*),''
    enddo
    close(unit=1)

    u_exact=0.0d0
    do i1=1,N+1
        do i2=1,N+1
            u_exact(i1,i2)=exact_solution2D(xcp(i2),ycp(i1),k1)
        enddo
    enddo

    open(unit=2,file='u_exact')
    do i1=1,N+1
        do i2=1,N+1
            write(2,*),xcp(i1),ycp(i2),u_exact(i1,i2)
        enddo
        write(2,*),''
    enddo
    close(unit=2)
    !figure
    !surf(xcp,ycp,u_dis)
    !title('u discrete', 'FontWeight','bold')
    !if(k1==2)
    !    zlim([0 0.08])
    !else
    !    zlim([0 0.025])
    !end

    !figure
    !surf(xcp,ycp,u_exact)
    !title('u exact', 'FontWeight','bold')
    !%% Odd indices figures: u discrete
    !%% Even indices figures: u exact

    !%%
    norm_l2(j)=0.0d0
    norm_l1(j)=0.0d0

    do i=1,N+1
        do t=1,N+1
        norm_l2(j)=norm_l2(j)+(u_dis(i,t)-u_exact(i,t))**2*h*k
        norm_l1(j)=norm_l1(j)+abs(u_dis(i,t)-u_exact(i,t))*h*k
        enddo
    enddo
    norm_l2(j)=sqrt(norm_l2(j))
    print*,'Mesh size:',N,'x',N
    print*,'L2 norm = ',norm_l2(j)
    print*,'L1 norm = ',norm_l1(j)

!    %%
!    %%
!    norm_h1(j)=0;
!    %norm_h1(j)=norm_h1(j)+(((u_dis(1)-u_exact(1))))^2;
!     for i=1:N
!         for t=1:N
!             if (i==1)&&(t==1)
!                 norm_h1(j)=norm_h1(j)+4*(((u_dis(i,t)-u_exact(i,t))))^2;
!             else
!                 norm_h1(j)=norm_h1(j)+(abs(u_dis(i+1,t)-u_exact(i+1,t))-abs(u_dis(i,t)-u_exact(i,t)))^2+(abs(u_dis(i,t+1)-u_exact(i,t+1))-abs(u_dis(i,t)-u_exact(i,t)))^2;
!             end
!         end
!     end
!     %norm_h1(j)=norm_h1(j)+(((u_dis(N)-u_exact(N))-(u_dis(N-1)-u_exact(N-1))))^2;
!     norm_h1(j)=(norm_h1(j))^(1/2);
!
!    %%

!figure
!%plot( log(M), -log(norm_l2), 'red', log(M), 2*log(M),'black', log(M), -log(norm_h1), 'cyan', log(M), (1/2)*log(M), 'green')
!%ylim([0 20])
!%legend('|| ||_{L^2}','2*Log(N)','|| ||_{H^1_0}','(1/2)*Log(N)')
!%legend('Location','northwest')
!plot( log(M), -log(norm_l2))
!hold on
!plot( log(M), 2*log(M),'r')
!xlabel("Log(N)")
!ylabel("-Log(err)")
    deallocate(x)
    deallocate(xcp)
    deallocate(y)
    deallocate(ycp)
    deallocate(B)
    deallocate(B_temp)
    deallocate(B_temp_inv)
    deallocate(A)
    deallocate(C)
    deallocate(D)
    deallocate(F)
    deallocate(u)
    deallocate(u_dis)
    deallocate(u_exact)
enddo

open(unit=3,file='L2 norm')
do j=1,4
    write(3,*),log(M(j)),log(norm_l2(j))
enddo
close(unit=3)

open(unit=4,file='L1 norm')
do j=1,4
    write(4,*),log(M(j)),log(norm_l1(j))
enddo
close(unit=4)

end program

function functionf2D(x,y,k)
   implicit none
   real(8)::functionf2D,y,x
   integer(4)::k

if (k==1) then
    functionf2D =2
endif

if (k==2) then
    functionf2D =2*y*(1-y)+2*x*(1-x)
endif

if (k==3) then
    functionf2D =2*x**2*(2*y - 2)*(x - 1) + 2*x**2*y*(x - 1) + 4*x*y*(y - 1)**2 + 2*y*(x - 1)*(y - 1)**2
endif

!if (k==4)
!    f=2*y^2*(2*x - 2)*(y - 1) + 4*x*y*(x - 1)^2 + 2*x*y^2*(y - 1) + 2*x*(x - 1)^2*(y - 1);
!end

if (k==5) then
    functionf2D =0.0d0
endif

end function functionf2D

function exact_solution2D(x,y,k)
    implicit none
    real(8)::exact_solution2D,y,x
    integer(4)::k
if (k==1) then
    exact_solution2D =x*(1-x)
endif
!% Truong added 1-11-2016
!% Neu truyen k=1, ham chi phu thuoc x nen nghiem xap xi ra sai

if(k==2)then
    exact_solution2D=x*(1-x)*y*(1-y)
endif

if (k==3) then
    exact_solution2D =x**2*(1-x)*y*(1-y)**2
endif

!if (k==4) then
!    u_ex=x*(1-x)^2*y^2*(1-y)
!endif

if (k==5) then
    exact_solution2D = sin(x)*exp(y)
endif
end function exact_solution2D


subroutine interp1(x,v,n1,n2,C,NZ)
   implicit none;
 integer*4 i,n1,n2;
   double precision,dimension(n1):: x,v;
 double precision,dimension(n2):: C,NZ;
 double precision,dimension(1)::NZ1,P1,P2,P3,V1,V2,V3;
 do i=1,n2
 P1=x(minloc(abs(x-C(i)))-1);P2=x(minloc(abs(x-C(i))));
 P3=x(minloc(abs(x-C(i)))+1);
 V1=v(minloc(abs(x-C(i)))-1);V2=v(minloc(abs(x-C(i))));
 V3=v(minloc(abs(x-C(i)))+1);
if (C(i)==x(i)) then
NZ1(1)=v(i);
else
 if (C(i)>P2(1)) then
  NZ1(1)=V2(1)+(V3(1)-V2(1))/(P3(1)-P2(1))*(C(i)-P2(1));
  else
  NZ1(1)=V2(1)+(V2(1)-V1(1))/(P2(1)-P1(1))*(C(i)-P2(1));
  end if;
end if;
  NZ(i)=NZ1(1);
  end do ;

   return;
end subroutine


text:
!> @ingroup insertion_sort_asc
!>
!> @brief
!> Run an insertion sort on a vector of double precision values.
!>
!> @details
!> Run an ascending on a vector of double precision values to sort them in the specified order.
!>
!> @param vec
!> The vector double precision values to be sorted.
!>
!> @author Arin R. Bratt
!> @date 11/07/2014
pure subroutine insertion_sort_asc_dble(vec)
  ! Arguments
  double precision, intent(inout) :: vec(:)

  ! Variables
  double precision :: double_precision_value
  double precision :: swap_value
  integer :: i
  integer :: j

  ! Sort to the first element
  if (vec(size(vec,1)) > vec(1)) then
     double_precision_value = vec(1)
  else
     double_precision_value = vec(size(vec,1))
     vec(size(vec,1)) = vec(1)
  end if

  do i = size(vec,1) - 1,2,-1
     if (double_precision_value > vec(i)) then
        ! Swap double precision values double_precision_value and vec(i).
        swap_value = double_precision_value
        double_precision_value = vec(i)
        vec(i) = swap_value
     end if
  end do

  ! Run regualr insertion sort
  do j = 3,size(vec,1)
     swap_value = vec(j)
     i = j - 1
     if (vec(i) > swap_value) then
        vec(j) = vec(i)
        i = i - 1
        do
           if (.not. vec(i) > swap_value) exit
           vec(i + 1) = vec(i)
           i = i - 1
        end do
        vec(i + 1) = swap_value
     end if
  end do
end subroutine insertion_sort_asc_dble

max_stars_repo_name: arinrb/bsort

max_issues_repo_name: arinrb/bsort

max_forks_repo_name: arinrb/bsort


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine complexfrequencys(inpc,textpart,nmethod,
     &  mei,iperturb,istep,istat,n,iline,ipol,inl,
     &  ipoinp,inp,ithermal,xboun,nboun,ipoinpc,mcs,cs,cyclicsymmetry,
     &  ier)
!
!     reading the input deck: *COMPLEX FREQUENCY
!
      implicit none
!
      character*1 inpc(*)
      character*132 textpart(16)
!
      integer nmethod,mei(4),istep,istat,iperturb(*),i,nboun,
     &  n,key,iline,ipol,inl,ipoinp(2,*),inp(3,*),nev,ithermal(*),
     &  ipoinpc(0:*),mcs,cyclicsymmetry,ier
!
      real*8 xboun(*),cs(17,*)
!
      mei(4)=0
!
      if(istep.lt.1) then
         write(*,*) '*ERROR reading *COMPLEX FREQUENCY:'
         write(*,*) '       *COMPLEX FREQUENCY can only be used'
         write(*,*) '       within a STEP'
         ier=1
         return
      endif
!
!     no heat transfer analysis
!
      if(ithermal(1).gt.1) then
         ithermal(1)=1
      endif
!
!     check for cyclic symmetry
!
      if((mcs.ne.0).and.(cs(2,1).ge.0.d0)) then
         cyclicsymmetry=1
      endif
!
      nmethod=0
      do i=2,n
         if(textpart(i)(1:8).eq.'CORIOLIS') then
            nmethod=6
         elseif(textpart(i)(1:7).eq.'FLUTTER') then
            nmethod=7
         elseif(textpart(i)(1:11).eq.'STORAGE=YES') then
            write(*,*) '*WARNING reading *COMPLEX FREQUENCY:'
            write(*,*) '         for this keyword'
            write(*,*) '         STORAGE=YES is deactivated'
            write(*,*) '         in the CalculiX code'

c            mei(4)=1
         else
            write(*,*)
     &           '*WARNING reading *COMPLEX FREQUENCY:'
            write(*,*) '         parameter not recognized:'
            write(*,*) '         ',
     &           textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*COMPLEX FREQUENCY%")
         endif
      enddo
      if(nmethod.eq.0) then
         write(*,*)
     &        '*ERROR reading *COMPLEX FREQUENCY:'
         write(*,*) '       either parameter CORIOLIS'
         write(*,*) '       or parameter FLUTTER is required'
         call inputerror(inpc,ipoinpc,iline,
     &        "*COMPLEX FREQUENCY%",ier)
         return
      endif
!
      if(iperturb(1).gt.1) iperturb(1)=0
      iperturb(2)=0
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
      if((istat.lt.0).or.(key.eq.1)) then
         write(*,*) '*ERROR reading *COMPLEX FREQUENCY:'
         write(*,*) '       definition not complete'
         write(*,*) '  '
         call inputerror(inpc,ipoinpc,iline,
     &        "*COMPLEX FREQUENCY%",ier)
         return
      endif
      read(textpart(1)(1:10),'(i10)',iostat=istat) nev
      if(istat.gt.0) then
         call inputerror(inpc,ipoinpc,iline,
     &        "*COMPLEX FREQUENCY%",ier)
         return
      endif
      if(nev.le.0) then
         write(*,*) '*ERROR reading *COMPLEX FREQUENCY:'
         write(*,*) '       less than 1 eigenvalue requested'
         ier=1
         return
      endif
!
      mei(1)=nev
!
!     removing nonzero boundary conditions
!
      do i=1,nboun
         xboun(i)=0.d0
      enddo
!
!     correction for cyclic symmetric structures:
!     if the present step was not preceded by a frequency step
!     no nodal diameter has been selected. To make sure that
!     mastructcs is called instead of mastruct a fictitious
!     minimum nodal diameter is stored
!
      if((cyclicsymmetry.eq.1).and.(mcs.ne.0).and.(cs(2,1)<0.d0))
     &       cs(2,1)=0.d0
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &     ipoinp,inp,ipoinpc)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text: select FLATTEN(t.a.b.c) AS c from dfs.`/drill/testdata/json_kvgenflatten/drill3562` t;


text:
MODULE mo2

  USE kind, ONLY: i8, dp

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: arr2, structype, struc, &
            alloc_arr2, alloc_strucarr, dealloc_arr2, dealloc_strucarr

  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: arr2

  TYPE structype
     REAL(dp), DIMENSION(:,:), ALLOCATABLE :: arr
  END TYPE structype
  TYPE(structype) :: struc

CONTAINS

  ! -------------------------------------------------------
  SUBROUTINE alloc_arr2(nx,ny)

    IMPLICIT NONE

    INTEGER(i8), INTENT(IN) :: nx, ny

    if (.not. allocated(arr2)) allocate(arr2(nx,ny))

  END SUBROUTINE alloc_arr2

  ! -------------------------------------------------------
  SUBROUTINE alloc_strucarr(nx,ny)

    IMPLICIT NONE

    INTEGER(i8), INTENT(IN) :: nx, ny

    if (.not. allocated(struc%arr)) allocate(struc%arr(nx,ny))

  END SUBROUTINE alloc_strucarr

  ! -------------------------------------------------------
  SUBROUTINE dealloc_arr2()

    IMPLICIT NONE

    if (allocated(arr2)) deallocate(arr2)

  END SUBROUTINE dealloc_arr2

  ! -------------------------------------------------------
  SUBROUTINE dealloc_strucarr()

    IMPLICIT NONE

    if (allocated(struc%arr)) deallocate(struc%arr)

  END SUBROUTINE dealloc_strucarr

END MODULE mo2

max_stars_repo_licenses: MIT Intel Unlicense NetCDF

max_issues_repo_licenses: MIT Intel Unlicense NetCDF

max_forks_repo_licenses: MIT Intel Unlicense NetCDF


text: #redirect Calgene LLC


text:
program main
    implicit none

    ! raw data information
    integer :: xdim  ! x-dimension of raw array
    integer :: ydim  ! y-dimension of raw array
    integer :: zdim  ! z-dimension of raw array
    integer :: xdelta ! number of points in each tile
    integer :: ydelta ! number of points in each tile
    real    :: known_lat ! the lat of bottom-left point
    real    :: known_lon ! the lon of bottom-left point
    real    :: dx ! the delta x in deg
    real    :: dy ! the delta y in deg

    integer :: signed   ! 0=unsigned data, 1=signed data
    integer :: endian   ! 0=big endian, 1=little endian
    integer :: urban    ! the id of urban
    integer :: wordsize ! number of bytes to use for each array element
    real    :: scale    ! value to multiply array elements by before truncation to integers

    ! new data
    integer :: newnx   ! x-dimension of raw array
    integer :: newny   ! y-dimension of raw array

    ! raw data
    real,allocatable,dimension(:) :: lats ! lats of raw data
    real,allocatable,dimension(:) :: lons ! lons of raw data
    real,allocatable,dimension(:,:,:) :: modis ! raw data

    real,allocatable,dimension(:,:) :: newlats ! lats of raw data
    real,allocatable,dimension(:,:) :: newlons ! lons of raw data
    real,allocatable,dimension(:,:) :: newdata ! raw data

    ! landindex
    integer :: oldindex = 13
    integer :: newindex = 13

    ! file
    character(len=200) :: oldfile
    character(len=200) :: newfile
    character(len=200) :: outfile

    character(len=128) :: oldpath
    character(len=128) :: newpath
    character(len=128) :: outpath

    ! local
    integer :: i,j,k
    integer :: status, strlen
    real    :: xleft, yleft
    real    :: maxlat,maxlon,minlat,minlon
    integer :: xstart, xend, ystart, yend
    character(len=23) :: fname
    ! function
    character(len=5),external :: integerToCharacter
    integer, external :: system

    ! Declare the local variables
    integer                      :: fnamelen    ! The len of the name of input file.
    integer                      :: ierror      ! Flag for open file error
    integer                      :: inputunit=9 ! Open unit
    character(len=100)           :: inputname   ! The name of input file
    character(len=200)           :: errmsg      ! error message

    namelist /description/ xdim, ydim, zdim, xdelta, ydelta, known_lat, known_lon, &
                    dx, dy, signed, endian, urban, wordsize, scale, newnx, newny, &
                    oldindex, newindex

    namelist /path/ oldpath, newfile, outpath

 ! write(*,*) "=========================================================="
 ! write(*,*) "          Copyright (C) 2018, Linhong Xiao"
 ! write(*,*) "         updateLanduse version 2.0.1 2018.05.02         "
 ! write(*,*) "=========================================================="

 !  read the input file name
 if (command_argument_count() == 0) then
     write(*,*) '-------------------parameter file name--------------------'
     write(*,*) 'Notice: The input file used is input.up in local folder. '
     write(*,*) '----------------------------------------------------------'
     write(*,*) '                                                             '
     inputname = 'input.up'
 elseif (command_argument_count() > 0) then
     call get_command_argument(1,inputname,fnamelen,ierror)
     ! write(*,*) fnamelen
     if (fnamelen < 100) then
        write(*,*) '-------------------parameter file name--------------------'
        write(*,*) 'Notice: The input file used is ' // trim(inputname)
        write(*,*) '----------------------------------------------------------'
        write(*,*) '                                                          '
     else
        write(*,*) '****************************************************'
        write(*,*) '*** Error: The name of input file is too LONG!   ***'
        write(*,*) '*** The length of the input name is:',fnamelen, '***'
        write(*,*) '*** Please keep the LEN less than 100 !          ***'
        write(*,*) '****************************************************'
        stop
     end if
 end if

 ! open input file and deteminte if exists
  open(inputunit,file=inputname,form='formatted',iostat=ierror,status='old',iomsg=errmsg)
  if(ierror /= 0)then
      write(*,*) errmsg
      write(*,*) inputname
      stop
  end if

  read(inputunit,NML=description,iostat=ierror,iomsg=errmsg )
   if(ierror/=0) write(*,*) errmsg
   if(ierror/=0) stop

  read(inputunit,NML=path,iostat=ierror,iomsg=errmsg )
   if(ierror/=0) write(*,*) errmsg
   if(ierror/=0) stop
  close(inputunit)

    allocate( lats(ydelta),lons(xdelta),modis(xdelta, ydelta, zdim) )
    allocate( newlats(newnx, newny),newlons(newnx, newny),newdata(newnx, newny) )
    call readnew(newfile, newnx, newny, newlats, newlons, newdata)

    ! handle the core operations
    do i=1,xdim/xdelta
        xleft=known_lon+(i-1)*xdelta*dx
        call createlist(xleft, dx, xdelta, lons)

        do j=1,ydim/ydelta
           yleft=known_lat+(j-1)*ydelta*dy
           call createlist(yleft, dy, ydelta, lats)

           xstart = (i-1)*xdelta+1
           xend   = i*xdelta
           ystart = (j-1)*ydelta+1
           yend   = j*ydelta
           fname  = integerToCharacter(xstart)//"-"//integerToCharacter(xend)// &
                    "."//integerToCharacter(ystart)//"-"//integerToCharacter(yend)
           ! read the raw data
           write(*,*) fname
           oldfile = trim(oldpath)//"/"//trim(fname)
           call countLen(oldfile,strlen)
           call read_geogrid(oldfile, strlen, modis, xdelta, ydelta, zdim,&
                              signed, endian, scale, wordsize, status)
                ! char * fname,        /* The name of the file to read from */
                ! int * len,           /* The length of the filename */
                ! float * rarray,      /* The array to be filled */
                ! int * nx,            /* x-dimension of the array */
                ! int * ny,            /* y-dimension of the array */
                ! int * nz,            /* z-dimension of the array */
                ! int * isigned,       /* 0=unsigned data, 1=signed data */
                ! int * endian,        /* 0=big endian, 1=little endian */
                ! float * scalefactor, /* value to multiply array elements by before truncation to integers */
                ! int * wordsize,      /* number of bytes to use for each array element */
                ! int * status)
           !write(*,"(F10.2)") maxval(modis)
           !stop
           ! do some change
           call match(xdelta, ydelta, lons, lats, modis, newnx, newny, newlons, newlats, newdata, oldindex, newindex)

           ! output the new
           outfile = trim(outpath)//"/"//trim(fname)
           !write(*,*) outfile
           status = system( "mkdir -p "// outpath )
           call countLen(outfile,strlen)
           call write_geogrid(modis, xdelta, ydelta, zdim, signed, endian, scale, wordsize, outfile, strlen)
           !stop
        end do
    end do
    status = system( "cp "//trim(oldpath)//"/index"//" "//outpath )
end program

!=====================================================================================!

character(len=*) function integerToCharacter(ii)
  ! usage:
    ! character(len=5),external :: integerToCharacter
    ! write(*,*) integerToCharacter(2)
    implicit none
    integer, intent(in) :: ii
    write(integerToCharacter,"(I5.5)") ii

end function integerToCharacter

subroutine countLen(string, strlen)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: countLen
!
! Purpose: This routine receives a fortran string, and returns the number of
!    characters in the string before the first "space" is encountered. It
!    considers ascii characters 33 to 126 to be valid characters, and ascii
!    0 to 32, and 127 to be "space" characters.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   implicit none

   ! Arguments
   character (len=*), intent(in) :: string
   integer, intent(out) :: strlen

   ! Local variables
   integer :: i, len_str, aval
   logical :: space

   space = .false.
   i = 1
   len_str = len(string)
   strlen = len_str
   do while ((i .le. len_str) .and. (.not. space))
      aval = ichar(string(i:i))
      if ((aval .lt. 33) .or. (aval .gt. 126)) then
         strlen = i - 1
         space = .true.
      end if
      i = i + 1
   end do
end subroutine countLen

subroutine createlist(start, delta, num, array)
  implicit none
  ! parameter
  real, intent(in)    :: start
  real, intent(in)    :: delta
  integer, intent(in) :: num
  real, intent(out) ,dimension(num) :: array

  ! local
  integer :: i
  do i =1, num
    array(i) = start + (i-1)*delta
  end do

end subroutine createlist

subroutine readnew(fname, nx, ny, lats, lons, array)
  implicit none
  ! calling
  character(len=*),intent(in) :: fname
  integer, intent(in)         :: nx
  integer, intent(in)         :: ny
  real, intent(out) ,dimension(nx, ny) :: lats
  real, intent(out) ,dimension(nx, ny) :: lons
  real, intent(out) ,dimension(nx, ny) :: array

  ! local
  character(200)    :: errormg ! the message of fail opening
  integer           :: ierror  ! Flag for allocate
  open(100,file=fname, form='unformatted', iostat=ierror, iomsg=errormg)
       read(100) lats
       read(100) lons
       read(100) array
  ! Check if writing data is successfule.
  if(ierror /= 0) then
    write(*,*) errormg
    stop
  end if
  close(100)
end subroutine readnew

subroutine match(nx, ny, lons, lats, modis, newnx, newny, newlons, newlats, newdata, oldindex, newindex)

  implicit none
  ! calling
  integer, intent(in)   :: nx
  integer, intent(in)   :: ny
  real, intent(in) ,dimension(nx) :: lons
  real, intent(in) ,dimension(ny) :: lats
  real, intent(inout) ,dimension(nx,ny) :: modis

  integer, intent(in)   :: newnx
  integer, intent(in)   :: newny
  real, intent(in) ,dimension(newnx, newny) :: newlons
  real, intent(in) ,dimension(newnx, newny) :: newlats
  real, intent(in) ,dimension(newnx, newny) :: newdata

  integer, intent(in) :: oldindex
  integer, intent(in) :: newindex

  ! local
  integer :: i,j
  integer :: ix
  integer :: iy

  do i=1,newnx
    do j=1,newny
        if ( newlons(i,j)>lons(1) .and. newlons(i,j)<lons(nx) .and. &
             newlats(i,j)>lats(1) .and. newlats(i,j)<lats(ny) ) then
        if ( int(newdata(i,j)) ==  newindex) then
           call ind(nx, ny, lons, lats, newlons(i,j), newlats(i,j), ix, iy)
           if (ix /=0 .and. iy /=0) then
               modis(ix,iy) = oldindex
           end if
        end if
        end if
    end do
  end do

end subroutine match

subroutine ind(nx, ny, lons, lats, lon, lat, ix, iy)
  implicit none
  ! calling
  integer, intent(in)   :: nx
  integer, intent(in)   :: ny
  real, intent(in), dimension(nx) :: lons
  real, intent(in), dimension(ny) :: lats

  real, intent(in)  ::  lon
  real, intent(in)  ::  lat

  integer, intent(out) :: ix
  integer, intent(out) :: iy

  ! local
   integer :: i

   ix=0
   iy=0
   do i=1,nx-1
      if ( lons(i)<= lon .and. lons(i+1)> lon  ) then
        ix=i
        exit
      end if
   end do

   do i=1,ny-1
      if ( lats(i)<= lat .and. lats(i+1)> lat  ) then
        iy=i
        exit
      end if
   end do
end subroutine  ind

max_stars_repo_name: XiaoLinhong/uplandsue

max_issues_repo_name: XiaoLinhong/uplandsue

max_forks_repo_name: XiaoLinhong/uplandsue


!
!     ------------------------------------------------------------------
!     W R T M A T
!     ------------------------------------------------------------------
!
      SUBROUTINE WRTMAT(A, NROW, NCOL, NMROW, NMCOL)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:49:11  11/20/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: NROW
      INTEGER , INTENT(IN) :: NCOL
      INTEGER , INTENT(IN) :: NMROW
      INTEGER , INTENT(IN) :: NMCOL
      REAL(DOUBLE) , INTENT(IN) :: A(NMROW,NMCOL)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, J
!-----------------------------------------------
!
      DO I = 1, NROW
         WRITE (6, 1010) I, (A(I,J),J=1,NCOL)
!mrg 1010 FORMAT(1H0,I3,2X,4(1X,E14.8),/,(1H ,5X,4(1X,E14.8)))
 1010    FORMAT('0',I5,2X,4(1X,E14.8),/,(' ',7X,4(1X,E14.8)))
      END DO
      RETURN
      END SUBROUTINE WRTMAT


text:
!-----------------------------------------------------------------------
!     Module:        beams3d_follow
!     Authors:       S. Lazerson (lazerson@pppl.gov), M. McMillan (matthew.mcmillan@my.wheaton.edu)
!     Date:          06/20/2012
!     Description:   This subroutine follows the particles through
!                    the grid.  The general ODE which must be solved
!                    can be written:
!                        dX
!                       ----   = fpart(t,q,qdot)
!                        dt
!                    where X=(R,phi,Z).  This allows particle trajectories
!                    to be parametrized in terms of time.
!
!-----------------------------------------------------------------------
SUBROUTINE beams3d_follow
    !-----------------------------------------------------------------------
    !     Libraries
    !-----------------------------------------------------------------------
    USE stel_kinds, ONLY: rprec
    USE beams3d_runtime
    USE beams3d_lines
    USE beams3d_grid, ONLY: tmin, tmax, delta_t, BR_spl, BZ_spl, BPHI_spl, &
                            MODB_spl, S_spl, U_spl, TE_spl, NE_spl, TI_spl, &
                            TE_spl, TI_spl, wall_load, wall_shine, &
                            plasma_mass, plasma_Zavg, plasma_Zmean, therm_factor
    USE mpi_params ! MPI
    USE beams3d_physics_mod
    USE beams3d_write_par
    USE safe_open_mod, ONLY: safe_open
    USE wall_mod, ONLY: wall_free, ihit_array, nface
    USE mpi_inc
    !-----------------------------------------------------------------------
    !     Local Variables
    !          status       MPI stats indicator
    !          mystart      Starting fieldline for each processor
    !          mypace       Number of fieldlines each processor should follow
    !          i            Dummy index
    !          j            Dummy index
    !          sender       Dummy index
    !          ier          Error flag
    !          l            Dummy index
    !          neqs_nag     Number of ODEs to solve (not limited to NAG routines)
    !          l2           Index helper
    !          itol         LSODE tolerance flag
    !          itask        LSODE flag (overshoot and interpolate)
    !          istate       LSODE restart flag
    !-----------------------------------------------------------------------
    IMPLICIT NONE
#if defined(MPI_OPT)
    INTEGER :: status(MPI_STATUS_size) !mpi stuff
    INTEGER :: mystart, mypace, i, j, sender
    INTEGER,ALLOCATABLE :: mnum(:), moffsets(:)
    INTEGER :: MPI_COMM_LOCAL
#endif
    INTEGER :: ier, l, neqs_nag, l2, itol, itask, &
               istate, iopt, lrw, liw, mf, out, iunit
    INTEGER, ALLOCATABLE :: iwork(:), itemp(:,:)
    REAL :: dist
    REAL(rprec) :: tf_max, vel_max, dt_out
    DOUBLE PRECISION, ALLOCATABLE :: w(:), q(:), t_last(:)
    DOUBLE PRECISION :: tf_nag, eps_temp, t_nag, t1_nag, &
                        tol_nag, rtol
    DOUBLE PRECISION :: atol(4), rwork(84)
    DOUBLE PRECISION :: rkh_work(4, 2)
    DOUBLE PRECISION :: qdot1(4)
    CHARACTER*1 :: relab

    DOUBLE PRECISION, PARAMETER :: electron_mass = 9.10938356D-31 !m_e
    DOUBLE PRECISION, PARAMETER :: sqrt_pi       = 1.7724538509   !pi^(1/2)
    DOUBLE PRECISION, PARAMETER :: e_charge      = 1.60217662E-19 !e_c

    !-----------------------------------------------------------------------
    !     External Functions
    !          fpart_nag            RHS of ODE integrator (for NAG)    for BEAMS3D
    !          D02CJF               NAG ODE Solver
    !          D02CJW               NAG Dummy function
    !          jacobian_lsode       Jacobian function (for LSODE, not currently utilized)
    !-----------------------------------------------------------------------
    EXTERNAL D02CJF, D02CJW, fpart_nag, D02CJX, out_beams3d_nag
    EXTERNAL fpart_lsode, jacobian_lsode
    EXTERNAL fpart_rkh68
    !-----------------------------------------------------------------------
    !     Begin Subroutine
    !-----------------------------------------------------------------------
    ! Initializations
    ier = 0
    tol_nag = follow_tol
    neqs_nag = 4
    relab = "M"
    mf = 10

    ! Calc max time to follow particles
    i = MAXLOC(ABS(t_end),1)
    tf_max = t_end(i)

    ! Calculate timestep for integration
    vel_max = MAX(MAXVAL(ABS(vll_start)),1E6)
    dt = SIGN(MAX(lendt_m/vel_max,1D-9),tf_max)

    ! Calculate number of integration timesteps per output timestep
    ndt_max = MAX(CEILING(tf_max/(dt*NPOINC)),1)

    ! Adjust dt to match ndt_max
    nsteps = ndt_max * NPOINC
    dt = tf_max/(ndt_max*NPOINC)
    dt_out = tf_max/NPOINC

    ! Break up the work
    CALL MPI_CALC_MYRANGE(MPI_COMM_BEAMS, 1, nparticles, mystart, myend)

    ! Save mystart and myend
    mystart_save = mystart
    myend_save = myend

    IF (lhitonly) THEN
        npoinc = 2
    END IF

#if defined(MPI_OPT)
      IF (ALLOCATED(mnum)) DEALLOCATE(mnum)
      IF (ALLOCATED(moffsets)) DEALLOCATE(moffsets)
      ALLOCATE(mnum(nprocs_beams), moffsets(nprocs_beams))
      CALL MPI_ALLGATHER((myend-mystart+1)*(npoinc+1),1,MPI_INTEGER,mnum,1,MPI_INTEGER,MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_ALLGATHER((mystart-1)*(npoinc+1),1,MPI_INTEGER,moffsets,1,MPI_INTEGER,MPI_COMM_BEAMS,ierr_mpi)
#endif

    ! Deallocations
    IF (ALLOCATED(q)) DEALLOCATE(q)
    IF (ALLOCATED(w)) DEALLOCATE(w)
    IF (ALLOCATED(R_lines)) DEALLOCATE(R_lines)
    IF (ALLOCATED(Z_lines)) DEALLOCATE(Z_lines)
    IF (ALLOCATED(PHI_lines)) DEALLOCATE(PHI_lines)
    IF (ALLOCATED(vll_lines)) DEALLOCATE(vll_lines)
    IF (ALLOCATED(S_lines)) DEALLOCATE(S_lines)
    IF (ALLOCATED(U_lines)) DEALLOCATE(U_lines)
    IF (ALLOCATED(B_lines)) DEALLOCATE(B_lines)
    IF (ALLOCATED(moment_lines)) DEALLOCATE(moment_lines)
    IF (ALLOCATED(neut_lines)) DEALLOCATE(neut_lines)

    ! Output some stuff
    IF (lverb) THEN
       WRITE(6, '(A)') '----- FOLLOWING PARTICLE TRAJECTORIES -----'
       WRITE(6, '(A,A)')          '      Method: ', TRIM(int_type)
       WRITE(6, '(A,I9)')         '   Particles: ', nparticles
       WRITE(6, '(A,I9,A,EN12.3)') '       Steps: ', nsteps, '   Delta-t: ', dt
       WRITE(6, '(A,I9,A,EN12.3)') '      NPOINC: ', npoinc, '    dt_out: ', dt_out
       SELECT CASE(TRIM(int_type))
          CASE("NAG")
             WRITE(6, '(A,EN12.3,A,A1)') '         Tol: ', follow_tol, '  Type: ', relab
          CASE("LSODE")
             WRITE(6, '(A,EN12.3,A,I2)') '         Tol: ', follow_tol, '  Type: ', mf
       END SELECT
       WRITE(6, '(5X,A,I3,A)', ADVANCE = 'no') 'Trajectory Calculation [', 0, ']%'
       CALL FLUSH(6)
    END IF
    ! Allocations
    ALLOCATE(q(neqs_nag), STAT = ier)
    IF (ier /= 0) CALL handle_err(ALLOC_ERR, 'Q', ier)
    ALLOCATE(R_lines(0:npoinc, mystart:myend), Z_lines(0:npoinc, mystart:myend), &
             PHI_lines(0:npoinc, mystart:myend), vll_lines(0:npoinc, mystart:myend), moment_lines(0:npoinc, mystart:myend), &
             neut_lines(0:npoinc, mystart:myend), S_lines(0:npoinc, mystart:myend), U_lines(0:npoinc, mystart:myend), &
              B_lines(0:npoinc, mystart:myend), STAT = ier)
    IF (ier /= 0) CALL handle_err(ALLOC_ERR, 'R_LINES, PHI_LINES, Z_LINES', ier)
    ALLOCATE(t_last(mystart:myend), STAT = ier)
    IF (ier /= 0) CALL handle_err(ALLOC_ERR, 't_last', ier)

    ! Initializations
    R_lines = 0.0; Z_lines = 0.0; PHI_lines = -1.0
    vll_lines = 0.0; moment_lines = 0.0
    S_lines = 1.5; U_lines = 0.0; B_lines = -1.0
    t_last = 0.0
    R_lines(0, mystart:myend) = R_start(mystart:myend)
    Z_lines(0, mystart:myend) = Z_start(mystart:myend)
    PHI_lines(0, mystart:myend) = phi_start(mystart:myend)
    vll_lines(0, mystart:myend) = vll_start(mystart:myend)
    moment_lines(0, mystart:myend) = mu_start(mystart:myend)
    neut_lines(0, mystart:myend) = .FALSE.
    IF (lbeam) neut_lines(0, mystart:myend) = .TRUE.

    ! Some helpers
    fact_vsound = 1.5*sqrt(e_charge/plasma_mass)*therm_factor
    fact_crit = SQRT(2*e_charge/plasma_mass)*(0.75*sqrt_pi*sqrt(plasma_mass/electron_mass))**(1.0/3.0) ! Wesson pg 226 5.4.9
    fact_kick = pi2*2*SQRT(pi*1E-7*plasma_mass)*E_kick*freq_kick

    ! Follow Trajectories
    IF (mystart <= nparticles) THEN
        SELECT CASE (TRIM(int_type))
            CASE ("NAG")
#if defined(NAG)
                ALLOCATE(w(neqs_nag * 21 + 28), STAT = ier)
                IF (ier /= 0) CALL handle_err(ALLOC_ERR, 'W', ier)
                DO l = mystart, myend
                    ltherm = .false.
                    lneut  = .false.
                    ndt    = ndt_max-1 ! we only record first point
                    q(1) = R_start(l)
                    q(2) = phi_start(l)
                    q(3) = Z_start(l)
                    q(4) = vll_start(l)
                    t_nag = 0.0
                    tf_nag = 0.0
                    mycharge = charge(l)
                    myZ = Zatom(l)
                    mymass = mass(l)
                    moment = mu_start(l)
                    my_end = t_end(l)
                    fact_pa   = plasma_mass*plasma_Zavg/(mymass*plasma_Zmean)
                    fact_coul = myZ*plasma_Zavg*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)
                    myv_neut(:) = v_neut(:,myline)
                    IF (lbeam) lneut = .TRUE.
                    CALL out_beams3d_nag(tf_nag,q)
                    IF (lbeam) THEN
                       lcollision = .FALSE.
                       ! Follow into plasma
                       CALL beams3d_follow_neut(t_nag,q)
                       mytdex = 1; ndt =1
                       tf_nag = t_nag
                       CALL out_beams3d_nag(tf_nag,q)
                       IF (tf_nag > t_end(l)) CYCLE  ! Detect end shinethrough particle
                       ! Ionize
                       CALL beams3d_ionize(tf_nag,q)
                       mytdex = 2; ndt =1
                       CALL out_beams3d_nag(tf_nag,q)
                       ltherm = .FALSE.
                       lcollision = .TRUE.
                       mytdex = 3; ndt =1
                       tf_nag = tf_nag - dt  ! Because out advances t by dt
                       ! Adjust timestep timestep
                       !CALL beams3d_calc_dt(q,moment,mymass,dt)
                    END IF
                    IF (ldepo) CYCLE
                    DO ! Must do it this way becasue lbeam changes q(4) values
                       CALL D02CJF(t_nag,tf_nag,neqs_nag,q,fpart_nag,tol_nag,relab,out_beams3d_nag,D02CJW,w,ier)
                       IF (ier < 0) CALL handle_err(D02CJF_ERR, 'beams3d_follow', ier)
                       t_last(l) = tf_nag ! Save the value here in case out_beams3d changes it
                       CALL out_beams3d_nag(tf_nag,q)
                       IF (ABS(tf_nag) > ABS(my_end)) EXIT
                    END DO
                END DO
#else
                ier = -1
                CALL handle_err(NAG_ERR, 'beams3d_follow', ier)
#endif
            CASE ("RKH68")
                ier = 0
                DO l = mystart, myend
                    myline = l
                    mytdex = 0
                    ndt    = ndt_max-1 ! we only record first point
                    ltherm = .false.
                    lneut  = .false.
                    q(1) = R_start(l)
                    q(2) = phi_start(l)
                    q(3) = Z_start(l)
                    q(4) = vll_start(l)
                    xlast = q(1)*cos(q(2))
                    ylast = q(1)*sin(q(2))
                    zlast = q(3)
                    t_nag = 0.0
                    tf_nag = 0.0
                    mycharge = charge(l)
                    myZ = Zatom(l)
                    mymass = mass(l)
                    mybeam = Beam(l)
                    moment = mu_start(l)
                    my_end = t_end(l)
                    fact_pa   = plasma_mass*plasma_Zavg/(mymass*plasma_Zmean)
                    fact_coul = myZ*plasma_Zavg*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)
                    myv_neut(:) = v_neut(:,myline)
                    IF (lbeam) lneut = .TRUE.
                    CALL out_beams3d_nag(tf_nag,q)
                    IF (lbeam) THEN
                       lcollision = .FALSE.
                       ! Follow into plasma
                       CALL beams3d_follow_neut(t_nag,q)
                       mytdex = 1; ndt =1
                       tf_nag = t_nag
                       CALL out_beams3d_nag(tf_nag,q)
                       IF (tf_nag > t_end(l)) CYCLE  ! Detect end shinethrough particle
                       ! Ionize
                       CALL beams3d_ionize(tf_nag,q)
                       mytdex = 2; ndt =1
                       CALL out_beams3d_nag(tf_nag,q)
                       ltherm = .FALSE.
                       lcollision = .TRUE.
                       mytdex = 3; ndt =1
                       tf_nag = tf_nag - dt  ! Because out advances t by dt
                       ! Adjust timestep timestep
                       !CALL beams3d_calc_dt(q,moment,mymass,dt)
                    END IF
                    IF (ldepo) CYCLE
                    DO
                        CALL drkhvg(t_nag, q, neqs_nag, dt, 2, fpart_rkh68, rkh_work, iopt, ier)
                        IF (ier < 0) CALL handle_err(RKH68_ERR, 'beams3d_follow', ier)
                        q(1)=rkh_work(1,2)
                        q(2)=rkh_work(2,2)
                        q(3)=rkh_work(3,2)
                        q(4)=rkh_work(4,2)
                        t_nag = t_nag+dt
                        tf_nag = tf_nag+dt
                        t_last(l) = tf_nag ! Save the value here in case out_beams3d changes it
                        CALL out_beams3d_nag(tf_nag,q)
                        IF ((istate == -1) .or. (istate ==-2) .or. (ABS(tf_nag) > ABS(my_end)) ) EXIT
                    END DO
                END DO
            CASE ("LSODE","DLSODE")
                ! Open a file
                iunit = myworkid+400
                ier = 0
                IF (ldebug) THEN
                   !CALL safe_open(iunit,ier,'lsode_err_'//TRIM(id_string),'replace','formatted')
                   iunit = 6
                   CALL xsetun(iunit)
                   iunit = 1
                   CALL xsetf(iunit)
                ELSE
                   iunit = 0
                   CALL xsetf(iunit)
                END IF
                ! Allocate helpers
                lrw = 20 + 16 * neqs_nag
                liw = 20
                ALLOCATE(w(lrw))
                ALLOCATE(iwork(liw))
                ier = 0
                DO l = mystart, myend
                    ! Setup LSODE parameters
                    iopt = 0 ! No optional output
                    w = 0; iwork = 0; itask = 1; istate = 1;
                    itol = 2; rtol = follow_tol; atol(:) = follow_tol
                    myline = l
                    mytdex = 0
                    ndt    = ndt_max-1 ! we only record first point
                    ! Initialize the calculation
                    ltherm = .false.
                    lneut  = .false.
                    q(1) = R_start(l)
                    q(2) = phi_start(l)
                    q(3) = Z_start(l)
                    q(4) = vll_start(l)
                    xlast = q(1)*cos(q(2))
                    ylast = q(1)*sin(q(2))
                    zlast = q(3)
                    t_nag = 0.0
                    tf_nag = 0.0
                    mycharge = charge(l)
                    myZ = Zatom(l)
                    mymass = mass(l)
                    mybeam = Beam(l)
                    moment = mu_start(l)
                    my_end = t_end(l)
                    fact_pa   = plasma_mass*plasma_Zavg/(mymass*plasma_Zmean)
                    fact_coul = myZ*plasma_Zavg*(mymass+plasma_mass)/(mymass*plasma_mass*6.02214076208E+26)
                    myv_neut(:) = v_neut(:,myline)
                    ! Setup timestep
                    !CALL beams3d_calc_dt(q,moment,mymass,dt)
                    ! Begin handling particle.
                    IF (lbeam) lneut = .TRUE.
                    CALL out_beams3d_nag(tf_nag,q)
                    IF (lbeam) THEN
                       lcollision = .FALSE.
                       ! Follow into plasma
                       CALL beams3d_follow_neut(t_nag,q)
                       mytdex = 1
                       tf_nag = t_nag
                       CALL out_beams3d_nag(tf_nag,q)
                       IF (tf_nag > t_end(l)) CYCLE  ! Detect end shinethrough particle
                       ! Ionize
                       CALL beams3d_ionize(tf_nag,q)
                       mytdex = 2
                       CALL out_beams3d_nag(tf_nag,q)
                       ltherm = .FALSE.
                       lcollision = .TRUE.
                       mytdex = 3; ndt =1
                       tf_nag = tf_nag - dt  ! Because out advances t by dt
                       ! Adjust timestep timestep
                       !CALL beams3d_calc_dt(q,moment,mymass,dt)
                    END IF
                    IF (ldepo) CYCLE
                    DO
                        IF (lcollision) istate = 1
                        CALL FLUSH(6)
                        CALL DLSODE(fpart_lsode, neqs_nag, q, t_nag, tf_nag, itol, rtol, atol, itask, istate, &
                                   iopt, w, lrw, iwork, liw, jacobian_lsode, mf)
                        IF ((istate == -3) .or. (istate == -4)) THEN
                           ! BIG  DEBUG MESSAGE
                           CALL FLUSH(6)
                           WRITE(6,*) '------------------'
                           WRITE(6,*) '     ',myworkid, l, myline,mf
                           WRITE(6,*) '     ',myworkid, neqs_nag, t_nag, tf_nag
                           WRITE(6,*) '     ',myworkid, q, moment
                           WRITE(6,*) '     ',myworkid, itol, rtol, atol
                           WRITE(6,*) '     ',myworkid, itask,istate,iopt
                           WRITE(6,*) '     ',myworkid, w(5),w(6),w(7)
                           WRITE(6,*) '     ',myworkid, iwork(5),iwork(6),iwork(7)
                           WRITE(6,*) '     ',myworkid, w(11:14)
                           WRITE(6,*) '     ',myworkid, iwork(11:18)
                           WRITE(6,*) '------------------'
                           CALL FLUSH(6)
                           CALL handle_err(LSODE_ERR, 'beams3d_follow', istate)
                        END IF
                        iwork(11) = 0; iwork(12) = 0; iwork(13) = 0
                        t_last(l) = tf_nag ! Save the value here in case out_beams3d changes it
                        CALL out_beams3d_nag(tf_nag,q)
                        IF ((istate == -1) .or. (istate ==-2) .or. (ABS(tf_nag) > ABS(my_end)) ) EXIT
                    END DO
                END DO
                IF (ldebug) CLOSE(iunit)
             CASE ('DEBUG')
                DO l = 0, npoinc
                   R_lines(l,mystart:myend) = REAL(l)
                END DO
                B_lines(0:npoinc,mystart:myend) = REAL(myworkid)
        END SELECT
    END IF

    IF (lverb) THEN
       CALL backspace_out(6,36)
       CALL FLUSH(6)
       WRITE(6,'(36X)',ADVANCE='no')
       CALL FLUSH(6)
       CALL backspace_out(6,36)
       WRITE(6,*)
       CALL FLUSH(6)
    END IF

    ! Check for crash
    istate = 0
    CALL handle_err(MPI_CHECK,'beams3d_follow',istate)

    !Deallocations
    IF (ALLOCATED(q)) DEALLOCATE(q)
    IF (ALLOCATED(w)) DEALLOCATE(w)
    IF (ALLOCATED(iwork)) DEALLOCATE(iwork)

    ! Fix U_lines
    CALL beams3d_fix_poloidal

    ! First reduce the cumulative arrays over shared memory groups then allreduce between shared memeory groups
#if defined(MPI_OPT)
    IF (myid_sharmem == master) THEN
       CALL MPI_REDUCE(MPI_IN_PLACE,   end_state,     nparticles,          MPI_INTEGER, MPI_MAX, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(MPI_IN_PLACE, epower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(MPI_IN_PLACE, ipower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(MPI_IN_PLACE,   ndot_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
    ELSE
       CALL MPI_REDUCE(end_state,     end_state,     nparticles,          MPI_INTEGER, MPI_MAX, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(epower_prof, epower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(ipower_prof, ipower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
       CALL MPI_REDUCE(ndot_prof,     ndot_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_SHARMEM, ierr_mpi)
    END IF

    i = MPI_UNDEFINED
    IF (myid_sharmem == master) i = 0
    CALL MPI_COMM_SPLIT( MPI_COMM_BEAMS,i,myworkid,MPI_COMM_LOCAL,ierr_mpi)
    IF (myid_sharmem == master) THEN
       CALL MPI_ALLREDUCE(MPI_IN_PLACE,   end_state,     nparticles,          MPI_INTEGER, MPI_MAX, MPI_COMM_LOCAL, ierr_mpi)
       CALL MPI_ALLREDUCE(MPI_IN_PLACE, epower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_LOCAL, ierr_mpi)
       CALL MPI_ALLREDUCE(MPI_IN_PLACE, ipower_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_LOCAL, ierr_mpi)
       CALL MPI_ALLREDUCE(MPI_IN_PLACE,   ndot_prof, nbeams*ns_prof1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_LOCAL, ierr_mpi)
       ! This only works becasue of how FORTRAN orders things.
       DO l = 1, ns_prof5
          CALL MPI_ALLREDUCE(MPI_IN_PLACE, dist5d_prof(:,:,:,:,:,l), nbeams*ns_prof1*ns_prof2*ns_prof3*ns_prof4, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_LOCAL, ierr_mpi)
       END DO
       IF (ASSOCIATED(ihit_array)) THEN
          CALL MPI_ALLREDUCE(MPI_IN_PLACE,ihit_array,nface,MPI_INTEGER,MPI_SUM,MPI_COMM_LOCAL,ierr_mpi)
       END IF
       IF (ASSOCIATED(wall_load)) THEN
          CALL MPI_ALLREDUCE(MPI_IN_PLACE,wall_load,nface*nbeams,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LOCAL,ierr_mpi)
       END IF
       IF (ASSOCIATED(wall_shine)) THEN
          CALL MPI_ALLREDUCE(MPI_IN_PLACE,wall_shine,nface*nbeams,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LOCAL,ierr_mpi)
       END IF
       CALL MPI_COMM_FREE(MPI_COMM_LOCAL,ierr_mpi)
    END IF
    CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)

    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'R_lines', DBLVAR=R_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,    'PHI_lines', DBLVAR=PHI_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'Z_lines', DBLVAR=Z_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,    'vll_lines', DBLVAR=vll_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend, 'moment_lines', DBLVAR=moment_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'S_lines', DBLVAR=S_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'U_lines', DBLVAR=U_lines)
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'B_lines', DBLVAR=B_lines)
    CALL beams3d_write1d_parhdf5(         1, nparticles, mystart, myend,      't_end',   DBLVAR=t_last,FILENAME='beams3d_'//TRIM(id_string))
    ALLOCATE(itemp(0:npoinc,mystart:myend))
    itemp = 0; WHERE(neut_lines) itemp=1;
    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,   'neut_lines', INTVAR=itemp)
    DEALLOCATE(itemp)
    IF (ALLOCATED(mnum)) DEALLOCATE(mnum)
    IF (ALLOCATED(moffsets)) DEALLOCATE(moffsets)
    IF (ALLOCATED(t_last)) DEALLOCATE(t_last)
    CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
    IF (ierr_mpi /= 0) CALL handle_err(MPI_BARRIER_ERR, 'beams3d_follow', ierr_mpi)
#endif

    RETURN
    !-----------------------------------------------------------------------
    !     End Subroutine
    !-----------------------------------------------------------------------
END SUBROUTINE beams3d_follow

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   Subroutine that builds G-matrix row for given raypath(s)          c
c   id = 1 for P; id = 2 for S; id = 3 for Qp/Qs from t*(S) and Qp    c
c   id = 4 for bulk attenuation t*(P) and t*(S)                       c
c   Add one point: ipt => piercing point into the Q model             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE getgrowQ(id,maxpts,xarc,yarc,zarc,npts,ipt1,
     &                      grow,ifskip)
      parameter (nxmax = 45, nymax = 45, nzmax = 30)  ! max nodes

      integer maxpts, pt, ipt, ifskip, casep      ! max # internal points in raypath
      integer nnx, nny, nnz, nnxQ, nnyQ, nnzQ, nnzQ1D, nparm, lps
      real x(nxmax), z(nzmax), y(nymax)
      real xq(nxmax), zq(nzmax), yq(nymax), zq1D(nzmax), attin1D(nzmax)
      real*8 sn(2,nzmax,nymax,nxmax)
      real*8 attin(nzmax,nymax,nxmax), Qout(2,nzmax)     ! 1/Qp input for Qp/Qs
      real*8 xarc(0:maxpts+2), yarc(0:maxpts+2), zarc(0:maxpts+2)
      real*8 savexarc(0:maxpts+2), saveyarc(0:maxpts+2)
      real*8 savezarc(0:maxpts+2)
      real*8 xp1,xp2,xp3,xp4,xp5,yp1,yp2,yp3,yp4,yp5,zp1,zp2,zp3,zp4,zp5
      real*8 wc(8), wd(8), slc, wtc, snQ(8), snQ2(8)     ! weighting coefficients
      real*8 xlen, ylen, zlen, sslen, Qx, Qy, Qz
      real grow(nxmax*nymax*nzmax+nzmax) ! G-matrix row for given raypaths
      real sslength, xlength, ylength, zlength, alllen

      common /a/ x, y, z, xq, yq, zq
      common /b/ ipolar,nnx,nny,nnz,re,pi,r2d,nnxQ,nnyQ,nnzQ,nnzQ1D
      common /c/ sn
      common /d/ attin, zq1D, attin1D

      ifskip = 0    ! If skip this event because most part of ray is outside
c     check if raypath is in the velocity and Q model space
      do pt = 0, npts+1
        if ((xarc(pt) .lt. x(1)).or.(xarc(pt) .gt. x(nnx)).or.
     &      (zarc(pt) .lt. z(1)).or.(zarc(pt) .gt. z(nnz)).or.
     &      (yarc(pt) .lt. y(1)).or.(yarc(pt) .gt. y(nny))) then
          write (*,*) "ERROR: ray is not in the velocity model space"
          write (*,*) pt, xarc(pt), yarc(pt), zarc(pt)
c          stop
          ifskip = 1
          goto 237
        endif
      enddo

c     Find the first point of the ray within the Q model
      ipt=0
      do i = 0, npts+2
        if ((xarc(i) .le. xq(nnxQ)).and.(xarc(i) .ge. xq(1)).and.
     &      (yarc(i) .le. yq(nnyQ)).and.(yarc(i) .ge. yq(1)).and.
     &      (zarc(i) .le. zq(nnzQ)).and.(zarc(i) .ge. zq(1))) then
          ipt = i
          goto 236
        endif
        ipt = npts+2
      enddo
236   continue
      if (ipt.eq.npts+2) then
        ifskip = 1
        goto 237
      endif
c      write(*,'(6f10.4)') xq(1),xq(nnxQ),yq(1),yq(nnyQ),zq(1),zq(nnzQ)
c      write(*,*) ipt,npts
c      write(*,300) "Earthquake:  ",xarc(0),yarc(0),zarc(0)
c      write(*,300) "First point: ",xarc(ipt),yarc(ipt),zarc(ipt)
c      write(*,300) "Point before:",xarc(ipt-1),yarc(ipt-1),zarc(ipt-1)
      do i = 0, npts+2
        savexarc(i) = xarc(i)
        saveyarc(i) = yarc(i)
        savezarc(i) = zarc(i)
      enddo

ccccc Initiating the row of G matrix
      nparm = nnxQ*nnyQ*nnzQ+nnzQ1D
c      write(*,*) "nparm = ",nparm
      do i = 1, nparm
        grow(i) = 0.
      enddo

ccccc If the whole ray is within the Q model space, jump to buiding G
      if (ipt.eq.0) then
c        ipt = 1
        goto 235
      endif

c      alllen=0.
cccc  -- build G-matrix row (grow) for arc segment in 1D Q model-- cccc
c     find and insert the piercing point
      casep=0
      if (zarc(ipt-1).gt.zq(nnzQ)) then   ! Piercing from bottom
        zp1 = zq(nnzQ)
        xp1 = xarc(ipt)-(zarc(ipt)-zq(nnzQ))/(zarc(ipt)-zarc(ipt-1))
     &                  *(xarc(ipt)-xarc(ipt-1))
        yp1 = yarc(ipt)-(zarc(ipt)-zq(nnzQ))/(zarc(ipt)-zarc(ipt-1))
     &                  *(yarc(ipt)-yarc(ipt-1))
        if ((xp1 .ge. xq(1)).and.(xp1 .le. xq(nnxQ)).and.
     &            (yp1 .ge. yq(1)).and.(yp1 .le. yq(nnyQ))) casep=1
      endif
      if (xarc(ipt-1).lt.xq(1)) then   ! Piercing from left
        xp2 = zq(1)
        zp2 = zarc(ipt)-(xarc(ipt)-xq(1))/(xarc(ipt)-xarc(ipt-1))
     &                  *(zarc(ipt)-zarc(ipt-1))
        yp2 = yarc(ipt)-(xarc(ipt)-xq(1))/(xarc(ipt)-xarc(ipt-1))
     &                  *(yarc(ipt)-yarc(ipt-1))
        if ((zp2 .ge. zq(1)).and.(zp2 .le. zq(nnzQ)).and.
     &            (yp2 .ge. yq(1)).and.(yp2 .le. yq(nnyQ))) casep=2
      endif
      if (xarc(ipt-1).gt.xq(nnxQ)) then   ! Piercing from right
        xp3 = xq(nnxQ)
        zp3 = zarc(ipt)-(xarc(ipt)-xq(nnxQ))/(xarc(ipt)-xarc(ipt-1))
     &                  *(zarc(ipt)-zarc(ipt-1))
        yp3 = yarc(ipt)-(xarc(ipt)-xq(nnxQ))/(xarc(ipt)-xarc(ipt-1))
     &                  *(yarc(ipt)-yarc(ipt-1))
        if ((zp3 .ge. zq(1)).and.(zp3 .le. zq(nnzQ)).and.
     &            (yp3 .ge. yq(1)).and.(yp3 .le. yq(nnyQ))) casep=3
      endif
      if (yarc(ipt-1).lt.yq(1)) then   ! Piercing from front
        yp4 = yq(1)
        zp4 = zarc(ipt)-(yarc(ipt)-yq(1))/(yarc(ipt)-yarc(ipt-1))
     &                  *(zarc(ipt)-zarc(ipt-1))
        xp4 = xarc(ipt)-(yarc(ipt)-yq(1))/(yarc(ipt)-yarc(ipt-1))
     &                  *(xarc(ipt)-xarc(ipt-1))
        if ((zp4 .ge. zq(1)).and.(zp4 .le. zq(nnzQ)).and.
     &            (xp4 .ge. xq(1)).and.(xp4 .le. xq(nnxQ))) casep=4
      endif
      if (yarc(ipt-1).gt.yq(nnyQ)) then   ! Piercing from back
        yp5 = yq(nnyQ)
        zp5 = zarc(ipt)-(yarc(ipt)-yq(nnyQ))/(yarc(ipt)-yarc(ipt-1))
     &                  *(zarc(ipt)-zarc(ipt-1))
        xp5 = xarc(ipt)-(yarc(ipt)-yq(nnyQ))/(yarc(ipt)-yarc(ipt-1))
     &                  *(xarc(ipt)-xarc(ipt-1))
        if ((zp5 .ge. zq(1)).and.(zp5 .le. zq(nnzQ)).and.
     &            (xp5 .ge. xq(1)).and.(xp5 .le. xq(nnxQ))) casep=5
      endif
c      write(*,*) "casep = ",casep
      npts = npts + 1
      ipt1 = ipt1 + 1
      do i = ipt+1, npts+2
        xarc(i) = savexarc(i-1)
        yarc(i) = saveyarc(i-1)
        zarc(i) = savezarc(i-1)
      enddo
      if (casep.eq.0) then
        write(*,*) "Warning: fail to find the piercing point"
        write(*,200) xarc(ipt),yarc(ipt),zarc(ipt)
        write(*,200) xarc(ipt-1),yarc(ipt-1),zarc(ipt-1)
        ifskip = 1
        goto 237
      elseif (casep.eq.1) then
        xarc(ipt) = xp1
        yarc(ipt) = yp1
        zarc(ipt) = zp1
      elseif (casep.eq.2) then
        xarc(ipt) = xp2
        yarc(ipt) = yp2
        zarc(ipt) = zp2
      elseif (casep.eq.3) then
        xarc(ipt) = xp3
        yarc(ipt) = yp3
        zarc(ipt) = zp3
      elseif (casep.eq.4) then
        xarc(ipt) = xp4
        yarc(ipt) = yp4
        zarc(ipt) = zp4
      elseif (casep.eq.5) then
        xarc(ipt) = xp5
        yarc(ipt) = yp5
        zarc(ipt) = zp5
      endif
c      write(*,*) xarc(ipt),yarc(ipt),zarc(ipt),casep,ipt
c      write(*,*) "One ray"

      if (zarc(ipt).lt.zq(nnzQ)) then
c        write(*,*) "Skip outside event piercing too shallow"
c        write(*,*) zarc(ipt-1), zarc(ipt), zarc(ipt+1)
        ifskip = 1
        goto 237
      endif
      x0 = xarc(0)
      y0 = yarc(0)
      z0 = zarc(0)
      x1 = xarc(1)
      y1 = yarc(1)
      z1 = zarc(1)
      pt = 1    ! End of the segment
c      write (*,*) (zarc(i),i=0,ipt)
      do while (pt .le. ipt)
        sign = +1.            ! ray moving upwards
        if (z1 .gt. z0) sign = -1.
        if (x1 .eq. x0) then
          xslope = 100000.
          yxslope = 100000.
        else
          xslope = (z1-z0)/(x1-x0)    ! positive means up & left or down & right
          yxslope = (y1-y0)/(x1-x0)
        endif
        if (y1 .eq. y0) then
          yslope = 100000.
          xyslope = 100000.
        else
          yslope = (z1-z0)/(y1-y0)
          xyslope = (x1-x0)/(y1-y0)
        endif
        b = z1 - xslope*x1
        c = z1 - yslope*y1
        do i = 1, nnx   ! i0: node nearest 0 in x direction
          if (x(i).lt.x0 .and. x(i+1).ge.x0 .and. xslope.ge.0.) i0 = i
          if (x(i).le.x0 .and. x(i+1).gt.x0 .and. xslope.lt.0.) i0 = i
        enddo
        if (x0 .eq. x(i0)   .and. sign*xslope .gt. 0.) i0 = i0-1
        if (x0 .eq. x(i0+1) .and. sign*xslope .lt. 0.) i0 = i0+1
        if (xslope*sign .ge. 0.) xi = x(i0)
        if (xslope*sign .lt. 0.) xi = x(i0+1)
        if (x1 .eq. x0) then
          lenz0=10000
          do ii = 1,nnz
            if (((z(ii)-z0)*(z(ii)-z1) .le. 0) .and.
     &            ((abs(z(ii)-z0)) .ge. (abs(z(ii)-z1))) .and.
     &            (abs(z(ii)-z0) .lt. lenz0)) then
              lenz0=abs(z(ii)-z0)
              zi = z(ii) ! between z1 and z0 but close to z0
            endif
          enddo
          leny0=10000
          do ii = 1,nny
            if (((y(ii)-y0)*(y(ii)-y1) .le. 0) .and.
     &            ((abs(y(ii)-y0)) .ge. (abs(y(ii)-y1))) .and.
     &            (abs(y(ii)-y0) .lt. leny0)) then
              leny0=abs(y(ii)-y0)
              yi = y(ii) ! between y1 and y0 but close to y0
            endif
          enddo
        else
          zi = xslope*xi + b
          yi = y0 + yxslope*(xi-x0)
        endif
        do j = 1, nny   ! j0: node nearest 0 in y direction
          if (y(j).lt.y0 .and. y(j+1).ge.y0 .and. yslope.ge.0.) j0 = j
          if (y(j).le.y0 .and. y(j+1).gt.y0 .and. yslope.lt.0.) j0 = j
        enddo
        if (y0 .eq. y(j0)   .and. sign*yslope .gt. 0.) j0 = j0-1
        if (y0 .eq. y(j0+1) .and. sign*yslope .lt. 0.) j0 = j0+1
        if (yslope*sign .ge. 0.) yj = y(j0)
        if (yslope*sign .lt. 0.) yj = y(j0+1)
        if (y1 .eq. y0) then
          lenz0=10000
          do ii = 1,nnz
            if (((z(ii)-z0)*(z(ii)-z1) .le. 0) .and.
     &            ((abs(z(ii)-z0)) .ge. (abs(z(ii)-z1))) .and.
     &            (abs(z(ii)-z0) .lt. lenz0)) then
              lenz0=abs(z(ii)-z0)
              zj = z(ii) ! between z1 and z0 but close to z0
            endif
          enddo
          lenx0=10000
          do ii = 1,nny
            if (((x(ii)-x0)*(x(ii)-x1) .le. 0) .and.
     &            ((abs(x(ii)-x0)) .ge. (abs(x(ii)-x1))) .and.
     &            (abs(x(ii)-x0) .lt. lenx0)) then
              leny0=abs(x(ii)-x0)
              xj = x(ii) ! between x1 and x0 but close to x0
            endif
          enddo
        else
          zj = yslope*yj + c
          xj = x0 + xyslope*(yj-y0)
        endif
        do k = 1, nnz   ! k0: node nearest 0 in z direction
          if (z(k) .lt. z0 .and. z(k+1) .ge. z0) k0 = k
        enddo
c        do k = 1, nnzQ
c          if (zq(k) .lt. z0 .and. zq(k+1) .ge. z0) kq0 = k
c        enddo
        if (z0 .eq. z(k0+1) .and. sign .lt. 0.) k0 = k0+1
c     Determine k1D from 1 to nnzQ1D
c        if (k0 .le. nnzQ) then
c          k1D = nnzQ+1
c        else
c          k1D = k0
c        endif
        k1D = 1
        do k = 1, nnzQ1D
          if (zq1D(k) .lt. z0 .and. zq1D(k+1) .ge. z0) k1D = k
        enddo
        if (z0 .eq. zq1D(k1D+1) .and. sign .lt. 0.) k1D = k1D + 1
c        print *,z0,zq1D(k1D),zq1D(k1D+1),k1D
c        if (z0 .eq. zq(kq0+1) .and. sign .lt. 0.) kq0 = kq0+1
        if (sign .gt. 0.) zk = z(k0)
        if (sign .lt. 0.) zk = z(k0+1)
        if (z1 .eq. z0) then
          lenx0=10000
          do ii = 1,nny
            if (((x(ii)-x0)*(x(ii)-x1) .le. 0) .and.
     &            ((abs(x(ii)-x0)) .ge. (abs(x(ii)-x1))) .and.
     &            (abs(x(ii)-x0) .lt. lenx0)) then
              leny0=abs(x(ii)-x0)
              xk = x(ii) ! between x1 and x0 but close to x0
            endif
          enddo
          leny0=10000
          do ii = 1,nny
            if (((y(ii)-y0)*(y(ii)-y1) .le. 0) .and.
     &            ((abs(y(ii)-y0)) .ge. (abs(y(ii)-y1))) .and.
     &            (abs(y(ii)-y0) .lt. leny0)) then
              leny0=abs(y(ii)-y0)
              yk = y(ii) ! between y1 and y0 but close to y0
            endif
          enddo
        else
          xk = (zk - b)/xslope
          yk = y0 + yxslope*(xk-x0)
        endif
        case1length = (x1-x0)**2+(y1-y0)**2+(z1-z0)**2 ! dist between 1 and 0
        xlength = (xi-x0)**2+(yi-y0)**2+(zi-z0)**2  ! dist between i0 and 0
        ylength = (xj-x0)**2+(yj-y0)**2+(zj-z0)**2  ! dist between j0 and 0
        zlength = (xk-x0)**2+(yk-y0)**2+(zk-z0)**2  ! dist between k0 and 0
        icase = 1       ! subsegment ends at arc node
        sslength = case1length
        if (xlength .le. sslength) then   ! subsegment ends at x(i)
          icase = 2
          sslength = xlength
        endif
        if (ylength .le. sslength) then    ! subsegment ends at y(j)
          icase = 3
          sslength = ylength
        endif
        if (zlength .le. sslength) then    ! subsegment ends at z(k)
          icase = 4
        endif
c        write(*,*) icase, sslength, xlength, ylength, zlength
        if (icase .eq. 1) then            ! set subsegment end location
          xss = x1
          yss = y1
          zss = z1
        else if (icase .eq. 2) then
          xss = xi
          yss = yi
          zss = zi
        else if (icase .eq. 3) then
          xss = xj
          yss = yj
          zss = zj
        else if (icase .eq. 4) then
          xss = xk
          yss = yk
          zss = zk
        else
          write (*,*) "icase out of bounds, (", icase, ")"
          write (*,*) 'x1,y1,z1,xi,yi,zi,xj,yj,zj,xk,yk,zk'
          write (*,*) x1,y1,z1,xi,yi,zi,xj,yj,zj,xk,yk,zk
          stop
        endif
        xlen = xss-x0
        ylen = yss-y0
        zlen = zss-z0
        if (ipolar .eq. 1) then     ! x,y distances decr w/ depth
          xlen = xlen*(1.-(zss+z0)/(2.*re))
          ylen = ylen*(1.-(zss+z0)/(2.*re))
        endif
        sslen = sqrt((xlen**2) + (ylen**2) + (zlen**2))
c        write(*,*) "Length of the subsegment is ",sslen,icase
c        alllen=alllen+sslen
c        write(*,*) x0,y0,z0,sn(id,k0,j0,i0)
ccccc Adding grow elements corresponding to 1D Q model from 1 to nnzQ1D
        if ((id .eq. 1) .or. (id .eq. 2)) then          ! Qp or Qs
          grow(nnxQ*nnyQ*nnzQ+k1D) = grow(nnxQ*nnyQ*nnzQ+k1D)+
     &                                    sslen*sn(id,k0,j0,i0)
        elseif (id .eq. 3) then    ! Qp/Qs
          grow(nnxQ*nnyQ*nnzQ+k1D) = grow(nnxQ*nnyQ*nnzQ+k1D)+
     &                                sslen*sn(2,k0,j0,i0)*attin1D(k1D)
C       elseif (id .eq. 4) then    ! Qk
C         grow(nnxQ*nnyQ*nnzQ+k1D) = grow(nnxQ*nnyQ*nnzQ+k1D)+
C    &                  1-3./4.*(sn(2,k0,j0,i0)**2)/(sn(1,k0,j0,i0)**2)
        endif
        x0 = xss
        y0 = yss
        z0 = zss
        if (z0 .eq. zarc(ipt)) then
          pt = pt+1
        elseif (icase .eq. 1) then      ! Move to next subsegment for case 1
          pt = pt + 1
          x1 = xarc(pt)
          y1 = yarc(pt)
          z1 = zarc(pt)
        endif
c        write(*,*) nnxQ*nnyQ*nnzQ+k0-nnzQ, grow(nnxQ*nnyQ*nnzQ+k0-nnzQ)
      enddo
c      write(*,*) "t* outside = ", tstarout
c      write(*,*) grow(nnxQ*nnyQ*nnzQ+1)
cccc  -- build G-matrix row (grow) for arc segment in 1D Q model-- cccc

cccc  -- build G-matrix row (grow) for arc segment in 3D Q model-- cccc
235   continue
c      write(*,*) 'npts = ', npts
      x0 = xarc(ipt)
      y0 = yarc(ipt)
      z0 = zarc(ipt)
      x1 = xarc(ipt+1)
      y1 = yarc(ipt+1)
      z1 = zarc(ipt+1)
c      write(*,*) x0,y0,z0
c      write(*,*) x1,y1,z1
c      write(*,'(6f10.4)') xq(1),xq(nnxQ),yq(1),yq(nnyQ),zq(1),zq(nnzQ)
c      write(*,*) "Adding G-matrix", ipt1,ipt
cccc  Start segments from bottom
      pt = ipt
      do while (pt .le. ipt1)
c        write(*,*) pt, x1, x0, y1, y0, z1, z0
        sign = +1.            ! ray moving upwards
        if (z1 .gt. z0) sign = -1.
        if (x1 .eq. x0) then
          xslope = 100000.
          yxslope = 100000.
        else
          xslope = (z1-z0)/(x1-x0)    ! positive means up & left or down & right
          yxslope = (y1-y0)/(x1-x0)
        endif
        if (y1 .eq. y0) then
          yslope = 100000.
          xyslope = 100000.
        else
          yslope = (z1-z0)/(y1-y0)
          xyslope = (x1-x0)/(y1-y0)
        endif
c        write(*,*) "yxslope, xslope, xyslope, yslope",
c     &            yxslope, xslope, xyslope, yslope
         ! sign*xslope => ray moving right-left (+left, -right)
         ! sign*yslope => ray moving front-back (+front, -back)

        b = z1 - xslope*x1
        c = z1 - yslope*y1

c        do i = 1, nnx   ! i0: node nearest 0 in x direction
c          if (x(i).lt.x0 .and. x(i+1).ge.x0 .and. xslope.ge.0.) i0 = i
c          if (x(i).le.x0 .and. x(i+1).gt.x0 .and. xslope.lt.0.) i0 = i
c        enddo
        iq0=0
        do i = 1, nnxQ
          if (xq(i).lt.x0 .and. xq(i+1).ge.x0 .and. xslope.ge.0.) iq0=i
          if (xq(i).le.x0 .and. xq(i+1).gt.x0 .and. xslope.lt.0.) iq0=i
        enddo
c        if (x0 .eq. x(i0)   .and. sign*xslope .gt. 0.) i0 = i0-1
c        if (x0 .eq. x(i0+1) .and. sign*xslope .lt. 0.) i0 = i0+1
        if (x0 .eq. xq(iq0)   .and. sign*xslope .gt. 0.) iq0 = iq0-1
        if (x0 .eq. xq(iq0+1) .and. sign*xslope .lt. 0.) iq0 = iq0+1
        if (xslope*sign .ge. 0.) xi = xq(iq0)
        if (xslope*sign .lt. 0.) xi = xq(iq0+1)
        if (x1 .eq. x0) then
          lenz0=10000
          do ii = 1,nnzQ
            if (((zq(ii)-z0)*(zq(ii)-z1) .le. 0) .and.
     &            ((abs(zq(ii)-z0)) .ge. (abs(zq(ii)-z1))) .and.
     &            (abs(zq(ii)-z0) .lt. lenz0)) then
              lenz0=abs(zq(ii)-z0)
              zi = zq(ii) ! between z1 and z0 but close to z0
            endif
          enddo
          leny0=10000
          do ii = 1,nnyQ
            if (((yq(ii)-y0)*(yq(ii)-y1) .le. 0) .and.
     &            ((abs(yq(ii)-y0)) .ge. (abs(yq(ii)-y1))) .and.
     &            (abs(yq(ii)-y0) .lt. leny0)) then
              leny0=abs(yq(ii)-y0)
              yi = yq(ii) ! between y1 and y0 but close to y0
            endif
          enddo
        else
          zi = xslope*xi + b
          yi = y0 + yxslope*(xi-x0)
        endif

c        do j = 1, nny   ! j0: node nearest 0 in y direction
c          if (y(j).lt.y0 .and. y(j+1).ge.y0 .and. yslope.ge.0.) j0 = j
c          if (y(j).le.y0 .and. y(j+1).gt.y0 .and. yslope.lt.0.) j0 = j
c        enddo
        jq0=0
        do j = 1, nnyQ
          if (yq(j).lt.y0 .and. yq(j+1).ge.y0 .and. yslope.ge.0.) jq0=j
          if (yq(j).le.y0 .and. yq(j+1).gt.y0 .and. yslope.lt.0.) jq0=j
        enddo
c        if (y0 .eq. y(j0)   .and. sign*yslope .gt. 0.) j0 = j0-1
c        if (y0 .eq. y(j0+1) .and. sign*yslope .lt. 0.) j0 = j0+1
        if (y0 .eq. yq(jq0)   .and. sign*yslope .gt. 0.) jq0 = jq0-1
        if (y0 .eq. yq(jq0+1) .and. sign*yslope .lt. 0.) jq0 = jq0+1
        if (yslope*sign .ge. 0.) yj = yq(jq0)
        if (yslope*sign .lt. 0.) yj = yq(jq0+1)
        if (y1 .eq. y0) then
          lenz0=10000
          do ii = 1,nnzQ
            if (((zq(ii)-z0)*(zq(ii)-z1) .le. 0) .and.
     &            ((abs(zq(ii)-z0)) .ge. (abs(zq(ii)-z1))) .and.
     &            (abs(zq(ii)-z0) .lt. lenz0)) then
              lenz0=abs(zq(ii)-z0)
              zj = zq(ii) ! between z1 and z0 but close to z0
            endif
          enddo
          lenx0=10000
          do ii = 1,nnyQ
            if (((xq(ii)-x0)*(xq(ii)-x1) .le. 0) .and.
     &            ((abs(xq(ii)-x0)) .ge. (abs(xq(ii)-x1))) .and.
     &            (abs(xq(ii)-x0) .lt. lenx0)) then
              leny0=abs(xq(ii)-x0)
              xj = xq(ii) ! between x1 and x0 but close to x0
            endif
          enddo
        else
          zj = yslope*yj + c
          xj = x0 + xyslope*(yj-y0)
        endif

c        do k = 1, nnz   ! k0: node nearest 0 in z direction
c          if (z(k) .lt. z0 .and. z(k+1) .ge. z0) k0 = k
c        enddo
        kq0 = 0
        do k = 1, nnzQ
          if (zq(k) .lt. z0 .and. zq(k+1) .ge. z0) kq0 = k
        enddo
c        if (z0 .eq. z(k0+1) .and. sign .lt. 0.) k0 = k0+1
        if (z0 .eq. zq(kq0+1) .and. sign .lt. 0.) kq0 = kq0+1
        if (kq0 .gt. 0) then
          if (sign .gt. 0.) zk = zq(kq0)
          if (sign .lt. 0.) zk = zq(kq0+1)
        else
          zk = z1
        endif
        if (z1 .eq. z0) then
          lenx0=10000
          do ii = 1,nnyQ
            if (((xq(ii)-x0)*(xq(ii)-x1) .le. 0) .and.
     &            ((abs(xq(ii)-x0)) .ge. (abs(xq(ii)-x1))) .and.
     &            (abs(xq(ii)-x0) .lt. lenx0)) then
              leny0=abs(xq(ii)-x0)
              xk = xq(ii) ! between x1 and x0 but close to x0
            endif
          enddo
          leny0=10000
          do ii = 1,nnyQ
            if (((yq(ii)-y0)*(yq(ii)-y1) .le. 0) .and.
     &            ((abs(yq(ii)-y0)) .ge. (abs(yq(ii)-y1))) .and.
     &            (abs(yq(ii)-y0) .lt. leny0)) then
              leny0=abs(yq(ii)-y0)
              yk = yq(ii) ! between y1 and y0 but close to y0
            endif
          enddo
        else
          xk = (zk - b)/xslope
          yk = y0 + yxslope*(xk-x0)
        endif

        case1length = (x1-x0)**2+(y1-y0)**2+(z1-z0)**2 ! dist between 1 and 0
        xlength = (xi-x0)**2+(yi-y0)**2+(zi-z0)**2  ! dist between i0 and 0
        ylength = (xj-x0)**2+(yj-y0)**2+(zj-z0)**2  ! dist between j0 and 0
        zlength = (xk-x0)**2+(yk-y0)**2+(zk-z0)**2  ! dist between k0 and 0

        icase = 1       ! subsegment ends at arc node
        sslength = case1length
        if (xlength .lt. sslength) then   ! subsegment ends at x(i)
          icase = 2
          sslength = xlength
        endif
        if (ylength .lt. sslength) then    ! subsegment ends at y(j)
          icase = 3
          sslength = ylength
        endif
        if (z0 .le. zq(1)) then
          sslength = 0.
        endif
        if (zlength .lt. sslength) then    ! subsegment ends at z(k)
          icase = 4
        endif
c        write(*,*) icase, sslength, xlength, ylength, zlength

        if (icase .eq. 1) then            ! set subsegment end location
          xss = x1
          yss = y1
          zss = z1
        else if (icase .eq. 2) then
          xss = xi
          yss = yi
          zss = zi
        else if (icase .eq. 3) then
          xss = xj
          yss = yj
          zss = zj
        else if (icase .eq. 4) then
          xss = xk
          yss = yk
          zss = zk
        else
          write (*,*) "icase out of bounds, (", icase, ")"
          write (*,*) 'x1,y1,z1,xi,yi,zi,xj,yj,zj,xk,yk,zk'
          write (*,*) x1,y1,z1,xi,yi,zi,xj,yj,zj,xk,yk,zk
          stop
        endif
c        write(*,*) x1,y1,z1,xss,yss,zss
c        write(*,*) icase, sslength, xss, yss, zss
c        write(*,*) "start"
c        write(*,*) x0,y0,z0,xlength,ylength
c        write(*,*) xss,yss,zss,case1length,zlength

c     Find slowness at each Q node with trilinear interpolation
cccc          010--------110
cccc          /|         /|
cccc        /  |       /  |
cccc     000--------100   |
cccc      |    |     |    |
cccc      |   011--- | --111
cccc      |   /      |   /
cccc      | /        | /
cccc     001--------101
        iw = 0
        if ((id .eq. 1) .or. (id .eq. 2)) then          ! Qp or Qs
          lps = id
        elseif (id .eq. 3) then    ! Qp/Qs
          lps = 2
        elseif (id .eq. 4) then    ! Qk
          lps = 1
        endif
        do jy = jq0, jq0+1
          do kz = kq0, kq0+1
            do ix = iq0, iq0+1
              iw = iw + 1
              Qx=xq(ix)
              Qy=yq(jy)
              do i = 1, nnx
                if (x(i).lt.Qx .and. x(i+1).ge.Qx) iv = i
              enddo
              do j = 1, nny
                if (y(j).lt.Qy .and. y(j+1).ge.Qy) jv = j
              enddo
              if ((ix.eq.0).or.(jy.eq.0).or.(kz.eq.0)) then
                snQ(iw)=0.
              else
                Qz=zq(kz)
                do k = 1, nnz
                  if (z(k).lt.Qz .and. z(k+1).ge.Qz) kv = k
                enddo
                sn00=sn(lps,kv,jv,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                    (sn(lps,kv,jv,iv+1)-sn(lps,kv,jv,iv))
                sn01=sn(lps,kv+1,jv,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                    (sn(lps,kv+1,jv,iv+1)-sn(lps,kv+1,jv,iv))
                sn10=sn(lps,kv,jv+1,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                    (sn(lps,kv,jv+1,iv+1)-sn(lps,kv,jv+1,iv))
                sn11=sn(lps,kv+1,jv+1,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                    (sn(lps,kv+1,jv+1,iv+1)-sn(lps,kv+1,jv+1,iv))
                sn0=sn00+(Qy-y(jv))/(y(jv+1)-y(jv))*(sn10-sn00)
                sn1=sn01+(Qy-y(jv))/(y(jv+1)-y(jv))*(sn11-sn01)
                snQ(iw)=sn0+(Qz-z(kv))/(z(kv+1)-z(kv))*(sn1-sn0)
              endif
            enddo
          enddo
        enddo
        sn000 = snQ(1)
        sn100 = snQ(2)
        sn001 = snQ(3)
        sn101 = snQ(4)
        sn010 = snQ(5)
        sn110 = snQ(6)
        sn011 = snQ(7)
        sn111 = snQ(8)

        if (id .eq. 4) then    ! Qk
          do jy = jq0, jq0+1
            do kz = kq0, kq0+1
              do ix = iq0, iq0+1
                iw = iw + 1
                Qx=xq(ix)
                Qy=yq(jy)
                do i = 1, nnx
                  if (x(i).lt.Qx .and. x(i+1).ge.Qx) iv = i
                enddo
                do j = 1, nny
                  if (y(j).lt.Qy .and. y(j+1).ge.Qy) jv = j
                enddo
                if ((ix.eq.0).or.(jy.eq.0).or.(kz.eq.0)) then
                  snQ2(iw)=0.
                else
                  Qz=zq(kz)
                  do k = 1, nnz
                    if (z(k).lt.Qz .and. z(k+1).ge.Qz) kv = k
                  enddo
                  sn00=sn(2,kv,jv,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                (sn(2,kv,jv,iv+1)-sn(2,kv,jv,iv))
                  sn01=sn(2,kv+1,jv,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                (sn(2,kv+1,jv,iv+1)-sn(2,kv+1,jv,iv))
                  sn10=sn(2,kv,jv+1,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                (sn(2,kv,jv+1,iv+1)-sn(2,kv,jv+1,iv))
                  sn11=sn(2,kv+1,jv+1,iv)+(Qx-x(iv))/(x(iv+1)-x(iv))*
     &                (sn(2,kv+1,jv+1,iv+1)-sn(2,kv+1,jv+1,iv))
                  sn0=sn00+(Qy-y(jv))/(y(jv+1)-y(jv))*(sn10-sn00)
                  sn1=sn01+(Qy-y(jv))/(y(jv+1)-y(jv))*(sn11-sn01)
                  snQ2(iw)=sn0+(Qz-z(kv))/(z(kv+1)-z(kv))*(sn1-sn0)
                endif
              enddo
            enddo
          enddo
          sn000 = sn000/snQ2(1)
          sn100 = sn100/snQ2(2)
          sn001 = sn001/snQ2(3)
          sn101 = sn101/snQ2(4)
          sn010 = sn010/snQ2(5)
          sn110 = sn110/snQ2(6)
          sn011 = sn011/snQ2(7)
          sn111 = sn111/snQ2(8)
        endif
c        write(*,*) sn000,iq0,jq0,kq0

c     Trapezoidal weighting: c - to begin of segment, d - to end of segment
        iw = 0
        slc = 0.
        sld = 0.
        do jy = j0, j0+1
          do kz = k0, k0+1
            do ix = i0, i0+1
              iw = iw + 1
              if ((ix.eq.0).or.(jy.eq.0).or.(kz.eq.0)) then
                wc(iw)=0.
                wd(iw)=0.
              else
                wc(iw) = sqrt(((x(ix)-x0)*(1.-z(kz)/re))**2 +
     &            ((y(jy)-y0)*(1.-z(kz)/re))**2 + (z(kz)-z0)**2)
                wd(iw) =  sqrt(((x(ix)-xss)*(1.-z(kz)/re))**2 +
     &            ((y(jy)-yss)*(1.-z(kz)/re))**2 + (z(kz)-zss)**2)
              endif
              slc = slc + wc(iw)                 ! sum of weighting lengths
              sld = sld + wd(iw)                 ! sum of weighting lengths
            enddo
          enddo
        enddo
        c000 = ((wc(2)+wc(3)+wc(4)+wc(5)+wc(6)+wc(7)+wc(8))/7.)/slc
        c100 = ((wc(1)+wc(3)+wc(4)+wc(5)+wc(6)+wc(7)+wc(8))/7.)/slc
        c001 = ((wc(1)+wc(2)+wc(4)+wc(5)+wc(6)+wc(7)+wc(8))/7.)/slc
        c101 = ((wc(1)+wc(2)+wc(3)+wc(5)+wc(6)+wc(7)+wc(8))/7.)/slc
        c010 = ((wc(1)+wc(2)+wc(3)+wc(4)+wc(6)+wc(7)+wc(8))/7.)/slc
        c110 = ((wc(1)+wc(2)+wc(3)+wc(4)+wc(5)+wc(7)+wc(8))/7.)/slc
        c011 = ((wc(1)+wc(2)+wc(3)+wc(4)+wc(5)+wc(6)+wc(8))/7.)/slc
        c111 = ((wc(1)+wc(2)+wc(3)+wc(4)+wc(5)+wc(6)+wc(7))/7.)/slc
        d000 = ((wd(2)+wd(3)+wd(4)+wd(5)+wd(6)+wd(7)+wd(8))/7.)/sld
        d100 = ((wd(1)+wd(3)+wd(4)+wd(5)+wd(6)+wd(7)+wd(8))/7.)/sld
        d001 = ((wd(1)+wd(2)+wd(4)+wd(5)+wd(6)+wd(7)+wd(8))/7.)/sld
        d101 = ((wd(1)+wd(2)+wd(3)+wd(5)+wd(6)+wd(7)+wd(8))/7.)/sld
        d010 = ((wd(1)+wd(2)+wd(3)+wd(4)+wd(6)+wd(7)+wd(8))/7.)/sld
        d110 = ((wd(1)+wd(2)+wd(3)+wd(4)+wd(5)+wd(7)+wd(8))/7.)/sld
        d011 = ((wd(1)+wd(2)+wd(3)+wd(4)+wd(5)+wd(6)+wd(8))/7.)/sld
        d111 = ((wd(1)+wd(2)+wd(3)+wd(4)+wd(5)+wd(6)+wd(7))/7.)/sld
        if (abs(1.-(c000+c001+c100+c101+c010+c011+c110+c111)).gt..0001)
     &          write(*,*) "WARNING: cumulative wc = ",
     &          c000+c001+c100+c101+c010+c011+c110+c111

        if (abs(1.-(d000+d001+d100+d101+d010+d011+d110+d111)).gt..0001)
     &          write(*,*) "WARNING: cumulative wd = ",
     &          d000+d001+d100+d101+d010+d011+d110+d111
        xlen = xss-x0
        ylen = yss-y0
        zlen = zss-z0
        if (ipolar .eq. 1) then     ! x,y distances decr w/ depth
          xlen = xlen*(1.-(zss+z0)/(2.*re))
          ylen = ylen*(1.-(zss+z0)/(2.*re))
        endif
        sslen = sqrt((xlen**2) + (ylen**2) + (zlen**2))
        wtc = 0.5
c        write(*,*) "Length of the subsegment is ",sslen,icase
c        write(*,*) "wtc*(d110+c110) = ", wtc*(d110+c110)
c        alllen=alllen+sslen

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      Qp: grow(element)=grow(element)+.5*(wgt)*distn*Pslo             C
C      Qs: grow(element)=grow(element)+.5*(wgt)*distn*Sslo             C
C   Qp/Qs: grow(element)=grow(element)+.5*(wgt)*distn*Sslo/Qp          C
C      Qk: grow(element)=grow(element)+.5*(wgt)*(1-3/4*(Pslo/Sslo)**2) C
C           wgt = some "trapezoidal weighting" function                C
C           distn = pathlength in that block                           C
C           ?slo = P or S slowness in that block                       C
C           element = pt-th block, counting in the order of z, x, y    C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        if ((id .eq. 1) .or. (id .eq. 2)) then          ! Qp or Qs
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0)+wtc*(d000+c000)*
     &        sslen*sn000
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1)+wtc*(d001+c001)
     &        *sslen*sn001
          grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) + wtc*(d010+c010)*
     &        sslen*sn010
          grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) + wtc*(d011+c011)*
     &        sslen*sn011
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d100+c100)*
     &        sslen*sn100
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d101+c101)*
     &        sslen*sn101
          grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d110+c110)*
     &        sslen*sn110
          grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d111+c111)*
     &        sslen*sn111
        else if (id .eq. 3) then    ! Qp/Qs
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0)+wtc*(d000+c000)*
     &        sslen*sn000*attin(kq0,jq0,iq0)
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1)+wtc*(d001+c001)
     &        *sslen*sn001*attin(kq0+1,jq0,iq0)
          grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) + wtc*(d010+c010)*
     &        sslen*sn010*attin(kq0,jq0,iq0+1)
          grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) + wtc*(d011+c011)*
     &        sslen*sn011*attin(kq0+1,jq0,iq0+1)
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d100+c100)*
     &        sslen*sn100*attin(kq0,jq0+1,iq0)
          grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
     &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d101+c101)*
     &        sslen*sn101*attin(kq0+1,jq0+1,iq0)
          grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d110+c110)*
     &        sslen*sn110*attin(kq0,jq0+1,iq0+1)
          grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
     &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d111+c111)*
     &        sslen*sn111*attin(kq0+1,jq0+1,iq0+1)
C       else if (id .eq. 4) then    ! Qp/Qs: sn = Vs/Vp
C         grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
C    &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0)+wtc*(d000+c000)*
C    &        (1-3./4.*sn000*sn000)
C         grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
C    &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1)+wtc*(d001+c001)
C    &        *(1-3./4.*sn001*sn001)
C         grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) =
C    &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0) + wtc*(d010+c010)*
C    &        (1-3./4.*sn010*sn010)
C         grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) =
C    &        grow(jq0*nnzQ+nnyQ*nnzQ*(iq0-1)+kq0+1) + wtc*(d011+c011)*
C    &        (1-3./4.*sn011*sn011)
C         grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) =
C    &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d100+c100)*
C    &        (1-3./4.*sn100*sn100)
C         grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
C    &        grow((jq0-1)*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d101+c101)*
C    &        (1-3./4.*sn101*sn101)
C         grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) =
C    &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0) + wtc*(d110+c110)*
C    &        (1-3./4.*sn110*sn110)
C         grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) =
C    &        grow(jq0*nnzQ+nnyQ*nnzQ*iq0+kq0+1) + wtc*(d111+c111)*
C    &        (1-3./4.*sn111*sn111)
        endif

        x0 = xss
        y0 = yss
        z0 = zss
c        write(*,*) x0, y0, z0
        if (z0 .eq. zarc(ipt1)) then
          pt = pt+1
        elseif (icase .eq. 1) then      ! Move to next subsegment for case 1
          pt = pt + 1
          x1 = xarc(pt)
          y1 = yarc(pt)
          z1 = zarc(pt)
        endif
c        write(*,*) x0,y0,z0,pt,icase
c        write(*,*) "pt, npts+1: ", pt, npts+1
      enddo       ! end of loop over each raypath node (pt -> ipt1)
c      write(*,*) "alllen2 = :",alllen
237   continue
c      write(*,400) ifskip,npts+2,pt-1,ipt,ipt1,
c     &                  zarc(pt-1),zarc(0),zarc(ipt),zarc(ipt1)
c      write(*,*) "Finish adding row of G matrix"
200   format(3(f10.4,1x))
300   format(a15,3(f10.4,1x))
400   format(5(i4,1x),4(f10.4,1x))
      END

max_stars_repo_name: swei-seismo/Qtomo

max_issues_repo_name: swei-seismo/Qtomo

max_forks_repo_name: swei-seismo/Qtomo