      subroutine pest_soil_tot

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates the total amount of pesticide in the soil

      use pesticide_data_module
      use hru_module, only : hru, ihru
      use soil_module
      use constituent_mass_module
      use output_ls_pesticide_module
      use plant_module

      implicit none

      integer :: j        !none          |HRU number
      integer :: k        !none          |sequential pesticide number
      integer :: ly       !none          |soil layer

      j = ihru

      if (cs_db%num_pests == 0) return

      do k = 1, cs_db%num_pests
        hpestb_d(j)%pest(k)%plant = cs_pl(j)%pest(k)    !store for each plant in future
        hpestb_d(j)%pest(k)%soil = 0.
        do ly = 1, soil(j)%nly
          hpestb_d(j)%pest(k)%soil = hpestb_d(j)%pest(k)%soil + cs_soil(j)%ly(ly)%pest(k)
        end do
      end do

      return
      end subroutine pest_soil_tot


text:
MODULE neo_parameters
  USE neo_precision
  USE nrtype
END MODULE neo_parameters

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
C     To debug sc_bounded_normalization

      program bounded07
      integer a, b

      if (a.gt.0.and.b.gt.0) then
         if (a+b.gt.0) then
            print *, 'test redondant a+b>0'
            print *, 'test redondant a+b>0'
         endif
      endif

      end

ext: f


text:
      subroutine dcopy_wrap(x,y,nx)
cf2py intent(inplace) y
cf2py intent(hide) nx
cf2py threadsafe
      double precision x(nx), y(nx)
      external dcopy
!       dcopy(n,dx,incx,dy,incy)

      CALL dcopy(nx,x,1,y,1)
      RETURN
      END

      SUBROUTINE dtrmm_wrap(M,N,A,B,SIDE,TRANSA,UPLO,ALPHA)
cf2py intent(inplace)::B
cf2py integer intent(hide)::M
cf2py integer intent(hide)::N
cf2py threadsafe

*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(M,M),B(M,N)

      EXTERNAL DTRMM
! DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*     X = alpha*op( A )*B,   or   X = alpha*B*op( A ),

      DIAG = 'N'
      LDA = M
      LDB = M

      CALL DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
      RETURN
      END

      SUBROUTINE dtrsm_wrap(M,N,A,B,SIDE,TRANSA,UPLO,ALPHA)
cf2py intent(inplace)::B
cf2py integer intent(hide)::M
cf2py integer intent(hide)::N
cf2py threadsafe

*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(M,M),B(M,N)

      EXTERNAL DTRSM
! DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,

      DIAG = 'N'
      LDA = M
      LDB = M

      CALL DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
      RETURN
      END

      SUBROUTINE dpotrf_wrap(A,n,info)

cf2py double precision dimension(n,n), intent(inplace) :: A
cf2py integer depend(A),intent(hide)::n=shape(A,0)
cf2py integer intent(out)::info
cf2py threadsafe

      DOUBLE PRECISION A(n,n)
      INTEGER n, info, i, j

      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization

!     C <- cholesky(C)
      call DPOTRF( 'U', n, A, n, info )
      do i=2,n
        do j=1,i-1
          A(i,j)=0.0D0
        enddo
      enddo

      return
      END

      subroutine dchdc_wrap(a,p,work,piv,info)

cf2py intent(inplace) a
cf2py intent(hide) p
cf2py intent(hide) work
cf2py intent(out) piv
cf2py intent(out) info
cf2py threadsafe

      integer lda,p,piv(p),job,info
      double precision a(p,p),work(p)

c
c     dchdc computes the cholesky decomposition of a positive definite
c     matrix.  a pivoting option allows the user to estimate the
c     condition of a positive definite matrix or determine the rank
c     of a positive semidefinite matrix.
c
c     on entry
c
c         a      double precision(lda,p).
c                a contains the matrix whose decomposition is to
c                be computed.  onlt the upper half of a need be stored.
c                the lower part of the array a is not referenced.
c
c         lda    integer.
c                lda is the leading dimension of the array a.
c
c         Since f2py will make sure that a is column-major contiguous,
c         lda is set to p in this wrapper.
c
c         p      integer.
c                p is the order of the matrix.
c
c         work   double precision.
c                work is a work array.
c
c         piv   integer(p).
c                piv contains integers that control the selection
c                of the pivot elements, if pivoting has been requested.
c                each diagonal element a(k,k)
c                is placed in one of three classes according to the
c                value of piv(k).
c
c                   if piv(k) .gt. 0, then x(k) is an initial
c                                      element.
c
c                   if piv(k) .eq. 0, then x(k) is a free element.
c
c                   if piv(k) .lt. 0, then x(k) is a final element.
c
c                before the decomposition is computed, initial elements
c                are moved by symmetric row and column interchanges to
c                the beginning of the array a and final
c                elements to the end.  both initial and final elements
c                are frozen in place during the computation and only
c                free elements are moved.  at the k-th stage of the
c                reduction, if a(k,k) is occupied by a free element
c                it is interchanged with the largest free element
c                a(l,l) with l .ge. k.  piv is not referenced if
c                job .eq. 0.
c
c        job     integer.
c                job is an integer that initiates column pivoting.
c                if job .eq. 0, no pivoting is done.
c                if job .ne. 0, pivoting is done.
c
c        In this wrapper, job=1.
c
c     on return
c
c         a      a contains in its upper half the cholesky factor
c                of the matrix a as it has been permuted by pivoting.
c
c         piv   piv(j) contains the index of the diagonal element
c                of a that was moved into the j-th position,
c                provided pivoting was requested.
c
c         info   contains the index of the last positive diagonal
c                element of the cholesky factor.
c
c     for positive definite matrices info = p is the normal return.
c     for pivoting with positive semidefinite matrices info will
c     in general be less than p.  however, info may be greater than
c     the rank of a, since rounding error can cause an otherwise zero
c     element to be positive. indefinite systems will always cause
c     info to be less than p.
c
c     linpack. this version dated 08/14/78 .
c     j.j. dongarra and g.w. stewart, argonne national laboratory and
c     university of maryland.
c
c
c     blas daxpy,dswap
c     fortran dsqrt
c
c     internal variables
c
      integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl
      double precision temp
      double precision maxdia
      logical swapk,negk
c
      pl = 1
      pu = 0
      info = p
      job=1
      lda=p
      if (job .eq. 0) go to 160
c
c        pivoting has been requested. rearrange the
c        the elements according to piv.
c
         do 70 k = 1, p
            swapk = piv(k) .gt. 0
            negk = piv(k) .lt. 0
            piv(k) = k
            if (negk) piv(k) = -piv(k)
            if (.not.swapk) go to 60
               if (k .eq. pl) go to 50
                  call dswap(pl-1,a(1,k),1,a(1,pl),1)
                  temp = a(k,k)
                  a(k,k) = a(pl,pl)
                  a(pl,pl) = temp
                  plp1 = pl + 1
                  if (p .lt. plp1) go to 40
                  do 30 j = plp1, p
                     if (j .ge. k) go to 10
                        temp = a(pl,j)
                        a(pl,j) = a(j,k)
                        a(j,k) = temp
                     go to 20
   10                continue
                     if (j .eq. k) go to 20
                        temp = a(k,j)
                        a(k,j) = a(pl,j)
                        a(pl,j) = temp
   20                continue
   30             continue
   40             continue
                  piv(k) = piv(pl)
                  piv(pl) = k
   50          continue
               pl = pl + 1
   60       continue
   70    continue
         pu = p
         if (p .lt. pl) go to 150
         do 140 kb = pl, p
            k = p - kb + pl
            if (piv(k) .ge. 0) go to 130
               piv(k) = -piv(k)
               if (pu .eq. k) go to 120
                  call dswap(k-1,a(1,k),1,a(1,pu),1)
                  temp = a(k,k)
                  a(k,k) = a(pu,pu)
                  a(pu,pu) = temp
                  kp1 = k + 1
                  if (p .lt. kp1) go to 110
                  do 100 j = kp1, p
                     if (j .ge. pu) go to 80
                        temp = a(k,j)
                        a(k,j) = a(j,pu)
                        a(j,pu) = temp
                     go to 90
   80                continue
                     if (j .eq. pu) go to 90
                        temp = a(k,j)
                        a(k,j) = a(pu,j)
                        a(pu,j) = temp
   90                continue
  100             continue
  110             continue
                  jt = piv(k)
                  piv(k) = piv(pu)
                  piv(pu) = jt
  120          continue
               pu = pu - 1
  130       continue
  140    continue
  150    continue
  160 continue
      do 270 k = 1, p
c
c        reduction loop.
c
         maxdia = a(k,k)
         kp1 = k + 1
         maxl = k
c
c        determine the pivot element.
c
         if (k .lt. pl .or. k .ge. pu) go to 190
            do 180 l = kp1, pu
               if (a(l,l) .le. maxdia) go to 170
                  maxdia = a(l,l)
                  maxl = l
  170          continue
  180       continue
  190    continue
c
c        quit if the pivot element is not positive.
c
         if (maxdia .gt. 0.0d0) go to 200
            info = k - 1
c     ......exit
            go to 280
  200    continue
         if (k .eq. maxl) go to 210
c
c           start the pivoting and update piv.
c
            km1 = k - 1
            call dswap(km1,a(1,k),1,a(1,maxl),1)
            a(maxl,maxl) = a(k,k)
            a(k,k) = maxdia
            jp = piv(maxl)
            piv(maxl) = piv(k)
            piv(k) = jp
  210    continue
c
c        reduction step. pivoting is contained across the rows.
c
         work(k) = dsqrt(a(k,k))
         a(k,k) = work(k)
         if (p .lt. kp1) go to 260
         do 250 j = kp1, p
            if (k .eq. maxl) go to 240
               if (j .ge. maxl) go to 220
                  temp = a(k,j)
                  a(k,j) = a(j,maxl)
                  a(j,maxl) = temp
               go to 230
  220          continue
               if (j .eq. maxl) go to 230
                  temp = a(k,j)
                  a(k,j) = a(maxl,j)
                  a(maxl,j) = temp
  230          continue
  240       continue
            a(k,j) = a(k,j)/work(k)
            work(j) = a(k,j)
            temp = -a(k,j)
            call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1)
  250    continue
  260    continue
  270 continue
  280 continue

C Zero the lower triangle
      do i=2,info
        do j=1,i-1
          a(i,j) = 0.0D0
        enddo
      enddo

C decrement the pivots.
      do i=1,p
        piv(i) = piv(i)-1
      enddo
      return
      end

      SUBROUTINE dpotrs_wrap(chol_fac, b, info, n, m, uplo)

cf2py double precision dimension(n,n), intent(in) :: chol_fac
cf2py integer depend(chol_fac), intent(hide)::n=shape(chol_fac,0)
cf2py integer depend(b), intent(hide)::m=shape(b,1)
cf2py optional character intent(in):: uplo='U'
cf2py integer intent(out)::info
cf2py double precision intent(inplace), dimension(n,m)::b
cf2py threadsafe

      DOUBLE PRECISION chol_fac(n,n), b(n,m)
      INTEGER n, info,m

      CHARACTER uplo

      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system

      call DPOTRS(uplo,n,m,chol_fac,n,b,n,info)

      return
      END

max_stars_repo_name: kyleabeauchamp/pymc

max_issues_repo_name: kyleabeauchamp/pymc

max_forks_repo_name: kyleabeauchamp/pymc


!
!  Copyright 2019-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.
!
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
module force_sub
  implicit none

contains

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

  subroutine calc_force(system,pp,fg,info,mg,stencil,poisson,srg,ppg,tpsi,ewald)
    use structures
    use math_constants,only : zi,pi
    use sendrecv_grid, only: s_sendrecv_grid, update_overlap_real8, update_overlap_complex8, dealloc_cache
    use communication, only: comm_summation
    use nonlocal_potential, only: calc_uVpsi_rdivided, calc_uVpsi
    use sym_vector_sub, only: sym_vector_xyz
    use sym_sub, only: use_symmetry
    use pseudo_pt_so_sub, only: SPIN_ORBIT_ON, calc_uVpsi_so
    use plusU_global, only: PLUS_U_ON, dm_mms_nla, U_eff
    use salmon_global, only: kion,cutoff_g,yn_periodic
    use code_optimization, only: force_omp_mode
    use timer
    implicit none
    type(s_dft_system)      ,intent(inout) :: system
    type(s_pp_info)         ,intent(in)    :: pp
    type(s_reciprocal_grid) ,intent(in)    :: fg
    type(s_parallel_info)   ,intent(in)    :: info
    type(s_rgrid)           ,intent(in)    :: mg
    type(s_stencil)         ,intent(in)    :: stencil
    type(s_poisson)         ,intent(in)    :: poisson
    type(s_sendrecv_grid)   ,intent(inout) :: srg
    type(s_pp_grid)         ,intent(in)    :: ppg
    type(s_orbital)         ,intent(inout) :: tpsi
    type(s_ewald_ion_ion)   ,intent(in)    :: ewald
    !
    integer :: ix,iy,iz,ia,nion,im,Nspin,ik_s,ik_e,io_s,io_e,nlma,ik,io,ispin,ilma,j
    integer :: m1,m2,jlma,n,l,Nproj_pairs,iprj,Nlma_ao
    real(8) :: kAc(3), rtmp,rtmp2(3),g(3),r(3),G2,Gd
    complex(8) :: rho_i, rho_e, egd, VG
    real(8),allocatable :: F_tmp(:,:),F_sum(:,:)
    real(8),allocatable :: dden(:,:,:,:)
    complex(8) :: w(3),duVpsi(3)
    complex(8),allocatable :: gtpsi(:,:,:,:),uVpsibox(:,:,:,:,:),uVpsibox2(:,:,:,:,:)
    complex(8),allocatable :: phipsibox(:,:),phipsibox2(:,:)
    complex(8),allocatable :: dphipsi_lma(:,:)
    complex(8) :: ddm_mms_nla(3), phipsi, dphipsi(3)
    complex(8),parameter :: zero=(0.0d0,0.0d0)
    complex(8),allocatable :: zF_tmp(:,:)
    complex(8) :: ztmp
    integer :: Norb,iorb,ilocal

    call timer_begin(LOG_CALC_ION_FORCE)

    nion = system%nion
    if(.not.allocated(system%Force)) allocate(system%Force(3,nion))
    allocate( F_tmp(3,nion), F_sum(3,nion) )
    if( PLUS_U_ON ) then
      allocate( zF_tmp(3,nion) )
      zF_tmp=zero
    end if

    if(info%im_s/=1 .or. info%im_e/=1) stop "error: calc_force_periodic" !??????
    im = 1

    Nspin = system%Nspin
    ik_s = info%ik_s
    ik_e = info%ik_e
    io_s = info%io_s
    io_e = info%io_e
    Norb = system%Nspin*info%numo*info%numk

    Nlma = ppg%Nlma

  ! Ewald sum
    call timer_begin(LOG_CALC_FORCE_ION_ION)
    call force_ewald_rspace(system%Force,F_tmp,system,info,ewald,pp,nion,info%icomm_r)
    call timer_end(LOG_CALC_FORCE_ION_ION)

    call timer_begin(LOG_CALC_FORCE_FOURIER)
    F_tmp   = 0d0

    select case(yn_periodic)
    case('n')
      allocate( dden(3,mg%is_array(1):mg%ie_array(1) &
                      ,mg%is_array(2):mg%ie_array(2) &
                      ,mg%is_array(3):mg%ie_array(3)))
      dden = 0d0
    case('y')

    allocate(dden(1,1,1,1)); dden=0d0 !! dummy allocation to avoid crush: not use in isolated system

    ! Fourier part (local part, etc)
      !$omp parallel do collapse(2) private(ix,iy,iz,ia,r,g,G2,Gd,rho_i,rho_e,rtmp,egd,VG) reduction(+:F_tmp)
      do iz=mg%is(3),mg%ie(3)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        if(fg%if_Gzero(ix,iy,iz)) cycle
        g(1) = fg%vec_G(1,ix,iy,iz)
        g(2) = fg%vec_G(2,ix,iy,iz)
        g(3) = fg%vec_G(3,ix,iy,iz)
        G2 = sum(g(:)**2)
        if(G2 .gt. cutoff_g**2) cycle   !xxx

        rho_i = ppg%zrhoG_ion(ix,iy,iz)
        rho_e = poisson%zrhoG_ele(ix,iy,iz)

        !OCL swp
        do ia=info%ia_s,info%ia_e
          r = system%Rion(1:3,ia)
          Gd = sum(g(:)*r(:))
          egd = exp(zI*Gd)
          rtmp = pp%Zps(Kion(ia))* fg%coef(ix,iy,iz) * fg%exp_ewald(ix,iy,iz)
          VG = ppg%zVG_ion(ix,iy,iz,kion(ia)) - fg%coef(ix,iy,iz) * pp%zps(kion(ia))
          F_tmp(:,ia) = F_tmp(:,ia) + g(:)* ( rtmp * aimag(rho_i*egd) + aimag(egd*rho_e*conjg(VG)) )
        end do
      end do
      end do
      end do
      !$omp end parallel do
    end select
    call timer_end(LOG_CALC_FORCE_FOURIER)

    call timer_begin(LOG_CALC_FORCE_ELEC_ION)
    if(allocated(tpsi%rwf)) then
      allocate(tpsi%zwf(mg%is_array(1):mg%ie_array(1) &
                       ,mg%is_array(2):mg%ie_array(2) &
                       ,mg%is_array(3):mg%ie_array(3) &
                       ,nspin,info%io_s:info%io_e,info%ik_s:info%ik_e,info%im_s:info%im_e))
      tpsi%zwf = cmplx(tpsi%rwf)
    end if

    if( SPIN_ORBIT_ON )then
      call calc_uVpsi_so(nspin,info,ppg,tpsi,uVpsibox2)
    else
    ! uVpsibox2 = < uV | exp(ikr) | psi >
      if (info%if_divide_rspace) then
        call calc_uVpsi_rdivided(nspin,info,ppg,tpsi,uVpsibox,uVpsibox2)
      else
        call calc_uVpsi(nspin,info,ppg,tpsi,uVpsibox2)
      end if
    end if

    if( PLUS_U_ON )then
      Nlma_ao = size(ppg%ia_tbl_ao)
      allocate( phipsibox(Nlma_ao,Norb)  ); phipsibox=zero
      allocate( phipsibox2(Nlma_ao,Norb) ); phipsibox2=zero
      iorb = 0
      do ik=ik_s,ik_e
      do io=io_s,io_e
      do ispin=1,Nspin
        iorb = iorb + 1
        do ilma=1,Nlma_ao
          ia = ppg%ia_tbl_ao(ilma)
          phipsi = 0.0d0
          do j=1,ppg%mps_ao(ia)
            ix = ppg%jxyz_ao(1,j,ia)
            iy = ppg%jxyz_ao(2,j,ia)
            iz = ppg%jxyz_ao(3,j,ia)
            phipsi = phipsi + conjg(ppg%zekr_phi_ao(j,ilma,ik)) &
                            * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
          end do
          phipsi = phipsi * system%Hvol
          phipsibox(ilma,iorb) = phipsi
        end do
      end do
      end do
      end do
      call comm_summation(phipsibox,phipsibox2,Nlma_ao*Norb,info%icomm_r)
    end if
    call timer_end(LOG_CALC_FORCE_ELEC_ION)

    if(info%if_divide_rspace) then
       if(tpsi%update_zwf_overlap) call update_overlap_complex8(srg,mg,tpsi%zwf)
    end if

    call timer_begin(LOG_CALC_FORCE_ELEC_ION)
    kAc = 0d0

!$omp parallel default(none) &
!$omp   firstprivate(kAc) &
!$omp   private(ik,io,ispin,rtmp,rtmp2,ilocal,ilma,ia,duVpsi,j,ix,iy,iz,w,dphipsi_lma,dphipsi,Nproj_pairs) &
!$omp   private(jlma,l,n,m1,m2,ddm_mms_nla,gtpsi) &
!$omp   shared(im,ik_s,ik_e,io_s,io_e,nspin,tpsi,mg,stencil,system,ppg,uVpsibox2,yn_periodic) &
!$omp   shared(PLUS_U_ON,Nlma_ao,phipsibox2,iorb,zF_tmp,U_eff,dm_mms_nla) &
!$omp   shared(SPIN_ORBIT_ON,ztmp) &
!$omp   shared(dden) &
!$omp   reduction(+:F_tmp) &
!$omp   if(force_omp_mode)

    if (.not. allocated(gtpsi)) then
      allocate(gtpsi(3,mg%is_array(1):mg%ie_array(1) &
                      ,mg%is_array(2):mg%ie_array(2) &
                      ,mg%is_array(3):mg%ie_array(3)))
    end if

!$omp do collapse(2) reduction(+:dden)
    do ik=ik_s,ik_e
    do io=io_s,io_e
    do ispin=1,Nspin

       ! gtpsi = (nabla) psi
       !call timer_begin(LOG_CALC_FORCE_GTPSI)
       call calc_gradient_psi(tpsi%zwf(:,:,:,ispin,io,ik,im),gtpsi,mg%is_array,mg%ie_array,mg%is,mg%ie &
            ,mg%idx,mg%idy,mg%idz,stencil%coef_nab,system%rmatrix_B)
       !call timer_end(LOG_CALC_FORCE_GTPSI)

       if(yn_periodic=='n') then
         rtmp = 2d0 * system%rocc(io,ik,ispin) * system%wtk(ik) * system%Hvol
!$omp parallel do collapse(2) private(iz,iy,ix,w,rtmp2)
         do iz=mg%is(3),mg%ie(3)
         do iy=mg%is(2),mg%ie(2)
         do ix=mg%is(1),mg%ie(1)
           w = conjg(gtpsi(:,ix,iy,iz)) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
           rtmp2(:) = rtmp * dble(w(:))
           dden(:,ix,iy,iz) = dden(:,ix,iy,iz) + rtmp2(:)
         enddo
         enddo
         enddo
!$omp end parallel do
       end if

       ! nonlocal part
       !call timer_begin(LOG_CALC_FORCE_NONLOCAL)
       if(yn_periodic=='y') kAc(1:3) = system%vec_k(1:3,ik) + system%vec_Ac(1:3)
       rtmp = 2d0 * system%rocc(io,ik,ispin) * system%wtk(ik) * system%Hvol

       if( SPIN_ORBIT_ON )then
         do ilma=1,size(ppg%ia_tbl_so)
           ia = ppg%ia_tbl_so(ilma)
           duVpsi=zero
           do j=1,ppg%mps(ia)
             ix = ppg%jxyz(1,j,ia)
             iy = ppg%jxyz(2,j,ia)
             iz = ppg%jxyz(3,j,ia)
             w(1) = gtpsi(1,ix,iy,iz) + zI*kAc(1) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             w(2) = gtpsi(2,ix,iy,iz) + zI*kAc(2) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             w(3) = gtpsi(3,ix,iy,iz) + zI*kAc(3) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             duVpsi(1) = duVpsi(1) + conjg( ppg%zekr_uV_so(j,ilma,ik,ispin,1) )*w(1)
             duVpsi(2) = duVpsi(2) + conjg( ppg%zekr_uV_so(j,ilma,ik,ispin,1) )*w(2)
             duVpsi(3) = duVpsi(3) + conjg( ppg%zekr_uV_so(j,ilma,ik,ispin,1) )*w(3)
           end do
           ztmp=uVpsibox2(1,ilma,io,ik,im)+uVpsibox2(2,ilma,io,ik,im)
           F_tmp(1,ia) = F_tmp(1,ia) - rtmp*dble( conjg(duVpsi(1)) * ztmp )
           F_tmp(2,ia) = F_tmp(2,ia) - rtmp*dble( conjg(duVpsi(2)) * ztmp )
           F_tmp(3,ia) = F_tmp(3,ia) - rtmp*dble( conjg(duVpsi(3)) * ztmp )
         end do
       else
!$omp parallel do private(ilocal,ilma,ia,duVpsi,j,ix,iy,iz,w) reduction(+:F_tmp)
       do ilocal=1,ppg%ilocal_nlma
          ilma=ppg%ilocal_nlma2ilma(ilocal)
          ia  =ppg%ilocal_nlma2ia  (ilocal)
          duVpsi = 0d0
!OCL swp
!OCL swp_freg_rate(115)
          do j=1,ppg%mps(ia)
             ix = ppg%jxyz(1,j,ia)
             iy = ppg%jxyz(2,j,ia)
             iz = ppg%jxyz(3,j,ia)
             w(1) = gtpsi(1,ix,iy,iz) + zI* kAc(1) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             w(2) = gtpsi(2,ix,iy,iz) + zI* kAc(2) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             w(3) = gtpsi(3,ix,iy,iz) + zI* kAc(3) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
             duVpsi(1) = duVpsi(1) + conjg(ppg%zekr_uV(j,ilma,ik)) * w(1) ! < uV | exp(ikr) (nabla) | psi >
             duVpsi(2) = duVpsi(2) + conjg(ppg%zekr_uV(j,ilma,ik)) * w(2) ! < uV | exp(ikr) (nabla) | psi >
             duVpsi(3) = duVpsi(3) + conjg(ppg%zekr_uV(j,ilma,ik)) * w(3) ! < uV | exp(ikr) (nabla) | psi >
          end do
          F_tmp(1,ia) = F_tmp(1,ia) - rtmp * dble( conjg(duVpsi(1)) * uVpsibox2(ispin,io,ik,im,ilma) )
          F_tmp(2,ia) = F_tmp(2,ia) - rtmp * dble( conjg(duVpsi(2)) * uVpsibox2(ispin,io,ik,im,ilma) )
          F_tmp(3,ia) = F_tmp(3,ia) - rtmp * dble( conjg(duVpsi(3)) * uVpsibox2(ispin,io,ik,im,ilma) )
       end do
!$omp end parallel do
       end if
       !call timer_end(LOG_CALC_FORCE_NONLOCAL)

       if( PLUS_U_ON )then
          if( .not.allocated(dphipsi_lma) )then
             allocate( dphipsi_lma(3,Nlma_ao) ); dphipsi_lma=zero
          end if
          do ilma=1,Nlma_ao
             ia = ppg%ia_tbl_ao(ilma)
             dphipsi = zero
             do j=1,ppg%mps_ao(ia)
                ix = ppg%jxyz_ao(1,j,ia)
                iy = ppg%jxyz_ao(2,j,ia)
                iz = ppg%jxyz_ao(3,j,ia)
                w  = gtpsi(:,ix,iy,iz) + zI* kAc(:) * tpsi%zwf(ix,iy,iz,ispin,io,ik,im)
                dphipsi(:) = dphipsi(:) + conjg(ppg%zekr_phi_ao(j,ilma,ik)) * w(:)
             end do
             dphipsi_lma(:,ilma) = dphipsi(:) * system%Hvol
          end do
       end if
       if( PLUS_U_ON )then
          Nproj_pairs = size(ppg%proj_pairs_ao,2)
          do iprj=1,Nproj_pairs
             ilma=ppg%proj_pairs_ao(1,iprj)
             jlma=ppg%proj_pairs_ao(2,iprj)
             ia = ppg%proj_pairs_info_ao(1,iprj)
             l  = ppg%proj_pairs_info_ao(2,iprj)
             n  = ppg%proj_pairs_info_ao(3,iprj)
             m1 = ppg%proj_pairs_info_ao(4,iprj)
             m2 = ppg%proj_pairs_info_ao(5,iprj)
             ddm_mms_nla(:)= & ! ddm_mms_nla(m1,m2,ispin,n,l,ia)=
                  system%rocc(io,ispin,ik)*system%wtk(ik) &
                  *( dphipsi_lma(:,ilma)*conjg(phipsibox2(jlma,iorb)) &
                  + phipsibox2(ilma,iorb)*conjg(dphipsi_lma(:,jlma)) )
             if( m1 == m2 )then
                zF_tmp(:,ia) = zF_tmp(:,ia) &
                     - 0.5d0*U_eff(n,l,ia)*( 1.0d0 - 2.0d0*dm_mms_nla(m1,m2,ispin,n,l,ia) ) &
                     * ddm_mms_nla(:)
             else
                zF_tmp(:,ia) = zF_tmp(:,ia) &
                     - 0.5d0*U_eff(n,l,ia)*( -2.0d0*dm_mms_nla(m1,m2,ispin,n,l,ia) ) &
                     * ddm_mms_nla(:)
             end if
          end do !iprj
       end if

    end do !ispin
    end do !io
    end do !ik
!$omp end do

    if (allocated(gtpsi)) deallocate(gtpsi)

!$omp end parallel
    call timer_end(LOG_CALC_FORCE_ELEC_ION)

    if(yn_periodic=='n') then
    ! local part (based on density gradient)
!$omp parallel do private(iz,iy,ix,ia)
      do ia=1,nion
        do iz=mg%is(3),mg%ie(3)
        do iy=mg%is(2),mg%ie(2)
        do ix=mg%is(1),mg%ie(1)
          F_tmp(:,ia) = F_tmp(:,ia) - dden(:,ix,iy,iz) * ppg%Vpsl_ion(ix,iy,iz,ia)
        end do
        end do
        end do
      end do
!$omp end parallel do
      deallocate(dden)
    end if

    !do ia=1,nion
    !  write(*,'(1x,i4,2f20.10)') ia,real(zF_tmp(1,ia)),aimag(zF_tmp(1,ia))
    !  write(*,'(1x,4x,2f20.10)')    real(zF_tmp(2,ia)),aimag(zF_tmp(2,ia))
    !  write(*,'(1x,4x,2f20.10)')    real(zF_tmp(3,ia)),aimag(zF_tmp(3,ia))
    !end do

    call comm_summation(F_tmp,F_sum,3*nion,info%icomm_rko)

!$omp parallel do private(ia)
    do ia=1,nion
      system%Force(:,ia) = system%Force(:,ia) + F_sum(:,ia)
    end do
!$omp end parallel do
!
    if (use_symmetry) then
      do ia=1,nion
        call sym_vector_xyz( system%Force(:,ia) )
      end do
    end if

    if(allocated(tpsi%rwf)) deallocate(tpsi%zwf)
    if(allocated(gtpsi)) deallocate(gtpsi)
    if(allocated(uVpsibox)) deallocate(uVpsibox)
    deallocate(F_tmp,F_sum,uVpsibox2)
    if( PLUS_U_ON )then
      deallocate( phipsibox, phipsibox2 )
      deallocate( dphipsi_lma )
      deallocate( zF_tmp )
    end if

    call timer_end(LOG_CALC_ION_FORCE)
    return
  end subroutine calc_force

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

  subroutine force_ewald_rspace(F_sum,F_tmp,system,info,ewald,pp,nion,comm)
    use structures
    use math_constants,only : pi,zi
    use salmon_math
    use salmon_global, only: yn_periodic,kion,aEwald,cutoff_r
    use communication, only: comm_summation
    implicit none
    type(s_dft_system),intent(in) :: system
    type(s_parallel_info),intent(in) :: info
    type(s_ewald_ion_ion),intent(in) :: ewald
    type(s_pp_info)   ,intent(in) :: pp
    integer  ,intent(in) :: nion,comm
    real(8)              :: F_tmp(3,nion),F_tmp_l(3,nion),F_sum(3,nion)
    !
    integer :: ix,iy,iz,iia,ia,ib,ipair
    real(8) :: rr,rab(3),r(3)

    select case(yn_periodic)
    case('n')

      F_sum = 0d0
      do ia=1,nion
        do ib=1,nion
          if(ia==ib) cycle
          rab(:) = system%Rion(:,ia) - system%Rion(:,ib)
          rr = sum(rab(:)**2)
          F_sum(:,ia) = F_sum(:,ia) + pp%Zps(Kion(ia))*pp%Zps(Kion(ib)) * rab(:)/sqrt(rr)**3
        end do
      end do

    case('y')

      F_sum = 0d0
      F_tmp = 0d0
      if(ewald%yn_bookkeep=='y') then

         F_tmp_l = 0d0
!$omp parallel do private(iia,ia,ipair,ix,iy,iz,ib,r,rab,rr)
         do iia= 1,info%nion_mg
        !do ia= system%nion_s, system%nion_e
        !do ia=1,nion
            ia = info%ia_mg(iia)
            do ipair = 1,ewald%npair_bk(iia)
               ix = ewald%bk(1,ipair,iia)
               iy = ewald%bk(2,ipair,iia)
               iz = ewald%bk(3,ipair,iia)
               ib = ewald%bk(4,ipair,iia)
              !if (ix**2+iy**2+iz**2 == 0 .and. ia == ib) cycle
               r(1) = ix*system%primitive_a(1,1) &
                    + iy*system%primitive_a(1,2) &
                    + iz*system%primitive_a(1,3)
               r(2) = ix*system%primitive_a(2,1) &
                    + iy*system%primitive_a(2,2) &
                    + iz*system%primitive_a(2,3)
               r(3) = ix*system%primitive_a(3,1) &
                    + iy*system%primitive_a(3,2) &
                    + iz*system%primitive_a(3,3)
               rab(1) = system%Rion(1,ia)-r(1) - system%Rion(1,ib)
               rab(2) = system%Rion(2,ia)-r(2) - system%Rion(2,ib)
               rab(3) = system%Rion(3,ia)-r(3) - system%Rion(3,ib)
               rr = sum(rab(:)**2)
               if(rr .gt. cutoff_r**2) then
                 !!!
               else
                 F_tmp_l(:,ia) = F_tmp_l(:,ia)  &
                               - pp%Zps(Kion(ia))*pp%Zps(Kion(ib))*rab(:)/sqrt(rr)*(-erfc_salmon(sqrt(aEwald*rr))/rr &
                               -2*sqrt(aEwald/(rr*Pi))*exp(-aEwald*rr))
               end if

            end do  !ipair
         end do     !ia
!$omp end parallel do
         call comm_summation(F_tmp_l,F_tmp,3*nion,comm)

      else

         stop  "must use book-keeping method for periodic condition"

      endif
      F_sum = F_sum + F_tmp

    end select
  end subroutine force_ewald_rspace
end module force_sub
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130


text:
c-----------------------------------------------------------------------
      program nekbone

      include 'SIZE'
      include 'TOTAL'
      include 'mpif.h'

      parameter (lxyz = lx1*ly1*lz1)
      parameter (lt=lxyz*lelt)

      real ah(lx1*lx1),bh(lx1),ch(lx1*lx1),dh(lx1*lx1)
     $    ,zpts(2*lx1),wght(2*lx1)

      real x(lt),f(lt),r(lt),w(lt),p(lt),z(lt),c(lt)
      real g(6,lt)
      real mfloplist(1024), avmflop
      integer icount

      logical ifbrick
      integer iel0,ielN,ielD   ! element range per proc.
      integer nx0,nxN,nxD      ! poly. order range
      integer npx,npy,npz      ! poly. order range
      integer mx,my,mz         ! poly. order range
      integer numthreads, omp_get_max_threads

      call iniproc(mpi_comm_world)    ! has nekmpi common block
      call read_param(ifbrick,iel0,ielN,ielD,nx0,nxN,nxD,
     +                npx,npy,npz,mx,my,mz)

      numthreads = 1
#ifdef _OPENMP
      numthreads= omp_get_max_threads()
#endif

      if (nid.eq.0) then
        write(*,*) "Max number of threads: ", numthreads
      end if

c     GET PLATFORM CHARACTERISTICS
c     iverbose = 1
c     call platform_timer(iverbose)   ! iverbose=0 or 1

      icount = 0

c     SET UP and RUN NEKBONE
      do nx1=nx0,nxN,nxD
         call init_dim
         do nelt=iel0,ielN,ielD
           call init_mesh(ifbrick,npx,npy,npz,mx,my,mz)
           call proxy_setupds    (gsh)     ! Has nekmpi common block
           call set_multiplicity (c)       ! Inverse of counting matrix

           call proxy_setup(ah,bh,ch,dh,zpts,wght,g)

           niter = 100
           n     = nx1*ny1*nz1*nelt

           call set_f(f,c,n)
           call cg(x,f,g,c,r,w,p,z,n,niter,flop_cg)

           call nekgsync()

           call set_timer_flop_cnt(0)
           call cg(x,f,g,c,r,w,p,z,n,niter,flop_cg)
           call set_timer_flop_cnt(1)

           call gs_free(gsh)

           icount = icount + 1
           mfloplist(icount)= mflops*np
         enddo
      enddo

      avmflop = 0.0
      do i = 1, icount
        avmflop = avmflop + mfloplist(i)
      end do

      if (icount .ne. 0) then
        avmflop = avmflop/icount
      end if

      if (nid .eq. 0) then
        write(6,1) avmflop
      end if
    1 format('Av MFlops = ', 1pe12.4)

c     TEST BANDWIDTH BISECTION CAPACITY
c     call xfer(np,cr_h)

      call exitt0

      end
c--------------------------------------------------------------
      subroutine set_f(f,c,n)
      real f(n),c(n)

      do i=1,n
         arg  = 1.e9*(i*i)
         arg  = 1.e9*cos(arg)
         f(i) = sin(arg)
      enddo

      call dssum(f)
      call col2 (f,c,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine init_dim

C     Transfer array dimensions to common

      include 'SIZE'
      include 'INPUT'

      ny1=nx1
      nz1=nx1

      ndim=ldim

      return
      end
c-----------------------------------------------------------------------
      subroutine init_mesh(ifbrick,npx,npy,npz,mx,my,mz)
      include 'SIZE'
      include 'TOTAL'
      logical ifbrick
      integer e,eg,offs

      if(.not.ifbrick) then   ! A 1-D array of elements of length P*lelt
  10     continue
         nelx = nelt*np
         nely = 1
         nelz = 1

         do e=1,nelt
            eg = e + nid*nelt
            lglel(e) = eg
         enddo
      else              ! A 3-D block of elements
        if (npx*npy*npz .ne. np) then
          call cubic(npx,npy,npz,np)  !xyz distribution of total proc
        end if
        if (mx*my*mz .ne. nelt) then
          call cubic(mx,my,mz,nelt)   !xyz distribution of elements per proc
        end if

c       if(mx.eq.nelt) goto 10

        nelx = mx*npx
        nely = my*npy
        nelz = mz*npz

        e = 1
        offs = (mod(nid,npx)*mx) + npx*(my*mx)*(mod(nid/npx,npy))
     $      + (npx*npy)*(mx*my*mz)*(nid/(npx*npy))
        do k = 0,mz-1
        do j = 0,my-1
        do i = 0,mx-1
           eg = offs+i+(j*nelx)+(k*nelx*nely)+1
           lglel(e) = eg
           e        = e+1
        enddo
        enddo
        enddo
      endif

      if (nid.eq.0) then
        write(6,*) "Processes: npx= ", npx, " npy= ", npy, " npz= ", npz
        write(6,*) "Local Elements: mx= ", mx, " my= ", my, " mz= ", mz
        write(6,*) "Elements: nelx= ", nelx, " nely= ", nely,
     &             " nelz= ", nelz
      end if

      return
      end
c-----------------------------------------------------------------------
      subroutine cubic(mx,my,mz,np)

      mx = np
      my = 1
      mz = 1
      ratio = np

      iroot3 = np**(1./3.) + 0.000001
      do i= iroot3,1,-1
        iz = i
        myx = np/iz
        nrem = np-myx*iz

        if (nrem.eq.0) then
          iroot2 = myx**(1./2.) + 0.000001
          do j=iroot2,1,-1
            iy = j
            ix = myx/iy
            nrem = myx-ix*iy
            if (nrem.eq.0) goto 20
          enddo
   20     continue

          if (ix < iy) then
            it = ix
            ix = iy
            iy = it
          end if

          if (ix < iz) then
            it = ix
            ix = iz
            iz = it
          end if

          if (iy < iz) then
            it = iy
            iy = iz
            iz = it
          end if

          if ( REAL(ix)/iz < ratio) then
            ratio = REAL(ix)/iz
            mx = ix
            my = iy
            mz = iz
          end if

        end if
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine set_multiplicity (c)       ! Inverse of counting matrix
      include 'SIZE'
      include 'TOTAL'

      real c(1)

      n = nx1*ny1*nz1*nelt

      call rone(c,n)
      call gs_op(gsh,c,1,1,0)  ! Gather-scatter operation  ! w   = QQ  w

      do i=1,n
         c(i) = 1./c(i)
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine set_timer_flop_cnt(iset)
      include 'SIZE'
      include 'TOTAL'
      include 'TIMER'

      integer i, numThrd, totThd
      integer omp_get_max_threads
      real tmp1(8), tmp2(8), tmp3(8)

      real time0,time1
      save time0,time1

      if (iset.eq.0) then
         flop_a  = 0
         flop_cg = 0

         do i = 1, tmax
           tzc(i)     = 0
           tglsc3(i)  = 0
           tax(i)     = 0
           tadd2sx(i) = 0
           tadd2s2(i) = 0
           tgsop(i)   = 0
           taxe(i)    = 0
         end do

         time0   = dnekclock()
      else
        time1   = dnekclock()-time0
        if (time1.gt.0) mflops = (flop_a+flop_cg)/(1.e6*time1)

        if (nid.eq.0) then
          write(6,1) nelt,np,nx1, nelt*np
          write(6,2) mflops*np, mflops
          write(6,3) flop_a,flop_cg,time1
        end if

    1   format('nelt = ', i7, ', np = ', i9, ', nx1 = ', i7,
     &         ', elements =', i10 )
    2   format('Tot MFlops = ', 1pe12.4, ', MFlops = ', e12.4)
    3   format('Ax FOp = ', 1pe12.4, ', CG FOp = ', e12.4,
     &         ', Solve Time = ', e12.4)

#ifdef TIMERS
        numThrd = 1
#ifdef _OPENMP
        numThrd = omp_get_max_threads()
#endif
        totThd = numThrd*np

        do i = 1,8
          tmp2(i) = 0.0
        end do

        tmp2(1)= time1
        do i = 1, numThrd
          tmp2(2)= tmp2(2) + tzc(i)
          tmp2(3)= tmp2(3) + tglsc3(i)
          tmp2(4)= tmp2(4) + tadd2sx(i)
          tmp2(5)= tmp2(5) + tax(i)
          tmp2(6)= tmp2(6) + tadd2s2(i)
          tmp2(7)= tmp2(7) + tgsop(i)
          tmp2(8)= tmp2(8) + taxe(i)
        end do

        call gop(tmp2, tmp3, '+  ', 8)

        if (nid.eq.0) then
          write(6,4) "av time: ", tmp2(1)/np, tmp2(2)/totThd,
     &               tmp2(3)/totThd, tmp2(4)/totThd, tmp2(5)/totThd
          write(6,5) "av time: ", tmp2(5)/totThd, tmp2(6)/totThd,
     &               tmp2(7)/totThd, tmp2(8)/totThd
        endif

        tmp2(1)= time1
        tmp2(2)= tzc(1)
        tmp2(3)= tglsc3(1)
        tmp2(4)= tadd2sx(1)
        tmp2(5)= tax(1)
        tmp2(6)= tadd2s2(1)
        tmp2(7)= tgsop(1)
        tmp2(8)= taxe(1)

        do i = 2, numThrd
          if (tzc(i) < tmp2(2)) tmp2(2)= tzc(i)
          if (tglsc3(i) < tmp2(3)) tmp2(3)= tglsc3(i)
          if (tadd2sx(i) < tmp2(4)) tmp2(4)= tadd2sx(i)
          if (tax(i) < tmp2(5)) tmp2(5)= tax(i)
          if (tadd2s2(i) < tmp2(6)) tmp2(6)= tadd2s2(i)
          if (tgsop(i) < tmp2(7)) tmp2(7)= tgsop(i)
          if (taxe(i) < tmp2(8)) tmp2(8)= taxe(i)
        end do

        call gop(tmp2, tmp3, 'm  ', 8)

        if (nid.eq.0) then
          write(6,4) "min time: ", tmp2(1), tmp2(2), tmp2(3),
     &               tmp2(4), tmp2(5)
          write(6,5) "min time: ", tmp2(5), tmp2(6), tmp2(7),
     &               tmp2(8)
        endif

        tmp2(1)= time1
        tmp2(2)= tzc(1)
        tmp2(3)= tglsc3(1)
        tmp2(4)= tadd2sx(1)
        tmp2(5)= tax(1)
        tmp2(6)= tadd2s2(1)
        tmp2(7)= tgsop(1)
        tmp2(8)= taxe(1)

        do i = 2, numThrd
          if (tzc(i) > tmp2(2)) tmp2(2)= tzc(i)
          if (tglsc3(i) > tmp2(3)) tmp2(3)= tglsc3(i)
          if (tadd2sx(i) > tmp2(4)) tmp2(4)= tadd2sx(i)
          if (tax(i) > tmp2(5)) tmp2(5)= tax(i)
          if (tadd2s2(i) > tmp2(6)) tmp2(6)= tadd2s2(i)
          if (tgsop(i) > tmp2(7)) tmp2(7)= tgsop(i)
          if (taxe(i) > tmp2(8)) tmp2(8)= taxe(i)
        end do

        call gop(tmp2, tmp3, 'M  ', 8)

        if (nid.eq.0) then
          write(6,4) "max time: ", tmp2(1), tmp2(2), tmp2(3),
     &               tmp2(4), tmp2(5)
          write(6,5) "max time: ", tmp2(5), tmp2(6), tmp2(7),
     &               tmp2(8)
        endif

    4   format(A, ' cg= ', 1pe12.4, ', zcm= ', e12.4,
     &         ', glsc3= ', e12.4, ', add2sx= ', e12.4,
     &         ', ax= ', e12.4)
    5   format(A, ' ax= ', 1pe12.4, ', add2s2= ', e12.4,
     &         ', gsop= ', e12.4, ', axe= ', e12.4)
#endif
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine xfer(np,gsh)
      include 'SIZE'
      parameter(npts_max = lx1*ly1*lz1*lelt)

      real buffer(2,npts_max)
      integer ikey(npts_max)

      nbuf = 800
      npts = 1
      do itest=1,200
         npoints = npts*np

         call load_points(buffer,nppp,npoints,npts,nbuf)
         iend   = mod1(npoints,nbuf)
         istart = 1
         if(nid.ne.0)istart = iend+(nid-1)*nbuf+1
         do i = 1,nppp
            icount=istart+(i-1)
            ikey(i)=mod(icount,np)
         enddo

         call nekgsync
         time0 = dnekclock()
         do loop=1,50
            call crystal_tuple_transfer(gsh,nppp,npts_max,
     $                ikey,1,ifake,0,buffer,2,1)
         enddo
         time1 = dnekclock()
         etime = (time1-time0)/50

         if (nid.eq.0) write(6,1) np,npts,npoints,etime
   1     format(2i7,i10,1p1e12.4,' bandwidth' )
         npts = 1.02*(npts+1)
         if (npts.gt.npts_max) goto 100
      enddo
 100  continue

      return
      end
c-----------------------------------------------------------------------
      subroutine load_points(buffer,nppp,npoints,npts,nbuf)
      include 'SIZE'
      include 'PARALLEL'

      real buffer(2,nbuf)

      nppp=0
      if(nbuf.gt.npts) then
       npass = 1+npoints/nbuf

       do ipass = 1,npass
          if(nid.eq.ipass.and.ipass.ne.npass) then
            do i = 1,nbuf
             buffer(1,i)=i
             buffer(2,i)=nid
            enddo
            nppp=nbuf
          elseif (npass.eq.ipass.and.nid.eq.0) then
            mbuf=mod1(npoints,nbuf)
            do i=1,mbuf
               buffer(1,i)=i
               buffer(2,i)=nid
            enddo
            nppp=mbuf
          endif
       enddo
      else
       do i = 1,npts
          buffer(1,i)=i
          buffer(2,i)=nid
       enddo
       nppp=npts
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine read_param(ifbrick,iel0,ielN,ielD,nx0,nxN,nxD,
     +                      npx,npy,npz,mx,my,mz)
      include 'SIZE'
      logical ifbrick
      integer iel0,ielN,ielD,nx0,nxN,nxD,npx,npy,npz,mx,my,mz

      !open .rea
      if(nid.eq.0) then
         open(unit=9,file='data.rea',status='old')
         read(9,*,err=100) ifbrick
         read(9,*,err=100) iel0,ielN,ielD
         read(9,*,err=100) nx0,nxN,nxD
         read(9,*,err=100) npx,npy,npz
         read(9,*,err=100) mx,my,mz
         close(9)
      endif
      call bcast(ifbrick,4)
      call bcast(iel0,4)
      call bcast(ielN,4)
      call bcast(ielD,4)
c     nx0=lx1
c     nxN=lx1
      call bcast(nx0,4)
      call bcast(nxN,4)
      call bcast(nxD,4)
      call bcast(npx,4)
      call bcast(npy,4)
      call bcast(npz,4)
      call bcast(mx,4)
      call bcast(my,4)
      call bcast(mz,4)
      if(iel0.gt.ielN.or.nx0.gt.nxN) goto 200

      return

  100 continue
      write(6,*) "ERROR READING data.rea....ABORT"
      call exitt0

  200 continue
      write(6,*) "ERROR data.rea :: iel0 > ielN or nx0 > nxN :: ABORT"
      call exitt0

      return
      end
c-----------------------------------------------------------------------

max_stars_repo_name: HPCToolkit/hpctest

max_issues_repo_name: HPCToolkit/hpctest

max_forks_repo_name: HPCToolkit/hpctest


text:
	SUBROUTINE PC_INTH  ( hght, adata, bdata, nparms, intflg,
     +			      angflg, jhght, outdat, iret )
C************************************************************************
C* PC_INTH								*
C*									*
C* This subroutine interpolates between two levels of data.  The	*
C* data are interpolated with respect to the height.  Height is 	*
C* found in the input arrays in the location given by jhght.		*
C* If errors are encountered no output data are changed.  Therefore,	*
C* data should be set to RMISSD before calling this subroutine.		*
C*									*
C* PC_INTH  ( HGHT, ADATA, BDATA, NPARMS, INTFLG, ANGFLG, JHGHT,	*
C*            OUTDAT, IRET )						*
C*									*
C* Input parameters:							*
C*	HGHT		REAL		Pressure to be interpolated	*
C*	ADATA (NPARMS)	REAL		Data at first level		*
C*	BDATA (NPARMS)	REAL		Data at second level		*
C*	NPARMS		INTEGER		Size of data arrays		*
C*	INTFLG (NPARMS)	LOGICAL		Interpolation flags		*
C*	ANGFLG (NPARMS)	LOGICAL		Angle interpolation flags	*
C*	JHGHT		INTEGER		Location of HGHT		*
C*									*
C* Output parameters:							*
C*	OUTDAT (NPARMS)	REAL		Data interpolated to HGHT	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-24 = vlev not between levels	*
C*					-28 = invalid pressure 		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 3/91	Add interpolation using hght		*
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
	REAL		adata (*), bdata (*), outdat (*)
	LOGICAL		intflg (*), angflg (*)
C*
	LOGICAL		betwen
C*
        INCLUDE         'ERMISS.FNC'
C-------------------------------------------------------------------------
C*	Check that HGHT is between two pressures input.
C
	hght1 = adata (jhght)
	hght2 = bdata (jhght)
	betwen = ( ( ( hght1 .lt. hght2 ) .and. ( hght1 .lt. hght ) .and.
     +		     ( hght .lt. hght2 ) )  .or.
     +		   ( ( hght2 .lt. hght1 ) .and. ( hght2 .lt. hght ) .and.
     +		     ( hght .lt. hght1 ) ) )
C
C*	Test for errors - can't interpolate or a non-positive pressure.
C
	IF  ( .not. betwen )  THEN
	    iret = -24
	  ELSE
	    iret = 0
C
C*	    Set up interpolation.  Check for interpolation of angles
C*	    through 360 degrees.
C
	    rmult = ( hght - adata (jhght) ) /
     +			( bdata (jhght) - adata (jhght) )
	    outdat (jhght) = hght
	    DO  i = 1, nparms
	      IF  ( i .ne. jhght )  THEN
		IF  ( intflg (i) )  THEN
		    IF  ( angflg (i) )  THEN
			ang1 = AMOD ( adata(i), 360. )
			ang2 = AMOD ( bdata(i), 360. )
			IF  ( ABS (ang1 - ang2) .gt. 180.)  THEN
			    IF  ( ang1 .lt. ang2 )  THEN
				ang1 = ang1 + 360.
			      ELSE
				ang2 = ang2 + 360.
			    END IF
			END IF
			ang = ang1 + (ang2 - ang1) * rmult
			outdat (i) = AMOD ( ang, 360. )
		      ELSE
		        outdat (i) = adata (i) +
     +				      ( bdata(i) - adata (i) ) * rmult
		    END IF
		    IF  ( ERMISS ( adata(i) ) .or.
     +				ERMISS ( bdata(i) ) )
     +					outdat (i) = RMISSD
		END IF
	      END IF
	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!PR fortran/32238
! { dg-do compile }
! { dg-final { cleanup-modules "bug_test" } }

module bug_test

contains
  subroutine bug(c)

  implicit none

  integer, parameter :: fp = selected_real_kind(13)
  complex(kind=fp)              :: c(:,:)
  where( abs( aimag( c ) ) < 1.e-10_fp )                             &
  &    c = cmplx( real( c , fp ) , 0._fp , fp )
  where( abs( real( c , fp ) ) < 1.e-10_fp )                         &
  &    c = cmplx( 0._fp , aimag( c ) , fp )

  return
  end subroutine bug

end module bug_test

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      MODULE vmec_params
      USE stel_kinds, ONLY: rprec, dp
      USE vparams, ONLY: mpold
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER, PARAMETER :: meven = 0, modd = 1
      INTEGER, PARAMETER :: ndamp = 10
      INTEGER, PARAMETER :: ns4 = 25

      INTEGER, PRIVATE :: ink
      INTEGER, PARAMETER, DIMENSION(0:mpold) ::
     1  jmin1 = (/ 1,1,(2,ink=2,mpold) /),        !starting js(m) values where R,Z are non-zero
     2  jmin2 = (/ 1,2,(2,ink=2,mpold) /),        !starting js(m) values for which R,Z are evolved
     3  jlam  = (/ 2,2,(2,ink=2,mpold) /)         !starting js(m) values for which Lambda is evolved

!  Besure to update werror in fileout.f when adding more error flags.
      INTEGER, PARAMETER :: norm_term_flag=0, bad_jacobian_flag=1,
     1                      more_iter_flag=2,
     2                      jac75_flag=4, input_error_flag=5,
     3                      phiedge_error_flag=7,
     4                      ns_error_flag=8,
     5                      misc_error_flag=9,
     6                      successful_term_flag=11, !ftol force criterion has been met
     7                      bsub_bad_js1_flag=12,
     8                      r01_bad_value_flag=13,
     9                      arz_bad_value_flag=14
      INTEGER, PARAMETER :: restart_flag=1, readin_flag=2,
     1                      timestep_flag=4,output_flag=8,
     2                      cleanup_flag=16, reset_jacdt_flag=32

      REAL(rprec), PARAMETER :: pdamp = 0.05_dp
      CHARACTER(LEN=*), PARAMETER :: version_ = '9.0'
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: ntmax, rcc, rss, rsc, rcs, zsc, zcs, zcc, zss
      INTEGER :: mnyq, nnyq
      INTEGER, ALLOCATABLE :: uminus(:)
      REAL(rprec), ALLOCATABLE :: mscale(:), nscale(:)
      REAL(rprec) :: signgs, lamscale
!-----------------------------------------------
!
!     VERSION INFORMATION HISTORY
!
!     8.00
!          a) added lforbal logical to fbal module to control whether to compute the flux-averaged
!             force balance equation printed in the threed1 file. This requires a modification of
!             the m=1,n=0 forces for R,Z in tomnsps subroutine. This works well, generally, and
!             yields an improved <EQUIF> in threed1 file. However, there have been some cases where
!             this non-variational departure fails to converge.
!          b) added "bias" to iotas in bcovar when determining preconditioner for very small iota
!             values. Seems to need this for improving convergence of preconditioned current-hole cases.
!             Eliminated in v8.20.
!          c) added "totzsps,a_hess" to recompute r,z,l rapidly for only those modes that are jogged during
!             Hessian calculation. NOTE: need to do this for lasym=true case, as well, eventually
!     8.20  (August, 2004)
!          a) removed 2-pt tie at origin for m=1 modes, to preserve tri-diagonal structure of Hessian.
!             This is needed for preconditioning, which assumes block-tridi structure of equations
!          b) fixed problem with free-boundary preconditioner, namely, ctor can not be extrapolated
!             at edge when computing preconditioner, because this breaks tri-diagonal structure
!          c) added new variables to input file to control preconditioning:
!             1) PRECON_TYPE: = 'default', default tri-di (block size = 1)
!                             = 'cg',      block tri-di, conjugate-gradient time-stepper
!                             = 'gmres',   "          ", gmres time-stepper
!                             = 'tfqmr',   "          ", transpose free qmr
!             2) PREC2D_THRESHOLD: value of (unpreconditioned) forces at which block (2D) preconditioner
!                                  is turned on (=0 block preconditioner never turned on); recommended
!                                  (default) value ~ 1.E-10, or smaller, if convergence is poor
!             3) LFORBAL: logical variable (default = .true.); when true, the force balance
!                         used in the threed1 file is used to evolve the m=1 R,Z components. This
!                         is a non-variational departure from the force equations for these modes,
!                         but generally does not have an unfavorable impact on convergence.
!          d) added new internal variable, ICTRL_PREC2D, to precon2d module. Replaces previous lprec2d
!             and iequi>1 variables.
!          e) removed lsweep_fast option from module precon2d. This slows the computation of the Hessian
!             by about 2/3, but is more accurate (includes pdamp, liota, lforbal correctly)
!          f) removed lflam logicals from bcovar and tomnsps, since now we must compute dFlam/dR,Z by
!             jogging
!          g) removed Compute_Hess_Flam_RZ from lamblks; this is now computed via jogging
!             (also removed Get_dFlam_dRZ, FFT2Hessian, Forbal_avg, GetGsqrtVar supporting routines)
!          h) removed internal liota logic, used to push iota profile rather than solving for it. Had
!             needed this for symmetric Hessian (lsweep_fast=true option), but no longer required. Also,
!             it was not implemented entirely correctly for lforbal=true case
!          i) for lasym m=1 constraint rsc = zcc, changed xc array so that R+ = .5*(rsc + zcc) is stored at
!             xc(rsc,m=1) and R- = .5*(rsc - zcc) is stored at xc(zcc,m=1). In residue, gcz(R-) == gcz(zcc)
!             is zeroed by "evolving" gcr(zcc) = azd*[xc(zcc)-xcint], and gcr(rsc) => .5*[gcr(rsc) + gcz(zcc)]
!             is evolved. In totzspa, the real rsc and zcc are computed from the internal representations
!             (check convert call, too) by calling a new routine convert_asym (also called from wrout before
!             writing out xc info). In readin, the original R+,R- are stored, so that for external "jogs",
!             there will be no change in forces. All these changes are needed to obtain an invertible Hessian.
!          j) added m=1 constraint for 3D case (similar to (i)), rss(n) = zcs(n), for n != 0. Imposed this
!             on forces by adding routine constrain_m1 in residue. Added convert_sym routine to totzsp to convert
!             from internal xc representation TO internal one.
!          k) Decreased exponent on pdamp factor r2 (in bcovar) from 2 to 1, to give better conditioning
!             especially for current hole cases
!          l) Eliminated iotas bias for determining preconditioner, previously added in v8.00 for stabilizing
!             current hole cases (not needed with corrected preconditioner)
!     8.30  (October, 2004)
!          a) Implemented flags for "reverse-communication" mode of vmec
!     8.40 a) Converted the m=1 constraints for 3D and asym back to old way; did not always
!             converge well with the new constraints introduced in 8.20 (i-j)
!     8.45  (December, 2005)
!          a) Added the lconm1 logical. If = True, new constraint; if = False, old m=1 constraint used
!          b) Added "perturbation" computation for lasym=TRUE case (totzspa_hess)
!     8.46  (June, 2009)
!          a) Added LRFP logical to allow easy switching of profiles between Stellarator/tokamak (PHIP=1, LRFP=F)
!             and RFP (CHIP=1, LRFP=T). When LRFP=T, AI coefficients are expansion of q = 1/iota. Added lrfp to
!             LIBSTELL/vmec_input module.
!     8.47  (July, 2010)
!          a) Rewrote magnetic field representation so that phip*lambda = new internal lambda. This greatly improves
!             the conditioning of the lambda equations which otherwise become singular at the RFP reversal point
!     8.48  (March 2012 - JDH)
!          a) Accumulated small changes from SPH & JDH
!          b) Modifications from J Geiger, March 2012
!             - to be able to get additional main iterations if the force tolerance is
!                  not met. Parameter MAX_MAIN_ITERATIONS specifies how many main iteration
!                  cycles should be run.
!             - to get a full output in the threed1-file if the force tolerance is not
!                  met. Specify the logical LFULL3D1OUT to true for this.
!             - if vmec2000 is compiled with netcdf, you can still get the ascii-output
!                  if you specify the logical LWOUTTXT as true.
!             - you get the output for diagno 1.0 and 1.5 if the logical LDIAGNO set true.
!             - you get a rather old fort.8-output if you specify LOLDOUT as true.
!
!             If none of these new variables is set, the behavior of vmec2000 is as
!             expected from the version without the changes.
!     8.49  (June, 2012)
!          a) Fixed bug in bcovar when averaging half-grid covariant components onto full grid: did not
!             zero components at (phantom) js=1 point, so edge force averages were incorrect
!          b) Added lamscale factor to scale lambda in contravariant B-components. Modify wrout to maintain
!             old-style lambda output
!          c) Set lbsubs=F in jxbforce by default to capture current sheets
!          d) Added lmove_axis INPUT logical (=T by default) so user can control whether or not the magnetic
!             axis can be initially shifted to improve the initial force residuals. It is important NOT to move
!             the helical axis for RFP equilibria requiring a helical seed (set l_moveaxis=F for this case!)
!
!     8.50   (Jan, 2013)
!          a) Improved scaling of lambda forces with respect to lamscale
!          b) Fixed fnorm1 scaling (remove hs dependence)
!          c) Added lgiveup logical (M. Drevlak/J. Geiger)
!-----------------------------------------------
      END MODULE vmec_params


text:
! { dg-do compile }
! Tests the patch for PRs 25084, 20852, 25085 and 25086, all of
! which involve assumed character length functions.
! This test checks the things that should not emit errors.
!
! Contributed by Paul Thomas  <pault@gcc.gnu.org>
!
function is_OK (ch)                ! { dg-warning "is obsolescent in fortran 95" }
  character(*) is_OK, ch           ! OK in an external function
  is_OK = ch
end function is_OK

! The warning occurs twice for the next line; for 'more_OK' and for 'fcn';
function more_OK (ch, fcn)         ! { dg-warning "is obsolescent in fortran 95" }
  character(*) more_OK, ch
  character (*), external :: fcn   ! OK as a dummy argument
  more_OK = fcn (ch)
end function more_OK

  character(4) :: answer
  character(4), external :: is_OK, more_OK

  answer = is_OK ("isOK")          ! LEN defined in calling scope
  print *, answer

  answer = more_OK ("okay", is_OK) ! Actual arg has defined LEN
  print *, answer

  answer = also_OK ("OKOK")
  print *, answer

contains
  function also_OK (ch)
    character(4) also_OK
    character(*) ch
    also_OK = is_OK (ch)            ! LEN obtained by host association
  end function also_OK
END

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
	SUBROUTINE HQCOMP ( icolr, color, ired, igreen, iblue, xname,
     +			    iret )
C************************************************************************
C* HQCOMP - XWP								*
C*									*
C* This subroutine returns the red, green, and blue components of a	*
C* color.  The color components are defined in the range 0 - 255.  If	*
C* the color was defined by name, the GEMPAK color name is also		*
C* returned.  The X Window System color name is returned if it is	*
C* available.								*
C*									*
C* HQCOMP  ( ICOLR, COLOR, IRED, IGREEN, IBLUE, XNAME, IRET )		*
C*									*
C* Input parameters:							*
C* 	ICOLR		INTEGER		Color number			*
C*									*
C* Output parameters:							*
C*	COLOR		CHAR*		GEMPAK color name		*
C*	IRED		INTEGER		Red color component		*
C*	IGREEN		INTEGER		Green color component		*
C*	IBLUE		INTEGER		Blue color component		*
C*	XNAME		CHAR*		X Window System color name	*
C*	IRET		INTEGER 	Return code			*
C**									*
C* Log:									*
C* L. Williams/EAI	 3/96						*
C* S. Jacobs/NCEP	10/96	Added checks for sub-devices		*
C* S. Maxwell/GSC        6/97   Documentation changes                   *
C************************************************************************
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'ERROR.PRM'
C*
	CHARACTER*(*) 	color, xname
C-----------------------------------------------------------------------
	iret = NORMAL
C
C*	Retrieve the values from structure.
C
	color(1:1) = ' '
	xname(1:1) = ' '
C
C*	Check the requested device.
C
	IF  ( ( ddev .eq. 'XW' ) .or. ( ddev .eq. 'XWP' ) )  THEN
	    CALL XQCOMP ( icolr, color, ired, igreen, iblue, xname,
     +			  ilen, ilxn, iret )
	  ELSE IF  ( ddev .eq. 'PS' )  THEN
	    CALL PQCOMP ( icolr, color, ired, igreen, iblue, xname,
     +			  ilen, ilxn, iret )
	END IF
C
C*	Clear the unused space of Color name and X color name
C
 	color(ilen+1: ) = ' '
 	xname (ilxn+1: ) = ' '
C*
	RETURN
	END


text:
subroutine readnc(lunmd     ,error     ,keyw      ,newkw     ,nlook     , &
                & record    ,cvar      ,cdefau    ,lenc      ,nrrec     , &
                & ntrec     ,lundia    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: readnc.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/general/readnc.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Reads NLOOK character values with a fixed length
!              from ONE record of MD-file.
!              nlook may NOT be 0
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer , pointer :: ifis
    integer , pointer :: itis
!
! Global variables
!
    integer                                 :: lenc   !!  Help var. (length of var. cvar to be looked for in the MD-file)
    integer                                 :: lundia !  Description and declaration in inout.igs
    integer                                 :: lunmd  !  Description and declaration in inout.igs
    integer                                 :: nlook  !!  Help var.: nr. of data to look for in the MD-file
    integer                                 :: nrrec  !!  Pointer to the record number in the MD-file
    integer                                 :: ntrec  !!  Help. var to keep track of NRREC
    logical                                 :: error  !!  Flag=TRUE if an error is encountered
    logical                                 :: newkw  !!  Logical var. specifying whether a
                                                      !!  new recnam should be read from the
                                                      !!  MD-file or just new data in the
                                                      !!  continuation line
    character(*)               , intent(in) :: cdefau
    character(*)                            :: record !!  Record read from either the MD-file or from the attribute file
    character(*), dimension(*)              :: cvar   !!  Character variable to read from a record
    character(6)                            :: keyw   !!  Name of record to look for in the
                                                      !!  MD-file (usually KEYWRD or RECNAM)
!
!
! Local variables
!
    integer     :: ibeg       ! Begin position in the RECORD from where the search for data/record is started
    integer     :: iend       ! Last position in the RECORD when the searched data/record is finished
    integer     :: iendm
    integer     :: iocond
    integer     :: lenc0
    integer     :: lkw        ! Length of char. str (usually the KEYWRD or RECNAM)
    integer     :: lrec       ! Help var. containing the length of RECORD
    integer     :: n
    integer     :: nn
    integer     :: nnlook     ! Help var.: nr. of data really read from the MD-file
    logical     :: found      ! If FOUND = TRUE then recnam in the MD-file was found
!
!
!! executable statements -------------------------------------------------------
!
    !
    !
    ifis  => gdp%gdrdpara%ifis
    itis  => gdp%gdrdpara%itis
    !
    lrec   = len(record)
    lkw    = index(keyw, ' ') - 1
    if (lkw<=0) lkw = 6
    found  = .false.
    error  = .false.
    lenc0  = lenc
    nnlook = nlook
    !
    !-----look for keyword in record
    !
    ibeg = itis
    call search(lunmd     ,error     ,newkw     ,nrrec     ,found     , &
              & ntrec     ,record    ,ibeg      ,keyw      ,lkw       , &
              & 'NO'      )
    !
    !-----if not found then check nlook
    !     nlook > 0, then wrong text, else no text but ok
    !
    if (.not.found) then
       call prterr(lundia    ,'U100'    ,keyw      )
       !
       error = .true.
       goto 999
    endif
    !
    !-----read value if ier <> 0 then no value not allowed =>
    !     because txtmrk creates string between ##, the initialisation of
    !     iendm should defined 1 position back
    !
    n = 0
    iendm = ifis - 1 - 1
    !-->
   10 continue
    n      = n + 1
    nnlook = n
    ibeg   = iendm + 2
    !
    !--------calculate part of mdfrec with text to search
    !
    call txtmrk(record    ,lrec      ,ibeg      ,iend      ,error     )
    iendm = iend
    !
    !--------read cvar, if blank set to default value
    !
    if (.not.error) then
       if (record(ibeg:iend)/=' ') then
          call txtstr(record    ,ibeg      ,iend      ,lenc      ,cvar(n)   , &
                    & error     )
          lenc = lenc0
       else
          cvar(n) = cdefau
       endif
    elseif (nlook==1) then
       error = .true.
    else
       error = .false.
       read (lunmd, '(a300)', iostat = iocond) record
       !
       !--------------iocond = 0, read continuation record
       !                          ensure the record is all blanks before = sign
       !                          and reset value of n (=> n < nlook per def.)
       !                     < 0, end of file encountered (not ok)
       !                     > 0, error occurred (not ok)
       !
       if (iocond==0) then
          nrrec = nrrec + 1
          if (record(:itis)/=' ') then
             do nn = n, nlook
                cvar(nn) = cdefau
             enddo
             nnlook = n - 1
             n = nlook
          else
             iendm = ifis - 1 - 1
             n = n - 1
          endif
       else
          error = .true.
       endif
    endif
    !
    if (error) then
       call prterr(lundia    ,'U036'    ,keyw      )
       !
       goto 999
    endif
    if (n<nlook) goto 10
    ! <--
    !
  999 continue
    nlook = nnlook
end subroutine readnc

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!
!     -------------------------------------------------------------
!      I T L S
!     -------------------------------------------------------------
!
!     Written by G. Gaigalas,                                      *
!     Vanderbilt University,  Nashville             October 1996   *
!
      INTEGER FUNCTION ITLS (IK, ID, KG, KGG, BD, IBT, BT, KG1, KGG1, ITP, ITG&
         , IQ)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE CONSTS_C
      USE RIBOLS_C
      USE RIBOLSF_C
      USE RIBOF_C
      USE RIBOLS3_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  10:51:45  11/16/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE ittk_I
      USE mes_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: KG
      INTEGER , INTENT(IN) :: KGG
      INTEGER  :: KG1
      INTEGER  :: KGG1
      INTEGER , INTENT(OUT) :: ITP
      INTEGER , INTENT(OUT) :: ITG
      INTEGER , INTENT(IN) :: IQ
      INTEGER  :: IK(7)
      INTEGER  :: ID(7)
      INTEGER , INTENT(OUT) :: IBT(7)
      REAL(DOUBLE) , INTENT(IN) :: BD(3)
      REAL(DOUBLE) , INTENT(OUT) :: BT(3)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: ITK, ITD, ITP1, ITG1
!-----------------------------------------------
      ITLS = 0
      IF (ID(3) > 9) RETURN
      KG1 = 2*KG
      KGG1 = 2*KGG
      IF (ITTK(ID(5),IK(5),KG1) == 0) RETURN
      IF (ITTK(ID(6),IK(6),KGG1) == 0) RETURN
      ITK = IK(1)
      ITD = ID(1)
      IF (ID(3) < 3) THEN
         ITP1 = IMPTLS(ITK)
         ITP = IMPTLS(ITD)
         IF (ITP1 /= ITP) RETURN
         ITG1 = IMGTLS(ITK)
         ITG = IMGTLS(ITD)
      ELSE IF (ID(3) == 3) THEN
         IF (ITK > 300) THEN
            IF (ITD < 300) CALL MES (51)
            IF (ID(4) > 2) CALL MES (11)
            IF (IK(4) > 2) CALL MES (11)
            ITK = ITK - 300
            ITD = ITD - 300
            ITP1 = IMPTLSF(ITK)
            ITP = IMPTLSF(ITD)
            IF (ITP1 /= ITP) RETURN
            ITG1 = IMGTLSF(ITK)
            ITG = IMGTLSF(ITD)
         ELSE
            IF (ITD > 300) CALL MES (51)
            ITP1 = IMPTF(ITK)
            ITP = IMPTF(ITD)
            IF (ITP1 /= ITP) RETURN
            ITG1 = IMGTF(ITK)
            ITG = IMGTF(ITD)
         ENDIF
      ELSE
         IF (ID(4) > 2) CALL MES (11)
         IF (IK(4) > 2) CALL MES (11)
         ITP1 = IMPTLS3(ITK)
         ITP = IMPTLS3(ITD)
         IF (ITP1 /= ITP) RETURN
         ITG1 = IMGTLS3(ITK)
         ITG = IMGTLS3(ITD)
      ENDIF
      IF (ITG1 /= ITG) RETURN
      ITLS = 1
      IBT(2) = ID(2)
      IBT(3) = ID(3)
      IBT(4) = ID(4) + IQ
      BT(3) = BD(3) + HALF*DBLE(IQ)
      RETURN
      END FUNCTION ITLS


text:
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!  bt_data module
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------

      module bt_data

!---------------------------------------------------------------------
! The following include file is generated automatically by the
! "setparams" utility. It defines
!      maxcells:      the square root of the maximum number of processors
!      problem_size:  12, 64, 102, 162 (for class T, A, B, C)
!      dt_default:    default time step for this problem size if no
!                     config file
!      niter_default: default number of iterations for this problem size
!---------------------------------------------------------------------

      include 'npbparams.h'

      integer           aa, bb, cc, BLOCK_SIZE
      parameter        (aa=1, bb=2, cc=3, BLOCK_SIZE=5)

      integer           grid_points(3)
      double precision  elapsed_time
!!$acc declare create(grid_points)
!!$omp declare target to(grid_points)

      double precision  tx1, tx2, tx3, ty1, ty2, ty3, tz1, tz2, tz3,  &
     &                  dx1, dx2, dx3, dx4, dx5, dy1, dy2, dy3, dy4,  &
     &                  dy5, dz1, dz2, dz3, dz4, dz5, dssp, dt,  &
     &                  ce(5,13), dxmax, dymax, dzmax, xxcon1, xxcon2,  &
     &                  xxcon3, xxcon4, xxcon5, dx1tx1, dx2tx1, dx3tx1,  &
     &                  dx4tx1, dx5tx1, yycon1, yycon2, yycon3, yycon4,  &
     &                  yycon5, dy1ty1, dy2ty1, dy3ty1, dy4ty1, dy5ty1,  &
     &                  zzcon1, zzcon2, zzcon3, zzcon4, zzcon5, dz1tz1,  &
     &                  dz2tz1, dz3tz1, dz4tz1, dz5tz1, dnxm1, dnym1,  &
     &                  dnzm1, c1c2, c1c5, c3c4, c1345, conz1, c1, c2,  &
     &                  c3, c4, c5, c4dssp, c5dssp, dtdssp, dttx1,  &
     &                  dttx2, dtty1, dtty2, dttz1, dttz2, c2dttx1,  &
     &                  c2dtty1, c2dttz1, comz1, comz4, comz5, comz6,  &
     &                  c3c4tx3, c3c4ty3, c3c4tz3, c2iv, con43, con16

      integer IMAX, JMAX, KMAX
      parameter (IMAX=problem_size,JMAX=problem_size,KMAX=problem_size)

!
!   field arrays
!
      double precision, allocatable ::  &
     &   us      (:, :, :),  &
     &   vs      (:, :, :),  &
     &   ws      (:, :, :),  &
     &   qs      (:, :, :),  &
     &   rho_i   (:, :, :),  &
     &   square  (:, :, :),  &
     &   forcing (:, :, :, :),  &
     &   u       (:, :, :, :),  &
     &   rhs     (:, :, :, :)
!!$acc declare create(u, us, vs, ws, qs, rho_i,  &
!!$acc&               square, rhs, forcing)
!!$omp declare target to(u, us, vs, ws, qs, rho_i,  &
!!$omp&               square, rhs, forcing)
      double precision cuf(0:problem_size),   q  (0:problem_size),  &
     &                 ue (0:problem_size,5), buf(0:problem_size,5)

!
!-----------------------------------------------------------------------
!   Timer constants
!-----------------------------------------------------------------------
      integer t_rhsx, t_rhsy, t_rhsz, t_xsolve, t_ysolve, t_zsolve,  &
     &        t_rdis1, t_solsub, t_add, t_rhs, t_last, t_total
      parameter (t_total = 1)
      parameter (t_rhsx = 2)
      parameter (t_rhsy = 3)
      parameter (t_rhsz = 4)
      parameter (t_rhs = 5)
      parameter (t_xsolve = 6)
      parameter (t_ysolve = 7)
      parameter (t_zsolve = 8)
      parameter (t_rdis1 = 9)
      parameter (t_solsub = 10)
      parameter (t_add = 11)
      parameter (t_last = 11)

      logical timeron

      end module bt_data

!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!  work_lhs module
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------

      module work_lhs

      use bt_data, only : problem_size

      integer :: bsize
!!$acc declare create(bsize)
!!$omp declare target to(bsize)

      integer, parameter :: ndim = problem_size - 1
      integer, parameter :: ndim1 = ndim - 1

      double precision, allocatable ::  &
     &   fjac(:, :, :, :, :),  &
!     >   njac(:, :, :, :, :),
     &   lhsa(:, :, :, :, :),  &
     &   lhsb(:, :, :, :, :),  &
     &   lhsc(:, :, :, :, :),  &
     &   rhsx(:, :,    :, :)
!!$acc declare create(fjac, lhsa, lhsb, lhsc, rhsx)
!!$omp declare target to(fjac, lhsa, lhsb, lhsc, rhsx)
      external matvec_sub, matmul_sub, binvcrhs, binvrhs, lhsinit
!!$acc routine(matvec_sub) vector
!!$acc routine(matmul_sub) vector
!!$acc routine(binvcrhs) vector
!!$acc routine(binvrhs) vector
!!$acc routine(lhsinit) vector
!!$omp declare target (matvec_sub)
!!$omp declare target (matmul_sub)
!!$omp declare target (binvcrhs)
!!$omp declare target (binvrhs)
!!$omp declare target (lhsinit)

      end module work_lhs

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

      subroutine lhsinit(ni, k)

!---------------------------------------------------------------------
!---------------------------------------------------------------------
      use work_lhs, only : lhsa, lhsb, lhsc, bsize
      implicit none
!!$acc routine vector
!!$omp declare target

      integer ni, k

      integer i, m, n, j

!---------------------------------------------------------------------
!     zero the whole left hand side for starters
!     set all diagonal values to 1. This is overkill, but convenient
!---------------------------------------------------------------------
!!$acc loop seq
      do i = 0, ni, ni
!!$acc loop vector
!!!$omp target teams distribute parallel do collapse(3)
         do j = 1, bsize
!!$acc loop seq
            do m = 1, 5
!!$acc loop seq
               do n = 1, 5
                  lhsa(j,m,n,i,k) = 0.0d0
                  lhsb(j,m,n,i,k) = 0.0d0
                  lhsc(j,m,n,i,k) = 0.0d0
               end do
               lhsb(j,m,m,i,k) = 1.0d0
            end do
         end do
      end do

      return
      end subroutine

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

      subroutine alloc_space

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

!---------------------------------------------------------------------
! allocate space dynamically for data arrays
!---------------------------------------------------------------------

      use bt_data
      implicit none

      integer ios

      integer IMAXP, JMAXP
      parameter (IMAXP=IMAX/2*2,JMAXP=JMAX/2*2)

!
!   to improve cache performance, grid dimensions padded by 1
!   for even number sizes only.
!
      allocate (  &
     &   us      (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   vs      (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   ws      (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   qs      (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   rho_i   (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   square  (0:IMAXP, 0:JMAXP, 0:KMAX-1),  &
     &   forcing (0:IMAXP, 0:JMAXP, 0:KMAX-1,5),  &
     &   u       (0:IMAXP, 0:JMAXP, 0:KMAX-1,5),  &
     &   rhs     (0:IMAXP, 0:JMAXP, 0:KMAX-1,5),  &
     &         stat = ios)

      if (ios .ne. 0) then
         write(*,*) 'Error encountered in allocating space'
         stop
      endif

      return
      end

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

      subroutine alloc_lhs_space()

!---------------------------------------------------------------------
!---------------------------------------------------------------------
      use work_lhs
      implicit none

      integer ios

      allocate (  &
     &   fjac(ndim1, 5, 5, 0:ndim, ndim-1),  &
!     >   njac(ndim1, 5, 5, 0:ndim, ndim-1),
     &   lhsa(ndim1, 5, 5, 0:ndim, ndim-1),  &
     &   lhsb(ndim1, 5, 5, 0:ndim, ndim-1),  &
     &   lhsc(ndim1, 5, 5, 0:ndim, ndim-1),  &
     &   rhsx(ndim1, 5,    0:ndim, ndim-1),  &
     &         stat = ios)

      if (ios .ne. 0) then
         write(*,*) 'Error encountered in allocating space'
         stop
      endif

      return
      end

      program repro
      end program repro

max_stars_repo_name: aperdeus-amd/aomp

max_issues_repo_name: aperdeus-amd/aomp

max_forks_repo_name: aperdeus-amd/aomp


module tests
    use assert
    implicit none
contains
    subroutine test_settings_create
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1

        print *, '[  RUN     ] settings.create'
        s1 = YMMSL_Settings_create()
        call assert_true(YMMSL_Settings_empty(s1))
        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.create'
    end subroutine test_settings_create

    subroutine test_settings_equals
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1, s2

        print *, '[  RUN     ] settings.equals'
        s1 = YMMSL_Settings_create()
        s2 = YMMSL_Settings_create()
        call assert_true(YMMSL_Settings_equals(s1, s2))

        call YMMSL_Settings_set(s1, 'key', 'value')
        call assert_false(YMMSL_Settings_equals(s1, s2))

        call YMMSL_Settings_set(s2, 'key', 'value')
        call assert_true(YMMSL_Settings_equals(s1, s2))

        call YMMSL_Settings_free(s1)
        call YMMSL_Settings_free(s2)
        print *, '[       OK ] settings.equals'
    end subroutine test_settings_equals

    subroutine test_settings_size
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1

        print *, '[  RUN     ] settings.size'
        s1 = YMMSL_Settings_create()

        call assert_eq_size(YMMSL_Settings_size(s1), 0_YMMSL_size)

        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.size'
    end subroutine test_settings_size

    subroutine test_settings_set_get_as
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1
        real(YMMSL_real8), dimension(2) :: ra1, ra2
        real(YMMSL_real8), dimension(3, 2) :: ra3, ra4

        print *, '[  RUN     ] settings.set_get_as'
        s1 = YMMSL_Settings_create()
        ra1 = (/1.0d0, 2.0d0/)
        ra3 = reshape((/1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0/), (/3, 2/))

        call YMMSL_Settings_set(s1, 'key1', 'value1')
        call YMMSL_Settings_set(s1, 'key2', 242424242_YMMSL_int4)
        call YMMSL_Settings_set(s1, 'key3', 42424242424242_YMMSL_int8)
        call YMMSL_Settings_set(s1, 'key4', .false.)
        call YMMSL_Settings_set(s1, 'key5', 13.13d0)
        call YMMSL_Settings_set(s1, 'key6', ra1)
        call YMMSL_Settings_set(s1, 'key7', ra3)

        call assert_true(YMMSL_Settings_contains(s1, 'key1'))
        call assert_true(YMMSL_Settings_contains(s1, 'key2'))
        call assert_true(YMMSL_Settings_contains(s1, 'key3'))
        call assert_true(YMMSL_Settings_contains(s1, 'key4'))
        call assert_true(YMMSL_Settings_contains(s1, 'key5'))
        call assert_true(YMMSL_Settings_contains(s1, 'key6'))
        call assert_true(YMMSL_Settings_contains(s1, 'key7'))
        call assert_false(YMMSL_Settings_contains(s1, 'nokey'))

        call assert_eq_character(YMMSL_Settings_get_as_character(s1, 'key1'), 'value1')
        call assert_eq_int4(YMMSL_Settings_get_as_int4(s1, 'key2'), 242424242_YMMSL_int4)
        call assert_eq_int8(YMMSL_Settings_get_as_int8(s1, 'key2'), 242424242_YMMSL_int8)
        call assert_eq_int8(YMMSL_Settings_get_as_int8(s1, 'key3'), 42424242424242_YMMSL_int8)
        call assert_eq_logical(YMMSL_Settings_get_as_logical(s1, 'key4'), .false.)
        call assert_eq_real8(YMMSL_Settings_get_as_real8(s1, 'key5'), 13.13d0)
        call assert_eq_real8(YMMSL_Settings_get_as_real8(s1, 'key2'), 242424242d0)
        call YMMSL_Settings_get_as_real8array(s1, 'key6', ra2)
        call assert_eq_real8array(ra2, ra1)
        call YMMSL_Settings_get_as_real8array2(s1, 'key7', ra4)
        call assert_eq_real8array2(ra4, ra3)

        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.set_get_as'
    end subroutine test_settings_set_get_as

    subroutine test_settings_is_a
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1
        real(YMMSL_real8), dimension(2) :: ra1
        real(YMMSL_real8), dimension(3, 2) :: ra2
        logical :: l1
        integer :: err_code

        print *, '[  RUN     ] settings.is_a'
        s1 = YMMSL_Settings_create()
        ra1 = (/1.0d0, 2.0d0/)
        ra2 = reshape((/1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0/), (/3, 2/))

        call YMMSL_Settings_set(s1, 'key1', 'value1')
        call YMMSL_Settings_set(s1, 'key2', 242424242_YMMSL_int4)
        call YMMSL_Settings_set(s1, 'key3', 42424242424242_YMMSL_int8)
        call YMMSL_Settings_set(s1, 'key4', .false.)
        call YMMSL_Settings_set(s1, 'key5', 13.13d0)
        call YMMSL_Settings_set(s1, 'key6', ra1)
        call YMMSL_Settings_set(s1, 'key7', ra2)

        call assert_true(YMMSL_Settings_is_a_character(s1, 'key1'))
        call assert_true(YMMSL_Settings_is_a_int4(s1, 'key2'))
        call assert_true(YMMSL_Settings_is_a_int8(s1, 'key2'))
        call assert_false(YMMSL_Settings_is_a_int4(s1, 'key3'))
        call assert_true(YMMSL_Settings_is_a_int8(s1, 'key3'))
        call assert_true(YMMSL_Settings_is_a_logical(s1, 'key4'))
        call assert_true(YMMSL_Settings_is_a_real8(s1, 'key5'))
        call assert_true(YMMSL_Settings_is_a_real8array(s1, 'key6'))
        call assert_true(YMMSL_Settings_is_a_real8array2(s1, 'key7'))

        call assert_false(YMMSL_Settings_is_a_int8(s1, 'key1'))
        call assert_false(YMMSL_Settings_is_a_logical(s1, 'key1'))
        call assert_false(YMMSL_Settings_is_a_real8(s1, 'key1'))
        call assert_false(YMMSL_Settings_is_a_real8array(s1, 'key1'))
        call assert_false(YMMSL_Settings_is_a_real8array2(s1, 'key1'))

        l1 = YMMSL_Settings_is_a_logical(s1, 'nokey', err_code)
        call assert_eq_integer(err_code, YMMSL_out_of_range)

        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.is_a'
    end subroutine test_settings_is_a

    subroutine test_settings_erase
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1

        print *, '[  RUN     ] settings.erase'
        s1 = YMMSL_Settings_create()
        call YMMSL_Settings_set(s1, 'key', 'value')
        call assert_true(YMMSL_Settings_contains(s1, 'key'))
        call assert_eq_size(YMMSL_Settings_erase(s1, 'key'), 1_YMMSL_size)
        call assert_false(YMMSL_Settings_contains(s1, 'key'))
        call assert_eq_size(YMMSL_Settings_erase(s1, 'key'), 0_YMMSL_size)
        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.erase'
    end subroutine test_settings_erase

    subroutine test_settings_clear
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1

        print *, '[  RUN     ] settings.clear'
        s1 = YMMSL_Settings_create()
        call YMMSL_Settings_set(s1, 'key1', 'value')
        call YMMSL_Settings_set(s1, 'key2', 42_YMMSL_int8)

        call assert_true(YMMSL_Settings_contains(s1, 'key1'))
        call assert_true(YMMSL_Settings_contains(s1, 'key2'))

        call YMMSL_Settings_clear(s1)

        call assert_false(YMMSL_Settings_contains(s1, 'key1'))
        call assert_false(YMMSL_Settings_contains(s1, 'key2'))
        call assert_eq_size(YMMSL_Settings_size(s1), 0_YMMSL_size)

        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.clear'
    end subroutine test_settings_clear

    subroutine test_settings_key
        use ymmsl
        implicit none

        type(YMMSL_Settings) :: s1
        character(len=:), allocatable :: c1, c2, err_msg
        integer :: err_code

        print *, '[  RUN     ] settings.key'
        s1 = YMMSL_Settings_create()
        call YMMSL_Settings_set(s1, 'key1', 'value')
        call YMMSL_Settings_set(s1, 'key2', 42_YMMSL_int8)

        c1 = YMMSL_Settings_key(s1, 1_YMMSL_size, err_code)
        call assert_eq_integer(err_code, YMMSL_success)
        c2 = YMMSL_Settings_key(s1, 2_YMMSL_size, err_code)
        call assert_eq_integer(err_code, YMMSL_success)

        call assert_true((c1 .eq. 'key1') .or. (c1 .eq. 'key2'))
        call assert_true((c2 .eq. 'key1') .or. (c2 .eq. 'key2'))
        call assert_true(c1 .ne. c2)

        c1 = YMMSL_Settings_key(s1, 0_YMMSL_size, err_code, err_msg)
        call assert_eq_integer(err_code, YMMSL_out_of_range)

        c1 = YMMSL_Settings_key(s1, 3_YMMSL_size, err_code)
        call assert_eq_integer(err_code, YMMSL_out_of_range)

        deallocate(c1)
        deallocate(c2)
        call YMMSL_Settings_free(s1)
        print *, '[       OK ] settings.key'
    end subroutine test_settings_key
end module tests

program test_settings
    use tests
    implicit none

    print *, ''
    print *, '[==========] Fortran API Settings'

        call test_settings_create
        call test_settings_equals
        call test_settings_size
        call test_settings_set_get_as
        call test_settings_erase
        call test_settings_clear
        call test_settings_is_a
        call test_settings_key

    print *, '[==========] Fortran API Settings'
    print *, '[  PASSED  ] Fortran API Settings'
    print *, ''
end program test_settings


! --- initialize parameters

      subroutine init_pars()

      use cdat

      implicit real*8 (a-h,o-z)

      open(5,file='IN')
      read(5,*) nbx,nby
      read(5,*) kmax,dt,kout
!      read(5,*) xmin,xmax
      read(5,*) amx,amy
      read(5,*) ax0,ay0,ax,ay,beta
      read(5,*) qx0,qy0
      read(5,*) px0,py0
      read(5,*) idum1
      read(5,*) order
      close(5)

!      ipot = 1 ! 1. HO 2. Morse

      end subroutine