text:
!
! written by Jeff Hammond, Intel Labs, April 2014
! see https://github.com/jeffhammond/nwchem-tce-triples-kernels/blob/master/src/new_sort6.f for history
!
      subroutine new_sort_6(u,s,a,b,c,d,e,f,i,j,k,l,m,n,z)
      implicit none
      integer a,b,c,d,e,f
      integer i,j,k,l,m,n
      double precision s(a*b*c*d*e*f)
      double precision u(a*b*c*d*e*f)
      double precision z ! factor
      integer id(6),jd(6)
      integer ia,ib
      integer j1,j2,j3,j4,j5,j6
      integer jdm,jdl,jdk,jdj
      integer iax,ibx
      jd(1) = a
      jd(2) = b
      jd(3) = c
      jd(4) = d
      jd(5) = e
      jd(6) = f
      if (n.eq.6) then
      jdm=jd(m)
      jdl=jd(l)
      jdk=jd(k)
      jdj=jd(j)
!$omp  parallel do shared(s,u) collapse(3)
!$omp& private(j1,j2,j3,j4,j5,j6,iax,ibx,ia,ib,id)
!$omp& firstprivate(a,b,c,d,e,f,i,j,k,l,m,n,z,jd)
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j4 = 1,d
          do j5 = 1,e
           id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           !id(6) = j6 ! 123456 are only used once and we know n=6, hence ijklm must not be 6
           iax = f*(j5-1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1)))))
           ibx = f*(id(m)-1+jdm*(id(l)-1+jdl*(id(k)-1+jdk
     &                *(id(j)-1+jdj*(id(i)-1)))))
!dec$ vector always nontemporal
!dir$ simd
           do j6 = 1,f
            ia = j6+iax
            ib = j6+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
!$omp end parallel do
      else if (n.eq.5) then
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j4 = 1,d
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           !id(5) = j5
           id(6) = j6
           iax = -1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1))))
           ibx = e*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j5 = 1,e
            ia = j6+f*(j5+iax)
            ib = j5+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.4) then
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           id(3) = j3
           !id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+d*(j3-1+c*(j2-1+b*(j1-1)))
           ibx = d*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j4 = 1,d
            ia = j6+f*(j5-1+e*(j4+iax))
            ib = j4+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.3) then
      do j1 = 1,a
       do j2 = 1,b
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           !id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+c*(j2-1+b*(j1-1))
           ibx = c*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j3 = 1,c
            ia = j6+f*(j5-1+e*(j4-1+d*(j3+iax)))
            ib = j3+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.2) then
      do j1 = 1,a
       do j3 = 1,c
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           !id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+b*(j1-1)
           ibx = b*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j2 = 1,b
            ia = j6+f*(j5-1+e*(j4-1+d*(j3-1+c*(j2+iax))))
            ib = j2+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else ! if (n.eq.1) then
      do j2 = 1,b
       do j3 = 1,c
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           !id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           ibx = a*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j1 = 1,a
            ia = j6+f*(j5-1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1)))))
            ib = j1+ibx
            s(ib) = u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      endif
      return
      end

      subroutine new_sortacc_6(u,s,a,b,c,d,e,f,i,j,k,l,m,n,z)
      implicit none
      integer a,b,c,d,e,f
      integer i,j,k,l,m,n
      double precision s(a*b*c*d*e*f)
      double precision u(a*b*c*d*e*f)
      double precision z
      integer id(6),jd(6)
      integer ia,ib
      integer j1,j2,j3,j4,j5,j6
      integer jdm,jdl,jdk,jdj
      integer iax,ibx
      jd(1) = a
      jd(2) = b
      jd(3) = c
      jd(4) = d
      jd(5) = e
      jd(6) = f
      if (n.eq.6) then
      jdm=jd(m)
      jdl=jd(l)
      jdk=jd(k)
      jdj=jd(j)
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j4 = 1,d
          do j5 = 1,e
           id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           !id(6) = j6 ! 123456 are only used once and we know n=6, hence ijklm must not be 6
           iax = f*(j5-1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1)))))
           ibx = f*(id(m)-1+jdm*(id(l)-1+jdl*(id(k)-1+jdk
     &                *(id(j)-1+jdj*(id(i)-1)))))
           do j6 = 1,f
            ia = j6+iax
            ib = j6+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.5) then
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j4 = 1,d
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           !id(5) = j5
           id(6) = j6
           iax = -1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1))))
           ibx = e*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j5 = 1,e
            ia = j6+f*(j5+iax)
            ib = j5+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.4) then
      do j1 = 1,a
       do j2 = 1,b
        do j3 = 1,c
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           id(3) = j3
           !id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+d*(j3-1+c*(j2-1+b*(j1-1)))
           ibx = d*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j4 = 1,d
            ia = j6+f*(j5-1+e*(j4+iax))
            ib = j4+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.3) then
      do j1 = 1,a
       do j2 = 1,b
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           id(2) = j2
           !id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+c*(j2-1+b*(j1-1))
           ibx = c*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j3 = 1,c
            ia = j6+f*(j5-1+e*(j4-1+d*(j3+iax)))
            ib = j3+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else if (n.eq.2) then
      do j1 = 1,a
       do j3 = 1,c
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           id(1) = j1
           !id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           iax = -1+b*(j1-1)
           ibx = b*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j2 = 1,b
            ia = j6+f*(j5-1+e*(j4-1+d*(j3-1+c*(j2+iax))))
            ib = j2+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      else ! if (n.eq.1) then
      do j2 = 1,b
       do j3 = 1,c
        do j4 = 1,d
         do j5 = 1,e
          do j6 = 1,f
           !id(1) = j1
           id(2) = j2
           id(3) = j3
           id(4) = j4
           id(5) = j5
           id(6) = j6
           ibx = a*(id(m)-1+jd(m)*(id(l)-1+jd(l)*(id(k)-1
     &           +jd(k)*(id(j)-1+jd(j)*(id(i)-1)))))
           do j1 = 1,a
            ia = j6+f*(j5-1+e*(j4-1+d*(j3-1+c*(j2-1+b*(j1-1)))))
            ib = j1+ibx
            s(ib) = s(ib) + u(ia) * z
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
      endif
      return
      end

max_stars_repo_name: tpops/nwchem-tce-triples-kernels

max_issues_repo_name: tpops/nwchem-tce-triples-kernels

max_forks_repo_name: tpops/nwchem-tce-triples-kernels


text:
module foo_module
  interface
    module subroutine foo()
    end subroutine
  end interface
end module

submodule(foo_module) foo_submodule
contains
  module procedure foo
  contains
    subroutine bar()
    end subroutine
  end procedure
end submodule

max_stars_repo_name: rouson/bug-reports

max_issues_repo_name: rouson/bug-reports

max_forks_repo_name: rouson/bug-reports


text:
!----------------------------------------------------------------------------!
!                                                                            !
!  Fortran Coarray MicroBenchmark Suite - Version 1.1                        !
!                                                                            !
!  David Henty, EPCC; d.henty@epcc.ed.ac.uk                                  !
!                                                                            !
!  Copyright 2014 the University of Edinburgh                                !
!                                                                            !
!  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 cafhalo

  use cafparams
  use cafcore
  use cafclock

  implicit none

contains

subroutine cafhaloswap(dimx, dimy, dimz, codimx, codimy, codimz, &
                         cafsynctype, cafhalotype, docheck)

  integer :: dimx, dimy, dimz, codimx, codimy, codimz

  double precision, allocatable, dimension(:,:,:), codimension[:,:,:] :: x

  integer :: coordx, coordy, coordz
  integer :: upcoordx, upcoordy, upcoordz
  integer :: dncoordx, dncoordy, dncoordz
  integer :: upimagex, upimagey, upimagez
  integer :: dnimagex, dnimagey, dnimagez

  double precision :: time1, time2, time

  integer :: i, ndata, nextent, irep, nrep, trialnrep, nsync
  integer :: cafsynctype, cafhalotype

  double precision :: targettime

  integer, parameter :: ndims = 3

  logical :: finished, oktime, docheck, active, lcheck, gcheck
  logical :: hcheck(2*ndims)

  double precision :: myhaloval

  integer, dimension(2*ndims) :: neighbours, neighsync
  double precision, dimension(2*ndims) :: neighhaloval

  integer, dimension(ndims) :: imagecoord

  if (codimx*codimy*codimz /= numimages) then

     if (thisimage == 1) then
        write(*,*) 'cafhaloswap: error - invalid codims for numimages = ', numimages
     end if

     stop

  end if

  allocate(x(0:dimx+1, 0:dimy+1, 0:dimz+1)[codimx, codimy, *])

  targettime = halotargettime

  trialnrep = 1

  ndata = 2*((dimy*dimz) + (dimx*dimz) + (dimx*dimy))

  imagecoord(1:ndims) = this_image(x)

  coordx = imagecoord(1)
  coordy = imagecoord(2)
  coordz = imagecoord(3)

  upcoordx = mod(coordx, codimx) + 1
  upcoordy = mod(coordy, codimy) + 1
  upcoordz = mod(coordz, codimz) + 1

  dncoordx = mod(coordx+codimx-2, codimx) + 1
  dncoordy = mod(coordy+codimy-2, codimy) + 1
  dncoordz = mod(coordz+codimz-2, codimz) + 1

  neighbours(1) = image_index(x, [upcoordx, coordy, coordz])
  neighbours(2) = image_index(x, [coordx, upcoordy, coordz])
  neighbours(3) = image_index(x, [coordx, coordy, upcoordz])
  neighbours(4) = image_index(x, [dncoordx, coordy, coordz])
  neighbours(5) = image_index(x, [coordx, dncoordy, coordz])
  neighbours(6) = image_index(x, [coordx, coordy, dncoordz])

!     write(*,*) 'cafhaloswap: image ', thisimage, 'neighbours = ', (neighbours(i), i = 1, 2*ndims)

  call getneighsync(neighbours, neighsync, nsync)

!     write(*,*) 'cafhaloswap: image ', thisimage, 'neighsync  = ', (neighsync(i), i = 1, nsync)

  active = .true.

  if (thisimage == 1) then

     write(*,*) "cafhaloswap: synchronisation: ", cafsyncname(cafsynctype)
     write(*,*) "cafhaloswap: transfer mode:   ", cafhaloname(cafhalotype)

     if (docheck) then
        write(*,*) 'cafhaloswap: verifying data'
     else
        write(*,*) 'cafhaloswap: NOT verifying data'
     end if

  end if

  oktime = .false.

  ! Because boundary data is copied to multiple halos (eg corner data
  ! is replicated even without diagonal neighbours) we have to set all
  ! the halos to the same value to allow for (simple) verification

  ! Multiply by a big number so we can make these (mostly) unique per
  ! iteration by adding on irep

  myhaloval = thisimage*1000

  neighhaloval(1) = neighbours(1)*1000
  neighhaloval(2) = neighbours(2)*1000
  neighhaloval(3) = neighbours(3)*1000
  neighhaloval(4) = neighbours(4)*1000
  neighhaloval(5) = neighbours(5)*1000
  neighhaloval(6) = neighbours(6)*1000

  do while (.not. oktime)

     lcheck = .true.
     gcheck = .true.
     hcheck(:) = .true.

     x(:,:,:) = 0.0

     x(dimx  , 1:dimy, 1:dimz) = dble(myhaloval)
     x(1:dimx, dimy,   1:dimz) = dble(myhaloval)
     x(1:dimx, 1:dimy, dimz  ) = dble(myhaloval)

     x(1,      1:dimy, 1:dimz) = dble(myhaloval)
     x(1:dimx, 1,      1:dimz) = dble(myhaloval)
     x(1:dimx, 1:dimy, 1     ) = dble(myhaloval)

     sync all

     time1 = caftime()

     do irep = 1, trialnrep

        ! Now set boundaries to iteration-dependent value if we
        ! are verifyiing data

        call cafset2d(x(dimx  , 1:dimy, 1:dimz), 1, dimy*dimz, dimy*dimz, &
                      dble(myhaloval+irep), docheck)

        call cafset2d(x(1:dimx, dimy,   1:dimz), 1, dimx*dimz, dimx*dimz, &
                      dble(myhaloval+irep), docheck)

        call cafset2d(x(1:dimx, 1:dimy, dimz  ), 1, dimx*dimy, dimx*dimy, &
                      dble(myhaloval+irep), docheck)

        call cafset2d(x(1,      1:dimy, 1:dimz), 1, dimy*dimz, dimy*dimz, &
                      dble(myhaloval+irep), docheck)

        call cafset2d(x(1:dimx, 1,      1:dimz), 1, dimx*dimz, dimx*dimz, &
                      dble(myhaloval+irep), docheck)

        call cafset2d(x(1:dimx, 1:dimy, 1     ), 1, dimx*dimy, dimx*dimy, &
                      dble(myhaloval+irep), docheck)

        ! Synchronise so we know we can change halo values on other images
        ! ie wait until halo data has been processed (put) or ensure boundary
        ! data is ready on other images (get)

        call cafdosync(cafsynctype, active, neighsync(1:nsync))

        select case (cafhalotype)

        case(cafmodehalopp)

! !DIR$ PGAS DEFER_SYNC
           x(0,      1:dimy, 1:dimz)[upcoordx, coordy, coordz] = &
                x(dimx,   1:dimy, 1:dimz)

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, 0,      1:dimz)[coordx, upcoordy, coordz] = &
                x(1:dimx, dimy,   1:dimz)

! !DIR$ PGAS DEFER_SYNC
           x(1:dimz, 1:dimy, 0     )[coordx, coordy, upcoordz] = &
                x(1:dimx, 1:dimy, dimz  )

           ! Send bottom boundaries downwards

! !DIR$ PGAS DEFER_SYNC
           x(dimx+1, 1:dimy, 1:dimz)[dncoordx, coordy, coordz] = &
                x(1,      1:dimy, 1:dimz)

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, dimy+1, 1:dimz)[coordx, dncoordy, coordz] = &
                x(1:dimx, 1,      1:dimz)

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, 1:dimy, dimz+1)[coordx, coordy, dncoordz] = &
                x(1:dimx, 1:dimy, 1)

        case(cafmodehalogg)

           ! Get bottom halo from downward neighbours

! !DIR$ PGAS DEFER_SYNC
           x(0,      1:dimy, 1:dimz) = &
                x(dimx,   1:dimy, 1:dimz)[dncoordx, coordy, coordz]

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, 0,      1:dimz) = &
                x(1:dimx, dimy,   1:dimz)[coordx, dncoordy, coordz]

! !DIR$ PGAS DEFER_SYNC
           x(1:dimz, 1:dimy, 0     ) = &
                x(1:dimx, 1:dimy, dimz  )[coordx, coordy, dncoordz]

           ! Get top halo from upward neighbour

! !DIR$ PGAS DEFER_SYNC
           x(dimx+1, 1:dimy, 1:dimz) = &
                x(1,      1:dimy, 1:dimz)[upcoordx, coordy, coordz]

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, dimy+1, 1:dimz) = &
                x(1:dimx, 1,      1:dimz)[coordx, upcoordy, coordz]

! !DIR$ PGAS DEFER_SYNC
           x(1:dimx, 1:dimy, dimz+1) = &
                x(1:dimx, 1:dimy, 1     )[coordx, coordy, upcoordz]

        case default

           write(*,*) 'cafhaloswap: invalid mode ', cafhaloname(cafhalotype)
           stop

        end select

        ! Now wait for all halo data to come in before proceeding (put)
        ! or wait for other images to read boundary data before altering (get)

        call cafdosync(cafsynctype, active, neighsync(1:nsync))

        hcheck(1) = cafcheck2d(x(dimx+1, 1:dimy, 1:dimz), 1, dimy*dimz, &
                               dimy*dimz, dble(neighhaloval(1)+irep), docheck)

        hcheck(2) = cafcheck2d(x(1:dimx, dimy+1, 1:dimz), 1, dimx*dimz, &
                               dimx*dimz, dble(neighhaloval(2)+irep), docheck)

        hcheck(3) = cafcheck2d(x(1:dimx, 1:dimy, dimz+1), 1, dimx*dimy, &
                               dimx*dimy, dble(neighhaloval(3)+irep), docheck)

        hcheck(4) = cafcheck2d(x(0, 1:dimy, 1:dimz), 1, dimy*dimz, &
                               dimy*dimz, dble(neighhaloval(4)+irep), docheck)

        hcheck(5) = cafcheck2d(x(1:dimx, 0, 1:dimz), 1, dimx*dimz,  &
                               dimx*dimz, dble(neighhaloval(5)+irep), docheck)

        hcheck(6) = cafcheck2d(x(1:dimx, 1:dimy, 0), 1, dimx*dimy, &
                               dimx*dimy, dble(neighhaloval(6)+irep), docheck)

        do i = 1, 2*ndims
           lcheck = lcheck .and. hcheck(i)
        end do

        if (docheck) then

           gcheck = gand(lcheck)

           if (gcheck .eqv. .false.) then

              exit

           end if

        end if

     end do

     sync all

     time2 = caftime()
     time = time2 - time1

     call dbcast(time, 1)

     nrep = trialnrep
     oktime = cafchecktime(trialnrep, time, targettime)

     if (gcheck .eqv. .false.) then

        exit

     end if

  end do

  if (thisimage == 1) then

     write(*,*) 'cafhaloswap: ndata, nrep, time = ', &
          ndata, nrep, time

     write(*,*) 'cafhaloswap: latency = ', time/dble(nrep), ' secs'
     write(*,*) 'cafhaloswap: bwidth  = ', &
          (dble(nrep)*dble(ndata)*dble(dblesize)) / &
          (dble(1024*1024)*time), ' MB/s'
     write(*,*)

     if (docheck) then

        if (gcheck) then

           write(*,*) 'cafhaloswap: All results validated'
           write(*,*)

        else

           write(*,*) 'cafhaloswap: ERROR, test FAILED to verify'
           write(*,*)

        end if

     end if

  end if

  deallocate(x)

end subroutine cafhaloswap

end module cafhalo

max_stars_repo_name: davidhenty/cafbench

max_issues_repo_name: davidhenty/cafbench

max_forks_repo_name: davidhenty/cafbench


Module quadratures
  Use math_constants
	Implicit None

Contains
  !This subroutine select the appropiate subroutine for whatever quadrature in terms
  !of the parameters of each quadrature; look inside the function for the parameters
  !order for each quadrature
	! "quad" accepts:
  !                 "CLCU" -> Clenshaw Curtis quadrature
  ! "params_i" -> array with the needed integer parameters
  ! "params_i" -> array with the needed floating point parameters
  ! "qp" -> quadrature points array
  ! "qw" -> quadrature weights array
	Subroutine quad_select(quad,params,qp,qw)
		Implicit None
		character(len=4), intent(in) :: quad !quadrature selector variable
    real * 8, dimension(:), intent(in) :: params !quadrature floating parameters
		real * 8, dimension(:), intent(out) :: qp, qw !quadrature points
    !real * 8, dimension(:), optional, intent(out) :: qw !quadrature weights
    integer :: np
    real * 8 :: xmin,xmax,dr
    xmin = params(1)
    xmax = params(2)
    np   = nint( params(3) )
    select case(quad)
    !Clenshaw Curtis quadrature
    case("CLCU")
      !INTEGER parameters needed order:
      !   Number of points "np"
      !Floating point parameters needed order:
      !   Domain minimum value "xmin"
      !   Domain maximum value "xmax"
      call Clenshaw_Curtis_quadrature(np,qp,qw)
      call rescale_OW(np,qp,qw,-1.d0,1.d0,xmin,xmax)
    !Rectangles quadratures
    case("RECT")
        !INTEGER parameters needed order:
        !   Number of points "np"
        !Floating point parameters needed order:
        !   Domain minimum value "xmin"
        !   Domain maximum value "xmax"
        dr = (xmax-xmin) / dble(np)
        qw = dr
        call rectangles(qp,xmin,dr)
    end select
    return
	End Subroutine

	Subroutine Clenshaw_Curtis_quadrature(order,x,w)
		Implicit None
    Integer, intent(in) :: order
    Real * 8, dimension(order), intent(out) :: x, w
	  Real * 8 :: b
	  Integer :: i
	  Integer :: j
	  Real * 8 :: theta
	  If ( order < 1 ) Then
	    write ( *, '(a)' ) ' '
	    write ( *, '(a)' ) 'CLENSHAW_CURTIS_COMPUTE - Fatal error!'
	    write ( *, '(a,i8)' ) '  Illegal value of ORDER = ', order
	    stop
	  End If
	  If ( order == 1 ) Then
	    x(1) = 0.0D+00
	    w(1) = 2.0D+00
	    return
	  End If
		Do i = 1, order
	    x(i) = cos ( dble ( order - i ) * pi / dble ( order - 1 ) )
	  End Do
	  x(1) = -1.0D+00
	  If ( mod ( order, 2 ) == 1 ) Then
	    x((order+1)/2) = 0.0D+00
	  End If
	  x(order) = +1.0D+00
	  Do i = 1, order
	    theta = dble( i - 1 ) * pi / dble( order - 1 )
	    w(i) = 1.0D+00
	    Do j = 1, ( order - 1 ) / 2
	      If ( 2 * j == ( order - 1 ) ) Then
	        b = 1.0D+00
	      Else
	        b = 2.0D+00
	      End If
	      w(i) = w(i) - b * cos ( 2.0D+00 * dble( j ) * theta ) / dble( 4 * j * j - 1 )
	    End Do
	  End Do
	  w(1) = w(1) / dble( order - 1 )
	  w(2:order-1) = 2.0D+00 * w(2:order-1) / dble( order - 1)
	  w(order) = w(order) / dble( order - 1 )
	  return
	End Subroutine

  Subroutine rectangles(x,x_min,dr)
		Implicit None
	  Real * 8, intent(out), dimension(:) :: x
    Real * 8, intent(in) :: dr,x_min
    integer :: i1
	  Do i1=1, size(x)
      x(i1) = x_min + i1 * dr
    End Do
	  return
	End Subroutine

  !This subroutines rescales the abscissas and weights values from a range of
  !x ϵ [xa, xb] to y ϵ [ya, yb] and obtains the new weights in the new arrays
  !y, wy
  Subroutine rescale_NA(n,x,wx,xa,xb,y,wy,ya,yb)
  	implicit none
  	integer, intent(in) :: n
  	real * 8, intent(in) :: xa,xb,ya,yb
  	real * 8, dimension(n), intent(in) :: x,wx
  	real * 8, dimension(n), intent(out) :: y,wy
  	y(1:n) = ya  + (( yb - ya ) * (x(1:n)-xa) / ( xb - xa ))
    wy(1:n) = ( yb - ya ) * wx(1:n) / ( xb - xa)
  End Subroutine

  !This subroutines rescales the abscissas and weights values from a range of
  !x ϵ [xa, xb] to y ϵ [ya, yb] and obtains the new wy in the arrays x,wx
  Subroutine rescale_OW(n,x,wx,xa,xb,ya,yb)
  	implicit none
  	integer, intent(in) :: n
  	real * 8, intent(in) :: xa,xb,ya,yb
  	real * 8, dimension(n), intent(inout) :: x,wx
  	x(1:n) = ( ( ya + yb ) + ( yb - ya ) * x(1:n) ) / ( xb - xa )
    wx(1:n) = ( yb - ya ) * wx(1:n) / ( xb - xa)
  End Subroutine

  !This subroutine writes down the information of a given quadrature used
	Subroutine quad_writting(quad,i_unit,params)
    Implicit None
    Real * 8, dimension(:), intent(in) :: params
    Character(len=4), intent(in) :: quad
    Integer, intent(in) :: i_unit
    Select case (quad)
			case("CLCU")
				Write(i_unit,*) "				Clenshaw-Curtis"
				Write(i_unit,*) "quadrature minimum value    = ", params(1)
				Write(i_unit,*) "quadrature maximum value    = ", params(2)
        Write(i_unit,*) "quadrature number of points = ", params(3)
      case("RECT")
        Write(i_unit,*) "				Rectangles"
        Write(i_unit,*) "quadrature minimum value    = ", params(1)
        Write(i_unit,*) "quadrature differential     = ", params(2)
			case Default
				Write(i_unit,*) "				"//quad
        Write(i_unit,*) "quadrature minimum value    = ", params(1)
				Write(i_unit,*) "quadrature maximum value    = ", params(2)
        Write(i_unit,*) "quadrature number of points = ", params(3)
		End Select
		Write(i_unit,*) "###############################"
  End Subroutine

  !This subroutine generates the quadrature points and weights of
  !up to Simpson 3/8 quadrature given an order greater than 0
  Subroutine Simpson_3_8_quad(qp,qw,order,qpmin,qpmax)
    Implicit None
    Real * 8, intent(out), dimension(:) :: qp, qw !Quadrature points and quadrature weights
    Real * 8, intent(in) :: qpmin, qpmax
    Integer, intent(in) :: order !Order of the quadrature
    Real * 8 :: dx
    Integer :: i1, i2, ratio, mod_ratio, imax
		dx = (qpmax-qpmin) / dble(order)
    !Calculating the quadrature points
    Do i1 = 1, order
      qp(i1) = qpmin + (i1 * dx)
    End Do
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    If (order<1) Then
      Print *, 'Error at "Simpson_3_8_quad", not a valid order '
    Else If (order == 1) Then
      qw(:) = dx
      return
    Else If (order==2) Then
      qw(:) =  0.5d0 * dx
			return
		Else If (order==3) Then
      qw(1) = 1.d0
      qw(2) = 4.d0
      qw(3) = 1.d0
      qw(:) = qw(:) * dx / 3.d0
			return
		Else If (order==4) Then
      qw(1) = 1.d0
      qw(2) = 3.d0
      qw(3) = 3.d0
      qw(4) = 1.d0
      qw(:) = qw(:) * dx * 0.375d0
			return
		Else If (order==5) Then
      qw(1) = 1.d0
      qw(2) = 4.d0
      qw(3) = 2.d0
      qw(4) = 4.d0
      qw(5) = 1.d0
      qw(:) = qw(:) * dx / 3.d0
			return
		Else If (order>5) Then
      ratio = order/3
			mod_ratio = mod(order,3)
      qw(1) = 1.d0
      Do i1=1, ratio
        i2 = (i1 * 3) - 1
        qw(i2)   = 3.d0
        qw(i2+1) = 3.d0
        if(i2+2 < order) qw(i2+2) = 2.d0
      End Do
      qw(:) = qw(:) * 3.d0 * dx / 8.d0
			If	(mod_ratio == 0) Then
				qw(order-2) = ((3.d0/8.d0) + (1.d0/3.d0)) * dx
        qw(order-1) = 4.d0 * dx / 3.d0
        qw(order)   = 1.d0 * dx / 3.d0
        return
      Else If (mod_ratio==1) Then
        qw(order-1) = (3.d0 / 8.d0 + 0.5d0) * dx
        qw(order)   =  0.5d0 * dx
        return
      Else If (mod_ratio==2) Then
				qw(order) = (3.d0/8.d0) * dx
        return
      End If
    End If
  End Subroutine

End Module


text:
      FUNCTION beta(z,w)
      REAL beta,w,z
CU    USES gammln
      REAL gammln
      beta=exp(gammln(z)+gammln(w)-gammln(z+w))
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/beta.for


text:
      subroutine write_to_file(outID, prn_pick, prn, elev,
     . azimuth, tod, s1,s2, bad_point,l1,l2,grange,
     .  blockIIRM,L2Ctracking,s5)
      implicit none
      include 'local.inc'
      real*8 s1, s2, tod, elev, azimuth, x,y, l1, l2,w1,w2,
     .  grange, s5
      logical bad_point
      integer prn_pick,outID,L2Ctracking,prn
      integer blockIIRM(maxsat)
c     author: Kristine Larson
c     16jul15 added s5
c     this supports old option that printed out L1 and L2 phase
      w1 = c/f1 ! L1 wavelength in meters
      w2 = c/f2 !L2 wavelength in meters
c     no longer writing out reflection point. columns 5 and 6 are zero
      x = 0
      y = 0

      if (prn_pick.eq.99.and.(elev.gt.5
     .  .and.elev.lt.30)) then
          write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 3f7.2)' )
     .    prn, elev, azimuth, tod, x,y, s1, s2, s5
      elseif ( prn_pick.eq.50 .and.  elev.lt.10 ) then
        write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 2f7.2)' )
     .    prn, elev, azimuth, tod, 0.d0,0.d0, s1, s2
c     all data above 5 degrees
      elseif ( prn_pick.eq.88 .and.  elev.gt.5 ) then
        write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 3f7.2)' )
     .    prn, elev, azimuth, tod, x,y, s1, s2,s5
c     everything < 30
      elseif (  prn_pick.eq.66 .and. elev.lt.30 ) then
        write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 3f7.2)' )
     .    prn, elev, azimuth, tod, 0.d0, 0.d0, s1, s2, s5
c       L2C data only - LLI indicator has to be a zero
c       assumes if no phase data, then it is not a good snr value
      elseif ((prn_pick.eq.77).and.
     .     (blockIIRM(prn).eq.1.and.elev.ge.5))then
        if (.not.bad_point.and.(L2Ctracking.eq.0) ) then
          write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 3f7.2)' )
     .        prn, elev, azimuth, tod, x,y,s1, s2,s5
        endif
c     if you only requested a single satellite
      elseif (prn.eq.prn_pick)  then
        write(outID,'(i3,  2f10.4, f10.0, 2f7.2, 3f7.2)' )
     .    prn, elev, azimuth, tod, x,y, s1, s2, s5
      endif
      end

ext: f


text:
program test_open
use stdlib_experimental_io, only: open
use stdlib_experimental_error, only: assert
implicit none

character(:), allocatable :: filename
integer :: io, u, a(3)

! Text file
filename = get_outpath() // "/io_open.dat"

! Test mode "w"
u = open(filename, "w")
write(u, *) 1, 2, 3
close(u)

! Test mode "r"
u = open(filename, "r")
read(u, *) a
call assert(all(a == [1, 2, 3]))
close(u)

! Test mode "a"
u = open(filename, "a")
write(u, *) 4, 5, 6
close(u)
u = open(filename, "r")
read(u, *) a
call assert(all(a == [1, 2, 3]))
read(u, *) a
call assert(all(a == [4, 5, 6]))
close(u)

! Stream file
filename = get_outpath() // "/io_open.stream"

! Test mode "w"
u = open(filename, "wb")
write(u) 1, 2, 3
close(u)

! Test mode "r"
u = open(filename, "rb")
read(u) a
call assert(all(a == [1, 2, 3]))
close(u)

! Test mode "a"
u = open(filename, "ab")
write(u) 4, 5, 6
close(u)
u = open(filename, "rb")
read(u) a
call assert(all(a == [1, 2, 3]))
read(u) a
call assert(all(a == [4, 5, 6]))
close(u)

!0 and non-0 open
filename = get_outpath() // "/io_open.stream"

u = open(filename, "rb", io)
call assert(io == 0)
if (io == 0) close(u)

u = open(filename, "ab", io)
call assert(io == 0)
if (io == 0) close(u)

filename = get_outpath() // "/does_not_exist.error"

u = open(filename, "a", io)
call assert(io /= 0)

u = open(filename, "r", io)
call assert(io /= 0)

contains

    function get_outpath() result(outpath)
    integer :: ierr
    character(256) :: argv
    character(:), allocatable :: outpath

    call get_command_argument(1, argv, status=ierr)
    if (ierr==0) then
        outpath = trim(argv)
    else
        outpath = '.'
    endif
    end function get_outpath

end program

max_stars_repo_name: nncarlson/stdlib

max_issues_repo_name: nncarlson/stdlib

max_forks_repo_name: nncarlson/stdlib


text:
! { dg-do run }
! { dg-options "-std=legacy" }
!
! PR 27784 - Different strings should compare unequal even if they
!            have CHAR(0) in them.

      program main
      character*3 str1, str2
      call setval(str1, str2)
      if (str1 == str2) call abort
      end

      subroutine setval(str1, str2)
      character*3 str1, str2
      str1 = 'a' // CHAR(0) // 'a'
      str2 = 'a' // CHAR(0) // 'c'
      end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!! https://github.com/flang-compiler/flang/issues/244
!!
!! This is a derived type-valued version of flang-20170924c.f90
!!
!! The following example should terminate normally, but instead
!! exits at the STOP 3 statement.
!!
!! $ flang flang-bug-20190503b.f90
!! $ ./a.out
!!   3

type point
  real a, b
end type

type foo
  class(*), allocatable :: val
end type
type(foo) :: x, y

allocate(x%val, source=point(1.0,2.0))
y = x

select type (p => x%val)
type is (point)
  if (p%a /= 1.0 .or. p%b /= 2.0) stop 1
class default
  stop 2
end select

select type (p => y%val)
type is (point)
  if (p%a /= 1.0 .or. p%b /= 2.0) stop 3
class default
  stop 4
end select

end

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
! Contains routines to evaluate XC energy and potential for radial charge distributions
module radial_xc

  use xc_f90_types_m
  use xc_f90_lib_m

  implicit none

  ! Numerical flag choosing functional type
  integer :: flag_functional_type
  character(len=120)  :: functional_description  ! DRB lengthened to contain LibXC names
  integer, parameter :: functional_lda_pz81        = 1
  integer, parameter :: functional_lda_gth96       = 2
  integer, parameter :: functional_lda_pw92        = 3    ! PRB 45, 13244 (1992) + PRL 45, 566 (1980)
  integer, parameter :: functional_xalpha          = 4    ! Slater/Dirac exchange only  ; no correlation
  integer, parameter :: functional_hartree_fock    = 10   ! Hartree-Fock exact exchange ; no correlation

  integer, parameter :: functional_gga_pbe96       = 101  ! Standard PBE
  integer, parameter :: functional_gga_pbe96_rev98 = 102  ! revPBE (PBE + Zhang-Yang 1998)
  integer, parameter :: functional_gga_pbe96_r99   = 103  ! RPBE   (PBE + Hammer-Hansen-Norskov 1999)
  integer, parameter :: functional_gga_pbe96_wc    = 104  ! WC   (Wu-Cohen 2006)

  integer, parameter :: functional_hyb_pbe0        = 201  ! PBE0   (hybrid PBE with exx_alpha=0.25)

  ! LibXC variables
  integer :: n_xc_terms
  integer, dimension(2) :: i_xc_family
  type(xc_f90_pointer_t), dimension(2) :: xc_func, xc_info
  logical :: flag_use_libxc

contains

  subroutine init_xc

    use global_module, ONLY : nspin, flag_dft_d2, iprint
    use GenComms, ONLY : inode, ionode, cq_abort
    use numbers

    implicit none

    ! Local variables
    integer :: vmajor, vminor, vmicro, i, j
    integer, dimension(2) :: xcpart
    character(len=120) :: name, kind, family, ref
    type(xc_f90_pointer_t) :: temp_xc_func
    type(xc_f90_pointer_t) :: temp_xc_info

    ! Test for LibXC or CQ
    if(flag_functional_type<0) then
       ! --------------------------
       ! LibXC functional specified
       ! --------------------------
       flag_use_libxc = .true.
       call xc_f90_version(vmajor, vminor, vmicro)
       if(inode==ionode.AND.iprint>0) then
          if(vmajor>2) then
             write(*,'("LibXC version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
          else
             write(*,'("LibXC version: ",I1,".",I1)') vmajor, vminor
          end if
       end if
       ! Identify the functional
       if(-flag_functional_type<1000) then ! Only exchange OR combined exchange-correlation
          n_xc_terms = 1
          xcpart(1) = -flag_functional_type
       else ! Separate the two parts
          n_xc_terms = 2
          ! Make exchange first, correlation second for consistency
          i = floor(-flag_functional_type/1000.0_double)
          ! Temporary init to find exchange or correlation
          if(nspin==1) then
             call xc_f90_func_init(temp_xc_func, temp_xc_info, i, XC_UNPOLARIZED)
          else if(nspin==2) then
             call xc_f90_func_init(temp_xc_func, temp_xc_info, i, XC_POLARIZED)
          end if
          select case(xc_f90_info_kind(temp_xc_info))
          case(XC_EXCHANGE)
             xcpart(1) = i
             xcpart(2) = -flag_functional_type - xcpart(1)*1000
          case(XC_CORRELATION)
             xcpart(2) = i
             xcpart(1) = -flag_functional_type - xcpart(2)*1000
          end select
          call xc_f90_func_end(temp_xc_func)
       end if
       ! Now initialise and output
       do i=1,n_xc_terms
          if(nspin==1) then
             call xc_f90_func_init(xc_func(i), xc_info(i), xcpart(i), XC_UNPOLARIZED)
          else if(nspin==2) then
             call xc_f90_func_init(xc_func(i), xc_info(i), xcpart(i), XC_POLARIZED)
          end if
          ! Consistent threshold with Conquest
          !if(vmajor>2) call xc_f90_func_set_dens_threshold(xc_func(i),RD_ERR)
          call xc_f90_info_name(xc_info(i), name)
          i_xc_family(i) = xc_f90_info_family(xc_info(i))
          if(inode==ionode) then
             select case(xc_f90_info_kind(xc_info(i)))
             case (XC_EXCHANGE)
                write(kind, '(a)') 'an exchange functional'
             case (XC_CORRELATION)
                write(kind, '(a)') 'a correlation functional'
             case (XC_EXCHANGE_CORRELATION)
                write(kind, '(a)') 'an exchange-correlation functional'
             case (XC_KINETIC)
                write(kind, '(a)') 'a kinetic energy functional'
             case default
                write(kind, '(a)') 'of unknown kind'
             end select

             select case (i_xc_family(i))
             case (XC_FAMILY_LDA);
                write(family,'(a)') "LDA"
             case (XC_FAMILY_GGA);
                write(family,'(a)') "GGA"
             case (XC_FAMILY_HYB_GGA);
                write(family,'(a)') "Hybrid GGA"
             case (XC_FAMILY_MGGA);
                write(family,'(a)') "MGGA"
             case (XC_FAMILY_HYB_MGGA);
                write(family,'(a)') "Hybrid MGGA"
             case default;
                write(family,'(a)') "unknown"
             end select
             functional_description = trim(name)
             if(iprint>2) then
                if(vmajor>2) then
                   write(*,'("The functional ", a, " is ", a, ", it belongs to the ", a, &
                        " family and is defined in the reference(s):")') &
                        trim(name), trim(kind), trim(family)
                   j = 0
                   call xc_f90_info_refs(xc_info(i), j, ref)
                   do while(j >= 0)
                      write(*, '(a,i1,2a)') '[', j, '] ', trim(ref)
                      call xc_f90_info_refs(xc_info(i), j, ref)
                   end do
                else
                   write(*,'("The functional ", a, " is ", a, ", and it belongs to the ", a, &
                        " family")') &
                        trim(name), trim(kind), trim(family)
                end if
             else if(iprint>0) then
                write(*,'(2x,"Using the ",a," functional ",a)') trim(family),trim(name)
             else
                write(*,fmt='(2x,"Using functional ",a)') trim(name)
             end if
          end if
       end do
    else
       ! -----------------------------
       ! Conquest functional specified
       ! -----------------------------
       flag_use_libxc = .false.
       if(nspin==2) then ! Check for spin-compatible functionals
          if ( flag_functional_type == functional_lda_pz81       .or. &
               flag_functional_type == functional_lda_gth96     ) then
             if (inode == ionode) &
                  write (*,'(/,a,/)') &
                  '*** WARNING: the chosen xc-functional is not &
                  &implemented for spin polarised calculation, &
                  &reverting to LDA-PW92. ***'
             flag_functional_type = functional_lda_pw92
          end if
       end if
       select case(flag_functional_type)
       case (functional_lda_pz81)
          functional_description = 'LDA PZ81'
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case (functional_lda_gth96)
          functional_description = 'LDA GTH96'
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case (functional_lda_pw92)
          functional_description = 'LSDA PW92'
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case (functional_gga_pbe96)
          functional_description = 'GGA PBE96'
       case (functional_gga_pbe96_rev98)            ! This is PBE with the parameter correction
          functional_description = 'GGA revPBE98'   !   in Zhang & Yang, PRL 80:4, 890 (1998)
       case (functional_gga_pbe96_r99)              ! This is PBE with the functional form redefinition
          functional_description = 'GGA RPBE99'     !   in Hammer et al., PRB 59:11, 7413-7421 (1999)
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case (functional_gga_pbe96_wc)               ! Wu-Cohen nonempirical GGA functional
          functional_description = 'GGA WC'         !   in Wu and Cohen, PRB 73. 235116, (2006)
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case (functional_hyb_pbe0)                   ! This is PB0E with the functional form redefinition
          functional_description = 'hyb PBE0'
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       case default
          functional_description = 'LSDA PW92'
          if(flag_dft_d2) call cq_abort("DFT-D2 only compatible with PBE and rPBE")
       end select
       if(inode==ionode) write(*,'(/10x, "The functional used will be ", a15)') functional_description
    end if ! if selecting LibXC or CQ
  end subroutine init_xc
  !!***

  subroutine get_vxc(n_tot,rr,rho,vxc,exc)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort
    use mesh, ONLY: make_derivatives

    implicit none

    ! Passed variables
    integer :: n_tot
    real(double), dimension(n_tot) :: rho, rr,vxc
    real(double), OPTIONAL :: exc
    real(double), allocatable, dimension(:) :: exc_array

    integer :: i, n
    logical :: flag_energy, flag_libxc ! This will be in global later
    real(double), dimension(:), allocatable :: drho_dr, sigma, vrho, vsigma, loc_rho, d2term

    flag_energy = .false.
    ! Necessary for LibXC
    allocate(exc_array(n_tot))
    exc_array = zero
    if(PRESENT(exc)) then
       flag_energy = .true.
       exc = zero
    end if
    vxc = zero
    ! Choose LibXC or Conquest
    if(flag_use_libxc) then
       allocate(drho_dr(n_tot),sigma(n_tot),vrho(n_tot),vsigma(n_tot),loc_rho(n_tot),d2term(n_tot))
       ! Make derivatives
       loc_rho = rho/(four*pi)
       call make_derivatives(loc_rho, drho_dr, sigma)
       do n=1,n_xc_terms
          ! Call routine
          select case (i_xc_family(n))
          case(XC_FAMILY_LDA)
             call xc_f90_lda_exc_vxc(xc_func(n),n_tot,loc_rho(1),exc_array(1),vrho(1))
             vxc = vxc + vrho
          case(XC_FAMILY_GGA)
             call xc_f90_gga_exc_vxc(xc_func(n),n_tot,loc_rho(1),sigma(1),exc_array(1),vrho(1),vsigma(1))
             vxc = vxc + vrho
             d2term = zero
             vsigma = vsigma*two*drho_dr
             call make_derivatives(vsigma,d2term)!drho_dr) ! Re-use variable - we only need sigma
             vxc = vxc - (two*vsigma/rr + d2term)!drho_dr) ! Radial part of div.(df/dg)
          end select
          if(PRESENT(exc)) then
             do i = 1, n_tot
                exc = exc + exc_array(i)
             end do
             ! Potentially also find Exc correction
          end if
       end do
       deallocate(drho_dr,sigma,vrho,vsigma)
    else
       select case(flag_functional_type)
       case(functional_lda_pz81)
          if(flag_energy) then
             call vxc_pz_ca(n_tot, rr, rho, vxc, exc_array)
          else
             call vxc_pz_ca(n_tot, rr, rho, vxc)
          end if
       case(functional_gga_pbe96)
          if(flag_energy) then
             call vxc_gga_pbe(n_tot, rr, rho, vxc, exc_array)
          else
             call vxc_gga_pbe(n_tot, rr, rho, vxc)
          end if
       end select
       if(PRESENT(exc)) then
          exc = zero
          do i = 1, n_tot
             exc = exc + exc_array(i)
          end do
          ! Potentially also find Exc correction
       end if
    end if
    deallocate(exc_array)
  end subroutine get_vxc

  subroutine vxc_pz_ca(n_tot,rr,rho,vxc,exc)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    implicit none

    ! Passed variables
    integer :: n_tot
    real(double), dimension(n_tot) :: rho, rr,vxc
    real(double), OPTIONAL, dimension(n_tot) :: exc

    ! Local variables
    integer :: i
    real(double) :: dx_sq_over_twelve, qtot, V0
    real(double), parameter :: prefac_rs = (three/fourpi)**third
    real(double), parameter :: prefac_vx = (three_halves/pi)**(two*third)
    real(double), dimension(:), allocatable :: y

    real(double) :: denominator, e_correlation, e_exchange, ln_rs, &
                    numerator, rcp_rs, rh, rs, rs_ln_rs, sq_rs

    real(double), parameter :: alpha  = -0.45817_double
    real(double), parameter :: beta_1 =  1.0529_double
    real(double), parameter :: beta_2 =  0.3334_double
    real(double), parameter :: gamma  = -0.1423_double
    real(double), parameter :: p =  0.0311_double
    real(double), parameter :: q = -0.048_double
    real(double), parameter :: r =  0.0020_double
    real(double), parameter :: s = -0.0116_double

    do i=1,n_tot
       !  rho is 4.pi.rho
       rh = rho(i)/fourpi
       if(rh>1e-8_double) then
          rs = prefac_rs*rh**(-third)
          if(rs>one) then
             sq_rs = sqrt(rs)
             denominator = (one + beta_1 * sq_rs + beta_2 * rs)
             vxc(i) = -prefac_vx/rs + gamma*(one + seven_sixths*beta_1*sq_rs + &
                  four_thirds*beta_2*rs)/(denominator*denominator)
             if(present(exc)) exc(i) = - 0.75_double*prefac_vx/rs + gamma/denominator
          else
             ln_rs = log(rs)
             vxc(i) = -prefac_vx/rs + p*ln_rs + (q-third*p) + two*third*r*rs*ln_rs + third*(two*s - r)*rs
             if(present(exc)) exc(i) = - 0.75_double*prefac_vx/rs + p*ln_rs + q + r*rs*ln_rs + s*rs
          end if
       else
          vxc(i) = zero
          if(present(exc)) exc(i) = zero
       end if
    end do
  end subroutine vxc_pz_ca

  ! It's not necessary to calculate the full set of arrays of derivatives etc for this simple GGA
  ! implementation - we can simply do it point by point (as it's one dimensional).
  subroutine vxc_gga_pbe(n_tot,rr,rho,vxc,exc)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort
    use mesh, ONLY: drdi, alpha, drdi_squared

    implicit none

    ! Passed variables
    integer :: n_tot
    real(double), dimension(n_tot) :: rho, rr,vxc
    real(double), dimension(n_tot), OPTIONAL :: exc

    ! Local variables
    ! rho/4pi, differentials wrt grid (n) and radius (r) and laplacian
    real(double) :: rho_sc, drho_dn, d2rho_dn2, drho_dr, d2rho_dr2, lapl_rho
    real(double), parameter :: thirty = 30.0_double
    real(double), parameter :: prefac_rs = (three/fourpi)**third
    real(double), parameter :: prefac_kf = (three *pi*pi)**third
    ! Parameters for GGA routines
    real(double) :: kF, s, t, u, v, rs, zeta, k_s
    ! GGA exchange energy and potential
    real(double) :: ex, vx, ec, vc, h, dvc

    integer :: i

    ! Work on main part first (where we can use FD) and add ends later
    do i=3,n_tot-2
       !  rho is 4.pi.rho
       rho_sc = rho(i)/fourpi
       ! Five point finite differences for differential wrt mesh points (not r)
       ! First derivative - extra factor of 4pi to scale rho
       drho_dn   = ( rho(i-2) -   eight*rho(i-1) + eight*rho(i+1) - rho(i+2))/(twelve*fourpi)
       ! Second derivative - extra factor of 4pi to scale rho
       d2rho_dn2 = (-rho(i-2) + sixteen*rho(i-1) -thirty*rho(i) + sixteen*rho(i+1) - rho(i+2))/(twelve*fourpi)
       ! Now differential wrt r
       drho_dr   = drho_dn/(drdi(i))
       ! NB the alpha here is from the radial mesh (NOT the alpha defined in the PW92 LDA correlation)
       d2rho_dr2 = (d2rho_dn2 - alpha*drho_dn)/drdi_squared(i)
       lapl_rho  = (d2rho_dn2 + alpha*drho_dn)/drdi_squared(i)
       if(rho_sc>1e-10_double) then ! Hamann's tolerance - check
          kF = prefac_kF*rho_sc**third
          s = abs(drho_dr)/(two * kF * rho_sc)
          u = abs(drho_dr)*d2rho_dr2/(rho_sc*rho_sc*eight*kF*kF*kF)
          v = lapl_rho/(rho_sc*four*kF*kF)
          call pbe_exch(rho_sc,s,u,v,ex,vx)
          rs = prefac_rs*rho_sc**(-third)
          k_s = sqrt(four*kF/pi) ! NB a0 = 1 as we use atomic units
          zeta = zero
          t = abs(drho_dr)/(two*k_s*rho_sc) ! Set phi = one as zeta = zero
          ! Redefine u and v with k_s instead of k_F
          u = abs(drho_dr)*d2rho_dr2/(rho_sc*rho_sc*eight*k_s*k_s*k_s)
          v = lapl_rho/(rho_sc*four*k_s*k_s)
          call pbe_corr(rs,t,u,v,ec,vc,h,dvc) ! No spin
       else
          vx = zero
          vc = zero
          dvc = zero
          ex = zero
          ec = zero
          h = zero
       end if
       vxc(i) = vx + vc + dvc
       if(PRESENT(exc)) exc(i) = ex + ec + h
    end do
    ! Assume that the end points are essentially constant
    vxc(1) = vxc(3)
    vxc(2) = vxc(3)
    vxc(n_tot-1) = vxc(n_tot-2)
    vxc(n_tot)   = vxc(n_tot-2)
    if(PRESENT(exc)) then
       exc(1) = exc(3)
       exc(2) = exc(3)
       exc(n_tot-1) = exc(n_tot-2)
       exc(n_tot)   = exc(n_tot-2)
    end if
    return
  end subroutine vxc_gga_pbe

  ! From the original PBE routines
  subroutine pbe_exch(rho, s, u, v, ex, vx)

    use datatypes
    use numbers

    implicit none

    ! Passed variables
    real(double) :: rho, s, u, v, ex, vx

    ! Precalculated constants
    real(double), parameter :: prefac_ex = -three_quarters * (3/pi)**third
    real(double), parameter :: mu = 0.2195149727645171_double ! beta (pi^2)/3 (Eq 12 in paper)
    real(double), parameter :: kappa = 0.804_double

    ! Local variables
    real(double) :: ex_lda, ul, p0, Fx, dFds, d2Fds2

    ul = mu/kappa
    ! ------
    ! Energy
    ! ------
    ! LDA exchange
    ex_lda = prefac_ex * rho**third
    ! PBE enhancement
    p0 = one + ul * s * s
    Fx = one + kappa - kappa / p0
    ex = ex_lda * Fx
    ! ---------
    ! Potential
    ! ---------
    ! Derivatives of Fx wrt s: dFds = (1/s)dF/ds; d2Fds2 = d(dFds)/ds
    dFds = two * mu / (p0 * p0)
    d2Fds2 = -four * ul * s * dFds / p0
    vx = ex_lda * (four * third * Fx - (u - four * third * s * s * s)*d2Fds2 - v * dFds)
  end subroutine pbe_exch

  subroutine pbe_corr(rs, t, u, v, ec_unif, vc_unif, h, dvc)

    use datatypes
    use numbers

    implicit none

    ! Passed variables
    real(double) :: rs, t, u, v, ec_unif, vc_unif, h, dvc ! Seitz radius !

    ! Local variables
    real(double) :: sq_rs, prefac_ec_lda, log_fac_ec_lda, denominator
    real(double) :: d_prefac_ec_lda, d_log_fac_ec_lda, d_denom_drs
    real(double) :: t2, t4, t6
    real(double) :: F1, F2, A, F3, F4, F6, ec_rs
    real(double) :: A2, bg, bec, F5, F8, F9, h_A, h_rs, fact0, fact1, h_BT, h_RST, h_T, fact2, fact3, h_TT
    ! Precalculated constants
    real(double), parameter :: prefac_ex = -three_quarters * (3/pi)**third
    real(double), parameter :: k02 = -0.062182_double       ! -2*A
    real(double), parameter :: k03 = -0.0132882934_double   ! -2*A*alpha1
    real(double), parameter :: k04 = 0.4723158174_double    ! 2*A*beta1
    real(double), parameter :: k05 = 0.2230841432_double    ! 2*A*beta2
    real(double), parameter :: k06 = 0.1018665524_double    ! 2*A*beta3
    real(double), parameter :: k07 = 0.03065199508_double   ! 2*A*beta4
    real(double), parameter :: k08 = -0.008858862267_double ! 2*k03/3
    real(double), parameter :: k09 = 0.0787193029_double    ! k04/6
    real(double), parameter :: k10 = 0.074361381067_double  ! k05/3
    real(double), parameter :: k11 = 0.0509332762_double    ! k06/2
    real(double), parameter :: k12 = 0.0204346633867_double ! 2*k07/3
    real(double), parameter :: beta = 0.066725_double
    real(double), parameter :: gamma = 0.031091_double
    real(double), parameter :: beta_gamma = 2.146119_double    ! beta/gamma

    sq_rs = sqrt(rs)

    ! Local correlation energy Eq 10 PRB 45, 13244
    prefac_ec_lda = k02 + k03*rs
    denominator = sq_rs * ( k04 + sq_rs * ( k05 + sq_rs * ( k06 + k07 * sq_rs)))
    if (denominator > zero) then
       log_fac_ec_lda = log( one + one/denominator )
    else
       log_fac_ec_lda = 0
    end if
    ec_unif = prefac_ec_lda * log_fac_ec_lda

    ! Non-local correlation energy
    ! Original Conquest routines adapted from PRL 77 3865 and PRB 45 13244
    F1 = ec_unif / gamma ! -PON
    F2 = exp(-F1)
    ! This is the variable A from the PBE paper Eq 8, NOT the constant in PW92 !
    A = beta_gamma / (F2 - one + RD_ERR) ! B in Burke's code
    t2 = t*t
    t4 = t2*t2
    F3 = t2 + A * t4 ! Q5 = 1+A*F3 in Burke; Q4 * t2 = F3 in Burke
    F6 = one + A * F3 !Q5
    F4 = beta_gamma * F3 / F6 !(one + A * F3)
    ! Eq 7 with phi = 1, atomic units
    h = gamma * log(one + F4)

    ! Local correlation potential
    ! ECRS - dec_drs for LDA from PW92 Appendix A
    ! rtrs == sq_rs
    ! Q0 == prefac_ec_lda
    ! Q1 == denominator
    ! Q2 == log_fac_ec_lda
    ! Q3 = A*(B1/rtrs + two*B2 + rtrs*(three*B3 + four*B4*rtrs))
    !    = half*sq_rs*(2A*B1 + 2*2A*B2*sq_rs + 3*2A*B3*rs + 4*2A*B4*rs*sq_rs  )/rs
    !    = half*(k04 + sq_rs*(two*k05 + sq_rs*(three*k06 + four*k07*sq_rs)))/sq_rs
    ! dec_drs = -2A * A1 * Q2 - Q0 * Q3 /(Q1*(one + Q1))
    d_denom_drs = half * (k04 + sq_rs * ( two * k05 + sq_rs * ( three * k06 + four * k07 * sq_rs)))/sq_rs
    ec_rs = k03 * log_fac_ec_lda - prefac_ec_lda * d_denom_drs / (denominator * (one + denominator))
    vc_unif = ec_unif - rs * third * ec_rs

    ! Non-local correlation potential
    t6 = t2*t4
    ! In Burke's routines,
    !   delt = beta/gamma which is beta_gamma here
    !   G is one if non-spin; F is zero if non-spin
    !   Q4 = 1 + A * t2
    !   GZ = zero if non-spin, as are FZ and ECZET
    A2 = A*A
    bg = -three * A2 * ec_unif * F2/beta
    bec = A2*F2/beta
    F5 = one + A * t2 !Q4
    F8 = F6 * F6 + beta_gamma * F5 * F6 * t2! Q8
    F9 = one + two * A * t2     ! Q9
    h_A = -beta * A * t6 * (two + A * t2) / F8
    h_rs = -third * rs * h_A * bec * ec_rs
    fact0 = two*beta_gamma - six * A ! FACT0
    fact1 = F6 * F9 + F5 * F9 * F9 !      FACT1 = Q5*Q9+Q4*Q9*Q9
    h_BT = two * beta * t4 * (F5 * F6 * fact0 - beta_gamma * fact1)/(F8 * F8)
    h_RST = third * rs * t2 * h_BT * bec * ec_rs ! RSTHRD*T2*hBT*BEC*ECRS
    h_T = two * beta * F9/F8 !2.d0*BET*G3*Q9/Q8
    fact2 = F5 * F6 + A * t2 * (F5 * F9 + F6) !FACT2 = Q4*Q5+B*T2*(Q4*Q9+Q5)
    fact3 = two * A * F6 * F9 + beta_gamma * fact2 !FACT3 = 2.D0*B*Q5*Q9+DELT*FACT2
    h_TT = four * beta * t * (two * A - F9*fact3/F8)/F8 !hTT = 4.D0*BET*G3*T*(2.D0*B/Q8-(Q9*FACT3/Q8)/Q8)
    ! Derivative wrt n (I'm pretty sure)
    dvc = H + h_rs + h_RST + t2 * h_t / six + seven * t2 * t * h_TT / six !COMM = H+HRS+HRST+T2*HT/6.D0+7.D0*T2*T*HTT/6.D0
    ! Derivative wrt grad n (again, I'm pretty sure)
    dvc = dvc - u * h_TT - v * h_T !COMM = COMM-PREF*ZET-UU*HTT-VV*HT-WW*(HZT-FACT5)

    ! These are Conquest lines which I will adapt and use to tidy the mess above later !
    ! Correlation
    !   (derivative of rho * e_correlation)
    ! NOTE: d_prefac_ec_lda is actually the derivative of rho*prefac_ec_lda
    !       d_log_fac_ec_lda is rho times the derivative of log_fac_ec_lda
    !d_prefac_ec_lda  = k02 + k08*rs
    !if (sq_rs > zero) then
    !   d_log_fac_ec_lda = &
    !        sq_rs * (k09 + sq_rs*(k10 + sq_rs*(k11 + k12 * sq_rs))) / &
    !        (denominator * (1 + denominator))
    !else
    !   d_log_fac_ec_lda = 0
    !end if
    !vc = d_prefac_ec_lda * log_fac_ec_lda + prefac_ec_lda * d_log_fac_ec_lda
    !if (present(drhoEps_c)) then
    !   drs_drho = - (third * rs / rho_tot_r)
    !   dkF_drho = third * kF / rho_tot_r
    !   dks_drho = half * ks * dkF_drho / kF
    !   deps_c_unif_drho = (Vc_unif(1) - eps_c_unif) / rho_tot_r
    !   dt_drho = (- t) * (dks_drho / ks + one / rho_tot_r)
    !   ! add RD_ERR to denominator to avoid div by 0
    !   dF1_drho = F1 * (deps_c_unif_drho / (eps_c_unif + RD_ERR))
    !   dF2_Drho = (- F2) * dF1_drho
    !   ! add RD_ERR to denominator to avoid div by 0
    !   dA_drho = A * dF2_drho / (one - F2 + RD_ERR)
    !   dF3_drho = (two * t + four * A * t**3) * dt_drho + dA_drho * t**4
    !   dF4_drho = F4 * (dF3_drho / F3 - &
    !        (dA_drho * F3 + A * DF3_drho) / (one + A * F3))
    !
    !   dH_drho = gamma * DF4_drho / (one + F4)
    !
    !   ! get drhoEps_c / drho(spin)
    !   drhoEps_c(0,1) = Vc_unif(1) + H + rho_tot_r * dH_drho
    !   do ii = 1, 3
    !      dt_dgrho = (t / mod_grho_tot_r) * grho_tot_r(ii) / mod_grho_tot_r
    !      dF3_dgrho = dt_dgrho * (two * t + four * A * t**3)
    !      dF4_dgrho = F4 * dF3_dgrho * (one / F3 - A / (one + A * F3))
    !      dH_dgrho = gamma * dF4_dgrho / (one + F4)
    !      drhoEps_c(ii,1) = rho_tot_r * dH_dgrho
    !   end do
    !end if ! present(drhoEps_c)

  end subroutine pbe_corr

end module radial_xc

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
What kind of information, or services should be added?  Seems like a topic that probably affects almost all aspects of life to me...
  Parenting
  Family Life
  Legal Services
  Counseling and Psychological Services

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!# MODULE: FUNCTION FMT
MODULE FUN_FMT

!!## PURPOSE
!! The formatting function <FMT> returns <VS> (a varying string)
!! for input/output of a single numeric entity, which could be used
!! in any FORMAT specificer (like that in a <READ> or <WRITE> statement).

!!## DETAILS
!! The format spec is of the form  <Pre><F><y><.z><Post> where
!!  @ <F> is a format string (f,e,a,i,g, etc.),
!!  @ <y> is the field specifier,
!!  @ <z> is the number of significant digits,
!!  @ <Post> is a post-append string, and
!!  @ <Pre> is a pre-append string.

!!## OWNER
! [waw]: william.wieselquist@gmail.com

!!## EXTERNAL KINDS
USE KND_IntrinsicTypes   ,ONLY: KIND_L1,KIND_L2,KIND_L4,& !!((01-A-KND_IntrinsicTypes.f90))
                  KIND_I1,KIND_I2,KIND_I4,KIND_I8,&
                  KIND_Rsp,KIND_Rdp,&
                  KIND_Csp,KIND_Cdp,&
                  KIND_S,KIND_Sfmt

!!## EXTERNAL PARAMETERS
USE PAR_IntrinsicLengths,ONLY: LEN_L1,LEN_L2,LEN_L4,&     !!((02-A-PAR_IntrinsicLengths.f90))
                  LEN_I1,LEN_I2,LEN_I4,LEN_I8,&
                  LEN_Rsp,LEN_Rdp,&
                  LEN_Csp,LEN_Cdp

!!#### FORTRAN STANDARDS MODULES
USE ISO_varying_string                                    !!((03-A-ISO_varying_string.f90))

!!## GLOBAL BINARY OPERATORS
USE BOP_sEQ                                               !!((03-A-BOP_sEQ.f90))

!!## EXTERNAL PROCEDURES
USE FUN_STR                                               !!((05-B-FUN_STR.f90))

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PRIVATE

!!## PROCEDURE OVERLOADING
INTERFACE FMT
 MODULE PROCEDURE FMT_Sen
 MODULE PROCEDURE FMT_Rsp
 MODULE PROCEDURE FMT_Rdp
 MODULE PROCEDURE FMT_Csp
 MODULE PROCEDURE FMT_Cdp
 MODULE PROCEDURE FMT_I1
 MODULE PROCEDURE FMT_I2
 MODULE PROCEDURE FMT_I4
 MODULE PROCEDURE FMT_I8
 MODULE PROCEDURE FMT_L1
 MODULE PROCEDURE FMT_L2
 MODULE PROCEDURE FMT_L4
END INTERFACE

INTERFACE FMT_
 MODULE PROCEDURE FMT_null
END INTERFACE

!!## PUBLIC ACCESS LIST
PUBLIC :: FMT,FMT_
PUBLIC :: FMT_Sen
PUBLIC :: FMT_Rsp
PUBLIC :: FMT_Rdp
PUBLIC :: FMT_Csp
PUBLIC :: FMT_Cdp
PUBLIC :: FMT_I1
PUBLIC :: FMT_I2
PUBLIC :: FMT_I4
PUBLIC :: FMT_I8
PUBLIC :: FMT_L1
PUBLIC :: FMT_L2
PUBLIC :: FMT_L4

!!## MODULE PROCEDURES
CONTAINS

!!### FUNCTION: FMT_NULL
FUNCTION FMT_null( F , y , z , Post , Pre ) RESULT(VS)
INCLUDE "06-B-FUN_FMT_null.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_null.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_Sen
FUNCTION FMT_Sen( X , y , Post , Pre ) RESULT(VS)
INCLUDE "06-B-FUN_FMT_S.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_S.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_Rsp
FUNCTION FMT_Rsp( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rsp             !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_R=LEN_Rsp

INCLUDE "06-B-FUN_FMT_R.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_R.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_Rdp
FUNCTION FMT_Rdp( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_R=>KIND_Rdp             !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_R=LEN_Rdp

INCLUDE "06-B-FUN_FMT_R.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_R.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_Csp
FUNCTION FMT_Csp( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Csp             !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_C=LEN_Csp

INCLUDE "06-B-FUN_FMT_C.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_C.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_Cdp
FUNCTION FMT_Cdp( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_C=>KIND_Cdp             !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_C=LEN_Cdp

INCLUDE "06-B-FUN_FMT_C.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_C.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_I1
FUNCTION FMT_I1( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I1              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_I=LEN_I1

INCLUDE "06-B-FUN_FMT_I.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_I.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_I2
FUNCTION FMT_I2( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I2              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_I=LEN_I2

INCLUDE "06-B-FUN_FMT_I.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_I.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_I4
FUNCTION FMT_I4( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I4              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_I=LEN_I4

INCLUDE "06-B-FUN_FMT_I.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_I.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_I8
FUNCTION FMT_I8( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_I=>KIND_I8              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_I=LEN_I8

INCLUDE "06-B-FUN_FMT_I.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_I.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_L1
FUNCTION FMT_L1( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L1              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_L=LEN_L1
INCLUDE "06-B-FUN_FMT_L.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_L.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_L2
FUNCTION FMT_L2( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L2              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_L=LEN_L2
INCLUDE "06-B-FUN_FMT_L.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_L.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

!!### FUNCTION: FMT_L4
FUNCTION FMT_L4( X , F , y , z , Post , Pre ) RESULT(VS)
!!#### LOCAL KINDS/LENGTHS
USE KND_IntrinsicTypes,ONLY: KIND_L=>KIND_L4              !!((01-A-KND_IntrinsicTypes.f90))
INTEGER,PARAMETER :: LEN_L=LEN_L4

INCLUDE "06-B-FUN_FMT_L.f90.hdr"
INCLUDE "06-B-FUN_FMT.f90.hdr"
!!--begin--
INCLUDE "06-B-FUN_FMT_L.f90.bdy"
INCLUDE "06-B-FUN_FMT.f90.bdy"
!!--end--
END FUNCTION

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


! { dg-do run }
!
! Test the fix for PR61147.
!
! Contributed by Thomas Clune  <Thomas.L.Clune@nasa.gov>
!
module B_mod

   type :: B
      character(:), allocatable :: string
   end type B

contains

   function toPointer(this) result(ptr)
      character(:), pointer :: ptr
      class (B), intent(in), target :: this

         ptr => this%string

   end function toPointer

end module B_mod

program main
   use B_mod

   type (B) :: obj
   character(:), pointer :: p

   obj%string = 'foo'
   p => toPointer(obj)

   If (len (p) .ne. 3) call abort
   If (p .ne. "foo") call abort

end program main


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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 registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!            DELWAQ - Deltares WAter Quality program
!
!                     Version 4.xx - june 2009 - DLL version
!                     Version 4.33 - july 1998
!                     Version 4.32 - june 1998
!
!     INFORMATION   : Deltares
!                     L. Postma,
!                     Rotterdamseweg 185,
!                     P.O. Box 177,
!                     2600 MH Delft,
!                     Netherlands.
!                     telephone (31) 88-3358273
!                     telefax   (31) 88-3358582
!
!     FUNCTION            : MAIN module for DELWAQ2 , dimensioning
!                           of the work array's.
!
!     SUBROUTINES CALLED  : DELWQ2, performs the simulation
!
!     PARAMETERS          :
!
!     NAME    KIND     LENGTH     FUNCT.  DESCRIPTION
!     ----    -----    ------     ------- -----------
!     ACTION  INTEGER  1          INPUT   Action to be taken
!     ARGC    INTEGER  1          INPUT   Number of simulated command-line arguments
!     ARGV    INTEGER  1          INPUT   Simulated command-line arguments
!
!     ITOTA   INTEGER  1          INPUT   length of real workarray
!     ITOTI   INTEGER  1          INPUT   length of integer workarray
!     ITOTC   INTEGER  1          INPUT   length of character workarray
!
!
!      PARAMETER (ITOTA=0       ,ITOTI=0       ,ITOTC=0       )

      subroutine delwaq2( argc, argv, errorcode )
      !DEC$ ATTRIBUTES DLLEXPORT::delwaq2

      use delwaq2_data
      implicit none
      include 'actions.inc'

      integer, intent(in)                           :: argc
      character(len=*), dimension(argc), intent(in) :: argv
      integer, intent(out)                          :: errorcode

      type(delwaq_data)                             :: dlwqd

      dlwqd%set_timer = .true.

      call dlwqmain( ACTION_FULLCOMPUTATION, argc, argv, dlwqd )

      ! Delwaq2_lib should never use a stop, but must be modified to return an error code instead (0 = normal end)
      ! Currently a return from the delwaq2_lib assumes a normal end.
      errorcode = 0
      return

      end subroutine delwaq2

      subroutine dlwqmain(ACTION, ARGC, ARGV, DLWQD)

      !DEC$ ATTRIBUTES DLLEXPORT::dlwqmain

      USE DELWAQ2
      USE DELWAQ2_DATA
      USE DHCOMMAND

      IMPLICIT NONE

      INTEGER, INTENT(IN)                           :: ACTION
      INTEGER, INTENT(IN)                           :: ARGC
      CHARACTER(LEN=*), DIMENSION(ARGC), INTENT(IN) :: ARGV
      TYPE(DELWAQ_DATA)                             :: DLWQD

      CHARACTER*20  RUNDAT
!
      LOGICAL                                       :: INIT        ! Do not save!
      INTEGER                                       :: LUNREP

      INTEGER, SAVE                                 :: ITOTA
      INTEGER, SAVE                                 :: ITOTI
      INTEGER, SAVE                                 :: ITOTC

      INCLUDE 'sysn.inc'
      INCLUDE 'sysi.inc'
      INCLUDE 'actions.inc'
      INCLUDE 'fsm-fix.i'

!
!     Initial step ...
!
      INIT = .FALSE.
      IF ( ACTION == ACTION_INITIALISATION  .OR.
     &     ACTION == ACTION_FULLCOMPUTATION      ) THEN
          INIT = .TRUE.

          call dhstore_command( argv )

          CALL AVUNDF
!
          ITOTA=0
          ITOTI=0
          ITOTC=0

          NULLIFY( DLWQD%RBUF  )
          NULLIFY( DLWQD%IBUF  )
          NULLIFY( DLWQD%CHBUF )

          ALLOCATE( DLWQD%RBUF(0) )
          ALLOCATE( DLWQD%IBUF(0) )
          ALLOCATE( DLWQD%CHBUF(0) )

      ENDIF

!     Computation step is always done

      CALL DELWQ2 ( DLWQD%RBUF, DLWQD%IBUF, DLWQD%CHBUF, ITOTA, ITOTI,
     &              ITOTC, INIT, ACTION, DLWQD )

!
!     Finalise - only if the full computation was done
!
      IF ( ACTION == ACTION_FULLCOMPUTATION      ) THEN

          CALL GETMLU(LUNREP)
          WRITE ( * , * )
          WRITE ( * , * ) ' SIMULATION ENDED '
          WRITE ( LUNREP , * )
          WRITE ( LUNREP , '(A)' ) ' Simulation ended normal'
          CALL DATTIM(RUNDAT)
          WRITE (LUNREP,'(2A)') ' Execution stop : ',RUNDAT
!
          close(lunrep)

          RETURN

      ENDIF

      IF ( ACTION == ACTION_FINALISATION      ) THEN

          CALL GETMLU(LUNREP)
          WRITE ( * , * )
          WRITE ( * , * ) ' SIMULATION ENDED '
          WRITE ( * , * )
          WRITE ( LUNREP , * )
          WRITE ( LUNREP , '(A)' ) ' Simulation ended normal'
          CALL DATTIM(RUNDAT)
          WRITE (LUNREP,'(2A)') ' Execution stop : ',RUNDAT
!
          close(lunrep)

      ENDIF

      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C$Procedure      READLN ( Read a text line from a logical unit )

      SUBROUTINE READLN ( UNIT, LINE, EOF )

C$ Abstract
C
C     This routine will read a single line of text from the Fortran
C     logical unit UNIT, reporting the end of file if it occurs.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     ASCII
C     TEXT
C     FILES
C
C$ Declarations
C
      INTEGER            UNIT
      CHARACTER*(*)      LINE
      LOGICAL            EOF
C
C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C      UNIT      I    The Fortran unit number to use for input.
C      LINE      O    The line read from the file.
C      EOF       O    A logical flag indicating the end of file.
C
C$ Detailed_Input
C
C     UNIT     The Fortran unit number for the input. This may
C              be either the unit number for the terminal, or the
C              unit number of a previously opened text file.
C
C$ Detailed_Output
C
C     LINE     On output, this will contain the next text line
C              encountered when reading from UNIT.
C
C              If the length of the character string LINE is shorter
C              than the length of the current line in the text file, the
C              line is truncated on the right by the Fortran READ
C              statement, filling LINE with the first LEN(LINE)
C              characters from the current line in the file.
C
C              If an error or the end of file occurs during the
C              attempt to read from UNIT, the value of this variable
C              is not guaranteed.
C
C     EOF      On output, this variable will be set to .TRUE. if the
C              end of file ( IOSTAT < 0 ) is encountered during the
C              attempt to read from unit UNIT. Otherwise, this
C              variable will be set to .FALSE..
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)   If an error occurs while attempting to read from the text
C          file attached to UNIT, the error SPICE(FILEREADFAILED) will
C          be signalled.
C
C     This routine only checks in with the error handler in the event
C     that an error occurred. (Discovery check in)
C
C$ Files
C
C     None.
C
C$ Particulars
C
C      This routine will read a single line, a text record, from the
C      logical unit UNIT. UNIT may be the terminal, or it may be a
C      logical unit number obtained from a Fortran OPEN or INQUIRE
C      statement. This routine will set a logical flag, EOF, on output
C      if the end of the file is encountered during the read attempt.
C
C$ Examples
C
C      CALL READLN ( UNIT, LINE, EOF )
C
C      IF ( EOF ) THEN
C         < The end of file, deal with it appropriately >
C      END IF
C
C      You now have a line of text from unit UNIT.
C
C$ Restrictions
C
C      None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer (JPL)
C
C$ Version
C
C-    SPICELIB     1.0.0, 20-DEC-1995 (KRG)
C
C        The routine graduated
C
C-    Beta Version 1.0.1, 22-NOV-1994 (KRG)
C
C        Cleaned up the comments a little bit. No code changes.
C
C-    Beta Version 1.0.0, 17-DEC-1992 (KRG)
C
C-&

C$ Index_Entries
C
C      read a text line from a logical unit
C
C-&

C
C     Local variables
C
      INTEGER               IOSTAT
C
C     Standard SPICE error handling.
C
C
C     Read in the next line from the text file attached to UNIT.
C
      READ ( UNIT, FMT='(A)', IOSTAT=IOSTAT ) LINE
C
C     Check to see if we got a read error, and signal it if we did.
C
      IF ( IOSTAT .GT. 0  ) THEN

         CALL CHKIN  ( 'READLN'                                  )
         CALL SETMSG ( 'Error reading from file: #. IOSTAT = #.' )
         CALL ERRFNM ( '#', UNIT                                 )
         CALL ERRINT ( '#', IOSTAT                               )
         CALL SIGERR ( 'SPICE(FILEREADFAILED)'                   )
         CALL CHKOUT ( 'READLN'                                  )
         RETURN

      END IF
C
C     Check to see if we got the end of file, and set the logical
C     flag EOF if we did.
C
      IF ( IOSTAT .LT. 0 ) THEN

         EOF = .TRUE.

      ELSE

         EOF = .FALSE.

      END IF

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      program archt2archv
      use mod_plot  ! HYCOM plot array interface
      use mod_za    ! HYCOM I/O interface
      implicit none
c
      integer, parameter   :: in_max =  9999    !maximum numder of input files
      integer, parameter   :: ib_off = 10000    !I/O unit offset for .B files
      integer, parameter   :: ia_off = 20000    !I/O unit offset for .A files
      real*4,  parameter   :: spval  = 2.0**100 !data void marker
c
c     convert a partial (tiled) archive to a full region archive file.
c
      character*80         :: cline
      character*256        :: flnm_in,flnm_out
      integer              :: i,id,ios,irec,j,k,l,n,ni,no
      real                 :: hmina,hmaxa
      integer, allocatable :: m_out(:,:)
      real,    allocatable :: a_out(:,:)
      real*4,  allocatable :: a_in(:)
c
      integer, dimension(in_max)
     &                     :: in_i0,in_ii,in_j0,in_jj
      integer              :: in_n,in_size
      logical              :: lend
c
      call xcspmd
      call zaiost
      lp=6
c
      allocate( m_out(idm,jdm) )
      allocate( a_out(idm,jdm) )
      a_out(:,:) = spval
c
c --- single output file.
c ---   'flnm_out  ' = name of output archive
      read (*,'(a)') flnm_out
      write (lp,'(2a)') 'output file: ',trim(flnm_out)
      call zhflsh(lp)
      no = 15
      l  = len_trim(flnm_out)
      open (unit=  no,
     &      file=  flnm_out(1:l-2)//'.b',
     &      form=  'formatted',
     &      status='new',
     &      action='write')
      call zaiopf(flnm_out(1:l-2)//'.a','new', no)
c
c --- many input files
c
      in_size = (idm*jdm)/1000 + idm+jdm  !guess the input array size
      write(lp,*) '1st in_size = ',in_size
      call zhflsh(lp)
      allocate( a_in(in_size) )
c
      n = 0
      do
c ---   'flnm_in   ' = name of input archive
        read (*,'(a)',iostat=ios) flnm_in
        if     (ios.ne.0) then
          exit !do
        endif
        n  = n + 1
        if     (n.gt.in_max) then
          write(lp,*)
          write(lp,*) 'error more than',in_max,' input files'
          write(lp,*)
          stop
        endif
        ni = n+ib_off
        l  = len_trim(flnm_in)
        open (unit=ni,
     &        file=flnm_in(1:l-2)//'.B',form='formatted',
     &        status='old',action='read')
c
c       process the archive header
c
        if     (n.eq.1) then
          do k= 1,7
            read( ni,'(a)') cline
            write(no,'(a)') cline
            write(lp,'(a)') cline
          enddo
c
          read( ni,*) in_i0(n)  !i1
          in_i0(n) = in_i0(n) - 1
          read( ni,*) in_j0(n)  !j1
          in_j0(n) = in_j0(n) - 1
          read( ni,*) in_ii(n)
          read( ni,*) in_jj(n)
c
          write(no,'(i5,4x,a)') idm,"'idm   ' = longitudinal array size"
          write(lp,'(i5,4x,a)') idm,"'idm   ' = longitudinal array size"
          write(no,'(i5,4x,a)') jdm,"'jdm   ' = latitudinal  array size"
          write(lp,'(i5,4x,a)') jdm,"'jdm   ' = latitudinal  array size"
c
          read( ni,'(a)') cline  ! field ...
          write(no,'(a)') cline
          write(lp,'(a)') cline
c
          call zhflsh(no)
          call zhflsh(lp)
        else
          do k= 1,7
            read( ni,'(a)') cline
          enddo
          read( ni,*) in_i0(n)  !i1
          in_i0(n) = in_i0(n) - 1
          read( ni,*) in_j0(n)  !j1
          in_j0(n) = in_j0(n) - 1
          read( ni,*) in_ii(n)
          read( ni,*) in_jj(n)
          read( ni,'(a)') cline  ! field ...
        endif
        if     (in_ii(n)*in_jj(n).gt.in_size) then
          in_size = in_ii(n)*in_jj(n) + idm+jdm
          write(lp,*) 'new in_size = ',in_size
          call zhflsh(lp)
          deallocate( a_in )
            allocate( a_in(in_size) )
        endif
        inquire( iolength=irec ) a_in(1:in_ii(n)*in_jj(n))
        open(unit=  n+ia_off,
     &       file=  flnm_in(1:l-2)//'.A',
     &       form=  'unformatted',
     &       status='old',
     +       access='direct',
     &       recl=  irec,
     &       iostat=ios)
        if     (ios.ne.0) then
          write(lp,*)
          write(lp,*) 'Error: can''t open ',flnm_in(1:l-2)//'.A'
          write(lp,*) 'ios   = ',ios
          write(lp,*) 'unit  = ',n+ia_off
          write(lp,*) 'nrecl = ',irec
          write(lp,*)
          stop
        endif
      enddo
      in_n = n
c
c --- any input?
c
      if     (in_n.eq.0) then
          write(lp,*)
          write(lp,*) 'Error: no input files'
          write(lp,*)
          stop
      endif
c
c     loop through all 2-d fields in archive file.
c
      lend = .false.
      do irec= 1,huge(irec)
        do n= 1,in_n
          read( n+ib_off,'(a)',iostat=ios) cline
          if     (ios.ne.0) then
            if     (n.eq.1) then
              lend = .true.
              exit !do n
            else
              write(lp,*)
              write(lp,*) 'error reading .B - n,irec = ',n,irec
              write(lp,*)
              stop
            endif
          endif
          call zaiordd(a_in,in_ii(n)*in_jj(n),n+ia_off,irec,ios)
          if     (ios.ne.0) then
            write(lp,*)
            write(lp,*) 'error reading .A - n,irec = ',n,irec
            write(lp,*) trim(cline)
            write(lp,*)
            stop
          endif
          do j= 1,in_jj(n)
            do i= 1,in_ii(n)
              id = mod(in_i0(n)+i-1+2*idm,idm)+1  ! allow periodic wrap
              a_out(id,in_j0(n)+j) = a_in(i+(j-1)*in_ii(n))
            enddo !i
          enddo !j
        enddo !n
c
        if     (lend) then
          exit !do irec
        endif
c
        call zaiowr(a_out,m_out,.false., hmina,hmaxa, no, .false.)
        write(no,'(a,1p2e16.7)') cline(1:41),hmina,hmaxa
        call zhflsh(no)
        write(lp,'(a,1p2e16.7)') cline(1:41),hmina,hmaxa
        call zhflsh(lp)
      enddo !irec
c
      end program archt2archv

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


module adjustl_test
    use iso_varying_string, only: adjustl, char
    use veggies, only: &
            input_t, &
            result_t, &
            string_input_t, &
            test_item_t, &
            assert_equals, &
            describe, &
            fail, &
            it, &
            ASCII_STRING_GENERATOR

    implicit none
    private
    public :: test_adjustl
contains
    function test_adjustl() result(tests)
        type(test_item_t) :: tests

        tests = describe( &
                "Sec 3.4.1: ADJUSTL", &
                [ it( &
                        "works the same for characters and strings", &
                        ASCII_STRING_GENERATOR, &
                        check_adjustl) &
                ])
    end function

    pure function check_adjustl(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        select type (input)
        type is (string_input_t)
            result_ = assert_equals( &
                    adjustl(char(input%input())), &
                    adjustl(input%input()))
        class default
            result_ = fail("Expected to get a string_input_t")
        end select
    end function
end module


text:
       module calc_maxcomp

        use polsynchemis, only: polsynchth

       contains

!         subroutine calc_maxjutt_subroutine(nu,ncgs,bcgs,tcgs,usim,incang,maxjutt_args,ktotal)
         subroutine calc_maxcomp_subroutine(nu,ncgs,bcgs,tcgs,incang,maxjutt_args,ktotal)
           double precision, dimension(:), intent(in) :: nu
           double precision, dimension(:), intent(in) :: ncgs,bcgs,tcgs,incang
           double precision, dimension(:), intent(in) :: maxjutt_args
!6 weights, 1 delta, 1 selection option
!           double precision, dimension(:), intent(in) :: usim
!2,3,5,7,11,13,17
           double precision, dimension(size(maxjutt_args) - 2) :: weights_arr, i_arr, delta_arr
           double precision, dimension(size(ncgs)) :: tcgs_min
!           double precision, dimension(size(ncgs)) :: theta_min
           double precision :: delta, a, dwsum

           double precision, dimension(size(ncgs),11) :: ktemp
           double precision, dimension(size(ncgs),11),intent(out) :: ktotal

           integer :: isel

!           ones = 1d0
!           write(6,*) 'maxjutt: ',size(maxjutt_args),size(weights_arr)
           a = 3d0
!Note derivation assumes relativistic a = 3.0
           delta = maxjutt_args(1) !first argument is delta
           selection = maxjutt_args(2) !second argument selects the decomposition
           weights_arr = maxjutt_args(3:)/SUM(maxjutt_args(3:)) !normalize
!           write(6,*) 'maxjutt weights: ',weights_arr,delta
           i_arr = (/ (I, I=0,size(weights_arr)-1) /) !indices
           delta_arr = delta**i_arr !all deltas to be multiplied by thetamin
           dwsum = SUM(weights_arr * delta_arr)
           tcgs_min = tcgs/dwsum

!pseudocode
!ktotal = whatever
!convert units to cgs?
!for i in (1,size(weights_arr)):
!    ktemp  = polsynchth(nu,weights_arr(i)*ncgs,bcgs,thetamin*delta_arr(i),e%incang,Kth)
!    ktotal = ktotal + ktemp
!total absorption
!partial absorption from selection
!selection = 6 returns the top maxwellian

!return K to emis.f90
           ktotal = 0d0
           do i=1,size(weights_arr)
              ktemp = 0d0
!              write(6,*) 'calc_maxjutt call polsynchth: ',weights_arr(i),ncgs(1),bcgs(1),tcgs_min(1)*delta_arr(1),incang(1)
              call polsynchth(nu,weights_arr(i)*ncgs,bcgs,tcgs_min*delta_arr(i),incang,ktemp)
              ktotal = ktotal + ktemp
           enddo

           isel = INT(selection)
           if(isel.gt.0) then
              if(isel.le.(size(weights_arr))) then
                 ktemp = 0d0
                 call polsynchth(nu,weights_arr(isel)*ncgs,bcgs,tcgs_min*delta_arr(isel),incang,ktemp)
                 ktotal(:,1:4) = ktemp(:,1:4)
              endif
           endif

!           write(6,*) 'maxjutt loop: ',tcgs_min,delta_arr,i,size(weights_arr)
         end subroutine calc_maxcomp_subroutine

         end module calc_maxcomp

max_stars_repo_name: HerculesJack/grtrans

max_issues_repo_name: HerculesJack/grtrans

max_forks_repo_name: HerculesJack/grtrans