text:
!   netJina_def.f
!
!   defines data layout for the JINA reaclib database
!   Edward Brown, Michigan State University
!
!   Requires installation of MESA (mesa.sourceforge.net) v 6022 for utils and
!   const modules
!

module netJina_def
    use const_def, only: sp,dp

    ! data storage parameter for nuclib
    ! each isotope has a name and "provenance" -- a reference to where the
    ! mass, spin, partition function are described.
    integer, parameter :: iso_name_length = 5
    integer, parameter :: provenance_length = 6
    integer, parameter :: reaction_reference_length = 4
    ! maximum number of nuclides in nucchem database
    integer, parameter :: max_nnuclib=10000
    ! no. entries in partition fcn table
    integer, parameter :: npfcn = 24

    ! maximum number of individual rates in reaclib.
    integer, parameter :: max_nreaclib=120000
    ! data storage parameters for reaclib
    integer, parameter :: max_species_per_reaction=6
    integer, parameter :: ncoefficients=7
    integer, parameter :: nchapters=11
    integer, parameter :: max_terms_per_rate = 12
    ! we keep track of reactions using a unique "handle" generated from the
    ! isotope names
    integer, parameter :: max_id_length = 42
    ! max. Z in the nucchem database
    integer, parameter :: max_element_z = 112
    ! used to output a formatted reaction string
    integer, parameter :: length_reaction_string = 45

    ! for starlib
    ! maximum number of individual rates in starlib
    integer, parameter :: max_nstarlib = 80000
    ! data storage parameters for starlib
    integer, parameter :: number_starlib_temps = 60
    integer, parameter :: starlib_max_species_per_reaction = 6

    ! table of elements; note that hydrogen is referred to as 'p' in the table
    ! and neutrons are referred to as 'n'
    character(len=iso_name_length), dimension(0:max_element_z) ::  &
    & element_name = [character(len=iso_name_length) :: &
    & 'neut','h','he','li','be','b','c','n','o','f','ne', 'na','mg','al','si', &
    & 'p', 's', 'cl', 'ar', 'k', 'ca', 'sc', 'ti', 'v','cr','mn', 'fe', 'co', &
    & 'ni','cu','zn','ga','ge','as','se','br','kr','rb','sr','y','zr', &
    & 'nb','mo','tc','ru','rh','pd','ag','cd','in','sn','sb','te','i','xe', &
    & 'cs','ba','la','ce','pr','nd','pm','sm','eu','gd','tb','dy','ho','er', &
    & 'tm','yb', 'lu','hf','ta','w','re','os','ir','pt','au','hg', &
    & 'tl','pb','bi','po','at','rn','fr','ra','ac','th','pa','u','np','pu', &
    & 'am','cm','bk','cf','es','fm','md','no','lr','rf','db','sg','bh','hs', &
    & 'mt','ds','rg','cn' ]

    ! synonyms of hydrogen isotopes, index by mass number
    character(len=iso_name_length), dimension(3) ::  &
    & h_isotopes = [character(len=iso_name_length) :: 'p','d','t']

    ! aluminum-26 isomers
    character(len=iso_name_length), dimension(2:3) ::  &
    & al_isomers = [character(len=iso_name_length) :: 'al-6','al*6']

    ! storage for nuclib database
    type nuclib_data
        integer :: Nnuclides
        ! for each nuclide, store its name, provenance, mass (real), charge and
        ! neutron numbers, ground-state spin, mass excess (MeV) and partition
        ! function table
        character(len=iso_name_length), dimension(:), allocatable :: name
        character(len=provenance_length),dimension(:), allocatable :: provenance
        real(dp), dimension(:), allocatable :: A
        integer, dimension(:), allocatable  :: Z
        integer, dimension(:), allocatable :: N
        real(dp), dimension(:), allocatable :: spin
        real(dp), dimension(:), allocatable :: mass_excess
        real(dp), dimension(:,:), allocatable :: pfcn
    end type nuclib_data

    ! temperatures, partition function evaluation
    real(dp), dimension(npfcn) :: pfcn_T9

    ! flags for reaction channels, bdat-style
    integer, parameter :: i_pg = 1
    integer, parameter :: i_an = i_pg+1
    integer, parameter :: i_ag = i_an+1
    integer, parameter :: i_ap = i_ag+1
    integer, parameter :: i_ng = i_ap+1
    integer, parameter :: i_np = i_ng+1
    integer, parameter :: i_gp = i_np+1
    integer, parameter :: i_na = i_gp+1
    integer, parameter :: i_ga = i_na+1
    integer, parameter :: i_pa = i_ga+1
    integer, parameter :: i_gn = i_pa+1
    integer, parameter :: i_pn = i_gn+1
    integer, parameter :: N_bdat_channels = i_pn

    integer,dimension(N_bdat_channels) ::  &
    & bdat_dZ = [1,2,2,1,0,-1,-1,-2,-2,-1,0,1],  &
    & bdat_dN = [0,1,2,2,1,1,0,-1,-2,-2,-1,-1]

    ! flags for reaction types: r_<number in>_<number out> = chapter id
    integer, parameter :: &
    &   r_one_one   = 1, &
    &   r_one_two   = 2, &
    &   r_one_three = 3, &
    &   r_two_one   = 4, &
    &   r_two_two   = 5, &
    &   r_two_three = 6, &
    &   r_two_four  = 7, &
    &   r_three_one = 8, &
    &   r_three_two = 9, &
    &   r_four_two  = 10,&
    &   r_one_four  = 11

    ! Nin(chapter), Nout(chapter) give number of nuclides on the entrance and
    ! exit channels for that chapter
    integer, dimension(nchapters) :: nJ_Nin = (/1,1,1,2,2,2,2,3,3,4,1/)
    integer, dimension(nchapters) :: nJ_Nout = (/1,2,3,1,2,3,4,1,2,2,4/)

    ! storage containers for rate data
    type rate_data
        integer :: Nentries
        integer,dimension(:),allocatable :: chapter
        character(len=iso_name_length),dimension(:,:),allocatable :: species
        character(len=reaction_reference_length),dimension(:),allocatable :: label
        character,dimension(:),allocatable :: reaction_flag
        character,dimension(:),allocatable :: reverse_flag
        real(dp),dimension(:),allocatable :: Qvalue
    end type rate_data

    ! storage container for reaclib file
    type, extends(rate_data) :: reaclib_data
        real(dp),dimension(:,:),allocatable :: coefficients
        integer,dimension(:),allocatable :: N_rate_terms
    end type reaclib_data

    ! storage container for starlib file
    type, extends(rate_data) :: starlib_data
        real(dp),dimension(:,:), allocatable :: T9
        real(dp),dimension(:,:), allocatable :: rate
        real(dp),dimension(:,:), allocatable :: uncertainty
    end type starlib_data

end module netJina_def

max_stars_repo_name: nworbde/netJina

max_issues_repo_name: nworbde/netJina

max_forks_repo_name: nworbde/netJina


text:
Julie Hooper was elected ASUCD Vice President on February 23, 2007.

Julie was an Business and Finance Commission ASUCD Business and Finance Commissioner, the Chief of Staff for the ASUCD Executive Office and a Cal Aggie Camp counselor. She is also the Philanthropy chair of Alpha Phi,the organizer of the King of Hearts philanthropy event as well as the ASUCD Blood Drive. She also volunteers for the ASUCD Lobby Corps and is an active LEAD LEADer.  She is an official candidate for Vice President in the Winter 2007 ASUCD Election under LEAD.
http://youtube.com/watch?vsunBnGCJ_c Kareem and Julies Winter 2007 video candidate statement.

20070225 02:23:06 nbsp Im honored that my photo of Julie is used on her profile! Im so happy and excited for the new presidential duo. Users/MissyWhitney

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 02.01.02.01/03                        *
C  *    TEST TITLE : Deletion of a non-existent structure  *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C PHIGS parameter types

      INTEGER NOSTID

      LOGICAL CSSEQ

      CALL INITGL ('02.01.02.01/03')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)

C creation of structures
      CALL POPST (100)
      CALL PEXST (105)
      CALL PEXST (110)
      CALL PCLST

      CALL SETMSG ('3', 'Deleting a non-existent structure ' //
     1             'should have no effect on the contents '  //
     2             'of the CSS.')

      NOSTID = 25
      CALL PDST (NOSTID)

      CALL IFPF (CSSEQ ('100,1, 105,0/ 100,2, 110,0/'))

      CALL ENDIT
      END


! diffusion_test.f90
! Generate test data for diffusion.f90
PROGRAM diffusion_test

  !------------------------------------------------------------------------------------------------!
  ! This software was written in 2016/17                                                           !
  ! by Michael P. Allen <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk>                        !
  ! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("the authors"),                             !
  ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"),     !
  ! published by Oxford University Press ("the publishers").                                       !
  !                                                                                !
  ! LICENCE                                                                                !
  ! Creative Commons CC0 Public Domain Dedication.                                                 !
  ! To the extent possible under law, the authors have dedicated all copyright and related         !
  ! and neighboring rights to this software to the PUBLIC domain worldwide.                        !
  ! This software is distributed without any warranty.                                             !
  ! You should have received a copy of the CC0 Public Domain Dedication along with this software.  !
  ! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.                               !
  !                                                                                !
  ! DISCLAIMER                                                                                !
  ! The authors and publishers make no warranties about the software, and disclaim liability       !
  ! for all uses of the software, to the fullest extent permitted by applicable law.               !
  ! The authors and publishers do not recommend use of this software for any purpose.              !
  ! It is made freely available, solely to clarify points made in the text. When using or citing   !
  ! the software, you should not imply endorsement by the authors or publishers.                   !
  !------------------------------------------------------------------------------------------------!

  ! Generates random configuration of atoms (positions, velocities)
  ! Cubic periodic boundary conditions
  ! Conducts Brownian dynamics without atomic interactions
  ! Resulting vacf should be exponential with supplied decay rate gamma
  ! Long-time mean square displacement should be diffusive with D = kT/gamma

  ! Reads several variables and options from standard input using a namelist nml
  ! Leave namelist empty to accept supplied defaults

  ! Input configuration, output configuration, all calculations, and results
  ! are given in simulation units defined by the model
  ! We assume mass m=1 throughout

  USE, INTRINSIC :: iso_fortran_env,  ONLY : input_unit, output_unit, error_unit, iostat_end, iostat_eor
  USE               config_io_module, ONLY : write_cnf_atoms
  USE               maths_module,     ONLY : random_normals

  IMPLICIT NONE

  ! Most important variables
  INTEGER :: n           ! Number of atoms
  REAL    :: box         ! Box length
  REAL    :: dt          ! Time step
  REAL    :: temperature ! Temperature (specified)
  REAL    :: gamma       ! Friction coefficient

  REAL, DIMENSION(:,:), ALLOCATABLE :: r    ! Positions (3,n)
  REAL, DIMENSION(:,:), ALLOCATABLE :: v    ! Velocities (3,n)
  REAL, DIMENSION(:,:), ALLOCATABLE :: zeta ! Random numbers (3,n)

  INTEGER :: blk, stp, nstep, nblock, ioerr

  CHARACTER(len=4), PARAMETER :: cnf_prefix = 'cnf.'
  CHARACTER(len=3)            :: sav_tag    = 'sav' ! May be overwritten with block number

  NAMELIST /nml/ n, nblock, nstep, dt, gamma, temperature, box

  WRITE ( unit=output_unit, fmt='(a)' ) 'diffusion_test'
  WRITE ( unit=output_unit, fmt='(a)' ) 'Brownian dynamics without interactions, constant-NVT ensemble'
  WRITE ( unit=output_unit, fmt='(a)' ) 'Particle mass m=1 throughout'

  ! Set sensible default run parameters for testing
  n           = 250
  nblock      = 999
  nstep       = 25
  dt          = 0.002
  gamma       = 1.0
  temperature = 1.0
  box         = 1.0

  ! Read run parameters from namelist
  ! Comment out, or replace, this section if you don't like namelists
  READ ( unit=input_unit, nml=nml, iostat=ioerr )
  IF ( ioerr /= 0 ) THEN
     WRITE ( unit=error_unit, fmt='(a,i15)') 'Error reading namelist nml from standard input', ioerr
     IF ( ioerr == iostat_eor ) WRITE ( unit=error_unit, fmt='(a)') 'End of record'
     IF ( ioerr == iostat_end ) WRITE ( unit=error_unit, fmt='(a)') 'End of file'
     STOP 'Error in diffusion_test'
  END IF

  ! Write out run parameters
  WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Number of atoms',           n
  WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Number of blocks',          nblock
  WRITE ( unit=output_unit, fmt='(a,t40,i15)'   ) 'Number of steps per block', nstep
  WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Time step',                 dt
  WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Friction coefficient',      gamma
  WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Temperature',               temperature
  WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Ideal diffusion coefft',    temperature / gamma

  ALLOCATE ( r(3,n), v(3,n), zeta(3,n) )

  ! Set random positions
  CALL RANDOM_SEED ()
  CALL RANDOM_NUMBER ( r )
  r = r - 0.5 ! Now in range (-1/2,1/2)
  r = r * box ! Now in range (-box/2,box/2)

  ! Set random velocities
  CALL random_normals ( 0.0, SQRT(temperature), v )

  sav_tag = '000' ! Initial configuration output file
  CALL write_cnf_atoms ( cnf_prefix//sav_tag, n, box, r, v ) ! Save configuration

  DO blk = 1, nblock ! Begin loop over blocks

     DO stp = 1, nstep ! Begin loop over steps

        CALL a_propagator ( dt/2.0 ) ! A drift half-step
        CALL o_propagator ( dt )     ! O random velocities and friction step
        CALL a_propagator ( dt/2.0 ) ! A drift half-step

     END DO ! End loop over steps

     r(:,:) = r(:,:) - box * ANINT ( r(:,:) / box )             ! Periodic boundaries
     IF ( nblock < 1000 ) WRITE(sav_tag,'(i3.3)') blk           ! Number configuration by block
     CALL write_cnf_atoms ( cnf_prefix//sav_tag, n, box, r, v ) ! Save configuration

  END DO ! End loop over blocks

  DEALLOCATE ( r, v, zeta )

  CALL exact ! Write out exact results for comparison

CONTAINS

  SUBROUTINE a_propagator ( t ) ! A propagator (drift)
    IMPLICIT NONE
    REAL, INTENT(in) :: t ! Time over which to propagate (typically dt/2)

    r(:,:) = r(:,:) + t * v(:,:) ! Positions

  END SUBROUTINE a_propagator

  SUBROUTINE o_propagator ( t ) ! O propagator (friction and random contributions)
    IMPLICIT NONE
    REAL, INTENT(in) :: t ! Time over which to propagate (typically dt)

    REAL            :: x, c
    REAL, PARAMETER :: c1 = 2.0, c2 = -2.0, c3 = 4.0/3.0, c4 = -2.0/3.0

    x = gamma * t
    IF ( x > 0.0001 ) THEN
       c = 1-EXP(-2.0*x)
    ELSE
       c = x * ( c1 + x * ( c2 + x * ( c3 + x * c4 )) ) ! Taylor expansion for low x
    END IF
    c = SQRT ( c )

    CALL random_normals ( 0.0, SQRT(temperature), zeta ) ! Random momenta
    v = EXP(-x) * v + c * zeta

  END SUBROUTINE o_propagator

  SUBROUTINE exact
    IMPLICIT NONE

    ! Writes out exact vacf, rvcf and msd for the Langevin equation
    ! for comparison with results of diffusion program

    INTEGER :: unit
    REAL    :: t, vacf, rvcf, msd

    WRITE ( unit=output_unit, fmt='(a)' ) 'Exact results output to diffusion_exact.out'

    OPEN ( newunit=unit, file='diffusion_exact.out', status='replace', iostat=ioerr )
    IF ( ioerr /= 0 ) THEN
       WRITE ( unit=error_unit, fmt='(a,i15)') 'Error opening file', ioerr
       STOP 'Error in diffusion'
    END IF

    DO blk = 0, nblock/2 ! Loop up to half the run length
       t    = (blk*nstep)*dt                                                    ! Time advances block by block
       vacf = 3.0*temperature * EXP(-gamma*t)                                   ! Velocity autocorrelation function
       rvcf = 3.0*temperature * ( 1.0 - EXP(-gamma*t) ) / gamma                 ! Velocity-displacement correlation
       msd  = 6.0*temperature * ( t - ( 1.0 - EXP(-gamma*t) ) / gamma ) / gamma ! Mean-square displacement
       WRITE ( unit=unit, fmt='(f15.6,3f15.8)' ) t, vacf, rvcf, msd
    END DO ! End loop up to half the run length

    CLOSE(unit=unit)

  END SUBROUTINE exact

END PROGRAM diffusion_test


text:
	SUBROUTINE GSLPAT  ( ilpat, iret )
C************************************************************************
C* GSLPAT								*
C* 									*
C* This subroutine sets the current line pattern to the given 		*
C* pattern.  The current line type number will be set to zero.  	*
C* The values specified in the line pattern designate the length 	*
C* of alternating on and off segments.  For example, a long-dash 	*
C* short-dash line pattern could be specified by LPAT values of:	*
C*									*
C*		10 5 5 5 0 0 0 0					*
C*									*
C* This line pattern will display a line for 10 units, space for 5 	*
C* units, a line for 5 units, and space for 5 units and repeat this	*
C* pattern.  Dots are specified using negative numbers.  The		*
C* absolute value of the number corresponds to the space in which	*
C* the dot will be centered.						*
C*									*
C* GSLPAT  ( ILPAT, IRET )						*
C*									*
C* Input parameters:							*
C*	ILPAT (8)	INTEGER		Line pattern values		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 8/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* S. Schotz/GSC	 2/90	Flush line before setting		*
C* S. Schotz/GSC	 8/90	Add lltyp to DSLPAT parameter list	*
C* M. Linda/GSC		 6/96	Removed FLLINE				*
C************************************************************************
	INCLUDE 	'ERROR.PRM'
	INCLUDE 	'DEVSET.CMN'
	INCLUDE 	'DEVCHR.CMN'
C*
	INTEGER 	 ilpat (*)
C------------------------------------------------------------------------
C* 	If device has not been set, return an error.
C
	IF  ( ddev .eq. ' ' )  THEN
	    iret = NDVICE
	  ELSE
C
C*	    Send values to device.
C
	    CALL DSLPAT  ( ilpat, lltyp, iret )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE QUADS (NEL,ITYPE,NINT,THIC,YM,PR,XX,S,IOUT)            QUA00001
C                                                                       QUA00002
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00003
C .                                                                   . QUA00004
C .   P R O G R A M                                                   . QUA00005
C .        TO CALCULATE ISOPARAMETRIC QUADRILATERAL ELEMENT STIFFNESS . QUA00006
C .        MATRIX FOR AXISYMMETRIC, PLANE STRESS, AND PLANE STRAIN    . QUA00007
C .        CONDITIONS                                                 . QUA00008
C .                                                                   . QUA00009
C .  - - INPUT VARIABLES - -                                          . QUA00010
C .        NEL       = NUMBER OF ELEMENT                              . QUA00011
C .        ITYPE     = ELEMENT TYPE                                   . QUA00012
C .                        EQ.0 = AXISYMMETRIC                        . QUA00013
C .                        EQ.1 = PLANE STRAIN                        . QUA00014
C .                        EQ.2 = PLANE STRESS                        . QUA00015
C .        NINT      = GAUSS NUMERICAL INTEGRATION ORDER              . QUA00016
C .        THIC      = THICKNESS OF ELEMENT                           . QUA00017
C .        YM        = YOUNG'S MODULUS                                . QUA00018
C .        PR        = POISSON'S RATIO                                . QUA00019
C .        XX(2,4)   = ELEMENT NODE COORDINATES                       . QUA00020
C .        S(8,8)    = STORAGE FOR STIFFNESS MATRIX                   . QUA00021
C .        IOUT      = UNIT NUMBER USED FOR OUTPUT                    . QUA00022
C .                                                                   . QUA00023
C .  - - OUTPUT - -                                                   . QUA00024
C .        S(8,8)    = CALCULATED STIFFNESS MATRIX                    . QUA00025
C .                                                                   . QUA00026
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00027
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               QUA00028
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00029
C .   THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CRAY     . QUA00030
C .   EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES,      . QUA00031
C .   ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR     . QUA00032
C .   SINGLE PRECISION ARITHMETIC.                                    . QUA00033
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00034
      DIMENSION D(4,4),B(4,8),XX(2,4),S(8,8),XG(4,4),WGT(4,4),DB(4)     QUA00035
C                                                                       QUA00036
C     MATRIX XG STORES GAUSS - LEGENDRE SAMPLING POINTS                 QUA00037
C                                                                       QUA00038
      DATA XG/   0.D0,   0.D0,   0.D0,   0.D0,   -.5773502691896D0,     QUA00039
     1 .5773502691896D0,   0.D0,   0.D0,   -.7745966692415D0,   0.D0,   QUA00040
     2 .7745966692415D0,   0.D0,   -.8611363115941D0,                   QUA00041
     3 -.3399810435849D0,   .3399810435849D0,   .8611363115941D0 /      QUA00042
C                                                                       QUA00043
C     MATRIX WGT STORES GAUSS - LEGENDRE WEIGHTING FACTORS              QUA00044
C                                                                       QUA00045
      DATA WGT /  2.D0,   0.D0,   0.D0,   0.D0,   1.D0,   1.D0,         QUA00046
     1 0.D0,   0.D0,   .5555555555556D0,   .8888888888889D0,            QUA00047
     2 .5555555555556D0,   0.D0,   .3478548451375D0,   .6521451548625D0,QUA00048
     3 .6521451548625D0,   .3478548451375D0 /                           QUA00049
C                                                                       QUA00050
C     O B T A I N  S T R E S S - S T R A I N  L A W                     QUA00051
C                                                                       QUA00052
      F=YM/(1.+PR)                                                      QUA00053
      G=F*PR/(1.-2.*PR)                                                 QUA00054
      H=F + G                                                           QUA00055
C                                                                       QUA00056
C     PLANE STRAIN ANALYSIS                                             QUA00057
C                                                                       QUA00058
      D(1,1)=H                                                          QUA00059
      D(1,2)=G                                                          QUA00060
      D(1,3)=0.                                                         QUA00061
      D(2,1)=G                                                          QUA00062
      D(2,2)=H                                                          QUA00063
      D(2,3)=0.                                                         QUA00064
      D(3,1)=0.                                                         QUA00065
      D(3,2)=0.                                                         QUA00066
      D(3,3)=F/2.                                                       QUA00067
      IF (ITYPE.EQ.1) THEN                                              QUA00068
      THIC=1.                                                           QUA00069
      GO TO 20                                                          QUA00070
      ENDIF                                                             QUA00071
C                                                                       QUA00072
C     AXISYMMETRIC ANALYSIS                                             QUA00073
C                                                                       QUA00074
      D(1,4)=G                                                          QUA00075
      D(2,4)=G                                                          QUA00076
      D(3,4)=0.                                                         QUA00077
      D(4,1)=G                                                          QUA00078
      D(4,2)=G                                                          QUA00079
      D(4,3)=0.                                                         QUA00080
      D(4,4)=H                                                          QUA00081
      IF (ITYPE.EQ.0) GO TO 20                                          QUA00082
C                                                                       QUA00083
C     FOR PLANE STRESS ANALYSIS CONDENSE STRESS-STRAIN MATRIX           QUA00084
C                                                                       QUA00085
      DO 10 I=1,3                                                       QUA00086
      A=D(I,4)/D(4,4)                                                   QUA00087
      DO 10 J=I,3                                                       QUA00088
      D(I,J)=D(I,J) - D(4,J)*A                                          QUA00089
   10 D(J,I)=D(I,J)                                                     QUA00090
C                                                                       QUA00091
C     C A L C U L A T E  E L E M E N T  S T I F F N E S S               QUA00092
C                                                                       QUA00093
   20 DO 30 I=1,8                                                       QUA00094
      DO 30 J=1,8                                                       QUA00095
   30 S(I,J)=0.                                                         QUA00096
      IST=3                                                             QUA00097
      IF (ITYPE.EQ.0) IST=4                                             QUA00098
      DO 80 LX=1,NINT                                                   QUA00099
      RI=XG(LX,NINT)                                                    QUA00100
      DO 80 LY=1,NINT                                                   QUA00101
      SI=XG(LY,NINT)                                                    QUA00102
C                                                                       QUA00103
C     EVALUATE DERIVATIVE OPERATOR B AND THE JACOBIAN DETERMINANT DET   QUA00104
C                                                                       QUA00105
      CALL STDM (XX,B,DET,RI,SI,XBAR,NEL,ITYPE,IOUT)                    QUA00106
C                                                                       QUA00107
C     ADD CONTRIBUTION TO ELEMENT STIFFNESS                             QUA00108
C                                                                       QUA00109
      IF (ITYPE.GT.0) XBAR=THIC                                         QUA00110
      WT=WGT(LX,NINT)*WGT(LY,NINT)*XBAR*DET                             QUA00111
      DO 70 J=1,8                                                       QUA00112
      DO 40 K=1,IST                                                     QUA00113
      DB(K)=0.0                                                         QUA00114
      DO 40 L=1,IST                                                     QUA00115
   40 DB(K)=DB(K) + D(K,L)*B(L,J)                                       QUA00116
      DO 60 I=J,8                                                       QUA00117
      STIFF=0.0                                                         QUA00118
      DO 50 L=1,IST                                                     QUA00119
   50 STIFF=STIFF + B(L,I)*DB(L)                                        QUA00120
   60 S(I,J)=S(I,J) + STIFF*WT                                          QUA00121
   70 CONTINUE                                                          QUA00122
   80 CONTINUE                                                          QUA00123
C                                                                       QUA00124
      DO 90 J=1,8                                                       QUA00125
      DO 90 I=J,8                                                       QUA00126
   90 S(J,I)=S(I,J)                                                     QUA00127
C                                                                       QUA00128
      RETURN                                                            QUA00129
C                                                                       QUA00130
      END                                                               QUA00131
      SUBROUTINE STDM (XX,B,DET,R,S,XBAR,NEL,ITYPE,IOUT)                QUA00132
C                                                                       QUA00133
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00134
C .                                                                   . QUA00135
C .   P R O G R A M                                                   . QUA00136
C .     TO EVALUATE THE STRAIN-DISPLACEMENT TRANSFORMATION MATRIX B   . QUA00137
C .     AT POINT (R,S) FOR A QUADRILATERAL ELEMENT                    . QUA00138
C .                                                                   . QUA00139
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . QUA00140
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               QUA00141
      DIMENSION XX(2,4),B(4,8),H(4),P(2,4),XJ(2,2),XJI(2,2)             QUA00142
C                                                                       QUA00143
      RP = 1.0 + R                                                      QUA00144
      SP = 1.0 + S                                                      QUA00145
      RM = 1.0 - R                                                      QUA00146
      SM = 1.0 - S                                                      QUA00147
C                                                                       QUA00148
C     INTERPOLATION FUNCTIONS                                           QUA00149
C                                                                       QUA00150
      H(1) = 0.25* RP* SP                                               QUA00151
      H(2) = 0.25* RM* SP                                               QUA00152
      H(3) = 0.25* RM* SM                                               QUA00153
      H(4) = 0.25* RP* SM                                               QUA00154
C                                                                       QUA00155
C     NATURAL COORDINATE DERIVATIVES OF THE INTERPOLATION FUNCTIONS     QUA00156
C                                                                       QUA00157
C        1. WITH RESPECT TO R                                           QUA00158
C                                                                       QUA00159
      P(1,1) = 0.25* SP                                                 QUA00160
      P(1,2) = - P(1,1)                                                 QUA00161
      P(1,3) = - 0.25* SM                                               QUA00162
      P(1,4) = - P(1,3)                                                 QUA00163
C                                                                       QUA00164
C        2. WITH RESPECT TO S                                           QUA00165
C                                                                       QUA00166
      P(2,1) = 0.25* RP                                                 QUA00167
      P(2,2) = 0.25* RM                                                 QUA00168
      P(2,3) = - P(2,2)                                                 QUA00169
      P(2,4) = - P(2,1)                                                 QUA00170
C                                                                       QUA00171
C     EVALUATE THE JACOBIAN MATRIX AT POINT (R,S)                       QUA00172
C                                                                       QUA00173
   10 DO 30 I=1,2                                                       QUA00174
      DO 30 J=1,2                                                       QUA00175
      DUM = 0.0                                                         QUA00176
      DO 20 K=1,4                                                       QUA00177
   20 DUM=DUM + P(I,K)*XX(J,K)                                          QUA00178
   30 XJ(I,J)=DUM                                                       QUA00179
C                                                                       QUA00180
C     COMPUTE THE DETERMINANT OF THE JACOBIAN MATRIX AT POINT (R,S)     QUA00181
C                                                                       QUA00182
      DET = XJ(1,1)* XJ(2,2) - XJ(2,1)* XJ(1,2)                         QUA00183
      IF (DET.GT.0.00000001) GO TO 40                                   QUA00184
      WRITE (IOUT,2000) NEL                                             QUA00185
      GO TO 800                                                         QUA00186
C                                                                       QUA00187
C     COMPUTE INVERSE OF THE JACOBIAN MATRIX                            QUA00188
C                                                                       QUA00189
   40 DUM=1./DET                                                        QUA00190
      XJI(1,1) = XJ(2,2)* DUM                                           QUA00191
      XJI(1,2) =-XJ(1,2)* DUM                                           QUA00192
      XJI(2,1) =-XJ(2,1)* DUM                                           QUA00193
      XJI(2,2) = XJ(1,1)* DUM                                           QUA00194
C                                                                       QUA00195
C     EVALUATE GLOBAL DERIVATIVE OPERATOR B                             QUA00196
C                                                                       QUA00197
      K2=0                                                              QUA00198
      DO 60 K=1,4                                                       QUA00199
      K2=K2 + 2                                                         QUA00200
      B(1,K2-1) = 0.                                                    QUA00201
      B(1,K2  ) = 0.                                                    QUA00202
      B(2,K2-1) = 0.                                                    QUA00203
      B(2,K2  ) = 0.                                                    QUA00204
      DO 50 I=1,2                                                       QUA00205
      B(1,K2-1) = B(1,K2-1) + XJI(1,I) * P(I,K)                         QUA00206
   50 B(2,K2  ) = B(2,K2  ) + XJI(2,I) * P(I,K)                         QUA00207
      B(3,K2  ) = B(1,K2-1)                                             QUA00208
   60 B(3,K2-1) = B(2,K2  )                                             QUA00209
C                                                                       QUA00210
C     IN CASE OF PLANE STRAIN OR PLANE STRESS ANALYSIS DO NOT INCLUDE   QUA00211
C     THE NORMAL STRAIN COMPONENT                                       QUA00212
C                                                                       QUA00213
      IF (ITYPE.GT.0) GO TO 900                                         QUA00214
C                                                                       QUA00215
C     COMPUTE THE RADIUS AT POINT (R,S)                                 QUA00216
C                                                                       QUA00217
      XBAR=0.0                                                          QUA00218
      DO 70 K=1,4                                                       QUA00219
   70 XBAR=XBAR + H(K)*XX(1,K)                                          QUA00220
C                                                                       QUA00221
C     EVALUATE THE HOOP STRAIN-DISPLACEMENT RELATION                    QUA00222
C                                                                       QUA00223
      IF (XBAR.GT.0.00000001) GO TO 90                                  QUA00224
C                                                                       QUA00225
C     FOR THE CASE OF ZERO RADIUS EQUATE RADIAL TO HOOP STRAIN          QUA00226
C                                                                       QUA00227
      DO 80 K=1,8                                                       QUA00228
   80 B(4,K)=B(1,K)                                                     QUA00229
      GO TO 900                                                         QUA00230
C                                                                       QUA00231
C     NON-ZERO RADIUS                                                   QUA00232
C                                                                       QUA00233
   90 DUM=1./XBAR                                                       QUA00234
      K2=0                                                              QUA00235
      DO 100 K=1,4                                                      QUA00236
      K2=K2 + 2                                                         QUA00237
      B(4,K2  ) = 0.                                                    QUA00238
  100 B(4,K2-1) = H(K)*DUM                                              QUA00239
      GO TO 900                                                         QUA00240
C                                                                       QUA00241
  800 STOP                                                              QUA00242
  900 RETURN                                                            QUA00243
C                                                                       QUA00244
 2000 FORMAT (//,' *** ERROR *** ',                                     QUA00245
     1    ' ZERO OR NEGATIVE JACOBIAN DETERMINANT FOR ELEMENT (',I8,')')QUA00246
C                                                                       QUA00247
      END                                                               QUA00248

ext: for


FUNCTION bessk1_s(x)
  use mo_kind
  use mo_nrutil, ONLY : assert,poly
  use mo_nr, ONLY : bessi1
  IMPLICIT NONE
  REAL(SP), INTENT(IN) :: x
  REAL(SP) :: bessk1_s
  REAL(sp) :: y
  REAL(sp), DIMENSION(7) :: p = (/1.0_sp,0.15443144_sp,&
       -0.67278579_sp,-0.18156897_sp,-0.1919402e-1_sp,&
       -0.110404e-2_sp,-0.4686e-4_sp/)
  REAL(sp), DIMENSION(7) :: q = (/1.25331414_sp,0.23498619_sp,&
       -0.3655620e-1_sp,0.1504268e-1_sp,-0.780353e-2_sp,&
       0.325614e-2_sp,-0.68245e-3_sp/)
  call assert(x > 0.0, 'bessk1_s arg')
  if (x <= 2.0) then
     y=x*x/4.0_sp
     bessk1_s=(log(x/2.0_sp)*bessi1(x))+(1.0_sp/x)*poly(y,p)
  else
     y=2.0_sp/x
     bessk1_s=(exp(-x)/sqrt(x))*poly(y,q)
  end if
END FUNCTION bessk1_s

FUNCTION bessk1_v(x)
  use mo_kind
  use mo_nrutil, ONLY : assert,poly
  use mo_nr, ONLY : bessi1
  IMPLICIT NONE
  REAL(SP), DIMENSION(:), INTENT(IN) :: x
  REAL(SP), DIMENSION(size(x)) :: bessk1_v
  REAL(sp), DIMENSION(size(x)) :: y
  LOGICAL(LGT), DIMENSION(size(x)) :: mask
  REAL(sp), DIMENSION(7) :: p = (/1.0_sp,0.15443144_sp,&
       -0.67278579_sp,-0.18156897_sp,-0.1919402e-1_sp,&
       -0.110404e-2_sp,-0.4686e-4_sp/)
  REAL(sp), DIMENSION(7) :: q = (/1.25331414_sp,0.23498619_sp,&
       -0.3655620e-1_sp,0.1504268e-1_sp,-0.780353e-2_sp,&
       0.325614e-2_sp,-0.68245e-3_sp/)
  call assert(all(x > 0.0), 'bessk1_v arg')
  mask = (x <= 2.0)
  where (mask)
     y=x*x/4.0_sp
     bessk1_v=(log(x/2.0_sp)*bessi1(x))+(1.0_sp/x)*poly(y,p,mask)
  elsewhere
     y=2.0_sp/x
     bessk1_v=(exp(-x)/sqrt(x))*poly(y,q,.not. mask)
  end where
END FUNCTION bessk1_v

FUNCTION dbessk1_s(x)
  use mo_kind
  use mo_nrutil, ONLY : assert,poly
  use mo_nr, ONLY : bessi1
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: x
  REAL(DP) :: dbessk1_s
  REAL(DP) :: y
  REAL(DP), DIMENSION(7) :: p = (/1.0_dp,0.15443144_dp,&
       -0.67278579_dp,-0.18156897_dp,-0.1919402e-1_dp,&
       -0.110404e-2_dp,-0.4686e-4_dp/)
  REAL(DP), DIMENSION(7) :: q = (/1.25331414_dp,0.23498619_dp,&
       -0.3655620e-1_dp,0.1504268e-1_dp,-0.780353e-2_dp,&
       0.325614e-2_dp,-0.68245e-3_dp/)
  call assert(x > 0.0, 'dbessk1_s arg')
  if (x <= 2.0) then
     y=x*x/4.0_dp
     dbessk1_s=(log(x/2.0_dp)*bessi1(x))+(1.0_dp/x)*poly(y,p)
  else
     y=2.0_dp/x
     dbessk1_s=(exp(-x)/sqrt(x))*poly(y,q)
  end if
END FUNCTION dbessk1_s

FUNCTION dbessk1_v(x)
  use mo_kind
  use mo_nrutil, ONLY : assert,poly
  use mo_nr, ONLY : bessi1
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), INTENT(IN) :: x
  REAL(DP), DIMENSION(size(x)) :: dbessk1_v
  REAL(DP), DIMENSION(size(x)) :: y
  LOGICAL(LGT), DIMENSION(size(x)) :: mask
  REAL(DP), DIMENSION(7) :: p = (/1.0_dp,0.15443144_dp,&
       -0.67278579_dp,-0.18156897_dp,-0.1919402e-1_dp,&
       -0.110404e-2_dp,-0.4686e-4_dp/)
  REAL(DP), DIMENSION(7) :: q = (/1.25331414_dp,0.23498619_dp,&
       -0.3655620e-1_dp,0.1504268e-1_dp,-0.780353e-2_dp,&
       0.325614e-2_dp,-0.68245e-3_dp/)
  call assert(all(x > 0.0), 'dbessk1_v arg')
  mask = (x <= 2.0)
  where (mask)
     y=x*x/4.0_dp
     dbessk1_v=(log(x/2.0_dp)*bessi1(x))+(1.0_dp/x)*poly(y,p,mask)
  elsewhere
     y=2.0_dp/x
     dbessk1_v=(exp(-x)/sqrt(x))*poly(y,q,.not. mask)
  end where
END FUNCTION dbessk1_v


text:
      subroutine trsurf(x,y,z,n0,xv,yv,nodes,net0,
     1  iscr,scrn,idim1,idim2,net,ibound)
C...01/09/78
      common /gftrg1/job,xloc,yloc,step,dop,itrgl(3)
      dimension ibound(2,1),x(1),y(1),z(1),net0(1),net(1),nodes(1),
     1 iscr(1),scrn(idim1,idim2,2)
      job=0
      k=iabs(nodes(n0+1))-1
      do 7 i=1,k
    7 net(i)=net0(i)
      nb=0
      do 11 ir=1,n0
      k1=iabs(nodes(ir))
      k2=iabs(nodes(ir+1))-1
      do 10 j=k1,k2
      l=net(j)
      if(l.ge.0) goto 10
      l=iabs(l)
      r=cosin1(y(ir)-y(l),x(l)-x(ir),x(l)-xv,y(l)-yv)
      if(r.gt.0) goto 11
      nb=nb+1
      ibound(1,nb)=l
      ibound(2,nb)=ir
      go to 11
   10 continue
   11 continue
   12 k=iabs(nodes(n0+1))-1
      call nextrg(ibound,x,y,net,nodes,nb,xv,yv)
      call trdraw(iscr,scrn,idim1,idim2,x,y,z)
      if(nb.gt.0) goto 12
      do 15 i=1,n0
   15 nodes(i)=iabs(nodes(i))
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


SUBROUTINE CHPEVX_F95( A, W, UPLO, Z, VL, VU, IL, IU, &
                       M, IFAIL, ABSTOL, INFO )
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS ..
   USE LA_PRECISION, ONLY: WP => SP
   USE LA_AUXMOD, ONLY: ERINFO, LSAME
   USE F77_LAPACK, ONLY: LAMCH_F77 => SLAMCH
   USE F77_LAPACK, ONLY: HPEVX_F77 => LA_HPEVX
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
!  .. CHARACTER ARGUMENTS ..
   CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
!  .. SCALAR ARGUMENTS ..
   INTEGER, INTENT(IN), OPTIONAL :: IL, IU
   INTEGER, INTENT(OUT), OPTIONAL :: INFO, M
   REAL(WP), INTENT(IN), OPTIONAL :: ABSTOL, VL, VU
!  .. ARRAY ARGUMENTS ..
   INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IFAIL(:)
   COMPLEX(WP), INTENT(OUT), OPTIONAL, TARGET :: Z(:,:)
   COMPLEX(WP), INTENT(INOUT) :: A(:)
   REAL(WP), INTENT(OUT) :: W(:)
!----------------------------------------------------------------------
!
! Purpose
! =======
!
!      LA_SPEVX / LA_HPEVX compute selected eigenvalues and, optionally,
! the corresponding eigenvectors of a real symmetric/complex hermitian
! matrix A in packed storage. Eigenvalues and eigenvectors can be
! selected by specifying either a range of values or a range of indices
! for the desired eigenvalues.
!
! =========
!
!        SUBROUTINE LA_SPEVX / LA_HPEVX( AP, W, UPLO=uplo, Z=z, &
!                 VL=vl, VU=vu, IL=il, IU=iu, M=m, IFAIL=ifail, &
!                 ABSTOL=abstol, INFO=info )
!          <type>(<wp>), INTENT(INOUT) :: AP(:)
!          REAL(<wp>), INTENT(OUT) :: W(:)
!          CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
!          <type>(<wp>), INTENT(OUT), OPTIONAL :: Z(:,:)
!          REAL(<wp>), INTENT(IN), OPTIONAL :: VL, VU
!          INTEGER, INTENT(IN), OPTIONAL :: IL, IU
!          INTEGER, INTENT(OUT), OPTIONAL :: M
!          INTEGER, INTENT(OUT), OPTIONAL :: IFAIL(:)
!          REAL(<wp>), INTENT(IN), OPTIONAL :: ABSTOL
!          INTEGER, INTENT(OUT), OPTIONAL :: INFO
!        where
!          <type> ::= REAL | COMPLEX
!          <wp> ::= KIND(1.0) | KIND(1.0D0)
!
! Arguments
! =========
!
!  AP      (input/output) REAL or COMPLEX array, shape (:) with size(AP)=
!          n*(n+1)/2, where n is the order of A.
!          On entry, the upper or lower triangle of matrix A in packed
!          storage. The elements are stored columnwise as follows:
!          if UPLO = 'U', AP(i+(j-1)*j/2)=A(i,j) for 1<=i<=j<=n;
!          if UPLO = 'L', AP(i+(j-1)*(2*n-j)/2)=A(i,j) for 1<=j<=i<=n.
!          On exit, AP is overwritten by values generated during the
!          reduction of A to a tridiagonal matrix T . If UPLO = 'U', the
!          diagonal and first superdiagonal of T overwrite the correspond-
!          ing diagonals of A. If UPLO = 'L', the diagonal and first
!          subdiagonal of T overwrite the corresponding diagonals of A.
!  W       (output) REAL array, shape (:) with size(W) = n.
!          The eigenvalues in ascending order.
!  UPLO    Optional (input) CHARACTER(LEN=1).
!          = 'U': Upper triangle of A is stored;
!          = 'L': Lower triangle of A is stored.
!          Default value: 'U'.
! Z        Optional (output) REAL or COMPLEX array, shape (:,:) with
!          size(Z,1) = n and size(Z,2) = M.
!          The first M columns of Z contain the orthonormal eigenvectors of
! 	   the matrix A corresponding to the selected eigenvalues, with the
! 	   i-th column of Z containing the eigenvector associated with
!          the eigenvalue in W(i) . If an eigenvector fails to converge,
! 	   then that column of Z contains the latest approximation to the
!    	   eigenvector, and the index of the eigenvector is returned in
! 	   IFAIL.
!          Note: The user must ensure that at least M columns are supplied
! 	   in the array Z. When the exact value of M is not known in
! 	   advance, an upper bound must be used. In all cases M <= n.
! VL,VU    Optional (input) REAL.
!          The lower and upper bounds of the interval to be searched for
! 	   eigenvalues. VL < VU.
!          Default values: VL = -HUGE(<wp>) and VU = HUGE(<wp>), where
! 	   <wp> ::= KIND(1.0) | KIND(1.0D0).
!          Note: Neither VL nor VU may be present if IL and/or IU is
! 	   present.
! IL,IU    Optional (input) INTEGER.
!          The indices of the smallest and largest eigenvalues to be
!          returned. The IL-th through IU-th
!          eigenvalues will be found. 1<=IL<=IU<=size(A,1).
!          Default values: IL = 1 and IU = size(A,1).
!          Note: Neither IL nor IU may be present if VL and/or VU is
!          present.
!          Note: All eigenvalues are calculated if none of the arguments
! 	   VL, VU, IL and IU are present.
! M        Optional (output) INTEGER.
!          The total number of eigenvalues found. 0<=M<=size(A,1).
!          Note: If IL and IU are present then M = IU - IL + 1.
! IFAIL    Optional (output) INTEGER array, shape (:) with
!          size(IFAIL) = n.
!          If INFO = 0, the first M elements of IFAIL are zero.
!          If INFO > 0, then IFAIL contains the indices of the
! 	   eigenvectors that failed to converge.
!          Note: If Z is present then IFAIL should also be present.
! ABSTOL   Optional (input) REAL.
!          The absolute error tolerance for the eigenvalues. An
! 	   approximate eigenvalue is accepted as converged when it is
! 	   determined to lie in an interval [a,b] of width less than or
! 	   equal to ABSTOL+EPSILON(1.0_<wp>) * max(|a|,|b|),
!          where <wp> is the working precision. If ABSTOL<=0, then
! 	   EPSILON(1.0_<wp>)*||T||1 will be used in its place, where
! 	   ||T||1 is the l1 norm of the tridiagonal matrix obtained by
! 	   reducing A to tridiagonal form. Eigenvalues will be computed
! 	   most accurately when ABSTOL is set to twice the underflow
! 	   threshold 2*LA_LAMCH(1.0_<wp>, 'Safe minimum'), not zero.
!          Default value: 0.0_<wp>.
!          Note: If this routine returns with INFO > 0, then some
! 	   eigenvectors did not converge. Try setting ABSTOL to
! 	   2*LA_LAMCH(1.0_<wp>, 'Safe minimum').
! INFO     Optional (output) INTEGER.
!          = 0: successful exit.
!          < 0: if INFO = -i, the i-th argument had an illegal value.
!          > 0: if INFO = i, then i eigenvectors failed to converge. Their
! 	   indices are stored in array IFAIL.
!          If INFO is not present and an error occurs, then the program is
! 	   terminated with an error message.
!----------------------------------------------------------------------
!  .. LOCAL PARAMETERS ..
   CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_HPEVX'
   CHARACTER(LEN=1) :: LJOBZ, LUPLO, LRANGE
!  .. LOCAL SCALARS ..
   INTEGER :: N, LINFO, LD, LIL, LIU, LM, ISTAT, &
              SIFAIL, S1Z, S2Z, NN
   INTEGER, TARGET :: ISTAT1(1)
   COMPLEX(WP), TARGET :: LLZ(1,1)
   REAL(WP) :: LABSTOL, LVL, LVU
   COMPLEX(WP) :: WW
!  .. LOCAL ARRAYS ..
   INTEGER, POINTER :: IWORK(:), LIFAIL(:)
   COMPLEX(WP), POINTER :: WORK(:)
   REAL(WP), POINTER :: RWORK(:)
!  .. INTRINSIC FUNCTIONS ..
   INTRINSIC HUGE, PRESENT, SIZE
!  .. EXECUTABLE STATEMENTS ..
   LINFO = 0; ISTAT = 0; NN = SIZE(A)
   WW = (-1+SQRT(1+8*REAL(NN,WP)))*0.5; N = INT(WW);  LD = MAX(1,N)
   IF( PRESENT(IFAIL) )THEN; SIFAIL = SIZE(IFAIL); ELSE; SIFAIL = N; END IF
   IF( PRESENT(UPLO) ) THEN; LUPLO = UPLO; ELSE; LUPLO = 'U'; END IF
   IF( PRESENT(VL) )THEN; LVL = VL; ELSE; LVL = -HUGE(LVL); ENDIF
   IF( PRESENT(VU) )THEN; LVU = VU; ELSE; LVU = HUGE(LVU); ENDIF
   IF( PRESENT(IL) )THEN; LIL = IL; ELSE; LIL = 1; ENDIF
   IF( PRESENT(IU) )THEN; LIU = IU; ELSE; LIU = N; ENDIF
   IF( PRESENT(Z) )THEN; S1Z = SIZE(Z,1); S2Z = SIZE(Z,2)
    ELSE; S1Z = 1; S2Z = 1; ENDIF
!  .. TEST THE ARGUMENTS
   IF( NN < 0 .OR. AIMAG(WW) /= 0 .OR. REAL(N,WP) /= REAL(WW) ) THEN; LINFO = -1
   ELSE IF( SIZE( W ) /= N )THEN; LINFO = -2
   ELSE IF( .NOT.LSAME(LUPLO,'U') .AND. .NOT.LSAME(LUPLO,'L') )THEN; LINFO = -3
   ELSE IF( PRESENT(Z) .AND. ( S1Z /= LD .OR. S2Z /= N ) )THEN; LINFO = -4
   ELSE IF( LVU < LVL )THEN; LINFO = -5
   ELSE IF( (PRESENT(VL) .OR. PRESENT(VU)) .AND. &
     &      (PRESENT(IL) .OR. PRESENT(IU)) )THEN; LINFO = -6
   ELSE IF(( LIU < LIL .OR. LIL < 1) .AND. N > 0 )THEN; LINFO = -7
   ELSE IF( N < LIU )THEN; LINFO = -8
   ELSE IF( SIFAIL /= N .OR. PRESENT(IFAIL).AND..NOT.PRESENT(Z) )THEN; LINFO = -10
   ELSE IF( N > 0 )THEN
      IF( PRESENT(VL) .OR. PRESENT(VU) )THEN; LRANGE = 'V'; LM = N
      ELSE IF( PRESENT(IL) .OR. PRESENT(IU) )THEN; LRANGE = 'I'; LM = LIU-LIL+1
      ELSE; LRANGE = 'A'; LM = N; END IF
      IF( PRESENT(Z) ) THEN; LJOBZ = 'V'
         IF( PRESENT(IFAIL) )THEN; LIFAIL => IFAIL
         ELSE; ALLOCATE( LIFAIL(N), STAT=ISTAT ); END IF
      ELSE; LJOBZ = 'N'; LIFAIL => ISTAT1; ENDIF
!     .. DETERMINE THE WORKSPACE
      IF( ISTAT == 0 ) THEN
         ALLOCATE(IWORK(MAX(1,5*N)), RWORK(MAX(1,7*N)), WORK(MAX(1,2*N)), STAT=ISTAT)
         IF( ISTAT == 0 )THEN
            IF( PRESENT(ABSTOL) )THEN; LABSTOL = ABSTOL
            ELSE; LABSTOL = 2*LAMCH_F77('Safe minimum'); ENDIF
	    IF (PRESENT (Z)) THEN
                CALL HPEVX_F77( LJOBZ, LRANGE, LUPLO, N, A, LVL, LVU, &
     &                       LIL, LIU, LABSTOL, LM, W, Z, S1Z, WORK, &
     &                       RWORK, IWORK, LIFAIL, LINFO )
            ELSE
	       CALL HPEVX_F77( LJOBZ, LRANGE, LUPLO, N, A, LVL, LVU, &
     &                       LIL, LIU, LABSTOL, LM, W, LLZ, S1Z, WORK, &
     &                       RWORK, IWORK, LIFAIL, LINFO )
            ENDIF
            IF( PRESENT(M) ) M = LM
            W(LM+1:N) = 0.0_WP
         ELSE; LINFO = -100; END IF
      END IF
      IF( PRESENT(Z) .AND. .NOT.PRESENT(IFAIL) ) DEALLOCATE( LIFAIL, STAT=ISTAT1(1) )
      DEALLOCATE(IWORK, RWORK, WORK, STAT=ISTAT1(1))
   END IF
   CALL ERINFO(LINFO,SRNAME,INFO,ISTAT)
END SUBROUTINE CHPEVX_F95


text:
      SUBROUTINE DTGFND(X,Y, TRIANG,NT,Q,INDTRI,TRI,S,MODE)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1997-06-15 DTGFND CLL In args: Removed NP. Changed meaning of MODE.
c>> 1996-03-30 DTGFND  Krogh  Removed Fortran 90 comments.
C>> 1996-02-02 DTGFND CLL
C>> 1995-09-26 DTGFND CLL Editing for inclusion into MATH77.
C
c     C.L.LAWSON, JPL, 1976 DEC 3
c         LOOK UP POINT Q IN TRIANGULAR GRID.
c
c   X(1:NP), Y(1:NP) [in]  (x,y) coordinates of vertices of the
c                          triangular grid.
c
c         {NP is the number of vertices in the triangular grid, but is
c          not explicitly used in this subroutine.}
c
c   TRIANG(1:6*NT) [in]  Array of integer pointers defining the
c            connectivity of the triangular grid.
c
c   NT [in]  No. of triangles in the triangular grid.
c
c   Q(1:2) [in]  The (x,y) coordinates of the point for which this
c             subr will attempt to find an enclosing triangle.
c
c   INDTRI [in]  DESIGNATES TRIANGLE AT WHICH SEARCH WILL BEGIN.
c              ON RETURN INDTRI IS THE INDEX OF THE LAST TRIANGLE
c              TESTED.  THIS IS THE TRIANGLE CONTAINING Q IF MODE=0.
c              IF MODE = 1, 2, or 3, THIS IS A BOUNDARY TRIANGLE AND Q
c              IS OUTSIDE A BOUNDARY EDGE OF THIS TRIANGLE.
C
c   TRI(1:7) [out]  INTEGER ARRAY OF POINTERS ASSOCIATED
c              WITH TRIANGLE INDTRI.
c
c   S(1:3, 1:3) [out]  If MODE=0 on return, S(,) CONTAINS
c              S(:,1) = UNNORMALIZED BARYCENTRIC COORDINATES OF Q.
c              S(:,2) = X COORDINATES OF EDGE VECTORS.
c              S(:,3) = Y COORDINATES OF EDGE VECTORS.
C
c   MODE [out]
c      = 0  MEANS OK.  Q IS INTERIOR OR ALMOST SO.
c      = 1, 2, or 3.   Means Q is exterior by more than the built-in
c           tolerance. The value of MODE identifies the edge of the
c           triangle relative to which Q is outside.  TRI(MODE) will be
c           zero, indicating a boundary edge.  Q is outside this edge.
c           The vertices at the ends of this edge, in counterclockwise
c           order are TRI(3+MODE) and TRI(4+MODE).
c
c      = -1  BAD.  SUBR IS CYCLING IN THE SEARCH.  THIS SHOULD NEVER
c              HAPPEN.
C
c     ------------------------------------------------------------------
c         Details of the contents of TRI(1:7) and S(1:3, 1:3).
c
c  For descriptive convenience we regard the subscript of P() and the
c  1st subsubscript of S(,) as always being reduced modulo 3 to 1, 2,
c  or 3.  Also for convenience we shall write P(i) to mean the vertex
c  indexed by P(i).
c
c  TRI() contains integer pointers defining the triangle indexed by
c  INDTRI.  The indices of the vertices of this triangle, in counter-
c  clockwise order are P(1), P(2), and P(3) which may be obtained as
c  P(i) = TRI(3+i) for i = 1, 2, and 3.  TRI(7) contains the same
c  value as TRI(4).  The triangle adjacent to this triangle across
c  the edge from P(3+i) to P(4+i) is indexed by TRI(i) for i = 1, 2,
c  and 3.  If there is no adjacent triangle across this edge then
c  TRI(i) = 0.
c
c  The unnormalized barycentric coordinate that is zero along the edge
c  from P(i) to P(i+1) and has a positive value at P(i+2) is stored in
c  S(i,1).  The (x,y) coordinates of the vector from P(i) to P(i+1) are
c  stored in (S(i,2), S(i,3)).
c     ------------------------------------------------------------------
c--D replaces "?": ?TGFND, ?TGGET
c     ------------------------------------------------------------------
      integer icount, indtri, iskip, j, j1, j2, j2save, MODE, NT
      integer TRI(7), TRIANG(1:*)
      double precision p1, p2, Q(2), S(3,3), X(*), Y(*)
c     ------------------------------------------------------------------
      IF( INDTRI .LE. 0 .OR. INDTRI .GT. NT) INDTRI=1
      MODE = -1
      ISKIP  = 0
      J2SAVE = 0
      DO 50 ICOUNT = 1, NT
         CALL DTGGET(INDTRI,TRI,TRIANG)
         TRI(7)=TRI(4)

            DO 10 J=1,3
               J1= TRI(J+3)
               IF ( J1 .EQ. J2SAVE) THEN
                  ISKIP=J
               ELSE
                  J2 = TRI(J+4)
                  S(J,2) = X(J2)-X(J1)
                  P2 = S(J,2) * (Q(2)-Y(J1))
                  S(J,3) = Y(J2)-Y(J1)
                  P1 = S(J,3) * (Q(1)-X(J1))
                  S(J,1) = -P1+P2
C
                  IF (S(J,1) .LT. -1.0D-4*(ABS(P1)+ABS(P2))) THEN
C
c                        Q IS OUTSIDE TRIANGLE INDTRI BY MORE THAN
c                        TOLERANCE.  MOVE TO NEIGHBORING TRIANGLE IF
c                        THERE IS ONE.
C
                     IF (TRI(J) .NE. 0) THEN
                        J2SAVE = J2
                        INDTRI = TRI(J)
                        go to 40
                     ELSE
                        MODE = J
                        go to 80
                     END IF
                  END IF
               END IF
   10       continue
c                                  Q IS INSIDE OR ALMOST INSIDE
c                                  THIS TRIANGLE
            MODE = 0
            go to 90
   40    continue
   50 continue
   80 continue
   90 continue
c                             COMPUTE SKIPPED INFO FOR USE IN
c                             INTERPOLATION SUBR.
      IF (ISKIP .NE. 0) THEN
         J1=TRI(ISKIP+3)
         J2=TRI(ISKIP+4)
         S(ISKIP,2)= X(J2)-X(J1)
         S(ISKIP,3)= Y(J2)-Y(J1)
         S(ISKIP,1)= -S(ISKIP,3)*(Q(1)-X(J1))
     *               +S(ISKIP,2)*(Q(2)-Y(J1))
      END IF
      RETURN
      END


text:
      SUBROUTINE CSD(Z,P,L,M,IWIN,DELT,STATE)
C
      REAL*8 TITLE(10)
      COMPLEX Z(1)
C
      COMMON/DEBUG1/ IAFM(5)
C
      COMMON/INPFFT/ICSD
C
      DIMENSION STORE(1024)
      DIMENSION F(10),AMP(10),STATE(7)
      DIMENSION P(1)
C
      DATA TWOPI /6.283185/
C
      DATA TITLE/'FIRST   ','SECOND  ','THIRD   ','FOURTH  ','FIFTH   ',
     1           'SIXTH   ','SEVENTH ','EIGTH   ','NINTH   ','TENTH   '/
C
      CALL FINDN(L,M,N,NS)
C
      L1=L+1
      LL=L+L
      IF(IAFM(4).EQ.0) GO TO 1975
    2 FORMAT(9H RECORD X)
      WRITE(6,2)
C     READ IN A SEQUENCE X (OR Y).
    1 FORMAT(1PE14.6,1PE13.6,1PE13.6,1PE13.6,1PE13.6,1PE13.6,1PE13.6
     &              ,1PE13.6,1PE13.6,1PE13.6)
      WRITE(6,1) (P(I),I=1,L)
 1975 CONTINUE
C     1. MAKE X INTO ZERO MEAN SEQUENCE AND STORE NEW X.
      CALL ZEROMN(P,L,SA)
      DO 88 I=1,L
   88 STORE(I)=P(I)
C     2. ADD ZERCES TO X TIL SEQUENCE LENGTH OF X BECOMES N.
      DO 99 I=L1,N
   99 P(I)=0.
      CALL AUTPOW(Z,P,L,N,SIG,M,NS,IWIN,DELT,DELF)
      CALL COMMAX(Z,NS,F,AMP,ICSD)
C
      DO 73 I=1,ICSD
      F(I)=F(I)*DELF
      WRITE(6,600) TITLE(I),F(I)
  600 FORMAT('0',10X,'*****  ',A8,'FREQUENCY = ',1PE13.6,' CPS  *****')
      F(I)=F(I)*TWOPI
   73 CONTINUE
C
      DO 50 J=1,ICSD
      FI=F(J)*DELT
      NT=TWOPI/FI
C     ARRANGE ZERO MEAN X AND Y ALTERNATELY,I.E., X1,Y1,X2,Y2,...,SL,YL.
      DO 44 I=1,L
      II=I+I
      P(II-1)=STORE(I)
      P(II)=SIN(FI*(I-1))
   44 CONTINUE
      LL1=LL+1
      N2=N+N
C     ADD ZEROES TO XY SEQUENCE TIL SEQUENCE LENGTH BECOMES 2*N.
      DO 111 I=LL1,N2
  111 P(I)=0.
C     COMPUTE CROSSCORRELATION
C
      CALL CROSS(Z,L,N)
      CALL PLTC(Z,1,NT,1,1,0,0)
      CALL COMMAX(Z,NT,TAW,AMPL,1)
C
      IF(TAW .EQ. 0.0) AMPL=REAL(Z(1))
C
      TAW=TAW*FI*360.0E0/TWOPI
C
      WRITE(6,601) TITLE(J),TAW,TITLE(J),AMPL
  601 FORMAT('0',20X,'PHASE ANGLE FOR ',A8,'MODE =',F8.4,' DEGREES',/
     1          ,20X,'AMPLITUDE   FOR ',A8,'MODE =',1PE13.6,/)
C
   50 CONTINUE
      RETURN
      END


text:
module declare

   implicit none

   integer(kind=8)::count_rate,cr,timestart,time0,timerun
   real(kind=8)::rate

   integer(kind=8)::nstep,jstep,jskip

   integer jobid,fanum,falen,myp2num,myo2num,myp2len,myo2len,nmemb,nfa,nbondfa,nmyp2,nmyo2,nbondmyp2,nbondmyo2,maxfalen
   integer nf,nb,jx,ja,ja1,ja2,jf,n,length,jdir,nm,jm,n1,n2,nstart,jfile,junit,junit1,junit2
   integer j,nl,jl,nltem,fanummax,nwalrad,nwalrep,nwalsurf,nwalsurf_p,nwalsurf_x,nlayer
   integer nrand,jrand,jreset,neinum,neinum5,ntotal,nbinds,nframe,jforce1,jforce2,jforce3
   integer jprint,nrforce,jrforce,crlknum1,crlknum2,crlknum,nlk,lklen,nbondlk,nnode,nabound,ntether,ntethertem
   integer days,hours,mins,secs,nthreads,jf1,jf2,j1,j2
   integer npair5_myp2ac,npair5_lkac,npair5_myp2lk,npair5_ac2,npair5_myp2,npair5_lk2,npair5_mb2,npair5_mymy
   integer npair_myp2ac,npair_lkac,npair_myp2lk,npair_ac2,npair_myp2,npair_lk2,npair_mymy
   integer npair5_myo2ac,npair5_myo2lk,npair5_myo2,npair_myo2ac,npair_myo2lk,npair_myo2,npair_mb2,nnei
   integer nhoop,nwall,nw,jp,nmeho,nbondwal,jm0,npick,nphi
   integer nxsol,nphisol,jh,jupdate,tension1,tension2,crlkorient,j_dep,cofilin,jmyoturn,jxldel,jxlturnover
   integer natom,natom_a1,natom_a2,natom_a2m,nact1,nact2,nfamax1,nfamax2,fanummax1,fanummax2
   integer jstart,nfahalf,nfafull,astart0,nfa1,nfa2,fanum1,fanum2,jbreak,nact,nmono
   integer natp,jsignal,crlknummax,crlknumactive,jget,nturnover,nmyoturn

   integer,allocatable,dimension(:)::astart,alen,adir,mydir,filid,my2mem,mbhooplen,waltyp,a2mem,boundtyp
   integer,allocatable,dimension(:,:)::bondfa,bondmyp2,myp2body,myp2head,myp2typ,bondmyo2,myo2body,myo2head,myo2typ
   integer,allocatable,dimension(:,:)::fa2myp2,fa2myo2,fa2lk,bondlk,apar,wposid,bondwal,mbhoop
   integer,allocatable,dimension(:,:,:)::myp2nei,myo2nei,lknei,myp2nei5,myo2nei5,lknei5
   integer,allocatable,dimension(:,:)::pair5_myp2ac,pair5_myo2ac,pair5_lkac,pair5_myp2lk,pair5_myo2lk,pair5_ac2
   integer,allocatable,dimension(:,:)::pair5_myp2,pair5_myo2,pair5_lk2,pair5_mb2,pair_mymy,pair5_mymy
   integer,allocatable,dimension(:,:)::pair_myp2ac,pair_myo2ac,pair_lkac,pair_myp2lk,pair_myo2lk
   integer,allocatable,dimension(:,:)::pair_ac2,pair_myp2,pair_myo2,pair_lk2
   integer,allocatable,dimension(:,:)::pair_mb2,pairpart
   integer,allocatable,dimension(:,:)::jfasol,jmyp2sol,jmyo2sol,jlksol,jmbsol,jsursol,nsurf
   integer,allocatable,dimension(:)::fadist,myp2dist,myo2dist,lkdist,fa1stbound,apos,lktyp,lkstart,mark
!   integer,allocatable,dimension(:)::fa1stboundtem,alentem,astarttem

   double precision walrad,mbrad,rrad,rwid,rthick,mbrad2,invmbrad,mwid,kwall,lsqueez,wallbead
   double precision l_lk1,l_lk2,l_lk,k_lk,klk_thet,thet_lk,phi_node_size,l_mem,k_mem,k_orient,kblock
   double precision pi,k_a,l_a,invl_a,ka_thet,thet_fa,kmh_thet,kmb_thet,thet_mh1,thet_mh2,thet_mb,l_mh,l_mb
   double precision xh,yh,zh,xgap,prob,xboundmin,xboundmax,xrmin,xrmax,shift
   double precision ltether,ktether,pres,fmemb,fturgor,df,mbradnew,p_tether
   double precision rvdw,kvdw,r_on,r_off,r_on2,r_off2,fvdwmax,krad
   double precision kpair,l_pair,rpair,thetpair,thet2by2,cos_t_2,kmemb!_thet,thet_memb
   double precision rvdw_mb,r_on_mb,r_off_mb,r_on2_mb,r_off2_mb,fvdwmax_mb
   double precision dphi,x0,y0,z0,xdel,r,phi0,phi,thick,rad,dx,dy,dz,dytem,dztem,dist
   double precision dx1,dy1,dz1,dx2,dy2,dz2,dx3,dy3,dz3,dx4,dy4,dz4
   double precision sinp,cosp,rxy,rxz,ryy,ryz,rzy,rzz
   double precision xmin,xmax,wthick,d2,dist2,fmax2,runtime,oldtime,printtime,dt,gap,kmy2mem,lmy2mem
   double precision dxsol,dphisol,wallrate,ringthick,aconc1,aconc2,lkconc,invmyp2len2,invmyo2len2

   double precision kbind,lbind,p1_hydr,p1_hydr0,p2_bind,p3_pi_rele,p4_adp_rele,p5_ubind,p_scale,p_dep,p_rate
   double precision delta,invdelta,beta,gam,fmaxfa2,fmaxmyp22,fmaxmyo22,fmaxlk2,plk_bind,fmaxmb2
   double precision plk_ubind1,plk_ubind2,plk_off,k_scale,k_scale_lk,plk_remove,pmyturn,plk_turn
   double precision halftime

   real,allocatable,dimension(:)::xwalsurf,ywalsurf,zwalsurf
   real,allocatable,dimension(:)::xwalrad,ywalrad,zwalrad,xwalrep,ywalrep,zwalrep

   double precision,allocatable,dimension(:)::xfa,yfa,zfa,xmyp2,ymyp2,zmyp2,xmyo2,ymyo2,zmyo2
   double precision,allocatable,dimension(:)::rwhoop,xwhoop
   double precision,allocatable,dimension(:)::xmemb,ymemb,zmemb
   double precision,allocatable,dimension(:,:)::xwall,ywall,zwall,rwall,rmemb,rmembmax
   double precision,allocatable,dimension(:)::xlk,ylk,zlk
   double precision,allocatable,dimension(:)::fxfa,fyfa,fzfa,fxfarep,fyfarep,fzfarep
   double precision,allocatable,dimension(:)::fxanglfa,fyanglfa,fzanglfa
   double precision,allocatable,dimension(:)::fxmyp2,fymyp2,fzmyp2,fxmyo2,fymyo2,fzmyo2
   double precision,allocatable,dimension(:)::fxmyp2rep,fymyp2rep,fzmyp2rep,fxmyo2rep,fymyo2rep,fzmyo2rep
   double precision,allocatable,dimension(:)::fxanglmyp2,fyanglmyp2,fzanglmyp2,fxanglmyo2,fyanglmyo2,fzanglmyo2
   double precision,allocatable,dimension(:)::fxlk,fylk,fzlk,fxlkrep,fylkrep,fzlkrep
   double precision,allocatable,dimension(:)::fxangllk,fyangllk,fzangllk
   double precision,allocatable,dimension(:)::rands,phi_node
   double precision,allocatable,dimension(:,:)::rxfa,ryfa,rzfa,rxmyp2,rymyp2,rzmyp2,rxmyo2,rymyo2,rzmyo2,rxlk,rylk,rzlk
   double precision,allocatable,dimension(:,:)::rxmemb,rymemb,rzmemb
!   double precision,allocatable,dimension(:)::fymemfa,fzmemfa,fymemmyo,fzmemmyo,fymemlk,fzmemlk
   double precision,allocatable,dimension(:)::fxmemb,fymemb,fzmemb
!   double precision,allocatable,dimension(:)::fxoll_vdw,fyoll_vdw,fzoll_vdw
!   double precision,allocatable,dimension(:)::fxill_vdw,fyill_vdw,fzill_vdw
   double precision,allocatable,dimension(:)::fxmemb_vdw,fymemb_vdw,fzmemb_vdw

!   double precision,allocatable,dimension(:)::fxfa_vdw,fyfa_vdw,fzfa_vdw
!   double precision,allocatable,dimension(:)::fxmyp2_vdw,fymyp2_vdw,fzmyp2_vdw
!   double precision,allocatable,dimension(:)::fxlk_vdw,fylk_vdw,fzlk_vdw
!   double precision,allocatable,dimension(:)::fxmemb_cst,fymemb_cst,fzmemb_cst
!   double precision,allocatable,dimension(:)::fxfa_cst,fyfa_cst,fzfa_cst
!   double precision,allocatable,dimension(:)::fxmyp2_cst,fymyp2_cst,fzmyp2_cst
!   double precision,allocatable,dimension(:)::fxlk_cst,fylk_cst,fzlk_cst
!   double precision,allocatable,dimension(:,:)::rmbsol,rwlsol

   double precision,allocatable,dimension(:,:)::xnorwall,ynorwall,znorwall,xnorsurf,ynorsurf,znorsurf
   double precision,allocatable,dimension(:,:)::xsurf,ysurf,zsurf
   real fa1,fa2,fa3,fa4,fa5

   character(len=91)::a91

   real timecheck,myorate

end module

max_stars_repo_name: nguyenthlam/ActoMyosin

max_issues_repo_name: nguyenthlam/ActoMyosin

max_forks_repo_name: nguyenthlam/ActoMyosin