text:
c******************************************************
      SUBROUTINE EGAUSS_2(lidd,nlidd,iin,nmm,nmma,ff,dens,wden)
      IMPLICIT NONE
      include "tomo.h"
      include "linea.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      real*8 swgt(3),wgt,tresh(3),s,s1,sigmw(3),x1(3),wgta
      real*4 wden,ff(NRAZ+NRAZA*(NAZIPL-1)),dens(NRAZA)
      integer*4 i,i1,i2,i3,iin,ilo,iou,ish,ish1,j,k,nmm,nmma
      integer*4 ip,nlidd,lidd(NRAZ+NRAZA*(NAZIPL-1))
c END TYPE definitions ++++++++++++++++++++++++++++++++++
C******************************************************
c*** exp(-s**2/2/sigma**2) < 10**(-pow), tresh=COS(s)
C******************************************************
      wden=0.0
      do i=2,ianiz+1
      swgt(i)=0.0d0
cxx   sigmw(i)=((alsi(i,4)-alsi(i,3))*(azmax-dazi(iin))/azmax+alsi(i,3))/R
      sigmw(i)=alsi(i,3)/R
      tresh(i)=DCOS(sigmw(i)*DSQRT(2.d0*DBLE(alsi(i,5))*DLOG(10.0d0)))
      enddo
      iou=ioutra(iin)
      i1=icrpnta(iou)/1000000
      i3=icrpnta(iou)-i1*1000000
      i2=i3/1000
      i3=i3-i2*1000
      s1=dra(i1)*dra(i1)+dra(i2)*dra(i2)+dra(i3)*dra(i3)
      s=DSQRT(s1)
      x1(1)=dra(i1)/s
      x1(2)=dra(i2)/s
      x1(3)=dra(i3)/s
      do 20 i=1,nalla
      s1=0.0d0
      ilo=ilocra(i)
      do 10 j=1,3
   10 s1=s1+x1(j)*dpxa(j,i)
      if(s1.gt.1.d0)s1=1.d0
      do 40 k=1,ianiz
      if(s1.lt.tresh(k+1)) goto 40
      s=DACOS(s1)
      wgta=s*s/sigmw(k+1)/sigmw(k+1)/2.0d0
      wgt=0.0d0
      if(wgta.lt.30.0d0) then
      wgt=DEXP(wgta)
      wgt=dpxa(4,i)/wgt
      endif
      swgt(k+1)=swgt(k+1)+wgt
      if(ilo.ne.0) then
      ff(ilo+(2*k-1)*nmma+nmm)=wgt
      if(k.eq.1) wden=wden+dens(ilo)*wgt
      endif
   40 continue
   20 continue
      do 50 k=1,ianiz
      ish=(2*k-1)*nmma+nmm
      do 30 j=1,nmma
   30 ff(j+ish)=ff(j+ish)*alsi(k+1,1)/swgt(k+1)
   50 ff(iin+ish)=(ff(iin+ish)/alsi(k+1,1)-1.0)*alsi(k+1,1)
      wden=wden/SNGL(swgt(2))
      do 70 k=1,ianiz
      ish=(2*k-1)*nmma+nmm
      ish1=ish-nmma
      do 60 i=1,nmma
   60 ff(i+ish1)=ff(i+ish)
   70 continue
      do j=nmm+1,nmm+2*ianiz*nmma
        if(ff(j).ne.0.0) then
          nlidd=nlidd+1
          ip=(j-nmm-1)/nmma
          lidd(nlidd)=j*10+ip
        endif
      enddo
      return
      end
c******************************************************
      SUBROUTINE E2PNTS_2(e,ngr_out,grid)
      IMPLICIT NONE
      include "linea.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      real*8    grid(4,3),e(3)
      real*8    dmax,de
      integer*4 ic,icol,index,ishift,ival,ndp1,ndm1
      integer*4 ix,iy,iz,j,k,nplane,ngr_out(4),ne(3)
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      equivalence (ne(1),ix),(ne(2),iy),(ne(3),iz)
      dmax=DMAX1(DABS(e(1)),DABS(e(2)),DABS(e(3)))
      nplane=0
      do 30 j=1,3
      if(nplane.eq.0.and.DABS(e(j)).eq.dmax) then
      nplane=j
      de=e(j)/dmax
      if(de.lt.-.5) ne(j)=1
      if(de.gt.0.5)ne(j)=nda+1
      goto 30
      endif
      de=e(j)/dmax
      icol=nda
      do 10 k=2,nda
      if(de.lt.dra(k)) then
      icol=k-1
      goto 20
      endif
   10 continue
   20 ne(j)=icol
   30 continue
      if(nplane.eq.0) then
      write(*,*)'NO CELL FOUND'
      stop
      endif
      ic=0
c****** Get integer coordinates *************************
      ishift=1000000
      do 40 j=1,3
      if(j.ne.nplane) then
      if(ic.lt.1) then
      ic=ic+1
      ngr_out(ic)=1000000*ne(1)+1000*ne(2)+ne(3)
      ic=ic+1
      ngr_out(ic)=1000000*ne(1)+1000*ne(2)+ne(3)+ishift
      else
      ic=ic+1
      ngr_out(ic)=ngr_out(ic-2)+ishift
      ic=ic+1
      ngr_out(ic)=ngr_out(ic-2)+ishift
      endif
      endif
      ishift=ishift/1000
   40 continue
      do 50 j=1,4
      ival=ngr_out(j)
c****** Searching by analitic expression ******************
      ix=ival/1000000
      iz=ival-ix*1000000
      iy=iz/1000
      iz=iz-iy*1000
      ndp1=nda+1
      ndm1=nda-1
      if(ix.eq.1) then
      index=(iy-1)*ndp1+iz
      else if(ix.eq.ndp1) then
      index=ndp1*ndp1+4*ndm1*nda+(iy-1)*ndp1+iz
      else
      index=(ix-2)*4*nda+ndp1*ndp1
      if(iy.eq.1) then
      index=index+iz
      else if(iy.eq.ndp1) then
      index=index+nda+1+2*ndm1+iz
      else
      index=index+2*(iy-2)+nda+2+iz/nda
      endif
      endif
      ngr_out(j)=index
c****** END Searching by analitic expression ******************
      de=DSQRT(dra(ix)*dra(ix)+dra(iy)*dra(iy)+dra(iz)*dra(iz))
      grid(j,1)=dra(ix)/de
      grid(j,2)=dra(iy)/de
   50 grid(j,3)=dra(iz)/de
      return
      end
c***************************************************************
      SUBROUTINE TRIAW_2 (ngr_out,grid,e,wei,ierr)
      IMPLICIT NONE
      include "linea.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      real*8    grid(4,3),e(3),wei(4)
      real*8    x1(3),x2(3),x3(3),u(3),w(3),uwn(3),x0(3),dnorm
      real*8    dlam(3),du,dw,duw,deu,dew
      integer*4    ngr_out(4),i,j,ierr
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      ierr=0
      do 9 j=0,1
      do 1 i=1,4
    1 wei(i)=0.0d0
      do 2 i=1,3
      x1(i)=grid(1+j,i)
      x2(i)=grid(2+j,i)
      x3(i)=grid(3+j,i)
      u(i)=x2(i)-x1(i)
    2 w(i)=x3(i)-x1(i)
c***** Intersection vector e with plane within x1, x2 and x3 ****
      CALL VECT(u,w,uwn)
      CALL NORM(uwn,uwn)
      dnorm=(uwn(1)*x1(1)+uwn(2)*x1(2)+uwn(3)*x1(3))/
     *(uwn(1)*e(1)+uwn(2)*e(2)+uwn(3)*e(3))
      do 3 i=1,3
    3 x0(i)=e(i)*dnorm-x1(i)
      du=u(1)*u(1)+u(2)*u(2)+u(3)*u(3)
      dw=w(1)*w(1)+w(2)*w(2)+w(3)*w(3)
      duw=u(1)*w(1)+u(2)*w(2)+u(3)*w(3)
      dnorm=du*dw-duw*duw
      deu=x0(1)*u(1)+x0(2)*u(2)+x0(3)*u(3)
      dew=x0(1)*w(1)+x0(2)*w(2)+x0(3)*w(3)
      dlam(2)=(deu*dw-dew*duw)/dnorm
      dlam(3)=(dew*du-deu*duw)/dnorm
      dlam(1)=1.d0-dlam(2)-dlam(3)
      do 5 i=1,3
      if(ilocra(ngr_out(i+j)).eq.0) ierr=1
    5 wei(i+j)=dlam(i)
      return
    9 continue
      end
c***************************************************************
      SUBROUTINE TRIAI_2 (grid,e,fun,funr)
      IMPLICIT NONE
      include "linea.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      real*8    grid(4,3),e(3),fun(4),funr
      real*8    x1(3),x2(3),x3(3),u(3),w(3),uwn(3),x0(3),dnorm
      real*8    dlam(3),du,dw,duw,deu,dew,ngr_out(4)
      integer*4 i,j,ierr,isu
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      do 9 j=0,1
      do 2 i=1,3
      x1(i)=grid(1+j,i)
      x2(i)=grid(2+j,i)
      x3(i)=grid(3+j,i)
      u(i)=x2(i)-x1(i)
    2 w(i)=x3(i)-x1(i)
c***** Intersection vector e with plane within x1, x2 and x3 ****
      CALL VECT(u,w,uwn)
      CALL NORM(uwn,uwn)
      dnorm=(uwn(1)*x1(1)+uwn(2)*x1(2)+uwn(3)*x1(3))/
     *(uwn(1)*e(1)+uwn(2)*e(2)+uwn(3)*e(3))
      do 3 i=1,3
    3 x0(i)=e(i)*dnorm-x1(i)
      du=u(1)*u(1)+u(2)*u(2)+u(3)*u(3)
      dw=w(1)*w(1)+w(2)*w(2)+w(3)*w(3)
      duw=u(1)*w(1)+u(2)*w(2)+u(3)*w(3)
      dnorm=du*dw-duw*duw
      deu=x0(1)*u(1)+x0(2)*u(2)+x0(3)*u(3)
      dew=x0(1)*w(1)+x0(2)*w(2)+x0(3)*w(3)
      dlam(2)=(deu*dw-dew*duw)/dnorm
      dlam(3)=(dew*du-deu*duw)/dnorm
      dlam(1)=1.d0-dlam(2)-dlam(3)
      funr=0.0d0
      do 5 i=1,3
      isu=ngr_out(i+j)
      if(ilocra(isu).eq.0) ierr=1
    5 funr=funr+dlam(i)*fun(i+j)
      return
    9 continue
      end
c***************************************************************
      SUBROUTINE EINTPOL_2 (xk,yk,nmm,nmma,fsl,value,nval,ierr)
      IMPLICIT NONE
      include "linea.h"
c TYPE definitions ++++++++++++++++++++++++++++++++++++++
      integer*4 nmm,nmma,ierr,j,mo,nclo,nval,ngr_out(4)
      real*4    fsl(nmm),value(5)
      real*8    grid(4,3),w(4),e(3),scell,xk,yk
c END TYPE definitions ++++++++++++++++++++++++++++++++++
      ierr=0
      e(1)=DSIN(xk)*DCOS(yk)
      e(2)=DSIN(xk)*DSIN(yk)
      e(3)=DCOS(xk)
      CALL E2PNTS_2(e,ngr_out,grid)
      do 20 mo=1,nval
      if(nval.gt.1.and.mo.eq.1) goto 20
      do 10 j=1,4
      nclo=ilocra(ngr_out(j))
      if(nclo.lt.1) then
      ierr=1
      return
      endif
      if(nval.eq.1) then
          w(j)=fsl(nclo)
      else
          w(j)=fsl(nmm+nmma*(mo-2)+nclo)
      endif
   10 continue
      CALL TRIAI_2(grid,e,w,scell)
      value(mo)=scell
   20 continue
      return
      end
c******************************************************

max_stars_repo_name: NoiseCIEI/RayTomo

max_issues_repo_name: NoiseCIEI/RayTomo

max_forks_repo_name: NoiseCIEI/RayTomo


module pressure_resistance_flow_c
implicit none
private

contains

!
!###################################################################################
!
!*calculate_stats:* Calculates statistics for the feto-placental circulation model
  subroutine calculate_stats_c(FLOW_GEN_FILE, filename_len, image_voxel_size) bind(C, name="calculate_stats_c")
    use iso_c_binding, only: c_ptr
    use utils_c, only: strncpy
    use other_consts, only: MAX_FILENAME_LEN
    use arrays, only: dp
    use pressure_resistance_flow, only: calculate_stats
    implicit none

    integer,intent(in) :: filename_len
    type(c_ptr), value, intent(in) :: FLOW_GEN_FILE
    real(dp),intent(in) :: image_voxel_size
    character(len=MAX_FILENAME_LEN) :: filename_f

    call strncpy(filename_f, FLOW_GEN_FILE, filename_len)

#if defined _WIN32 && defined __INTEL_COMPILER
    call so_calculate_stats(filename_f, image_voxel_size)
#else
    call calculate_stats(filename_f, image_voxel_size)
#endif

  end subroutine calculate_stats_c

!
!###################################################################################
!

!!!###################################################################################

subroutine evaluate_prq_c(mesh_type,mesh_type_len,bc_type,bc_type_len,inlet_flow, &
               inlet_pressure,outlet_pressure) bind(C, name="evaluate_prq_c")

use iso_c_binding, only: c_ptr
use utils_c, only: strncpy
use other_consts, only: MAX_STRING_LEN
use arrays, only: dp
use pressure_resistance_flow, only: evaluate_prq
implicit none

type(c_ptr), value, intent(in) :: mesh_type,bc_type
integer,intent(in) :: mesh_type_len,bc_type_len
character(len=MAX_STRING_LEN) :: mesh_type_f,bc_type_f
real(dp),intent(in) :: inlet_flow,inlet_pressure,outlet_pressure

call strncpy(mesh_type_f, mesh_type, mesh_type_len)
call strncpy(bc_type_f, bc_type, bc_type_len)

#if defined _WIN32 && defined __INTEL_COMPILER
call so_evaluate_prq(mesh_type_f,bc_type_f,inlet_flow,inlet_pressure,outlet_pressure)
#else
call evaluate_prq(mesh_type_f,bc_type_f,inlet_flow,inlet_pressure,outlet_pressure)
#endif

end subroutine evaluate_prq_c

!###################################################################################
end module pressure_resistance_flow_c


text:
!
! wavy - A spectral ocean wave modeling and development framework
! Copyright (c) 2017, Wavebit Scientific LLC
! All rights reserved.
!
! Licensed under the BSD-3 clause license. See LICENSE for details.
!
program test_advection_2d

use mod_const
use mod_nondimensional
use mod_precision,only:rk => realkind
use mod_spectrum,only:spectrum_type
use mod_domain,only:domain_type
use mod_grid,only:grid_type
use mod_spectral_shapes
use mod_source_functions,only:sin_DCCM2012,sds_DCCM2012,sdt_DCCM2012,snl_DCCM2012
use mod_time_integration,only:backward_euler,forward_euler,exact_exponential
use mod_advection
use mod_utility,only:tile,range,ones
use netcdf

implicit none

type(spectrum_type),dimension(:,:),allocatable :: spec
!type(spectrum_type) :: spectrum
type(domain_type) :: domain,tend
!type(grid_type) :: grid
type(domain_type) :: sin_dom,sds_dom,sdt_dom,sbf_dom,snl_dom

real :: t0,t1

real(kind=rk) :: wspd
real(kind=rk) :: fetch
real(kind=rk),parameter :: grav = 9.8_rk
integer :: i,m,n,idm,nfreqs,jp,ip, j

real(kind=rk),dimension(:),allocatable :: freq
real(kind=rk),dimension(:,:),allocatable :: f
real(kind=rk),dimension(:,:),allocatable :: cg
real(kind=rk),dimension(:),allocatable :: dx
real(kind=rk),dimension(:,:),allocatable :: adv

real, dimension(:,:,:), allocatable :: swh

integer :: varid, ncid, xid, yid, ncid2
integer :: stat
integer :: swhid
integer :: xdimid,ydimid,tdimid

idm = 20

write(*,*)
write(*,*)'Initialize omnidirectional spectrum with JONSWAP shape;'
write(*,*)

wspd = 3e1_rk
fetch = 1e5_rk

! Initialize domain
domain = domain_type(grid = grid_type([1,1],[50,20],&
  dx = 1e3_rk*tile(ones(50,WAVY_INT),20),&
  dy = 1e3_rk*tile(ones(50,WAVY_INT),20)),&
  spectrum = spectrum_type(0.0313_rk,2._rk,1.1_rk,32,1e3_rk))

sin_dom = domain
sds_dom = domain
snl_dom = domain
sdt_dom = domain
sbf_dom = domain

spec = domain % getSpectrum()
tend = domain

! Initialize spectrum arrays

!do concurrent(i = 1:idm)
  spec(5,5) = donelanHamiltonHuiDirectionalSpectrum(f = domain % getFrequency(),&
                              theta = domain % getDirections(),&
                              wspd  = wspd,&
                              fpeak = jonswapPeakFrequency(5._rk,fetch,grav),&
                              theta_mean = -1.17_rk,&
                              grav  = grav)

!enddo
domain = spec

call cpu_time(t0)
call domain % writeJSON('domain.json',minify=.true.)
call cpu_time(t1)
write(*,*)'domain.json elapsed',t1-t0,'seconds'

write(*,fmt='(a)')'   wspd      Hs      Tp       Tm1      Tm2      mss'&
                //'      m0(f)    m1(f)    m2(f)'
write(*,fmt='(a)')'----------------------------------------------------------'&
                //'----------------------'

do n = 1,1000

  ip = 1
  jp = 5

  spec = domain % getSpectrum()

  ! wind input
  sin_dom = sin_DCCM2012(domain % getSpectrum(),wspd,wdir=1.17_rk,&
                         input_height=10._rk,ustar=1e-2_rk,vonkarman=0.4_rk)

  ! wave dissipation
  sds_dom = sds_DCCM2012(domain % getSpectrum(),sds_coefficient=42._rk,&
                         sds_power=2.5_rk,mss_coefficient=120._rk)

  ! dissipation from turbulence
  sdt_dom = sdt_DCCM2012(domain % getSpectrum(),1e-2_rk,1e-2_rk)

  ! non-linear wave-wave transfer
  snl_dom = snl_DCCM2012(domain % getSpectrum(),sds_dom % getSpectrum(),5._rk)

  !domain = forward_euler(domain,(sin_dom-sds_dom-sdt_dom)*domain+snl_dom,dt=1._rk)
  domain = exact_exponential(domain,sin_dom-sds_dom-sdt_dom,dt=10._rk)
  !domain = forward_euler(domain,snl_dom,dt=10._rk)

  write(*,fmt='(9(f8.4,1x))')wspd,spec(ip,jp) % significantWaveHeight(),&
    1./spec(ip,jp) % peakFrequency(),spec(ip,jp) % meanPeriod(),&
    spec(ip,jp) % meanPeriodZeroCrossing(),spec(ip,jp) % meanSquareSlope(),&
    spec(ip,jp) % frequencyMoment(0),spec(ip,jp) % frequencyMoment(1),&
    spec(ip,jp) % frequencyMoment(2)

  !domain = forward_euler(domain,&
  !  domain % advect(advectUpwind1stOrder1dRank1,1,&
  !  WAVY_OMNIDIRECTIONAL),dt=30._rk)

  !write(*,fmt='(20(f6.4,1x))')domain % significantWaveHeight()
  !write(*,fmt='(20(f7.4,1x))')domain % meanPeriod()

  domain = forward_euler(domain,&
    domain % advect(advectUpwind1stOrder2dRank2,[1,1]),dt=10._rk)

!call domain % writeJSON('domain.json',minify=.false.)
!
!
call nc_check(nf90_create('output.nc',nf90_clobber,ncid2))

call nc_check(nf90_def_dim(ncid2,'x',idm,xdimid))
call nc_check(nf90_def_dim(ncid2,'y',idm,ydimid))
stat = nf90_def_var(ncid2,'swh',nf90_float,[xdimid,ydimid],swhid)
call nc_check(nf90_enddef(ncid2))

stat = nf90_put_att(ncid2,swhid,name='description',values='Significant wave height')
stat = nf90_put_att(ncid2,swhid,name='units',values='m')

stat = nf90_put_var(ncid2,swhid,domain % significantWaveHeight(),start=[1,1],count=[idm,idm])
stat = nf90_close(ncid2)
enddo

contains
  SUBROUTINE nc_check(stat,varb)
!======================================================================!
!                                                                      !
! DESCRIPTION: Checks for NetCDF errors and if any, print and abort    !
!                                                                      !
!======================================================================!

INTEGER,INTENT(IN) :: stat
CHARACTER,INTENT(IN), optional :: VARB
!======================================================================!

IF(stat /= nf90_noerr)THEN
if (present(varb))then
  WRITE(*,*)(VARB)
endif
WRITE(*,*)'Error in NetCDF I/O'
WRITE(*,*)TRIM(nf90_strerror(stat))
STOP
ENDIF

ENDSUBROUTINE nc_check

endprogram test_advection_2d

max_stars_repo_name: dksasaki/wavy

max_issues_repo_name: dksasaki/wavy

max_forks_repo_name: dksasaki/wavy


program main
  implicit none
  real(8),parameter :: pi=4d0*atan(1d0)
  complex(8),parameter :: zI=(0d0,1d0)
  integer,parameter :: Nt = 50000, Nw = 6000
!  real(8),parameter :: ww_max = 120d0/27.211d0, dw = ww_max/Nw
  real(8),parameter :: ww_max = 300d0, dw = ww_max/Nw
  real(8)  :: Tprop = 100d0/0.02418d0
  real(8) :: tt(0:Nt),pz(0:Nt),weight(0:Nt)
  integer :: it,iw
  real(8) :: f1,f2,f3,f4,f5,f6,f7,f8,xx,ww,dt
  complex(8) :: zpzw

  do it = 0,Nt
    read(*,*)tt(it),f1,f2,f3,f4,f5,f6,f7,f8,pz(it)
  end do
  dt = tt(1) - tt(0)
  Tprop = tt(Nt)

  weight = 0d0
  do it = 0,Nt
!    if(tt(it) < 0.5d0*Tprop)then
!      weight(it) = 1d0
!    else
!      xx = (tt(it)-0.5d0*Tprop)/(tt(Nt)-0.5d0*Tprop)
      xx = (tt(it)-tt(0))/(tt(Nt)-tt(0))
!      weight(it) = 1d0 - 3d0*xx**2 + 2d0*xx**3
!      weight(it) = exp(-1d0/27.2d0*(tt(it)-Tprop))
!    end if
  end do
  weight = 1d0

  open(21,file="zpzw.out")
  do iw = 0,Nw
    ww = dw*iw
    zpzw = 0d0
    do it = 0,Nt
      zpzw = zpzw +exp(zI*ww*tt(it))*weight(it)*pz(it)
    end do
    zpzw = zpzw*dt
    write(21,"(999e26.16e3)")ww,abs(zpzw)**2
  end do
  close(21)

end program main


text:
Photographers Photographer

   Davis has a wide variety of resources for Photography

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module mod_inparam

use mod_types
use mod_cmdline

implicit none

type type_inparam
  integer(kind=ik):: benchmark_id !Auswahl der Laplacian-Implementierung (0,2)
  integer(kind=ik):: num_threads !Anzahl der Threads
  integer(kind=ik):: num_numas !Anzahl der Numa-Nodes
  integer(kind=ik):: numa_cores !Anzahl der Cores pro Numa-Nodes
  integer(kind=ik):: has_numas !Ob der compute node Numa-NOdes hat
  integer(kind=ik):: size_xx !Gebietsgroeßen zum Testen
  integer(kind=ik):: size_yy !Gebietsgroeßen zum Testen
  integer(kind=ik):: size_zz !Gebietsgroeßen zum Testen
  integer(kind=ik):: min_size_xx !Gebietsgroeßen zum Testen
  integer(kind=ik):: max_size_xx !Gebietsgroeßen zum Testen
  integer(kind=ik):: min_size_yy !Gebietsgroeßen zum Testen
  integer(kind=ik):: max_size_yy !Gebietsgroeßen zum Testen
  integer(kind=ik):: min_size_zz !Gebietsgroeßen zum Testen
  integer(kind=ik):: max_size_zz !Gebietsgroeßen zum Testen
  integer(kind=ik):: step_size_xx !Gebietsgroeßen zum Testen
  integer(kind=ik):: step_size_yy !Gebietsgroeßen zum Testen
  integer(kind=ik):: step_size_zz !Gebietsgroeßen zum Testen
  real(kind=rk):: min_time !Angabe für die minimale Dauer eines Testes (Groeßenordnung)
  integer(kind=ik):: blk_size_xx !Bloeckgroesse in X-Richtung
  integer(kind=ik):: blk_size_yy !Bloeckgroesse in Y-Richtung
  integer(kind=ik):: blk_size_zz !Bloeckgroesse in Z-Richtung
  integer(kind=ik):: min_blk_size_xx !Bloeckgroesse in X-Richtung
  integer(kind=ik):: max_blk_size_xx !Bloeckgroesse in X-Richtung
  integer(kind=ik):: step_blk_size_xx !Bloeckgroesse in X-Richtung
  integer(kind=ik):: min_blk_size_yy !Bloeckgroesse in Y-Richtung
  integer(kind=ik):: max_blk_size_yy !Bloeckgroesse in Y-Richtung
  integer(kind=ik):: step_blk_size_yy !Bloeckgroesse in Y-Richtung
  integer(kind=ik):: min_blk_size_zz !Bloeckgroesse in Z-Richtung
  integer(kind=ik):: max_blk_size_zz !Bloeckgroesse in Z-Richtung
  integer(kind=ik):: step_blk_size_zz !Bloeckgroesse in Z-Richtung
  integer(kind=ik):: num_repetitions !Anzahl der Wiederholungen
  integer(kind=ik):: check_solution !0 - no check, >0 chek each solution
  integer(kind=ik):: clear_cache !0 - Verdränge Keine Daten aus dem cache! 1-Verdränge Daten aus dem Cache
  integer(kind=ik):: deep_level !0 - Alle JPEGS aus dem Verzeichnis auslesen; > 0 Anzahl der JPEGS zu verarbeiten
  integer(kind=ik):: check_laplacian !0 - Arbeite mit JPEGs, sonst uu = xx*yy*yy+zz*zz*zz->dd=2.0_rk*xx+6.0*zz
  real(kind=rk)   :: scale_coeff!default 1 and is used to calculate the step h
  integer(kind=ik):: verbosity !Verbosity
  integer(kind=ik):: jobid
  character(len=1024) :: filepath !Dateiname für die Performance-Statistik
  character(len=1024) :: jpeg_topdir !Verzeichnis mit JPEG-Dateien
  character(len=1024) :: png_outdir !Verzeichnis für neue JPEG-Dateien
  character(len=1024) :: profile_filepath !Dateiname mit den Anwendungspahen als Timestamps gespeichert
  integer(kind=ik) :: has_profile !Ob die Anwendungsphasen, die Timestamps, gespeichert werden
end type type_inparam

type(type_inparam)   :: input_param
contains

subroutine input_param_print()

  write(*,'(A,I0)') "Eingabeparameter:"
  write(*,'(A,I0)') "verbosity:", input_param%verbosity
  write(*,'(A,I0)') "jobid:", input_param%jobid
  write(*,'(A,I0)') "benchmark_id:", input_param%benchmark_id
  write(*,'(A,I0)') "num_threads:", input_param%num_threads
  write(*,'(A,I0)') "size_xx:", input_param%size_xx
  write(*,'(A,I0)') "size_yy:", input_param%size_yy
  write(*,'(A,I0)') "size_zz:", input_param%size_zz
  write(*,'(A,I0)') "min_size_xx:", input_param%min_size_xx
  write(*,'(A,I0)') "max_size_xx:", input_param%max_size_xx
  write(*,'(A,I0)') "min_size_yy:", input_param%min_size_yy
  write(*,'(A,I0)') "max_size_yy:", input_param%max_size_yy
  write(*,'(A,I0)') "min_size_kk:", input_param%min_size_zz
  write(*,'(A,I0)') "max_size_kk:", input_param%max_size_zz
  write(*,'(A,I0)') "step_size_xx:", input_param%step_size_xx
  write(*,'(A,I0)') "step_size_yy:", input_param%step_size_yy
  write(*,'(A,I0)') "step_size_zz:", input_param%step_size_zz
  write(*,'(A,E13.6)') "min_time:", input_param%min_time
  write(*,'(A,I0)') "blk_size_xx:", input_param%blk_size_xx
  write(*,'(A,I0)') "blk_size_yy:", input_param%blk_size_yy
  write(*,'(A,I0)') "blk_size_zz:", input_param%blk_size_zz
  write(*,'(A,I0)') "min_blk_size_xx:", input_param%min_blk_size_xx
  write(*,'(A,I0)') "max_blk_size_xx:", input_param%max_blk_size_xx
  write(*,'(A,I0)') "step_blk_size_xx:", input_param%step_blk_size_xx
  write(*,'(A,I0)') "min_blk_size_yy:", input_param%min_blk_size_yy
  write(*,'(A,I0)') "max_blk_size_yy:", input_param%max_blk_size_yy
  write(*,'(A,I0)') "step_blk_size_yy:", input_param%step_blk_size_yy
  write(*,'(A,I0)') "min_blk_size_zz:", input_param%min_blk_size_zz
  write(*,'(A,I0)') "max_blk_size_zz:", input_param%max_blk_size_zz
  write(*,'(A,I0)') "step_blk_size_zz:", input_param%step_blk_size_zz
  write(*,'(A,I0)') "check_solution:", input_param%check_solution
  write(*,'(A,I0)') "num_repetitions:", input_param%num_repetitions
  write(*,'(A,I0)') "clear_cache:", input_param%clear_cache
  write(*,'(A,A)') "filepath:", trim(input_param%filepath)
  write(*,'(A,A)') "profile_filepath:", trim(input_param%profile_filepath)
  write(*,'(A,A)') "jpeg_topdir:", trim(input_param%jpeg_topdir)
  write(*,'(A,A)') "png_outdir:", trim(input_param%png_outdir)
  write(*,'(A,I0)') "deep_level:", input_param%deep_level
  write(*,'(A,I0)') "check_laplacian:", input_param%check_laplacian
  write(*,'(A,I0)') "num_numas:", input_param%num_numas
  write(*,'(A,I0)') "numa_cores:", input_param%numa_cores
  write(*,'(A,I0)') "has_numas:", input_param%has_numas
  write(*,'(A,E13.6)') "scale_coeff:", input_param%scale_coeff
  write(*,'(A,I0)') "has_profile:", input_param%has_profile

end subroutine input_param_print

subroutine input_param_read()
  character(len=2048) :: message
  character(len=2048) :: cmd
  integer(kind=ik) :: help

  call input_param_def()
  call get_command(cmd)
  message='Eingabeparameter:' &
    &//' -verbosity <int>' &
    &//' -jobid <int>' &
    &//' -benchmark_id <int> ' &
    &//' -num_threads <int> ' &
    &//' -min_size_x <int>' &
    &//' -max_size_x <int>' &
    &//' -min_size_y <int>' &
    &//' -max_size_y <int>' &
    &//' -min_size_z <int>' &
    &//' -max_size_z <int>' &
    &//' -size_x <int>' &
    &//' -size_y <int>' &
    &//' -size_z <int>' &
    &//' -step_size_x <int>' &
    &//' -step_size_y <int>' &
    &//' -step_size_z <int>' &
    &//' -min_time <real>' &
    &//' -blk_size_x <int>' &
    &//' -blk_size_y <int>' &
    &//' -blk_size_z <int>' &
    &//' -min_blk_size_x <int>' &
    &//' -max_blk_size_x <int>' &
    &//' -step_blk_size_x <int>' &
    &//' -min_blk_size_y <int>' &
    &//' -max_blk_size_y <int>' &
    &//' -step_blk_size_y <int>' &
    &//' -min_blk_size_z <int>' &
    &//' -max_blk_size_z <int>' &
    &//' -step_blk_size_z <int>' &
    &//' -num_repetitions <int>' &
    &//' -filepath <filename>' &
    &//' -check_solution <int>' &
    &//' -clear_cache <int>' &
    &//' -jpeg_topdir <filename>' &
    &//' -png_outdir <filename>' &
    &//' -profile_filepath <filename>' &
    &//' -deep_level <int>' &
    &//' -check_laplacian <int>' &
    &//' -scale_coeff <real>' &
    &//' -num_numas <int>' &
    &//' -numa_cores <int>' &
    &//' -has_numas <int>' &
    &//' -has_profile <int>' &
    &//' -help <int>'
    if(extract_command_parameter(cmd,'-help',stop_on_error=.false.,&
                value=help, syntax=message) ==0) then
      if(help .GE. 0) then
        write (*,'(A)') trim(message)
        if(help .GT. 0) then
          stop
        endif
      end if
  endif
  if(extract_command_parameter(cmd,'-verbosity',stop_on_error=.false.,&
                value=input_param%verbosity, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-jobid',stop_on_error=.false.,&
                value=input_param%jobid, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-benchmark_id',stop_on_error=.false.,&
                value=input_param%benchmark_id, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-num_threads',stop_on_error=.false.,&
                value=input_param%num_threads, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_time',stop_on_error=.false.,&
                value=input_param%min_time, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-size_x',stop_on_error=.false.,&
                value=input_param%size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-size_y',stop_on_error=.false.,&
                value=input_param%size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-size_z',stop_on_error=.false.,&
                value=input_param%size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_size_x',stop_on_error=.false.,&
                value=input_param%min_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_size_x',stop_on_error=.false.,&
                value=input_param%max_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_size_y',stop_on_error=.false.,&
                value=input_param%min_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_size_y',stop_on_error=.false.,&
                value=input_param%max_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_size_z',stop_on_error=.false.,&
                value=input_param%min_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_size_z',stop_on_error=.false.,&
                value=input_param%max_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_size_x',stop_on_error=.false.,&
                value=input_param%step_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_size_y',stop_on_error=.false.,&
                value=input_param%step_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_size_z',stop_on_error=.false.,&
                value=input_param%step_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-blk_size_x',stop_on_error=.false.,&
                value=input_param%blk_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-blk_size_y',stop_on_error=.false.,&
                value=input_param%blk_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-blk_size_z',stop_on_error=.false.,&
                value=input_param%blk_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_blk_size_x',stop_on_error=.false.,&
                value=input_param%min_blk_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_blk_size_x',stop_on_error=.false.,&
                value=input_param%max_blk_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_blk_size_x',stop_on_error=.false.,&
                value=input_param%step_blk_size_xx, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_blk_size_y',stop_on_error=.false.,&
                value=input_param%min_blk_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_blk_size_y',stop_on_error=.false.,&
                value=input_param%max_blk_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_blk_size_y',stop_on_error=.false.,&
                value=input_param%step_blk_size_yy, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-min_blk_size_z',stop_on_error=.false.,&
                value=input_param%min_blk_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-max_blk_size_z',stop_on_error=.false.,&
                value=input_param%max_blk_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-step_blk_size_z',stop_on_error=.false.,&
                value=input_param%step_blk_size_zz, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-check_solution',stop_on_error=.false.,&
                value=input_param%check_solution, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-clear_cache',stop_on_error=.false.,&
                value=input_param%clear_cache, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-filepath',stop_on_error=.false.,&
                value=input_param%filepath, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-jpeg_topdir',stop_on_error=.false.,&
                value=input_param%jpeg_topdir, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-png_outdir',stop_on_error=.false.,&
                value=input_param%png_outdir, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-deep_level',stop_on_error=.false.,&
                value=input_param%deep_level, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-check_laplacian',stop_on_error=.false.,&
                value=input_param%check_laplacian, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-scale_coeff',stop_on_error=.false.,&
                value=input_param%scale_coeff, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-has_numas',stop_on_error=.false.,&
                value=input_param%has_numas, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-num_numas',stop_on_error=.false.,&
                value=input_param%num_numas, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-numa_cores',stop_on_error=.false.,&
                value=input_param%numa_cores, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-num_repetitions',stop_on_error=.false.,&
                value=input_param%num_repetitions, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-profile_filepath',stop_on_error=.false.,&
                value=input_param%profile_filepath, syntax=message) /=0) then
  endif
  if(extract_command_parameter(cmd,'-has_profile',stop_on_error=.false.,&
                value=input_param%has_profile, syntax=message) /=0) then
  endif
  if(input_param%size_xx .gt. 0) then
    input_param%min_size_xx=input_param%size_xx
    input_param%max_size_xx=input_param%size_xx
    input_param%step_size_xx=input_param%size_xx
  endif
  if(input_param%size_yy .gt. 0) then
    input_param%min_size_yy=input_param%size_yy
    input_param%max_size_yy=input_param%size_yy
    input_param%step_size_yy=input_param%size_yy
  endif
  if(input_param%size_zz .gt. 0) then
    input_param%min_size_zz=input_param%size_zz
    input_param%max_size_zz=input_param%size_zz
    input_param%step_size_zz=input_param%size_zz
  endif

  if(input_param%blk_size_xx .gt. 0) then
    input_param%min_blk_size_xx=input_param%blk_size_xx
    input_param%max_blk_size_xx=input_param%blk_size_xx
    input_param%step_blk_size_xx=input_param%blk_size_xx
  endif
  if(input_param%blk_size_yy .gt. 0) then
    input_param%min_blk_size_yy=input_param%blk_size_yy
    input_param%max_blk_size_yy=input_param%blk_size_yy
    input_param%step_blk_size_yy=input_param%blk_size_yy
  endif
  if(input_param%blk_size_zz .gt. 0) then
    input_param%min_blk_size_zz=input_param%blk_size_zz
    input_param%max_blk_size_zz=input_param%blk_size_zz
    input_param%step_blk_size_zz=input_param%blk_size_zz
  endif
  if(input_param%check_laplacian .gt. 0) then
    input_param%min_blk_size_zz=input_param%blk_size_zz
    input_param%max_blk_size_zz=input_param%blk_size_zz
    input_param%step_blk_size_zz=input_param%blk_size_zz
  endif
  call input_param_check()

end subroutine input_param_read

subroutine input_param_def()

  input_param%benchmark_id=0
  input_param%num_threads=1
  input_param%num_numas=2
  input_param%numa_cores=12
  input_param%has_numas=0
  input_param%size_xx=0
  input_param%size_yy=0
  input_param%size_zz=0
  input_param%min_size_xx=1920
  input_param%min_size_yy=1920
  input_param%min_size_zz=180
  input_param%max_size_xx=1920
  input_param%max_size_yy=1920
  input_param%max_size_zz=180
  input_param%step_size_xx=64
  input_param%step_size_yy=64
  input_param%step_size_zz=64
  input_param%min_time=1.0
  input_param%blk_size_xx=0
  input_param%blk_size_yy=0
  input_param%blk_size_zz=0
  input_param%min_blk_size_xx=128
  input_param%max_blk_size_xx=128
  input_param%step_blk_size_xx=128
  input_param%min_blk_size_yy=128
  input_param%max_blk_size_yy=128
  input_param%step_blk_size_yy=128
  input_param%min_blk_size_zz=180
  input_param%max_blk_size_zz=180
  input_param%step_blk_size_zz=180
  input_param%num_repetitions=0
  input_param%check_solution=1
  input_param%clear_cache=1
  input_param%deep_level=0
  input_param%check_laplacian=0
  input_param%scale_coeff=1.0
  input_param%has_profile=0
  write(input_param%filepath,'(A,I0,4(A,I0,A,I0),A,I0,A)') &
      & "./bench_",input_param%benchmark_id, &
      & "_size_",input_param%min_size_xx,"-",input_param%max_size_xx, &
      & "_size_",input_param%min_size_yy,"-",input_param%max_size_yy, &
      & "_size_",input_param%min_size_zz,"-",input_param%max_size_zz, &
      & "_block_",input_param%blk_size_xx,"x",input_param%blk_size_yy, &
      & "x",input_param%blk_size_zz, &
      & ".csv"
  write(input_param%jpeg_topdir,'(A)') "/nas_home/hpcdkhab/example_cache_blocking/fotos/"
  write(input_param%png_outdir,'(A)') "/nas_home/hpcdkhab/example_cache_blocking/new_fotos/"
  write(input_param%png_outdir,'(A)') "/nas_home/hpcdkhab/example_cache_blocking/data/profile.csv"
end subroutine input_param_def

subroutine input_param_check()

  if(input_param%benchmark_id .gt. 4) then
    write(*,'(A,I0)') "Error: check benchmark_id (0..4):",&
      & input_param%benchmark_id
    stop
  endif
  if(input_param%min_size_xx .gt. input_param%max_size_xx) then
    write(*,'(A,I0)') "Error: check min_size_x and max_size_x (max_size must be greater than min_size)"
    write(*,'(2(A,I0))') "min_size_x: ", input_param%min_size_xx, "; max_size_x: ", input_param%max_size_xx
    stop
  endif
  if(input_param%min_size_yy .gt. input_param%max_size_yy) then
    write(*,'(A,I0)') "Error: check min_size_y and max_size_y (max_size must be greater than min_size)"
    write(*,'(2(A,I0))') "min_size_y: ", input_param%min_size_yy, "; max_size_y: ", input_param%max_size_yy
    stop
  endif
  if(input_param%min_size_zz .gt. input_param%max_size_zz) then
    write(*,'(A,I0)') "Error: check min_size_z and max_size_z (max_size must be greater than min_size)"
    write(*,'(2(A,I0))') "min_size_z: ", input_param%min_size_zz, "; max_size_z: ", input_param%max_size_zz
    stop
  endif
  if(input_param%max_blk_size_xx .gt. input_param%min_size_xx) then
    write(*,'(A,I0)') "Error: block_size_x (block_size must be equal or less  than min_size)"
    stop
  endif
  if(input_param%max_blk_size_yy .gt. input_param%min_size_yy) then
    write(*,'(A,I0)') "Error: block_size_y (block_size must be equal or less  than min_size)"
    stop
  endif
  if(input_param%max_blk_size_zz .gt. input_param%min_size_zz) then
    write(*,'(A,I0)') "Error: block_size_z (block_size must be equal or less  than min_size)"
    stop
  endif
  if(input_param%has_numas .gt. 0) then
    if(input_param%num_threads .gt. input_param%num_numas*input_param%numa_cores) then
      write(*,'(A,I0)') "Error: num_threads greater than number of cores=num_numas*numa_cores"
      stop
    endif
  endif
  if(input_param%min_time .lt. 1.e-6) then
    write(*,'(A,I0)') "Error: check min_time (too small)"
    stop
  endif

end subroutine input_param_check

end module mod_inparam


text:
      subroutine qqb_z2jet_v(p,msqv)
      implicit none
************************************************************************
*     Author: R.K. Ellis                                               *
*     September 2001.                                                  *
*                                                                      *
*     Calculate the virtual matrix element squared and                 *
*     subtraction terms for the process                                *
*                                                                      *
*     q(-p1) + qbar(-p2) --> Z + j(p5) + j(p6)                         *
*                            |                                         *
*                            --> e^-(p3) + e^+(p4)                     *
*                                                                      *
*     where the partons are either q(p5) and qbar(p6) [Qflag = .true.] *
*                               or g(p5) and g(p6)    [Gflag = .true.] *
*                                                                      *
*     The value of COLOURCHOICE determines which colour structures     *
*     are included in the terms for the QQGG piece                     *
*                                                                      *
************************************************************************
*                                                                      *
*     By setting compare=.true. it is possible to test the squared     *
*     tree amplitudes in this routine with those from qqb_z2jet.f      *
*      (checked on 10/22/01 by JMC)                                    *
*                                                                      *
************************************************************************

      include 'constants.f'
      include 'qcdcouple.f'
      include 'ewcouple.f'
      include 'ewcharge.f'
      include 'zcouple.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'zprods_com.f'
      include 'epinv.f'
      include 'scheme.f'
      include 'flags.f'
      include 'msq_cs.f'
      include 'lc.f'
      include 'first.f'
      double precision msq(-nf:nf,-nf:nf),msqv(-nf:nf,-nf:nf),
     . mmsq_qqb(2,2),mmsq_qbq(2,2),mmsq_gq(2,2),mmsq_qg(2,2),
     . mmsq_qbg(2,2),mmsq_gqb(2,2),mmsq_gg(2,2),
     . p(mxpart,4),pswap(mxpart,4),fac
      double complex mmsq_qqb_ax(2,2),mmsq_qbq_ax(2,2),
     . mmsq_gq_ax(2,2),mmsq_qg_ax(2,2),mmsq_qbg_ax(2,2),
     . mmsq_gqb_ax(2,2),mmsq_gg_ax(2,2),
     . mmsq_qqb_vec(2,2),mmsq_qbq_vec(2,2),
     . mmsq_gq_vec(2,2),mmsq_qg_vec(2,2),mmsq_qbg_vec(2,2),
     . mmsq_gqb_vec(2,2),mmsq_gg_vec(2,2)
      double complex atreez,a61z,a62z,a63z,prop,vcouple(2),
     . tamp,lamp,tampup,lampup,tampdo,lampdo,tamps,lamps,lampx,lampsx
      integer nu,j,k,cs,polq,polb,polz
      integer nup,ndo,rvcolourchoice
      integer, parameter :: jj(-nf:nf)=(/-1,-2,-1,-2,-1,0,1,2,1,2,1/)
      double precision subuv(0:2)
      double precision faclo,v2(2),vQ(nf,2),idfac
      double precision mqq(0:2,fn:nf,fn:nf)
      double precision msqx(0:2,-nf:nf,-nf:nf,-nf:nf,-nf:nf)
      double precision msqx_cs(0:2,-nf:nf,-nf:nf)
      double complex atreez_526143(2,2,2),atreez_251643(2,2,2),
     .               atreez_625143(2,2,2),atreez_261543(2,2,2),
     .               atreez_162543(2,2,2),atreez_615243(2,2,2),
     .               atreez_152643(2,2,2),atreez_516243(2,2,2),
     .               atreez_562143(2,2,2),atreez_651243(2,2,2),
     .               atreez_265143(2,2,2),atreez_621543(2,2,2),
     .               atreez_561243(2,2,2),atreez_652143(2,2,2),
     .               atreez_165243(2,2,2),atreez_612543(2,2,2),
     .               a61z_526143(2,2,2),a61z_251643(2,2,2),
     .               a61z_625143(2,2,2),a61z_261543(2,2,2),
     .               a61z_162543(2,2,2),a61z_615243(2,2,2),
     .               a61z_152643(2,2,2),a61z_516243(2,2,2),
     .               a61z_562143(2,2,2),a61z_651243(2,2,2),
     .               a61z_265143(2,2,2),a61z_621543(2,2,2),
     .               a61z_561243(2,2,2),a61z_652143(2,2,2),
     .               a61z_165243(2,2,2),a61z_612543(2,2,2),
     .               a62z_625143(2,2,2),a62z_261543(2,2,2),
     .               a62z_526143(2,2,2),a62z_251643(2,2,2),
     .               a62z_152643(2,2,2),a62z_516243(2,2,2),
     .               a62z_162543(2,2,2),a62z_615243(2,2,2),
     .               a62z_265143(2,2,2),a62z_621543(2,2,2),
     .               a62z_562143(2,2,2),a62z_651243(2,2,2),
     .               a62z_165243(2,2,2),a62z_612543(2,2,2),
     .               a62z_561243(2,2,2),a62z_652143(2,2,2),
     .               a63z_526143(2,2,2),a63z_625143(2,2,2),
     .               a63z_162543(2,2,2),a63z_152643(2,2,2),
     .               a63z_562143(2,2,2),a63z_265143(2,2,2),
     .               a63z_561243(2,2,2),a63z_165243(2,2,2)

      logical compare,checkvector,checkaxial
      common/mqq/mqq
      common/rvcolourchoice/rvcolourchoice
      parameter (nup=2,ndo=nf-nup)
      parameter (compare=.false.)
c--- These parameters allow for a point-by-point comparison of the
c--- vector and axial pieces with MadLoop. They should normally
c--- both be set to .false.
      parameter(checkvector=.false.,checkaxial=.false.)
!$omp threadprivate(/mqq/,/rvcolourchoice/)

      if (Qflag .and. Gflag) then
        write(6,*) 'Both Qflag and Gflag cannot be true'
        write(6,*) 'They are set in file options.DAT'
        write(6,*) 'Failed in qqb_z2jet_v.f'
        stop
      endif

      if (first) then
        first=.false.
        if ((Gflag) .or. (QandGflag)) then
          write(*,*) 'Using QQGG (VIRTUAL) matrix elements'
          write(*,*) '[LC is     N   ]'
          write(*,*) '[SLC is   1/N  ]'
          write(*,*) '[SSLC is 1/N**3]'
        endif
        if ((Qflag) .or. (QandGflag)) then
          write(*,*) 'Using QQBQQB (VIRTUAL) matrix elements'
          write(*,*) '[LC is   1 ]'
          write(*,*) '[SLC is 1/N]'
        endif
        if     (rvcolourchoice .eq. 1) then
          write(*,*) 'Leading colour only in VIRTUAL'
        elseif (rvcolourchoice .eq. 2) then
          write(*,*) 'Sub-leading colour only in VIRTUAL'
        elseif (rvcolourchoice .eq. 3) then
          write(*,*) 'Sub-sub-leading colour only in VIRTUAL'
        elseif (rvcolourchoice .eq. 0) then
          write(*,*) 'Total of all colour structures in VIRTUAL'
        else
          write(*,*) 'Bad colourchoice'
          stop
        endif
      endif

c--- get special phase space point for MadLoop comparison
      if (checkvector .or. checkaxial) then
        call ps_check(p,1)
      endif

      scheme='dred'
c--- calculate the lowest order matrix element and fill the
c--- common block twopij with s_{ij}
      call qqb_z2jetx(p,msq,mqq,msqx,msqx_cs)

c--- write out ug Born amplitude when checking
      if (checkvector .or. checkaxial) then
        write(6,*) 'Madloop check: ug Born',msq(2,0)
      endif

      do j=-nf,nf
      do k=-nf,nf
      msqv(j,k)=0d0
      enddo
      enddo

      prop=s(3,4)/dcmplx(s(3,4)-zmass**2,zmass*zwidth)

      v2(1)=l1
      v2(2)=r1

      do j=1,nf
        vQ(j,1)=L(j)
        vQ(j,2)=R(j)
      enddo

c---- DEBUG: check vector piece
      if (checkvector) then
        prop=czip
      endif

************************************************************************
*     Endpoint contributions from QQGG matrix elements                 *
*     Loop matrix elements are also initialized here                   *
************************************************************************
      if (Gflag) then
c----UV counterterm contains the finite renormalization to arrive
c----at MS bar scheme.
      if     (colourchoice .eq. 1) then
        subuv(1)=2d0*xn*(epinv*(11d0-2d0*dble(nf)/xn)-1d0)/6d0
        subuv(2)=subuv(1)
      elseif (colourchoice .eq. 2) then
        subuv(0)=2d0*xn*(epinv*(11d0-2d0*dble(nf)/xn)-1d0)/6d0
      elseif (colourchoice .eq. 3) then
c--- all zero already
      elseif (colourchoice .eq. 0) then
        subuv(1)=2d0*xn*(epinv*(11d0-2d0*dble(nf)/xn)-1d0)/6d0
        subuv(2)=subuv(1)
        subuv(0)=subuv(1)
      endif

c--- transfer lowest order matrix elements
c--- NB: this breaks the routine if Qflag = Gflag = .true.

      do cs=0,2
        do j=-nf,nf
        do k=-nf,nf
        msq_cs(cs,j,k)=msqx_cs(cs,j,k)
        enddo
        enddo
      enddo

c--- Now calculate the relevant lowest-order matrix elements
c--- for each possible initial state from the QQGG contribution

c---  calculate the qqb terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p2)+g(p6)+g(p5)+qb(p1)+lbar(p4)+l(p3)
      do nu=1,4
      pswap(1,nu)=p(2,nu)
      pswap(2,nu)=p(6,nu)
      pswap(3,nu)=p(5,nu)
      pswap(4,nu)=p(1,nu)
      pswap(5,nu)=p(4,nu)
      pswap(6,nu)=p(3,nu)
      enddo
      call spinoru(6,pswap,za,zb)
      call xzqqgg_v(mmsq_qqb,mmsq_qqb_vec,mmsq_qqb_ax)

c--- obtain qbq from qqb by symmetry
      do polq=1,2
      do polz=1,2
        mmsq_qbq(polq,polz)=mmsq_qqb(3-polq,polz)
        mmsq_qbq_vec(polq,polz)=mmsq_qqb_vec(3-polq,polz)
        mmsq_qbq_ax(polq,polz)=-mmsq_qqb_ax(3-polq,polz)
      enddo
      enddo
c---  calculate the qbq terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p1)+g(p6)+g(p5)+qb(p2)+lbar(p4)+l(p3)
c      do nu=1,4
c      pswap(1,nu)=p(1,nu)
c      pswap(2,nu)=p(6,nu)
c      pswap(3,nu)=p(5,nu)
c      pswap(4,nu)=p(2,nu)
c      pswap(5,nu)=p(4,nu)
c      pswap(6,nu)=p(3,nu)
c      enddo
c      call spinoru(6,pswap,za,zb)
c      call xzqqgg_v(mmsq_qbq,mmsq_qbq_vec,mmsq_qbq_ax)
c      do polq=1,2
c      do polz=1,2
c        write(6,*) mmsq_qbq(polq,polz)-mmsq_qqb(3-polq,polz)
c        write(6,*) mmsq_qbq_vec(polq,polz)-mmsq_qqb_vec(3-polq,polz)
c        write(6,*) mmsq_qbq_ax(polq,polz)+mmsq_qqb_ax(3-polq,polz)
c      enddo
c      enddo
c      pause

c---  calculate the gq terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p5)+g(p6)+g(p1)+qb(p2)+lbar(p4)+l(p3)
      do nu=1,4
      pswap(1,nu)=p(5,nu)
      pswap(2,nu)=p(6,nu)
      pswap(3,nu)=p(1,nu)
      pswap(4,nu)=p(2,nu)
      pswap(5,nu)=p(4,nu)
      pswap(6,nu)=p(3,nu)
      enddo
      call spinoru(6,pswap,za,zb)
      call xzqqgg_v(mmsq_gq,mmsq_gq_vec,mmsq_gq_ax)

c--- obtain gqb from gq by symmetry
      do polq=1,2
      do polz=1,2
        mmsq_gqb(polq,polz)=mmsq_gq(3-polq,polz)
        mmsq_gqb_vec(polq,polz)=mmsq_gq_vec(3-polq,polz)
        mmsq_gqb_ax(polq,polz)=-mmsq_gq_ax(3-polq,polz)
      enddo
      enddo
c---  calculate the gqb terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p2)+g(p6)+g(p1)+qb(p5)+lbar(p4)+l(p3)
c      do nu=1,4
c      pswap(1,nu)=p(2,nu)
c      pswap(2,nu)=p(6,nu)
c      pswap(3,nu)=p(1,nu)
c      pswap(4,nu)=p(5,nu)
c      pswap(5,nu)=p(4,nu)
c      pswap(6,nu)=p(3,nu)
c      enddo
c      call spinoru(6,pswap,za,zb)
c      call xzqqgg_v(mmsq_gqb,mmsq_gqb_vec,mmsq_gqb_ax)

c---  calculate the qg terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p5)+g(p6)+g(p2)+qb(p1)+lbar(p4)+l(p3)
      do nu=1,4
      pswap(1,nu)=p(5,nu)
      pswap(2,nu)=p(6,nu)
      pswap(3,nu)=p(2,nu)
      pswap(4,nu)=p(1,nu)
      pswap(5,nu)=p(4,nu)
      pswap(6,nu)=p(3,nu)
      enddo
      call spinoru(6,pswap,za,zb)
      call xzqqgg_v(mmsq_qg,mmsq_qg_vec,mmsq_qg_ax)

c--- obtain qbg from qg by symmetry
      do polq=1,2
      do polz=1,2
        mmsq_qbg(polq,polz)=mmsq_qg(3-polq,polz)
        mmsq_qbg_vec(polq,polz)=mmsq_qg_vec(3-polq,polz)
        mmsq_qbg_ax(polq,polz)=-mmsq_qg_ax(3-polq,polz)
      enddo
      enddo
c---  calculate the qbg terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p1)+g(p6)+g(p2)+qb(p5)+lbar(p4)+l(p3)
c      do nu=1,4
c      pswap(1,nu)=p(1,nu)
c      pswap(2,nu)=p(6,nu)
c      pswap(3,nu)=p(2,nu)
c      pswap(4,nu)=p(5,nu)
c      pswap(5,nu)=p(4,nu)
c      pswap(6,nu)=p(3,nu)
c      enddo
c      call spinoru(6,pswap,za,zb)
c      call xzqqgg_v(mmsq_qbg,mmsq_qbg_vec,mmsq_qbg_ax)

c--- calculate the gg terms
c-BDKW  0--> q(p1)+g(p2)+g(p3)+qb(p4)+lbar(p5)+l(p6)
CALL    0--> q(p6)+g(p1)+g(p2)+qb(p5)+lbar(p4)+l(p3)
      do nu=1,4
      pswap(1,nu)=p(6,nu)
      pswap(2,nu)=p(1,nu)
      pswap(3,nu)=p(2,nu)
      pswap(4,nu)=p(5,nu)
      pswap(5,nu)=p(4,nu)
      pswap(6,nu)=p(3,nu)
      enddo
      call spinoru(6,pswap,za,zb)
      call xzqqgg_v(mmsq_gg,mmsq_gg_vec,mmsq_gg_ax)
      endif

************************************************************************
*     Endpoint contributions from QQQQ matrix elements                 *
************************************************************************
      if (Qflag) then
c--- UV counter-term is already included in a6routine.f
      subuv(1)=0d0
      subuv(2)=subuv(1)
      subuv(0)=subuv(1)

c--- transfer lowest order matrix elements
c--- NB: this breaks the routine if Qflag = Gflag = .true.

      do cs=0,2
        do j=-nf,nf
        do k=-nf,nf
        msq_cs(cs,j,k)=mqq(cs,j,k)
        enddo
        enddo
      enddo

      endif

c--- Add VIRTUAL terms

************************************************************************
*     Include loop contributions from QQGG matrix elements             *
************************************************************************
      if (Gflag) then

c--- compute correct vector-like coupling for diagrams with Z coupled to a loop
      vcouple(1)=czip
      vcouple(2)=czip
      do j=1,nf
      do polz=1,2
      vcouple(polz)=vcouple(polz)
     & +Q(j)*q1+0.5d0*(vQ(j,1)+vQ(j,2))*v2(polz)*prop
      enddo
      enddo

      do j=-nf,nf
      do k=-nf,nf
      do polq=1,2
      do polz=1,2

c--- quark-antiquark
      if ((j .gt. 0) .and. (k .lt. 0)) then
        if (j .eq. -k)
     .  msqv(j,k)=msqv(j,k)+half*(aveqq/avegg)*(mmsq_qqb(polq,polz)*(
     .          cdabs(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_qqb_vec(polq,polz)
     .        *dconjg(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_qqb_ax(polq,polz)
     .        *dconjg(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

c--- antiquark-quark
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
        if (j .eq. -k)
     .  msqv(j,k)=msqv(j,k)+half*(aveqq/avegg)*(mmsq_qbq(polq,polz)*(
     .          cdabs(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_qbq_vec(polq,polz)
     .        *dconjg(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_qbq_ax(polq,polz)
     .        *dconjg(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

c--- quark-gluon
      elseif ((j .gt. 0) .and. (k .eq. 0)) then

************************** BEGIN MADLOOP CHECKING CODE ************************
        if (checkvector) then
c--- MadLoop check: vector couplings only - no Z coupling to leptons and quarks
        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(
     .                   +dble(mmsq_qg_vec(polq,polz)
     .          *dconjg(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)*vcouple(polz))
     .                   )
        endif

        if (checkaxial) then
c--- MadLoop check: axial couplings only - photon contribution removed
c------- full Z in Born (no photon)
        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(
     .              +dble(mmsq_qg_ax(polq,polz)
     .        *dconjg(Q(j)*q1*0+vQ(j,polq)*v2(polz)*prop)
     .         *(v2(polz)*prop)/sin2w))

c------- axial Z everywehere (Born and virt, quarks and leptons)
c        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(
c     .              +dble(mmsq_qg_ax(polq,polz)
c     .         *dconjg(Q(j)*q1*0+
c     .           (vQ(j,polq)-vQ(j,3-polq))*(v2(polz)-v2(3-polz))*prop)
c     .         *((v2(polz)-v2(3-polz))*prop)/sin2w))/8d0
      endif
*************************** END MADLOOP CHECKING CODE *************************

      if ((checkvector.eqv..false.).and.(checkaxial.eqv..false.)) then
c--- normal case
      msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(mmsq_qg(polq,polz)*(
     .        cdabs(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_qg_vec(polq,polz)
     .        *dconjg(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_qg_ax(polq,polz)
     .        *dconjg(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

        endif

c--- antiquark-gluon
      elseif ((j .lt. 0) .and. (k .eq. 0)) then

        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(mmsq_qbg(polq,polz)*(
     .          cdabs(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_qbg_vec(polq,polz)
     .        *dconjg(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_qbg_ax(polq,polz)
     .        *dconjg(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

c--- gluon-quark
      elseif ((j .eq. 0) .and. (k .gt. 0)) then
        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(mmsq_gq(polq,polz)*(
     .          cdabs(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_gq_vec(polq,polz)
     .        *dconjg(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_gq_ax(polq,polz)
     .        *dconjg(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

c--- gluon-antiquark
      elseif ((j .eq. 0) .and. (k .lt. 0)) then
        msqv(j,k)=msqv(j,k)+(aveqg/avegg)*(mmsq_gqb(polq,polz)*(
     .          cdabs(Q(-k)*q1+vQ(-k,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_gqb_vec(polq,polz)
     .        *dconjg(Q(-k)*q1+vQ(-k,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_gqb_ax(polq,polz)
     .        *dconjg(Q(-k)*q1+vQ(-k,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))

c--- gluon-gluon
      elseif ((j .eq. 0) .and. (k .eq. 0)) then
        msqv(j,k)=msqv(j,k)+dfloat(ndo)*(mmsq_gg(polq,polz)*(
     .             cdabs(Q(1)*q1+vQ(1,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_gg_vec(polq,polz)
     .        *dconjg(Q(1)*q1+vQ(1,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_gg_ax(polq,polz)
     .        *dconjg(Q(1)*q1+vQ(1,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))
        msqv(j,k)=msqv(j,k)+dfloat(nup)*(mmsq_gg(polq,polz)*(
     .          cdabs(Q(2)*q1+vQ(2,polq)*v2(polz)*prop)**2)
     .                     +dble(mmsq_gg_vec(polq,polz)
     .        *dconjg(Q(2)*q1+vQ(2,polq)*v2(polz)*prop)*vcouple(polz))
     .                 +dble(mmsq_gg_ax(polq,polz)
     .        *dconjg(Q(2)*q1+vQ(2,polq)*v2(polz)*prop)
     .          *(v2(polz)*prop)/sin2w))
      endif
      enddo
      enddo
      enddo
      enddo

      endif

************************************************************************
*     Include loop contributions from QQQQ matrix elements             *
************************************************************************
      if (Qflag) then

      call spinoru(6,p,za,zb)

      faclo=4d0*V*aveqq*esq**2*gsq**2
      fac=faclo*xn*0.5d0*ason2pi

c--- Set-up the desired amplitudes, in order to minimize number of calls
      do polq=1,2
      do polz=1,2
      do polb=1,2
c--- atreez
      atreez_526143(polq,polb,polz)=
     .       atreez(polq,polb,polz,5,2,6,1,4,3,za,zb)
      atreez_251643(polq,polb,polz)=
     .       atreez(polq,polb,polz,2,5,1,6,4,3,za,zb)
      atreez_625143(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,2,5,1,4,3,za,zb)
      atreez_261543(polq,polb,polz)=
     .       atreez(polq,polb,polz,2,6,1,5,4,3,za,zb)
      atreez_162543(polq,polb,polz)=
     .       atreez(polq,polb,polz,1,6,2,5,4,3,za,zb)
      atreez_615243(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,1,5,2,4,3,za,zb)
      atreez_152643(polq,polb,polz)=
     .       atreez(polq,polb,polz,1,5,2,6,4,3,za,zb)
      atreez_516243(polq,polb,polz)=
     .       atreez(polq,polb,polz,5,1,6,2,4,3,za,zb)
      atreez_562143(polq,polb,polz)=
     .       atreez(polq,polb,polz,5,6,2,1,4,3,za,zb)
      atreez_651243(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,5,1,2,4,3,za,zb)
      atreez_265143(polq,polb,polz)=
     .       atreez(polq,polb,polz,2,6,5,1,4,3,za,zb)
      atreez_621543(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,2,1,5,4,3,za,zb)
      atreez_561243(polq,polb,polz)=
     .       atreez(polq,polb,polz,5,6,1,2,4,3,za,zb)
      atreez_652143(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,5,2,1,4,3,za,zb)
      atreez_165243(polq,polb,polz)=
     .       atreez(polq,polb,polz,1,6,5,2,4,3,za,zb)
      atreez_612543(polq,polb,polz)=
     .       atreez(polq,polb,polz,6,1,2,5,4,3,za,zb)
c--- a61z
      a61z_526143(polq,polb,polz)=a61z(polq,polb,polz,5,2,6,1,4,3,za,zb)
      a61z_251643(polq,polb,polz)=a61z(polq,polb,polz,2,5,1,6,4,3,za,zb)
      a61z_625143(polq,polb,polz)=a61z(polq,polb,polz,6,2,5,1,4,3,za,zb)
      a61z_261543(polq,polb,polz)=a61z(polq,polb,polz,2,6,1,5,4,3,za,zb)
      a61z_162543(polq,polb,polz)=a61z(polq,polb,polz,1,6,2,5,4,3,za,zb)
      a61z_615243(polq,polb,polz)=a61z(polq,polb,polz,6,1,5,2,4,3,za,zb)
      a61z_152643(polq,polb,polz)=a61z(polq,polb,polz,1,5,2,6,4,3,za,zb)
      a61z_516243(polq,polb,polz)=a61z(polq,polb,polz,5,1,6,2,4,3,za,zb)
      a61z_562143(polq,polb,polz)=a61z(polq,polb,polz,5,6,2,1,4,3,za,zb)
      a61z_651243(polq,polb,polz)=a61z(polq,polb,polz,6,5,1,2,4,3,za,zb)
      a61z_265143(polq,polb,polz)=a61z(polq,polb,polz,2,6,5,1,4,3,za,zb)
      a61z_621543(polq,polb,polz)=a61z(polq,polb,polz,6,2,1,5,4,3,za,zb)
      a61z_561243(polq,polb,polz)=a61z(polq,polb,polz,5,6,1,2,4,3,za,zb)
      a61z_652143(polq,polb,polz)=a61z(polq,polb,polz,6,5,2,1,4,3,za,zb)
      a61z_165243(polq,polb,polz)=a61z(polq,polb,polz,1,6,5,2,4,3,za,zb)
      a61z_612543(polq,polb,polz)=a61z(polq,polb,polz,6,1,2,5,4,3,za,zb)
c--- a62z
      a62z_625143(polq,polb,polz)=a62z(polq,polb,polz,6,2,5,1,4,3,za,zb)
      a62z_261543(polq,polb,polz)=a62z(polq,polb,polz,2,6,1,5,4,3,za,zb)
      a62z_526143(polq,polb,polz)=a62z(polq,polb,polz,5,2,6,1,4,3,za,zb)
      a62z_251643(polq,polb,polz)=a62z(polq,polb,polz,2,5,1,6,4,3,za,zb)
      a62z_152643(polq,polb,polz)=a62z(polq,polb,polz,1,5,2,6,4,3,za,zb)
      a62z_516243(polq,polb,polz)=a62z(polq,polb,polz,5,1,6,2,4,3,za,zb)
      a62z_162543(polq,polb,polz)=a62z(polq,polb,polz,1,6,2,5,4,3,za,zb)
      a62z_615243(polq,polb,polz)=a62z(polq,polb,polz,6,1,5,2,4,3,za,zb)
      a62z_265143(polq,polb,polz)=a62z(polq,polb,polz,2,6,5,1,4,3,za,zb)
      a62z_621543(polq,polb,polz)=a62z(polq,polb,polz,6,2,1,5,4,3,za,zb)
      a62z_562143(polq,polb,polz)=a62z(polq,polb,polz,5,6,2,1,4,3,za,zb)
      a62z_651243(polq,polb,polz)=a62z(polq,polb,polz,6,5,1,2,4,3,za,zb)
      a62z_165243(polq,polb,polz)=a62z(polq,polb,polz,1,6,5,2,4,3,za,zb)
      a62z_612543(polq,polb,polz)=a62z(polq,polb,polz,6,1,2,5,4,3,za,zb)
      a62z_561243(polq,polb,polz)=a62z(polq,polb,polz,5,6,1,2,4,3,za,zb)
      a62z_652143(polq,polb,polz)=a62z(polq,polb,polz,6,5,2,1,4,3,za,zb)
c--- a63z
      a63z_526143(polq,polb,polz)=a63z(polq,polb,polz,5,2,6,1,4,3,za,zb)
      a63z_625143(polq,polb,polz)=a63z(polq,polb,polz,6,2,5,1,4,3,za,zb)
      a63z_162543(polq,polb,polz)=a63z(polq,polb,polz,1,6,2,5,4,3,za,zb)
      a63z_152643(polq,polb,polz)=a63z(polq,polb,polz,1,5,2,6,4,3,za,zb)
      a63z_562143(polq,polb,polz)=a63z(polq,polb,polz,5,6,2,1,4,3,za,zb)
      a63z_265143(polq,polb,polz)=a63z(polq,polb,polz,2,6,5,1,4,3,za,zb)
      a63z_561243(polq,polb,polz)=a63z(polq,polb,polz,5,6,1,2,4,3,za,zb)
      a63z_165243(polq,polb,polz)=a63z(polq,polb,polz,1,6,5,2,4,3,za,zb)
      enddo
      enddo
      enddo

      do j=-nf,nf
      do k=-nf,nf

c---Desired formula =(Att*(A61+A61o+(A63-A63s/n+A62s+A62os)/n)
c                     +Atts*(A61s+A61os+(A63s-A63/n +A62+A62o)/n))
c                   =tamp*(lamp+lampx)
c                   +tamps*(lamps+lampsx)
c   where lamp,lampsx are the pieces without the "s" (5<->6 swap)
c     and lamps,lampx are the pieces that have the "s"
      do polq=1,2
      do polz=1,2
      do polb=1,2
        tamp=czip
        lamp=czip
        lampx=czip
        tamps=czip
        lamps=czip
        lampsx=czip
c--- Q-Q
        if ((j .gt. 0) .and. (k .gt. 0)) then
          idfac=1d0
          tamp=atreez_526143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -atreez_251643(3-polb,3-polq,polz)
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop)
          lamp=a61z_526143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_251643(3-polb,3-polq,polz)
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop)
     .        +a63z_526143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          if (j .eq. k) then
c--- statistical factor
          idfac=0.5d0
          tamps=-(atreez_625143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .         -atreez_261543(3-polb,3-polq,polz)
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop))
          lampx=-(
     .        -a63z_625143(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_625143(polq,polb,polz)/xn
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a62z_261543(3-polb,3-polq,polz)/xn
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop))
          lamps=-(a61z_625143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_261543(3-polb,3-polq,polz)
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop)
     .        +a63z_625143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop)
          lampsx=
     .        -a63z_526143(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_526143(polq,polb,polz)/xn
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a62z_251643(3-polb,3-polq,polz)/xn
     .         *(Q(k)*q1+vQ(k,polb)*v2(polz)*prop)
          endif
          if (compare) then
            msqv(j,k)=msqv(j,k)
     .      +idfac*faclo*dble(tamp*dconjg(tamp))
     .      +idfac*faclo*dble(tamps*dconjg(tamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +idfac*faclo*dble(tamp*dconjg(tamps))*(-2d0/xn)
          else
            msqv(j,k)=msqv(j,k)
     .      +idfac*fac*2d0*dble(tamp*dconjg(lamp))
     .      +idfac*fac*2d0*dble(tamps*dconjg(lamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +idfac*fac*2d0*dble(tamp*dconjg(lampx))
     .      +idfac*fac*2d0*dble(tamps*dconjg(lampsx))
          endif
c--- Qbar-Qbar
        elseif ((j .lt. 0) .and. (k .lt. 0)) then
          idfac=1d0
          tamp=atreez_162543(polq,polb,polz)
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .        -atreez_615243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
          lamp=a61z_162543(polq,polb,polz)
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .        -a61z_615243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
     .        +a63z_162543(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          if (j .eq. k) then
c--- statistical factor
          idfac=0.5d0
          tamps=-(atreez_152643(polq,polb,polz)
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .        -atreez_516243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop))
          lampx=-(
     .        -a63z_152643(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .         +a62z_152643(polq,polb,polz)/xn
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .        -a62z_516243(3-polb,3-polq,polz)/xn
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop))
          lamps=-(a61z_152643(polq,polb,polz)
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .        -a61z_516243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
     .        +a63z_152643(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop)
          lampsx=
     .         -a63z_162543(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .         +a62z_162543(polq,polb,polz)/xn
     .         *(Q(-j)*q1+vQ(-j,polq)*v2(polz)*prop)
     .         -a62z_615243(3-polb,3-polq,polz)/xn
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
          endif
          if (compare) then
            msqv(j,k)=msqv(j,k)
     .      +idfac*faclo*dble(tamp*dconjg(tamp))
     .      +idfac*faclo*dble(tamps*dconjg(tamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +idfac*faclo*dble(tamp*dconjg(tamps))*(-2d0/xn)
          else
            msqv(j,k)=msqv(j,k)
     .      +idfac*fac*2d0*dble(tamp*dconjg(lamp))
     .      +idfac*fac*2d0*dble(tamps*dconjg(lamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +idfac*fac*2d0*dble(tamp*dconjg(lampx))
     .      +idfac*fac*2d0*dble(tamps*dconjg(lampsx))
          endif
c--- Q-Qbar
        elseif ((j .gt. 0) .and. (k .lt. 0)) then
          tamp=atreez_562143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -atreez_651243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
          lamp=a61z_562143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_651243(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
     .        +a63z_562143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          if (j .eq. -k) then
          tamps=-(atreez_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -atreez_621543(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop))
          lampx=-(
     .        -a63z_265143(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_265143(polq,polb,polz)/xn
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a62z_621543(3-polb,3-polq,polz)/xn
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop))
          lamps=-(a61z_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_621543(3-polb,3-polq,polz)
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
     .        +a63z_265143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop)
          lampsx=
     .        -a63z_562143(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_562143(polq,polb,polz)/xn
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a62z_651243(3-polb,3-polq,polz)/xn
     .         *(Q(-k)*q1+vQ(-k,polb)*v2(polz)*prop)
          tampdo=atreez_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -atreez_621543(3-polb,3-polq,polz)
     .         *(Q(1)*q1+vQ(1,polb)*v2(polz)*prop)
          lampdo=a61z_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_621543(3-polb,3-polq,polz)
     .         *(Q(1)*q1+vQ(1,polb)*v2(polz)*prop)
     .        +a63z_265143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          tampup=atreez_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -atreez_621543(3-polb,3-polq,polz)
     .         *(Q(2)*q1+vQ(2,polb)*v2(polz)*prop)
          lampup=a61z_265143(polq,polb,polz)
     .         *(Q(j)*q1+vQ(j,polq)*v2(polz)*prop)
     .        -a61z_621543(3-polb,3-polq,polz)
     .         *(Q(2)*q1+vQ(2,polb)*v2(polz)*prop)
     .        +a63z_265143(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          endif
          if (compare) then
            msqv(j,k)=msqv(j,k)
     .      +faclo*dble(tamp*dconjg(tamp))
     .      +faclo*dble(tamps*dconjg(tamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +faclo*dble(tamp*dconjg(tamps))*(-2d0/xn)
          else
            msqv(j,k)=msqv(j,k)+fac*2d0*(
     .        +dble(tamp*dconjg(lamp))
     .        +dble(tamps*dconjg(lamps)))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)+fac*2d0*(
     .        +dble(tamp*dconjg(lampx))
     .        +dble(tamps*dconjg(lampsx)))
          endif
c--- add additional annihilation diagrams if necessary
          if (j .eq. -k) then
            if (jj(j) .eq. 1) then
              if (compare) then
                msqv(j,k)=msqv(j,k)+faclo*(
     .          +dble(tampup*dconjg(tampup))*dfloat(nup)
     .          +dble(tampdo*dconjg(tampdo))*dfloat(ndo-1))
              else
                msqv(j,k)=msqv(j,k)+fac*2d0*(
     .          +dble(tampup*dconjg(lampup))*dfloat(nup)
     .          +dble(tampdo*dconjg(lampdo))*dfloat(ndo-1))
              endif
            elseif (jj(j) .eq. 2) then
              if (compare) then
                msqv(j,k)=msqv(j,k)+faclo*(
     .          +dble(tampup*dconjg(tampup))*dfloat(nup-1)
     .          +dble(tampdo*dconjg(tampdo))*dfloat(ndo))
              else
                msqv(j,k)=msqv(j,k)+fac*2d0*(
     .          +dble(tampup*dconjg(lampup))*dfloat(nup-1)
     .          +dble(tampdo*dconjg(lampdo))*dfloat(ndo))
              endif
            endif
          endif
c--- Qbar-Q
        elseif ((j .lt. 0) .and. (k .gt. 0)) then
          tamp=atreez_561243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -atreez_652143(3-polb,3-polq,polz)
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop)
          lamp=a61z_561243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a61z_652143(3-polb,3-polq,polz)
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop)
     .        +a63z_561243(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          if (j .eq. -k) then
          tamps=-(atreez_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -atreez_612543(3-polb,3-polq,polz)
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop))
          lampx=-(
     .        -a63z_165243(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_165243(polq,polb,polz)/xn
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a62z_612543(3-polb,3-polq,polz)/xn
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop))
          lamps=-(a61z_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a61z_612543(3-polb,3-polq,polz)
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop)
     .        +a63z_165243(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop)
          lampsx=
     .        -a63z_561243(polq,polb,polz)/xn**2
     .         *v2(polz)/sin2w*prop
     .        +a62z_561243(polq,polb,polz)/xn
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a62z_652143(3-polb,3-polq,polz)/xn
     .         *(Q(-j)*q1+vQ(-j,polb)*v2(polz)*prop)
          tampdo=atreez_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -atreez_612543(3-polb,3-polq,polz)
     .         *(Q(1)*q1+vQ(1,polb)*v2(polz)*prop)
          lampdo=a61z_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a61z_612543(3-polb,3-polq,polz)
     .         *(Q(1)*q1+vQ(1,polb)*v2(polz)*prop)
     .        +a63z_165243(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          tampup=atreez_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -atreez_612543(3-polb,3-polq,polz)
     .         *(Q(2)*q1+vQ(2,polb)*v2(polz)*prop)
          lampup=a61z_165243(polq,polb,polz)
     .         *(Q(k)*q1+vQ(k,polq)*v2(polz)*prop)
     .        -a61z_612543(3-polb,3-polq,polz)
     .         *(Q(2)*q1+vQ(2,polb)*v2(polz)*prop)
     .        +a63z_165243(polq,polb,polz)/xn
     .         *v2(polz)/sin2w*prop
          endif
          if (compare) then
            msqv(j,k)=msqv(j,k)
     .      +faclo*dble(tamp*dconjg(tamp))
     .      +faclo*dble(tamps*dconjg(tamps))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)
     .      +faclo*dble(tamp*dconjg(tamps))*(-2d0/xn)
          else
            msqv(j,k)=msqv(j,k)+fac*2d0*(
     .        +dble(tamp*dconjg(lamp))
     .        +dble(tamps*dconjg(lamps)))
            if (polq .eq. polb) msqv(j,k)=msqv(j,k)+fac*2d0*(
     .        +dble(tamp*dconjg(lampx))
     .        +dble(tamps*dconjg(lampsx)))
          endif
c--- add additional annihilation diagrams if necessary
          if (j .eq. -k) then
            if (jj(j) .eq. -1) then
              if (compare) then
                msqv(j,k)=msqv(j,k)+faclo*(
     .          +dble(tampup*dconjg(tampup))*dfloat(nup)
     .          +dble(tampdo*dconjg(tampdo))*dfloat(ndo-1))
              else
                msqv(j,k)=msqv(j,k)+fac*2d0*(
     .          +dble(tampup*dconjg(lampup))*dfloat(nup)
     .          +dble(tampdo*dconjg(lampdo))*dfloat(ndo-1))
              endif
            elseif (jj(j) .eq. -2) then
              if (compare) then
                msqv(j,k)=msqv(j,k)+faclo*(
     .          +dble(tampup*dconjg(tampup))*dfloat(nup-1)
     .          +dble(tampdo*dconjg(tampdo))*dfloat(ndo))
              else
                msqv(j,k)=msqv(j,k)+fac*2d0*(
     .          +dble(tampup*dconjg(lampup))*dfloat(nup-1)
     .          +dble(tampdo*dconjg(lampdo))*dfloat(ndo))
              endif
            endif
          endif
        endif
      enddo
      enddo
      enddo

      enddo
      enddo

      endif

c--- write out ug virtual amplitude when checking
      if (checkvector .or. checkaxial) then
        write(6,*) 'Madloop check: ug Virt',msqv(2,0)
c      pause
      endif

************************************************************************
*     UV contributions are included here                               *
C     This is the correction to put the answer in UV renormalized      *
C     dred scheme with msbar coupling                                  *
************************************************************************
      do j=-nf,nf
      do k=-nf,nf

      do cs=0,2
      msqv(j,k)=msqv(j,k)-ason2pi*subuv(cs)*msq_cs(cs,j,k)
      enddo

      enddo
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
    subroutine terrestrial_dep(time)
!   ********************************
    implicit none
    include 'combine_foam.inc'

    ept=0.023
    fodt=0.00d+12*clo
    return
    end


text:
C******************************************************************************
C FILE: omp_bug3.f
C DESCRIPTION:
C   Run time bug
C AUTHOR: Blaise Barney  01/09/04
C LAST REVISED: 06/28/05
C******************************************************************************

      PROGRAM BUG3

      INTEGER N, I, NTHREADS, TID, SECTION, OMP_GET_NUM_THREADS,
     +        OMP_GET_THREAD_NUM
      PARAMETER (N=50)
      REAL A(N), B(N), C(N)

C     Some initializations
      DO I = 1, N
        A(I) = I * 1.0
        B(I) = A(I)
      ENDDO

!$OMP PARALLEL PRIVATE(C,I,TID,SECTION)
      TID = OMP_GET_THREAD_NUM()
      IF (TID .EQ. 0) THEN
        NTHREADS = OMP_GET_NUM_THREADS()
        PRINT *, 'Number of threads = ', NTHREADS
      END IF

C     Use barriers for clean output
!$OMP BARRIER
      PRINT *, 'Thread ', TID, ' starting...'
!$OMP BARRIER

!$OMP SECTIONS
!$OMP SECTION
      SECTION = 1
      DO I = 1, N
         C(I) = A(I) * B(I)
      ENDDO
      CALL PRINT_RESULTS(C, TID, SECTION)

!$OMP SECTION
      SECTION = 2
      DO I = 1, N
         C(I) = A(I) + B(I)
      ENDDO
      CALL PRINT_RESULTS(C, TID, SECTION)

!$OMP END SECTIONS

C     Use barrier for clean output
!$OMP BARRIER
      PRINT *, 'Thread', tid, ' exiting...'

!$OMP END PARALLEL

      END

      SUBROUTINE PRINT_RESULTS(C, TID, SECTION)
      INTEGER TID, SECTION, N, I, J
      PARAMETER (N=50)
      REAL C(N)

      J = 1
C     Use critical for clean output
!$OMP CRITICAL
      PRINT *, ' '
      PRINT *, 'Thread', TID, ' did section', SECTION
      DO I=1, N
        WRITE(*,100) C(I)
 100    FORMAT(E12.6,$)
        J = J + 1
        IF (J .EQ. 6) THEN
          PRINT *, ' '
          J = 1
        END IF
      END DO
      PRINT *, ' '
!$OMP END CRITICAL

!$OMP BARRIER
      PRINT *,'Thread', TID, ' done and synchronized'

      END SUBROUTINE PRINT_RESULTS


    Program dgecon_example

!     DGECON Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_interfaces, Only: dgecon, dgetrf, dlange
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
      Character (1), Parameter :: norm = '1'
!     .. Local Scalars ..
      Real (Kind=dp) :: anorm, rcond
      Integer :: i, info, lda, n
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), work(:)
      Integer, Allocatable :: ipiv(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon
!     .. Executable Statements ..
      Write (nout, *) 'DGECON Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      Allocate (a(lda,n), work(4*n), ipiv(n), iwork(n))

!     Read A from data file

      Read (nin, *)(a(i,1:n), i=1, n)

!     Compute norm of A

      anorm = dlange(norm, n, n, a, lda, work)

!     Factorize A
      Call dgetrf(n, n, a, lda, ipiv, info)

      Write (nout, *)
      If (info==0) Then

!       Estimate condition number

        Call dgecon(norm, n, a, lda, anorm, rcond, work, iwork, info)

        If (rcond>=epsilon(1.0E0_dp)) Then
          Write (nout, 100) 'Estimate of condition number =', 1.0E0_dp/rcond
        Else
          Write (nout, *) 'A is singular to working precision'
        End If
      Else
        Write (nout, *) 'The factor U is singular'
      End If

100   Format (1X, A, 1P, E10.2)
    End Program


text:
    module RangesTests
    use RangeUtils
    implicit none
    contains

    function RunRangesTests() result (fails)
    integer fails
    Type(TRanges) R
    double precision, allocatable :: tmp(:)
    double precision, pointer :: P(:)

    fails = 0
    !Combine set of ranges with specified approx maximum spacing in each
    call R%Add(1.d0,5.d0,4) !1, 2,3, 4, 5
    call R%Add_delta(0.5d0, 2.d0,0.5d0)
    call R%Add_delta(2.d0, 4.d0, 1.1d0) !Should be ignored as wider spacing
    call R%GetArray()
    if (any(R%points /= [0.5d0, 1.d0, 1.5d0, 2.d0, 3.d0, 4.d0, 5.d0])) then
        fails = fails + 1
        print *, 'Error in Ranges', R%Points
    end if
    !Check very close points merged
    call R%Add_delta(1.39999d0, 1.800001d0, 0.1d0)
    call R%GetArray()
    if (any(R%points==1.5d0)) then
        fails=fails+1
        print *, 'Error in ranges merge', R%points
    end if
    call R%Free()
    call R%Add(1d-3, 10.d0, 4, isLog=.true.)
    call R%Add_delta(-2.d0, 3.d0, 1.d0)
    call R%GetArray()
    if (R%Points(1)/=-2.d0 .or. abs(R%Points(4)-1d-2)>1d-10) then
        fails = fails + 1
        print *, 'Error in neg range', R%Points
    end if
    call R%Free()
    call R%Add(1d-3, 10.d0, 4, isLog=.true.)
    call R%Add_delta(1.d0, 15.d0, 1.d0)
    call R%GetArray()
    if (abs(R%points(5)-2.d0) > 1d-10 .or. abs(R%Points(2)-0.01d0)>1d-10) then
        fails = fails + 1
        print *, 'Error in log range', R%Points
    end if
    !Test ordering invariance (not generally true, but is here)
    call R%Free()
    call R%Add(1d-3, 10.d0, 4, isLog=.true.)
    call R%Add_delta(0.d0, 3.d0, 1.d0)
    call R%GetArray()
    allocate(tmp, source = R%Points)
    call R%Free()
    call R%Add_delta(0.d0, 3.d0, 1.d0)
    call R%Add(1d-3, 10.d0, 4, isLog=.true.)
    call R%GetArray()
    if (size(tmp) /= size(R%Points)) then
        fails = fails + 1
        print *, 'Error in size on range add ordering'
        print *, 'arr 1', tmp
        print *, 'arr 2 ', R%Points
    elseif (any(abs(tmp-R%points)>1d-7)) then
        fails = fails + 1
        print *, 'Error in range add ordering'
        print *, 'arr 1', tmp
        print *, 'arr 2 ', R%Points
    end if

    call R%Free()
    call R%Add_delta(1.d0, 5.d0, 1.d0)
    call R%Add(1d-3, 4.d0, 4, isLog=.true.)
    call R%GetArray()
    if (abs(R%points(5)-1)>1d-10) then
        fails = fails + 1
        print *, 'error in log left of linear'
        print *, R%points
    end if
    deallocate(tmp)
    allocate(tmp, source = R%Array())
    call R%Free()
    call R%Add(1d-3, 4.d0, 4, isLog=.true.)
    call R%Add_delta(1.d0, 5.d0, 1.d0)
    if (any(abs(R%Array() - tmp)>1d-10)) then
        fails = fails + 1
        print *, 'error in linear right of log'
        print *, R%points
    end if
    call R%Free()
    call R%Add(1d-3, 4.d0, 4, isLog=.true.)
    call R%Add_delta(3.d0, 7.d0, 1.d0)
    call R%Add_delta(2.d0, 4.3d0, 0.1d0)
    P => R%Array()
    if (abs(P(size(P)-1)-6.1d0)>1d-10) then
        fails = fails + 1
        print *, 'error in intersecting arrays'
        print *, R%Array()
    end if

    if (fails==0) print *, 'Ranges OK'
    end function

    end module

max_stars_repo_name: LBJ-Wade/forutils

max_issues_repo_name: LBJ-Wade/forutils

max_forks_repo_name: LBJ-Wade/forutils


text:
      program bs

C     Standard (?) extension found in Zebulon
C     SUn f77 and GNU g77 do not produce code with same behavior
C     Non-ANSI extension only recognized by xlf, which produces
C     code with same behavior as SUN f77

c      print *, 'Hello!'
      print *, '---\'\a\\\n''---'
      end


text:
*DECK STWAY
      SUBROUTINE STWAY (U, V, YHP, INOUT, STOWA)
C***BEGIN PROLOGUE  STWAY
C***SUBSIDIARY
C***PURPOSE  Subsidiary to BVSUP
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (STWAY-S, DSTWAY-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C  This subroutine stores (recalls) integration data in the event
C  that a restart is needed (the homogeneous solution vectors become
C  too dependent to continue)
C
C***SEE ALSO  BVSUP
C***ROUTINES CALLED  STOR1
C***COMMON BLOCKS    ML15TO, ML18JR, ML8SZ
C***REVISION HISTORY  (YYMMDD)
C   750601  DATE WRITTEN
C   890921  Realigned order of variables in certain COMMON blocks.
C           (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C***END PROLOGUE  STWAY
C
      DIMENSION U(*),V(*),YHP(*),STOWA(*)
C
      COMMON /ML8SZ/ C,XSAV,IGOFX,INHOMO,IVP,NCOMP,NFC
      COMMON /ML15TO/ PX,PWCND,TND,X,XBEG,XEND,XOT,XOP,INFO(15),ISTKOP,
     1                KNSWOT,KOP,LOTJP,MNSWOT,NSWOT
      COMMON /ML18JR/ AE,RE,TOL,NXPTS,NIC,NOPG,MXNON,NDISK,NTAPE,NEQ,
     1                INDPVT,INTEG,NPS,NTP,NEQIVP,NUMORT,NFCC,
     2                ICOCO
C
C***FIRST EXECUTABLE STATEMENT  STWAY
      IF (INOUT .EQ. 1) GO TO 100
C
C     SAVE IN STOWA ARRAY AND ISTKOP
C
      KS=NFC*NCOMP
      CALL STOR1(STOWA,U,STOWA(KS+1),V,1,0,0)
      KS=KS+NCOMP
      IF (NEQIVP .EQ. 0) GO TO 50
      DO 25 J=1,NEQIVP
      KSJ=KS+J
   25 STOWA(KSJ)=YHP(KSJ)
   50 KS=KS+NEQIVP
      STOWA(KS+1)=X
      ISTKOP=KOP
      IF (XOP .EQ. X) ISTKOP=KOP+1
      RETURN
C
C     RECALL FROM STOWA ARRAY AND ISTKOP
C
  100 KS=NFC*NCOMP
      CALL STOR1(YHP,STOWA,YHP(KS+1),STOWA(KS+1),1,0,0)
      KS=KS+NCOMP
      IF (NEQIVP .EQ. 0) GO TO 150
      DO 125 J=1,NEQIVP
      KSJ=KS+J
  125 YHP(KSJ)=STOWA(KSJ)
  150 KS=KS+NEQIVP
      X=STOWA(KS+1)
      INFO(1)=0
      KO=KOP-ISTKOP
      KOP=ISTKOP
      IF (NDISK .EQ. 0  .OR.  KO .EQ. 0) RETURN
      DO 175 K=1,KO
  175 BACKSPACE NTAPE
      RETURN
      END


text:
      subroutine testWBFVVApartComb(j1,j2,j7,j8,partSwapOk)
      implicit none
      include 'constants.f'
      include 'pid_pdg.f'
      integer j1,j2,j7,j8
      logical partSwapOk

      partSwapOk=.true.
      if (
     & (pid_pdg(1).eq.21) .or.
     & (pid_pdg(2).eq.21) .or.
     & (pid_pdg(7).eq.21) .or.
     & (pid_pdg(8).eq.21) .or.
     & (
     & (j1.eq.7 .or. j2.eq.7) .and. (
     &      (pid_pdg(7).ge.1 .and. pid_pdg(7).le.5)
     & .or. (pid_pdg(7).ge.11 .and. pid_pdg(7).le.16)
     & )
     & ) .or. (
     & (j1.eq.8 .or. j2.eq.8) .and. (
     &      (pid_pdg(8).ge.1 .and. pid_pdg(8).le.5)
     & .or. (pid_pdg(8).ge.11 .and. pid_pdg(8).le.16)
     & )
     & ) .or. (
     & (j7.eq.7 .or. j8.eq.7) .and. (
     &      (-pid_pdg(7).ge.1 .and. -pid_pdg(7).le.5)
     & .or. (-pid_pdg(7).ge.11 .and. -pid_pdg(7).le.16)
     & )
     & ) .or. (
     & (j7.eq.8 .or. j8.eq.8) .and. (
     &      (-pid_pdg(8).ge.1 .and. -pid_pdg(8).le.5)
     & .or. (-pid_pdg(8).ge.11 .and. -pid_pdg(8).le.16)
     & )
     & ) .or. (! Tests from the opposite direction
     & (j7.eq.1 .or. j8.eq.1) .and. (
     &      (pid_pdg(1).ge.1 .and. pid_pdg(1).le.5)
     & )
     & ) .or. (
     & (j7.eq.2 .or. j8.eq.2) .and. (
     &      (pid_pdg(2).ge.1 .and. pid_pdg(2).le.5)
     & )
     & ) .or. (
     & (j1.eq.1 .or. j2.eq.1) .and. (
     &      (-pid_pdg(1).ge.1 .and. -pid_pdg(1).le.5)
     & )
     & ) .or. (
     & (j1.eq.2 .or. j2.eq.2) .and. (
     &      (-pid_pdg(2).ge.1 .and. -pid_pdg(2).le.5)
     & )
     & ) ) then
      partSwapOk = .false.
      endif

      return
      end


text:
*----------------------------------------------------------------------*
      subroutine symidx_ab(idxlist_out,
     &     idxlist_in,nlist,mel,
     &     op_info,str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
*
*     given a list of indices (within bounds of ME-list mel)
*     return the corresponding indices if AB-symmetrization is
*     carried out incl. sign due to phase changes
*
*     nov 2008, andreas
*
*----------------------------------------------------------------------*

      implicit none

      include 'stdunit.h'
      include 'opdim.h'
      include 'def_orbinf.h'
      include 'def_graph.h'
      include 'mdef_operator_info.h'
      include 'def_strinf.h'
      include 'def_strmapinf.h'
      include 'ifc_memman.h'

      integer, parameter ::
     &     ntest = 00

      type(orbinf), intent(in) ::
     &     orb_info
      type(strinf), intent(in) ::
     &     str_info
      type(strmapinf), intent(in) ::
     &     strmap_info
      type(operator_info), intent(in) ::
     &     op_info
      type(me_list), intent(inout) ::
     &     mel
      integer, intent(in) ::
     &     nlist, idxlist_in(nlist)
      integer, intent(out) ::
     &     idxlist_out(nlist)

      integer ::
     &     nocc_cls, njoined, ioff_blk,
     &     ifree, nblk, nbuff, iblk,
     &     idx, ngam, ioff, lenblk
      real(8) ::
     &     fac_rel, fac_pre

      type(operator), pointer ::
     &     op

      real(8), external ::
     &     ddot
      logical, external ::
     &     list_in_bounds

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'symidx_ab')
        write(lulog,*) 'MEL:  ',trim(mel%label)
        write(lulog,*) 'idxlist(IN):'
        write(lulog,'(1x,5i6,x,5i6)') idxlist_in(1:nlist)
      end if

      idxlist_out(1:nlist) = 0

      op  => mel%op

      nocc_cls = op%n_occ_cls
      njoined  = op%njoined

      ! Number of irreps in symmetry group.
      ngam = orb_info%nsym

      ! Loop over occupation classes.
      iocc_loop: do iblk = 1, nocc_cls

        if(op%formal_blk(iblk)) cycle

        ioff_blk = (iblk-1)*njoined

        ioff = mel%off_op_occ(iblk)
        lenblk = mel%len_op_occ(iblk)

        if (.not.list_in_bounds(idxlist_in,nlist,ioff+1,ioff+lenblk))
     &       cycle

        if (ntest.ge.100)
     &       write(lulog,*) 'iblk, ioff, lenblk: ',iblk, ioff, lenblk

        call symidx_ab_blk(idxlist_out,
     &                  idxlist_in, nlist,
     &                  mel,iblk,
     &                  str_info,strmap_info,orb_info)

      end do iocc_loop

      if (ntest.ge.100) then
        write(lulog,*) 'idxlist(OUT):'
        write(lulog,'(1x,5i6,x,5i6)') idxlist_out(1:nlist)
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public