text:
      SUBROUTINE SSGETD (ELID,TI,GRIDS)
C
C     THIS ROUTINE (CALLED BY -EDTL-) READS ELEMENT TEMPERATURE
C     DATA FROM A PRE-POSITIONED RECORD
C
C     ELID   = ID OF ELEMENT FOR WHICH DATA IS DESIRED
C     TI     = BUFFER DATA IS TO BE RETURNED IN
C     GRIDS  = 0 IF EL-TEMP FORMAT DATA IS TO BE RETURNED
C            = NO. OF GRID POINTS IF GRID POINT DATA IS TO BE RETURNED.
C     ELTYPE = ELEMENT TYPE TO WHICH -ELID- BELONGS
C     OLDEL  = ELEMENT TYPE CURRENTLY BEING WORKED ON (INITIALLY 0)
C     EORFLG = .TRUE. WHEN ALL DATA HAS BEEN EXHAUSTED IN RECORD
C     ENDID  = .TRUE. WHEN ALL DATA HAS BEEN EXHAUSTED WITHIN AN ELEMENT
C              TYPE.
C     BUFFLG = NOT USED
C     ITEMP  = TEMPERATURE LOAD SET ID
C     IDEFT  = NOT USED
C     IDEFM  = NOT USED
C     RECORD = .TRUE. IF A RECORD OF DATA IS INITIALLY AVAILABLE
C     DEFALT = THE DEFALT TEMPERATURE VALUE OR -1 IF IT DOES NOT EXIST
C     AVRAGE = THE AVERAGE ELEMENT TEMPERATURE
C
      LOGICAL         EORFLG   ,ENDID    ,BUFFLG   ,RECORD
      INTEGER         TI(7)    ,GRIDS    ,ELID     ,ELTYPE   ,OLDEL   ,
     1                NAME(2)  ,GPTT     ,DEFALT
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /SYSTEM/ DUM      ,IOUT
      COMMON /SSGETT/ ELTYPE   ,OLDEL    ,EORFLG   ,ENDID    ,BUFFLG  ,
     1                ITEMP    ,IDEFT    ,IDEFM    ,RECORD
      COMMON /LOADX / DUMMY(9) ,GPTT
      COMMON /FPT   / DEFALT
      DATA    NAME  / 4HSSGE,4HTD   /    ,MAXWDS   / 15 /
C
      IF (ITEMP .NE. 0) GO TO 20
      DO 10 I = 1,MAXWDS
   10 TI(I) = 0
      RETURN
C
   20 IF (.NOT.RECORD .OR. EORFLG) GO TO 80
   30 IF (ELTYPE .NE. OLDEL) GO TO 150
      IF (ENDID) GO TO 80
C
C     HERE WHEN ELTYPE IS AT HAND AND END OF THIS TYPE DATA
C     HAS NOT YET BEEN REACHED.  READ AN ELEMENT ID
C
   40 CALL READ (*200,*210,GPTT,ID,1,0,FLAG)
      IF (ID) 50,80,50
   50 IF (IABS(ID) .EQ. ELID) IF (ID) 90,90,70
      IF (ID) 40,40,60
   60 CALL READ (*200,*210,GPTT,TI,NWORDS,0,FLAG)
      GO TO 40
C
C     MATCH ON ELEMNT ID MADE, AND IT WAS WITH DATA.
C     IF QUAD4 OR TRIA3 ELEMENT, SET THE TI(7) FLAG FOR TLQD4D/S (QAUD4)
C     OR TLTR3D/S (TRIA3)
C
   70 CALL READ (*200,*210,GPTT,TI,NWORDS,0,FLAG)
      IF (ELTYPE.NE.64 .AND. ELTYPE.NE.83) RETURN
      TI(7) = 13
      IF (TI(6) .NE. 1) TI(7) = 2
      RETURN
C
C     NO MORE DATA FOR THIS ELEMENT TYPE
C
   80 ENDID = .TRUE.
C
C     NO DATA FOR ELEMENT ID DESIRED, THUS USE DEFALT
C
   90 IF (DEFALT .EQ. -1) GO TO 130
      IF (GRIDS  .GT.  0) GO TO 110
      DO 100 I = 2,MAXWDS
  100 TI(I) = 0
      TI(1) = DEFALT
      IF (ELTYPE .EQ. 34) TI(2) = DEFALT
      RETURN
  110 DO 120 I = 1,GRIDS
  120 TI(I) = DEFALT
      TI(GRIDS+1) = DEFALT
      RETURN
C
C     NO TEMP DATA OR DEFALT
C
  130 WRITE  (IOUT,140) UFM,ELID,ITEMP
  140 FORMAT (A23,' 4017. THERE IS NO TEMPERATURE DATA FOR ELEMENT',I9,
     1        ' IN SET',I9)
      CALL MESAGE (-61,0,0)
C
C     LOOK FOR MATCH ON ELTYPE (FIRST SKIP ANY UNUSED ELEMENT DATA)
C
  150 IF (ENDID) GO TO 180
  160 CALL READ (*200,*210,GPTT,ID,1,0,FLAG)
      IF (ID) 160,180,170
  170 CALL READ (*200,*210,GPTT,TI,NWORDS,0,FLAG)
      GO TO 160
C
C     READ ELTYPE AND COUNT
C
  180 CALL READ (*200,*190,GPTT,TI,2,0,FLAG)
      OLDEL  = TI(1)
      NWORDS = TI(2)
      ENDID  = .FALSE.
      GO TO 30
C
C     END OF RECORD HIT
C
  190 EORFLG = .TRUE.
      GO TO 80
C
  200 CALL MESAGE (-2,GPTT,NAME)
  210 CALL MESAGE (-3,GPTT,NAME)
      RETURN
      END

ext: f


module intwmod
!$$$ module documentation block
!           .      .    .                                       .
! module:   intwmod    module for intw and its tangent linear intw_tl
!   prgmmr:
!
! abstract: module for intw and its tangent linear intw_tl
!
! program history log:
!   2005-05-16  Yanqiu zhu - wrap intw and its tangent linear intw_tl into one module
!   2005-11-16  Derber - remove interfaces
!   2008-11-26  Todling - remove intw_tl; add interface back
!   2009-08-13  lueken - update documentation
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - implemented obs adjoint test
!   2014-04-12       su - add non linear qc from Purser's scheme
!   2015-02-26       su - add njqc as an option to chose new non linear qc
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!
! subroutines included:
!   sub intw_
!
! variable definitions:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

use m_obsNode, only: obsNode
use m_wNode, only: wNode
use m_wNode, only: wNode_typecast
use m_wNode, only: wNode_nextcast
use m_obsdiagNode, only: obsdiagNode_set
implicit none

PRIVATE
PUBLIC intw

interface intw; module procedure &
          intw_
end interface

contains

subroutine intw_(whead,rval,sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    intw        apply nonlin qc obs operator for winds
!   prgmmr: derber           org: np23                date: 1991-02-26
!
! abstract: apply observation operator and adjoint for winds with
!             nonlinear qc operator
!
! program history log:
!   1991-02-26  derber
!   1997-12-12  weiyu yang
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-08-02  treadon - add only to module use, add intent in/out
!   2004-10-09  parrish - add nonlinear qc option
!   2005-03-01  parrish - nonlinear qc change to account for inflated obs error
!   2005-04-11  treadon - merge intw and intw_qc into single routine
!   2005-08-02  derber  - modify for variational qc parameters for each ob
!   2005-09-28  derber  - consolidate location and weight arrays
!   2005-10-21  su      - modify for variational qc
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!                       - unify NL qc
!   2007-03-19  tremolet - binning of observations
!   2007-06-05  tremolet - use observation diagnostics structure
!   2007-07-09  tremolet - observation sensitivity
!   2008-01-04  tremolet - Don't apply H^T if l_do_adjoint is false
!   2008-11-28  todling  - turn FOTO optional; changed ptr%time handle
!   2010-03-13  todling  - update to use gsi_bundle; update interface
!   2012-09-14  Syed RH Rizvi, NCAR/NESL/MMM/DAS  - introduced ladtest_obs
!   2014-04-12       su - add non linear qc from Purser's scheme
!   2014-12-03  derber  - modify so that use of obsdiags can be turned off
!   2015-12-21  yang    - Parrish's correction to the previous code in new varqc.
!
!   input argument list:
!     whead    - obs type pointer to obs structure
!     su       - u increment in grid space
!     sv       - v increment in grid space
!     ru
!     rv
!
!   output argument list:
!     ru       - u results from observation operator
!     rv       - v results from observation operator
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: half,one,tiny_r_kind,cg_term,r3600,two
  use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
  use qcmod, only: nlnqc_iter,varqc_iter,njqc,vqc
  use jfunc, only: jiter
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use gsi_4dvar, only: ladtest_obs
  implicit none

! Declare passed variables
  class(obsNode), pointer,intent(in   ) :: whead
  type(gsi_bundle),       intent(in   ) :: sval
  type(gsi_bundle),       intent(inout) :: rval

! Declare local variables
  integer(i_kind) i1,i2,i3,i4,i5,i6,i7,i8,ier,istatus
! real(r_kind) penalty
  real(r_kind) valu,valv,w1,w2,w3,w4,w5,w6,w7,w8
  real(r_kind) cg_w,p0,gradu,gradv,wnotgross,wgross,term,w_pg
  real(r_kind),pointer,dimension(:) :: su,sv
  real(r_kind),pointer,dimension(:) :: ru,rv
  type(wNode), pointer :: wptr

!  If no w data return
  if(.not. associated(whead))return

  ier=0
  call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
  call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
  if(ier/=0)return

  !wptr => whead
  wptr => wNode_typecast(whead)
  do while(associated(wptr))
     i1=wptr%ij(1)
     i2=wptr%ij(2)
     i3=wptr%ij(3)
     i4=wptr%ij(4)
     i5=wptr%ij(5)
     i6=wptr%ij(6)
     i7=wptr%ij(7)
     i8=wptr%ij(8)
     w1=wptr%wij(1)
     w2=wptr%wij(2)
     w3=wptr%wij(3)
     w4=wptr%wij(4)
     w5=wptr%wij(5)
     w6=wptr%wij(6)
     w7=wptr%wij(7)
     w8=wptr%wij(8)

!    Forward model
     valu=w1* su(i1)+w2* su(i2)+w3* su(i3)+w4* su(i4)+&
          w5* su(i5)+w6* su(i6)+w7* su(i7)+w8* su(i8)
     valv=w1* sv(i1)+w2* sv(i2)+w3* sv(i3)+w4* sv(i4)+&
          w5* sv(i5)+w6* sv(i6)+w7* sv(i7)+w8* sv(i8)

     if(luse_obsdiag)then
        if (lsaveobsens) then
           gradu = valu*wptr%raterr2*wptr%err2
           gradv = valv*wptr%raterr2*wptr%err2
           !-- wptr%diagu%obssen(jiter) = gradu
           !-- wptr%diagv%obssen(jiter) = gradv
           call obsdiagNode_set(wptr%diagu,jiter=jiter,obssen=gradu)
           call obsdiagNode_set(wptr%diagv,jiter=jiter,obssen=gradv)
        else
           if (wptr%luse) then
              !-- wptr%diagu%tldepart(jiter)=valu
              !-- wptr%diagv%tldepart(jiter)=valv
              call obsdiagNode_set(wptr%diagu,jiter=jiter,tldepart=valu)
              call obsdiagNode_set(wptr%diagv,jiter=jiter,tldepart=valv)
           endif
        endif
     endif

     if (l_do_adjoint) then
        if (.not. lsaveobsens) then
           if( .not. ladtest_obs) then
              valu=valu-wptr%ures
              valv=valv-wptr%vres
           end if

!          gradient of nonlinear operator

           if (vqc .and. nlnqc_iter .and. wptr%pg > tiny_r_kind .and.  &
                                wptr%b  > tiny_r_kind) then
              w_pg=wptr%pg*varqc_iter
              cg_w=cg_term/wptr%b
              wnotgross= one-w_pg
              wgross =w_pg*cg_w/wnotgross                ! wgross is gama in Enderson
              p0=wgross/(wgross+                      &  ! p0 is P in Enderson
              exp(-half*wptr%err2*(valu**2+valv**2)))
              term=one-p0                                !  term is Wqc in Enderson
              valu = valu*term
              valv = valv*term
           endif
           if (njqc .and. wptr%jb  > tiny_r_kind .and. wptr%jb <10.0_r_kind) then
              valu=sqrt(two*wptr%jb)*tanh(sqrt(wptr%err2)*valu/sqrt(two*wptr%jb))
              valv=sqrt(two*wptr%jb)*tanh(sqrt(wptr%err2)*valv/sqrt(two*wptr%jb))
              gradu = valu*wptr%raterr2*sqrt(wptr%err2)
              gradv = valv*wptr%raterr2*sqrt(wptr%err2)
           else
              gradu = valu*wptr%raterr2*wptr%err2
              gradv = valv*wptr%raterr2*wptr%err2
           endif
           if( ladtest_obs) then
              gradu = valu
              gradv = valv
           end if
        endif

!       Adjoint
        ru(i1)=ru(i1)+w1*gradu
        ru(i2)=ru(i2)+w2*gradu
        ru(i3)=ru(i3)+w3*gradu
        ru(i4)=ru(i4)+w4*gradu
        ru(i5)=ru(i5)+w5*gradu
        ru(i6)=ru(i6)+w6*gradu
        ru(i7)=ru(i7)+w7*gradu
        ru(i8)=ru(i8)+w8*gradu

        rv(i1)=rv(i1)+w1*gradv
        rv(i2)=rv(i2)+w2*gradv
        rv(i3)=rv(i3)+w3*gradv
        rv(i4)=rv(i4)+w4*gradv
        rv(i5)=rv(i5)+w5*gradv
        rv(i6)=rv(i6)+w6*gradv
        rv(i7)=rv(i7)+w7*gradv
        rv(i8)=rv(i8)+w8*gradv

     endif

     !wptr => wptr%llpoint
     wptr => wNode_nextcast(wptr)

  end do
  return
end subroutine intw_

end module intwmod


text:
       integer*1 one
       integer*2 two
       parameter (one=1)
       parameter (two=2)
       select case (I)
       case (one)
       case (two)
       end select
       end

max_stars_repo_name: ArrogantWombatics/openbsd-src

max_issues_repo_name: ArrogantWombatics/openbsd-src

max_forks_repo_name: ArrogantWombatics/openbsd-src


text:
SUBROUTINE auxval
!
!   Set auxiliary values
!
  USE basic
  USE fields
  USE tracers
  USE numerics
  IMPLICIT none
!
  INCLUDE 'mpif.h'
!   Local vars and arrays
  INTEGER :: k1,j2
  DOUBLE PRECISION :: length
!
  INTEGER :: ierr, myid, master, numprocs
!    Preparing for parallelization
  CALL mpi_comm_size(MPI_COMM_WORLD,numprocs,ierr)
  CALL mpi_comm_rank(MPI_COMM_WORLD,myid,ierr)
  master = 0
!________________________________________________________________________________
  IF (myid == master) THEN
      WRITE(*,'(a/)') '=== Set auxiliary values ==='
  END IF
!________________________________________________________________________________
!
! Pay attention that stagbin is a multiple of numprocs (number of processors)
   stagbin = stagbin/numprocs
   ntracer = stagpart*stagbin
   IF (all_bins_same_pot == .FALSE.) THEN
	numtpts = staglength+(stagbin-1)*stagdist+1
   ELSE
	numtpts = staglength
   END IF

   length = fieldtimestep*staglength/tracerdtfactor

   IF (length > nrun) THEN
      STOP 'length is greater than nrun...stopping'
   END IF

   IF (Bfield_type=='slab') THEN
      Rcurv = 0.
      E_basis = 'cart'
      IF (Bx > 0 .OR. Bz > 0) THEN
         STOP 'in slab field, only abs(By) can be larger than zero...stopping'
      END IF
   ELSE
      E_basis = 'tilt'
   END IF

   IF (Bratio > 0) THEN
      theta=atan(Bratio)
      sineoftheta = sin(theta)
      cosineoftheta = cos(theta)
   ELSE
      theta = 0.
      sineoftheta = 0.
      cosineoftheta = 1.
   END IF

   tracerdeltat=dt*tracerdtfactor

   IF (Lr .NE. 0.) THEN
     kr = 2*pi/Lr
   ELSE
     kz = 0.
   END IF

   IF (Lz .NE. 0.) THEN
     kz = 2*pi/Lz
   ELSE
     kz = 0.
   END IF

   dr = Lr/numrpts
   dz = Lz/numzpts

   IF (ODE_solver=='boris' .OR. ODE_solver=='rk4_3d') THEN
        ndim_tracer=3
   ELSE IF (ODE_solver=='sympl') THEN
        ndim_tracer=2
   ELSE
       write(*,*)'Error: Invalid ODE solver'
   END IF
!
! Allocate memory for global arrays
!
  CALL memory('i')
!

  IF (Efield_type=='kicks' .AND. readkicks) THEN
     IF (myid == master)THEN
	OPEN(UNIT = 7,FILE = "./kickrph.Gauss")
	DO j2 = 1,1000
	   DO k1 = 1,323
	      READ(7,*) kicksteps(j2,k1)
	   END DO
	END DO
     END IF
     CALL MPI_BCAST(kicksteps,1000*323,MPI_INTEGER,master,MPI_COMM_WORLD,ierr)
  END IF

END SUBROUTINE auxval

max_stars_repo_name: alexbovet/fastpart

max_issues_repo_name: alexbovet/fastpart

max_forks_repo_name: alexbovet/fastpart


text:
      PROGRAM xchint
C     driver for routine chint
      INTEGER NVAL
      REAL PIO2
      PARAMETER(NVAL=40, PIO2=1.5707963)
      INTEGER i,mval
      REAL a,b,chebev,fint,func,x,c(NVAL),cint(NVAL)
      EXTERNAL fint,func
      a=-PIO2
      b=PIO2
      call chebft(a,b,c,NVAL,func)
C     test integral
10    write(*,*) 'How many terms in Chebyshev evaluation?'
      write(*,'(1x,a,i2,a)') 'Enter n between 6 and ',NVAL,
     *     '. Enter n=0 to END.'
      read(*,*) mval
      if ((mval.le.0).or.(mval.gt.NVAL)) goto 20
      call chint(a,b,c,cint,mval)
      write(*,'(1x,t10,a,t19,a,t29,a)') 'X','Actual','Cheby. Integ.'
      do 11 i=-8,8,1
        x=i*PIO2/10.0
        write(*,'(1x,3f12.6)') x,fint(x)-fint(-PIO2),
     *       chebev(a,b,cint,mval,x)
11    continue
      goto 10
20    END
      REAL FUNCTION func(x)
      REAL x
      func=(x**2)*(x**2-2.0)*sin(x)
      END
      REAL FUNCTION fint(x)
C     integral of FUNC
      REAL x
      fint=4.0*x*((x**2)-7.0)*sin(x)-((x**4)-14.0*(x**2)+28.0)*cos(x)
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xchint.for


text:
      subroutine cdfhwwcuts(p,maxparts,passed)
c--- given the event momenta in p (with maxparts entries),
c--- performs H->WW search cuts along the lines of CDF Note 9887;
c--- a point that fails the cuts returns passed=.false.
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'plabel.f'
      integer ilept(2),ile,ineut(2),inu,i,j,maxparts
      logical passed
      double precision pt,etarap,p(mxpart,4),etvec(2),
     . missinget,m_ll,R,pt1,pt2,eta1,eta2,phi,phimin
c--- pt and eta thresholds for trigger lepton
      double precision trigpt,trigeta
c--- pt and eta thresholds for second lepton (should be looser)
      double precision scndpt,scndeta
c--- minimum threshold for dilepton invariant mass
      double precision m_llcut
c--- minimum threshold for missing Et "star" or "spec" (see defn below)
      double precision missingetcut
c--- size of cone for lepton isolation
      double precision Risol
c--- fractional pt threshold for lepton isolation
      double precision threshisol
      double precision mtmax,mtmin,mt45,m45,et_vec(4),etmiss,pttwo
      data trigpt,trigeta/20d0,0.8d0/
      data scndpt,scndeta/10d0,1.1d0/
      data m_llcut/16d0/
      data missingetcut/25d0/
      data Risol/0.4d0/
      data threshisol/0.1d0/

c--- default behaviour is that the cuts have been passed
      passed=.true.

c--- initialize all variables
      etvec(1)=0d0
      etvec(2)=0d0

************************** START BASIC CUTS ****************************

      ile=0
      inu=0
c--- identify the two leptons and neutrinos
      do j=3,maxparts
        if (     (plabel(j).eq.'el') .or. (plabel(j).eq.'ea')
     .      .or. (plabel(j).eq.'ml') .or. (plabel(j).eq.'ma')) then
          ile=ile+1
          ilept(ile)=j
        elseif  ((plabel(j).eq.'nl') .or. (plabel(j).eq.'na')) then
          inu=inu+1
          ineut(inu)=j
        endif
      enddo

c--- print warning if there are not 2 leptons
      if (ile .ne. 2) then
        write(6,*) 'WARNING: cdfhwwcuts.f must be upgraded, <2 leptons'
      stop
      endif

c--- compute lepton pt and rapidity
      pt1=pt(ilept(1),p)
      pt2=pt(ilept(2),p)
      eta1=abs(etarap(ilept(1),p))
      eta2=abs(etarap(ilept(2),p))

c--- CUT on basic acceptance cuts for the leptons
      if ((  (pt1 .gt. trigpt) .and. (eta1 .lt. trigeta)
     .  .and.(pt2 .gt. scndpt) .and. (eta2 .lt. scndeta)) .or.
     .    (  (pt2 .gt. trigpt) .and. (eta2 .lt. trigeta)
     .  .and.(pt1 .gt. scndpt) .and. (eta1 .lt. scndeta)) ) then
        continue ! cut passed
      else
        passed=.false.
      return
      endif

c--- compute invariant mass of leptons
      m_ll=dsqrt((p(ilept(1),4)+p(ilept(2),4))**2
     .          -(p(ilept(1),1)+p(ilept(2),1))**2
     .          -(p(ilept(1),2)+p(ilept(2),2))**2
     .          -(p(ilept(1),3)+p(ilept(2),3))**2)

c--- CUT on dilepton invariant mass
      if (m_ll .lt. m_llcut) then
        passed=.false.
      return
      endif

c--- CUT on the isolation of the leptons
      do i=1,2
      do j=7,maxparts ! jets correspond to indices 7,...,maxparts
        if (R(p,ilept(i),j) .lt. Risol) then
        if (pt(j,p) .gt. pt(ilept(i),p)*threshisol) then
          passed=.false.
          return
        endif
      endif
      enddo
      enddo

c--- print warning if there are not 2 neutrinos
      if (inu .ne. 2) then
        write(6,*) 'WARNING: cdfhwwcuts.f must be upgraded, <2 neuts'
      stop
      endif

c--- form the missing et vector
      do j=1,2
        etvec(j)=etvec(j)+p(ineut(1),j)+p(ineut(2),j)
      enddo
      missinget=dsqrt(etvec(1)**2+etvec(2)**2)

c--- find the object (charged lepton or jet) that is closest in
c--- azimuthal angle to missing Et vector
      phimin=10d0
      do i=3,maxparts
        if ((i .ne. ineut(1)) .and. (i .ne. ineut(2))) then
          phi= (p(i,1)*etvec(1)+p(i,2)*etvec(2))
     .         /dsqrt((p(i,1)**2+p(i,2)**2)*missinget**2)
          if (phi .gt. +0.9999999D0) phi=+1D0
          if (phi .lt. -0.9999999D0) phi=-1D0
          phi=dacos(phi)
        endif
        if (phi .lt. phimin) phimin=phi
      enddo

c--- convert missing Et into missing Et "star" (or "spec")
      if (phi .lt. pi/2d0) then
        missinget=missinget*dsin(phi)
      endif

c--- CUT on missing Et
      if (missinget .lt. missingetcut) then
        passed=.false.
      return
      endif

! OPTIONAL mt cut
      m45=0d0
      do i=1,4
         if(i.ne.4) then
            m45=m45-(p(4,i)+p(5,i))**2
         else
            m45=m45+(p(4,i)+p(5,i))**2
         endif
      enddo

      mtmin=0d0*hmass
      mtmax=hmass
      mt45=0d0
      mt45=(dsqrt(dsqrt(pttwo(4,5,p)**2+m45)+etmiss(p,et_vec))**2)
!      write(6,*) mt45
 !     mt45=dsqrt(max(0d0,dsqrt(mt45-ptsq)))
 !     write(6,*) mt45
 !     pause
!------ Mt cuts
      if((mt45.lt.mtmin).or.(mt45.gt.mtmax)) then
         passed=.false.
         return
      endif

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
module probdata_module

  use network, only: nspec, network_species_index
  use eos_type_module, only: eos_t, eos_input_rt
  use eos_module, only: eos
  use bl_constants_module, only: ZERO, THIRD, HALF, ONE, TWO, THREE, M_PI, FOUR
  use fundamental_constants_module, only: Gconst, M_solar, AU
  use initial_model_module, only: initial_model

  ! Initial stellar properties
  ! Note that the envelope mass is included within the total mass of the star

  double precision, save :: mass_P = ONE
  double precision, save :: mass_S = ONE
  double precision, save :: central_density_P = -ONE
  double precision, save :: central_density_S = -ONE
  double precision, save :: stellar_temp = 1.0d7
  double precision, save :: primary_envelope_mass, secondary_envelope_mass
  double precision, save :: primary_envelope_comp(nspec), secondary_envelope_comp(nspec)

  ! Ambient medium

  double precision, save :: ambient_density = 1.0d-4
  double precision, save :: ambient_temp = 1.0d7
  double precision, save :: ambient_comp(nspec)

  ! Smallest allowed velocity on the grid

  double precision, save :: smallu = ZERO

  ! Parameters for interpolation from 1D model to 3D model:

  ! Number of sub-grid-scale zones to use

  integer, save :: nsub = 1

  ! Default to interpolation that preserves temperature; otherwise, use pressure

  logical, save :: interp_temp = .true.

  ! Method for determining the initial problem setup.
  !
  ! 0 = Collision; distance determined by a multiple of the secondary WD radius
  ! 1 = Keplerian orbit; distance determined by the rotation period
  ! 2 = Keplerian orbit; distance set so that the secondary exactly fills its Roche lobe radius
  ! 3 = Problem 2 with an initial relaxation step
  ! 4 = Free-fall; distance determined by a multiple of the secondary WD radius
  ! 5 = Tidal disruption event; distance determined by a multiple of the WD tidal radius

  integer, save :: problem = 2

  ! If we're automatically determining the initial distance based on the Roche lobe
  ! radii for problem 3, this is the sizing factor we use. We set the default value
  ! to be a large enough distance so that the system is close to stable.

  double precision, save :: roche_radius_factor = 2.0d0

  ! Collision parameters

  ! For a collision, number of (secondary) WD radii to
  ! separate the WDs by.

  double precision, save :: collision_separation = 4.0d0

  ! For a collision, the impact parameter measured in
  ! units of the primary's initial radius.

  double precision, save :: collision_impact_parameter = 0.0d0

  ! TDE parameters

  ! For a TDE, number of WD tidal radii to separate the WD and BH.

  double precision, save :: tde_separation = 8.0d0

  ! For a TDE, the parameter beta: the ratio of the tidal radius to
  ! the Schwarzschild radius of the BH.

  double precision, save :: tde_beta = 6.0d0

  ! For a TDE, should we give the star an initial kick of velocity
  ! corresponding to its parabolic orbit? By default we will, but
  ! this option exists so we can test for HSE.

  integer, save :: tde_initial_velocity = 1

  double precision, save :: tde_tidal_radius
  double precision, save :: tde_schwarzschild_radius
  double precision, save :: tde_pericenter_radius

  ! Binary orbit properties

  double precision, save :: r_P_initial, r_S_initial, a_P_initial, a_S_initial, a
  double precision, save :: v_P_r, v_S_r, v_P_phi, v_S_phi, v_P, v_S
  double precision, save :: center_P_initial(3), center_S_initial(3)
  double precision, save :: orbital_eccentricity = 0.0d0
  double precision, save :: orbital_angle = 0.0d0

  ! Axis is in orbital plane; we measure angle with respect to this axis. Normally the x axis.

  integer, save :: axis_1 = 1

  ! Perpendicular axis in the orbital plane. Normally the y axis.

  integer, save :: axis_2 = 2

  ! Perpendicular to both other axes. Normally the z axis and also the rotation axis.

  integer, save :: axis_3 = 3

  ! Location of the physical center of the problem, as a fraction of domain size

  double precision, save :: center_fracx = HALF
  double precision, save :: center_fracy = HALF
  double precision, save :: center_fracz = HALF

  ! Bulk system motion

  double precision, save :: bulk_velx = ZERO
  double precision, save :: bulk_vely = ZERO
  double precision, save :: bulk_velz = ZERO

  ! Whether we're doing an initialization or a restart

  integer, save :: init

  ! Are we doing a single star simulation?

  logical, save :: single_star

  ! 1D initial models

  type (initial_model) :: model_P, model_S

  ! For the grid spacing for our model, we'll use
  ! 6.25 km. No simulation we do is likely to have a resolution
  ! higher than that inside the stars (it represents
  ! three jumps by a factor of four compared to our
  ! normal coarse grid resolution). By using 4096
  ! grid points, the size of the 1D domain will be 2.56e9 cm,
  ! which is larger than any reasonable mass white dwarf.

  double precision, save :: initial_model_dx = 6.25d5
  integer, save          :: initial_model_npts = 4096

  ! initial_model_mass_tol is tolerance used for getting the total WD mass
  ! equal to the desired mass. It can be reasonably small, since there
  ! will always be a central density value that can give the desired
  ! WD mass on the grid we use.

  double precision, save :: initial_model_mass_tol = 1.d-6

  ! hse_tol is the tolerance used when iterating over a zone to force
  ! it into HSE by adjusting the current density (and possibly
  ! temperature).  hse_tol should be very small (~ 1.e-10).

  double precision, save :: initial_model_hse_tol = 1.d-10

  ! Composition properties of initial models.
  ! We follow the prescription of Dan et al. 2012 for determining
  ! the composition of the WDs. In this approach, below a certain
  ! mass there are pure He WDs; just above that are hybrid WDs
  ! with pure CO cores and a He mantle; above that are pure CO WDs
  ! with slightly more oxygen than carbon; and above that are
  ! ONeMg WDs. All masses are in solar masses.

  double precision, save :: max_he_wd_mass = 0.45d0
  double precision, save :: max_hybrid_wd_mass = 0.6d0
  double precision, save :: hybrid_wd_he_shell_mass = 0.1d0
  double precision, save :: max_co_wd_mass = 1.05d0
  double precision, save :: co_wd_he_shell_mass = 0.0d0

  double precision, save :: hybrid_wd_c_frac = 0.50d0
  double precision, save :: hybrid_wd_o_frac = 0.50d0

  double precision, save :: co_wd_c_frac = 0.40d0
  double precision, save :: co_wd_o_frac = 0.60d0

  double precision, save :: onemg_wd_o_frac  = 0.60d0
  double precision, save :: onemg_wd_ne_frac = 0.35d0
  double precision, save :: onemg_wd_mg_frac = 0.05d0

  ! Tagging criteria

  integer,          save :: max_stellar_tagging_level = 20
  integer,          save :: max_temperature_tagging_level = 20
  integer,          save :: max_center_tagging_level = 20
  double precision, save :: stellar_density_threshold = 1.0d0
  double precision, save :: temperature_tagging_threshold = 5.0d8
  double precision, save :: center_tagging_radius = 0.0d0
  double precision, save :: max_tagging_radius = 0.75d0
  double precision, save :: roche_tagging_factor = 2.0d0

  ! Stores the center of mass location of the stars throughout the run

  double precision, save :: com_P(3), com_S(3)
  double precision, save :: vel_P(3), vel_S(3)

  ! Stores the effective Roche radii

  double precision, save :: roche_rad_P, roche_rad_S

  ! Relaxation parameters for problem 3

  double precision, save :: relaxation_damping_factor = 1.0d-1
  double precision, save :: relaxation_density_cutoff = 1.0d3
  integer,          save :: relaxation_is_done = 0

  ! Radial damping parameters for problem 3

  double precision, save :: radial_damping_factor = 1.0d3
  double precision, save :: initial_radial_velocity_factor = -1.0d-3

  ! Distance (in kpc) used for calculation of the gravitational wave amplitude
  ! (this wil be calculated along all three coordinate axes).

  double precision, save :: gw_dist = 10.0d0

  ! Current value of the dynamical timescale for each star

  double precision, save :: t_ff_P, t_ff_S

  ! Global extrema

  double precision, save :: T_global_max, rho_global_max, ts_te_global_max

  namelist /fortin/ &
       mass_P, mass_S, &
       central_density_P, central_density_S, &
       nsub, &
       roche_radius_factor, &
       problem, &
       collision_separation, &
       collision_impact_parameter, &
       tde_separation, &
       tde_beta, &
       tde_initial_velocity, &
       interp_temp, &
       relaxation_damping_factor, &
       relaxation_density_cutoff, &
       initial_radial_velocity_factor, &
       radial_damping_factor, &
       ambient_density, &
       stellar_temp, &
       ambient_temp, &
       max_he_wd_mass, &
       max_hybrid_wd_mass, hybrid_wd_he_shell_mass, &
       max_co_wd_mass, &
       co_wd_he_shell_mass, &
       hybrid_wd_c_frac, hybrid_wd_o_frac, &
       co_wd_c_frac, co_wd_o_frac, &
       onemg_wd_o_frac, onemg_wd_ne_frac, onemg_wd_mg_frac, &
       orbital_eccentricity, orbital_angle, &
       axis_1, axis_2, axis_3, &
       max_stellar_tagging_level, &
       max_temperature_tagging_level, &
       max_center_tagging_level, &
       stellar_density_threshold, &
       temperature_tagging_threshold, &
       center_tagging_radius, &
       max_tagging_radius, &
       roche_tagging_factor, &
       bulk_velx, bulk_vely, bulk_velz, &
       smallu, &
       center_fracx, center_fracy, center_fracz, &
       initial_model_dx, &
       initial_model_npts, &
       initial_model_mass_tol, &
       initial_model_hse_tol, &
       gw_dist

  ! Stores whether we assert that the simulation has completed.

  logical, save :: jobIsDone = .false.
  logical, save :: signalJobIsNotDone = .false.

  ! Auxiliary data for determining whether the job is done.

  integer, parameter :: num_previous_ener_timesteps = 5
  double precision :: total_ener_array(num_previous_ener_timesteps)

end module probdata_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
subroutine rayleigh (a, u, xlam, eps, nmaxit, ier)

  implicit none

  real, dimension(:,:), intent(in) :: a
  real, dimension(:), intent(out) :: u
  real, intent(inout) :: xlam
  real, intent(in) :: eps
  integer, intent(in) :: nmaxit
  integer, intent(inout) :: ier

  real,dimension(:), allocatable :: v

  integer :: n, iter, i, j
  real :: xn , s, xlam1, xlam2

  n = size(a(1,:))

  allocate(v(n))

  u = 0
  u(1) = 1

  do iter= 1, nmaxit
      do i =1, n
        v(i) = 0
        do j = 1,n
          v(i) = v(i) + a(i,j)*u(j)
        end do
      end do
      s = 0
      xn = 0
      do i =1,n
          s = s + u(i)*v(i)
          xn = xn + v(i)*v(i)
      end do
      xn = sqrt(xn)
      do i = 1, n
        u(i) = v(i)/xn
      end do
      xlam2 = s
      if(iter .ne. 1) then
        if(abs(xlam2 - xlam1)/(abs(xlam1) + 1)<= eps) then
          xlam = xlam2
          ier = 1
          return !salimos
        end if
      end if
      xlam1 = xlam2
  end do
  xlam = xlam2
  ier = 2

end subroutine

max_stars_repo_name: JP-Amboage/Matrix-Numerical-Analysis

max_issues_repo_name: JP-Amboage/Matrix-Numerical-Analysis

max_forks_repo_name: JP-Amboage/Matrix-Numerical-Analysis


text:
program watch_point
      implicit none
      integer, parameter :: x_max = 4, y_max = 5, t_max = 100
      real, parameter :: prob = 0.7, boundary  = 0.0
      real, dimension(x_max, y_max) :: temp
      integer :: t

      call init_temp(temp, boundary)
      call show_temp(temp)
      do t = 1, t_max
          call update_temp(temp)
          call perturb_temp(temp, prob)
      end do
      call show_temp(temp)

contains

      subroutine init_temp(temp, boundary)
          implicit none
          real, dimension(:, :), intent(inout) :: temp
          real, intent(in) :: boundary
          integer :: i
          call random_number(temp)
          do i = 1, size(temp, 1)
              temp(i, 1) = boundary
              temp(i, shape(temp, 2)) = boundary
          end do
          do i = 1, size(temp, 2)
              temp(1, i) = boundary
              temp(shape(temp, 1), i) = boundary
          end do
      end subroutine init_temp

      subroutine update_temp(temp)
          implicit none
          real, dimension(:, :), intent(inout) :: temp
          integer :: i, j
          do j = 2, size(temp, 2) - 1
              do i = 2, size(temp, 1) - 1
                  temp(i, j) = 0.25*(temp(i - 1, j) + &
                                     temp(i + 1, j) + &
                                     temp(i,  j - 1) + &
                                     temp(i, j + 1))
              end do
          end do
      end subroutine update_temp

      subroutine perturb_temp(temp, prob)
          implicit none
          real, dimension(:, :), intent(inout) :: temp
          real, intent(in) :: prob
          integer :: i, j
          real :: r
          call random_number(r)
          if (r > prob) then
              call random_number(r)
              i = 1 + int(r*size(temp, 1))
              call random_number(r)
              j = 1 + int(r*size(temp, 2))
              call random_number(r)
              temp(i, j) = r
          end if
      end subroutine perturb_temp

      subroutine show_temp(temp)
          implicit none
          real, dimension(:, :), intent(inout) :: temp
          integer :: i
          character(len=32) :: fmt_str
          write (fmt_str, '(A, I0, A)') '(', size(temp, 2), 'ES15.5)'
          do i = 1, size(temp, 1)
              print fmt_str, temp(i, :)
          end do
      end subroutine show_temp

end program watch_point

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     This routines is for solving nonlinear structural mechanics
!     problem (pure displacement-based formulation).
!
!--------------------------------------------------------------------

      SUBROUTINE CONSTRUCT_dSOLID(lM, Ag, Yg, Dg)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      TYPE(mshType), INTENT(IN) :: lM
      REAL(KIND=RKIND), INTENT(IN) :: Ag(tDof,tnNo), Yg(tDof,tnNo),
     2   Dg(tDof,tnNo)

      INTEGER(KIND=IKIND) a, e, g, Ac, eNoN, cPhys, iFn, nFn
      REAL(KIND=RKIND) w, Jac, ksix(nsd,nsd)

      INTEGER(KIND=IKIND), ALLOCATABLE :: ptr(:)
      REAL(KIND=RKIND), ALLOCATABLE :: xl(:,:), al(:,:), yl(:,:),
     2   dl(:,:), bfl(:,:), fN(:,:), pS0l(:,:), pSl(:), ya_l(:), N(:),
     3   Nx(:,:), lR(:,:), lK(:,:,:)

      eNoN = lM%eNoN
      nFn  = lM%nFn
      IF (nFn .EQ. 0) nFn = 1

!     STRUCT: dof = nsd
      ALLOCATE(ptr(eNoN), xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN),
     2   dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN),
     3   pSl(nsymd), ya_l(eNoN), N(eNoN), Nx(nsd,eNoN), lR(dof,eNoN),
     4   lK(dof*dof,eNoN,eNoN))

!     Loop over all elements of mesh
      DO e=1, lM%nEl
!        Update domain and proceed if domain phys and eqn phys match
         cDmn  = DOMAIN(lM, cEq, e)
         cPhys = eq(cEq)%dmn(cDmn)%phys
         IF (cPhys .NE. phys_struct) CYCLE

!        Update shape functions for NURBS
         IF (lM%eType .EQ. eType_NRB) CALL NRBNNX(lM, e)

!        Create local copies
         fN   = 0._RKIND
         pS0l = 0._RKIND
         ya_l = 0._RKIND
         DO a=1, eNoN
            Ac = lM%IEN(a,e)
            ptr(a)   = Ac
            xl(:,a)  = x(:,Ac)
            al(:,a)  = Ag(:,Ac)
            yl(:,a)  = Yg(:,Ac)
            dl(:,a)  = Dg(:,Ac)
            bfl(:,a) = Bf(:,Ac)
            IF (ALLOCATED(lM%fN)) THEN
               DO iFn=1, nFn
                  fN(:,iFn) = lM%fN((iFn-1)*nsd+1:iFn*nsd,e)
               END DO
            END IF
            IF (ALLOCATED(pS0)) pS0l(:,a) = pS0(:,Ac)
            IF (cem%cpld) ya_l(a) = cem%Ya(Ac)
         END DO

!        Gauss integration
         lR = 0._RKIND
         lK = 0._RKIND
         DO g=1, lM%nG
            IF (g.EQ.1 .OR. .NOT.lM%lShpF) THEN
               CALL GNN(eNoN, nsd, lM%Nx(:,:,g), xl, Nx, Jac, ksix)
               IF (ISZERO(Jac)) err = "Jac < 0 @ element "//e
            END IF
            w = lM%w(g) * Jac
            N = lM%N(:,g)

            pSl = 0._RKIND
            IF (nsd .EQ. 3) THEN
               CALL STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
     2            pS0l, pSl, ya_l, lR, lK)

            ELSE IF (nsd .EQ. 2) THEN
               CALL STRUCT2D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
     2            pS0l, pSl, ya_l, lR, lK)

            END IF

!           Prestress
            IF (pstEq) THEN
               DO a=1, eNoN
                  Ac = ptr(a)
                  pSn(:,Ac) = pSn(:,Ac) + w*N(a)*pSl(:)
                  pSa(Ac)   = pSa(Ac)   + w*N(a)
               END DO
            END IF
         END DO ! g: loop

!        Assembly
#ifdef WITH_TRILINOS
         IF (eq(cEq)%assmTLS) THEN
            CALL TRILINOS_DOASSEM(eNoN, ptr, lK, lR)
         ELSE
#endif
            CALL DOASSEM(eNoN, ptr, lK, lR)
#ifdef WITH_TRILINOS
         END IF
#endif
      END DO ! e: loop

      DEALLOCATE(ptr, xl, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, N, Nx,
     2   lR, lK)

      RETURN
      END SUBROUTINE CONSTRUCT_dSOLID
!####################################################################
      SUBROUTINE STRUCT3D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
     2   pS0l, pSl, ya_l, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN, nFn
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN),
     2   al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(3,eNoN),
     3   fN(3,nFn), pS0l(6,eNoN), ya_l(eNoN)
      REAL(KIND=RKIND), INTENT(OUT) :: pSl(6)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      INTEGER(KIND=IKIND) :: a, b, i, j, k
      REAL(KIND=RKIND) :: rho, dmp, T1, amd, afl, ya_g, fb(3), ud(3),
     2   NxSNx, BmDBm, F(3,3), S(3,3), P(3,3), Dm(6,6), DBm(6,3),
     3   Bm(6,3,eNoN), S0(3,3)

!     Define parameters
      rho     = eq(cEq)%dmn(cDmn)%prop(solid_density)
      dmp     = eq(cEq)%dmn(cDmn)%prop(damping)
      fb(1)   = eq(cEq)%dmn(cDmn)%prop(f_x)
      fb(2)   = eq(cEq)%dmn(cDmn)%prop(f_y)
      fb(3)   = eq(cEq)%dmn(cDmn)%prop(f_z)
      amd     = eq(cEq)%am*rho + eq(cEq)%af*eq(cEq)%gam*dt*dmp
      afl     = eq(cEq)%af*eq(cEq)%beta*dt*dt
      i       = eq(cEq)%s
      j       = i + 1
      k       = j + 1

!     Inertia, body force and deformation tensor (F)
      ud     = -rho*fb
      F      = 0._RKIND
      F(1,1) = 1._RKIND
      F(2,2) = 1._RKIND
      F(3,3) = 1._RKIND
      S0     = 0._RKIND
      ya_g   = 0._RKIND
      DO a=1, eNoN
         ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a))
         ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a))
         ud(3) = ud(3) + N(a)*(rho*(al(k,a)-bfl(3,a)) + dmp*yl(k,a))

         F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
         F(1,2) = F(1,2) + Nx(2,a)*dl(i,a)
         F(1,3) = F(1,3) + Nx(3,a)*dl(i,a)
         F(2,1) = F(2,1) + Nx(1,a)*dl(j,a)
         F(2,2) = F(2,2) + Nx(2,a)*dl(j,a)
         F(2,3) = F(2,3) + Nx(3,a)*dl(j,a)
         F(3,1) = F(3,1) + Nx(1,a)*dl(k,a)
         F(3,2) = F(3,2) + Nx(2,a)*dl(k,a)
         F(3,3) = F(3,3) + Nx(3,a)*dl(k,a)

         S0(1,1) = S0(1,1) + N(a)*pS0l(1,a)
         S0(2,2) = S0(2,2) + N(a)*pS0l(2,a)
         S0(3,3) = S0(3,3) + N(a)*pS0l(3,a)
         S0(1,2) = S0(1,2) + N(a)*pS0l(4,a)
         S0(2,3) = S0(2,3) + N(a)*pS0l(5,a)
         S0(3,1) = S0(3,1) + N(a)*pS0l(6,a)

         ya_g    = ya_g + N(a)*ya_l(a)
      END DO
      S0(2,1) = S0(1,2)
      S0(3,2) = S0(2,3)
      S0(1,3) = S0(3,1)

!     2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in
!     Voigt notationa (Dm)
      CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm)

!     Prestress
      pSl(1) = S(1,1)
      pSl(2) = S(2,2)
      pSl(3) = S(3,3)
      pSl(4) = S(1,2)
      pSl(5) = S(2,3)
      pSl(6) = S(3,1)
      S = S + S0

!     1st Piola-Kirchhoff tensor (P)
      P = MATMUL(F, S)

      DO a=1, eNoN
         Bm(1,1,a) = Nx(1,a)*F(1,1)
         Bm(1,2,a) = Nx(1,a)*F(2,1)
         Bm(1,3,a) = Nx(1,a)*F(3,1)

         Bm(2,1,a) = Nx(2,a)*F(1,2)
         Bm(2,2,a) = Nx(2,a)*F(2,2)
         Bm(2,3,a) = Nx(2,a)*F(3,2)

         Bm(3,1,a) = Nx(3,a)*F(1,3)
         Bm(3,2,a) = Nx(3,a)*F(2,3)
         Bm(3,3,a) = Nx(3,a)*F(3,3)

         Bm(4,1,a) = (Nx(1,a)*F(1,2) + F(1,1)*Nx(2,a))
         Bm(4,2,a) = (Nx(1,a)*F(2,2) + F(2,1)*Nx(2,a))
         Bm(4,3,a) = (Nx(1,a)*F(3,2) + F(3,1)*Nx(2,a))

         Bm(5,1,a) = (Nx(2,a)*F(1,3) + F(1,2)*Nx(3,a))
         Bm(5,2,a) = (Nx(2,a)*F(2,3) + F(2,2)*Nx(3,a))
         Bm(5,3,a) = (Nx(2,a)*F(3,3) + F(3,2)*Nx(3,a))

         Bm(6,1,a) = (Nx(3,a)*F(1,1) + F(1,3)*Nx(1,a))
         Bm(6,2,a) = (Nx(3,a)*F(2,1) + F(2,3)*Nx(1,a))
         Bm(6,3,a) = (Nx(3,a)*F(3,1) + F(3,3)*Nx(1,a))
      END DO

!     Local residue and tangent matrices
      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*ud(1) + Nx(1,a)*P(1,1) +
     2      Nx(2,a)*P(1,2) + Nx(3,a)*P(1,3))
         lR(2,a) = lR(2,a) + w*(N(a)*ud(2) + Nx(1,a)*P(2,1) +
     2      Nx(2,a)*P(2,2) + Nx(3,a)*P(2,3))
         lR(3,a) = lR(3,a) + w*(N(a)*ud(3) + Nx(1,a)*P(3,1) +
     2      Nx(2,a)*P(3,2) + Nx(3,a)*P(3,3))

         DO b=1, eNoN
!           Geometric stiffness
            NxSNx = Nx(1,a)*S(1,1)*Nx(1,b) + Nx(2,a)*S(2,1)*Nx(1,b) +
     2              Nx(3,a)*S(3,1)*Nx(1,b) + Nx(1,a)*S(1,2)*Nx(2,b) +
     3              Nx(2,a)*S(2,2)*Nx(2,b) + Nx(3,a)*S(3,2)*Nx(2,b) +
     4              Nx(1,a)*S(1,3)*Nx(3,b) + Nx(2,a)*S(2,3)*Nx(3,b) +
     5              Nx(3,a)*S(3,3)*Nx(3,b)
            T1 = amd*N(a)*N(b) + afl*NxSNx

!           Material Stiffness (Bt*D*B)
            DBm = MATMUL(Dm, Bm(:,:,b))

            BmDBm = Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1) +
     2              Bm(3,1,a)*DBm(3,1) + Bm(4,1,a)*DBm(4,1) +
     2              Bm(5,1,a)*DBm(5,1) + Bm(6,1,a)*DBm(6,1)
            lK(1,a,b) = lK(1,a,b) + w*(T1 + afl*BmDBm)

            BmDBm = Bm(1,1,a)*DBm(1,2) + Bm(2,1,a)*DBm(2,2) +
     2              Bm(3,1,a)*DBm(3,2) + Bm(4,1,a)*DBm(4,2) +
     2              Bm(5,1,a)*DBm(5,2) + Bm(6,1,a)*DBm(6,2)
            lK(2,a,b) = lK(2,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,1,a)*DBm(1,3) + Bm(2,1,a)*DBm(2,3) +
     2              Bm(3,1,a)*DBm(3,3) + Bm(4,1,a)*DBm(4,3) +
     2              Bm(5,1,a)*DBm(5,3) + Bm(6,1,a)*DBm(6,3)
            lK(3,a,b) = lK(3,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,2,a)*DBm(1,1) + Bm(2,2,a)*DBm(2,1) +
     2              Bm(3,2,a)*DBm(3,1) + Bm(4,2,a)*DBm(4,1) +
     2              Bm(5,2,a)*DBm(5,1) + Bm(6,2,a)*DBm(6,1)
            lK(dof+1,a,b) = lK(dof+1,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,2,a)*DBm(1,2) + Bm(2,2,a)*DBm(2,2) +
     2              Bm(3,2,a)*DBm(3,2) + Bm(4,2,a)*DBm(4,2) +
     2              Bm(5,2,a)*DBm(5,2) + Bm(6,2,a)*DBm(6,2)
            lK(dof+2,a,b) = lK(dof+2,a,b) + w*(T1 + afl*BmDBm)

            BmDBm = Bm(1,2,a)*DBm(1,3) + Bm(2,2,a)*DBm(2,3) +
     2              Bm(3,2,a)*DBm(3,3) + Bm(4,2,a)*DBm(4,3) +
     2              Bm(5,2,a)*DBm(5,3) + Bm(6,2,a)*DBm(6,3)
            lK(dof+3,a,b) = lK(dof+3,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,3,a)*DBm(1,1) + Bm(2,3,a)*DBm(2,1) +
     2              Bm(3,3,a)*DBm(3,1) + Bm(4,3,a)*DBm(4,1) +
     2              Bm(5,3,a)*DBm(5,1) + Bm(6,3,a)*DBm(6,1)
            lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,3,a)*DBm(1,2) + Bm(2,3,a)*DBm(2,2) +
     2              Bm(3,3,a)*DBm(3,2) + Bm(4,3,a)*DBm(4,2) +
     2              Bm(5,3,a)*DBm(5,2) + Bm(6,3,a)*DBm(6,2)
            lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,3,a)*DBm(1,3) + Bm(2,3,a)*DBm(2,3) +
     2              Bm(3,3,a)*DBm(3,3) + Bm(4,3,a)*DBm(4,3) +
     2              Bm(5,3,a)*DBm(5,3) + Bm(6,3,a)*DBm(6,3)
            lK(2*dof+3,a,b) = lK(2*dof+3,a,b) + w*(T1 + afl*BmDBm)
         END DO
      END DO

      RETURN
      END SUBROUTINE STRUCT3D
!####################################################################
      SUBROUTINE STRUCT2D(eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN,
     2   pS0l, pSl, ya_l, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN, nFn
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(2,eNoN),
     2   al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(2,eNoN),
     3   fN(2,nFn), pS0l(3,eNoN), ya_l(eNoN)
      REAL(KIND=RKIND), INTENT(OUT) :: pSl(3)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      INTEGER(KIND=IKIND) :: a, b, i, j
      REAL(KIND=RKIND) :: rho, dmp, T1, amd, afl, ya_g, fb(2), ud(2),
     2   NxSNx, BmDBm, F(2,2), S(2,2), P(2,2), Dm(3,3), DBm(3,2),
     3   Bm(3,2,eNoN), S0(2,2)

!     Define parameters
      rho     = eq(cEq)%dmn(cDmn)%prop(solid_density)
      dmp     = eq(cEq)%dmn(cDmn)%prop(damping)
      fb(1)   = eq(cEq)%dmn(cDmn)%prop(f_x)
      fb(2)   = eq(cEq)%dmn(cDmn)%prop(f_y)
      amd     = eq(cEq)%am*rho + eq(cEq)%af*eq(cEq)%gam*dt*dmp
      afl     = eq(cEq)%af*eq(cEq)%beta*dt*dt
      i       = eq(cEq)%s
      j       = i + 1

!     Inertia, body force and deformation tensor (F)
      ud     = -rho*fb
      F      = 0._RKIND
      F(1,1) = 1._RKIND
      F(2,2) = 1._RKIND
      S0     = 0._RKIND
      ya_g   = 0._RKIND
      DO a=1, eNoN
         ud(1) = ud(1) + N(a)*(rho*(al(i,a)-bfl(1,a)) + dmp*yl(i,a))
         ud(2) = ud(2) + N(a)*(rho*(al(j,a)-bfl(2,a)) + dmp*yl(j,a))

         F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
         F(1,2) = F(1,2) + Nx(2,a)*dl(i,a)
         F(2,1) = F(2,1) + Nx(1,a)*dl(j,a)
         F(2,2) = F(2,2) + Nx(2,a)*dl(j,a)

         S0(1,1) = S0(1,1) + N(a)*pS0l(1,a)
         S0(2,2) = S0(2,2) + N(a)*pS0l(2,a)
         S0(1,2) = S0(1,2) + N(a)*pS0l(3,a)

         ya_g    = ya_g + N(a)*ya_l(a)
      END DO
      S0(2,1) = S0(1,2)

!     2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in
!     Voigt notation
      CALL GETPK2CC(eq(cEq)%dmn(cDmn), F, nFn, fN, ya_g, S, Dm)

!     Prestress
      pSl(1) = S(1,1)
      pSl(2) = S(2,2)
      pSl(3) = S(1,2)
      S = S + S0

!     1st Piola-Kirchhoff tensor (P)
      P = MATMUL(F, S)

      DO a=1, eNoN
         Bm(1,1,a) = Nx(1,a)*F(1,1)
         Bm(1,2,a) = Nx(1,a)*F(2,1)

         Bm(2,1,a) = Nx(2,a)*F(1,2)
         Bm(2,2,a) = Nx(2,a)*F(2,2)

         Bm(3,1,a) = (Nx(1,a)*F(1,2) + F(1,1)*Nx(2,a))
         Bm(3,2,a) = (Nx(1,a)*F(2,2) + F(2,1)*Nx(2,a))
      END DO

!     Local residue and tangent matrices
      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*ud(1) + Nx(1,a)*P(1,1) +
     2      Nx(2,a)*P(1,2))
         lR(2,a) = lR(2,a) + w*(N(a)*ud(2) + Nx(1,a)*P(2,1) +
     2      Nx(2,a)*P(2,2))

         DO b=1, eNoN
!           Geometric stiffness
            NxSNx = Nx(1,a)*S(1,1)*Nx(1,b) + Nx(2,a)*S(2,1)*Nx(1,b) +
     2              Nx(1,a)*S(1,2)*Nx(2,b) + Nx(2,a)*S(2,2)*Nx(2,b)
            T1 = amd*N(a)*N(b) + afl*NxSNx

!           Material stiffness (Bt*D*B)
            DBm(1,1) = Dm(1,1)*Bm(1,1,b) + Dm(1,2)*Bm(2,1,b) +
     2         Dm(1,3)*Bm(3,1,b)
            DBm(1,2) = Dm(1,1)*Bm(1,2,b) + Dm(1,2)*Bm(2,2,b) +
     2         Dm(1,3)*Bm(3,2,b)

            DBm(2,1) = Dm(2,1)*Bm(1,1,b) + Dm(2,2)*Bm(2,1,b) +
     2         Dm(2,3)*Bm(3,1,b)
            DBm(2,2) = Dm(2,1)*Bm(1,2,b) + Dm(2,2)*Bm(2,2,b) +
     2         Dm(2,3)*Bm(3,2,b)

            DBm(3,1) = Dm(3,1)*Bm(1,1,b) + Dm(3,2)*Bm(2,1,b) +
     2         Dm(3,3)*Bm(3,1,b)
            DBm(3,2) = Dm(3,1)*Bm(1,2,b) + Dm(3,2)*Bm(2,2,b) +
     2         Dm(3,3)*Bm(3,2,b)

            BmDBm = Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1) +
     2         Bm(3,1,a)*DBm(3,1)
            lK(1,a,b) = lK(1,a,b) + w*(T1 + afl*BmDBm)

            BmDBm = Bm(1,1,a)*DBm(1,2) + Bm(2,1,a)*DBm(2,2) +
     2         Bm(3,1,a)*DBm(3,2)
            lK(2,a,b) = lK(2,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,2,a)*DBm(1,1) + Bm(2,2,a)*DBm(2,1) +
     2         Bm(3,2,a)*DBm(3,1)
            lK(dof+1,a,b) = lK(dof+1,a,b) + w*afl*BmDBm

            BmDBm = Bm(1,2,a)*DBm(1,2) + Bm(2,2,a)*DBm(2,2) +
     2         Bm(3,2,a)*DBm(3,2)
            lK(dof+2,a,b) = lK(dof+2,a,b) + w*(T1 + afl*BmDBm)
         END DO
      END DO

      RETURN
      END SUBROUTINE STRUCT2D
!####################################################################
      SUBROUTINE BSTRUCT3D(eNoN, w, N, Nx, dl, hl, nV, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN), hl(eNoN),
     2   dl(tDof,eNoN), nV(3)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      INTEGER(KIND=IKIND) :: i, j, k, a, b
      REAL(KIND=RKIND) :: af, h, Jac, wl, Ku, F(3,3), Fi(3,3), nFi(3),
     2   NxFi(3,eNoN)

      af     = eq(cEq)%af*eq(cEq)%beta*dt*dt
      i      = eq(cEq)%s
      j      = i + 1
      k      = j + 1

      h      = 0._RKIND
      F      = 0._RKIND
      F(1,1) = 1._RKIND
      F(2,2) = 1._RKIND
      F(3,3) = 1._RKIND
      DO a=1, eNoN
         h       = h      + N(a)*hl(a)
         F(1,1)  = F(1,1) + Nx(1,a)*dl(i,a)
         F(1,2)  = F(1,2) + Nx(2,a)*dl(i,a)
         F(1,3)  = F(1,3) + Nx(3,a)*dl(i,a)
         F(2,1)  = F(2,1) + Nx(1,a)*dl(j,a)
         F(2,2)  = F(2,2) + Nx(2,a)*dl(j,a)
         F(2,3)  = F(2,3) + Nx(3,a)*dl(j,a)
         F(3,1)  = F(3,1) + Nx(1,a)*dl(k,a)
         F(3,2)  = F(3,2) + Nx(2,a)*dl(k,a)
         F(3,3)  = F(3,3) + Nx(3,a)*dl(k,a)
      END DO
      Jac = MAT_DET(F, 3)
      Fi  = MAT_INV(F, 3)

      nFi(1) = nV(1)*Fi(1,1) + nV(2)*Fi(2,1) + nV(3)*Fi(3,1)
      nFi(2) = nV(1)*Fi(1,2) + nV(2)*Fi(2,2) + nV(3)*Fi(3,2)
      nFi(3) = nV(1)*Fi(1,3) + nV(2)*Fi(2,3) + nV(3)*Fi(3,3)

      DO a=1, eNoN
         NxFi(1,a) = Nx(1,a)*Fi(1,1) + Nx(2,a)*Fi(2,1) + Nx(3,a)*Fi(3,1)
         NxFi(2,a) = Nx(1,a)*Fi(1,2) + Nx(2,a)*Fi(2,2) + Nx(3,a)*Fi(3,2)
         NxFi(3,a) = Nx(1,a)*Fi(1,3) + Nx(2,a)*Fi(2,3) + Nx(3,a)*Fi(3,3)
      END DO

      wl = w*Jac*h
      DO a=1, eNoN
         lR(1,a) = lR(1,a) - wl*N(a)*nFi(1)
         lR(2,a) = lR(2,a) - wl*N(a)*nFi(2)
         lR(3,a) = lR(3,a) - wl*N(a)*nFi(3)

         DO b=1, eNoN
            Ku = wl*af*N(a)*(nFi(2)*NxFi(1,b) - nFi(1)*NxFi(2,b))
            lK(2,a,b)     = lK(2,a,b)     + Ku
            lK(dof+1,a,b) = lK(dof+1,a,b) - Ku

            Ku = wl*af*N(a)*(nFi(3)*NxFi(1,b) - nFi(1)*NxFi(3,b))
            lK(3,a,b)       = lK(3,a,b)       + Ku
            lK(2*dof+1,a,b) = lK(2*dof+1,a,b) - Ku

            Ku = wl*af*N(a)*(nFi(3)*NxFi(2,b) - nFi(2)*NxFi(3,b))
            lK(dof+3,a,b)   = lK(dof+3,a,b)   + Ku
            lK(2*dof+2,a,b) = lK(2*dof+2,a,b) - Ku
         END DO
      END DO

      RETURN
      END SUBROUTINE BSTRUCT3D
!--------------------------------------------------------------------
      SUBROUTINE BSTRUCT2D(eNoN, w, N, Nx, dl, hl, nV, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: eNoN
      REAL(KIND=RKIND), INTENT(IN) :: w, N(eNoN), Nx(2,eNoN), hl(eNoN),
     2   dl(tDof,eNoN), nV(2)
      REAL(KIND=RKIND), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      INTEGER(KIND=IKIND) :: i, j, a, b
      REAL(KIND=RKIND) :: af, h, Jac, wl, Ku, F(2,2), Fi(2,2), nFi(2),
     2   NxFi(2,eNoN)

      af     = eq(cEq)%af*eq(cEq)%beta*dt*dt
      i      = eq(cEq)%s
      j      = i + 1

      h      = 0._RKIND
      F      = 0._RKIND
      F(1,1) = 1._RKIND
      F(2,2) = 1._RKIND
      DO a=1, eNoN
         h       = h      + N(a)*hl(a)
         F(1,1)  = F(1,1) + Nx(1,a)*dl(i,a)
         F(1,2)  = F(1,2) + Nx(2,a)*dl(i,a)
         F(2,1)  = F(2,1) + Nx(1,a)*dl(j,a)
         F(2,2)  = F(2,2) + Nx(2,a)*dl(j,a)
      END DO
      Jac = MAT_DET(F, 2)
      Fi  = MAT_INV(F, 2)

      nFi(1) = nV(1)*Fi(1,1) + nV(2)*Fi(2,1)
      nFi(2) = nV(1)*Fi(1,2) + nV(2)*Fi(2,2)

      DO a=1, eNoN
         NxFi(1,a) = Nx(1,a)*Fi(1,1) + Nx(2,a)*Fi(2,1)
         NxFi(2,a) = Nx(1,a)*Fi(1,2) + Nx(2,a)*Fi(2,2)
      END DO

      wl = w*Jac*h
      DO a=1, eNoN
         lR(1,a) = lR(1,a) - wl*N(a)*nFi(1)
         lR(2,a) = lR(2,a) - wl*N(a)*nFi(2)

         DO b=1, eNoN
            Ku = wl*af*N(a)*(nFi(2)*NxFi(1,b) - nFi(1)*NxFi(2,b))
            lK(2,a,b)     = lK(2,a,b)     + Ku
            lK(dof+1,a,b) = lK(dof+1,a,b) - Ku
         END DO
      END DO

      RETURN
      END SUBROUTINE BSTRUCT2D
!####################################################################

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


text:
	SUBROUTINE POLY_INIT ( iret )
C************************************************************************
C* POLY_INIT 								*
C*									*
C* This subroutine initializes polygon common block.			*
C*									*
C* POLY_INIT ( IRET )							*
C*									*
C* Input parameters:							*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					 0 = normal return		*
C*									*
C**									*
C* Log:									*
C* T. Lee/SAIC		02/08						*
C************************************************************************
	INCLUDE		'gpolyg.cmn'
C-----------------------------------------------------------------------
	iret = 0
C
C*	Default warning categories.  Set WSCA to a high number because small
C*	craft advisory has not been implemented operationally.  Once in
C*	operation, its value is defined by poly_parm.tbl.  Current value
C*	will be overwritten.
C
	wsca  = 9999.
	wgale = 32.5
	wstrm = 47.5
	whurr = 62.5
C
C*	Initialize wave/wind information.
C
	DO nc = 1, NCATEG
	    ngroup ( nc ) = 0
	    DO ng = 1, MAXGRP
		mpoint ( ng, nc ) = 0
		npoint ( ng, nc ) = 0
		numpts ( ng, nc ) = 0
		DO np = 1, MAXPTS
		    pgrlat ( np, ng, nc ) = RMISSD
		    pgrlon ( np, ng, nc ) = RMISSD
		    kgridu ( np, ng, nc ) = IMISSD
		    kgridv ( np, ng, nc ) = IMISSD
		    isnake ( np, ng, nc ) = IMISSD
		    jsnake ( np, ng, nc ) = IMISSD
		END DO
C
		DO kp = 1, 2
		    imargn ( kp, ng, nc ) = IMISSD
		    jmargn ( kp, ng, nc ) = IMISSD
		END DO
	    END DO
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
c                                                   02/05/93 - ktname0.f
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.... Initialize array ktname with the known data type names.

      subroutine ktname0
      use m_die,  only : die
      use m_stdio,only : stdout

c.......................................................................
c.... Table ktname to hold data type names.

	use m_output, only : output_ObsErrCov
	use m_stdio
	implicit none

      include     "ktmax.h"
      include     "kttabl.h"

c.......................................................................
c.... Data type names.

c     DATA ktname /
c    1 'u_SeaLev' , 'v_SeaLev' , 'p_SeaLev' , 'u_UprAir' , 'v_UprAir',
c    2 'H_UprAir' , 'q_UprAir' /

	character*7 myname
	parameter  (myname='ktname0')

	integer istat

	integer ikt,jkt

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	call rdkttbl(RC_kt,ktmax,ktname,ktunit,ktdesc,ktmvar,istat)
	if(istat.ne.0) call die(myname,'rdkttbl()',istat)

c	..Echo input in its original form
c	=================================
	if(output_ObsErrCov) then
	  write(stdout,'(/a)') RC_kt
	  do ikt=1,ktmax
	    write(stdout,'(i2,4(1x,a),1x,7l2)') ikt,
     &	      ktname(ikt),ktunit(ikt),ktdesc(ikt),'$',
     &	      (ktmvar(ikt,jkt),jkt=1,ikt)
	  end do
	  write(stdout,'(a)') '::'
	endif

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	end
c.

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
*     TG01DD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          LMAX, NMAX, PMAX
      PARAMETER        ( LMAX = 20, NMAX = 20, PMAX = 20)
      INTEGER          LDA, LDC, LDE, LDZ
      PARAMETER        ( LDA = LMAX, LDC = PMAX,
     $                   LDE = LMAX, LDZ = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MIN(LMAX,NMAX)+MAX(LMAX,NMAX,PMAX) )
*     .. Local Scalars ..
      CHARACTER*1      COMPZ
      INTEGER          I, INFO, J, L, N, P
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), C(LDC,NMAX),
     $                 DWORK(LDWORK), E(LDE,NMAX), Z(LDZ,NMAX)
*     .. External Subroutines ..
      EXTERNAL         TG01DD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) L, N, P
      COMPZ = 'I'
      IF ( L.LT.0 .OR. L.GT.LMAX ) THEN
         WRITE ( NOUT, FMT = 99992 ) L
      ELSE
         IF( N.LT.0 .OR. N.GT.NMAX ) THEN
            WRITE ( NOUT, FMT = 99991 ) N
         ELSE
            READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,L )
            READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,L )
            IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
               WRITE ( NOUT, FMT = 99990 ) P
            ELSE
               READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
*              Find the transformed descriptor system pair
*              (A-lambda E,B).
               CALL TG01DD( COMPZ, L, N, P, A, LDA, E, LDE, C, LDC,
     $                      Z, LDZ, DWORK, LDWORK, INFO )
*
               IF( INFO.NE.0 ) THEN
                  WRITE ( NOUT, FMT = 99998 ) INFO
               ELSE
                  WRITE ( NOUT, FMT = 99997 )
                  DO 10 I = 1, L
                     WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
   10             CONTINUE
                  WRITE ( NOUT, FMT = 99996 )
                  DO 20 I = 1, L
                     WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N )
   20             CONTINUE
                  WRITE ( NOUT, FMT = 99994 )
                  DO 30 I = 1, P
                     WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N )
   30             CONTINUE
                  WRITE ( NOUT, FMT = 99993 )
                  DO 40 I = 1, N
                     WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N )
   40             CONTINUE
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' TG01DD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01DD = ',I2)
99997 FORMAT (/' The transformed state dynamics matrix A*Z is ')
99996 FORMAT (/' The transformed descriptor matrix E*Z is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (/' The transformed input/state matrix C*Z is ')
99993 FORMAT (/' The right transformation matrix Z is ')
99992 FORMAT (/' L is out of range.',/' L = ',I5)
99991 FORMAT (/' N is out of range.',/' N = ',I5)
99990 FORMAT (/' P is out of range.',/' P = ',I5)
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


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

C=======================================================================
      SUBROUTINE PLTISP(SV,CV)
      REAL DEVCAP(23)
      REAL DEFOUT(7)
      COMMON /STATUS/DEVCAP,DEFOUT
      REAL DEVP(5)
      COMMON /DEVICE/DEVP
      REAL COLP(3)
      REAL PALETT(3,16)
      COMMON /COLOR/COLP,PALETT
      REAL TEXTP(40)
      COMMON /TEXT/TEXTP
      REAL VECTP(5)
      REAL XCUR
      REAL YCUR
      COMMON /VECTRC/VECTP,XCUR,YCUR
      INTEGER IDEX(200,2)
      INTEGER NVECT(200,2)
      REAL XSIZE(200,2)
      REAL YSIZE(200,2)
      REAL X0(2300,2)
      REAL Y0(2300,2)
      REAL X1(2300,2)
      REAL Y1(2300,2)
      COMMON /FONT/IDEX,NVECT,XSIZE,YSIZE,X0,Y0,X1,Y1
      REAL GRAPHP(100)
      COMMON /GRAPH/GRAPHP
      COMMON /MAPPAR/MAPP(11)
      REAL MAPP
      COMMON /STORAG/MEMORY(1000)

      IF (SV.LT.0.) THEN
         SV = 0.
      END IF

      IF (SV.GT.1.) THEN
         SV = 1.
      END IF

      CV = COLP(2) + REAL(NINT(SV*COLP(3)))
      IF (CV.GT.DEVCAP(4)-1.) THEN
         CV = DEVCAP(4) - 1.
      END IF

      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
program main
  call sggsvd()
  call dggsvd()
  call cggsvd()
  call zggsvd()
end program

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
      PROGRAM PARTONSG
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      character*5 flav(6)

      COMMON/IOUT/IPRINT
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/RESULT/S1,S2,S3,S4
      COMMON/IFLAVOR/IFLAVOR

      data flav/'qg   ','gg   ','qqb  ','qqbp ','qq   ','qqp  '/

      EXTERNAL BORN, SOVI, HARD, SOVISCA, HARDSCA

      CALL RSTART(12,34,56,78)

C --- PARAMETER TO CHANGE  -------------------------------------

C***  1:Q G   2:G G   3:Q QBAR   4:QP QBAR   5:Q Q   6:Q QP

      DO 100 IFLAVOR = 2,2

C***  THE MASSES:

      MS = 280.D0
      MG = 200.D0
      MT = 175.D0

C***  START VALUE IN LOG(ETA), NUMBERS OF STEPS, STEP WIDTH IN LOG(ETA)

      ETALOG = -3D0
      INUM = 31
      STEP = 0.2D0

      ETALOG = -0.54D0
      INUM = 20
      STEP = 0.01D0

C***  THE NUMBER OF VEGAS CALLS ( DEFAULT = 500 )

      IBORN = 1000
      ISOVI = 500
      IHARD = 500

C***  PRINT VEGAS STATISTICS OR NOT

      IPRINT = 0

C --- CONSTANTS DON'T CHANGE PLEASE ------------------------------

      PI = 4.D0*DATAN(1.D0)
      CONV = 389379660.D0
      ALPHAS = 0.1D0
      SCA = 4.D0

C -------------------------------------------------------------

c$$$      IF (IFLAVOR.EQ.1) OPEN(10,FILE="partonsg_qg.dat")
c$$$      IF (IFLAVOR.EQ.2) OPEN(10,FILE="partonsg_gg.dat")
c$$$      IF (IFLAVOR.EQ.3) OPEN(10,FILE="partonsg_qb.dat")
c$$$      IF (IFLAVOR.EQ.4) OPEN(10,FILE="partonsg_qbp.dat")
c$$$      IF (IFLAVOR.EQ.5) OPEN(10,FILE="partonsg_qq.dat")
c$$$      IF (IFLAVOR.EQ.6) OPEN(10,FILE="partonsg_qp.dat")

      PRINT 20, MS,MG,MT
      PRINT 30
c$$$      WRITE(10, 20) MS,MG,MT
c$$$      WRITE(10, 30)

 20   FORMAT(' MS = ',G10.4,' MG= ',G10.4,' MT = ',G10.4)
 30   FORMAT('  ETA        BORN          SOVI       REL     ',
     +   '   HARD       REL       SCALE')

C --- INTEGRATION BY VEGAS ------------------------------------

      DO 100  I = 1,INUM

        ETA = 10.D0**(ETALOG)
        ETALOG = ETALOG +STEP

        S = (MG +MS)**2*(ETA + 1.D0)
        DEL = 1.D-6*MG*MS*(1 -(MS+MG)**2/S)

        IF (IFLAVOR.EQ.1) THEN
           CALL INTEG(BORN,1,IBORN,3,5*IBORN,5,1.D-4,RES)
           RESBORN = (MS+MG)**2/4.D0/ALPHAS**2/CONV*S1
           DEVBORN = (MS+MG)**2/4.D0/ALPHAS**2/CONV*S2

           CALL INTEG(SOVI,1,ISOVI,3,5*ISOVI,5,5.D-4,RES)
           RESSOVI = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S1
           DEVSOVI = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S2
           RELSOVI = DABS(DEVSOVI/RESSOVI)

           CALL INTEG(SOVISCA,1,ISOVI,3,5*ISOVI,5,5.D-4,RES)
           RESSOVISCA=(MS+MG)**2/4D0/ALPHAS**3/4D0/PI/CONV/DLOG(SCA)*S1
           DEVSOVISCA=(MS+MG)**2/4D0/ALPHAS**3/4D0/PI/CONV/DLOG(SCA)*S2

           CALL INTEG(HARD,3,ISOVI,3,5*ISOVI,5,5.D-4,RES)
           RESHARD = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S1
           DEVHARD = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S2
           RELHARD = DABS(DEVHARD/RESHARD)

           CALL INTEG(HARDSCA,2,ISOVI,3,5*ISOVI,5,5.D-4,RES)
          RESHARDSCA=(MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S1
          DEVHARDSCA=(MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S2

        ELSE

           CALL INTEG(HARD,3,IHARD,5,5*IHARD,10,5.D-4,RES)
           RESHARD = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S1
           DEVHARD = (MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV*S2
           RELHARD = DABS(DEVHARD/RESHARD)

           CALL INTEG(HARDSCA,2,ISOVI,3,5*ISOVI,5,5.D-4,RES)
          RESHARDSCA=(MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S1
          DEVHARDSCA=(MS+MG)**2/4.D0/ALPHAS**3/4.D0/PI/CONV/DLOG(SCA)*S2

          RESBORN = 0D0
          RESSOVI = 0D0
          RELSOVI = 0D0
          RESSOVISCA = 0D0

      ENDIF

      RESSCALE = RESSOVISCA +RESHARDSCA

      PRINT 50,ETA,RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE
c$$$      WRITE(10,50)ETA,RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE
c$$$      PRINT 51,flav(iflavor),ETA,
c$$$     +     RESBORN,RESSOVI,RELSOVI,RESHARD,RELHARD, RESSCALE

 51   FORMAT(A5,G9.3,' ',G11.5,'  ',G11.5,' ',F5.4,'  ',G11.5,' ',
     +     F5.4,'  ',G11.5)

 50   FORMAT(G10.4,' ',G11.5,'  ',G11.5,' ',G8.2,'  ',G11.5,' ',
     +     G8.2,'  ',G11.5)
 100  CONTINUE

      END
C ======================================================================

      DOUBLE PRECISION FUNCTION BORN(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      XLAM = DSQRT((S -MG**2 -MS**2)**2 -4*MG**2*MS**2)
      TG = -S/2.D0*(1.D0 +(MG**2 -MS**2)/S -(2*X -1.D0)*XLAM/S)
      IF (IFLAVOR.EQ.1) BORN = XLAM*DSGQGB(ALPHAS,S,TG,MS,MG)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION SOVI(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      XLAM = DSQRT((S -MG**2 -MS**2)**2 -4*MG**2*MS**2)
      TG = -S/2.D0*(1.D0 +(MG**2 -MS**2)/S -(2*X -1.D0)*XLAM/S)
      IF (IFLAVOR.EQ.1) SOVI = XLAM*DSGQGV(ALPHAS,S,TG,MS,MG,MT)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION HARD(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(3)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      Y = VAR(2)
      Z = VAR(3)
      IF (IFLAVOR.EQ.1) THEN
         XLAM = DSQRT((S -MG**2 -MS**2)**2 -4*MG**2*MS**2)
         TG = -(S +MG**2 -MS**2 -DEL)/2.D0
     +        + (X -0.5D0) *DSQRT((S+MG**2-MS**2-DEL)**2 -4*MG**2*S)
         S4MAX = S +TG + MG**2 -MS**2 +MG**2*S/TG
         S4 = ( S4MAX -DEL)*Y +DEL
         PREF = S4MAX*XLAM
      ELSE
C***     EPSG = ( GAMMA(GLUINO)/MASS(GLUINO) )**2
C***     EPSS = ( GAMMA(SQUARK)/MASS(SQUARK) )**2
         EPSS = 1.D-10
         EPSG = 1.D-10
         M2 = MG**2 -MS**2
         S4MAX0 = S + M2 -2.D0 * DSQRT(S*MG**2)
         GAMG = DSQRT(EPSG)*MG**2
         TYMIN = DATAN(-M2/GAMG)
         TYMAX = DATAN((S4MAX0-M2)/GAMG)
         TY = (TYMAX -TYMIN) * X + TYMIN
         PREFTY = (TYMAX -TYMIN)/GAMG
         S40 = GAMG*DTAN(TY) + M2
         S4G0 = S40 -M2

         PRETG0 = DSQRT((S - S4G0)**2 - 4*S*MG**2)
         TG0 = - (S -S4G0)/2.D0 +(Y -0.5D0)*PRETG0
         PRETG1 = 1
         IF(S**2-4*MG**2*S.GT.0.D0)THEN
          PRETG1 = DSQRT(S**2 - 4*S*MG**2)
         ENDIF
         TG1 = - S/2.D0 +(Y -0.5D0)*PRETG1
         PREF0 = PREFTY * PRETG0
         PREF1 = PREFTY * PRETG1

         S3GMAX2 = S - M2 - 2*DSQRT(S*MS**2)
         GAMS = DSQRT(EPSS)*MS**2
         TZMIN = DATAN(+M2/GAMS)
         TZMAX = DATAN((S3GMAX2+M2)/GAMS)
         TZ = (TZMAX -TZMIN) * X + TZMIN
         PREFTZ = (TZMAX -TZMIN)/GAMS
         S32 = GAMS*DTAN(TZ)
         S33 = 0D0
         S3G2 = S32 -M2
         S3G3 = S33 -M2

         PRES42 = S3G2/(S32+MS**2)*DSQRT((S-S32)**2 -4*MS**2*S)
         PRES43 = 1
         IF((S-S33)**2 -4*MS**2*S.GT.0.D0)THEN
          PRES43 = S3G3/(S33+MS**2)*DSQRT((S-S33)**2 -4*MS**2*S)
         ENDIF
         S42 = S3G2/2D0/(S32+MS**2)*(S -S32 -2*MS**2)
     +        +PRES42 *(Y -0.5D0)
         S43 = S3G3/2D0/(S33+MS**2)*(S -S33 -2*MS**2)
     +        +PRES43 *(Y -0.5D0)
         S4G2 = S42 -M2
         S4G3 = S43 -M2

         PRETG2 = DSQRT( (S-S4G2)**2 -4*MG**2*S )
         PRETG3 = 1
         IF((S-S4G3)**2 -4*MG**2*S.GT.0.D0)THEN
          PRETG3 = DSQRT( (S-S4G3)**2 -4*MG**2*S )
         ENDIF
         TG2 = -(S-S4G2)/2.D0 +(Z -0.5D0)*PRETG2
         TG3 = -(S-S4G3)/2.D0 +(Z -0.5D0)*PRETG3

         PREX2 = S42/2D0/(S42+MS**2)*DSQRT((S-S4G2)**2-4*MG**2*S)
         PREX3 = 1
         IF((S-S4G3)**2-4*MG**2*S.GT.0.D0)THEN
          PREX3 = S43/2D0/(S43+MS**2)*DSQRT((S-S4G3)**2-4*MG**2*S)
         ENDIF

         PREF2 = PRES42 *PRETG2 /PREX2 * PREFTZ
         PREF3 = PRES43 *PRETG3 /PREX3 * PREFTZ

         S4G02 = S4G0**2 + EPSG*MG**4
      END IF

      IF (IFLAVOR.EQ.1) HARD = PREF* (
     +     + DSGQGH(ALPHAS,S,TG,S4,MS,MG)
     +     + DSGQGD(ALPHAS,S,TG,S4,MS,MG,DEL,S4MAX)/(S4MAX -DEL) )

      IF (IFLAVOR.EQ.2) THEN
         HARD1 = 0D0
     +        + PREF0 *DSGGGH(ALPHAS,S,TG0,S40,MS,MG,EPSS,EPSG) *S4G02
         IF ((S.GT.4*MG**2).AND.(MG.GT.MS)) THEN
            HARD2 = 0D0
     +           + PREF0 *DSGGGS(ALPHAS,S,TG0,S40,MS,MG)
     +           - PREF1 *DSGGGS(ALPHAS,S,TG1,M2 ,MS,MG)
         ELSE
            HARD2 = 0D0
     +           + PREF0 *DSGGGS(ALPHAS,S,TG0,S40,MS,MG)
         END IF
         IF ((S.GT.4*MS**2).AND.(MS.GT.MG)) THEN
            HARD3 = 0D0
     +           + PREF2*DSGGGT(ALPHAS,S,TG2,S42,S32,MS,MG)
     +           - PREF3*DSGGGT(ALPHAS,S,TG3,S43,0D0,MS,MG)
         ELSE
            HARD3 = 0D0
     +           + PREF2*DSGGGT(ALPHAS,S,TG2,S42,S32,MS,MG)
         END IF
         HARD = HARD1 + HARD2 +HARD3
      END IF

      IF (IFLAVOR.EQ.3) THEN
         HARD1 =
     +        PREF0 *DSGQBH(ALPHAS,S,TG0,S40,MS,MG,1,EPSS,EPSG)*S4G02
         IF ((S.GT.4*MG**2).AND.(MG.GT.MS)) THEN
            HARD2 =
     +           + PREF0 *DSGQBS(ALPHAS,S,TG0,S40,MS,MG,1)
     +           - PREF1 *DSGQBS(ALPHAS,S,TG1,M2 ,MS,MG,1)
         ELSE
            HARD2 =
     +           + PREF0 *DSGQBS(ALPHAS,S,TG0,S40,MS,MG,1)
         END IF
         IF ((S.GT.4*MS**2).AND.(MS.GT.MG)) THEN
            HARD3 =
     +           + PREF2*DSGQBT(ALPHAS,S,TG2,S42,S32,MS,MG,1)
     +           - PREF3*DSGQBT(ALPHAS,S,TG3,S43,0D0,MS,MG,1)
         ELSE
            HARD3 =
     +           + PREF2*DSGQBT(ALPHAS,S,TG2,S42,S32,MS,MG,1)
         END IF
         HARD = HARD1 + HARD2 +HARD3
      END IF

      IF (IFLAVOR.EQ.4) THEN
         HARD1 =
     +        PREF0* DSGQBH(ALPHAS,S,TG0,S40,MS,MG,0,EPSS,EPSG)*S4G02
         HARD2 = 0.D0
         IF ((S.GT.4*MS**2).AND.(MS.GT.MG)) THEN
            HARD3 =
     +           + PREF2*DSGQBT(ALPHAS,S,TG2,S42,S32,MS,MG,0)
     +           - PREF3*DSGQBT(ALPHAS,S,TG3,S43,0D0,MS,MG,0)
         ELSE
            HARD3 =
     +           + PREF2*DSGQBT(ALPHAS,S,TG2,S42,S32,MS,MG,0)
         END IF
         HARD = HARD1 + HARD2 +HARD3
      END IF

      IF (IFLAVOR.EQ.5) THEN
         HARD1 =
     +        PREF0* DSGQQH(ALPHAS,S,TG0,S40,MS,MG,1,EPSS,EPSG) *S4G02
         HARD2 = 0.D0
         IF ((S.GT.4*MS**2).AND.(MS.GT.MG)) THEN
            HARD3 =
     +           + PREF2*DSGQQT(ALPHAS,S,TG2,S42,S32,MS,MG,1)
     +           - PREF3*DSGQQT(ALPHAS,S,TG3,S43,0D0,MS,MG,1)
         ELSE
            HARD3 =
     +           + PREF2*DSGQQT(ALPHAS,S,TG2,S42,S32,MS,MG,1)
         END IF
         HARD = HARD1 + HARD2 +HARD3
      END IF

      IF (IFLAVOR.EQ.6) THEN
         HARD1 =
     +        PREF0* DSGQQH(ALPHAS,S,TG0,S40,MS,MG,0,EPSS,EPSG) *S4G02
         HARD2 = 0.D0
         IF ((S.GT.4*MS**2).AND.(MS.GT.MG)) THEN
            HARD3 =
     +           + PREF2*DSGQQT(ALPHAS,S,TG2,S42,S32,MS,MG,0)
     +           - PREF3*DSGQQT(ALPHAS,S,TG3,S43,0D0,MS,MG,0)
         ELSE
            HARD3 =
     +           + PREF2*DSGQQT(ALPHAS,S,TG2,S42,S32,MS,MG,0)
         END IF
         HARD = HARD1 + HARD2 +HARD3
      END IF

      IF (.NOT.(DABS(HARD).LT.1.D40)) THEN
         HARD = 0.D0
         PRINT *,'NAN'
      END IF

      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION SOVISCA(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(1)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      XLAM = DSQRT((S -MG**2 -MS**2)**2 -4*MG**2*MS**2)
      TG = -S/2.D0*(1.D0 +(MG**2 -MS**2)/S -(2*X -1.D0)*XLAM/S)
      IF (IFLAVOR.EQ.1) SOVISCA = XLAM*DSGQG1(ALPHAS,S,TG,MS,MG,SCA)
      RETURN
      END

C ======================================================================

      DOUBLE PRECISION FUNCTION HARDSCA(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      IMPLICIT COMPLEX*16 (K)
      DIMENSION VAR(2)
      COMMON/CONST/S,ALPHAS,MS,MG,MT,SCA
      COMMON/DELTA/DEL
      COMMON/IFLAVOR/IFLAVOR
      X = VAR(1)
      Y = VAR(2)
      XLAM = DSQRT((S -MG**2 -MS**2)**2 -4*MG**2*MS**2)
      TG = -(S +MG**2 -MS**2 -DEL)/2.D0
     +     + (X -0.5D0) *DSQRT((S+MG**2-MS**2-DEL)**2 -4*MG**2*S)
      S4MAX = S +TG + MG**2 -MS**2 +MG**2*S/TG
      S4 = ( S4MAX -DEL)*Y +DEL
      PREF = S4MAX*XLAM
      IF (IFLAVOR.EQ.1) HARDSCA = PREF* (
     +     + DSGQG2(ALPHAS,S,TG,S4,MS,MG,DEL,S4MAX,SCA)/(S4MAX -DEL)
     +     + DSGQG3(ALPHAS,S,TG,S4,MS,MG,SCA) )
      IF (IFLAVOR.EQ.2) HARDSCA = PREF* (
     +     + DSGGG3(ALPHAS,S,TG,S4,MS,MG,SCA) )
      IF (IFLAVOR.EQ.3) HARDSCA = PREF* (
     +     + DSGQB3(ALPHAS,S,TG,S4,MS,MG,SCA,1) )
      IF (IFLAVOR.EQ.4) HARDSCA = PREF* (
     +     + DSGQB3(ALPHAS,S,TG,S4,MS,MG,SCA,0) )
      IF (IFLAVOR.EQ.5) HARDSCA = PREF* (
     +     + DSGQQ3(ALPHAS,S,TG,S4,MS,MG,SCA,1) )
      IF (IFLAVOR.EQ.6) HARDSCA = PREF* (
     +     + DSGQQ3(ALPHAS,S,TG,S4,MS,MG,SCA,0) )

      RETURN
      END

C ======================================================================

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE GAUSS_BL (ISYTYP, VARFAC, LENF, NLNKF, LINKF, HIDEF,
     &   XN, YN, ZN, ISVOK, LVARF, *)
C=======================================================================

C   --*** GAUSS *** (DETOUR) Plot element symbols at gauss points
C   --   Written by Amy Gilkey, revised 03/14/88
C   --   D. P. Flanagan, 12/08/83
C   --
C   --GAUSS drives the user symbol interface for element variables.
C   --It processes each element by element block, computing scaling factors
C   --and element information, then calling the user symbol routine.
C   --Only elements of selected element blocks are drawn.
C   --
C   --Parameters:
C   --   ISYTYP - IN - the symbol type (as in MODTYP of /DETOPT/)
C   --   VARFAC - IN - the face variable values
C   --   LENF - IN - the cumulative face counts by element block
C   --   NLNKF - IN - the number of nodes per face
C   --   LINKF - IN - the connectivity for all faces
C   --   HIDEF(i) - IN - true iff face i is hidden (3D only)
C   --   XN, YN, ZN - IN - the nodal coordinates
C   --   ISVOK - IN - ISVOK(i) is true iff the gauss variables are defined
C   --      for element block i
C   --   * - return statement if the cancel function is active
C   --
C   --Common Variables:
C   --   Uses NELBLK of /DBNUMS/
C   --   Uses IS3DIM of /D3NUMS/
C   --   Uses VECSCL of /ETCOPT/
C   --   Uses DTW, VWSCL of /DEVDAT/

      PARAMETER (KSCHSZ=2)

      include 'dbnums.blk'
      COMMON /D3NUMS/ IS3DIM, NNPSUR, NUMNPF, LLNSET
      LOGICAL IS3DIM
      COMMON /ETCOPT/ VECSCL
      COMMON /DEVDAT/ DTW, VWSCL

      CHARACTER*(*) ISYTYP
      INTEGER LENF(0:NELBLK)
      INTEGER NLNKF(NELBLK)
      INTEGER LINKF(*)
      REAL VARFAC(LVARF,4)
      LOGICAL HIDEF(*)
      REAL XN(*), YN(*), ZN(*)
      LOGICAL ISVOK(NELBLK)

      LOGICAL GRABRT
      LOGICAL PLTSTT, PLTGTT, LDUM
      REAL XGAUSS(4), YGAUSS(4), ZGAUSS(4)

      SSCL = .25 * VECSCL * DTW

C   --Set symbol size, must be reset before exit
      LDUM = PLTGTT (KSCHSZ, SYMSZ)
      LDUM = PLTSTT (KSCHSZ, VECSCL*SYMSZ)

      DO 120 IELB = 1, NELBLK
         IF (ISVOK(IELB) .AND. (NLNKF(IELB) .EQ. 4)) THEN

C         --Set the symbol color

            CALL GRCOLR (IELB)

            DO 110 IFAC = LENF(IELB-1)+1, LENF(IELB)
               IF (IS3DIM) THEN
                  IF (HIDEF(IFAC)) GOTO 110
               END IF

C            --Compute gauss coordinates

               IXL = IDBLNK (IELB, IFAC, LENF, NLNKF)
               CALL GAUSSF (IS3DIM, NLNKF(IELB), LINKF(IXL), XN, YN, ZN,
     &            XGAUSS, YGAUSS, ZGAUSS)

C            --Call symbol routine each gauss point

               IF (GRABRT ()) GOTO 130
               DO 100 I = 1, 4
                  CALL USRSYM (ISYTYP, IS3DIM,
     &               XGAUSS(I), YGAUSS(I), ZGAUSS(I), VARFAC(IFAC,I),
     &               SSCL)
  100          CONTINUE
  110       CONTINUE

            CALL PLTFLU
         END IF
  120 CONTINUE

C   --Reset symbol size
      LDUM = PLTSTT (KSCHSZ, SYMSZ)

      RETURN

  130 CONTINUE
C   --Reset symbol size
      LDUM = PLTSTT (KSCHSZ, SYMSZ)

      RETURN 1
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
!
! JACOBI
! Jacobi diagonalizer with sorted output.  Same calling sequence as
! EISPACK routine, but must specify nrot!
!
      SUBROUTINE jacobi (a, n, np, d, v, nrot)
      IMPLICIT CHARACTER (A-Z)
!
      INTEGER n, np, nrot
      DOUBLEPRECISION a (np, n)
      DOUBLEPRECISION d (n)
      DOUBLEPRECISION v (np, n)
!
      DOUBLEPRECISION onorm, dnorm
      DOUBLEPRECISION b, dma, q, t, c, s
      DOUBLEPRECISION atemp, vtemp, dtemp
      INTEGER i, j, k, l
!
      DO 10000 j = 1, n
        DO 10010 i = 1, n
          v (i, j) = 0.0D0
10010     CONTINUE
        v (j, j) = 1.0D0
        d (j) = a (j, j)
10000   CONTINUE
!
      DO 20000 l = 1, nrot
        dnorm = 0.0D0
        onorm = 0.0D0
        DO 20100 j = 1, n
          dnorm = dnorm + ABS (d (j))
          DO 20110 i = 1, j - 1
            onorm = onorm + ABS (a (i, j))
20110       CONTINUE
20100     CONTINUE
        IF (onorm / dnorm .LE. 0.0D0) GOTO 19999
        DO 21000 j = 2, n
          DO 21010 i = 1, j - 1
            b = a (i, j)
            IF (ABS (b) .GT. 0.0D0) THEN
              dma = d (j) - d (i)
              IF (ABS (dma) + ABS (b) .LE. ABS (dma)) THEN
                t = b / dma
              ELSE
                q = 0.5D0 * dma / b
                t = SIGN (1.0D0 / (ABS (q) + SQRT (1.0D0 + q * q)), q)
              ENDIF
              c = 1.0D0 / SQRT (t * t + 1.0D0)
              s = t * c
              a (i, j) = 0.0D0
              DO 21110 k = 1, i - 1
                atemp    = c * a (k, i) - s * a (k, j)
                a (k, j) = s * a (k, i) + c * a (k, j)
                a (k, i) = atemp
21110           CONTINUE
              DO 21120 k = i + 1, j - 1
                atemp    = c * a (i, k) - s * a (k, j)
                a (k, j) = s * a (i, k) + c * a (k, j)
                a (i, k) = atemp
21120           CONTINUE
              DO 21130 k = j + 1, n
                atemp    = c * a (i, k) - s * a (j, k)
                a (j, k) = s * a (i, k) + c * a (j, k)
                a (i, k) = atemp
21130           CONTINUE
              DO 21140 k = 1, n
                vtemp    = c * v (k, i) - s * v (k, j)
                v (k, j) = s * v (k, i) + c * v (k, j)
                v (k, i) = vtemp
21140           CONTINUE
              dtemp = c * c * d (i) + s * s * d (j) -&
                      2.0D0 * c * s * b
              d (j) = s * s * d (i) + c * c * d (j) +&
                      2.0D0 * c * s * b
              d (i) = dtemp
              ENDIF
21010       CONTINUE
21000     CONTINUE
20000   CONTINUE
19999 CONTINUE
      nrot = l
!
      DO 30000 j = 1, n - 1
        k = j
        dtemp = d (k)
        DO 30100 i = j + 1, n
          IF (d (i) .LT. dtemp) THEN
            k = i
            dtemp = d (k)
            ENDIF
30100     CONTINUE
        IF (k .GT. j) THEN
          d (k) = d (j)
          d (j) = dtemp
          DO 30200 i = 1, n
            dtemp    = v (i, k)
            v (i, k) = v (i, j)
            v (i, j) = dtemp
30200       CONTINUE
          ENDIF
30000   CONTINUE
!
RETURN
END subroutine jacobi

max_stars_repo_name: donshen/packmol

max_issues_repo_name: donshen/packmol

max_forks_repo_name: donshen/packmol