      MODULE CFAST
!=======================================================================
!     LATEST CORRECTION BY
!
!     PURPOSE
!            Module defining the the Faster structures.
!
!     PROGRAMMED BY:  bjorn.melhus@akersolutions.com
!     CREATED......:  10.05.2021
!=======================================================================
      CHARACTER*64      :: CIDENT='<description>' ! User defined ID of the analyses
      INTEGER           :: NPARA=0            ! Number of parameters, i.e. number of openfast runs
      LOGICAL           :: LANAL=.FALSE.      ! Logical for performing the Openfast analyses
      LOGICAL           :: LSIMUL=.TRUE.      ! Logical for running the analyses, default=.true.
                                              !  (LSIMUL=.FALSE. and LANAL=.TRUE. will generate the files, but omit the simulation)
      INTEGER           :: MOOR_LINES=0       ! Number of mooring lines
      INTEGER           :: NSEGMENTS=0        ! Number of segments on each mooring line
      INTEGER           :: MAX_MASTER=3       ! Number of masters (will be increased if needed)
      INTEGER           :: MAX_SLAVE=4        ! Number of slaves (will be increased if needed)
      INTEGER           :: MAX_PAR=30         ! Number of assigned parameters = number of openfast analyses (will be increased if needed)
      INTEGER           :: MAX_VAR=4          ! Number of variables variables; WaveTp, WaveTP, etc.
      INTEGER           :: NLINES_OF_OUTUT_FILE
      INTEGER,PARAMETER :: MAX_CHARACTERS=2500     ! Max characters in one line of the output file
      INTEGER,PARAMETER :: NUMCOMM=6               ! Initial number of comment lines of the output file
      LOGICAL           :: CLEAN_UP_MESS=.FALSE.   ! Logical for clean up, i.e. delete all slave files etc.
!
      TYPE TMASTER
        INTEGER                    :: NMASTER= 0            ! Number of master files
        INTEGER                    :: NSLAVE = 0            ! Number of slaves for given master
        INTEGER , allocatable      :: MSLAVES(:)            ! Array of slave numbers
        CHARACTER*128, allocatable :: CMASTER(:)            ! Array of masters MAX_MASTER
        CHARACTER*128              :: CMAESTRO(2)           ! Meastros for Openfast and Turbsim
        CHARACTER*128, allocatable :: CSLAVE(:,:)           ! Array of slaves for given master
      END TYPE TMASTER
!
      TYPE TSLAVE
        INTEGER                    :: NSLAVE = 0            ! Number of slaves
        INTEGER                    :: ITYPE=1               ! ITYPE=1 -> Openfast,  ITYPE=2 -> Turbsim
        INTEGER                    :: NSWAP=0               ! Number of files where strings are to be swapped
        INTEGER,allocatable        :: NVAR(:)               ! Number of variables for given slave
        CHARACTER*128,allocatable  :: CSLAVE(:)             ! Slave id; MAXSLAVE
        CHARACTER*128,allocatable  :: CSWAP_STR_FROM(:)     ! String to be swappeed  ; MAX_SLAVE
        CHARACTER*128,allocatable  :: CSWAP_STR_TO(:,:)     ! String to be swappeed to in npar files ; NPAR*MAXSLAVE
        CHARACTER*128,allocatable  :: CSWAP(:)              ! Files to be swapped from  ; MAX_SLAVE
        CHARACTER*14, allocatable  :: CPAR(:,:,:)           ! parameters ; MAX_PAR*MAX_VAR*MAX_SLAVE
        CHARACTER*14, allocatable  :: CVAR(:,:)             ! variables ; MAX_VAR*MAX_SLAVE
        LOGICAL , allocatable      :: CONST(:,:)            ! Constant for all runs?  MAX_VAR*MAX_SLAVE
      END TYPE TSLAVE
!
      TYPE TRES
         CHARACTER*14, allocatable  :: CHAN(:)              ! Channel ID's
         CHARACTER*14, allocatable  :: CUNI(:)              ! Channel units
         INTEGER                    :: NROW                 ! Number of rows = Number of time instances
         INTEGER                    :: NCOL                 ! Number of columns = Number of channels
       END TYPE TRES
!
      TYPE (TMASTER) ::  FAST
      TYPE (TSLAVE)  ::  SLAVE
      TYPE (TRES)    ::  RES
!
     END MODULE


text:
      subroutine readops

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine reads data from the HRU/subbasin management input file
!!    (.mgt). This file contains data related to management practices used in
!!    the HRU/subbasin.

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name       |units            |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!                                 |urban.dat
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name            |units          |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!                                    |daily
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    eof         |none          |end of file flag (=-1 if eof, else = 0)
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    day         |none          |day operation occurs
!!    mgt_op      |none          |operation code number
!!                               |0 end of rotation year
!!                               |1 plant/beginning of growing season
!!    mgt1i       |none          |first management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt2i       |none          |second management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt3i       |none          |third management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt4        |none          |fourth management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt5        |none          |fifth management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt6        |none          |sixth management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt7        |none          |seventh management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt8        |none          |eighth management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    mgt9        |none          |ninth management parameter out of .mgt
!!                               |file (definition changes depending on
!!                               |mgt_op)
!!    titldum     |NA            |title line from input dataset
!!    vfscon(:)   |none          |Fraction of the total runoff from the entire field
!!                               |entering the most concentrated 10% of the VFS.
!!    vfsratio(:) |none          |Field area/VFS area ratio
!!    vfsch(:)    |none          |Fraction of flow entering the most concentrated 10% of the VFS.
!!	                           |which is fully channelized
!!	grwat_n(:)      |none      |Mannings's n for grassed waterway
!!	grwat_i(:)      |none      |Flag for the simulation of grass waterways
!!                                 | gwat_i = 0 inactive
!!                                 |        = 1 active
!!	grwat_l(:)      |km	   |Length of grass Waterway
!!	grwat_w(:)      |m         |Average width of grassed waterway
!!	grwat_d(:)      |m         |Depth of grassed waterway from top of bank
!!                                 |  to bottom
!!	grwat_s(:)      |m         |Average slope of grassed waterway channel
!!	grwat_spcon(:)  |none      |Linear parameter for calculating sediment
!!                                 |  in grassed waterways
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    SWAT: Jdt

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      character (len=80) :: titldum
	integer :: eof
	integer :: mon, day, mgt_op, mgt2i, mgt1i
	real :: mgt3,mgt4,mgt5,mgt6,mgt7,mgt8,mgt9

	read (111,5000,end=999) titldum

!!      read scheduled operations

          iops = 0

          do
          mon = 0
          day = 0
          year = 0.
          mgt_op = 0
          mgt1i = 0
          mgt2i = 0
          mgt3 = 0.
          mgt4 = 0.
          mgt5 = 0.
          mgt6 = 0
          mgt7 = 0.
          mgt8 = 0.
		mgt9 = 0.

          read (111,5200,iostat=eof) mon, day, iyear, mgt_op, mgt1i,
     & mgt2i, mgt3, mgt4, mgt5, mgt6, mgt7, mgt8, mgt9
	    if (eof < 0) exit
          iops = iops + 1
          iopday(iops,ihru) = Jdt (ndays,day,mon)
          iopyr(iops,ihru) = iyear
          mgt_ops(iops,ihru) = mgt_op

          select case (mgt_op)

          case (0)  !! end of operations file

          case (1)  !! terracing
             terr_p(iops,ihru) = mgt4
             terr_cn(iops,ihru) = mgt5
             terr_sl(iops,ihru) = mgt6

          case (2)  !! tile drainage
             drain_d(iops,ihru) = mgt4
             drain_t(iops,ihru) = mgt5
             drain_g(iops,ihru) = mgt6
             if (mgt7 < 1.e-6) mgt7 = 6000.0
             drain_idep(iops,ihru) = mgt7

          case (3)  !! contouring
             cont_cn(iops,ihru) = mgt4
             cont_p(iops,ihru) =  mgt5

          case (4)  !! filter
		     filter_i(iops,ihru) = mgt1i  !! on off flag
		     filter_ratio(iops,ihru) = mgt3
		     filter_con(iops,ihru) = mgt4
		     filter_ch(iops,ihru) = mgt5

          case (5)  !! strip cropping
             strip_n(iops,ihru) = mgt4
             strip_cn(iops,ihru) = mgt5
             strip_c(iops,ihru) = mgt6
             strip_p(iops,ihru) = mgt7

          case (6)  !! fire
             fire_cn(iops,ihru) = mgt4

          case (7)  !! grass waterways
             gwati(iops,ihru) = mgt1i
             gwatn(iops,ihru) = mgt3
             gwatspcon(iops,ihru) = mgt4
             gwatd(iops,ihru) = mgt5
             gwatw(iops,ihru) = mgt6
             gwatl(iops,ihru) = mgt7
             gwats(iops,ihru) = mgt8
!! Set defaults
!! Mannings via Fiener, 2006
	if (gwatn(iops,ihru) <=0.) gwatn(iops,ihru) = 0.35
!! length based on one side of a square HRU
	if (gwatl(iops,ihru) <=0.) gwatl(iops,ihru) = hru_km(ihru)**.5
!! default to a workable depth based on with and 8:1 sideslopes
	if (gwatd(iops,ihru) <= 0.) then
		gwatd(iops,ihru) = 3. / 64. * gwatw(iops,ihru)
	end if
!! Default to 3/4 HRU slope
	if (gwats(iops,ihru) <=0.) gwats(iops,ihru) = hru_slp(ihru)*.75
!! default sed coeff to 0.005
	if (gwatspcon(iops,ihru)<= 0.) gwatspcon(iops,ihru) = 0.005

          case (8) !! plant parameter update
             cropno_upd(iops,ihru) = mgt1i
             hi_upd(iops,ihru) = mgt4
             laimx_upd(iops,ihru) = mgt5

	!!	case (9) !! Residue Managment  mjw
	!! Force residue to a minimum value regardless of tillage.		mjw
	!!	so_res_flag(iops,ihru) = mgt1i !!mjw
	!!	so_res(iops,ihru) = mgt4 !!mjw

		case (10) !! Generic Conservation Practice  mjw
	!! Get user defined removal eff and use these		mjw
		ro_bmp_flag(iops,ihru) = mgt1i !! Flag to turn on or off user BMP
		ro_bmp_sed(iops,ihru) = mgt3 !! Sediment
		ro_bmp_pp(iops,ihru) = mgt4 !! Particulate P
		ro_bmp_sp(iops,ihru) = mgt5 !! Soluble P
		ro_bmp_pn(iops,ihru) = mgt7 !! Particulate N1123
		ro_bmp_sn(iops,ihru) = mgt8 !! Soluble N
		ro_bmp_bac(iops,ihru) = mgt9 !! Bacteria
          end select
          end do

      close (111)

  999 return
 5000 format (a)
 5200 format (1x,i2,1x,i2,5x,i4,1x,i2,1x,i4,1x,i3,1x,f6.2,1x,f12.5,1x,
     &        f6.2,1x,f11.5,1x,f8.2,1x,f6.2,1x,2f5.2)
      end

ext: f


text:
cc      MODULE FEMUTIL
cc      CONTAINS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C11111111122222222223333333333444444444455555555556666666666777777777777
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C         --G E N E R A L  F E M  S U B R O U T I N E S--              C
C                                                                      C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C STRAIN DISPLACEMENT INTERPOLATORS                                    C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM9:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 9-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE STDM9(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,XBAR)

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)

      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(2,NNE),AUX1(2,NNE),
     1          XJ(2,2),XJI(2,2)

C     Initialize arrays

      XBAR=ONE

C     Shape functions derivatives w.r.t natural coordinates

      CALL SFDER9(IELT,NDOFEL,NNE,R,S,P)

C     Computes the Jacobian operator and its determinant at point (r,s)

      CALL JACOPER(NNE,XJ,XX,P)
      DDET=XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1)
C     Computes the inverse of the Jacobiam operator
C
      CALL JACINVE(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      AUX1=MATMUL(XJI,P)
C
C     Assembles B matrix.
C
      DO I=1,NNE
        II=2*(I-1)+1
        B(1,II)=AUX1(1,I)
        B(2,II+1)=AUX1(2,I)

        B(3,II)=AUX1(2,I)
        B(3,II+1)=AUX1(1,I)

        B(1,II+1)=0.0D0
        B(2,II  )=0.0D0

      END DO
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM8:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 8-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE STDM8(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,XBAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)
C
      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(2,NNE),AUX1(2,NNE),
     1          XJ(2,2),XJI(2,2)
C
C     Initialize arrays
C
      XBAR=ONE
C
C     Shape functions derivatives w.r.t natural coordinates
C
      CALL SFDER8(IELT,NDOFEL,NNE,R,S,P)
C
C     Computes the Jacobian operator and its determinant at point (r,s)
C
      CALL JACOPER(NNE,XJ,XX,P)
      DDET=XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1)
C
C     Computes the inverse of the Jacobiam operator
C
      CALL JACINVE(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      AUX1=MATMUL(XJI,P)
C
C     Assembles B matrix for a
C     Cosserat element.
C
      DO I=1,NNE
        II=2*(I-1)+1
        B(1,II)=AUX1(1,I)
        B(2,II+1)=AUX1(2,I)
        B(3,II)=AUX1(2,I)
        B(3,II+1)=AUX1(1,I)

        B(1,II+1)=0.0D0
        B(2,II  )=0.0D0

      END DO
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM4:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 4-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE STDM4(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,XBAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)
C
      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(2,NNE),AUX1(2,NNE),
     1          XJ(2,2),XJI(2,2)
C
C     Initialize arrays
C
      XBAR=ONE
C
C     Shape functions derivatives w.r.t natural coordinates
C
      CALL SFDER4(IELT,NDOFEL,NNE,R,S,P)
C
C     Computes the Jacobian operator and its determinant at point (r,s)
C
      CALL JACOPER(NNE,XJ,XX,P)
      DDET=XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1)
C
C     Computes the inverse of the Jacobian operator
C
      CALL JACINVE(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      AUX1=MATMUL(XJI,P)
C
C     Assembles B matrix.
C
      DO I=1,NNE
        II=2*(I-1)+1
        B(1,II)=AUX1(1,I)
        B(2,II+1)=AUX1(2,I)
        B(3,II)=AUX1(2,I)
        B(3,II+1)=AUX1(1,I)

        B(1,II+1)=0.0D0
        B(2,II  )=0.0D0
      END DO
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM6:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 6-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE STDM6(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,XBAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)
C
      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(2,NNE),AUX1(2,NNE),
     1          XJ(2,2),XJI(2,2)
C
C     Initialize arrays
C
      XBAR=ONE
C
C     Shape functions derivatives w.r.t natural coordinates
C
      CALL SFDER6(IELT,NDOFEL,NNE,R,S,P)
C
C     Computes the Jacobian operator and its determinant at point (r,s)
C
      CALL JACOPER(NNE,XJ,XX,P)
      DDET=XJ(1,1)*XJ(2,2)-XJ(1,2)*XJ(2,1)
C
C     Computes the inverse of the Jacobiam operator
C
      CALL JACINVE(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      AUX1=MATMUL(XJI,P)
C
C     Assembles B matrix.
C
      DO I=1,NNE
        II=2*(I-1)+1
        B(1,II)=AUX1(1,I)
        B(2,II+1)=AUX1(2,I)
        B(3,II)=AUX1(2,I)
        B(3,II+1)=AUX1(1,I)

        B(1,II+1)=0.0D0
        B(2,II  )=0.0D0

      END DO
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM8:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 8-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE STDM3D8(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,T,
     1                   XBAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)
C
      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(3,NNE),AUX1(3,NNE),
     1          XJ(3,3),XJI(3,3)
C
C     Initialize arrays
C
      XBAR=ONE
C
C     Shape functions derivatives w.r.t natural coordinates
C
ccccccccESTA RUTINA TAMBIÉN TIENE PROBLEMAS
C     CALL SFDER3D8(IELT,NDOFEL,NNE,R,S,T,P)
      CALL SFDER3D8(NDOFEL,NNE,R,S,T,P)
C
C     Computes the Jacobian operator and its determinant at point (r,s,t)
C
      CALL JACOPER3D(NNE,XJ,XX,P)
      DD1= XJ(1,1)*(XJ(2,2)*XJ(3,3)-XJ(2,3)*XJ(3,2))
      DD2=-XJ(1,2)*(XJ(2,1)*XJ(3,3)-XJ(2,3)*XJ(3,1))
      DD3= XJ(1,3)*(XJ(2,1)*XJ(3,2)-XJ(2,2)*XJ(3,1))
      DDET= DD1+DD2+DD3
C
C     Computes the inverse of the Jacobiam operator
C
      CALL JACINVE3D(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      AUX1=MATMUL(XJI,P)
C
C     Assembles B matrix for a
C     three dimensional element.
C
      DO I=1,NNE
        II=3*(I-1)+1
        B(1,II  )=AUX1(1,I)
        B(1,II+1)=0.0D0
        B(1,II+2)=0.0D0
C
        B(2,II  )=0.0D0
        B(2,II+1)=AUX1(2,I)
        B(2,II+2)=0.0D0
C
        B(3,II  )=0.0D0
        B(3,II+1)=0.0D0
        B(3,II+2)=AUX1(3,I)
C
        B(4,II  )=AUX1(2,I)
        B(4,II+1)=AUX1(1,I)
        B(4,II+2)=0.0D0
C
        B(5,II  )=AUX1(3,I)
        B(5,II+1)=0.0D0
        B(5,II+2)=AUX1(1,I)
C
        B(6,II  )=0.0D0
        B(6,II+1)=AUX1(3,I)
        B(6,II+2)=AUX1(2,I)
C
      END DO
C
      RETURN
C
      END
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE STDM8:GENERATES THE STRAIN-DISPLACEMENT MATRIX B      C
C     AND JACOBIAN DETERMINANT DDET AT THE POINT r ,s                  C
C                                                 i  j                 C
C     FOR AN 8-NODED 2D ELEMENT-PLANE STRAIN                           C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE STDM3D10(NNE,NDOFEL,NTENS,XX,MCRD,B,DDET,R,S,T,
     1                   XBAR)
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (ONE=1.D0,TWO=2.D0,THREE=3.D0,HALF=0.5D0)
C
      DIMENSION XX(MCRD,NNE),B(NTENS,NDOFEL),P(3,NNE),AUX1(3,NNE),
     1          XJ(3,3),XJI(3,3)
C
C     Initialize arrays
C
      XBAR=ONE
C
C     Shape functions derivatives w.r.t natural coordinates
C
      CALL SFDER3D10(NDOFEL,NNE,R,S,T,P)
C
C     Computes the Jacobian operator and its determinant at point (r,s,t)
C
      CALL JACOPER3D(NNE,XJ,XX,P)
      DD1= XJ(1,1)*(XJ(2,2)*XJ(3,3)-XJ(2,3)*XJ(3,2))
      DD2=-XJ(1,2)*(XJ(2,1)*XJ(3,3)-XJ(2,3)*XJ(3,1))
      DD3= XJ(1,3)*(XJ(2,1)*XJ(3,2)-XJ(2,2)*XJ(3,1))
      DDET= DD1+DD2+DD3
C
C     Computes the inverse of the Jacobiam operator
C
      CALL JACINVE3D(XJ,DDET,XJI)
C
C     Jacobian Inverse times Shape Functions derivatives w.r.t natural coordinates
C     to produce shape function derivatives w.r.t x,y coordinates.
C
      CALL MMULT(XJI,3,3,P,3,NNE,AUX1)
C
C     Assembles B matrix for a
C     three dimensional element.
C
      DO I=1,NNE
        II=3*(I-1)+1
        B(1,II  )=AUX1(1,I)
        B(1,II+1)=0.0D0
        B(1,II+2)=0.0D0
C
        B(2,II  )=0.0D0
        B(2,II+1)=AUX1(2,I)
        B(2,II+2)=0.0D0
C
        B(3,II  )=0.0D0
        B(3,II+1)=0.0D0
        B(3,II+2)=AUX1(3,I)
C
        B(4,II  )=AUX1(2,I)
        B(4,II+1)=AUX1(1,I)
        B(4,II+2)=0.0D0
C
        B(5,II  )=AUX1(3,I)
        B(5,II+1)=0.0D0
        B(5,II+2)=AUX1(1,I)
C
        B(6,II  )=0.0D0
        B(6,II+1)=AUX1(3,I)
        B(6,II+2)=AUX1(2,I)
C
      END DO
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER3D8:GENERATES THE SHAPE FUNCTION DERIVATIVES     C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s , t                  C
C                                             i  j   k                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C     SUBROUTINE SFDER3D8(IELT,NDOFEL,NNE,R,S,T,P)
      SUBROUTINE SFDER3D8(NDOFEL,NNE,R,S,T,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0 , QUARTER=0.25D0 , EIGTH =0.125D0)
C
      DIMENSION P(3,NNE)
C

      RP = ONE+R
      SP = ONE+S
      TP = ONE+T
      RM = ONE-R
      SM = ONE-S
      TM = ONE-T

C
C     3D LINEAR 8-NODED ELEMENT
C
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,1)= EIGTH*SP*TP
      P(1,2)=-EIGTH*SP*TP
      P(1,3)=-EIGTH*SM*TP
      P(1,4)= EIGTH*SM*TP
      P(1,5)= EIGTH*SP*TM
      P(1,6)=-EIGTH*SP*TM
      P(1,7)=-EIGTH*SM*TM
      P(1,8)= EIGTH*SM*TM
C
C     w.r.t.s
C
      P(2,1)= EIGTH*RP*TP
      P(2,2)= EIGTH*RM*TP
      P(2,3)=-EIGTH*RM*TP
      P(2,4)=-EIGTH*RP*TP
      P(2,5)= EIGTH*RP*TM
      P(2,6)= EIGTH*RM*TM
      P(2,7)=-EIGTH*RM*TM
      P(2,8)=-EIGTH*RP*TM
C
C     w.r.t.t
C
      P(3,1)= EIGTH*RP*SP
      P(3,2)= EIGTH*RM*SP
      P(3,3)= EIGTH*RM*SM
      P(3,4)= EIGTH*RP*SM
      P(3,5)=-EIGTH*RP*SP
      P(3,6)=-EIGTH*RM*SP
      P(3,7)=-EIGTH*RM*SM
      P(3,8)=-EIGTH*RP*SM

C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER3D8:GENERATES THE SHAPE FUNCTION DERIVATIVES     C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s , t                  C
C                                             i  j   k                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SFDER3D10(NDOFEL,NNE,R,S,T,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0 , QUARTER=0.25D0 , EIGTH =0.125D0)
C
      DIMENSION P(3,NNE)
C

      RP = ONE+R
      SP = ONE+S
      TP = ONE+T
      RM = ONE-R
      SM = ONE-S
      TM = ONE-T

C
C     3D LINEAR 8-NODED ELEMENT
C
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,1)= EIGTH*SP*TP
      P(1,2)=-EIGTH*SP*TP
      P(1,3)=-EIGTH*SM*TP
      P(1,4)= EIGTH*SM*TP
      P(1,5)= EIGTH*SP*TM
      P(1,6)=-EIGTH*SP*TM
      P(1,7)=-EIGTH*SM*TM
      P(1,8)= EIGTH*SM*TM
C
C     w.r.t.s
C
      P(2,1)= EIGTH*RP*TP
      P(2,2)= EIGTH*RM*TP
      P(2,3)=-EIGTH*RM*TP
      P(2,4)=-EIGTH*RP*TP
      P(2,5)= EIGTH*RP*TM
      P(2,6)= EIGTH*RM*TM
      P(2,7)=-EIGTH*RM*TM
      P(2,8)=-EIGTH*RP*TM
C
C     w.r.t.t
C
      P(3,1)= EIGTH*RP*SP
      P(3,2)= EIGTH*RM*SP
      P(3,3)= EIGTH*RM*SM
      P(3,4)= EIGTH*RP*SM
      P(3,5)=-EIGTH*RP*SP
      P(3,6)=-EIGTH*RM*SP
      P(3,7)=-EIGTH*RM*SM
      P(3,8)=-EIGTH*RP*SM

C
      RETURN
C
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE KJACOPER                                                C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE JACOPER3D(NNE,XJA,XCORD,RDER)

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER(ZERO=0.0D0)

      DIMENSION XJA(3,3),XCORD(3,NNE),RDER(3,NNE)
C
      DUM=ZERO
      DO K=1,3
        DO J=1,3
          DO I=1,NNE
            DUM=DUM+RDER(K,I)*XCORD(J,I)
          END DO
          XJA(K,J)=DUM
          DUM=ZERO
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      3X3 Jacobian inverse                                            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE JACINVE3D(XJA,DD,XJAI)

      IMPLICIT REAL*8 (A-H,O-Z)

      DIMENSION XJA(3,3),XJAI(3,3),XADJ(3,3),XADJT(3,3)
C
      XADJT(1 , 1) =  XJA(2,2)*XJA(3,3)-XJA(2,3)*XJA(3,2)
      XADJT(1 , 2) =  XJA(2,1)*XJA(3,3)-XJA(2,3)*XJA(3,1)
      XADJT(1 , 3) =  XJA(2,1)*XJA(3,2)-XJA(2,2)*XJA(3,1)
C
      XADJT(2 , 1) =  XJA(1,2)*XJA(3,3)-XJA(1,3)*XJA(3,2)
      XADJT(2 , 2) =  XJA(1,1)*XJA(3,3)-XJA(1,3)*XJA(3,1)
      XADJT(2 , 3) =  XJA(1,1)*XJA(3,2)-XJA(1,2)*XJA(3,1)
C
      XADJT(3 , 1) =  XJA(1,2)*XJA(2,3)-XJA(1,3)*XJA(2,2)
      XADJT(3 , 2) =  XJA(1,1)*XJA(2,3)-XJA(1,3)*XJA(2,1)
      XADJT(3 , 3) =  XJA(1,1)*XJA(2,2)-XJA(1,2)*XJA(2,1)
C
      CALL MTRAN(XADJT,3,3,XADJ)
C
      DO J=1,3
        DO K=1,3
          COFA=((-1)**(J+K))*XADJ(J,K)
          XJAI(J,K)=COFA/DD
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE AMASS3D8(MASS,NDOFEL,RII,SII)                        C
C      Computes the mass matrix for an 8-noded element                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS3D8(AMASS,NDOFEL,RII,SII,TII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=3)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(8),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL CLEAR(SF,NTENS,NDOFEL)
      CALL CLEAR(SFT,NDOFEL,NTENS)
      CALL SHAPEF3D8(SN,RII,SII,TII)
C
      DO I=1,8
        KK=3*(I-1)+1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(3,KK+2)=SN(I)
C
      END DO

      CALL MTRAN(SF,NTENS,NDOFEL,SFT)
      CALL MMULT(SFT,NDOFEL,NTENS,SF,NTENS,NDOFEL,AUX1)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=AUX1(I,J)
        END DO
      END DO
C
      RETURN

      END
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE AMASS3D8(MASS,NDOFEL,RII,SII)                        C
C      Computes the mass matrix for an 8-noded element                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS3D10(AMASS,NDOFEL,RII,SII,TII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=3)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(8),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL CLEAR(SF,NTENS,NDOFEL)
      CALL CLEAR(SFT,NDOFEL,NTENS)
      CALL SHAPEF3D8(SN,RII,SII,TII)
C
      DO I=1,8
        KK=3*(I-1)+1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(3,KK+2)=SN(I)
C
      END DO

      CALL MTRAN(SF,NTENS,NDOFEL,SFT)
      CALL MMULT(SFT,NDOFEL,NTENS,SF,NTENS,NDOFEL,AUX1)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=AUX1(I,J)
        END DO
      END DO
C
      RETURN

      END

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF8(SN,R,S)                                      C
C      8-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF3D8(SN,R,S,T)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,EIGTH =0.125D0)
C
      DIMENSION SN(8)
C
      CALL CLEARV(SN,8)
C
      RP =ONE+R
      RM =ONE-R
      SP =ONE+S
      SM =ONE-S
      TP =ONE+T
      TM =ONE-T
C
      SN(1)= EIGTH*RP*SP*TP
      SN(2)= EIGTH*RM*SP*TP
      SN(3)= EIGTH*RM*SM*TP
      SN(4)= EIGTH*RP*SM*TP
      SN(5)= EIGTH*RP*SP*TM
      SN(6)= EIGTH*RM*SP*TM
      SN(7)= EIGTH*RM*SM*TM
      SN(8)= EIGTH*RP*SM*TM
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF8(SN,R,S)                                      C
C      8-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF3D10(SN,R,S,T)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,EIGTH =0.125D0)
C
      DIMENSION SN(8)
C
      CALL CLEARV(SN,8)
C
      RP =ONE+R
      RM =ONE-R
      SP =ONE+S
      SM =ONE-S
      TP =ONE+T
      TM =ONE-T
C
      SN(1)= EIGTH*RP*SP*TP
      SN(2)= EIGTH*RM*SP*TP
      SN(3)= EIGTH*RM*SM*TP
      SN(4)= EIGTH*RP*SM*TP
      SN(5)= EIGTH*RP*SP*TM
      SN(6)= EIGTH*RM*SP*TM
      SN(7)= EIGTH*RM*SM*TM
      SN(8)= EIGTH*RP*SM*TM
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SHAPE FUNCTIONS AND THEIR DERIVATIVES                            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER9:GENERATES THE SHAPE FUNCTION DERIVATIVES       C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s                      C
C                                             i  j                     C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SFDER9(IELT,NDOFEL,NNE,R,S,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0,TWO=2.D0, HALF=0.5D0,QUARTER=0.25D0,FOUR=4.D0)
C
      DIMENSION P(2,NNE)
C
      RP= ONE+R
      SP= ONE+S
      RM= ONE-R
      SM= ONE-S
      RMS=ONE-R**TWO
      SMS=ONE-S**TWO
C
C     9-NODED ELEMENT
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,9)=-TWO*R*SMS
      P(1,8)=-HALF*SMS-HALF*P(1,9)
      P(1,7)=-R*SP-HALF*P(1,9)
      P(1,6)= HALF*SMS-HALF*P(1,9)
      P(1,5)=-R*SM-HALF*P(1,9)
      P(1,4)=-QUARTER*SP-HALF*P(1,7)-HALF*P(1,8)-QUARTER*P(1,9)
      P(1,3)= QUARTER*SP-HALF*P(1,7)-HALF*P(1,6)-QUARTER*P(1,9)
      P(1,2)= QUARTER*SM-HALF*P(1,5)-HALF*P(1,6)-QUARTER*P(1,9)
      P(1,1)=-QUARTER*SM-HALF*P(1,8)-HALF*P(1,5)-QUARTER*P(1,9)
C
C        w.r.t.s
C
      P(2,9)=-TWO*S*RMS
      P(2,8)=-S*RM-HALF*P(2,9)
      P(2,7)= HALF*RMS-HALF*P(2,9)
      P(2,6)=-S*RP-HALF*P(2,9)
      P(2,5)=-HALF*RMS-HALF*P(2,9)
      P(2,4)= QUARTER*RM-HALF*P(2,7)-HALF*P(2,8)-QUARTER*P(2,9)
      P(2,3)= QUARTER*RP-HALF*P(2,7)-HALF*P(2,6)-QUARTER*P(2,9)
      P(2,2)=-QUARTER*RP-HALF*P(2,5)-HALF*P(2,6)-QUARTER*P(2,9)
      P(2,1)=-QUARTER*RM-HALF*P(2,8)-HALF*P(2,5)-QUARTER*P(2,9)
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER8:GENERATES THE SHAPE FUNCTION DERIVATIVES       C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s                      C
C                                             i  j                     C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SFDER8(IELT,NDOFEL,NNE,R,S,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0,TWO=2.D0, HALF=0.5D0,QUARTER=0.25D0,FOUR=4.D0)
C
      DIMENSION P(2,NNE)
C
      RP= ONE+R
      SP= ONE+S
      RM= ONE-R
      SM= ONE-S
      RMS=ONE-R**TWO
      SMS=ONE-S**TWO
C
C     8-NODED ELEMENT
C
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,6)= HALF*SMS
      P(1,5)=-R*SM
      P(1,8)=-HALF*SMS
      P(1,7)=-R*SP
      P(1,2)= QUARTER*SM-HALF*P(1,5)-HALF*P(1,6)
      P(1,1)=-QUARTER*SM-HALF*P(1,8)-HALF*P(1,5)
      P(1,4)=-QUARTER*SP-HALF*P(1,7)-HALF*P(1,8)
      P(1,3)= QUARTER*SP-HALF*P(1,7)-HALF*P(1,6)
C
C     w.r.t.s
C
      P(2,6)=-S*RP
      P(2,5)=-HALF*RMS
      P(2,8)=-S*RM
      P(2,7)= HALF*RMS
      P(2,2)=-QUARTER*RP-HALF*P(2,5)-HALF*P(2,6)
      P(2,1)=-QUARTER*RM-HALF*P(2,8)-HALF*P(2,5)
      P(2,4)= QUARTER*RM-HALF*P(2,7)-HALF*P(2,8)
      P(2,3)= QUARTER*RP-HALF*P(2,7)-HALF*P(2,6)
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER4:GENERATES THE SHAPE FUNCTION DERIVATIVES       C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s                      C
C                                             i  j                     C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SFDER4(IELT,NDOFEL,NNE,R,S,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0,TWO=2.D0,HALF=0.5D0,QUARTER=0.25D0,FOUR=4.D0)
C
      DIMENSION P(2,NNE)
C
      RP= ONE+R
      SP= ONE+S
      RM= ONE-R
      SM= ONE-S
      RMS=ONE-R**TWO
      SMS=ONE-S**TWO
C
C     4-NODED ELEMENT
C
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,2)= QUARTER*SM
      P(1,1)=-QUARTER*SM
      P(1,4)=-QUARTER*SP
      P(1,3)= QUARTER*SP
C
C     w.r.t.s
C
      P(2,2)=-QUARTER*RP
      P(2,1)=-QUARTER*RM
      P(2,4)= QUARTER*RM
      P(2,3)= QUARTER*RP
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C     SUBROUTINE SFDER6:GENERATES THE SHAPE FUNCTION DERIVATIVES       C
C     ACCORDING TO ELEMENT TYPE AT THE POINT r ,s                      C
C                                             i  j                     C
C     B    =STRAIN-DISPLACEMENT MATRIX                                 C
C     DDET =JACOBIAN DETERMINANT                                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SFDER6(IELT,NDOFEL,NNE,R,S,P)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER(ONE=1.D0,TWO=2.D0,THREE=3.0D0, HALF=0.5D0,
     1          QUARTER=0.25D0,FOUR=4.D0,EIGHT=8.0D0, ZERO=0.0D0)
C
      DIMENSION P(2,NNE)
C
C     6-NODED ELEMENT
C
C     Derivatives w.r.t the natural coordinates
C     w.r.t.r
C
      P(1,6)=-FOUR*S
      P(1,5)= FOUR*S
      P(1,4)= FOUR-EIGHT*R-FOUR*S
      P(1,3)= ZERO
      P(1,2)= FOUR*R-ONE
      P(1,1)=-THREE+FOUR*R+FOUR*S
C
C     w.r.t.s
C
      P(2,6)= FOUR-FOUR*R-EIGHT*S
      P(2,5)= FOUR*R
      P(2,4)=-FOUR*R
      P(2,3)= FOUR*S-ONE
      P(2,2)= ZERO
      P(2,1)=-THREE+FOUR*R+FOUR*S
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF9(SN,R,S)                                      C
C      9-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF9(SN,R,S)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,HALF=0.5D0,QUART=0.25D0)
C
      DIMENSION SN(9)
C
      CALL CLEARV(SN,9)
C
      RP =ONE+R
      RM =ONE-R
      RMS=ONE-R*R
      SP =ONE+S
      SM =ONE-S
      SMS=ONE-S*S
C
      SN(9)=RMS*SMS
      SN(8)=HALF*SMS*RM-HALF*SN(9)
      SN(7)=HALF*RMS*SP-HALF*SN(9)
      SN(6)=HALF*SMS*RP-HALF*SN(9)
      SN(5)=HALF*RMS*SM-HALF*SN(9)
      SN(1)=QUART*RM*SM-HALF*SN(8)-HALF*SN(5)-QUART*SN(9)
      SN(2)=QUART*RP*SM-HALF*SN(6)-HALF*SN(5)-QUART*SN(9)
      SN(3)=QUART*RP*SP-HALF*SN(6)-HALF*SN(7)-QUART*SN(9)
      SN(4)=QUART*RM*SP-HALF*SN(8)-HALF*SN(7)-QUART*SN(9)
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF8(SN,R,S)                                      C
C      8-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF8(SN,R,S)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,HALF=0.5D0,QUART=0.25D0)
C
      DIMENSION SN(8)
C
      CALL CLEARV(SN,8)
C
      RP =ONE+R
      RM =ONE-R
      RMS=ONE-R*R
      SP =ONE+S
      SM =ONE-S
      SMS=ONE-S*S
C
      SN(8)=HALF*SMS*RM
      SN(7)=HALF*RMS*SP
      SN(6)=HALF*SMS*RP
      SN(5)=HALF*RMS*SM
      SN(1)=QUART*RM*SM-HALF*SN(8)-HALF*SN(5)
      SN(2)=QUART*RP*SM-HALF*SN(6)-HALF*SN(5)
      SN(3)=QUART*RP*SP-HALF*SN(6)-HALF*SN(7)
      SN(4)=QUART*RM*SP-HALF*SN(8)-HALF*SN(7)
C
      RETURN
C
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF4(SN,R,S)                                      C
C      4-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF4(SN,R,S)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,HALF=0.5D0,QUART=0.25D0)
C
      DIMENSION SN(4)
C
      CALL CLEARV(SN,4)
C
      RP =ONE+R
      RM =ONE-R
      RMS=ONE-R*R
      SP =ONE+S
      SM =ONE-S
      SMS=ONE-S*S
C
      SN(1)=QUART*RM*SM
      SN(2)=QUART*RP*SM
      SN(3)=QUART*RP*SP
      SN(4)=QUART*RM*SP
C
      RETURN
C
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF4DASH(SN,R,S)                                  C
C      4-noded element shape function at point R,S                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SHAPEF4DASH(SMN,R,S, KINC)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (ONE=1.0D0,HALF=0.5D0,QUART=0.25D0)
C
      DIMENSION SN(4), SMN(12,3)

      RP =ONE+R
      RM =ONE-R
      SP =ONE+S
      SM =ONE-S

      SN(1)=QUART*RP*SP
      SN(2)=QUART*RM*SP
      SN(3)=QUART*RM*SM
      SN(4)=QUART*RP*SM

	  SMN=0.0D0

	  SMN( 1,1)= SN(1)
	  SMN( 4,1)= SN(2)
	  SMN( 7,1)= SN(3)
	  SMN(10,1)= SN(4)

	  SMN( 2,2)= SN(1)
	  SMN( 5,2)= SN(2)
	  SMN( 8,2)= SN(3)
	  SMN(11,2)= SN(4)

	  SMN( 3,3)= SN(1)
	  SMN( 6,3)= SN(2)
	  SMN( 9,3)= SN(3)
	  SMN(12,3)= SN(4)

      RETURN
C
      END

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE KJACOPER                                                C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE JACOPER(NNE,XJA,XCORD,RDER)

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER(ZERO=0.0D0)

      DIMENSION XJA(2,2),XCORD(2,NNE),RDER(2,NNE)

      DUM=ZERO
      DO K=1,2
        DO J=1,2
          DO I=1,NNE
            DUM=DUM+RDER(K,I)*XCORD(J,I)
          END DO
          XJA(K,J)=DUM
          DUM=ZERO
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      2X2 Jacobian inverse                                            C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE JACINVE(XJA,DD,XJAI)

      IMPLICIT REAL*8 (A-H,O-Z)

      DIMENSION XJA(2,2),XJAI(2,2),XADJ(2,2)

      XADJ(1,1)=XJA(2,2)
      XADJ(1,2)=XJA(1,2)
      XADJ(2,1)=XJA(2,1)
      XADJ(2,2)=XJA(1,1)

      DO J=1,2
        DO K=1,2
          COFA=((-1)**(J+K))*XADJ(J,K)
          XJAI(J,K)=COFA/DD
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE AMASS9(MASS,NDOFEL,RII,SII)                          C
C      Computes the mass matrix for an 8-noded element                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS9(AMASS,NDOFEL,RII,SII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=2)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(9),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL SHAPEF9(SN,RII,SII)

      DO I=1,9
        KK=2*I-1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(2,  KK)=0.0D0
        SF(1,KK+1)=0.0D0
      END DO

      SFT=TRANSPOSE(SF)
      AUX1=MATMUL(SFT,SF)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=AUX1(I,J)
        END DO
      END DO

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE CMASS8(MASS,NDOFEL,RII,SII)                          C
C      Computes the mass matrix for an 8-noded element                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS8(AMASS,NDOFEL,RII,SII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=2)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(8),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL SHAPEF8(SN,RII,SII)

      DO I=1,8
        KK=2*I-1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(2,  KK)=0.0D0
        SF(1,KK+1)=0.0D0
      END DO

      SFT=TRANSPOSE(SF)
      AUX1=MATMUL(SFT,SF)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=AUX1(I,J)
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE CMASS6(MASS,NDOFEL,RII,SII)                          C
C      Computes the mass matrix for an 6-noded triangular element      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS6(AMASS,NDOFEL,RII,SII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=2)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(6),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL SHAPEF6(SN,RII,SII)

      DO I=1,6
        KK=2*I-1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(2,  KK)=0.0D0
        SF(1,KK+1)=0.0D0
      END DO

      SFT=TRANSPOSE(SF)
      AUX1=MATMUL(SFT,SF)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=DCMPLX(AUX1(I,J))
        END DO
      END DO

      RETURN

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE CMASSA(MASS,NDOFEL,RII,SII)                          C
C      Computes the mass matrix for an 4-noded element                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE AMASS4(AMASS,NDOFEL,RII,SII)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NTENS=2)

      DIMENSION SF(NTENS,NDOFEL),SFT(NDOFEL,NTENS),SN(4),
     1          AUX1(NDOFEL,NDOFEL),AMASS(NDOFEL,NDOFEL)

      CALL SHAPEF4(SN,RII,SII)

      DO I=1,4
        KK=2*I-1
        SF(1,  KK)=SN(I)
        SF(2,KK+1)=SN(I)
        SF(2,  KK)=0.0D0
        SF(1,KK+1)=0.0D0
      END DO

      SFT=TRANSPOSE(SF)
      AUX1=MATMUL(SFT,SF)

      DO I=1,NDOFEL
        DO J=1,NDOFEL
          AMASS(I,J)=AUX1(I,J)
        END DO
      END DO

      RETURN

      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012
C                                                                      C
C      SUBROUTINE SHAPEF6(SN,R,S)                                      C
C      6-noded triangular element shape function at point R,S          C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE SHAPEF6(SN,R,S)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ONE=1.0D0,HALF=0.5D0,FOUR=4.0D0, TWO=2.0D0)

      DIMENSION SN(6)

C     CALL CLEARV(SN,6)

      OMRMS=ONE-R-S

      SN(6)=FOUR*S*OMRMS
      SN(5)=FOUR*R*S
      SN(4)=FOUR*R*OMRMS
      SN(3)=TWO*S*(S-HALF)
      SN(2)=TWO*R*(R-HALF)
      SN(1)=TWO*(HALF-R-S)*OMRMS

      RETURN

      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE SHAPE1D                                                 C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE SHAPE1D(SM,ETA)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ONE=1.0D0,HALF=0.5D0)

      DIMENSION SM(6,2)

      EP=ONE+ETA
      EM=ONE-ETA
      EMS=ONE-ETA**2

      S11=HALF*EM-HALF*EMS
      S51=HALF*EP-HALF*EMS
      S31=EMS

      SM(1,1)=S11
      SM(2,2)=SM(1,1)
      SM(3,1)=S31
      SM(4,2)=SM(3,1)
      SM(5,1)=S51
      SM(6,2)=SM(5,1)

      SM(1,2)=0.0D0
      SM(2,1)=0.0D0
      SM(3,2)=0.0D0
      SM(4,1)=0.0D0
      SM(5,2)=0.0D0
      SM(6,1)=0.0D0

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE SHAPE1DM                                                C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE SHAPE1DM(SM,ETA)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ONE=1.0D0,HALF=0.5D0)

      DIMENSION SM(6,2)

      EP=ONE+ETA
      EM=ONE-ETA
      EMS=ONE-ETA**2

      S1=HALF*EM-HALF*EMS
      S2=HALF*EP-HALF*EMS
      S3=EMS

      SM(1,1)=S1
      SM(2,2)=SM(1,1)
      SM(3,1)=S2
      SM(4,2)=SM(3,1)
      SM(5,1)=S3
      SM(6,2)=SM(5,1)

      SM(1,2)=0.0D0
      SM(2,1)=0.0D0
      SM(3,2)=0.0D0
      SM(4,1)=0.0D0
      SM(5,2)=0.0D0
      SM(6,1)=0.0D0

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE SHAPE1D                                                 C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE SHAPEF31D(SM,ETA)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (ONE=1.0D0,HALF=0.5D0)

      DIMENSION SM(3)

      EP=ONE+ETA
      EM=ONE-ETA
      EMS=ONE-ETA**2

      SM(1)=HALF*EM-HALF*EMS
      SM(2)=HALF*EP-HALF*EMS
      SM(3)=EMS

      RETURN

      END

C23456789012345678901234567890123456789012345678901234567890123456789012
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   SUBROUTINE SHAPEMAT8                                               C
C                                                                      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890123456789012

      SUBROUTINE SHAPEMAT8(SNN,SN)

      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION SNN(2,16),SN(8)

      DO II=1,8
        SNN(1,2*II-1)=SN(II)
        SNN(1,2*II  )=0.0
        SNN(2,2*II-1)=0.0
        SNN(2,2*II  )=SN(II)
      END DO

      RETURN

      END

cc      END MODULE FEMUTIL

ext: for


text:
! Test that -print-target-triple prints correct triple

!-----------------------------------------
! RUN LINES
!-----------------------------------------
! RUN: %flang -print-target-triple 2>&1 --target=aarch64-linux-gnu \
! RUN:   | FileCheck --check-prefix=AARCH64 %s

! RUN: %flang -print-target-triple 2>&1 --target=x86_64-linux-gnu \
! RUN:   | FileCheck --check-prefix=X86_64 %s

!-----------------
! EXPECTED OUTPUT
!-----------------
! X86_64: x86_64-unknown-linux-gnu
! AARCH64: aarch64-unknown-linux-gnu

max_stars_repo_name: lxbndr/llvm-project

max_issues_repo_name: lxbndr/llvm-project

max_forks_repo_name: lxbndr/llvm-project


text:
C    @(#)kmpuvov.f	20.3 2/13/96
	integer function kmpuvov (p, q)
        integer p, q
c
c	This function compares vltsrt(p) with vltsrt(q) using fields
c       <area_name><bus_name><base_kv>
c
        include 'ipfinc/parametr.inc'
        include 'ipfinc/blank.inc'
        include 'ipfinc/bus.inc'
        include 'ipfinc/area.inc'
        include 'ipfinc/arcntl.inc'
        include 'ipfinc/sortuvov.inc'

        if (p .eq. q) then
           kmpuvov = 0
        else
           m = vltsrt(p)
           n = vltsrt(q)
           if (inp2alf(m) .gt. ntot_alf .or. inp2alf(n) .gt. ntot_alf)
     &        then
              kmpuvov = kompr (bus(m), bus(n), junk)
              if (kmpuvov .eq. 0) then
                 kmpuvov = 100.0 * (base(m) - base(n))
              endif
              if (kmpuvov .eq. 0) then
                 kmpuvov = kompr (arcnam(jarzn(m)), arcnam(jarzn(n)),
     &                            junk)
              endif
           else
              kmpuvov = kompr (arcnam(jarzn(m)), arcnam(jarzn(n)), junk)
              if (kmpuvov .eq. 0) then
                 kmpuvov = kompr (bus(m), bus(n), junk)
              endif
              if (kmpuvov .eq. 0) then
                 kmpuvov = 100.0 * (base(m) - base(n))
              endif
           endif
        endif
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C 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 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 NTESS 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 $Log: getvar.f,v $
C Revision 1.2  2009/03/25 12:36:44  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.1  1994/04/07 20:01:52  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.2  1990/12/14  08:51:15  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE GETVAR (A, IVAR, IELBLK, INSTEP, LENVAR, VAR)
C=======================================================================

C   --*** GETVAR *** (BLOT) Read variable
C   --   Written by Amy Gilkey - revised 05/17/88
C   --
C   --GETVAR returns the values for the requested variable for the
C   --requested time step.  It either reads the values from the sequential
C   --database file and writes them to a direct access scratch file or it
C   --reads the values from the direct access scratch file.
C   --
C   --This routine uses MDFIND to find the following dynamic memory arrays:
C   --   NUMELB - the number of elements per element block
C   --   ISEVOK - the element block variable truth table;
C   --      variable i of block j exists iff ISEVOK(j,i)
C   --   WHOTIM - true iff whole (versus history) time step
C   --
C   --Parameters:
C   --   A - IN - the dynamic memory base array
C   --   IVAR - IN - the variable index
C   --      (0 to initialize random file only)
C   --   IELBLK - IN - the element block number, <=0 for all
C   --      (for element blocks only)
C   --   INSTEP - IN - the time step number
C   --      = +n to read time step n
C   --      = -n to transfer time step n to random file only
C   --      =  0 to transfer all time steps to random file
C   --   LENVAR - IN - the length of VAR
C   --   VAR - OUT - the variable values (indeterminate if INSTEP <= 0)
C   --
C   --Common Variables:
C   --   Uses NDB of /DBASE/
C   --   Uses NUMNP, NUMEL, NELBLK, NVARHI, NVARGL, NVARNP, NVAREL,
C   --      NSTEPS, NSTEPW of /DBNUMS/
C   --
C   --Database is rewound upon the first entry of this routine; upon
C   --exit a flag is set to keep track of the database position; the
C   --database should not be moved between calls to this routine.
C   --
C   --A scratch random file is created and read and written in this routine.
C   --It is connected to unit 90.

      DIMENSION A(*)
      REAL VAR(*)

      IF (IVAR .LE. 0) RETURN
      CALL RNDVAR (A, A, A, IVAR, IELBLK, INSTEP, LENVAR, VAR)

      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
C    @(#)file_sts.f	20.3 11/11/97
C****************************************************************
C
C   File: file_sts.f
C   Purpose: Routine to process /GET_DATA, TYPE=FILE_WRT commands
C            Returned value:
C
C            status = 0 : file is writable or does not exsit
C                     1 : file exists and is not writable
C
C   Author: Walt Powell  Date: 10 Nov 1992
C   Called by: p_gtdata.f
C   Re-written by Jay Coleman  30-SEP-93
C
C****************************************************************
        integer function file_sts (in_buffer, out_buffer)
        character in_buffer *(*)
        character out_buffer *(*)

        include 'ipfinc/prt.inc'
        include 'ipfinc/miscfile.inc'

        integer index, lastch

        character word(20)*80, capital*80, null*1, linefeed*1,
     &            filename*80, formatted*10, form*1
        logical yes, found
        integer o2, iwrd, nwrd, last, status

        null = char(0)
        linefeed = char(10)
        file_sts = 0            ! Initialize to "does not exist"
        out_buffer(1:1) = null
        yes = .false.
        last = index (in_buffer, null)
        if (last .eq. 0) last = 120

        call uscan (in_buffer(1:last), word, nwrd, '=' ,
     &              ', ' // linefeed // null )

c       Align "iwrd" after "FILE ="

        iwrd = 1
        do while (iwrd .le. nwrd)
           word(iwrd) = capital(word(iwrd))
           if (word(iwrd)(1:5) .eq. 'FILE '  .and.
     &         word(iwrd+1) .eq. '='               ) then
              iwrd = iwrd + 2
              filename = word(iwrd)
              last = lastch (filename)
              inquire (file = filename(1:last), exist = yes,
     &                 formatted = formatted)
              if (formatted .eq. 'FORMATTED') then
                 form = 'F'
              else
                 form = 'U'
              endif
              if (yes) then
                 call close_file(lunscr1)
                 status = open_file (lunscr1, filename(1:last), form,
     &                               'U', iostat)
                 if (iostat .eq. 0) then
                    file_sts = 0   ! File exists and is writable
                 else
                    file_sts = 1   ! File exists and is not writable
                 endif
                 call close_file(lunscr1)
                 return
              endif
              return
           endif
           iwrd = iwrd + 1
        enddo

        write (errbuf(1), 100)
  100   format (' Command parsing error, missing "FILE =" keyword ')
        call prterx ('W', 1)
        return
        end


    MODULE cdf_binomial_mod
! ______________________________________________________________________
!
!                             cdf_binomial_mod
!                             *=*=*=*=*=*=*=*=
!
!  -  SUBROUTINE CDF_BINOMIAL(WHICH, CUM, CCUM, S, N, PR, CPR, STATUS,
!                             CHECK_INPUT)
!  -  REAL (dpkind) FUNCTION CUM_BINOMIAL(S, N, PR, CPR, STATUS, CHECK_INPUT)
!  -  REAL (dpkind)FUNCTION CCUM_BINOMIAL(S, N, PR, CPR, STATUS, CHECK_INPUT)
!  -  REAL (dpkind) FUNCTION INV_BINOMIAL(CUM, CCUM, N, PR, CPR, STATUS,
!                                         CHECK_INPUT)
!
!                             The Distribution
!                             ================
!
!   The density of the binomial distribution provides the probability of S
!   successes in N independent trials, each with probability of success
!   PR.  The density is proportional to:
!
!                                S       N-S
!                              PR  (1-PR)
!
!     The binomial is extended to non-integer values via the connection
!     between the cumulative binomial and the incomplete beta.
!
!                                Arguments
!                                =========
!
!  - INTEGER, INTENT(IN) :: WHICH. Integer indicating which of the next
!    four arguments is to be calculated.
!  Input Range: [ 1:4 ]
!     - 1 CUM and CCUM
!     - 2 S
!     - 3 N
!     - 4 PR and CPR
!  - REAL (dpkind), OPTIONAL :: CUM. The CDF of the binomial distribution,
!    i.e., the probability of 0 to S successes in N trials.
!  Range: [ 0:1-10^-10 ]
!  - REAL (dpkind), OPTIONAL :: CCUM. One minus the CDF of the binomial
!    distribution.
!  Range: [ 10^-10:1 ]
!  - REAL (dpkind) :: S. The upper limit of summation of the binomial
!    density. Note that S must be less than or equal to N.
!  Range: [ 0:10^10 ]
!  - REAL (dpkind) :: N. The number of independent trials generating the
!    binomial density. N must be greater than or equal to S.
!  Range: [ 0:10^10 ]
!  - REAL (dpkind), OPTIONAL :: PR. The probability of success in each
!    independent trial.
!  Range: [ 0:1 ]
!  - REAL (dpkind), OPTIONAL :: CPR. One minus the probability of success
!    in each independent trial; the probability of failure in each trial.
!  Range: [ 0:1 ]
!  - INTEGER, OPTIONAL, INTENT(OUT) :: STATUS. Return code.
!        0 problem solved successfully
!       -1 WHICH outside input range
!       -2 CUM outside range
!       -3 CCUM outside range
!       -4 S outside range
!       -5 N outside range
!       -6 PR outside range
!       -7 CPR outside range
!        3 CUM + CCUM is not nearly one
!        4 PR + CPR is not nearly one
!        5 S not between 0 and N.
!       10 cdf_beta (in local_cum_binomial) has no answer.
!      -50 Answer (if any) is BELOW the LOWER search bound
!       50 Answer (if any) is ABOVE the UPPER search bound
!  - LOGICAL, INTENT(IN), OPTIONAL :: CHECK_INPUT. If PRESENT and
!       .TRUE. input argument values are not checked for validity.
!
!  NOTE: CUM and CCUM and also PR and CPR must add to (nearly) one.
! ______________________________________________________________________
! .. Use Statements ..
      USE biomath_constants_mod
! ..
! .. Implicit None Statement ..
      IMPLICIT NONE
! ..
! .. Default Accessibility ..
      PRIVATE
! ..
! .. Public Statements ..
      PUBLIC :: ccum_binomial, cdf_binomial, cum_binomial, inv_binomial
! ..
    CONTAINS

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

      FUNCTION ccum_binomial(s,n,pr,cpr,status,check_input)
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Scalar Arguments ..
        REAL (dpkind), OPTIONAL :: cpr, pr
        REAL (dpkind) :: n, s
        INTEGER, OPTIONAL, INTENT (OUT) :: status
        LOGICAL, OPTIONAL, INTENT (IN) :: check_input
! ..
! .. Function Return Value ..
        REAL (dpkind) :: ccum_binomial
! ..
! .. Executable Statements ..

        CALL cdf_binomial(which=1,ccum=ccum_binomial,s=s,n=n,pr=pr, &
          cpr=cpr,status=status,check_input=check_input)

        RETURN

      END FUNCTION ccum_binomial

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

      SUBROUTINE cdf_binomial(which,cum,ccum,s,n,pr,cpr,status, &
          check_input)
! .. Use Statements ..
        USE cdf_aux_mod
        USE zero_finder
! ..
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Scalar Arguments ..
        REAL (dpkind), OPTIONAL :: ccum, cpr, cum, pr
        REAL (dpkind) :: n, s
        INTEGER, OPTIONAL, INTENT (OUT) :: status
        INTEGER, INTENT (IN) :: which
        LOGICAL, OPTIONAL, INTENT (IN) :: check_input
! ..
! .. Local Structures ..
        TYPE (zf_locals) :: local
! ..
! .. Local Arrays ..
        REAL (dpkind) :: params(6)
! ..
! .. Local Scalars ..
        REAL (dpkind) :: fx, local_ccum, local_cpr, local_cum, local_pr, &
          try_ccum, try_cpr, try_cum, try_pr
        INTEGER :: zf_status
        LOGICAL :: has_status, local_check_input, match_cum, vary_pr
! ..
! .. Intrinsic Functions ..
        INTRINSIC MAX, PRESENT
! ..
! .. Executable Statements ..

        has_status = PRESENT(status)

        ! status = 0 means NO error

        IF (has_status) THEN
           status = 0
        END IF

        ! Check that at least one of cum, ccum is present in the calling list

        CALL check_complements(cum,ccum,the_binomial%name,'cum', 'ccum', &
             local_cum,local_ccum,set_values=(which/=1), bad_status=3, &
             status=status)

        IF (has_status) THEN
           IF (status /= 0) THEN
              RETURN
           END IF
        END IF

        ! Check that at least one of pr, cpr is present in the calling list

        CALL check_complements(pr,cpr,the_binomial%name,'pr','cpr', &
             local_pr,local_cpr,set_values=(which/=4), bad_status=4, &
             status=status)

        IF (has_status) THEN
           IF (status /= 0) THEN
              RETURN
           END IF
        END IF

        !==========
        ! Check that probabilities add to one
        !==========

        IF (which /= 4) THEN
           IF ( .NOT. add_to_one(local_pr,local_cpr,the_binomial%name, &
                'pr','cpr',4,status)) RETURN
        END IF

        params(1) = local_cum
        params(2) = local_ccum
        params(3) = s
        params(4) = n
        params(5) = local_pr
        params(6) = local_cpr

        IF (PRESENT(check_input)) THEN
          local_check_input = check_input
        ELSE
          local_check_input = .TRUE.
        END IF

        IF (local_check_input) THEN

          CALL validate_parameters(the_binomial,which,params,status)

          IF (has_status) THEN
             IF (status /= 0) THEN
                RETURN
             END IF
          END IF

          IF (which/=2 .AND. which/=3) THEN

            ! Assure that s is between 0 and n

            IF ( .NOT. in_range(value=s, lo=zero, hi=n, routine_name= &
              the_binomial%name, arg_name='s', bad_status=5,status=status)) &

              RETURN
          END IF

        END IF

        !==========
        ! Calculate answers
        !==========

        IF (which>1) match_cum = (local_cum <= half)

        SELECT CASE (which)

        CASE (1)
           ! Solve for the cumulative binomial (cdf)

           CALL local_cum_binomial(s,n,local_pr,local_cpr,local_cum, &
                local_ccum, status)

           IF (has_status) THEN
              IF (status /= 0) THEN
                 ! Problem can NOT be solved because cdf_beta has NO solution.
                 RETURN
              END IF
           END IF

           IF (PRESENT(cum)) cum = local_cum
           IF (PRESENT(ccum)) ccum = local_ccum

           RETURN

        CASE (2)
           ! Solve for s (number of successes)

          s = half * n
          zf_status = 0

! DMS
! Clearly s <= n, so the UPPER bound for s should be n
! We can not call cdf_set_zero() as in the other cases, because
! the_binomial structures contains PARAMETERS which can NOT be modified.
! The values are set like in  cdf_set_zero()

          CALL set_zero_finder( &
               low_limit=the_binomial%parameters(3)%low_bound, &
               hi_limit=n, abs_step=half, rel_step=half, &
               step_multiplier=five, abs_tol=1.0E-50_dpkind, &
               rel_tol=1.0E-8_dpkind, local=local)

          DO
!            CALL rc_interval_zf(zf_status, s, fx, local)
             CALL rc_step_zf(zf_status, s, fx, local)

            IF (zf_status /= 1) EXIT

            CALL local_cum_binomial(s, n, local_pr, local_cpr, try_cum, &
                                    try_ccum, status)

            IF (has_status) THEN
               IF (status /= 0) THEN
                  RETURN
               END IF
            END IF

            IF (match_cum) THEN
              fx = try_cum - local_cum
            ELSE
              fx = local_ccum - try_ccum
            END IF

          END DO

        CASE (3)
          ! Solve for n (sample size)

          n = two * MAX(s,one)
          zf_status = 0

          CALL cdf_set_zero_finder(the_binomial,4,local)

          DO
            CALL rc_step_zf(zf_status,n,fx,local)

            IF (zf_status/=1) EXIT

            CALL local_cum_binomial(s,n,local_pr,local_cpr,try_cum, &
                                    try_ccum, status)

            IF (has_status) THEN
               IF (status /= 0) THEN
                  RETURN
               END IF
            END IF

            IF (match_cum) THEN
              fx = try_cum - local_cum
            ELSE
              fx = try_ccum - local_ccum
            END IF
          END DO

        CASE (4)
           ! Solve for the binomial probability (pr) or
           ! its complement (cpr)

          zf_status = 0

!!! A difficult  case -- we  will try to  match the lesser of  cum and
!!! ccum If the answer <= 1/2 we will vary pr else we will vary cpr

! Decide whether to vary pr or cpr

          CALL local_cum_binomial(s,n,half,half,try_cum,try_ccum, status)

          IF (has_status) THEN
             IF (status /= 0) THEN
                RETURN
             END IF
          END IF

          vary_pr = (try_cum<=local_cum)

          IF (match_cum .AND. vary_pr) THEN

             ! We are trying to match cum varying pr

            CALL cdf_set_zero_finder(the_dummy_binomial,1,local)

            DO
              CALL rc_interval_zf(zf_status,try_pr,fx,local)

              try_cpr = one - try_pr

              IF (zf_status/=1) EXIT

              CALL local_cum_binomial(s,n,try_pr,try_cpr,try_cum, &
                                      try_ccum, status)

              IF (has_status) THEN
                 IF (status /= 0) THEN
                    RETURN
                 END IF
              END IF

              fx = local_cum - try_cum
            END DO

          ELSE IF (match_cum .AND. .NOT. vary_pr) THEN
            ! Try to match cum varying cpr

            CALL cdf_set_zero_finder(the_dummy_binomial,2,local)

            DO

              CALL rc_interval_zf(zf_status,try_cpr,fx,local)

              try_pr = one - try_cpr

              IF (zf_status/=1) EXIT

              CALL local_cum_binomial(s,n,try_pr,try_cpr,try_cum, &
                                      try_ccum, status)

              IF (has_status) THEN
                 IF (status /= 0) THEN
                    RETURN
                 END IF
              END IF

              fx = local_cum - try_cum
            END DO

          ELSE IF ( .NOT. match_cum .AND. vary_pr) THEN
! Try to match ccum varying pr

            CALL cdf_set_zero_finder(the_dummy_binomial,1,local)

            DO

              CALL rc_interval_zf(zf_status,try_pr,fx,local)

              try_cpr = one - try_pr

              IF (zf_status/=1) EXIT

              CALL local_cum_binomial(s,n,try_pr,try_cpr,try_cum, &
                                      try_ccum, status)

              IF (has_status) THEN
                 IF (status /= 0) THEN
                    RETURN
                 END IF
              END IF

              fx = local_ccum - try_ccum
            END DO

          ELSE

            ! Try to match ccum varying cpr

            CALL cdf_set_zero_finder(the_dummy_binomial,2,local)

            DO

              CALL rc_interval_zf(zf_status,try_cpr,fx,local)

              try_pr = one - try_cpr

              IF (zf_status/=1) EXIT

              CALL local_cum_binomial(s,n,try_pr,try_cpr,try_cum, &
                                      try_ccum, status)

              IF (has_status) THEN
                 IF (status /= 0) THEN
                    RETURN
                 END IF
              END IF

              fx = local_ccum - try_ccum
            END DO

          END IF

          IF (PRESENT(pr)) pr = try_pr

          IF (PRESENT(cpr)) cpr = try_cpr

        END SELECT

        IF (has_status) THEN
           ! Set the status of the zero finder
           CALL cdf_finalize_status( local, status )
        END IF

        RETURN

      END SUBROUTINE cdf_binomial

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

      FUNCTION cum_binomial(s,n,pr,cpr,status,check_input)
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Scalar Arguments ..
        REAL (dpkind), OPTIONAL :: cpr, pr
        REAL (dpkind) :: n, s
        INTEGER, OPTIONAL, INTENT (OUT) :: status
        LOGICAL, OPTIONAL, INTENT (IN) :: check_input
! ..
! .. Function Return Value ..
        REAL (dpkind) :: cum_binomial
! ..
! .. Executable Statements ..

        CALL cdf_binomial(which=1,cum=cum_binomial,s=s,n=n,pr=pr,cpr=cpr, &
                          status=status,check_input=check_input)

        RETURN

      END FUNCTION cum_binomial

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

      FUNCTION inv_binomial(cum,ccum,n,pr,cpr,status,check_input)
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Scalar Arguments ..
        REAL (dpkind), OPTIONAL :: ccum, cpr, cum, pr
        REAL (dpkind) :: n
        INTEGER, OPTIONAL, INTENT (OUT) :: status
        LOGICAL, OPTIONAL, INTENT (IN) :: check_input
! ..
! .. Function Return Value ..
        REAL (dpkind) :: inv_binomial
! ..
! .. Executable Statements ..

        CALL cdf_binomial(2,cum,ccum,inv_binomial,n,pr,cpr,status, &
                          check_input)

        RETURN

      END FUNCTION inv_binomial

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

      SUBROUTINE local_cum_binomial(s,xn,pr,ompr,cum,ccum, status)
!----------------------------------------------------------------------
!                    CUMulative BINomial distribution
!
!                              Function
!
!     Returns the probability   of 0  to  S  successes in  XN   binomial
!     trials, each of which has a probability of success, PBIN.
!                              Arguments
!     S --> The upper limit of cumulation of the binomial distribution.
!
!     XN --> The number of binomial trials.
!
!     PR --> The probability of success in each binomial trial.
!
!     OMPR --> 1 - PBIN
!
!     CUM <-- Cumulative binomial distribution.
!
!     CCUM <-- Complement of Cumulative binomial distribution.
!
!     STATUS <-- status of the computation.
!                Set to 10 if cdf_beta has no answer.
!
!                              Method
!
!     Formula  26.5.24    of   Abramowitz  and    Stegun,  Handbook   of
!     Mathematical   Functions (1966) is   used  to reduce the  binomial
!     distribution  to  the  cumulative    beta distribution.
!----------------------------------------------------------------------
! .. Use Statements ..
        USE cdf_beta_mod
! ..
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Scalar Arguments ..
        REAL (dpkind), INTENT (OUT) :: ccum, cum
        REAL (dpkind), INTENT (IN) :: ompr, pr, s, xn
        INTEGER, OPTIONAL, INTENT(OUT) :: status
! ..
! .. Local Scalars
        INTEGER :: beta_status
! ..
! .. Intrinsic Functions ..
        INTRINSIC PRESENT
! ..
! .. Executable Statements ..

        IF (s < xn) THEN
          CALL cdf_beta(1,ccum,cum,pr,ompr,s+one,xn-s, status=beta_status)

          IF (beta_status /= 0) THEN
             ! cdf_beta has NO answer.

             IF (PRESENT(status)) THEN
                status = 10
             ELSE
                WRITE (*,*) 'Error in local_cum_binomial call to cdf_beta'
                WRITE (*,*) 'Status: ', beta_status

                STOP 'Error in local_cum_binomial call to cdf_beta'
             END IF
          END IF
        ELSE
          cum = one
          ccum = zero
        END IF

        RETURN

      END SUBROUTINE local_cum_binomial

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

    END MODULE cdf_binomial_mod


text:
	SUBROUTINE G2T_RANGE ( ktype, mval, irange, iret )
C************************************************************************
C* G2T_RANGE 								*
C*									*
C* This subroutine returns the range for a give wave or wind value.	*
C*									*
C* G2T_RANGE ( KTYPE, MVAL, IRANGE, IRET )				*
C*									*
C* Input parameters:							*
C*	KTYPE		INTEGER		Data type			*
C*					 1 = Wave			*
C*					 2 = Wind			*
C*	MVAL		INTEGER		Maximum wave height/wind	*
C*									*
C* Output parameters:							*
C*	IRANGE		INTEGER		Maxium range			*
C*	IRET		INTEGER		Return code			*
C*					 0 = normal return		*
C**									*
C* Log:									*
C* T. Lee/SAIC		 9/06						*
C* T. Lee/SAIC		11/07	Use g2t_parm.tbl to get range		*
C************************************************************************
	INCLUDE		'goftxt.cmn'
	LOGICAL		contin
	INCLUDE 	'EQUAL.FNC'
C-----------------------------------------------------------------------
	iret = 0
C
	IF ( IRMISS ( mval ) )  THEN
	    irange = 0
	    RETURN
	END IF
C
C*	Wave heights.
C
	IF  ( ktype .eq. 1 )  THEN
	    IF ( ncnter .eq. 2 )  THEN
		irange = NINT ( mval / 3. )
	      ELSE
		kw = 1
		contin = .true.
		DO WHILE ( contin )
		    IF ( mval .ge. mnwgap ( kw, 1  )  .and.
     +			 mval .lt. mnwgap ( kw, 2  ) )  THEN
			contin = .false.
			irange = mnwgap ( kw, 3 )
		    END IF
		    kw = kw + 1
		    IF ( kw .gt. llwgap ) contin = .false.
		END DO
	    END IF
	  ELSE IF ( ktype .eq. 2 )  THEN
	    ks = 1
	    contin = .true.
	    DO WHILE ( contin )
		IF ( mval .gt. mnsgap ( ks, 1  )  .and.
     +		     mval .le. mnsgap ( ks, 2  ) )  THEN
		    contin = .false.
		    irange = mnsgap ( ks, 3 )
		END IF
		ks = ks + 1
		IF ( ks .gt. llsgap ) contin = .false.
	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Campus Recycling has seen a major decrease in revenue thanks to the administrations decision to stop baling cardboard.  (see http://www.californiaaggie.com/media/storage/paper981/news/2005/01/19/FrontPage/Campus.Scraps.OnSite.Cardboard.Baling1319230.shtml aggie story)

UC Davis used to collect, sort (to remove contamination), and bale all cardboard collected on campus. This was important for campus recycling because we could separate out contamination and sell the clean bales of cardboard directly to the mill. This allowed us to receive higher rates for the material.  BR
Many campus Environmental Organizations are upset about the decision to stop baling because cardboard is now shipped loose to Sacramento  a process that costs us much more fuel for transportation and a sharp decrease in payment. Occasionally the loads end up in the Landfill due to contamination or loads that do not fill the truck (these loads are then mixed with trash and sent together to the Landfill).

The baler is now for sale at the Bargain Barn.

It has been suggested that ASUCD should buy the baler and turn it into a source of revenue.

The cardboard is hauled into Sacramento by the same truck that goes to collect garbage from the Med Center bins. There isnt an increased fuel usage. Think of it as the truck being reused which is one of the four Rs.  Users/JasonAller

Well, when the cardboard is shipped loose the trucks actually have to make more runs to sac to drop the cardboard off... more trips   more fuel ~~ Users/MikeyNolan

20050221 20:59:41 nbsp What are the economics of the cardboard bailer? Was it losing money? Users/JackHaskel

20050221 21:36:47 nbsp Jason makes a good point  the big problem with recycling cardboard at UCD is medwaste contamination from the Med Center and Vet school. Blood and urine were sorted out by R4 staff (full timers) and the shipping info from these sources of contamination was collected so R4 could kindly ask those contaminators to stop (a major function of the group is eductaion). If cardboard is contaminated (by whats in the truck or dumpsters) then it is thrown out. Now, no one will see what goes to the recycler before it gets there  entire truck loads can be landfilled if contamination is significant.  Users/MikeSiminitus
20050221 21:39:20 nbsp In terms of economics  baling cardboard allowed us to store the bales until prices were best for selling the stuff. Cardboard prices fluctuate quite a bit. Save the lengthy costbenifit for the admin types. Is the landfill generating revinue? Why should recycling need to be selfsufficent when solid waste is not? Users/MikeSiminitus

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE DS3 (R1,R2,R3,R4,R5,R6, ARRAY)

C  Set an array with six real numbers

      REAL R1,R2,R3,R4,R5,R6,ARRAY(6)

      ARRAY(1) = R1
      ARRAY(2) = R2
      ARRAY(3) = R3
      ARRAY(4) = R4
      ARRAY(5) = R5
      ARRAY(6) = R6

      END