text:
!! Fixed in 6.1 build 6144 (report 99312)
!!
!! Bad intermediate C code is being generated for this example
!! resulting in errors from the C compiler.
!!
!! $ nagfor -c nag-bug-20171110.f90
!! NAG Fortran Compiler Release 6.1(Tozai) Build 6143
!! Warning: nag-bug-20171110.f90, line 42: Unused dummy variable THIS
!! Warning: nag-bug-20171110.f90, line 46: Unused dummy variable THIS
!! Warning: nag-bug-20171110.f90, line 46: INTENT(OUT) dummy argument THIS never set
!! [NAG Fortran Compiler normal termination, 3 warnings]
!! nag-bug-20171110.f90: In function 'fubar_MP_init':
!! nag-bug-20171110.f90:44:14: error: 'this_' is a pointer; did you mean to use '->'?
!!    subroutine init(this)
!!               ^
!!               ->
!! nag-bug-20171110.f90: In function 'fubar_MP__zfinal0':
!! nag-bug-20171110.f90:21:11: error: 'v_' is a pointer; did you mean to use '->'?
!!  module fubar
!!            ^
!!            ->

module fubar

  type x
  contains
    final :: x_final
  end type

  type y
    type(x) :: xval
  end type

  type z
    type(y) :: yval(2)
  contains
    procedure :: init
  end type

contains

  elemental subroutine x_final(this)
    type(x), intent(inout) :: this
  end subroutine

  subroutine init(this)
    class(z), intent(out) :: this
  end subroutine

end module

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:
http://wheel.ucdavis.edu/ The Wheel is an educational technology blog published by Dr. Andy Jones and Mary Stewart (PhD student, Education). The blog is part of Dr. Andy’s initiative to reach out to faculty through his role as academic associate director at Academic Technology Services.

Topics include:
  http://wheel.ucdavis.edu/category/edtechinnovations/ Educational Technology Innovations
  http://wheel.ucdavis.edu/category/dolce/ Discussing Online Learning and Collaborative Education
  http://wheel.ucdavis.edu/category/socialmedia/ Social Media
  http://wheel.ucdavis.edu/category/inpersonbriefings/ SmartSite Support
  http://wheel.ucdavis.edu/category/audiovideobriefings/ Audio & Video
  http://wheel.ucdavis.edu/category/profiles/ Meet the Staff
  http://wheel.ucdavis.edu/category/facultyspotlights/ Faculty Spotlights

Our goal is to provide information that will help you integrate technology into your classroom without reinventing the wheel. This includes posting information about innovations in the field of educational technology, as well as updates about the monthly Discussing Online Learning and Collaborative Education meetings (first Friday of the month, 12pm, 1310 Surge III).

We also host monthly Faculty Panels (third Friday of the month, 12pm, 1310 Surge III), which are informal faculty discussions about teaching with specific technologies, e.g., Blogging or YouTube.

Visit http://wheel.ucdavis.edu/ The Wheel for more details about our upcoming events!

You can also follow us on Twitter https://twitter.com/ucdaviswheel @ucdaviswheel.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      integer qq,qg,gq,gg
      parameter (qq=1,qg=2,gq=3,gg=4)
      logical qqproc,qgproc,gqproc,ggproc
      common/dipproc/qqproc,qgproc,gqproc,ggproc
!$omp threadprivate(/dipproc/)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


SUBROUTINE CloseVariables
! By Allan P. Engsig-Karup.
USE GlobalVariables
USE MGlevels
IMPLICIT NONE
DEALLOCATE(Wavefield%E,Wavefield%W,Wavefield%P,FineGrid%h,tmp2D)
IF (FineGrid%Nx>1) THEN
	DEALLOCATE(Wavefield%Ex,Wavefield%Exx,Wavefield%Px)
	DEALLOCATE(FineGrid%hx,FineGrid%hxx)
!	DEALLOCATE(FineGrid%DiffStencils%StencilX)
!	DEALLOCATE(FullRankStencils%StencilX)
END IF
IF (FineGrid%Ny>1) THEN
	DEALLOCATE(Wavefield%Ey,Wavefield%Eyy,Wavefield%Py)
	DEALLOCATE(FineGrid%hy,FineGrid%hyy)
!	DEALLOCATE(FineGrid%DiffStencils%StencilY)
!	DEALLOCATE(FullRankStencils%StencilY)
END IF
!DEALLOCATE(FullRankStencils%StencilZ)
!DEALLOCATE(FullRankStencils%Indexes)
!DEALLOCATE(FullRankStencils%Indexesnew)
!DEALLOCATE(FullRankStencils%StencilXZorYZ)
DEALLOCATE(rhsE,rhsP,PHI,RHS)
!DEALLOCATE(FineGrid%dsigma)
IF (Precond==3) THEN
	DEALLOCATE(arrLevels) ! From Multigrid part
ENDIF
DEALLOCATE(COLSCA, ROWSCA,workspace)
! Deallocate grid information (FIXME)
IF (relaxONOFF==1) THEN
	DEALLOCATE(RelaxZones)
ENDIF
! GD: New variable definition, SWENSE part
IF(swenseONOFF==1) THEN !We deallocate the incident wavefields values
    DEALLOCATE(Wavefield%E_I,Wavefield%Et_I,Wavefield%P_I_s,Wavefield%Pz_I_s,Wavefield%Pt_I_s,Wavefield%E_I_bp,Wavefield%Pz_I_bp)
    IF (FineGrid%Nx>1) THEN
    	DEALLOCATE(Wavefield%Ex_I, Wavefield%Exx_I, Wavefield%Px_I_s, Wavefield%Ex_I_bp,Wavefield%Px_I_bp)
    ENDIF
    IF (FineGrid%Ny>1) THEN
    	DEALLOCATE(Wavefield%Ey_I, Wavefield%Eyy_I, Wavefield%Py_I_s, Wavefield%Ey_I_bp, Wavefield%Py_I_bp)
	ENDIF
! Deallocate the surfacic source terms
! FIXME: Define even if non-SWENSE ?
! APEK: moved lines below inside loop as these are only allocated in SWENSE?
IF (FineGrid%Nx>1) THEN
    DEALLOCATE(Wavefield%SourceEx, Wavefield%SourcePx)
ENDIF
IF (FineGrid%Ny>1) THEN
    DEALLOCATE(Wavefield%SourceEy, Wavefield%SourcePy)
ENDIF
ENDIF
END SUBROUTINE CloseVariables


text:
subroutine intpol(kmax      ,kmxdt     ,kmxt      ,h0        ,zeta      , &
                & dps       ,siglim    ,barcli    ,riav      ,rilz      , &
                & xk        ,u0        ,dudz      ,d2u0      ,ustbx     , &
                & ustwix    ,v0        ,dvdz      ,d2v0      ,ustby     , &
                & ustwiy    ,sig       ,thick     ,rich      ,bruvai    , &
                & rtur0     ,vicww     ,scale     ,bvav      ,bvmx      , &
                & xkh       ,kbed      ,ktop      ,umean     ,vmean     , &
                & utg       ,vtg       ,d2u       ,d2v       ,bv2       , &
                & tke       ,eps       ,edvis     ,zlw       ,tlw       , &
                & zlo       ,useref    ,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: intpol.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/taylor/intpol.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Interpolation on vertical profiles of:
!             -horizontal velocity,
!             -bouyancy frequency and
!             -turbulence parameters.
!              Estimate vertical curvature of horizontal
!              velocity profile.
! Method used: linear interpolation and matching to
!              logarithmic layers.
!
!!--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
    !
    real(fp)               , pointer :: ricrit
    real(fp)               , pointer :: bvmin
    real(fp)               , pointer :: epsuko
    real(fp)               , pointer :: epsbv2
    real(fp)               , pointer :: clw
    real(fp)               , pointer :: ckw
    real(fp)               , pointer :: viscof
    real(fp)               , pointer :: ag
    real(fp)               , pointer :: z0
    real(fp)               , pointer :: vonkar
!
! Global variables
!
    integer                                     :: kbed   !!  KTG-value of upper level stratified
                                                          !!  layer
    integer                       , intent(in)  :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer                       , intent(in)  :: kmxdt  !  Description and declaration in dimens.igs
    integer                       , intent(in)  :: kmxt   !  Description and declaration in dimens.igs
    integer                                     :: ktop   !!  KTG-value of lower level stratified
                                                          !!  layer
    logical                       , intent(out) :: barcli !!  Switch for applying IWE model
    logical                       , intent(in)  :: useref
    real(fp)                                    :: bvav   !!  Mean bouyancy frequency
    real(fp)                      , intent(out) :: bvmx   !!  Maximal buoyancy frequency
    real(prec)                    , intent(in)  :: dps    !  Description and declaration in esm_alloc_real.f90
    real(fp)                                    :: h0     !!  Water depth
    real(fp)                                    :: riav   !!  Mean value of all RICH with RICH >
                                                          !!  RICRIT
    real(fp)                                    :: rilz   !!  Vertical interval with RICH > RICRIT
    real(fp)                                    :: scale  !!  Length scale for non-dimensionalising
    real(fp)                      , intent(in)  :: siglim !  Description and declaration in iwepar.igs
    real(fp)                                    :: umean  !  Description and declaration in esm_alloc_real.f90
    real(fp)                                    :: ustbx
    real(fp)                                    :: ustby
    real(fp)                      , intent(in)  :: ustwix !!  X-component of wind-shear stress
                                                          !!  velocity
    real(fp)                      , intent(in)  :: ustwiy !!  Y-component of wind-shear stress
                                                          !!  velocity
    real(fp)                                    :: vmean  !  Description and declaration in esm_alloc_real.f90
    real(fp)                      , intent(in)  :: xk     !!  Magnitude of horizontal wave number
                                                          !!  vector, in [m^-1]
    real(fp)                      , intent(out) :: xkh    !!  XK normalized by SCALE
    real(fp)                      , intent(in)  :: zeta   !!  Watrface elevation
    real(fp)                                    :: zlo    !!  Ozmidov length scale, depth-averaged
                                                          !!  over stratified layer
    real(fp), dimension(0:kmax)   , intent(in)  :: bruvai !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)                 :: d2u0   !!  Dummay array, curvature of vertical
                                                          !!  profile of u-velocity component
                                                          !!  component on hydrodynamic grid
    real(fp), dimension(0:kmax)                 :: d2v0   !!  Dummy array, curvature of vertical
                                                          !!  profile of v-velocity component
                                                          !!  component on hydrodynamic grid
    real(fp), dimension(0:kmax)   , intent(in)  :: dudz   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)   , intent(in)  :: dvdz   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)   , intent(in)  :: rich   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)   , intent(in)  :: vicww  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax, 2), intent(in)  :: rtur0  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmxdt)                :: bv2    !!  Interpolated bouyancy frequency,
                                                          !!   squared
    real(fp), dimension(0:kmxdt)                :: d2u    !!  D2U0, interpolated
    real(fp), dimension(0:kmxdt)                :: d2v    !!  D2V0, interpolated
    real(fp), dimension(0:kmxdt)                :: edvis  !!  Non-dimensional eddy viscosity, on
                                                          !!  interpolated grid
    real(fp), dimension(0:kmxdt)                :: eps    !  Description and declaration in numeco.igs
    real(fp), dimension(0:kmxdt)                :: tke    !!  Non-dimensional interpolated TKE=
                                                          !!  RTUR0(K,1)
    real(fp), dimension(0:kmxdt)  , intent(out) :: tlw    !!  Non-dimensional integral time scale
                                                          !!  of vertical turbulent motions
    real(fp), dimension(0:kmxdt)                :: utg    !!  Non-dimensional interpolated u-vel.
                                                          !!  component, relative to vel. at KREF
    real(fp), dimension(0:kmxdt)                :: vtg    !!  As UTG for v-velocity component
    real(fp), dimension(0:kmxdt)                :: zlw    !!  Non-dimensional integral length scale
                                                          !!  of turbulence in z/w-direction
    real(fp), dimension(kmax)     , intent(in)  :: sig    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)     , intent(in)  :: thick  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)     , intent(in)  :: u0     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)     , intent(in)  :: v0     !  Description and declaration in esm_alloc_real.f90
!
! Local variables
!
    integer  :: k
    integer  :: kd
    integer  :: kdpm   ! Sigma-level pointer
    integer  :: kdpmd
    integer  :: kref   ! Reference level for velocity, on interpolated grid
    integer  :: ksum
    integer  :: ktg    ! Pointer for interpolation arrays
    real(fp) :: bv2av  ! Mean buoyancy frequency squared
    real(fp) :: bv2lvl ! Definition of minimal buoyancy frequency, squared, in stratified layer
    real(fp) :: bv2mx  ! Maximal buoyancy frequency squared
    real(fp) :: diff
    real(fp) :: dsig
    real(fp) :: dz     ! Thickness of sigma layer
    real(fp) :: dztg   ! Layer thickness of interpolation grid
    real(fp) :: epsktg
    real(fp) :: epsscl ! Scaling factor for dissipation [m^2/s^3]
    real(fp) :: ratio
    real(fp) :: rmsw   ! Rms of turbulent vertical velocity
    real(fp) :: sumbv2
    real(fp) :: tkektg
    real(fp) :: tkescl ! Scaling factor for turbulent kinetic energy [m^2/s^2]
    real(fp) :: uref   ! U-velocity component as reference level KREF
    real(fp) :: velscl ! Scaling factor for velocity [m/s]
    real(fp) :: vis3
    real(fp) :: vismol ! Non-dimensional kinematic viscosity of water
    real(fp) :: visscl ! Scaling factor for viscosity [m^2/s]
    real(fp) :: vref
    real(fp) :: w0     ! Weighting factor for linear inter- polation
    real(fp) :: w1     ! Weighting factor for linear inter- polation
    real(fp) :: xlam
    real(fp) :: z0min
    real(fp) :: z1
    real(fp) :: zbv2   ! Weight factor for centre of gravity determination of velocity reference level KREF
    real(fp) :: zk     ! Vertical coordinate relative to reference plane distance ZETA below water surface
    real(fp) :: zlk    ! Kolmogoroff length scale
    real(fp) :: zsur   ! Distance between (u,v) velocity level and water surface
    real(fp) :: ztg    ! Distance, scaled by water depth, from water surface to w-level in inter- polation grid
    real(fp) :: zw     ! Distance, scaled by  water depth, from water surface to w-level in hydrodynamic grid
!
!! executable statements -------------------------------------------------------
!
    ag       => gdp%gdphysco%ag
    z0       => gdp%gdphysco%z0
    vonkar   => gdp%gdphysco%vonkar
    ricrit   => gdp%gdiwearr%ricrit
    bvmin    => gdp%gdiwearr%bvmin
    epsuko   => gdp%gdiwearr%epsuko
    epsbv2   => gdp%gdiwearr%epsbv2
    clw      => gdp%gdiwearr%clw
    ckw      => gdp%gdiwearr%ckw
    viscof   => gdp%gdiwearr%viscof
    !
    h0 = zeta + real(dps,fp)
    rilz = 0.0
    riav = 0.0
    do k = 1, kmax - 1
       if (rich(k)>ricrit) then
          dsig = 0.5*(thick(k + 1) + thick(k))
          rilz = rilz + dsig
          riav = riav + rich(k)*dsig
       endif
    enddo
    if (riav>0.0) riav = riav/rilz
    !
    if (rilz<siglim) then
       barcli = .false.
       rilz = rilz*h0
       goto 9999
    else
       barcli = .true.
       rilz = rilz*h0
    endif
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Interpolate vertical gradients and turbulence parameters,
    ! defined in w-points of DPM, to TG-grid:
    !
    dztg = 1./real(kmxt,fp)
    kdpm = 1
    zw = 1. - thick(1)
    !
    do ktg = 0, kmxt
       ztg = 1. - ktg*dztg
       !-->
   30  continue
       if (ztg<zw .and. kdpm<kmax) then
          kdpm = kdpm + 1
          zw = zw - thick(kdpm)
          goto 30
       ! <--
       endif
       w0 = max(0.0_fp, (ztg - zw)/thick(kdpm))
       w1 = 1. - w0
       kdpmd = kdpm - 1
       bv2(ktg) = w0*bruvai(kdpmd) + w1*bruvai(kdpm)
       tke(ktg) = w0*rtur0(kdpmd, 1) + w1*rtur0(kdpm, 1)
       eps(ktg) = w0*rtur0(kdpmd, 2) + w1*rtur0(kdpm, 2)
       edvis(ktg) = w0*vicww(kdpmd) + w1*vicww(kdpm)
    enddo
    !
    ! Depth-averaged velocity and vertical curvature of (u,v),
    ! defined in (u,v)-points of DPM grid:
    !
    umean = 0.0
    vmean = 0.0
    do k = 1, kmax
       kd = k - 1
       dz = thick(k)*h0
       d2u0(k) = (dudz(kd) - dudz(k))/dz
       d2v0(k) = (dvdz(kd) - dvdz(k))/dz
       umean = umean + thick(k)*u0(k)
       vmean = vmean + thick(k)*v0(k)
    enddo
    !
    ! Interpolate velocity and its curvature of hydrodynamic module, defined
    ! at (u,v) level, to TG-grid:
    !
    dztg = h0/kmxt
    kdpm = 1
    !
    do ktg = 1, kmxt - 1
       ztg = zeta - ktg*dztg
       ! -->
   80  continue
       zk = zeta + h0*sig(kdpm)
       if (ztg<=zk .and. kdpm<kmax) then
          kdpm = kdpm + 1
          goto 80
       ! <--
       endif
       if (kdpm==1) then
          ! log layer near surface:
          zsur = zeta - ztg
          ratio = log(zsur/(zeta - zk))
          utg(ktg) = u0(1) + ratio*ustwix/vonkar
          vtg(ktg) = v0(1) + ratio*ustwiy/vonkar
          d2u(ktg) = d2u0(1)
          d2v(ktg) = d2v0(1)
       elseif (ztg<zk) then
          ! log layer near bed:
          z0min = 1.E-6
          z1 = max(z0min, z0)
          ratio = log((real(dps,fp) + ztg + z1)/z1)/log((real(dps,fp) + zk + z1)/z1)
          utg(ktg) = ratio*u0(kmax)
          vtg(ktg) = ratio*v0(kmax)
          d2u(ktg) = d2u0(kmax)
          d2v(ktg) = d2v0(kmax)
       else
          ! internal levels:
          kdpmd = kdpm - 1
          w0 = (ztg - zk)/(0.5*h0*(thick(kdpm) + thick(kdpmd)))
          w1 = 1. - w0
          utg(ktg) = w0*u0(kdpmd) + w1*u0(kdpm)
          vtg(ktg) = w0*v0(kdpmd) + w1*v0(kdpm)
          d2u(ktg) = w0*d2u0(kdpmd) + w1*d2u0(kdpm)
          d2v(ktg) = w0*d2v0(kdpmd) + w1*d2v0(kdpm)
       endif
    enddo
    ! extrapolations to surface:
    utg(0) = 2*utg(1) - utg(2)
    vtg(0) = 2*vtg(1) - vtg(2)
    d2u(0) = d2u(1)
    d2v(0) = d2v(1)
    ! extrapolations to bed:
    utg(kmxt) = 2*utg(kmxt - 1) - utg(kmxt - 2)
    vtg(kmxt) = 2*vtg(kmxt - 1) - vtg(kmxt - 2)
    d2u(kmxt) = d2u(kmxt - 1)
    d2v(kmxt) = d2v(kmxt - 1)
    !
    ! Normalisation of buoyancy and velocity profiles for
    ! application in subroutine TG.
    !
    sumbv2 = 0.
    zbv2 = 0.
    ksum = 0
    do ktg = 0, kmxt
       if (bv2(ktg)>0) then
          ksum = ksum + 1
          sumbv2 = sumbv2 + bv2(ktg)
          zbv2 = zbv2 + ktg*bv2(ktg)
       endif
    enddo
    !
    ! Reference level for velocity is ZBVV2 and it equals
    ! centre-of-gravity of stably-stratified part of water
    ! column:
    !
    sumbv2 = max(sumbv2, 1.E-12_fp)
    kref = nint(zbv2/sumbv2)
    bv2av = sumbv2/ksum
    bvav = max(bvmin, sqrt(bv2av))
    bv2av = bvav**2
    !
    ! Normalisation of buoyancy freq. profile, velocity vector, its
    ! derivatives and turbulence state parameters. Length SCALE depends
    ! on water depth or horizontal IW wave length:
    !
    bv2mx = 0.
    do ktg = 0, kmxt
       bv2(ktg) = bv2(ktg)/bv2av
       if (bv2mx<bv2(ktg)) then
          bv2mx = bv2(ktg)
       endif
    enddo
    bvmx = sqrt(bv2mx)
    !
    xlam = 1./xk
    scale = min(h0, xlam)
    xkh = xk*scale
    !
    tkescl = (bvav*scale)**2
    epsscl = (bvav**3)*(scale**2)
    visscl = bvav*(scale**2)
    velscl = bvav*scale
    !
    if (useref) then
       uref = utg(kref)
       vref = vtg(kref)
    else
       uref = 0.0
       vref = 0.0
    endif
    !
    do ktg = 0, kmxt
       d2u(ktg) = d2u(ktg)*scale/bvav
       d2v(ktg) = d2v(ktg)*scale/bvav
       utg(ktg) = (utg(ktg) - uref)/velscl
       vtg(ktg) = (vtg(ktg) - vref)/velscl
       tke(ktg) = tke(ktg)/tkescl
       eps(ktg) = eps(ktg)/epsscl
       edvis(ktg) = edvis(ktg)/visscl
    enddo
    !
    ! Check for neutrally-stratified or unstably-stratified top layer:
    !
    bv2lvl = epsbv2*bv2mx
    ktg = -1
    ! -->
  120 continue
    ktg = ktg + 1
    diff = bv2(ktg) - bv2lvl
    if (diff<0 .and. ktg<kmxt) goto 120
    ! <--
    if (ktg<kmxt) then
       ktop = ktg
    else
       ktop = 0
    endif
    !
    ! Check for neutrally-stratified or unstably-stratified bed layer:
    !
    ktg = kmxt + 1
    ! -->
  130 continue
    ktg = ktg - 1
    diff = bv2(ktg) - bv2lvl
    if (diff<0 .and. ktg>0) goto 130
    ! <--
    if (ktg>0) then
       kbed = ktg
    else
       kbed = kmxt
    endif
    ! safety first:
    ktop = min(ktop, kbed)
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Determine turbulence parameters:
    !
    vismol = viscof/visscl
    vis3 = vismol**3
    zlo = 0.0
    !
    do ktg = 0, kmxt
       epsktg = max(epsuko, eps(ktg))
       tkektg = max(epsuko, tke(ktg))
       zlw(ktg) = clw*tkektg**1.5/epsktg
       !
       ! ZLW is truncated at Kolmogoroff length scale ZLK:
       !
       zlk = (vis3/epsktg)**0.25
       rmsw = ckw*sqrt(tkektg)
       tlw(ktg) = zlw(ktg)/rmsw
       zlw(ktg) = max(zlw(ktg), zlk)
       !
       ! Truncation ZLW for spuriously large ZLW without turbulence:
       !
       if (edvis(ktg)<vismol) zlw(ktg) = zlk
       if (ktop<=ktg .and. ktg<=kbed .and. bv2(ktg)>bvmin) then
          zlo = zlo + sqrt(epsktg/bv2(ktg)**1.5)
       endif
    enddo
    zlo = zlo/(kbed - ktop + 2)
    !
 9999 continue
end subroutine intpol

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 1998 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                          Spherepack                           *
!     *                                                               *
!     *       A Package of Fortran77 Subroutines and Programs         *
!     *                                                               *
!     *              for Modeling Geophysical Processes               *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *                  John Adams and Paul Swarztrauber             *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!    a program testing subroutines shpgi and shpg
!
!    shpg is the n**2 projection with complement, odd/even
!    factorization and zero truncation on a Gaussian distributed
!    grid. The projection is defined in the JCP paper "Generalized
!    discrete spherical harmonic transforms" by
!    Paul N. Swarztrauber and William F. Spotz
!    J. Comp. Phys., 159(2000) pp. 213-230.
!
!                     April 2002
!
program tshpg

    use, intrinsic :: ISO_Fortran_env, only: &
        stdout => OUTPUT_UNIT

    use spherepack

    ! Explicit typing only
    implicit none

    real(wp) :: discretization_error
    real(wp) :: g
    real(wp) :: ga
    real(wp) :: gb
    real(wp) :: gh
    integer(ip) :: i
    integer(ip) :: idimg
    integer(ip) :: idp
    integer(ip) :: ierror
    integer(ip) :: iprint
    integer(ip) :: isym
    integer(ip) :: iwshp
    integer(ip) :: j
    integer(ip) :: jdimg
    integer(ip) :: kdp
    integer(ip) :: liwshp
    integer(ip) :: lwork
    integer(ip) :: lwrk
    integer(ip) :: lwrk1
    integer(ip) :: lwsha
    integer(ip) :: lwshp
    integer(ip) :: lwshs
    integer(ip) :: mode
    integer(ip) :: mp1
    integer(ip) :: mtr
    integer(ip) :: mtrunc
    integer(ip) :: nlat
    integer(ip) :: nlon
    integer(ip) :: np1
    integer(ip) :: nt
    real(wp) :: sx
    real(wp) :: sy

    real(wp) :: toe
    real(wp) :: tusl
    real(wp) :: wrk1
    real(wp) :: wrk2
    real(wp) :: wshagc
    real(wp) :: wshp
    real(wp) :: wshsgc
    parameter (idp=8)
    parameter (kdp=idp+idp-2)
    parameter (lwshp=2*(idp+1)**2+kdp+20, &
        liwshp=4*(idp+1), lwrk=1.25*(idp+1)**2+7*idp+8)
    parameter (lwrk1=idp*kdp)
    parameter(lwork = 4*idp*(idp-1), &
        lwsha=idp*(4*idp+1)+idp+idp+15)
    real(wp), parameter :: ZERO = 0.0_wp
    real(wp) :: work(lwrk)
    dimension sx(idp, kdp), sy(idp, kdp), &
        wshp(lwshp), iwshp(liwshp), wrk1(lwrk1)
    dimension g(idp, kdp, 2), ga(idp, idp, 2), gb(idp, idp, 2), &
        gh(idp, kdp, 2), &
        wrk2(lwork), wshagc(lwsha), wshsgc(lwsha)

    iprint = 0
    nt = 1
    isym = 0
    mode = 0

    do nlat=6, 8
        do mtr=1, 2
            nlon = 2*(nlat-1)
            mtrunc = nlat-mtr
            mtrunc = min(mtrunc, nlat-1, nlon/2)
            idimg = idp
            jdimg = kdp
            call shagci(nlat, nlon, wshagc, ierror)
            call check_error(ierror)

            lwshs = lwsha
            call shsgci(nlat, nlon, wshsgc, ierror)
            call check_error(ierror)

            ! Initiate faster filter
            call shpgi(nlat, nlon, isym, mtrunc, wshp, lwshp, iwshp, liwshp, &
                work, lwrk, ierror)
            call check_error(ierror)

            if (iprint /= 0) write (stdout, '(3(a, i5))') &
                ' mode =' , nlon, '  nlat =', nlat, '  nlon =', nlon

            ! Initialize with pseudo random field
            do i=1, nlat
                do j=1, nlon
                    sx(i, j) =  cos(real(i*j, kind=wp))
                    g(i, j, 1) = sx(i, j)
                end do
            end do

            call shagc(nlat, nlon, mode, nt, g, idimg, jdimg, ga, gb, idimg, idimg, &
                wshagc, ierror)
            call check_error(ierror)

            if(mtrunc < nlat-1) then
                do np1=mtrunc+2, nlat
                    do mp1=1, np1
                        ga(mp1, np1, 1) = ZERO
                        gb(mp1, np1, 1) = ZERO
                    end do
                end do
            end if

            call shsgc(nlat, nlon, mode, nt, gh, idimg, jdimg, ga, gb, idimg, idimg, wshsgc, ierror)
            call check_error(ierror)

            call shpg(nlat, nlon, isym, mtrunc, sx, sy, idp, &
                wshp, lwshp, iwshp, liwshp, wrk1, lwrk1, ierror)
            call check_error(ierror)

            if (iprint > 0) write (stdout, '(/a/)') ' approx and exact solution'

            do j=1, nlon
                if(iprint > 0) write (stdout, 437) j, (sy(i, j), gh(i, j, 1), i=1, nlat)
437             format(' j=', i5, 1p4e15.6/(8x, 1p4e15.6))
            end do

            discretization_error = maxval(abs(sy(:nlat, :nlon)-gh(:nlat, :nlon, 1)))

            write (stdout, '(/2(a, i5)/)') &
                'case nlat =', nlat, ' and mtrunc =', mtrunc

            write (stdout, '(3(a, 1pe15.6/))') &
                ' error =', discretization_error, &
                ' tusl  =', tusl, &
                ' toe   =', toe
        end do
    end do
end program tshpg


text:
      SUBROUTINE CONCAL( NIN, NOUT, LERR,
     &			 CVECT, CENTRE,
     &			 ATMAX, ATNO, ATBRK, ATRES, ATCHN, ATRNO,
     &			 ATXYZ, ATVDW, ATBND,
     &                   SPDBSP, FPDBSP,
     &			 SHORTO,
     &			 ENDRAD, CUTSIZE, CONNR,
     &			 PERPVE, PERPVN, REQUIV, PEG_WRITEALL, IREC,
     &                   SAMPLE)
      IMPLICIT NONE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Oliver Smart. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Oliver Smart are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 2000 Oliver Smart *
C *                                                                  *
C ********************************************************************
C
C s/r to perform connolly type surface calc within hole
C
C Modification history:
C
C Date  Author          Modification
C 01/00 O.S.S.    	CONN option introduced
C 10/00 O.S.S.		PEG_WRITEALL option added for PEG effect calcs

C passed vbles *****************

C input/output stream numbers
      INTEGER			NIN, NOUT

C error indicator
      LOGICAL			LERR

C a vector in the direction of the channel
C must be a unit vector - return unchanged
      DOUBLE PRECISION          CVECT(3)

C The centre about which the connolly grid is to be
C done (s/r holcal normal supplies lowcen)
      DOUBLE PRECISION          CENTRE(3)

C maximum no. of atoms
      INTEGER                   ATMAX

C number of atoms read in from each file:
      INTEGER                   ATNO

C atom names found in Brookhaven file:
      CHARACTER*4               ATBRK(ATMAX)

C residue name in brook
      CHARACTER*3               ATRES(ATMAX)

C chain identifier in brookhaven pdb file
      CHARACTER*1               ATCHN(ATMAX)

C integer residue no.
      INTEGER                   ATRNO(ATMAX)

C co-ordinates
      DOUBLE PRECISION          ATXYZ( 3, ATMAX)

C vdw and bond radii of each atoms
      DOUBLE PRECISION          ATVDW(ATMAX), ATBND(ATMAX)

C pdb format output file 'atoms' are sphere centres
C default to none by setting FPDBSP to 'NONE'
      CHARACTER*200             FPDBSP
      INTEGER                   SPDBSP

C 11/95 introduce new variable to reduce text output -
C if oshort =
C 0: full text output
C 1: no in run progress reported
C 2: ditto plus no graph type output - only leaving
C    Minimum radius & conductance calculations.
C 3: All text output other than input card mirroring
C    and error messages turned off
      INTEGER                   SHORTO

C ENDRAD: Radius above which a point is considered to have escaped
C ENDRP3 also need to have a slightly higher value for the purposes of this s/r
      DOUBLE PRECISION          ENDRAD, ENDRP3

C cut off list control vble, -ve: no cutoff lists used in calculation
C otherwise this is the additional radius above the sphere radius
C to be used for cutoff list (see s/r holeen)
      DOUBLE PRECISION		CUTSIZE

C option CONN introduced 1/00 -
C CONNR(1) is the connolly surface probe radius (default 1.15Angs)
C          the option is turned off by setting this radius negative
C CONNR(2) is the grid size to be used in the calculation - default .25 angs
      DOUBLE PRECISION		CONNR(2)

C unit vectors in plane normal to cvect
C needed to a number of routines - originally h2dmap but now conn
C and addend
C these are normal - first is east direction (x on plane)
C                    and second is north     (y on plane)
C return unchanged!!!
      DOUBLE PRECISION		PERPVE(3), PERPVN(3)

C area of final assembly returned as a effective radius = sqrt(area/pi)
      DOUBLE PRECISION		REQUIV

C write all the records? (for PEG calcs)
      LOGICAL 			PEG_WRITEALL

C record number for sph write
      INTEGER                   IREC

! a sampling distance - the distance between planes
      DOUBLE PRECISION		SAMPLE

C internal variables ***********

C a new point
      DOUBLE PRECISION          NEWCEN(3), NEWENG

C atom list no. with smallest dist-vdw radius, 2nd smallest
C 07/06/94 as iat's may or may not be supplied with previously
C found numbers then we can use this to speed up the procedure
C Dec 97 add iat3
      INTEGER			IAT1, IAT2, IAT3

C 2nd/3rd smallest distance-vdw radius (of iat2/3)
      DOUBLE PRECISION		DAT2, DAT3

C arrays for connolly calc
C maximum size of the array STRMAX
      INTEGER			SCOMAX
      PARAMETER(		SCOMAX= 10000)
C number of entries stored
      INTEGER			SCONUM
C number of entries - that have been "checked"
C that is looking at the neighbouring eight points
C
C      INTEGER			SCOCHK
C the next entry to check
      INTEGER			SCONXT

C a stored entry - a point with a radius above the probe radius connr(1)
      DOUBLE PRECISION		SCOXYZ( 3, SCOMAX)
C the radius of the point
      DOUBLE PRECISION		SCORAD( SCOMAX)
C logical vble whether point needs to be checked\
C if point is on grid then it must be if
C  produced by the spike procedure it need no
C
C at the end of the procedure change meaning -
C it then becomes a store of whether point should be
C included in output/ area calculations
      LOGICAL			SCOTDO( SCOMAX)
C s/r coarea needs a store for the east/north coords of
C each point on grid - these are measured relative to the
C initial point
! SCOCIRCRAD: the radius of the point in the plane under consideration
! in the case of calls from concal this we equal scocrad
      DOUBLE PRECISION		SCOENC( 2, SCOMAX),
     &				SCOCIRCRAD( SCOMAX)

      ! vbles needed by coarea - not used here only required
      ! for peg calls to coarea by sph_process
      INTEGER			SCO_PEG_ACC(SCOMAX)
      LOGICAL			SCO_PEG_FLAG(SCOMAX)

CORIGC the status of the  surronding points
CORIGC if .false. then point in relevant direction has been accepted or not yet checked
CORIGC if .true. then last point in that direction has been rejected
CORIGC first variable is "East" and next is "North":
CORIGC i.e. point 0,0 is the point itself
CORIGC and directions:
CORIGC  (-1,+1) ( 0,+1) (+1,+1)
CORIGC  (-1, 0)  Cpoint (+1, 0)
CORIGC  (-1,-1) ( 0,-1) (+1,-1)
CORIG      LOGICAL			SCONEI( -1:1, -1:1, SCOMAX)
CORIG
CORIGC number of new points
CORIG INTEGER			SCONEW

C GRID SIZE, working variable - dist2
      DOUBLE PRECISION 		GRIDS, DIST2, DIFF1, DIFF2, DIFF3

C loop counts
      INTEGER			NCOUNT, ECOUNT, SCOUNT, JCOUNT, XCOUNT,
     &				CCOUNT

C test point
      DOUBLE PRECISION		TSTXYZ(3)

C variables for spike like routine
      DOUBLE PRECISION		VDISP(3), DELTA

C max radius indicator
      DOUBLE PRECISION		MAXRAD

C number of points eliminated, elimination distance
      INTEGER			NELIM
      DOUBLE PRECISION		DELIM2

! for writing .sph file correctly
      INTEGER			WRESNO
      DOUBLE PRECISION 		WREQUIV

C end of decs ******************

C bigger value of endrad
      ENDRP3 = ENDRAD + 3D0

      IF (SHORTO.LT.2) WRITE( NOUT, '(A)')
     & ' Connolly routine for this plane'
D	WRITE( NOUT, '(//A/ A,F8.3/ A,F8.3/ A,3F8.3)')
D    & 'D-line: Call to s/r concal',
D    &' D-line: CONNR(1) = ', CONNR(1),
D    &' D-line: CONNR(2) = ', CONNR(2),
D    &' D-line: CENTRE = ', CENTRE(1), CENTRE(2), CENTRE(3)

C find the pore radius of the initial point
      CALL HOLEEN( CENTRE, NEWENG, IAT1, IAT2, IAT3, DAT2, DAT3,
     &           ATMAX, ATNO, ATXYZ, ATVDW, CUTSIZE)

D      WRITE( NOUT, '(A,F8.3,A)')
D    & 'D-line pore radius of supplied CENTRE = ', -NEWENG, ' Angs'

C is the hole pore radius less than the probe radius - if so cannot do connolly
       IF (-NEWENG.LT.CONNR(1)) THEN
	   WRITE(NOUT, '(/ A/ A,F8.3,A,F8.3/ A)')
     &' Connolly routine (s/r concal)',
     &'   initial point probe radius=', -NEWENG,
     &              ' less than probe radius=',CONNR(1),
     &'   So using HOLE point for calcs........'

C simply return the hole radius as the effective
        REQUIV = -NEWENG

      ELSE
C initial point has a pore radius above the probe radius so can check
C initialize array - one point - the initial
        SCONUM = 1
C next point to check is the first
        SCONXT = 1
        SCOXYZ(1,SCONUM) = CENTRE(1)
        SCOXYZ(2,SCONUM) = CENTRE(2)
        SCOXYZ(3,SCONUM) = CENTRE(3)
        SCORAD(SCONUM) = -NEWENG
	SCOTDO(SCONUM) = .TRUE.
C gridsize - set outside this routine
C (n.b. neighbouring points must overlap)
        GRIDS = CONNR(2)

C jump back here to check next unchecked point
222    CONTINUE
C ignore points that have a radius above ENDRP3 - these should not spawn more
C and points that have been generated by spike procedure
        IF ( (SCORAD(SCONXT).LT.ENDRP3) .AND.
     &              SCOTDO(SCONXT)            ) THEN
C mark the point as checked
          SCOTDO(SCONXT) = .FALSE.
C check points around scoxyz(1,SCONXT), scoxyz(2,SCONXT), scoxyz(3,SCONXT)
C increment south, zero, north
          DO 15 NCOUNT = -1, 1
            DO 25 ECOUNT = -1, 1
C test point:
              TSTXYZ(1) = SCOXYZ(1,SCONXT)
              TSTXYZ(2) = SCOXYZ(2,SCONXT)
              TSTXYZ(3) = SCOXYZ(3,SCONXT)
C movement east by grids
              TSTXYZ(1) = TSTXYZ(1) + REAL(ECOUNT)*GRIDS*PERPVE(1)
              TSTXYZ(2) = TSTXYZ(2) + REAL(ECOUNT)*GRIDS*PERPVE(2)
              TSTXYZ(3) = TSTXYZ(3) + REAL(ECOUNT)*GRIDS*PERPVE(3)
C movement north by grids
              TSTXYZ(1) = TSTXYZ(1) + REAL(NCOUNT)*GRIDS*PERPVN(1)
              TSTXYZ(2) = TSTXYZ(2) + REAL(NCOUNT)*GRIDS*PERPVN(2)
              TSTXYZ(3) = TSTXYZ(3) + REAL(NCOUNT)*GRIDS*PERPVN(3)
C Q: have we already produced this point?
C look through existing list
              DO 67 SCOUNT = 1, SCONUM
C find distance**2 between the points
                DIFF1 = TSTXYZ(1)-SCOXYZ(1,SCOUNT)
                DIFF2 = TSTXYZ(2)-SCOXYZ(2,SCOUNT)
                DIFF3 = TSTXYZ(3)-SCOXYZ(3,SCOUNT)
                DIST2 = DIFF1*DIFF1 + DIFF2*DIFF2 + DIFF3+DIFF3
                IF (DIST2.LT.GRIDS/1000D0) THEN
C this point already considered and accepted
CORIG                  SCONEI( ECOUNT, NCOUNT, SCONXT) = .FALSE.
D                 WRITE( NOUT, '(A,3F8.3)')
D    & 'D-line Have already considered + accepted point= ',
D    &       TSTXYZ(1), TSTXYZ(2),TSTXYZ(3)
C jump to next point
                  GOTO 25
                ENDIF
67            CONTINUE

C find the radius of this point
              CALL HOLEEN( TSTXYZ, NEWENG, IAT1, IAT2, IAT3, DAT2, DAT3,
     &                   ATMAX, ATNO, ATXYZ, ATVDW, CUTSIZE)
C change sign to give +ve pore radius
              NEWENG = -NEWENG
C Q: is this an accepted point?
              IF (NEWENG.GE.CONNR(1)) THEN
C YES this is an accepted point (and its unique)
CORIGC first store information for point that produced this
CORIG                SCONEI( ECOUNT, NCOUNT, SCONXT) = .FALSE.

C now add point onto end of list to be considered
C array check here
                IF (SCONUM+1.GT.SCOMAX) THEN
                  WRITE( NOUT, '(A)')
     &' ERROR SCOMAX reached in s/r concal',
     &' this should not generally happen (unless grid set too small)!',
     &' please report to Oliver Smart o.s.smart@bham.ac.uk',
     &' with input and output files'
                  LERR = .TRUE.
                  GOTO 55555
                ENDIF
C store the extra point
                SCONUM = SCONUM + 1
                SCOXYZ(1,SCONUM) = TSTXYZ(1)
                SCOXYZ(2,SCONUM) = TSTXYZ(2)
                SCOXYZ(3,SCONUM) = TSTXYZ(3)
                SCORAD(SCONUM) = NEWENG
C this is a not spike so must be checked
		SCOTDO(SCONUM) = .TRUE.
D               WRITE( NOUT, '(A,3F8.3,A,F8.3,A)')
D    & 'D-line ACCEPT point= ', TSTXYZ(1), TSTXYZ(2),TSTXYZ(3),
D    &               ' radius= ', NEWENG, ' Angs'
C must finish as soon as a ENDRP3 point arrived at
                IF (NEWENG.GT.ENDRP3) THEN
		  IF (SHORTO.LT.2) WRITE( NOUT, '(A)')
     &' Escaped to outside at this level'
                  GOTO 333
		ENDIF

              ELSE
C NO this point is rejected (Q: is this an accepted point?)
CORIGC first store information for point that produced this
CORIG                SCONEI( ECOUNT, NCOUNT, SCONXT) = .TRUE.
D               WRITE( NOUT, '(A,3F8.3,A,F8.3,A)')
D    & 'D-line REJECT point= ', TSTXYZ(1), TSTXYZ(2),TSTXYZ(3),
D    &               ' radius= ', NEWENG, ' Angs'
C "spike" procedure - adjust point towards originating vector until
C it gets to desired radius
C vector starts from SCOXYZ(*,SCONXT) and goes to TSTXYZ(*)
                VDISP(1) = TSTXYZ(1)-SCOXYZ(1,SCONXT)
                VDISP(2) = TSTXYZ(2)-SCOXYZ(2,SCONXT)
                VDISP(3) = TSTXYZ(3)-SCOXYZ(3,SCONXT)
C make unit
                CALL DUVEC2( VDISP)
C DISP is the proportion displacement along this vector
C will adjust this to give a pore radius equal to the probe
C first estimate approx way along
                DELTA = 0.25D0*GRIDS
C next bit is analogous to procedure in hodotu.f for producing spikes
C do ten cycles of refinement to find sphere centre SEE f048 FOR METHOD
                DO 55 CCOUNT = 1, 100
C load up current point dum is the distance along line
C from the sphere which is current estimate of overlap point.
                  TSTXYZ(1) = SCOXYZ(1,SCONXT) + DELTA*VDISP(1)
                  TSTXYZ(2) = SCOXYZ(2,SCONXT) + DELTA*VDISP(2)
                  TSTXYZ(3) = SCOXYZ(3,SCONXT) + DELTA*VDISP(3)
C find the radius of this point
                  CALL HOLEEN( TSTXYZ, NEWENG, IAT1, IAT2, IAT3,
     &			       DAT2, DAT3,
     &                        ATMAX, ATNO, ATXYZ, ATVDW, CUTSIZE)
C make neweng=poreradius
                  NEWENG = -NEWENG
D               WRITE( NOUT, '(A,I5,AF8.3,A,F8.3,A)')
D    & 'D-line SPIKE cycle ',CCOUNT,' DELTA= ', DELTA,
D    &                    ' radius= ', NEWENG, ' Angs'
C take off the probe radius plus a very little bit
                  NEWENG =  NEWENG - (CONNR(1)+0.0001D0)
		  IF (ABS(NEWENG).LT.0.0005) GOTO 56
C if this number is greater than zero then we are too close to point SCONXT
C so add number on
                  DELTA = DELTA + NEWENG
C as we go thru this cycle the pore radius at the point should get closer to probe radius
55              CONTINUE
C at the end should have a point close to the probe radius
C store it if good enough
C very occasionally delta goes less than zero - this can produce big jumps
C so eliminate these points
56              CONTINUE
                IF ((ABS(NEWENG).LT.0.0005).AND.DELTA.GT.0.0) THEN
C check to see whether point is close to an already existing point
C look through existing list
                   DO 667 SCOUNT = 1, SCONUM
C find distance**2 between the points
                     DIFF1 = TSTXYZ(1)-SCOXYZ(1,SCOUNT)
                     DIFF2 = TSTXYZ(2)-SCOXYZ(2,SCOUNT)
                     DIFF3 = TSTXYZ(3)-SCOXYZ(3,SCOUNT)
                     DIST2 = DIFF1*DIFF1 + DIFF2*DIFF2 + DIFF3+DIFF3
C distance must be above 0.3 angs
                     IF (DIST2.LT.0.09D0) THEN
C this point already considered and accepted
CORIG                  SCONEI( ECOUNT, NCOUNT, SCONXT) = .FALSE.
D                 WRITE( NOUT, '(A,3F8.3)')
D    & 'D-line Point close (.3 angs) to already stored point '
C jump to next point
                       GOTO 25
                     ENDIF
667               CONTINUE

C array check here
                  IF (SCONUM+1.GT.SCOMAX) THEN
                    WRITE( NOUT, '(A)')
     &' ERROR SCOMAX reached in s/r concal',
     &' this should not generally happen (unless grid set too small)!',
     &' please report to Oliver Smart o.s.smart@bham.ac.uk',
     &' with input and output files'
                    LERR = .TRUE.
                    GOTO 55555
                  ENDIF
C store point
                  SCONUM = SCONUM + 1
                  SCOXYZ(1,SCONUM) = TSTXYZ(1)
                  SCOXYZ(2,SCONUM) = TSTXYZ(2)
                  SCOXYZ(3,SCONUM) = TSTXYZ(3)
C get back to the pore radius
                  SCORAD(SCONUM) = NEWENG + (CONNR(1)+0.0001D0)
C this is a spike - so should not produce any new point
		  SCOTDO(SCONUM) = .FALSE.
	        ENDIF
              ENDIF
25          CONTINUE
15        CONTINUE
C finished considering point SCONXT
        ENDIF

C jump back to consider next point if necessary
C at each stage the highest radius point propogates
C this is to keep number of points down when we get escape to infinity
C initialize maximum radius array
        MAXRAD = -1D10
	SCONXT = -1000
C go thru all the point
        DO 701 SCOUNT = 1, SCONUM
          IF (SCOTDO(SCOUNT)) THEN
	     IF (SCORAD(SCOUNT).GT.MAXRAD) THEN
	       SCONXT = SCOUNT
	       MAXRAD = SCORAD(SCOUNT)
	     ENDIF
	  ENDIF
701     CONTINUE

C do we have to jump back - are there any points to check?
        IF (SCONXT.GT.0) THEN
D         	WRITE( NOUT, '(A,I5,A,F8.3)')
D    & 'D-line Jump back next point SCONXT= ', SCONXT, ' MAXRAD= ', MAXRAD
          GOTO 222
        ENDIF

C Jump out here if ENDRP3 is hit
333     CONTINUE
C tell user vectors used to define rotation
        IF (SHORTO.LT.2) WRITE( NOUT, '(A,I5,A)')
     & ' Have stored total of ', SCONUM, ' points'

C eliminate unnecessary points.
C use SCOTDO as a flag for consideration
C set all to true
        DO 702 SCOUNT = 1, SCONUM
          SCOTDO(SCOUNT) = .TRUE.
702     CONTINUE
C if we have found an end (which will allways be the last point found)
C then eliminate all points whose centres lie within it
        IF (SCORAD(SCONUM).GT.ENDRP3) THEN
	  NELIM = 0
          DELIM2 = SCORAD(SCONUM)
	  DELIM2 = DELIM2*DELIM2
C go thru points - but not the last
          DO 703 SCOUNT = 1, SCONUM-1
C find distance**2 between the point and end point
            DIFF1 = SCOXYZ(1,SCONUM)-SCOXYZ(1,SCOUNT)
            DIFF2 = SCOXYZ(2,SCONUM)-SCOXYZ(2,SCOUNT)
            DIFF3 = SCOXYZ(3,SCONUM)-SCOXYZ(3,SCOUNT)
            DIST2 = DIFF1*DIFF1 + DIFF2*DIFF2 + DIFF3+DIFF3
C within limits?
            IF (DIST2.LT.DELIM2) THEN
	      SCOTDO(SCOUNT) = .FALSE.
	      NELIM = NELIM + 1
	    ENDIF
703       CONTINUE
          IF (SHORTO.LT.2) WRITE( NOUT, '(A,I5,A)')
     & ' Have eliminated ', NELIM,
     &    ' of these points as they are within the end found',
     & ' leaves ', SCONUM-NELIM, ' points '
	ENDIF

C 2nd cycle of elimination - points that lie well within original HOLE point
C that is with 1.5 probe radii
        IF (SCORAD(1).GT.2.5D0*CONNR(1)) THEN
	  NELIM = 0
          DELIM2 = SCORAD(1)-1.5D0*CONNR(1)
	  DELIM2 = DELIM2*DELIM2

C go thru points - but not the first
          DO 704 SCOUNT = 2, SCONUM
	    IF (.NOT.SCOTDO(SCOUNT)) GOTO 704
C find distance**2 between the point and original point
            DIFF1 = SCOXYZ(1,1)-SCOXYZ(1,SCOUNT)
            DIFF2 = SCOXYZ(2,1)-SCOXYZ(2,SCOUNT)
            DIFF3 = SCOXYZ(3,1)-SCOXYZ(3,SCOUNT)
            DIST2 = DIFF1*DIFF1 + DIFF2*DIFF2 + DIFF3+DIFF3
C within limits?
            IF (DIST2.LT.DELIM2) THEN
	      SCOTDO(SCOUNT) = .FALSE.
	      NELIM = NELIM + 1
	    ENDIF
704       CONTINUE
          IF (SHORTO.LT.2) WRITE( NOUT, '(A,I5,A)')
     & ' Have eliminated ', NELIM,
     &    ' points as they are well within original HOLE centre',
     & ' leaves ', SCONUM-NELIM, ' points '
	ENDIF

C now do the area calculation for this set
        CALL COAREA( REQUIV, NOUT, LERR, SHORTO, ENDRAD,
     &		     CENTRE, PERPVE, PERPVN,
     &               SCOMAX, SCONUM, SCOXYZ, SCORAD, SCOTDO,
     &               SCOENC, SCOCIRCRAD, CVECT, SAMPLE,
     &		     .FALSE., SCO_PEG_ACC, SCO_PEG_FLAG, -1)

C write out the records if we need to
        IF (FPDBSP.NE.'NONE') THEN
	  WRESNO = IREC ! first record is the centre supplied
	                ! to routine write this as a normal record
	  WREQUIV = REQUIV
	  IF (WREQUIV.GT.1000.) WREQUIV = 999.99
	  DO 4100 SCOUNT = 1, SCONUM

C PEG_WRITEALL option do not edit out records
C set the record number to -999 for all records
             IF (       SCOTDO(SCOUNT) .OR.
     &            (SCORAD(SCOUNT).GT.ENDRAD) .OR.
     &                    PEG_WRITEALL                ) THEN
	       WRITE( SPDBSP, '(A,I4,4X,3F8.3,2F6.2)')
     &   'ATOM      1  QSS SPH S', WRESNO,
     &       (SCOXYZ(XCOUNT,SCOUNT), XCOUNT= 1, 3),
     &       SCORAD(SCOUNT), WREQUIV
               IF (SCORAD(SCOUNT).GT.ENDRAD)
     &           WRITE( SPDBSP, '(A)') 'LAST-REC-END'
               WRESNO = -999 ! after first write mark as connolly type records
             ENDIF
4100      CONTINUE
	ENDIF

C end of if Q is the hole pore radius less than the probe radius?
      ENDIF

C      STOP 'temp stop in HOLCAL'

55555 RETURN
CORIGD       DO 5010 SCOUNT = 1, SCONUM
CORIGD          WRITE( NOUT, '(A,3F8.3,A,F8.3,A/ 3L1/ 3L1/ 3L1/)')
CORIGD    & ' xyz=', (SCOXYZ(XCOUNT,SCOUNT), XCOUNT = 1,3),
CORIGD    &   '  rad=', SCORAD(SCOUNT),
CORIGD    &   '   flags=',
CORIGD    &   SCONEI( -1,  1, SCOUNT),
CORIGD    &   SCONEI(  0,  1, SCOUNT),
CORIGD    &   SCONEI( +1,  1, SCOUNT),
CORIGD    &   SCONEI( -1,  0, SCOUNT),
CORIGD    &   SCONEI(  0,  0, SCOUNT),
CORIGD    &   SCONEI( +1,  0, SCOUNT),
CORIGD    &   SCONEI( -1, -1, SCOUNT),
CORIGD    &   SCONEI(  0, -1, SCOUNT),
CORIGD    &   SCONEI( +1, -1, SCOUNT)
CORIGD5010 CONTINUE

CORIGC original procedure - involved halving gridsize but failed
CORIGC just in case leave in commented CORIG
CORIGC have finished original gridsize - halve it and then
CORIGC check the half points
CORIG        IF (GRIDS.GT.CONNR(1)/1.999) THEN
CORIG
CORIGC halve the grid size
CORIG          GRIDS = GRIDS/2D0
CORIG
CORIGD         WRITE( NOUT, '(A,F10.4)')
CORIGD    & 'D-line GRIDS reduced to = ', GRIDS
CORIG
CORIGC go thru the existing list
CORIG          SCONEW = SCONUM
CORIG          DO 89 SCOUNT = 1, SCONUM
CORIGC look at each direction of stored point
CORIGC increment south, zero, north
CORIG            DO 157 NCOUNT = -1, 1
CORIG              DO 257 ECOUNT = -1, 1
CORIGC if the flag SCONEI( ECOUNT, NCOUNT, SCOUNT) is .true.
CORIGC then this point was rejected for previous grid size
CORIGC so check again
CORIG                IF (SCONEI( ECOUNT, NCOUNT, SCOUNT)) THEN
CORIGC test point for this direction with new grid
CORIG                  TSTXYZ(1) = SCOXYZ(1,SCOUNT)
CORIG                  TSTXYZ(2) = SCOXYZ(2,SCOUNT)
CORIG                  TSTXYZ(3) = SCOXYZ(3,SCOUNT)
CORIGC movement east by grids
CORIG                  TSTXYZ(1) = TSTXYZ(1) + REAL(ECOUNT)*GRIDS*PERPVE(1)
CORIG                  TSTXYZ(2) = TSTXYZ(2) + REAL(ECOUNT)*GRIDS*PERPVE(2)
CORIG                  TSTXYZ(3) = TSTXYZ(3) + REAL(ECOUNT)*GRIDS*PERPVE(3)
CORIGC movement north by grids
CORIG                  TSTXYZ(1) = TSTXYZ(1) + REAL(NCOUNT)*GRIDS*PERPVN(1)
CORIG                  TSTXYZ(2) = TSTXYZ(2) + REAL(NCOUNT)*GRIDS*PERPVN(2)
CORIG                  TSTXYZ(3) = TSTXYZ(3) + REAL(NCOUNT)*GRIDS*PERPVN(3)
CORIGC Q: have we already produced this point?
CORIGC look through existing list including new points
CORIG                  DO 677 JCOUNT = 1, SCONEW
CORIGC find distance**2 between the points
CORIG                    DIFF1 = TSTXYZ(1)-SCOXYZ(1,JCOUNT)
CORIG                    DIFF2 = TSTXYZ(2)-SCOXYZ(2,JCOUNT)
CORIG                    DIFF3 = TSTXYZ(3)-SCOXYZ(3,JCOUNT)
CORIG                    DIST2 = DIFF1*DIFF1 + DIFF2*DIFF2 + DIFF3+DIFF3
CORIG                    IF (DIST2.LT.GRIDS/1000D0) THEN
CORIGC this point already considered and accepted
CORIG                      SCONEI( ECOUNT, NCOUNT, SCOCHK) = .FALSE.
CORIGD                     WRITE( NOUT, '(A,3F8.3)')
CORIGD    & 'D-line 2nd loop Have already considered + accepted point= ',
CORIGD    &       TSTXYZ(1), TSTXYZ(2),TSTXYZ(3)
CORIGC jump to next point
CORIG                     GOTO 257
CORIG                   ENDIF
CORIG677              CONTINUE
CORIG
CORIGC find the radius of this point
CORIG                 CALL HOLEEN( TSTXYZ, NEWENG, IAT1, IAT2, IAT3,
CORIG     &			      DAT2, DAT3,
CORIG     &                        ATMAX, ATNO, ATXYZ, ATVDW, CUTSIZE)
CORIGC change sign to give +ve pore radius
CORIG                 NEWENG = -NEWENG
CORIGC Q: is this an accepted point?
CORIG                 IF (NEWENG.GE.CONNR(1)) THEN
CORIGC YES this is an accepted point (and its unique)
CORIGC first store information for point that produced this
CORIG                   SCONEI( ECOUNT, NCOUNT, SCOUNT) = .FALSE.
CORIGC now add point onto end of list to be considered
CORIGC array check here
CORIG                   IF (SCONEW+1.GT.SCOMAX) THEN
CORIG                     WRITE( NOUT, '(A)')
CORIG     &' ERROR SCOMAX reached in s/r concal',
CORIG     &' this should not generally happen!',
CORIG     &' please report to Oliver Smart o.s.smart@bham.ac.uk',
CORIG     &' with input and output files'
CORIG                     LERR = .TRUE.
CORIG                     GOTO 55555
CORIG                  ENDIF
CORIGC store the extra point
CORIG                  SCONEW = SCONEW + 1
CORIG                  SCOXYZ(1,SCONEW) = TSTXYZ(1)
CORIG                  SCOXYZ(2,SCONEW) = TSTXYZ(2)
CORIG                  SCOXYZ(3,SCONEW) = TSTXYZ(3)
CORIG                  SCORAD(SCONEW) = NEWENG
CORIGD                 WRITE( NOUT, '(A,3F8.3,A,F8.3,A)')
CORIGD    & 'D-line 2nd loop ACCEPT point= ', TSTXYZ(1), TSTXYZ(2),TSTXYZ(3),
CORIGD    &               ' radius= ', NEWENG, ' Angs'
CORIG                  ELSE
CORIGC NO this point is rejected (Q: is this an accepted point?)
CORIGC first store information for point that produced this
CORIG                    SCONEI( ECOUNT, NCOUNT, SCOUNT) = .TRUE.
CORIGD                   WRITE( NOUT, '(A,3F8.3,A,F8.3,A)')
CORIGD    & 'D-line REJECT point= ', TSTXYZ(1), TSTXYZ(2),TSTXYZ(3),
CORIGD    &               ' radius= ', NEWENG, ' Angs'
CORIG                  ENDIF
CORIG                ENDIF
CORIG257           CONTINUE
CORIG157        CONTINUE
CORIG89        CONTINUE
CORIGC finished adding extra points for half grids
CORIG          SCONUM = SCONEW
CORIG
CORIGC jump back to consider next point if necessary
CORIG          IF (SCOCHK.LT.SCONUM) THEN
CORIGD           WRITE( NOUT, '(A,I5,A,I5)')
CORIGD    & ' Jump back as SCOCHK= ', SCOCHK, ' SCONUM= ', SCONUM
CORIG            GOTO 222
CORIG          ENDIF
CORIG        ENDIF
      END


text:
      MODULE frame_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:15  03/09/06
      SUBROUTINE frame (FMAT, NUMAT, MODE)
      USE vast_kind_param,ONLY: DOUBLE
      REAL(DOUBLE), DIMENSION(*), INTENT(INOUT) :: FMAT
      INTEGER, INTENT(IN) :: NUMAT
      INTEGER, INTENT(IN) :: MODE
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      SUBROUTINE GWAVE (ERV,EZV,ERH,EZH,EPH)
C
C     GWAVE COMPUTES THE ELECTRIC FIELD, INCLUDING GROUND WAVE, OF A
C     CURRENT ELEMENT OVER A GROUND PLANE USING FORMULAS OF K.A. NORTON
C     (PROC. IRE, SEPT., 1937, PP.1203,1236.)
C
      COMPLEX FJ,TPJ,U2,U,RK1,RK2,T1,T2,T3,T4,P1,RV,OMR,W,F,Q1,RH,V,G,XR
     11,XR2,X1,X2,X3,X4,X5,X6,X7,EZV,ERV,EZH,ERH,EPH,XX1,XX2,ECON,FBAR
      COMMON /GWAV/ U,U2,XX1,XX2,R1,R2,ZMH,ZPH
      DIMENSION FJX(2), TPJX(2), ECONX(2)
      EQUIVALENCE (FJ,FJX), (TPJ,TPJX), (ECON,ECONX)
      DATA PI/3.141592654/,FJX/0.,1./,TPJX/0.,6.283185308/
      DATA ECONX/0.,-188.367/
      SPPP=ZMH/R1
      SPPP2=SPPP*SPPP
      CPPP2=1.-SPPP2
      IF (CPPP2.LT.1.E-20) CPPP2=1.E-20
      CPPP=SQRT(CPPP2)
      SPP=ZPH/R2
      SPP2=SPP*SPP
      CPP2=1.-SPP2
      IF (CPP2.LT.1.E-20) CPP2=1.E-20
      CPP=SQRT(CPP2)
      RK1=-TPJ*R1
      RK2=-TPJ*R2
      T1=1.-U2*CPP2
      T2=CSQRT(T1)
      T3=(1.-1./RK1)/RK1
      T4=(1.-1./RK2)/RK2
      P1=RK2*U2*T1/(2.*CPP2)
      RV=(SPP-U*T2)/(SPP+U*T2)
      OMR=1.-RV
      W=1./OMR
      W=(4.,0.)*P1*W*W
      F=FBAR(W)
      Q1=RK2*T1/(2.*U2*CPP2)
      RH=(T2-U*SPP)/(T2+U*SPP)
      V=1./(1.+RH)
      V=(4.,0.)*Q1*V*V
      G=FBAR(V)
      XR1=XX1/R1
      XR2=XX2/R2
      X1=CPPP2*XR1
      X2=RV*CPP2*XR2
      X3=OMR*CPP2*F*XR2
      X4=U*T2*SPP*2.*XR2/RK2
      X5=XR1*T3*(1.-3.*SPPP2)
      X6=XR2*T4*(1.-3.*SPP2)
      EZV=(X1+X2+X3-X4-X5-X6)*ECON
      X1=SPPP*CPPP*XR1
      X2=RV*SPP*CPP*XR2
      X3=CPP*OMR*U*T2*F*XR2
      X4=SPP*CPP*OMR*XR2/RK2
      X5=3.*SPPP*CPPP*T3*XR1
      X6=CPP*U*T2*OMR*XR2/RK2*.5
      X7=3.*SPP*CPP*T4*XR2
      ERV=-(X1+X2-X3+X4-X5+X6-X7)*ECON
      EZH=-(X1-X2+X3-X4-X5-X6+X7)*ECON
      X1=SPPP2*XR1
      X2=RV*SPP2*XR2
      X4=U2*T1*OMR*F*XR2
      X5=T3*(1.-3.*CPPP2)*XR1
      X6=T4*(1.-3.*CPP2)*(1.-U2*(1.+RV)-U2*OMR*F)*XR2
      X7=U2*CPP2*OMR*(1.-1./RK2)*(F*(U2*T1-SPP2-1./RK2)+1./RK2)*XR2
      ERH=(X1-X2-X4-X5+X6+X7)*ECON
      X1=XR1
      X2=RH*XR2
      X3=(RH+1.)*G*XR2
      X4=T3*XR1
      X5=T4*(1.-U2*(1.+RV)-U2*OMR*F)*XR2
      X6=.5*U2*OMR*(F*(U2*T1-SPP2-1./RK2)+1./RK2)*XR2/RK2
      EPH=-(X1-X2+X3-X4+X5+X6)*ECON
      RETURN
      END

max_stars_repo_name: ViperScience/ViperNEC

max_issues_repo_name: ViperScience/ViperNEC

max_forks_repo_name: ViperScience/ViperNEC


text:
      SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X)
C
C **  FROM NUMERICAL RECIPES
C CHANGE RECORD
C
      DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP)
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::TMP
      IF(.NOT.ALLOCATED(TMP))ALLOCATE(TMP(N))
C
      DO 12 J=1,N
        S=0.
        IF(W(J).NE.0.)THEN
          DO 11 I=1,M
            S=S+U(I,J)*B(I)
   11     CONTINUE
          S=S/W(J)
        ENDIF
        TMP(J)=S
   12 CONTINUE
      DO 14 J=1,N
        S=0.
        DO 13 JJ=1,N
          S=S+V(J,JJ)*TMP(JJ)
   13   CONTINUE
        X(J)=S
   14 CONTINUE
      DEALLOCATE(TMP)
      RETURN
      END


text:
      subroutine gen4mdkrad(r,p,pswt,*)
c--- this routine is an extension of gen4 to include the decay
c--- of one of the heavy particles, with radiation included in the decay

      implicit none
      include 'constants.f'
      include 'mxdim.f'
      include 'masses.f'
      include 'zerowidth.f'
      include 'limits.f'
      include 'phasemin.f'
      include 'breit.f'
      include 'x1x2.f'
      double precision r(mxdim)
      double precision p(mxpart,4),pswt,smin
      double precision p1(4),p2(4),p12(4),p8(4),p34568(4),
     & p7(4),p3458(4),p345(4),p6(4),wt12,wt34568,wt3458
      double precision p34(4),p5(4),p3(4),p4(4),wt345,wt34
      double precision mtbsq
      integer nu
      double precision xjac,p1ext(4),p2ext(4),wt0
      double precision tau,x1mx2,surd,lntaum
      common/pext/p1ext,p2ext
      parameter(wt0=1d0/twopi**2)
!$omp threadprivate(/pext/)

c--- this part is taken from gen4
      lntaum=dlog(taumin)
      tau=dexp(lntaum*(one-r(14)))
      xjac=-lntaum*tau

      x1mx2=two*r(15)-one
      surd=dsqrt(x1mx2**2+four*tau)

      xx(1)=half*(+x1mx2+surd)
      xx(2)=half*(-x1mx2+surd)

      pswt=xjac*two/surd

      if   ((xx(1) .gt. 1d0)  .or. (xx(2) .gt. 1d0)
     & .or. (xx(1) .lt. xmin) .or. (xx(2) .lt. xmin)) return 1

      do nu=1,4
      p1(nu)=xx(1)*p1ext(nu)
      p2(nu)=xx(2)*p2ext(nu)
      p12(nu)=-p1(nu)-p2(nu)
      enddo

c--- these must be set this way for this part
      mass3=mb ! check this
      n3=0

c--- taken from phase4
      mtbsq=(mt+mb)**2
      call phi1_2m(0d0,r(1),r(2),r(3),mtbsq,p12,p7,p34568,wt12,*99)
      call phi3m(r(4),r(5),p34568,p3458,p6,mt,mb,wt34568,*99)
c      r(6)=1d0-1d-4*r(6) ! debug - soft p8
      call phi1_2m(0d0,r(6),r(7),r(8),mb,p3458,p8,p345,wt3458,*99)
      pswt=pswt*wt0*wt12*wt34568*wt3458
c--- alternative
c      r(2)=r(2)*1d-5 ! DEBUG: s78 small
c      call phi1_2(r(1),r(2),r(3),r(4),p12,p3456,p78,wt12,*99)
c      call phi3m0(r(5),r(6),p78,p7,p8,wt78,*99)
c      call phi3m(r(7),r(8),p3456,p345,p6,mt,mb,wt3456,*99)
c      pswt=pswt*wt0*wt12*wt78*wt3456

c--- these must be set this way for this part
      mass3=wmass
      n3=1

c--- set up minimum invariant mass for the W in the top decay
      if (zerowidth) then
        smin=wmass**2
      else
        smin=wsqmin
      endif

c--- decay top -> b W
      call phi1_2m(mb,r(9),r(10),r(11),smin,p345,p5,p34,wt345,*99)
c--- NOTE: for now, generate massless b in top quark decay
c      call phi1_2m(zip,r(9),r(10),r(11),smin,p345,p5,p34,wt345,*99)
c--- decay W -> e n
      call phi3m0(r(12),r(13),p34,p3,p4,wt34,*99)

c--- compute new weight
      pswt=pswt/twopi**2*wt345*wt34*pi*mt*twidth

      do nu=1,4
      p(1,nu)=p1(nu)
      p(2,nu)=p2(nu)
      p(3,nu)=p3(nu)
      p(4,nu)=p4(nu)
      p(5,nu)=p5(nu)
      p(6,nu)=p6(nu)
      p(7,nu)=p7(nu)
      p(8,nu)=p8(nu)
      enddo

      return

   99 continue
      return 1

      end


text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module blip
! ------------------------------------------------------------------------------
! Code area 11: Basis functions
! ------------------------------------------------------------------------------

!!****h* Conquest/blip *
!!  NAME
!!   blip_module
!!  PURPOSE
!!   Stores various variables relating to blips
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   16/05/01
!!  MODIFICATION HISTORY
!!   28/05/2001 dave
!!    Added ROBODoc header
!!   18/03/2002 dave
!!    Added RCS id and log tags in header and static RCS id for object file
!!   04/07/2002 mike
!!   Added arrays region_single and region_double, defined as follows:
!!    the location of a blip function on the blip-grid is specified
!!    by the integer triplet (nx,ny,nz), with (0,0,0) at the centre
!!    of the support region. All the integers are in the range
!!  -bliparraysize <= nx,ny,nz <= bliparraysize. Now fix the
!!    magnitude of one of the integers to be na, with 0 <= na <= bliparraysize.
!!    Then the maximum magnitude of either of the other two integers
!!    is region_single(na). With one of the integers fixed at na
!!    with 0 <= na <= bliparraysize, and another of the integers fixed
!!    at nb, with 0 <= nb <= region_single(na), the maximum magnitude
!!    of the third integer is given by region_double(na,nb).
!!   11:30, 12/11/2004 dave
!!    Changed to use proper variables from maxima_module
!!   2006/07/12 16:56 dave
!!    Changing variables to dynamically allocatable
!!   2006/07/20 08:33 dave
!!    Removed use common, tidied - all variables now dynamically allocated, basis functions
!!    stored in derived types in support_spec_format
!!   2006/10/20 14:01 dave
!!    Included make_pre and set_blip_index
!!   2006/11/01 17:04 dave
!!    Added alpha, beta and gauss2blip
!!   2008/02/04 08:28 dave
!!    Changed for output to file not stdout
!!   2008/05/28 ast
!!    Added timers
!!   2011/11/15 07:59 dave
!!    Updating the derived type for blips
!!   2014/09/15 18:30 lat
!!    fixed call start/stop_timer to timer_module (not timer_stdlocks_module !)
!!  SOURCE
!!
module blip

  use datatypes
  use global_module,          only: io_lun
  use timer_module,           only: start_timer, stop_timer
  use timer_stdclocks_module, only: tmr_std_basis, tmr_std_allocation

  implicit none

  save

  real(double) :: alpha, beta

  type blip_data
     integer, pointer, dimension(:,:) :: blip_location ! dim( 3, n_blips)
     integer, pointer, dimension(:,:,:) :: blip_number ! dim(-bliparraysize:bliparraysize x3 )
     integer :: BlipArraySize
     integer :: NBlipsRegion
     integer :: OneArraySize, FullArraySize
     integer :: Extent
     real(double) :: SupportGridSpacing, RegionRadius
     real(double) :: BlipWidth
     real(double) :: FourOnBlipWidth
     ! These last two are temporary and will be removed once G2B is working
     integer, pointer, dimension(:) :: region_single   ! (0:bliparraysize)
     integer, pointer, dimension(:,:) :: region_double ! (0:bliparraysize,0:bliparraysize)
  end type blip_data

  type blip_precon
     integer :: size
     real(double), pointer, dimension(:,:) :: coeffs
  end type blip_precon

  type(blip_data), allocatable, dimension(:) :: blip_info
  type(blip_precon), allocatable, dimension(:) :: PreCond !(n_blips_region, n_blips_region )
  character(len=10) :: init_blip_flag

  integer :: blip_FFT_size, blip_FFT_off
  !integer, parameter :: nfft = 30
  !integer, dimension(nfft), parameter :: blip_FFT_grids = (/ &
  !     2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81/)
  integer, parameter :: nfft = 9
  integer, dimension(nfft), parameter :: blip_FFT_grids = (/ &
       3, 5, 9, 15, 25, 27, 45, 75, 81/)
!!***

contains

!!****f* blip/make_pre *
!!
!!  NAME
!!   make_pre
!!  USAGE
!!
!!  PURPOSE
!!   Calculates the preconditioning matrix for length-scale
!!    ill conditioning
!!  INPUTS
!!
!!
!!  USES
!!   datatypes, numbers, global_module, blip, dimens, GenComms, GenBlas
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   27/2/98
!!  MODIFICATION HISTORY
!!   10/05/01 DRB
!!    Converted to ROBODoc format and F90
!!   28/05/2001 dave
!!    Stripped subroutine calls
!!   08/06/2001 dave
!!    Added RCS Id and Log tags and GenComms for my_barrier
!!    Misc tidying
!!   11/06/2001 dave
!!    Changed dpot calls to generic ones through GenBlas
!!   25/06/2002 drb
!!    Tweaked headers, added USES line
!!   15:40, 2003/06/09 dave
!!    Added flag to switch preconditioning on/off
!!   2006/09/13 08:06 dave
!!    Changed to use species-dependent parameters from blip
!!   2006/10/20 13:57 dave
!!    Added memory registration and fixed small bug (stat not passed to allocate)
!!   2006/10/20 14:00 dave
!!    Included into blip
!!   2008/05/28 ast
!!    Added timers
!!  SOURCE
!!
  subroutine make_pre(inode, ionode)

    use datatypes
    use numbers
    use global_module, ONLY: iprint_basis, area_basis, flag_precondition_blips
    use GenComms, ONLY: my_barrier, cq_abort
    use GenBlas, ONLY: potrf, potri
    use species_module, ONLY: n_species
    use memory_module, ONLY: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables
    integer :: inode, ionode

    ! Local variables
    integer :: i, N, n_blip1, n_blip2, nx, n_blip3,  &
         nx1, nx2, ny, ny1, ny2, nz, nz1, nz2

    real(double) :: p0(0:3),y0(0:3),p2(0:3),asq, sum
    real(double), allocatable, dimension(:,:)  :: Kmat
    !real(double) :: Kmat(n_blips_region, n_blips_region)
    real(double), parameter :: k0=1.3229425_double
    integer :: n_blips,info,stat, spec

    call start_timer(tmr_std_basis)
    if(flag_precondition_blips) then
       stat=0
       call start_timer(tmr_std_allocation)
       allocate(PreCond(n_species),STAT=stat)
       if(stat/=0) call cq_abort('ERROR allocating PreCond: ',n_species)
       call stop_timer(tmr_std_allocation)
       do spec = 1,n_species
          call start_timer(tmr_std_allocation)
          allocate(PreCond(spec)%coeffs(blip_info(spec)%NBlipsRegion,blip_info(spec)%NBlipsRegion), &
               Kmat(blip_info(spec)%NBlipsRegion,blip_info(spec)%NBlipsRegion),STAT=stat)
          if(stat/=0) call cq_abort('ERROR allocating PreCond and Kmat: ',blip_info(spec)%NBlipsRegion)
          call reg_alloc_mem(area_basis,2*blip_info(spec)%NBlipsRegion*blip_info(spec)%NBlipsRegion,type_dbl)
          call stop_timer(tmr_std_allocation)
          PreCond(spec)%size = blip_info(spec)%NBlipsRegion
          PreCond(spec)%coeffs = zero
          Kmat = zero
          if(inode==ionode.AND.iprint_basis>0) write(io_lun,fmt='(6x,"Applying preconditioning")')
          asq = one/(k0*k0*blip_info(spec)%SupportGridSpacing*blip_info(spec)%SupportGridSpacing)
          ! Set up the coefficients for the blip multiples
          p0(0)=151.0_double/140.0_double
          p0(1)=1191.0_double/2240.0_double
          p0(2)=3.0_double/56.0_double
          p0(3)=1.0_double/2240.0_double
          y0(0)=3.0_double/2.0_double
          y0(1)=-9.0_double/32.0_double
          y0(2)=-9.0_double/20.0_double
          y0(3)=-3.0_double/160.0_double
          p2(0)=y0(0)/p0(0)
          p2(1)=y0(1)/p0(1)
          p2(2)=y0(2)/p0(2)
          p2(3)=y0(3)/p0(3)
          ! Build the approximate Hessian
          call my_barrier()
          do n_blip1 = 1, blip_info(spec)%NBlipsRegion
             do n_blip2 = 1, blip_info(spec)%NBlipsRegion

                nx1 = blip_info(spec)%blip_location( 1, n_blip1 )
                ny1 = blip_info(spec)%blip_location( 2, n_blip1 )
                nz1 = blip_info(spec)%blip_location( 3, n_blip1 )

                nx2 = blip_info(spec)%blip_location( 1, n_blip2 )
                ny2 = blip_info(spec)%blip_location( 2, n_blip2 )
                nz2 = blip_info(spec)%blip_location( 3, n_blip2 )

                nx = iabs(nx1-nx2)
                ny = iabs(ny1-ny2)
                nz = iabs(nz1-nz2)
                if((nx.le.3).AND.(ny.LE.3).AND.(nz.LE.3)) then
                   Kmat(n_blip1,n_blip2) = p0(nx)*p0(ny)*p0(nz)* &
                        (one+asq*(p2(nx)+p2(ny)+p2(nz)))
                else
                   Kmat(n_blip1,n_blip2) = zero
                endif
             end do
          end do
          ! Invert
          ! Cholesky decomposition
          call my_barrier()
          info = 0
          do n_blip1 = 1,blip_info(spec)%NBlipsRegion
             do n_blip2 = 1,blip_info(spec)%NBlipsRegion
                PreCond(spec)%coeffs(n_blip1,n_blip2) = Kmat(n_blip1,n_blip2)
             enddo
          enddo
          call potrf('U',blip_info(spec)%NBlipsRegion,Kmat,blip_info(spec)%NBlipsRegion,info)
          if(info/=0) call cq_abort("Problem with potrf in preconditioning: ",info)
          ! Inversion
          call my_barrier()
          info = 0
          call potri('U',blip_info(spec)%NBlipsRegion,Kmat,blip_info(spec)%NBlipsRegion,info)
          if(info/=0) call cq_abort("Problem with potri in preconditioning: ",info)
          ! Now make the upper triangle a full matrix
          call my_barrier()
          do n_blip1 = 1, blip_info(spec)%NBlipsRegion
             do n_blip2 = n_blip1, blip_info(spec)%NBlipsRegion
                Kmat(n_blip2,n_blip1) = Kmat(n_blip1,n_blip2)
             enddo
          enddo
          do n_blip1 = 1, blip_info(spec)%NBlipsRegion
             do n_blip2 = 1, blip_info(spec)%NBlipsRegion
                PreCond(spec)%coeffs(n_blip1,n_blip2) = Kmat(n_blip1,n_blip2)
             enddo
          enddo
          deallocate(Kmat)
          call reg_dealloc_mem(area_basis,blip_info(spec)%NBlipsRegion*blip_info(spec)%NBlipsRegion,type_dbl)
       end do
    else
       if(inode==ionode.AND.iprint_basis>0) &
            write(io_lun,fmt='(6x,"Not applying preconditioning")')
       !do n_blip1 = 1, blip_info(spec)%NBlipsRegion
       !   do n_blip2 = 1, blip_info(spec)%NBlipsRegion
       !      if(n_blip1.EQ.n_blip2) then
       !         PreCond(spec)%coeffs(n_blip1,n_blip2) = one
       !      else
       !         PreCond(spec)%coeffs(n_blip1,n_blip2) = zero
       !      endif
       !   enddo
       !enddo
    endif
    call stop_timer(tmr_std_basis)
    return
  end subroutine make_pre
!!***

!!****f* blip/set_blip_index *
!!
!!  NAME
!!   set_blip_index
!!  USAGE
!!
!!  PURPOSE
!!   Creates variables and index arrays characterising the set of
!!   blip functions wholly contained within the support region.
!!   Specifically:
!!   bliparraysize: the smallest integer such that for all blip
!!     functions wholly contained within the support region, these
!!     blip functions being centred on blip-grid sites (nx,ny,nz),
!!     the absolute values of nx, ny and nz never exceed bliparraysize.
!!   n_blips_region: total number of blip functions wholly contained
!!     within the support region.
!!   blip_number(nx,ny,nz): gives the sequence number of blip function
!!     at site (nx,ny,nz) in the enumeration in which nx varies
!!     most rapidly and nz least rapidly.
!!   blip_location(i,n): for blip function having sequence number n,
!    gives its blip-grid site (nx,ny,nz), with components nx, ny and
!!     nz corresponding to index i = 1, 2, 3 respectively.
!!   region_single(na), region_double(na,nb): defined as follows:
!!    the location of a blip function on the blip-grid is specified
!!    by the integer triplet (nx,ny,nz), with (0,0,0) at the centre
!!    of the support region. All the integers are in the range
!!  -bliparraysize <= nx,ny,nz <= bliparraysize. Now fix the
!!    magnitude of one of the integers to be na, with 0 <= na <= bliparraysize.
!!    Then the maximum magnitude of either of the other two integers
!!    is region_single(na). With one of the integers fixed at na
!!    with 0 <= na <= bliparraysize, and another of the integers fixed
!!    at nb, with 0 <= nb <= region_single(na), the maximum magnitude
!!    of the third integer is given by region_double(na,nb).
!!  INPUTS
!!   inode: id number of node running executable
!!   ionode: id number of input/output node
!!
!!  USES
!!
!!  AUTHOR
!!   E.H.Hernandez
!!  CREATION DATE
!!   14/04/95
!!  MODIFICATION HISTORY
!!   16/05/2001 dave
!!    Converted to F90, added ROBODoc header and reduced number
!!    of passed variables
!!   28/05/2001 dave
!!    Removed r_H as passed variables and included from dimens
!!   20/06/2001 dave
!!    Added cq_abort and RCS Id and Log tags
!!   25/06/2002 mike
!!    Code simplified, protection against array-bound overflow
!!    improved, comments added
!!   29/06/2002 mike
!!    Code added to create arrays region_single and region_double
!!   2006/10/20 14:00 dave
!!    Included into blip
!!   2008/05/28 ast
!!    Added timers
!!   2011/11/17 10:38 dave
!!    Changes for new blip data (removing allocation)
!!   2012/01/18 16:56 dave
!!    Initialisation for blip coefficient transfer
!!   2012/04/26 16:18 dave
!!    More changes for analytic blips
!!   2016/07/29 18:30 nakata
!!    Renamed supports_on_atom -> blips_on_atom
!!  SOURCE
!!
  subroutine set_blip_index(inode,ionode)

    use datatypes
    use GenComms, ONLY: cq_abort
    use global_module, ONLY: iprint_init, area_basis
    use numbers
    use support_spec_format, ONLY : blips_on_atom, support_gradient, support_elec_gradient, &
         coefficient_array,grad_coeff_array, elec_grad_coeff_array, coeff_array_size,&
         allocate_supp_coeff_array, associate_supp_coeff_array
    use species_module, ONLY: nsf_species, n_species
    use primary_module, ONLY: bundle
    use dimens, ONLY: RadiusSupport
    use memory_module, ONLY: reg_alloc_mem, type_dbl, type_int
    use fft_module, ONLY: set_SF_fft
    use matrix_data, ONLY: blip_trans

    implicit none

    ! Variables passed in argument list
    integer, intent(in) :: inode, ionode

    ! Local variables
    integer :: i, n_blips, na, nb, nx, ny, nz, stat, size, this_nsf, spec, maxsize, ni, np
    real(double) :: a2, b2, dx, dy, dz, r2_over_b2

    call start_timer(tmr_std_basis)
    if((inode == ionode).and.(iprint_init >= 2)) then
       write(unit=io_lun,fmt='(//10x,65("*")/25x,"REPORT FROM SET_BLIP_INDEX"/&
            &10x,65("*"))')
    end if

    maxsize=0
    do spec = 1, n_species
       ! It is required that the full width of the B-spline blip-function
       ! be exactly four times the blip-grid spacing. This condition
       ! is checked here, and the calculation stops if the condition
       ! is not satisfied:
       !blip_info(spec)%BlipWidth = four*blip_info(spec)%SupportGridSpacing
       if(inode==ionode.AND.iprint_init>2) write(io_lun,fmt='(10x,"Blip width: ",f20.12)') blip_info(spec)%BlipWidth
       blip_info(spec)%FourOnBlipWidth = four/blip_info(spec)%BlipWidth
       !if(abs(blip_width - four*support_grid_spacing) > RD_ERR) then
       !   call cq_abort('set_blip_index: blip width must be exactly four &
       !        &times blip-grid spacing')
       !end if

       ! check that the relation between region radius and blip-grid
       ! spacing is such that at least one blip function is
       ! wholly contained within the region (RadiusSupport(spec) is support-region radius,
       ! support_grid_spacing is blip-grid spacing):
       r2_over_b2 = (RadiusSupport(spec)*RadiusSupport(spec))/ &
            (blip_info(spec)%SupportGridSpacing*blip_info(spec)%SupportGridSpacing)
       if(r2_over_b2 < 12.0_double) then
          call cq_abort('set_blip_grid: support region too small to hold&
               & one blip function')
       end if

       ! For the meaning of the integer bliparraysize, see header. Here,
       ! bliparraysize is calculated by formula, printed, and
       ! tested against the relevant array bound:
       blip_info(spec)%BlipArraySize = sqrt(r2_over_b2 - 8.0_double) - 2.0_double
       if((inode == ionode).and.(iprint_init >= 2)) then
          write(unit=io_lun,fmt='(/10x," set_blip_index: bliparraysize:",i5/10x,&
               &"(This is the smallest integer such that for all blip functions",/10x,&
               &" wholly contained within the support region, these blip",/10x,&
               &" functions being centred on sites (nx,ny,nz), the absolute",/10x,&
               &" values of nx, ny and nz never exceed this integer)")') &
               &blip_info(spec)%BlipArraySize
       end if
       blip_info(spec)%OneArraySize = 1 + 2*(blip_info(spec)%BlipArraySize+4)
       blip_info(spec)%FullArraySize = blip_info(spec)%OneArraySize*blip_info(spec)%OneArraySize*&
            blip_info(spec)%OneArraySize
       if(blip_info(spec)%BlipArraySize>maxsize) maxsize = blip_info(spec)%BlipArraySize
       !if(bliparraysize > MAXBAS) then
       !   call cq_abort('set_blip_index: MAXBAS too small',bliparraysize,MAXBAS)
       !end if
       call start_timer(tmr_std_allocation)
       allocate(blip_info(spec)%region_single(0:blip_info(spec)%BlipArraySize),&
            blip_info(spec)%region_double(0:blip_info(spec)%BlipArraySize,0:blip_info(spec)%BlipArraySize),&
            STAT=stat)
       if(stat/=0) call cq_abort('Error allocating region_single and region_double: ',&
            blip_info(spec)%BlipArraySize)
       call reg_alloc_mem(area_basis,(blip_info(spec)%BlipArraySize+1)*(blip_info(spec)%BlipArraySize+2),type_int)
       call stop_timer(tmr_std_allocation)
       ! Calculate the index arrays region_single, region_double
       do na = 0, blip_info(spec)%BlipArraySize
          a2 = (na+2)*(na+2)
          blip_info(spec)%region_single(na) = sqrt(r2_over_b2 - a2 - 4.0_double) - 2.0_double
          do nb = 0, blip_info(spec)%region_single(na)
             b2 = (nb+2)*(nb+2) + a2
             blip_info(spec)%region_double(na,nb) = sqrt(r2_over_b2 - b2) - 2.0_double
          end do
       end do

       n_blips = 0
       do nz = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
          dz = 2 + abs(nz)
          do ny = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
             dy = 2 + abs(ny)
             do nx = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
                dx = 2 + abs(nx)
                if (dx*dx + dy*dy + dz*dz <= r2_over_b2) n_blips = n_blips + 1
             end do
          end do
       end do
       ! For meaning of the array blip_info(nx,ny,nz), see description
       ! at head of routine. The entire array
       ! is initialised to zero here to maintain compatibility with
       ! the original version of the code.

       call start_timer(tmr_std_allocation)
       allocate(blip_info(spec)%blip_number(-blip_info(spec)%BlipArraySize:blip_info(spec)%BlipArraySize,&
            -blip_info(spec)%BlipArraySize:blip_info(spec)%BlipArraySize,&
            -blip_info(spec)%BlipArraySize:blip_info(spec)%BlipArraySize),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating blip_number: ',blip_info(spec)%BlipArraySize)
       call reg_alloc_mem(area_basis,(2*blip_info(spec)%BlipArraySize+1)*(2*blip_info(spec)%BlipArraySize+1)*&
            (2*blip_info(spec)%BlipArraySize+1),type_int)
       call stop_timer(tmr_std_allocation)
       blip_info(spec)%blip_number = 0

       call start_timer(tmr_std_allocation)
       allocate(blip_info(spec)%blip_location(3,n_blips),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating blip_number and blip_location: ',n_blips)
       call reg_alloc_mem(area_basis,3*n_blips,type_int)
       call stop_timer(tmr_std_allocation)
       blip_info(spec)%NBlipsRegion = n_blips
       n_blips = 0
       do nz = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
          dz = 2 + abs(nz)
          do ny = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
             dy = 2 + abs(ny)
             do nx = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
                dx = 2 + abs(nx)
                if (dx*dx + dy*dy + dz*dz <= r2_over_b2) then
                   n_blips = n_blips + 1
                   blip_info(spec)%blip_number(nx,ny,nz) = n_blips
                   blip_info(spec)%blip_location(1,n_blips) = nx
                   blip_info(spec)%blip_location(2,n_blips) = ny
                   blip_info(spec)%blip_location(3,n_blips) = nz
                end if
             end do
          end do
       end do
       if((inode == ionode).and.(iprint_init >= 2)) then
          write(unit=io_lun,fmt='(/10x," set_blip_index: total number of blip &
               &functions wholly contained within support region:",i5)') &
               &blip_info(spec)%NBlipsRegion
       end if

       if((inode == ionode).and.(iprint_init > 2)) then
          write(unit=io_lun,fmt='(/10x," set_blip_index: index array region_single:"/)')
          do na = 0, blip_info(spec)%BlipArraySize
             write(unit=io_lun,fmt='(10x,i5,3x,i5)') na, blip_info(spec)%region_single(na)
          end do
          write(unit=io_lun,fmt='(/10x," set_blip_index: index array region_double:"/)')
          do na = 0, blip_info(spec)%BlipArraySize
             do nb = 0, blip_info(spec)%region_single(na)
                write(unit=io_lun,fmt='(10x,2i5,3x,i5)') na, nb, blip_info(spec)%region_double(na,nb)
             end do
          end do
          write(unit=io_lun,fmt='(/10x," set_blip_index: index array blip_number:"/)')
          do nz = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
             do ny = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
                do nx = -blip_info(spec)%BlipArraySize, blip_info(spec)%BlipArraySize
                   write(unit=io_lun,fmt='(10x,3i5,3x,i5)') nx, ny, nz, &
                        &blip_info(spec)%blip_number(nx,ny,nz)
                end do
             end do
          end do
          write(unit=io_lun,fmt='(/10x," set_blip_index: index array blip_location:"/)')
          do i = 1, blip_info(spec)%NBlipsRegion
             write(unit=io_lun,fmt='(10x,i5,3x,3i5)') i, blip_info(spec)%blip_location(1,i), &
                  &blip_info(spec)%blip_location(2,i), blip_info(spec)%blip_location(3,i)
          end do
       end if
    end do ! Loop over species
    blip_FFT_size = 2*(maxsize+4)+1
    do i=1,nfft
       if(blip_FFT_size==blip_FFT_grids(i)) then
          exit
       else if(blip_FFT_size<blip_FFT_grids(i)) then
          blip_FFT_size=blip_FFT_grids(i)
          exit
       end if
    end do
    blip_FFT_off = (blip_FFT_size-1)/2
    call set_SF_fft(blip_FFT_size)!2*(maxsize+4)+1)
    ! Set up data storage for atoms on this processor
    size = 0
    call start_timer(tmr_std_allocation)
    allocate(blips_on_atom(bundle%n_prim),support_gradient(bundle%n_prim),support_elec_gradient(bundle%n_prim))
    allocate(blip_trans%partst(bundle%groups_on_node),blip_trans%partlen(bundle%groups_on_node))
    call stop_timer(tmr_std_allocation)
    i=0
    !do i=1,bundle%n_prim
    blip_trans%partst = 0
    blip_trans%partlen = 0
    blip_trans%partst(1) = 1
    do np=1,bundle%groups_on_node
       blip_trans%partlen(np) = 0
       if(bundle%nm_nodgroup(np) > 0) then
          do ni=1,bundle%nm_nodgroup(np)
             i=i+1
             ! Check on species
             spec = bundle%species(i)
             this_nsf = nsf_species(spec)
             size = size + this_nsf*blip_info(spec)%NBlipsRegion
             blip_trans%partlen(np)=blip_trans%partlen(np)+ this_nsf*blip_info(spec)%NBlipsRegion
             blips_on_atom(i)%nsuppfuncs = this_nsf
             call start_timer(tmr_std_allocation)
             allocate(blips_on_atom(i)%supp_func(this_nsf),support_gradient(i)%supp_func(this_nsf),&
                  support_elec_gradient(i)%supp_func(this_nsf))
             call stop_timer(tmr_std_allocation)
             blips_on_atom(i)%supp_func(:)%ncoeffs = blip_info(spec)%NBlipsRegion
             support_gradient(i)%nsuppfuncs = this_nsf
             support_gradient(i)%supp_func(:)%ncoeffs = blip_info(spec)%NBlipsRegion
             support_elec_gradient(i)%nsuppfuncs = this_nsf
             support_elec_gradient(i)%supp_func(:)%ncoeffs = blip_info(spec)%NBlipsRegion
          end do
       end if
       if(np>1) then
          blip_trans%partst(np) = blip_trans%partst(np-1)+blip_trans%partlen(np-1)
       end if
    end do
    coeff_array_size = size
    call allocate_supp_coeff_array(size)
    call associate_supp_coeff_array(blips_on_atom,bundle%n_prim,coefficient_array,size)
    call associate_supp_coeff_array(support_gradient,bundle%n_prim,grad_coeff_array,size)
    call associate_supp_coeff_array(support_elec_gradient,bundle%n_prim,elec_grad_coeff_array,size)
    call setup_blip_transfer
    call stop_timer(tmr_std_basis)

    return

  end subroutine set_blip_index
!!***

!!****f* blip/setup_blip_transfer *
!!
!!  NAME
!!   setup_blip_transfer
!!  USAGE
!!
!!  PURPOSE
!!   Creates arrays needed for transfer of blip coefficients for analytic blip integrals
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   D. R. Bowler
!!  CREATION DATE
!!   2012/01/12
!!  MODIFICATION HISTORY
!!
!!  SOURCE
!!
  subroutine setup_blip_transfer

    use matrix_data, ONLY: halo, Srange, blip_trans
    use mult_module, ONLY: mult, LHL_SL_S
    use group_module, ONLY: parts
    use global_module, ONLY: species_glob, id_glob
    use species_module, ONLY: nsf_species

    implicit none

    ! Local variables
    integer :: i, j, kpart, jpart, partlen, part, iproc, thisproc, gcspartlab, ni, specj, jseq, j_in_halo, &
         maxps, maxpr

    blip_trans%nproc = mult(LHL_SL_S)%comms%inode
    if(blip_trans%nproc>0) then
       maxps = 0
       do i=1,blip_trans%nproc
          if(mult(LHL_SL_S)%comms%np_send(i)>maxps) maxps = mult(LHL_SL_S)%comms%np_send(i)
       end do
       allocate(blip_trans%np_send(blip_trans%nproc),blip_trans%pl_send(maxps,blip_trans%nproc+1), &
            blip_trans%neigh_pl(halo(Srange)%np_in_halo),blip_trans%ncomm(blip_trans%nproc), &
            blip_trans%len_recv(halo(Srange)%np_in_halo))
       blip_trans%npart_send = 0
       do i=1,blip_trans%nproc
          blip_trans%np_send(i) = mult(LHL_SL_S)%comms%np_send(i)
          blip_trans%npart_send = blip_trans%npart_send + blip_trans%np_send(i)
          blip_trans%ncomm(i) = mult(LHL_SL_S)%comms%ncomm(i)
          do kpart=1,blip_trans%np_send(i)
             blip_trans%pl_send(kpart,i) = mult(LHL_SL_S)%comms%pl_send(kpart,i)
             ! The length and start are blip_trans%partst and blip_trans%partlen, accessed by pl_send(kpart,i)
             ! and defined in blip_module
          end do
       end do
    else
       allocate(blip_trans%neigh_pl(halo(Srange)%np_in_halo), &
            blip_trans%len_recv(halo(Srange)%np_in_halo), blip_trans%pl_send(halo(Srange)%np_in_halo,1))
       do i=1,halo(Srange)%np_in_halo
          blip_trans%pl_send(i,1) = parts%i_cc2seq(halo(Srange)%lab_hcell(i))
       end do
    end if
    do jpart=1,halo(Srange)%np_in_halo
       partlen = 0
       blip_trans%neigh_pl(jpart) =  mult(LHL_SL_S)%comms%neigh_node_list(jpart)
       gcspartlab = halo(Srange)%i_hbeg(halo(Srange)%lab_hcover(jpart))
       do j=1,halo(Srange)%nh_part(jpart)
          j_in_halo = halo(Srange)%j_beg(jpart) + j-1
          jseq = halo(Srange)%j_seq(j_in_halo)
          specj = species_glob( id_glob( parts%icell_beg(halo(Srange)%lab_hcell(jpart))+jseq-1) )
          !partlen = partlen + halo(Srange)%ndimj(j_in_halo)*blip_info(specj)%NBlipsRegion
       end do
       ! Store partlen
       do j=1,parts%nm_group(halo(Srange)%lab_hcell(jpart))
          specj = species_glob( id_glob( parts%icell_beg(halo(Srange)%lab_hcell(jpart))+j-1) )
          partlen = partlen + nsf_species(specj)*blip_info(specj)%NBlipsRegion
       end do
       blip_trans%len_recv(jpart) = partlen
    end do
    return
  end subroutine setup_blip_transfer
!!***

!!****f* blip/gauss2blip *
!!
!!  NAME
!!   gauss2blip
!!  USAGE
!!
!!  PURPOSE
!!   Sets blip coefficients as values of a Gaussian function
!!   at the distance of each blip from the atom position. Can
!!   only be used if there are four supports on each atom,
!!   in which case the support have the form of teh s-function
!!   and three p-functions. The Gaussian exponents for s- and
!!   p-functions are alpha and beta (atomic units). This way
!!   of initiating blip coefficients is a relic of the very
!!   early history of Conquest, and its use is strongly
!!   discouraged, unless you know what you are doing. It
!!   continues to be included in Conquest mainly in order
!!   to maintain compatibility with earlier version.
!!
!!  INPUTS
!!
!!
!!  USES
!!
!!  AUTHOR
!!   E.H.Hernandez
!!  CREATION DATE
!!   11/2/95
!!  MODIFICATION HISTORY
!!   10/05/01 DRB
!!    Added ROBODoc header and updated to F90
!!   28/05/2001 dave
!!    Stripped subroutine call
!!   13/06/02 MJG
!!    Name of routine changed from set_support to gauss2blip
!!    and comments added. In addition, the number of
!!    supports on each atom is checked, and the code aborts
!!    if this number is not equal to 4. fortran coding tidied up.
!!   15:48, 04/02/2003 drb
!!    Changed output so that it's only seen for iprint>2 and changed format statement for grid spacing
!!   2006/09/13 07:49 dave
!!    Changed loop to be over primary set and added species dependence
!!   2006/11/01 17:06 dave
!!    Included in blip
!!   2008/03/03 18:41 dave
!!    Changed float to real (probably unnecessary)
!!   2008/05/28 ast
!!    Added timers
!!   2008/09/17 10:01 dave
!!    Added nsf=1 case
!!   2016/07/29 18:30 nakata
!!    Renamed supports_on_atom -> blips_on_atom
!!  SOURCE
!!
  subroutine gauss2blip

    use datatypes
    use numbers
    use global_module, ONLY: iprint_basis
    use GenComms, ONLY: cq_abort, inode, ionode
    use support_spec_format, ONLY: coefficient_array, blips_on_atom
    use primary_module, ONLY: bundle

    implicit none

    ! Local variables
    integer :: i, n_blip, nx, ny, nz, spec
    real(double) :: gauss_1, gauss_2, r2, x, y, z

    ! variables accessed by use:
    ! -------------------------------------------
    !   variable                  module
    ! -------------------------------------------
    !   blip_info(spec)%blip_location:  blip
    !   data_blip:                      blip
    !   iprint:                         global module
    !   blip_info(spec)%NBlipsRegion:             blip
    !   blip_info(spec)%SupportGridSpacing:       blip
    ! -------------------------------------------

    ! since the initialisation method provided in this routine
    ! can only be used if there are 4 support functions on
    ! every atom, this condition is checked. It is a fatal
    ! error if the condition is not satisfied, and the
    ! calculation is aborted.

    call start_timer(tmr_std_basis)
    if((inode == ionode).and.(iprint_basis >= 0)) then
       write(unit=io_lun,fmt='(/10x," gauss2blip: sbrt entered")')
    end if

    !if((inode == ionode) .AND. (iprint_basis >= 0)) then
    !   write(unit=io_lun,fmt='(10x," gauss2blip: support-grid spacing:",f20.12)') &
    !        &blip_info(spec)%SupportGridSpacing
    !end if

    ! data_blip is the array holding the blip coeffs for the support
    ! support functions of atoms in the primary set of my processor. T
    ! maintain consistency with earlier versions of the code, the
    ! entire array is initialised to zero before values are
    ! assigned to the blip coeffs.

    coefficient_array = zero

    ! loop over all atoms in my primary set, and over all blips in
    ! support region of each atom. At present, all species of atoms
    ! have the same support region and the same blip grid.
    do i = 1, bundle%n_prim
       spec = bundle%species(i)
       do n_blip = 1, blip_info(spec)%NBlipsRegion
          nx = blip_info(spec)%blip_location( 1, n_blip )
          ny = blip_info(spec)%blip_location( 2, n_blip )
          nz = blip_info(spec)%blip_location( 3, n_blip )
          x = real( nx, double ) * blip_info(spec)%SupportGridSpacing
          y = real( ny, double ) * blip_info(spec)%SupportGridSpacing
          z = real( nz, double ) * blip_info(spec)%SupportGridSpacing
          r2 = x * x + y * y + z * z
          gauss_1 = exp( - alpha * r2)
          gauss_2 = exp( - beta * r2)
          if(blips_on_atom(i)%nsuppfuncs>=4) then
             blips_on_atom(i)%supp_func(1)%coefficients(n_blip) = gauss_1
             blips_on_atom(i)%supp_func(2)%coefficients(n_blip) = x * gauss_2
             blips_on_atom(i)%supp_func(3)%coefficients(n_blip) = y * gauss_2
             blips_on_atom(i)%supp_func(4)%coefficients(n_blip) = z * gauss_2
             if(iprint_basis>2.AND.inode==ionode) then
                write(io_lun,fmt='(10x,"Blip Values: ",4i4,f20.12)') inode, i, 1, n_blip, &
                     blips_on_atom(i)%supp_func(1)%coefficients(n_blip)
                write(io_lun,fmt='(10x,"Blip Values: ",4i4,f20.12)') inode, i, 2, n_blip, &
                     blips_on_atom(i)%supp_func(2)%coefficients(n_blip)
                write(io_lun,fmt='(10x,"Blip Values: ",4i4,f20.12)') inode, i, 3, n_blip, &
                     blips_on_atom(i)%supp_func(3)%coefficients(n_blip)
                write(io_lun,fmt='(10x,"Blip Values: ",4i4,f20.12)') inode, i, 4, n_blip, &
                     blips_on_atom(i)%supp_func(4)%coefficients(n_blip)
             end if
          else if (blips_on_atom(i)%nsuppfuncs==1) then
             blips_on_atom(i)%supp_func(1)%coefficients(n_blip) = gauss_1
             if(iprint_basis>2.AND.inode==ionode) then
                write(io_lun,fmt='(10x,"Blip Values: ",4i4,f20.12)') inode, i, 1, n_blip, &
                     blips_on_atom(i)%supp_func(1)%coefficients(n_blip)
             end if
          else
             call cq_abort("With this version of Conquest, gaussian initialisation needs 1 or 4 support functions ",&
                  blips_on_atom(i)%nsuppfuncs)
          end if
       end do
    end do
    call stop_timer(tmr_std_basis)

    return
  end subroutine gauss2blip
!!***

end module blip

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
subroutine iterate()
  real :: x

  x = 100.0;
end subroutine iterate

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf