text:
!=======================================================================
!  SoilKi_init
!
!  Initalize inorganic soil K variables using extractable K from
!     soil file (EXK).
!-----------------------------------------------------------------------
!  REVISION HISTORY
!  10/31/2007 US / CHP Written
!-----------------------------------------------------------------------
!  Called by: SoilKi
!=======================================================================

      SUBROUTINE SoilKi_init (CONTROL, ISWPOT, SOILPROP,  !Input
     &    KFertIndex, Ki_Avail, SKi_Avail,                !Output
     &    SKi_Tot, SKiAvlProf, SKiTotProf)                !Output

!     ------------------------------------------------------------------
      USE ModuleDefs
      IMPLICIT  NONE
      SAVE
!     ------------------------------------------------------------------
      CHARACTER*1  ISWPOT
      CHARACTER*6, PARAMETER :: ERRKEY = 'SKINIT'
      CHARACTER*78 MSG(3)
      TYPE (ControlType) CONTROL

      INTEGER L
      INTEGER NLAYR

!     Soil properties:
      REAL SKiAvlProf, SKiTotProf
      REAL, DIMENSION(NL) :: Ki_Avail, Ki_Tot, SKi_Avail, SKi_Tot
      REAL, DIMENSION(NL) :: KFertIndex, KExchgIndex, KAvailIndex,KG2PPM
      TYPE (SoilType)    SOILPROP

!     ------------------------------------------------------------------
      IF (ISWPOT .NE. 'Y') RETURN

      NLAYR    = SOILPROP % NLAYR
      KG2PPM   = SOILPROP % KG2PPM

!     Initialize total soil K only once
      IF (CONTROL % RUN == 1 .OR. INDEX('QF',CONTROL % RNMODE) <= 0)THEN
        DO L = 1, NLAYR
          IF (SOILPROP % EXK(L) < 0.0) THEN
            MSG(1) = "Missing value found in soil file for potassium."
            MSG(2) = "Cannot run potassium model for this data."
            MSG(3) = "Model will stop."
            CALL WARNING(3,ERRKEY,MSG)
            CALL ERROR(ERRKEY,0," ",0)
          ENDIF
          Ki_Tot(L) = SOILPROP % EXK(L) * 391.
        ENDDO
      ENDIF

!     Update availability indices every season. May be important for
!       long-term sequence simulations.
      SKiAvlProf = 0.
      SKiTotProf = 0.

      DO L = 1, NLAYR
        IF (SOILPROP % OC(L) > 1.E-5) THEN
           KFertIndex(L) = 1.05 + 0.5 * LOG(SOILPROP % OC(L))
        ELSE
           KFertIndex(L) = 0.05   !Value of function at 1.E-6
        ENDIF

        IF (SOILPROP % CEC(L) > 1.E-5) THEN
         KExchgIndex(L) = MIN(1.0, 17.5* SOILPROP % EXK(L) /
     &                      SOILPROP % CEC(L))
        ELSE
          KExchgIndex(L) = 1.00
        ENDIF

        KAvailIndex(L) = MIN(KFertIndex(L), KExchgIndex(L), 1.0)

        Ki_Avail(L) = Ki_Tot(L) * KAvailIndex(L)
        SKi_Avail(L) = Ki_Avail(L) / KG2PPM(L)
        SKi_Tot(L)   = Ki_Tot(L)   / KG2PPM(L)
        SKiAvlProf = SKiAvlProf + SKi_Avail(L)
        SKiTotProf = SKiTotProf + SKi_Tot(L)
      ENDDO

!     ------------------------------------------------------------------
      RETURN
      END SUBROUTINE SoilKi_init

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

ext: for

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


!
!  Copyright 2020 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!=======================================================================
module poisson_isolated
  implicit none

! NOTE: Conjugate-Gradient (CG) method has a reproducibility problem
!       in SALMON calculation...
!       This problem is caused by `omp reduction` clause, which uses
!       atomic instructions as popular implementation.
!       We assure the computation order when merging results between
!       thread.

contains

!============================ Hartree potential (Solve Poisson equation)
subroutine poisson_isolated_cg(lg,mg,info,system,poisson,trho,tVh,srg_scalar,stencil)
  use inputoutput, only: threshold_cg
  use structures, only: s_rgrid,s_parallel_info,s_dft_system,s_poisson,s_sendrecv_grid,s_stencil
  use communication, only: comm_is_root, comm_summation
  use math_constants, only : pi
  use sendrecv_grid, only: update_overlap_real8
  use parallelization, only: get_nthreads, get_thread_id
  implicit none
  integer,parameter :: ndh=4
  type(s_rgrid),intent(in) :: lg
  type(s_rgrid),intent(in) :: mg
  type(s_parallel_info),intent(in) :: info
  type(s_dft_system),intent(in) :: system
  type(s_poisson),intent(inout) :: poisson
  real(8) :: trho(mg%is(1):mg%ie(1),    &
                  mg%is(2):mg%ie(2),      &
                  mg%is(3):mg%ie(3))
  real(8) :: tVh(mg%is(1):mg%ie(1),    &
                 mg%is(2):mg%ie(2),      &
                 mg%is(3):mg%ie(3))
  type(s_sendrecv_grid),intent(inout) :: srg_scalar
  type(s_stencil),intent(in) :: stencil
  integer,parameter :: maxiter=1000
  integer :: ix,iy,iz,iter
  real(8) :: sum1,sum2,ak,ck
  real(8) :: tottmp
  real(8) :: totbox
  real(8) :: rlap_wk(mg%is_array(1):mg%ie_array(1),    &
                     mg%is_array(2):mg%ie_array(2),      &
                     mg%is_array(3):mg%ie_array(3))
  real(8) :: zk(mg%is_array(1):mg%ie_array(1),   &
                mg%is_array(2):mg%ie_array(2),   &
                mg%is_array(3):mg%ie_array(3))
  real(8) :: pk(mg%is_array(1):mg%ie_array(1),   &
                mg%is_array(2):mg%ie_array(2),   &
                mg%is_array(3):mg%ie_array(3))

  integer :: tid
  real(8),allocatable :: reduce_work(:)
  allocate(reduce_work(0:get_nthreads()-1))

  call poisson_boundary(lg,mg,info,system,poisson,trho,pk)

!------------------------- C-G minimization

!$omp parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is_array(3),mg%ie_array(3)
  do iy=mg%is_array(2),mg%ie_array(2)
  do ix=mg%is_array(1),mg%ie_array(1)
    zk(ix,iy,iz)=0.d0
  end do
  end do
  end do

!$omp parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    pk(ix,iy,iz)=tVh(ix,iy,iz)
    zk(ix,iy,iz)=-4.d0*pi*trho(ix,iy,iz)
  end do
  end do
  end do
  if(info%if_divide_rspace) call update_overlap_real8(srg_scalar, mg, pk)
  call laplacian_poisson(mg,pk,rlap_wk,stencil%coef_lap0,stencil%coef_lap)

!$omp parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is_array(3),mg%ie_array(3)
  do iy=mg%is_array(2),mg%ie_array(2)
  do ix=mg%is_array(1),mg%ie_array(1)
    pk(ix,iy,iz)=0.d0
  end do
  end do
  end do

!$omp parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    zk(ix,iy,iz)=zk(ix,iy,iz)-rlap_wk(ix,iy,iz)
    pk(ix,iy,iz)=zk(ix,iy,iz)
  end do
  end do
  end do

  reduce_work=0d0
!$omp parallel private(iz,iy,ix,tid)
  tid=get_thread_id()
!$omp do collapse(2)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    reduce_work(tid)=reduce_work(tid)+zk(ix,iy,iz)**2*system%hvol
  end do
  end do
  end do
!$omp end do
!$omp end parallel
  sum1=sum(reduce_work)

  if(info%isize_r==1)then
  else
    call comm_summation(sum1,sum2,info%icomm_r)
    sum1=sum2
  end if

  iteration : do iter=1,maxiter

    if(info%if_divide_rspace) call update_overlap_real8(srg_scalar, mg, pk)
    call laplacian_poisson(mg,pk,rlap_wk,stencil%coef_lap0,stencil%coef_lap)

    reduce_work=0d0
!$omp parallel private(iz,iy,ix,tid)
    tid=get_thread_id()
!$omp do collapse(2)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      reduce_work(tid)=reduce_work(tid)+(zk(ix,iy,iz)*rlap_wk(ix,iy,iz))
    end do
    end do
    end do
!$omp end do
!$omp end parallel
    totbox=sum(reduce_work)

    if(info%isize_r==1)then
      tottmp=totbox
    else
      call comm_summation(totbox,tottmp,info%icomm_r)
    end if

    ak=sum1/tottmp/system%hvol

!$omp parallel do private(iz,iy,ix) firstprivate(ak) collapse(2)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
       tVh(ix,iy,iz)=tVh(ix,iy,iz)+ak*pk(ix,iy,iz)
       zk(ix,iy,iz)=zk(ix,iy,iz)-ak*rlap_wk(ix,iy,iz)
    end do
    end do
    end do

    reduce_work=0d0
!$omp parallel private(iz,iy,ix,tid)
    tid=get_thread_id()
!$omp do collapse(2)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      reduce_work(tid)=reduce_work(tid)+zk(ix,iy,iz)**2
    end do
    end do
    end do
!$omp end do
!$omp end parallel
    totbox=sum(reduce_work)

    if(info%isize_r==1)then
      tottmp=totbox
    else
      call comm_summation(totbox,tottmp,info%icomm_r)
    end if

    sum2=tottmp*system%hvol

    if ( abs(sum2) < threshold_cg*dble(lg%num(1)*lg%num(2)*lg%num(3)) ) exit

    ck=sum2/sum1 ; sum1=sum2

!$omp parallel do private(iz,iy,ix) firstprivate(ck) collapse(2)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      pk(ix,iy,iz)=zk(ix,iy,iz)+ck*pk(ix,iy,iz)
    end do
    end do
    end do

  end do iteration

  poisson%iterVh=iter
  if ( poisson%iterVh>maxiter .and. comm_is_root(info%id_r)) then
     write(*,*) "Warning:Vh iteration is not converged"
     write(*,'("||tVh(i)-tVh(i-1)||**2/(# of grids) = ",e15.8)') &
                                sum2/dble(lg%num(1)*lg%num(2)*lg%num(3))
  end if

  return

end subroutine poisson_isolated_cg

subroutine laplacian_poisson(mg,pk,rlap_wk,lap0,lapt)
  use structures, only: s_rgrid
  implicit none
  type(s_rgrid),intent(in) :: mg
  real(8),intent(in) :: pk(mg%is_array(1):mg%ie_array(1),  &
                           mg%is_array(2):mg%ie_array(2),  &
                           mg%is_array(3):mg%ie_array(3))
  real(8),intent(out) :: rlap_wk(mg%is_array(1):mg%ie_array(1),  &
                                 mg%is_array(2):mg%ie_array(2),  &
                                 mg%is_array(3):mg%ie_array(3))
  real(8),intent(in)  :: lap0,lapt(4,3)
  integer :: ix,iy,iz

!$omp parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    rlap_wk(ix,iy,iz)=-2.d0*lap0*pk(ix,iy,iz)+(  &
                      lapt(1,1)*(pk(ix+1,iy,iz) + pk(ix-1,iy,iz)) &
                     +lapt(2,1)*(pk(ix+2,iy,iz) + pk(ix-2,iy,iz)) &
                     +lapt(3,1)*(pk(ix+3,iy,iz) + pk(ix-3,iy,iz)) &
                     +lapt(4,1)*(pk(ix+4,iy,iz) + pk(ix-4,iy,iz)) &
                     +lapt(1,2)*(pk(ix,iy+1,iz) + pk(ix,iy-1,iz)) &
                     +lapt(2,2)*(pk(ix,iy+2,iz) + pk(ix,iy-2,iz)) &
                     +lapt(3,2)*(pk(ix,iy+3,iz) + pk(ix,iy-3,iz)) &
                     +lapt(4,2)*(pk(ix,iy+4,iz) + pk(ix,iy-4,iz)) &
                     +lapt(1,3)*(pk(ix,iy,iz+1) + pk(ix,iy,iz-1)) &
                     +lapt(2,3)*(pk(ix,iy,iz+2) + pk(ix,iy,iz-2)) &
                     +lapt(3,3)*(pk(ix,iy,iz+3) + pk(ix,iy,iz-3)) &
                     +lapt(4,3)*(pk(ix,iy,iz+4) + pk(ix,iy,iz-4)))
  end do
  end do
  end do

  return

end subroutine laplacian_poisson

subroutine poisson_boundary(lg,mg,info,system,poisson,trho,wk2)
  use inputoutput, only: natom,lmax_multipole,layout_multipole,natom
  use structures, only: s_rgrid,s_parallel_info,s_dft_system,s_poisson
  use salmon_math, only: ylm
  use communication, only: comm_summation
  use parallelization, only: get_nthreads, get_thread_id

  use omp_lib, only: omp_get_num_threads, omp_get_thread_num, omp_get_max_threads
  use misc_routines, only: ceiling_pow2

  implicit none
  integer,parameter :: ndh=4
  type(s_rgrid),intent(in) :: lg
  type(s_rgrid),intent(in) :: mg
  type(s_parallel_info),intent(in) :: info
  type(s_dft_system),intent(in) :: system
  type(s_poisson),intent(inout) :: poisson
  real(8) :: trho(mg%is(1):mg%ie(1),    &
                 mg%is(2):mg%ie(2),      &
                 mg%is(3):mg%ie(3))
  real(8) :: wk2(mg%is_array(1):mg%ie_array(1),    &
                 mg%is_array(2):mg%ie_array(2),      &
                 mg%is_array(3):mg%ie_array(3))
  integer,parameter :: maxiter=1000
  integer :: ii,jj,ix,iy,iz,lm,ll,m,icen  !,kk,pl,cl
  integer :: ixbox,iybox,izbox
  integer :: j,k
  integer :: istart(0:info%isize_r-1),iend(0:info%isize_r-1)
  integer :: icount
  integer,allocatable :: itrho(:)
  integer :: num_center
  real(8) :: ylm2(25)
  integer :: l2(25)
  real(8) :: xx,yy,zz,rr,sum1,xxxx,yyyy,zzzz,rrrr
  real(8),allocatable :: rholm(:,:),rholm2(:,:) !,rholm3(:,:)
 !integer :: tid
  real(8) :: center_trho2(3)
  real(8),allocatable :: center_trho(:,:)
  real(8),allocatable :: center_trho_nume_deno(:,:)
  real(8),allocatable :: center_trho_nume_deno2(:,:)
  real(8) :: xp2,yp2,zp2,xy,yz,xz
  real(8) :: deno(25)
  real(8) :: rinv, rbox, hvol
  real(8),allocatable :: Rion2(:,:)
  integer,allocatable :: ig_num(:)
  integer,allocatable :: ig(:,:,:)
  real(8),allocatable :: coordinate(:,:)
 !integer :: lmax_multipole_tmp !iwata
  integer :: comm_xyz(3),nproc_xyz(3),myrank_xyz(3)

  integer :: tid
  real(8),allocatable :: reduce_work(:)
  real(8),allocatable :: reduce_work4(:,:)
  allocate(reduce_work(0:get_nthreads()-1))
  allocate(reduce_work4(4,0:get_nthreads()-1))

  comm_xyz(1) = info%icomm_x
  comm_xyz(2) = info%icomm_y
  comm_xyz(3) = info%icomm_z
  nproc_xyz(1) = info%isize_x
  nproc_xyz(2) = info%isize_y
  nproc_xyz(3) = info%isize_z
  myrank_xyz(1) = info%id_x
  myrank_xyz(2) = info%id_y
  myrank_xyz(3) = info%id_z

  !------------------------- Boundary condition (multipole expansion)

  hvol=system%hvol
  if(allocated(poisson%ig_num)) then
    if(.not.allocated(ig_num)) then
      if(layout_multipole==2)then
        allocate(ig_num(natom))
      else if(layout_multipole==3)then
        allocate(ig_num(poisson%npole_partial))
      end if
    end if
    ig_num=poisson%ig_num
  end if
  if(layout_multipole/=1) then
    if(.not.allocated(ig))then
      allocate(ig(3,maxval(poisson%ig_num(:)),poisson%npole_partial))
    end if
    ig=poisson%ig
  end if

  if(.not.allocated(coordinate))then
    allocate(coordinate(minval(lg%is_overlap(1:3)):maxval(lg%ie_overlap(1:3)),3))
  end if
  do j=1,3
    coordinate(lg%is_overlap(j):lg%ie_overlap(j),j)=lg%coordinate(lg%is_overlap(j):lg%ie_overlap(j),j)
  end do

  if(.not.allocated(poisson%wkbound))then
    allocate(poisson%wkbound(lg%num(1)*lg%num(2)*lg%num(3)/minval(lg%num(1:3))*6*ndh))
    allocate(poisson%wkbound2(lg%num(1)*lg%num(2)*lg%num(3)/minval(lg%num(1:3))*6*ndh))
  end if

  select case( layout_multipole )

  case(1)

  num_center=1
  allocate (rholm((lmax_multipole+1)**2,1))
  allocate (rholm2((lmax_multipole+1)**2,1))
  allocate(itrho(1))
  allocate(center_trho(3,1))
  do jj=1,3
    center_trho(jj,1)=0.d0
    center_trho2(jj)=0.d0
  end do
  do lm=1,(lmax_multipole+1)**2
    rholm(lm,1)=0.d0
    rholm2(lm,1)=0.d0
  end do
  itrho(1)=1

  do ll=0,lmax_multipole
  do m=-ll,ll
    lm=ll*ll+ll+1+m
    reduce_work=0d0
!$omp parallel private(ix,iy,iz,xx,yy,zz,rr,xxxx,yyyy,zzzz,tid)
    tid=get_thread_id()
!$omp do collapse(2)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      xx=coordinate(ix,1)-center_trho(1,1)
      yy=coordinate(iy,2)-center_trho(2,1)
      zz=coordinate(iz,3)-center_trho(3,1)
      rr=sqrt(xx*xx+yy*yy+zz*zz)+1.d-50 ; xxxx=xx/rr ; yyyy=yy/rr ; zzzz=zz/rr
      reduce_work(tid)=reduce_work(tid)+rr**ll*ylm(xxxx,yyyy,zzzz,ll,m)*trho(ix,iy,iz)*hvol
    end do
    end do
    end do
!$omp end do
!$omp end parallel
    rholm2(lm,1)=sum(reduce_work)
  end do
  end do

  case(2)

  num_center=natom
  allocate (rholm((lmax_multipole+1)**2,natom))
  allocate (rholm2((lmax_multipole+1)**2,natom))
  allocate(itrho(natom))
  allocate(center_trho(3,natom))
  allocate(Rion2(3,natom))
  Rion2(:,:)=system%Rion(:,:)

  !$OMP parallel do private(icen, lm, jj)
  do icen=1,num_center
    do lm=1,(lmax_multipole+1)**2
      rholm(lm,icen)=0.d0
      rholm2(lm,icen)=0.d0
    end do
    do jj=1,3
      center_trho(jj,icen)=Rion2(jj,icen)
    end do
    itrho(icen)=1
  end do

  do icen=1,poisson%npole_partial
    do ll=0,lmax_multipole
    do m=-ll,ll
      lm=ll*ll+ll+1+m
      reduce_work=0d0
!$omp parallel private(jj,ix,iy,iz,xx,yy,zz,rr,xxxx,yyyy,zzzz,tid)
      tid=get_thread_id()
!$omp do
      do jj=1,ig_num(icen)
        ix=ig(1,jj,icen)
        iy=ig(2,jj,icen)
        iz=ig(3,jj,icen)
        xx=coordinate(ix,1)-Rion2(1,poisson%ipole_tbl(icen))
        yy=coordinate(iy,2)-Rion2(2,poisson%ipole_tbl(icen))
        zz=coordinate(iz,3)-Rion2(3,poisson%ipole_tbl(icen))
        rr=sqrt(xx*xx+yy*yy+zz*zz)+1.d-50 ; xxxx=xx/rr ; yyyy=yy/rr ; zzzz=zz/rr
        reduce_work(tid)=reduce_work(tid)+rr**ll*ylm(xxxx,yyyy,zzzz,ll,m)*trho(ix,iy,iz)*hvol
      end do
!$omp end do
!$omp end parallel
      rholm2(lm,poisson%ipole_tbl(icen))=sum(reduce_work)
    end do
    end do
  end do

  case(3)

  num_center=poisson%npole_total

  allocate (rholm((lmax_multipole+1)**2,num_center))
  allocate (rholm2((lmax_multipole+1)**2,num_center))
  allocate(itrho(num_center))
  allocate(center_trho(3,num_center))
  allocate(center_trho_nume_deno(4,num_center))
  allocate(center_trho_nume_deno2(4,num_center))

  !$OMP parallel do private(icen, jj, lm)
  do icen=1,num_center
    do jj=1,4
      center_trho_nume_deno2(jj,icen)=0.d0
    end do
    do lm=1,(lmax_multipole+1)**2
      rholm(lm,icen)=0.d0
      rholm2(lm,icen)=0.d0
    end do
  end do

  do ii=1,poisson%npole_partial
    reduce_work4=0d0
!$omp parallel private(jj,ixbox,iybox,izbox,xx,yy,zz,tid)
    tid=get_thread_id()
!$omp do
    do jj=1,ig_num(ii)
      ixbox=ig(1,jj,ii)
      iybox=ig(2,jj,ii)
      izbox=ig(3,jj,ii)
      xx=coordinate(ixbox,1)
      yy=coordinate(iybox,2)
      zz=coordinate(izbox,3)
      tid=get_thread_id()
      reduce_work4(1,tid)=reduce_work4(1,tid)+trho(ixbox,iybox,izbox)*xx
      reduce_work4(2,tid)=reduce_work4(2,tid)+trho(ixbox,iybox,izbox)*yy
      reduce_work4(3,tid)=reduce_work4(3,tid)+trho(ixbox,iybox,izbox)*zz
      reduce_work4(4,tid)=reduce_work4(4,tid)+trho(ixbox,iybox,izbox)
    end do
!$omp end do
!$omp end parallel
    center_trho_nume_deno2(1,poisson%ipole_tbl(ii))=sum(reduce_work4(1,:))
    center_trho_nume_deno2(2,poisson%ipole_tbl(ii))=sum(reduce_work4(2,:))
    center_trho_nume_deno2(3,poisson%ipole_tbl(ii))=sum(reduce_work4(3,:))
    center_trho_nume_deno2(4,poisson%ipole_tbl(ii))=sum(reduce_work4(4,:))
  end do

  call comm_summation(center_trho_nume_deno2,center_trho_nume_deno,4*poisson%npole_total,info%icomm_r)

  do ii=1,poisson%npole_total
    if(center_trho_nume_deno(4,ii)*hvol>=1.d-12)then
      itrho(ii)=1
      center_trho(1:3,ii)=center_trho_nume_deno(1:3,ii)/center_trho_nume_deno(4,ii)
    else
      itrho(ii)=0
      center_trho(1:3,ii)=0.d0
    end if
  end do

!  if(omp_get_max_threads() > 16) then
!  !$omp parallel shared(rholm3,lmax_multipole)
!  !$omp master
!      allocate(rholm3((lmax_multipole+1)**2,0:ceiling_pow2(omp_get_num_threads())-1))
!  !$omp end master
!  !$omp end parallel
!
!    lmax_multipole_tmp=lmax_multipole !iwata
!    rholm2=0.d0
!    rholm3=0.d0
!    do ii=1,poisson%npole_partial
!      pl=poisson%ipole_tbl(ii)
!      cl=ig_num(ii)
!      if(itrho(pl)==1)then
!  !$omp parallel default(none) &
!  !$omp          shared(ig,coordinate,center_trho,trho,rholm3) &
!  !$omp          private(tid,kk,jj,ll,lm,ixbox,iybox,izbox,xx,yy,zz,rr,rinv,xxxx,yyyy,zzzz) &
!  !$omp          firstprivate(ii,pl,cl,lmax_multipole_tmp,hvol) !iwata
!        tid=omp_get_thread_num()
!        rholm3(:,tid)=0.d0
!
!  !$omp do
!        do jj=1,cl
!          ixbox=ig(1,jj,ii)
!          iybox=ig(2,jj,ii)
!          izbox=ig(3,jj,ii)
!          xx=coordinate(ixbox,1)-center_trho(1,pl)
!          yy=coordinate(iybox,2)-center_trho(2,pl)
!          zz=coordinate(izbox,3)-center_trho(3,pl)
!          rr=sqrt(xx*xx+yy*yy+zz*zz)+1.d0-50.d0
!          rinv=1.0d0/rr
!          xxxx=xx*rinv
!          yyyy=yy*rinv
!          zzzz=zz*rinv
!          do ll=0,lmax_multipole_tmp !iwata
!          do m=-ll,ll
!            lm=ll*ll+ll+1+m
!            rholm3(lm,tid)=rholm3(lm,tid)+rr**ll*ylm(xxxx,yyyy,zzzz,ll,m)*trho(ixbox,iybox,izbox)*hvol
!          end do
!          end do
!        end do
!  !$omp end do
!
!        kk = ceiling_pow2(omp_get_num_threads())/2
!        do while(kk > 0)
!          if(tid < kk) then
!            rholm3(:,tid) = rholm3(:,tid) + rholm3(:,tid+kk)
!          end if
!          kk = kk/2
!  !$omp barrier
!        end do
!  !$omp end parallel
!      end if
!      rholm2(:,pl)=rholm3(:,0)
!    end do
!    deallocate(rholm3)
!  else
    rholm2=0.d0
    do ii=1,poisson%npole_partial
      if(itrho(poisson%ipole_tbl(ii))==1)then
        rholm=0.d0
        do ll=0,lmax_multipole
        do m=-ll,ll
          lm=ll*ll+ll+1+m
          reduce_work=0d0
!$omp parallel private(jj,ixbox,iybox,izbox,xx,yy,zz,rr,xxxx,yyyy,zzzz,tid)
          tid=get_thread_id()
!$omp do
          do jj=1,ig_num(ii)
            ixbox=ig(1,jj,ii)
            iybox=ig(2,jj,ii)
            izbox=ig(3,jj,ii)
            xx=coordinate(ixbox,1)-center_trho(1,poisson%ipole_tbl(ii))
            yy=coordinate(iybox,2)-center_trho(2,poisson%ipole_tbl(ii))
            zz=coordinate(izbox,3)-center_trho(3,poisson%ipole_tbl(ii))
            rr=sqrt(xx*xx+yy*yy+zz*zz)+1.d-50 ; xxxx=xx/rr ; yyyy=yy/rr ; zzzz=zz/rr
            reduce_work(tid)=reduce_work(tid)+rr**ll*ylm(xxxx,yyyy,zzzz,ll,m)*trho(ixbox,iybox,izbox)*hvol
          end do
!$omp end do
!$omp end parallel
          rholm2(lm,poisson%ipole_tbl(ii))=sum(reduce_work)
        end do
        end do
      end if
    end do
!  endif

  deallocate(center_trho_nume_deno)
  deallocate(center_trho_nume_deno2)

  end select

  if(info%isize_r==1)then
  !$OMP parallel do
    do icen=1,num_center
      rholm(:,icen)=rholm2(:,icen)
    end do
  else
    call comm_summation(rholm2,rholm,(lmax_multipole+1)**2*num_center,info%icomm_r)
  end if

  !$OMP parallel do private(iz,iy,ix) collapse(2)
  do iz=mg%is_array(3),mg%ie_array(3)
  do iy=mg%is_array(2),mg%ie_array(2)
  do ix=mg%is_array(1),mg%ie_array(1)
    wk2(ix,iy,iz)=0.d0
  end do
  end do
  end do

  do k=1,3

  !$OMP parallel do
    do jj=1,lg%num(1)*lg%num(2)*lg%num(3)/minval(lg%num(1:3))*6*ndh
      poisson%wkbound2(jj)=0.d0
    end do

    icount=mg%num(1)*mg%num(2)*mg%num(3)/mg%num(k)*2*ndh

  !$OMP parallel do
    do ii=0,nproc_xyz(k)-1
      istart(ii)=ii*icount/nproc_xyz(k)+1
      iend(ii)=(ii+1)*icount/nproc_xyz(k)
    end do

    do ll=0,lmax_multipole
      do lm=ll**2+1,(ll+1)**2
        l2(lm)=ll
      end do
    end do

  !$OMP parallel do &
  !$OMP private(xx,yy,zz,rr,xxxx,yyyy,zzzz,lm,ll,sum1,ylm2,rrrr,xp2,yp2,zp2,xy,yz,xz,rinv,rbox,deno,icen)
    do jj=istart(myrank_xyz(k)),iend(myrank_xyz(k))
      do icen=1,num_center
        if(itrho(icen)==1)then
          xx=coordinate(poisson%ig_bound(1,jj,k),1)-center_trho(1,icen)
          yy=coordinate(poisson%ig_bound(2,jj,k),2)-center_trho(2,icen)
          zz=coordinate(poisson%ig_bound(3,jj,k),3)-center_trho(3,icen)
          rr=sqrt(xx**2+yy**2+zz**2)+1.d-50
          rinv=1.d0/rr
  !        xxxx=xx/rr ; yyyy=yy/rr ; zzzz=zz/rr
          xx=xx*rinv ; yy=yy*rinv ; zz=zz*rinv

          xp2=xx**2
          yp2=yy**2
          zp2=zz**2
          xy=xx*yy
          yz=yy*zz
          xz=xx*zz

          rrrr=(xp2+yp2+zp2)**2

          ylm2(1)=1.d0
          ylm2(2)=yy
          ylm2(3)=zz
          ylm2(4)=xx
          ylm2(5)=sqrt(3.d0)*xy                      ! lm=5  (2 -2)
          ylm2(6)=sqrt(3.d0)*yz                      ! lm=6  (2 -1)
          ylm2(7)=(2*zp2-xp2-yp2)/2.d0               ! lm=7  (2 0)
          ylm2(8)=sqrt(3.d0)*xz                      ! lm=8  (2 1)
          ylm2(9)=sqrt(3.d0/4.d0)*(xp2-yp2)          ! lm=9  (2 2)
          ylm2(10)=sqrt(5.d0/8.d0)*yy*(3*xp2-yp2)    ! lm=10 (3 -3)
          ylm2(11)=sqrt(15.d0)*xx*yy*zz               ! lm=11 (3 -2)
          ylm2(12)=sqrt(3.d0/8.d0)*yy*(4*zp2-xp2-yp2)  ! lm=12 (3 -1)
          ylm2(13)=zz*(2*zp2-3*xp2-3*yp2)/2.d0         ! lm=13 (3 0)
          ylm2(14)=sqrt(3.d0/8.d0)*xx*(4*zp2-xp2-yp2)  ! lm=14 (3 1)
          ylm2(15)=sqrt(15.d0/4.d0)*zz*(xp2-yp2)       ! lm=15 (3 2)
          ylm2(16)=sqrt(5.d0/8.d0)*xx*(xp2-3*yp2)      ! lm=16 (3 3)

          ylm2(17)=rrrr*sqrt(35.d0)/2.d0*xy*(xp2-yp2)
          ylm2(18)=rrrr*sqrt(35.d0/8.d0)*yz*(3*xp2-yp2)
          ylm2(19)=rrrr*sqrt(5.d0)/2.d0*xy*(7*zp2-1.d0)
          ylm2(20)=rrrr*sqrt(5.d0/8.d0)*yz*(7*zp2-3.d0)
          ylm2(21)=rrrr*(35*zp2**2-30*zp2+3.d0)/8.d0
          ylm2(22)=rrrr*sqrt(5.d0/8.d0)*xz*(7*zp2-3.d0)
          ylm2(23)=rrrr*sqrt(5.d0)/4.d0*(7*zp2-1)*(xp2-yp2)
          ylm2(24)=rrrr*sqrt(35.d0/8.d0)*xz*(xp2-3*yp2)
          ylm2(25)=rrrr*sqrt(35.d0)/8.d0*(xp2**2+yp2**2-6*xp2*yp2)

          deno(1)=rinv
          rbox=rinv*rinv
          deno(2:4)=rbox
          rbox=rbox*rinv
          deno(5:9)=rbox
          rbox=rbox*rinv
          deno(10:16)=rbox
          rbox=rbox*rinv
          deno(17:25)=rbox

          sum1=0.d0
          do lm=1,(lmax_multipole+1)**2
            sum1=sum1+ylm2(lm)*deno(lm)*rholm(lm,icen)
          end do
          poisson%wkbound2(jj) = poisson%wkbound2(jj) + sum1
        end if
      end do
    end do

    if(info%isize_r==1)then
  !$OMP parallel do
      do jj=1,icount
        poisson%wkbound(jj)=poisson%wkbound2(jj)
      end do
    else
      call comm_summation( &
        poisson%wkbound2,              poisson%wkbound,              icount/2, comm_xyz(k), 0                    )
      call comm_summation( &
        poisson%wkbound2(icount/2+1:), poisson%wkbound(icount/2+1:), icount/2, comm_xyz(k), nproc_xyz(k)-1)
    end if

    if(myrank_xyz(k)==0) then
  !$OMP parallel do
      do jj=1,icount/2
        wk2(poisson%ig_bound(1,jj,k),poisson%ig_bound(2,jj,k),poisson%ig_bound(3,jj,k))=poisson%wkbound(jj)
      end do
    end if
    if(myrank_xyz(k)==nproc_xyz(k)-1) then
  !$OMP parallel do
      do jj=icount/2+1,icount
        wk2(poisson%ig_bound(1,jj,k),poisson%ig_bound(2,jj,k),poisson%ig_bound(3,jj,k))=poisson%wkbound(jj)
      end do
    end if

  end do

  deallocate(rholm,rholm2)
  deallocate(center_trho)

  return

end subroutine poisson_boundary

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

subroutine poisson_isolated_ft(lg,mg,info,fg,rho,Vh,poisson)
  use structures
  use communication, only: comm_summation, comm_bcast
  implicit none
  type(s_rgrid)          ,intent(in) :: lg
  type(s_rgrid)          ,intent(in) :: mg
  type(s_parallel_info)  ,intent(in) :: info
  type(s_reciprocal_grid),intent(in) :: fg
  type(s_scalar)         ,intent(in) :: rho
  type(s_scalar)                     :: Vh
  type(s_poisson)                    :: poisson
  !
  integer :: ix,iy,iz
  integer :: kx,ky,kz
  integer :: iix,iiy,iiz

  integer :: imgx_s,imgx_e
  integer :: imgy_s,imgy_e
  integer :: imgz_s,imgz_e
  integer :: ifgx_s,ifgx_e
  integer :: ifgy_s,ifgy_e
  integer :: ifgz_s,ifgz_e

  imgx_s = mg%is(1)-lg%is(1)+1
  imgx_e = mg%is(1)-lg%is(1)+mg%num(1)
  imgy_s = mg%is(2)-lg%is(2)+1
  imgy_e = mg%is(2)-lg%is(2)+mg%num(2)
  imgz_s = mg%is(3)-lg%is(3)+1
  imgz_e = mg%is(3)-lg%is(3)+mg%num(3)

  ifgx_s = (mg%is(1)-lg%is(1))*2+1
  ifgx_e = (mg%is(1)-lg%is(1))*2+mg%num(1)*2
  ifgy_s = (mg%is(2)-lg%is(2))*2+1
  ifgy_e = (mg%is(2)-lg%is(2))*2+mg%num(2)*2
  ifgz_s = (mg%is(3)-lg%is(3))*2+1
  ifgz_e = (mg%is(3)-lg%is(3))*2+mg%num(3)*2

  poisson%ff1z = 0d0

!$OMP parallel do private(iz,iy,ix,iix,iiy,iiz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    iix = ix-lg%is(1)+1
    iiy = iy-lg%is(2)+1
    iiz = iz-lg%is(3)+1
    poisson%ff1z(iix,iiy,iiz) = dcmplx(rho%f(ix,iy,iz))
  end do
  end do
  end do

  call comm_summation(poisson%ff1z,poisson%ff2z,mg%num(1)*mg%num(2)*2*lg%num(3),info%icomm_z)

  poisson%ff1y = 0d0

!$OMP parallel do private(kz,iy,ix)
  do kz = ifgz_s,ifgz_e
  do iy = imgy_s,imgy_e
  do ix = imgx_s,imgx_e
    poisson%ff1y(ix,iy,kz) = sum(fg%egzc(kz,:)*poisson%ff2z(ix,iy,:))
  end do
  end do
  end do

  call comm_summation(poisson%ff1y,poisson%ff2y,mg%num(1)*2*lg%num(2)*2*mg%num(3),info%icomm_y)

  poisson%ff1x = 0.d0

!$OMP parallel do private(kz,ky,ix)
  do kz = ifgz_s,ifgz_e
  do ky = ifgy_s,ifgy_e
  do ix = imgx_s,imgx_e
    poisson%ff1x(ix,ky,kz) = sum(fg%egyc(ky,:)*poisson%ff2y(ix,:,kz))
  end do
  end do
  end do

  call comm_summation(poisson%ff1x,poisson%ff2x,2*lg%num(1)*2*mg%num(2)*2*mg%num(3),info%icomm_x)

  poisson%ff1 = 0.d0

!$OMP parallel do private(kz,ky,kx)
  do kz = ifgz_s,ifgz_e
  do ky = ifgy_s,ifgy_e
  do kx = ifgx_s,ifgx_e
    poisson%ff1(kx,ky,kz) = sum(fg%egxc(kx,:)*poisson%ff2x(:,ky,kz))/dble(8*lg%num(1)*lg%num(2)*lg%num(3))
  end do
  end do
  end do

  call comm_summation(poisson%ff1,poisson%ff2,2*lg%num(1)*2*mg%num(2)*2*mg%num(3),info%icomm_x)

  poisson%ff3z = 0.d0

!$OMP parallel do private(kz,ky,kx)
  do kz = ifgz_s,ifgz_e
  do ky = ifgy_s,ifgy_e
  do kx = ifgx_s,ifgx_e
    poisson%zrhoG_ele(kx,ky,kz) = poisson%ff4z(kx,ky,kz)
    poisson%ff3z(kx,ky,kz) = fg%coef(kx,ky,kz)*poisson%ff2(kx,ky,kz)
  end do
  end do
  end do

  call comm_summation(poisson%ff3z,poisson%ff4z,2*mg%num(1)*2*mg%num(2)*2*lg%num(3),info%icomm_z)

  poisson%ff3y = 0.d0

!$OMP parallel do private(iz,ky,kx)
  do iz = imgz_s,imgz_e
  do ky = ifgy_s,ifgy_e
  do kx = ifgx_s,ifgx_e
    poisson%ff3y(kx,ky,iz) = sum(fg%egz(:,iz)*poisson%ff4z(kx,ky,:))
  end do
  end do
  end do

  call comm_summation(poisson%ff3y,poisson%ff4y,2*mg%num(1)*2*lg%num(2)*mg%num(3),info%icomm_y)

  poisson%ff3x = 0.d0

!$OMP parallel do private(iz,iy,kx)
  do iz = imgz_s,imgz_e
  do iy = imgy_s,imgy_e
  do kx = ifgx_s,ifgx_e
    poisson%ff3x(kx,iy,iz) = sum(fg%egy(:,iy)*poisson%ff4y(kx,:,iz))
  end do
  end do
  end do

  call comm_summation(poisson%ff3x,poisson%ff4x,2*lg%num(1)*mg%num(2)*mg%num(3),info%icomm_x)

!$OMP parallel do private(iz,iy,ix,iix,iiy,iiz)
  do iz = mg%is(3),mg%ie(3)
  do iy = mg%is(2),mg%ie(2)
  do ix = mg%is(1),mg%ie(1)
    iix = ix-lg%is(1)+1
    iiy = iy-lg%is(2)+1
    iiz = iz-lg%is(3)+1
    Vh%f(ix,iy,iz) = sum(fg%egx(:,iix)*poisson%ff4x(:,iiy,iiz))
  end do
  end do
  end do

end subroutine poisson_isolated_ft

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

subroutine poisson_isolated_ffte(lg,mg,info,fg,rho,Vh,poisson)
  use structures
  use communication, only: comm_summation, comm_bcast
  implicit none
  type(s_rgrid)          ,intent(in) :: lg
  type(s_rgrid)          ,intent(in) :: mg
  type(s_parallel_info)  ,intent(in) :: info
  type(s_reciprocal_grid),intent(in) :: fg
  type(s_scalar)         ,intent(in) :: rho
  type(s_scalar)                     :: Vh
  type(s_poisson)                    :: poisson
  !
  integer :: ix,iy,iz
  integer :: iiy,iiz,iix
  real(8) :: inv_lgnum3

  integer :: ishift_y,ishift_z
  integer :: nproc_ob

  inv_lgnum3=1.d0/(8*lg%num(1)*lg%num(2)*lg%num(3))

  poisson%a_ffte=0.d0
  poisson%b_ffte=0.d0

  nproc_ob  = info%nporbital

  if(mod(nproc_ob,4)==0)then

    if(info%id_o_isolated_ffte == 0) then
!$OMP parallel do private(iiz,iiy,ix) collapse(2)
      do iz=1,mg%num(3)
      do iy=1,mg%num(2)
        iiz=iz+mg%is(3)-1
        iiy=iy+mg%is(2)-1
        poisson%b_ffte(1:lg%num(1),iy,iz) = dcmplx(rho%f(lg%is(1):lg%ie(1),iiy,iiz))
      end do
      end do
    end if

    call comm_summation(poisson%b_ffte,poisson%a_ffte,size(poisson%a_ffte),info%icomm_x)

    CALL PZFFT3DV_MOD(poisson%a_ffte,poisson%b_ffte,2*lg%num(1),2*lg%num(2),2*lg%num(3),   &
                      info%isize_y_isolated_ffte,info%isize_z_isolated_ffte,-1, &
                      info%icomm_y_isolated_ffte,info%icomm_z_isolated_ffte)

    poisson%zrhoG_ele=0d0
    if(info%id_z_isolated_ffte>=info%isize_z_isolated_ffte/2)then
      ishift_z = lg%num(3)
    else
      ishift_z = 0
    end if
    if(info%id_y_isolated_ffte>=info%isize_y_isolated_ffte/2)then
      ishift_y = lg%num(2)
    else
      ishift_y = 0
    end if

!$omp parallel do collapse(2) default(none) &
!$omp             private(iz,iy,ix,iiy,iiz,iix) &
!$omp             shared(mg,lg,poisson,inv_lgnum3,fg,info,ishift_z,ishift_y)
    do iz=1,mg%num(3)
    do iy=1,mg%num(2)
      iiz=iz+mg%is(3)-lg%is(3)+ishift_z
      iiy=iy+mg%is(2)-lg%is(2)+ishift_y
      do ix=1,2*lg%num(1)
        iix=ix+lg%is(1)-1
        poisson%zrhoG_ele(iix,iiy,iiz) = poisson%b_ffte(ix,iy,iz)*inv_lgnum3
      end do
      do ix=1,2*lg%num(1)
        poisson%b_ffte(ix,iy,iz) = poisson%b_ffte(ix,iy,iz) * fg%coef(ix,iiy,iiz)
      end do
    end do
    end do
!$omp end parallel do

    CALL PZFFT3DV_MOD(poisson%b_ffte,poisson%a_ffte,2*lg%num(1),2*lg%num(2),2*lg%num(3), &
                      info%isize_y_isolated_ffte,info%isize_z_isolated_ffte,1, &
                      info%icomm_y_isolated_ffte,info%icomm_z_isolated_ffte)

    call comm_bcast(poisson%a_ffte,info%icomm_o_isolated_ffte,0)

!$OMP parallel do private(iiz,iiy) collapse(2)
    do iz=1,mg%num(3)
    do iy=1,mg%num(2)
      iiz=iz+mg%is(3)-1
      iiy=iy+mg%is(2)-1
      Vh%f(lg%is(1):lg%ie(1),iiy,iiz) = poisson%a_ffte(1:lg%num(1),iy,iz)
    end do
    end do

  else

!$OMP parallel do private(iiz,iiy) collapse(2)
    do iz=1+mg%is(3)-lg%is(3),mg%num(3)+mg%is(3)-lg%is(3)
    do iy=1+mg%is(2)-lg%is(2),mg%num(2)+mg%is(2)-lg%is(2)
    do ix=1+mg%is(1)-lg%is(1),mg%num(1)+mg%is(1)-lg%is(1)
      iiz=iz+lg%is(3)-1
      iiy=iy+lg%is(2)-1
      iix=ix+lg%is(1)-1
      poisson%b_ffte(ix,iy,iz) = dcmplx(rho%f(iix,iiy,iiz))
    end do
    end do
    end do

    call comm_summation(poisson%b_ffte,poisson%a_ffte,size(poisson%a_ffte),info%icomm_r)

    CALL PZFFT3DV_MOD(poisson%a_ffte,poisson%b_ffte,2*lg%num(1),2*lg%num(2),2*lg%num(3),   &
                      info%isize_y_isolated_ffte,info%isize_z_isolated_ffte,-1, &
                      info%icomm_y_isolated_ffte,info%icomm_z_isolated_ffte)

    poisson%zrhoG_ele=0d0

!$omp parallel do collapse(2) default(none) &
!$omp             private(iz,iy,ix) &
!$omp             shared(mg,lg,poisson,inv_lgnum3,fg,info,ishift_z,ishift_y)
    do iz=1,2*lg%num(3)
    do iy=1,2*lg%num(2)
      do ix=1,2*lg%num(1)
        poisson%zrhoG_ele(ix,iy,iz) = poisson%b_ffte(ix,iy,iz)*inv_lgnum3
      end do
      do ix=1,2*lg%num(1)
        poisson%b_ffte(ix,iy,iz) = poisson%b_ffte(ix,iy,iz) * fg%coef(ix,iy,iz)
      end do
    end do
    end do
!$omp end parallel do

    CALL PZFFT3DV_MOD(poisson%b_ffte,poisson%a_ffte,2*lg%num(1),2*lg%num(2),2*lg%num(3), &
                      info%isize_y_isolated_ffte,info%isize_z_isolated_ffte,1, &
                      info%icomm_y_isolated_ffte,info%icomm_z_isolated_ffte)

!$OMP parallel do private(iiz,iiy) collapse(2)
    do iz=1+mg%is(3)-lg%is(3),mg%num(3)+mg%is(3)-lg%is(3)
    do iy=1+mg%is(2)-lg%is(2),mg%num(2)+mg%is(2)-lg%is(2)
      iiz=iz+lg%is(3)-1
      iiy=iy+lg%is(2)-1
      Vh%f(mg%is(1):mg%ie(1),iiy,iiz) = poisson%a_ffte(1+mg%is(1)-lg%is(1):mg%num(1)+mg%is(1)-lg%is(1),iy,iz)
    end do
    end do

  end if

  return
end subroutine poisson_isolated_ffte

end module poisson_isolated


! { dg-do compile }
! This checks the fix for PR 26064
!
! Contributed by Sven Buijssen <sven.buijssen@math.uni-dortmund.de>
module ice
  implicit none
  contains

    subroutine foo()
    contains

      subroutine bar(baz)
        integer, optional :: baz
        if (present(baz)) then
        endif
      end subroutine bar
    end subroutine foo
end module


module boundary_conditions

    use interface
    use define_indexes
    use define_keys
    use define_globals
    use utility, only: util_print_warning
    use utility_interpolate
    use define_settings, only: setting
    use face, only: face_interpolate_bc
    use define_xsect_tables
    use xsect_tables
    use utility_crash

    implicit none

    private

    public :: bc_step
    public :: bc_update

contains
!%
!%==========================================================================
!% PUBLIC
!%==========================================================================
!%
    subroutine bc_update()
        !%------------------------------------------------------------------
        !% Description:
        !%------------------------------------------------------------------
        !% Declarations:
            integer :: ii
            logical :: isBConly = .true.
            character(64) :: subroutine_name = "bc_update"
        !%------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%boundary_conditions)  &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%------------------------------------------------------------------

        !% Gets boundary flow and head data from SWMM-C and stores in the BC structure
        call bc_step()
        if (crashI==1) return

        ! print *, 'in bcupdate AAA'
        ! print *, faceR(47:48,fr_Head_u)
        ! print *, faceR(47:48,fr_Head_d)

        if (N_flowBC > 0 .or. N_headBC > 0) then
            !% interpolate the BC in time
            !print *, 'in 223563 ',trim(subroutine_name)
            !write(*,'(8f9.2)') faceR(1:N_elem(1)+1,fr_Flowrate)
            !print *, faceR(1:N_elem(1)+1,fr_Flowrate)

            call bc_interpolate()

            ! print *, 'in bcupdate BBB'
            ! print *, faceR(47:48,fr_Head_u)
            ! print *, faceR(47:48,fr_Head_d)

            !print *, 'in 4789615 ',trim(subroutine_name)
            !write(*,'(8f9.2)') faceR(1:N_elem(1)+1,fr_Flowrate)

            call face_interpolate_bc(isBConly) ! broadcast interpolation to face & elem arrays

            ! print *, 'in bcupdate CCC'
            ! print *, faceR(47:48,fr_Head_u)
            ! print *, faceR(47:48,fr_Head_d)

            !print *, 'in 836673 ',trim(subroutine_name)
            !    write(*,'(8f9.2)') faceR(1:N_elem(1)+1,fr_Flowrate)
        end if

        !%------------------------------------------------------------------
        !% Closing
            if (setting%Debug%File%boundary_conditions) then
                print *, "INFLOW BC"
                print *, "BC times"
                do ii = 1, setting%BC%TimeSlotsStored
                    print *, BC%flowR_timeseries(:, ii, br_time)
                end do
                print *, "BC values"
                do ii = 1, setting%BC%TimeSlotsStored
                    print *, BC%flowR_timeseries(:, ii, br_value)
                end do
                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

                print *, "HEAD BC"
                print *, "BC times"
                do ii = 1, setting%BC%TimeSlotsStored
                    print *, BC%headR_timeseries(:, ii, br_time)
                end do
                print *, "BC values"
                do ii = 1, setting%BC%TimeSlotsStored
                    print *, BC%headR_timeseries(:, ii, br_value)
                end do

                write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
            end if

    end subroutine bc_update
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine bc_step()
        !%-----------------------------------------------------------------------------
        integer :: ii, nidx, lidx
        integer :: tstep_larger_than_resolution
        real(8) :: ttime, tnow, tend
        character(64) :: subroutine_name = "bc_step"
        !%-----------------------------------------------------------------------------
        !if (crashYN) return
        if (setting%Debug%File%boundary_conditions)  &
            write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

        tnow = setting%Time%Now
        tend = setting%Time%End

        !% --- flow BC
        if (N_flowBC > 0) then
            do ii = 1, N_flowBC
                nidx = BC%flowI(ii, bi_node_idx)
                if (BC%flowYN(ii,bYN_read_input_file)) then
                    if (BC%flowIdx(ii) == 0) then ! First fetch
                        call bc_fetch_flow(ii)
                    else
                        ttime = BC%flowR_timeseries(ii, BC%flowIdx(ii), br_time) ! Current time slot (upper bound of time interval)
                        if (tnow > ttime) then ! Needs update
                            if (BC%flowIdx(ii) == setting%BC%TimeSlotsStored) then
                                call bc_fetch_flow(ii)
                            else
                                tstep_larger_than_resolution = -1
                                do while((tnow > ttime) .and. (BC%flowIdx(ii) < setting%BC%TimeSlotsStored))
                                    BC%flowIdx(ii) = BC%flowIdx(ii) + 1
                                    ttime = BC%flowR_timeseries(ii, BC%flowIdx(ii), br_time)
                                    if ((BC%flowIdx(ii) == setting%BC%TimeSlotsStored) .and. (tnow > ttime)) then
                                        call bc_fetch_flow(ii)
                                        ttime = BC%flowR_timeseries(ii, BC%flowIdx(ii), br_time)
                                    end if
                                    tstep_larger_than_resolution = tstep_larger_than_resolution + 1
                                end do
                                if ((tstep_larger_than_resolution > 0) .and. setting%Output%Warning) then
                                    call util_print_warning("Warning (bc.f08): The flow boundary condition for node " &
                                    // node%Names(nidx)%str // " has a higher resolution than the time step")
                                end if
                            end if
                        end if
                    end if
                else
                    call util_print_warning("Error (bc.f08): The flow boundary condition for node " &
                    // node%Names(nidx)%str // " should always read from an input file")
                    !stop
                    call util_crashpoint( 87453)
                    !return
                end if
            end do
        end if

        !% --- Head BC (outfall)
        if (N_headBC > 0) then
            do ii = 1, N_headBC
                nidx = BC%headI(ii, bi_node_idx)
                if (BC%headYN(ii,bYN_read_input_file)) then
                    if (BC%headIdx(ii) == 0) then ! First fetch
                        call bc_fetch_head(ii)
                        !% HACK - we are assuming that outfalls can only have one link upstream
                        lidx = node%I(nidx, ni_Mlink_u1)
                        link%R(lidx, lr_InitialDnstreamDepth) = BC%headR_timeseries(ii, 1, br_value) - node%R(nidx,nr_Zbottom)
                    else
                        ttime = BC%headR_timeseries(ii, BC%headIdx(ii), br_time) ! Current time slot (upper bound of time interval)
                        if (tnow > ttime) then ! Needs update
                            if (BC%headIdx(ii) == setting%BC%TimeSlotsStored) then
                                call bc_fetch_head(ii)
                            else
                                tstep_larger_than_resolution = -1
                                do while((tnow > ttime) .and. (BC%headIdx(ii) < setting%BC%TimeSlotsStored))
                                    BC%headIdx(ii) = BC%headIdx(ii) + 1
                                    ttime = BC%headR_timeseries(ii, BC%headIdx(ii), br_time)
                                    if ((BC%headIdx(ii) == setting%BC%TimeSlotsStored) .and. (tnow > ttime)) then
                                        call bc_fetch_head(ii)
                                        ttime = BC%headR_timeseries(ii, BC%headIdx(ii), br_time)
                                    end if
                                    tstep_larger_than_resolution = tstep_larger_than_resolution + 1
                                end do
                                if ((tstep_larger_than_resolution > 0) .and. setting%Output%Warning) then
                                    call util_print_warning("Warning (bc.f08): The head boundary condition for node " &
                                    // node%Names(nidx)%str // " has a higher resolution than the time step")
                                end if
                            end if
                        end if
                    end if
                end if
            end do
        end if

        if (setting%Debug%File%boundary_conditions)  &
            write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine bc_step
!%
!%==========================================================================
!% PRIVATE
!%==========================================================================
!%
    subroutine bc_fetch_flow(bc_idx)
        !%-------------------------------------------------------------------
        !% Description
        !% To prevent having to read the inflow file at every time step
        !% we store a number of "slots" of data controlled by the setting
        !% parameter "setting.BC.TimeSlotsStored". This code reads in and
        !% stores from the inflows in the SWMM input files
        !%-------------------------------------------------------------------
            integer, intent(in) :: bc_idx
            integer             :: ii
            integer, pointer    ::  NN
            real(8)             :: new_inflow_time
            real(8)             :: new_inflow_value
            character(64)       :: subroutine_name = "bc_fetch_flow"
        !%-------------------------------------------------------------------
        !% Preliminaries
            !if (crashYN) return
            if (setting%Debug%File%boundary_conditions)  &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------
        !% Aliases
            NN => setting%BC%TimeSlotsStored
        !%-------------------------------------------------------------------

        !% get the first of the next set of data
        if (BC%flowIdx(bc_idx) == 0) then ! First fetch is for the simulation start time
            BC%flowR_timeseries(bc_idx, 1, br_time)  = setting%Time%Start
            BC%flowR_timeseries(bc_idx, 1, br_value) = interface_get_flowBC(bc_idx, setting%Time%Start)
        else ! last value becomes first in the new storred set
            BC%flowR_timeseries(bc_idx, 1, br_time)  = BC%flowR_timeseries(bc_idx, NN, br_time)
            BC%flowR_timeseries(bc_idx, 1, br_value) = BC%flowR_timeseries(bc_idx, NN, br_value)
        end if

        !% read in additional data
        new_inflow_time = setting%Time%Start
        do ii = 2, NN
            !print *, 'from interface ',interface_get_next_inflow_time(bc_idx, new_inflow_time)
            new_inflow_time = min(setting%Time%End, interface_get_next_inflow_time(bc_idx, new_inflow_time))
            BC%flowR_timeseries(bc_idx, ii, br_time) = new_inflow_time
            BC%flowR_timeseries(bc_idx, ii, br_value) = interface_get_flowBC(bc_idx, new_inflow_time)
            !print *,  ii, new_inflow_time
            if (new_inflow_time == setting%Time%End) exit
        end do
        BC%flowIdx(bc_idx) = 2

        if (setting%Debug%File%boundary_conditions) then
            do ii = 1, NN
                write(*, "(*(G0.4 : ','))") BC%flowR_timeseries(bc_idx, ii, :)
            end do
            write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        end if

    end subroutine bc_fetch_flow
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine bc_fetch_head(bc_idx)
        !%-----------------------------------------------------------------------------
        integer, intent(in) :: bc_idx
        integer             :: ii, NN
        real(8)             :: new_head_time
        real(8)             :: new_head_value, normDepth, critDepth
        integer, pointer    :: nodeIdx, faceIdx, elemUpIdx
        character(64)       :: subroutine_name = "bc_fetch_head"
        !%-----------------------------------------------------------------------------
        !if (crashYN) return
        if (setting%Debug%File%boundary_conditions)  &
            write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

        NN = setting%BC%TimeSlotsStored

        if (BC%headIdx(bc_idx) == 0) then ! First fetch

            BC%headR_timeseries(bc_idx, 1, br_time) = setting%Time%Start
            BC%headR_timeseries(bc_idx, 1, br_value) &
                 = interface_get_headBC(bc_idx, setting%Time%Start) &
                 - setting%Solver%ReferenceHead

            ! print *, ' '
            ! print *, 'in AAA ',trim(subroutine_name)
            ! print *, setting%Solver%ReferenceHead, BC%headR_timeseries(bc_idx, 1, br_value)

        else ! last value becomes first
            BC%headR_timeseries(bc_idx, 1, br_time) = BC%headR_timeseries(bc_idx, NN, br_time)
            BC%headR_timeseries(bc_idx, 1, br_value)  &
                = BC%headR_timeseries(bc_idx, NN, br_value)

                ! print *, ' '
                ! print *, 'in BBB ',trim(subroutine_name)
                ! print *, setting%Solver%ReferenceHead, BC%headR_timeseries(bc_idx, 1, br_value)
        end if

        do ii = 2, NN
            new_head_time = min(setting%Time%End, interface_get_next_head_time(bc_idx, setting%Time%Start))
            BC%headR_timeseries(bc_idx, ii, br_time) = new_head_time
            BC%headR_timeseries(bc_idx, ii, br_value) &
                 = interface_get_headBC(bc_idx, new_head_time) &
                 - setting%Solver%ReferenceHead

            ! print *, ' '
            ! print *, 'in CCC ',trim(subroutine_name)
            ! print *, setting%Solver%ReferenceHead, BC%headR_timeseries(bc_idx, 1, br_value)

            if (new_head_time == setting%Time%End) exit
        end do
        BC%headIdx(bc_idx) = 2

        if (setting%Debug%File%boundary_conditions) &
            write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
    end subroutine bc_fetch_head
!%
!%==========================================================================
!%==========================================================================
!%
    subroutine bc_interpolate()
        !%-------------------------------------------------------------------
        !% Description:
        !% This subroutine is for boundary condition interpolation.
        !% Base on the time passed from the time loop, we interpolate (linear interpolation for now)
        !% the boundary condition to get the corresponding value.
        !%-------------------------------------------------------------------
        !% Declarations:
            real(8) :: normDepth, critDepth
            real(8), pointer :: tnow
            integer :: ii,  lower_idx
            integer, pointer :: nodeIdx, faceIdx, elemUpIdx, upper_idx
            character(64) :: subroutine_name = 'bc_interpolate'
        !%-------------------------------------------------------------------
        !% Preliminaries:
            if (setting%Debug%File%boundary_conditions)  &
                write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"
        !%-------------------------------------------------------------------
        !% Aliases
            tnow => setting%Time%Now
        !%-------------------------------------------------------------------

        do ii=1, N_flowBC
            !slot_idx  => BC%flowIdx(ii)
            upper_idx => BC%flowIdx(ii)
            !upper_idx = slot_idx
            lower_idx = upper_idx - 1

            !% Find the closest index first, assign it to lower_idx for now
            if (BC%flowR_timeseries(ii, lower_idx, br_time) == tnow) then
                !% exact match -- take the existing BC data, no need to do the interpolation
                BC%flowRI(ii) = BC%flowR_timeseries(ii, lower_idx, br_value)

            else if (lower_idx > 0) then
                if ( BC%flowR_timeseries(ii, lower_idx, br_value) == BC%flowR_timeseries(ii, upper_idx, br_value)) then
                    !% constant value, no need to do the interpolation
                    BC%flowRI(ii) = BC%flowR_timeseries(ii, lower_idx, br_value)
                else
                    !% interpolation step
                    if (.not. setting%BC%disableInterpolationYN) then
                        BC%flowRI(ii) = util_interpolate_linear( &
                            tnow, &
                            BC%flowR_timeseries(ii, lower_idx, br_time), &
                            BC%flowR_timeseries(ii, upper_idx, br_time), &
                            BC%flowR_timeseries(ii, lower_idx, br_value), &
                            BC%flowR_timeseries(ii, upper_idx, br_value))
                    else
                        !% no interpolatoni -- take the upper index value
                        BC%flowRI(ii) = BC%flowR_timeseries(ii, upper_idx, br_value)
                    end if
                end if
            else
                !% lower_idx <= 0
                write(*,*), 'CODE ERROR? unexpected else in BC'
                !stop
                call util_crashpoint( 987098)
                !return
            end if
        end do

        do ii=1, N_headBC
            nodeIdx     => BC%headI(ii,bi_node_idx)
            faceIdx     => BC%headI(ii,bi_face_idx)
            elemUpIdx   => faceI(faceIdx,fi_Melem_uL)
            !slot_idx    => BC%headIdx(ii)
            upper_idx   => BC%headIdx(ii)
            lower_idx   =  upper_idx - 1

            !% --- prescribed head at outlet
            if ((BC%headI(ii,bi_subcategory) == BCH_fixed)   .or. &
                (BC%headI(ii,bi_subcategory) == BCH_tseries) .or. &
                (BC%headI(ii,bi_subcategory) == BCH_tidal)) then

                !% Find the closest index first, assign it to lower_idx just for now
                if (BC%headR_timeseries(ii, lower_idx, br_time) == tnow) then
                    !% exact match -- directly take the existing BC data, no need to do the interpolation,
                    BC%headRI(ii) = BC%headR_timeseries(ii, lower_idx, br_value)
                else if (lower_idx > 0) then
                    if ( BC%headR_timeseries(ii, lower_idx, br_value) == BC%headR_timeseries(ii, upper_idx, br_value)) then
                        !% constant value, no need to do the interpolation
                        BC%headRI(ii) = BC%headR_timeseries(ii, lower_idx, br_value)
                    else
                        !% do the interpolation -- NOTE the setting..disableInterpolationYN not available here
                        BC%headRI(ii) = util_interpolate_linear( &
                            tnow, &
                            BC%headR_timeseries(ii, lower_idx, br_time), &
                            BC%headR_timeseries(ii, upper_idx, br_time), &
                            BC%headR_timeseries(ii, lower_idx, br_value), &
                            BC%headR_timeseries(ii, upper_idx, br_value))
                    end if
                else
                    !% lower_idx <= 0
                    write(*,*), 'CODE ERROR? unexpected else in BC'
                    !stop
                    call util_crashpoint( 786985)
                    !return
                end if

            !% --- normal flow at outlet
            else if (BC%headI(ii,bi_subcategory) == BCH_normal) then

                !% for normal dnBC, if the connecting link has an offset,
                !% the depth in the node is zero
                if (link%R(node%I(nodeIdx,ni_Mlink_u1),lr_OutletOffset) > zeroR) then
                    BC%headRI(ii) = faceR(faceIdx,fr_Zbottom)
                else
                    if (elemI(elemUpIdx,ei_elementType) == CC) then
                        BC%headRI(ii) = elemR(elemUpIdx,er_Depth) + faceR(faceIdx,fr_Zbottom)
                    else
                        !% for normal dnBC, if the upstream link is not CC (i.e. weir, orifice etc)
                        !% the depth in the node is zero
                        BC%headRI(ii) =  faceR(faceIdx,fr_Zbottom)
                    end if
                end if

            !% --- free overflow at outlet
            else if (BC%headI(ii,bi_subcategory) == BCH_free) then

                !% for free dnBC, if the connecting link has an offset,
                !% the depth in the node is zero
                ! print*, link%R(node%I(nodeIdx,ni_Mlink_u1),lr_OutletOffset), 'outlet offset'
                if (link%R(node%I(nodeIdx,ni_Mlink_u1),lr_OutletOffset) > zeroR) then
                    BC%headRI(ii) = faceR(faceIdx,fr_Zbottom)
                else
                    if (elemI(elemUpIdx,ei_elementType) == CC) then
                        normDepth = elemR(elemUpIdx,er_Depth)
                        critDepth = bc_get_CC_critical_depth(elemUpIdx)
                        BC%headRI(ii) = min(critDepth,normDepth) + faceR(faceIdx,fr_Zbottom)
                    else
                        !% for free dnBC, if the upstream link is not CC (i.e. weir, orifice etc)
                        !% the depth in the node is zero
                        BC%headRI(ii) =  faceR(faceIdx,fr_Zbottom)
                    end if
                end if
            !% --- error in specifying outlet
            else
                call util_print_warning("Error (bc.f08): Unknown downstream boundary condition type at " &
                    // node%Names(nodeIdx)%str // " node")
                !stop
                call util_crashpoint( 86474)
                !return
            end if
        end do

        if (setting%Debug%File%boundary_conditions) &
            write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

    end subroutine bc_interpolate
!%
!%==========================================================================
!%==========================================================================
!%
function bc_get_CC_critical_depth(elemIdx) result (criticalDepth)
    !%-----------------------------------------------------------------------------
    !% Description:
    !% gets the critical depth of a CC element
    !% this code has been adapted from SWMM5-C code
    !%-----------------------------------------------------------------------------
        integer, intent(in)  :: elemIdx
        real(8)              :: criticalDepth
        integer, pointer     :: elemGeometry
        real(8), pointer     :: Flowrate, FullDepth, BreadthMax, grav, FullArea
        real(8), pointer     :: alpha, bottomWidth
        real(8)              :: Q2g, critDepthEstimate, ratio, sideSlope, epsilon_c, tc0
        character(64) :: subroutine_name = 'bc_get_CC_critical_depth'
    !%-----------------------------------------------------------------------------
        !if (crashYN) return
        if (setting%Debug%File%boundary_conditions)  &
            write(*,"(A,i5,A)") '*** enter ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

        elemGeometry => elemI(elemIdx,ei_geometryType)
        Flowrate     => elemR(elemIdx,er_Flowrate)
        FullDepth    => elemR(elemIdx,er_FullDepth)
        FullArea     => elemR(elemIdx,er_FullArea)
        BreadthMax   => elemR(elemIdx,er_BreadthMax)
        grav         => setting%constant%gravity
        alpha        => setting%constant%energy_correction_factor

        Q2g = (Flowrate**twoR) / grav

        if (Q2g == zeroR) then
            criticalDepth =  zeroR
        else
            select case (elemGeometry)

            case (rectangular, rectangular_closed)
                criticalDepth = (Q2g/BreadthMax**twoR)**onethirdR

            case (triangular)
                sideSlope = elemSGR(elemIdx,esgr_Triangular_Slope)
                criticalDepth = (twoR * Q2g / sideSlope ** twoR) ** onefifthR

            case (trapezoidal)
                !% use the average side slope
                sideSlope = onehalfR * (  elemSGR(elemIdx,esgr_Trapezoidal_LeftSlope)   &
                                        + elemSGR(elemIdx,esgr_Trapezoidal_RightSlope))
                bottomWidth => elemSGR(elemIdx,esgr_Trapezoidal_Breadth)
                !% Using approach of Vatakhah (2013)
                !% non-dimensional discharge (epsilon_c), eq. 14
                epsilon_c = fourR * sideSlope * ( alpha * (Flowrate**2) / (grav * (bottomWidth**5))  )**(onethirdR)
                !% non-dimensional critical flow estimat (tc0), eq. 23
                tc0 = (oneR + 1.161d0 * epsilon_c * (oneR + 0.666d0 * (epsilon_c**(1.041d0)) )**0.374d0 )**0.144d0
                !% critical depth, eq. 22
                criticalDepth = -onehalfR + onehalfR             &
                          * (                                    &
                                (fiveR * (tc0**6) + oneR       ) &
                              / ( sixR * (tc0**5) - epsilon_c  ) &
                            )**3

            case (circular)
                !% first estimate Critical Depth for an equivalent circular conduit
                critDepthEstimate = min(1.01*(Q2g/FullDepth)**onefourthR, FullDepth)

                !% find ratio of conduit area to equiv. circular area
                ratio = FullArea/(pi/fourR*(FullDepth**twoR))

                if ((ratio >= onehalfR) .and. (ratio <= twoR)) then
                    criticalDepth = bc_critDepth_enum (elemIdx, Flowrate, critDepthEstimate)
                else
                    criticalDepth = bc_critDepth_ridder (elemIdx, Flowrate, critDepthEstimate)
                end if

            case default
                print *, 'in ',trim(subroutine_name)
                print *, 'CODE ERROR: unknown geometry of # ',elemGeometry
                print *, 'which has key ',trim(reverseKey(elemGeometry))
                !stop
                call util_crashpoint( 389753)
                !return
            end select
        end if

        if (setting%Debug%File%boundary_conditions) &
            write(*,"(A,i5,A)") '*** leave ' // trim(subroutine_name) // " [Processor ", this_image(), "]"

end function bc_get_CC_critical_depth
!%
!%==========================================================================
!%==========================================================================
!%
function bc_critDepth_enum (elemIdx, Q, yEstimate) result (criticalDepth)
    !%-----------------------------------------------------------------------------
    !% Description:
    !% gets the critical depth of a CC by enumeration method
    !% this code has been adapted from SWMM5-C code
    !%-----------------------------------------------------------------------------
        integer, intent(in)  :: elemIdx
        real(8), intent(in)  :: Q, yEstimate
        real(8)              :: criticalDepth
        real(8), pointer     :: FullDepth
        real(8)              :: dY, Y, Q0, QC
        integer              :: i1, ii
        character(64)        :: subroutine_name = 'bc_critDepth_enum'
    !%-----------------------------------------------------------------------------

    FullDepth => elemR(elemIdx,er_FullDepth)

    !% find the critical depth from SWMM5 ennumeration method
    dY = FullDepth/25.0
    i1 = int(yEstimate/dY)
    Q0 = bc_get_critical_flow(elemIdx,i1*dY,zeroR)

    if (Q0 < Q) then
        criticalDepth = FullDepth
        do  ii = i1+1,25
            QC = bc_get_critical_flow(elemIdx,ii*dY,zeroR)
            if (QC > Q) then
                criticalDepth = ((Q-Q0)/(QC-Q0)+ real((ii-1),8))*dY
                return
            end if
            Q0 = QC
        end do
    else
        criticalDepth = zeroR
        do  ii = i1-1,0,-1
            QC = bc_get_critical_flow(elemIdx,ii*dY,zeroR)
            if (QC < Q) then
                criticalDepth = ((Q-QC)/(Q0-QC)+ real((ii),8))*dY
                return
            end if
            Q0 = QC
        end do
    end if

end function bc_critDepth_enum
!%
!%==========================================================================
!%==========================================================================
!%
function bc_critDepth_ridder (elemIdx, Q, yEstimate) result (criticalDepth)
    !%-----------------------------------------------------------------------------
    !% Description:
    !% gets the critical depth of a CC by ridder method
    !% this code has been adapted from SWMM5-C code
    !%-----------------------------------------------------------------------------
        integer, intent(in)  :: elemIdx
        real(8), intent(in)  :: Q, yEstimate
        real(8)              :: criticalDepth, Y1, Y2
        real(8), pointer     :: FullDepth
        real(8)              :: Q0, Q1, Q2, Tol
        character(64)        :: subroutine_name = 'bc_critDepth_ridder'
    !%-----------------------------------------------------------------------------

    FullDepth => elemR(elemIdx,er_FullDepth)

    Y1 = zeroR
    Y2 = 0.99d0*FullDepth

    !% check if critical flow at full depth < target flow
    Q2 = bc_get_critical_flow(elemIdx,Y2,zeroR)
    if (Q2 < Q) then
        criticalDepth = FullDepth

    !% otherwise evaluate critical flow at initial estimate
    !% and 1/2 of fullDepth
    else
        Q0 = bc_get_critical_flow(elemIdx,yEstimate,zeroR)
        Q1 = bc_get_critical_flow(elemIdx, onehalfR*FullDepth,zeroR)

        if (Q0 > Q) then
            Y2 = yEstimate
            if (Q1 < Q) Y1 = onehalfR*FullDepth
        else
            Y1 = yEstimate
            if (Q1 > Q) Y2 = onehalfR*FullDepth
        end if
        Tol = 0.001
        criticalDepth = bc_findRoot_ridder_for_critDepth(elemIdx,Q,Y1,Y2,Tol)
    end if

end function bc_critDepth_ridder
!%
!%==========================================================================
!%==========================================================================
!%
function bc_findRoot_ridder_for_critDepth(elemIdx,Q,Y1,Y2,Tol) result (criticalDepth)
    !%-----------------------------------------------------------------------------
    !% Description:
    !% find root for the critical depth of a CC by ridder method
    !% this code has directly been adapted from SWMM5-C code
    !%-----------------------------------------------------------------------------
        integer, intent(in)  :: elemIdx
        real(8), intent(in)  :: Q, Y1, Y2, Tol
        real(8)              :: criticalDepth
        integer              :: ii, MaxIter
        real(8)              :: ans, fhi, flo, fm, fnew, s, Yhi, Ylo, Ym, Ynew

        character(64)        :: subroutine_name = 'bc_findRoot_ridder_for_critDepth'
    !%-----------------------------------------------------------------------------
        MaxIter = 60
        flo = bc_get_critical_flow(elemIdx,Y1,Q)
        fhi = bc_get_critical_flow(elemIdx,Y2,Q)

        if (flo == zeroR) then
            criticalDepth = Y1
            return
        end if

        if (fhi == zeroR) then
            criticalDepth = Y2
            return
        end if

        ans = onehalfR*(Y1+Y2)

        if ((flo > zeroR .and. fhi < zeroR) .or. (flo < zeroR .and. fhi > zeroR)) then
            Ylo = Y1
            Yhi = Y2
            do ii = 1,MaxIter
                Ym = onehalfR*(Ylo+Yhi)
                fm = bc_get_critical_flow(elemIdx,Ym,Q)
                s = sqrt(fm*fm - flo*fhi)
                if (s == zeroR) then
                    criticalDepth = ans
                    return
                end if

                if (flo >= fhi) then
                    Ynew = Ym + (Ym-Ylo)*fm/s
                else
                    Ynew = Ym + (Ym-Ylo)*(-oneI*fm/s)
                endif

                if(abs(Ynew-ans) <= Tol) then
                    criticalDepth = Ynew
                    return
                end if
                ans = Ynew
                fnew = bc_get_critical_flow(elemIdx,Ynew,Q)

                if (sign(fm,fnew) /= fm) then
                    Ylo = Ym
                    flo = fm
                    Yhi = ans
                    fhi = fnew
                else if (sign(flo,fnew) /= flo) then
                    Yhi = ans
                    fhi = fnew
                else if (sign(fhi,fnew) /= fhi) then
                    Ylo = ans
                    flo = fnew
                else
                    criticalDepth = ans
                    return
                end if

                if ( abs(Yhi - Ylo) <= Tol) then
                    criticalDepth = ans
                    return
                end if
            end do
            criticalDepth = ans
        end if

end function bc_findRoot_ridder_for_critDepth
!%
!%==========================================================================
!%==========================================================================
!%
function bc_get_critical_flow (elemIdx, Depth, Flow_0) result (Flow)
    !%-----------------------------------------------------------------------------
    !% Description:
    !% gets the critical flow of a CC element for a given depth
    !%-----------------------------------------------------------------------------
        integer, intent(in)  :: elemIdx
        real(8), intent(in)  :: Depth, Flow_0
        real(8)              :: Flow
        integer, pointer     :: elemGeometry
        real(8), pointer     :: fullDepth, fullArea, grav
        real(8)              :: YoverYfull, Area, Topwidth
        character(64) :: subroutine_name = 'bc_get_critical_flow'
    !%-----------------------------------------------------------------------------

    elemGeometry => elemI(elemIdx,ei_geometryType)
    fullDepth    => elemR(elemIdx,er_FullDepth)
    fullArea     => elemR(elemIdx,er_FullArea)
    grav         => setting%constant%gravity

    select case (elemGeometry)

    case (circular)
        YoverYfull = Depth/fullDepth
        Area     = fullArea * xsect_table_lookup_singular (YoverYfull, ACirc ) !% 20220506brh removed NACirc
        Topwidth = max(fullDepth * xsect_table_lookup_singular (YoverYfull, TCirc), &   !% 20220506brh removed NTirc
                    setting%ZeroValue%Topwidth)

        Flow  = Area * sqrt(grav*Area/Topwidth) - Flow_0

    case default
        print *, 'in ',trim(subroutine_name)
        print *, 'CODE ERROR: unknown geometry of # ',elemGeometry
        print *, 'which has key ',trim(reverseKey(elemGeometry))
        !stop
        call util_crashpoint( 84198)
        !return
     end select

end function bc_get_critical_flow
!%
!%==========================================================================
!% END OF MODULE
!%+=========================================================================
end module boundary_conditions


text:
program ASA

implicit none
real,allocatable,dimension(:,:)::A,L,U,C,inv
real,allocatable,dimension(:)::sol,B,D,X
real:: det
integer::n,m,f,i

call matriz(A,n,n) !ingresar una matriz nxm
call matcof(A,n,C,inv)  !matriz de cofactores

allocate(B(n))

write(*,*) 'hoa'

do i=1,n
read(*,*) B(i)
end do
call solcof(B,A,n,D)

 write(*,*) 'yea'
do i=1,n
 write(*,*) D(i)
end do

 write(*,*) 'yea'

call kramer(A,n,B,X)
do i=1,n
write(*,*) X(i)
end do
!-----------------------------------------------------------------------
contains
! ----------------------------------------------------------------------
subroutine  matlu(A,n,m,L,U)

! esta subrutina come una matriz A con n filas y m columnas y devuelve una matriz L y una matriz U

implicit none
real,allocatable,dimension(:,:)::A,L,U,D
integer:: k,h,j,n,m

allocate(L(n,m),U(n,m),D(n,m))
D=A
U=A
do h= 1,m
 do j= 1,n
! j filas, h columnas
if (j>h) then
   do k=1,m
  	U(j,k)=A(j,k)-(A(j,h)/A(h,h))*A(h,k)
   end do
   L(j,h)=(A(j,h)*(1/A(h,h)))
else if (j==h) then
 L(j,h)=1
else
  L(j,h)=0
end if
 end do
A=U
end do
A=D
return
end subroutine
!--------------------------------------------------------------

subroutine matdet(A,n,det)
! esta matriz calcula determinante de una matriz A nxn

implicit none
real,allocatable,dimension(:,:)::A,L,U
integer:: i,n
real::det

if (n==1) then
 det=A(1,1)
else
call matlu(A,n,n,L,U)
det=1
do i=1,n
det=U(i,i)*det
end do
end if
return
end subroutine

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

subroutine matcof(A,n,C,inv)
! esta subrutina come una matriz A con n filas y m columnas y devuelve una matriz C de cofactores

implicit none
real,allocatable,dimension(:,:)::A,C,L,U,C2,inv
integer:: k,h,j,n,i,z
real::det
allocate(C(n,n),C2(n-1,n-1))

!h=f,k=c,i=f,j=c
do h=1,n
do k=1,n

do i=1,n
do j=1,n

!ME ARMO UNA MATRIZ SIN ESA FILA Y COLUMNA Y LE HAGO EL DETERMINANTE
if (j/=k.and.i/=h) then

   if (i<h.and.j<k) then
      C2(i,j)=A(i,j)
   else if (i<h.and.j>k) then
      C2(i,j-1)=A(i,j)
   else if (i>h.and.j<k) then
      C2(i-1,j)=A(i,j)
   else if (i>h.and.j>k) then
      C2(i-1,j-1)=A(i,j)
      else
   end if

end if

end do
end do
  call matdet(C2,(n-1),det)
  C(h,k)=det*((-1)**(h+k))
end do

end do

call matdet(A,n,det)

allocate(inv(n,n))

do i=1,n
do j=1,n
 inv(i,j)=C(j,i)/det
end do
end do

deallocate(C2)
return

end subroutine
!--------------------------------------------------------------
subroutine matriz(mat,fila,colum)
!esta subrutina te genera una matriz mat de fila y columna ingresada
real,allocatable,dimension(:,:)::mat
integer:: fila,colum,fila2,colum2,i

Write(*,*) "Ingrese las dimenciones de la matriz"
Write(*,*) "Filas"
read(*,*) fila
Write(*,*) "Columnas"
read(*,*) colum

allocate(mat(fila,colum))

write(*,*) "Ingrese los elementos de la matriz"
fila2=1
colum2=1

do while(fila2<=fila)
 do while(fila2<=fila.and.colum2<=colum)
 Write(*,*)"Ingrese el elemento de la matriz",fila2,",",colum2
 read(*,*)mat(fila2,colum2)
 colum2=colum2+1
end do
fila2=fila2+1
colum2=1
enddo
Write(*,*)"La matriz ingresada fue :"
do i=1,fila,1
write(*,*) mat(i,:)
end do
return
end subroutine

!---------------------------------------------
!error en la division. reveer
subroutine solveLU(L,U,n,B,solu)
! esta subrutina come una matriz L y U de la descomposicion LU, un vector B
! que es el resultado, donde todas tiene n filas y devuelve el vector solucion

implicit none
real,allocatable,dimension(:,:)::L,U
real,allocatable,dimension(:)::solu,B,sol2
integer:: j,i,n,k
real::suma,suma2

allocate(sol2(n),solu(n))
sol2(:)=1.0d0
solu(:)=1.0d0

do i=1,n

if (i==1) then
 sol2(1)=B(1)/L(1,1)
else
suma=0
do k=1,(i-1)
suma=L(i,k)*sol2(k)+suma
end do
sol2(i)=(1./(L(i,i)))*(B(i)-suma)
end if

end do

solu(n)=sol2(n)/U(n,n)

do i=1,(n-1)
suma2=0
do k=(n-i+1),n
suma2=U(n-i,k)*solu(k)+suma2
end do
solu(n-i)=(1./(U(n-i,n-i)))*(sol2(n-i)-suma2)
end do

deallocate (sol2)
return
end subroutine
!------------------------------
subroutine solcof(B,A,n,D)

implicit none
real,allocatable,dimension(:,:)::A,inv,C
real,allocatable,dimension(:)::D,B
integer:: n

call matcof(A,n,C,inv)
allocate(D(n))
D=matmul(inv,B)

return
end subroutine
!----------------------------
subroutine kramer(A,n,B,X)
implicit none
real, allocatable, dimension(:,:)::A,G,F
real, allocatable, dimension(:)::B,X
integer:: n,k,j,i
real:: det,deter

call matdet(A,n,det)

allocate(X(n))

do k=1,n
allocate(F(n,n))

do i=1,n
do j=1,n

if (k==j) then
F(i,j)=B(i)
else
F(i,j)=A(i,j)
end if
end do
end do
deter=0
call matdet(F,n,deter)
X(k)=deter/det
deallocate(F)

end do
return
end subroutine

!--------------------------------
end program

max_stars_repo_name: patricio-c/fortran

max_issues_repo_name: patricio-c/fortran

max_forks_repo_name: patricio-c/fortran


text: #redirect Mother & Baby Source

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!-----------------------------------------------------------------------
!     Module:        beams3d_randomize_particles
!     Authors:       S. Lazerson (samuel.lazerson@ipp.mpg.de)
!     Date:          05/07/2020
!     Description:   This subroutine radomizes the particles so that
!                    no processor follows only slow or fast particles.
!-----------------------------------------------------------------------
      SUBROUTINE beams3d_randomize_particles
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE beams3d_lines, ONLY: nparticles
      USE mpi_params
      USE mpi_inc
      USE beams3d_runtime

!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      INTEGER :: i, d1, d2, b_temp
      INTEGER, DIMENSION(:), ALLOCATABLE :: randomindex
      REAL(rprec) :: r_temp, phi_temp, z_temp, vll_temp, w_temp, &
                     za_temp, m_temp, c_temp, mu_temp, t_temp
      REAL(rprec), DIMENSION(3) :: v_temp
      REAL, DIMENSION(:), ALLOCATABLE :: randomnumbers

!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      IF (myworkid == master) THEN
         ALLOCATE(randomindex(nparticles))
         ALLOCATE(randomnumbers(nparticles))
         CALL RANDOM_NUMBER(randomnumbers)
         randomindex = NINT(randomnumbers*nparticles)
         WHERE(randomindex < 1) randomindex = 1
         DEALLOCATE(randomnumbers)
         DO i = 1, nparticles
            d1 = i
            d2 = randomindex(i)
            ! Save swap
            r_temp = R_start(d2)
            phi_temp = phi_start(d2)
            z_temp   = Z_start(d2)
            vll_temp = vll_start(d2)
            w_temp   = weight(d2)
            za_temp  = Zatom(d2)
            m_temp   = mass(d2)
            c_temp   = charge(d2)
            mu_temp  = mu_start(d2)
            t_temp   = t_end(d2)
            v_temp   = v_neut(:,d2)
            b_temp   = Beam(d2)
            !swap
            R_start(d2)   = R_start(d1)
            phi_start(d2) = phi_start(d1)
            Z_start(d2)   = Z_start(d1)
            vll_start(d2) = vll_start(d1)
            weight(d2)    = weight(d1)
            Zatom(d2)     = Zatom(d1)
            mass(d2)      = mass(d1)
            charge(d2)    = charge(d1)
            mu_start(d2)  = mu_start(d1)
            t_end(d2)     = t_end(d1)
            v_neut(:,d2)  = v_neut(:,d1)
            Beam(d2)      = Beam(d1)
            ! Finish
            R_start(d1)   = r_temp
            phi_start(d1) = phi_temp
            Z_start(d1)   = z_temp
            vll_start(d1) = vll_temp
            weight(d1)    = w_temp
            Zatom(d1)     = za_temp
            mass(d1)      = m_temp
            charge(d1)    = c_temp
            mu_start(d1)  = mu_temp
            t_end(d1)     = t_temp
            v_neut(:,d1)  = v_temp
            Beam(d1)      = b_temp
         END DO
         DEALLOCATE(randomindex)
      END IF

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(R_start,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(phi_start,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(Z_start,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(vll_start,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(weight,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(Zatom,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(mass,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(charge,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(mu_start,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(t_end,nparticles,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(v_neut,nparticles*3,MPI_REAL8, master, MPI_COMM_BEAMS,ierr_mpi)
      CALL MPI_BCAST(beam,nparticles,MPI_INTEGER, master, MPI_COMM_BEAMS,ierr_mpi)
#endif

      RETURN
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE beams3d_randomize_particles

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20080911 21:39:23 nbsp Welcome to the Wiki Howdy Brendan, and Welcome to the Wiki!  My names Evan, pleased to meet you!  Thanks for the recent automotive reviews.  If you can add anything to any of the entries, dont be afraid to edit the directly... everything on the wiki was written by people like you. Users/JabberWokky Evan JabberWokky Edwards