text:
module kernel_flux_mod

use constants_mod

use mesh_types_mod, only: quad_t

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

contains

subroutine kernel_fill_facevars(quad)

    use mesh_types_mod, only: quad_t

    type(quad_t), intent(inout) :: quad

    integer :: i,j,h

    integer, parameter :: M = N_NODES-1
    integer, parameter :: N = N_NODES

    !! Convention: increaseing indices denotes further distance from interface
    associate(fvstate => quad%pld%hydro%state, fvars => quad%aux%hydro%facevars)
    do h = 1,N_VARS; do j = 1,N_NODES; do i = 1,N_NODES;
        fvars(i,j,h,2,SOU) = fvstate(M,i,j,h)
        fvars(i,j,h,1,SOU) = fvstate(N,i,j,h)

        fvars(i,j,h,1,NOR) = fvstate(1,i,j,h)
        fvars(i,j,h,2,NOR) = fvstate(2,i,j,h)

        fvars(i,j,h,2,EAS) = fvstate(i,M,j,h)
        fvars(i,j,h,1,EAS) = fvstate(i,N,j,h)

        fvars(i,j,h,1,WES) = fvstate(i,1,j,h)
        fvars(i,j,h,2,WES) = fvstate(i,2,j,h)

        fvars(i,j,h,2,BAC) = fvstate(i,j,M,h)
        fvars(i,j,h,1,BAC) = fvstate(i,j,N,h)

        fvars(i,j,h,1,FRO) = fvstate(i,j,1,h)
        fvars(i,j,h,2,FRO) = fvstate(i,j,2,h)
    end do; end do; end do;
    end associate

end subroutine

subroutine kernel_fill_faceflux(side)

    use mesh_types_mod, only: quad_t, side_t
    use mesh_const_mod, only: MESH_HALF, MESH_FACES

    use riemann_inner_mod, only: xriemann_inner => xriemann
    use riemann_inner_mod, only: yriemann_inner => yriemann
    use riemann_inner_mod, only: zriemann_inner => zriemann

# if !PP_MESH_PERIODIC
    use riemann_north_mod, only: riemann_north => xriemann
    use riemann_south_mod, only: riemann_south => xriemann

    use riemann_west_mod, only: riemann_west => yriemann
    use riemann_east_mod, only: riemann_east => yriemann

    use riemann_front_mod, only: riemann_front => zriemann
    use riemann_back_mod, only: riemann_back => zriemann
# endif

    use setup_boundary_mod, only: setup_side_boundary

    use kernel_share_mod, only: rfaceMat,rfaceTam
    use kernel_share_mod, only: cfaceMat,cfaceTam

    use equations_mod, only: cons2prim, prim2cons
    use equations_const_mod, only: posvars
    use utils_mod, only: vmax

    use kernel_config_mod, only: lofactor
    use math_mod, only: minmod3

# if PP_NON_CONSERVATIVE_FLUX
    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_xFlux
    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_yFlux
    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_zFlux
# endif

    integer, parameter :: N = N_NODES
    integer, parameter :: V = N_VARS

    type(side_t), intent(inout) :: side

    real(dp) :: recovars(N_NODES,N_NODES,N_VARS,2,MESH_HALF,2)
    real(dp) :: facevars(N_NODES,N_NODES,N_VARS,  MESH_HALF,2)
    real(dp) :: flexvars(N_NODES,N_NODES,N_VARS,  MESH_HALF,2)
    real(dp) :: fluxvars(N_NODES,N_NODES,N_VARS,  MESH_HALF,2)

    real(dp) :: ocervars(N_NODES,N_NODES,N_VARS,2,MESH_HALF,2)

    integer :: i,j,k,h,f, iside, iquad,nquads

    !! Reconstruct face data for involved orders for conforming and non-conforming faces.
    do iside = 1,2
        !! non-conforming face
        if (side%nquads(iside) > 0 .and. side%nquads(iside) < side%nquads(3-iside)) then
            do iquad = 1,MESH_HALF
                do j = 1,2; do h = 1,N_VARS;
                    recovars(:,:,h,j,iquad,iside) &
                        = matmul(rfaceMat(:,:,iquad),&
                          matmul(side%quads(1,iside)%p%aux%hydro%facevars(:,:,h,j,side%faceid(iside)),rfaceTam(:,:,iquad)))
                end do; end do;
                recovars(:,:,:,2,iquad,iside) = 2.0_dp/3.0_dp*recovars(:,:,:,1,iquad,iside) + 1.0_dp/3.0_dp*recovars(:,:,:,2,iquad,iside)
            end do

        !! conforming face
        else
            do iquad = 1,side%nquads(iside)
                recovars(:,:,:,:,iquad,iside) = side%quads(iquad,iside)%p%aux%hydro%facevars(:,:,:,:,side%faceid(iside))
            end do
        end if
    end do

    nquads = maxval(side%nquads)

    !! Mirror inner states at boundaries.
    if (side%boundary > 0) then
        recovars(:,:,:,:,:,side%outer) = recovars(:,:,:,:,:,side%inner)
    end if

    do iside = 1,2
        do iquad = 1,side%nquads(iside)
            side%quads(iquad,iside)%p%aux%hydro%halovars(:,:,:,side%faceid(iside)) = recovars(:,:,:,1,iquad,3-iside)
        end do
    end do

    ocervars = recovars

    !! Transform to primitive variables.
    do iside = 1,2; do iquad = 1,nquads; do k = 1,2; do j = 1,N_NODES; do i = 1,N_NODES;
        recovars(i,j,:,k,iquad,iside) = cons2prim(recovars(i,j,:,k,iquad,iside))
    end do; end do; end do; end do; end do;

    !! Reconstruct via 2nd order monotonized minmod limiter.
    do iquad = 1,nquads
        do h = 1,N_VARS
            facevars(:,:,h,iquad,side%plus) = recovars(:,:,h,1,iquad,side%plus) + 0.25_dp*minmod3(&
                4.0_dp*(recovars(:,:,h,1,iquad,side%plus )-recovars(:,:,h,2,iquad,side%plus)), &
                       (recovars(:,:,h,1,iquad,side%minus)-recovars(:,:,h,2,iquad,side%plus)), &
                4.0_dp*(recovars(:,:,h,1,iquad,side%minus)-recovars(:,:,h,1,iquad,side%plus)))
        end do
    end do

    do iquad = 1,nquads
        do h = 1,N_VARS
            facevars(:,:,h,iquad,side%minus) = recovars(:,:,h,1,iquad,side%minus) - 0.25_dp*minmod3(&
                4.0_dp*(recovars(:,:,h,1,iquad,side%minus)-recovars(:,:,h,1,iquad,side%plus)), &
                       (recovars(:,:,h,2,iquad,side%minus)-recovars(:,:,h,1,iquad,side%plus)), &
                4.0_dp*(recovars(:,:,h,2,iquad,side%minus)-recovars(:,:,h,1,iquad,side%minus)))
        end do
    end do

    !! Reconstruct via 2nd order monotonized minmod limiter.
    do iside = 1,2; do iquad = 1,nquads; do h = 1,size(posvars);
        facevars(:,:,posvars(h),iquad,iside) = vmax(lofactor*recovars(:,:,posvars(h),1,iquad,iside),facevars(:,:,posvars(h),iquad,iside))
    end do; end do; end do;

    !! Transform to conservative variables.
    do iside = 1,2; do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
        facevars(i,j,:,iquad,iside) = prim2cons(facevars(i,j,:,iquad,iside))
    end do; end do; end do; end do;

    !! Process boundary conditions.
    if (side%boundary > 0) then
        call setup_side_boundary(side,facevars(:,:,:,1,side%inner),facevars(:,:,:,1,side%outer))
    end if

    !! Calculate face fluxes.
    if (side%boundary > 0) then
# if !PP_MESH_PERIODIC
        select case(side%boundary)
            case (NOR)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_north(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

            case (SOU)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_south(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

            case (WES)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_west(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

            case (EAS)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_east(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

           case (FRO)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_front(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

           case (BAC)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = riemann_back(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;
        end select
# endif
    else
        select case(side%direction)
            case (X_DIR)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = xriemann_inner(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

            case (Y_DIR)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = yriemann_inner(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;

            case (Z_DIR)
                do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                    flexvars(i,j,:,iquad,1) = zriemann_inner(facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
                end do; end do; end do;
        end select
    end if

    flexvars(:,:,:,:,2) = flexvars(:,:,:,:,1)

# if PP_NON_CONSERVATIVE_FLUX
    select case(side%direction)
        case (X_DIR)
            do iside = 1,2; do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                flexvars(i,j,:,iquad,iside) = flexvars(i,j,:,iquad,iside) + Non_Conservative_two_point_xFlux(&
                    ocervars(i,j,:,1,iquad,iside),facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
            end do; end do; end do; end do;

        case (Y_DIR)
            do iside = 1,2; do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                flexvars(i,j,:,iquad,iside) = flexvars(i,j,:,iquad,iside) + Non_Conservative_two_point_yFlux(&
                    ocervars(i,j,:,1,iquad,iside),facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
            end do; end do; end do; end do;

        case (Z_DIR)
            do iside = 1,2; do iquad = 1,nquads; do j = 1,N_NODES; do i = 1,N_NODES;
                flexvars(i,j,:,iquad,iside) = flexvars(i,j,:,iquad,iside) + Non_Conservative_two_point_zFlux(&
                    ocervars(i,j,:,1,iquad,iside),facevars(i,j,:,iquad,side%plus),facevars(i,j,:,iquad,side%minus))
            end do; end do; end do; end do;
    end select
# endif

    !! Project fluxes back to block faces.
    do iside = 1,2
        if (side%nquads(iside) > 0 .and. side%nquads(iside) < side%nquads(3-iside)) then
            do iquad = 1,MESH_HALF
                forall (h = 1:N_VARS)
                    fluxvars(:,:,h,iquad,iside) = matmul(cfaceMat(:,:,iquad),matmul(flexvars(:,:,h,iquad,iside),cfaceTam(:,:,iquad)))
                end forall
            end do

            do iquad = 2,MESH_HALF
                fluxvars(:,:,:,1,iside) = fluxvars(:,:,:,1,iside) + fluxvars(:,:,:,iquad,iside)
            end do
        else
            do iquad = 1,side%nquads(iside)
                fluxvars(:,:,:,iquad,iside) = flexvars(:,:,:,iquad,iside)
            end do
        end if
    end do

    !! Copy facefluxes to respective quads.
    do iside = 1,2
        do iquad = 1,side%nquads(iside)
            side%quads(iquad,iside)%p%aux%hydro%faceflux(:,:,:,side%faceid(iside)) = fluxvars(:,:,:,iquad,iside)
        end do
    end do

end subroutine

subroutine kernel_calc_flux(quad,fvrhs)

    use mesh_types_mod, only: quad_t
    use riemann_inner_mod, only: xriemann
    use riemann_inner_mod, only: yriemann
    use riemann_inner_mod, only: zriemann

    use equations_mod, only: cons2prim, prim2cons
    use equations_const_mod, only: posvars
    use utils_mod, only: vmax

    use kernel_config_mod, only: lofactor
    use math_mod, only: minmod3

# if PP_NON_CONSERVATIVE_FLUX
    use equations_share_mod, only: GLM_alpha
    use equations_const_mod, only: GLMP_VAR

    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_xFlux
    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_yFlux
    use Non_Conservative_flux_mod, only: Non_Conservative_two_point_zFlux
# endif

    integer, parameter :: R = N_NODES-1
    integer, parameter :: N = N_NODES
    integer, parameter :: P = N_NODES+1
    integer, parameter :: V = N_VARS

    type(quad_t), intent(inout) :: quad
    real(dp), intent(out)       :: fvrhs(N_NODES,N_NODES,N_NODES,N_VARS)

    !! lower thresholds
    real(dp) :: xlob(0:P,N,N,size(posvars))
    real(dp) :: ylob(N,0:P,N,size(posvars))
    real(dp) :: zlob(N,N,0:P,size(posvars))

    !! stencils for reconstruction
    real(dp) :: xrec(0:P,N,N,N_VARS)
    real(dp) :: yrec(N,0:P,N,N_VARS)
    real(dp) :: zrec(N,N,0:P,N_VARS)

    !! inner interfaces
    real(dp) :: xuM(R,N,N,N_VARS)
    real(dp) :: xuP(R,N,N,N_VARS)

    real(dp) :: yuM(N,R,N,N_VARS)
    real(dp) :: yuP(N,R,N,N_VARS)

    real(dp) :: zuM(N,N,R,N_VARS)
    real(dp) :: zuP(N,N,R,N_VARS)

    !! interface fluxes
    real(dp) :: xrm(0:N_NODES,1:N_NODES,1:N_NODES,N_VARS)
    real(dp) :: yrm(1:N_NODES,0:N_NODES,1:N_NODES,N_VARS)
    real(dp) :: zrm(1:N_NODES,1:N_NODES,0:N_NODES,N_VARS)

    integer :: i,j,k,h

    real(dp) :: sdx,sdy,sdz

    sdx = REAL(N_NODES,dp)/quad%delta(1)
    sdy = REAL(N_NODES,dp)/quad%delta(2)
    sdz = REAL(N_NODES,dp)/quad%delta(3)

    !! ---------------------------------------------------------------------- !!
    !! Interpolating inner interfaces. Monotonized centrel reconstruction.

    do h = 1,V; do k = 1,N; do j = 1,N; do i = 1,N;
        xrec(i,j,k,h) = quad%pld%hydro%state(i,j,k,h)
        yrec(i,j,k,h) = quad%pld%hydro%state(i,j,k,h)
        zrec(i,j,k,h) = quad%pld%hydro%state(i,j,k,h)
    end do; end do; end do; end do;

    xrec(0,:,:,:) = quad%aux%hydro%halovars(:,:,:,NOR)
    xrec(P,:,:,:) = quad%aux%hydro%halovars(:,:,:,SOU)
    yrec(:,0,:,:) = quad%aux%hydro%halovars(:,:,:,WES)
    yrec(:,P,:,:) = quad%aux%hydro%halovars(:,:,:,EAS)
    zrec(:,:,0,:) = quad%aux%hydro%halovars(:,:,:,FRO)
    zrec(:,:,P,:) = quad%aux%hydro%halovars(:,:,:,BAC)

    !! Transform to primitive variables.
    do k = 1,N; do j = 1,N; do i = 0,P;
        xrec(i,j,k,:) = cons2prim(xrec(i,j,k,:))
    end do; end do; end do;

    do k = 1,N; do j = 0,P; do i = 1,N;
        yrec(i,j,k,:) = cons2prim(yrec(i,j,k,:))
    end do; end do; end do;

    do k = 0,P; do j = 1,N; do i = 1,N;
        zrec(i,j,k,:) = cons2prim(zrec(i,j,k,:))
    end do; end do; end do;

    !! Reconstruct interfaces.
    do h = 1,V; do k = 1,N; do j = 1,N; do i = 1,R;
        xuP(i,j,k,h) = xrec(i  ,j,k,h) + 0.25_dp*minmod3(&
            4.0_dp*(xrec(i  ,j,k,h) - xrec(i-1,j,k,h)), &
                   (xrec(i+1,j,k,h) - xrec(i-1,j,k,h)), &
            4.0_dp*(xrec(i+1,j,k,h) - xrec(i  ,j,k,h)))

        xuM(i,j,k,h) = xrec(i+1,j,k,h) - 0.25_dp*minmod3(&
            4.0_dp*(xrec(i+1,j,k,h) - xrec(i  ,j,k,h)), &
                   (xrec(i+2,j,k,h) - xrec(i  ,j,k,h)), &
            4.0_dp*(xrec(i+2,j,k,h) - xrec(i+1,j,k,h)))
    end do; end do; end do; end do;

    do h = 1,V; do k = 1,N; do j = 1,R; do i = 1,N;
        yuP(i,j,k,h) = yrec(i,j  ,k,h) + 0.25_dp*minmod3(&
            4.0_dp*(yrec(i,j  ,k,h) - yrec(i,j-1,k,h)), &
                   (yrec(i,j+1,k,h) - yrec(i,j-1,k,h)), &
            4.0_dp*(yrec(i,j+1,k,h) - yrec(i,j  ,k,h)))

        yuM(i,j,k,h) = yrec(i,j+1,k,h) - 0.25_dp*minmod3(&
            4.0_dp*(yrec(i,j+1,k,h) - yrec(i,j  ,k,h)), &
                   (yrec(i,j+2,k,h) - yrec(i,j  ,k,h)), &
            4.0_dp*(yrec(i,j+2,k,h) - yrec(i,j+1,k,h)))
    end do; end do; end do; end do;

    do h = 1,V; do k = 1,R; do j = 1,N; do i = 1,N;
        zuP(i,j,k,h) = zrec(i,j,k  ,h) + 0.25_dp*minmod3(&
            4.0_dp*(zrec(i,j,k  ,h) - zrec(i,j,k-1,h)), &
                   (zrec(i,j,k+1,h) - zrec(i,j,k-1,h)), &
            4.0_dp*(zrec(i,j,k+1,h) - zrec(i,j,k  ,h)))

        zuM(i,j,k,h) = zrec(i,j,k+1,h) - 0.25_dp*minmod3(&
            4.0_dp*(zrec(i,j,k+1,h) - zrec(i,j,k  ,h)), &
                   (zrec(i,j,k+2,h) - zrec(i,j,k  ,h)), &
            4.0_dp*(zrec(i,j,k+2,h) - zrec(i,j,k+1,h)))
    end do; end do; end do; end do;

    !! Limit interface states to lower threshold.
    do h = 1,size(posvars);
        xlob(:,:,:,h) = lofactor*xrec(:,:,:,posvars(h))
        ylob(:,:,:,h) = lofactor*yrec(:,:,:,posvars(h))
        zlob(:,:,:,h) = lofactor*zrec(:,:,:,posvars(h))
    end do

    forall (h = 1:size(posvars), k = 1:N, j = 1:N, i = 1:R)
        xuP(i,j,k,posvars(h)) = vmax(xlob(i,  j,k,h),xuP(i,j,k,posvars(h)))
        xuM(i,j,k,posvars(h)) = vmax(xlob(i+1,j,k,h),xuM(i,j,k,posvars(h)))
    end forall

    forall (h = 1:size(posvars), k = 1:N, j = 1:R, i = 1:N)
        yuP(i,j,k,posvars(h)) = vmax(ylob(i,j  ,k,h),yuP(i,j,k,posvars(h)))
        yuM(i,j,k,posvars(h)) = vmax(ylob(i,j+1,k,h),yuM(i,j,k,posvars(h)))
    end forall

    forall (h = 1:size(posvars), k = 1:R, j = 1:N, i = 1:N)
        zuP(i,j,k,posvars(h)) = vmax(zlob(i,j,k  ,h),zuP(i,j,k,posvars(h)))
        zuM(i,j,k,posvars(h)) = vmax(zlob(i,j,k+1,h),zuM(i,j,k,posvars(h)))
    end forall

    !! Convert back to conservative variables.
    forall (k = 1:N, j = 1:N, i = 1:R)
        xuP(i,j,k,:) = prim2cons(xuP(i,j,k,:))
        xuM(i,j,k,:) = prim2cons(xuM(i,j,k,:))
    end forall

    forall (k = 1:N, j = 1:R, i = 1:N)
        yuP(i,j,k,:) = prim2cons(yuP(i,j,k,:))
        yuM(i,j,k,:) = prim2cons(yuM(i,j,k,:))
    end forall

    forall (k = 1:R, j = 1:N, i = 1:N)
        zuP(i,j,k,:) = prim2cons(zuP(i,j,k,:))
        zuM(i,j,k,:) = prim2cons(zuM(i,j,k,:))
    end forall

    !! ---------------------------------------------------------------------- !!
    !! Solve Riemann problems.

    forall (k = 1:N, j = 1:N, i = 1:R) xrm(i,j,k,:) = xriemann(xuP(i,j,k,:),xuM(i,j,k,:))
    forall (k = 1:N, j = 1:R, i = 1:N) yrm(i,j,k,:) = yriemann(yuP(i,j,k,:),yuM(i,j,k,:))
    forall (k = 1:R, j = 1:N, i = 1:N) zrm(i,j,k,:) = zriemann(zuP(i,j,k,:),zuM(i,j,k,:))

    !! ---------------------------------------------------------------------- !!
    !! Append block face fluxes.

    xrm(0,:,:,:) = quad%aux%hydro%faceflux(:,:,:,NOR)
    xrm(N,:,:,:) = quad%aux%hydro%faceflux(:,:,:,SOU)

    yrm(:,0,:,:) = quad%aux%hydro%faceflux(:,:,:,WES)
    yrm(:,N,:,:) = quad%aux%hydro%faceflux(:,:,:,EAS)

    zrm(:,:,0,:) = quad%aux%hydro%faceflux(:,:,:,FRO)
    zrm(:,:,N,:) = quad%aux%hydro%faceflux(:,:,:,BAC)

    !! ---------------------------------------------------------------------- !!
    !! Apply FV operator.

    xrm = sdx*xrm
    yrm = sdy*yrm
    zrm = sdz*zrm

    forall (h = 1:V, k = 1:N, j = 1:N, i = 1:N)
        fvrhs(i,j,k,h) = xrm(i-1,j,k,h)-xrm(i,j,k,h) + yrm(i,j-1,k,h)-yrm(i,j,k,h) + zrm(i,j,k-1,h)-zrm(i,j,k,h)
    end forall

# if PP_NON_CONSERVATIVE_FLUX
    associate(u => quad%pld%hydro%state)
    do k = 1,N; do j = 1,N; do i = 1,R;
        fvrhs(i  ,j,k,:) = fvrhs(i  ,j,k,:) - sdx*Non_Conservative_two_point_xFlux(u(i  ,j,k,:),xuP(i,j,k,:),xuM(i,j,k,:))
        fvrhs(i+1,j,k,:) = fvrhs(i+1,j,k,:) + sdx*Non_Conservative_two_point_xFlux(u(i+1,j,k,:),xuP(i,j,k,:),xuM(i,j,k,:))
    end do; end do; end do;

    do k = 1,N; do j = 1,R; do i = 1,N;
        fvrhs(i,j  ,k,:) = fvrhs(i,j  ,k,:) - sdy*Non_Conservative_two_point_yFlux(u(i,j  ,k,:),yuP(i,j,k,:),yuM(i,j,k,:))
        fvrhs(i,j+1,k,:) = fvrhs(i,j+1,k,:) + sdy*Non_Conservative_two_point_yFlux(u(i,j+1,k,:),yuP(i,j,k,:),yuM(i,j,k,:))
    end do; end do; end do;

    do k = 1,R; do j = 1,N; do i = 1,N;
        fvrhs(i,j,k  ,:) = fvrhs(i,j,k  ,:) - sdz*Non_Conservative_two_point_zFlux(u(i,j,k  ,:),zuP(i,j,k,:),zuM(i,j,k,:))
        fvrhs(i,j,k+1,:) = fvrhs(i,j,k+1,:) + sdz*Non_Conservative_two_point_zFlux(u(i,j,k+1,:),zuP(i,j,k,:),zuM(i,j,k,:))
    end do; end do; end do;
    end associate

    !! GLM damping term.
    fvrhs(:,:,:,GLMP_VAR) = fvrhs(:,:,:,GLMP_VAR) - GLM_alpha*quad%pld%hydro%state(:,:,:,GLMP_VAR)
# endif

end subroutine

end module

max_stars_repo_name: jmark/nemo

max_issues_repo_name: jmark/nemo

max_forks_repo_name: jmark/nemo


text:
 ... because you have nothing else better to do.

Describe SeanRiddle here.

Sean went to Bishop ODowd High School ... the poor basterd.

Sean is a Computer Science and Engineering CSE major (for now) ... poor bastard.

Sean lives in West Davis ... poor bastard.

Seans Senior Quote was, Dont drink and derive. ... poor bastard.

Seans http://www.xanga.com/home.aspx?userJustifiedByMargins Xanga screen name is JustifiedByMargins ... poor bastard.

Seans favorite reason for doing anything is, Shut up, I know it wont work! ... poor bastard.

In effect, the best two word description of Sean would quite possibly be ... a great guy.

... poor bastard.

Q.E.D.

There are many who believe that Sean Riddle is the secret identity of Genghis Khan.  Riddle by day, Khan by night.  Sean is plotting to end world hunger and establish world peace.

But evidence . . .

See, totally the same person!  Although not identical, you cannot deny the similarities!  There you have it, proof the Sean is Khan

SSSEEEAAANNNNN!!!!!!!!!!!!!!!!!!!!!!!

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program driver2
   use KindParamModule, only : IntKind, RealKind
   use MathParamModule, only : PI, THREE, THIRD, PI4
   use MPPModule, only : initMPP, endMPP
   use MPPModule, only : syncAllPEs
   use MPPModule, only : MyPE, NumPEs, AllPEs, AnyPE
   use MPPModule, only : bcastMessage, sendMessage, recvMessage
   use MPPModule, only : packMessage, unpackMessage, sendPackage, recvPackage
   use PolyhedraModule, only : initPolyhedra, genPolyhedron, endPolyhedra
   use PolyhedraModule, only : getNeighborIndex, getNumPlanes, getPlaneArea
   use PolyhedraModule, only : getInscrSphRadius, getOutscrSphRadius, getVolume
   use PolyhedraModule, only : getWignerSeitzRadius, getBoxVolume
   use PolyhedraModule, only : printPolyhedron, getNeighborDistance
   implicit none
!
   character (len=70) :: text
   character (len=70) :: path
   character (len=30) :: info_table
!
   integer (kind=IntKind) :: num_atoms
   integer (kind=IntKind) :: i, j, k, lmax, n
   integer (kind=IntKind) :: mtasa
   integer (kind=IntKind) :: ngaussr, ngaussq, np
   integer (kind=IntKind), allocatable :: alist(:)
   integer (kind=IntKind), allocatable :: anum(:)
!
   real (kind=RealKind) :: bravais(3,3)
   real (kind=RealKind) :: rins, rins0
   real (kind=RealKind) :: rcut, r1, r2, r3, r4, rmt, rws, rcirc
   real (kind=RealKind), allocatable :: position(:,:)
   real (kind=RealKind), allocatable :: rad(:)
   real (kind=RealKind), allocatable :: area(:), dist(:)
   real (kind=RealKind) :: omegmt, omegws, omegbra, v
!
   call initMPP()
!
   if (MyPE == 0) then
      do i=1,6
         read(5,'(a)')text
      enddo
      read(5,'(a)')path
      do i=1,13
         read(5,'(a)')text
      enddo
      read(5,*)num_atoms
      read(5,'(a)')text
      read(5,*)i,j,mtasa,rins
      do i=1,5
         read(5,'(a)')text
      enddo
      read(5,*)ngaussr, ngaussq
      read(5,'(a)')text
      read(5,'(a)')info_table
!     ----------------------------------------------------------------
      call packMessage(num_atoms)
      call packMessage(mtasa)
      call packMessage(ngaussr)
      call packMessage(ngaussq)
      call packMessage(rins)
      call sendPackage(1000,AllPEs)
!     ----------------------------------------------------------------
   else
!     ----------------------------------------------------------------
      call recvPackage(1000,0)
      call unpackMessage(num_atoms)
      call unpackMessage(mtasa)
      call unpackMessage(ngaussr)
      call unpackMessage(ngaussq)
      call unpackMessage(rins)
!     ----------------------------------------------------------------
   endif
!
   allocate( position(3,num_atoms) )
   allocate( rad(num_atoms) )
!
   if (MyPE == 0) then
      allocate( anum(num_atoms) )
      open(unit=10,file=trim(path)//info_table,status='old')
      read(10,'(a)')text
      read(10,'(a)')text
      read(10,'(a)')text
      read(10,*)bravais(1:3,1)
      read(10,*)bravais(1:3,2)
      read(10,*)bravais(1:3,3)
!     ----------------------------------------------------------------
      call bcastMessage(bravais,3,3,0)
!     ----------------------------------------------------------------
      do i=1,5
         read(10,'(a)')text
      enddo
      do i=1,num_atoms
         read(10,*)anum(i),k,position(1,i),position(2,i),position(3,i),lmax, &
                   rcut,r1,r2,r3,r4,rad(i)
      enddo
!     ----------------------------------------------------------------
      call bcastMessage(position,3,num_atoms,0)
      call bcastMessage(rad,num_atoms,0)
!     ----------------------------------------------------------------
      close(unit=10)
   else
!     ----------------------------------------------------------------
      call bcastMessage(bravais,3,3,0)
      call bcastMessage(position,3,num_atoms,0)
      call bcastMessage(rad,num_atoms,0)
!     ----------------------------------------------------------------
   endif
!
   n = 0
   do i=1,num_atoms
      if (mod(i-1,NumPEs) == MyPE) then
         n = n+1
      endif
   enddo
!
   call initPolyhedra(n,bravais,'none',-1)
!
   if (MyPE == 0) then
      open(unit=11,file=trim(path)//'voronoi.dat',status='unknown')
      open(unit=12,file=trim(path)//'neighbor.dat',status='unknown')
      write(11,'(80(''=''))')
      write(11,'(a)')            &
      ' Atom    No. Surf     Vcell       Vmt         Rmt       Rcirc       Rws'
      write(11,'(80(''-''))')
   endif
!
!  -------------------------------------------------------------------
   call syncAllPEs()
!  -------------------------------------------------------------------
!
   rins0 = rins
!  -------------------------------------------------------------------
   omegbra = getBoxVolume()
!  -------------------------------------------------------------------
   v = 0.0d0
   n = 0
   do i=1,num_atoms
      if (mod(i-1,NumPEs) == MyPE) then
         rins = rins0
         n = n+1
!        -------------------------------------------------------------
         call genPolyhedron(n,i,num_atoms,position,rad)
!        -------------------------------------------------------------
         np = getNumPlanes(n)
         rmt = getInscrSphRadius(n)
         rws = getWignerSeitzRadius(n)
         rcirc = getOutscrSphRadius(n)
         omegws = getVolume(n)
!        -------------------------------------------------------------
         allocate(alist(np), area(np), dist(np))
         do j=1,np
            alist(j) = getNeighborIndex(n,j)
            area(j)  = getPlaneArea(n,j)
            dist(j)  = getNeighborDistance(n,j)
         enddo
         omegmt = PI4*THIRD*rmt**3
!
         if (mtasa >= 1) then
            if(rins > 0.0d0) then
!              if rins is preset in the input then force all rws the same
               rws=(THREE*omegbra/(PI4*num_atoms))**THIRD
               omegws=PI4*rws**3/THREE
            else
!              otherwise rins is the real rmt
               rins=rmt
            endif
            rmt = rws
            omegmt = omegws
            if (mtasa == 1) then
               rins=rws
            endif
         else if (rins <= 0.0d0) then
!           otherwise rins is the real rmt
            rins=rmt
         endif
!
         if (MyPE /= 0) then
            call packMessage(omegws)
            call packMessage(omegmt)
            call packMessage(rmt)
            call packMessage(rcirc)
            call packMessage(rws)
            call sendPackage(3000+i,0)
            call sendMessage(np,4000+i,0)
            call sendMessage(alist,np,5000+i,0)
            call sendMessage(area,np,6000+i,0)
            call sendMessage(dist,np,7000+i,0)
         endif
      else if (MyPE == 0) then
         j = mod(i-1,NumPEs)
         call recvPackage(3000+i,j)
         call unpackMessage(omegws)
         call unpackMessage(omegmt)
         call unpackMessage(rmt)
         call unpackMessage(rcirc)
         call unpackMessage(rws)
         call recvMessage(np,4000+i,j)
         allocate( alist(np), area(np), dist(np) )
         call recvMessage(alist,np,5000+i,j)
         call recvMessage(area,np,6000+i,j)
         call recvMessage(dist,np,7000+i,j)
      endif
      if (MyPE == 0) then
         v = v + omegws
         write(11,'(i5,4x,i5,4x,5(1x,f10.5))')i,np,omegws,omegmt,rmt,rcirc,rws
         write(12,'(3i5))')i,anum(i),np
         write(12,'(3(i5,2f10.5,a))')(alist(j),area(j),dist(j),';',j=1,np)
      endif
      if (allocated(alist)) then
         deallocate( alist )
      endif
      if (allocated(area)) then
         deallocate( area )
      endif
      if (allocated(dist)) then
         deallocate( dist )
      endif
   enddo
   if (MyPE == 0) then
      close(unit=11)
      close(unit=12)
      deallocate( anum )
      print *,'v = ',v
      print *,'omegbra = ',omegbra
      print *,'Volume difference = ',abs(v-omegbra)
   endif
!
   deallocate( position, rad )
   call endPolyhedra()
   call endMPP()
end program driver2

max_forks_repo_name: hkershaw-brown/MuST


text:
! Copyright (C) 2015-2017, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

module dvmlm_wrapper
  use kinds
  use model, only: etot
  ! wrapper to dvmlm subroutine
  implicit none
  private

  public :: dvmlm_minimize

contains

  subroutine dvmlm_min(prm,grd,dim,accuracy,ndim,mstep,task,wa)
    integer :: dim
    real(kflt), intent(inout) :: prm(dim)
    real(kflt), intent(inout) :: grd(dim)
    integer,       intent(in) :: accuracy,ndim,mstep
    character(60), intent(inout) :: task
    real(kflt),    intent(in) :: wa(:)
    integer       :: isave(5)
    real(kflt) :: dsave(24)
    real(kflt) :: f
    real(kflt) :: frtol
    real(kflt) :: fatol
    real(kflt) :: fmin
    external dvmlm

    ! set prms for minimization
    ! frtol: desired relative error
    ! fatol: desired absolute error
    select case(accuracy)
    case(0)
       ! low accuracy
       frtol = 1.0e-2_kflt
       fatol = 1.0e-3_kflt
    case(1)
       ! moderate accuracy
       frtol = 1.0e-4_kflt
       fatol = 1.0e-6_kflt
    case(2)
       ! high accuracy
       frtol = 1.0e-8_kflt
       fatol = 1.0e-12_kflt
    end select

    fmin = -1.9e30_kflt

    f = - etot
    call dvmlm(ndim,prm,f,grd,frtol,fatol,fmin,task,mstep,&
         wa(1),wa(ndim*mstep+1),wa(2*ndim*mstep+1),&
         isave,dsave,wa(2*ndim*mstep+mstep+1),wa(2*ndim*mstep+mstep+ndim+1))

  end subroutine dvmlm_min

  subroutine dvmlm_minimize(nv,ns,nd,dim,data_samples,w,prm,grd,accuracy,iter,totiter)
    use model, only: update_gradient
    integer, intent(in) :: nv,ns,nd,dim
    integer, intent(in) :: data_samples(nv,ns)
    real(kflt), intent(in) :: w(nd)
    real(kflt), intent(inout) :: prm(dim)
    real(kflt), intent(inout) :: grd(dim)
    integer, intent(in) :: accuracy
    integer, intent(out) :: iter,totiter
    integer :: err
    integer :: ndim,mstep
    character(60) :: task
    real(kflt), allocatable :: wa(:)
    integer :: lwa

    iter = 0
    totiter = 0
    task = 'START'

    ndim = size(prm)
    ! this is the number of steps for hessian approximation
    mstep = 100
    lwa = 2*ndim*mstep + 2*ndim + mstep
    allocate(wa(lwa),stat=err)

    call update_gradient(nv,ns,nd,data_samples,w,prm(:ns),prm(ns+1:),grd(:ns),grd(ns+1:))
    do
       if(totiter > 100) then
          write(0,*) 'warning: totiter > 100'
          flush(0)
       end if
       call dvmlm_min(prm,grd,size(prm),accuracy,ndim,mstep,task,wa)
       if(task(1:2) == 'FG') then
          ! update etot and gradient for line search
          totiter = totiter + 1
          call update_gradient(nv,ns,nd,data_samples,w,prm(:ns),prm(ns+1:),grd(:ns),grd(ns+1:))
       elseif(task(1:4) == 'NEWX') then
          ! start new line search
          iter = iter + 1
       elseif(task(1:4) == 'WARN') then
          write(0,*) 'warning ', iter
          flush(0)
       elseif(task(1:4) == 'CONV') then
          ! compute final values for likelihood
          call update_gradient(nv,ns,nd,data_samples,w,prm(:ns),prm(ns+1:),grd(:ns),grd(ns+1:))
          exit
       end if
    end do

    deallocate(wa)

  end subroutine dvmlm_minimize

end module dvmlm_wrapper

max_stars_repo_name: simomarsili/fmpl

max_issues_repo_name: simomarsili/fmpl

max_forks_repo_name: simomarsili/fmpl


text:
module Mod_nsm_PorousMedia
   use typre
   use Mod_Mesh
   use Mod_Element
   use Mod_PointerSetter
   use Mod_nsm_BaseElmope
   implicit none
   private
   public SetPointersPorousMedia

   type, extends(PointerSetter) :: SPPorousMedia
contains
      procedure :: SpecificSet => SpecificSetPorousMedia
   end type
   type(SPPorousMedia) :: SetPointersPorousMedia

 contains

   subroutine SpecificSetPorousMedia(d)
      implicit none
      class(SPPorousMedia) :: d

      if (a%kfl_PorousMedia == 1) then
         call ConcatenateProcedures(ProcHook_InGaussElmats,PorousBox)
      endif
   end subroutine

!----------------------------------------------------------------------------
   subroutine PorousBox
      implicit none
      real(rp) :: gpcod(3)
      real(rp) ::  midPoint

      midPoint = a%PME%minX + a%PME%Dis       !Divided to 2 box

      call e%interpg(e%ndime,e%elcod,gpcod)
      if (a%PME%minX <= gpcod(1) .and. gpcod(1) < midPoint  .and. a%PME%minY <= gpcod(2) .and. gpcod(2) <= a%PME%maxY ) then      !Box 1
         call elm_porous(e,dvol,a%PME%perm1,wrmat1)
      elseif(midPoint <= gpcod(1) .and. gpcod(1)<= a%PME%maxX .and. a%PME%minY <= gpcod(2) .and. gpcod(2) <= a%PME%maxY ) then    !Box 2
         call elm_porous(e,dvol,a%PME%perm2,wrmat1)
      endif
   end subroutine

   subroutine elm_porous(e,dvol,perm,myelmat)
         implicit none
         type(FiniteElement)     :: e
         real(rp), intent(in)    :: dvol
         real(rp), intent(inout) :: myelmat(e%mnode,e%mnode)
         integer(ip) :: inode,jnode
         real(rp)    :: tmp, perm,invperm,acvis

         invperm = 1.0_rp/perm
         do inode=1,e%pnode
            do jnode=1,e%pnode
               myelmat(inode,jnode) = e%shape(inode,e%igaus)*e%shape(jnode,e%igaus)*invperm*dvol+ myelmat(inode,jnode)
            end do
         end do
   end subroutine

 end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
SUBROUTINE TIMELOG(N,DAYJUL,OUTDIR1)

  ! *** WRITE OUT MODEL TIME STEP AND SUN/PC SYSTEM CLOCK TIME TO TIME.LOG FILE

  USE GLOBAL,ONLY:IK8
  IMPLICIT NONE
  INTEGER(IK8) :: N
  DOUBLE PRECISION, INTENT(IN) :: DAYJUL
  CHARACTER*8,      INTENT(IN) :: OUTDIR1
  CHARACTER*8 MRMDATE,MRMTIME*10

  ! *** WRITE OUT MODEL TIME STEP AND SYSTEM CLOCK TIME TO TIME.LOG
  CALL DATE_AND_TIME(MRMDATE,MRMTIME)

  OPEN(9,FILE=OUTDIR1//'TIME.LOG',POSITION='APPEND')
  WRITE(9,100)N,DAYJUL,MRMDATE,MRMTIME
  CLOSE(9)

  100 FORMAT(' ','N =',I12,5X,'TIMEDAY = ',F12.4,5X,'DATE = ',A8,5X,'TIME = ',A10)

END SUBROUTINE

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
      module comgrid_sg
!***********************************************************************
!***********************************************************************
!     Module containing grid information variables for subgrid
!     calculations.
!
!     Identifier      Type     Description
!
!     n_grid_points_sg  integer  number of points (cells) in the grid
!     lda_sg            integer  input flag denoting the format of the
!                                stor file: 1 - ASCII, 2 - UNFORMATTED
!     x_sg, y_sg, z_sg  real*8   coordinates of each grid cell (center of
!                                cell), arrays of size n_grid_points_sg
!     sx1_sg            real*8   volume of each cell (m^3), array of size
!                                n_grid_points_sg
!     izonef_sg         integer  array denoting the zone number associated
!                                with each grid cell, array of size
!                                n_grid_points_sg
!***********************************************************************

      integer n_grid_points_sg
      integer lda_sg
      real*8,  allocatable :: x_sg(:)
      real*8,  allocatable :: y_sg(:)
      real*8,  allocatable :: z_sg(:)
      real*8,  allocatable :: sx1_sg(:)
      integer, allocatable :: izonef_sg(:)

      end module comgrid_sg

max_stars_repo_name: lanl/PLUMECALC

max_issues_repo_name: lanl/PLUMECALC

max_forks_repo_name: lanl/PLUMECALC


! https://github.com/MichaelSiehl/Atomic_Subroutines--How_the_Parallel_Codes_may_look_like--Part_1

module Main_Sub
!
contains
!
!**********
!
subroutine Entry_Main_Sub
  !
  use OOOPimma_admImageManager
  !
  implicit none
  !
  call OOOPimma_Start (OOOPimmaImageManager_1) ! start the ImageManager on all images
  !
end subroutine Entry_Main_Sub
!
!**********
!
end module Main_Sub


text:
      SUBROUTINE DRWREC (XL,XH, YL,YH)

C  DRWREC draws a rectangle as specified by x-low,x-high,y-low,y-high

      REAL     XL,XH, YL,YH, XA(5),YA(5)

      XA(1) = XL
      XA(2) = XL
      XA(3) = XH
      XA(4) = XH
      XA(5) = XL

      YA(1) = YL
      YA(2) = YH
      YA(3) = YH
      YA(4) = YL
      YA(5) = YL
      CALL PPL (5,XA,YA)

      END


text:
*
*     -------------------------------------------------------------
*      J F A Z E
*     -------------------------------------------------------------
*                                                                  *
*     DETERMINATE THE PHASE FACTOR WHICH APPEAR FROM PERMUTATION   *
*     OF OPERATORS OF SECOND QUANTIZATION                          *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University,  Nashville            February 1994   *
*
      FUNCTION JFAZE(I1,I2,I3,I4)
      JFAZE=1
      IF(I1.GT.I2)JFAZE=-JFAZE
      IF(I1.GT.I3)JFAZE=-JFAZE
      IF(I1.GT.I4)JFAZE=-JFAZE
      IF(I2.GT.I3)JFAZE=-JFAZE
      IF(I2.GT.I4)JFAZE=-JFAZE
      IF(I3.GT.I4)JFAZE=-JFAZE
      RETURN
      END


text:
	SUBROUTINE GD_WPGDW  ( iacss, grid, igx, igy, ighdr, gdattm1,
     +			       gdattm2, level1, level2, ivcord, parm,
     +			       rewrit, ipktyp, nbits, iret )
C************************************************************************
C* GD_WPGDW								*
C*									*
C* This subroutine packs an input grid of real values and writes it	*
C* to a grid file.  IPKTYP should be one of the following parameter	*
C* names from GEMPRM.PRM:						*
C*									*
C*         MDGNON        No grid packing				*
C*         MDGGRB        Pack in GEMPAK GRIB format given nbits		*
C*         MDGDEC        Pack in GEMPAK GRIB format given precision	*
C*         MDGDIF        Pack in GEMPAK DIF format given nbits		*
C*									*
C* If the packing type is MDGNON, the real data will be stored as if	*
C* GD_WDAT were called.  If MDGGRB or MDGDIF is specified, the		*
C* number of bits given in NBITS will be used to store the data.	*
C* For packing type MDGDEC, NBITS is the precision.  The grid data	*
C* is multiplied by 10 ** NBITS and rounded to the nearest integer.	*
C* The actual number of bits used to store the data is the minimum	*
C* number required to store the resulting integers.			*
C*									*
C* GD_WPGDW  ( IACSS, GRID, IGX, IGY, IGHDR, GDATTM1, GDATTM2, LEVEL1,	*
C*	       LEVEL2, IVCORD,PARM, REWRIT, IPKTYP, NBITS, IRET )	*
C*									*
C* Input parameters:							*
C*	IACSS 		INTEGER		Grid access number		*
C*	GRID (IGX,IGY)	REAL		Grid data			*
C*	IGX		INTEGER		Number of horizontal points	*
C*	IGY		INTEGER		Number of vertical points 	*
C*	IGHDR (IHDRSZ)	INTEGER		Grid header			*
C*	GDATTM1		CHAR*20		GEMPAK times			*
C*	GDATTM2		CHAR*20		GEMPAK times			*
C*	LEVEL1		INTEGER		Vertical levels			*
C*	LEVEL2		INTEGER		Vertical levels			*
C*	IVCORD		INTEGER		Vertical coordinate		*
C*				  	   0 = NONE			*
C*				  	   1 = PRES			*
C*					   2 = THTA			*
C*					   3 = HGHT			*
C*	PARM		CHAR*12		Parameter name			*
C*	REWRIT		LOGICAL		Flag to replace existing grid	*
C*	IPKTYP		INTEGER		Packing type			*
C*	NBITS		INTEGER		Number of bits / precision	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = file not open		*
C*					 -5 = no write access		*
C*					 -6 = read/ write error		*
C*					 -9 = invalid grid size		*
C*					-10 = grid already exists	*
C*					-11 = grid file is full		*
C**									*
C* Log:									*
C* R. Tian/SAIC          3/06						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		grid (*)
	CHARACTER*(*)	gdattm1, gdattm2, parm
	INTEGER		ighdr (*)
	LOGICAL		rewrit
C*
	CHARACTER	gdattm (2)*20
	INTEGER		level (2)
C------------------------------------------------------------------------
	gdattm (1) = gdattm1
	gdattm (2) = gdattm2
	level (1) = level1
	level (2) = level2
	CALL GD_WPGD ( iacss, grid, igx, igy, ighdr, gdattm,
     +		       level, ivcord, parm, rewrit, ipktyp, nbits, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
cc Copyright (C) 2010-2011: Leslie Greengard, Zydrunas Gimbustas
cc and Manas Rachh
cc Contact: greengard@cims.nyu.edu
cc
cc This program is free software; you can redistribute it and/or modify
cc it under the terms of the GNU General Public License as published by
cc the Free Software Foundation; either version 2 of the License, or
cc (at your option) any later version.  This program is distributed in
cc the hope that it will be useful, but WITHOUT ANY WARRANTY; without
cc even the implied warranty of MERCHANTABILITY or FITNESS FOR A
cc PARTICULAR PURPOSE.  See the GNU General Public License for more
cc details. You should have received a copy of the GNU General Public
cc License along with this program;
cc if not, see <http://www.gnu.org/licenses/>.

c
c   Helmholtz FMM in R^2: evaluate all pairwise particle
c   interactions (ignoring self-interaction)
c   and interactions with targets.
c
c   We use H_0(kr)*(i/4) for the Green's function.
c   Self-interactions are not included
c
c   h2d: charge and dipstr are complex valued, x in \R^2
c
c   \phi(x_i) = (i/4)\sum_{j\ne i} charge_j H^{(1)}_0(k |x_i - x_j|)
c   + dipstr_j (dipvec_j \dot (x_i - x_j)) H^{(1)}_1(k |x_i - x_j|*
c                                          k/|x_i-x_j|
c
c

c
c
c
c
c
c

      subroutine hfmm2d_s_c_p(eps,zk,ns,sources,
     1            charge,pot,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge
cf2py  intent(out) pot,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      complex *16 charge(ns)

      complex *16 pot(ns)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 1
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_s_c_g(eps,zk,ns,sources,
     1            charge,pot,grad,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge
cf2py  intent(out) pot,grad,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      complex *16 charge(ns)

      complex *16 pot(ns),grad(2,ns)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      nt = 0

      ifpgh = 2
      ifpghtarg = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_s_c_h(eps,zk,ns,sources,
     1            charge,pot,grad,hess,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge
cf2py  intent(out) pot,grad,hess,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      complex *16 charge(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 3
      ifpghtarg = 0

      nd = 1

      nt = 0

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_s_d_p(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec
cf2py  intent(out) pot,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns)

c
cc     temporary variables
c
      complex *16 charge
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 1
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_s_d_g(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,grad,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec
cf2py  intent(out) pot,grad,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns),grad(2,ns)

c
cc     temporary variables
c
      complex *16 charge
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 2
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_s_d_h(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,grad,hess,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec
cf2py  intent(out) pot,grad,hess,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)

c
cc     temporary variables
c
      complex *16 charge
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 3
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_s_cd_p(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,charge
cf2py  intent(out) pot,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns)

c
cc     temporary variables
c
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 1
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c

      subroutine hfmm2d_s_cd_g(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,grad,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,charge
cf2py  intent(out) pot,grad,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns),grad(2,ns)

c
cc     temporary variables
c
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 2
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_s_cd_h(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,grad,hess,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,charge
cf2py  intent(out) pot,grad,hess,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,ier
      real *8 sources(2,ns)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)

c
cc     temporary variables
c
      integer nt
      real *8 targ(2)
      complex *16 pottarg(1)
      complex *16 gradtarg(2)
      complex *16 hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 3
      ifpghtarg = 0

      nt = 0

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

      subroutine hfmm2d_t_c_p(eps,zk,ns,sources,
     1            charge,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      complex *16 pot(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 0
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_t_c_g(eps,zk,ns,sources,
     1            charge,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pottarg,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 0
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_t_c_h(eps,zk,ns,sources,
     1            charge,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pottarg,gradtarg,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 0
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_t_d_p(eps,zk,ns,sources,
     1            dipstr,dipvec,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 charge
      complex *16 pot(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_t_d_g(eps,zk,ns,sources,
     1            dipstr,dipvec,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 charge
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_t_d_h(eps,zk,ns,sources,
     1            dipstr,dipvec,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,gradtarg,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      complex *16 charge
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_t_cd_p(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 pot(1)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c

      subroutine hfmm2d_t_cd_g(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_t_cd_h(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pottarg,gradtarg,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd
      complex *16 pot(1)
      complex *16 grad(2)
      complex *16 hess(3)

      ifcharge = 1
      ifdipole = 1

      ifpgh = 0
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c
c
c
c
c
c

      subroutine hfmm2d_st_c_p(eps,zk,ns,sources,
     1            charge,pot,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pot,pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pot(ns)
      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 1
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_st_c_g(eps,zk,ns,sources,
     1            charge,pot,grad,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pot(ns),grad(2,ns)
      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 2
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_st_c_h(eps,zk,ns,sources,
     1            charge,pot,grad,hess,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,hess,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      complex *16 charge(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)
      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      complex *16 dipstr
      real *8 dipvec(2)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 0

      ifpgh = 3
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_st_d_p(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns)
      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 charge
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 1
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------

      subroutine hfmm2d_st_d_g(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,grad,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns),grad(2,ns)
      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 charge
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 2
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_st_d_h(eps,zk,ns,sources,
     1            dipstr,dipvec,pot,grad,hess,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,hess,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 dipstr(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)
      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      complex *16 charge
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 0
      ifdipole = 1

      ifpgh = 3
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end

c-------------------------------

      subroutine hfmm2d_st_cd_p(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,nt,targ,pottarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   pottarg(nt)   : potential at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns)
      complex *16 pottarg(nt)

c
cc     temporary variables
c
      complex *16 grad(2),gradtarg(2)
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg

      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 1
      ifpghtarg = 1

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c

      subroutine hfmm2d_st_cd_g(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,grad,nt,targ,pottarg,gradtarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns),grad(2,ns)
      complex *16 pottarg(nt),gradtarg(2,nt)

c
cc     temporary variables
c
      complex *16 hess(3),hesstarg(3)
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 2
      ifpghtarg = 2

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end
c------------------------------
c
c
c
c
c
      subroutine hfmm2d_st_cd_h(eps,zk,ns,sources,charge,
     1            dipstr,dipvec,pot,grad,hess,nt,targ,pottarg,
     2            gradtarg,hesstarg,ier)
cf2py  intent(in) eps
cf2py  intent(in) zk
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,nt,targ
cf2py  intent(out) pot,pottarg,grad,gradtarg,hess,hesstarg,ier

c----------------------------------------------
c   INPUT PARAMETERS:
c   eps           : FMM precision requested
c   zk            : Helmholtz parameter
c   ns            : number of sources
c   sources(2,ns) : source locations
c   charge(ns)    : charge strengths
c   dipstr(ns)    : dipole strengths
c   dipvec(2,ns)  : dipole orientation vectors
c   nt            : number of targets
c   targ(2,nt)    : target locations
c
c   OUTPUT PARAMETERS
c   pot(ns)       : potential at the source locations
c   grad(2,ns)    : gradients at the source locations
c   hess(3,ns)    : hessian at the source locations
c   pottarg(nt)   : potential at the target locations
c   gradtarg(2,nt): gradient at the target locations
c   hesstarg(3,nt): hessian at the target locations
c

      implicit none
c
cc      calling sequence variables
c
      real *8 eps
      complex *16 zk
      integer ns,nt,ier
      real *8 sources(2,ns),targ(2,nt)
      real *8 dipvec(2,ns)
      complex *16 charge(ns),dipstr(ns)

      complex *16 pot(ns),grad(2,ns),hess(3,ns)
      complex *16 pottarg(nt),gradtarg(2,nt),hesstarg(3,nt)

c
cc     temporary variables
c
      integer ifcharge,ifdipole,iper
      integer ifpgh,ifpghtarg
      integer nd

      ifcharge = 1
      ifdipole = 1

      ifpgh = 3
      ifpghtarg = 3

      nd = 1

      call hfmm2d(nd,eps,zk,ns,sources,ifcharge,charge,
     1            ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess,
     2            nt,targ,ifpghtarg,pottarg,gradtarg,
     3            hesstarg,ier)
      return
      end


text:
      SUBROUTINE PARTN2 (CP,RP,CORE,BUF)
C
C     THIS IS AN INITIALIZATION ROUTINE FOR PARTN1 AND MERGE1.
C     IT CALLS PARTN3 TO BUILD THE BIT STRINGS FROM THE PARTITIONING
C     VECTORS -CP- AND -RP- AND SETS DEFAULT OPTIONS FOR -CP- AND  -RP-
C     BASED ON -SYM-.
C
C
      LOGICAL         CPNULL  ,RPNULL  ,CPHERE  ,RPHERE
      INTEGER         CP      ,RP      ,CORE    ,BUF(4)  ,SUBR(2) ,
     1                CPSIZE  ,RPSIZE  ,CPONES  ,RPONES  ,Z       ,
     2                SYM     ,TYPE    ,FORM    ,SYSBUF  ,OUTPT   ,
     3                CPCOL   ,RPCOL
      CHARACTER       UFM*23  ,UWM*25  ,UIM*29  ,SFM*25  ,SWM*27
      COMMON /XMSSG / UFM     ,UWM     ,UIM     ,SFM     ,SWM
      COMMON /SYSTEM/ SYSBUF  ,OUTPT
      COMMON /PRTMRG/ CPSIZE  ,RPSIZE  ,CPONES  ,RPONES  ,CPNULL  ,
     1                RPNULL  ,CPHERE  ,RPHERE  ,ICP     ,NCP     ,
     2                IRP     ,NRP
      COMMON /ZZZZZZ/ Z(1)
      COMMON /BLANK / SYM     ,TYPE    ,FORM(4) ,CPCOL   ,RPCOL   ,
     1                IREQCL
      DATA    SUBR  / 4HPART  ,4HN2    /
C
C
C              I             I             I                I
C       SYM    I  RP-PURGED  I  CP-PURGED  I NEITHER-PURGED I
C     ---------+-------------+-------------+----------------+--------
C              I             I             I                I
C       .LT.0  I  RP IS SET  I  CP IS SET  I  RP MUST HAVE  I
C              I  = TO CP    I  = TO RP    I  SAME ONES-S   I
C              I             I             I  COUNT AS CP   I
C     ---------+-------------+-------------+----------------+--------
C              I             I             I                I
C       .GE.0  I  RP IS SET  I  CP IS SET  I  USE CP AND RP I
C              I  TO ALL 0   I  TO ALL 0   I                I
C              I             I             I                I
C
C     IN ALL CASES, RESULTANT -CP- AND -RP- DIMENSIONS MUST BE
C     COMPATIBLE TO THOSE OF  -A-
C
C
C     CONVERT COLUMN PARTITIONING VECTOR TO BIT STRING.
C
      ICP = 1
      IREQCL = CPCOL
      CALL PARTN3 (CP,CPSIZE,CPONES,ICP,NCP,CPHERE,BUF,CORE)
      CPCOL = IREQCL
      IF (CPHERE) GO TO 10
      IRP = 1
      GO TO 20
   10 IRP = NCP + 1
C
C     CONVERT ROW PARTITIONING VECTOR TO BIT STRING.
C
   20 IREQCL = RPCOL
      CALL PARTN3 (RP,RPSIZE,RPONES,IRP,NRP,RPHERE,BUF,CORE)
      RPCOL = IREQCL
C
C     BRANCH ON SYMMETRIC OR  NON-SYMMETRIC DMAP VARIABLE SYM
C
      CPNULL = .FALSE.
      RPNULL = .FALSE.
      IF (SYM) 30,140,140
C
C     DMAP USER CLAIMS SYMMETRIC INPUT AND OUTPUT
C
   30 IF (CPHERE) GO TO 70
C
C     -CP- IS PURGED.  CHECK FOR -RP- PURGED (ERROR), AND IF NOT SET
C     -CP- BITS EQUAL TO -RP- BITS
C
      IF (RPHERE) GO TO 60
C
C     BOTH -RP- AND -CP- PURGED AND -A- IS NOT PURGED (ERROR).
C
   40 WRITE  (OUTPT,50) SFM
   50 FORMAT (A25,' 2170, BOTH THE ROW AND COLUMN PARTITIONING VECTORS',
     1       ' ARE PURGED AND ONLY ONE MAY BE.')
      CALL MESAGE (-61,0,SUBR)
C
C     SET CP-ONES = RP-ONES BY SIMPLE EQUIVALENCE OF CORE SPACE
C
   60 ICP = IRP
      NCP = NRP
      CPONES = RPONES
      CPSIZE = RPSIZE
      GO TO 170
C
C     -CP- IS NOT PURGED.  IF -RP- IS PURGED IT IS SET EQUAL TO -CP-.
C
   70 IF (RPHERE) GO TO 80
      IRP = ICP
      NRP = NCP
      RPONES = CPONES
      RPSIZE = CPSIZE
      GO TO 170
C
C     BOTH -RP- AND -CP- ARE PRESENT AND SINCE USER HAS SPECIFIED A
C     SYMMETRIC OUTPUT PARTITION IS DESIRED THE NUMBER OF
C     NON-ZEROS IN-CP- MUST EQUAL THE NUMBER OF NON-ZEROS IN -RP- FOR NO
C     ERROR HERE.
C
   80 IF (CPONES.EQ.RPONES .AND. CPSIZE.EQ.RPSIZE) GO TO 100
      WRITE  (OUTPT,90) SWM,CP,RP
   90 FORMAT (A27,' 2171, SYM FLAG INDICATES TO THE PARTITION OR MERGE',
     1       ' MODULE THAT A SYMMETRIC MATRIX IS TO BE', /5X,
     2       ' OUTPUT.  THE PARTITIONING VECTORS',2I4,' HOWEVER DO NOT',
     3       ' CONTAIN AN IDENTICAL NUMBER OF ZEROS AND NON-ZEROS.')
C
C     CHECK FOR ORDER OF ONES IN ROW AND COLUMN PARTITIONING VECTOR.
C
  100 IF (CPSIZE .NE. RPSIZE) GO TO 170
      J = IRP
      DO 110 I = ICP,NCP
      IF (Z(I) .NE. Z(J)) GO TO 120
      J = J + 1
  110 CONTINUE
      GO TO 170
C
C     ROW AND COLUMN PARTITIONING VECTORS DO NOT HAVE SAME ORDER.
C
  120 WRITE  (OUTPT,130) SWM
  130 FORMAT (A27,' 2172, ROW AND COLUMN PARTITIONING VECTORS DO NOT ',
     1       'HAVE IDENTICAL ORDERING OF ZERO', /5X,' AND NON-ZERO ',
     2       'ELEMENTS, AND SYM FLAG INDICATES THAT A SYMMETRIC ',
     3       'PARTITION OR MERGE IS TO BE PERFORMED.')
      GO TO 170
C
C     DMAP USER DOES NOT REQUIRE SYMMETRY
C
  140 IF (CPHERE) GO TO 160
C
C     -CP- IS PURGED.  THUS -RP- MUST BE PRESENT FOR NO ERROR.
C
      IF (RPHERE) GO TO 150
      GO TO 40
C
C     SET CP-ONES EQUAL TO 0 AND CPSIZE = 0
C
  150 CPNULL = .TRUE.
      CPSIZE = 0
      CPONES = 0
      GO TO 170
C
C     -CP- NOT PURGED.  IF -RP- IS PURGED SET IT NULL.
C
  160 IF (RPHERE) GO TO 170
      RPNULL = .TRUE.
      NRP    = IRP - 1
      RPSIZE = 0
      RPONES = 0
  170 RETURN
      END


program main
    implicit none

    integer :: x(1:6) = (/ 3, 4, 5, 6, 7, 8 /)
    integer :: i
    integer, dimension(size(x)) :: res
    res = [(square(x(i)), i = 1, size(x))]

    write(*, *) res

    contains
    function square(x) result(res)
        implicit none
        integer, intent(in) :: x
        integer :: res
        !integer, intent(out) :: res
        res = x * x
    end function square

end program main


    program demo_besj0
    use, intrinsic :: iso_fortran_env, only : real_kinds, &
    & real32, real64, real128
       implicit none
       real(kind=real64) :: x = 0.0_real64
       x = bessel_j0(x)
       write(*,*)x
    end program demo_besj0


text:
      SUBROUTINE SHFDI ( INTL,DI )
C
      CHARACTER DI*5,TUNITS(5)*1
      DIMENSION IZINTL(5)
      INTEGER*4 IZINTL, INTL                                            M
C
      DATA IZINTL /1,60,1440,43200,525600 /
      DATA TUNITS /'N','H','D','M','Y' /
C
      DI='DI'
      N = 1
   5  CONTINUE
      IF ( N.LE.4.AND.INTL.GE.IZINTL(N+1) ) THEN
          N=N+1
          GO TO 5
      ENDIF
      DI(3:3) = TUNITS(N)
      NINT = INTL/IZINTL(N)
      IF ( NINT.LT.10 ) THEN
          WRITE ( DI(4:4),'(I1)') NINT
      ELSE
          WRITE ( DI(4:5),'(I2)') NINT
      ENDIF
      RETURN
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
      subroutine parse_targets(tgt_info,wlist)

      implicit none

      include 'stdunit.h'
      include 'mdef_target_info.h'
      include 'def_word_list.h'
      include 'par_parse_targets.h'

      integer, parameter ::
     &     ntest = 00
      character(len=13), parameter ::
     &     i_am = 'parse_targets'

      type(target_info), intent(inout) ::
     &    tgt_info
      type(word_list) ::
     &    wlist

      character(len=256) ::
     &     word, cur_target, cur_rule
      character ::
     &     sep
      integer ::
     &     mode, error, line, col

      logical, external ::
     &     advance_word_list_entry,
     &     pt_handle_target, pt_handle_depend, pt_handle_rules

      if (ntest.gt.0) call write_title(lulog,wst_dbg_subr,i_am)

      ! safty exit
      if (.not.associated(wlist%head))
     &      call quit(1,i_am,'buggy wlist')

      if (ntest.ge.100) then
        write(lulog,*) 'word list:'
        call print_word_list(lulog,wlist)
      end if

      ! set "current" pointer to head
      wlist%current => wlist%head

      mode = scan_for_targets
      error = no_error
      cur_target = '- '
      cur_rule = '- '
      main_loop: do

        ! read present entry
        call get_word_list_entry(word,sep,wlist)

        if (ntest.ge.100) then
          write(lulog,*) 'word = ',trim(word),'   sep = ',sep
        end if

        if (mode.ne.scan_for_arguments) then
          select case(trim(word))
            ! nothing
            case('')
              ! only newline or close_block is an acceptable separator
              if (sep.ne.'E'.and.sep.ne.')') then
                error = error_unexp_sep
                exit main_loop
              end if
              if (sep.eq.')') then
                if (mode.ne.scan_for_rules) then
                  error = error_unexp_eob
                  exit main_loop
                end if
                ! current target definition has ended
                if (.not.advance_word_list_entry(wlist,'u')) then
                  call quit(1,i_am,'unexpected error for u-branch')
                end if
                mode = scan_for_targets
                cur_target = '- '
              end if
            ! a keyword
            !case('if')
            case('target')
              if (ntest.ge.100) write(lulog,*) 'handle case: target'
              error = pt_handle_target(tgt_info,cur_target,wlist,mode)
              if (error.ne.no_error) exit main_loop
              cycle main_loop
            case('depend')
              if (ntest.ge.100) write(lulog,*) 'handle case: depend'
              error = pt_handle_depend(tgt_info,cur_target,wlist,mode)
              if (error.ne.no_error) exit main_loop
              cycle main_loop
            case default
              ! check for rules
              if (ntest.ge.100) write(lulog,*) 'handle case: rules'
              error = pt_handle_rules(tgt_info,cur_target,wlist,mode)
              if (error.ne.no_error) exit main_loop
              cycle main_loop
          end select

        else

        end if

        if (.not.advance_word_list_entry(wlist,' ')) exit main_loop

      end do main_loop

      if (error.eq.no_error) return

      ! error handling:
      call get_word_list_entry(word,sep,wlist)
      call get_word_list_position(line,col,wlist)
      write(lulog,'(/,x,a,i4,a,i3)')
     &       'Input error near line ',line,' column ',col

      select case(error)
      case(error_unexp_sep)
        write(lulog,*) 'Unexpected separator "'//sep//'"'
      case(error_exp_bob)
        write(lulog,*) 'I expected the beginning of a block: ('
      case(error_unexp_eob)
        write(lulog,*) 'Unexpected end of block'
      case(error_here_no_tgt)
        write(lulog,*) 'You cannot define a target in this context'
      case(error_here_no_dpd)
        write(lulog,*) 'You cannot define a dependency in this context'
      case(error_req_nlsc)
        write(lulog,*) 'I expected one of <newline> or ; '
      case(error_exp_copa)
        write(lulog,*) 'I expected one of , or ) '
      case(error_label_too_long)
        write(lulog,*) 'The label is too long'
      case(error_unknown_comkey)
        write(lulog,*) 'Unknown command or keyword: "'//trim(word)//'"'
      case(error_here_no_rule)
        write(lulog,*) 'You cannot define a rule in this context'
      case(error_unknown_arglab)
        write(lulog,*) 'Unknown argument label: "'//trim(word)//'"'
      case(error_exp_eq)
        write(lulog,*) 'Expected = sign'
      case(error_exp_arg)
        write(lulog,*) 'Expected value for argument here'
      case(error_no_str_arr)
        write(lulog,*) 'A string argument cannot be an array'
      case(error_read_log)
        write(lulog,*) 'Error reading logical value from: "'
     &                  //trim(word)//'"'
      case(error_read_int)
        write(lulog,*) 'Error reading integer value from: "'
     &                  //trim(word)//'"'
      case(error_read_rl8)
        write(lulog,*) 'Error reading real value from: "'
     &                  //trim(word)//'"'
      case(error_arg_not_set)
        write(lulog,*) 'One or more required arguments are missing!'
      case default
        write(lulog,*) 'error code = ',error
      end select

      call quit(0,i_am,'Input error')
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


! { dg-do compile }
! { dg-options "-fopenmp -fdump-tree-omplower" }
  !$omp  parallel
call bar
	!$omp end parallel
  !$omp 	 	p&
!$omp&arallel
call bar
!$omp e&
!$omp&ndparallel
!$omp  &
!$omp  &  &
!$omp pa&
!$omp rallel
call bar
!$omp end parallel
end

! { dg-final { scan-tree-dump-times "pragma omp parallel" 3 "omplower" } }