text:
C$Procedure ZZPDPLTC (Planetodetic coordinates, point latitude check)

      LOGICAL FUNCTION ZZPDPLTC ( RE, F, P, LAT )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines. Users should not call this routine directly due to the
C     volatile nature of this routine.
C
C     Indicate whether a given point on a planetodetic latitude cone
C     has the correct latitude sign.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     COORDINATES
C     GEOMETRY
C     PLANETODETIC
C     LATITUDE
C     MATH
C
C$ Declarations

      IMPLICIT NONE

      DOUBLE PRECISION      RE
      DOUBLE PRECISION      F
      DOUBLE PRECISION      P   ( 3 )
      DOUBLE PRECISION      LAT

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     RE         I   Equatorial radius.
C     F          I   Flattening coefficient.
C     P          I   Three-dimensional point.
C     LAT        I   Planetodetic latitude.
C
C     The function returns .TRUE. if the sign of the planetodetic
C     latitude of the input point matches that of the input
C     planetodetic latitude.
C
C$ Detailed_Input
C
C     RE,
C     F          are, respectively, the equatorial radius
C                and flattening coefficient of a biaxial
C                spheroid.
C
C                The polar radius RP of the spheroid is
C
C                   RP = RE * ( 1 - F )
C
C                RP may be less than, equal to, or greater than RE.
C
C
C     P          is a point (equivalently, a vector) in
C                three-dimensional space. P is expressed in Cartesian
C                coordinates. P must lie on the planetodetic
C                latitude cone corresponding to LAT.
C
C                The units of P must be consistent with those of RE.
C
C
C     LAT        is a planetodetic latitude value defining a cone.
C                Units are radius.
C
C
C$ Detailed_Output
C
C     The function returns .TRUE. if any of the following
C     are true:
C
C        - The input spheroid is prolate or spherical.
C
C        - The input latitude is zero.
C
C        - The input point is determined to have planetodetic
C          latitude of the same sign as the input latitude.
C
C    Otherwise the function returns .FALSE.
C
C    See Particulars below for details.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If either the equatorial radius or flattening coefficient is
C         invalid, the error SPICE(VALUEOUTOFRANGE) will be signaled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine supports computation of the intersection between
C     a ray and a planetodetic volume element. The routine is used
C     to determine whether the intersection of the ray and latitude
C     cone is actually on the element's boundary.
C
C     This function serves as a "macro" that executes a logical test
C     that would be awkward to perform in-line.
C
C     For a given reference spheroid, all points having a given
C     planetodetic latitude lie on a cone. The axis of the cone
C     coincides with the Z axis of the coordinate system.
C
C     The possibility that a point on a given latitude can have a Z
C     coordinate with the wrong sign exists for oblate reference
C     spheroids. For these shapes, the vertex of a positive latitude
C     cone has a negative Z component, and the vertex of a negative
C     latitude cone has a positive Z component.
C
C     The purpose of the function is to determine, for a point that
C     lies on a given planetodetic latitude cone, whether that point is
C     on the correct side of the X-Y plane: that is, the side
C     corresponding to the sign of the input latitude value.
C
C     This check is not as simple as checking the sign of the Z
C     component of the input point. For values of LAT that are non-zero
C     but have very small magnitude, points that are nominally on the
C     corresponding latitude cone may have Z components of the wrong
C     sign due to round-off errors that occurred in the process of
C     computing those points. For such cases, the input point is
C     checked by comparing its distance from the Z axis to the X
C     intercept of a line normal to the reference spheroid at a point
C     having planetodetic latitude LAT and longitude 0. This check can
C     be performed accurately even when the Z component of the input
C     point consists of noise.
C
C     This routine does not test the input point to determine whether
C     it lies on the latitude cone defined by the input argument LAT.
C     The caller must ensure that this is the case.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     This routine does not test the input point to determine whether
C     it lies on the latitude cone defined by the input argument LAT.
C     The caller must ensure that this is the case.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 13-JAN-2017 (NJB)
C
C-&

C$ Index_Entries
C
C     is point on planetodetic latitude boundary
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
      DOUBLE PRECISION      LIMIT
      PARAMETER           ( LIMIT = 1.D-2 )

C
C     Local variables
C
      DOUBLE PRECISION      A
      DOUBLE PRECISION      B
      DOUBLE PRECISION      R
      DOUBLE PRECISION      R2
      DOUBLE PRECISION      XXPT
      DOUBLE PRECISION      YXPT

C
C     Give the function a default value.
C
      ZZPDPLTC = .FALSE.

      IF ( RETURN() ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'ZZPDPLTC' )

C
C     The equatorial radius must be greater than zero.
C
      IF ( RE .LE. 0.0D0 ) THEN

          CALL SETMSG ( 'Equatorial radius was *.' )
          CALL ERRDP  ( '*', RE                    )
          CALL SIGERR ( 'SPICE(VALUEOUTOFRANGE)'   )
          CALL CHKOUT ( 'ZZPDPLTC'                 )
          RETURN

      END IF

C
C     If the flattening coefficient is greater than one, the polar
C     radius computed below is negative. If it's equal to one, the
C     polar radius is zero. Either case is a problem, so signal an
C     error and check out.
C
      IF ( F .GE. 1 ) THEN

          CALL SETMSG ( 'Flattening coefficient was *.'  )
          CALL ERRDP  ( '*', F                           )
          CALL SIGERR ( 'SPICE(VALUEOUTOFRANGE)'         )
          CALL CHKOUT ( 'ZZPDPLTC'                       )
          RETURN

      END IF

C
C     The input point is assumed to be on the cone
C     corresponding to the input latitude.
C
C     If the reference spheroid is prolate or spherical, there's
C     nothing to do: the point is automatically on the correct side of
C     the X-Y plane.
C
      IF ( F .LE. 0.D0 ) THEN

         ZZPDPLTC = .TRUE.

      ELSE
C
C        This is the oblate case.
C
C
C        If the point is on the "correct" side of the X-Y plane---
C        that is, its Z component as the same sign as LAT, the
C        point is considered to have the correct latitude.
C
C        If the point is on the X-Y plane, or if LAT is zero, the point
C        is considered to have the indicated latitude. We condense
C        these cases by requiring that
C
C              LAT * P(3) >= 0
C
C           rather than
C
C              LAT * P(3) > 0
C
C
         IF (  ( P(3) * LAT )  .GE.  0.D0  ) THEN

            ZZPDPLTC = .TRUE.

         ELSE IF ( ABS(LAT) .GE. LIMIT ) THEN
C
C           Ideally, the input point is considered to have the given
C           latitude if the point is on the side of the X-Y plane
C           corresponding to the sign of the input latitude. The
C           problem with this criterion is that it can't be applied
C           correctly when LAT has very small magnitude.
C
C           If the magnitude of LAT is above the limit, it's ok to
C           use the sign of P(3) to determine whether the point
C           has the given latitude.
C
C           The point has the indicated latitude if both LAT and P(3)
C           have the same sign. In this case, we know they have the
C           opposite sign.
C
            ZZPDPLTC = .FALSE.

         ELSE
C
C           At this point we know LAT is non-zero, so the cone
C           corresponding to LAT has its vertex on the opposite side of
C           the X-Y plane from any point having latitude LAT. So it's
C           possible for a point to be on the cone but not have the
C           correct latitude.
C
C           We're in the special case where the point's Z component
C           has the opposite sign as LAT, and the magnitude of LAT
C           is below the limit. We don't automatically reject the
C           point in this case: we'll accept it if it is far enough
C           from the Z axis to be outside the portion of the latitude
C           cone on the wrong side of the X-Y plane.
C
            A = RE
            B = A * ( 1.D0 - F )

C
C           Compute the intercepts of a normal vector of a point
C           at latitude LAT, longitude 0, with the X and Y axes.
C
            CALL ZZELNAXX ( A, B, LAT, XXPT, YXPT )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZPDPLTC' )
               RETURN
            END IF

C
C           We check the point's distance from the Z axis. This can be
C           done accurately even when the Z component of P consists of
C           noise.
C
C           The point is considered to have the correct latitude when
C           it is farther from the Z axis than the intercept on the X
C           axis of a normal line passing through a point having
C           latitude LAT and longitude 0. Ideally, a point that is on
C           the latitude cone and that satisfies this criterion must be
C           on the correct side of the X-Y plane.
C

            R2       = ( P(1) * P(1) ) +  ( P(2) * P(2) )

            R        = SQRT(  MAX( R2, 0.D0 )  )

            ZZPDPLTC = R .GE. XXPT

         END IF

      END IF

      CALL CHKOUT ( 'ZZPDPLTC' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
Filipino Martial Arts refers to many forms of fighting systems that originated in, or were refined in, the Philippines. The most wellknown, traditional branch of Filipino martial arts is Eskrima, an umbrella term primarily used for several forms of weaponfighting. Although the aspect of Eskrima most often portrayed in the media is rattanstickfighting, training with other weaponry such as knives, swords, and longer rods is also common. Some Eskrima schools also emphasize emptyhanded techniques.

Instructors that have a Davis presence include those associated with Baciarinis Martial Arts and Doce Pares Eskrima Filipino Martial Arts Academy.

Baciarinis started out as Pallens Martial Arts. The Eskrima taught there follows the style known as Senkotiros, continuing the lineage of Grandmaster Max Pallen. Pallen stripped down Eskrima to the core essentials, creating Senkotiros, which, in his belief, made the form more practical and effective.

Doce Pares is notable for teaching not only Eskrima but also YawYan Ardigma (Filipino Kickboxing). Being that Baciarinis began including a wide range of martial art practices in recent years, Doce Pares is probably the only school to be 100% based on Filipino techniques.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*> \brief \b ZSYSWAPR
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZSYSWAPR + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZSYSWAPR( UPLO, N, A, LDA, I1, I2)
*
*       .. Scalar Arguments ..
*       CHARACTER        UPLO
*       INTEGER          I1, I2, LDA, N
*       ..
*       .. Array Arguments ..
*       COMPLEX*16       A( LDA, N )
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZSYSWAPR applies an elementary permutation on the rows and the columns of
*> a symmetric matrix.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          Specifies whether the details of the factorization are stored
*>          as an upper or lower triangular matrix.
*>          = 'U':  Upper triangular, form is A = U*D*U**T;
*>          = 'L':  Lower triangular, form is A = L*D*L**T.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA,N)
*>          On entry, the NB diagonal matrix D and the multipliers
*>          used to obtain the factor U or L as computed by ZSYTRF.
*>
*>          On exit, if INFO = 0, the (symmetric) inverse of the original
*>          matrix.  If UPLO = 'U', the upper triangular part of the
*>          inverse is formed and the part of A below the diagonal is not
*>          referenced; if UPLO = 'L' the lower triangular part of the
*>          inverse is formed and the part of A above the diagonal is
*>          not referenced.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in] I1
*> \verbatim
*>          I1 is INTEGER
*>          Index of the first row to swap
*> \endverbatim
*>
*> \param[in] I2
*> \verbatim
*>          I2 is INTEGER
*>          Index of the second row to swap
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complex16SYauxiliary
*
*  =====================================================================
      SUBROUTINE ZSYSWAPR( UPLO, N, A, LDA, I1, I2)
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER        UPLO
      INTEGER          I1, I2, LDA, N
*     ..
*     .. Array Arguments ..
      COMPLEX*16       A( LDA, N )
*
*  =====================================================================
*
*     ..
*     .. Local Scalars ..
      LOGICAL            UPPER
      INTEGER            I
      COMPLEX*16         TMP
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           ZSWAP
*     ..
*     .. Executable Statements ..
*
      UPPER = LSAME( UPLO, 'U' )
      IF (UPPER) THEN
*
*         UPPER
*         first swap
*          - swap column I1 and I2 from I1 to I1-1
         CALL ZSWAP( I1-1, A(1,I1), 1, A(1,I2), 1 )
*
*          second swap :
*          - swap A(I1,I1) and A(I2,I2)
*          - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1
         TMP=A(I1,I1)
         A(I1,I1)=A(I2,I2)
         A(I2,I2)=TMP
*
         DO I=1,I2-I1-1
            TMP=A(I1,I1+I)
            A(I1,I1+I)=A(I1+I,I2)
            A(I1+I,I2)=TMP
         END DO
*
*          third swap
*          - swap row I1 and I2 from I2+1 to N
         DO I=I2+1,N
            TMP=A(I1,I)
            A(I1,I)=A(I2,I)
            A(I2,I)=TMP
         END DO
*
        ELSE
*
*         LOWER
*         first swap
*          - swap row I1 and I2 from I1 to I1-1
         CALL ZSWAP( I1-1, A(I1,1), LDA, A(I2,1), LDA )
*
*         second swap :
*          - swap A(I1,I1) and A(I2,I2)
*          - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1
          TMP=A(I1,I1)
          A(I1,I1)=A(I2,I2)
          A(I2,I2)=TMP
*
          DO I=1,I2-I1-1
             TMP=A(I1+I,I1)
             A(I1+I,I1)=A(I2,I1+I)
             A(I2,I1+I)=TMP
          END DO
*
*         third swap
*          - swap col I1 and I2 from I2+1 to N
          DO I=I2+1,N
             TMP=A(I,I1)
             A(I,I1)=A(I,I2)
             A(I,I2)=TMP
          END DO
*
      ENDIF
      END SUBROUTINE ZSYSWAPR

max_stars_repo_name: sergey-v-kuznetsov/lapack

max_issues_repo_name: sergey-v-kuznetsov/lapack

max_forks_repo_name: sergey-v-kuznetsov/lapack


text:
MODULE load_mod
 USE param_mod,ONLY: mnam
 USE ctrl_mod,ONLY: nmdof
 USE esets_mod,ONLY: nnode,snames,add_name,search_name!snhead,sntail,n_names
 USE history_mod,ONLY: contact
 USE shape_mod
 USE quadrature_mod
 IMPLICIT NONE
 !Derived type for load pressure
 TYPE press_load
  REAL(kind=8) :: pload              !load magnitud (pressure)
  REAL(kind=8) :: thres1,thres2      !define zone for pressure
 END TYPE press_load
 !Derived type for load contact (stent)
 TYPE contc_load
  REAL(kind=8) :: cnor,ctan       !norm and tang elastic coeff
  REAL(kind=8) :: radiu              !radius of stent-graft
  REAL(kind=8) :: erfun              !factor for error function
  REAL(kind=8) :: thres1,thres2      !define zone for contact
  REAL(kind=8) :: dz,dr      !force direction
  REAL(kind=8) :: friction   !friction coefficient
 END TYPE contc_load
 !list of elements in a load set
 TYPE setload
  CHARACTER(len=mnam) :: sname              !load set name
  !pressure load
  TYPE(press_load) :: press         !pointer to pressure
  !contact load (stent)
  TYPE(contc_load) :: contc         !pointer to stent-graft
  TYPE(setload),POINTER :: next
 END TYPE setload
 !list of load sets
 TYPE list_load
   INTEGER(kind=4) :: nsload           !number of load sets
   TYPE(setload),POINTER :: head,tail  !first and last sets
 END TYPE list_load

 TYPE(list_load) :: loadlist

 CONTAINS
!******************************************************************************
!**                  Initialization of set load list                         **
!******************************************************************************
 SUBROUTINE ini_list_load(list)
  TYPE(list_load),INTENT(INOUT) :: list
   NULLIFY(list%head,list%tail)
   list%nsload = 0
 END SUBROUTINE ini_list_load
!******************************************************************************
!**               Add load set to the end of loadlist                        **
!******************************************************************************
 SUBROUTINE add_ldset_to_list(list,new)
  TYPE(list_load),INTENT(INOUT) :: list
  TYPE(setload),POINTER :: new
   IF (list%nsload==0) THEN
     list%head => new
     list%tail => new
     nullify(list%tail%next)
   ELSE
     list%tail%next => new
     nullify(new%next)
     list%tail => new
   END IF
   list%nsload = list%nsload + 1
 END SUBROUTINE add_ldset_to_list
!******************************************************************************
!**                 Initialization of load set (ldset)                       **
!******************************************************************************
 SUBROUTINE ini_load_set(list)
  TYPE(setload),POINTER :: list
   ALLOCATE(list)
   list%sname = ''
   CALL new_pressure(list%press)
   CALL new_contact(list%contc)
   NULLIFY(list%next)
 END SUBROUTINE ini_load_set
!******************************************************************************
!**              Add pressure to the end of load set                         **
!******************************************************************************
! SUBROUTINE add_press_to_ldset(new,head,tail)
!  TYPE(press_load),POINTER :: new,head,tail
!  !check if a list is empty
!  IF (.NOT.ASSOCIATED(head)) THEN
!     head => new
!     tail => new
!     nullify(tail%next)
!  ELSE
!     tail%next => new
!     nullify(new%next)
!     tail => new
!  END IF
!  !RETURN
! END SUBROUTINE add_press_to_ldset
!******************************************************************************
!**                 Initialization of pressure load                          **
!******************************************************************************
 SUBROUTINE new_pressure(new)
  TYPE(press_load),INTENT(INOUT) :: new
   new%pload = 0.d0        !initialize pressure
   new%thres1 = 0.d0       !first point for pressure zone
   new%thres2 = 0.d0       !second point for pressure zone
 END SUBROUTINE new_pressure
!******************************************************************************
!**              Add contact to the end of load set                          **
!******************************************************************************
! SUBROUTINE add_contc_to_ldset(new,head,tail)
!  TYPE(contc_load),POINTER :: new,head,tail
!  !check if a list is empty
!  IF (.NOT.ASSOCIATED(head)) THEN
!     head => new
!     tail => new
!     nullify(tail%next)
!  ELSE
!     tail%next => new
!     nullify(new%next)
!     tail => new
!  END IF
!  !RETURN
! END SUBROUTINE add_contc_to_ldset
!******************************************************************************
!**                 Initialization of contact load                           **
!******************************************************************************
 SUBROUTINE new_contact(new)
  TYPE(contc_load),INTENT(INOUT) :: new
   new%cnor = 0.d0        !normal elastic coefficient
   new%radiu = 0.d0        !initialize radius of stent
   new%erfun = 0.d0        !initialize factor for error function
   new%thres1 = 0.d0       !first point for stent zone
   new%thres2 = 0.d0       !second point for stent zone
   new%dz = 0.d0       !force in z direction
   new%dr = 0.d0       !force in r direction
   new%friction = 0.d0  !friction coefficient
   new%ctan = 0.d0      !tangential elastic coefficient
 END SUBROUTINE new_contact
!******************************************************************************
!**               Elemental Stiffnes matrix from loads                       **
!******************************************************************************
 SUBROUTINE loads_elem(time,coor,corp,cosd,disp,nodes,K_load,R_load)
  IMPLICIT NONE
  !--Dummy arguments
  INTEGER,INTENT(IN) :: nodes(nnode)
  REAL(kind=8),INTENT(IN) :: time(2),coor(2,nnode),corp(2,nnode),cosd(2),disp(3,nnode)
  REAL(kind=8),INTENT(OUT) :: K_load(3,3,nnode,nnode),R_load(3,nnode)
  !--Local arguments
  INTEGER(kind=4) :: isload,i1,i2,i3,i4
  REAL(kind=8) :: pressure,Cn,radius,delta,oversiz,factor,Ct,mu
  REAL(kind=8) :: pthres1,pthres2,cthres1,cthres2
  REAL(kind=8) :: pressp1,pressp2,VariLoadp
  REAL(kind=8) :: pressc1,pressc2,VariLoadc,direction(2)
  REAL(kind=8) :: tang1,tang2,VariTang1(2),VariTang2(2)
  REAL(kind=8) :: K_pres(2,2,nnode,nnode),R_pres(2,nnode)
  REAL(kind=8) :: K_cntc(2,2,nnode,nnode),R_cntc(2,nnode)
  REAL(kind=8) :: K_tang(2,2,nnode,nnode),R_tang(2,nnode)
  REAL(kind=8) :: coorc(2,nnode) !coorc(dof,node)
  REAL(kind=8) :: z1,z2,r1,r2,l,coseno,seno,z1p,z2p,r1p,r2p
  REAL(kind=8) :: delta_sliding(2)
  REAL(kind=8) :: time_thres
  TYPE(press_load) :: press
  TYPE(contc_load) :: stent
  TYPE(setload),POINTER :: ldset

!COORDINATES DEFINITION
  coorc(1,1) = coor(1,1) + disp(1,1)*cosd(2) - disp(2,1)*cosd(1)
  coorc(1,2) = coor(1,2) + disp(1,2)*cosd(2) - disp(2,2)*cosd(1)
  coorc(2,1) = coor(2,1) + disp(1,1)*cosd(1) + disp(2,1)*cosd(2)
  coorc(2,2) = coor(2,2) + disp(1,2)*cosd(1) + disp(2,2)*cosd(2)
  !current coordinates by node
  z1 = coorc(1,1)
  z2 = coorc(1,2)
  r1 = coorc(2,1)
  r2 = coorc(2,2)
  l = SQRT((r2-r1)**2.d0+(z2-z1)**2.d0)
  coseno = (z2-z1)/l
  seno = (r2-r1)/l

  !initialize matrices and vectors
  K_load = 0.d0
  R_load = 0.d0
  K_pres = 0.d0
  R_pres = 0.d0
  K_cntc = 0.d0
  R_cntc = 0.d0
  K_tang = 0.d0
  R_tang = 0.d0
  pressp1 = 0.d0
  pressp2 = 0.d0
  pressc1 = 0.d0
  pressc2 = 0.d0

  ldset => loadlist%head          !pointer to a load set in loadlist
  DO isload=1,loadlist%nsload     !CYCLE in loadlist (list of load sets)
  !PRESSURE DEFINITION
    press = ldset%press           !pressure load
      pressure = press%pload      !pressure value for this zone
      pthres1 = press%thres1
      pthres2 = press%thres2
    IF (pressure.GT.0.d0) THEN
  !PRESSURE BLOCK FOR LOADS SETS
      !check node 1 from element is inside the zone pressure
      IF ((coorc(1,1).GE.pthres1).AND.(coorc(1,1).LE.pthres2)) THEN
        pressp1 = pressure
      ELSE
        pressp1 = 0.d0
      END IF
      !check node 2 from element is inside the zone pressure
      IF ((coorc(1,2).GE.pthres1).AND.(coorc(1,2).LE.pthres2)) THEN
        pressp2 = pressure
      ELSE
        pressp2 = 0.d0
      END IF
    END IF

  !STENT DEFINITION
    time_thres = 40.d0/30.d0
    stent = ldset%contc           !contact load
      Cn = stent%cnor     !stent stiffness in this zone
      radius = stent%radiu        !stent radius in this zone
      factor = stent%erfun        !error function factor
      cthres1 = stent%thres1
      cthres2 = stent%thres2
      direction(1) = stent%dz
      direction(2) = stent%dr
      mu = stent%friction
      Ct = stent%ctan
    IF (Cn.GT.0.d0) THEN
  !CONTACT-SPRING BLOCK FOR STENT IN LOADS SETS
      !check node 1 from element is inside the zone stent-graft
      IF ((coorc(1,1).GE.cthres1).AND.(coorc(1,1).LE.cthres2)) THEN
        oversiz = radius-coor(2,1)
        delta = radius-coorc(2,1)
        IF (delta.GE.0.d0) THEN
          pressc1 = Cn*delta
        ELSE
          pressc1 = 0.d0
        END IF
        !STENT AGAINST PRESSURE
        IF ((coorc(1,1).GE.pthres1).AND.(coorc(1,1).LE.pthres2)) &
                & pressp1 = 0.5d0*pressure*(ERF(factor*5.d0*(delta/oversiz &
                & - 0.5d0/factor)) + 1.d0) &
                & - 0.5d0*pressure*(ERF(factor*5.d0*(delta/oversiz &
                & - 0.5d0/factor)) - 1.d0)*EXP(-(time(1)-1.d0)/time_thres)
      ELSE
        pressc1 = 0.d0
      END IF
      !check node 2 from element is inside the zone stent-graft
      IF ((coorc(1,2).GE.cthres1).AND.(coorc(1,2).LE.cthres2)) THEN
        oversiz = radius-coor(2,2)
        delta = radius-coorc(2,2)
        IF (delta.GE.0.d0) THEN
          pressc2 = Cn*delta
        ELSE
          pressc2 = 0.d0
        END IF
        !STENT AGAINST PRESSURE
        IF ((coorc(1,2).GE.pthres1).AND.(coorc(1,2).LE.pthres2)) &
                & pressp2 = 0.5d0*pressure*(ERF(factor*5.d0*(delta/oversiz &
                & - 0.5d0/factor))+1.d0) &
                & - 0.5d0*pressure*(ERF(factor*5.d0*(delta/oversiz &
                & - 0.5d0/factor)) - 1.d0)*EXP(-(time(1)-1.d0)/time_thres)
      ELSE
        pressc2 = 0.d0
      END IF
    END IF

    !GET LOAD FOR DIFERENT CASES
      !follower load for pressure
    IF (pressure.GT.0.d0) THEN
      VariLoadp = 0.d0
      CALL follower_load(pressp1,pressp2,VariLoadp,coorc,K_pres,R_pres)
    END IF
      !surface load for stent
    IF (Cn.GT.0.d0) THEN
      !previous coordinates by node
      z1p = corp(1,1)
      z2p = corp(1,2)
      r1p = corp(2,1)
      r2p = corp(2,2)
      !sliding
      delta_sliding(1) = (z1 - z1p)*coseno + (r1 - r1p)*seno
      delta_sliding(2) = (z2 - z2p)*coseno + (r2 - r2p)*seno
      !contact force
      VariLoadc = -Cn
      CALL surface_load(pressc1,pressc2,direction,VariLoadc,coorc,&
      & K_cntc,R_cntc)
      !node 1
      CALL slip_condition(mu,Ct,Cn,pressc1,coseno,seno,nodes(1),&
      & delta_sliding(1),tang1,VariTang1)
      !node 2
      CALL slip_condition(mu,Ct,Cn,pressc2,coseno,seno,nodes(2),&
      & delta_sliding(2),tang2,VariTang2)
      !tangential force
      CALL tangential_load(tang1,tang2,VariTang1,VariTang2,coorc,K_tang,R_tang)
    END IF
    ldset => ldset%next
  END DO

  DO i1=1,2
    DO i2=1,2
      DO i3=1,2
        DO i4=1,2
          K_load(i1,i2,i3,i4) = K_pres(i1,i2,i3,i4) + K_cntc(i1,i2,i3,i4) &
          & + K_tang(i1,i2,i3,i4)
        END DO
      END DO
    END DO
  END DO

  DO i1=1,2
    DO i3=1,2
      R_load(i1,i3) = R_pres(i1,i3) + R_cntc(i1,i3) + R_tang(i1,i3)
    END DO
  END DO
 END SUBROUTINE loads_elem
!******************************************************************************
!**                        FOLLOWER LOAD                                     **
!******************************************************************************
SUBROUTINE follower_load(press1,press2,DpDw,coorc,K_foll,R_foll)
 IMPLICIT NONE
 !**dummy arguments
 REAL(kind=8),INTENT(IN) :: press1,press2,coorc(2,2),DpDw
 REAL(kind=8),INTENT(OUT) :: K_foll(2,2,nnode,nnode),R_foll(2,nnode)
 !**local arguments
 REAL(kind=8) :: p(2),z(2),r(2)
 REAL(kind=8) :: K_DpDw(2,2,nnode,nnode)

 !pressure by node
 p(1) = press1
 p(2) = press2
 !coordinates by node
 z(1) = coorc(1,1)
 z(2) = coorc(1,2)
 r(1) = coorc(2,1)
 r(2) = coorc(2,2)

 CALL gauscheb(ngauss,gpoin,gcoef)   !gauss points and coefficients
 K_foll = 0.d0    !Stiffness matrix for constant face-load
 K_DpDw = 0.d0    !Stiffness matrix for Dp (variational load)
 R_foll = 0.d0    !Force vector from C*w
 DO igaus=1,ngauss          !gauss quadrature
 !shape functions
 CALL shape_1d_select(gpoin(igaus),Nfunc,dNds)
!----------------------------
 K_foll_g = 0.d0   !Stiffness matrix for constant face-load
 K_DpDw_g = 0.d0   !Stiffness matrix for Dp (variational load)
 R_foll_g = 0.d0   !Force vector from C*w
 !This cycle is to build stiffness matrix and residual vector by gauss point
 DO i1=1,nnode
   DO i2=1,nnode
     DO i3=1,nnode
       DO i4=1,nnode
       !Stiffness matrix for constant face-load
         K_foll_g(1,1,i1,i2) = 0.d0
         K_foll_g(1,2,i1,i2) = -Nfunc(i1)*p(i3)*Nfunc(i3)*r(i4)*(Nfunc(i4)*dNds(i2)&
                             & + dNds(i4)*Nfunc(i2)) + K_foll_g(1,2,i1,i2)
         K_foll_g(2,1,i1,i2) = Nfunc(i1)*p(i3)*Nfunc(i3)*r(i4)*Nfunc(i4)*dNds(i2)&
                             & + K_foll_g(2,1,i1,i2)
         K_foll_g(2,2,i1,i2) = Nfunc(i1)*p(i3)*Nfunc(i3)*z(i4)*dNds(i4)*Nfunc(i2)&
                             & + K_foll_g(2,2,i1,i2)
       !-----------------------
       !Stiffness matrix for Dp (variational load)
         K_DpDw_g(1,1,i1,i2) = 0.d0
         K_DpDw_g(1,2,i1,i2) = -Nfunc(i1)*Nfunc(i2)*r(i3)*Nfunc(i3)*r(i4)*dNds(i4)
         K_DpDw_g(2,1,i1,i2) = 0.d0
         K_DpDw_g(2,2,i1,i2) = Nfunc(i1)*Nfunc(i2)*r(i3)*Nfunc(i3)*z(i4)*dNds(i4)
       !-----------------------
        !Force vector from C*w
         R_foll_g(1,i1) = -Nfunc(i1)*p(i2)*Nfunc(i2)*r(i3)&
                        &*Nfunc(i3)*r(i4)*dNds(i4) + R_foll(1,i1)
         R_foll_g(2,i1) = Nfunc(i1)*p(i2)*Nfunc(i2)*r(i3)&
                        &*Nfunc(i3)*z(i4)*dNds(i4) + R_foll(2,i1)
       END DO
     END DO
   END DO
 END DO
 !Finally Gauss integration
 K_foll = gcoef(igaus)*K_foll_g + K_foll
 K_Dpdw = gcoef(igaus)*K_DpDw_g + K_DpDw
 R_foll = gcoef(igaus)*R_foll_g + R_foll
!---------------------------
 END DO                  !gauss quadrature

 !Stiffness matrix for constant face-load
! K_foll(1,1,1,1) = 0.d0
! K_foll(1,2,1,1) = 2.d0*(r1*(3.d0*p1+p2)-r2*p1)
! K_foll(2,1,1,1) = -r1*(3.d0*p1+p2)-r2*(p1+p2)
! K_foll(2,2,1,1) = (z2-z1)*(3.d0*p1+p2)

! K_foll(1,1,1,2) = 0.d0
! K_foll(1,2,1,2) = -2.d0*(r1*p1+r2*(p1+p2))
! K_foll(2,1,1,2) = r1*(3.d0*p1+p2)+r2*(p1+p2)
! K_foll(2,2,1,2) = (z2-z1)*(p1+p2)

! K_foll(1,1,2,1) = 0.d0
! K_foll(1,2,2,1) = 2.d0*(r1*(p1+p2)+r2*p2)
! K_foll(2,1,2,1) = -r1*(p1+p2)-r2*(p1+3.d0*p2)
! K_foll(2,2,2,1) = (z2-z1)*(p1+p2)

! K_foll(1,1,2,2) = 0.d0
! K_foll(1,2,2,2) = 2.d0*(r1*p2-r2*(p1+3.d0*p2))
! K_foll(2,1,2,2) = r1*(p1+p2)+r2*(p1+3.d0*p2)
! K_foll(2,2,2,2) = (z2-z1)*(p1+3.d0*p2)

 !Stiffness matrix for Dp (variational load)
! K_DpDw(1,1,1,1) = 0.d0
! K_Dpdw(1,2,1,1) = -(r2-r1)*(3.d0*r1+r2)
! K_DpDw(2,1,1,1) = 0.d0
! K_DpDw(2,2,1,1) = (z2-z1)*(3.d0*r1+r2)

! K_DpDw(1,1,1,2) = 0.d0
! K_DpDw(1,2,1,2) = -(r2-r1)*(r1+r2)
! K_DpDw(2,1,1,2) = 0.d0
! K_DpDw(2,2,1,2) = (z2-z1)*(r1+r2)

! K_DpDw(1,1,2,1) = 0.d0
! K_DpDw(1,2,2,1) = -(r2-r1)*(r1+r2)
! K_DpDw(2,1,2,1) = 0.d0
! K_DpDw(2,2,2,1) = (z2-z1)*(r1+r2)

! K_DpDw(1,1,2,2) = 0.d0
! K_DpDw(1,2,2,2) = -(r2-r1)*(r1+3.d0*r2)
! K_DpDw(2,1,2,2) = 0.d0
! K_DpDw(2,2,2,2) = (z2-z1)*(r1+3.d0*r2)

 !Force vector from C*w
! R_foll(1,1) = (r1-r2)*(3.d0*p1*r1+p1*r2+p2*r1+p2*r2)
! R_foll(2,1) = (z2-z1)*(3.d0*p1*r1+p1*r2+p2*r1+p2*r2)
! R_foll(1,2) = (r1-r2)*(p1*r1+p1*r2+p2*r1+3.d0*p2*r2)
! R_foll(2,2) = (z2-z1)*(p1*r1+p1*r2+p2*r1+3.d0*p2*r2)

 !set pressure to vector and matrix
 R_foll = R_foll/12.d0
 K_foll = (K_foll + DpDw*K_DpDw)/12.d0
END SUBROUTINE follower_load
!******************************************************************************
!**                        SURFACE LOAD                                     **
!******************************************************************************
SUBROUTINE surface_load(p1,p2,vec,DpDw,coorc,K_surf,R_surf)
 IMPLICIT NONE
 !**dummy arguments
 REAL(kind=8),INTENT(IN) :: p1,p2,vec(2),coorc(2,2),DpDw
 REAL(kind=8),INTENT(OUT) :: K_surf(2,2,2,2),R_surf(2,2)
 !**local arguments
 REAL(kind=8) :: factor1,factor2
 REAL(kind=8) :: z1,z2,r1,r2,l,v1,v2
 REAL(kind=8) :: K_DpDw(2,2,2,2)

 !direction vector
 v1 = vec(1)
 v2 = vec(2)
 !coordinates by node
 z1 = coorc(1,1)
 z2 = coorc(1,2)
 r1 = coorc(2,1)
 r2 = coorc(2,2)
 l = SQRT((r2-r1)**2.d0+(z2-z1)**2.d0)

 !Stiffness matrix for constant surf-load
 factor1 = (3.d0*p1*r1 + p1*r2 + p2*r1 + p2*r2)/l
 factor2 = (p1*r1 + p1*r2 + p2*r1 + 3.d0*p2*r2)/l

 K_surf(1,1,1,1) = 0.d0 - v1*(z2-z1)*factor1
 K_surf(1,2,1,1) = v1*l*(3.d0*p1+p2) - v1*(r2-r1)*factor1
 K_surf(2,1,1,1) = 0.d0 - v2*(z2-z1)*factor1
 K_surf(2,2,1,1) = v2*l*(3.d0*p1+p2) - v2*(r2-r1)*factor1

 K_surf(1,1,1,2) = 0.d0 + v1*(z2-z1)*factor1
 K_surf(1,2,1,2) = v1*l*(p1+p2) + v1*(r2-r1)*factor1
 K_surf(2,1,1,2) = 0.d0 + v2*(z2-z1)*factor1
 K_surf(2,2,1,2) = v2*l*(p1+p2) + v2*(r2-r1)*factor1

 K_surf(1,1,2,1) = 0.d0 - v1*(z2-z1)*factor2
 K_surf(1,2,2,1) = v1*l*(p1+p2) - v1*(r2-r1)*factor2
 K_surf(2,1,2,1) = 0.d0 - v2*(z2-z1)*factor2
 K_surf(2,2,2,1) = v2*l*(p1+p2) - v2*(r2-r1)*factor2

 K_surf(1,1,2,2) = 0.d0 + v1*(z2-z1)*factor2
 K_surf(1,2,2,2) = v1*l*(p1+3.d0*p2) + v1*(r2-r1)*factor2
 K_surf(2,1,2,2) = 0.d0 + v2*(z2-z1)*factor2
 K_surf(2,2,2,2) = v2*l*(p1+3.d0*p2) + v2*(r2-r1)*factor2

 !Stiffness matrix for Dp (variational load)
 K_DpDw(1,1,1,1) = 0.d0
 K_Dpdw(1,2,1,1) = v1*(3.d0*r1+r2)
 K_DpDw(2,1,1,1) = 0.d0
 K_DpDw(2,2,1,1) = v2*(3.d0*r1+r2)

 K_DpDw(1,1,1,2) = 0.d0
 K_DpDw(1,2,1,2) = v1*(r1+r2)
 K_DpDw(2,1,1,2) = 0.d0
 K_DpDw(2,2,1,2) = v2*(r1+r2)

 K_DpDw(1,1,2,1) = 0.d0
 K_DpDw(1,2,2,1) = v1*(r1+r2)
 K_DpDw(2,1,2,1) = 0.d0
 K_DpDw(2,2,2,1) = v2*(r1+r2)

 K_DpDw(1,1,2,2) = 0.d0
 K_DpDw(1,2,2,2) = v1*(r1+3.d0*r2)
 K_DpDw(2,1,2,2) = 0.d0
 K_DpDw(2,2,2,2) = v2*(r1+3.d0*r2)

 !Force vector from C*w
 R_surf(1,1) = v1*l**2.d0*factor1
 R_surf(2,1) = v2*l**2.d0*factor1
 R_surf(1,2) = v1*l**2.d0*factor2
 R_surf(2,2) = v2*l**2.d0*factor2

 !set pressure to vector and matrix
 R_surf = R_surf/12.d0
 K_surf = (K_surf + DpDw*l*K_DpDw)/12.d0
END SUBROUTINE surface_load
!******************************************************************************
!**                 TANGENTIAL CONTACT LOAD                                  **
!******************************************************************************
SUBROUTINE tangential_load(t1,t2,Dt1,Dt2,coorc,K_tang,R_tang)
 IMPLICIT NONE
 !**dummy arguments
 REAL(kind=8),INTENT(IN) :: t1,t2,coorc(2,2),Dt1(2),Dt2(2)
 REAL(kind=8),INTENT(OUT) :: K_tang(2,2,2,2),R_tang(2,2)
 !**local arguments
 REAL(kind=8) :: factor1,factor2
 REAL(kind=8) :: z1,z2,r1,r2
 REAL(kind=8) :: K_Dtan(2,2,2,2)

 !coordinates by node
 z1 = coorc(1,1)
 z2 = coorc(1,2)
 r1 = coorc(2,1)
 r2 = coorc(2,2)

 !Stiffness matrix for constant surf-load
 factor1 = (3.d0*t1*r1 + t1*r2 + t2*r1 + t2*r2)
 factor2 = (t1*r1 + t1*r2 + t2*r1 + 3.d0*t2*r2)

 K_tang(1,1,1,1) = 0.d0 - factor1
 K_tang(1,2,1,1) = (z2-z1)*(3.d0*t1+t2)
 K_tang(2,1,1,1) = 0.d0
 K_tang(2,2,1,1) = (r2-r1)*(3.d0*t1+t2) - factor1

 K_tang(1,1,1,2) = 0.d0 + factor1
 K_tang(1,2,1,2) = (z2-z1)*(t1+t2)
 K_tang(2,1,1,2) = 0.d0
 K_tang(2,2,1,2) = (r2-r1)*(t1+t2) + factor1

 K_tang(1,1,2,1) = 0.d0 - factor2
 K_tang(1,2,2,1) = (z2-z1)*(t1+t2)
 K_tang(2,1,2,1) = 0.d0
 K_tang(2,2,2,1) = (r2-r1)*(t1+t2) - factor2

 K_tang(1,1,2,2) = 0.d0 + factor2
 K_tang(1,2,2,2) = (z2-z1)*(t1+3.d0*t2)
 K_tang(2,1,2,2) = 0.d0
 K_tang(2,2,2,2) = (r2-r1)*(t1+3.d0*t2) + factor2

 !Stiffness matrix for DT (variational load) - SLIDING
 K_Dtan(1,1,1,1) = (z2-z1)*(3.d0*r1+r2)*Dt1(1)
 K_Dtan(1,2,1,1) = (z2-z1)*(3.d0*r1+r2)*Dt1(2)
 K_Dtan(2,1,1,1) = (r2-r1)*(3.d0*r1+r2)*Dt1(1)
 K_Dtan(2,2,1,1) = (r2-r1)*(3.d0*r1+r2)*Dt1(2)

 K_Dtan(1,1,1,2) = (z2-z1)*(r1+r2)*Dt2(1)
 K_Dtan(1,2,1,2) = (z2-z1)*(r1+r2)*Dt2(2)
 K_Dtan(2,1,1,2) = (r2-r1)*(r1+r2)*Dt2(1)
 K_Dtan(2,2,1,2) = (r2-r1)*(r1+r2)*Dt2(2)

 K_Dtan(1,1,2,1) = (z2-z1)*(r1+r2)*Dt1(1)
 K_Dtan(1,2,2,1) = (z2-z1)*(r1+r2)*Dt1(2)
 K_Dtan(2,1,2,1) = (r2-r1)*(r1+r2)*Dt1(1)
 K_Dtan(2,2,2,1) = (r2-r1)*(r1+r2)*Dt1(2)

 K_Dtan(1,1,2,2) = (z2-z1)*(r1+3.d0*r2)*Dt2(1)
 K_Dtan(1,2,2,2) = (z2-z1)*(r1+3.d0*r2)*Dt2(2)
 K_Dtan(2,1,2,2) = (r2-r1)*(r1+3.d0*r2)*Dt2(1)
 K_Dtan(2,2,2,2) = (r2-r1)*(r1+3.d0*r2)*Dt2(2)

 !Force vector from C*w
 R_tang(1,1) = (z2-z1)*factor1
 R_tang(2,1) = (r2-r1)*factor1
 R_tang(1,2) = (z2-z1)*factor2
 R_tang(2,2) = (r2-r1)*factor2

 !set pressure to vector and matrix
 R_tang = R_tang/12.d0
 K_tang = (K_tang + K_Dtan)/12.d0
END SUBROUTINE tangential_load
!******************************************************************************
!**                          SLIP CONDITION                                  **
!******************************************************************************
SUBROUTINE slip_condition(friction,Ct,Cn,pres,coseno,seno,nodes,delta_slid,&
  & tang,Dtang)
 IMPLICIT NONE
 !**dummy arguments
 INTEGER,INTENT(IN) :: nodes
 REAL(kind=8),INTENT(IN) :: friction,Ct,Cn,pres,coseno,seno,delta_slid
 REAL(kind=8),INTENT(OUT) :: Dtang(2),tang
 !**local arguments
 REAL(kind=8) :: pn,force
 REAL(kind=8) :: delta_stick,stick,sliding

 !get sliding and stick from history
 sliding = contact(nodes,1)
 stick = contact(nodes,2)
 !friction threshold force
 pn = pres*coseno            !normal pressure to element
 force = friction*pn               !threshold friction force
 !total sliding
 sliding = sliding + delta_slid     !add delta to previous sliding
 !trial Force
 tang = Ct*(sliding - stick)     !get tangential force with actual sliding
 !Slip condition
 IF (tang.LE.force) THEN           !NO SLIDING
   tang = tang                !trial force is the force
   stick = stick              !no evolution of stick-plastic displacement
   Dtang(1) = Ct*coseno       !variational for tangential force
   Dtang(2) = Ct*seno         !variational for tangential force
 ELSE                              !SLIDING
   delta_stick = (tang - force)/Ct
   tang = force                  !threhold force is the force
   stick = stick + delta_stick   !evolution of stick-plastic displacement
   Dtang(1) = 0.d0                   !variational of tangential force
   Dtang(2) = -friction*Cn*coseno    !variational of tangential force
 END IF
 !update history
 contact(nodes,1) = sliding
 contact(nodes,2) = stick
END SUBROUTINE slip_condition

END MODULE load_mod

max_stars_repo_name: jdlaubrie/shell-elem

max_issues_repo_name: jdlaubrie/shell-elem

max_forks_repo_name: jdlaubrie/shell-elem


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE BART ( OPT, L, AREA, I1, I2, JTOR, SCOEFF, K1, K2, I12, E, G, ALPHA, TREF )

! Calculates, for 1-D Timoshenko beam element

!  1) PTE       = element thermal load vectors         , if OPT(2) = 'Y'
!  2) SEi, STEi = element stress data recovery matrices, if OPT(3) = 'Y'
!  3) KE        = element linear stiffness matrix      , if OPT(6) = 'N' (i.e. always calc KE linear unless OPT(6) = 'Y')
!  4) KED       = element differen stiff matrix        , if OPT(6) = 'Y'

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  FATAL_ERR, NTSUB, BLNK_SUB_NAM, SOL_NAME
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BAR1_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO, ONE, TWO, THREE, FOUR, FIVE, SIX, TEN, TWELVE
      USE DEBUG_PARAMETERS
      USE PARAMS, ONLY                :  EPSIL, ART_KED, ART_ROT_KED, ART_TRAN_KED
      USE NONLINEAR_PARAMS, ONLY      :  LOAD_ISTEP
      USE MODEL_STUF, ONLY            :  ELDOF, DOFPIN, DT, EID, NUM_EMG_FATAL_ERRS, KE, KED, PEL, PTE, SE1, SE2, STE1, STE2, TYPE,&
                                         UEL
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG

      USE BAR1_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BART'
      CHARACTER(1*BYTE), INTENT(IN)   :: OPT(6)            ! 'Y'/'N' flags for whether to calc certain elem matrices

      INTEGER(LONG)                   :: I,J               ! DO loop induces
      INTEGER(LONG)                   :: IERROR            ! Local error indicator
      INTEGER(LONG)                   :: NUM_PFLAG_DOFS    ! The number of pin flagged DOF's for this element
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BAR1_BEGEND

      REAL(DOUBLE) , INTENT(IN)       :: ALPHA             ! Coefficient of thermal expansion
      REAL(DOUBLE) , INTENT(IN)       :: AREA              ! Cross-sectional area
      REAL(DOUBLE) , INTENT(IN)       :: E                 ! Youngs modulus
      REAL(DOUBLE) , INTENT(IN)       :: G                 ! Shear modulus
      REAL(DOUBLE) , INTENT(IN)       :: I1                ! Bending inertia in plane 1
      REAL(DOUBLE) , INTENT(IN)       :: I12               ! Product of inertia
      REAL(DOUBLE) , INTENT(IN)       :: I2                ! Bending inertia in plane 2
      REAL(DOUBLE) , INTENT(IN)       :: JTOR              ! Torsional constant
      REAL(DOUBLE) , INTENT(IN)       :: K1                ! Shear constant for plane 1 (used in K1*AREA*G)
      REAL(DOUBLE) , INTENT(IN)       :: K2                ! Shear constant for plane 2 (used in K1*AREA*G)
      REAL(DOUBLE) , INTENT(IN)       :: L                 ! Elem length
      REAL(DOUBLE) , INTENT(IN)       :: SCOEFF            ! Stress recovery coeff for torsion
      REAL(DOUBLE) , INTENT(IN)       :: TREF              ! Element reference temperature
      REAL(DOUBLE)                    :: ABAR(6,5)         ! Thermal displacements per unit TPRIME
      REAL(DOUBLE)                    :: B1(3,6)           ! Intermediate array used in calculating SE1 stress matrix
      REAL(DOUBLE)                    :: B2(3,6)           ! Intermediate array used in calculating SE2 stress matrix
      REAL(DOUBLE)                    :: BT1(3,5)          ! Intermediate array used in calculating STE1 thermal stress effects
      REAL(DOUBLE)                    :: BT2(3,5)          ! Intermediate array used in calculating STE2 thermal stress effects
      REAL(DOUBLE)                    :: BTA(6,5)          ! = KAA*ABAR. Each col gives thermal loads at end a per unit TPRIME value
      REAL(DOUBLE)                    :: BTB(6,5)          ! = KBA*ABAR. Each col gives thermal loads at end a per unit TPRIME value
      REAL(DOUBLE)                    :: C01               ! Intermediate variable used in calc terms for the stiff matrix, KE
      REAL(DOUBLE)                    :: DELTA1            ! I2/DEN
      REAL(DOUBLE)                    :: DELTA2            ! I1/DEN
      REAL(DOUBLE)                    :: DELTA12           ! I12/DEN
      REAL(DOUBLE)                    :: DEN               ! I1*I2 - I12**2. If DEN = 0 cannot calc stress rec matrices (SEi, STEi)
      REAL(DOUBLE)                    :: DUM1(3,NTSUB)     ! Intermediate matrix
      REAL(DOUBLE)                    :: DUM2(3,NTSUB)     ! Intermediate matrix
      REAL(DOUBLE)                    :: EPS1              ! A small number to compare for real zero
      REAL(DOUBLE)                    :: G1                ! Intermediate variable used in calc terms for the stiff matrix, KE
      REAL(DOUBLE)                    :: G2                ! Intermediate variable used in calc terms for the stiff matrix, KE
      REAL(DOUBLE)                    :: C1,C2,C3,C4,C5,C6,C7,C8,C9
      REAL(DOUBLE)                    :: KAA(6,6)          ! Upper left  6 x 6 partition of KE
      REAL(DOUBLE)                    :: KBA(6,6)          ! Lower left  6 x 6 partition of KE
      REAL(DOUBLE)                    :: KAB(6,6)          ! Upper right 6 x 6 partition of KE
      REAL(DOUBLE)                    :: PTA(6,NTSUB)      ! Thermal loads at end a of BAR (rows 1- 6 of PTE)
      REAL(DOUBLE)                    :: PTb(6,NTSUB)      ! Thermal loads at end b of BAR (rows 7-12 of PTE)
      REAL(DOUBLE)                    :: RG                ! Intermediate variable used in calc terms for the stiff matrix, KE
      REAL(DOUBLE)                    :: S11(3,6)          ! Intermediate matrix used in calculating SE1 stress matrix
      REAL(DOUBLE)                    :: S12(3,6)          ! Intermediate matrix used in calculating SE1 stress matrix
      REAL(DOUBLE)                    :: S21(3,6)          ! Intermediate matrix used in calculating SE2 stress matrix
      REAL(DOUBLE)                    :: S22(3,6)          ! Intermediate matrix used in calculating SE2 stress matrix
      REAL(DOUBLE)                    :: TBAR              ! Average elem temperature

! The following are used for the differential stiffness matrix calc. See NASTRAN Prog's Manual (COSMIC 1972) page 4.87.30

      REAL(DOUBLE)                    :: M1a               ! Bend mom, plane 1, end a (also Maz, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: M2a               ! Bend mom, plane 2, end a (also May, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: M1b               ! Bend mom, plane 1, end b (also Mbz, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: M2b               ! Bend mom, plane 2, end b (also Mby, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: V1                ! Plane 1 shear (also Vy, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: V2                ! Plane 2 shear (also Vz, in NASTRAN Prog Man differ stiff)
      REAL(DOUBLE)                    :: Fx                ! Axial force (also Fx, in NASTRAN Prog Man differ stiff)

      INTRINSIC DABS

       REAL(DOUBLE)                    :: TPRIME(5,NTSUB)   ! Matrix where each col has the 5 temperature/gradients for the BAR elem

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      EPS1 = EPSIL(1)

! Calculate element stiffness matrix KE(12,12), regardless of OPT(4) since the calculation of PTE and/or SEi, STEi
! uses values from KE. Initialize KE here to make sure it starts as zero for the BAR element

      DO I=1,12
         DO J=1,12
            KE(I,J) = ZERO
         ENDDO
      ENDDO

      G1 = K1*AREA*G
      G2 = K2*AREA*G

      C1 =  AREA*E/L

      C2 =  G1/L
      C3 =  G1/2                                           ! NOTE: TESSLER HAS OPPOSITE SIGN BUT NEGATIVE, NOT POSITIVE, WORKS
      C4 =  E*I2/L + G1*L/4
      C5 = -E*I2/L + G1*L/4

      C6 =  G2/L
      C7 =  G2/2
      C8 =  E*I1/L + G2*L/4
      C9 = -E*I1/L + G2*L/4

      RG =  G*JTOR/L

      DEN    = I1*I2 - I12*I12                             ! I1*I2 > I12^2 was checked when PBAR's were read in subr LOADB
      DELTA1 = I2/DEN
      DELTA2 = I1/DEN
      DELTA12= 0.D0

      IF (DEBUG(203) > 0) CALL DEBUG_BART ( 1 )

! Generate KE

      KE( 1, 1) = C1
      KE( 1, 7) =-C1

      KE( 2, 2) = C6
      KE( 2, 6) = C7
      KE( 2, 8) =-C6
      KE( 2,12) = C7

      KE( 3, 3) = C2
      KE( 3, 5) =-C3
      KE( 3, 9) =-C2
      KE( 3,11) =-C3

      KE( 4, 4) = RG
      KE( 4,10) =-RG

      KE( 5, 5) = C4
      KE( 5, 9) = C3
      KE( 5,11) = C5

      KE( 6, 6) = C8
      KE( 6, 8) =-C7
      KE( 6,12) = C9

      KE( 7, 7) = C1

      KE( 8, 8) = C6
      KE( 8,12) =-C7

      KE( 9, 9) = C2
      KE( 9,11) = C3

      KE(10,10) = RG

      KE(11,11) = C4

      KE(12,12) = C8

      DO I=2,12                                            ! Set lower triangular partition of KE using symmetry
         DO J=1,I-1
            KE(I,J) = KE(J,I)
         ENDDO
      ENDDO

! Process Pin Flags. NUM_PFLAG_DOFS is a count of the total number of DOF pin flagged. DOFPIN(i) generated in ELMDAT is an
! integer array of the DOF numbers of the pin flagged DOF.

      NUM_PFLAG_DOFS = 0
      DO I=1,12
         IF (DOFPIN(I) > 0) THEN
            NUM_PFLAG_DOFS = NUM_PFLAG_DOFS + 1
         ENDIF
      ENDDO
      IF (NUM_PFLAG_DOFS /= 0) THEN
         CALL PINFLG ( NUM_PFLAG_DOFS )
      ENDIF

      DO I=1,6                                             ! Upper left partition of KE is KAA. Need this for calc'ing PTE, SEi
         DO J=1,6
            KAA(I,J) = KE(I,J)
         ENDDO
      ENDDO

      DO I=1,6                                             ! Upper right partition of KE is KAB. Need this for calculating SEi
         DO J=7,12
            KAB(I,J-6) = KE(I,J)
         ENDDO
      ENDDO

      DO I=7,12                                            ! Lower left partition of KE is KBA. Need this for calculating PTE
         DO J=1,6
            KBA(I-6,J) = KE(I,J)
         ENDDO
      ENDDO

! The following matrices are needed if either OPT(2) or OPT(3) or OPT(6) is called for

      IF ((OPT(2) == 'Y') .OR. (OPT(3) == 'Y') .OR. (OPT(6) == 'Y')) THEN
         IF (NTSUB > 0) THEN

            DO J=1,NTSUB
               TBAR        = (DT(1,J) + DT(2,J))/TWO
               TPRIME(1,J) = TBAR - TREF
               TPRIME(2,J) = DT(3,J)
               TPRIME(3,J) = DT(4,J)
               TPRIME(4,J) = DT(5,J)
               TPRIME(5,J) = DT(6,J)
            ENDDO

         ENDIF

      ENDIF

! **********************************************************************************************************************************
! Determine element thermal loads.

!     IF ((OPT(2) == 'Y') .OR. (OPT(6) == 'Y')) THEN

         IF (NTSUB > 0) THEN

            DO I=1,6
               DO J=1,5
                  ABAR(I,J) = ZERO
               ENDDO
            ENDDO

            ABAR(1,1) =  ONE
            ABAR(2,2) =  DELTA1*I1*L/SIX
            ABAR(2,3) =  DELTA1*I1*L/THREE
            ABAR(2,4) = -DELTA12*I2*L/SIX
            ABAR(2,5) = -DELTA12*I2*L/THREE
            ABAR(3,2) = -DELTA12*I1*L/SIX
            ABAR(3,3) = -DELTA12*I1*L/THREE
            ABAR(3,4) =  DELTA2*I2*L/SIX
            ABAR(3,5) =  DELTA2*I2*L/THREE
            ABAR(5,2) = -DELTA12*I1/TWO
            ABAR(5,3) = -DELTA12*I1/TWO
            ABAR(5,4) =  DELTA2*I2/TWO
            ABAR(5,5) =  DELTA2*I2/TWO
            ABAR(6,2) = -DELTA1*I1/TWO
            ABAR(6,3) = -DELTA1*I1/TWO
            ABAR(6,4) =  DELTA12*I2/TWO
            ABAR(6,5) =  DELTA12*I2/TWO

            DO I=1,6
               DO J=1,5
                  ABAR(I,J) = -ALPHA*L*ABAR(I,J)
               ENDDO
            ENDDO

            CALL MATMULT_FFF ( KAA, ABAR, 6, 6, 5, BTA )
            CALL MATMULT_FFF ( KBA, ABAR, 6, 6, 5, BTB )

            CALL MATMULT_FFF ( BTA, TPRIME, 6, 5, NTSUB, PTA )
            CALL MATMULT_FFF ( BTB, TPRIME, 6, 5, NTSUB, PTB )
            DO I=1,6
               DO J=1,NTSUB
                  PTE(I,J)   = PTA(I,J)
                  PTE(I+6,J) = PTB(I,J)
               ENDDO
            ENDDO

         ENDIF

!     ENDIF

! **********************************************************************************************************************************
! Calculate SE matrices for stress data recovery.

!     IF ((OPT(3) == 'Y') .OR. (OPT(6) == 'Y')) THEN

         DO I=1,3
            DO J=1,6
               B1(I,J) = ZERO
               B2(I,J) = ZERO
            ENDDO
         ENDDO

         IF (DABS(AREA) > EPS1) THEN                       ! Can't calculate axial stress if the area = 0
            B1(1,1) = -ONE/AREA
         ENDIF

         B1(2,5) = -DELTA12
         B1(2,6) = -DELTA1

         B1(3,5) =  DELTA2
         B1(3,6) =  DELTA12

         B2(1,2) =  DELTA1*L
         B2(1,3) = -DELTA12*L
         B2(1,5) = -DELTA12
         B2(1,6) = -DELTA1

         B2(2,2) = -DELTA12*L
         B2(2,3) =  DELTA2*L
         B2(2,5) =  DELTA2
         B2(2,6) =  DELTA12

         IF (DABS(JTOR) > EPS1) THEN                       ! Can't calc torsional stress if torsional constant = 0
            B2(3,4) = -SCOEFF/JTOR
         ENDIF

         CALL MATMULT_FFF ( B1, KAA, 3, 6, 6, S11 )
         CALL MATMULT_FFF ( B1, KAB, 3, 6, 6, S12 )
         CALL MATMULT_FFF ( B2, KAA, 3, 6, 6, S21 )
         CALL MATMULT_FFF ( B2, KAB, 3, 6, 6, S22 )

         DO I=1,3
            DO J=1,6
               SE1(I,J,1) = S11(I,J)
               SE2(I,J,1) = S21(I,J)
            ENDDO
            DO J=7,12
               SE1(I,J,1) = S12(I,J-6)
               SE2(I,J,1) = S22(I,J-6)
            ENDDO
         ENDDO

         IF (NTSUB > 0) THEN                               ! Calculate STEi for thermal stress data recovery

            DO I=1,3
               DO J=1,5
                  BT1(I,J) = ZERO
                  BT2(I,J) = ZERO
               ENDDO
            ENDDO

            CALL MATMULT_FFF ( S11, ABAR, 3, 6, 5, BT1 )
            CALL MATMULT_FFF ( S21, ABAR, 3, 6, 5, BT2 )

            CALL MATMULT_FFF ( BT1, TPRIME, 3, 5, NTSUB, DUM1 )
            CALL MATMULT_FFF ( BT2, TPRIME, 3, 5, NTSUB, DUM2 )
            DO I=1,3
               DO J=1,NTSUB
                  STE1(I,J,1) = DUM1(I,J)
                  STE2(I,J,1) = DUM2(I,J)
               ENDDO
            ENDDO

         ENDIF

!     ENDIF

! **********************************************************************************************************************************
! Calculate linear differential stiffness matrix

      IF ((OPT(6) == 'Y') .AND. (LOAD_ISTEP > 1)) THEN

         CALL ELMDIS

                                                           ! Calc BAR forces
         CALL CALC_ELEM_NODE_FORCES
         M1a = -PEL(6)                                     ! M1a (bending moment, plane 1, end a for BAR) - NASTRAN Maz
         M2a =  PEL(5)                                     ! M2a (bending moment, plane 2, end a for BAR) - NASTRAN May
         M1b = -PEL(6) + PEL(2)*L                          ! M1b (bending moment, plane 1, end b for BAR) - NASTRAN Mbz
         M2b =  PEL(5) + PEL(3)*L                          ! M2b (bending moment, plane 2, end b for BAR) - NASTRAN Mby
         V1  = -PEL(2)                                     ! V1  (plane 1 shear for BAR)                  - NASTRAN Vy
         V2  = -PEL(3)                                     ! V2  (plane 2 shear for BAR)                  - NASTRAN Vz
         Fx  = -PEL(1)                                     ! Fx  (axial force for BAR or ROD)             - NASTRAN Fx

         IF (ART_KED == 'Y') THEN
            KED( 1, 1) = ART_TRAN_KED
            KED( 4, 4) = ART_ROT_KED
            KED( 7, 7) = ART_TRAN_KED
            KED(10,10) = ART_ROT_KED
         ENDIF

         C01 = Fx/L

         KED( 2, 2) =  (SIX/FIVE)*Fx/L
         KED( 2, 4) =  M2b/L
         KED( 2, 6) =  Fx/TEN
         KED( 2, 8) = -KED( 2, 2)
         KED( 2,10) =  M2a/L
         KED( 2,12) =  KED( 2, 6)

         KED( 3, 3) =  KED( 2, 2)
         KED( 3, 4) =  M1b/L
         KED( 3, 5) = -KED( 2, 6)
         KED( 3, 9) = -KED( 2, 2)
         KED( 3,10) =  M1a/L
         KED( 3,11) = -KED( 2, 6)

         KED( 4, 4) =  JTOR*Fx/(L*L*AREA)
         KED( 4, 5) = -V1*L/SIX
         KED( 4, 6) = -V2*L/SIX
         KED( 4, 8) = -M2b/L
         KED( 4, 9) = -M1b/L
         KED( 4,10) = -KED( 4, 4)
         KED( 4,11) = -KED( 4, 5)
         KED( 4,12) = -KED( 4, 6)

         KED( 5, 5) =  FOUR*Fx*L/(THREE*TEN)
         KED( 5, 9) =  KED( 2, 6)
         KED( 5,10) = -KED( 4, 5)
         KED( 5,11) = -Fx*L/(THREE*TEN)

         KED( 6, 6) =  KED( 5, 5)
         KED( 6, 8) = -KED( 2, 6)
         KED( 6,10) = -KED( 4, 6)
         KED( 6,12) =  KED( 5,11)

         KED( 8, 8) =  KED( 2, 2)
         KED( 8,10) = -KED( 2,10)
         KED( 8,12) = -KED( 2, 6)

         KED( 9, 9) =  KED( 2, 2)
         KED( 9,10) = -KED( 3,10)
         KED( 9,11) =  KED( 2, 6)

         KED(10,10) =  KED( 4, 4)
         KED(10,11) =  KED( 4, 5)
         KED(10,12) =  KED( 4, 6)

         KED(11,11) =  KED( 5, 5)

         KED(12,12) =  KED( 5, 5)

         DO I=2,12
            DO J=1,I-1
               KED(I,J) = KED(J,I)
            ENDDO
         ENDDO

         IF (DEBUG(178) >= 1) THEN
            WRITE(F06,2001) TRIM(SUBR_NAME), OPT(6), LOAD_ISTEP, TYPE, EID
            DO I=1,12
               WRITE(F06,2101) (KED(I,J),J=1,12)
            ENDDO
            WRITE(F06,*)
         ENDIF

      ENDIF

! If this is DIFFEREN or NLSTATIC, add linear and differential stiffness matrices

      IF ((SOL_NAME(1:8) == 'DIFFEREN') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN

         DO I=1,ELDOF
            DO J=1,ELDOF
               KE(I,J) = KE(I,J) + KED(I,J)
            ENDDO
         ENDDO

         IF (DEBUG(178) >= 1) THEN
            IF ((OPT(6) == 'Y') .AND. (LOAD_ISTEP > 1)) THEN
               WRITE(F06,2002) TRIM(SUBR_NAME), OPT(6), LOAD_ISTEP, TYPE, EID
            ELSE
               WRITE(F06,2003) TRIM(SUBR_NAME), OPT(6), LOAD_ISTEP, TYPE, EID
            ENDIF
            DO I=1,12
               WRITE(F06,2101) (KE(I,J),J=1,12)
            ENDDO
            WRITE(F06,*)
         ENDIF

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1918 FORMAT(' *ERROR  1918: NEGATIVE SHEAR FACTOR, K1 OR K2, ON ',A,' ELEMENT ',I8,' IS CAUSING A SINGULAR STIFFNESS')

 2001 FORMAT(' In ', A, ' with OPT(6) = ', A, ', and LOAD_ISTEP = ', I8, ': Differential Stiffness Matrix for ', A,                &
             ' element number ', I8)

 2002 FORMAT(' In ', A, ' with OPT(6) = ', A, ', and LOAD_ISTEP = ', I8, ': Linear + Differential Stiffness Matrix for ', A,       &
             ' element number ', I8)

 2003 FORMAT(' In ', A, ' with OPT(6) = ', A, ', and LOAD_ISTEP = ', I8, ': Linear Stiffness Matrix for ', A,' element number ', I8)

 2101 FORMAT(12(1ES14.6))

! ##################################################################################################################################

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE DEBUG_BART (WHAT)

      USE PENTIUM_II_KIND
      USE IOUNT1, ONLY                :  ERR, F04, F06

      IMPLICIT NONE

      INTEGER(LONG), INTENT(IN)       :: WHAT

! **********************************************************************************************************************************
      IF (WHAT == 1) THEN

         WRITE(F06,*)
         WRITE(F06,1997)
         WRITE(F06,'(A,I8)') 'In subr BART with BAR element ',EID
         WRITE(F06,*) '-------------------------------------'
         WRITE(F06,1998) 'L       = ',L
         WRITE(F06,1998) 'AREA    = ',AREA
         WRITE(F06,1998) 'I1      = ',I1
         WRITE(F06,1998) 'I2      = ',I2
         WRITE(F06,1998) 'I12     = ',I12
         WRITE(F06,1998) 'JTOR    = ',JTOR
         WRITE(F06,1998) 'K1      = ',K1
         WRITE(F06,1998) 'K2      = ',K2
         WRITE(F06,1998) 'E       = ',E
         WRITE(F06,1998) 'G       = ',G
         WRITE(F06,1998) 'SCOEFF  = ',SCOEFF
         WRITE(F06,1998) 'ALPHA   = ',ALPHA
         WRITE(F06,1998) 'TREF    = ',TREF
         WRITE(F06,1998) 'G1      = ',G1
         WRITE(F06,1998) 'G2      = ',G2
         WRITE(F06,1998) 'C1      = ',C1
         WRITE(F06,1998) 'C2      = ',C2
         WRITE(F06,1998) 'C3      = ',C3
         WRITE(F06,1998) 'C4      = ',C4
         WRITE(F06,1998) 'C5      = ',C5
         WRITE(F06,1998) 'C6      = ',C6
         WRITE(F06,1998) 'C7      = ',C7
         WRITE(F06,1998) 'C8      = ',C8
         WRITE(F06,1998) 'C9      = ',C9
         WRITE(F06,1998) 'RG      = ',RG
         WRITE(F06,1998) 'DEN     = ',DEN
         WRITE(F06,1998) 'DELTA1  = ',DELTA1
         WRITE(F06,1998) 'DELTA2  = ',DELTA2
         WRITE(F06,1998) 'DELTA12 = ',DELTA12

      ELSE IF (WHAT == 2) THEN

      ENDIF

! **********************************************************************************************************************************
 1997 FORMAT('************************************************************')

 1998 FORMAT(A, 1ES14.6)

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

      END SUBROUTINE DEBUG_BART

      END SUBROUTINE BART


text:
module init_data_module

contains

  subroutine init_data_3d(tlo,thi,lo,hi,cons,ng,dx,phlo,phhi) bind(c,name='init_data_3d')

    use variables_module, only : irho, imx,imy,imz,iene,iry1,ncons, iH2, iO2, iN2
    use chemistry_module, only : nspecies, patm

    implicit none

    integer,          intent(in   ) :: tlo(3),thi(3),lo(3),hi(3),ng
    double precision, intent(in   ) :: dx(3),phlo(3),phhi(3)
    double precision, intent(inout) :: cons(-ng+lo(1):hi(1)+ng,-ng+lo(2):hi(2)+ng,-ng+lo(3):hi(3)+ng,ncons)

    double precision, parameter :: Pi = 3.141592653589793238462643383279502884197d0
    double precision, parameter :: rfire = 0.01d0
    double precision, parameter :: Lperiodic = 0.1d0

    integer          :: i,j,k,n
    double precision :: x, y, z, r

    double precision Xt(nspecies), Yt(nspecies)
    double precision rhot,u1t,u2t,u3t,Tt,et,Pt
    integer :: iwrk
    double precision :: rwrk, expfac, kx, ky, kz

    kx = 2.d0*Pi/(phhi(1) - phlo(1))
    ky = 2.d0*Pi/(phhi(2) - phlo(2))
    kz = 2.d0*Pi/(phhi(3) - phlo(3))

    do k=tlo(3),thi(3)
       z = mod(dx(3)*(k+0.5d0), Lperiodic) - 0.5d0*Lperiodic
       do j=tlo(2),thi(2)
          y = mod(dx(2)*(j+0.5d0), Lperiodic) - 0.5d0*Lperiodic
          do i=tlo(1),thi(1)
             x = mod(dx(1)*(i+0.5d0), Lperiodic) - 0.5d0*Lperiodic

             r = sqrt(x**2+y**2+z**2)

             Pt = Patm
             Tt = 300.0d0

             Xt = 0.0d0
             Xt(iH2) = 0.10d0
             Xt(iO2) = 0.25d0

             expfac = exp(-(r / rfire)**2)
             Pt      = Pt      + 0.1d0*patm * expfac
             Tt      = Tt      + 1100.0d0   * expfac
             Xt(iH2) = Xt(iH2) + 0.025d0    * expfac
             Xt(iO2) = Xt(iO2) - 0.050d0    * expfac

             Xt(iN2) = 1.0d0 - Xt(iH2) - Xt(iO2)

             u1t =  sin(kx*x)*cos(ky*y)*cos(kz*z) * 300.d0
             u2t = -cos(kx*x)*sin(ky*y)*cos(kz*z) * 300.d0
             u3t = 0.d0

             CALL CKXTY (Xt, IWRK, RWRK, Yt)
             CALL CKRHOY(Pt,Tt,Yt,IWRK,RWRK,rhot)
             call CKUBMS(Tt,Yt,IWRK,RWRK,et)

             cons(i,j,k,irho) = rhot
             cons(i,j,k,imx)  = rhot*u1t
             cons(i,j,k,imy)  = rhot*u2t
             cons(i,j,k,imz)  = rhot*u3t
             cons(i,j,k,iene) = rhot*(et + 0.5d0*(u1t**2 + u2t**2 + u3t**2))

             do n=1,nspecies
                cons(i,j,k,iry1-1+n) = Yt(n)*rhot
             end do

          enddo
       enddo
    enddo

  end subroutine init_data_3d

end module init_data_module

max_stars_repo_name: ylunalin/amrex

max_issues_repo_name: ylunalin/amrex

max_forks_repo_name: ylunalin/amrex


text:
Fernanda Guimarães Rosa

a.k.a Fer, Fefa, Fezoca

journalist & web worker

Brazilian, from Campinas, São Paulo

living in Davis, California, USA

Married with Children

one husband, one son

Two Cats cats: Senhor Misty Gray & Roux

Love to: talk, dance, walk, swim, cook, write

Good things: movies, music, books, Photography pictures, magazines, Thrift Stores thrift stores

Food blog in Portuguese

http://www.chucrutecomsalsicha.com


text:
  !! Local Data
  INTEGER :: i, j
  INTEGER :: total_values

  !! There can't be more than one entry per row
  CALL ConstructTripletList(triplet_list, this%local_rows)

  total_values = 0
  !! Find local identity values
  row_iter: DO j = 1, this%local_rows
     column_iter: DO i = 1, this%local_columns
        IF (j + this%start_row - 1 .EQ. i + this%start_column - 1 .AND. &
             & j+this%start_row-1 .LE. this%actual_matrix_dimension) THEN
           total_values = total_values + 1
           triplet_list%data(total_values)%index_column = i
           triplet_list%data(total_values)%index_row = j
           triplet_list%data(total_values)%point_value = 1.0
        END IF
     END DO column_iter
  END DO row_iter

  !! Finish constructing
  CALL ConstructTripletList(unsorted_triplet_list, total_values)
  unsorted_triplet_list%data = triplet_list%data(:total_values)
  CALL SortTripletList(unsorted_triplet_list,this%local_columns,&
       & this%local_rows, sorted_triplet_list)
  CALL ConstructMatrixFromTripletList(local_matrix, sorted_triplet_list, &
       & this%local_rows, this%local_columns)

  CALL SplitMatrixToLocalBlocks(this, local_matrix)

  CALL DestructMatrix(local_matrix)
  CALL DestructTripletList(triplet_list)
  CALL DestructTripletList(unsorted_triplet_list)
  CALL DestructTripletList(sorted_triplet_list)

max_stars_repo_name: shoshijak/NTPoly

max_issues_repo_name: omni-compiler/NTPoly

max_forks_repo_name: omni-compiler/NTPoly


text:
      recursive subroutine dqawf(f,a,omega,integr,epsabs,result,
     *   abserr,neval,ier,limlst,lst,leniw,maxp1,lenw,iwork,work)
c***begin prologue  dqawf
c***date written   800101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a3a1
c***keywords  automatic integrator, special-purpose,fourier
c             integral, integration between zeros with dqawoe,
c             convergence acceleration with dqelg
c***author  piessens,robert ,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math & progr. div. - k.u.leuven
c***purpose  the routine calculates an approximation result to a given
c            fourier integral i=integral of f(x)*w(x) over (a,infinity)
c            where w(x) = cos(omega*x) or w(x) = sin(omega*x).
c            hopefully satisfying following claim for accuracy
c            abs(i-result).le.epsabs.
c***description
c
c        computation of fourier integrals
c        standard fortran subroutine
c        double precision version
c
c
c        parameters
c         on entry
c            f      - double precision
c                     function subprogram defining the integrand
c                     function f(x). the actual name for f needs to be
c                     declared e x t e r n a l in the driver program.
c
c            a      - double precision
c                     lower limit of integration
c
c            omega  - double precision
c                     parameter in the integrand weight function
c
c            integr - integer
c                     indicates which of the weight functions is used
c                     integr = 1      w(x) = cos(omega*x)
c                     integr = 2      w(x) = sin(omega*x)
c                     if integr.ne.1.and.integr.ne.2, the routine
c                     will end with ier = 6.
c
c            epsabs - double precision
c                     absolute accuracy requested, epsabs.gt.0.
c                     if epsabs.le.0, the routine will end with ier = 6.
c
c         on return
c            result - double precision
c                     approximation to the integral
c
c            abserr - double precision
c                     estimate of the modulus of the absolute error,
c                     which should equal or exceed abs(i-result)
c
c            neval  - integer
c                     number of integrand evaluations
c
c            ier    - integer
c                     ier = 0 normal and reliable termination of the
c                             routine. it is assumed that the requested
c                             accuracy has been achieved.
c                     ier.gt.0 abnormal termination of the routine.
c                             the estimates for integral and error are
c                             less reliable. it is assumed that the
c                             requested accuracy has not been achieved.
c            error messages
c                    if omega.ne.0
c                     ier = 1 maximum number of cycles allowed
c                             has been achieved, i.e. of subintervals
c                             (a+(k-1)c,a+kc) where
c                             c = (2*int(abs(omega))+1)*pi/abs(omega),
c                             for k = 1, 2, ..., lst.
c                             one can allow more cycles by increasing
c                             the value of limlst (and taking the
c                             according dimension adjustments into
c                             account). examine the array iwork which
c                             contains the error flags on the cycles, in
c                             order to look for eventual local
c                             integration difficulties.
c                             if the position of a local difficulty
c                             can be determined (e.g. singularity,
c                             discontinuity within the interval) one
c                             will probably gain from splitting up the
c                             interval at this point and calling
c                             appropriate integrators on the subranges.
c                         = 4 the extrapolation table constructed for
c                             convergence accelaration of the series
c                             formed by the integral contributions over
c                             the cycles, does not converge to within
c                             the requested accuracy.
c                             as in the case of ier = 1, it is advised
c                             to examine the array iwork which contains
c                             the error flags on the cycles.
c                         = 6 the input is invalid because
c                             (integr.ne.1 and integr.ne.2) or
c                              epsabs.le.0 or limlst.lt.1 or
c                              leniw.lt.(limlst+2) or maxp1.lt.1 or
c                              lenw.lt.(leniw*2+maxp1*25).
c                              result, abserr, neval, lst are set to
c                              zero.
c                         = 7 bad integrand behaviour occurs within
c                             one or more of the cycles. location and
c                             type of the difficulty involved can be
c                             determined from the first lst elements of
c                             vector iwork.  here lst is the number of
c                             cycles actually needed (see below).
c                             iwork(k) = 1 the maximum number of
c                                          subdivisions (=(leniw-limlst)
c                                          /2) has been achieved on the
c                                          k th cycle.
c                                      = 2 occurrence of roundoff error
c                                          is detected and prevents the
c                                          tolerance imposed on the k th
c                                          cycle, from being achieved
c                                          on this cycle.
c                                      = 3 extremely bad integrand
c                                          behaviour occurs at some
c                                          points of the k th cycle.
c                                      = 4 the integration procedure
c                                          over the k th cycle does
c                                          not converge (to within the
c                                          required accuracy) due to
c                                          roundoff in the extrapolation
c                                          procedure invoked on this
c                                          cycle. it is assumed that the
c                                          result on this interval is
c                                          the best which can be
c                                          obtained.
c                                      = 5 the integral over the k th
c                                          cycle is probably divergent
c                                          or slowly convergent. it must
c                                          be noted that divergence can
c                                          occur with any other value of
c                                          iwork(k).
c                    if omega = 0 and integr = 1,
c                    the integral is calculated by means of dqagie,
c                    and ier = iwork(1) (with meaning as described
c                    for iwork(k),k = 1).
c
c         dimensioning parameters
c            limlst - integer
c                     limlst gives an upper bound on the number of
c                     cycles, limlst.ge.3.
c                     if limlst.lt.3, the routine will end with ier = 6.
c
c            lst    - integer
c                     on return, lst indicates the number of cycles
c                     actually needed for the integration.
c                     if omega = 0, then lst is set to 1.
c
c            leniw  - integer
c                     dimensioning parameter for iwork. on entry,
c                     (leniw-limlst)/2 equals the maximum number of
c                     subintervals allowed in the partition of each
c                     cycle, leniw.ge.(limlst+2).
c                     if leniw.lt.(limlst+2), the routine will end with
c                     ier = 6.
c
c            maxp1  - integer
c                     maxp1 gives an upper bound on the number of
c                     chebyshev moments which can be stored, i.e. for
c                     the intervals of lengths abs(b-a)*2**(-l),
c                     l = 0,1, ..., maxp1-2, maxp1.ge.1.
c                     if maxp1.lt.1, the routine will end with ier = 6.
c            lenw   - integer
c                     dimensioning parameter for work
c                     lenw must be at least leniw*2+maxp1*25.
c                     if lenw.lt.(leniw*2+maxp1*25), the routine will
c                     end with ier = 6.
c
c         work arrays
c            iwork  - integer
c                     vector of dimension at least leniw
c                     on return, iwork(k) for k = 1, 2, ..., lst
c                     contain the error flags on the cycles.
c
c            work   - double precision
c                     vector of dimension at least
c                     on return,
c                     work(1), ..., work(lst) contain the integral
c                      approximations over the cycles,
c                     work(limlst+1), ..., work(limlst+lst) contain
c                      the error extimates over the cycles.
c                     further elements of work have no specific
c                     meaning for the user.
c
c***references  (none)
c***routines called  dqawfe,xerror
c***end prologue  dqawf
c
       double precision a,abserr,epsabs,f,omega,result,work
       integer ier,integr,iwork,last,leniw,lenw,limit,limlst,ll2,lvl,
     *  lst,l1,l2,l3,l4,l5,l6,maxp1,neval
c
       dimension iwork(leniw),work(lenw)
c
       external f
c
c         check validity of limlst, leniw, maxp1 and lenw.
c
c***first executable statement  dqawf
      ier = 6
      neval = 0
      last = 0
      result = 0.0d+00
      abserr = 0.0d+00
      if(limlst.lt.3.or.leniw.lt.(limlst+2).or.maxp1.lt.1.or.lenw.lt.
     *   (leniw*2+maxp1*25)) go to 10
c
c         prepare call for dqawfe
c
      limit = (leniw-limlst)/2
      l1 = limlst+1
      l2 = limlst+l1
      l3 = limit+l2
      l4 = limit+l3
      l5 = limit+l4
      l6 = limit+l5
      ll2 = limit+l1
      call dqawfe(f,a,omega,integr,epsabs,limlst,limit,maxp1,result,
     *  abserr,neval,ier,work(1),work(l1),iwork(1),lst,work(l2),
     *  work(l3),work(l4),work(l5),iwork(l1),iwork(ll2),work(l6))
c
c         call error handler if necessary
c
      lvl = 0
10    if(ier.eq.6) lvl = 1
      if(ier.ne.0) call xerror('abnormal return from dqawf',26,ier,lvl)
      return
      end

max_stars_repo_name: jcampbell/scipy

max_issues_repo_name: jcampbell/scipy

max_forks_repo_name: jcampbell/scipy


text:
      subroutine rtpart(iplane,rtmasy)
c
c     + + + PURPOSE + + +
c     Partitions total root mass into layers.
c     Called from PTGRA & PTGRP.
c     Author(s):
c     Reference in User Guide:
c
c     Changes:
c
c     Version: This module recoded from WEPP version 91.50.
c     Date recoded: 12/05/91.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pmxpln.inc'
      include 'pmxnsl.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer iplane
      real rtmasy
c
c     + + + ARGUMENT DEFINITIONS + + +
c     iplane - current OFE
c     rtmasy - root mass yesterday
c
c     + + + COMMON BLOCKS + + +
      include 'ccrpout.inc'
c        read: rtd(mxplan),rtmass(mxplan)
c      modify: rtm15(mxplan),rtm30(mxplan),rtm60(mxplan)
c
c     + + + LOCAL VARIABLES + + +
      real delt
c
c     + + + LOCAL DEFINITIONS + + +
c     delt   - increase in rootmass from yesterday to today
c
c     + + + END SPECIFICATIONS + + +
c
c
      if (rtd(iplane).le.0.15) then
        rtm15(iplane) = rtmass(iplane)
        rtm30(iplane) = 0.0
        rtm60(iplane) = 0.0
c
      else if (rtd(iplane).gt.0.15.and.rtd(iplane).le.0.30) then
        delt = rtmass(iplane) - rtmasy
        rtm15(iplane) = rtm15(iplane) + delt * 0.6
        rtm30(iplane) = rtm30(iplane) + delt * 0.4
        rtm60(iplane) = 0.0
c
      else if (rtd(iplane).gt.0.30.and.rtd(iplane).le.0.60) then
        delt = rtmass(iplane) - rtmasy
        rtm15(iplane) = rtm15(iplane) + delt * 0.45
        rtm30(iplane) = rtm30(iplane) + delt * 0.30
        rtm60(iplane) = rtm60(iplane) + delt * 0.25
c
      else
        delt = rtmass(iplane) - rtmasy
        rtm15(iplane) = rtm15(iplane) + delt * 0.42
        rtm30(iplane) = rtm30(iplane) + delt * 0.28
        rtm60(iplane) = rtm60(iplane) + delt * 0.20
c
c     rtm60(iplane)=rtm60(iplane)+delt*0.20
c     -- XXX -- Shouldn't the above fractions add to '1.00' ? -- CRM -- 12/5/91
c     -  NOTE - Changed 0.20 to 0.30 for rtm60    dcf  5/25/93
c     Changed back - after discussions with Arnold and Alberts - the extra
c     10% is assumed to go to the root mass below 60 centimeters  dcf  11/18/93
c
      end if
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
      subroutine mxm(a,n1,b,n2,c,n3)

#if defined(XSMM_DISPATCH)
      USE :: LIBXSMM
#endif

#define LIBXSMM_DMM1(N, a, b, c)  LIBXSMM_DMM1_str(N, a, b, c)
#define LIBXSMM_DMM1_str(N, a, b, c)  libxsmm_dmm_##N##x##N##_##N##_##N(a, b, c)
#define LIBXSMM_DMM2(N, a, b, c)   LIBXSMM_DMM2_str(N, a, b, c)
#define LIBXSMM_DMM2_str(N, a, b, c)  libxsmm_dmm_##N##_##N##x##N##_##N(a, b, c)
#define LIBXSMM_DMM3(N, a, b, c)  LIBXSMM_DMM3_str(N, a, b, c)
#define LIBXSMM_DMM3_str(N, a, b, c)  libxsmm_dmm_##N##_##N##_##N(a, b, c)

c
c     Compute matrix-matrix product C = A*B
c     for contiguously packed matrices A,B, and C.
c
#if defined (MKL)
#     include  "mkl_direct_call.fi"
#endif
      real a(n1,n2),b(n2,n3),c(n1,n3)
      real alpha, beta
c
      include 'SIZE'
      include 'TOTAL'
c
      integer aligned
      integer K10_mxm
      integer init, prevn2

#if defined(XSMM_DISPATCH)
      TYPE(LIBXSMM_DMMFUNCTION) :: xmm1,xmm2,xmm3
#endif

      data init /0/, prevn2 /0/
      save init, prevn2
#if defined(XSMM_DISPATCH)
      save xsmm1, xsmm2, xsmm3
#endif

c     write(*,*) "in", init, prevn2, LOC(xsmm1), LOC(xsmm2), LOC(xsmm3)

#if defined (MKL)
      alpha = 1.0
      beta = 0.0
      call dgemm('N','N',n1,n3,n2,alpha,A,n1,B,n2,beta,C,n1)
#elif defined (BLAS_MXM)
      alpha = 1.0
      beta = 0.0
      call dgemm('N','N',n1,n3,n2,alpha,a,n1,b,n2,beta,c,n1)
#elif defined (BG)
      call bg_aligned3(a,b,c,aligned)
      if (n2.eq.2) then
         call mxm44_2(a,n1,b,n2,c,n3)
      else if ((aligned.eq.1) .and.
     $         (n1.ge.8) .and. (n2.ge.8) .and. (n3.ge.8) .and.
     $         (modulo(n1,2).eq.0) .and. (modulo(n2,2).eq.0) ) then
         if (modulo(n3,4).eq.0) then
            call bg_mxm44(a,n1,b,n2,c,n3)
         else
            call bg_mxm44_uneven(a,n1,b,n2,c,n3)
         endif
      else if((aligned.eq.1) .and.
     $        (modulo(n1,6).eq.0) .and. (modulo(n3,6).eq.0) .and.
     $        (n2.ge.4) .and. (modulo(n2,2).eq.0) ) then
         call bg_mxm3(a,n1,b,n2,c,n3)
      else
         call mxm44_0(a,n1,b,n2,c,n3)
      endif
#elif defined (K10_MXM)
      ! fow now only supported for lx1=8
      ! tuned for AMD K10
      ierr = K10_mxm(a,n1,b,n2,c,n3)
      if (ierr.gt.0) call mxmf2(a,n1,b,n2,c,n3)
#elif defined (XSMM_DISPATCH)
      if (init == 0) then
        CALL libxsmm_init()
        init = 1
        write(*,*) "initializing libxsmm"
      end if

      if (prevn2 /= n2) then
        prevn2 = n2

        CALL libxsmm_dispatch(xmm1, n2, n2, n2*n2, alpha=1D0, beta=0D0)
c       write(*,*) "initialized xmm1"
        IF (.NOT. libxsmm_available(xmm1)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF

        CALL libxsmm_dispatch(xmm2, n2, n2, n2, alpha=1D0, beta=0D0)
c       write(*,*) "initialized xmm2"
        IF (.NOT. libxsmm_available(xmm2)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF

        CALL libxsmm_dispatch(xmm3, n2*n2, n2, n2, alpha=1D0, beta=0D0)
c       write(*,*) "initialized xmm3"
        IF (.NOT. libxsmm_available(xmm3)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF
      end if

      if (n1 .eq. n2*n2) then
c       write(*,*) "call to xmm3", n1, n2, n3
        IF (.NOT. libxsmm_available(xmm3)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF
        call libxsmm_call(xmm3, C_LOC(a), C_LOC(b), C_LOC(c))
c       call libxsmm_dmm_256_16_16(a, b, c)
      else if (n3 .eq. n2*n2) then
c       write(*,*) "call to xmm1", n1, n2, n3
        IF (.NOT. libxsmm_available(xmm1)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF
        call libxsmm_call(xmm1, C_LOC(a), C_LOC(b), C_LOC(c))
c       call libxsmm_dmm_16_256_16(a, b, c)
      else
c       write(*,*) "call to xmm2", n1, n2, n3
        IF (.NOT. libxsmm_available(xmm2)) THEN
          write(*,*) "  ** Error: unable to dispatch libxsmm call"
          STOP
        END IF
        call libxsmm_call(xmm2, C_LOC(a), C_LOC(b), C_LOC(c))
c       call libxsmm_dmm_16_16_16(a, b, c)
      end if
#elif defined (XSMM_FIXED)
      if (n2 .eq. NPOLY) then
        if (n1 .eq. n2*n2) then
          call LIBXSMM_DMM1(NPOLY, a, b, c)
        else if (n3 .eq. n2*n2) then
          call LIBXSMM_DMM2(NPOLY, a, b, c)
        else
          call LIBXSMM_DMM3(NPOLY, a, b, c)
        end if
      else
        write(*,*) "Invalid matrix size"
        stop
      end if
#elif defined (XSMM)
      alpha = 1.0
      beta = 0.0
      CALL libxsmm_dgemm('N','N',n1,n3,n2,alpha,A,n1,B,n2,beta,C,n1)
#elif defined (MXMBASIC)
      do j=1,n3
        do i=1,n1
          c(i,j) = 0.0
          do k=1,n2
            c(i,j) = c(i,j) + a(i,k)*b(k,j)
          enddo
        enddo
      enddo
#else
      call mxmf2(a,n1,b,n2,c,n3)
#endif

c     write(*,*) "out", init, prevn2, xsmm1, xsmm2, xsmm3

      return
      end


text:
	PROGRAM ANTENNAZ
C
C	CALCULATES RESISTANCE, CAPACITANCE AND SHOT NOISE ON WIND ANTENNAS
C
C	X ANTENNA
	PHYSLEN = 5000.
	EFFLEN = .5*PHYSLEN
	ICS = 1
	TPHI = 1.3		! eV
	TPHI = .7		! eV
	DIA = .038
	DENSE = 7.
	TE = 12.
	RATIO = 1.1
	DENSE = 2.
	DO N = 1,30
	  DENSE = DENSE*RATIO
	  FLUX = DENSE*SQRT(TE)
	  CALL ANTENNA(EFFLEN,PHYSLEN,DIA,ICS,DENSE,TE,TPHI,ANTCURR,
     1		RES,CAP)
	  PRINT*,ANTCURR,RES,CAP
	  WRITE(66,*) DENSE,FLUX,RES
	ENDDO
	STOP
	END


text:
County Fair Mall is a Malls mall in wiki:Woodland. Yolobus lines http://www.yolobus.com/42.html 42 and http://www.yolobus.com/242.html 242 will take you right there.

A few of the major businesses:

  Avon
  Burlington Coat Factory
  Sees Candy
  Chase bank formerly WaMu (ATM is drive up and walk up)
  CitiWear
  County Fair Mall Theater Cinemark Movies 5
  Claires
  J & G Gift & Watch Repair
  JCPenney
  Matthews Mattress
  Pac Sun
  Payless ShoeSource
  RadioShack
  Sassy

There was an Arcade called Timeout. It closed sometime in 2010. It replaced Fun and Games, which replaced Tilt.

A Cindys Cinnamon Rolls managed to stay open from the opening of the mall to very recently. It served fare much like Cinnabon.

A new General Manager was hired in May 2008 and the mall underwent many upgrades and opened some new stores. The Daily Democrat reported that tenancy increased to 85 percent, up from 55 percent in 2005.  In December 2005, the mall was purchased by Encinobased County Fair Fashion Mall LLC for $15.1 million. It was previously owned by Santa Anabased Triple Net Properties LLC. (Source: http://www.sacbee.com/content/community_news/yolo/story/13981716p14815402c.html Sac Bee 12/15/2005).

At many times this mall has been referred to as the Dirt Mall amongst local youth because of plants growing inside, with the dirt sometimes tracked around by kids.

A brand new Food Court was built, but a Sees Candy and a Debs Deli look to be the only businesses selling food these days.

The Target at this shopping center closed and moved to a new store at the Woodland Gateway retail center in July 2008. The old store was small and featured Target decor from the 80s. Burlington Coat Factory has taken the place of the 88,000 square foot former Target as of March 2009. Mervyns filed for bankruptcy and closed in late 2008. Washington Mutual was sold to JPMorgan Chase and was rebranded to Chase in March of 2009. Gottschalks filed for bankruptcy and was sold to liquidators and the store eventually closed, essentially cutting foot traffic in the mall down to basically nothing.

As of Fall, 2007, there is an alternative high school smack dab in the middle of the mall, courtesy of Einstein Education Center. Its funded, in part, by a grant from the Bill and Melinda Gates Foundation.  E^2^C served 150 students the first year (20072008), with a full capacity of 400 students at the end of the second year.  There have been fundraisers and benefit events, such as Sacramento Kings player http://www.mikebibbyinwoodland.com/index2.php Mike Bibby coming to a meet and greet.

Also see Shopping Centers.

The mall also hosts a http://pi.cdfa.ca.gov/scripts/cfm/pdf_app_view.asp?id174 Farmers Market in their parking lot in the summer and fall  Tuesdays, 5PM7PM and Saturdays, 9AM  12PM.  Theres a http://insidewoodland.com/farmersmarket.cfm story and some photos at the http://www.insidewoodland.com InsideWoodland web site.

The Mall bus stop is also a transfer point where you can take a Yolobus to get to the new Woodland Target, Best Buy, and Costco.

20060327 14:34:00 nbsp Right now is a good time to shop here. There are a lot of sales going on  I found the shoes I wanted at a decent price at, of all places, Mervyns. One store is planning to close for the summer and wants to clear out the inventory.

If you require the whole real mall experience  you know: no parking, too many people, a real food court  then try Arden Fair. We found County Fair Mall to be quiet and almost restful. And the pizza is pretty decent. users/JeffreyNonken

20060719 15:16:34 nbsp Good things: They have an orange julius which I love, a $3 photobooth and the amazing shiatsu massage chairs.  users/CatStarr

20070322 23:47:02 nbsp This is by far the most depressing mall I have ever been in. The ceiling is too low, the vacancies are prolific, and the remaining stores are mostly that kind of startup looking to fill an unexploited niche, usually with names like Bobs Shoehorn Shack.  users/JoePomidor

20070520 10:54:55 nbsp pretty lame excuse for a mall..i work here, and its depressing..seriously, what kind of mall closes at 7pm on a saturday night?! Users/DonaldJaye

20070527 14:27:38 nbsp The Washington Mutual has a drivethru ATM. In fact, that is the only ATM they have. You have to walk up to the Drive Up ATM if you have no car. Users/Jedron
  Just the same way you walk up to a normal ATM if you dont have a car?? You just cant please stupid people.  Just so were clear, i just called you stupid, stupid. Users/piratechef

20070705 17:32:39 nbsp added some stores, info about Target Users/PatrickSing

20071002 11:51:52 nbsp If you are looking for clothes at nonDavis prices, I recommend a trip to Woodland just so you can visit this mall.  I usually go into Citiwear, Sassy, and Target.  Oh yeah, and I got my bed here too (at Matthews Mattress) for a good price... I guess they have discontinued models, which translates to good savings. Users/CalamityJanie

20090504 16:17:56 nbsp a good place to shop for clothes if you cant go far. Users/thtly

20091130 11:57:19 nbsp I wrote a little script for a short movie called Ghost Mall. This place was the inspiration. When I was a kid I once saw a guy with a miniature zoos worth of exotic snakes showing them off to all the little chillins. They also had a pizza version of the pepsi challenge featuring all the Woodland pizza places. Now its a nice quiet place to study and feel lonely. Users/GreatRyan


text:
 \ \033[<fgColor>;<bgColor>m see bgColor
 variable escLeading $5b1b escLeading ! \ $5b1b escLeading ! \ esc [
 variable fgColor 0 , $3733 fgColor ! \ fgWhite default
 variable bgColor 0 , $3034 fgColor ! \ bgBlack default
 : fgColor: create bl word cell+ @ , does> r> r> drop @ fgColor ! escLeading print fgColor print [char] m emit ;
 : bgColor: create bl word cell+ @ , does> r> r> drop @ bgColor ! escLeading print bgColor print [char] m emit ;
 fgColor: fgBlack 30
 fgColor: fgRed 31
 fgColor: fgGreen 32
 fgColor: fgYellow 33
 fgColor: fgBlue 34
 fgColor: fgMagenta 35
 fgColor: fgCyan 36
 fgColor: fgWhite 37
 fgColor: fgBrightBlack 30;1
 fgColor: fgBrightRed 31;1
 fgColor: fgBrightGreen 32;1
 fgColor: fgBrightYellow 33;1
 fgColor: fgBrightBlue 34;1
 fgColor: fgBrightMagenta 35;1
 fgColor: fgBrightCyan 36;1
 fgColor: fgBrightWhite 37;1
 bgColor: bgBlack 40
 bgColor: bgRed 41
 bgColor: bgGreen 42
 bgColor: bgYellow 43
 bgColor: bgBlue 44
 bgColor: bgMagenta 45
 bgColor: bgCyan 46
 bgColor: bgWhite 47
 bgColor: bgBrightBlack 100
 bgColor: bgBrightRed 101
 bgColor: bgBrightGreen 102
 bgColor: bgBrightYellow 103
 bgColor: bgBrightBlue 104
 bgColor: bgBrightMagenta 105
 bgColor: bgBrightCyan 106
 bgColor: bgBrightWhite 107