text:
!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2021  National Renewable Energy Laboratory
!
!    This file is part of ROSCO.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!**********************************************************************************************************************************
MODULE SysSubs

    USE ROSCO_Types
    USE Constants

    IMPLICIT NONE

    INTEGER, PARAMETER            :: ConRecL     = 120                               ! The record length for console output.
    INTEGER, PARAMETER            :: CU          = 6                                 ! The I/O unit for the console.  Unit 6 causes ADAMS to crash.
    INTEGER, PARAMETER            :: MaxWrScrLen = 98                                ! The maximum number of characters allowed to be written to a line in WrScr
    LOGICAL, PARAMETER            :: KBInputOK   = .TRUE.                            ! A flag to tell the program that keyboard input is allowed in the environment.
    CHARACTER(*),  PARAMETER      :: NewLine     = ACHAR(10)                         ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}]
    CHARACTER(*),  PARAMETER      :: OS_Desc     = 'Intel Fortran for Linux'         ! Description of the language/OS
    CHARACTER( 1), PARAMETER      :: PathSep     = '/'                               ! The path separator.
    CHARACTER( 1), PARAMETER      :: SwChar      = '-'                               ! The switch character for command-line options.
    CHARACTER(11), PARAMETER      :: UnfForm     = 'UNFORMATTED'                     ! The string to specify unformatted I/O files.

! This module contains routines with system-specific logic and references, including all references to the console unit, CU.
! It also contains standard (but not system-specific) routines it uses.
! SysGnuLinux.f90 is specifically for the GNU Fortran (gfortran) compiler on Linux and macOS.

    CONTAINS
!=======================================================================
    SUBROUTINE LoadDynamicLib ( DLL, ErrStat, ErrMsg )

        ! This SUBROUTINE is used to dynamically load a DLL.

        TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be loaded.
        INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
        CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

#ifdef USE_DLL_INTERFACE

     !bjj: these are values I found on the web; I have no idea if they actually work...
     !bjj: hopefully we can find them pre-defined in a header somewhere
        INTEGER(C_INT), PARAMETER :: RTLD_LAZY=1            ! "Perform lazy binding. Only resolve symbols as the code that references them is executed. If the symbol is never referenced, then it is never resolved. (Lazy binding is only performed for function references; references to variables are always immediately bound when the library is loaded.) "
        INTEGER(C_INT), PARAMETER :: RTLD_NOW=2             ! "If this value is specified, or the environment variable LD_BIND_NOW is set to a nonempty string, all undefined symbols in the library are resolved before dlopen() returns. If this cannot be done, an error is returned."
        INTEGER(C_INT), PARAMETER :: RTLD_GLOBAL=256        ! "The symbols defined by this library will be made available for symbol resolution of subsequently loaded libraries"
        INTEGER(C_INT), PARAMETER :: RTLD_LOCAL=0           ! "This is the converse of RTLD_GLOBAL, and the default if neither flag is specified. Symbols defined in this library are not made available to resolve references in subsequently loaded libraries."

        INTERFACE !linux API routines
           !bjj see http://linux.die.net/man/3/dlopen
           !    and https://developer.apple.com/library/mac/documentation/Darwin/Reference/ManPages/man3/dlopen.3.html

           FUNCTION dlOpen(filename,mode) BIND(C,NAME="dlopen")
           ! void *dlopen(const char *filename, int mode);
              USE ISO_C_BINDING
              IMPLICIT NONE
              TYPE(C_PTR)                   :: dlOpen
              CHARACTER(C_CHAR), INTENT(IN) :: filename(*)
              INTEGER(C_INT), VALUE         :: mode
           END FUNCTION

        END INTERFACE

        ErrStat = ErrID_None
        ErrMsg = ''

        ! Load the DLL and get the file address:

        DLL%FileAddrX = dlOpen( TRIM(DLL%FileName)//C_NULL_CHAR, RTLD_LAZY )  !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end)

        IF( .NOT. C_ASSOCIATED(DLL%FileAddrX) ) THEN
           ErrStat = ErrID_Fatal
           WRITE(ErrMsg,'(I2)') BITS_IN_ADDR
           ErrMsg  = 'The dynamic library '//TRIM(DLL%FileName)//' could not be loaded. Check that the file '// &
                     'exists in the specified location and that it is compiled for '//TRIM(ErrMsg)//'-bit applications.'
           RETURN
        END IF

        ! Get the procedure address:

        CALL LoadDynamicLibProc ( DLL, ErrStat, ErrMsg )
#else

        ErrStat = ErrID_Fatal
        ErrMsg = ' LoadDynamicLib: Not compiled with -DUSE_DLL_INTERFACE for '//TRIM(OS_Desc)

#endif

        RETURN
     END SUBROUTINE LoadDynamicLib
     !=======================================================================
     SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg )

        ! This SUBROUTINE is used to dynamically load a procedure from a DLL.

        TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be loaded.
        INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
        CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
        INTEGER(IntKi)                            :: i

#ifdef USE_DLL_INTERFACE

        INTERFACE !linux API routines

           !bjj see http://linux.die.net/man/3/dlsym
           !    and https://developer.apple.com/library/mac/documentation/Darwin/Reference/ManPages/man3/dlsym.3.html

           FUNCTION dlSym(handle,name) BIND(C,NAME="dlsym")
           ! void *dlsym(void *handle, const char *name);
              USE ISO_C_BINDING
              IMPLICIT NONE
              TYPE(C_FUNPTR)                :: dlSym ! A function pointer
              TYPE(C_PTR), VALUE            :: handle
              CHARACTER(C_CHAR), INTENT(IN) :: name(*)
           END FUNCTION

        END INTERFACE

        ErrStat = ErrID_None
        ErrMsg = ''

        do i=1,NWTC_MAX_DLL_PROC
           if ( len_trim( DLL%ProcName(i) ) > 0 ) then

              DLL%ProcAddr(i) = dlSym( DLL%FileAddrX, TRIM(DLL%ProcName(i))//C_NULL_CHAR )  !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end)

              IF(.NOT. C_ASSOCIATED(DLL%ProcAddr(i))) THEN
                 ErrStat = ErrID_Fatal + i - 1
                 ErrMsg  = 'The procedure '//TRIM(DLL%ProcName(i))//' in file '//TRIM(DLL%FileName)//' could not be loaded.'
                 RETURN
              END IF

           end if
        end do

#else

        ErrStat = ErrID_Fatal
        ErrMsg = ' LoadDynamicLibProc: Not compiled with -DUSE_DLL_INTERFACE for '//TRIM(OS_Desc)

#endif

        RETURN
     END SUBROUTINE LoadDynamicLibProc
     !=======================================================================
     SUBROUTINE FreeDynamicLib ( DLL, ErrStat, ErrMsg )

           ! This SUBROUTINE is used to free a dynamically loaded DLL (loaded in LoadDynamicLib).

           TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be freed.
           INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
           CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
           INTEGER(C_INT)                            :: Success     ! Whether or not the call to dlClose was successful
           INTEGER(C_INT), PARAMETER                 :: TRUE  = 0

#ifdef USE_DLL_INTERFACE
     !bjj: note that this is not tested.

        INTERFACE !linux API routine
           !bjj see http://linux.die.net/man/3/dlclose
           !    and https://developer.apple.com/library/mac/documentation/Darwin/Reference/ManPages/man3/dlclose.3.html

           FUNCTION dlClose(handle) BIND(C,NAME="dlclose")
           ! int dlclose(void *handle);
              USE ISO_C_BINDING
              IMPLICIT NONE
              INTEGER(C_INT)       :: dlClose
              TYPE(C_PTR), VALUE   :: handle
           END FUNCTION

        END INTERFACE

        ! Close the library:

        IF( .NOT. C_ASSOCIATED(DLL%FileAddrX) ) RETURN
        Success = dlClose( DLL%FileAddrX ) !The function dlclose() returns 0 on success, and nonzero on error.

        IF ( Success /= TRUE ) THEN !bjj: note that this is not the same as LOGICAL .TRUE.
           ErrStat = ErrID_Fatal
           ErrMsg  = 'The dynamic library could not be freed.'
           RETURN
        ELSE
           ErrStat = ErrID_None
           ErrMsg = ''
           DLL%FileAddrX = C_NULL_PTR
        END IF

#else

           ErrStat = ErrID_Fatal
           ErrMsg = ' FreeDynamicLib: Not compiled with -DUSE_DLL_INTERFACE for '//TRIM(OS_Desc)

#endif

           RETURN
        END SUBROUTINE FreeDynamicLib
     !=======================================================================

END MODULE SysSubs

max_stars_repo_name: ptrbortolotti/ROSCO

max_issues_repo_name: ptrbortolotti/ROSCO

max_forks_repo_name: mvanv/ROSCO


! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

! Testing 7.6 enum

  ! OK
  enum, bind(C)
    enumerator :: red, green
    enumerator blue, pink
    enumerator yellow
    enumerator :: purple = 2
  end enum

  integer(yellow) anint4

  enum, bind(C)
    enumerator :: square, cicrle
    !ERROR: 'square' is already declared in this scoping unit
    enumerator square
  end enum

  dimension :: apple(4)
  real :: peach

  enum, bind(C)
    !ERROR: 'apple' is already declared in this scoping unit
    enumerator :: apple
    enumerator :: pear
    !ERROR: 'peach' is already declared in this scoping unit
    enumerator :: peach
    !ERROR: 'red' is already declared in this scoping unit
    enumerator :: red
  end enum

  enum, bind(C)
    !ERROR: Enumerator value could not be computed from the given expression
    !ERROR: Must be a constant value
    enumerator :: wrong = 0/0
  end enum

end


text:
!-------------------------------------------------------------------------------
!
!> \file    GetOutputChemPot.f90
!> \brief   Get specific thermodynamic output.
!> \author  M.H.A. Piro
!> \date    Sept. 16, 2015
!
!
! Revisions:
! ==========
!
!   Date            Programmer          Description of change
!   ----            ----------          ---------------------
!   09/16/2015      M.H.A. Piro         Original code
!
!
! Purpose:
! ========
!
!> \details The purpose of this subroutine is to get specific thermodynamic
!! output from an equilibrium calculation.
!
!
!
! Pertinent variables:
! ====================
!
!> \param[inout]  cElementNameRequest   A three letter character string
!!                                       representing the name of a chemical
!!                                       element.
!> \param[out]    dElementChemPot       A double real scalar representing the
!!                                       chemical potential of this element.
!> \param[out]    INFO                  An integer scalar indicating a successful
!!                                       exit (== 0) or an error (/= 0).
!
!-------------------------------------------------------------------------------

subroutine GetOutputChemPot(cElementNameRequest, dElementChemPot, INFO)

    USE ModuleThermo
    USE ModuleThermoIO

    implicit none

    integer,      intent(out)   :: INFO
    integer                     :: i, j
    real(8),      intent(out)   :: dElementChemPot
    character(3), intent(inout) :: cElementNameRequest

    ! Initialize variables:
    INFO            = 0
    dElementChemPot = 0D0

    ! Only proceed if Thermochimica solved successfully:
    if (INFOThermo == 0) then

        ! Remove trailing blanks:
        cElementNameRequest = TRIM(cElementNameRequest)

        ! Loop through elements to find the one corresponding to the element
        ! being requested:
        j = 0
        LOOP_A: do i = 1, nElements
            if (cElementNameRequest == cElementName(i)) then
                j = i
                exit LOOP_A
            end if
        end do LOOP_A

        ! Check to make sure that the element was found:
        if (j /= 0) then
            ! The element was found in the list.  Store the chemical potential
            ! of this element and convert back to units of J/g-at:
            dElementChemPot = dElementPotential(j) * dIdealConstant * dTemperature
        else
            ! This element was not found.  Report an error:
            INFO = 1
        end if

    else
        ! Record an error with INFO if INFOThermo /= 0.
        INFO = -1
    end if

    return

end subroutine GetOutputChemPot

max_stars_repo_name: ORNL-CEES/thermochimica

max_issues_repo_name: nikolasscuro/thermochimica

max_forks_repo_name: nikolasscuro/thermochimica


text:
      module MOD_PLOTT

      contains

      SUBROUTINE PLOTT (ND,R,T,MODHEAD,JOBNUM)

C***  DIRECT TRANSFER OF THE TEMPERATURE STRATIFICATION T(R) VERSUS LOG(R/R*-1)

      use MOD_ISRCHFLT
      use MOD_ISRCHFLE
      use MOD_PLOTANF
      implicit real*8(a-h,o-z)

C***  TRANSFER OF "YMAX" FROM SUBR. DECSTAR
      COMMON / COMPLOT / YMAX
      DIMENSION X(100),Y(100)
      DIMENSION R(ND),T(ND)
      CHARACTER HEAD1*60,HEAD2*60,MODHEAD*104

C***  INITIALIZATION
      KANAL=2
      OPEN (KANAL,FILE='PLOT')
C      CALL JSYMSET (2LG2,'TRANSFER')
      write (6,*) 'T-STRATIFICATION TO BE ROUTED'

C***  RADIUS RANGE:   2. > LOG(R/R*-1) > -3.5
      LMIN= ISRCHFLT(ND,R(1),1,101.)
      LMAX= ISRCHFLE(ND,R(1),1,1.000316228) - 1
      NDL=LMAX-LMIN+1
      IF (NDL .GT. 100) THEN
         write (6,*) 'T-PLOT IMPOSSIBLE'
         RETURN
         ENDIF
      DO 5 L=1,NDL
      X(L)=LOG10(R(LMAX+1-L)-1.d0)
      Y(L)=T(LMAX+1-L)/1.d3
    5 CONTINUE
C***  PLOT PARAMETER
      XMIN=-3.5
      XMAX=2.
      YMIN=0.
      IF (YMAX .EQ. 0.) YMAX=10.*(INT(Y(1)/10.)+1.d0)
      XSCALE=4.
      YSCALE=15./YMAX
      XTICK=0.5
      YTICK=10.
      XABST=1.
      YABST=10.
      HEAD1=' WR TEMPERATURE STRATIFICATION T(R) VERSUS LOG(R/R*-1) '
c      ENCODE(60,9,HEAD2) JOBNUM
      write (head2,9) jobnum
    9 FORMAT (19X,1HJ,I3)
      HEAD2(1:18)=MODHEAD(15:32)
      HEAD2(24:56)=' TEMPERATURE STRATIFICATION T(R)'
      CALL PLOTANF (KANAL,HEAD1,HEAD2
     $ ,60H      LOG (R/R* - 1)
     $ ,60H      T(R) / KK
     $ ,XSCALE,XMIN,XMAX,XTICK,XABST,.0
     $ ,YSCALE,YMIN,YMAX,YTICK,YABST,.0
     $ ,X,Y,NDL,5)

      RETURN
      END subroutine
      end module

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


      subroutine mgt_killop (jj, iplant)

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine performs the kill operation

      use basin_module
      use organic_mineral_mass_module
      use hru_module, only : hru, ihru, ipl
      use soil_module
      use plant_module
      use constituent_mass_module

      implicit none

      integer :: j                     !none           |HRU number
      integer :: k                     !none           |counter
      integer, intent (in) :: jj       !none           |counter
      integer, intent (in) :: iplant   !               |plant number xwalked from hlt_db()%plant and plants.plt
      integer :: ly                    !none           |soil layer

      j = jj
      ipl = iplant

      !! update root fractions in each layer
      call pl_rootfr

      !! allocate dead roots, N, P to soil layers
	  do ly = 1, soil(j)%nly
	      soil1(j)%tot(ly) = soil(j)%ly(ly)%rtfr * pl_mass(j)%root_com + soil1(j)%tot(ly)
      end do

      !! add above ground mass to residue pool
      rsd1(j)%tot(1) = pl_mass(j)%ab_gr(ipl) + rsd1(j)%tot(1)

      !! zero all plant mass
      pl_mass(j)%tot(ipl) = plt_mass_z
      pl_mass(j)%ab_gr(ipl) = plt_mass_z
      pl_mass(j)%leaf(ipl) = plt_mass_z
      pl_mass(j)%stem(ipl) = plt_mass_z
      pl_mass(j)%seed(ipl) = plt_mass_z
      pl_mass(j)%root(ipl) = plt_mass_z

      do k = 1, cs_db%num_pests
        cs_soil(j)%ly(1)%pest(k) = cs_soil(j)%ly(1)%pest(k) !+ cs_pl(j)%pest(k)
        cs_pl(j)%pest(k) = 0.
      end do

	  !! reset plant variables
      pcom(j)%plg(ipl) = plgz
      pcom(j)%plm(ipl) = plmz
      pcom(j)%plstr(ipl) = plstrz
      !! can't reset entire plcur - harv_num can't be zero'd
      pcom(j)%plcur(ipl)%gro = "n"
      pcom(j)%plcur(ipl)%idorm = "n"
      pcom(j)%plcur(ipl)%phuacc = 0.
      pcom(j)%plcur(ipl)%curyr_mat = 1

      return
      end subroutine mgt_killop


text:
!**************************************************************
!* AceGen    6.702 Windows (4 May 16)                         *
!*           Co. J. Korelc  2013           18 Mar 17 18:58:10 *
!**************************************************************
! User     : Full professional version
! Notebook : dRdFFunction
! Evaluation time                 : 20 s    Mode  : Optimal
! Number of formulae              : 572     Method: Automatic
! Subroutine                      : dRdF1 size: 8737
! Total size of Mathematica  code : 8737 subexpressions
! Total size of Fortran code      : 23060 bytes

!******************* S U B R O U T I N E **********************
SUBROUTINE dRdF1(v,x,props,statev,Fnew,dRdF)
USE SMSUtility
IMPLICIT NONE
LOGICAL b126,b165
DOUBLE PRECISION v(829),x(19),props(9),statev(19),Fnew(9),dRdF(19,9)
v(808)=props(4)*x(1)
v(757)=-(x(10)*x(2))
v(756)=x(5)*x(7)
v(755)=-(x(4)*x(9))
v(754)=x(6)*x(7)
v(753)=x(5)*x(6)-x(3)*x(8)
v(752)=x(10)*x(9)
v(751)=-(x(3)*x(7))
v(750)=-(x(2)*x(6))+x(8)*x(9)
v(749)=x(10)*x(8)
v(748)=-(x(4)*x(5))
v(747)=x(2)*x(3)-x(5)*x(9)
v(746)=-(x(7)*x(8))
v(745)=x(2)*x(4)
v(744)=-(x(10)*x(6))
v(743)=x(3)*x(4)
v(742)=x(12)*x(13)-x(15)*x(19)
v(741)=-(x(13)*x(14))+x(17)*x(19)
v(740)=x(14)*x(16)-x(11)*x(19)
v(739)=-(x(12)*x(16))+x(18)*x(19)
v(738)=x(15)*x(16)-x(13)*x(18)
v(737)=-(x(11)*x(15))+x(17)*x(18)
v(736)=x(11)*x(12)-x(14)*x(18)
v(735)=x(11)*x(13)-x(16)*x(17)
v(734)=x(14)*x(15)-x(12)*x(17)
v(733)=1d0/(v(750)*x(10)+v(747)*x(4)+v(753)*x(7))
v(732)=-statev(1)+x(1)
v(57)=v(733)*(v(743)+v(744))
v(59)=v(733)*(v(745)+v(746))
v(60)=v(733)*v(747)
v(61)=v(733)*(v(748)+v(749))
v(62)=v(733)*v(750)
v(63)=v(733)*(v(751)+v(752))
v(64)=v(733)*v(753)
v(65)=v(733)*(v(754)+v(755))
v(66)=v(733)*(v(756)+v(757))
v(67)=Fnew(1)*v(57)+Fnew(7)*v(63)+Fnew(4)*v(65)
v(758)=2d0*v(67)
v(199)=v(63)*v(758)
v(208)=-v(199)/3d0
v(196)=v(65)*v(758)
v(205)=-v(196)/3d0
v(193)=v(57)*v(758)
v(202)=-v(193)/3d0
v(68)=Fnew(2)*v(59)+Fnew(8)*v(61)+Fnew(5)*v(66)
v(759)=2d0*v(68)
v(218)=v(61)*v(759)
v(227)=-v(218)/3d0
v(215)=v(66)*v(759)
v(224)=-v(215)/3d0
v(212)=v(59)*v(759)
v(221)=-v(212)/3d0
v(69)=Fnew(3)*v(60)+Fnew(9)*v(62)+Fnew(6)*v(64)
v(760)=2d0*v(69)
v(237)=v(62)*v(760)
v(246)=-v(237)/3d0
v(234)=v(64)*v(760)
v(243)=-v(234)/3d0
v(231)=v(60)*v(760)
v(240)=-v(231)/3d0
v(70)=Fnew(4)*v(59)+Fnew(1)*v(61)+Fnew(7)*v(66)
v(761)=2d0*v(70)
v(253)=v(66)*v(67)+v(63)*v(70)
v(250)=v(59)*v(67)+v(65)*v(70)
v(247)=v(61)*v(67)+v(57)*v(70)
v(217)=v(66)*v(761)
v(226)=-v(217)/3d0
v(214)=v(59)*v(761)
v(223)=-v(214)/3d0
v(211)=v(61)*v(761)
v(220)=-v(211)/3d0
v(71)=Fnew(5)*v(60)+Fnew(2)*v(62)+Fnew(8)*v(64)
v(762)=2d0*v(71)
v(263)=v(64)*v(68)+v(61)*v(71)
v(260)=v(60)*v(68)+v(66)*v(71)
v(257)=v(62)*v(68)+v(59)*v(71)
v(236)=v(64)*v(762)
v(245)=-v(236)/3d0
v(233)=v(60)*v(762)
v(242)=-v(233)/3d0
v(230)=v(62)*v(762)
v(239)=-v(230)/3d0
v(72)=Fnew(6)*v(57)+Fnew(3)*v(63)+Fnew(9)*v(65)
v(763)=2d0*v(72)
v(273)=v(65)*v(69)+v(62)*v(72)
v(270)=v(57)*v(69)+v(64)*v(72)
v(267)=v(63)*v(69)+v(60)*v(72)
v(201)=v(65)*v(763)
v(210)=-v(201)/3d0
v(198)=v(57)*v(763)
v(207)=-v(198)/3d0
v(195)=v(63)*v(763)
v(204)=-v(195)/3d0
v(73)=Fnew(7)*v(60)+Fnew(4)*v(62)+Fnew(1)*v(64)
v(764)=2d0*v(73)
v(271)=v(60)*v(67)+v(63)*v(73)
v(268)=v(62)*v(67)+v(65)*v(73)
v(265)=v(64)*v(67)+v(57)*v(73)
v(262)=v(60)*v(70)+v(66)*v(73)
v(259)=v(62)*v(70)+v(59)*v(73)
v(256)=v(64)*v(70)+v(61)*v(73)
v(235)=v(60)*v(764)
v(244)=-v(235)/3d0
v(232)=v(62)*v(764)
v(241)=-v(232)/3d0
v(229)=v(64)*v(764)
v(238)=-v(229)/3d0
v(74)=Fnew(8)*v(57)+Fnew(5)*v(63)+Fnew(2)*v(65)
v(765)=2d0*v(74)
v(272)=v(57)*v(71)+v(64)*v(74)
v(269)=v(63)*v(71)+v(60)*v(74)
v(266)=v(65)*v(71)+v(62)*v(74)
v(254)=v(57)*v(68)+v(61)*v(74)
v(251)=v(63)*v(68)+v(66)*v(74)
v(248)=v(65)*v(68)+v(59)*v(74)
v(200)=v(57)*v(765)
v(209)=-v(200)/3d0
v(197)=v(63)*v(765)
v(206)=-v(197)/3d0
v(194)=v(65)*v(765)
v(203)=-v(194)/3d0
v(75)=Fnew(9)*v(59)+Fnew(6)*v(61)+Fnew(3)*v(66)
v(766)=2d0*v(75)
v(264)=v(59)*v(69)+v(62)*v(75)
v(261)=v(61)*v(69)+v(64)*v(75)
v(258)=v(66)*v(69)+v(60)*v(75)
v(255)=v(59)*v(72)+v(65)*v(75)
v(252)=v(61)*v(72)+v(57)*v(75)
v(249)=v(66)*v(72)+v(63)*v(75)
v(219)=v(59)*v(766)
v(228)=-v(219)/3d0
v(216)=v(61)*v(766)
v(225)=-v(216)/3d0
v(213)=v(66)*v(766)
v(222)=-v(213)/3d0
v(76)=(v(67)*v(67))+(v(72)*v(72))+(v(74)*v(74))
v(88)=-v(76)/3d0
v(77)=(v(68)*v(68))+(v(70)*v(70))+(v(75)*v(75))
v(89)=-v(77)/3d0
v(78)=(v(69)*v(69))+(v(71)*v(71))+(v(73)*v(73))
v(333)=(2d0/3d0)*v(78)+v(88)+v(89)
v(84)=-v(78)/3d0
v(323)=(2d0/3d0)*v(77)+v(84)+v(88)
v(313)=(2d0/3d0)*v(76)+v(84)+v(89)
v(79)=v(67)*v(70)+v(68)*v(74)+v(72)*v(75)
v(767)=2d0*v(79)
v(279)=-(v(767)*v(78))
v(276)=(v(79)*v(79))
v(771)=-v(276)+v(76)*v(77)
v(80)=v(68)*v(71)+v(70)*v(73)+v(69)*v(75)
v(282)=v(767)*v(80)
v(280)=(-2d0)*v(76)*v(80)
v(274)=(v(80)*v(80))
v(769)=-v(274)+v(77)*v(78)
v(81)=v(69)*v(72)+v(67)*v(73)+v(71)*v(74)
v(768)=2d0*v(81)
v(283)=-(v(768)*v(77))
v(774)=v(282)+v(283)
v(281)=v(768)*v(79)
v(773)=v(280)+v(281)
v(277)=v(768)*v(80)
v(772)=v(277)+v(279)
v(275)=(v(81)*v(81))
v(770)=-v(275)+v(76)*v(78)
v(291)=v(201)*v(769)+v(219)*v(770)+v(237)*v(771)+v(255)*v(772)+v(264)*v(773)+v(273)*v(774)
v(290)=v(200)*v(769)+v(218)*v(770)+v(236)*v(771)+v(254)*v(772)+v(263)*v(773)+v(272)*v(774)
v(289)=v(199)*v(769)+v(217)*v(770)+v(235)*v(771)+v(253)*v(772)+v(262)*v(773)+v(271)*v(774)
v(288)=v(198)*v(769)+v(216)*v(770)+v(234)*v(771)+v(252)*v(772)+v(261)*v(773)+v(270)*v(774)
v(287)=v(197)*v(769)+v(215)*v(770)+v(233)*v(771)+v(251)*v(772)+v(260)*v(773)+v(269)*v(774)
v(286)=v(196)*v(769)+v(214)*v(770)+v(232)*v(771)+v(250)*v(772)+v(259)*v(773)+v(268)*v(774)
v(285)=v(195)*v(769)+v(213)*v(770)+v(231)*v(771)+v(249)*v(772)+v(258)*v(773)+v(267)*v(774)
v(284)=v(194)*v(769)+v(212)*v(770)+v(230)*v(771)+v(248)*v(772)+v(257)*v(773)+v(266)*v(774)
v(278)=v(193)*v(769)+v(211)*v(770)+v(229)*v(771)+v(247)*v(772)+v(256)*v(773)+v(265)*v(774)
v(82)=-(v(274)*v(76))-v(275)*v(77)+v(771)*v(78)+v(277)*v(79)
v(304)=1d0/v(82)**0.13333333333333333d1
v(775)=-v(304)/3d0
v(312)=v(291)*v(775)
v(311)=v(290)*v(775)
v(310)=v(289)*v(775)
v(309)=v(288)*v(775)
v(308)=v(287)*v(775)
v(307)=v(286)*v(775)
v(306)=v(285)*v(775)
v(305)=v(284)*v(775)
v(303)=v(278)*v(775)
v(292)=sqrt(v(82))
v(86)=props(2)*(-v(292)+v(82))
v(85)=1d0/v(82)**0.3333333333333333d0
v(806)=props(1)*v(85)
v(803)=props(1)*(v(303)*v(313)+((2d0/3d0)*v(193)+v(220)+v(238))*v(85))
v(802)=props(1)*(v(303)*v(333)+(v(202)+v(220)+(2d0/3d0)*v(229))*v(85))
v(801)=props(1)*(v(303)*v(323)+(v(202)+(2d0/3d0)*v(211)+v(238))*v(85))
v(800)=props(1)*(v(305)*v(313)+((2d0/3d0)*v(194)+v(221)+v(239))*v(85))
v(799)=props(1)*(v(305)*v(333)+(v(203)+v(221)+(2d0/3d0)*v(230))*v(85))
v(798)=props(1)*(v(305)*v(323)+(v(203)+(2d0/3d0)*v(212)+v(239))*v(85))
v(797)=props(1)*(v(306)*v(313)+((2d0/3d0)*v(195)+v(222)+v(240))*v(85))
v(796)=props(1)*(v(306)*v(333)+(v(204)+v(222)+(2d0/3d0)*v(231))*v(85))
v(795)=props(1)*(v(306)*v(323)+(v(204)+(2d0/3d0)*v(213)+v(240))*v(85))
v(794)=props(1)*(v(307)*v(313)+((2d0/3d0)*v(196)+v(223)+v(241))*v(85))
v(793)=props(1)*(v(307)*v(333)+(v(205)+v(223)+(2d0/3d0)*v(232))*v(85))
v(792)=props(1)*(v(307)*v(323)+(v(205)+(2d0/3d0)*v(214)+v(241))*v(85))
v(791)=props(1)*(v(308)*v(313)+((2d0/3d0)*v(197)+v(224)+v(242))*v(85))
v(790)=props(1)*(v(308)*v(333)+(v(206)+v(224)+(2d0/3d0)*v(233))*v(85))
v(789)=props(1)*(v(308)*v(323)+(v(206)+(2d0/3d0)*v(215)+v(242))*v(85))
v(788)=props(1)*(v(309)*v(313)+((2d0/3d0)*v(198)+v(225)+v(243))*v(85))
v(787)=props(1)*(v(309)*v(333)+(v(207)+v(225)+(2d0/3d0)*v(234))*v(85))
v(786)=props(1)*(v(309)*v(323)+(v(207)+(2d0/3d0)*v(216)+v(243))*v(85))
v(785)=props(1)*(v(310)*v(313)+((2d0/3d0)*v(199)+v(226)+v(244))*v(85))
v(784)=props(1)*(v(310)*v(333)+(v(208)+v(226)+(2d0/3d0)*v(235))*v(85))
v(783)=props(1)*(v(310)*v(323)+(v(208)+(2d0/3d0)*v(217)+v(244))*v(85))
v(782)=props(1)*(v(311)*v(313)+((2d0/3d0)*v(200)+v(227)+v(245))*v(85))
v(781)=props(1)*(v(311)*v(333)+(v(209)+v(227)+(2d0/3d0)*v(236))*v(85))
v(780)=props(1)*(v(311)*v(323)+(v(209)+(2d0/3d0)*v(218)+v(245))*v(85))
v(779)=props(1)*(v(312)*v(313)+((2d0/3d0)*v(201)+v(228)+v(246))*v(85))
v(778)=props(1)*(v(312)*v(333)+(v(210)+v(228)+(2d0/3d0)*v(237))*v(85))
v(777)=props(1)*(v(312)*v(323)+(v(210)+(2d0/3d0)*v(219)+v(246))*v(85))
v(369)=props(1)*(v(312)*v(81)+v(273)*v(85))
v(368)=props(1)*(v(311)*v(81)+v(272)*v(85))
v(367)=props(1)*(v(310)*v(81)+v(271)*v(85))
v(366)=props(1)*(v(309)*v(81)+v(270)*v(85))
v(365)=props(1)*(v(308)*v(81)+v(269)*v(85))
v(364)=props(1)*(v(307)*v(81)+v(268)*v(85))
v(363)=props(1)*(v(306)*v(81)+v(267)*v(85))
v(362)=props(1)*(v(305)*v(81)+v(266)*v(85))
v(361)=props(1)*(v(303)*v(81)+v(265)*v(85))
v(360)=props(1)*(v(312)*v(80)+v(264)*v(85))
v(359)=props(1)*(v(311)*v(80)+v(263)*v(85))
v(358)=props(1)*(v(310)*v(80)+v(262)*v(85))
v(357)=props(1)*(v(309)*v(80)+v(261)*v(85))
v(356)=props(1)*(v(308)*v(80)+v(260)*v(85))
v(355)=props(1)*(v(307)*v(80)+v(259)*v(85))
v(354)=props(1)*(v(306)*v(80)+v(258)*v(85))
v(353)=props(1)*(v(305)*v(80)+v(257)*v(85))
v(352)=props(1)*(v(303)*v(80)+v(256)*v(85))
v(351)=props(1)*(v(312)*v(79)+v(255)*v(85))
v(350)=props(1)*(v(311)*v(79)+v(254)*v(85))
v(349)=props(1)*(v(310)*v(79)+v(253)*v(85))
v(348)=props(1)*(v(309)*v(79)+v(252)*v(85))
v(347)=props(1)*(v(308)*v(79)+v(251)*v(85))
v(346)=props(1)*(v(307)*v(79)+v(250)*v(85))
v(345)=props(1)*(v(306)*v(79)+v(249)*v(85))
v(344)=props(1)*(v(305)*v(79)+v(248)*v(85))
v(343)=props(1)*(v(303)*v(79)+v(247)*v(85))
v(94)=1d0/(v(736)*x(13)+v(734)*x(16)+v(737)*x(19))**2
v(100)=1d0/v(94)**0.3333333333333333d0
v(804)=props(7)*v(100)
v(805)=v(804)*v(94)
v(95)=((v(734)*v(734))+(v(736)*v(736))+(v(737)*v(737)))*v(94)
v(99)=-v(95)/3d0
v(96)=((v(735)*v(735))+(v(740)*v(740))+(v(741)*v(741)))*v(94)
v(102)=-v(96)/3d0
v(97)=((v(738)*v(738))+(v(739)*v(739))+(v(742)*v(742)))*v(94)
v(103)=-v(97)/3d0
v(98)=v(804)*(v(102)+(2d0/3d0)*v(97)+v(99))
v(101)=v(804)*(v(103)+(2d0/3d0)*v(96)+v(99))
v(104)=v(804)*(v(102)+v(103)+(2d0/3d0)*v(95))
v(112)=(v(735)*v(738)+v(739)*v(740)+v(741)*v(742))*v(805)
v(117)=(v(735)*v(737)+v(736)*v(740)+v(734)*v(741))*v(805)
v(119)=(v(737)*v(738)+v(736)*v(739)+v(734)*v(742))*v(805)
v(120)=v(313)*v(806)+v(86)-v(98)
v(121)=-v(101)+v(323)*v(806)+v(86)
v(122)=-v(104)+v(333)*v(806)+v(86)
v(123)=-v(112)+v(79)*v(806)
v(819)=4d0*v(123)
v(124)=-v(117)+v(80)*v(806)
v(820)=4d0*v(124)
v(125)=-v(119)+v(806)*v(81)
v(821)=4d0*v(125)
IF(dabs(props(5)).lt.0.1d-11) THEN
 v(807)=v(808)
 v(127)=v(807)
ELSE
 v(809)=1d0/props(5)
 v(127)=v(809)*(1d0-dexp(-(props(5)*v(808))))
ENDIF
v(135)=1d0/(props(3)+v(127))
v(810)=0.15d1*v(135)
v(495)=v(369)*v(810)
v(504)=-(v(495)*v(732))
v(494)=v(368)*v(810)
v(503)=-(v(494)*v(732))
v(493)=v(367)*v(810)
v(502)=-(v(493)*v(732))
v(492)=v(366)*v(810)
v(501)=-(v(492)*v(732))
v(491)=v(365)*v(810)
v(500)=-(v(491)*v(732))
v(490)=v(364)*v(810)
v(499)=-(v(490)*v(732))
v(489)=v(363)*v(810)
v(498)=-(v(489)*v(732))
v(488)=v(362)*v(810)
v(497)=-(v(488)*v(732))
v(487)=v(361)*v(810)
v(496)=-(v(487)*v(732))
v(477)=v(360)*v(810)
v(486)=-(v(477)*v(732))
v(476)=v(359)*v(810)
v(485)=-(v(476)*v(732))
v(475)=v(358)*v(810)
v(484)=-(v(475)*v(732))
v(474)=v(357)*v(810)
v(483)=-(v(474)*v(732))
v(473)=v(356)*v(810)
v(482)=-(v(473)*v(732))
v(472)=v(355)*v(810)
v(481)=-(v(472)*v(732))
v(471)=v(354)*v(810)
v(480)=-(v(471)*v(732))
v(470)=v(353)*v(810)
v(479)=-(v(470)*v(732))
v(469)=v(352)*v(810)
v(478)=-(v(469)*v(732))
v(459)=v(351)*v(810)
v(468)=-(v(459)*v(732))
v(458)=v(350)*v(810)
v(467)=-(v(458)*v(732))
v(457)=v(349)*v(810)
v(466)=-(v(457)*v(732))
v(456)=v(348)*v(810)
v(465)=-(v(456)*v(732))
v(455)=v(347)*v(810)
v(464)=-(v(455)*v(732))
v(454)=v(346)*v(810)
v(463)=-(v(454)*v(732))
v(453)=v(345)*v(810)
v(462)=-(v(453)*v(732))
v(452)=v(344)*v(810)
v(461)=-(v(452)*v(732))
v(451)=v(343)*v(810)
v(460)=-(v(451)*v(732))
v(450)=v(778)*v(810)
v(449)=v(781)*v(810)
v(448)=v(784)*v(810)
v(447)=v(787)*v(810)
v(446)=v(790)*v(810)
v(445)=v(793)*v(810)
v(444)=v(796)*v(810)
v(443)=v(799)*v(810)
v(442)=v(802)*v(810)
v(441)=v(777)*v(810)
v(440)=v(780)*v(810)
v(439)=v(783)*v(810)
v(438)=v(786)*v(810)
v(437)=v(789)*v(810)
v(436)=v(792)*v(810)
v(435)=v(795)*v(810)
v(434)=v(798)*v(810)
v(433)=v(801)*v(810)
v(432)=v(779)*v(810)
v(431)=v(782)*v(810)
v(430)=v(785)*v(810)
v(429)=v(788)*v(810)
v(428)=v(791)*v(810)
v(427)=v(794)*v(810)
v(426)=v(797)*v(810)
v(425)=v(800)*v(810)
v(424)=v(803)*v(810)
v(128)=-v(121)/3d0
v(129)=-v(122)/3d0
v(132)=(2d0/3d0)*v(120)+v(128)+v(129)
v(824)=2d0*v(132)
v(130)=-v(120)/3d0
v(137)=(2d0/3d0)*v(122)+v(128)+v(130)
v(823)=2d0*v(137)
v(134)=(2d0/3d0)*v(121)+v(129)+v(130)
v(822)=2d0*v(134)
v(507)=1d0/sqrt(0.15d1*(2d0*v(123)**2+2d0*v(124)**2+2d0*v(125)**2+v(132)**2+v(134)**2+v(137)**2))
v(818)=0.75d0*v(507)
v(133)=v(132)*v(810)
v(136)=v(134)*v(810)
v(138)=v(137)*v(810)
v(139)=v(123)*v(810)
v(140)=v(124)*v(810)
v(141)=v(125)*v(810)
IF(sqrt(0.15d1*(2d0*v(112)**2+2d0*v(117)**2+2d0*v(119)**2+((2d0/3d0)*v(101)-v(104)/3d0-v(98)/3d0)**2+(-v(101)/3d0+(2d0&
&/3d0)*v(104)-v(98)/3d0)**2+(-v(101)/3d0-v(104)/3d0+(2d0/3d0)*v(98))**2)).lt.0.1d-11) THEN
 v(516)=0d0
 v(517)=0d0
 v(518)=0d0
 v(519)=0d0
 v(520)=0d0
 v(521)=0d0
 v(522)=0d0
 v(523)=0d0
 v(524)=0d0
 v(525)=0d0
 v(526)=0d0
 v(527)=0d0
 v(528)=0d0
 v(529)=0d0
 v(530)=0d0
 v(531)=0d0
 v(532)=0d0
 v(533)=0d0
 v(534)=0d0
 v(535)=0d0
 v(536)=0d0
 v(537)=0d0
 v(538)=0d0
 v(539)=0d0
 v(540)=0d0
 v(541)=0d0
 v(542)=0d0
 v(543)=0d0
 v(544)=0d0
 v(545)=0d0
 v(546)=0d0
 v(547)=0d0
 v(548)=0d0
 v(549)=0d0
 v(550)=0d0
 v(551)=0d0
 v(552)=0d0
 v(553)=0d0
 v(554)=0d0
 v(555)=0d0
 v(556)=0d0
 v(557)=0d0
 v(558)=0d0
 v(559)=0d0
 v(560)=0d0
 v(561)=0d0
 v(562)=0d0
 v(563)=0d0
 v(564)=0d0
 v(565)=0d0
 v(566)=0d0
 v(567)=0d0
 v(568)=0d0
 v(569)=0d0
 v(570)=0d0
 v(571)=0d0
 v(572)=0d0
 v(573)=0d0
 v(574)=0d0
 v(575)=0d0
 v(576)=0d0
 v(577)=0d0
 v(578)=0d0
 v(579)=0d0
 v(580)=0d0
 v(581)=0d0
 v(582)=0d0
 v(583)=0d0
 v(584)=0d0
 v(585)=0d0
 v(586)=0d0
 v(587)=0d0
 v(588)=0d0
 v(589)=0d0
 v(590)=0d0
 v(591)=0d0
 v(592)=0d0
 v(593)=0d0
 v(594)=0d0
 v(595)=0d0
 v(596)=0d0
ELSE
 v(815)=props(8)**(1d0+props(9))
 v(814)=1d0-props(6)
 v(813)=2d0*v(119)
 v(812)=2d0*v(117)
 v(811)=2d0*v(112)
 v(605)=v(101)*v(441)+v(104)*v(450)+v(459)*v(811)+v(477)*v(812)+v(495)*v(813)+v(432)*v(98)
 v(604)=v(101)*v(440)+v(104)*v(449)+v(458)*v(811)+v(476)*v(812)+v(494)*v(813)+v(431)*v(98)
 v(603)=v(101)*v(439)+v(104)*v(448)+v(457)*v(811)+v(475)*v(812)+v(493)*v(813)+v(430)*v(98)
 v(602)=v(101)*v(438)+v(104)*v(447)+v(456)*v(811)+v(474)*v(812)+v(492)*v(813)+v(429)*v(98)
 v(601)=v(101)*v(437)+v(104)*v(446)+v(455)*v(811)+v(473)*v(812)+v(491)*v(813)+v(428)*v(98)
 v(600)=v(101)*v(436)+v(104)*v(445)+v(454)*v(811)+v(472)*v(812)+v(490)*v(813)+v(427)*v(98)
 v(599)=v(101)*v(435)+v(104)*v(444)+v(453)*v(811)+v(471)*v(812)+v(489)*v(813)+v(426)*v(98)
 v(598)=v(101)*v(434)+v(104)*v(443)+v(452)*v(811)+v(470)*v(812)+v(488)*v(813)+v(425)*v(98)
 v(597)=v(101)*v(433)+v(104)*v(442)+v(451)*v(811)+v(469)*v(812)+v(487)*v(813)+v(424)*v(98)
 v(175)=sqrt(0.15d1*(v(101)**2+v(104)**2+2d0*v(112)**2+2d0*v(117)**2+2d0*v(119)**2+v(98)**2))
 v(606)=(v(175)**props(9)*v(815))/2d0
 v(178)=(0.15d1*props(6))/v(175)
 v(622)=v(119)*v(178)+v(141)*v(814)
 v(621)=v(117)*v(178)+v(140)*v(814)
 v(620)=v(112)*v(178)+v(139)*v(814)
 v(619)=v(104)*v(178)+v(138)*v(814)
 v(618)=v(101)*v(178)+v(136)*v(814)
 v(617)=v(133)*v(814)+v(178)*v(98)
 v(176)=v(101)*v(136)+v(104)*v(138)+v(139)*v(811)+v(140)*v(812)+v(141)*v(813)+v(133)*v(98)
 v(608)=dsign(1.d0,v(176))
 v(816)=v(606)*(1d0+v(608))
 v(616)=v(605)*v(816)
 v(615)=v(604)*v(816)
 v(614)=v(603)*v(816)
 v(613)=v(602)*v(816)
 v(612)=v(601)*v(816)
 v(611)=v(600)*v(816)
 v(610)=v(599)*v(816)
 v(609)=v(598)*v(816)
 v(607)=v(597)*v(816)
 v(177)=v(606)*(v(176)+dabs(v(176)))
 v(817)=v(177)*v(814)
 v(516)=v(607)*v(617)+v(424)*v(817)
 v(517)=v(609)*v(617)+v(425)*v(817)
 v(518)=v(610)*v(617)+v(426)*v(817)
 v(519)=v(611)*v(617)+v(427)*v(817)
 v(520)=v(612)*v(617)+v(428)*v(817)
 v(521)=v(613)*v(617)+v(429)*v(817)
 v(522)=v(614)*v(617)+v(430)*v(817)
 v(523)=v(615)*v(617)+v(431)*v(817)
 v(524)=v(616)*v(617)+v(432)*v(817)
 v(525)=v(607)*v(618)+v(433)*v(817)
 v(526)=v(609)*v(618)+v(434)*v(817)
 v(527)=v(610)*v(618)+v(435)*v(817)
 v(528)=v(611)*v(618)+v(436)*v(817)
 v(529)=v(612)*v(618)+v(437)*v(817)
 v(530)=v(613)*v(618)+v(438)*v(817)
 v(531)=v(614)*v(618)+v(439)*v(817)
 v(532)=v(615)*v(618)+v(440)*v(817)
 v(533)=v(616)*v(618)+v(441)*v(817)
 v(534)=v(607)*v(619)+v(442)*v(817)
 v(535)=v(609)*v(619)+v(443)*v(817)
 v(536)=v(610)*v(619)+v(444)*v(817)
 v(537)=v(611)*v(619)+v(445)*v(817)
 v(538)=v(612)*v(619)+v(446)*v(817)
 v(539)=v(613)*v(619)+v(447)*v(817)
 v(540)=v(614)*v(619)+v(448)*v(817)
 v(541)=v(615)*v(619)+v(449)*v(817)
 v(542)=v(616)*v(619)+v(450)*v(817)
 v(543)=v(607)*v(620)+v(451)*v(817)
 v(544)=v(609)*v(620)+v(452)*v(817)
 v(545)=v(610)*v(620)+v(453)*v(817)
 v(546)=v(611)*v(620)+v(454)*v(817)
 v(547)=v(612)*v(620)+v(455)*v(817)
 v(548)=v(613)*v(620)+v(456)*v(817)
 v(549)=v(614)*v(620)+v(457)*v(817)
 v(550)=v(615)*v(620)+v(458)*v(817)
 v(551)=v(616)*v(620)+v(459)*v(817)
 v(552)=v(607)*v(621)+v(469)*v(817)
 v(553)=v(609)*v(621)+v(470)*v(817)
 v(554)=v(610)*v(621)+v(471)*v(817)
 v(555)=v(611)*v(621)+v(472)*v(817)
 v(556)=v(612)*v(621)+v(473)*v(817)
 v(557)=v(613)*v(621)+v(474)*v(817)
 v(558)=v(614)*v(621)+v(475)*v(817)
 v(559)=v(615)*v(621)+v(476)*v(817)
 v(560)=v(616)*v(621)+v(477)*v(817)
 v(561)=v(607)*v(622)+v(487)*v(817)
 v(562)=v(609)*v(622)+v(488)*v(817)
 v(563)=v(610)*v(622)+v(489)*v(817)
 v(564)=v(611)*v(622)+v(490)*v(817)
 v(565)=v(612)*v(622)+v(491)*v(817)
 v(566)=v(613)*v(622)+v(492)*v(817)
 v(567)=v(614)*v(622)+v(493)*v(817)
 v(568)=v(615)*v(622)+v(494)*v(817)
 v(569)=v(616)*v(622)+v(495)*v(817)
 v(570)=v(561)
 v(571)=v(562)
 v(572)=v(563)
 v(573)=v(564)
 v(574)=v(565)
 v(575)=v(566)
 v(576)=v(567)
 v(577)=v(568)
 v(578)=v(569)
 v(579)=v(543)
 v(580)=v(544)
 v(581)=v(545)
 v(582)=v(546)
 v(583)=v(547)
 v(584)=v(548)
 v(585)=v(549)
 v(586)=v(550)
 v(587)=v(551)
 v(588)=v(552)
 v(589)=v(553)
 v(590)=v(554)
 v(591)=v(555)
 v(592)=v(556)
 v(593)=v(557)
 v(594)=v(558)
 v(595)=v(559)
 v(596)=v(560)
ENDIF
dRdF(1,1)=v(818)*(v(343)*v(819)+v(352)*v(820)+v(361)*v(821)+v(801)*v(822)+v(802)*v(823)+v(803)*v(824))
dRdF(1,2)=v(818)*(v(344)*v(819)+v(353)*v(820)+v(362)*v(821)+v(798)*v(822)+v(799)*v(823)+v(800)*v(824))
dRdF(1,3)=v(818)*(v(345)*v(819)+v(354)*v(820)+v(363)*v(821)+v(795)*v(822)+v(796)*v(823)+v(797)*v(824))
dRdF(1,4)=v(818)*(v(346)*v(819)+v(355)*v(820)+v(364)*v(821)+v(792)*v(822)+v(793)*v(823)+v(794)*v(824))
dRdF(1,5)=v(818)*(v(347)*v(819)+v(356)*v(820)+v(365)*v(821)+v(789)*v(822)+v(790)*v(823)+v(791)*v(824))
dRdF(1,6)=v(818)*(v(348)*v(819)+v(357)*v(820)+v(366)*v(821)+v(786)*v(822)+v(787)*v(823)+v(788)*v(824))
dRdF(1,7)=v(818)*(v(349)*v(819)+v(358)*v(820)+v(367)*v(821)+v(783)*v(822)+v(784)*v(823)+v(785)*v(824))
dRdF(1,8)=v(818)*(v(350)*v(819)+v(359)*v(820)+v(368)*v(821)+v(780)*v(822)+v(781)*v(823)+v(782)*v(824))
dRdF(1,9)=v(818)*(v(351)*v(819)+v(360)*v(820)+v(369)*v(821)+v(777)*v(822)+v(778)*v(823)+v(779)*v(824))
dRdF(2,1)=-(v(424)*v(732))
dRdF(2,2)=-(v(425)*v(732))
dRdF(2,3)=-(v(426)*v(732))
dRdF(2,4)=-(v(427)*v(732))
dRdF(2,5)=-(v(428)*v(732))
dRdF(2,6)=-(v(429)*v(732))
dRdF(2,7)=-(v(430)*v(732))
dRdF(2,8)=-(v(431)*v(732))
dRdF(2,9)=-(v(432)*v(732))
dRdF(3,1)=-(v(433)*v(732))
dRdF(3,2)=-(v(434)*v(732))
dRdF(3,3)=-(v(435)*v(732))
dRdF(3,4)=-(v(436)*v(732))
dRdF(3,5)=-(v(437)*v(732))
dRdF(3,6)=-(v(438)*v(732))
dRdF(3,7)=-(v(439)*v(732))
dRdF(3,8)=-(v(440)*v(732))
dRdF(3,9)=-(v(441)*v(732))
dRdF(4,1)=-(v(442)*v(732))
dRdF(4,2)=-(v(443)*v(732))
dRdF(4,3)=-(v(444)*v(732))
dRdF(4,4)=-(v(445)*v(732))
dRdF(4,5)=-(v(446)*v(732))
dRdF(4,6)=-(v(447)*v(732))
dRdF(4,7)=-(v(448)*v(732))
dRdF(4,8)=-(v(449)*v(732))
dRdF(4,9)=-(v(450)*v(732))
dRdF(5,1)=v(460)
dRdF(5,2)=v(461)
dRdF(5,3)=v(462)
dRdF(5,4)=v(463)
dRdF(5,5)=v(464)
dRdF(5,6)=v(465)
dRdF(5,7)=v(466)
dRdF(5,8)=v(467)
dRdF(5,9)=v(468)
dRdF(6,1)=v(478)
dRdF(6,2)=v(479)
dRdF(6,3)=v(480)
dRdF(6,4)=v(481)
dRdF(6,5)=v(482)
dRdF(6,6)=v(483)
dRdF(6,7)=v(484)
dRdF(6,8)=v(485)
dRdF(6,9)=v(486)
dRdF(7,1)=v(496)
dRdF(7,2)=v(497)
dRdF(7,3)=v(498)
dRdF(7,4)=v(499)
dRdF(7,5)=v(500)
dRdF(7,6)=v(501)
dRdF(7,7)=v(502)
dRdF(7,8)=v(503)
dRdF(7,9)=v(504)
dRdF(8,1)=v(496)
dRdF(8,2)=v(497)
dRdF(8,3)=v(498)
dRdF(8,4)=v(499)
dRdF(8,5)=v(500)
dRdF(8,6)=v(501)
dRdF(8,7)=v(502)
dRdF(8,8)=v(503)
dRdF(8,9)=v(504)
dRdF(9,1)=v(460)
dRdF(9,2)=v(461)
dRdF(9,3)=v(462)
dRdF(9,4)=v(463)
dRdF(9,5)=v(464)
dRdF(9,6)=v(465)
dRdF(9,7)=v(466)
dRdF(9,8)=v(467)
dRdF(9,9)=v(468)
dRdF(10,1)=v(478)
dRdF(10,2)=v(479)
dRdF(10,3)=v(480)
dRdF(10,4)=v(481)
dRdF(10,5)=v(482)
dRdF(10,6)=v(483)
dRdF(10,7)=v(484)
dRdF(10,8)=v(485)
dRdF(10,9)=v(486)
dRdF(11,1)=-((-v(424)+v(516))*v(732))
dRdF(11,2)=-((-v(425)+v(517))*v(732))
dRdF(11,3)=-((-v(426)+v(518))*v(732))
dRdF(11,4)=-((-v(427)+v(519))*v(732))
dRdF(11,5)=-((-v(428)+v(520))*v(732))
dRdF(11,6)=-((-v(429)+v(521))*v(732))
dRdF(11,7)=-((-v(430)+v(522))*v(732))
dRdF(11,8)=-((-v(431)+v(523))*v(732))
dRdF(11,9)=-((-v(432)+v(524))*v(732))
dRdF(12,1)=-((-v(433)+v(525))*v(732))
dRdF(12,2)=-((-v(434)+v(526))*v(732))
dRdF(12,3)=-((-v(435)+v(527))*v(732))
dRdF(12,4)=-((-v(436)+v(528))*v(732))
dRdF(12,5)=-((-v(437)+v(529))*v(732))
dRdF(12,6)=-((-v(438)+v(530))*v(732))
dRdF(12,7)=-((-v(439)+v(531))*v(732))
dRdF(12,8)=-((-v(440)+v(532))*v(732))
dRdF(12,9)=-((-v(441)+v(533))*v(732))
dRdF(13,1)=-((-v(442)+v(534))*v(732))
dRdF(13,2)=-((-v(443)+v(535))*v(732))
dRdF(13,3)=-((-v(444)+v(536))*v(732))
dRdF(13,4)=-((-v(445)+v(537))*v(732))
dRdF(13,5)=-((-v(446)+v(538))*v(732))
dRdF(13,6)=-((-v(447)+v(539))*v(732))
dRdF(13,7)=-((-v(448)+v(540))*v(732))
dRdF(13,8)=-((-v(449)+v(541))*v(732))
dRdF(13,9)=-((-v(450)+v(542))*v(732))
dRdF(14,1)=-((-v(451)+v(543))*v(732))
dRdF(14,2)=-((-v(452)+v(544))*v(732))
dRdF(14,3)=-((-v(453)+v(545))*v(732))
dRdF(14,4)=-((-v(454)+v(546))*v(732))
dRdF(14,5)=-((-v(455)+v(547))*v(732))
dRdF(14,6)=-((-v(456)+v(548))*v(732))
dRdF(14,7)=-((-v(457)+v(549))*v(732))
dRdF(14,8)=-((-v(458)+v(550))*v(732))
dRdF(14,9)=-((-v(459)+v(551))*v(732))
dRdF(15,1)=-((-v(469)+v(552))*v(732))
dRdF(15,2)=-((-v(470)+v(553))*v(732))
dRdF(15,3)=-((-v(471)+v(554))*v(732))
dRdF(15,4)=-((-v(472)+v(555))*v(732))
dRdF(15,5)=-((-v(473)+v(556))*v(732))
dRdF(15,6)=-((-v(474)+v(557))*v(732))
dRdF(15,7)=-((-v(475)+v(558))*v(732))
dRdF(15,8)=-((-v(476)+v(559))*v(732))
dRdF(15,9)=-((-v(477)+v(560))*v(732))
dRdF(16,1)=-((-v(487)+v(561))*v(732))
dRdF(16,2)=-((-v(488)+v(562))*v(732))
dRdF(16,3)=-((-v(489)+v(563))*v(732))
dRdF(16,4)=-((-v(490)+v(564))*v(732))
dRdF(16,5)=-((-v(491)+v(565))*v(732))
dRdF(16,6)=-((-v(492)+v(566))*v(732))
dRdF(16,7)=-((-v(493)+v(567))*v(732))
dRdF(16,8)=-((-v(494)+v(568))*v(732))
dRdF(16,9)=-((-v(495)+v(569))*v(732))
dRdF(17,1)=-((-v(487)+v(570))*v(732))
dRdF(17,2)=-((-v(488)+v(571))*v(732))
dRdF(17,3)=-((-v(489)+v(572))*v(732))
dRdF(17,4)=-((-v(490)+v(573))*v(732))
dRdF(17,5)=-((-v(491)+v(574))*v(732))
dRdF(17,6)=-((-v(492)+v(575))*v(732))
dRdF(17,7)=-((-v(493)+v(576))*v(732))
dRdF(17,8)=-((-v(494)+v(577))*v(732))
dRdF(17,9)=-((-v(495)+v(578))*v(732))
dRdF(18,1)=-((-v(451)+v(579))*v(732))
dRdF(18,2)=-((-v(452)+v(580))*v(732))
dRdF(18,3)=-((-v(453)+v(581))*v(732))
dRdF(18,4)=-((-v(454)+v(582))*v(732))
dRdF(18,5)=-((-v(455)+v(583))*v(732))
dRdF(18,6)=-((-v(456)+v(584))*v(732))
dRdF(18,7)=-((-v(457)+v(585))*v(732))
dRdF(18,8)=-((-v(458)+v(586))*v(732))
dRdF(18,9)=-((-v(459)+v(587))*v(732))
dRdF(19,1)=-((-v(469)+v(588))*v(732))
dRdF(19,2)=-((-v(470)+v(589))*v(732))
dRdF(19,3)=-((-v(471)+v(590))*v(732))
dRdF(19,4)=-((-v(472)+v(591))*v(732))
dRdF(19,5)=-((-v(473)+v(592))*v(732))
dRdF(19,6)=-((-v(474)+v(593))*v(732))
dRdF(19,7)=-((-v(475)+v(594))*v(732))
dRdF(19,8)=-((-v(476)+v(595))*v(732))
dRdF(19,9)=-((-v(477)+v(596))*v(732))
END

max_stars_repo_name: kengwit/MaterialModels

max_issues_repo_name: kengwit/MaterialModels

max_forks_repo_name: kengwit/MaterialModels


program multi_dimensional_bruteforce_mc

      !Author : Anantha Rao
      ! Reg no : 20181044
      ! 6D (x1,x2,x3,x4,x5,x6) Integral
      implicit none
      integer:: n,i,j
      real*8 :: x(6), y, func, int_mc, var, sigma, p(2)
      ! p : Dummy argument for random number generator (2 values)
      real*8 :: length, volume, sqrt2, gauss_dev

      length=5.0d0 ! use one side of each dimension (0,5)
      volume=acos(-1.0d0)**3  ! volume of 6D space
      sqrt2=1.0d0/sqrt(2.0d0)

      open(unit=1, file="20181044_multi_impsampling_mc.dat")
      n=1

      ! Initiate variables

print *,"Welcome to this program.&
& This program solves mutidimensional integral using Monte-Carlo rule with sampling"

    7 int_mc=0.0d0
      var=0.0d0
      sigma=0.0d0
      write(*,10) n
      10 format("Computig for n=",i10)

      do i=1,n
        do j=1,6
            call random_number(p)
            x(j) = gauss_dev(p)*sqrt2
        end do
        int_mc=int_mc+func(x)
        sigma=sigma+func(x)*func(x)
      end do

      int_mc = int_mc/real(n)
      sigma=sigma/real(n)
      var=sigma-int_mc*int_mc

      int_mc = volume*int_mc
      sigma=volume*sqrt(var/real(n))

      write(1,*) n, "", int_mc, "", sigma

      ! begin automation
      n=n*10
      if (n .lt. 1000000000 ) goto 7
      print*,"Done! Output stored in file 20181044_multi_impsampling_mc.dat"
      ! end automation

end program

real*8 function func(x)
      implicit none
      real*8::x(6), xy, a
      a=0.5d0

      xy= (x(1)-x(4))**2 + (x(2)-x(5))**2 + (x(3)-x(6))**2
      func=exp(- a*xy)

end function

real*8 function gauss_dev(x)
    implicit none
    real*8:: fact, sqr, p, x1, x2, x(2)

    7 call random_number(p)
    x1 = 2.0d0*p - 1.0d0
    call random_number(p)
    x2 = 2.0d0*p - 1.0d0
    sqr = x1*x1 + x2*x2

    if (sqr .ge. 1.0d0 .or. sqr .eq. 0.0d0) goto 7

    fact=sqrt(-2.0d0*log(sqr)/sqr)
    gauss_dev=x2*fact
    end function


! { dg-do run }
use omp_lib
  call test_parallel
  call test_do
  call test_sections
  call test_single

contains
  subroutine test_parallel
    integer :: a, b, c, e, f, g, i, j
    integer, dimension (20) :: d
    logical :: h
    a = 6
    b = 8
    c = 11
    d(:) = -1
    e = 13
    f = 24
    g = 27
    h = .false.
    i = 1
    j = 16
!$omp para&
!$omp&llel &
!$omp if (a .eq. 6) private (b, c) shared (d) private (e) &
  !$omp firstprivate(f) num_threads (a - 1) first&
!$ompprivate(g)default (shared) reduction (.or. : h) &
!$omp reduction(*:i)
    if (i .ne. 1) h = .true.
    i = 2
    if (f .ne. 24) h = .true.
    if (g .ne. 27) h = .true.
    e = 7
    b = omp_get_thread_num ()
    if (b .eq. 0) j = 24
    f = b
    g = f
    c = omp_get_num_threads ()
    if (c .gt. a - 1 .or. c .le. 0) h = .true.
    if (b .ge. c) h = .true.
    d(b + 1) = c
    if (f .ne. g .or. f .ne. b) h = .true.
!$omp endparallel
    if (h) call abort
    if (a .ne. 6) call abort
    if (j .ne. 24) call abort
    if (d(1) .eq. -1) call abort
    e = 1
    do g = 1, d(1)
      if (d(g) .ne. d(1)) call abort
      e = e * 2
    end do
    if (e .ne. i) call abort
  end subroutine test_parallel

  subroutine test_do_orphan
    integer :: k, l
!$omp parallel do private (l)
    do 600 k = 1, 16, 2
600   l = k
  end subroutine test_do_orphan

  subroutine test_do
    integer :: i, j, k, l, n
    integer, dimension (64) :: d
    logical :: m

    j = 16
    d(:) = -1
    m = .true.
    n = 24
!$omp parallel num_threads (4) shared (i, k, d) private (l) &
!$omp&reduction (.and. : m)
    if (omp_get_thread_num () .eq. 0) then
      k = omp_get_num_threads ()
    end if
    call test_do_orphan
!$omp do schedule (static) firstprivate (n)
    do 200 i = 1, j
      if (i .eq. 1 .and. n .ne. 24) call abort
      n = i
200   d(n) = omp_get_thread_num ()
!$omp enddo nowait

!$omp do lastprivate (i) schedule (static, 5)
    do 201 i = j + 1, 2 * j
201   d(i) = omp_get_thread_num () + 1024
    ! Implied omp end do here

    if (i .ne. 33) m = .false.

!$omp do private (j) schedule (dynamic)
    do i = 33, 48
      d(i) = omp_get_thread_num () + 2048
    end do
!$omp end do nowait

!$omp do schedule (runtime)
    do i = 49, 4 * j
      d(i) = omp_get_thread_num () + 4096
    end do
    ! Implied omp end do here
!$omp end parallel
    if (.not. m) call abort

    j = 0
    do i = 1, 64
      if (d(i) .lt. j .or. d(i) .ge. j + k) call abort
      if (i .eq. 16) j = 1024
      if (i .eq. 32) j = 2048
      if (i .eq. 48) j = 4096
    end do
  end subroutine test_do

  subroutine test_sections
    integer :: i, j, k, l, m, n
    i = 9
    j = 10
    k = 11
    l = 0
    m = 0
    n = 30
    call omp_set_dynamic (.false.)
    call omp_set_num_threads (4)
!$omp parallel num_threads (4)
!$omp sections private (i) firstprivate (j, k) lastprivate (j) &
!$omp& reduction (+ : l, m)
!$omp section
    i = 24
    if (j .ne. 10 .or. k .ne. 11 .or. m .ne. 0) l = 1
    m = m + 4
!$omp section
    i = 25
    if (j .ne. 10 .or. k .ne. 11) l = 1
    m = m + 6
!$omp section
    i = 26
    if (j .ne. 10 .or. k .ne. 11) l = 1
    m = m + 8
!$omp section
    i = 27
    if (j .ne. 10 .or. k .ne. 11) l = 1
    m = m + 10
    j = 271
!$omp end sections nowait
!$omp sections lastprivate (n)
!$omp section
    n = 6
!$omp section
    n = 7
!$omp endsections
!$omp end parallel
    if (j .ne. 271 .or. l .ne. 0) call abort
    if (m .ne. 4 + 6 + 8 + 10) call abort
    if (n .ne. 7) call abort
  end subroutine test_sections

  subroutine test_single
    integer :: i, j, k, l
    logical :: m
    i = 200
    j = 300
    k = 400
    l = 500
    m = .false.
!$omp parallel num_threads (4), private (i, j), reduction (.or. : m)
    i = omp_get_thread_num ()
    j = omp_get_thread_num ()
!$omp single private (k)
    k = 64
!$omp end single nowait
!$omp single private (k) firstprivate (l)
    if (i .ne. omp_get_thread_num () .or. i .ne. j) then
      j = -1
    else
      j = -2
    end if
    if (l .ne. 500) j = -1
    l = 265
!$omp end single copyprivate (j)
    if (i .ne. omp_get_thread_num () .or. j .ne. -2) m = .true.
!$omp endparallel
    if (m) call abort
  end subroutine test_single
end


text:
module mod1

end module mod1

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
Im Mike Molina.

Not related to Users/JohnMolina, but they are buddies.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


PROGRAM  three
    IMPLICIT NONE
    REAL:: x, xe, beta, De, MorseVx, HarmonicVx
    CHARACTER*10:: filename
    INTEGER:: i, outunit
    xe = 0.74
    beta = 1.0
    De = 109.4 * 1.59 * (10**(-3.0))

    WRITE(filename, fmt = "(a, i0)") "output.", 1
    OPEN(UNIT = outunit, FILE = filename, FORM = "formatted")
    x = -4.2
    DO i = 1, 100
        x = x + 0.20
        MorseVx = De * (1 - EXP(-beta * (x - xe))) * (1 - EXP(-beta * (x - xe)))
        HarmonicVx = De * beta * beta * (x-xe) * (x-xe)
        WRITE(outunit, * ) x, MorseVx, HarmonicVx
    ENDDO
    CLOSE(outunit)
END


text:
c-----------------------------------------------------------------------
      program trans

c     this program will take a nek all-ascii rea file and
c     extract geometry, bcs, and curve data to a new binary re2 file, plus
c     an ascii rea file for just the parameters

      character*4 ans
      real*4      buf(10)
      real*8      buf2(30)
      real*8      rgroup,rcurve,rbc

      character*80 file,fout,fbout,string
      character*1  file1(80),fout1(80),fbout1(80),string1(80)
      equivalence (file1,file)
      equivalence (fout1,fout)
      equivalence (fbout1,fbout)
      equivalence (string,string1)

      parameter(lelt=MAXNEL)
      common /array/ x(8),y(8),z(8),bc(5,6,lelt),curve(6,8)
      common /arrai/ nel,ncurve
      common /arrac/ cbc
      character*3 cbc(6,lelt)
      integer e,f
      character*1 ccurve(4)
      logical ifflow,ifheat

      character*80 hdr
      real*4 test
      data   test  / 6.54321 /

      integer ibc(6,lelt)
      real*8  bc8(5)

      write(6,*) 'Input .rea name:'
      call blank(file,80)
      read(5,80) file
      len = ltrunc(file,80)
   80 format(a80)
   81 format(80a1)

      write(6,*) 'Input .rea/.re2 output name:'
      call blank(fout,80)
      read(5,80) fout
      fbout = fout
      lou = ltrunc(fout,80)

      write(6,*)
      write(6,'(A)') 'Start converting ...'
      write(6,*)

      call chcopy(file1(len+1),'.rea',4)
      call chcopy(fout1(lou+1),'.rea',4)
      call chcopy(fbout1(lou+1),'.re2' // char(0),5)

      open(unit=10, file=file)
      open(unit=11, file=fout)
      call byte_open(fbout,ierr)

      call scanout(string,'MESH DATA',9,10,11)
      call lineout(11,string,80)
      write(6,'(A,A)') 'read: ',string

      read (10,*)   nel, ndim, nelv
      write(11,11) -nel, ndim, nelv
   11 format(i12,2x,i1,2x,i12,5x,'NELT,NDIM,NELV')

      if(nel.ge.lelt) then
        write(6,*) 'Abort: number of elements too large',nel
        write(6,*) 'change MAXNEL and recompile'
        stop
      endif

      call blank(hdr,80)
      write(hdr,1) nel,ndim,nelv
    1 format('#v002',i9,i3,i9,' this is the hdr')
      call byte_write(hdr,20,ierr)   ! assumes byte_open() already issued
      call byte_write(test,1,ierr)   ! write the endian discriminator

C MESH
      igroup = 0
      do ie=1,nel
         read(10,*) ! for now skip igroup (not used at the moment)
         rgroup=igroup
         call byte_write(rgroup,2,ierr)

         if (ndim.eq.3) then
            read (10,*)   (x(k),k=1,4)
            read (10,*)   (y(k),k=1,4)
            read (10,*)   (z(k),k=1,4)
            read (10,*)   (x(k),k=5,8)
            read (10,*)   (y(k),k=5,8)
            read (10,*)   (z(k),k=5,8)

            call copy48(buf2(1) ,x,8)
            call copy48(buf2(9) ,y,8)
            call copy48(buf2(17),z,8)
            call byte_write(buf2(1) ,16,ierr)
            call byte_write(buf2(9) ,16,ierr)
            call byte_write(buf2(17),16,ierr)
         else
            read (10,*)   (x(k),k=1,4)
            read (10,*)   (y(k),k=1,4)

            call copy48(buf2(1) ,x,4)
            call copy48(buf2(5) ,y,4)
            call byte_write(buf2(1),8,ierr)
            call byte_write(buf2(5),8,ierr)
         endif
      enddo
      write(6,'(i12,1x,i12,1x,i1,A)') nel,nelv,ndim,
     &     ' nelt/nelv/ndim'

c CURVED SIDES
      read (10,80) string !  curved sides
      write(6,'(A,A)') 'read: ',string
      read(10,*) ncurve
      rcurve=ncurve
      call byte_write(rcurve,2,ierr)

      call blank(ccurve,4)
      if (ncurve.gt.0) then
         do icurve=1,ncurve
c           if (mod(icurve,10000).eq.0) write(6,*) icurve,' curve'
            if (nel.lt.1000) then
               read(10,60) f,e,(buf(k),k=1,5),ccurve(1)
            elseif (nel.lt.1 000 000) then
               read(10,61) f,e,(buf(k),k=1,5),ccurve(1)
            else
               read(10,62) f,e,(buf(k),k=1,5),ccurve(1)
            endif
            buf2(1)=e
            buf2(2)=f
            call copy48(buf2(3),buf(1),5)
            call blank(buf2(8),8)
            call chcopy(buf2(8),ccurve,1)
            call byte_write(buf2,16,ierr)
         enddo
   60    format(i3,i3 ,5g14.6,1x,a1)
   61    format(i2,i6 ,5g14.6,1x,a1)
   62    format(i2,i12,5g14.6,1x,a1)

          write(6,*) ncurve,' Number of curved sides'
      endif

C BOUNDARY CONDITIONS
      read (10,80) string ! ***** BOUNDARY CONDITIONS *****
      write(6,'(A,A)') 'read: ',string

      nface = 2*ndim
      do kpass = 1,50   ! for now, at most 50 fields
         read (10,80) string
         if (indx1(string,'BOUN',4).ne.0) then ! we might have bcs
            if (indx1(string,'NO ',3).eq.0)  then ! we have bcs, read and count
               write(6,'(A,A)') 'read: ',string
               nbc = 0
               nelb = nelv
               if(kpass.eq.2) nelb=nel     ! only ifield2 is a T mesh
               do e=1,nelb
               do f=1,nface
                  if (nel.lt.1 000 000) then
                     read(10,20) cbc(f,e),(bc(k,f,e),k=1,5)
                  else
                     read(10,21) cbc(f,e),(bc8(k),k=1,5)
                     do ii = 1,5
                        bc(ii,f,e)=bc8(ii)
                     enddo
                     ibc(f,e) = bc8(1)
                  endif
                  if (cbc(f,e).ne.'E  ') nbc = nbc+1
               enddo
               enddo
   20          format(1x,a3,6x,5g14.6)
   21          format(1x,a3,12x,5g18.11)

               write(6,*) kpass,nbc,' Number of bcs'
               rbc=nbc
               call byte_write(rbc,2,ierr)

               do e=1,nelb
               do f=1,nface
                  if (cbc(f,e).ne.'E  ') then
                     buf2(1)=e
                     buf2(2)=f
                     call copy48(buf2(3),bc(1,f,e),5)
                     call blank      (buf2(8),8)
                     call chcopy     (buf2(8),cbc(f,e),3)
                     if(nelb.ge.1000000) buf2(3)=ibc(f,e)
                     call byte_write (buf2,16,ierr)
                  endif
               enddo
               enddo

            endif
         else
            goto 51
         endif

      enddo
   51 continue

      len = ltrunc(string,80)
      write(11,81) (string1(k),k=1,len)
c     write(6 ,81) (string1(k),k=1,len)

      call scanout(string,'end',3,10,11)

      close (unit=10)
      close (unit=11)
      call byte_close (fbout,ierr)

  999 continue
      write(6,*) 'done'
      stop

 9999 continue
      write(6,*) 'could not find file "indat". abort.'
      stop
      end
c-----------------------------------------------------------------------
      subroutine exitt

      stop

      return
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYEVOL
C...Handles intertwined pT-ordered spacelike initial-state parton
C...and multiple interactions.

      SUBROUTINE PYEVOL(MODE,PT2MAX,PT2MIN)
C...Mode = -1 : Initialize first time. Determine MAX and MIN scales.
C...MODE =  0 : (Re-)initialize ISR/MI evolution.
C...Mode =  1 : Evolve event from PT2MAX to PT2MIN.

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)

C...External
      EXTERNAL PYALPS
      DOUBLE PRECISION PYALPS
C...Parameter statement for maximum size of showers.
      PARAMETER (MAXNUR=1000)
C...Commonblocks.
      include 'inc/pypart'
      include 'inc/pyjets'
      include 'inc/pydat1'
      include 'inc/pydat2'
      include 'inc/pypars'
      include 'inc/pyint1'
      include 'inc/pyint2'
      include 'inc/pyint3'
      include 'inc/pyintm'
      include 'inc/pyctag'
      include 'inc/pyismx'
      include 'inc/pyisjn'
C...Max size of hard system = HEPEUP size
      INTEGER MAXNUP
      PARAMETER (MAXNUP=500)
C...Local arrays and saved variables.
      DIMENSION VINTSV(11:80),KSAV(MAXNUP,5),PSAV(MAXNUP,5),
     &     VSAV(MAXNUP,5),SHAT(240)
      SAVE NSAV,NPARTS,M15SV,M16SV,M21SV,M22SV,VINTSV,SHAT,ISUBHD,ALAM3
     &     ,PSAV,KSAV,VSAV

C----------------------------------------------------------------------
C...MODE=-1: Pre-initialization. Store info on hard scattering etc,
C...done only once per event, while MODE=0 is repeated each time the
C...evolution needs to be restarted.
      IF (MODE.EQ.-1) THEN
        ISUBHD=MINT(1)
        NSAV=N
        NPARTS=NPART
C...Store hard scattering variables
        M15SV=MINT(15)
        M16SV=MINT(16)
        M21SV=MINT(21)
        M22SV=MINT(22)
        DO 100 J=11,80
          VINTSV(J)=VINT(J)
  100   CONTINUE
        DO 120 J=1,5
          DO 110 IS=1,NSAV-MINT(84)
            I=IS+MINT(84)
            PSAV(IS,J)=P(I,J)
            KSAV(IS,J)=K(I,J)
            VSAV(IS,J)=V(I,J)
  110     CONTINUE
  120   CONTINUE

C...Set shat for hardest scattering
        SHAT(1)=VINT(44)
        IF(ISET(ISUBHD).GE.3.AND.ISET(ISUBHD).LE.5) SHAT(1)=VINT(26)
     &       *VINT(2)

C...Compute 3-Flavour Lambda_QCD (sets absolute lowest PT scale below)
        RMC=PMAS(4,1)
        RMB=PMAS(5,1)
        ALAM4=PARP(61)
        IF(MSTU(112).LT.4) ALAM4=PARP(61)*(PARP(61)/RMC)**(2D0/25D0)
        IF(MSTU(112).GT.4) ALAM4=PARP(61)*(RMB/PARP(61))**(2D0/25D0)
        ALAM3=ALAM4*(RMC/ALAM4)**(2D0/27D0)

C----------------------------------------------------------------------
C...MODE= 0: Initialize ISR/MI evolution, i.e. begin from hardest
C...interaction initiators, with no previous evolution. Check the input
C...PT2MAX and PT2MIN and impose extra constraints on minimum PT2 (e.g.
C...must be larger than Lambda_QCD) and maximum PT2 (e.g. must be
C...smaller than the CM energy / 2.)
      ELSEIF (MODE.EQ.0) THEN
C...Reset counters and switches
        N=NSAV
        NPART=NPARTS
        MINT(30)=0
        MINT(31)=1
        MINT(36)=1
C...Reset hard scattering variables
        MINT(1)=ISUBHD
        DO 130 J=11,80
          VINT(J)=VINTSV(J)
  130   CONTINUE
        DO 150 J=1,5
          DO 140 IS=1,NSAV-MINT(84)
            I=IS+MINT(84)
            P(I,J)=PSAV(IS,J)
            K(I,J)=KSAV(IS,J)
            V(I,J)=VSAV(IS,J)
            P(MINT(83)+4+IS,J)=PSAV(IS,J)
            V(MINT(83)+4+IS,J)=VSAV(IS,J)
  140     CONTINUE
  150   CONTINUE
C...Reset statistics on activity in event.
        DO 160 J=351,359
          MINT(J)=0
          VINT(J)=0D0
  160   CONTINUE
C...Reset extra companion reweighting factor
        VINT(140)=1D0

C...We do not generate MI for soft process (ISUB=95), but the
C...initialization must be done regardless, for later purposes.
        MINT(36)=1

C...Initialize multiple interactions.
        CALL PYPTMI(-1,PTDUM1,PTDUM2,PTDUM3,IDUM)
        IF(MINT(51).NE.0) RETURN

C...Decide whether quarks in hard scattering were valence or sea
        PT2HD=VINT(54)
        DO 170 JS=1,2
          MINT(30)=JS
          CALL PYPTMI(2,PT2HD,PTDUM2,PTDUM3,IDUM)
          IF(MINT(51).NE.0) RETURN
  170   CONTINUE

C...Set lower cutoff for PT2 iteration and colour interference PT2 scale
        VINT(18)=0D0
        PT2MIN=MAX(PT2MIN,(1.1D0*ALAM3)**2)
        IF (MSTP(70).EQ.2) THEN
C...VINT(18) is freezeout scale of alpha_s: alpha_eff(0) = alpha_s(VINT(18))
          VINT(18)=(PARP(82)*(VINT(1)/PARP(89))**PARP(90))**2
        ELSEIF (MSTP(70).EQ.3) THEN
C...MSTP(70) = 3 : Derive VINT(18) from alpha_eff(Lambda3) = PARP(73)
          ALPHA0 = MAX(1D-6,PARP(73))
          Q20 = ALAM3**2/PARP(64)
          IF (MSTP(64).EQ.3) Q20 = Q20 * 1.661**2
          VINT(18) = Q20 * (EXP(12*PARU(1)/27D0/ALPHA0)-1D0)
        ENDIF
C...Also store PT2MIN in VINT(17).
        VINT(17)=PT2MIN

C...Set FS masses zero now.
        VINT(63)=0D0
        VINT(64)=0D0

C...Initialize IS showers with VINT(56) as max scale.
        PT2ISR=VINT(56)
        PT20=PT2MIN
        IF (MSTP(70).EQ.0) THEN
          PT20=MAX(PT2MIN,PARP(62)**2)
        ELSEIF (MSTP(70).EQ.1) THEN
          PT20=MAX(PT2MIN,(PARP(81)*(VINT(1)/PARP(89))**PARP(90))**2)
        ENDIF
        CALL PYPTIS(-1,PT2ISR,PT20,PT2DUM,IFAIL)
        IF(MINT(51).NE.0) RETURN

        RETURN

C----------------------------------------------------------------------
C...MODE= 1: Evolve event from PTMAX to PTMIN.
      ELSEIF (MODE.EQ.1) THEN

C...Skip if no phase space.
        IF (PT2MAX.LE.PT2MIN) GOTO 330

C...Starting pT2 max scale (to be udpated successively).
        PT2CMX=PT2MAX

C...Evolve two sides of the event to find which branches at highest pT.
  200   JSMX=-1
        MIMX=0
        PT2MX=0D0

C...Loop over current shower initiators.
        IF (MSTP(61).GE.1) THEN
          DO 230 MI=1,MINT(31)
            IF (MI.GE.2.AND.MSTP(84).LE.0) GOTO 230
            ISUB=96
            IF (MI.EQ.1) ISUB=ISUBHD
            MINT(1)=ISUB
            MINT(36)=MI
C...Set up shat, initiator x values, and x remaining in BR.
            VINT(44)=SHAT(MI)
            VINT(141)=XMI(1,MI)
            VINT(142)=XMI(2,MI)
            VINT(143)=1D0
            VINT(144)=1D0
            DO 210 JI=1,MINT(31)
              IF (JI.EQ.MINT(36)) GOTO 210
              VINT(143)=VINT(143)-XMI(1,JI)
              VINT(144)=VINT(144)-XMI(2,JI)
  210       CONTINUE
C...Loop over sides.
C...Generate trial branchings for this interaction. The hardest
C...branching so far is automatically updated if necessary in /PYISMX/.
            DO 220 JS=1,2
              MINT(30)=JS
              PT20=PT2MIN
              IF (MSTP(70).EQ.0) THEN
                PT20=MAX(PT2MIN,PARP(62)**2)
              ELSEIF (MSTP(70).EQ.1) THEN
                PT20=MAX(PT2MIN,
     &              (PARP(81)*(VINT(1)/PARP(89))**PARP(90))**2)
              ENDIF
              CALL PYPTIS(0,PT2CMX,PT20,PT2NEW,IFAIL)
              IF (MINT(51).NE.0) RETURN
  220       CONTINUE
  230     CONTINUE
        ENDIF

C...Generate trial additional interaction.
        MINT(36)=MINT(31)+1
        IF (MOD(MSTP(81),10).GE.1) THEN
          MINT(1)=96
C...Set up X remaining in BR.
          VINT(143)=1D0
          VINT(144)=1D0
          DO 250 JI=1,MINT(31)
            VINT(143)=VINT(143)-XMI(1,JI)
            VINT(144)=VINT(144)-XMI(2,JI)
  250     CONTINUE
C...Generate trial interaction
          CALL PYPTMI(0,PT2CMX,PT2MIN,PT2NEW,IFAIL)
          IF (MINT(51).EQ.1) RETURN
        ENDIF

C...And the winner is:
        IF (PT2MX.LT.PT2MIN) THEN
          GOTO 330
        ELSEIF (JSMX.EQ.0) THEN
C...Accept additional interaction (may still fail).
          CALL PYPTMI(1,PT2NEW,PT2MIN,PT2DUM,IFAIL)
          IF(MINT(51).NE.0) RETURN
          IF (IFAIL.EQ.0) THEN
            SHAT(MINT(36))=VINT(44)
C...Decide on flavours (valence/sea/companion).
            DO 270 JS=1,2
              MINT(30)=JS
              CALL PYPTMI(2,PT2NEW,PT2MIN,PT2DUM,IFAIL)
              IF(MINT(51).NE.0) RETURN
  270       CONTINUE
          ENDIF
        ELSEIF (JSMX.EQ.1.OR.JSMX.EQ.2) THEN
C...Reconstruct kinematics of acceptable ISR branching.
C...Set up shat, initiator x values, and x remaining in BR.
          MINT(30)=JSMX
          MINT(36)=MIMX
          VINT(44)=SHAT(MINT(36))
          VINT(141)=XMI(1,MINT(36))
          VINT(142)=XMI(2,MINT(36))
          VINT(143)=1D0
          VINT(144)=1D0
          DO 280 JI=1,MINT(31)
            IF (JI.EQ.MINT(36)) GOTO 280
            VINT(143)=VINT(143)-XMI(1,JI)
            VINT(144)=VINT(144)-XMI(2,JI)
  280     CONTINUE
          PT2NEW=PT2MX
          CALL PYPTIS(1,PT2NEW,PT2DM1,PT2DM2,IFAIL)
          IF (MINT(51).EQ.1) RETURN
        ELSEIF (JSMX.EQ.3.OR.JSMX.EQ.4) THEN
C...Bookeep joining. Cannot (yet) be constructed kinematically.
          MINT(354)=MINT(354)+1
          VINT(354)=VINT(354)+SQRT(PT2MX)
          IF (MINT(354).EQ.1) VINT(359)=SQRT(PT2MX)
          MJOIND(JSMX-2,MJN1MX)=MJN2MX
          MJOIND(JSMX-2,MJN2MX)=MJN1MX
        ENDIF

C...Update PT2 iteration scale.
        PT2CMX=PT2MX

C...Loop back to continue evolution.
        IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
          CALL PYERRM(11,'(PYEVOL:) no more memory left in PYJETS')
        ELSE
          IF (JSMX.GE.0.AND.PT2CMX.GE.PT2MIN) GOTO 200
        ENDIF

C----------------------------------------------------------------------
C...MODE= 2: (Re-)store user information on hardest interaction etc.
      ELSEIF (MODE.EQ.2) THEN

C...Revert to "ordinary" meanings of some parameters.
        DO 310 JS=1,2
          MINT(12+JS)=K(IMI(JS,1,1),2)
          VINT(140+JS)=XMI(JS,1)
          IF(MINT(18+JS).EQ.1) VINT(140+JS)=VINT(154+JS)*XMI(JS,1)
          VINT(142+JS)=1D0
          DO 300 MI=1,MINT(31)
            VINT(142+JS)=VINT(142+JS)-XMI(JS,MI)
  300     CONTINUE
  310   CONTINUE

C...Restore saved quantities for hardest interaction.
        MINT(1)=ISUBHD
        MINT(15)=M15SV
        MINT(16)=M16SV
        MINT(21)=M21SV
        MINT(22)=M22SV
        DO 320 J=11,80
          VINT(J)=VINTSV(J)
  320   CONTINUE

      ENDIF

  330 RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      subroutine setupmom(p)
      implicit none
      include 'constants.f'
      include 'masses.f'
      double precision al,be,ro,si,th,ph,p(mxpart,4),E1,E2,E7,y

C----setup p
      al=1.04225677d0
      be=0.5564225677d0
      ro=0.86454322d0
      si=1.04225677d0
      th=0.63425677d0
      ph=0.145563425677d0
      zmass=91.188d0
      p(:,:)=0d0

      E1=343.43678d0
      E2=436.68843678d0
      y=0.5d0*(1d0+cos(th))
      E7=-(zmass**2-zmass*E1-zmass*E2+E1*E2)/(zmass-E1*y-E2*(1d0-y))

      p(1,4)=-E1
      p(1,1)=0d0
      p(1,2)=0d0
      p(1,3)=+E1

      p(2,4)=-E2
      p(2,1)=0d0
      p(2,2)=0d0
      p(2,3)=-E2

      p(3,4)=+0.5d0*zmass
      p(3,1)=+0.5d0*zmass*sin(al)*cos(ro)
      p(3,2)=+0.5d0*zmass*sin(al)*sin(ro)
      p(3,3)=+0.5d0*zmass*cos(al)

      p(4,4)=+0.5d0*zmass
      p(4,1)=-0.5d0*zmass*sin(al)*cos(ro)
      p(4,2)=-0.5d0*zmass*sin(al)*sin(ro)
      p(4,3)=-0.5d0*zmass*cos(al)

      p(5,4)=+0.5d0*zmass
      p(5,1)=+0.5d0*zmass*sin(be)*cos(si)
      p(5,2)=+0.5d0*zmass*sin(be)*sin(si)
      p(5,3)=+0.5d0*zmass*cos(be)

      p(6,4)=+0.5d0*zmass
      p(6,1)=-0.5d0*zmass*sin(be)*cos(si)
      p(6,2)=-0.5d0*zmass*sin(be)*sin(si)
      p(6,3)=-0.5d0*zmass*cos(be)

      p(7,4)=E7
      p(7,1)=E7*sin(th)*sin(ph)
      p(7,2)=E7*sin(th)*cos(ph)
      p(7,3)=E7*cos(th)

      p(8,:)=-p(1,:)-p(2,:)-p(3,:)-p(4,:)-p(5,:)-p(6,:)-p(7,:)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!  Description: set string key with spaces
!
!
USE eccodes

 INTEGER :: IGRIBH
 CHARACTER (LEN=17) :: CLNOMA1, CLNOMA2

 CALL codes_grib_new_from_samples(IGRIBH, "regular_ll_sfc_grib2")

 ! set centre to MeteoFrance and use their local definition
 CALL codes_set (IGRIBH, 'centre', 85)
 CALL codes_set (IGRIBH, 'grib2LocalSectionPresent', 1)
 CALL codes_set (IGRIBH, 'grib2LocalSectionNumber', 1)

 CLNOMA1 = 'SUNSHI. DURATION'
 CALL codes_set_string(IGRIBH, 'CLNOMA', CLNOMA1)
 CALL codes_get_string(IGRIBH, 'CLNOMA', CLNOMA2)

 PRINT *, " CLNOMA1 = ", CLNOMA1
 PRINT *, " CLNOMA2 = ", CLNOMA2

 IF (CLNOMA1 /= CLNOMA2) STOP 1

END

max_stars_repo_name: bbakernoaa/eccodes

max_issues_repo_name: bbakernoaa/eccodes

max_forks_repo_name: bbakernoaa/eccodes


text:
************************************************************************
* This file contains the OMEGA parametrisation routines for the
* computational domain:
*   OPARX - Obtain the X-coordinate of a point from its parameter value
*   OPARX - Obtain the Y-coordinate of a point from its parameter value
*   OTMAX - Obtain the maximum parameter value on a boundary component
*
* The implementation here is used as the standard implementation in case
* of that the domain is read in by RDPARM from an OMEGA-style file.
************************************************************************
*   parxypre.f:
*     - contains the parametrization functions:
*           OPARX, OPARY, OTMAX, OTNBC
*       for the preprocessing variant ('rdparm')
*     - the COMMON-block /TDATA/ is used frequently and has to be
*       initialized by the modul 'initxx.f'
************************************************************************

************************************************************************
* Get X-coordinate of parameter value.
*
* This routine returns the X-coordinate of a parameter value.
*
* In:
*   T      - parameter value
*   IBCT   - number of boundary component
*
* Return:
*   X-coordinate of parameter value T on boundary component IBCT.
************************************************************************

      DOUBLE PRECISION FUNCTION OPARX(T,IBCT)

      IMPLICIT NONE

      INCLUDE 'cmem.inc'
      INCLUDE 'cparqdata.inc'

C parameters

      DOUBLE PRECISION T
      INTEGER IBCT

C local variables

      INTEGER NCOMP,ICPTR,ICOMP,IICOMP,ITYP,NSPLIN,NPAR,IPPTR
      INTEGER KXPAR,KYPAR
      DOUBLE PRECISION TDIFF,X1,XDIFF,XM,R,PHI,PHI1,PHI2,PHI3,PHI4

      INTEGER ISPLIN,KK
      DOUBLE PRECISION TTDIF,TT,C1,C2,C3,C4

C=======================================================================
C *** IBCT dependent quantities
      NCOMP=KWORK(L(LNCOMP)+IBCT-1)
      ICPTR=KWORK(L(LICPTR)+IBCT-1)
C
C *** determine boundary component ICOMP of part IBCT and
C     the ICOMP-dependent quantities
C
      ICOMP=T+1.d0
      tdiff=T-ICOMP+1.d0
      if(tdiff.lt.0d0 .or. tdiff.ge.1d0) then
        write(*,*)'error in OPARX (userpre): conflict with ICOMP'
        stop
      endif
      iicomp=ICPTR+ICOMP-1
      ITYP=KWORK(L(LITYP)+iicomp)
      NSPLIN=KWORK(L(LNSPLN)+iicomp)
      NPAR=KWORK(L(LNPAR)+iicomp)
      IPPTR=KWORK(L(LIPPTR)+iicomp)
C
C *** pointer for the LXPAR/LYPAR-arrays
      KXPAR=L(LXPAR)+IPPTR
      KYPAR=L(LYPAR)+IPPTR
C=======================================================================
C   ITYP=1:  line
C=======================================================================
      if(ITYP.eq.1) then
C
        if(NSPLIN.ne.1 .or. NPAR.ne.2) then
          write(*,*)'error in OPARX (userpre): conflict at ITYP=1'
          stop
        endif
        x1=DWORK(KXPAR)
        xdiff=DWORK(KXPAR+1)
        OPARX=x1 + xdiff*tdiff
C=======================================================================
C   ITYP=2:  part of a circle
C=======================================================================
      else if(ITYP.eq.2) then
C
        if(NSPLIN.ne.1 .or. NPAR.ne.3) then
          write(*,*)'error in OPARX (userpre): conflict at ITYP=2'
          stop
        endif
        xm=DWORK(KXPAR)
        r=DWORK(KXPAR+1)
        phi1=DWORK(KXPAR+2)
        phi2=DWORK(KYPAR+2)
        phi=phi1 + tdiff*(phi2-phi1)
        OPARX=xm + r*COS(phi)
C=======================================================================
C   ITYP=3:  spline
C=======================================================================
      else if(ITYP.eq.3) then
C
C *** determine "number of the subspline -1"=ISPLIN and tt=t_tilde
        ttdif=tdiff*NSPLIN
        ISPLIN=ttdif
        tt=ttdif-ISPLIN
        if(tt.lt.0d0 .or. tt.ge.1d0 .or. ISPLIN.gt.NSPLIN-1) then
          write(*,*)'error in OPARX (userpre): conflict with ISPLIN'
          stop
        endif
C
C *** determine cubic spline-function values at 'tt'
        phi1=(2d0*tt*tt - 3d0*tt)*tt + 1d0
        phi2=(-2d0*tt + 3d0)*tt*tt
        phi3=(tt*tt - 2d0*tt + 1d0)*tt
        phi4=(tt - 1d0)*tt*tt
C
C *** determine parameter constants
        kk=KXPAR + 4*ISPLIN
        c1=DWORK(kk)
        c2=DWORK(kk+1)
        c3=DWORK(kk+2)
        c4=DWORK(kk+3)
C
        OPARX=c1*phi1+c2*phi2+c3*phi3+c4*phi4
C=======================================================================
      else
        write(*,*)'error in OPARX (userpre): wrong ITYP-value'
        stop
      endif
C
      END

************************************************************************
* Get Y-coordinate of parameter value.
*
* This routine returns the Y-coordinate of a parameter value.
*
* In:
*   T      - parameter value
*   IBCT   - number of boundary component
*
* Return:
*   Y-coordinate of parameter value T on boundary component IBCT.
************************************************************************

      DOUBLE PRECISION FUNCTION OPARY(T,IBCT)

      IMPLICIT NONE

      INCLUDE 'cmem.inc'
      INCLUDE 'cparqdata.inc'

C parameters

      DOUBLE PRECISION T
      INTEGER IBCT

C local variables

      INTEGER NCOMP,ICPTR,ICOMP,IICOMP,ITYP,NSPLIN,NPAR,IPPTR
      INTEGER KXPAR,KYPAR
      DOUBLE PRECISION TDIFF,R,PHI,PHI1,PHI2,PHI3,PHI4

      INTEGER ISPLIN,KK
      DOUBLE PRECISION TTDIF,TT,C1,C2,C3,C4,Y1,YDIFF,YM
C=======================================================================
C *** IBCT dependent quantities
      NCOMP=KWORK(L(LNCOMP)+IBCT-1)
      ICPTR=KWORK(L(LICPTR)+IBCT-1)
C
C *** determine boundary component ICOMP of part IBCT and
C     the ICOMP-dependent quantities
C
      ICOMP=T+1.d0
      tdiff=T-ICOMP+1.d0
      if(tdiff.lt.0d0 .or. tdiff.ge.1d0) then
        write(*,*)'error in OPARY (userpre): conflict with ICOMP'
        stop
      endif
      iicomp=ICPTR+ICOMP-1
      ITYP=KWORK(L(LITYP)+iicomp)
      NSPLIN=KWORK(L(LNSPLN)+iicomp)
      NPAR=KWORK(L(LNPAR)+iicomp)
      IPPTR=KWORK(L(LIPPTR)+iicomp)
C
C *** pointer for the LXPAR/LYPAR-arrays
      KXPAR=L(LXPAR)+IPPTR
      KYPAR=L(LYPAR)+IPPTR
C=======================================================================
C   ITYP=1:  line
C=======================================================================
      if(ITYP.eq.1) then
C
        if(NSPLIN.ne.1 .or. NPAR.ne.2) then
          write(*,*)'error in OPARY (userpre): conflict at ITYP=1'
          stop
        endif
        y1=DWORK(KYPAR)
        ydiff=DWORK(KYPAR+1)
        OPARY=y1 + ydiff*tdiff
C=======================================================================
C   ITYP=2:  part of a circle
C=======================================================================
      else if(ITYP.eq.2) then
C
        if(NSPLIN.ne.1 .or. NPAR.ne.3) then
          write(*,*)'error in OPARY (userpre): conflict at ITYP=2'
          stop
        endif
        ym=DWORK(KYPAR)
        r=DWORK(KXPAR+1)
        phi1=DWORK(KXPAR+2)
        phi2=DWORK(KYPAR+2)
        phi=phi1 + tdiff*(phi2-phi1)
        OPARY=ym + r*SIN(phi)
C=======================================================================
C   ITYP=3:  spline
C=======================================================================
      else if(ITYP.eq.3) then
C
C *** determine "number of the subspline -1"=ISPLIN and tt=t_tilde
        ttdif=tdiff*NSPLIN
        ISPLIN=ttdif
        tt=ttdif-ISPLIN
        if(tt.lt.0d0 .or. tt.ge.1d0 .or. ISPLIN.gt.NSPLIN-1) then
          write(*,*)'error in OPARY (userpre): conflict with ISPLIN'
          stop
        endif
C
C *** determine cubic spline-function values at 'tt'
        phi1=(2d0*tt*tt - 3d0*tt)*tt + 1d0
        phi2=(-2d0*tt + 3d0)*tt*tt
        phi3=(tt*tt - 2d0*tt + 1d0)*tt
        phi4=(tt - 1d0)*tt*tt
C
C *** determine parameter constants
        kk=KYPAR + 4*ISPLIN
        c1=DWORK(kk)
        c2=DWORK(kk+1)
        c3=DWORK(kk+2)
        c4=DWORK(kk+3)
C
        OPARY=c1*phi1+c2*phi2+c3*phi3+c4*phi4
C=======================================================================
      else
        write(*,*)'error in OPARY (userpre): wrong ITYP-value'
        stop
      endif
C
      END

************************************************************************
* Get maximum parameter value
*
* This routine returns the maximum parameter value on a boundary
* component
*
* In:
*   IBCT   - number of boundary component
*
* Return:
*   Maximum parameter value on that boundary component.
************************************************************************

      DOUBLE PRECISION FUNCTION OTMAX(IBCT)

      IMPLICIT NONE

      INCLUDE 'cmem.inc'
      INCLUDE 'cparqdata.inc'

C parameters

      INTEGER IBCT

C local variables

      INTEGER NCOMP

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

      NCOMP=KWORK(L(LNCOMP)+IBCT-1)

      OTMAX=NCOMP

      END

************************************************************************
* Get number of real boundary components
*
* This routine returns the number of real boundary components in the
* parametrization.
*
* In:
*   -
*
* Return:
*   Number of boundary components
************************************************************************

      INTEGER FUNCTION OTNBC ()

      IMPLICIT NONE

      INCLUDE 'cparqdata.inc'

      OTNBC = NBCT

      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
Animal Science is one of the many majors offered at UCD.  It is affiliated with the Department of Animal Science and the College of Ag College of Agriculture and Environmental Science.  The fact that UC Davis is a research campus makes this one of the best places for Animal Science majors since the campus houses the animals that the students are studying.  Resources like the Cole Facility and the Vet Med Teaching Hospital (VMTH) provide students with learning experiences that are difficult to find elsewhere.   Fistulated Cow Fistulated Cows, http://www.news.ucdavis.edu/dateline/051702/dl_calf.html cloned cows and http://en.wikipedia.org/wiki/Geep geeps are just some of the more unorthodox beasts you will encounter in animal science courses. Cows Cow clones have been created on campus.

Also, Davis is one of few schools in the country that incorporates small/companion animals into their program.

The department has one of the largest undergraduates undergraduate enrollments at UC Davis, with approximately 850 students currently enrolled.

Many Animal Science majors at UC Davis are pursuing acceptance to the School of Veterinary Medicine.

Animal Science majors must each choose one of the seven specializations:

  Aquatic Animals
  Avian Sciences
  Companion and Captive Animals
  Equine Science
  Laboratory Animals
  Livestock and Dairy
  Poultry