! Test the behaviour of the driver when input is missing or is invalid. Note
! that with the compiler driver (flang-new), the input _has_ to be specified.
! Indeed, the driver decides what "job/command" to create based on the input
! file's extension. No input file means that it doesn't know what to do
! (compile?  preprocess? link?). The frontend driver (flang-new -fc1) simply
! assumes that "no explicit input == read from stdin"

! REQUIRES: new-flang-driver

!--------------------------
! FLANG DRIVER (flang-new)
!--------------------------
! RUN: not %flang  2>&1 | FileCheck %s --check-prefix=FLANG-NO-FILE
! RUN: not %flang %t 2>&1 | FileCheck %s --check-prefix=FLANG-NONEXISTENT-FILE

!-----------------------------------------
! FLANG FRONTEND DRIVER (flang-new -fc1)
!-----------------------------------------
! RUN: not %flang_fc1 %t 2>&1  | FileCheck %s --check-prefix=FLANG-FC1-NONEXISTENT-FILE
! RUN: not %flang_fc1 %S 2>&1  | FileCheck %s --check-prefix=FLANG-FC1-DIR

!-----------------------
! EXPECTED OUTPUT
!-----------------------
! FLANG-NO-FILE: flang-new: error: no input files

! FLANG-NONEXISTENT-FILE: flang-new: error: no such file or directory: {{.*}}
! FLANG-NONEXISTENT-FILE: flang-new: error: no input files

! FLANG-FC1-NONEXISTENT-FILE: error: {{.*}} does not exist
! FLANG-FC1-DIR: error: {{.*}} is not a regular file


text:
           program demo_spread
           implicit none
             integer :: a = 1, b(2) = [ 1, 2 ]
             write(*,*) spread(a, 1, 2)            ! "1 1"
             write(*,*) spread(b, 1, 2)            ! "1 1 2 2"
           end program demo_spread

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


text:
module generation_rate_interface
  !! Encapsulate mass and energy generation rates to facilitate functional programming:
  !! because the energy generation rate depends on the mass generation rate, calculating
  !! both in one function avoids redundancy. Functions can have only one result so a
  !! generation_rate_t object captures both results in one object.
  use kind_parameters, only : rkind
  implicit none

  private

  type, public :: generation_rate_t
    private
    real(rkind) burn_rate_
    real(rkind) e_dot_gen_
    real(rkind) m_dot_gen_
  contains
    procedure :: burn_rate
    procedure :: e_dot_gen
    procedure :: m_dot_gen
  end type

  interface generation_rate_t
    module procedure new_generation_rate
  end interface

  interface

    pure module function new_generation_rate(burn_rate, mass_generation_rate, enthalpy_generation_rate)
      implicit none
      real(rkind), intent(in) :: burn_rate, mass_generation_rate, enthalpy_generation_rate
      type(generation_rate_t) new_generation_rate
    end function

    pure module function burn_rate(this)
      implicit none
      class(generation_rate_t), intent(in) :: this
      real(rkind) burn_rate
    end function

    pure module function m_dot_gen(this)
      implicit none
      class(generation_rate_t), intent(in) :: this
      real(rkind) m_dot_gen
    end function

    pure module function e_dot_gen(this)
      implicit none
      class(generation_rate_t), intent(in) :: this
      real(rkind) e_dot_gen
    end function

  end interface

end module generation_rate_interface

max_stars_repo_name: sourceryinstitute/rocket-science

max_issues_repo_name: sourceryinstitute/rocket-science

max_forks_repo_name: sourceryinstitute/rocket-science


text:
\ STEP? can help if you don't have +LOOP -- an 31jan2022
(*
        limit index DO ... delta +LOOP
   can be coded as
        limit index BEGIN ... delta STEP? UNTIL 2drop
        \ Loop parameters on data stack, I changes to DUP
   or
        limit index BEGIN 2>r ... 2r> delta STEP? UNTIL 2drop
        \ Loop parameters on return stack, I changes to R@
*)

\ ----- Code using WITHIN
: STEP?             ( limit index delta -- limit index' stop? )
    s>d >r          \ "     "     "             r: deltasign
    2dup +          \ "     "     delta idx'
    2over within >r \ "     "     delta         r: deltasign stop?
    + 2r> xor ;     \ "     idx'  stop?         r: -
\ ----- End

\ ----- Code without WITHIN
: STEP?             ( limit index delta -- limit index' stop? )
    s>d over 2>r +  \ limit index'                  r: deltasign delta
    2dup - negate   \ "     "      idx'-lim
    dup r> -        \ "     "      "        idx-lim r: deltasign
    u<              \ "     "      stop?
    r> xor ;        \ delta<0? -> invert stop?      r: -
\ ----- End

(*
\ With DO+LOOP
: T1+ ( x y -- ) do i .      1 +loop ;
: T+  ( x y -- ) do i .      4 +loop ;
: TT+ ( x y -- ) do i u.  2000 +loop ;  \ For 16b forth!
: T1- ( x y -- ) do i .     -1 +loop ;
: T-  ( x y -- ) do i .     -4 +loop ;
: TT- ( x y -- ) do i u. -2000 +loop ;  \ For 16b forth!
*)

\ With STEP?
hex
: T1+ ( x y -- ) begin dup .      1 step? until 2drop ;
: T+  ( x y -- ) begin dup .      4 step? until 2drop ;
: TT+ ( x y -- ) begin dup u.  2000 step? until 2drop ;
: T1- ( x y -- ) begin dup .     -1 step? until 2drop ;
: T-  ( x y -- ) begin dup .     -4 step? until 2drop ;
: TT- ( x y -- ) begin dup .  -2000 step? until 2drop ;

\ ---------- tests
 4 -4    t1+
 25 -4    t+
 24 -4    t+
 23 -4    t+
 F000 0  tt+
 0 0     tt+
\ ----------
 -4 4    t1-
 -25 4    t-
 -24 4    t-
 -23 4    t-
 -F000 0 tt-
 0 0     tt-
 \ ----------

(* Results
@)\ ---------- tests  OK.0
@) 4 -4    t1+ -4 -3 -2 -1 0 1 2 3  OK.0
@) 25 -4    t+ -4 0 4 8 C 10 14 18 1C 20 24  OK.0
@) 24 -4    t+ -4 0 4 8 C 10 14 18 1C 20  OK.0
@) 23 -4    t+ -4 0 4 8 C 10 14 18 1C 20  OK.0
@) F000 0  tt+ 0 2000 4000 6000 8000 A000 C000 E000  OK.0
@) 0 0     tt+ 0 2000 4000 6000 8000 A000 C000 E000  OK.0
@)\ ----------  OK.0
@) -4 4    t1- 4 3 2 1 0 -1 -2 -3 -4  OK.0
@) -25 4    t- 4 0 -4 -8 -C -10 -14 -18 -1C -20 -24  OK.0
@) -24 4    t- 4 0 -4 -8 -C -10 -14 -18 -1C -20 -24  OK.0
@) -23 4    t- 4 0 -4 -8 -C -10 -14 -18 -1C -20  OK.0
@) -F000 0 tt- 0 -2000 -4000 -6000 -8000 6000 4000 2000  OK.0
@) 0 0     tt- 0  OK.0
@) \ ----------  OK.0
*)

\ <><>

ext: f

max_stars_repo_name: embeddingforth/embeddinForth

max_issues_repo_name: embeddingforth/embeddinForth

max_forks_repo_name: embeddingforth/embeddinForth


subroutine control2model(xhat,sval,bval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    control2model
!   prgmmr: tremolet
!
! abstract:  Converts control variable to physical space
!
! program history log:
!   2007-04-13  tremolet - initial code
!   2007-04-27  tremolet - multiply by sqrt(B) (from ckerror D. Parrish)
!   2008-12-04  todling  - update interface to ckgcov; add tsen/p3d
!   2008-12-29  todling  - add call to strong balance contraint
!   2009-04-21  derber   - modify call to getuv to getuv(*,0)
!   2009-08-14  lueken   - update documentation
!   2010-03-15  zhu      - make changes to ckbcov, add assign_cs2array
!   2010-04-28  todling  - update to use gsi_bundle
!   2010-05-31  todling  - better consistency checks; add co/co2
!                        - ready to bypass analysis of (any) meteorological fields
!   2010-06-15  todling  - generalized handling of chemistry
!   2011-05-15  auligne/todling - generalized cloud handling
!   2011-07-12  zhu      - add cw_to_hydro and cwhydromod
!   2011-12-14  mkim     - changed clouds4crtm to clouds in metguess
!   2012-06-25  parrish  - modify wbundle by adding motley variables to control vector
!                          so will be in form expected by new version of ckgcov which uses general_grid2sub.
!   2012-10-09  wgu      - update interface to tbalance (fut2ps)
!   2013-10-25  todling  - nullify work pointers
!   2013-10-28  todling  - rename p3d to prse
!   2013-05-23   zhu     - add ntclen and predt for aircraft temperature bias correction
!
!   input argument list:
!     xhat - Control variable
!     sval - State variable
!     bval - Bias predictors
!
!   output argument list:
!     sval - State variable
!     bval - Bias predictors
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use control_vectors, only: control_vector
use control_vectors, only: cvars3d,cvars2d,cvarsmd
use control_vectors, only: nc2d,mvars
use bias_predictors, only: predictors
use gsi_4dvar, only: nsubwin, l4dvar, lsqrtb
use gridmod, only: lat2,lon2,nsig
use jfunc, only: nsclen,npclen,ntclen
use berror, only: varprd,fpsproj,fut2ps
use balmod, only: balance
use cwhydromod, only: cw2hydro_tl
use gsi_bundlemod, only: gsi_bundlecreate
use gsi_bundlemod, only: gsi_gridcreate
use gsi_bundlemod, only: gsi_grid
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: gsi_bundlegetvar
use gsi_bundlemod, only: gsi_bundledestroy
use gsi_chemguess_mod, only: gsi_chemguess_get
use gsi_metguess_mod, only: gsi_metguess_get
use mpeu_util, only: getindex
use constants, only: max_varname_length
implicit none

! Declare passed variables
type(control_vector), intent(in   ) :: xhat
type(gsi_bundle)    , intent(inout) :: sval(nsubwin)
type(predictors)    , intent(inout) :: bval

! Declare local variables
character(len=*),parameter:: myname ='control2model'
real(r_kind),dimension(lat2,lon2,nsig) :: workst,workvp,workrh
type(gsi_bundle) :: wbundle
type(gsi_grid)   :: grid
integer(i_kind) :: ii,jj,i,ic,id,ngases,nclouds,istatus
character(len=10),allocatable,dimension(:) :: gases
character(len=max_varname_length),allocatable,dimension(:) :: clouds

! Note: The following does not aim to get all variables in
!       the state and control vectors, but rather the ones
!       explicitly needed by this routine.
! Declare required local state variables
logical :: ls_u,ls_v,ls_prse,ls_q,ls_tsen,ls_tv,ls_ps,ls_ql,ls_qi
integer(i_kind), parameter :: nsvars = 9
integer(i_kind) :: isps(nsvars)
character(len=4), parameter :: mysvars(nsvars) = (/  &  ! vars from SV needed here
        'u   ', 'v   ', 'prse', 'q   ', 'tsen',  'tv  ', 'ps  ','ql  ', 'qi  ' /)
character(len=max_varname_length),allocatable,dimension(:) :: cvars2dpm  ! names of 2d fields including
                                                                         !  motley vars (if any)
real(r_kind),pointer,dimension(:,:)   :: sv_ps=>NULL()
real(r_kind),pointer,dimension(:,:)   :: sv_sst=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_u=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_v=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_prse=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_q=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_tsen=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_tv=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_oz=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: sv_rank3=>NULL()
real(r_kind),pointer,dimension(:,:)   :: sv_rank2=>NULL()

logical :: do_balance,do_getprs_tl,do_normal_rh_to_q,do_tv_to_tsen,do_getuv,cw_to_hydro

!******************************************************************************

if (.not.lsqrtb) then
   write(6,*)trim(myname),': assumes lsqrtb'
   call stop2(104)
end if
if (nsubwin/=1 .and. .not.l4dvar) then
   write(6,*)trim(myname),': error 3dvar',nsubwin
   call stop2(105)
end if

!  create extended list of 2d variable names to include motley variables.
!     NOTE: if mvars=0, there are no motley variables so cvars2dpm=cvars2d.
   allocate(cvars2dpm(nc2d+mvars))
   do i=1,nc2d
      cvars2dpm(i)=cvars2d(i)
   end do
   do i=1,mvars
      cvars2dpm(nc2d+i)=cvarsmd(i)
   end do
   call gsi_gridcreate(grid,lat2,lon2,nsig)

! Inquire about chemistry
call gsi_chemguess_get('dim',ngases,istatus)
if (ngases>0) then
    allocate(gases(ngases))
    call gsi_chemguess_get('gsinames',gases,istatus)
endif

! Since each internal vector of xhat has the same structure, pointers are
! the same independent of the subwindow jj
call gsi_bundlegetpointer (sval(1),mysvars,isps,istatus)
ls_u  =isps(1)>0; ls_v   =isps(2)>0; ls_prse=isps(3)>0
ls_q  =isps(4)>0; ls_tsen=isps(5)>0; ls_tv =isps(6)>0
ls_ps =isps(7)>0; ls_ql  =isps(8)>0; ls_qi =isps(9)>0

! Define what to do depending on what's in SV and
! what's explictly needed in this routine
do_balance       =ls_tv.and.ls_ps
do_getprs_tl     =ls_ps.and.ls_tv .and.ls_prse
do_normal_rh_to_q=ls_tv.and.ls_prse.and.ls_q
do_tv_to_tsen    =ls_tv.and.ls_q  .and.ls_tsen
do_getuv         =ls_u .and.ls_v

! Inquire about cloud-variables
cw_to_hydro=.false.
call gsi_metguess_get ('clouds::3d',nclouds,istatus)
if (nclouds>0) then
   allocate(clouds(nclouds))
   call gsi_metguess_get ('clouds::3d',clouds,istatus)
   if (getindex(cvars3d,'cw')>0 .and. ls_ql .and. ls_qi) cw_to_hydro=.true.
end if

! Loop over control steps
do jj=1,nsubwin

!  create an internal structure w/ the same vars as those in the control vector, including motley vars
   call gsi_bundlecreate (wbundle,grid,'control2model work',istatus,names2d=cvars2dpm,names3d=cvars3d)
   if(istatus/=0) then
      write(6,*) trim(myname), ': trouble creating work bundle'
      call stop2(999)
   endif

!  Multiply by sqrt of background error (ckerror)
!  -----------------------------------------------------------------------------
!  Apply sqrt of variance, as well as vertical & horizontal parts of background
!  error

   call ckgcov(xhat%step(jj)%values(:),wbundle,size(xhat%step(jj)%values(:)))

!  Get pointers to required state variables
   call gsi_bundlegetpointer (sval(jj),'u'   ,sv_u,   istatus)
   call gsi_bundlegetpointer (sval(jj),'v'   ,sv_v,   istatus)
   call gsi_bundlegetpointer (sval(jj),'ps'  ,sv_ps,  istatus)
   call gsi_bundlegetpointer (sval(jj),'prse',sv_prse,istatus)
   call gsi_bundlegetpointer (sval(jj),'tv'  ,sv_tv,  istatus)
   call gsi_bundlegetpointer (sval(jj),'tsen',sv_tsen,istatus)
   call gsi_bundlegetpointer (sval(jj),'q'   ,sv_q ,  istatus)
   call gsi_bundlegetpointer (sval(jj),'oz'  ,sv_oz , istatus)
   call gsi_bundlegetpointer (sval(jj),'sst' ,sv_sst, istatus)

!  Copy variables from CV to SV
   call gsi_bundlegetvar ( wbundle, 'sf' , workst, istatus )
   call gsi_bundlegetvar ( wbundle, 'vp' , workvp, istatus )
   call gsi_bundlegetvar ( wbundle, 'q'  , workrh, istatus )
   call gsi_bundlegetvar ( wbundle, 't'  , sv_tv,  istatus )
   call gsi_bundlegetvar ( wbundle, 'oz' , sv_oz,  istatus )
   call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps,  istatus )
   call gsi_bundlegetvar ( wbundle, 'sst', sv_sst, istatus )

!  Same one-to-one map for chemistry-vars; take care of them together
   do ic=1,ngases
      id=getindex(cvars3d,gases(ic))
      if (id>0) then
         call gsi_bundlegetpointer (sval(jj),gases(ic),sv_rank3,istatus)
         call gsi_bundlegetvar     (wbundle, gases(ic),sv_rank3,istatus)
      endif
      id=getindex(cvars2d,gases(ic))
      if (id>0) then
         call gsi_bundlegetpointer (sval(jj),gases(ic),sv_rank2,istatus)
         call gsi_bundlegetvar     (wbundle, gases(ic),sv_rank2,istatus)
      endif
   enddo

!  Balance equation
   if(do_balance) call balance(sv_tv,sv_ps,workst,workvp,fpsproj,fut2ps)

!  -----------------------------------------------------------------------------

!  Get 3d pressure
   if(do_getprs_tl) call getprs_tl(sv_ps,sv_tv,sv_prse)

!  Convert input normalized RH to q
   if(do_normal_rh_to_q) call normal_rh_to_q(workrh,sv_tv,sv_prse,sv_q)

!  Calculate sensible temperature
   if(do_tv_to_tsen) call tv_to_tsen(sv_tv,sv_q,sv_tsen)

!  Convert streamfunction and velocity potential to u,v
   if(do_getuv) call getuv(sv_u,sv_v,workst,workvp,0)

   if (nclouds>0) then
      if (cw_to_hydro) then
!        Case when cw is generated from hydrometeors
         call cw2hydro_tl(sval(jj),wbundle,clouds,nclouds)
      else
!        Case when cloud-vars map one-to-one, take care of them together
         do ic=1,nclouds
            id=getindex(cvars3d,clouds(ic))
            if (id>0) then
                call gsi_bundlegetpointer (sval(jj),clouds(ic),sv_rank3,istatus)
                call gsi_bundlegetvar     (wbundle, clouds(ic),sv_rank3,istatus)
            endif
         enddo
      end if
   end if

!  destroy temporary bundle
   call gsi_bundledestroy(wbundle,istatus)
   if(istatus/=0) then
      write(6,*) trim(myname), ': trouble destroying work bundle'
      call stop2(999)
   endif

end do

! Bias correction terms
do ii=1,nsclen
   bval%predr(ii)=xhat%predr(ii)*sqrt(varprd(ii))
enddo

do ii=1,npclen
   bval%predp(ii)=xhat%predp(ii)*sqrt(varprd(nsclen+ii))
enddo

if (ntclen>0) then
   do ii=1,ntclen
      bval%predt(ii)=xhat%predt(ii)*sqrt(varprd(nsclen+npclen+ii))
   enddo
end if

! Clean up
if (ngases>0) then
    deallocate(gases)
endif

if (nclouds>0) deallocate(clouds)

return
end subroutine control2model


text:
C*************************************************
      subroutine CXSOIL(N1)
C***********************************************************************
C
C   THIS ROUTINE CALCULATES THE complex SOIL PROPERTIES AND TRANSFER
C   FUNCTIONS FOR THE LAYERS
C
C        N1      = NUMBER OF SOIL LAYERS
C        BL      = RATIO OF CRITICAL DAMPING
C        GL      = SHEAR MODULUS
C        R       = DENSITY
C        G       = complex SHEAR MODULUS
C        V       = complex SHEAR WAVE VELOCITY
C        PLUS    = complex TRANSFER FUNCTION
C        MINUS   = complex TRANSFER FUNCTION
C
C   CODED BY PER B SCHNABEL OCT 1971
C
C***********************************************************************
C
      complex G, V, PLUS, MINUS, MU
      character*6 IDNT
      common /SOILA/ IDNT(6),BL(51),GL(51),FACT(51),H(51),R(51),BF(51)
      common /CSOIL/ G(51), V(51), PLUS(51), MINUS(51)
C
      N = N1 + 1
      do 1 I = 1,N
      GIMAG=2.*BL(I)*GL(I)*SQRT(1.-BL(I)*BL(I))
      GREAL=GL(I)*(1.-2.*BL(I)*BL(I))
      G(I)=CMPLX(GREAL,GIMAG)
      V(I) = CSQRT(G(I)/R(I))
    1 continue
      do 2 I = 1,N1
      J = I + 1
      MU = CSQRT(R(I)/R(J)*G(I)/G(J))
      PLUS(I) = (1. + MU)/2.
      MINUS(I)= (1. - MU)/2.
    2 continue
      return
      end


Module Stop3BodyDecays
!---------------------------------------------------------------
! 23.11: adding decay modes as given in Djouadi and Mambrini
!        PRD63, 115005 (2001) hep-ph/0011364
!---------------------------------------------------------------
! load modules
Use Control
Use DecayFunctions, only: ScalarToTwoFermions
Use Mathematics
Use StandardModel, Only: CKM, mW, mW2, mZ2, mf_u, mf_u2, mf_d, mf_d2  &
  & , mf_l, mf_l2, mZ
Use ThreeBodyPhaseSpaceS
! load modules

! private variables
Integer, private :: istosle
Real(dp), Private :: mL2, mb2, mSl2, mst2, gammc(2), gammc2(2), gammsb(2)   &
         & , gammsb2(2), gamtmt, gamtmt2, amc2(2), SumM(3), amN2, msbot2(2) &
         & , mt2
Complex(dp), Private :: CoeffSl(3,5), ca(3,9), cb(2,8), cc(2,2,7), cd(3,6) &
                      & , ce(3,2)

! private routines
Private :: Ftilde
! private variables

Contains

 Real(dp) Function Ftilde(muChi, muPhi, muSq)
 implicit none
  Real(dp), Intent(in) :: muChi, muPhi, muSq

  Real(dp) :: parts(2), rp, rm

  parts(1) = 1._dp + muPhi - MuSq
  parts(2) = Sqrt( (1._dp - muPhi - muSq)**2 - 4._dp * muPhi * muSq)
  rp = 0.5_dp * (parts(1) + parts(2) )
  rm = 0.5_dp * (parts(1) - parts(2) )

  FTilde = - f(1._dp) + f(rp) + f(rm) &
       & + Log(muPhi) * Log((muChi-muSq)/(muChi-1._dp) )

 contains
  real(dp) function f(x)
  implicit none
   Real(dp), intent(in) :: x
    f = Log(muChi - 1._dp + x) * Log((1._dp - muChi)/(muSq-muChi))  &
      & + Li2((muChi-muSq)/(muchi - 1._dp + x) )                    &
      & - Li2((muChi-1._dp)/(muchi - 1._dp + x) )
  end function
 End Function Ftilde

 Real(dp) Function Sq2toSq1Chi1(muChi, muSq)
 Implicit None
  Real(dp), Intent(in) :: muChi(2), muSq

  Real(dp) :: ProdChi

  If (muChi(1).Eq.muChi(2)) Then
    Sq2toSq1Chi1 = (muSq - 1._dp) * (5._dp - 6._dp * muChi(1)                 &
              &                    + 5._dp * muSq- 2._dp * muSq / muChi(1) ) &
              &  - 2._dp * (muSq/muchi(1))**2 * Log( muSq )                  &
              &  + 2._dp * (muSq-muChi(1)) * (muChi(1)-1._dp)                &
              &    * (muSq + muSq*muChi(1) + muChi(1) - 3._dp*muChi(1)**2 )  &
              &    * Log((muChi(1)-1._dp)/(muChi(1)-muSq)) / muChi(1)**2

  Else
   prodChi = muChi(1) * muChi(2)
   Sq2toSq1Chi1 = (muSq - 1._dp) * (3._dp +3._dp * muSq -2._dp * Sum(muChi))  &
              & - 2._dp * muSq**2/ProdChi * Log( muSq )                       &
              &  + ( (muSq-muChi(1))**2 * (muChi(1)-1._dp)**2                 &
              &            * Log((muChi(1)-muSq)/(muChi(1)-1._dp)) / muchi(1) &
              &    -muChi(1) * (muSq-muChi(2))**2 * (muChi(2)-1._dp)**2       &
              &        * Log((muChi(2)-muSq)/(muChi(2)-1._dp)) / muchi(2)     &
              &    ) / (0.5_dp * (muChi(2)-muChi(1)))
  End If

  Sq2toSq1Chi1 = 0.25_dp * Sq2toSq1Chi1

 End Function Sq2toSq1Chi1

 Real(dp) Function Sq2toSq1Chi2(muChi, muSq)
 Implicit None
  Real(dp), Intent(in) :: muChi(2), muSq

  Real(dp) :: prodChi

  If (muChi(1).Eq.muChi(2)) Then
   Sq2toSq1Chi2 = 0.5_dp * (muSq - 1._dp)                                   &
              &   * (muchi(1) - 2._dp * muchi(1)**2 - 2._dp * muSq          &
              &     + muSq * muChi(1) )                                     &
              & + muSq * (1._dp + muSq - muSq/muChi(1) ) * Log( muSq )      &
              & + (muSq- muChi(1)) * (muChi(1)-1._dp) * (muSq-muChi(1)**2)  &
              &   * Log( (muChi(1)-1._dp) / (muChi(1)-muSq) ) / muChi(1)

  Else
   prodChi = muChi(1) * muChi(2)
   Sq2toSq1Chi2 = muSq * (1._dp + muSq - 0.5_dp *muSq * Sum(muChi) / ProdChi) &
             &         * Log(muSq)                                            &
             &  + 0.5_dp * (muSq - 1._dp) * (muSq + ProdChi)                  &
             &  + (muChi(2) * (muSq-muChi(1))**2 * (muChi(1)-1._dp)**2        &
             &    * Log((muChi(1)-muSq)/(muChi(1)-1._dp)) / muchi(1)          &
             &    -muChi(1) * (muSq-muChi(2))**2 * (muChi(2)-1._dp)**2        &
             &    * Log((muChi(2)-muSq)/(muChi(2)-1._dp)) / muchi(2)          &
             &    ) / (2._dp * (muChi(2)-muChi(1)) )

  End If

  Sq2toSq1Chi2 = Sq2toSq1Chi2 / (muChi(1)*muChi(2))

 End Function Sq2toSq1Chi2

 Subroutine Sq2ToSq1ff(mSq, C_Sq2Sq1Z, Lf, Rf, mS0, c_Sq2Sq1S0, c_FFS0 &
           & , mP0, c_Sq2Sq1P0, c_FFP0, mC, C_CFSfp_L, C_CFSfp_R, mN   &
           & , C_FNSf_L, C_FNSf_R, mGlu, C_FGSf_L, C_FGSf_R, Nc        &
           & , SameGeneration, width)
 !--------------------------------------------------------------------------
 ! calculates the decay sq_2 -> sq_1 f \bar{f}
 ! uses the formulas given in A. Djouadi and Y. Mambrini, hep-ph/0011364
 ! 23.11: implementing gauge boson exchange
 !--------------------------------------------------------------------------
 implicit none
  Real(dp), Intent(in) :: mSq(2), Lf, Rf, mS0(2), mP0(2), c_FFS0(2)  &
       &  , c_FFP0(2), mC(2), mN(4), mGlu, Nc
  Complex(dp), Intent(in) :: C_Sq2Sq1Z, c_Sq2Sq1S0(2), c_Sq2Sq1P0(2)     &
       & , C_CFSfp_L(2,2), C_CFSfp_R(2,2), C_FNSf_L(4,2), C_FNSf_R(4,2)  &
       & , C_FGSf_L(2), C_FGSf_R(2)
  logical :: SameGeneration(2)
  Real(dp), intent(out) :: width

  Integer :: i1, i2, i_count
  Real(dp) :: parts(16), muV, muSq, muPhi(2)

  width = 0._dp
  parts = 0._dp
  !------------------
  ! Z boson exchange
  !------------------
  muV = mZ / mSq(2)
  muSq = mSq(1) / mSq(2)
  parts(1) = 2._dp * (Rf**2 + Lf**2 ) * Abs(C_Sq2Sq1Z)**2 &
         &         * Sq2ToSq1VV(muV, muSq)
  !-------------------------------------------------------------------------
  ! h^0 h^0 exchange, note that Higgs-F-F coupling is defined differently
  !-------------------------------------------------------------------------
  muPhi = mS0(1) / mSq(2)
  parts(2) = 4._dp * c_FFS0(1)**2 * Abs(c_Sq2Sq1S0(1))**2    &
           & * Sq2ToSq1PhiPhi(muPhi, muSq)
  !-------------------------------------------------------------------------
  ! h^0 H^0 exchange, note that Higgs-F-F coupling is defined differently
  !-------------------------------------------------------------------------
  muPhi(1) = mS0(1) / mSq(2)
  muPhi(2) = mS0(2) / mSq(2)
  parts(3) = 8._dp * c_FFS0(1) * c_FFS0(2)                    &
           & * Real(c_Sq2Sq1S0(1)*Conjg(c_Sq2Sq1S0(2)),dp)**2    &
           & * Sq2ToSq1PhiPhi(muPhi, muSq)
  !-------------------------------------------------------------------------
  ! H^0 H^0 exchange, note that Higgs-F-F coupling is defined differently
  !-------------------------------------------------------------------------
  muPhi = mS0(2) / mSq(2)
  parts(4) = 4._dp * c_FFS0(2)**2 * Abs(c_Sq2Sq1S0(2))**2    &
           & * Sq2ToSq1PhiPhi(muPhi, muSq)
  !-------------------------------------------------------------------------
  ! A^0 A^0 exchange, note that Higgs-F-F coupling is defined differently
  !-------------------------------------------------------------------------
  muPhi = mP0(2) / mSq(2)
  parts(5) = 4._dp * c_FFP0(2)**2 * Abs(c_Sq2Sq1P0(2))**2    &
           & * Sq2ToSq1PhiPhi(muPhi, muSq)
  !--------------------------------
  ! chargino contribution
  !--------------------------------
  If (SameGeneration(1)) then
   muPhi = mC(1) / mSq(2)
   parts(6) = ( Abs( C_CFSfp_L(1,1) ) * Abs( C_CFSfp_L(1,2) )     &
          &   + Abs( C_CFSfp_R(1,1) ) * Abs( C_CFSfp_R(1,2) )  )  &
          &  * Sq2toSq1Chi1(muPhi, muSq)                          &
          & + Real( C_CFSfp_L(1,2) * C_CFSfp_L(1,1)               &
          &         * Conjg( C_CFSfp_R(1,2) * C_CFSfp_R(1,1) ) ,dp)  &
          &  * Sq2toSq1Chi2(muPhi, muSq)
   muPhi = mC / mSq(2)
   parts(7) = ( Real( C_CFSfp_L(1,1) * C_CFSfp_L(2,2)                &
          &          * Conjg( C_CFSfp_L(1,2) * C_CFSfp_L(2,1) ),dp )    &
          &   + Real( C_CFSfp_R(1,1) * C_CFSfp_R(2,2)                &
          &          * Conjg( C_CFSfp_R(1,2) * C_CFSfp_R(2,1) ),dp )  ) &
          & * Sq2toSq1Chi1(muPhi, muSq)                              &
          & + ( Real( C_CFSfp_L(1,2) * C_CFSfp_L(2,1)                &
          &         * Conjg( C_CFSfp_R(1,1) * C_CFSfp_R(2,2) ),dp )     &
          &   +  Real( C_CFSfp_R(1,2) * C_CFSfp_R(2,1)               &
          &         * Conjg( C_CFSfp_L(1,1) * C_CFSfp_L(2,2) ),dp) )   &
          &  * Sq2toSq1Chi2(muPhi, muSq)
   parts(7) = 2._dp * parts(7)
   muPhi = mC(2) / mSq(2)
   parts(8) = ( Abs( C_CFSfp_L(2,1) ) * Abs( C_CFSfp_L(2,2) )     &
          &   + Abs( C_CFSfp_R(2,1) ) * Abs( C_CFSfp_R(2,2) )  )  &
          & * Sq2toSq1Chi1(muPhi, muSq)                           &
          & + Real( C_CFSfp_L(2,2) * C_CFSfp_L(2,1)               &
          &         * Conjg( C_CFSfp_R(2,2) * C_CFSfp_R(2,1) ),dp )  &
          &  * Sq2toSq1Chi2(muPhi, muSq)
  end if
  !--------------------------------
  ! neutralino contribution
  !--------------------------------
  If (SameGeneration(2)) then
   i_count = 5
   Do i1=1,4
    muPhi = mN(i1) / mSq(2)
    parts(i_count) = ( Abs( C_FNSf_L(i1,1) ) * Abs( C_FNSf_L(i1,2) )     &
                 &   + Abs( C_FNSf_R(i1,1) ) * Abs( C_FNSf_R(i1,2) )  )  &
                 & * Sq2toSq1Chi1(muPhi, muSq)                           &
                 & + Real( C_FNSf_L(i1,2) * C_FNSf_L(i1,1)               &
                 &         * Conjg( C_FNSf_R(i1,2) * C_FNSf_R(i1,1) ),dp )  &
                 &  * Sq2toSq1Chi2(muPhi, muSq)
    i_count = i_count + 1
    Do i2=i1+1,4
     muPhi(1) = mN(i1) / mSq(2)
     muPhi(2) = mN(i2) / mSq(2)
     parts(i_count) = ( Real( C_FNSf_L(i1,1) * C_FNSf_L(i2,2)        &
          &          * Conjg( C_FNSf_L(i1,2) * C_FNSf_L(i2,1) ) )    &
          &   + Real( C_FNSf_R(i1,1) * C_FNSf_R(i2,2)                &
          &          * Conjg( C_FNSf_R(i1,2) * C_FNSf_R(i2,1) ),dp)) &
          & * Sq2toSq1Chi1(muPhi, muSq)                              &
          & + ( Real( C_FNSf_L(i1,2) * C_FNSf_L(i2,1)                &
          &         * Conjg( C_FNSf_R(i1,1) * C_FNSf_R(i2,2) ),dp )     &
          &   +  Real( C_FNSf_R(i1,2) * C_FNSf_R(i2,1)               &
          &         * Conjg( C_FNSf_L(i1,1) * C_FNSf_L(i2,2) ),dp ) )   &
          &  * Sq2toSq1Chi2(muPhi, muSq)
     parts(i_count) = 2._dp * parts(i_count)
     i_count = i_count + 1
    end do
   end do
  end if
  !----------------------
  ! gluino contribution
  !----------------------
  If (SameGeneration(2)) then
   muPhi = mGlu / mSq(2)
   parts(i_count) = ( Abs( C_FGSf_L(1) ) * Abs( C_FGSf_L(2) )     &
                &   + Abs( C_FGSf_R(1) ) * Abs( C_FGSf_R(2) )  )  &
                & * Sq2toSq1Chi1(muPhi, muSq)                           &
                & + Real( C_FGSf_L(2) * C_FGSf_L(1)               &
                &         * Conjg( C_FGSf_R(2) * C_FGSf_R(1) ),dp )  &
                &  * Sq2toSq1Chi2(muPhi, muSq)
  end if

  width = oo128pi3 * mSq(2) * Nc * Sum(parts)

 end subroutine Sq2ToSq1ff

 Real(dp) Function Sq2ToSq1PhiChi(muChi, muPhi, muSq)
 !-------------------------------------------------------------------
 ! calculates the gaugino - scalar interference for sfermion decays
 ! input: muChi .... m_gaugino / m_sfermion_2
 !        muPhi .... m_scalar / m_sfermion_2
 !        muSq ..... m_sfermion_1 /  m_sfermion_2
 ! uses: function Ftilde
 ! written by Werner Porod, 09.12.02
 !-------------------------------------------------------------------
 implicit none
  Real(dp), Intent(in) :: muChi, muPhi, muSq

  Sq2ToSq1PhiChi = muSq - 1._dp + muSq * Log(muSq) / muChi  &
        & + (muSQ-muchi)*(muChi-1._dp)/muChi                &
        &    *  Log( (muChi-1._dp) / (muChi-muSq) )         &
        & + muPhi * Ftilde(muChi, muPhi, muSq)

 End Function Sq2ToSq1PhiChi

 Real(dp) Function Sq2ToSq1PhiPhi(muPhi, muSq)
 implicit none
  Real(dp), intent(in) :: muPhi(2), muSq

  Real(dp) :: lam1, sq_lam1, lam2, sq_lam2

  if (muPhi(1).eq.muPhi(2)) then
   lam1 = -1._dp + 2._dp * (muPhi(1) + muSq) - (muSq - muPhi(1))**2
   sq_lam1 = Sqrt(abs(lam1))
   Sq2ToSq1PhiPhi = 2._dp * (muSq - 1._dp)                                 &
                & + 0.5_dp * (1._dp + muSq - Sum(muPhi)) * Log(muSq)       &
                & - ( muPhi(1) * (1._dp+muSq-muPhi(1)) + lam1)             &
                &   * ATan( (1-muSq) * sq_lam1                             &
                &         / (muPhi(1)*(1._dp - muPhi(1) + muSq) - lam1) )  &
                &   / sq_lam1
  else
   lam1 = -1._dp + 2._dp * (muPhi(1) + muSq) - (muSq - muPhi(1))**2
   sq_lam1 = Sqrt(abs(lam1))
   lam2 = -1._dp + 2._dp * (muPhi(2) + muSq) - (muSq - muPhi(2))**2
   sq_lam2 = Sqrt(abs(lam2))
   Sq2ToSq1PhiPhi = muSq - 1._dp                                           &
                & + 0.5_dp * (1._dp + muSq - Sum(muPhi)) * Log(muSq)       &
                & + ( muPhi(1) * sq_lam1                                   &
                &   * ATan( (1-muSq) * sq_lam1                             &
                &         / (muPhi(1)*(1._dp - muPhi(1) + muSq) - lam1) )  &
                &   - muPhi(2) * sq_lam2                                   &
                &   * ATan( (1-muSq) * sq_lam2                             &
                &         / (muPhi(2)*(1._dp - muPhi(2) + muSq) - lam2) )  &
                &   ) / (muPhi(2) - muPhi(1))

  end if

 end Function Sq2ToSq1PhiPhi

 Real(dp) Function Sq2ToSq1VChi(muChi, muV, muSq)
 !-------------------------------------------------------------------
 ! calculates the gaugino - vector boson interference for sfermion decays
 ! input: muChi .... m_gaugino / m_sfermion_2
 !        muV ...... m_V / m_sfermion_2
 !        muSq ..... m_sfermion_1 /  m_sfermion_2
 ! uses: function Ftilde
 ! written by Werner Porod, 09.12.02
 !-------------------------------------------------------------------
 implicit none
  Real(dp), Intent(in) :: muChi, muV, muSq

  real(dp) :: lam, sqlam

  lam = -1._dp + 2._dp * (muV + muSq) - (muV - muSq)**2
  sqlam = Sqrt(lam)
  Sq2ToSq1VChi = (muSq - 1._dp) * ( 0.25_dp * (1._dp+muSQ+2._dp*muV)          &
        &                         + (1._dp+muSq-2._dp*(muChi+muV))*Log(muV))  &
        & - (muSQ-muchi)*(muChi-1._dp) *  Log( (muChi-1._dp) / (muChi-muSq) ) &
        & - 0.25_dp * (lam + 2._dp * muchi * (1._dp - muV + muSq)             &
        &             - 6._dp * muSq) * Log(muSq)                             &
        & + 0.5_dp * sqlam * (1._dp + muSq - 2._dp + muChi - muV)             &
        &     * Atan( (1._dp - muSq) * sqlam / (lam + muV*(muV-1._dp-muSq)) ) &
        & + muV * Ftilde(muChi, muV, muSq)

 End Function Sq2ToSq1VChi

 real(dp) function Sq2ToSq1VV(muV, muSq)
 implicit none
  Real(dp), intent(in) :: muV, muSq
  Real(dp) :: lam, sq_lam

  lam = -1._dp + 2._dp * (muV + muSq) - (muSq - muV)**2
  sq_lam = Sqrt(abs(lam))
  Sq2ToSq1VV = (1-muSq) * ( 5._dp * (1._dp + muSq) - 4._dp * muV  &
           &              + 2._dp * lam / muV) / 3._dp            &
           & + (lam - 2._dp * muSq) * Log( muSq )                 &
           & + 2._dp * (1._dp - muV + muSq) * sq_lam              &
           &         * Atan( (1-muSq) * sq_lam                    &
           &               / (muV*(1._dp - muV + muSq) - lam) )
  Sq2ToSq1VV = 0.25_dp * Sq2ToSq1VV

 end function Sq2ToSq1VV

 Real(dp) Function stbWchi(s)
 !-----------------------------------------------------
 ! Integrand for stop_1 -> b + W+ + neutralino_1
 ! The following functions are called:
 !   -dambda
 !   -I2dsmg1tmg2
 !    I2dtm1g1
 !    I2dtm12g12
 !    I2tdsmg1tmg2
 !    I2tdtm1g1
 !    I2tdtm12g12
 !    I2t2dtm1g1
 !    I2t2dtm12g12
 ! written by Werner Porod
 ! 21.01.02: portation to f90
 !-----------------------------------------------------
  Implicit None

  Real(dp), Intent(in) :: s

  Integer j, k
  Real(dp) :: s2, s3, ws, tmin, tmax, vwert, uwert, fakt
  Complex(dp) :: intI(4,3)

  Iname = Iname + 1
  NameOfUnit(Iname) = "StBWChi"

  s2 = s*s
  s3 = s2*s

  ws = kappa(s, mst2, mw2) * kappa(s, amN2, mb2) / (2._dp*s)
 !-------------------------
 ! Range of t-Integration
 !-------------------------
  tmin = 0.5_dp * (sumM(1)-s) - (mst2-mW2)*(amN2-mb2) / (2._dp*s) - ws
  tmax = tmin + 2._dp * ws

  stbWchi = 0._dp
 !--------------------------------------------------
 ! t-Integrals for chargino-exchange,
 ! chargino exchange: j=1 ... chargino_1 chargino_1
 !                    j=2 ... chargino_2 chargino_2
 !                    j=3 ... chargino_1 chargino_2
 !--------------------------------------------------
  Do j=1,2
   vwert = sumM(j+1)-s
   intI(j,1) = I2dtm1g1(tmin,tmax,vwert,gammc2(j))
   intI(j,2) = I2tdtm1g1(tmin,tmax,vwert,gammc2(j))
   intI(j,3) = I2t2dtm1g1(tmin,tmax,vwert,gammc2(j))
  End Do

  vwert = sumM(2)-s
  uwert = sumM(3)-s
  intI(3,1) = I2dtm12g12(tmin,tmax,vwert,gammc(1),uwert,gammc(2))
  intI(3,2) = I2tdtm12g12(tmin,tmax,vwert,gammc(1),uwert,gammc(2))
  intI(3,3) = I2t2dtm12g12(tmin,tmax,vwert,gammc(1),uwert,gammc(2))

!Write(40,*) " "
!Write(40,123) s,ws,tmin,tmax,vwert
!123 Format(5f15.6)
  Do j = 1,3
   stbWchi = stbWchi                                                      &
         & + ( ca(j,1) + ca(j,2)*s + ca(j,3)*s2 + ca(j,4)*s3) * intI(j,1) &
         & + ( ca(j,5) + ca(j,6)*s + ca(j,7)*s2) * intI(j,2)              &
         & + ( ca(j,8) + ca(j,9)*s) * intI(j,3)
  End Do
 !-------------------
 ! chargino_j top interference term
 !-------------------
  Do j=1,2
   uwert = sumM(j+1) - s
   intI(j,1) = I2dtm12g12(tmin,tmax,uwert,-gammc(j),mt2,gamtmt)
   intI(j,2) = I2tdtm12g12(tmin,tmax,uwert,-gammc(j),mt2,gamtmt)
   intI(j,3) = I2t2dtm12g12(tmin,tmax,uwert,-gammc(j),mt2,gamtmt)

 !--------------------------------------------------
 ! neagitive sign because chargino is u-channel
 !--------------------------------------------------
   stbWchi = stbWchi - ( cb(j,1) + cb(j,2)*s + cb(j,3)*s2) * intI(j,1)   &
           &         - ( cb(j,4) + cb(j,5)*s + cb(j,6)*s2) * intI(j,2)   &
           &         - ( cb(j,7) + cb(j,8)*s) * intI(j,3)
  End Do
 !-------------------
 ! chargino_j sbottom_k interference term
 !-------------------
  Do j=1,2
   uwert = sumM(j+1) - s
   Do k=1,2
    intI(j,1) = I2dsmg1tmg2(s,tmin,tmax,msbot2(k),gammsb(k),uwert,-gammc(j))
    intI(j,2) = I2tdsmg1tmg2(s,tmin,tmax,msbot2(k),gammsb(k),uwert,-gammc(j))

 !--------------------------------------------------
 ! negative sign because chargino is u-channel
 !--------------------------------------------------
    stbWchi = stbWchi                                                       &
          & - (cc(k,j,1)+cc(k,j,2)*s+cc(k,j,3)*s2+cc(k,j,4)*s3) * intI(j,1) &
          & - ( cc(k,j,5) + cc(k,j,6)*s + cc(k,j,7)*s2) * intI(j,2)
   End Do
  End Do
 !--------------------------------------------------
 ! t-Integrals for top-top term,
 !--------------------------------------------------
  intI(1,1) = I2dtm1g1(tmin,tmax,mt2,gamtmt2)
  intI(1,2) = I2tdtm1g1(tmin,tmax,mt2,gamtmt2)
  intI(1,3) = I2t2dtm1g1(tmin,tmax,mt2,gamtmt2)
!j=1
!Write(40,*) j,Real(inti(j,:))
!Write(*,*) "s,tmin,tmax,mt2,gamtmt2,intI(1,:)"
!Write(42,123) s,tmin,tmax,mt2,gamtmt2,real(intI(1,:))
!123 format(4e16.7)
!write(42,*) " "
!stop 99
  stbWchi = stbWchi + (cd(1,1) + cd(1,2)*s ) * intI(1,1)  &
          &         + (cd(1,3) + cd(1,4)*s ) * intI(1,2)  &
          &         + (cd(1,5) + cd(1,6)*s ) * intI(1,3)
 !-------------------
 ! top-sbottom_k interference term
 !-------------------
  Do k=1,2
   intI(k,1) = I2dsmg1tmg2(s,tmin,tmax,msbot2(k),gammsb(k),mt2,gamtmt)
   intI(k,2) = I2tdsmg1tmg2(s,tmin,tmax,msbot2(k),gammsb(k),mt2,gamtmt)

   stbWchi = stbWchi + (cd(k+1,1) + cd(k+1,2)*s + cd(k+1,3)*s2) * intI(k,1) &
           &         + (cd(k+1,4) + cd(k+1,5)*s + cd(k+1,6)*s2) * intI(k,2)
  End Do
 !--------------------------------------------------
 ! t-Integrals for sbottom_k-sbottom_j terms,
 ! difference if decay-width of sbottom_k is known
 !--------------------------------------------------
 ! k=j :
 !-------
  fakt = kappa(s,mst2,mW2)**2 / mW2

  Do k=1,2
   If (gammsb(k).Eq.0.0_dp) Then
    intI(k,1) = fakt * (ce(k,1) + ce(k,2)*s) / (s-msbot2(k))**2
   Else
    intI(k,1) = fakt * (ce(k,1) + ce(k,2)*s) / ((s-msbot2(k))**2 + gammsb2(k))
   End If
  End Do
 !-----------
 ! k=1, j=2 :
 !-----------
  If ((gammsb(1).Eq.0._dp).Or.(gammsb(2).Eq.0)) Then
   intI(3,1) = fakt * (ce(3,1) + ce(3,2)*s) / ((s-msbot2(1))*(s-msbot2(2)))
  Else
   intI(3,1) = fakt * (ce(3,1) + ce(3,2)*s)                                 &
           &    * ( (s-msbot2(1)) * (s-msbot2(2)) + gammsb(1)*gammsb(2) )   &
           &    / ( ((s-msbot2(1))**2 + gammsb(1)**2)                       &
           &        * ((s-msbot2(2))**2 + gammsb(2)**2) )
  End If

  stbWchi = stbWchi + 2._dp * (intI(1,1) + intI(2,1) + intI(3,1) ) * ws

  Iname = Iname - 1

 End Function stbWchi

 Subroutine StopDecays3(id_l, id_nu, n_su, n_sd, n_sle, n_snu              &
        & , id_d, id_W, n_n, n_c, Sup, Sdown, Chi0, g, c_UNSu_L, c_UNSu_R  &
        & , c_DNSd_L, c_DNSd_R, c_SdSuW, Sneut, ChiPm, c_CLSn_L, c_CLSn_R  &
        & , c_CDSu_L, c_CDSu_R, c_CNW_L, c_CNW_R, Slept, c_CNuSl_R, prec   &
        & , Check_Real_States)
 !-----------------------------------------------------------------------
 ! calculates the 3-body stop decays
 ! written by Werner Porod, 08.01.02
 ! 22.09.03: including the possiblity to check for real intermediates states
 ! 20.09.10: adapting to new variable types
 !-----------------------------------------------------------------------
 Implicit None
  Integer, Intent(in) :: id_l(:), id_nu(:), n_su, n_sd, n_sle  &
     & , n_snu, id_d(:), id_W(:), n_n, n_c
  Real(dp), Intent(in) :: prec, g
  Complex(dp), Intent(in) :: c_CLSn_L(:,:,:), c_CLSn_R(:,:,:), c_CDSu_L(:,:,:) &
     & , c_CDSu_R (:,:,:), c_CNuSl_R(:,:,:), c_UNSu_L(:,:,:), c_UNSu_R(:,:,:)  &
     & , c_CNW_L(:,:), c_CNW_R(:,:), c_DNSd_L(:,:,:), c_DNSd_R(:,:,:)          &
     & , c_SdSuW(:,:)
  Logical, Intent(in) :: Check_Real_States
  Type(particle2), Intent(in) :: Sdown(:)
  Type(particle23), Intent(in) :: Chi0(:), ChiPm(:), Slept(:), Sneut(:)
  Type(particle23), Intent(inout) :: Sup(:)

  Integer :: i1, i2, i3, i, j, k, i_c
  Real(dp) ::  mStop(n_su), mStop2(n_su), mSbottom2(n_sd) &
     & , mN(n_n), mN2(n_n), mSneut2(n_snu), mC(n_c), mC2(n_c), g_C(n_c)  &
     & , mSlept2(n_sle), g_sb(n_sd)
  Real(dp) :: kl2kt2, kl2lt2 &
     & , kll2, ktl2, mbmcha, mtaumc, mbmtau, mbmcha1                          &
     & , mbmcha2, smin, smax, Abslt11Sq, Abskt11Sq, AbslL11Sq, mSl4, f_stop   &
     & , qlqr, ql2qr2, kl, k2l2, lqlkqr, kqllqr,  mbchp1, mchp1chi, mbchi     &
     & , mbt, mb4, mchi2, mchi4, mchit, mchp1chp2, mst4   &
     & , mb2chi2, mb2st2, mchi2st2, oomw, mbomw, mbtomw, mbchiomw, mchitomw   &
     & , mstomw, mchiomw, kt1l1CQliCQri, kt1Cl1QliCQri, mbchp2, mchp2chi      &
     & , mchp1t, mbchp1omw, mchp1chiomw, ft2, ht2, fht, mSb(2)
  Complex(dp) :: lt11, kt11  &
     & , lt12C, kt12C, lL11C, kL11C, lL12, kL12, mcha1mcha2, kt11C, lt11C   &
     & , l1l2Ql2Qr1, k1l2Ql2Qr1, k2l1Ql2Qr1, k1k2Ql2Qr1, l1l2Ql1Ql2         &
     & , k1l2Ql1Ql2, k2l1Ql1Ql2, k1k2Ql1Ql2, fltQl1, fktQl1, fltQr1, fktQr1 &
     & , hltQl1, hktQl1, hltQr1, hktQr1, hbltQl1, fbltQl1, hbktQl1, fbktQl1 &
     & , ffb1, fhb1, hfb1, hhb1
  Logical :: calc

  Iname = Iname + 1
  NameOfUnit(Iname) = "StopDecays3"

  mStop = Sup%m
  mStop2 = Sup%m2
  mSbottom2 = Sdown%m2
  g_sb = Sdown%g
  mSneut2 = Sneut%m2
  mSlept2 = Slept%m2
  mN = Chi0%m
  mN2 = Chi0%m2
  mC = ChiPm%m
  mC2 = ChiPm%m2
  g_C = ChiPm%g

  Do i1=1,6
   Sup(i1)%gi3 = 0._dp
   Sup(i1)%bi3 = 0._dp
  End Do

  amc2 = mC2
  MSbot2 = mSbottom2(5:6)
  If (Check_Real_States) Then
   GamMSB2 = 0._dp
   gammc = 0._dp
   gamtmt2 = 0._dp
  Else
   GamMSB2 = mSbottom2(5:6) * g_sb(5:6)**2
   gammc = mC * g_C
   gamtmt2 = 2737.37816_dp
  End If
  gammc2 = gammc**2
  gamtmt = Sqrt(gamtmt2)
  GamMSB = Sqrt( GamMSB2 )
!  Write(*,*) "hier",mc,g_c,gammc

  i_c = 1
  If (GenerationMixing) Then
  Else
   f_stop = 2._dp * oo512pi3 / mStop(5)**3
  End If
  !---------------------------------------------
  ! decays into d-quarks, sneutrinos and leptons
  !---------------------------------------------
  istosle = 1 ! needed for integration routine

  If (GenerationMixing) Then
  Else
   Sup(5)%gi3 = 0._dp

   lt11 = c_CDSu_R(1,3,5)
   kt11 = c_CDSu_L(1,3,5)
   lt12C = Conjg(c_CDSu_R(2,3,5))
   kt12C = Conjg(c_CDSu_L(2,3,5) )
   mst2 = mStop2(5)
   mcha1mcha2 = mC(1) * mC(2)
   mb2 = mf_d2(3)
   mbmcha1 = mf_d(3) * mC(1)
   mbmcha2 = mf_d(3) * mC(2)

   Do i1=1,3 ! lepton generation
    mSl2 = mSneut2(i1)
    !-----------------------------------------
    ! checking kinematics
    !-----------------------------------------
    If (mStop(5).Lt.(Sqrt(mSl2)+mf_d(3)+mf_l(i1))) Cycle

    mbmtau = mf_d(3) * mf_l(i1)
    coeffsl = 0._dp
    Do i2=1,2 ! number of charginos
     kl2kt2 = 0.5_dp * (Abs(c_CLSn_L(i2,i1,i1))**2 * Abs(c_CDSu_L(i2,3,5))**2 &
          &            + Abs(c_CLSn_R(i2,i1,i1))**2 * Abs(c_CDSu_R(i2,3,5))**2)
     kl2lt2 = 0.5_dp * (Abs(c_CLSn_L(i2,i1,i1))**2 * Abs(c_CDSu_R(i2,3,5))**2 &
          &            + Abs(c_CLSn_R(i2,i1,i1))**2 * Abs(c_CDSu_L(i2,3,5))**2)
     kll2 = 2._dp * Real( Conjg(c_CDSu_L(i2,3,5)) * c_CDSu_R(i2,3,5),dp )   &
        & * ( Abs(c_CLSn_L(i2,i1,i1))**2 + Abs(c_CLSn_R(i2,i1,i1))**2 )
     ktl2 = 2._dp * Real(Conjg(c_CLSn_L(i2,i1,i1)) * c_CLSn_R(i2,i1,i1)) &
        &   * (Abs(c_CDSu_L(i2,3,5))**2 + Abs(c_CDSu_R(i2,3,5))**2)

     mbmcha = mf_d(3) * mC(i2)
     mtaumc = mf_l(i1) * mC(i2)

     CoeffSl(i2,1) = - kll2 * mbmcha - kl2lt2 * mC2(i2) - ktl2 * mtaumc     &
        &  + kl2kt2 * (mstop2(5) + mSneut2(i1) - mf_d2(3)-mf_l2(i1))        &
        &  - 4._dp * Real( Conjg(c_CDSu_L(i2,3,5) * c_CLSn_R(i2,i1,i1))     &
        &                * c_CDSu_R(i2,3,5) * c_CLSn_L(i2,i1,i1),dp ) * mbmtau

     CoeffSl(i2,2) = kl2lt2 * mC2(i2) * (mf_d2(3)-mstop2(5))  &
        &                  * (mSneut2(i1)-mf_l2(i1))

     CoeffSl(i2,3) = kl2kt2 * (mstop2(5)-mf_d2(3))*(mf_l2(i1)-mSneut2(i1))    &
        &    + kl2lt2 * mC2(i2) * (mstop2(5)+mSneut2(i1)-mf_d2(3)-mf_l2(i1))  &
        &    + ktl2 * mtaumc * (mstop2(5)-mf_d2(3))                           &
        &    + kll2 * mbmcha * (mSneut2(i1)-mf_l2(i1))                        &
        &    - 4._dp * Real( Conjg(c_CDSu_L(i2,3,5) * c_CLSn_L(i2,i1,i1))     &
        &         * c_CDSu_R(i2,3,5) * c_CLSn_R(i2,i1,i1),dp ) *mbmtau *mC2(i2)

     CoeffSl(i2,4) = - kl2kt2
    End Do

    ! this results from a Mathematica output
    lL11C = Conjg(c_CLSn_R(1,i1,i1))
    kL11C = Conjg(c_CLSn_L(1,i1,i1))
    lL12 = c_CLSn_R(2,i1,i1)
    kL12 = c_CLSn_L(2,i1,i1)
    mL2 = mf_l2(i1)

    CoeffSl(3,1) =  - kL11C*kL12*kt11*kt12C*mb2 -                 &
     &  lL11C*lL12*lt11*lt12C*mb2 -                               &
     &  2*kL11C*kL12*kt12C*lt11*mbmcha1 -                         &
     &  2*kt11*lL11C*lL12*lt12C*mbmcha1 -                         &
     &  2*kt12C*lL11C*lL12*lt11*mbmcha2 -                         &
     &  2*kL11C*kL12*kt11*lt12C*mbmcha2 -                         &
     &  kt11*kt12C*lL11C*lL12*mcha1mcha2 -                        &
     &  kL11C*kL12*lt11*lt12C*mcha1mcha2 +                        &
     &  kL11C*kL12*kt11*kt12C*mL2 + lL11C*lL12*lt11*lt12C*mL2 +   &
     &  kL11C*kL12*kt11*kt12C*mSl2 +                              &
     &  lL11C*lL12*lt11*lt12C*mSl2 +                              &
     &  kL11C*kL12*kt11*kt12C*mst2 + lL11C*lL12*lt11*lt12C*mst2

    CoeffSl(3,2) = -(kt11*kt12C*lL11C*lL12 + kL11C*kL12*lt11*lt12C) &
      &    * mcha1mcha2 * (mL2 - mSl2) * (mb2 - mst2)

    CoeffSl(3,3) = kt11*(2*lt12C*(lL11C*lL12*mbmcha1 + kL11C*kL12*mbmcha2)*  &
     &      mSl2 + kt12C*(-(lL11C*lL12*mb2*mcha1mcha2) +                     &
     &        kL11C*kL12*mb2*mL2 - lL11C*lL12*mcha1mcha2*mL2 +               &
     &        kL11C*kL12*mb2*mSl2 +                                          &
     &        lL11C*lL12*mcha1mcha2*mSl2 +                                   &
     &        lL11C*lL12*mcha1mcha2*mst2 -                                   &
     &        kL11C*kL12*mL2*mst2 - kL11C*kL12*mSl2*mst2)) +                 &
     &  lt11*(kL11C*kL12*(-(lt12C*mb2*mcha1mcha2) -                          &
     &        lt12C*mcha1mcha2*mL2 + 2*kt12C*mbmcha1*mSl2 +                  &
     &        lt12C*mcha1mcha2*mSl2 + lt12C*mcha1mcha2*mst2) +               &
     &     lL11C*lL12*(lt12C*mb2*mL2 + lt12C*mb2*mSl2 +                      &
     &        2*kt12C*mbmcha2*mSl2 - lt12C*mL2*mst2 -                        &
     &        lt12C*mSl2*mst2))

    CoeffSl(3,4) = -(kL11C*kL12*kt11*kt12C + lL11C*lL12*lt11*lt12C)

    smin = (Sqrt(msneut2(i1)) + mf_l(i1))**2
    smax = (mstop(5) - mf_d(3))**2

    calc = .True.
    If (Check_Real_States) Then
     If (mStop(5).Gt.(mC(2)+mf_d(3))) Then  ! all intermediate states are real
      calc = .False.
     Else If (mStop(5).Gt.(mC(1)+mf_d(3))) Then ! only chi^+_2 contributes
      CoeffSl(1,:) = 0._dp
      CoeffSl(3,:) = 0._dp
     End If
    End If

    If (calc) Then
      Sup(5)%gi3(i_c) = f_stop * Dgauss(stoslep,smin,smax,prec)
      Sup(5)%id3(i_c,1) = Sneut(i1)%id
      Sup(5)%id3(i_c,2) = id_l(i1) + 1
      Sup(5)%id3(i_c,3) = id_d(3)
      i_c = i_c +1
    End If

   End Do
  End If

  !---------------------------------------------
  ! decays into d-quarks, neutrinos and sleptons
  !---------------------------------------------
  istosle = 2 ! needed for integration routine

  If (GenerationMixing) Then
  Else
   mb2 = mf_d2(3)
   mst2 = mstop2(5)

   mst2 = mstop2(5)
   mcha1mcha2 = mC(1) * mC(2)
   mb2 = mf_d2(3)
   mbmcha1 = mf_d(3) * mC(1)
   mbmcha2 = mf_d(3) * mC(2)

   Do i1=1,3 ! lepton generation
    Do i2=1,2 ! number of sleptons
     mSl2 = mSlept2(2*(i1-1)+i2)
     !-----------------------------------------
     ! checking kinematics
     !-----------------------------------------
     If (mStop(5).Lt.(Sqrt(mSl2)+mf_d(3))) Cycle
     mSl4 = mSl2**2

     Do i3=1,2 ! number of charginos
      lt11 = c_CDSu_R(i3,3,5)
      kt11 = c_CDSu_L(i3,3,5)
      lt11C = Conjg(c_CDSu_R(i3,3,5))
      kt11C = Conjg(c_CDSu_L(i3,3,5))
      Abslt11Sq = lt11 * lt11C
      Abskt11Sq = kt11 * kt11C
      mbmcha1 = mf_d(3) * mC(i3)
      AbslL11Sq = Abs(c_CNuSl_R(i3,i1,2*(i1-1)+i2))**2

      CoeffSl(i3,1) = 0.5_dp * Abslt11Sq * mC2(i3)                  &
                  &          * (mst2 + 2._dp * mSl2 - mb2)          &
                  & + Abskt11Sq*mSl2*(mb2 - mst2 - 0.5_dp * mSl2)   &
                  & + 4._dp*Real(kt11C*lt11,dp)*mbmcha1*mSl2

      CoeffSl(i3,2) = 0.5_dp * Abslt11Sq*mC2(i3)*mSl4*(mst2 - mb2)

      CoeffSl(i3,3) = - 0.5_dp * mSl2*(mSl2*                 &
       &        (2*(kt11C*lt11 + kt11*lt11C)*mbmcha1 +                  &
       &          Abskt11Sq*(mb2 - mst2)) +                             &
       &       Abslt11Sq*mC2(i3)*(-2*mb2 + mSl2 + 2*mst2))

      CoeffSl(i3,4) = - 0.5_dp * (4*Real(kt11C*lt11,dp)*mbmcha1        &
       &       + Abslt11Sq*mC2(i3) + Abskt11Sq*(mb2 - 2*mSl2 - mst2))

      CoeffSl(i3,5) = - 0.5_dp * Abskt11Sq
      coeffSl(i3,:) = AbslL11Sq * coeffSl(i3,:)
     End Do

     lt11 = c_CDSu_R(1,3,5)
     kt11 = c_CDSu_L(1,3,5)
     lt11C = Conjg(c_CDSu_R(1,3,5))
     kt11C = Conjg(c_CDSu_L(1,3,5))
     lt12C = Conjg(c_CDSu_R(2,3,5))
     kt12C = Conjg(c_CDSu_L(2,3,5) )
     Abslt11Sq = lt11 * lt11C
     Abskt11Sq = kt11 * kt11C
     lL11C = Conjg(c_CNuSl_R(1,i1,2*(i1-1)+i2))
     lL12 = c_CNuSl_R(2,i1,2*(i1-1)+i2)

     CoeffSl(3,1) = -(lL11C*lL12*lt11*lt12C*mb2*mcha1mcha2) +  &
       &  2*kt11*kt12C*lL11C*lL12*mb2*mSl2 +                   &
       &  4*kt12C*lL11C*lL12*lt11*mbmcha1*mSl2 +               &
       &  4*kt11*lL11C*lL12*lt12C*mbmcha2*mSl2 +               &
       &  2*lL11C*lL12*lt11*lt12C*mcha1mcha2*mSl2 -            &
       &  kt11*kt12C*lL11C*lL12*mSl4 +                         &
       &  lL11C*lL12*lt11*lt12C*mcha1mcha2*mst2 -              &
       &  2*kt11*kt12C*lL11C*lL12*mSl2*mst2

     CoeffSl(3,2) = lL11C*lL12*lt11*lt12C*mcha1mcha2*mSl4*(mst2 -mb2)

     CoeffSl(3,3) = -(lL11C*lL12*mSl2*                           &
       &    (kt11*mSl2*(2*lt12C*mbmcha2 + kt12C*(mb2 - mst2)) +  &
       &      lt11*(2*kt12C*mbmcha1*mSl2 +                       &
       &         lt12C*mcha1mcha2*(-2*mb2 + mSl2 + 2*mst2))))

     CoeffSl(3,4) = - lL11C*lL12*(lt11*                              &
       &       (2*kt12C*mbmcha1 + lt12C*mcha1mcha2) +               &
       &      kt11*(2*lt12C*mbmcha2 + kt12C*(mb2 - 2*mSl2 - mst2)) )

     CoeffSl(3,5) = - kt11*kt12C*lL11C*lL12

     smin = mSl2
     smax = (mstop(5) - mf_d(3))**2

     calc = .True.
     If (Check_Real_States) Then
      If (mStop(5).Gt.(mC(2)+mf_d(3))) Then  ! all intermediate states are real
       calc = .False.
      Else If (mStop(5).Gt.(mC(1)+mf_d(3))) Then ! only chi^+_2 contributes
       CoeffSl(1,:) = 0._dp
       CoeffSl(3,:) = 0._dp
      End If
     End If

     If (calc) Then
      Sup(5)%gi3(i_c) = f_stop * Dgauss(stoslep,smin,smax,prec)
      Sup(5)%id3(i_c,1) = Slept(2*(i1-1)+i2)%id + 1
      Sup(5)%id3(i_c,2) = id_nu(i1)
      Sup(5)%id3(i_c,3) = id_d(3)
      i_c = i_c +1
     End If
    End Do
   End Do
  End If

  !---------------------------
  ! stop_1 -> b + W + chi^0_j
  !---------------------------
  If (GenerationMixing) Then
  Else ! .not.GenerationMixing
   j = 1 ! lightest neutralino
   amN2 = mN2(j)

   If (mStop(5).Gt.(mW + mf_d(3) + Abs(mN(j)) ) ) Then ! kinematics
    !--------------------------------------------------
    ! sums that are needed for chargino-exchange terms
    !--------------------------------------------------
    sumM(1) = mf_d2(3) + mN2(j) + mW2 + mStop2(5)
    sumM(2) = sumM(1) - mC2(1)
    sumM(3) = sumM(1) - mC2(2)

    !--------------------------------------------------
    ! mass-combinations that are needed in the following
    !--------------------------------------------------
     mbchi = mf_d(3) * mN(j)
     mbt = mf_d(3) * mf_u(3)
     mb2 = mf_d2(3)
     mb4 = mf_d2(3)**2
     mchi2 = mN2(j)
     mchi4 = mchi2 * mchi2
     mchit = mN(j) * mf_u(3)
     mchp1chp2 = mC(1) * mC(2)
     mst2 = mStop2(5)
     mst4 = mst2 * mst2
     mt2 = mf_u2(3)

     mb2chi2 = mf_d2(3) * mchi2
     mb2st2 = mf_d2(3) * mst2
     mchi2st2 = mchi2 * mst2

     oomw = 1._dp / mW2
     mbomw = mf_d2(3) * oomw
     mbtomw = mbt * oomw
     mbchiomw = mbchi * oomw
     mchiomw = mchi2 * oomw
     mchitomw = mchit * oomw
     mstomw = mst2 * oomw
     !--------------------------------------------------
     ! coefficients for chargino-exchange terms
     !--------------------------------------------------
     ! 11,22
     !-------
     Do i=1,2
      qlqr = Real( c_CNW_L(i,j) * Conjg( c_CNW_R(i,j) ), dp )
      ql2qr2 = Abs( c_CNW_L(i,j) )**2 + Abs( c_CNW_R(i,j) )**2
      kl = Real(  c_CDSu_L(i,3,5) * Conjg( c_CDSu_R(i,3,5) ), dp )
      k2l2 = Abs( c_CDSu_L(i,3,5) )**2  + Abs( c_CDSu_R(i,3,5) )**2
      lqlkqr = Abs( c_CDSu_R(i,3,5) )**2 * Abs( c_CNW_R(i,j) )**2    &
           & + Abs( c_CDSu_L(i,3,5) )**2 * Abs( c_CNW_L(i,j) )**2
      kqllqr = Abs( c_CDSu_L(i,3,5) )**2 * Abs( c_CNW_R(i,j) )**2    &
           & + Abs( c_CDSu_R(i,3,5) )**2 * Abs( c_CNW_L(i,j) )**2
      kt1l1CQliCQri = Real( c_CNW_L(i,j) *  c_CDSu_L(i,3,5)             &
            &             * Conjg(  c_CNW_R(i,j) * c_CDSu_R(i,3,5)), dp )
      kt1Cl1QliCQri = Real( c_CNW_R(i,j) *  c_CDSu_L(i,3,5)             &
            &             * Conjg(  c_CNW_L(i,j) * c_CDSu_R(i,3,5)), dp )

      mbchp1 = mf_d(3) * mC(i)
      mchp1chi = mC(i) * mN(j)

      ca(i,1) = ( 6._dp*qlqr* k2l2 * mchp1chi * (2._dp * mb2 + mchi2 + mW2) &
            &   + 12._dp * kt1l1CQliCQri * mbchi * (mb2+mchi2+mst2+mW2)     &
            &   + 12._dp * kt1Cl1QliCQri * mbchi * mC2(i)                   &
            &   - 2._dp * kl * ql2qr2 * mbchp1                              &
            &           * ( 3._dp*mb2 + 2._dp*mchi2 + 3._dp*mst2            &
            &             + (mst2+mb2)**2 * oomw)                           &
            &  - lqlkqr * ( (4._dp*mb2+mW2) * (mchi2+mb2) + mchi4           &
            &             + mst2 * (2._dp*mW2 + 4._dp*mb2 + mchi2)          &
            &             + (mb2+mst2) * (mb4+mb2st2+mchi2st2) * oomw )     &
            &  - kqllqr * mC2(i) * ( mchi2 + mst2 + 2._dp * mb2             &
            &                      + (mb4+mb2st2)*oomw)  )

      ca(i,2) = - 12._dp * kt1l1CQliCQri * mbchi                            &
            & - 6._dp * k2l2 * qlqr * mchp1chi                              &
            & + 2._dp * kl* ql2qr2 * mbchp1 * (3._dp+2._dp*(mb2+mst2)*oomw) &
            & + lqlkqr * (6._dp * mb2 + 2._dp * mchi2 + 2._dp * mst2 + mW2  &
            &            + (3._dp * mb4 + mb2chi2 + 4.*mb2st2               &
            &            + 2._dp * mchi2st2 + mst4)*oomw )                  &
            & + kqllqr * mC2(i) * (2._dp + mb2 * oomw)

      ca(i,3) = - 2._dp * kl * ql2qr2 * mbchp1 * oomw         &
            & - lqlkqr * (2._dp + (mchi2+3._dp*mb2+2._dp*mst2) * oomw)

      ca(i,4) = lqlkqr * oomw

      ca(i,5) = -12._dp * kl * qlqr * mbchi                                 &
            & - 6._dp * k2l2 * qlqr * mchp1chi                              &
            & + 2._dp * kl * ql2qr2 * mbchp1 *(3._dp+2._dp*(mb2+mst2)*oomw) &
            & + lqlkqr * ( 6._dp*mb2 + 3._dp*mchi2 + 2._dp*mst2 + 2._dp*mW2 &
            &            + (2._dp*mb4 + 2._dp*mb2st2 + mchi2st2) * oomw)    &
            & + kqllqr * mC2(i) * (1._dp + (2._dp*mb2+mst2) * oomw)

      ca(i,6) = - 4.* kl * ql2qr2 * mbchp1 * oomw                          &
            & - lqlkqr * (4._dp + (mchi2 + 4.*mb2 + 2._dp*mst2 ) * oomw )  &
            & - kqllqr * mC2(i) * oomw

      ca(i,7) = 2._dp * lqlkqr * oomw

      ca(i,8) = -2._dp * kl * ql2qr2 * mbchp1 * oomw                    &
            & - lqlkqr * (2._dp + mb2*oomw) - kqllqr * mC2(i) * oomw

      ca(i,9) = lqlkqr * oomw
     End Do
     !-------
     ! 12
     !-------
     l1l2Ql2Qr1 = lt11 * lt12C * Conjg( c_CNW_R(2,j) ) * c_CNW_L(1,j) &
              & + kt11 * kt12C * c_CNW_R(1,j) * Conjg( c_CNW_L(2,j) )
     k1l2Ql2Qr1 = kt11 * lt12C * Conjg( c_CNW_R(2,j) ) *c_CNW_L(1,j)  &
              & + kt12C * lt11 * c_CNW_R(1,j) * Conjg( c_CNW_L(2,j) )
     k2l1Ql2Qr1 = kt12C * lt11 * Conjg( c_CNW_R(2,j) ) * c_CNW_L(1,j) &
              & + kt11 * lt12C * c_CNW_R(1,j) * Conjg( c_CNW_L(2,j) )
     k1k2Ql2Qr1 = kt11 * kt12C * Conjg( c_CNW_R(2,j) ) * c_CNW_L(1,j) &
              & + lt11 * lt12C * c_CNW_R(1,j) * Conjg( c_CNW_L(2,j) )
     l1l2Ql1Ql2 = lt11 * lt12C * c_CNW_R(1,j) * Conjg( c_CNW_R(2,j) ) &
              & + kt11 * kt12C * c_CNW_L(1,j) * Conjg( c_CNW_L(2,j) )
     k1l2Ql1Ql2 = kt11 * lt12C * c_CNW_R(1,j) * Conjg( c_CNW_R(2,j) ) &
              & + kt12C * lt11 * c_CNW_L(1,j) * Conjg( c_CNW_L(2,j) )
     k2l1Ql1Ql2 = kt12C * lt11 * c_CNW_R(1,j) * Conjg( c_CNW_R(2,j) ) &
              & + kt11 * lt12C * c_CNW_L(1,j) * Conjg( c_CNW_L(2,j) )
     k1k2Ql1Ql2 = kt11 * kt12C * c_CNW_R(1,j) * Conjg( c_CNW_R(2,j) ) &
              & + lt11 * lt12C * c_CNW_L(1,j) * Conjg( c_CNW_L(2,j) )

     mbchp1 = mf_d(3) * mC(1)
     mbchp2 = mf_d(3) * mC(2)
     mchp1chi = mC(1) * mN(j)
     mchp2chi = mC(2) * mN(j)

     ca(3,1) = 3._dp * l1l2Ql2Qr1 * mchp1chi * (2._dp*mb2+mchi2+mW2)   &
           & + 6._dp * k1l2Ql2Qr1 * mbchi * (mb2+mchi2+mst2+mW2)       &
           & + 6._dp * k2l1Ql2Qr1 * mbchi * mchp1chp2                  &
           & + 3._dp * k1k2Ql2Qr1 * mchp2chi * (2._dp*mb2+mchi2+mW2)   &
           & - l1l2Ql1Ql2 * ( (mchi2+mb2)*(mchi2+mst2+mW2+mb2)         &
           &                + mb2 * (3._dp*mb2+2._dp*mchi2+3._dp*mst2) &
           &                + 2._dp*mst2*mW2                           &
           &                + (mb2+mst2)*(mb4+mst2*(mchi2+mb2))*oomw ) &
           & - k1l2Ql1Ql2 * mbchp1 * (3._dp*mb2+2._dp*mchi2+3._dp*mst2 &
           &                         + (mb2+mst2)*(mb2+mst2)*oomw )    &
           & - k2l1Ql1Ql2 * mbchp2 * (3._dp*mb2+2._dp*mchi2+3._dp*mst2 &
           &                         + (mb2+mst2)*(mb2+mst2)*oomw )    &
           & - k1k2Ql1Ql2 * mchp1chp2 * (2._dp*mb2+mchi2+mst2          &
           &                            + (mb2+mst2)*mbomw )

     ca(3,2) = -3._dp * l1l2Ql2Qr1 * mchp1chi                            &
           & -6._dp * k1l2Ql2Qr1 * mbchi                                 &
           & -3._dp * k1k2Ql2Qr1 * mchp2chi                              &
           & + l1l2Ql1Ql2 * ( 6._dp*mb2 + 2._dp*mchi2 + 2._dp*mst2 + mW2 &
           &                + ((2._dp*mb2+mchi2)*(mb2+2._dp*mst2)        &
           &                   +mb4+mst4)*oomw )                         &
           & + k1l2Ql1Ql2 * mbchp1 * (3._dp+2._dp*(mb2+mst2)*oomw)       &
           & + k2l1Ql1Ql2 * mbchp2 * (3._dp+2._dp*(mb2+mst2)*oomw)       &
           & + k1k2Ql1Ql2 * mchp1chp2 * (2._dp+mbomw)

     ca(3,3) = -l1l2Ql1Ql2 * (2._dp+(mchi2+3._dp*mb2+2._dp*mst2)*oomw)  &
           & - ( k1l2Ql1Ql2 * mbchp1 + k2l1Ql1Ql2 * mbchp2 ) * oomw

     ca(3,4) = l1l2Ql1Ql2 * oomw

     ca(3,5) = - 3._dp * l1l2Ql2Qr1 * mchp1chi                            &
           & - 6._dp * k1l2Ql2Qr1 * mbchi                                 &
           & - 3._dp * k1k2Ql2Qr1 * mchp2chi                              &
           & + l1l2Ql1Ql2 * ( 6._dp*mb2+3._dp*mchi2+2._dp*mst2+2._dp*mW2  &
           &                + (2._dp*mb4+mst2*(mchi2+2._dp*mb2))*oomw )   &
           & + k1l2Ql1Ql2 * mbchp1 * (3._dp+2._dp*(mb2+mst2)*oomw )       &
           & + k2l1Ql1Ql2 * mbchp2 * (3._dp+2._dp*(mb2+mst2)*oomw )       &
           & + k1k2Ql1Ql2 * mchp1chp2 * (1._dp+(2._dp*mb2+mst2)*oomw )

     ca(3,6) = - l1l2Ql1Ql2 * (4.+(mchi2+4.*mb2+2._dp*mst2)*oomw)        &
           & - ( 2._dp*k1l2Ql1Ql2 * mbchp1 + 2._dp*k2l1Ql1Ql2 * mbchp2   &
           &   + k1k2Ql1Ql2 * mchp1chp2 ) * oomw

     ca(3,7) = 2._dp*l1l2Ql1Ql2 * oomw

     ca(3,8) = - l1l2Ql1Ql2 * (2._dp+mbomw)                    &
           & - ( k1l2Ql1Ql2 * mbchp1 + k2l1Ql1Ql2 * mbchp2     &
           &   + k1k2Ql1Ql2 * mchp1chp2 ) * oomw

     ca(3,9) = l1l2Ql1Ql2 * oomw

     ca(3,:) = 2._dp * ca(3,:)  ! factor because of interference

     !--------------------------------------------------
     ! coefficients for top-chargino_i terms
     !--------------------------------------------------
     Do i=1,2
      fltQl1 = Conjg( c_UNSu_R(3,j,5) ) * c_CDSu_R(i,3,5) * c_CNW_R(i,j)
      fktQl1 = Conjg( c_UNSu_R(3,j,5) ) * c_CDSu_L(i,3,5) * c_CNW_R(i,j)
      fltQr1 = Conjg( c_UNSu_R(3,j,5) ) * c_CDSu_R(i,3,5) * c_CNW_L(i,j)
      fktQr1 = Conjg( c_UNSu_R(3,j,5) ) * c_CDSu_L(i,3,5) * c_CNW_L(i,j)
      hltQl1 = Conjg( c_UNSu_L(3,j,5) ) * c_CDSu_R(i,3,5) * c_CNW_R(i,j)
      hktQl1 = Conjg( c_UNSu_L(3,j,5) ) * c_CDSu_L(i,3,5) * c_CNW_R(i,j)
      hltQr1 = Conjg( c_UNSu_L(3,j,5) ) * c_CDSu_R(i,3,5) * c_CNW_L(i,j)
      hktQr1 = Conjg( c_UNSu_L(3,j,5) ) * c_CDSu_L(i,3,5) * c_CNW_L(i,j)

      mbchp1 = mf_d(3) * mC(i)
      mchp1chi = mC(i) * mN(j)
      mchp1t = mC(i) * mf_u(3)

      cb(i,1) = fltQl1 * ( 2._dp * (mchi2-mst2) * (mW2+mchi2+mb2)         &
            &            - mb2 * mst2 * (1._dp+(mchi2-mb2-mst2)*oomw) )   &
            & + 3._dp * fktQl1 * mbchp1 * (mchi2-mst2)                    &
            & + 3._dp * hltQl1 * mchit * (2._dp*mb2+mchi2+mW2)            &
            & + 6._dp * hktQl1 * mchit * mbchp1                           &
            & + fltQr1 * mchp1chi * (2._dp*mW2-mb2*(1._dp+mbomw))         &
            & + fktQr1 * mbchi                                            &
            &    * (3._dp*mW2+mst2*(2._dp-mbomw)-mb2*(2._dp+mbomw)-mchi2) &
            & - hltQr1 * mchp1t                                           &
            &    * (mchi2+mst2*(1._dp+mbomw)+mb2*(2._dp+mbomw))           &
            & - hktQr1 * mbt                                              &
            &    * (2._dp*mchi2 + (mb2+mst2)*(3._dp+(mb2+mst2)*oomw))

      cb(i,2) = fltQl1 * (mst2-mchi2)*(2._dp-mbomw)               &
             & - 3._dp* hltQl1 * mchit                            &
             & + fktQr1 * mbchi * (mbomw-1._dp)                   &
             & + hltQr1 * mchp1t * (2._dp+mbomw)                  &
             & + hktQr1 * mbt * (3._dp+2._dp*(mb2+mst2)*oomw)

      cb(i,3) = - hktQr1 * mbt * oomw

      cb(i,4) = fltQl1 * (4.*mb2+2._dp*mW2+mchi2*mst2*oomw       &
            &              +mst2*(2._dp-mbomw) )                 &
            & + 3._dp * fktQl1 * mbchp1                          &
            & - 3._dp * hltQl1 * mchit                           &
            & - fltQr1 * mchp1chi * (1._dp-2._dp*mbomw)          &
            & + fktQr1 * mbchi * (1._dp+2._dp*mbomw+mst2*oomw)   &
            & + hltQr1 * mchp1t * (1._dp+2._dp*mbomw+mst2*oomw)  &
            & + hktQr1 * mbt * (3._dp+2._dp*(mb2+mst2)*oomw)

      cb(i,5) = - fltQl1 * (3._dp+mbomw+(mst2+mchi2)*oomw)    &
            & - ( fktQr1 * mbchi  + hltQr1 * mchp1t           &
            &   + 2._dp * hktQr1 * mbt ) * oomw

      cb(i,6) = fltQl1 * oomw

      cb(i,7) = - 2._dp*fltQl1                                &
            & - ( fltQr1 * mchp1chi  + fktQr1 * mbchi         &
            &   + hltQr1 * mchp1t  + hktQr1 * mbt ) * oomw

      cb(i,8) = fltQl1 * oomw
     End Do

     cb = - sqrt2 * g * cb  ! factor because of interference

     !--------------------------------------------------
     ! coefficients for sbottom_k-chargino_i terms
     !--------------------------------------------------
     Do i=1,2
      mbchp1 = mf_d(3) * mC(i)
      mchp1chi = mC(i) * mN(j)
      mbchp1omw = mbchp1 * oomw
      mchp1chiomw = mchp1chi * oomw

      Do k=1,2
       hbltQl1 = Conjg( c_DNSd_L(3,j,k+4) ) * c_CDSu_R(i,3,5) * c_CNW_R(i,j) &
             & + Conjg( c_DNSd_R(3,j,k+4) ) * c_CDSu_L(i,3,5) * c_CNW_L(i,j)
       fbltQl1 = Conjg( c_DNSd_R(3,j,k+4) ) * c_CDSu_R(i,3,5) * c_CNW_R(i,j) &
             & + Conjg( c_DNSd_L(3,j,k+4) ) * c_CDSu_L(i,3,5) * c_CNW_L(i,j)
       hbktQl1 = Conjg( c_DNSd_L(3,j,k+4) ) * c_CDSu_L(i,3,5) * c_CNW_R(i,j) &
             & + Conjg( c_DNSd_R(3,j,k+4) ) * c_CDSu_R(i,3,5) * c_CNW_L(i,j)
       fbktQl1 = Conjg( c_DNSd_R(3,j,k+4) ) * c_CDSu_L(i,3,5) * c_CNW_R(i,j) &
             & + Conjg( c_DNSd_L(3,j,k+4) ) * c_CDSu_R(i,3,5) * c_CNW_L(i,j)

       cc(k,i,1) = hbltQl1 *mbchi *(mb2-mst2-2._dp*mchi2+mst2*(mbomw+mstomw)) &
              & + fbltQl1 * ( mb2*mW2+mchi2*mst2*(mstomw-1._dp)               &
              &             + mb2*mst2*(mbomw-2._dp)+mb4-2._dp*mb2*mchi2      &
              &             + mb2*mst2*mstomw )                               &
              & + hbktQl1 *mchp1chi *(mb2-mW2-2._dp*mchi2+mst2*(1._dp+mbomw)) &
              & + fbktQl1 * mbchp1 *(mb2-2._dp*mchi2+mst2*(mbomw+mstomw-1._dp))

       cc(k,i,2) = - hbltQl1 * mbchi * (1._dp+mbomw+2._dp*mstomw)            &
              & - fbltQl1 * (mchi2*(1._dp+2._dp*mstomw)+mst2*mstomw          &
              &             +mW2+mb2*(3._dp+mbomw+3._dp*mstomw) )            &
              & - hbktQl1 * mchp1chi * (mbomw-1._dp)                         &
              & - fbktQl1 * mbchp1 * (1._dp+mbomw+2._dp*mstomw)

       cc(k,i,3) = hbltQl1 * mbchiomw + fbktQl1 * mbchp1omw  &
              & + fbltQl1 * (2._dp+2._dp*mbomw+mchiomw+2._dp*mstomw)

       cc(k,i,4) = - fbltQl1 * oomw

       cc(k,i,5) = ( hbktQl1 * mchp1chi + fbktQl1 * mbchp1            &
               &   + hbltQl1 * mbchi + fbltQl1 * mb2 ) * (1._dp-mstomw)

       cc(k,i,6) = hbltQl1 * mbchiomw + fbltQl1 * (1._dp+mbomw+mstomw) &
               & + hbktQl1 * mchp1chiomw + fbktQl1 * mbchp1omw

       cc(k,i,7) = cc(k,i,4)

       cc(k,i,:) = cc(k,i,:) * Conjg( C_SdSuW(k+4,5) )

      End Do
     End Do
     cc = -2._dp * cc

     !--------------------------------------------------
     ! coefficients for top-top term
     !--------------------------------------------------
     ft2 = Abs( c_UNSu_R(3,j,5) )**2
     ht2 = Abs( c_UNSu_L(3,j,5) )**2
     fht = Real( c_UNSu_R(3,j,5) * Conjg( c_UNSu_L(3,j,5) ),dp )

     cd(1,1) = ft2 * (mchi2-mst2) * (2._dp*mW2-mb2*(1._dp+mbomw))    &
           & - ht2 * mt2 * (mchi2+mst2+mb2*(2._dp+mbomw+mst2*oomw))  &
           & + 2._dp * fht * mchit * (2._dp*mW2-mb2*(1._dp+mbomw))

     cd(1,2) = ht2 * mt2 * (2._dp+mbomw)

     cd(1,3) = ft2 * (mb2+2._dp*mst2+2._dp*mW2+(2._dp*mchi2-mst2)*mbomw)  &
           & + ht2 * mt2 *(1._dp+(2._dp*mb2+mst2)*oomw)                   &
           & - 2._dp * fht * mchit * (1._dp-2._dp*mbomw)

     cd(1,4) = - ft2 * (2._dp+mbomw) - ht2 * mt2 * oomw

     cd(1,5) = - ft2 * (2._dp+mchi2*oomw) &
           & - ( 2._dp*fht * mchit  + ht2 * mt2 ) * oomw

     cd(1,6) = ft2 * oomw

     cd(1,:) = 0.5_dp *g**2 * cd(1,:)

     !--------------------------------------------------
     ! coefficients for top-sbottom_k terms
     !--------------------------------------------------
     Do i=1,2
      ffb1 = Conjg( c_DNSd_R(3,j,i+4) ) * c_UNSu_R(3,j,5)
      fhb1 = Conjg( c_DNSd_R(3,j,i+4) ) * c_UNSu_L(3,j,5)
      hfb1 = Conjg( c_DNSd_L(3,j,i+4) ) * c_UNSu_R(3,j,5)
      hhb1 = Conjg( c_DNSd_L(3,j,i+4) ) * c_UNSu_L(3,j,5)

      cd(1+i,1) = ffb1 * (mb2*(mchi2+mst2)+2._dp*mchi2*(mst2-mchi2-mW2)  &
              &          +(mchi2-mst2)*mb2*mstomw )                      &
              & + (fhb1 * mchit + hfb1 * mbchi)                          &
              &              * (mb2+mst2-2._dp*mchi2-mW2+mb2*mstomw)     &
              & + hhb1 * mbt * (mb2-2._dp*mchi2-mst2+(mb2+mst2)*mstomw)

      cd(1+i,2) = ffb1 * (mchi2-mst2)*(2._dp-mbomw)                     &
              & + (fhb1 * mchit + hfb1 * mbchi) * (1._dp-mbomw)         &
              & - hhb1 * mbt * (1._dp+mbomw+2._dp*mstomw)

      cd(1+i,3) = hhb1 * mbt * oomw

      cd(1+i,4) = ( ffb1*mchi2 + fhb1 * mchit + hfb1 * mbchi  &
              & + hhb1 * mbt) * (1._dp-mstomw)

      cd(1+i,5) = ffb1 * (1._dp+mchiomw+mstomw)                     &
              & + fhb1 * mchitomw + hfb1 * mbchiomw + hhb1 * mbtomw

      cd(1+i,6) = -ffb1 * oomw

      cd(1+i,:) =  sqrt2 * g * Conjg( c_SdSuW(i+4,5) ) * cd(1+i,:)

     End Do

     !--------------------------------------------------
     ! coefficients for sbottom_i-sbottom_k terms
     !--------------------------------------------------
     ! i=k :
     !--------
     Do i=1,2
      ft2 = Abs( c_SdSuW(i+4,5) * c_DNSd_R(3,j,i+4) )**2
      ht2 = Abs( c_SdSuW(i+4,5) * c_DNSd_L(3,j,i+4) )**2
      fht = Real( c_DNSd_R(3,j,i+4) * Conjg( c_DNSd_L(3,j,i+4) ),dp ) &
          & * Abs( c_SdSuW(i+4,5) )**2

      ce(i,1) = -( (ft2 +ht2 ) * (mb2+mchi2) + 4._dp * fht * mbchi )

      ce(i,2) = (ft2 +ht2 )
     End Do
     !-----------
     ! i=1, k=2 :
     !-----------
     ce(3,2) = 2._dp * ( c_DNSd_R(3,j,5) * Conjg( c_DNSd_R(3,j,6) )    &
           &     + c_DNSd_L(3,j,5) * Conjg( c_DNSd_L(3,j,6) ) )  &
           &    * c_SdSuW(5,5) * Conjg( c_SdSuW(6,5) )

     ce(3,1) = - ce(3,2) * (mb2+mchi2)                           &
           & - 4._dp * ( c_DNSd_R(3,j,5) * Conjg( c_DNSd_L(3,j,6) )    &
           &           + c_DNSd_L(3,j,5) * Conjg( c_DNSd_R(3,j,6) ) )  &
           &         * c_SdSuW(5,5) * Conjg( c_SdSuW(6,5) ) * mbchi

     smin = (mf_d(3) + mN(j))**2
     smax = (mstop(5) - mW)**2

     calc = .True.
     If (Check_Real_States) Then
      mSb = Sqrt(MSbot2)

      If (mStop(5).Gt.Max(mC(2)+mf_d(3),mf_u(3)+mN(j),mW+mSb(2))) Then
       calc = .False.   ! all intermediate states are real
      Else ! check now for intermediate real states
       If (mStop(5).Gt.(mC(2)+mf_d(3))) Then ! no chargino contribution
        ca = 0._dp
        cb = 0._dp
        cc = 0._dp
       Else If (mStop(5).Gt.(mC(1)+mf_d(3))) Then ! only chi^+_2 contributes
        ca(1,:) = 0._dp
        ca(3,:) = 0._dp
        cb(1,:) = 0._dp
        cc(1,:,:) = 0._dp
       End If

       If (mStop(5).Gt.(mN(j)+mf_u(3))) Then ! no top contribution
        cb = 0._dp
        cd = 0._dp
       End If

       If (mStop(5).Gt.(mSb(2)+mW)) Then ! no sbottom contribution
        cc = 0._dp
        cd(2:3,:) = 0._dp
        ce = 0._dp

       Else If (mStop(5).Gt.(mSb(1)+mW)) Then ! only sbottom_2 contributes
        cc(:,1,:) = 0._dp
        cd(2,:) = 0._dp
        ce(1,:) = 0._dp
        ce(3,:) = 0._dp
       End If

      End If
     End If

    If (calc)  Then
      Sup(5)%gi3(i_c) = f_stop * dgauss(stbWchi, smin, smax, prec)
      Sup(5)%id3(i_c,1) = Chi0(1)%id
      Sup(5)%id3(i_c,2) = id_W(1)
      Sup(5)%id3(i_c,3) = id_d(3)
      i_c = i_c +1
    End If
   End If ! kinematic check
  End If
  Sup(5)%g = Sum(Sup(5)%gi2) + Sum(Sup(5)%gi3)
  If (Sup(5)%g.Ne.0._dp) Then
   Sup(5)%bi2 = Sup(5)%gi2 / Sup(5)%g
   Sup(5)%bi3 = Sup(5)%gi3 / Sup(5)%g
  End If

  Iname = Iname - 1

 End Subroutine StopDecays3

 Real(dp) Function  stoslep(s)
 !----------------------------------------------------------------------
 ! Integrand for the decay stop_1 -> b + slepton + lepton'.
 ! written by Werner Porod
 ! 24.01.02:portation to f90
 !----------------------------------------------------------------------
 Implicit None
  Real(dp), Intent(in) :: s

  Integer :: i
  Real(dp) :: wsr, s2
  Complex(dp) :: inte(3)
 !-----------------------------------------------------
 ! the following variables from the commons are needed:
 !      integer istosle
 !  -coeffSl,istosle
 !   (coffecients for the integrals in the function stoslep) and
 !   istosle controls the kind of decay:
 !     istosle=1 : stop_1 -> b + sneutrino + electron
 !     istosle=2 : stop_1 -> b + tau-sneutrino + tau
 !     istosle=3 : stop_1 -> b + selectron_L + neutrino
 !     istosle=4 : stop_1 -> b + stau_1 + tau-neutrino
 !     istosle=5 : stop_1 -> b + stau_2 + tau-neutrino
 !  -amc2(i,j), i,j=1,2
 !  -mb2
 !  -amtau2
 !  -msne2,mtsneut2
 !  -mStop2(5)
 !-----------------------------------------------------

  Iname = Iname + 1
  NameOfUnit(Iname) = "StoSle"

  s2 = s*s

  If (istosle.Eq.1) Then ! sneutrino + lepton
   wsr = kappa(s,mst2,mb2) * kappa(s, mSl2, mL2)

  Elseif (istosle.Eq.2) Then
   wsr = kappa(s,mst2,mb2)

  Else
   Write (ErrCan,*) 'Error in function stoslep: istolse =',istosle
   Call TerminateProgram()
  End If

  inte(1) = wsr / ( (s-amc2(1))**2 + gammc2(1) )
  inte(2) = wsr / ( (s-amc2(2))**2 + gammc2(2) )
!  inte(3) = wsr * ( (s-amc2(1))*(s-amc2(2)) + gammc(1)*gammc(2) )         &
!        &  / ( ((s-amc2(1))**2 + gammc2(1)) *((s-amc2(2))**2 + gammc2(2)) )
  inte(3) = wsr / ( (s-amc2(1) + (0._dp,1._dp) * gammc(1) )         &
          &       * (s-amc2(2) - (0._dp,1._dp) * gammc(2) ) )

  stoslep = 0.
  If (istosle.Eq.1) Then
   Do i=1,3
    stoslep = stoslep                                                     &
          &  + Real( inte(i) * ( coeffSl(i,1) + coeffSl(i,2)/s2           &
          &                    + coeffSl(i,3)/s + coeffSl(i,4)*s),dp )
   Enddo
  Else
   Do i=1,3
    stoslep = stoslep                                                       &
          &  + Real( inte(i) * ( coeffSl(i,1) + coeffSl(i,2)/s2             &
          &        + coeffSl(i,3)/s + coeffSl(i,4) * s + coeffSl(i,5)*s2),dp )
   Enddo
  End If
  If (stoslep.Lt.0) Then
   Write (ErrCan,*) 'istosle,s,stoslep'
   Write (ErrCan,*) istosle,s,stoslep
   Write (ErrCan,*) 'inte(1),inte(2),inte(3)'
   Write (ErrCan,*) inte(1),inte(2),inte(3)
   Write (ErrCan,*) 'coeffSl(i,j)'
   Write (ErrCan,*) coeffSl(1,:)
   Write (ErrCan,*) coeffSl(2,:)
   Write (ErrCan,*) coeffSl(2,:)
  End If

  Iname = Iname - 1

 End Function  stoslep

End Module Stop3BodyDecays


text: David Dionisi is the Noteworthy Published Authors author of http://www.americanhiroshima.info/ American Hiroshima and a resident of Davis. He is the founder of Teach Peace.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)
      DOUBLE PRECISION A(*),LAMBDA,X(*),V(*),PAR(*)
      INTEGER IPAR(*),K
C
C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHO, RHOJAC.
C
C RETURN IN THE VECTOR  V  THE KTH COLUMN OF THE JACOBIAN
C MATRIX [D RHO/D LAMBDA, D RHO/DX] EVALUATED AT THE POINT
C (A, LAMBDA, X).
C
C THE FOLLOWING CODE IS SPECIFICALLY FOR THE POLYNOMIAL SYSTEM DRIVER
C  POLSYS , AND SHOULD BE USED VERBATUM WITH  POLSYS .  IF THE USER IS
C CALLING  FIXP??  OR   STEP??  DIRECTLY, HE MUST SUPPLY APPROPRIATE
C REPLACEMENT CODE HERE.
      INTEGER J,N,N2
      N=IPAR(1)
      N2=2*N
      IF (K .EQ. 1) THEN
C FORCE PREDICTED POINT TO HAVE  LAMBDA .GE. 0  .
        IF (LAMBDA .LT. 0.0) LAMBDA=0.0
        CALL HFUNP(A,LAMBDA,X,PAR,IPAR)
        DO 10 J=1,N2
          V(J)=PAR(IPAR(3 + (6-1)) + (J-1))
   10   CONTINUE
        RETURN
      ELSE
        DO 20 J=1,N2
          V(J)=PAR(IPAR(3 + (5-1)) + (J-1) + N2*(K-2))
   20   CONTINUE
      ENDIF
C
      RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
program main
  use my_mod
  real(real64) :: vf(2,25)
  vf = 1.0_real64
  call use_count_with_dim (vf)
end program

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
C
C	$Id: wmlgpl.f,v 1.7 2008-07-27 00:17:37 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE WMLGPL(N,X,Y)
C
C  Draws polylines either using GPL (IWDTYP=0) or WMDRFL (IWDTYP=1).
C
      DIMENSION X(N),Y(N)
C
      include 'wmcomn.h'
C
C  Save the current line attributes and reset.
C
      CALL GQPLCI(IER,ICOLD)
      CALL GSPLCI(ICOLOR)
C
C  Handle the squall lines (IFRONT = 5), tropical fronts (IFRONT = 6),
C  and convergence lines (IFRONT = 7) separately.
C
      IF (IFRONT .EQ. 5) THEN
        CALL GQLWSC(IER,OSCL)
        CALL GSLWSC(DLINWD)
        CALL DPGETI('DPS',IDPS)
        CALL DPGETR('WOG',OWOG)
        CALL DPGETR('WOS',OWOS)
        CALL DPGETR('LS1',OLS1)
C
        CALL DPSETR('WOG',.003)
        CALL DPSETR('WOS',.0035)
        CALL DPSETR('LS1',0.0)
        CALL DPSETC('DPT','___:F34X150Y150:7_:F34X150Y150:7___$$$$$$$$$$
     +$$$$$$')
        CALL DPCURV(X,Y,N)
C
        CALL GSLWSC(OSCL)
        CALL DPSETI('DPS',IDPS)
        CALL DPSETC('DPT',' ')
        CALL DPSETR('LS1',OLS1)
        CALL DPSETR('WOG',OWOG)
        CALL DPSETR('WOS',OWOS)
        RETURN
      ELSE IF (IFRONT .EQ. 6) THEN
        CALL GQLWSC(IER,OSCL)
        CALL GSLWSC(DLINWD)
        CALL DPGETI('DPS',IDPS)
        CALL DPGETI('DPT',IDPT)
        CALL DPGETR('WOG',OWOG)
        CALL DPGETR('WOS',OWOS)
C
        CALL DPSETI('DPS',-22)
        CALL DPSETR('WOG',.007)
        CALL DPSETR('WOS',.007)
        CALL DPSETI('DPT',4177920)
        CALL GSPLCI(ITRO1C)
        CALL DPCURV(X,Y,N)
        CALL DPSETI('DPT',2040)
        CALL GSPLCI(ITRO2C)
        CALL DPCURV(X,Y,N)
C
        CALL GSLWSC(OSCL)
        CALL DPSETI('DPS',IDPS)
        CALL DPSETI('DPT',IDPT)
        CALL DPSETR('WOG',OWOG)
        CALL DPSETR('WOS',OWOS)
        RETURN
      ELSE IF (IFRONT .EQ. 7) THEN
        CALL GQLWSC(IER,OSCL)
        CALL GSLWSC(DLINWD)
        CALL DPGETI('DPS',IDPS)
        CALL DPGETI('DPT',IDPT)
        CALL DPGETR('WOG',OWOG)
        CALL DPGETR('WOS',OWOS)
C
        CALL DPSETI('DPS',-5)
        CALL DPSETR('WOG',.0075)
        CALL DPSETR('WOS',.0075)
        CALL DPSETI('DPT',31)
        CALL DPCURV(X,Y,N)
C
        CALL GSLWSC(OSCL)
        CALL DPSETI('DPS',IDPS)
        CALL DPSETI('DPT',IDPT)
        CALL DPSETR('WOG',OWOG)
        CALL DPSETR('WOS',OWOS)
        RETURN
      ENDIF
C
      IF (IWDTYP .EQ. 1) THEN
        CALL WMDRFL(N,X,Y)
      ELSE
        CALL GPL(N,X,Y)
      ENDIF
C
C  Restore line attributes.
C
      CALL GSPLCI(ICOLD)
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
MODULE constants
!
!   Define some constants
!
  INTEGER, PARAMETER       :: db     = SELECTED_REAL_KIND(8)
!
  DOUBLE PRECISION, PARAMETER :: vlight = 299792458.0_db         ! c
  DOUBLE PRECISION, PARAMETER :: vacimp = 376.73031346177066_db  ! \mu_0*c
  DOUBLE PRECISION, PARAMETER :: eev    = 510998.89613320108_db  ! m_e*c^2/e
  DOUBLE PRECISION, PARAMETER :: pi     = 3.1415926535897931_db
END MODULE constants

max_stars_repo_name: alexbovet/fastpart

max_issues_repo_name: alexbovet/fastpart

max_forks_repo_name: alexbovet/fastpart


text:
Yolo Compassionate Friends is a local chapter of The Compassionate Friends.  It is made up of parents who have children that have died.  It also includes support for siblings and grandparents.

The chapter is currently seeking new leadership, and is inactive.  It used to  meet monthly at the conference room of Yolo Hospice on the third Wednesday of each month.

At Compassionate Friends meetings, people have a chance to share experiences, exchange stories, and provide support to other.  Sometimes they have a speaker; other times they discuss a common issue.  There are events shared with thousands other people around the nation and world including their annual candle lighting in December. They do not have outside professional leadership.  Each of them has a common experience.  They are available to provide support to people who have immediate needs.

Contact information for nearby TCF chapters is listed above.  Additional information about other grief resources and The Compassionate Friends national organization may be found on the Yolo Compassionate Friends website.


text:
      SUBROUTINE PLT_CNTR_S
     &          (mxlat,nrlat,mxlon,nrlon,
     &           amp1,std1,
     &           icomp,iplot,ntlist,ntl,
     &           ampt,stdt)

      dimension     amp1(mxlon,mxlat    ),std1(mxlon,mxlat    ),
     &              ampt(mxlon,mxlat,2,5),stdt(mxlon,mxlat,2,5),
     &              ntlist(5)

      ntl=1
      nt=iplot
      ntlist(1)=iplot

      do l=1,nrlat
         do k=1,nrlon
            amp1(k,l)=ampt(k,l,icomp,nt)
            std1(k,l)=stdt(k,l,icomp,nt)
         end do
      end do

      RETURN
      END      ! PLT_CNTR_S

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
Thille Hall was named for Albert J. Thille and was the collective name for Tercero buildings A through F. This group made up the northwestern portion of the Tercero area, bordered by Dairy Road and Tercero Hill Circle, until it was demolished in April 2012.

The university has been planning the expansion of Tercero for quite some time. First there were the Tercero South Residence Halls and the dining commons renovation. Then Tercero South Phase II. The plan after this was to demolish Leach Hall and build additional dorms there. However, in January of 2011, a seismic evaluation found the Pierce Hall buildings seismically deficient. The concrete structures employ steel connections to connect floors to supporting walls, which was consistent with thencurrent building codes but do not comply with current seismic standards. Under University of California seismic policy, the buildings were determined to be poor to very poor. http://www.universityofcalifornia.edu/regents/regmeet/jan11/gb2b.pdf (Source) Since Leach Hall was safe for the time being, student housing decided to demolish Pierce and Thille Halls instead and expand there. The buildings were demolished in early 2012. Construction on the new buildings is scheduled to begin in September 2012 and is expected to be completed by May 2014.

A note of interest: While physically Thille Hall was buildings AF, the Thille Resident Advisors RA Staff encompassed only buildings AE. F Building RAs were included in the Pierce RA Staff.

Gardner Action Dailey designed this complex, along with Pierce Hall. It was built in 1967.

 Buildings
  Tercero A Building A Building  Better known as Casa (short for Casa Cuauhtémoc), this building offers a crosscultural opportunity with emphasis on the exploration of Chicano/Latino thought and culture. (from http://www.housing.ucdavis.edu/Halls/New/programnew.htm#casa)
  Tercero B Building B Building  Literally the closest dorm to the dairy farm on campus.  Check out a campus map for yourself.
  Tercero C Building C Building
  Tercero D Building D Building  Hammarskjöld House (aka Hamm House) is a International Studies/Community Service themed program building. It houses many International Relations majors, though one doesnt have to have such a major to live there.
  Tercero E Building E Building
  Tercero F Building F Building  Home of the Communications ATP, though once again, one doesnt have to be a Communications major to live there.


text:
*     MB03QG EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 10 )
      INTEGER          LDA, LDE, LDU, LDV
      PARAMETER        ( LDA = NMAX, LDE = NMAX, LDU = NMAX, LDV = NMAX)
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = 8*NMAX + 16 )
*     .. Local Scalars ..
      CHARACTER*1      DICO, JOBU, JOBV, STDOM
      INTEGER          I, INFO, J, N, NDIM, NLOW, NSUP
      DOUBLE PRECISION ALPHA
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), BETA(NMAX), DWORK(LDWORK),
     $                 E(LDE,NMAX), U(LDU,NMAX), V(LDV,NMAX), WI(NMAX),
     $                 WR(NMAX)
      LOGICAL          BWORK(NMAX)
*     .. External Functions ..
      LOGICAL          DELCTG
*     .. External Subroutines ..
      EXTERNAL         DGGES, MB03QG
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, NLOW, NSUP, ALPHA, DICO, STDOM, JOBU,
     $                      JOBV
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99990 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
*        Compute Schur form, eigenvalues and Schur vectors.
         CALL DGGES( 'Vectors', 'Vectors', 'Not sorted', DELCTG, N,
     $               A, LDA, E, LDE, NDIM, WR, WI, BETA, U, LDU, V, LDV,
     $               DWORK, LDWORK, BWORK, INFO )
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
*           Block reordering.
            CALL MB03QG( DICO, STDOM, JOBU, JOBV, N, NLOW, NSUP, ALPHA,
     $                   A, LDA, E, LDE, U, LDU, V, LDV, NDIM, DWORK,
     $                   LDWORK, INFO )
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99997 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99996 ) NDIM
               WRITE ( NOUT, FMT = 99994 )
               DO 10 I = 1, N
                  WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
   10          CONTINUE
               WRITE ( NOUT, FMT = 99993 )
               DO 20 I = 1, N
                  WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N )
   20          CONTINUE
               WRITE ( NOUT, FMT = 99992 )
               DO 30 I = 1, N
                  WRITE ( NOUT, FMT = 99995 ) ( U(I,J), J = 1,N )
   30          CONTINUE
               WRITE ( NOUT, FMT = 99991 )
               DO 40 I = 1, N
                  WRITE ( NOUT, FMT = 99995 ) ( V(I,J), J = 1,N )
   40          CONTINUE
           END IF
         END IF
      END IF
*
      STOP
*
99999 FORMAT (' MB03QG EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from DGEES  = ',I2)
99997 FORMAT (' INFO on exit from MB03QG = ',I2)
99996 FORMAT (' The number of eigenvalues in the domain is ',I5)
99995 FORMAT (8X,20(1X,F8.4))
99994 FORMAT (/' The ordered Schur form matrix is ')
99993 FORMAT (/' The ordered triangular matrix is ')
99992 FORMAT (/' The transformation matrix U is ')
99991 FORMAT (/' The transformation matrix V is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine hyperelastics(inpc,textpart,elcon,nelcon,
     &  nmat,ntmat_,ncmat_,irstrt,istep,istat,n,iperturb,iline,ipol,
     &  inl,ipoinp,inp,ipoinpc)
!
!     reading the input deck: *HYPERELASTIC
!
      implicit none
!
      character*1 inpc(*)
      character*132 textpart(16)
!
      integer nelcon(2,*),nmat,ntmat,ntmat_,istep,istat,ipoinpc(0:*),
     &  n,key,i,j,k,ityp,iperturb(*),iend,jcoef(3,14),ncmat_,irstrt,
     &  iline,ipol,inl,ipoinp(2,*),inp(3,*)
!
      real*8 elcon(0:ncmat_,ntmat_,*),um
!
!     jcoef indicates for each hyperelastic model the position of
!     the compressibility coefficients in the field elcon (max. 3
!     positions per model)
!
      data jcoef /3,0,0,3,0,0,2,0,0,3,0,0,5,6,0,7,8,9,3,0,0,
     &            6,7,0,12,13,14,2,0,0,3,4,0,4,5,6,5,0,0,4,5,6/
!
      ntmat=0
      iperturb(1)=3
      iperturb(2)=1
      write(*,*) '*INFO reading *HYPERELASTIC: nonlinear geometric'
      write(*,*) '      effects are turned on'
      write(*,*)
!
      if((istep.gt.0).and.(irstrt.ge.0)) then
         write(*,*)
     &       '*ERROR reading *HYPERELASTIC: *HYPERELASTIC should be'
         write(*,*) '  placed before all step definitions'
         call exit(201)
      endif
!
      if(nmat.eq.0) then
         write(*,*)
     &        '*ERROR reading *HYPERELASTIC: *HYPERELASTIC should be'
         write(*,*) '  preceded by a *MATERIAL card'
         call exit(201)
      endif
!
      ityp=-7
!
      do i=2,n
         if(textpart(i)(1:12).eq.'ARRUDA-BOYCE') then
            ityp=-1
         elseif(textpart(i)(1:13).eq.'MOONEY-RIVLIN') then
            ityp=-2
         elseif(textpart(i)(1:8).eq.'NEOHOOKE') then
            ityp=-3
         elseif(textpart(i)(1:5).eq.'OGDEN') then
            ityp=-4
         elseif(textpart(i)(1:10).eq.'POLYNOMIAL') then
            ityp=-7
         elseif(textpart(i)(1:17).eq.'REDUCEDPOLYNOMIAL') then
            ityp=-10
         elseif(textpart(i)(1:11).eq.'VANDERWAALS') then
            ityp=-13
         elseif(textpart(i)(1:4).eq.'YEOH') then
            ityp=-14
         elseif(textpart(i)(1:2).eq.'N=') then
            if(textpart(i)(3:3).eq.'1') then
            elseif(textpart(i)(3:3).eq.'2') then
               if(ityp.eq.-4) then
                  ityp=-5
               elseif(ityp.eq.-7) then
                  ityp=-8
               elseif(ityp.eq.-10) then
                  ityp=-11
               else
                  write(*,*) '*WARNING reading *HYPERELASTIC: N=2 is not
     & applicable for this material type; '
                  call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
               endif
            elseif(textpart(i)(3:3).eq.'3') then
               if(ityp.eq.-4) then
                  ityp=-6
               elseif(ityp.eq.-7) then
                  ityp=-9
               elseif(ityp.eq.-10) then
                  ityp=-12
               else
                  write(*,*) '*WARNING reading *HYPERELASTIC: N=3 is not
     & applicable for this material type; '
                  call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
               endif
            else
               write(*,*) '*WARNING reading *HYPERELASTIC: only N=1, N=2
     &, or N=3 are allowed; '
               call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
            endif
         else
            write(*,*)
     &       '*WARNING reading *HYPERELASTIC: parameter not recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
         endif
      enddo
!
      nelcon(1,nmat)=ityp
!
      if((ityp.ne.-6).and.(ityp.ne.-9)) then
         if((ityp.eq.-3).or.(ityp.eq.-10)) then
            iend=2
         elseif((ityp.eq.-1).or.(ityp.eq.-2).or.(ityp.eq.-4).or.
     &      (ityp.eq.-7)) then
            iend=3
         elseif(ityp.eq.-11) then
            iend=4
         elseif(ityp.eq.-13) then
            iend=5
         elseif((ityp.eq.-5).or.(ityp.eq.-12).or.(ityp.eq.-14)) then
            iend=6
         elseif(ityp.eq.-8) then
            iend=7
         endif
         do
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) exit
            ntmat=ntmat+1
            nelcon(2,nmat)=ntmat
            if(ntmat.gt.ntmat_) then
               write(*,*)'*ERROR reading *HYPERELASTIC: increase ntmat_'
               call exit(201)
            endif
            do i=1,iend
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &                  elcon(i,ntmat,nmat)
               if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
            enddo
            read(textpart(iend+1)(1:20),'(f20.0)',iostat=istat)
     &                  elcon(0,ntmat,nmat)
            if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
         enddo
      else
         do
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) exit
            ntmat=ntmat+1
            nelcon(2,nmat)=ntmat
            if(ntmat.gt.ntmat_) then
               write(*,*)'*ERROR reading *HYPERELASTIC: increase ntmat_'
               call exit(201)
            endif
            do i=1,8
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &                     elcon(i,ntmat,nmat)
               if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
            enddo
!
            if(ityp.eq.-6) then
               iend=1
            else
               iend=4
            endif
            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &           ipoinp,inp,ipoinpc)
            if((istat.lt.0).or.(key.eq.1)) then
               write(*,*)
     &           '*ERROR reading *HYPERELASTIC: hyperelastic definition'
               write(*,*) '  is not complete. '
               call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
               call exit(201)
            endif
            do i=1,iend
               read(textpart(i)(1:20),'(f20.0)',iostat=istat)
     &                  elcon(8+i,ntmat,nmat)
               if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
            enddo
            read(textpart(iend+1)(1:20),'(f20.0)',iostat=istat)
     &                  elcon(0,ntmat,nmat)
            if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
     &"*HYPERELASTIC%")
         enddo
      endif
!
!     if any of the compressibility coefficients is zero (incompressible
!     material), it is replaced. The lowest order coefficient is replaced
!     such that it corresponds to a Poisson coeffient of 0.475, the
!     following ones are replaced by a power of the first one
!
      do j=1,ntmat
!
!        calculating the shear coefficient in the undeformed state
!
         if(ityp.eq.-1) then
            um=elcon(1,j,nmat)
         elseif(ityp.eq.-2) then
            um=2.d0*(elcon(1,j,nmat)+elcon(2,j,nmat))
         elseif(ityp.eq.-3) then
            um=2.d0*elcon(1,j,nmat)
         elseif(ityp.eq.-4) then
            um=elcon(1,j,nmat)
         elseif(ityp.eq.-5) then
            um=elcon(1,j,nmat)+elcon(3,j,nmat)
         elseif(ityp.eq.-6) then
            um=elcon(1,j,nmat)+elcon(3,j,nmat)+elcon(5,j,nmat)
         elseif((ityp.eq.-7).or.(ityp.eq.-8).or.(ityp.eq.-9)) then
            um=2.d0*(elcon(1,j,nmat)+elcon(2,j,nmat))
         elseif((ityp.eq.-10).or.(ityp.eq.-11).or.(ityp.eq.-12)) then
            um=2.d0*elcon(1,j,nmat)
         elseif(ityp.eq.-13) then
            um=elcon(1,j,nmat)
         elseif(ityp.eq.-14) then
            um=2.d0*elcon(1,j,nmat)
         endif
!
         do i=1,3
            k=jcoef(i,abs(ityp))
            if(k.eq.0) exit
            if(dabs(elcon(k,j,nmat)).lt.1.d-10) then
               elcon(k,j,nmat)=(0.1d0/um)**i
               write(*,*)
     &             '*WARNING reading *HYPERELASTIC: default value was'
               write(*,*) '         used for compressibility coefficient
     &s'
               write(*,100) i,elcon(k,j,nmat)
            endif
         enddo
      enddo
!
 100  format('   D',i1,' = ',e11.4)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
Subroutine pyspli(kf, kflin, kflch, kflsp)
  Dimension kfl(3)
  kfa = iabs(kf)
  kfs = isign(1, kf)
  kfl(1) = mod(kfa/1000, 10)
  kfl(2) = mod(kfa/100, 10)
  kfl(3) = mod(kfa/10, 10)
  kflr = kflin*kfs
  kflch = 0
  If (kfl(1)==0) Then
    kfl(2) = kfl(2)*(-1)**kfl(2)
    kfl(3) = -kfl(3)*(-1)**iabs(kfl(2))
    If (kflr==kfl(2)) Then
      kflsp = kfl(3)
    Else If (kflr==kfl(3)) Then
      kflsp = kfl(2)
    Else If (iabs(kflr)==21 .And. rlu(0)>0.5) Then
      kflsp = kfl(2)
      kflch = kfl(3)
    Else If (iabs(kflr)==21) Then
      kflsp = kfl(3)
      kflch = kfl(2)
    Else If (kflr*kfl(2)>0) Then
      Call lukfdi(-kflr, kfl(2), kfdump, kflch)
      kflsp = kfl(3)
    Else
      Call lukfdi(-kflr, kfl(3), kfdump, kflch)
      kflsp = kfl(2)
    End If
  Else
    nagr = 0
    Do j = 1, 3
      If (kflr==kfl(j)) nagr = nagr + 1
    End Do
    If (nagr>=1) Then
      ragr = 0.00001 + (nagr-0.00002)*rlu(0)
      iagr = 0
      Do j = 1, 3
        If (kflr==kfl(j)) ragr = ragr - 1.
        If (iagr==0 .And. ragr<=0.) iagr = j
      End Do
    Else
      iagr = int(1.00001+2.99998*rlu(0))
    End If
    id1 = 1
    If (iagr==1) id1 = 2
    If (iagr==1 .And. kfl(3)>kfl(2)) id1 = 3
    id2 = 6 - iagr - id1
    ksp = 3
    If (mod(kfa,10)==2 .And. kfl(1)==kfl(2)) Then
      If (iagr/=3 .And. rlu(0)>0.25) ksp = 1
    Else If (mod(kfa,10)==2 .And. kfl(2)>=kfl(3)) Then
      If (iagr/=1 .And. rlu(0)>0.25) ksp = 1
    Else If (mod(kfa,10)==2) Then
      If (iagr==1) ksp = 1
      If (iagr/=1 .And. rlu(0)>0.75) ksp = 1
    End If
    kflsp = 1000*kfl(id1) + 100*kfl(id2) + ksp
    If (kflin==21) Then
      kflch = kfl(iagr)
    Else If (nagr==0 .And. kflr>0) Then
      Call lukfdi(-kflr, kfl(iagr), kfdump, kflch)
    Else If (nagr==0) Then
      Call lukfdi(10000+kflsp, -kflr, kfdump, kflch)
      kflsp = kfl(iagr)
    End If
  End If
  kflch = kflch*kfs
  kflsp = kflsp*kfs
  Return
End Subroutine pyspli

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT