text:
! Copyright (c) 2018 Etienne Descamps
! All rights reserved.
!
! 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. Neither the name of the copyright holder nor the names of its contributors may be
!    used to endorse or promote products derived from this software without specific prior
!    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 COPYRIGHT HOLDER OR 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.

Module test_cstr_kmc_model
  Use ieee_arithmetic
  Use iso_c_binding
  Use mlf_intf
  Use mlf_rsc_array
  Use mlf_models
  Use mlf_step_algo
  Use mlf_fun_intf
  Use mlf_utils
  Use mlf_hybrid_kmc
  Use mlf_supervised_model
  Use mlf_ode45
  Use mlf_ode_class
  Use iso_fortran_env
  IMPLICIT NONE
  PRIVATE

  Type, Public, Extends(mlf_model_real_parameters) :: test_cstr_kmc_parameters
    real(c_double) :: Alpha, Beta, Delta, Kappa, Zeta, Volume
  Contains
    procedure :: set => test_cstr_kmc_parameters_set
    procedure :: getNParameters => test_cstr_kmc_parameters_getNParameters
  End Type test_cstr_kmc_parameters

  Type, Public, Extends(mlf_model_experiment) :: test_cstr_kmc_experiment
    real(c_double) :: a, b, c, d ! Element concentration
    real(c_double) :: T
  End Type test_cstr_kmc_experiment

  ! Values: X -> [c, d, T]
  ! T (with constaint T > 0) -> when T == 0 create A + increase T (beta/Volume)
  ! c, d
  ! c catalyse A -> B (alpha)
  ! reaction c -> d (kappa) decrease T (-delta)
  ! A catalyse d -> c (zeta) increase T (delta)
  Type, Public, Extends(mlf_hybrid_kmc_cstrModel) :: model_cstr_kmc
    integer(c_int64_t), pointer :: NIndiv(:) ! [A, B]
    real(c_double), pointer :: Alpha, Beta, Delta, Kappa, Zeta, Volume
  Contains
    procedure :: init => test_hybrid_init
    procedure :: setupModel => test_setupModel
    procedure :: applyAction => test_applyAction
    procedure :: funTransitionRates => test_funTransitionRates
    procedure :: evalOde => test_evalOde
    procedure :: m_getHMax => test_getHMax
    procedure :: m_updateCstr => test_updateCstr
    procedure :: m_reachCstr => test_reachCstr
    procedure :: m_getDerivatives => test_getDerivatives
  End Type model_cstr_kmc
Contains
  Subroutine test_cstr_kmc_parameters_getNParameters(this, nPar, nCstr)
    class(test_cstr_kmc_parameters), intent(in), target :: this
    integer, intent(out) :: nPar, nCstr
    nPar = 5
    nCstr = 0
  End Subroutine test_cstr_kmc_parameters_getNParameters

  Integer Function test_cstr_kmc_parameters_set(this, X, Cstr) Result(ret)
    class(test_cstr_kmc_parameters), intent(inout), target :: this
    real(c_double), intent(in) :: X(:)
    real(c_double), intent(out), optional :: Cstr(:)
    this%Alpha = X(1)
    this%Beta  = X(2)
    this%Delta = X(3)
    this%Kappa = X(4)
    this%Zeta  = X(5)
    ret = 0
  End Function test_cstr_kmc_parameters_set

  Integer Function test_setupModel(this, param, experiment) Result(info)
    class(model_cstr_kmc), intent(inout), target :: this
    class(mlf_model_parameters), intent(in) :: param
    class(mlf_model_experiment), intent(in) :: experiment
    real(c_double) :: X0(3)
    info = -1
    Select Type(param)
    Class is (test_cstr_kmc_parameters)
      this%Volume = param%Volume
      this%Alpha  = param%Alpha
      this%Beta   = param%Beta
      this%Delta  = param%Delta
      this%Kappa  = param%Kappa
      this%Zeta   = param%Zeta
    Class Default
      RETURN
    End Select
    Select Type(experiment)
    Class is (test_cstr_kmc_experiment)
      this%NIndiv(1) = INT(experiment%a*this%Volume, KIND=8)
      this%NIndiv(2) = INT(experiment%b*this%Volume, KIND=8)
      X0(1) = experiment%c
      X0(2) = experiment%d
      X0(3) = experiment%T
      info = this%initModel(X0)
    Class Default
      RETURN
    End Select
  End Function test_setupModel

  Integer Function test_hybrid_init(this, data_handler) Result(info)
    class(model_cstr_kmc), intent(inout), target :: this
    class(mlf_data_handler), intent(inout), optional :: data_handler
    type(mlf_step_numFields) :: numFields
    class(mlf_ode_algo), pointer :: ode
    integer(c_int64_t) :: NCat
    CALL numFields%initFields(nRPar = 6, nRsc = 1)
    If(info < 0) RETURN
    ALLOCATE(mlf_ode45_obj :: ode)
    If(.NOT. PRESENT(data_handler)) Then
      Select Type(ode)
      Class is (mlf_ode45_obj)
        info = ode%init(4_8, atoli = 1d-6, rtoli = 1d-6)
        If(info < 0) RETURN
      End Select
    Endif
    info = mlf_hybrid_kmc_h_init(this, numFields, ode, &
      numCstr = 1, NActions = 1, data_handler = data_handler)
    If(info < 0) RETURN
    NCat = 2
    info = this%add_i64array(numFields, NCat, this%NIndiv, C_CHAR_"NIndiv", &
      data_handler = data_handler)
    If(info < 0) RETURN
    CALL this%addRPar(numFields, this%Alpha, "Alpha")
    CALL this%addRPar(numFields, this%Beta, "Beta")
    CALL this%addRPar(numFields, this%Delta, "Delta")
    CALL this%addRPar(numFields, this%Kappa, "Kappa")
    CALL this%addRPar(numFields, this%Zeta, "Zeta")
    CALL this%addRPar(numFields, this%Volume, "Volume")
  End Function test_hybrid_init

  Integer Function test_funTransitionRates(this, t, X, F, Rates) Result(N)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(in) :: t
    real(c_double), intent(in), target :: X(:), F(:)
    real(c_double), intent(out), target :: Rates(:) ! [A->B, B->A]
    ASSOCIATE(A => REAL(this%NIndiv(1)), d => X(2), zeta => this%Zeta)
      Rates(1) = zeta*A*d
    END ASSOCIATE
    N = 1
  End Function test_funTransitionRates

  Integer Function test_evalOde(this, t, X, F) Result(info)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(in) :: t
    real(c_double), intent(in), target :: X(:) ! [c, d]
    real(c_double), intent(out), target :: F(:)
    real(c_double) :: Val
    ASSOCIATE(A => REAL(this%NIndiv(1)), B => REAL(this%NIndiv(2)), c => X(1), &
        d => X(2), kappa => this%Kappa, zeta => this%Zeta, V => this%Volume, delta => this%Delta)
      Val = zeta*d*B/V - kappa*c
      F(1) = Val
      F(2) = -Val
      F(3) = Val*delta
    END ASSOCIATE
    info = 0
  End Function test_evalOde

  Integer Function test_applyAction(this, id, t, X, F, Rate) Result(info)
    class(model_cstr_kmc), intent(inout), target :: this
    integer, intent(in) :: id
    real(c_double), intent(in) :: Rate
    real(c_double), intent(inout) :: t
    real(c_double), intent(inout), target :: X(:), F(:)
    this%NIndiv = this%NIndiv + [-1, +1]
    info = 0
  End Function test_applyAction

  Real(c_double) Function test_getHMax(this, t, X, F) Result(hMax)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(in) :: t
    real(c_double), intent(in), target :: X(:), F(:)
    hMax = HUGE(hMax)
    If(X(3) > 0 .AND. F(3) > 0) hMax = 2d0*this%kmc_alpha*X(3)/F(3)
  End Function test_getHMax

  Integer Function test_reachCstr(this, t, tMin, tMax, ids, X, F) Result(info)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(inout) :: t
    real(c_double), intent(in) :: tMin, tMax
    integer, intent(in) :: ids(:)
    real(c_double), intent(inout), target :: X(:), F(:)
    real(c_double) :: dt
    dt = -X(3)/F(3)
    dt = MIN(MAX(dt, tMin-t), tMax-t)
    t = t + dt
    X = X + dt*F
    X(3) = this%Beta/this%Volume
    this%NIndiv(1) = this%NIndiv(1)+1
    info = mlf_ODE_Continue
  End Function test_reachCstr

  Integer Function test_updateCstr(this, t, X0, X, F0, F, ids, hMax) Result(N)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(in) :: t
    real(c_double), intent(in), target :: X0(:), X(:), F0(:), F(:)
    integer, intent(out), target :: ids(:)
    real(c_double), intent(inout) :: hMax
    If(X(3) < 0) Then
      ids(1) = 1
      N = 1
    Else
      N = 0
    Endif
  End Function test_updateCstr

  Subroutine test_getDerivatives(this, ids, X0, X, K, C0, C, Q)
    class(model_cstr_kmc), intent(inout), target :: this
    real(c_double), intent(in), target :: K(:,:), X0(:), X(:)
    real(c_double), intent(out), target :: C0(:), C(:), Q(:,:)
    integer, intent(in), target :: ids(:)
    C0(1) = X0(3)
    C(1) = X(3)
    Q(1,:) = K(3,:)
  End Subroutine test_getDerivatives
End Module test_cstr_kmc_model

max_stars_repo_name: Etdescamps/MlFortran

max_issues_repo_name: Etdescamps/MlFortran

max_forks_repo_name: Etdescamps/MlFortran


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20081029 14:39:59 nbsp Welcome to the Wiki Howdy Steph and Welcome to the Wiki!     My names Evan, pleased to meet you!  Thanks for your review... what else do you like about Davis? Users/JabberWokky Evan JabberWokky Edwards

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c  Sets pa(i)=a(iperm(i)) i=1,2,...,n for logical a and pa.
c  The vectors a and pa can share the same storage locations

      subroutine permutl( a, iperm, n, pa )

      logical    a(n), pa(n), asave
      integer    iperm(n)

      if( n.gt.0 ) then

         do 100 i = 1, n
            iperm(i) = -iperm(i)
  100    continue

         kbeg = 0
  200    kbeg = kbeg + 1

         if( kbeg.le.n ) then

            if( iperm(kbeg).lt.0 ) then

               asave = a(kbeg)
               k = kbeg
  300          l = -iperm(k)

               if( l.eq.kbeg ) then
                  pa(k) = asave
                  iperm(k) = l
                  goto 200
               else
                  pa(k) = a(l)
                  iperm(k) = l
                  k = l
                  goto 300
               endif

            else
               goto 200

            endif

         endif

      endif

      return
      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
! as_mdt_def.for - Aligned Structure for Multiple Data Types
!
! low_byte_def.for must be included before including as_mdt_def.for
!	include		'low_byte_def.for'
!
$IF ABSOFT_FORTRAN
	integer*4	max_ints
	integer*4	max_floats
	integer*4	max_doubles
	integer*4	max_bytes
	parameter	(max_ints = 4096*2)
	parameter	(max_floats = max_ints)
	parameter	(max_doubles = max_ints/2)
	parameter	(max_bytes = max_ints*4)
$ELSE
	parameter	max_ints = 4096*2
	parameter	max_floats = max_ints
	parameter	max_doubles = max_ints/2
	parameter	max_bytes = max_ints*4
$ENDIF
	structure /multi_data_types/
	   union
	   map
	      integer*4	i4(max_ints)
	   end map
	   map
	      integer*2	i2(max_ints*2)
	   end map
	   map
	      real*4	r4(max_floats)
	   end map
	   map
	      real*8	r8(max_doubles)
	   end map
	   map
	      byte	b(max_bytes)
	   end map
	   map
	      character*(max_bytes)	c
	   end map
	   end union
	end structure

	structure /multi_type/
	   union
	   map
	      integer*4	i4
	      integer*4	i4_2nd
	   end map
	   map
	      integer*2	i2
	      integer*2	i2_2nd
	      integer*2	i2_3rd
	      integer*2	i2_4th
	   end map
	   map
	      real*4	r4
	      real*4	r4_2nd
	   end map
	   map
	      real*8	r8
	   end map
	   map
	      byte	b1,b2,b3,b4
	   end map
	   map
	      record /low_byte/ lb
	   end map
	   end union
	end structure

ext: for


text:
	PROGRAM COMPARE_FLUX3
C
C	COMPARE ACCURATE FLUXES FROM STUART AND CDF DATA
C
	include		'wind_examples:wind_tm_routine_def.for'
	include		'wind_examples:wind_return_code_def.for'
C
	CHARACTER*50 	FILENAME
	character*80	stream
	character*32	item
	character*4	event
	INTEGER*4 	DATE,YYYYMMDD,ANGLE
	integer*4	ok,okt,return_size,CDFCH
	real*4		etemp(100),EFLUX(100),dens3dp(100),swedens(100)
	real*8		scetcdf,SWETIME(100),TIME3DP(100),FLUXTIME(100)
	REAL*8		NPTIME,TDIFF,TRIAL
	data event /'CDF'/
	DATA TWOPI /6.28318531/
	DATA BOLTZK /1.3807E-16/
	DATA CLIGHT /2.9979E10/
	DATA EVME /.511E6/
	DATA GMME / 9.10939E-28/
c
	OPEN(UNIT=72,FILE='COMPAREFLUX.RESULTS',ACCESS='APPEND')
C
 10	write(6,6)
  6	format(1x,'Enter TM stream name: ',$)
c	read(5,5,err=10,end=200) iq, stream
	read(5,*,err=10,end=200) yyyymmdd
  5	format(q,a)
	print*,'input ',yyyymmdd
	write(stream,1001) yyyymmdd
 1001	format('wi_lz_wav_',i8.8,'_v*.dat')
	print*,'stream ',stream
	ok = w_channel_open(cdfch,stream)
	if (ok.ne.1) stop 'Cannot open cdf channel'
C
	scetcdf = 0.
	call w_channel_position(cdfch,scetcdf)
	print*,'cdf file starts at scetcdf',scetcdf
	ndayfl = scetcdf
c	now add nearly two hours, as stuart's files start at 0200
	scetcdf = scetcdf + 1.9d00/24.d00
	call w_channel_position(cdfch,scetcdf)
	print*,'cdf file set to',scetcdf
C
 1000 	FORMAT(A)
	IS = 0
C
 100	CONTINUE
	  ok = w_event(cdfch,'CDF')
	  if (ok.ne.1) type*, 'cannot get CDF event,ok=',ok
	  item = 'WIND_SWE_DENSITY_R4'
	  ok = w_item_r4(cdfch, item, swedens, 100, return_size)
	  item = 'WIND_SWE_SCET_R8'
	  ok = w_item_r8(cdfch, item, swetime, 100, return_size)
	  item = 'WIND_3DP_ION_DENSITY_R4'
	  ok = w_item_r4(cdfch, item, dens3DP, 100, return_size)
	  item = 'WIND_3DP_E_TEMP_R4'
	  ok = w_item_r4(cdfch, item, ETEMP, 100, return_size)
	  item = 'WIND_3DP_SCET_R8'
	  ok = w_item_r8(cdfch, item, time3dp, 100, return_size)
C
	do is = 1,50
	  write(66,*) time3dp(is),dens3dp(is),etemp(is)
	enddo
C	  if(is.gt.20) stop
C	  go to 100
C
c	FLUX = 1.68E7*FLUX
C
C
C	THIS SECTION READS IN STUART'S DATA AND FINDS THE CLOSEST
C	OF THE CDF DATA
C
	WRITE(FILENAME,123) YYYYMMDD
C 123	FORMAT('MONTHLY:[KELLOGG.EFLUX]FLUX_',I8.8,'.ASC')
 123	FORMAT('USER_A:[KELLOGG.WIND.PAPERS]FLUX_',I8.8,'.ASC')
	OPEN(UNIT=83,FILE=FILENAME,STATUS='OLD',READONLY)
C
	WRITE(FILENAME,124) YYYYMMDD
 124	FORMAT('USER_A:[KELLOGG.WIND.PAPERS]NP_',I8.8,'.ASC')
C	OPEN(UNIT=84,FILE=FILENAME,STATUS='OLD',READONLY)
C
	NSS = 0
	DO NS =  1,50
	  READ(83,*,END=140,ERR=150) NYR,MO,NDAY,NHR,MIN,NSEC,EFLUXT
	  NSS = NSS+1
	  FDAYF = NHR/24. + MIN/1440. + NSEC/86400.
C	  print*,'nyr,MO,NDAY,NSEC,eflux',nyr,MO,NDAY,NSEC,efluxt
     	  FLUXTIME(NSS) = DFLOAT(NDAYFL) + DBLE(FDAYF)
	  EFLUX(NSS) = EFLUXT
 150	  CONTINUE
	ENDDO
 140	continue
	print*,' lines read in',nss,'  last ',fdayf
C
C	SEARCH CDF DATA (NCDF = 1,100) FOR CLOSEST TIME
C
	DO NS = 1,NSS
	  TDIFF = 1.D00
	  DO NCDF = 1,100
C	  print*,'time3dp,flux',time3dp(ncdf),ns,fluxtime(ns)
	    IF(TIME3DP(NCDF).GT.FLUXTIME(NS)) GO TO 160
	    TRIAL = ABS(FLUXTIME(NS)-TIME3DP(NCDF))
	    IF(TRIAL.LT.TDIFF) THEN
		TDIFF = TRIAL
	        NBEST = NCDF
	    ENDIF
	  ENDDO
 160	  CONTINUE
C
C	  NOW NBEST IS THE CDF DATA NEARET STUART'S FLUX

	PRINT*,'ns,fluxtime,ncdf,time3dp ',ns,FLUXTIME(ns)
     1		,ncdf,time3dp(ncdf)
	  IF(DENS3DP(NBEST).LE.0.) GO TO 162
	  XTEMP = EFLUX(NS)/DENS3DP(NBEST)
	  EFFTEMP = TWOPI*EVME*(XTEMP/CLIGHT)**2	! IN EV
	  RATIO = 1.
	  IF(ETEMP(NBEST).GT.0.) RATIO = EFFTEMP/ETEMP(NBEST)
	  WRITE(72,1072) YYYYMMDD,FLUXTIME(NS),DENS3DP(NBEST),EFLUX(NS),
     1	    ETEMP(NBEST),EFFTEMP,RATIO
 1072	  FORMAT(I9,F12.5,F8.3,E12.4,E12.4,F10.4,E12.4)
  162	  CONTINUE
	ENDDO
C
C	READ(84,*,END=141,ERR=140) NYR,MO,NDAY,NHR,MIN,NSEC,DENSI
C	print*,'nyr,MO,NDAY,NSEC,DENSI',nyr,MO,NDAY,NSEC,DENSI
C	FDAYN = NHR/24. + MIN/1440. + NSEC/86400.
C	NPTIME = 365.D00*DFLOAT(NYR-1982) + 30.D00*DFLOAT(MO) +
C     1	  DFLOAT(NDAY) + DBLE(FDAYN)
C	PRINT*,FLUXTIME,NPTIME
C	IF(ABS(FDAYF-FDAYN).GT.1.E-7) STOP 'UNEQUAL FLUX AND NP TIMES'
C
 141	CONTINUE
	OK = W_CHANNEL_CLOSE(CDFCH)
	GO TO 10
C
 200	CONTINUE
	PRINT*,'END OF .COM INPUT'
	CLOSE(UNIT=83)
C	CLOSE(UNIT=84)
	CLOSE(UNIT=72)
	STOP
	END


text:
SUBROUTINE crust(crustf_rm,tillf,crustf,lmosf,lmosm)
!
IMPLICIT NONE
!
! Subroutine arguments
!
REAL :: crustf,crustf_rm,lmosf,lmosm,tillf
!
!     + + + purpose + + +
!
!     this subroutine destroys the surface crust after a tillage event.
!
!     + + + keywords + + +
!     crust, tillage (primary/secondary)

!     include 'p1werm.inc'
!
!     + + + argument declarations + + +
!
!
!     + + + argument definitions + + +
!
!     lmosf - fraction of crusted surface containing loose erodible material
!     lmosm - mass of loose erodible material on crusted portion of surface
!     crustf - current fraction of surface crusted (before & after operation)
!     crustf_rm - fraction of crust removed (0 <= crustf_rm <= 1)
!     tillf - fraction of the surface tilled (0 <= tillf <= 1)
!
!     + + + accessed common block variable definitions + + +
!
!
!     + + + parameters + + +
!
!     + + + local variables + + +
!
!     + + + local variable definitions + + +
!
!     + + + end specifications + + +
!
!
!     crf = cri * ( (1.0 - tillf) + (tillf * (1.0-crustf_rm)))

      ! determine fraction of surface that remains crusted
crustf = crustf*(1.0-tillf*crustf_rm)

! currently the crust function doesn't modify the loose erodible
! material variables on the crusted surface.  that could be changed
! in the future if it was deemed necessary.

! the following should be removed.  need to check soil and erosion
! first to make sure they aren't adversely affected. - lew
! 8/25/1999

!     check to see if the loose material on the surface is still there
!     if enough of the crust is removed set lmosf to zero (loose material)
!     this was done according to l. hagen

      ! just clear them out if it close to zero
      ! (lh shouldn't have erosion or soil submodels this sensitive)
IF (crustf.LT.0.01) THEN
  lmosf = 0.0
  lmosm = 0.0
END IF
!
END SUBROUTINE crust

max_stars_repo_path: Project Documents/Source Code Original/Crust.f90

max_stars_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_issues_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_forks_repo_name: USDA-ARS-WMSRU/upgm-standalone


text:
Best Sorority Forever is a philanthropic event presented by the brothers of Delta Lambda Phi.  Each assigned a DLP brother as its coach, the sororities at UC Davis select a representative to Greek Competitive Events compete a fun and sassy beauty pageant.

Winners
2012 6 sororities competing
1st place:Tri Delta
2nd place:Alpha Phi


module zfp
    use zFORp
    use, intrinsic :: iso_c_binding
    implicit none

contains
! Fortran convenience subroutines

    subroutine compress_pixels(pixels, compressed)
        use, intrinsic :: iso_c_binding
        implicit none

        ! inputs
        real(kind=4), dimension(:, :), contiguous, target, intent(in) :: pixels

        ! the output
        character(kind=c_char), allocatable, intent(out) :: compressed(:)

        ! internal variables
        integer :: nx, ny

        ! zfp_field
        type(zFORp_field) :: field

        ! bitstream
        character(kind=c_char), allocatable, target :: buffer(:)
        integer(kind=8) buffer_size_bytes, bitstream_offset_bytes
        type(zFORp_bitstream) :: bitstream, queried_bitstream

        ! zfp_stream
        type(zFORp_stream) :: stream
        real(kind=8) :: desired_rate, rate_result
        integer :: dims, wra
        integer :: zfp_type

        nx = ubound(pixels, 1)
        ny = ubound(pixels, 2)

        print *, 'nx=', nx, ', ny=', ny

        ! setup zfp_field
        zfp_type = zFORp_type_float
        field = zFORp_field_2d(c_loc(pixels), zfp_type, nx, ny)

        print *, 'got here#0'

        ! in the worst case assume no compression
        buffer_size_bytes = sizeof(pixels)
        print *, 'buffer_size_bytes=', buffer_size_bytes
        allocate (buffer(buffer_size_bytes))
        bitstream = zFORp_bitstream_stream_open(c_loc(buffer), buffer_size_bytes)
        print *, 'got here#1'

        ! setup zfp_stream
        stream = zFORp_stream_open(bitstream)
        print *, 'got here#2'

        desired_rate = 8.0
        dims = 2
        wra = 0
        zfp_type = zFORp_type_float
        rate_result = zFORp_stream_set_rate(stream, desired_rate, zfp_type, dims, wra)
        print *, 'got here#3'

        queried_bitstream = zFORp_stream_bit_stream(stream)
        print *, 'got here#4'

        ! compress
        bitstream_offset_bytes = zFORp_compress(stream, field)
        write (*, *) "After compression, bitstream offset at "
        write (*, *) bitstream_offset_bytes

        ! deallocations
        call zFORp_stream_close(stream)
        call zFORp_bitstream_stream_close(queried_bitstream)
        call zFORp_field_free(field)

        if (bitstream_offset_bytes .gt. 0) compressed = reshape(buffer, (/bitstream_offset_bytes/))
    end subroutine compress_pixels

end module zfp


text:
CCC   'station.f' RECORDS SEISMOGRAMS AND WAVEFIELDS

      subroutine getrec(coords,maxdim,rank,nxt,nyt,nzt,nz,recproc,nprec,nrec,
     +           nbgx,nedx,nskpx,nbgy,nedy,nskpy,nbgz,nedz,nskpz,casenumber)

CCC   DETERMINES STATION DISTRIBUTION AMONG PROCS

      use parstat
      integer :: i,j,k
      integer*8 :: r
      integer :: l
      integer :: nprec
      integer(i8) :: nrec
      integer :: nxt,nyt,nzt, nz, rank,recproc, casenumber
      integer :: nbgx,nedx,nskpx,nbgy,nedy,nskpy,nbgz,nedz,nskpz
      integer, dimension(maxdim) :: coords
      integer*8 :: Li, Lj, Lk

C     SEARCH OVER ALL RECEIVER NODES

C     DETERMINE IF NODE BELONGS TO PROC AND COMPUTE FILE POSITION

      l = 1
      r = -1

      do k=nbgz,nedz,nskpz
        if(k==0) then
           nzp = int((nz-1)/nzt);
        else
           nzp = int((nz-k)/nzt);
        end if

        if (coords(3)==nzp) then
          do j=nbgy,nedy,nskpy
            nyp = int((j-1)/nyt)

            if (coords(2)==nyp) then
              do i=nbgx,nedx,nskpx
                nxp = int((i-1)/nxt)

                if (coords(1)==nxp) then
                  a = mod(i,nxt)
                  b = ANINT(a/(a+1))
                  tprec(l,1) = a + (1-b) * nxt

                  a = mod(j,nyt)
                  b = ANINT(a/(a+1))
                  tprec(l,2) = a + (1-b) * nyt

                  if(k == 0) then
                     tprec(l,3) = nzt + 1;
                  else
                     a = mod(nz-k+1,nzt)
                     b = ANINT(a/(a+1))
                     tprec(l,3) = a + (1-b) * nzt
                  end if

C    TAG PROC AS STATION BEARING
c    Calculate r first,
                  Li = int((nedx-nbgx)/nskpx) + 1
                  Lj = int((nedy-nbgy)/nskpy) + 1
                  Lk = int((nedz-nbgz)/nskpz) + 1
                  r = Li*Lj*int((k-nbgz)/nskpz) + Li*int((j-nbgy)/nskpy) + int((i-nbgx)/nskpx)

                  tpmap(l) = r
                  recproc = rank
                  l=l+1
                endif
              end do
            endif
          end do
        endif
      end do

C     SET NUMBER OF STATION NODES IN PROC

      nprec=l-1
      nrec=(int((nedx-nbgx)/nskpx)+1) * (int((nedy-nbgy)/nskpy)+1) * (int((nedz-nbgz)/nskpz)+1)

C     ALLOCATE RECEIVER ARRAYS AND DEALLOCATE TMP ARRAYS
c      allocate(prec(nprec,maxdim))
c      allocate(pmap(nprec))

C     TRANSFER RECEIVER INFO TO PERMANENT ARRAYS
c      mm=1
CC     use casenumber to assign to different prec* and pmap* variables
CC     1: ; 2: original getrec2; 3: for rupt output; 4: for all other 19 variables in SGSN

      SELECT case (casenumber)
        case (1)
          allocate(prec(nprec,maxdim))
     	  allocate(pmap(nprec))

      	  do i=1,nprec
      	    prec(i,1) = tprec(i,1)
      	    prec(i,2) = tprec(i,2)
       	    prec(i,3) = tprec(i,3)

       	    pmap(i) = tpmap(i)
      	  end do

        case (2)
          allocate(prec2(nprec,maxdim))
  	  allocate(pmap2(nprec))

      	  do i=1,nprec
      	    prec2(i,1) = tprec(i,1)
	    prec2(i,2) = tprec(i,2)
       	    prec2(i,3) = tprec(i,3)

       	    pmap2(i) = tpmap(i)
          end do

        case (3)
          allocate(prec3(nprec,maxdim))
   	  allocate(pmap3(nprec))

      	  do i=1,nprec
            prec3(i,1) = tprec(i,1)
       	    prec3(i,2) = tprec(i,2)
       	    prec3(i,3) = tprec(i,3)

       	    pmap3(i) = tpmap(i)
      	  end do

	END SELECT

C     DEALLOCATE TEMPORARY RECEIVER ARRAYS

        deallocate(tprec)
        deallocate(tpmap)

        return
      end

c     wrtrec is not in use
      subroutine wrtrec(ii,ntiskp,fhx,fhy,fhz,nprec,nrec)

CCC   WRITE OUT SEISMOGRAMS IN SPECIFIED VOLUME
      use parstat

      !include 'mpif.h'
c      integer,parameter :: rtype=selected_int_kind( 18 )
      integer :: nprec
      integer :: nrec
c      integer(i8) :: nrec
c      integer(rtype) :: nprec,nrec
      integer :: ii,j, ntiskp,it
      integer :: err,filetype
      integer :: xn,yn,zn
      integer :: fhx,fhy,fhz

      integer (kind=MPI_OFFSET_KIND) :: disp
      integer :: status(MPI_STATUS_SIZE)
      real, dimension(nprec) :: bufx,bufy,bufz

      disp=0

      it = ntiskp*(ii/ntiskp)
      if(ii==it) then
        call MPI_TYPE_CREATE_INDEXED_BLOCK(nprec,1,pmap,MPI_REAL,
     +                                   filetype,err)

        call MPI_TYPE_COMMIT(filetype,err)

        call MPI_FILE_SET_VIEW(fhx,disp,MPI_REAL,filetype,'native',
     +                       MPI_INFO_NULL,err)

        call MPI_FILE_SET_VIEW(fhy,disp,MPI_REAL,filetype,'native',
     +                       MPI_INFO_NULL,err)

        call MPI_FILE_SET_VIEW(fhz,disp,MPI_REAL,filetype,'native',
     +                       MPI_INFO_NULL,err)

        do j=1,nprec

          xn = prec(j,1)
          yn = prec(j,2)
          zn = prec(j,3)

          bufx(j) = u1(xn,yn,zn)
          bufy(j) = v1(xn,yn,zn)
          bufz(j) = w1(xn,yn,zn)

        end do

        call MPI_FILE_WRITE_ALL(fhx,bufx,nprec,MPI_REAL,status,err)
        call MPI_FILE_WRITE_ALL(fhy,bufy,nprec,MPI_REAL,status,err)
        call MPI_FILE_WRITE_ALL(fhz,bufz,nprec,MPI_REAL,status,err)

        pmap = pmap + nrec

      end if
      return
      end


text:
program failed_allocate_test
    use, intrinsic :: iso_fortran_env, only : INT32, INT64
    implicit none
    integer(kind=INT64), parameter :: init_size = 5000, size_step = 5000
    integer(kind=INT32), dimension(:, :), allocatable :: A
    integer :: istat
    integer(kind=INT64) :: size
    real :: total_size

    size = init_size
    do while(.true.)
        allocate(A(size, size), stat=istat)
        total_size = 4.0*size**2/(1024**2)
        if (.not. allocated(A)) &
            print '(A)', 'allocated returns .false.'
        if (istat /= 0) then
            print '(A, F10.2, A)', 'array of size ', total_size, &
                    ' MB not allocated'
            exit
        else
            print '(A, F10.2, A)', 'array of size ', total_size, &
                    ' MB allocated'
        end if
        deallocate(A)
        size = size + size_step
    end do

end program failed_allocate_test

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
C-------GF 28.08.1980------------------------------------- P R O S
      INTEGER*2 FPRO,PROHIB
      INTEGER*2 PRO   (XXPROH),PROLNG,PROLEF(XXPROH)
      INTEGER*2 PROMON(XXPROH),PROSEM(XXPROH)
      COMMON /PROS/
     =  FPRO     ! 1ST FREE PRODUCTION
     = ,PROHIB   ! HIGH BOUND OF 'PROS'
     = ,PRO     ! -> NEXT PRODUCTION OF THE LEFT SIDE
     = ,PROLEF  ! SYMBOL ON THE LEFT SIDE
     = ,PROMON  ! -> FIRST MEMBER
     = ,PROSEM  ! -> SEMANTICS FOR THE PRODUCTION

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


module matrix_utils

  use iso_fortran_env, only : dp => REAL64

contains

  !===================================================================!
  ! Setup Toepliz matrix with supplied columns
  !===================================================================!

  subroutine toepliz(col, A)

    real(dp), allocatable :: A(:,:)
    real(dp) :: col(:)
    integer :: m, i, j

    m = size(A, 1)
    do j = 1, m
       do i = 1, m
          if (i.eq.j) then
             A(i,j) = col(1)
          else if (i+1.eq.j .or.i-1.eq.j) then
             A(i,j) = col(2)
             A(j,i) = col(2)
          else if (i+2.eq.j .or.i-2.eq.j) then
             A(i,j) = col(3)
             A(j,i) = col(3)
          else if (i+3.eq.j .or.i-3.eq.j) then
             A(i,j) = col(4)
             A(j,i) = col(4)
          else if (i+4.eq.j .or.i-4.eq.j) then
             A(i,j) = col(5)
             A(j,i) = col(5)
          end if
       end do
    end do

  end subroutine toepliz

  subroutine print(matrix)

    real(dp), intent(in) :: matrix(:,:)
    integer :: i, j, m, n

    m = size(matrix, 1)
    n = size(matrix, 2)

    do i = 1, m
       write(*,'(100g15.5)') ( matrix(i,j), j=1,n )
    end do
  end subroutine print

end module matrix_utils


text:
Adrienne McCarthy: n.,, 1. A recent graduate of UC Davis in the area of microbiology who now lives and works in the bay area.

She plays DDR with me at our awesome apartment, too. Users/BrentLaabs


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE UMSCHM (IA, NPER, KKK, LLL, NNN, ML, MS, NSPR, ISLIST,
     &   NINT, IFLINE, NLPS, ILLIST, LINKL, LINKS, MXNPER, MAXPRM, MAX3,
     &   MXND, X, Y, NID, NNPS, ANGLE, XN, YN, NUID, LXK, KXL, NXL, LXN,
     &   XSUB, YSUB, NIDSUB, INDX, IAVAIL, NAVAIL, CCW, ERR)
C***********************************************************************

C  UMSCHM - "U" MESH SCHEME; CALCULATE A "PENTAGON" MAPPED MESH
C           (3 RECTANGULAR SUBREGIONS)

C***********************************************************************

      DIMENSION IA(1)
      DIMENSION ISLIST(NSPR), NINT(ML), IFLINE(MS), NLPS(MS)
      DIMENSION ILLIST(MS*3), LINKL(2, ML), LINKS(2, MS)
      DIMENSION X(MXNPER), Y(MXNPER), NID(MXNPER*MAXPRM), NNPS(MAX3)
      DIMENSION ANGLE(MXNPER), XN(MXND), YN(MXND), NUID(MXND)
      DIMENSION LXK(4, MXND), KXL(2, 3*MXND), NXL(2, 3*MXND)
      DIMENSION LXN(4, MXND)
      DIMENSION XSUB(MXNPER), YSUB(MXNPER), NIDSUB(MXNPER), INDX(MXND)

      LOGICAL CCW, ERR, FINAL

C  SET UP THE PENTAGON DIVISIONS, AND FIND THE CENTER POINT

      CALL GETM5 (IA, ML, MS, MAX3, NSPR, ISLIST, NINT, IFLINE, NLPS,
     &   ILLIST, LINKL, LINKS, X, Y, NID, NNPS, ANGLE, NPER, M1A, M1B,
     &   M2, M3A, M3B, M4A, M4B, M5, MC, XCEN, YCEN, CCW, ERR)
      FINAL = .FALSE.

C  SET UP THE FIRST SUBREGION, AND SEND IT OFF TO BE GENERATED

      IF (.NOT.ERR) THEN
         CALL SUBPEN (NPER, NEWPER, X, Y, NID, XSUB, YSUB, NIDSUB,
     &      M1B + M2 + M3A, M4A, MC, M1A, 1, XCEN, YCEN)
         NNNOLD = NNN
         KKKOLD = KKK
         LLLOLD = LLL
         CALL RMESH (NEWPER, MXND, XSUB, YSUB, NIDSUB, XN, YN, NUID,
     &      LXK, KXL, NXL, LXN, M1B, M2, KKK, KKKOLD, NNN, NNNOLD, LLL,
     &      LLLOLD, IAVAIL, NAVAIL, ERR)
      END IF

C  SET UP THE SECOND SUBREGION, AND SEND IT OFF TO BE GENERATED

      IF (.NOT.ERR) THEN
         CALL SUBPEN (NPER, NEWPER, X, Y, NID, XSUB, YSUB, NIDSUB,
     &      M3B + M4A, M3B, M4A, M1A + M1B + M2 + M3A, 2, XCEN, YCEN)
         NNNOLD = NNN
         KKKOLD = KKK
         LLLOLD = LLL
         CALL RMESH (NEWPER, MXND, XSUB, YSUB, NIDSUB, XN, YN, NUID,
     &      LXK, KXL, NXL, LXN, M3B, M4A, KKK, KKKOLD, NNN, NNNOLD, LLL,
     &      LLLOLD, IAVAIL, NAVAIL, ERR)
         CALL FIXSUB (MXND, NNNOLD, NNN, LLLOLD, LLL, KKKOLD, KKK, XN,
     &      YN, NUID, LXK, KXL, NXL, LXN, INDX, IAVAIL, NAVAIL, FINAL)
      END IF

C  SET UP THE THIRD SUBREGION, AND SEND IT OFF TO BE GENERATED

      IF (.NOT.ERR) THEN
         CALL SUBPEN (NPER, NEWPER, X, Y, NID, XSUB, YSUB, NIDSUB,
     &      M4B + M5 + M1A, MC, M3B, M1A + M1B + M2 + M3A + M3B + M4A,
     &      3, XCEN, YCEN)
         NNNOLD = NNN
         KKKOLD = KKK
         LLLOLD = LLL
         CALL RMESH (NEWPER, MXND, XSUB, YSUB, NIDSUB, XN, YN, NUID,
     &      LXK, KXL, NXL, LXN, M4B, M5, KKK, KKKOLD, NNN, NNNOLD, LLL,
     &      LLLOLD, IAVAIL, NAVAIL, ERR)
         FINAL = .TRUE.
         CALL FIXSUB (MXND, NNNOLD, NNN, LLLOLD, LLL, KKKOLD, KKK, XN,
     &      YN, NUID, LXK, KXL, NXL, LXN, INDX, IAVAIL, NAVAIL, FINAL)
      END IF

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
program test_comm_split
  use libmpifx_module
  implicit none

  type(mpifx_comm) :: allproc, groupproc
  integer :: groupsize, mygroup

  call mpifx_init()
  call allproc%init()
  groupsize = allproc%size / 2
  mygroup = allproc%rank / groupsize
  call allproc%split(mygroup, allproc%rank, groupproc)
  write(*, "(3(A,1X,I0,1X))") "GLOBAL ID:", allproc%rank, "SUBGROUP", &
      & mygroup, "SUBGROUP ID", groupproc%rank
  call mpifx_finalize()

end program test_comm_split

max_stars_repo_name: terminationshock/mpifx

max_issues_repo_name: terminationshock/mpifx

max_forks_repo_name: terminationshock/mpifx


text:
Barnyard Theatre is a community Live Theatre theatre group dedicated to creating theatre from scratch. Each summer this group of theatre artists gathers in Davis, CA to develop a new script, and transform a working barn into a nontraditional theatre space. Situated in the Outskirts fields of Yolo County, the unique location of the Historic Schmeiser Barn supports Barnyard Theatre’s innovative staging and the creation of dynamic relationships between actors and audience members. The artists are dedicated to producing highquality theatre in a stressfree environment.

The Barnyard Theatre motto, Think Organic, reflects a dedication to theatre from scratch and a return to Yolo County’s agricultural roots in the Historic Schmeiser Barn.

Barnyard Theatres productions of original works include the premieres of the fulllength shows My Avisia Winger, The Sterling Affair, and A Flag Touched The Ground.

In 2006, Barnyard Theatre partnered with UC Davis Friends of Acoustic Music to present the first Barnyard Acoustic Music Festival in The Historic Schmeiser Barn. The Barnyard Acoustic Music Festival featured four separate evenings of acoustic music on the dark nights of the July/August theatre run, with a different theme highlighted each evening.

 Past Showings

Things that Fall from the Sky: an evening of short plays.
Directed by Lindsay Carpenter.

Featuring:
Unusual Epitaphs by Rob Rinow, Fly Me to the Moon by Eva Konstantopoulos, and Things That Fall From the Sky by Brenda Varda.

20070524 18:29:53 nbsp Theyre so good that you can sometimes forget that theyre just kids, especially when you remember that theyve done everything from writing the script to cleaning up after.  Bring your mosquito repellent and enjoy. Users/NotTires

20100618 10:48:07 nbsp This year they will be doing a series of short plays in a performance entitled Things that Fall from the Sky. Dont miss it!
 Users/EdwardNiemand


module NavierStokes2DM

  use UtilitiesM

  use BaseIntegrandM

  use IntegrandM

  implicit none

  private
  public :: NavierStokes2DDT, SetNavierStokes2D

  type, extends(IntegrandDT) :: NavierStokes2DDT
     real(rkind), dimension(:), allocatable :: lhs
   contains
     procedure :: t          => dnavierStokes2D_dt
     procedure :: add        => add_navierStokes2D
     procedure :: multiply   => multiply_navierStokes2D
     procedure :: assign     => assign_navierStokes2D
     procedure :: getState
  end type NavierStokes2DDT

  interface SetNavierStokes2D
     procedure constructor
  end interface SetNavierStokes2D

contains

  type(navierStokes2DDT) function constructor(initial_state, lhs&
       , this_strategy, step)
    class(BaseIntegrandDT)       , intent(in) :: this_strategy
    real(rkind), dimension(:), intent(in) :: initial_state
    real(rkind), dimension(:), intent(in) :: lhs
    integer(ikind), intent(in), optional  :: step
    constructor%state             = initial_state
    constructor%lhs               = lhs
    if (present(step)) constructor%step = step
    call constructor%set_quadrature(this_strategy)
  end function constructor

  function dNavierStokes2D_dt(this, dof) result(dState_dt)
    class(NavierStokes2DDT)  , intent(in) :: this
    class(IntegrandDT)    , allocatable   :: dState_dt
    type(NavierStokes2DDT), allocatable   :: local_dState_dt
    real(rkind), dimension(:), intent(in) :: dof
    allocate(local_dState_dt)
    call local_dState_dt%set_quadrature(this%get_quadrature())
    allocate(local_dState_dt%state(size(this%state)))
    local_dState_dt%state = this%lhs
    call move_alloc(local_dState_dt,dState_dt)
  end function dNavierStokes2D_dt

  function add_navierStokes2D(lhs,rhs) result(sum)
    class(NavierStokes2DDT), intent(in) :: lhs
    class(IntegrandDT)     , intent(in) :: rhs
    class(IntegrandDT), allocatable     :: sum
    type(NavierStokes2DDT), allocatable :: local_sum
    select type(rhs)
    class is (navierStokes2DDT)
       allocate(local_sum)
       call local_sum%set_quadrature(lhs%get_quadrature())
       local_sum%state = lhs%state + rhs%state
    class default
       stop 'add_navierStokes2D: unsupported class'
    end select
    call move_alloc(local_sum,sum)
  end function add_navierStokes2D

  function multiply_navierStokes2D(lhs,rhs) result(product)
    class(NavierStokes2DDT), intent(in)  :: lhs
    class(IntegrandDT)     , allocatable :: product
    type(NavierStokes2DDT) , allocatable :: local_product
    real(rkind)            , intent(in)  :: rhs
    allocate(local_product)
    call local_product%set_quadrature(lhs%get_quadrature())
    local_product%state = lhs%state * rhs
    call move_alloc(local_product,product)
  end function multiply_navierStokes2D

  subroutine assign_navierStokes2D(lhs,rhs)
    class(NavierStokes2DDT), intent(inout) :: lhs
    class(IntegrandDT)     , intent(in)    :: rhs
    select type(rhs)
    class is (NavierStokes2DDT)
       lhs%state = rhs%state
       lhs%step = rhs%step
       lhs%lhs = rhs%lhs
       call lhs%set_quadrature(rhs%get_quadrature())
    class default
       stop 'assign_navierStokes2D: unsupported class'
    end select
  end subroutine assign_navierStokes2D

  function getState(this) result(coordinates)
    class(NavierStokes2DDT)  , intent(in)  :: this
    real(rkind), dimension(:), allocatable :: coordinates
    coordinates = this%state
  end function getState

end module NavierStokes2DM


text:
*DECK R1MACH
      REAL FUNCTION R1MACH (I)
C***BEGIN PROLOGUE  R1MACH
C***PURPOSE  Return floating point machine dependent constants.
C***LIBRARY   SLATEC
C***CATEGORY  R1
C***TYPE      SINGLE PRECISION (R1MACH-S, D1MACH-D)
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  Fox, P. A., (Bell Labs)
C           Hall, A. D., (Bell Labs)
C           Schryer, N. L., (Bell Labs)
C***DESCRIPTION
C
C   R1MACH can be used to obtain machine-dependent parameters for the
C   local machine environment.  It is a function subprogram with one
C   (input) argument, and can be referenced as follows:
C
C        A = R1MACH(I)
C
C   where I=1,...,5.  The (output) value of A above is determined by
C   the (input) value of I.  The results for various values of I are
C   discussed below.
C
C   R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
C   R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C   R1MACH(3) = B**(-T), the smallest relative spacing.
C   R1MACH(4) = B**(1-T), the largest relative spacing.
C   R1MACH(5) = LOG10(B)
C
C   Assume single precision numbers are represented in the T-digit,
C   base-B form
C
C              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
C   EMIN .LE. E .LE. EMAX.
C
C   The values of B, T, EMIN and EMAX are provided in I1MACH as
C   follows:
C   I1MACH(10) = B, the base.
C   I1MACH(11) = T, the number of base-B digits.
C   I1MACH(12) = EMIN, the smallest exponent E.
C   I1MACH(13) = EMAX, the largest exponent E.
C
C   To alter this function for a particular environment, the desired
C   set of DATA statements should be activated by removing the C from
C   column 1.  Also, the values of R1MACH(1) - R1MACH(4) should be
C   checked for consistency with the local operating system.
C
C***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
C                 a portable library, ACM Transactions on Mathematical
C                 Software 4, 2 (June 1978), pp. 177-188.
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790101  DATE WRITTEN
C   890213  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900618  Added DEC RISC constants.  (WRB)
C   900723  Added IBM RS 6000 constants.  (WRB)
C   910710  Added HP 730 constants.  (SMR)
C   911114  Added Convex IEEE constants.  (WRB)
C   920121  Added SUN -r8 compiler option constants.  (WRB)
C   920229  Added Touchstone Delta i860 constants.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C   920625  Added CONVEX -p8 and -pd8 compiler option constants.
C           (BKS, WRB)
C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
      SAVE RMACH
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C

C     MACHINE CONSTANTS FOR Linux on IBM PC

c Machine variables were incorrect, replaced with decimal representation
c      DATA RMACH(1) / Z'00800000' /
c      DATA RMACH(2) / Z'7F7FFFFF' /
c      DATA RMACH(3) / Z'33800000' /
c      DATA RMACH(4) / Z'34000000' /
c      DATA RMACH(5) / Z'3E9A209B' /

       DATA RMACH(1) / 1.17549435e-38 /
       DATA RMACH(2) / 3.40280562e+38 /
       DATA RMACH(3) / 5.96046448e-08 /
       DATA RMACH(4) / 1.1920929e-07  /
       DATA RMACH(5) / 0.30103001     /

C
C***FIRST EXECUTABLE STATEMENT  R1MACH
      IF (I .LT. 1 .OR. I .GT. 5) CALL XERMSG ('SLATEC', 'R1MACH',
     +   'I OUT OF BOUNDS', 1, 2)
C
      R1MACH = RMACH(I)
      RETURN
C
      END

max_stars_repo_name: sfarrens/cosmostat

max_issues_repo_name: sfarrens/cosmostat

max_forks_repo_name: sfarrens/cosmostat


text:
module PRIMLINEULER_bc
    use mod_kinds,          only: rk,ik
    use mod_constants,      only: ONE,TWO,THREE,FOUR,FIVE,EIGHT,NINE,HALF,ZERO,PI

    use type_operator,      only: operator_t
    use type_chidg_worker,  only: chidg_worker_t
    use type_properties,    only: properties_t
    use DNAD_D

    use mod_primitive_linearized_euler
    implicit none
    private

    !>  BC flux for Linearized Euler equations.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   1/22/2018
    !!
    !!
    !----------------------------------------------------------------------------------
    type, extends(operator_t), public :: PRIMLINEULER_bc_t

    contains

        procedure   :: init
        procedure   :: compute

    end type PRIMLINEULER_bc_t
    !***********************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   1/22/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(PRIMLINEULER_bc_t),  intent(inout)   :: self

        ! Set operator name
        call self%set_name("PRIMLINEULER BC")

        ! Set operator type
        call self%set_operator_type("BC Advective Flux")

        ! Set operator equations
        call self%add_primary_field("Density"   )
        call self%add_primary_field("Velocity-1")
        call self%add_primary_field("Velocity-2")
        call self%add_primary_field("Velocity-3")
        call self%add_primary_field("Pressure"  )

    end subroutine init
    !********************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   1/23/2018
    !!
    !----------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(PRIMLINEULER_bc_t),   intent(inout)   :: self
        type(chidg_worker_t),       intent(inout)   :: worker
        class(properties_t),        intent(inout)   :: prop

        type(AD_D), allocatable, dimension(:)   ::  &
            rho_r, u_r, v_r, w_r, p_r,              &
            flux_1, flux_2, flux_3

        print*, 'Computing BC'

        !
        ! Zero primary matrix container to get just contribution from this integral
        !
        call worker%solverdata%lhs%clear()

        !
        ! Interpolate solution to quadrature nodes: REAL
        !
        rho_r = worker%get_field('Density'   , 'value', 'face interior')
        u_r   = worker%get_field('Velocity-1', 'value', 'face interior')
        v_r   = worker%get_field('Velocity-2', 'value', 'face interior')
        w_r   = worker%get_field('Velocity-3', 'value', 'face interior')
        p_r   = worker%get_field('Pressure'  , 'value', 'face interior')

        !===========================
        !        MASS FLUX
        !===========================
        flux_1 = rho_1_rho * rho_r  + &
!                 rho_1_u   * u_r    + &
!                 rho_1_v   * v_r    + &
!                 rho_1_w   * w_r    + &
                 rho_1_p   * p_r

        flux_2 = rho_2_rho * rho_r  + &
!                 rho_2_u   * u_r    + &
!                 rho_2_v   * v_r    + &
!                 rho_2_w   * w_r    + &
                 rho_2_p   * p_r

        flux_3 = rho_3_rho * rho_r  + &
!                 rho_3_u   * u_r    + &
!                 rho_3_v   * v_r    + &
!                 rho_3_w   * w_r    + &
                 rho_3_p   * p_r

        flux_2 = ZERO
        flux_3 = ZERO
        call worker%integrate_boundary_condition('Density','Advection',flux_1,flux_2,flux_3)

        !===========================
        !     X-MOMENTUM FLUX
        !===========================
        flux_1 = u_1_rho * rho_r  + &
!                 u_1_u   * u_r    + &
!                 u_1_v   * v_r    + &
!                 u_1_w   * w_r    + &
                 u_1_p   * p_r

        flux_2 = u_2_rho * rho_r  + &
!                 u_2_u   * u_r    + &
!                 u_2_v   * v_r    + &
!                 u_2_w   * w_r    + &
                 u_2_p   * p_r

        flux_3 = u_3_rho * rho_r  + &
!                 u_3_u   * u_r    + &
!                 u_3_v   * v_r    + &
!                 u_3_w   * w_r    + &
                 u_3_p   * p_r

        flux_2 = ZERO
        flux_3 = ZERO
        call worker%integrate_boundary_condition('Velocity-1','Advection',flux_1,flux_2,flux_3)

        !============================
        !     Y-MOMENTUM FLUX
        !============================
        flux_1 = v_1_rho * rho_r  + &
!                 v_1_u   * u_r    + &
!                 v_1_v   * v_r    + &
!                 v_1_w   * w_r    + &
                 v_1_p   * p_r

        flux_2 = v_2_rho * rho_r  + &
!                 v_2_u   * u_r    + &
!                 v_2_v   * v_r    + &
!                 v_2_w   * w_r    + &
                 v_2_p   * p_r

        flux_3 = v_3_rho * rho_r  + &
!                 v_3_u   * u_r    + &
!                 v_3_v   * v_r    + &
!                 v_3_w   * w_r    + &
                 v_3_p   * p_r

        flux_2 = ZERO
        flux_3 = ZERO
        call worker%integrate_boundary_condition('Velocity-2','Advection',flux_1,flux_2,flux_3)

        !============================
        !     Z-MOMENTUM FLUX
        !============================
        flux_1 = w_1_rho * rho_r  + &
!                 w_1_u   * u_r    + &
!                 w_1_v   * v_r    + &
!                 w_1_w   * w_r    + &
                 w_1_p   * p_r

        flux_2 = w_2_rho * rho_r  + &
!                 w_2_u   * u_r    + &
!                 w_2_v   * v_r    + &
!                 w_2_w   * w_r    + &
                 w_2_p   * p_r

        flux_3 = w_3_rho * rho_r  + &
!                 w_3_u   * u_r    + &
!                 w_3_v   * v_r    + &
!                 w_3_w   * w_r    + &
                 w_3_p   * p_r

        flux_2 = ZERO
        flux_3 = ZERO
        call worker%integrate_boundary_condition('Velocity-3','Advection',flux_1,flux_2,flux_3)

        !============================
        !       ENERGY FLUX
        !============================
        flux_1 = p_1_rho * rho_r  + &
!                 p_1_u   * u_r    + &
!                 p_1_v   * v_r    + &
!                 p_1_w   * w_r    + &
                 p_1_p   * p_r

        flux_2 = p_2_rho * rho_r  + &
!                 p_2_u   * u_r    + &
!                 p_2_v   * v_r    + &
!                 p_2_w   * w_r    + &
                 p_2_p   * p_r

        flux_3 = p_3_rho * rho_r  + &
!                 p_3_u   * u_r    + &
!                 p_3_v   * v_r    + &
!                 p_3_w   * w_r    + &
                 p_3_p   * p_r

        flux_2 = ZERO
        flux_3 = ZERO
        call worker%integrate_boundary_condition('Pressure','Advection',flux_1,flux_2,flux_3)

!        !
!        ! Store off matrix contributions
!        !
!        if (worker%iface == 1) then
!            worker%solverdata%A3a = worker%solverdata%lhs
!        else if (worker%iface == 2) then
!            worker%solverdata%A3b = worker%solverdata%lhs
!        end if

    end subroutine compute
    !******************************************************************************************************

end module PRIMLINEULER_bc

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
needs util/tasks

| This will run for about three seconds, then stop
ms@ variable, start
: idle yield idle ;
:: 'a emit  ms@ start @ - 3000 >if cr ." Good bye!" cr bye then ; >task
:: 'b emit  ; >task
task: 'c emit  ;

| Start the idle loop
idle

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
include "bug-2601829-mid.h"

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
      SUBROUTINE NA1 2 A8 (*,A,N,B,NOTUSE)
C
      INTEGER         A(1),     B(2),    CDC
      CHARACTER*1     C(1),     T(8)
      CHARACTER*8     D(1)
      CHARACTER*10    BLNK,     TEMP
      COMMON /XREADX/ NOUT
      COMMON /MACHIN/ MACH
      EQUIVALENCE     (T(1),TEMP)
      DATA            BLNK / '          '  /, CDC / 4  /
C
C     THESE ROUTTNES CONVERT N A1 BCD WORDS IN A, OR N A1 CHARACTERS IN
C     C TO AN 8-BYTE BCD WORD IN B (CDC ONLY), (OR TO TWO 4-BYTE BCD
C     WORDS IN B, ALL OTHER NON-CDC MACHINES), OR AN 8-CHARACTER WORD
C     IN D, LEFT ADJUSTED.
C     CALLING ROUTINE MUST NOT USE LOGICAL*1 FOR A-ARRAY.
C     (NO SYSTEM ENCODE/DECODE FUNCTIONS ARE USED)
C
C     ENTRY POINTS   NA1 2 A8  (BCD-BYTE  VERSION)
C                    NK1 2 K8  (CHARACTER VERSION)
C
C
C     WRITTEN BY G.CHAN/SPERRY IN AUG. 1985
C     PARTICULARLY FOR XREAD ROUTINE, IN SUPPORT OF ITS NEW FREE-FIELD
C     INPUT FORMAT.  THIS SUBROUTINE IS MACHINE INDEPENDENT
C
C     LAST REVISED  8/1988
C
      IF (N .GT. 8) GO TO 40
      TEMP = BLNK
      CALL B2K (A,TEMP,N)
      IF (MACH .NE. CDC) CALL KHRBC2 (TEMP,B(1))
CWKBD IF (MACH .EQ. CDC) B(1) = ISWAP(TEMP)
      RETURN
C
      ENTRY NK1 2 K8 (*,C,N,D,NOTUSE)
C     ===============================
C
      IF (N .GT. 8) GO TO 40
      TEMP = BLNK
      DO 30 I = 1,N
 30   T(I) = C(I)
      D(1) = TEMP
      RETURN
C
 40   WRITE  (NOUT,50) N
 50   FORMAT ('   N.GT.8/NA12A8',I6)
      J = NOTUSE
      RETURN 1
      END


text: #redirect ASUCD Housing Day


text:
      SUBROUTINE SB10RD( N, M, NP, NCON, NMEAS, GAMMA, A, LDA, B, LDB,
     $                   C, LDC, D, LDD, F, LDF, H, LDH, TU, LDTU, TY,
     $                   LDTY, X, LDX, Y, LDY, AK, LDAK, BK, LDBK, CK,
     $                   LDCK, DK, LDDK, IWORK, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute the matrices of an H-infinity (sub)optimal controller
C
C              | AK | BK |
C          K = |----|----|,
C              | CK | DK |
C
C     from the state feedback matrix F and output injection matrix H as
C     determined by the SLICOT Library routine SB10QD.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the system.  N >= 0.
C
C     M       (input) INTEGER
C             The column size of the matrix B.  M >= 0.
C
C     NP      (input) INTEGER
C             The row size of the matrix C.  NP >= 0.
C
C     NCON    (input) INTEGER
C             The number of control inputs (M2).  M >= NCON >= 0.
C             NP-NMEAS >= NCON.
C
C     NMEAS   (input) INTEGER
C             The number of measurements (NP2).  NP >= NMEAS >= 0.
C             M-NCON >= NMEAS.
C
C     GAMMA   (input) DOUBLE PRECISION
C             The value of gamma. It is assumed that gamma is
C             sufficiently large so that the controller is admissible.
C             GAMMA >= 0.
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
C             The leading N-by-N part of this array must contain the
C             system state matrix A.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= max(1,N).
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
C             The leading N-by-M part of this array must contain the
C             system input matrix B.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= max(1,N).
C
C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
C             The leading NP-by-N part of this array must contain the
C             system output matrix C.
C
C     LDC     INTEGER
C             The leading dimension of the array C.  LDC >= max(1,NP).
C
C     D       (input) DOUBLE PRECISION array, dimension (LDD,M)
C             The leading NP-by-M part of this array must contain the
C             system input/output matrix D.
C
C     LDD     INTEGER
C             The leading dimension of the array D.  LDD >= max(1,NP).
C
C     F       (input) DOUBLE PRECISION array, dimension (LDF,N)
C             The leading M-by-N part of this array must contain the
C             state feedback matrix F.
C
C     LDF     INTEGER
C             The leading dimension of the array F.  LDF >= max(1,M).
C
C     H       (input) DOUBLE PRECISION array, dimension (LDH,NP)
C             The leading N-by-NP part of this array must contain the
C             output injection matrix H.
C
C     LDH     INTEGER
C             The leading dimension of the array H.  LDH >= max(1,N).
C
C     TU      (input) DOUBLE PRECISION array, dimension (LDTU,M2)
C             The leading M2-by-M2 part of this array must contain the
C             control transformation matrix TU, as obtained by the
C             SLICOT Library routine SB10PD.
C
C     LDTU    INTEGER
C             The leading dimension of the array TU.  LDTU >= max(1,M2).
C
C     TY      (input) DOUBLE PRECISION array, dimension (LDTY,NP2)
C             The leading NP2-by-NP2 part of this array must contain the
C             measurement transformation matrix TY, as obtained by the
C             SLICOT Library routine SB10PD.
C
C     LDTY    INTEGER
C             The leading dimension of the array TY.
C             LDTY >= max(1,NP2).
C
C     X       (input) DOUBLE PRECISION array, dimension (LDX,N)
C             The leading N-by-N part of this array must contain the
C             matrix X, solution of the X-Riccati equation, as obtained
C             by the SLICOT Library routine SB10QD.
C
C     LDX     INTEGER
C             The leading dimension of the array X.  LDX >= max(1,N).
C
C     Y       (input) DOUBLE PRECISION array, dimension (LDY,N)
C             The leading N-by-N part of this array must contain the
C             matrix Y, solution of the Y-Riccati equation, as obtained
C             by the SLICOT Library routine SB10QD.
C
C     LDY     INTEGER
C             The leading dimension of the array Y.  LDY >= max(1,N).
C
C     AK      (output) DOUBLE PRECISION array, dimension (LDAK,N)
C             The leading N-by-N part of this array contains the
C             controller state matrix AK.
C
C     LDAK    INTEGER
C             The leading dimension of the array AK.  LDAK >= max(1,N).
C
C     BK      (output) DOUBLE PRECISION array, dimension (LDBK,NMEAS)
C             The leading N-by-NMEAS part of this array contains the
C             controller input matrix BK.
C
C     LDBK    INTEGER
C             The leading dimension of the array BK.  LDBK >= max(1,N).
C
C     CK      (output) DOUBLE PRECISION array, dimension (LDCK,N)
C             The leading NCON-by-N part of this array contains the
C             controller output matrix CK.
C
C     LDCK    INTEGER
C             The leading dimension of the array CK.
C             LDCK >= max(1,NCON).
C
C     DK      (output) DOUBLE PRECISION array, dimension (LDDK,NMEAS)
C             The leading NCON-by-NMEAS part of this array contains the
C             controller input/output matrix DK.
C
C     LDDK    INTEGER
C             The leading dimension of the array DK.
C             LDDK >= max(1,NCON).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (LIWORK), where
C             LIWORK = max(2*(max(NP,M)-M2-NP2,M2,N),NP2)
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) contains the optimal
C             LDWORK.
C
C     LDWORK  INTEGER
C             The dimension of the array DWORK.
C             LDWORK >= max(1, M2*NP2 + NP2*NP2 + M2*M2 +
C                           max(D1*D1 + max(2*D1, (D1+D2)*NP2),
C                               D2*D2 + max(2*D2, D2*M2), 3*N,
C                               N*(2*NP2 + M2) +
C                               max(2*N*M2, M2*NP2 +
C                                           max(M2*M2+3*M2, NP2*(2*NP2+
C                                                  M2+max(NP2,N))))))
C             where D1 = NP1 - M2, D2 = M1 - NP2,
C                  NP1 = NP - NP2, M1 = M - M2.
C             For good performance, LDWORK must generally be larger.
C             Denoting Q = max(M1,M2,NP1,NP2), an upper bound is
C             max( 1, Q*(3*Q + 3*N + max(2*N, 4*Q + max(Q, N)))).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 1:  if the controller is not admissible (too small value
C                   of gamma);
C             = 2:  if the determinant of Im2 + Tu*D11HAT*Ty*D22 is zero.
C
C     METHOD
C
C     The routine implements the Glover's and Doyle's formulas [1],[2].
C
C     REFERENCES
C
C     [1] Glover, K. and Doyle, J.C.
C         State-space formulae for all stabilizing controllers that
C         satisfy an Hinf norm bound and relations to risk sensitivity.
C         Systems and Control Letters, vol. 11, pp. 167-172, 1988.
C
C     [2] Balas, G.J., Doyle, J.C., Glover, K., Packard, A., and
C         Smith, R.
C         mu-Analysis and Synthesis Toolbox.
C         The MathWorks Inc., Natick, Mass., 1995.
C
C     NUMERICAL ASPECTS
C
C     The accuracy of the result depends on the condition numbers of the
C     input and output transformations.
C
C     CONTRIBUTORS
C
C     P.Hr. Petkov, D.W. Gu and M.M. Konstantinov, October 1998.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, May 1999,
C     Sept. 1999, Oct. 2001.
C
C     KEYWORDS
C
C     Algebraic Riccati equation, H-infinity optimal control, robust
C     control.
C
C  *********************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
C     ..
C     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDAK, LDB, LDBK, LDC, LDCK, LDD,
     $                   LDDK, LDF, LDH, LDTU, LDTY, LDWORK, LDX, LDY,
     $                   M, N, NCON, NMEAS, NP
      DOUBLE PRECISION   GAMMA
C     ..
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), AK( LDAK, * ), B( LDB, * ),
     $                   BK( LDBK, * ), C( LDC, * ), CK( LDCK, * ),
     $                   D( LDD, * ), DK( LDDK, * ), DWORK( * ),
     $                   F( LDF, * ), H( LDH, * ), TU( LDTU, * ),
     $                   TY( LDTY, * ), X( LDX, * ), Y( LDY, * )
C     ..
C     .. Local Scalars ..
      INTEGER            I, ID11, ID12, ID21, IJ, INFO2, IW1, IW2, IW3,
     $                   IW4, IWB, IWC, IWRK, J, LWAMAX, M1, M2, MINWRK,
     $                   ND1, ND2, NP1, NP2
      DOUBLE PRECISION   ANORM, EPS, RCOND
C     ..
C     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
      EXTERNAL           DLAMCH, DLANGE, DLANSY
C     ..
C     .. External Subroutines ..
      EXTERNAL           DGECON, DGEMM, DGETRF, DGETRI, DGETRS, DLACPY,
     $                   DLASET, DPOTRF, DSYCON, DSYRK, DSYTRF, DSYTRS,
     $                   DTRMM, MA02AD, MB01RX, XERBLA
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC          DBLE, INT, MAX
C     ..
C     .. Executable Statements ..
C
C     Decode and Test input parameters.
C
      M1  = M - NCON
      M2  = NCON
      NP1 = NP - NMEAS
      NP2 = NMEAS
C
      INFO = 0
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( NP.LT.0 ) THEN
         INFO = -3
      ELSE IF( NCON.LT.0 .OR. M1.LT.0 .OR. M2.GT.NP1 ) THEN
         INFO = -4
      ELSE IF( NMEAS.LT.0 .OR. NP1.LT.0 .OR. NP2.GT.M1 ) THEN
         INFO = -5
      ELSE IF( GAMMA.LT.ZERO ) THEN
         INFO = -6
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -10
      ELSE IF( LDC.LT.MAX( 1, NP ) ) THEN
         INFO = -12
      ELSE IF( LDD.LT.MAX( 1, NP ) ) THEN
         INFO = -14
      ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
         INFO = -16
      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
         INFO = -18
      ELSE IF( LDTU.LT.MAX( 1, M2 ) ) THEN
         INFO = -20
      ELSE IF( LDTY.LT.MAX( 1, NP2 ) ) THEN
         INFO = -22
      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
         INFO = -24
      ELSE IF( LDY.LT.MAX( 1, N ) ) THEN
         INFO = -26
      ELSE IF( LDAK.LT.MAX( 1, N ) ) THEN
         INFO = -28
      ELSE IF( LDBK.LT.MAX( 1, N ) ) THEN
         INFO = -30
      ELSE IF( LDCK.LT.MAX( 1, M2 ) ) THEN
         INFO = -32
      ELSE IF( LDDK.LT.MAX( 1, M2 ) ) THEN
         INFO = -34
      ELSE
C
C        Compute workspace.
C
         ND1 = NP1 - M2
         ND2 = M1 - NP2
         MINWRK = MAX( 1, M2*NP2 + NP2*NP2 + M2*M2 +
     $                 MAX( ND1*ND1 + MAX( 2*ND1, ( ND1 + ND2 )*NP2 ),
     $                      ND2*ND2 + MAX( 2*ND2, ND2*M2 ), 3*N,
     $                      N*( 2*NP2 + M2 ) +
     $                         MAX( 2*N*M2, M2*NP2 +
     $                                MAX( M2*M2 + 3*M2, NP2*( 2*NP2 +
     $                                     M2 + MAX( NP2, N ) ) ) ) ) )
         IF( LDWORK.LT.MINWRK )
     $      INFO = -37
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SB10RD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( N.EQ.0 .OR. M.EQ.0 .OR. NP.EQ.0 .OR. M1.EQ.0 .OR. M2.EQ.0
     $    .OR. NP1.EQ.0 .OR. NP2.EQ.0 ) THEN
         DWORK( 1 ) = ONE
         RETURN
      END IF
C
C     Get the machine precision.
C
      EPS = DLAMCH( 'Epsilon' )
C
C     Workspace usage.
C
      ID11 = 1
      ID21 = ID11 + M2*NP2
      ID12 = ID21 + NP2*NP2
      IW1  = ID12 + M2*M2
      IW2  = IW1  + ND1*ND1
      IW3  = IW2  + ND1*NP2
      IWRK = IW2
C
C     Set D11HAT := -D1122 .
C
      IJ = ID11
      DO 20 J = 1, NP2
         DO 10 I = 1, M2
            DWORK( IJ ) = -D( ND1+I, ND2+J )
            IJ = IJ + 1
   10    CONTINUE
   20 CONTINUE
C
C     Set D21HAT := Inp2 .
C
      CALL DLASET( 'Upper', NP2, NP2, ZERO, ONE, DWORK( ID21 ), NP2 )
C
C     Set D12HAT := Im2 .
C
      CALL DLASET( 'Lower', M2, M2, ZERO, ONE, DWORK( ID12 ), M2 )
C
C     Compute D11HAT, D21HAT, D12HAT .
C
      LWAMAX = 0
      IF( ND1.GT.0 ) THEN
         IF( ND2.EQ.0 ) THEN
C
C           Compute D21HAT'*D21HAT = Inp2 - D1112'*D1112/gamma^2 .
C
            CALL DSYRK( 'U', 'T', NP2, ND1, -ONE/GAMMA**2, D, LDD, ONE,
     $                  DWORK( ID21 ), NP2 )
         ELSE
C
C           Compute gdum = gamma^2*Ind1 - D1111*D1111' .
C
            CALL DLASET( 'U', ND1, ND1, ZERO, GAMMA**2, DWORK( IW1 ),
     $                   ND1 )
            CALL DSYRK( 'U', 'N', ND1, ND2, -ONE, D, LDD, ONE,
     $                  DWORK( IW1 ), ND1 )
            ANORM = DLANSY( 'I', 'U', ND1, DWORK( IW1 ), ND1,
     $                      DWORK( IWRK ) )
            CALL DSYTRF( 'U', ND1, DWORK( IW1 ), ND1, IWORK,
     $                   DWORK( IWRK ), LDWORK-IWRK+1, INFO2 )
            IF( INFO2.GT.0 ) THEN
               INFO = 1
               RETURN
            END IF
            LWAMAX = INT( DWORK( IWRK ) ) + IWRK - 1
            CALL DSYCON( 'U', ND1, DWORK( IW1 ), ND1, IWORK, ANORM,
     $                   RCOND, DWORK( IWRK ), IWORK( ND1+1 ), INFO2 )
C
C           Return if the matrix is singular to working precision.
C
            IF( RCOND.LT.EPS ) THEN
               INFO = 1
               RETURN
            END IF
C
C           Compute inv(gdum)*D1112 .
C
            CALL DLACPY( 'Full', ND1, NP2, D( 1, ND2+1 ), LDD,
     $                   DWORK( IW2 ), ND1 )
            CALL DSYTRS( 'U', ND1, NP2, DWORK( IW1 ), ND1, IWORK,
     $                   DWORK( IW2 ), ND1, INFO2 )
C
C           Compute D11HAT = -D1121*D1111'*inv(gdum)*D1112 - D1122 .
C
            CALL DGEMM( 'T', 'N', ND2, NP2, ND1, ONE, D, LDD,
     $                  DWORK( IW2 ), ND1, ZERO, DWORK( IW3 ), ND2 )
            CALL DGEMM( 'N', 'N', M2, NP2, ND2, -ONE, D( ND1+1, 1 ),
     $                  LDD, DWORK( IW3 ), ND2, ONE, DWORK( ID11 ), M2 )
C
C           Compute D21HAT'*D21HAT = Inp2 - D1112'*inv(gdum)*D1112 .
C
            CALL MB01RX( 'Left', 'Upper', 'Transpose', NP2, ND1, ONE,
     $                   -ONE, DWORK( ID21 ), NP2, D( 1, ND2+1 ), LDD,
     $                   DWORK( IW2 ), ND1, INFO2 )
C
            IW2  = IW1 + ND2*ND2
            IWRK = IW2
C
C           Compute gdum = gamma^2*Ind2 - D1111'*D1111 .
C
            CALL DLASET( 'L', ND2, ND2, ZERO, GAMMA**2, DWORK( IW1 ),
     $                   ND2 )
            CALL DSYRK( 'L', 'T', ND2, ND1, -ONE, D, LDD, ONE,
     $                  DWORK( IW1 ), ND2 )
            ANORM = DLANSY( 'I', 'L', ND2, DWORK( IW1 ), ND2,
     $                      DWORK( IWRK ) )
            CALL DSYTRF( 'L', ND2, DWORK( IW1 ), ND2, IWORK,
     $                   DWORK( IWRK ), LDWORK-IWRK+1, INFO2 )
            IF( INFO2.GT.0 ) THEN
               INFO = 1
               RETURN
            END IF
            LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
            CALL DSYCON( 'L', ND2, DWORK( IW1 ), ND2, IWORK, ANORM,
     $                   RCOND, DWORK( IWRK ), IWORK( ND2+1 ), INFO2 )
C
C           Return if the matrix is singular to working precision.
C
            IF( RCOND.LT.EPS ) THEN
               INFO = 1
               RETURN
            END IF
C
C           Compute inv(gdum)*D1121' .
C
            CALL MA02AD( 'Full', M2, ND2, D( ND1+1, 1 ), LDD,
     $                   DWORK( IW2 ), ND2 )
            CALL DSYTRS( 'L', ND2, M2, DWORK( IW1 ), ND2, IWORK,
     $                   DWORK( IW2 ), ND2, INFO2 )
C
C           Compute D12HAT*D12HAT' = Im2 - D1121*inv(gdum)*D1121' .
C
            CALL MB01RX( 'Left', 'Lower', 'NoTranspose', M2, ND2, ONE,
     $                   -ONE, DWORK( ID12 ), M2, D( ND1+1, 1 ), LDD,
     $                   DWORK( IW2 ), ND2, INFO2 )
         END IF
      ELSE
         IF( ND2.GT.0 ) THEN
C
C           Compute D12HAT*D12HAT' = Im2 - D1121*D1121'/gamma^2 .
C
            CALL DSYRK( 'L', 'N', M2, ND2, -ONE/GAMMA**2, D, LDD, ONE,
     $                  DWORK( ID12 ), M2 )
         END IF
      END IF
C
C     Compute D21HAT using Cholesky decomposition.
C
      CALL DPOTRF( 'U', NP2, DWORK( ID21 ), NP2, INFO2 )
      IF( INFO2.GT.0 ) THEN
         INFO = 1
         RETURN
      END IF
C
C     Compute D12HAT using Cholesky decomposition.
C
      CALL DPOTRF( 'L', M2, DWORK( ID12 ), M2, INFO2 )
      IF( INFO2.GT.0 ) THEN
         INFO = 1
         RETURN
      END IF
C             _
C     Compute Z = In - Y*X/gamma^2 and its LU factorization in AK .
C
      IWRK = IW1
      CALL DLASET( 'Full', N, N, ZERO, ONE, AK, LDAK )
      CALL DGEMM( 'N', 'N', N, N, N, -ONE/GAMMA**2, Y, LDY, X, LDX,
     $            ONE, AK, LDAK )
      ANORM = DLANGE( '1', N, N, AK, LDAK, DWORK( IWRK ) )
      CALL DGETRF( N, N, AK, LDAK, IWORK, INFO2 )
      IF( INFO2.GT.0 ) THEN
         INFO = 1
         RETURN
      END IF
      CALL DGECON( '1', N, AK, LDAK, ANORM, RCOND, DWORK( IWRK ),
     $             IWORK( N+1 ), INFO )
C
C     Return if the matrix is singular to working precision.
C
      IF( RCOND.LT.EPS ) THEN
         INFO = 1
         RETURN
      END IF
C
      IWB = IW1
      IWC = IWB + N*NP2
      IW1 = IWC + ( M2 + NP2 )*N
      IW2 = IW1 + N*M2
C
C     Compute C2' + F12' in BK .
C
      DO 40 J = 1, N
         DO 30 I = 1, NP2
            BK( J, I ) = C( NP1 + I, J ) + F( ND2 + I, J )
   30    CONTINUE
   40 CONTINUE
C                                                          _
C     Compute the transpose of (C2 + F12)*Z , with Z = inv(Z) .
C
      CALL DGETRS( 'Transpose', N, NP2, AK, LDAK, IWORK, BK, LDBK,
     $             INFO2 )
C
C     Compute the transpose of F2*Z .
C
      CALL MA02AD( 'Full', M2, N, F( M1+1, 1 ), LDF, DWORK( IW1 ), N )
      CALL DGETRS( 'Transpose', N, M2, AK, LDAK, IWORK, DWORK( IW1 ), N,
     $             INFO2 )
C
C     Compute the transpose of C1HAT = F2*Z - D11HAT*(C2 + F12)*Z .
C
      CALL DGEMM( 'N', 'T', N, M2, NP2, -ONE, BK, LDBK, DWORK( ID11 ),
     $            M2, ONE, DWORK( IW1 ), N )
C
C     Compute CHAT .
C
      CALL DGEMM( 'N', 'T', M2, N, M2, ONE, TU, LDTU, DWORK( IW1 ), N,
     $            ZERO, DWORK( IWC ), M2+NP2 )
      CALL MA02AD( 'Full', N, NP2, BK, LDBK, DWORK( IWC+M2 ), M2+NP2 )
      CALL DTRMM( 'L', 'U', 'N', 'N', NP2, N, -ONE, DWORK( ID21 ), NP2,
     $            DWORK( IWC+M2 ), M2+NP2 )
C
C     Compute B2 + H12 .
C
      IJ = IW2
      DO 60 J = 1, M2
         DO 50 I = 1, N
            DWORK( IJ ) = B( I, M1 + J ) + H( I, ND1 + J )
            IJ = IJ + 1
   50    CONTINUE
   60 CONTINUE
C
C     Compute A + HC in AK .
C
      CALL DLACPY( 'Full', N, N, A, LDA, AK, LDAK )
      CALL DGEMM( 'N', 'N', N, N, NP, ONE, H, LDH, C, LDC, ONE, AK,
     $            LDAK )
C
C     Compute AHAT = A + HC + (B2 + H12)*C1HAT in AK .
C
      CALL DGEMM( 'N', 'T', N, N, M2, ONE, DWORK( IW2 ), N,
     $            DWORK( IW1 ), N, ONE, AK, LDAK )
C
C     Compute B1HAT = -H2 + (B2 + H12)*D11HAT in BK .
C
      CALL DLACPY( 'Full', N, NP2, H( 1, NP1+1 ), LDH, BK, LDBK )
      CALL DGEMM( 'N', 'N', N, NP2, M2, ONE, DWORK( IW2 ), N,
     $            DWORK( ID11 ), M2, -ONE, BK, LDBK )
C
C     Compute the first block of BHAT, BHAT1 .
C
      CALL DGEMM( 'N', 'N', N, NP2, NP2, ONE, BK, LDBK, TY, LDTY, ZERO,
     $            DWORK( IWB ), N )
C
C     Compute Tu*D11HAT .
C
      CALL DGEMM( 'N', 'N', M2, NP2, M2, ONE, TU, LDTU, DWORK( ID11 ),
     $            M2, ZERO, DWORK( IW1 ), M2 )
C
C     Compute Tu*D11HAT*Ty in DK .
C
      CALL DGEMM( 'N', 'N', M2, NP2, NP2, ONE, DWORK( IW1 ), M2, TY,
     $            LDTY, ZERO, DK, LDDK )
C
C     Compute P = Im2 + Tu*D11HAT*Ty*D22 and its condition.
C
      IW2  = IW1 + M2*NP2
      IWRK = IW2 + M2*M2
      CALL DLASET( 'Full', M2, M2, ZERO, ONE, DWORK( IW2 ), M2 )
      CALL DGEMM( 'N', 'N', M2, M2, NP2, ONE, DK, LDDK,
     $            D( NP1+1, M1+1 ), LDD, ONE, DWORK( IW2 ), M2 )
      ANORM = DLANGE( '1', M2, M2, DWORK( IW2 ), M2, DWORK( IWRK ) )
      CALL DGETRF( M2, M2, DWORK( IW2 ), M2, IWORK, INFO2 )
      IF( INFO2.GT.0 ) THEN
         INFO = 2
         RETURN
      END IF
      CALL DGECON( '1', M2, DWORK( IW2 ), M2, ANORM, RCOND,
     $             DWORK( IWRK ), IWORK( M2+1 ), INFO2 )
C
C     Return if the matrix is singular to working precision.
C
      IF( RCOND.LT.EPS ) THEN
         INFO = 2
         RETURN
      END IF
C
C     Find the controller matrix CK, CK = inv(P)*CHAT(1:M2,:) .
C
      CALL DLACPY( 'Full', M2, N, DWORK( IWC ), M2+NP2, CK, LDCK )
      CALL DGETRS( 'NoTranspose', M2, N, DWORK( IW2 ), M2, IWORK, CK,
     $             LDCK, INFO2 )
C
C     Find the controller matrices AK, BK, and DK, exploiting the
C     special structure of the relations.
C
C     Compute Q = Inp2 + D22*Tu*D11HAT*Ty and its LU factorization.
C
      IW3  = IW2 + NP2*NP2
      IW4  = IW3 + NP2*M2
      IWRK = IW4 + NP2*NP2
      CALL DLASET( 'Full', NP2, NP2, ZERO, ONE, DWORK( IW2 ), NP2 )
      CALL DGEMM( 'N', 'N', NP2, NP2, M2, ONE, D( NP1+1, M1+1 ), LDD,
     $            DK, LDDK, ONE, DWORK( IW2 ), NP2 )
      CALL DGETRF( NP2, NP2, DWORK( IW2 ), NP2, IWORK, INFO2 )
      IF( INFO2.GT.0 ) THEN
         INFO = 2
         RETURN
      END IF
C
C     Compute A1 = inv(Q)*D22 and inv(Q) .
C
      CALL DLACPY( 'Full', NP2, M2, D( NP1+1, M1+1 ), LDD, DWORK( IW3 ),
     $             NP2 )
      CALL DGETRS( 'NoTranspose', NP2, M2, DWORK( IW2 ), NP2, IWORK,
     $             DWORK( IW3 ), NP2, INFO2 )
      CALL DGETRI( NP2, DWORK( IW2 ), NP2, IWORK, DWORK( IWRK ),
     $             LDWORK-IWRK+1, INFO2 )
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
C
C     Compute A2 = ( inv(Ty) - inv(Q)*inv(Ty) -
C                    A1*Tu*D11HAT )*inv(D21HAT) .
C
      CALL DLACPY( 'Full', NP2, NP2, TY, LDTY, DWORK( IW4 ), NP2 )
      CALL DGETRF( NP2, NP2, DWORK( IW4 ), NP2, IWORK, INFO2 )
      CALL DGETRI( NP2, DWORK( IW4 ), NP2, IWORK, DWORK( IWRK ),
     $             LDWORK-IWRK+1, INFO2 )
C
      CALL DLACPY( 'Full', NP2, NP2, DWORK( IW4 ), NP2, DWORK( IWRK ),
     $             NP2 )
      CALL DGEMM( 'N', 'N', NP2, NP2, NP2, -ONE, DWORK( IW2), NP2,
     $            DWORK( IWRK ), NP2, ONE, DWORK( IW4 ), NP2 )
      CALL DGEMM( 'N', 'N', NP2, NP2, M2, -ONE, DWORK( IW3), NP2,
     $            DWORK( IW1 ), M2, ONE, DWORK( IW4 ), NP2 )
      CALL DTRMM( 'R', 'U', 'N', 'N', NP2, NP2, ONE, DWORK( ID21 ), NP2,
     $            DWORK( IW4 ), NP2 )
C
C     Compute [ A1  A2 ]*CHAT .
C
      CALL DGEMM( 'N', 'N', NP2, N, M2+NP2, ONE, DWORK( IW3 ), NP2,
     $            DWORK( IWC ), M2+NP2, ZERO, DWORK( IWRK ), NP2 )
C
C     Compute AK := AHAT - BHAT1*[ A1  A2 ]*CHAT .
C
      CALL DGEMM( 'N', 'N', N, N, NP2, -ONE, DWORK( IWB ), N,
     $            DWORK( IWRK ), NP2, ONE, AK, LDAK )
C
C     Compute BK := BHAT1*inv(Q) .
C
      CALL DGEMM( 'N', 'N', N, NP2, NP2, ONE, DWORK( IWB ), N,
     $            DWORK( IW2 ), NP2, ZERO, BK, LDBK )
C
C     Compute DK := Tu*D11HAT*Ty*inv(Q) .
C
      CALL DGEMM( 'N', 'N', M2, NP2, NP2, ONE, DK, LDDK, DWORK( IW2 ),
     $            NP2, ZERO, DWORK( IW3 ), M2 )
      CALL DLACPY( 'Full', M2, NP2, DWORK( IW3 ), M2, DK, LDDK )
C
      DWORK( 1 ) = DBLE( LWAMAX )
      RETURN
C *** Last line of SB10RD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      subroutine qqb_w2jet(p,msq)
      implicit none
c--- matrix element squared and averaged over initial colours and spins
c     q(-p1) + qbar(-p2) --> W + f(p5) + f(p6)
c                            |
c                            --> nu(p3) + e^+(p4)
c     where the fermions are either q(p5) and qbar(p6) [Qflag = .true.]
c                                or g(p5) and g(p6)    [Gflag = .true.]
c--- all momenta are incoming
      include 'constants.f'
      include 'masses.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'ckm.f'
      include 'sprods_com.f'
      include 'zprods_com.f'
      include 'flags.f'
      include 'msq_cs.f'
      include 'lc.f'
      include 'part.f'
      include 'first.f'
      double precision msq(-nf:nf,-nf:nf),p(mxpart,4),
     .                 facgg,facqq,prop,Vfac
      double precision qqbWgg2,qbqWgg2,qgWqg2,qbgWqbg2,
     .                 gqbWqbg2,gqWqg2,ggWqbq2
      double precision qqbWgg2_cs(0:2),qbqWgg2_cs(0:2),qgWqg2_cs(0:2),
     .                 qbgWqbg2_cs(0:2),gqbWqbg2_cs(0:2),
     .                 gqWqg2_cs(0:2),ggWqbq2_cs(0:2)
      double precision
     . qqb_ijkk(0:2),qqb_ijii(0:2),qqb_ijjj(0:2),qqb_ijkj(0:2),
     . qqb_ijik(0:2),qqb_ijkl(0:2),qqb_iiij(0:2),qqb_iiji(0:2),
     . qbq_ijkk(0:2),qbq_ijii(0:2),qbq_ijjj(0:2),qbq_ijkj(0:2),
     . qbq_ijik(0:2),qbq_ijkl(0:2),qbq_iiij(0:2),qbq_iiji(0:2),
     . qq_iiji(0:2),qq_ijkj(0:2),qq_ijik(0:2),
     . qq_ijjj(0:2),qq_ijii(0:2),
     . qbqb_iiji(0:2),qbqb_ijkj(0:2),qbqb_ijik(0:2),
     . qbqb_ijjj(0:2),qbqb_ijii(0:2)
      double precision mqq(0:2,fn:nf,fn:nf)
      double complex qqb1(3),qqb2(3),qqb3(3),qqb4(3),
     .               qq1(4),qq2(4),qq3(4),qq4(4),
     .               qbq1(3),qbq2(3),qbq3(3),qbq4(3),
     .               qbqb1(4),qbqb2(4),qbqb3(4),qbqb4(4)
      integer rcolourchoice
      common/mqq/mqq
c--- we label the amplitudes by helicity (qqb1 ... qqb4)
c--- and by type of contribution qqb(1) ... qqb(n)
      integer i,j,k,n1,n2
!$omp threadprivate(/mqq/)

      if (first) then
      first=.false.
        if (Gflag) then
          write(*,*) 'Using QQGG matrix elements'
          write(*,*) '[LC is   N ]'
          write(*,*) '[SLC is 1/N]'
        endif
        if (Qflag) then
          write(*,*) 'Using QQBQQB matrix elements'
          write(*,*) '[LC is   1 ]'
          write(*,*) '[SLC is 1/N]'
        endif
        if (part .eq. 'lord') then
          if     (colourchoice .eq. 1) then
            write(*,*) 'Leading colour only'
          elseif (colourchoice .eq. 2) then
            write(*,*) 'Sub-leading colour only'
          elseif (colourchoice .eq. 0) then
            write(*,*) 'Total of both colour structures'
          else
            write(*,*) 'Bad colourchoice'
            stop
          endif
        else
          write(*,*) 'Calculating all colour structures in LO'
        endif
      endif

c--- if we're calculating the REAL or VIRT matrix elements, we
c--- need all the colour structures, but want to preserve
c--- the actual value of colourchoice
      if ((part .eq. 'real') .or. (part .eq. 'virt')) then
        rcolourchoice=colourchoice
        colourchoice=0
      endif

c--- initialize matrix elements
      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      enddo
      enddo

c--- set up spinors
      call spinoru(6,p,za,zb)
      prop=s(3,4)**2/((s(3,4)-wmass**2)**2+wmass**2*wwidth**2)
      facqq=4d0*V*gsq**2*(gwsq/2d0)**2*aveqq*prop
      facgg=V*xn/four*(gwsq/2d0)**2*gsq**2*prop

c--- calculate 2-quark, 2-gluon amplitudes
      if (Gflag) then
        call w2jetsq(1,2,3,4,5,6,za,zb,qqbWgg2)
        call storecs(qqbWgg2_cs)
        call w2jetsq(2,1,3,4,5,6,za,zb,qbqWgg2)
        call storecs(qbqWgg2_cs)
        call w2jetsq(1,5,3,4,2,6,za,zb,qgWqg2)
        call storecs(qgWqg2_cs)
        call w2jetsq(2,5,3,4,1,6,za,zb,gqWqg2)
        call storecs(gqWqg2_cs)
        call w2jetsq(5,1,3,4,2,6,za,zb,qbgWqbg2)
        call storecs(qbgWqbg2_cs)
        call w2jetsq(5,2,3,4,1,6,za,zb,gqbWqbg2)
        call storecs(gqbWqbg2_cs)
        call w2jetsq(5,6,3,4,1,2,za,zb,ggWqbq2)
        call storecs(ggWqbq2_cs)
        do i=0,2
          qqbWgg2_cs(i) = half*aveqq*facgg*qqbWgg2_cs(i)
          qbqWgg2_cs(i) = half*aveqq*facgg*qbqWgg2_cs(i)
          gqWqg2_cs(i)  = aveqg*facgg*gqWqg2_cs(i)
          qgWqg2_cs(i)  = aveqg*facgg*qgWqg2_cs(i)
          gqbWqbg2_cs(i)= aveqg*facgg*gqbWqbg2_cs(i)
          qbgWqbg2_cs(i)= aveqg*facgg*qbgWqbg2_cs(i)
          ggWqbq2_cs(i) = avegg*facgg*ggWqbq2_cs(i)
        enddo
        qqbWgg2 = qqbWgg2_cs(1) +qqbWgg2_cs(2) +qqbWgg2_cs(0)
        qbqWgg2 = qbqWgg2_cs(1) +qbqWgg2_cs(2) +qbqWgg2_cs(0)
        gqWqg2  = gqWqg2_cs(1)  +gqWqg2_cs(2)  +gqWqg2_cs(0)
        qgWqg2  = qgWqg2_cs(1)  +qgWqg2_cs(2)  +qgWqg2_cs(0)
        gqbWqbg2= gqbWqbg2_cs(1)+gqbWqbg2_cs(2)+gqbWqbg2_cs(0)
        qbgWqbg2= qbgWqbg2_cs(1)+qbgWqbg2_cs(2)+qbgWqbg2_cs(0)
        ggWqbq2 = ggWqbq2_cs(1) +ggWqbq2_cs(2) +ggWqbq2_cs(0)
      endif

c--- calculate four-quark amplitudes
      if (Qflag) then
c--- basic amplitudes - q qb --> W + g* (--> q qb) (amps 1 and 3)
c---                and q qb --> g* --> q (--> W q) qb (amps 2 and 4)
        call amp_q_QbQ_qb(1,2,5,6,qqb1(1),qqb2(1),qqb3(1),qqb4(1))
c--- crossed - q qb --> q qb ( --> W qb)
        call amp_q_QbQ_qb(1,5,2,6,qqb1(2),qqb2(2),qqb3(2),qqb4(2))
c--- crossed - q qb --> q ( --> W q) qb
        call amp_q_QbQ_qb(6,2,5,1,qqb1(3),qqb2(3),qqb3(3),qqb4(3))

c--- now the qb q amplitudes
        call amp_q_QbQ_qb(2,1,5,6,qbq1(1),qbq2(1),qbq3(1),qbq4(1))
c--- crossed - qb q --> qb q ( --> W q)
        call amp_q_QbQ_qb(2,5,1,6,qbq1(2),qbq2(2),qbq3(2),qbq4(2))
c--- crossed - qb q --> qb ( --> W qb) q
        call amp_q_QbQ_qb(6,1,5,2,qbq1(3),qbq2(3),qbq3(3),qbq4(3))

c--- crossed q q --> q ( --> W q) q
        call amp_q_QbQ_qb(1,5,6,2,qq1(1),qq2(1),qq3(1),qq4(1))
c--- crossed q q --> q q ( --> W q)
        call amp_q_QbQ_qb(2,5,6,1,qq1(2),qq2(2),qq3(2),qq4(2))
c--- crossed q q --> q q ( --> W q)
        call amp_q_QbQ_qb(2,6,5,1,qq1(3),qq2(3),qq3(3),qq4(3))
c--- crossed q q --> q ( --> W q) q
        call amp_q_QbQ_qb(1,6,5,2,qq1(4),qq2(4),qq3(4),qq4(4))

c--- crossed qb qb --> qb ( --> W qb) qb
        call amp_q_QbQ_qb(5,1,2,6,qbqb1(1),qbqb2(1),qbqb3(1),qbqb4(1))
c--- crossed qb qb --> qb qb ( --> W qb)
        call amp_q_QbQ_qb(5,2,1,6,qbqb1(2),qbqb2(2),qbqb3(2),qbqb4(2))
c--- crossed qb qb --> qb ( --> W qb) qb
        call amp_q_QbQ_qb(6,2,1,5,qbqb1(3),qbqb2(3),qbqb3(3),qbqb4(3))
c--- crossed qb qb --> qb qb ( --> W qb)
        call amp_q_QbQ_qb(6,1,2,5,qbqb1(4),qbqb2(4),qbqb3(4),qbqb4(4))

c--- now square these amplitudes separating into color structures
c   1) Amplitude
c   2) Amplitude with (5<-->6)
c   0) Interference between above
c
c--- q(i) qb(j) --> W + g* (--> q(k) qb(k)) with k != i,j
        qqb_ijkk(1)=abs(qqb1(1))**2+abs(qqb3(1))**2
        qqb_ijkk(2)=zip
        qqb_ijkk(0)=zip

        qbq_ijkk(1)=abs(qbq1(1))**2+abs(qbq3(1))**2
        qbq_ijkk(2)=zip
        qbq_ijkk(0)=zip
c--- q(i) qb(j) --> W + g* (--> q(i) qb(i)) i.e. k = i
        qqb_ijii(1)=abs(qqb1(1))**2+abs(qqb3(1))**2
        qqb_ijii(2)=abs(qqb2(2))**2+abs(qqb4(2))**2
        qqb_ijii(0)=+2d0/xn*dble(qqb1(1)*Dconjg(qqb2(2)))

        qbq_ijii(1)=abs(qbq1(1))**2+abs(qbq3(1))**2
        qbq_ijii(2)=abs(qbq2(2))**2+abs(qbq4(2))**2
        qbq_ijii(0)=+2d0/xn*dble(qbq1(1)*Dconjg(qbq2(2)))
c--- q(i) qb(j) --> W + g* (--> q(j) qb(j)) i.e. k = j
        qqb_ijjj(1)=abs(qqb1(1))**2+abs(qqb3(1))**2
        qqb_ijjj(2)=abs(qqb2(3))**2+abs(qqb4(3))**2
        qqb_ijjj(0)=2d0/xn*dble(qqb1(1)*Dconjg(qqb2(3)))

        qbq_ijjj(1)=abs(qbq1(1))**2+abs(qbq3(1))**2
        qbq_ijjj(2)=abs(qbq2(3))**2+abs(qbq4(3))**2
        qbq_ijjj(0)=2d0/xn*dble(qbq1(1)*Dconjg(qbq2(3)))
c--- q (i) qb(j) --> q(i) qb(j) ( --> W qb(k)) with k != i,j
        qqb_ijik(1)=zip
        qqb_ijik(2)=abs(qqb2(2))**2+abs(qqb4(2))**2
        qqb_ijik(0)=zip

        qbq_ijik(2)=abs(qbq2(2))**2+abs(qbq4(2))**2
        qbq_ijik(1)=zip
        qbq_ijik(0)=zip
c--- q (i) qb(j) --> q(i) ( --> W q(k)) qb(j) with k != i,j
        qqb_ijkj(1)=zip
        qqb_ijkj(2)=abs(qqb2(3))**2+abs(qqb4(3))**2
        qqb_ijkj(0)=zip

        qbq_ijkj(2)=abs(qbq2(3))**2+abs(qbq4(3))**2
        qbq_ijkj(1)=zip
        qbq_ijkj(0)=zip
c--- q(i) qb(j) --> g* --> q(l) (--> W q(k)) qb(l)
        qqb_ijkl(1)=abs(qqb2(1))**2+abs(qqb4(1))**2
        qqb_ijkl(2)=zip
        qqb_ijkl(0)=zip

        qbq_ijkl(1)=abs(qbq2(1))**2+abs(qbq4(1))**2
        qbq_ijkl(2)=zip
        qbq_ijkl(0)=zip
c--- q(i) qb(i) --> g* --> q(j) (--> W q(i)) qb(j)
        qqb_iiij(1)=abs(qqb2(1))**2+abs(qqb4(1))**2
        qqb_iiij(2)=abs(qqb1(3))**2+abs(qqb3(3))**2
        qqb_iiij(0)=2d0/xn*dble(qqb2(1)*Dconjg(qqb1(3)))

        qbq_iiij(1)=abs(qbq2(1))**2+abs(qbq4(1))**2
        qbq_iiij(2)=abs(qbq1(3))**2+abs(qbq3(3))**2
        qbq_iiij(0)=2d0/xn*dble(qbq2(1)*Dconjg(qbq1(3)))
c--- q(i) qb(i) --> g* --> q(j) qb(j) (--> W qb(i))
        qqb_iiji(1)=abs(qqb2(1))**2+abs(qqb4(1))**2
        qqb_iiji(2)=abs(qqb1(2))**2+abs(qqb3(2))**2
        qqb_iiji(0)=2d0/xn*dble(qqb2(1)*Dconjg(qqb1(2)))

        qbq_iiji(1)=abs(qbq2(1))**2+abs(qbq4(1))**2
        qbq_iiji(2)=abs(qbq1(2))**2+abs(qbq3(2))**2
        qbq_iiji(0)=2d0/xn*dble(qbq2(1)*Dconjg(qbq1(2)))

c--- q(i) q(i) --> q(i) ( --> W q(j) ) q(i)
        qq_iiji(1)=abs(qq1(1))**2+abs(qq3(1))**2
        qq_iiji(2)=abs(qq1(2))**2+abs(qq3(2))**2
        qq_iiji(0)=2d0/xn*dble(qq1(1)*Dconjg(qq1(2)))
c--- q(i) q(j) --> q(i) ( --> W q(k) ) q(j)
        qq_ijkj(1)=abs(qq1(1))**2+abs(qq3(1))**2
        qq_ijkj(2)=zip
        qq_ijkj(0)=zip
c--- q(i) q(j) --> q(i) q(j) ( --> W q(k) )
        qq_ijik(1)=abs(qq1(3))**2+abs(qq3(3))**2
        qq_ijik(2)=zip
        qq_ijik(0)=zip
c--- q(i) q(j) --> q(i) ( --> W q(j) ) q(j)
        qq_ijjj(1)=abs(qq1(1))**2+abs(qq3(1))**2
        qq_ijjj(2)=abs(qq1(4))**2+abs(qq3(4))**2
        qq_ijjj(0)=2d0/xn*dble(qq1(1)*Dconjg(qq1(4)))
c--- q(i) q(j) --> q(i) q(j) ( --> W q(i) )
        qq_ijii(1)=abs(qq1(3))**2+abs(qq3(3))**2
        qq_ijii(2)=abs(qq1(2))**2+abs(qq3(2))**2
        qq_ijii(0)=2d0/xn*dble(qq1(3)*Dconjg(qq1(2)))

c--- qb(i) qb(i) --> qb(i) ( --> W qb(j) ) qb(i)
        qbqb_iiji(1)=abs(qbqb1(1))**2+abs(qbqb3(1))**2
        qbqb_iiji(2)=abs(qbqb1(2))**2+abs(qbqb3(2))**2
        qbqb_iiji(0)=2d0/xn*dble(qbqb1(1)*Dconjg(qbqb1(2)))
c--- qb(i) qb(j) --> qb(i) ( --> W qb(k) ) qb(j)
        qbqb_ijkj(1)=abs(qbqb1(1))**2+abs(qbqb3(1))**2
        qbqb_ijkj(2)=zip
        qbqb_ijkj(0)=zip
c--- qb(i) qb(j) --> qb(i) qb(j) ( --> W qb(k) )
        qbqb_ijik(1)=abs(qbqb1(3))**2+abs(qbqb3(3))**2
        qbqb_ijik(2)=zip
        qbqb_ijik(0)=zip
c--- qb(i) qb(j) --> qb(i) ( --> W qb(j) ) qb(j)
        qbqb_ijjj(1)=abs(qbqb1(1))**2+abs(qbqb3(1))**2
        qbqb_ijjj(2)=abs(qbqb1(4))**2+abs(qbqb3(4))**2
        qbqb_ijjj(0)=2d0/xn*dble(qbqb1(1)*Dconjg(qbqb1(4)))
c--- qb(i) qb(j) --> qb(i) qb(j) ( --> W qb(i) )
        qbqb_ijii(2)=abs(qbqb1(2))**2+abs(qbqb3(2))**2
        qbqb_ijii(1)=abs(qbqb1(3))**2+abs(qbqb3(3))**2
        qbqb_ijii(0)=2d0/xn*dble(qbqb1(2)*Dconjg(qbqb1(3)))
      endif

c--- 4-quark contribution to matrix elements
      if (Qflag) then

      do j=-nf,nf
      do k=-nf,nf
      mqq(0,j,k)=0d0
      mqq(1,j,k)=0d0
      mqq(2,j,k)=0d0
      if ((j .gt. 0) .and. (k .lt. 0)) then
        if (j .ne. -k) then
c--- Q QBAR - different flavours
           mqq(0,j,k)=facqq*Vsq(j,k)*(qqb_ijii(0)+qqb_ijjj(0))
           mqq(1,j,k)=facqq*Vsq(j,k)*(
     .             (nf-2)*qqb_ijkk(1)+qqb_ijii(1)+qqb_ijjj(1))
           mqq(2,j,k)=facqq*Vsq(j,k)*(qqb_ijii(2)+qqb_ijjj(2))
     .            +facqq*(Vsum(j)-Vsq(j,k))*qqb_ijkj(2)
     .            +facqq*(Vsum(k)-Vsq(j,k))*qqb_ijik(2)
        else
c--- Q QBAR - same flavours
          Vfac=0d0
          do n1=1,nf
          do n2=-nf,-1
          if ((n1 .ne. j) .and. (n2 .ne. k)) then
            Vfac=Vfac+Vsq(n1,n2)
          endif
          enddo
          enddo
          mqq(0,j,k)=
     .             +facqq*Vsum(k)*qqb_iiij(0)
     .             +facqq*Vsum(j)*qqb_iiji(0)
          mqq(1,j,k)=
     .             +facqq*Vsum(k)*qqb_iiij(1)
     .             +facqq*Vsum(j)*qqb_iiji(1)
     .             +facqq*Vfac*qqb_ijkl(1)
          mqq(2,j,k)=
     .             +facqq*Vsum(k)*qqb_iiij(2)
     .             +facqq*Vsum(j)*qqb_iiji(2)
        endif
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
        if (j .ne. -k) then
c--- QBAR Q - different flavours
          mqq(0,j,k)=facqq*Vsq(j,k)*(
     .             +qbq_ijii(0)+qbq_ijjj(0))
          mqq(1,j,k)=facqq*Vsq(j,k)*(
     .             (nf-2)*qbq_ijkk(1)
     .             +qbq_ijii(1)+qbq_ijjj(1))
          mqq(2,j,k)=facqq*Vsq(j,k)*(qbq_ijii(2)+qbq_ijjj(2))
     .             +facqq*(Vsum(k)-Vsq(j,k))*qbq_ijkj(2)
     .             +facqq*(Vsum(j)-Vsq(j,k))*qbq_ijik(2)
        else
c--- QBAR Q - same flavours
          Vfac=0d0
          do n1=-nf,-1
          do n2=1,nf
          if ((n1 .ne. j) .and. (n2 .ne. k)) then
            Vfac=Vfac+Vsq(n1,n2)
          endif
          enddo
          enddo
          mqq(0,j,k)=
     .             +facqq*Vsum(j)*qbq_iiij(0)
     .             +facqq*Vsum(k)*qbq_iiji(0)
          mqq(1,j,k)=
     .             +facqq*Vsum(j)*qbq_iiij(1)
     .             +facqq*Vsum(k)*qbq_iiji(1)
     .             +facqq*Vfac*qbq_ijkl(1)
          mqq(2,j,k)=
     .             +facqq*Vsum(j)*qbq_iiij(2)
     .             +facqq*Vsum(k)*qbq_iiji(2)
        endif
      elseif ((j .gt. 0) .and. (k .gt. 0)) then
        if (j .ne. k) then
c--- Q Q - different flavours
          mqq(0,j,k)=
     .          +facqq*half*Vsq(j,-k)*qq_ijjj(0)
     .          +facqq*half*Vsq(k,-j)*qq_ijii(0)
          mqq(1,j,k)=
     .          +facqq*(Vsum(j)-Vsq(j,-k))*qq_ijkj(1)
     .          +facqq*(Vsum(k)-Vsq(k,-j))*qq_ijik(1)
     .          +facqq*half*Vsq(j,-k)*qq_ijjj(1)
     .          +facqq*half*Vsq(k,-j)*qq_ijii(1)
          mqq(2,j,k)=
     .          +facqq*half*Vsq(j,-k)*qq_ijjj(2)
     .          +facqq*half*Vsq(k,-j)*qq_ijii(2)
        else
c--- Q Q - same flavours
          mqq(0,j,k)=facqq*Vsum(j)*qq_iiji(0)
          mqq(1,j,k)=facqq*Vsum(j)*qq_iiji(1)
          mqq(2,j,k)=facqq*Vsum(j)*qq_iiji(2)
        endif
      elseif ((j .lt. 0) .and. (k .lt. 0)) then
        if (j .ne. k) then
c--- QBAR QBAR - different flavours
          mqq(0,j,k)=
     .+facqq*half*Vsq(j,-k)*qbqb_ijjj(0)
     .+facqq*half*Vsq(k,-j)*qbqb_ijii(0)
          mqq(1,j,k)=facqq*(Vsum(j)-Vsq(j,-k))*qbqb_ijkj(1)
     .              +facqq*(Vsum(k)-Vsq(k,-j))*qbqb_ijik(1)
     .+facqq*half*Vsq(j,-k)*qbqb_ijjj(1)
     .+facqq*half*Vsq(k,-j)*qbqb_ijii(1)
          mqq(2,j,k)=
     .+facqq*half*Vsq(j,-k)*qbqb_ijjj(2)
     .+facqq*half*Vsq(k,-j)*qbqb_ijii(2)
        else
c--- QBAR QBAR - same flavours
          mqq(0,j,k)=facqq*Vsum(j)*qbqb_iiji(0)
          mqq(1,j,k)=facqq*Vsum(j)*qbqb_iiji(1)
          mqq(2,j,k)=facqq*Vsum(j)*qbqb_iiji(2)
        endif
      endif
      if     (colourchoice .eq. 1) then
        mqq(0,j,k)=0d0
      elseif (colourchoice .eq. 2) then
        mqq(1,j,k)=0d0
        mqq(2,j,k)=0d0
      endif
      msq(j,k)=mqq(0,j,k)+mqq(1,j,k)+mqq(2,j,k)
      enddo
      enddo

      endif

c--- 2-quark, 2-gluon contribution to matrix elements
      if (Gflag) then

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

      do i=0,2
      msq_cs(i,j,k)=0d0
      enddo

      if     ((j .gt. 0) .and. (k .lt. 0)) then
          msq(j,k)=msq(j,k)+Vsq(j,k)*qqbWgg2
          do i=0,2
            msq_cs(i,j,k)=Vsq(j,k)*qqbWgg2_cs(i)
          enddo
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
          msq(j,k)=msq(j,k)+Vsq(j,k)*qbqWgg2
          do i=0,2
            msq_cs(i,j,k)=Vsq(j,k)*qbqWgg2_cs(i)
          enddo
      elseif ((j .gt. 0) .and. (k .eq. 0)) then
          msq(j,k)=msq(j,k)+
     &(Vsq(j,-1)+Vsq(j,-2)+Vsq(j,-3)+Vsq(j,-4)+Vsq(j,-5))*qgWqg2
          do i=0,2
            msq_cs(i,j,k)=
     &(Vsq(j,-1)+Vsq(j,-2)+Vsq(j,-3)+Vsq(j,-4)+Vsq(j,-5))*qgWqg2_cs(i)
          enddo
      elseif ((j .lt. 0) .and. (k .eq. 0)) then
          msq(j,k)=msq(j,k)+
     &(Vsq(j,+1)+Vsq(j,+2)+Vsq(j,+3)+Vsq(j,+4)+Vsq(j,+5))*qbgWqbg2
          do i=0,2
            msq_cs(i,j,k)=
     &(Vsq(j,+1)+Vsq(j,+2)+Vsq(j,+3)+Vsq(j,+4)+Vsq(j,+5))*qbgWqbg2_cs(i)
          enddo
      elseif ((j .eq. 0) .and. (k .gt. 0)) then
          msq(j,k)=msq(j,k)+
     &(Vsq(-1,k)+Vsq(-2,k)+Vsq(-3,k)+Vsq(-4,k)+Vsq(-5,k))*gqWqg2
          do i=0,2
            msq_cs(i,j,k)=
     &(Vsq(-1,k)+Vsq(-2,k)+Vsq(-3,k)+Vsq(-4,k)+Vsq(-5,k))*gqWqg2_cs(i)
          enddo
      elseif ((j .eq. 0) .and. (k .lt. 0)) then
          msq(j,k)=msq(j,k)+
     &(Vsq(+1,k)+Vsq(+2,k)+Vsq(+3,k)+Vsq(+4,k)+Vsq(+5,k))*gqbWqbg2
          do i=0,2
            msq_cs(i,j,k)=
     &(Vsq(+1,k)+Vsq(+2,k)+Vsq(+3,k)+Vsq(+4,k)+Vsq(+5,k))*gqbWqbg2_cs(i)
          enddo
      elseif ((j .eq. 0) .and. (k .eq. 0)) then
          Vfac=0d0
          do n1=1,nf
            do n2=-nf,-1
              Vfac=Vfac+Vsq(n1,n2)
            enddo
          enddo
          msq(j,k)=msq(j,k)+Vfac*ggWqbq2
          do i=0,2
            msq_cs(i,j,k)=Vfac*ggWqbq2_cs(i)
          enddo
      endif

      enddo
      enddo

      endif

c--- restore proper colourchoice if necessary
      if ((part .eq. 'real') .or. (part .eq. 'virt')) then
        colourchoice=rcolourchoice
      endif

      return
      end

      subroutine amp_q_QbQ_qb(i1,i2,i5,i6,amp1,amp2,amp3,amp4)
      implicit none
      integer i1,i2,i5,i6
      double complex aqqb_zbb_new,amp1,amp2,amp3,amp4
c--- Amplitudes for q(i1) + qb(i2) --> qb(i6) + q(i5) + W (-> 3+4)
c
c--- the form of this function is taken from the subroutine msqzbb
c--- in qqb_zbb.f. See there for comments.
c--- We return four amplitudes, in pairs corresponding to diagrams
c--- where the W couples to both quark lines

c--- quark i5 is left-handed
      amp1=+aqqb_zbb_new(i1,i6,i5,i2,3,4)
      amp2=-Dconjg(aqqb_zbb_new(i5,i2,i1,i6,4,3))

c--- quark i5 is right-handed
      amp3=-aqqb_zbb_new(i1,i5,i6,i2,3,4)
      amp4=-Dconjg(aqqb_zbb_new(i5,i1,i2,i6,4,3))

      return
      end

      subroutine storecs(mcs)
c-- this routine transfers the information on the colour structure
c-- for the W2jet matrix elements into separate arrays for each
c-- incoming parton case
      implicit none
      include 'mmsq_cs.f'
      integer i
      double precision mcs(0:2)

      do i=0,2
        mcs(i)=mmsq_cs(i,+1,+1)
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen