text:
      PROGRAM FDLDASHC
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C INVOKE DEMO DRIVER
C
      CALL LINEEX (IWKID)
C
C     DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
      STOP
      END
C
      SUBROUTINE LINEEX (IWKID)
C
C PURPOSE                To provide a simple demonstration of the
C                        different line drawing techniques.
C
C USAGE                  CALL LINEEX (IWKID)
C
C ARGUMENTS
C
C ON INPUT               IWKID
C                          Workstation id
C
C LANGUAGE               FORTRAN
C
C HISTORY                Written  by members of the
C                        Scientific Computing Division of NCAR,
C                        Boulder Colorado
C
C PORTABILITY            FORTRAN 77
C
C NOTE                   The call to GOPWK will have to be modified
C                        when using a non-NCAR GKS package.  The third
C                        argument must be the workstation type for WISS.
C
C
      CHARACTER STR*60
C
C  Data for the graphical objects.
C
      DATA    TWOPI/6.283185/
      DATA    X0P, Y0P, RP, NPTP/ 0.500, 0.5, 0.45 ,  16/
C
C Declare the constant for converting from degrees to radians.
C
      DATA DTR / .017453292519943 /
C
C Establish the viewport and window.
C
      CALL SET(.01,.99,0.01,.99,0.,1.,0.,1.,1)
C
C Turn buffering off
C
      CALL SETUSV('PB',2)
C
C Set up a color table
C
C White background
C
      CALL GSCR (IWKID,0,1.,1.,1.)
C
C Black foreground
C
      CALL GSCR (IWKID,1,0.,0.,0.)
C
C Red
C
      CALL GSCR (IWKID,2,1.,0.,0.)
C
C Green
C
      CALL GSCR (IWKID,3,0.,1.,0.)
C
C Blue
C
      CALL GSCR (IWKID,4,0.,0.,1.)
C
C  Create a polygonal fan
C
      DTHETA = TWOPI/NPTP
      DO 10 I=1,6
         ANG = DTHETA*REAL(I) + .19625
         XC = RP*COS(ANG)
         YC = RP*SIN(ANG)
C
C Set the line color
C
         CALL GSPLCI (MOD(I,4)+1)
C
C Set line width
C
         CALL GSLWSC(2.)

         IF (I .EQ. 1) THEN
C
C  Solid Line (the default)
C
            STR = '$$$$Solid$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL GSLN(1)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 2) THEN
C
C Dashed line
C
            STR = '$$$$Dashed$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL GSLN(2)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 3) THEN
C
C Dotted line
C
            STR = '$$$$Dotted$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL GSLN(3)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 4) THEN
C
C Dashed dotted line
C
            STR = '$$$$Dotted''Dashed$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL GSLN(4)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 5) THEN
C
C Don't do anything different here
C Color is changed at beginning of loop
C
            STR = '$$$$Color''Any$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 6) THEN
C
C Increase the line width
C
            STR = '$$$$Width''Any$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
            CALL DASHDC(STR,20,20)
            CALL GSLWSC(10.)
            CALL LINED(X0P,Y0P,X0P+XC,Y0P+YC)
         ELSEIF (I .EQ. 7) THEN
         ELSEIF (I .EQ. 8) THEN
         ENDIF
C
C Reset the line type to a thin black solid line
C
         CALL GSLN(1)
         CALL GSLWSC(2.)
         CALL GSPLCI (1)
 10   CONTINUE
C
C Create and Plot 2 Sine Curves
C
C Set line width
C
      CALL GSLWSC(4.)
C
C Set the dash pattern
C
      STR = '$''$Draw''Any''Curve$''$''$''$''$''$''$'''
      CALL DASHDC(STR,15,20)
C
C Move plotter pen
C
      CALL FRSTD(0.,.25)
C
C Compute the curve coordinates
C
      DO 20 I=1,360
         Y = (SIN(REAL(I)*(TWOPI/360.)) * .25) + .25
         CALL VECTD(REAL(I)/360.,Y)
 20   CONTINUE
C
C Set line width
C
      CALL GSLWSC(2.)
C
C  Set the dash pattern
C
      STR = '$''$$''$$$''$$$$ Any''Pattern'
      CALL DASHDC(STR,20,20)
C
C Set the line color to green
C
      CALL GSPLCI (3)
C
C Move plotter pen
C
      CALL FRSTD(0.,.125)
C
C Compute the curve coordinates
C
      DO 30 I=1,360
         Y = (SIN(REAL(I)*(TWOPI/360.)) * .125) + .125
         CALL VECTD(REAL(I)/720.,Y)
 30   CONTINUE
C
C Create and plot a spiral curve
C
      CALL SET (.4,.9,.0,.55,-1.,1.,-1.,1.,1)
      STR = '$$$$$$$$$$$$$$$$$$$$$$$$$$Any''Shape'
      CALL DASHDC(STR,80,20)
C
C Move plotter pen
C
      RAD=.001
      XCD=.25+.5*RAD*COS(0.)
      YCD=.25+.5*RAD*SIN(0.)
      CALL FRSTD(XCD,YCD)
C
C Set the line color to red
C
      CALL GSPLCI (2)
      DO 40 ING=1,1500
         RAD=REAL(ING)/1000.
         ANG=DTR*REAL(ING-1)
         XCD=.25+.5*RAD*COS(ANG)
         YCD=.25+.5*RAD*SIN(ANG)
         CALL VECTD(XCD, YCD)
 40   CONTINUE
      CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1)
C
C Set the line color to black
C
      CALL GSPLCI (1)
C
C  Create a background perimeter
C
      CALL FRSTPT( 0.0, 0.0)
      CALL VECTOR( 1.0, 0.0)
      CALL VECTOR( 1.0, 1.0)
      CALL VECTOR( 0.0, 1.0)
      CALL VECTOR( 0.0, 0.0)
C
C  Label the plot
C
C     CALL PLCHLQ(0.5,0.91,'',25.,0.,0.)
      CALL FRAME
C
      END


text:
MODULE VCA_BATH_FUNCTIONS
  USE SF_CONSTANTS, only: zero
  USE SF_IOTOOLS, only:free_unit,reg,file_length,txtfy
  USE SF_LINALG, only: eye,inv,zeye,inv_her
  USE VCA_INPUT_VARS
  USE VCA_VARS_GLOBAL
  USE VCA_BATH_SETUP
  USE VCA_AUX_FUNX
  implicit none

  private

  interface delta_bath_freq
    module procedure delta_bath_freq_main
  endinterface delta_bath_freq

  !##################################################################
  !
  !\DELTA HYBRIDIZATION FUNCTION MATSUBARA
  !
  !##################################################################
  !NORMAL
  interface delta_bath_mats
     module procedure delta_bath_mats_main
     module procedure delta_bath_mats_main_
  end interface delta_bath_mats

  !##################################################################
  !
  !\DELTA HYBRIDIZATION FUNCTION REAL
  !
  !##################################################################
  !NORMAL
  interface delta_bath_real
     module procedure delta_bath_real_main
     module procedure delta_bath_real_main_
  end interface delta_bath_real
  !

  !##################################################################
  !
  !NON-INTERACTING GREEN'S FUNCTION MATSUBARA
  !
  !##################################################################
  !NORMAL
  interface g0and_bath_mats
     module procedure g0and_bath_mats_main
     module procedure g0and_bath_mats_main_
  end interface g0and_bath_mats
  !

  !##################################################################
  !
  !INVERSE NON-INTERACTING GREEN'S FUNCTION MATSUBARA
  !
  !##################################################################
  !NORMAL
  interface invg0_bath_mats
     module procedure invg0_bath_mats_main
     module procedure invg0_bath_mats_main_
  end interface invg0_bath_mats

  !##################################################################
  !
  !NON-INTERACTING GREEN'S FUNCTION REAL-AXIS
  !
  !##################################################################
  !NORMAL
  interface g0and_bath_real
     module procedure g0and_bath_real_main
     module procedure g0and_bath_real_main_
  end interface g0and_bath_real
  !

  !##################################################################
  !
  !INVERSE NON-INTERACTING GREEN'S FUNCTION REAL-AXIS
  !
  !##################################################################
  !NORMAL
  interface invg0_bath_real
     module procedure invg0_bath_real_main
     module procedure invg0_bath_real_main_
  end interface invg0_bath_real
  !
  public :: delta_bath_freq

  public :: delta_bath_mats
  public :: g0and_bath_mats
  public :: invg0_bath_mats
  !
  public :: delta_bath_real
  public :: g0and_bath_real
  public :: invg0_bath_real

  integer :: ilat,jlat,iorb,jorb,ispin,jspin

contains

  !##################################################################
  !
  !     DELTA FUNCTIONS
  !     G0 FUNCTIONS
  !     G0^{-1} FUNCTIONS
  !     at a point x from type(effective_bath) :: vca_bath
  !
  !##################################################################
  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  Delta functions for specific frequency:
  !+-----------------------------------------------------------------------------+!
  !NORMAL:
  function delta_bath_freq_main(x,vca_bath_) result(Delta)
    complex(8),intent(in)                                         :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb)         :: Delta
    integer                                                       :: i,ih,L
    integer                                                       :: ibath
    integer                                                       :: io,jo
    real(8),dimension(Nbath)                                      :: eps,vps
    real(8),dimension(Nlat,Norb,Nbath)                            :: vops
    !
    Delta=zero
    !
    !
    if(.not.vca_bath_%status)return
    !\Delta_{aa} = \sum_k [ V_{a}(k) * V_{a}(k)/(iw_n - E_{a}(k)) ]
    do ilat=1,Nlat
       do iorb=1,Norb
          do ispin=1,Nspin
             eps = vca_bath_%e(ilat,ispin,iorb,1:Nbath)
             vps = vca_bath_%v(ilat,ispin,iorb,1:Nbath)
             Delta(ilat,ilat,ispin,ispin,iorb,iorb) = sum( vps(:)*vps(:)/(x - eps(:)+XMU) )
          enddo
       enddo
    enddo
  end function delta_bath_freq_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  Delta functions on the Matsubara axis:
  !+-----------------------------------------------------------------------------+!
  !NORMAL:
  function delta_bath_mats_main(x,vca_bath_) result(Delta)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: Delta
    integer                                                       :: i,ih,L
    integer                                                       :: ibath
    integer                                                       :: io,jo
    real(8),dimension(Nbath)                                      :: eps,vps
    real(8),dimension(Nlat,Norb,Nbath)                            :: vops
    !
    Delta=zero
    !
    L = size(x)
    !
    if(.not.vca_bath_%status)return
    !\Delta_{aa} = \sum_k [ V_{a}(k) * V_{a}(k)/(iw_n - E_{a}(k)) ]
    do ilat=1,Nlat
       do iorb=1,Norb
          do ispin=1,Nspin
             eps = vca_bath_%e(ilat,ispin,iorb,1:Nbath)
             vps = vca_bath_%v(ilat,ispin,iorb,1:Nbath)
             do i=1,L
                Delta(ilat,ilat,ispin,ispin,iorb,iorb,i) = sum( vps(:)*vps(:)/(x(i) - eps(:) + XMU) )
             enddo
          enddo
       enddo
    enddo
  end function delta_bath_mats_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  Delta functions on the Real axis:
  !+-----------------------------------------------------------------------------+!
  function delta_bath_real_main(x,vca_bath_) result(Delta)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: Delta
    integer                                                       :: i,ih,L
    integer                                                       :: ibath
    integer                                                       :: io,jo
    real(8),dimension(Nbath)                                      :: eps,vps
    real(8),dimension(Nlat,Norb,Nbath)                            :: vops
    !
    !
    Delta=zero
    !
    L = size(x)
    !
    if(.not.vca_bath_%status)return
    !\Delta_{aa} = \sum_k [ V_{a}(k) * V_{a}(k)/(w+i\h - E_{a}(k)) ]
    do ilat=1,Nlat
       do iorb=1,Norb
          do ispin=1,Nspin
             eps = vca_bath_%e(ilat,ispin,iorb,1:Nbath)
             vps = vca_bath_%v(ilat,ispin,iorb,1:Nbath)
             do i=1,L
                Delta(ilat,ilat,ispin,ispin,iorb,iorb,i) = sum( vps(:)*vps(:)/(x(i) - eps(:) +XMU) )
             enddo
          enddo
       enddo
    enddo
    !
  end function delta_bath_real_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  G0 non-interacting Green's functions on the Matsubara axis
  !+-----------------------------------------------------------------------------+!
  !NORMAL:
  function g0and_bath_mats_main(x,vca_bath_) result(G0and)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and,Delta
    integer                                                       :: iorb,jorb,ispin,jspin,io,jo,Nso,i,L
    real(8),dimension(size(x))                                    :: det
    complex(8),dimension(size(x))                                 :: fg,ff
    complex(8),dimension(:,:),allocatable                         :: fgorb,zeta
    !
    G0and = zero
    !
    L=size(x)
    !
    !
    Delta = delta_bath_mats(x,vca_bath_)
    do ilat=1,Nlat
       do ispin=1,Nspin
          do iorb=1,Norb
             fg(:) = x(:) + xmu - impHloc(ilat,ilat,ispin,ispin,iorb,iorb) - Delta(ilat,ilat,ispin,ispin,iorb,iorb,:)
             G0and(ilat,ilat,ispin,ispin,iorb,iorb,:) = one/fg(:)
          enddo
       enddo
    enddo
    !
  end function g0and_bath_mats_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  G0 non-interacting Green's functions on the real-axis
  !+-----------------------------------------------------------------------------+!
  !NORMAL:
  function g0and_bath_real_main(x,vca_bath_) result(G0and)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and,Delta
    integer                                                       :: iorb,jorb,ispin,jspin,io,jo,Nso,i,L
    complex(8),dimension(size(x))                                 :: det,fg,ff
    complex(8),dimension(:,:),allocatable                         :: fgorb,zeta
    !
    G0and = zero
    !
    L = size(x)
    !
    Delta = delta_bath_real(x,vca_bath_)
    do ilat=1,Nlat
       do ispin=1,Nspin
          do iorb=1,Norb
             fg(:)    = x(:) + xmu - impHloc(ilat,ilat,ispin,ispin,iorb,iorb) - Delta(ilat,ilat,ispin,ispin,iorb,iorb,:)
             G0and(ilat,ilat,ispin,ispin,iorb,iorb,:) = one/fg(:)
          enddo
       enddo
    enddo
    !
  end function g0and_bath_real_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  G0^{-1} non-interacting Green's functions on the Matsubara axis:
  !+-----------------------------------------------------------------------------+!
  !NORMAL:
  function invg0_bath_mats_main(x,vca_bath_) result(G0and)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and,Delta
    integer                                                       :: i,iorb,jorb,ispin,jspin,io,jo,Nso,L
    complex(8),dimension(:,:),allocatable                         :: zeta!,fgorb
    !
    G0and = zero
    !
    L=size(x)
    !
    Delta = delta_bath_mats(x,vca_bath_)
    do ilat=1,Nlat
       do ispin=1,Nspin
          do iorb=1,Norb
             G0and(ilat,ilat,ispin,ispin,iorb,iorb,:) = x(:) + xmu - impHloc(ilat,ilat,ispin,ispin,iorb,iorb) - Delta(ilat,ilat,ispin,ispin,iorb,iorb,:)
          enddo
       enddo
    enddo
  end function invg0_bath_mats_main

  !+-----------------------------------------------------------------------------+!
  !PURPOSE:  G0 non-interacting Green's functions on the real-axis:
  !+-----------------------------------------------------------------------------+!
  function invg0_bath_real_main(x,vca_bath_) result(G0and)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and,Delta
    integer                                                       :: iorb,jorb,ispin,jspin,io,jo,Nso,i,L
    complex(8),dimension(:,:),allocatable                         :: zeta!,fgorb
    !
    G0and = zero
    !
    L = size(x)
    !
    Delta = delta_bath_real(x,vca_bath_)
    do ilat=1,Nlat
       do ispin=1,Nspin
          do iorb=1,Norb
             G0and(ilat,ilat,ispin,ispin,iorb,iorb,:) =  x(:) + xmu - impHloc(ilat,ilat,ispin,ispin,iorb,iorb) - Delta(ilat,ilat,ispin,ispin,iorb,iorb,:)
          enddo
       enddo
    enddo
    !    !
  end function invg0_bath_real_main

  !##################################################################
  !
  !     DELTA FUNCTIONS
  !     G0 FUNCTIONS
  !     G0^{-1} FUNCTIONS
  !     at a point x from real(8),dimension(:) :: bath_array
  !
  !##################################################################
  function delta_bath_mats_main_(x,bath_) result(Delta)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: Delta
    real(8),dimension(:)                                          :: bath_
    logical                                                       :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "delta_bath_mats_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    Delta = delta_bath_mats_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function delta_bath_mats_main_

  function delta_bath_real_main_(x,bath_) result(Delta)
    complex(8),dimension(:),intent(in)                            :: x
    type(effective_bath)                                          :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: Delta
    real(8),dimension(:)                                          :: bath_
    logical                                                       :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "delta_bath_real_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    Delta = delta_bath_real_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function delta_bath_real_main_

  function g0and_bath_mats_main_(x,bath_) result(G0and)
    complex(8),dimension(:),intent(in)                  :: x
    type(effective_bath)                                :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and
    real(8),dimension(:)                                :: bath_
    logical                                             :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "g0and_bath_mats_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    G0and = g0and_bath_mats_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function g0and_bath_mats_main_

  function g0and_bath_real_main_(x,bath_) result(G0and)
    complex(8),dimension(:),intent(in)                  :: x
    type(effective_bath)                                :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and
    real(8),dimension(:)                                :: bath_
    logical                                             :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "g0and_bath_real_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    G0and = g0and_bath_real_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function g0and_bath_real_main_

  function invg0_bath_mats_main_(x,bath_) result(G0and)
    complex(8),dimension(:),intent(in)                  :: x
    type(effective_bath)                                :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and
    real(8),dimension(:)                                :: bath_
    logical                                             :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "invg0_bath_mats_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    G0and = invg0_bath_mats_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function invg0_bath_mats_main_

  function invg0_bath_real_main_(x,bath_) result(G0and)
    complex(8),dimension(:),intent(in)                  :: x
    type(effective_bath)                                :: vca_bath_
    complex(8),dimension(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(x)) :: G0and
    real(8),dimension(:)                                :: bath_
    logical                                             :: check
    check= check_bath_dimension(bath_)
    if(.not.check)stop "invg0_bath_real_main_ error: wrong bath dimensions"
    call vca_allocate_bath(vca_bath_)
    call vca_set_bath(bath_,vca_bath_)
    G0and = invg0_bath_real_main(x,vca_bath_)
    call vca_deallocate_bath(vca_bath_)
  end function invg0_bath_real_main_

END MODULE VCA_BATH_FUNCTIONS

max_stars_repo_name: QcmPlab/VCA

max_issues_repo_name: QcmPlab/VCA

max_forks_repo_name: QcmPlab/VCA


text:
      PROGRAM DRIVER

      USE GET_ENV_VARS

      USE CSQY_PARAMETERS
      USE CSQY_REFER_DATA

      IMPLICIT NONE

      CHARACTER(  5 )    :: SPLIT_OUTPUT
      CHARACTER( 13 )    :: TYPE_OUTPUTS  = 'SPLIT_OUTPUT'

      INTERFACE
       SUBROUTINE CONVERT_CASE ( BUFFER, UPPER )
           CHARACTER(LEN= *), INTENT( INOUT ) :: BUFFER
           LOGICAL,           INTENT( IN    ) :: UPPER
       END SUBROUTINE CONVERT_CASE
      END INTERFACE

!Bin 01 177.5-202.5 nm
!Bin 02 177.5-202.5 nm
!Bin 03 177.5-202.5 nm
!Bin 04 177.5-202.5 nm
!Bin 05 202.5-206.5 nm
!Bin 06 206.5-209.5 nm
!Bin 07 209.5-212.5 nm
!Bin 08 212.5-215.5 nm
!Bin 09 233.0-275.5 nm
!Bin 10 275.5-286.5 nm
!Bin 11 286.5-291.0 nm
!Bin 12 291.0-298.3 nm
!Bin 13 298.3-307.5 nm
!Bin 14 307.5-312.5 nm
!Bin 15 312.5-320.3 nm
!Bin 16 320.3-345.0 nm
!Bin 17 345.0-412.5 nm
!Bin 18 412.5-850.0 nm

         CALL VALUE_NAME ( TYPE_OUTPUTS,  SPLIT_OUTPUT)

         CALL CONVERT_CASE( SPLIT_OUTPUT, .TRUE.)

         IF( SPLIT_OUTPUT(1:1) .EQ. 'T' .OR. SPLIT_OUTPUT(1:1) .EQ. 'Y' )THEN
             SPLIT_OUTPUTS = .TRUE.
             WRITE(6,'(A)')'Environment Variable  SPLIT_OUTPUT set to '
     &       // TRIM( SPLIT_OUTPUT ) // '. CSQY and Optics Data in two'
     &      //  ' separate files'
         ELSE IF(  SPLIT_OUTPUT(1:1) .EQ. 'F' .OR. SPLIT_OUTPUT(1:1) .EQ. 'N' )THEN
             SPLIT_OUTPUTS = .FALSE.
             WRITE(6,'(A)')'Environment Variable  SPLIT_OUTPUT set to '
     &      // TRIM( SPLIT_OUTPUT ) // '. CSQY and Optics Data in one file'
         ELSE
             WRITE(6,' (A)')'Environment Variable  SPLIT_OUTPUT set to '
     &       // TRIM( SPLIT_OUTPUT ) // ' and must equal T, Y, F, or N.'
     &       // ' Using default value of T'
             SPLIT_OUTPUTS = .TRUE.
         END IF

         CALL PROCESS_CSQYS ( )

         IF(SPLIT_OUTPUTS)CALL WRT_OPTICS_DATA( )

         CLOSE( JTABLE_UNIT )

        WRITE(6, *)'NORMAL STOP'
      STOP
      END

ext: f

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
! moved here from amg4psblas-extension
!
!
!                             AMG4PSBLAS  Extensions
!
!    (C) Copyright 2019
!
!                        Salvatore Filippone  Cranfield University
!        Pasqua D'Ambra         IAC-CNR, Naples, IT
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
! File: amg_s_krm_solver_impl.f90
!
!  This is the implementation file corresponding to amg_s_krm_solver_mod.
!
!
subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)

  use psb_base_mod
  use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_bld

  Implicit None

  ! Arguments
  type(psb_sspmat_type), intent(in), target        :: a
  Type(psb_desc_type), Intent(inout)                  :: desc_a
  class(amg_s_krm_solver_type), intent(inout)         :: sv
  integer(psb_ipk_), intent(out)                      :: info
  type(psb_sspmat_type), intent(in), target, optional :: b
  class(psb_s_base_sparse_mat), intent(in), optional  :: amold
  class(psb_s_base_vect_type), intent(in), optional   :: vmold
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
  integer(psb_lpk_) :: lnr
  integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level
  type(psb_ctxt_type) :: ctxt, l_ctxt
  character(len=20)   :: name='@Z@_krm_solver_bld', ch_err

  info=psb_success_
  call psb_erractionsave(err_act)
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()
  ctxt       = desc_a%get_context()
  call psb_info(ctxt, me, np)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' start'

  call sv%prec%free(info)

  if (sv%global) then
    call sv%prec%init(ctxt,sv%kprec,info)
    if (sv%i_sub_solve>0) then
      call sv%prec%set('sub_solve',sv%i_sub_solve,info)
    else
      call sv%prec%set('sub_solve',sv%sub_solve,info)
    end if
    call sv%prec%set('sub_fillin',sv%fillin,info)
    call sv%prec%hierarchy_build(a,desc_a,info)
    call sv%prec%smoothers_build(a,desc_a,info,amold=amold,vmold=vmold)
    sv%a => a
  else
    call psb_init(l_ctxt,np=1_psb_ipk_,basectxt=ctxt,ids=(/me/))
    n_row = desc_a%get_local_rows()
    lnr = n_row
    call psb_cdall(l_ctxt,sv%desc_local,info,mg=lnr,repl=.true.)
    call psb_cdasb(sv%desc_local,info)
    call sv%prec%init(l_ctxt,sv%kprec,info)
    !This is a gigantic kludge, to be revised
    if (sv%i_sub_solve>0) then
      call sv%prec%set('sub_solve',sv%i_sub_solve,info)
    else
      call sv%prec%set('sub_solve',sv%sub_solve,info)
    end if
    call sv%prec%set('sub_fillin',sv%fillin,info)
    call a%csclip(sv%a_local,info,jmax=n_row)
    call sv%prec%hierarchy_build(sv%a_local,sv%desc_local,info)
    call sv%prec%smoothers_build(sv%a_local,sv%desc_local,info,amold=amold,vmold=vmold)
    call psb_geall(sv%x_local,sv%desc_local,info)
    call sv%x_local%zero()
    call psb_geasb(sv%x_local,sv%desc_local,info,mold=vmold)
    call psb_geall(sv%z_local,sv%desc_local,info)
    call sv%z_local%zero()
    call psb_geasb(sv%z_local,sv%desc_local,info,mold=vmold)
  end if

  if(info /= psb_success_) then
    info=psb_err_from_subroutine_
    ch_err='inner precbld'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' end'

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return
end subroutine amg_s_krm_solver_bld

subroutine amg_s_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
     & trans,work,wv,info,init,initu)

  use psb_base_mod
  use psb_krylov_mod
  use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_apply_vect

  Implicit None
  type(psb_desc_type), intent(in)             :: desc_data
  class(amg_s_krm_solver_type), intent(inout) :: sv
  type(psb_s_vect_type),intent(inout)         :: x
  type(psb_s_vect_type),intent(inout)         :: y
  real(psb_spk_),intent(in)                    :: alpha,beta
  character(len=1),intent(in)                   :: trans
  real(psb_spk_),target, intent(inout)         :: work(:)
  type(psb_s_vect_type),intent(inout)         :: wv(:)
  integer(psb_ipk_), intent(out)                :: info
  character, intent(in), optional                :: init
  type(psb_s_vect_type),intent(inout), optional   :: initu

  type(psb_s_vect_type)   :: z
  integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level
  type(psb_ctxt_type) :: ctxt
  character(len=20) :: name='@Z@_krm_solver_apply_v', ch_err

  info=psb_success_
  call psb_erractionsave(err_act)

  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()
  ctxt       = desc_data%get_context()
  call psb_info(ctxt, me, np)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' start'
  if (sv%global) then
    call psb_geall(z,desc_data,info)
    call z%zero()
    call psb_geasb(z,desc_data,info,mold=x%v)

    if (debug_level >= psb_debug_outer_) &
         & write(debug_unit,*) me,' ',trim(name),' norms ',&
         & psb_genrm2(x,desc_data,info),psb_genrm2(y,desc_data,info),&
         & psb_genrm2(z,desc_data,info)

    call psb_krylov(sv%method,sv%a,sv%prec,x,z,sv%eps,&
         & desc_data,info,itmax=sv%itmax,itrace=sv%itrace,&
         & istop=sv%istopc,irst=sv%irst)
!!$  call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work)
    call psb_geaxpby(alpha,z,beta,y,desc_data,info)
  else
    call psb_geaxpby(sone,x,szero,sv%x_local,sv%desc_local,info)
    call psb_krylov(sv%method,sv%a_local,sv%prec,sv%x_local,sv%z_local,sv%eps,&
         & sv%desc_local,info,itmax=sv%itmax,itrace=sv%itrace,&
         & istop=sv%istopc,irst=sv%irst)
    call psb_geaxpby(alpha,sv%z_local,beta,y,sv%desc_local,info)
  end if

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return
end subroutine amg_s_krm_solver_apply_vect

subroutine amg_s_krm_solver_apply(alpha,sv,x,beta,y,desc_data,&
     & trans,work,info,init,initu)
  use psb_base_mod
  use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_apply
  implicit none
  type(psb_desc_type), intent(in)      :: desc_data
  class(amg_s_krm_solver_type), intent(inout) :: sv
  real(psb_spk_),intent(inout)         :: x(:)
  real(psb_spk_),intent(inout)         :: y(:)
  real(psb_spk_),intent(in)            :: alpha,beta
  character(len=1),intent(in)           :: trans
  real(psb_spk_),target, intent(inout) :: work(:)
  integer(psb_ipk_), intent(out)        :: info
  character, intent(in), optional       :: init
  real(psb_spk_),intent(inout), optional :: initu(:)
  real(psb_spk_), allocatable    :: z(:)
  integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level
  type(psb_ctxt_type) :: ctxt
  character(len=20) :: name='@Z@_krm_solver_apply', ch_err

  info=psb_success_
  call psb_erractionsave(err_act)
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()
  ctxt       = desc_data%get_context()
  call psb_info(ctxt, me, np)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),' start'

  call psb_geasb(z,desc_data,info,scratch=.true.)

  call sv%prec%apply(x,z,desc_data,info,trans=trans,work=work)

  call psb_geaxpby(alpha,z,beta,y,desc_data,info)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return
end subroutine amg_s_krm_solver_apply


            program demo_string_to_values
             use M_strings, only : string_to_values
             implicit none
             character(len=80)  :: s=' 10 20e3;3.45 -400.3e-2;1234; 5678 '
             integer,parameter  :: isz=10
             real               :: array(isz)
             integer            :: inums, ierr, ii

             call string_to_values(s,10,array,inums,' ;',ierr)
             call reportit()

             call string_to_values('10;2.3;3.1416',isz,array,inums,' ;',ierr)
             call reportit()

             contains
                subroutine reportit()
                   write(*,*)'string_to_values:'
                   write(*,*)'input string.............',trim(s)
                   write(*,*)'number of values found...',inums
                   write(*,*)'values...................',(array(ii),ii=1,inums)
                end subroutine reportit
            end program demo_string_to_values


text:
program persist
    implicit none
    ! Explicit types for blas calls
    integer, parameter :: i32 = 4
    integer, parameter :: i64 = 8
    integer, parameter :: f32 = kind(1.e0)
    integer, parameter :: f64 = kind(1.d0)
    real(f64) :: ddot

    ! Iteration variables
    integer(i32) i, j

    ! Physical constants
    integer(i32),  parameter :: shape = 43_i32
    real(f64),     parameter :: scale = 1.91163_f64

    ! Specified constants
    real(f64)    :: kon, koff, sigma
    integer(i32) :: m, n, d

    ! Program variables
    real(f64),    allocatable :: A(:,:), x(:), y(:)
    integer(i32), allocatable :: ipiv(:)
    real(f64)                 :: z

    ! Parse command-line input
    character(len=128) :: arg
    if (command_argument_count() < 4) then
        print *, "Usage: ./persist mass minsize beta gamma"
        call exit()
    end if
    call get_command_argument(1, arg)
    read(arg,*) m
    call get_command_argument(2, arg)
    read(arg,*) n
    call get_command_argument(3, arg)
    read(arg,*) kon
    call get_command_argument(4, arg)
    read(arg,*) koff

    ! Program constants
    d = m-n+1
    sigma = 2.0_f64 * kon / koff

    ! Allocate and initialize the memory for the problem
    allocate(A(n+2, n:m), x(n:m), y(n:m), ipiv(n:m))
    ! super-diagonals
    A(2:n, n:m) = 1.0_f64
    ! diagonal
    A(n+1, n:m) = (/ (-sigma*(m-i) - (i-1), i=n, m) /)
    ! sub-diagonal
    A(n+2, n:m) = (/ (sigma*(m-i), i=n, m) /)
    ! summation vector
    x(n:(2*n-1)) = 0.0_f64
    x((2*n):m)   = (/ (i-2*n+1, i=2*n, m) /)
    ! unit vector
    y(n) = 1.0_f64
    y((n+1):m) = 0.0_f64

    ! LU-factorize the matrix operator
    call dgbtrf(d, d, 1, n-1, A, n+2, ipiv, i)

    ! update y <- inv(A)*y
    call dgbtrs('N', d, 1, n-1, 1, A, n+2, ipiv, y, d, i)

    z = ddot(d, x, 1, y, 1)

    ! update A again (maybe could update factors directly if too slow)
    ! super-diagonals
    A(2:n, n:m) = -koff*scale*1.0_f64
    ! diagonal
    A(n+1, n:m) = (/ (1 + (koff*scale)*(sigma*(m-i) + (i-1)), i=n, m) /)
    ! sub-diagonal
    A(n+2, n:m) = (/ (-(koff*scale)*sigma*(m-i), i=n, m) /)

    ! LU-factorize the matrix operator
    call dgbtrf(d, d, 1, n-1, A, n+2, ipiv, i)

    do j = 1, shape
        call dgbtrs('N', d, 1, n-1, 1, A, n+2, ipiv, y, d, i)
    end do

    print *, m, n, kon, koff, ddot(d, x, 1, y, 1)-z

    ! Cleanup the problem
    deallocate(A, x, y, ipiv)

end program persist

max_stars_repo_name: jasondark/dissertation

max_issues_repo_name: jasondark/dissertation

max_forks_repo_name: jasondark/dissertation


subroutine init()
use uio_bulk_module
implicit none
call uio_init(progrm='pybold')
end subroutine init
!
subroutine pyread(uio_data_reg, modelfile, parfile, imodel)
use uio_reader_module
use rhd_dq_module
implicit none
external uio_data_reg
character*(*),              intent(in)  :: modelfile
character*(*),              intent(in)  :: parfile
integer,                    intent(in)  :: imodel
integer                                 :: imodel0, nmodel
!
record => uio_data_reg
call uio_lbl_count(modelfile, startlabel='dataset', endlabel='enddataset', &
                   nlabel=nmodel)
if (trim(parfile) /= '' .and. nmodel>0) then
        use_dq = .True.
else
        use_dq = .False.
endif
if (imodel<0) then
        imodel0=nmodel+imodel+1
else
        imodel0=imodel
endif
!print*, 'nmodel: ', nmodel
if (use_dq) then
        call rhd_dq_init(parfile, modelfile, nmodel, imodel0, 6)
endif
if (nmodel>0) then
        call uio_struct_rd(modelfile, startlabel='dataset', &
                           endlabel='enddataset', nlabel=imodel0, &
                           substartlabel='box', subendlabel='endbox', &
                           headread=.true.)
else
        call uio_struct_rd(modelfile, endlabel='enddataset')
endif
if (use_dq) then
        call rhd_dq_end_init()
endif
end subroutine pyread
!
subroutine open_model(modelfile, parfile)
use uio_reader_module
implicit none
character*(*),              intent(in)  :: modelfile
character*(*),              intent(in)  :: parfile
!
r_modelfile = modelfile
r_parfile = parfile
call uio_openrd(r_nc, r_modelfile)
call uio_lbl_count(r_modelfile, r_nc, &
                   startlabel='dataset', endlabel='enddataset', &
                   nlabel=r_nmodel, maxlabel=1)
if (trim(r_parfile) /= '' .and. r_nmodel>0) then
        use_dq = .True.
else
        use_dq = .False.
endif
call uio_chpos(r_nc, 'rewind')
r_imodel = 0
if (r_nmodel>0) r_nmodel = -1
end subroutine open_model
!
subroutine count_model()
use uio_reader_module
implicit none
integer                                 :: pos
!
call uio_chpos(r_nc, 'get', posout=pos)
call uio_chpos(r_nc, 'rewind')
call uio_lbl_count(r_modelfile, r_nc, &
                   startlabel='dataset', endlabel='enddataset', nlabel=r_nmodel)
call uio_chpos(r_nc, 'goto', pos)
end subroutine count_model
!
subroutine close_model()
use uio_reader_module
implicit none
!
call uio_closrd(r_nc)
end subroutine close_model
!
subroutine read_header(uio_data_reg)
use uio_reader_module
implicit none
external uio_data_reg
integer                                 :: pos
!
record => uio_data_reg
if (r_nmodel/=0) then
        call uio_chpos(r_nc, 'get', posout=pos)
        call uio_chpos(r_nc, 'rewind')
        call uio_struct_rd(r_modelfile, channel=r_nc, headread=.true.)
        call uio_chpos(r_nc, 'goto', pos)
endif
end subroutine read_header
!
subroutine read_next(uio_data_reg)
use uio_reader_module
implicit none
external uio_data_reg
!
record => uio_data_reg
if (r_imodel < r_nmodel .or. r_nmodel < 0) then
        r_imodel = r_imodel+1
        if (use_dq) then
                call rhd_dq_init(r_parfile, r_modelfile, r_nmodel, r_imodel, 6)
        endif
        if (r_nmodel/=0) then
                call uio_struct_rd(r_modelfile, channel=r_nc, &
                                   startlabel='dataset', &
                                   endlabel='enddataset', &
                                   substartlabel='box', &
                                   subendlabel='endbox')
        else
                !print*, 'Inside else !!!'
                call uio_struct_rd(r_modelfile, channel=r_nc, &
                                   endlabel='enddataset')
        endif
        if (use_dq .and. r_imodel >=0) then
                call rhd_dq_end_init()
        endif
endif
end subroutine read_next
!
subroutine get_dq(uio_data_reg, dq)
use uio_reader_module
use rhd_dq_module
implicit none
external uio_data_reg
character*(*), intent(in)               :: dq
real, pointer, dimension(:,:,:)         :: tmp3D
!
record => uio_data_reg
select case (dq)
  case ('T')
          call rhd_dq_vars(T=tmp3D)
          label='T'
          call record('link')
  case ('P')
          call rhd_dq_vars(P=tmp3D)
          label='P'
          call record('link')
  case ('s')
          call rhd_dq_vars(s=tmp3D)
          label='s'
          call record('link')
  case ('j1')
          call rhd_dq_vars(j1=tmp3D)
          label='j1'
          call record('link')
  case ('j2')
          call rhd_dq_vars(j2=tmp3D)
          label='j2'
          call record('link')
  case ('j3')
          call rhd_dq_vars(j3=tmp3D)
          label='j3'
          call record('link')
  case ('jabs')
          call rhd_dq_vars(jabs=tmp3D)
          label='jabs'
          call record('link')
  case ('kappa')
          call rhd_dq_vars(kappa=tmp3D)
          label='kappa'
          call record('link')
  case ('Bc1')
          call rhd_dq_vars(Bc1=tmp3D)
          label='Bc1'
          call record('link')
  case ('Bc2')
          call rhd_dq_vars(Bc2=tmp3D)
          label='Bc2'
          call record('link')
  case ('Bc3')
          call rhd_dq_vars(Bc3=tmp3D)
          label='Bc3'
          call record('link')
  case ('divB')
          call rhd_dq_vars(divB=tmp3D)
          label='divB'
          call record('link')
  case ('vabs')
          call rhd_dq_vars(vabs=tmp3D)
          label='vabs'
          call record('link')
  case ('vh')
          call rhd_dq_vars(vh=tmp3D)
          label='vh'
          call record('link')
  case ('ekin')
          call rhd_dq_vars(ekin=tmp3D)
          label='ekin'
          call record('link')
  case ('plin')
          call rhd_dq_vars(plin=tmp3D)
          label='plin'
          call record('link')
  case ('vmflux')
          call rhd_dq_vars(vmflux=tmp3D)
          label='vmflux'
          call record('link')
  case ('g1')
          call rhd_dq_vars(g1=tmp3D)
          label='g1'
          call record('link')
  case ('g3')
          call rhd_dq_vars(g3=tmp3D)
          label='g3'
          call record('link')
  case ('cs')
          call rhd_dq_vars(cs=tmp3D)
          label='cs'
          call record('link')
  case ('mach')
          call rhd_dq_vars(mach=tmp3D)
          label='mach'
          call record('link')
  case ('mmu')
          call rhd_dq_vars(mmu=tmp3D)
          label='mmu'
          call record('link')
  case ('tau')
          call rhd_dq_vars(tau=tmp3D)
          label='tau'
          call record('link')
  case ('Babs')
          call rhd_dq_vars(Babs=tmp3D)
          label='Babs'
          call record('link')
  case ('Bh')
          call rhd_dq_vars(Bh=tmp3D)
          label='Bh'
          call record('link')
  case ('B2')
          call rhd_dq_vars(B2=tmp3D)
          label='B2'
          call record('link')
  case ('emag')
          call rhd_dq_vars(emag=tmp3D)
          label='emag'
          call record('link')
  case ('cA')
          call rhd_dq_vars(cA=tmp3D)
          label='cA'
          call record('link')
  case ('beta')
          call rhd_dq_vars(beta=tmp3D)
          label='beta'
          call record('link')
  case ('csca')
          call rhd_dq_vars(csca=tmp3D)
          label='csca'
          call record('link')
end select
end subroutine get_dq
!
subroutine write_model(modelfile, xb1, xb2, xb3, xc1, xc2, xc3, &
                       v1, v2, v3, rho, ei, Bb_flag, Bb1, Bb2, Bb3, &
                       B1_unit, B2_unit, B3_unit, dtime, itime, time, &
                       time_db, description, history, version, &
                       time_out_mean_last, time_out_full_last, &
                       m1, n1, m2, n2, m3, n3)
!use const_module
use rhd_gl_module
use rhd_box_module
use rhd_sub_module
!use rhd_io_module
!use tabinter_coeff_module
use uio_base_module
use rhd_prop_module
implicit none
!
! --- I/O parameters ---
character*(*),              intent(in)  :: modelfile
integer,                    intent(in)  :: itime, m1, n1, m2, n2, m3, n3
real, target, dimension(m1:n1+1), &
                         intent(inout)  :: xb1
real, target, dimension(m2:n2+1), &
                         intent(inout)  :: xb2
real, target, dimension(m3:n3+1), &
                         intent(inout)  :: xb3
real, target, dimension(m1:n1), &
                         intent(inout)  :: xc1
real, target, dimension(m2:n2), &
                         intent(inout)  :: xc2
real, target, dimension(m3:n3), &
                         intent(inout)  :: xc3
real, target, &
       dimension(m1:n1,m2:n2,m3:n3), &
                         intent(inout)  :: v1, v2, v3, rho, ei
real, target, &
     dimension(m1:n1+1,m2:n2,m3:n3), &
                         intent(inout)  :: Bb1
real, target, &
     dimension(m1:n1,m2:n2+1,m3:n3), &
                         intent(inout)  :: Bb2
real, target, &
     dimension(m1:n1,m2:n2,m3:n3+1), &
                         intent(inout)  :: Bb3
real,                       intent(in)  :: dtime, time, &
                                           time_out_mean_last,time_out_full_last
real(kind=kind(0.0D+00)),   intent(in)  :: time_db
character, dimension(:,:), &
                            intent(in)  :: history
character(len=80),          intent(in)  :: version, description, &
                                           B1_unit, B2_unit, B3_unit
logical,                    intent(in)  :: Bb_flag
!
! --- Local variables ---
type(box_type)                          :: box
character(len=80)                       :: date_end
character(len=80), dimension(nhismax)   :: history_in, history_end
character(len=80), dimension(1)         :: description_end
integer                                 :: ncbox, nhis_in, nhis_end, nchars, &
                                           i, j
!
nhis_in=size(history,1)
nchars=size(history,2)
do i=1,nhismax
  history_in(i)=repeat(' ',80)
  history_end(i)=repeat(' ',80)
end do
do i=1,nhis_in
  do j=1,nchars
    history_in(i)(j:j)=history(i,j)
  end do
end do
!print*, modelfile, xb1, xb2, xb3, xc1, xc2, xc3, &
!                       v1, v2, v3, rho, ei, Bb_flag, Bb1, Bb2, Bb3, &
!                       itime, time, time_db, history_in, &
!                       time_out_mean_last, time_out_full_last, &
!                       m1, n1, m2, n2, m3, n3, shape(history)
!
call rhd_box_Init(box)
!
box%m1=m1
box%n1=n1
box%m2=m2
box%n2=n2
box%m3=m3
box%n3=n3
!
box%itime=itime
box%time=time
box%time_db=time_db
!
box%xb1=>xb1(m1:n1+1)
box%xb2=>xb2(m2:n2+1)
box%xb3=>xb3(m3:n3+1)
box%xc1=>xc1(m1:n1)
box%xc2=>xc2(m2:n2)
box%xc3=>xc3(m3:n3)
!
box%rho=>rho(m1:n1,m2:n2,m3:n3)
box%ei=>ei(m1:n1,m2:n2,m3:n3)
box%v1=>v1(m1:n1,m2:n2,m3:n3)
box%v2=>v2(m1:n1,m2:n2,m3:n3)
box%v3=>v3(m1:n1,m2:n2,m3:n3)
!
prop%Bb1_IsInFile_flag=.false.
prop%Bb2_IsInFile_flag=.false.
prop%Bb3_IsInFile_flag=.false.
if(Bb_flag) then
  box%Bb1=>Bb1(m1:n1+1,m2:n2,m3:n3)
  box%Bb2=>Bb2(m1:n1,m2:n2+1,m3:n3)
  box%Bb3=>Bb3(m1:n1,m2:n2,m3:n3+1)
  prop%Bb1_IsInFile_flag=.true.
  prop%Bb2_IsInFile_flag=.true.
  prop%Bb3_IsInFile_flag=.true.
  prop%B1_unit=B1_unit
  prop%B2_unit=B2_unit
  prop%B3_unit=B3_unit
endif
!
! --- History ---
call uio_dattim(date_end)
history_end(1:nhis_in)=history_in(1:nhis_in)
if (nhis_in > (nhismax-1)) then
  history_end(nhismax/2-1)='...'
  history_end(nhismax/2:nhismax-3)=history_in(nhis_in-(nhismax-3)+nhismax/2:nhis_in)
  nhis_end=nhismax-1
else
  nhis_end=nhis_in
endif
nhis_end=nhis_end+1
history_end(nhis_end)='Box saved with pybold: ' // trim(date_end)
!
if (trim(description) == '') then
  description_end=(/'Pybold-generated RHD-simulation model snapshot'/)
else
  description_end(1) = description
endif
!print*, 'open,model,box,endmo,close', ncbox, dtime, &
!                  time_out_full_last, &
!                  time_out_mean_last, modelfile, &
!                  'unformatted', 'ieee_4', 1, &
!                  description, nhis_end, history_end, &
!                  trim(version)
call rhd_box_Write('open,model,box,endmo,close', ncbox, box=box, dtime=dtime, &
                  time_out_full_last=time_out_full_last, &
                  time_out_mean_last=time_out_mean_last, file=modelfile, &
                  form='unformatted', conv='ieee_4', ndes=1, &
                  description=description_end, nhis=nhis_end, &
                  history=history_end, version=trim(version))
!
! --- Open file for writing in uio-form ---
!call rhd_box_WrData('open', gl%nc_p, ncbox, file=modelfile, &
!                   form='unformatted', conv='ieee_4', ndes=1, &
!                   description=description, nhis=nhis_end, &
!                   history=history_end, version=trim(version))
! --- Write model header ---
!call rhd_box_WrData('model', gl%nc_p, ncbox, time_unit=time_unit, &
!                   itime=itime, time=time, time_db=time_db, dtime=dtime, &
!                   time_out_full_last=time_out_full_last, &
!                   time_out_mean_last=time_out_mean_last)
! --- Write modelbox ---
!call rhd_box_WrData('box', gl%nc_p, ncbox, &
!                   time_unit=time_unit, &
!                   x1_unit=x1_unit, x2_unit=x2_unit, x3_unit=x3_unit, &
!                   rho_unit=rho_unit, v_unit=v_unit, e_unit=e_unit, &
!                   box_id=box_id, &
!                   m1=m1, m2=m2, m3=m3, n1=n1, n2=n2, n3=n3, &
!                   itime=itime, time=time, time_db=time_db, dtime=dtime, &
!                   xc1=xc1, xc2=xc2, xc3=xc3, &
!                   xb1=xb1, xb2=xb2, xb3=xb3, &
!                   rho=rho, ei=ei, v1=v1, v2=v2, v3=v3, &
!                   Bb1=Bb1, Bb2=Bb2, Bb3=Bb3)
! --- Close file ---
!call rhd_box_WrData('endmo', gl%nc_p, ncbox)
! --- Close file ---
!call rhd_box_WrData('close', gl%nc_p, ncbox)
end subroutine write_model


text:
!=======================================================================!
! Copyright (c) Intel Corporation - All rights reserved.                !
! This file is part of the LIBXSMM library.                             !
!                                                                       !
! For information on the license, see the LICENSE file.                 !
! Further information: https://github.com/hfp/libxsmm/                  !
! SPDX-License-Identifier: BSD-3-Clause                                 !
!=======================================================================!
! Hans Pabst (Intel Corp.), Alexander Heinecke (Intel Corp.)
!=======================================================================!

      PROGRAM smm
        USE :: LIBXSMM, libxsmm_mmcall => libxsmm_dmmcall_abc
        !$ USE omp_lib
        IMPLICIT NONE

        INTEGER, PARAMETER :: T = KIND(0D0)

        REAL(T), ALLOCATABLE :: a(:,:,:), b(:,:,:), c(:,:,:)
        REAL(T), ALLOCATABLE :: sumxsmm(:,:), sumblas(:,:)
        REAL(T), ALLOCATABLE, SAVE :: sumtmp(:,:)
        !DIR$ ATTRIBUTES ALIGN:64 :: a, b, c, sumxsmm, sumtmp
        !$OMP THREADPRIVATE(sumtmp)
        TYPE(LIBXSMM_DMMFUNCTION) :: xmm
        INTEGER(8) :: i, r, s, size0, size1, size2, repetitions, start
        TYPE(LIBXSMM_MATDIFF_INFO) :: diff, max_diff
        INTEGER(LIBXSMM_BLASINT_KIND) :: m, n, k
        DOUBLE PRECISION :: duration, scale
        CHARACTER(32) :: argv
        INTEGER :: argc

        argc = COMMAND_ARGUMENT_COUNT()
        IF (1 <= argc) THEN
          CALL GET_COMMAND_ARGUMENT(1, argv)
          READ(argv, "(I32)") m
        ELSE
          m = 23
        END IF
        IF (3 <= argc) THEN
          CALL GET_COMMAND_ARGUMENT(3, argv)
          READ(argv, "(I32)") k
        ELSE
          k = m
        END IF
        IF (2 <= argc) THEN
          CALL GET_COMMAND_ARGUMENT(2, argv)
          READ(argv, "(I32)") n
        ELSE
          n = k
        END IF
        IF (4 <= argc) THEN
          CALL GET_COMMAND_ARGUMENT(4, argv)
          READ(argv, "(I32)") size1
        ELSE
          size1 = 0
        END IF
        IF (5 <= argc) THEN
          CALL GET_COMMAND_ARGUMENT(5, argv)
          READ(argv, "(I32)") size2
        ELSE
          size2 = 0 ! 1 repetition by default
        END IF

        ! Initialize LIBXSMM
        CALL libxsmm_init()

        ! JIT-compile stream(a,b)-kernel
        CALL libxsmm_dispatch(xmm, m, n, k,                             &
     &    alpha=REAL(1,T), beta=REAL(1,T))

        ! workload is about 2 GByte in memory by default
        size0 = (m * k + k * n + m * n) * T ! size of stream(a,b,c)-element in Byte
        size1 = MERGE(2048_8, MERGE(size1, ISHFT(ABS(size0 * size1)     &
     &            + ISHFT(1, 20) - 1, -20), 0.LE.size1), 0.EQ.size1)
        size2 = ISHFT(MERGE(MAX(size2, size1), ISHFT(ABS(size2) * size0 &
     &            + ISHFT(1, 20) - 1, -20), 0.LE.size2), 20) / size0
        s = ISHFT(size1, 20) / size0
        repetitions = size2 / s
        scale = 1D0 / s
        duration = 0

        CALL libxsmm_matdiff_clear(max_diff)
        ALLOCATE(a(m,k,s))
        ALLOCATE(b(k,n,s))

        ! Initialize a, b
        !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(s, a, b, scale)
        DO i = 1, s
          CALL init(42, a(:,:,i), scale, i - 1)
          CALL init(24, b(:,:,i), scale, i - 1)
        END DO
        !$OMP END PARALLEL DO

        WRITE(*, "(3(A,I0),A,I0,A,I0,A,I0)")                            &
     &    "m=", m, " n=", n, " k=", k, " elements=", UBOUND(a, 3),      &
     &    " size=", size1, " MB repetitions=", repetitions

        ! compute reference solution and warmup BLAS library
        ALLOCATE(sumblas(m,n))
        sumblas(:,:) = 0
        !$OMP PARALLEL REDUCTION(+:sumblas) PRIVATE(i, r)               &
        !$OMP   DEFAULT(NONE) SHARED(m, n, k, a, b, repetitions)
        ALLOCATE(sumtmp(m,n))
        sumtmp(:,:) = 0
        DO r = 1, repetitions
          !$OMP DO
          DO i = LBOUND(a, 3), UBOUND(a, 3)
            ! PGI: cannot deduce generic procedure (libxsmm_blas_gemm)
            CALL libxsmm_blas_dgemm(m=m, n=n, k=k,                      &
     &              a=a(:,:,i), b=b(:,:,i), c=sumtmp)
          END DO
        END DO
        sumblas(:,:) = sumblas(:,:) + sumtmp(:UBOUND(sumblas,1),:)
        ! Deallocate thread-local arrays
        DEALLOCATE(sumtmp)
        !$OMP END PARALLEL

        WRITE(*, "(A)") "Streamed (A,B)... (BLAS)"
        ALLOCATE(sumxsmm(m,n))
        sumxsmm(:,:) = 0
        !$OMP PARALLEL REDUCTION(+:sumxsmm) PRIVATE(i, r, start)        &
        !$OMP   DEFAULT(NONE)                                           &
        !$OMP   SHARED(m, n, k, a, b, duration, repetitions)
        ALLOCATE(sumtmp(m,n))
        sumtmp(:,:) = 0
        !$OMP MASTER
        start = libxsmm_timer_tick()
        !$OMP END MASTER
        !$OMP BARRIER
        DO r = 1, repetitions
          !$OMP DO
          DO i = LBOUND(a, 3), UBOUND(a, 3)
            ! PGI: cannot deduce generic procedure (libxsmm_blas_gemm)
            CALL libxsmm_blas_dgemm(m=m, n=n, k=k,                      &
     &              a=a(:,:,i), b=b(:,:,i), c=sumtmp)
          END DO
        END DO
        !$OMP BARRIER
        !$OMP MASTER
        duration = libxsmm_timer_duration(start, libxsmm_timer_tick())
        !$OMP END MASTER
        sumxsmm(:,:) = sumxsmm(:,:) + sumtmp(:UBOUND(sumxsmm,1),:)
        ! Deallocate thread-local arrays
        DEALLOCATE(sumtmp)
        !$OMP END PARALLEL
        CALL performance(duration, m, n, k, size2, m * k + k * n)

        WRITE(*, "(A)") "Streamed (A,B)... (auto-dispatched)"
        sumxsmm(:,:) = 0
        !$OMP PARALLEL REDUCTION(+:sumxsmm) PRIVATE(i, r, start)        &
        !$OMP   DEFAULT(NONE)                                           &
        !$OMP   SHARED(m, n, k, a, b, duration, repetitions)
        ALLOCATE(sumtmp(m,n))
        sumtmp(:,:) = 0
        !$OMP MASTER
        start = libxsmm_timer_tick()
        !$OMP END MASTER
        !$OMP BARRIER
        DO r = 1, repetitions
          !$OMP DO
          DO i = LBOUND(a, 3), UBOUND(a, 3)
            ! PGI: cannot deduce generic procedure (libxsmm_gemm)
            CALL libxsmm_dgemm(m=m, n=n, k=k,                           &
     &              a=a(:,:,i), b=b(:,:,i), c=sumtmp)
          END DO
        END DO
        !$OMP BARRIER
        !$OMP MASTER
        duration = libxsmm_timer_duration(start, libxsmm_timer_tick())
        !$OMP END MASTER
        sumxsmm(:,:) = sumxsmm(:,:) + sumtmp(:UBOUND(sumxsmm,1),:)
        ! Deallocate thread-local arrays
        DEALLOCATE(sumtmp)
        !$OMP END PARALLEL
        CALL performance(duration, m, n, k, size2, m * k + k * n)
        CALL libxsmm_matdiff(diff, LIBXSMM_DATATYPE_F64, m, n,          &
     &    libxsmm_ptr(sumblas), libxsmm_ptr(sumxsmm))
        WRITE(*, "(1A,A,F10.1)") CHAR(9), "diff:      ", diff%l2_abs
        CALL libxsmm_matdiff_reduce(max_diff, diff)

        IF (libxsmm_available(xmm)) THEN
          sumxsmm(:,:) = 0
          WRITE(*, "(A)") "Streamed (A,B)... (specialized)"
          !$OMP PARALLEL REDUCTION(+:sumxsmm) PRIVATE(i, r, start)
            !DEFAULT(NONE) SHARED(m, n, a, b, duration, repetitions, xmm)
          ALLOCATE(sumtmp(m,n))
          sumtmp(:,:) = 0
          !$OMP MASTER
          start = libxsmm_timer_tick()
          !$OMP END MASTER
          !$OMP BARRIER
          DO r = 1, repetitions
            !$OMP DO
            DO i = LBOUND(a, 3), UBOUND(a, 3)
              CALL libxsmm_mmcall(xmm, a(:,:,i), b(:,:,i), sumtmp)
            END DO
          END DO
          !$OMP BARRIER
          !$OMP MASTER
          duration = libxsmm_timer_duration(start, libxsmm_timer_tick())
          !$OMP END MASTER
          sumxsmm(:,:) = sumxsmm(:,:) + sumtmp(:UBOUND(sumxsmm,1),:)
          ! Deallocate thread-local arrays
          DEALLOCATE(sumtmp)
          !$OMP END PARALLEL
          CALL performance(duration, m, n, k, size2, m * k + k * n)
          CALL libxsmm_matdiff(diff, LIBXSMM_DATATYPE_F64, m, n,        &
     &      libxsmm_ptr(sumblas), libxsmm_ptr(sumxsmm))
          WRITE(*, "(1A,A,F10.1)") CHAR(9), "diff:      ", diff%l2_abs
          CALL libxsmm_matdiff_reduce(max_diff, diff)
        END IF

        WRITE(*, "(A)") "Streamed (A,B,C/Beta=0)... (specialized)"
        ALLOCATE(c(m,n,s))
        start = libxsmm_timer_tick()
        DO r = 1, repetitions
          !$OMP PARALLEL DEFAULT(NONE) SHARED(a, b, c)
          CALL smm_stream_abc_beta0(a, b, c)
          !$OMP END PARALLEL
        END DO
        duration = libxsmm_timer_duration(start, libxsmm_timer_tick())
        CALL performance(duration, m, n, k, size2,                      &
     &    m * k + k * n + m * n * 2) ! 2: assume RFO

        WRITE(*, "(A)") "Streamed (A,B,C/Beta=1)... (specialized)"
        start = libxsmm_timer_tick()
        DO r = 1, repetitions
          !$OMP PARALLEL DEFAULT(NONE) SHARED(a, b, c)
          CALL smm_stream_abc_beta1(a, b, c)
          !$OMP END PARALLEL
        END DO
        duration = libxsmm_timer_duration(start, libxsmm_timer_tick())
        CALL performance(duration, m, n, k, size2,                      &
     &    m * k + k * n + m * n * 2) ! 2: assume RFO

        ! Deallocate global arrays
        DEALLOCATE(a, b, c)
        DEALLOCATE(sumxsmm)
        DEALLOCATE(sumblas)

        ! finalize LIBXSMM
        CALL libxsmm_finalize()

        IF (1.LT.(max_diff%l2_rel)) STOP 1

      CONTAINS
        PURE SUBROUTINE init(seed, matrix, scale, n)
          INTEGER, INTENT(IN) :: seed
          REAL(T), INTENT(OUT) :: matrix(:,:)
          REAL(8), INTENT(IN) :: scale
          INTEGER(8), INTENT(IN), OPTIONAL :: n
          INTEGER(8) :: minval, addval, maxval
          INTEGER :: ld, i, j
          REAL(8) :: val, norm
          ld = UBOUND(matrix, 1) - LBOUND(matrix, 1) + 1
          minval = MERGE(n, 0_8, PRESENT(n)) + seed
          addval = (UBOUND(matrix, 1) - LBOUND(matrix, 1)) * ld         &
     &           + (UBOUND(matrix, 2) - LBOUND(matrix, 2))
          maxval = MAX(ABS(minval), addval)
          norm = MERGE(scale / maxval, scale, 0.NE.maxval)
          DO j = LBOUND(matrix, 2), UBOUND(matrix, 2)
            DO i = LBOUND(matrix, 1),                                   &
     &             LBOUND(matrix, 1) + UBOUND(matrix, 1) - 1
              val = (i - LBOUND(matrix, 1)) * ld                        &
     &            + (j - LBOUND(matrix, 2)) + minval
              matrix(i,j) = norm * (val - 0.5D0 * addval)
            END DO
          END DO
        END SUBROUTINE

        SUBROUTINE disp(matrix, ld, format)
          REAL(T), INTENT(IN) :: matrix(:,:)
          INTEGER, INTENT(IN), OPTIONAL :: ld
          CHARACTER(*), INTENT(IN), OPTIONAL :: format
          CHARACTER(32) :: fmt
          INTEGER :: i0, i1, i, j
          IF (.NOT.PRESENT(format)) THEN
            fmt = "(16F20.0)"
          ELSE
            WRITE(fmt, "('(16',A,')')") format
          END IF
          i0 = LBOUND(matrix, 1)
          i1 = MIN(                                                     &
     &          MERGE(i0 + ld - 1, UBOUND(matrix, 1), PRESENT(ld)),     &
     &          UBOUND(matrix, 1))
          DO i = i0, i1
            DO j = LBOUND(matrix, 2), UBOUND(matrix, 2)
              WRITE(*, fmt, advance='NO') matrix(i,j)
            END DO
            WRITE(*, *)
          END DO
        END SUBROUTINE

        SUBROUTINE performance(duration, m, n, k, s, c)
          INTEGER(LIBXSMM_BLASINT_KIND), INTENT(IN) :: m, n, k, c
          INTEGER(8), INTENT(IN) :: s
          REAL(T), INTENT(IN) :: duration
          IF (0.LT.duration) THEN
            WRITE(*, "(1A,A,F10.1,A)") CHAR(9), "performance:",         &
     &        2D0 * s * m * n * k * 1D-9 / duration, " GFLOPS/s"
            WRITE(*, "(1A,A,F10.1,A)") CHAR(9), "bandwidth:  ",         &
     &        s * c * T / (duration * ISHFT(1_8, 30)), " GB/s"
          END IF
          WRITE(*, "(1A,A,F10.1,A)") CHAR(9), "duration:   ",           &
     &        1D3 * duration, " ms"
        END SUBROUTINE

        SUBROUTINE smm_stream_abc_beta0(a, b, c)
          REAL(T), INTENT(IN)  :: a(:,:,:), b(:,:,:)
          REAL(T), INTENT(OUT) :: c(:,:,:)
          TYPE(LIBXSMM_DMMFUNCTION) :: xmm
          INTEGER(8) :: i
          CALL libxsmm_dispatch(xmm, SIZE(c,1), SIZE(c,2), SIZE(b,1),   &
     &      flags=MERGE(                                                &
     &        LIBXSMM_GEMM_FLAG_ALIGN_C_NTS_HINT_BETA_0,                &
     &        LIBXSMM_GEMM_FLAG_BETA_0,                                 &
     &        libxsmm_aligned(libxsmm_ptr(c(1,1,1)),                    &
     &        SIZE(c,1) * SIZE(c,2) * T)))
          !$OMP DO
          DO i = LBOUND(c, 3, 8), UBOUND(c, 3, 8)
            CALL libxsmm_mmcall(xmm, a(:,:,i), b(:,:,i), c(:,:,i))
          END DO
        END SUBROUTINE

        SUBROUTINE smm_stream_abc_beta1(a, b, c)
          REAL(T), INTENT(IN)    :: a(:,:,:), b(:,:,:)
          REAL(T), INTENT(INOUT) :: c(:,:,:)
          TYPE(LIBXSMM_DMMFUNCTION) :: xmm
          INTEGER(8) :: i
          CALL libxsmm_dispatch(xmm, SIZE(c,1), SIZE(c,2), SIZE(b,1))
          !$OMP DO
          DO i = LBOUND(c, 3, 8), UBOUND(c, 3, 8)
            CALL libxsmm_mmcall(xmm, a(:,:,i), b(:,:,i), c(:,:,i))
          END DO
        END SUBROUTINE
      END PROGRAM

max_stars_repo_name: schtomerintel/libxsmm

max_issues_repo_name: schtomerintel/libxsmm

max_forks_repo_name: schtomerintel/libxsmm


text:
program test_ewald2

! This test compares lattice energy 'eew' calculated from the Ewald summation
! against the energy calculated using the Madelung constant for various lattice
! constants L. The agreement is essentially to machine precision.

! Similar to test_ewald, but here we the diagonal Na atom is moved by 3/8
! towards the Cl atom.

use types, only: dp
use constants, only: ang2bohr, kJmol2Ha
use ewald_sums, only: ewald, ewald2, fred2fcart
use utils, only: assert, init_random
implicit none

integer :: natom, ntypat
! Various NaCl lattice constants in A
real(dp), parameter :: Llist(5) = [5.6402_dp, 5.5_dp, 4.5_dp, 6.5_dp, 10._dp]
! The components of the correct force and stress tensor (multiplied by L)
real(dp), parameter :: fcorrect0 = 10.180280846982683_dp
real(dp), parameter :: stress1 = -6.6903565769534099_dp
real(dp), parameter :: stress2 = -5.7797035916557675_dp

real(dp) :: ucvol
integer, allocatable :: typat(:)
real(dp) :: gmet(3, 3), rmet(3, 3), rprim(3, 3), gprim(3, 3), stress(6)
real(dp) :: eew
real(dp), allocatable :: xred(:, :), zion(:), grewtn(:, :), fcart(:, :)
real(dp) :: L, alpha, E_ewald, E_madelung
integer :: i, j

! Madelung constant for NaCl, where the diagonal Na atom is moved by 3/8
! towards the Cl atom
alpha = 2.5088837163575413_dp

! Conventional cell:
natom = 8
ntypat = 2
allocate(xred(3, natom), zion(ntypat), grewtn(3, natom), typat(natom), &
    fcart(3, natom))
! Cl^-
xred(:, 1) = [0, 0, 0]
xred(:, 2) = [1, 1, 0] / 2._dp
xred(:, 3) = [1, 0, 1] / 2._dp
xred(:, 4) = [0, 1, 1] / 2._dp
! Na^+
xred(:, 5) = [1, 1, 1] / 2._dp
xred(:, 6) = [1, 0, 0] / 2._dp
xred(:, 7) = [0, 1, 0] / 2._dp
xred(:, 8) = [0, 0, 1] / 2._dp
xred(:, 5:8) = xred(:, 5:8) - spread([3, 3, 3] / 8._dp, 2, 4)
xred(:, 6:8) = xred(:, 6:8) - floor(xred(:, 6:8))
typat = [1, 1, 1, 1, 2, 2, 2, 2]
zion = [-1, +1]

do i = 1, size(Llist)
    L = Llist(i) * ang2bohr
    print *, L

    rmet = 0
    rmet(1, 1) = L**2
    rmet(2, 2) = L**2
    rmet(3, 3) = L**2

    ! gmet = inv(rmet)
    gmet = 0
    gmet(1, 1) = 1/L**2
    gmet(2, 2) = 1/L**2
    gmet(3, 3) = 1/L**2

    ! ucvol = sqrt(det(rmet))
    ucvol = L**3

    ! Reciprocal primitive vectors (without 2*pi) in cartesian coordinates.
    ! gmet = matmul(transpose(gprim), gprim)
    gprim = 0
    gprim(1, 1) = 1 / L
    gprim(2, 2) = 1 / L
    gprim(3, 3) = 1 / L

    ! Real space primitive vectors
    rprim = 0
    rprim(1, 1) = L
    rprim(2, 2) = L
    rprim(3, 3) = L

    call ewald(eew,gmet,grewtn,natom,ntypat,rmet,typat,ucvol,xred,zion)

    E_ewald = eew / (natom/ntypat)
    E_madelung = -2*alpha/L

    call ewald2(gmet,natom,ntypat,rmet,rprim,gprim,stress,typat,ucvol, &
                xred,zion)

    print *, "a =", L/ang2bohr*100, "pm"
    print *, "Madelung:", E_madelung / kJmol2Ha, "kJ/mol"
    print *, "Ewald:   ", E_ewald / kJmol2Ha, "kJ/mol"
    print *, "error:   ", abs(E_ewald - E_madelung), "a.u."
    call assert(abs(E_ewald - E_madelung) < 1e-14_dp)
    stress = -stress * ucvol
    call assert(all(abs(stress - [stress1, stress1, stress1, &
        stress2, stress2, stress2]/L) < 1e-10_dp))

    call fred2fcart(fcart, grewtn, gprim)
    do j = 1, 4
        call assert(all(abs(fcart(:, j) - (-fcorrect0/L**2)) < 1e-10_dp))
    end do
    do j = 5, 8
        call assert(all(abs(fcart(:, j) - (+fcorrect0/L**2)) < 1e-10_dp))
    end do
end do
deallocate(xred, zion, grewtn, typat, fcart)

print *, "--------"

! Primitive cell (FCC lattice)
natom = 2
ntypat = 2
allocate(xred(3, natom), zion(ntypat), grewtn(3, natom), typat(natom), &
    fcart(3, natom))
! Cl^-
xred(:, 1) = [0, 0, 0]
! Na^+
xred(:, 2) = [1, 1, 1] / 2._dp
xred(:, 2) = xred(:, 2) - [3, 3, 3] / 8._dp
typat = [1, 2]
zion = [-1._dp, 1._dp]

do i = 1, size(Llist)
    L = Llist(i) * ang2bohr

    rmet(1, :) = [2, 1, 1]
    rmet(2, :) = [1, 2, 1]
    rmet(3, :) = [1, 1, 2]
    rmet = rmet * L**2 / 4

    ! gmet = inv(rmet)
    gmet(1, :) = [ 3, -1, -1]
    gmet(2, :) = [-1,  3, -1]
    gmet(3, :) = [-1, -1,  3]
    gmet = gmet / L**2

    ! ucvol = sqrt(det(rmet))
    ucvol = L**3 / 4

    ! Reciprocal primitive vectors (without 2*pi) in cartesian coordinates.
    ! gmet = matmul(transpose(gprim), gprim)
    gprim(:, 1) = [ 1,  1, -1] / L
    gprim(:, 2) = [-1,  1,  1] / L
    gprim(:, 3) = [ 1, -1,  1] / L

    rprim(:, 1) = [1, 0, 1] * L / 2
    rprim(:, 2) = [1, 1, 0] * L / 2
    rprim(:, 3) = [0, 1, 1] * L / 2

    call ewald(eew,gmet,grewtn,natom,ntypat,rmet,typat,ucvol,xred,zion)

    E_ewald = eew / (natom/ntypat)
    E_madelung = -2*alpha/L

    call ewald2(gmet,natom,ntypat,rmet,rprim,gprim,stress,typat,ucvol, &
                xred,zion)

    print *, "a =", L/ang2bohr*100, "pm"
    print *, "Madelung:", E_madelung / kJmol2Ha, "kJ/mol"
    print *, "Ewald:   ", E_ewald / kJmol2Ha, "kJ/mol"
    print *, "error:   ", abs(E_ewald - E_madelung), "a.u."
    call assert(abs(E_ewald - E_madelung) < 1e-14_dp)
    call fred2fcart(fcart, grewtn, gprim)
    call assert(all(abs(fcart(:, 1) - (-fcorrect0/L**2)) < 1e-10_dp))
    call assert(all(abs(fcart(:, 2) - (+fcorrect0/L**2)) < 1e-10_dp))
    stress = -stress * L**3
    call assert(all(abs(stress - [stress1, stress1, stress1, &
        stress2, stress2, stress2]/L) < 1e-10_dp))
end do
deallocate(xred, zion, grewtn, typat, fcart)
end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
      SUBROUTINE PHO_POMSCA(Ii,Mspom,Mhpom,Msreg,Ival1,Ival2,Mspar1,
     &                      Mspar2,Mhpar1,Mhpar2,Irej)
C**********************************************************************
C
C     parton orientated formulation of soft and hard inelastic events
C
C
C     input:    II        particle combiantion (1..4)
C               MSPOM     number of soft pomerons
C               MHPOM     number of semihard pomerons
C               MSREG     number of soft reggeons
C
C     output:   IVAL1,2   0 no valence quark engaged
C                         otherwise:  position of valence quark engaged
C                         neg.number: gluon connected to valence quark
C                                     by color flow
C               MSPAR1,2  number of realized soft partons
C               MHPAR1,2  number of realized hard partons
C               IREJ      1 failure
C                         0 success
C
C**********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION ahs , alnhs , DT_RNDM , emin , p1 , p2 , pd1 ,
     &                 ptmax , ptmx , q2h , sah , sas , TINY , valpro ,
     &                 wg1 , wgx , xhmax1 , xhmax2 , xisr1 , xisr2
      DOUBLE PRECISION xmax1 , xmax2 , xmaxh1 , xmaxh2 , xmaxp1 ,
     &                 xmaxp2 , xmaxx1 , xmaxx2 , xscut , xss1 , xss2 ,
     &                 xsss1 , xsss2 , xtmp1 , xtmp2 , xx , z1difs ,
     &                 z2difs
      INTEGER i , i1 , i2 , i3 , i4 , ifl1 , ifl2 , ihard , Ii , inmax ,
     &        iptm , iqua1 , iqua2 , Irej , itry , Ival1 , Ival2 ,
     &        ivg1 , ivg2 , ivglu1
      INTEGER ivglu2 , ivq1 , ivq2 , j , k , l , line , mhard , mhcha ,
     &        Mhpar1 , Mhpar2 , Mhpom , mscha , msdiff , msg1 , msg2 ,
     &        msm1 , msm2 , Mspar1 , Mspar2
      INTEGER Mspom , Msreg , n , nn , ntry
      SAVE

      PARAMETER (TINY=1.D-30)

C  input/output channels
      INCLUDE 'inc/poinou'
C  event debugging information
      INCLUDE 'inc/podebg'
C  model switches and parameters
      INCLUDE 'inc/pomdls'
C  current beam selection
      INCLUDE 'inc/pobeam'
C  general process information
      INCLUDE 'inc/poprcs'
C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
      INCLUDE 'inc/pohdfl'
C  event weights and generated cross section
      INCLUDE 'inc/powght'
C  hard cross sections and MC selection weights
      INCLUDE 'inc/pohrcs'
C  hard scattering parameters used for most recent hard interaction
      INCLUDE 'inc/pohapa'
C  data of c.m. system of Pomeron / Reggeon exchange
      INCLUDE 'inc/popcms'
C  obsolete cut-off information
      INCLUDE 'inc/pocut1'
C  some hadron information, will be deleted in future versions
      INCLUDE 'inc/pohdrn'
C  data on most recent hard scattering
      INCLUDE 'inc/pockin'
C  light-cone x fractions and c.m. momenta of soft cut string ends
      INCLUDE 'inc/posoft'
C  hard scattering data
      INCLUDE 'inc/pohslt'
C  table of particle indices for recursive PHOJET calls
      INCLUDE 'inc/porecu'
C  internal rejection counters
      INCLUDE 'inc/poloop'

      DIMENSION p1(4) , p2(4) , pd1(-6:6)

      IF ( LPRi.GT.4 .AND. IDEb(24).GT.20 ) WRITE (LO,'(1X,A,3I5)')
     &      'PHO_POMSCA: MSPOM,MHPOM,MSREG' , Mspom , Mhpom , Msreg

      itry = 0
      ntry = 10
      Irej = 0
      inmax = 10
      mhard = Mhpom

C  phase space limitation (single hard valence-valence quark scattering)
      IF ( Mhpom.GT.0 ) THEN
         emin = 2.D0*PTWant + 0.2D0
         IF ( ECMp.LT.emin ) THEN
            IF ( LPRi.GT.4 .AND. IDEb(24).GE.1 )
     &            WRITE (LO,'(1X,A,1P,3E10.3)') 'PHO_POMSCA: ' ,
     &           'kin. rejection (1) (Ecm,Ptcut,Emin)' , ECMp , PTWant ,
     &           emin
            Irej = 50
            IFAil(6) = IFAil(6) + 1
            RETURN
         END IF
      END IF

      sas = PARmdl(160+Ii)/ECMp
      sah = 2.D0*PTWant/ECMp
      AS = sas**2
      AH = sah**2

C  save energy for leading particle effect
      xmaxp1 = 1.D0
      IF ( (IPAmdl(13).GT.0) .AND. (IHFls(1).NE.0) ) xmaxp1 = 1.D0 -
     &     PARmdl(165)*XPSub
      xmaxp2 = 1.D0
      IF ( (IPAmdl(13).GT.0) .AND. (IHFls(2).NE.0) ) xmaxp2 = 1.D0 -
     &     PARmdl(165)*XTSub

C
C  main loop to select hard and soft parton kinematics
C -----------------------------------------------------
      IFAil(31) = IFAil(31) + mhard
 100  Irej = 0
      ihard = 0
      LSC1hd = 0
      itry = itry + 1
      IF ( itry.GT.1 ) IFAil(5) = IFAil(5) + 1
      IF ( itry.GE.ntry ) THEN
         Irej = 1
         GOTO 500
      END IF
      line = 0
      LSCahd = 0
      IF ( (IPAmdl(13).GT.0) .AND. (IPOix3.EQ.0) .AND. (IPRoce.EQ.1) )
     &     THEN
         xss1 = MAX(0.D0,1.D0-XPSub)
         xss2 = MAX(0.D0,1.D0-XTSub)
      ELSE
         xss1 = 0.D0
         xss2 = 0.D0
      END IF

C  partons needed to construct soft/hard interactions
 200  Mspar1 = 2*Mspom + Msreg + Mhpom
      Mspar2 = Mspar1
      Mhpar1 = Mhpom
      Mhpar2 = Mhpom

C  number of strings
      mscha = 2*Mspom + Msreg
      mhcha = 2*Mhpom

      KSOft = mscha
      KHArd = mhcha

C  check actual phase space limit
      xx = sas*DBLE(mscha) + sah*DBLE(mhcha)/2.D0
      IF ( xx.GE.1.D0 ) THEN
         IF ( LPRi.GT.4 .AND. IDEb(24).GE.3 )
     &         WRITE (LO,'(1X,2A,/1X,4I3,1P4E12.4)')
     &         'PHO_POMSCA: internal kin. rejection ' ,
     &        '(MSpom,MHpom,MSchain,MHchain,Ecm,AS,AH,XX):' , Mspom ,
     &        Mhpom , mscha , mhcha , ECMp , AS , AH , xx
         IF ( Mspom+Msreg+Mhpom.GT.1 ) THEN
            IF ( Msreg.GT.0 ) THEN
               Msreg = Msreg - 1
            ELSE IF ( Mspom.GT.0 ) THEN
               Mspom = Mspom - 1
            ELSE IF ( Mhpom.GT.1 ) THEN
               Mhpom = Mhpom - 1
            END IF
            GOTO 200
         END IF
         IF ( LPRi.GT.4 .AND. IDEb(24).GE.1 )
     &         WRITE (LO,'(1X,A,1P2E10.3)')
     &         'PHO_POMSCA: kin. rejection (2) (Ecm,Ptcut)' , ECMp ,
     &        PTWant
         Irej = 50
         IFAil(6) = IFAil(6) + 1
         RETURN
      END IF

      xmaxx1 = MAX(TINY,1.D0-MIN(Mspar1,1)*AS-MIN(Mhpar1,1)*AH)
      xmaxx2 = MAX(TINY,1.D0-MIN(Mspar2,1)*AS-MIN(Mhpar2,1)*AH)

C  very low energy phase space restriction
      IF ( mhard.GT.0 ) THEN
         IF ( (xmaxx1*xmaxx2.LE.AH) ) THEN
            IF ( LPRi.GT.4 .AND. IDEb(24).GE.1 )
     &            WRITE (LO,'(1X,A,1P2E10.3)')
     &            'PHO_POMSCA: kin. rejection (3) (Ecm,Ptcut)' , ECMp ,
     &           PTWant
            Irej = 50
            IFAil(6) = IFAil(6) + 1
            RETURN
         END IF
      END IF

      AS = MAX(AS,PSOmin/PCMp)
      ALNs = LOG(AS)
      ALNh = LOG(AH)
      Z1Max = LOG(xmaxx1)
      Z2Max = LOG(xmaxx2)
      Z1Dif = Z1Max + Z2Max - ALNh
      Z2Dif = Z1Dif
      ptmax = 0.D0
C
C  select hard parton momenta
C ------------------- begin of inner loop -------------------
      IF ( IPOix3.EQ.0 ) IPOwgc(4+Ii) = 0

      IF ( mhard.GT.MSCAHD ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,2I3)') 'PHO_POMSCA: ' ,
     &        'no space left in /POHSLT/ (MHARD,MSCAHD):' , mhard ,
     &        MSCAHD
         Irej = 1
         RETURN
      END IF

      DO nn = 1 , mhard
C
C  generate one resolved hard scattering
C
C  high-pt option
         IF ( (nn.EQ.1) .AND. (Ii.EQ.1) .AND. (HSWcut(4+Ii).GT.PTWant) )
     &        THEN
            CALL PHO_HARINT(-1,ECMp,PVIrtp(1),PVIrtp(2),-1,MAX_PRO_2,1,
     &                      4,Mspom+Mhpom)
            xscut = HSIg(9,IDXmpar)
            ahs = AH
            alnhs = ALNh
            z1difs = Z1Dif
            z2difs = Z2Dif
            AH = (2.D0*PTWant/ECMp)**2
            ALNh = LOG(AH)
            Z1Dif = Z1Max + Z2Max - ALNh
            Z2Dif = Z1Dif
            IF ( (Z1Dif.LE.0.01D0) .OR. (Z2Dif.LE.0.01D0) ) THEN
               IF ( LPRi.GT.4 .AND. IDEb(24).GE.1 )
     &               WRITE (LO,'(1X,2A,/1X,1P4E12.3)')
     &               'PHO_POMSCA: kin.rejection, high-pt option ' ,
     &              '(Z1/2max,ALNH,Z1dif):' , Z1Max , Z2Max , ALNh ,
     &              Z1Dif
               Irej = 5
               RETURN
            END IF
            CALL PHO_HARSCA(2,Ii)
            CALL PHO_HARINT(1,ECMp,PVIrtp(1),PVIrtp(2),-1,MAX_PRO_2,1,4,
     &                      Mspom+Mhpom)
            AH = ahs
            ALNh = alnhs
            Z1Dif = z1difs
            Z2Dif = z2difs
            IPOwgc(4+Ii) = IPOwgc(4+Ii) + 1
            HSWght(4+Ii) = xscut/HSIg(9,IDXmpar)*DBLE(mhard)
C  minimum bias option
         ELSE
            CALL PHO_HARSCA(2,Ii)
         END IF

C  fill /POHSLT/
         LSIdx(nn) = nn
         LSCahd = nn
         XHD(nn,1) = X1
         XHD(nn,2) = X2
         X0Hd(nn,1) = X1
         X0Hd(nn,2) = X2
         VHD(nn) = V
         ETAhd(nn,1) = ETAc
         ETAhd(nn,2) = ETAd
         PTHd(nn) = PT
         NPRohd(nn) = MSPr
         Q2Sca(nn,1) = QQPd
         Q2Sca(nn,2) = QQPd
         PDFva(nn,1) = PDF1(IA)
         PDFva(nn,2) = PDF2(IB)
         NINhd(nn,1) = IA
         NINhd(nn,2) = IB
         N0Inhd(nn,1) = IA
         N0Inhd(nn,2) = IB
         NIVal(nn,1) = IV1
         NIVal(nn,2) = IV2
         N0Ival(nn,1) = IV1
         N0Ival(nn,2) = IV2
         NOUthd(nn,1) = IC
         NOUthd(nn,2) = ID
         NBRahd(nn,1) = IDPdg1
         NBRahd(nn,2) = IDPdg2
         i3 = 8*(nn-1)
         i4 = 8*(nn-1) + 4
         DO i = 1 , 4
            PPH(i3+i,1) = PHI1(i)
            PPH(i3+i,2) = PHI2(i)
            PPH(i4+i,1) = PHO1(i)
            PPH(i4+i,2) = PHO2(i)
         END DO

      END DO

C  sort according to pt-hat
      DO nn = 1 , mhard
         ptmx = PTHd(LSIdx(nn))
         iptm = nn
         DO i = nn + 1 , mhard
            IF ( PTHd(LSIdx(i)).GT.ptmx ) THEN
               iptm = i
               ptmx = PTHd(LSIdx(i))
            END IF
         END DO
         IF ( iptm.NE.nn ) CALL PHO_SWAPI(LSIdx(nn),LSIdx(iptm))
      END DO
      iptm = LSIdx(1)

C  copy partons, generate ISR
      DO l = 1 , mhard
         nn = LSIdx(l)
         xsss1 = xss1 + XHD(nn,1)
         xsss2 = xss2 + XHD(nn,2)
C  debug output
         IF ( LPRi.GT.4 .AND. IDEb(24).GE.10 )
     &         WRITE (LO,'(1X,A,3I4,1P,3E11.3)')
     &         'PHO_POMSCA: NR,LSIDX,MSPR,X1,X2,PT' , l , nn ,
     &        NPRohd(nn) , XHD(nn,1) , XHD(nn,2) , PTHd(nn)
C  check phase space
         IF ( (xsss1.GT.xmaxx1) .OR. (xsss2.GT.xmaxx2) .OR.
     &        ((1.D0-xsss1)*(1.D0-xsss2).LT.AS) ) THEN
            IF ( ihard.EQ.0 ) THEN
               IF ( ISWmdl(2).NE.1 ) GOTO 100
               Mhpom = 0
               Mspom = 1
               Msreg = 0
            END IF
            GOTO 300
         END IF

C  reweight according to photon virtuality
         IF ( IPAmdl(115).GE.1 ) THEN
            QQPd = Q2Sca(nn,1)
            wgx = 1.D0
            IF ( IDPdg1.EQ.22 ) THEN
               IF ( IPAmdl(115).EQ.1 ) THEN
                  IF ( QQPd.LT.PVIrtp(1)+PARmdl(144) ) THEN
                     wg1 = 0.D0
                  ELSE
                     wg1 = LOG(QQPd/(PVIrtp(1)+PARmdl(144)))
     &                     /LOG(QQPd/PARmdl(144))
                  END IF
                  IF ( NINhd(nn,1).EQ.0 ) wg1 = wg1*wg1
               ELSE IF ( IPAmdl(115).EQ.2 ) THEN
                  CALL PHO_PDF(1,X0Hd(nn,1),QQPd,PVIrtp(1),pd1)
                  wg1 = pd1(NINhd(nn,1))/PDFva(nn,1)
               END IF
               wgx = wg1
            END IF
            QQPd = Q2Sca(nn,2)
            IF ( IDPdg2.EQ.22 ) THEN
               IF ( IPAmdl(115).EQ.1 ) THEN
                  IF ( QQPd.LT.PVIrtp(2)+PARmdl(144) ) THEN
                     wg1 = 0.D0
                  ELSE
                     wg1 = LOG(QQPd/(PVIrtp(2)+PARmdl(144)))
     &                     /LOG(QQPd/PARmdl(144))
                  END IF
                  IF ( NINhd(nn,2).EQ.0 ) wg1 = wg1*wg1
               ELSE IF ( IPAmdl(115).EQ.2 ) THEN
                  CALL PHO_PDF(2,X0Hd(nn,2),QQPd,PVIrtp(2),pd1)
                  wg1 = pd1(NINhd(nn,2))/PDFva(nn,2)
               END IF
               wgx = wgx*wg1
            END IF

            IF ( IDEb(24).GE.25 .AND. LPRi.GT.4 )
     &            WRITE (LO,'(1X,2A,/5X,I10,I3,1P6E10.3)')
     &            'PHO_POMSCA: ' ,
     &           ' re-weight with (EVE, MSPR, X1/2, Q2, PV1/2, W1/W2)' ,
     &           KEVent , MSPr , X0Hd(nn,1) , X0Hd(nn,2) , QQPd ,
     &           PVIrtp , wgx

            IF ( wgx.LT.DT_RNDM(wgx) ) THEN
               IF ( l.NE.1 ) GOTO 300
               Irej = 50
               RETURN
            END IF

            IF ( LPRi.GT.4 .AND. wgx.GT.1.D0 )
     &            WRITE (LO,'(1X,2A,/5X,I10,I3,1P6E10.3)')
     &            'PHO_POMSCA: ' ,
     &           'weight >1 (EVE, MSPR, X1/2, Q2, PV1/2, W1/W2)' ,
     &           KEVent , MSPr , X0Hd(nn,1) , X0Hd(nn,2) , QQPd ,
     &           PVIrtp , wgx

         END IF

C  generate ISR
         IF ( (ISWmdl(8).GE.2) .AND. ((IPAmdl(101).NE.1) .OR. (l.EQ.1))
     &        ) THEN
            IF ( IPAmdl(109).EQ.1 ) THEN
               q2h = PARmdl(93)*PTHd(nn)**2
            ELSE
               q2h = -PARmdl(93)*VHD(nn)*XHD(nn,1)*XHD(nn,2)*ECMp*ECMp
            END IF
            xhmax1 = 1.D0 - xsss1 - mscha*AS + XHD(nn,1)
            xhmax2 = 1.D0 - xsss2 - mscha*AS + XHD(nn,2)
            i3 = 8*nn - 4
            DO j = 1 , 4
               p1(j) = PPH(i3+j,1)
               p2(j) = PPH(i3+j,2)
            END DO
            IF ( IDEb(24).GE.10 .AND. LPRi.GT.4 )
     &            WRITE (LO,'(1X,A,/5X,2I3,1P,3E12.4)')
     &            'PHO_POMSCA: generate ISR for (L,NN,X1,X2,Q2H)' , l ,
     &           nn , XHD(nn,1) , XHD(nn,2) , q2h
            j = nn
            IF ( l.EQ.1 ) j = -nn
            CALL PHO_HARISR(j,p1,p2,NOUthd(nn,1),NOUthd(nn,2),
     &                      N0Inhd(nn,1),N0Inhd(nn,2),N0Ival(nn,1),
     &                      N0Ival(nn,2),q2h,X0Hd(nn,1),X0Hd(nn,2),
     &                      xhmax1,xhmax2,ifl1,ifl2,NIVal(nn,1),
     &                      NIVal(nn,2),xisr1,xisr2,Irej)
            xsss1 = xsss1 + xisr1 - XHD(nn,1)
            xsss2 = xsss2 + xisr2 - XHD(nn,2)
            NINhd(nn,1) = ifl1
            NINhd(nn,2) = ifl2
            XHD(nn,1) = xisr1
            XHD(nn,2) = xisr2
         END IF

C  check phase space
         IF ( (xsss1.GT.xmaxx1) .OR. (xsss2.GT.xmaxx2) .OR.
     &        ((1.D0-xsss1)*(1.D0-xsss2).LT.AS) ) THEN
            IF ( ihard.EQ.0 ) THEN
               IF ( ISWmdl(2).NE.1 ) GOTO 100
               Mhpom = 0
               Mspom = 1
               Msreg = 0
            END IF
            GOTO 300
         END IF

C  leave energy for leading particle effect
         IF ( (ihard.GT.0) .AND.
     &        ((xsss1.GT.xmaxp1) .OR. (xsss2.GT.xmaxp2)) ) GOTO 300

C  hard scattering accepted
         ihard = ihard + 1
         xss1 = xsss1
         xss2 = xsss2
         IFAil(31) = IFAil(31) - 1

      END DO

C ------------------- end of inner (hard) loop -------------------

 300  Mhpom = ihard
      Mhpar1 = ihard
      Mhpar2 = ihard

C  count valences involved in hard scattering
      Ival1 = 0
      Ival2 = 0
      DO l = 1 , ihard
         nn = LSIdx(l)
         IF ( (NIVal(nn,1).NE.0) .AND. (Ival1.EQ.0) ) Ival1 = nn
         IF ( (NIVal(nn,2).NE.0) .AND. (Ival2.EQ.0) ) Ival2 = nn
      END DO

      iqua1 = 0
      iqua2 = 0
      ivglu1 = 0
      ivglu2 = 0
      DO l = 1 , ihard
         nn = LSIdx(l)

C  photon, pomeron valences
         IF ( (Ival1.EQ.0) .AND. (NINhd(nn,1).NE.0) ) THEN
            IF ( (IDPdg1.EQ.22) .OR. (IDPdg1.EQ.990) ) THEN
               NIVal(nn,1) = 1
               Ival1 = nn
            END IF
         END IF
         IF ( (Ival2.EQ.0) .AND. (NINhd(nn,2).NE.0) ) THEN
            IF ( (IDPdg2.EQ.22) .OR. (IDPdg2.EQ.990) ) THEN
               NIVal(nn,2) = 1
               Ival2 = nn
            END IF
         END IF

C  total number of quarks
         IF ( NINhd(nn,1).NE.0 ) THEN
            iqua1 = iqua1 + 1
         ELSE IF ( ivglu1.EQ.0 ) THEN
            ivglu1 = nn
         END IF
         IF ( NINhd(nn,2).NE.0 ) THEN
            iqua2 = iqua2 + 1
         ELSE IF ( ivglu2.EQ.0 ) THEN
            ivglu2 = nn
         END IF
      END DO

C  gluons emitted by valence quarks
      valpro = 1.D0
      IF ( Ii.EQ.1 ) valpro = VALprg(1)
      ivq1 = 1
      ivg1 = 0
      Ival1 = MAX(Ival1,0)
      IF ( Ival1.EQ.0 ) THEN
         ivq1 = 0
         IF ( (ivglu1.NE.0) .AND. (DT_RNDM(xss1).LT.valpro) ) THEN
            Ival1 = -ivglu1
            ivg1 = 1
         END IF
      END IF
      valpro = 1.D0
      IF ( Ii.EQ.1 ) valpro = VALprg(2)
      ivq2 = 1
      ivg2 = 0
      Ival2 = MAX(Ival2,0)
      IF ( Ival2.EQ.0 ) THEN
         ivq2 = 0
         IF ( (ivglu2.NE.0) .AND. (DT_RNDM(xss2).LT.valpro) ) THEN
            Ival2 = -ivglu2
            ivg2 = 1
         END IF
      END IF
      Mspom = MAX(0,Mspom-iqua1-iqua2)
C  debug output
      IF ( LPRi.GT.4 .AND. IDEb(24).GE.5 ) WRITE (LO,'(1X,A,6I4)')
     &      'PHO_POMSCA: IVAL1/2,IQUA1/2,IVGLU1/2' , Ival1 , Ival2 ,
     &     iqua1 , iqua2 , ivglu1 , ivglu2

C  select soft X values
C  number of soft/remnant quarks
 400  IF ( Mspom.EQ.0 ) THEN
         IF ( IPAmdl(18).EQ.0 ) THEN
            Mspar1 = 2 + 2*Mhpom + Msreg - iqua1 - 2*ivq1 - 2*ivg1
            Mspar2 = 2 + 2*Mhpom + Msreg - iqua2 - 2*ivq2 - 2*ivg2
         ELSE
            Mspar1 = 2 + Msreg + iqua1 - 2*ivq1
            Mspar2 = 2 + Msreg + iqua2 - 2*ivq2
         END IF
      ELSE IF ( IPAmdl(18).EQ.0 ) THEN
         Mspar1 = 2*Mspom + Msreg + 2*Mhpom - iqua1
         Mspar2 = 2*Mspom + Msreg + 2*Mhpom - iqua2
      ELSE
         Mspar1 = 2*Mspom + Msreg + iqua1 + 2*ivg1
         Mspar2 = 2*Mspom + Msreg + iqua2 + 2*ivg2
      END IF
C  debug output
      IF ( LPRi.GT.4 .AND. IDEb(24).GE.15 ) WRITE (LO,'(1X,A,9I3)')
     &      'PHO_POMSCA: MSP,MSR,MHP,IVQ1/2,IVG1/2,MSPAR1/2' , Mspom ,
     &     Msreg , Mhpom , ivq1 , ivq2 , ivg1 , ivg2 , Mspar1 , Mspar2

      xmax1 = 1.D0 - MAX(Mspar1-1,0)*AS - xss1
      xmax2 = 1.D0 - MAX(Mspar2-1,0)*AS - xss2
      i1 = ivq1
      i2 = ivq2
      IF ( Ival1.LE.0 ) i1 = 0
      IF ( Ival2.LE.0 ) i2 = 0
      IF ( (ivq1+ivg1)*(ivq2+ivg2).NE.0 ) THEN
         msdiff = 2*Mspom
      ELSE
         msdiff = 2*MAX(0,Mspom-1)
      END IF
      msg1 = Mspar1
      msg2 = Mspar2
      msm1 = Mspar1 - msdiff
      msm2 = Mspar2 - msdiff
      xmaxh1 = MIN(xmax1,PARmdl(44))
      xmaxh2 = MIN(xmax2,PARmdl(44))
      CALL PHO_SOFTXX(NPOsp(1),NPOsp(2),msg1,msg2,i1,i2,msm1,msm2,xss1,
     &                xss2,xmaxh1,xmaxh2,XS1,XS2,Irej)

C  correct for proper simulation of high pt tail
      IF ( Irej.NE.0 ) THEN
         IF ( LPRi.GT.4 .AND. IDEb(48).GE.2 ) WRITE (LO,'(1X,A,4I4)')
     &         'PHO_STDPAR: rejection (PHO_SOFTXX): MSPOM,MHPOM,I1,I2' ,
     &        Mspom , Mhpom , i1 , i2
         IF ( Mspom*Mhpom.GT.0 ) THEN
            Mspom = Mspom - 1
            GOTO 400
         ELSE IF ( Mspom.GT.1 ) THEN
            Mspom = Mspom - 1
            GOTO 400
         ELSE IF ( Mhpom.GT.1 ) THEN
            ihard = ihard - 1
            IF ( (IPAmdl(13).GT.0) .AND. (IPOix3.EQ.0) .AND.
     &           (IPRoce.EQ.1) ) THEN
               xss1 = MAX(0.D0,1.D0-XPSub)
               xss2 = MAX(0.D0,1.D0-XTSub)
            ELSE
               xss1 = 0.D0
               xss2 = 0.D0
            END IF
            DO k = 1 , ihard
               i = LSIdx(k)
               xss1 = xss1 + XHD(i,1)
               xss2 = xss2 + XHD(i,2)
            END DO
            GOTO 300
         END IF
         Irej = 4
         GOTO 500
      END IF
C  accepted
      Mspom = Mspom - (Mspar1-msg1)/2
      Mspar1 = msg1
      Mspar2 = msg2
C  ------------ kinematics sampled ---------------
C  debug output
      IF ( IDEb(24).GE.10 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,I3)')
     &         'PHO_POMSCA: soft x values, ITRY' , itry
         DO i = 2 , MAX(Mspar1,Mspar2)
            IF ( LPRi.GT.4 ) WRITE (LO,'(10X,I3,2E12.3)') i , XS1(i) ,
     &           XS2(i)
         END DO
      END IF
      IF ( (1.D0-xss1)*(1.D0-xss2).LT.AS ) GOTO 100

C  end of loop
      XS1(1) = 1.D0 - xss1
      XS2(1) = 1.D0 - xss2

C  process counting
      DO n = 1 , LSCahd
         MH_acc_1(NPRohd(n),Ii,IDXmpar) = MH_acc_1(NPRohd(n),Ii,IDXmpar)
     &      + 1
      END DO

C  soft particle momenta

      IF ( MAX(Mspar1,Mspar2).GT.MAXSOF ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,3I4)')
     &         'PHO_POMSCA: no space left in ' ,
     &        '/POSOFT/ (MSPAR1/2,MAXSOF):' , Mspar1 , Mspar2 , MAXSOF
         Irej = 1
         RETURN
      END IF

      DO i = 1 , Mspar1
         PSOft1(1,i) = 0.D0
         PSOft1(2,i) = 0.D0
         PSOft1(3,i) = XS1(i)*ECMp/2.D0
         PSOft1(4,i) = XS1(i)*ECMp/2.D0
      END DO
      DO i = 1 , Mspar2
         PSOft2(1,i) = 0.D0
         PSOft2(2,i) = 0.D0
         PSOft2(3,i) = -XS2(i)*ECMp/2.D0
         PSOft2(4,i) = XS2(i)*ECMp/2.D0
      END DO

      KSOft = MAX(Mspar1,Mspar2)
      KHArd = MAX(Mhpar1,Mhpar2)
      KSPom = Mspom
      KSReg = Msreg
      KHPom = Mhpom

C  debug output
      IF ( IDEb(24).GE.10 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,A,2I3,2I5)')
     &         'PHO_POMSCA: accepted IVAL1,IVAL2,ITRY,NTRY' , Ival1 ,
     &        Ival2 , itry , ntry
         IF ( Mspar1+Mspar2.GT.0 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A)')
     &            'soft x particle1   particle2:'
            xtmp1 = 0.D0
            xtmp2 = 0.D0
            DO i = 1 , MAX(Mspar1,Mspar2)
               IF ( i.LE.MIN(Mspar1,Mspar2) ) THEN
                  IF ( LPRi.GT.4 ) WRITE (LO,'(10X,I3,2E13.4)') i ,
     &                 XS1(i) , XS2(i)
                  xtmp1 = xtmp1 + XS1(i)
                  xtmp2 = xtmp2 + XS2(i)
               ELSE IF ( i.LE.Mspar1 ) THEN
                  IF ( LPRi.GT.4 ) WRITE (LO,'(10X,I3,2E13.4)') i ,
     &                 XS1(i) , 0.D0
                  xtmp1 = xtmp1 + XS1(i)
               ELSE IF ( i.LE.Mspar2 ) THEN
                  IF ( LPRi.GT.4 ) WRITE (LO,'(10X,I3,2E13.4)') i ,
     &                 0.D0 , XS2(i)
                  xtmp2 = xtmp2 + XS2(i)
               END IF
            END DO
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A,2E13.4)')
     &            'sum X1/2 (soft):' , xtmp1 , xtmp2
         END IF
         IF ( Mhpar1.GT.0 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A)')
     &        'NR  IDX  MSPR hard X / hard X ISR / flavor particle 1,2:'
            DO k = 1 , Mhpar1
               i = LSIdx(k)
               IF ( LPRi.GT.4 ) WRITE (LO,'(5X,3I3,4E12.3,2I3)') k , i ,
     &              NPRohd(i) , X0Hd(i,1) , X0Hd(i,2) , XHD(i,1) ,
     &              XHD(i,2) , NINhd(i,1) , NINhd(i,2)
               xtmp1 = xtmp1 + XHD(i,1)
               xtmp2 = xtmp2 + XHD(i,2)
            END DO
            IF ( LPRi.GT.4 ) WRITE (LO,'(1X,A,2E13.4)')
     &            'sum X1/2 (soft+hard):' , xtmp1 , xtmp2
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A)')
     &            'hard momenta  particle1:'
            DO k = 1 , Mhpar1
               i = LSIdx(k)
               i3 = 8*i - 4
               IF ( LPRi.GT.4 ) WRITE (LO,'(5X,2I3,1P,4E12.3,I5)') k ,
     &              i , (PPH(i3+l,1),l=1,4) , NOUthd(i,1)
            END DO
            IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A)')
     &            'hard momenta  particle2:'
            DO k = 1 , Mhpar2
               i = LSIdx(k)
               i3 = 8*i - 4
               IF ( LPRi.GT.4 ) WRITE (LO,'(5X,2I3,1P,4E12.3,I5)') k ,
     &              i , (PPH(i3+l,2),l=1,4) , NOUthd(i,2)
            END DO
         END IF
      END IF
      RETURN

C  event rejected, print debug information
 500  IFAil(4) = IFAil(4) + 1
      IF ( IDEb(24).GE.2 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,/,10X,7I5)')
     &        'PHO_POMSCA: ' ,
     &        'rejection (MSPOM,MHPOM,IHARD,MHARD,ITRY,NTRY,IREJ)' ,
     &        Mspom , Mhpom , ihard , mhard , itry , ntry , Irej
         IF ( LPRi.GT.4 ) WRITE (LO,'(5X,A,I4,1P,2E12.4)')
     &         'IP,Ecm,PTcut:' , Ii , ECMp , PTWant
         IF ( IDEb(24).GE.5 ) THEN
            CALL PHO_PREVNT(0)
         ELSE
            CALL PHO_PREVNT(-1)
         END IF
      END IF

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET