text:
      SUBROUTINE FA1
C
C     FA1 IS THE DRIVER FOR PART ONE OF FLUTTER ANALYSIS
C
      INTEGER         SYSBUF,OUT,BUFF,BUFF1,NS(2),FLOOP,TSTART,
     1                KHH,BHH,MHH,QHHL,CASECC,FLIST,FSAVE,KXHH,MXHH,
     2                BXHH,SCR1,REC0(8),FLUT(10),IMETH(2),FMETHD,SMETH,
     3                TRL(10),AERO(2),FLFACT(2),FLUTER(2),
     4                SR,SM,SK,PR,PM,PK,SL,IBLOCK(12),METHOD(4)
      REAL            BLOCK(12),REC(8),KFREQ,RHO
      DIMENSION       DLT(3),Z(1)
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /SYSTEM/ SYSBUF,OUT
      COMMON /BLANK / FLOOP,TSTART,ICEAD
      COMMON /OUTPUT/ HDG(96)
      COMMON /PACKX / ITI,ITO,IJ,NN,INCR1
      COMMON /ZZZZZZ/ IZ(1)
      COMMON /UNPAKX/ IOUT,II,JJ,INCR
      EQUIVALENCE     (REC0(1),REC(1))
      EQUIVALENCE     (BLOCK(1),IBLOCK(1))
      EQUIVALENCE     (IZ(1),Z(1))
      DATA  KHH /101/,   BHH /102/,  MHH /103/, QHHL /104/, CASECC /105/
     1,   FLIST /106/, FSAVE /201/, KXHH /202/, BXHH /203/,   MXHH /204/
      DATA  SCR1/301/,    NS /4HFA1 ,4H    /
      DATA  IMETH   /  4HS   ,4HL   /
      DATA  NMD /4  /, METHOD/4HK   ,4HKE  ,4HPK  ,4HINV /
      DATA  TRL     /  90,1006,7*0,6 /
      DATA  AERO    /  3202,32 /,   FLFACT /4102,41/,   FLUTER /3902,39/
C
      DO 5 I = 1,12
    5 IBLOCK(I) = 0
      NCORE = KORSZ(IZ)
      BUFF  = NCORE - SYSBUF - 1
      BUFF1 = BUFF  - SYSBUF
      IF (FLOOP .NE. 0) GO TO 200
C
C     FIRST TIME THROUGH FIND FMETHOD ON CASECC
C
      IFILE = CASECC
      CALL GOPEN (CASECC,IZ(BUFF+1),0)
      CALL READ (*530,*10,CASECC,IZ,BUFF,1,NWR)
   10 LCC = NWR
      CALL CLOSE (CASECC,1)
C
C     GET DATA FOR REC0 OF FSAVE
C
      CALL FNAME (FSAVE,REC0)
      IFILE = FLIST
      CALL PRELOC (*480,IZ(BUFF+1),FLIST)
      CALL LOCATE (*470,IZ(BUFF+1),AERO,IDUM)
      CALL READ (*530,*530,FLIST,REC0(4),4,1,NWR)
      REC(6) = REC(6)*0.5
      CALL LOCATE (*15,IZ(BUFF+1),FLFACT,IDUM)
      CALL READ (*530,*20,FLIST,IZ(LCC+1),BUFF,1,NWR)
      GO TO 400
   15 NWR = 0
   20 LFL = NWR + LCC
      CALL LOCATE (*450,IZ(BUFF+1),FLUTER,IDUM)
   30 CALL READ (*530,*450,FLIST,FLUT,10,0,NWR)
      I165 = 165
      IF (FLUT(1) .NE. IZ(I165)) GO TO 30
      CALL CLOSE (FLIST,1)
      REC0(8) = FLUT(9)
      IEP = FLUT(10)
      DO 40 I = 1,NMD
      IF (FLUT(2) .EQ. METHOD(I)) GO TO 50
   40 CONTINUE
      GO TO 490
   50 REC0(3) = I
      FMETHD = I
      GO TO (60,60,61,490), I
   60 REC0(4) = 0
      IF (FLUT(7) .EQ. IMETH(1)) REC0(4) = 1
      IF (FLUT(7) .EQ. IMETH(2)) REC0(4) = 2
      IF (REC0(4) .EQ. 0) GO TO 430
      SMETH = REC0(4)
      GO TO 65
C
C     PK METHOD HAS LINEAR SPLINE ONLY
C
   61 REC0(4) = 2
      SMETH = 2
   65 CONTINUE
C
C     BUILD RECORDS 0,1,2,3 OF SAVE
C
      IFILE = FSAVE
      CALL OPEN (*480,FSAVE,IZ(BUFF+1),1)
      CALL WRITE (FSAVE,REC0,8,1)
      BREF = REC(6)
      RREF = REC(7)
      NEIW = REC0(8)
C
C     BUILD M,K,RHO LIST FOR FLUTTER LOOP
C
      SR = 0
      SM = 0
      SK = 0
      I  = LCC
      IF (I .EQ. LFL) GO TO 410
   70 I  = I + 1
      IF (IZ(I) .EQ. FLUT(4)) SR = I
      IF (IZ(I) .EQ. FLUT(5)) SM = I
      IF (IZ(I) .EQ. FLUT(6)) SK = I
   80 I = I + 1
      IF (I .GE. LFL) GO TO 90
      IF (IZ(I) .EQ. -1) GO TO 70
      GO TO 80
   90 IF (SR.EQ.0 .OR. SM.EQ.0 .OR. SK.EQ.0) GO TO 410
      NRHO = 0
      PR = SR
   95 PR = PR + 1
      IF (IZ(PR) .EQ. -1) GO TO 97
      NRHO = NRHO + 1
      GO TO 95
   97 NLOOPS = 0
      IF (FMETHD .NE. 3) GO TO 105
C
C     J.PETKAS/LOCKHEED      3/91
C     19 LINES OF OLD CODE FOR BUILDING ELEMENTS OF FSAVE FOR PK METHOD
C     WERE IN ERROR, AND ARE NOW REPLACED BY NEXT 29 NEW LINES
C
      PM = SM
  101 PM = PM + 1
      IF (IZ(PM) .EQ. -1) GO TO 130
      DLT(1) = Z(PM)
C
C     CENTER LOOP ON RHO
C
      PR = SR
  102 PR = PR + 1
      IF (IZ(PR) .EQ. -1) GO TO 101
      DLT(3) = Z(PR)
C
C     INNER LOOP ON VELOCITY
C
      PK = SK
  103 PK = PK + 1
      IF (IZ(PK) .EQ. -1) GO TO 102
      DLT(2) = Z(PK)
      NLOOPS = NLOOPS + 1
      CALL WRITE (FSAVE,DLT,3,0)
      GO TO 103
C
C     ALGORITHM FOR BUILDING ELEMENTS OF FSAVE FOR K AND KE METHODS
C
  105 CONTINUE
C
C     OUTER LOOP ON MACH NUMBER
C
      PM = SM
  107 PM = PM + 1
      IF (IZ(PM) .EQ. -1) GO TO 130
      DLT(1) = Z(PM)
C
C     CENTER LOOP ON KFREQ
C
      PK = SK
  110 PK = PK + 1
      IF (IZ(PK) .EQ. -1) GO TO 107
      DLT(2) = Z(PK)
C
C     INNER LOOP ON RHO
C
      PR = SR
  120 PR = PR + 1
      IF (IZ(PR) .EQ. -1) GO TO 110
      DLT(3) = Z(PR)
      NLOOPS = NLOOPS + 1
      CALL WRITE (FSAVE,DLT,3,0)
      GO TO 120
  130 CALL WRITE (FSAVE,0,0,1)
C
C     PICK UP M AND K FROM QHHL
C
      IFILE = QHHL
      CALL OPEN (*480,QHHL,IZ(BUFF1+1),0)
      CALL READ (*530,*140,QHHL,IZ(LCC+1),BUFF1,1,NWR)
      GO TO 400
  140 LFL = NWR + LCC
      SL  = LCC + 5
      CALL CLOSE (QHHL,1)
      REC0(1) = QHHL
      CALL RDTRL (REC0)
      NP  = MIN0(IZ(SL-1),REC0(2)/REC0(3))
      LFL = MIN0(LFL,2*NP+SL-1)
      NP  = LFL - SL + 1
      CALL WRITE (FSAVE,IZ(SL),NP,1)
      NP  = NP/2
C
C     WRITE CASECC RECORD AND TRAILER
C
      CALL WRITE (FSAVE,IZ(1),LCC,1)
      CALL CLOSE (FSAVE,1)
      REC0(1) = FSAVE
      REC0(2) = FLOOP
      REC0(3) = NLOOPS
      REC0(4) = NP
      REC0(5) = LCC
      REC0(6) = 0
      REC0(7) = NRHO
      CALL WRTTRL (REC0)
      GO TO 210
  200 IFILE = FSAVE
      CALL OPEN (*480,FSAVE,IZ(BUFF+1),0)
      CALL READ (*530,*530,FSAVE,IZ(1),8,1,NWR)
      CALL CLOSE (FSAVE,1)
      IZX   = 0
      FMETHD= IZ(IZX+3)
      SMETH = IZ(IZX+4)
      BREF  =  Z(IZX+6)
      RREF  =  Z(IZX+7)
      NEIW  = IZ(IZX+8)
  210 REC0(1) = FSAVE
      CALL RDTRL (REC0)
C
C     START OF LOOPING BUMP LOOP COUNTER SET TIME AND GO
C
      FLOOP = FLOOP + 1
      NLOOPS= REC0(3)
      CALL KLOCK (TSTART)
      GO TO (220,230,240,490), FMETHD
C
C     K METHOD BUILD PROPER QHH ON SCR1
C
  220 CALL FA1K (SMETH,KFREQ,RHO,SCR1,0)
      GO TO 300
C
C     KE METHOD DO INCORE EIGNVALUE EXTRACTION
C
  230 REC0(1) = BHH
      CALL RDTRL (REC0)
      IF (REC0(1).GT.0 .AND. REC0(7).GT.0) GO TO 510
      REC0(1) = KHH
      CALL RDTRL (REC0)
      ICO = REC0(2)*REC0(2)*4 + 4
  235 CALL FA1K (SMETH,KFREQ,RHO,SCR1,ICO)
      CALL FA1KE (SCR1,KFREQ,BREF,RHO,RREF,FLOOP,NLOOPS)
      IF (FLOOP .GE. NLOOPS) GO TO 350
      FLOOP = FLOOP + 1
      GO TO 235
C
C     PK METHOD  LINEAR INTERPOLATION  AND INCORE LOOP FOR
C     EIGENVALUE CONVERGENCE
C
  240 CALL FA1PKI (FSAVE,QHHL)
      CALL FA1PKE (KHH,BHH,MHH,BXHH,FSAVE,NLOOPS,BREF,RREF,NEIW,IEP)
      IF (FLOOP .GE. NLOOPS) GO TO 250
      FLOOP = FLOOP  + 1
      GO TO 240
C
C     PHID  - KXHH   CLAMAD - BXHH
C
  250 IBUF = BUFF1 - SYSBUF
      TRL(1) = SCR1
      CALL RDTRL (TRL)
      IF (TRL(2) .EQ. 0) GO TO 350
      CALL OPEN (*350,SCR1,Z(IBUF),0)
      CALL READ (*290,*255,SCR1,REC,6,1,NWR)
  255 CALL READ (*290,*260,SCR1,Z,IBUF,1,NWR)
  260 NN = NWR/2
      CALL GOPEN (KXHH,Z(BUFF),1)
      CALL GOPEN (BXHH,Z(BUFF1),1)
      CALL WRITE (BXHH,TRL(1),50,0)
      CALL WRITE (BXHH,HDG,96,1)
      TRL(1) = KXHH
      TRL(2) = 0
      TRL(3) = NN
      TRL(4) = 2
      TRL(5) = 3
      ITI = 3
      ITO = 3
      IJ  = 1
      INCR1 = 1
  265 CALL WRITE (BXHH,REC,6,0)
      CALL PACK  (Z,KXHH,TRL)
      CALL READ  (*280,*270,SCR1,REC,6,1,NWR)
  270 CALL READ  (*280,*265,SCR1,Z,IBUF,1,NWR)
  280 CALL WRITE (BXHH,0,0,1)
      CALL CLOSE (BXHH,1)
      CALL CLOSE (KXHH,1)
      CALL WRTTRL (TRL)
      TRL(1) = BXHH
      TRL(2) = 1006
      TRL(7) = 0
      CALL WRTTRL (TRL)
  290 CALL CLOSE (SCR1,1)
      GO TO 350
C
C     COPY KHH TO KXHH
C
  300 CALL GOPEN (KHH,IZ(BUFF+1),0)
      CALL GOPEN (KXHH,IZ(BUFF1+1),1)
      REC0(1) = KHH
      CALL RDTRL (REC0)
      REC0(1) = KXHH
      IOUT = REC0(5)
      INCR = 1
      I = REC0(2)
      REC0(2) = 0
      REC0(6) = 0
      REC0(7) = 0
      CALL CYCT2B (KHH,KXHH,I,Z,REC0)
      CALL CLOSE  (KHH,1)
      CALL CLOSE  (KXHH,1)
      CALL WRTTRL (REC0)
C
C     BUILD BXHH = (K/B)BHH
C
      REC0(1) = BHH
      CALL RDTRL (REC0)
      IF (REC0(1) .LE. 0) GO TO 310
      IBLOCK(2) = 1
      BLOCK(3) = KFREQ/BREF
      CALL SSG2C (BHH,0,BXHH,0,BLOCK(2))
  310 CONTINUE
C
C                2  2
C     MXHH  =  (K /B ) MHH  + (RHO*RREF/2.0) QHH
C
      IBLOCK(2) = 1
      BLOCK (3) = (KFREQ*KFREQ)/(BREF*BREF)
      IBLOCK(8) = 1
      BLOCK (9) = RHO*RREF/2.0
      CALL SSG2C (MHH,SCR1,MXHH,0,BLOCK(2))
C
C     THE END
C
  350 CONTINUE
      REC0(1) = FSAVE
      CALL RDTRL (REC0)
      REC0(2) = FLOOP
      CALL WRTTRL (REC0)
      IF (FLOOP .EQ. NLOOPS) FLOOP = -1
      ICEAD = 1
      IF (FMETHD .EQ. 2) ICEAD = -1
      IF (FMETHD .EQ. 3) ICEAD = -1
      GO TO 600
C
C     ERROR MESSAGES
C
  400 CALL MESAGE (-8,0,NS)
  410 WRITE  (OUT,420) UFM,FLUT(4),FLUT(5),FLUT(6)
  420 FORMAT (A23,', ONE OR MORE OF THE FOLLOWING FLFACT SETS WERE NOT',
     1       ' FOUND - ',3I9)
      GO TO 540
  430 WRITE  (OUT,440) UFM,FLUT(7)
  440 FORMAT (A23,' 2267, INTERPOLATION METHOD ',A4,' UNKNOWN')
      GO TO 540
  450 I165 = 165
      WRITE  (OUT,460) UFM,IZ(I165)
  460 FORMAT (A23,' 2268, FMETHOD SET',I9,' NOT FOUND')
      GO TO 540
  470 CALL MESAGE (-7,0,NS)
  480 CALL MESAGE (-1,IFILE,NS)
  490 WRITE  (OUT,500) UFM,FLUT(2)
  500 FORMAT (A23,' 2269, FLUTTER METHOD ',A4,' NOT IMPLEMENTED')
      GO TO 540
  510 WRITE  (OUT,520) UFM,FLUT(2)
  520 FORMAT (A23,', FLUTTER METHOD ',A4,' NOT IMPLEMENTED WITH B ',
     1       'MATRIX')
      GO TO 540
  530 CALL MESAGE (-3,IFILE,NS)
  540 CALL MESAGE (-61,0,NS)
C
  600 RETURN
      END


text:
Crickets feats this year

Chased Koko 20 ft for no reason
Got moved to protected contact
Got out of his cage a bazillion times
Convinced everyone that if he were human he would undoubtedly be a murderer

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
An apiary, or Bees bee yard, is a place where apiarists, or bee keepers, keep lots of Honeybees bee hives.

Apiaries in Davis
  Harry H. Laidlaw Jr. Honey Bee Research Facility in UCDs west campus area.
  The Domes
  The Tricoops
  Davis Bee Charmers
  Bee Charmers Apiary
  Sola Bee Farms


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> \brief
!> This module contains the classes and procedures for setting up the `delayedRejectionScaleFactorVec` attribute of samplers of class [ParaDRAM_type](@ref paradram_mod::paradram_type).
!> For more information, see the description of this attribute in the body of the module.
!> \author Amir Shahmoradi

module SpecDRAM_DelayedRejectionScaleFactorVec_mod

    use Constants_mod, only: RK
    implicit none

    character(*), parameter         :: MODULE_NAME = "@SpecDRAM_DelayedRejectionScaleFactorVec_mod"

    real(RK), allocatable           :: DelayedRejectionScaleFactorVec(:) ! namelist input

    type                            :: DelayedRejectionScaleFactorVec_type
        real(RK), allocatable       :: Val(:)
        real(RK)                    :: def
        real(RK)                    :: null
        character(:), allocatable   :: desc
    contains
        procedure, pass             :: set, checkForSanity, nullifyNameListVar
    end type DelayedRejectionScaleFactorVec_type

    interface DelayedRejectionScaleFactorVec_type
        module procedure            :: construct
    end interface DelayedRejectionScaleFactorVec_type

    private :: construct, set, checkForSanity, nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function construct(nd,methodName) result(self)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: construct
#endif
        use Decoration_mod, only: TAB
        use Constants_mod, only: IK, RK, NULL_RK
        use String_mod, only: num2str
        implicit none
        integer(IK), intent(in)                     :: nd
        character(*), intent(in)                    :: methodName
        type(DelayedRejectionScaleFactorVec_type)   :: self
        self%def  = 0.5_RK**(1._RK/real(nd,kind=RK))  ! This gives a half volume to the covariance
        self%null = NULL_RK
        self%desc = &
        "delayedRejectionScaleFactorVec is a real-valued positive vector of length (1:delayedRejectionCount) by which &
        &the covariance matrix of the proposal distribution of "//methodName//" sampler is scaled when the Delayed Rejection (DR) &
        &scheme is activated (by setting delayedRejectionCount>0). At each ith stage of the DR process, &
        &the proposal distribution from the last stage is scaled by the factor delayedRejectionScaleFactorVec(i). &
        &Missing elements of the delayedRejectionScaleFactorVec in the input to "//methodName//" will be set to the default value. &
        &The default value at all stages is 0.5^(1/ndim) = "//num2str(self%def)//", which reduces the &
        &volume of the covariance matrix of the proposal from the last DR stage by one half. The variable ndim represents the &
        &number of dimensions of the Domain of the objective function."
    end function construct

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine nullifyNameListVar(self)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: nullifyNameListVar
#endif
        use SpecDRAM_DelayedRejectionCount_mod, only: MAX_DELAYED_REJECTION_COUNT
        use Constants_mod, only: IK
        implicit none
        class(DelayedRejectionScaleFactorVec_type), intent(in) :: self
        if (allocated(DelayedRejectionScaleFactorVec)) deallocate(DelayedRejectionScaleFactorVec)
        allocate(DelayedRejectionScaleFactorVec(MAX_DELAYED_REJECTION_COUNT), source = self%null)
    end subroutine nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pure subroutine set(self,delayedRejectionCount, DelayedRejectionScaleFactorVec)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: set
#endif
        use Constants_mod, only: IK, RK
        implicit none
        class(DelayedRejectionScaleFactorVec_type), intent(inout)   :: self
        integer(IK) , intent(in)                                    :: delayedRejectionCount
        real(RK)    , intent(in), optional                          :: DelayedRejectionScaleFactorVec(:) ! This must remain colon
        if (present(DelayedRejectionScaleFactorVec)) then
            self%Val = pack( DelayedRejectionScaleFactorVec, mask = DelayedRejectionScaleFactorVec /= self%null )
            if (size(self%Val) == 0) then
                deallocate(self%Val)
                allocate(self%Val(delayedRejectionCount), source = self%def)
            end if
        else
            if (allocated(self%Val)) then
                if (size(self%Val) == 0) then
                    deallocate(self%Val)
                    allocate(self%Val(delayedRejectionCount), source = self%def)
                end if
            end if
        end if
    end subroutine set

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine checkForSanity(self,Err,methodName,delayedRejectionCount)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: checkForSanity
#endif
        use Constants_mod, only: IK, RK
        use Err_mod, only: Err_type
        use String_mod, only: num2str
        implicit none
        class(DelayedRejectionScaleFactorVec_type), intent(in)  :: self
        type(Err_type), intent(inout)                           :: Err
        character(*), intent(in)                                :: methodName
        integer(IK), intent(in)                                 :: delayedRejectionCount
        character(*), parameter                                 :: PROCEDURE_NAME = "@checkForSanity()"
        integer(IK)                                             :: i
        if ( size(self%Val)/=delayedRejectionCount ) then
            Err%occurred = .true.
            Err%msg =   Err%msg // &
                        MODULE_NAME//PROCEDURE_NAME//": Error occurred. The length of the vector delayedRejectionScaleFactorVec (" // &
                        num2str(size(self%Val)) // ") is not equal to delayedRejectionCount = " // &
                        num2str(delayedRejectionCount) // ". If you are not sure how to set the values of &
                        &delayedRejectionScaleFactorVec, drop it from the input. " &
                        // methodName // " will automatically set the appropriate value for delayedRejectionScaleFactorVec.\n\n"
        end if
        do i = 1,size(self%Val)
            if ( self%Val(i)<=0._RK ) then
                Err%occurred = .true.
                Err%msg =   Err%msg // &
                            MODULE_NAME // PROCEDURE_NAME // ": Error occurred. The input value for the element " // &
                            num2str(i) // " of the variable delayedRejectionScaleFactorVec cannot be smaller than or equal to 0.\n\n"
            end if
        end do
        deallocate(DelayedRejectionScaleFactorVec)
    end subroutine checkForSanity

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module SpecDRAM_DelayedRejectionScaleFactorVec_mod ! LCOV_EXCL_LINE

max_stars_repo_name: cdslaborg/paramonte

max_issues_repo_name: cdslaborg/paramonte

max_forks_repo_name: cdslaborg/paramonte


text:
      subroutine cubic (xk0, wp, alph, rad, qplus, qminus)

c     input:  xk0, wp, alph
c     output: rad, qplus, qminus

      implicit double precision (a-h, o-z)
      complex*16 s1,s13
      parameter (three = 3)
      parameter (third = 1/three)

c     this subroutine finds the roots of the equation
c     4xk0 * q^3  +  (alph-4xk0^2) * q^2  +  wp^2 = 0
c     see abramowitz and stegun pg 17 for formulae.

      a2 = (alph / (4*xk0**2)  -  1) * xk0
      a0 = wp**2 / (4*xk0)
      a1 = 0
      q = a1/3 - a2**2/9
      r = (a1*a2 - 3*a0)/6  -  a2**3/27
      rad = q**3 + r**2
      if (rad .gt. 0) then
         qplus = 0
         qminus = 0
         return
      endif

      s13 = dcmplx (r, sqrt(-rad))
      s1 = s13 ** third
      qz1 = 2*s1 - a2/3
c     qz2 = -(s1 + sqrt(three)*dimag(s1) + a2/3)
      qz3 = -(s1 - sqrt(three)*dimag(s1) + a2/3)
      qplus = qz1
      qminus = qz3

      return
      end


text:
! This test code to demonstrate that the function declaration following a "contains" statement are significant.
! It works with gfortran, but fails with OFP.

    PROGRAM  Cubes
       IMPLICIT   NONE

  ! contains statement without an associated function causes parser to fail
    CONTAINS

    END PROGRAM  Cubes

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
great_rsd(k1,x1).
great_rsd(p1,g1).
great_rsd(e1,u1).
great_rsd(v1,l1).
great_rsd(f1,k1).
great_rsd(w1,hh1).
great_rsd(e1,z1).
great_rsd(t1,aa1).
great_rsd(t1,jj1).
great_rsd(e1,jj1).
great_rsd(e1,m1).
great_rsd(w1,k1).
great_rsd(l1,x1).
great_rsd(f1,h1).
great_rsd(p1,q1).
great_rsd(ff1,m1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
 SUBROUTINE elem10(TASK, nelms, elsnam, istop, flag2)

 ! Tasks Routine for element RIGID - HEAT

 USE ctrl_db, ONLY: ndime, ndofn, npoin, lumped
 USE outp_db, ONLY: sumat, iwrit
 USE ele10_db
 USE kinc_db
 USE npo_db

 IMPLICIT NONE
 CHARACTER(len=*),INTENT(IN):: TASK

 ! optional parameters

 LOGICAL, OPTIONAL :: flag2
 CHARACTER (len=*), OPTIONAL  :: elsnam
 INTEGER (kind=4), OPTIONAL  :: istop, nelms(:)

 CHARACTER(len=mnam) :: sname

 INTEGER (kind=4) nelem,nnode,nmast,ntype

 TYPE (ele10_set), POINTER :: elset, anter

 IF ( .NOT.ASSOCIATED (head) ) RETURN

 NULLIFY (anter)
 elset => head

 DO

   CALL comm10(1,nelem,ntype,nnode,nmast,sname,elset)

   SELECT CASE (TRIM(task))

   CASE ('GAUSSC')
     IF( elset%heat ) CALL gaus10(nelem, ntype, elset%ngaus, nnode, elset%lnods, coord,    &
                   elset%cartd, elset%dvolu, elset%shape, istop)

   CASE ('LUMASS')
     IF( ntype == 0 ) THEN   !defined by particles
       IF( elset%tmass /= 0d0 )THEN
       IF( nmast /= 0 ) CALL luma10p (ndime,nnode,elset%lnods,emass,sumat,elset%tmass)

       ELSE
         !nothing, Concentrated masses must be used
       END IF
     ELSE
       IF( nmast /= 0 ) CALL luma10 (nelem,ntype,ndime,nnode,elset%lnods,elset%matno,coord,  &
                                     emass,iwrit,sumat, elset%tmass, lumped,ymass,ifpre)
     END IF

   CASE ('NEW','NSTRA1','NSTRA2')
     CALL acvd10 (nnode,ndime,nelem,ndumn,elset%lnods,nmast,    &
                  ndofn,ifpre,elset%heat)

   CASE ('UPDLON')
     CALL updl10 (elset,oldlb)

   CASE ('RIGBDY')
     !compute center of mass and inertia tensor in principal directions
     IF( ntype == 0 ) THEN   !defined by particles
       CALL rigbdy(nmast,elset%lnods(:,1),nnode)
     ELSE
       IF( nmast /= 0 ) CALL rigb10(nelem,elset%lnods,nmast,nnode)
     END IF

   CASE ('DUMPIN')
     CALL dump10 (sname, nelem, nnode, elset%rbnod, nmast, ntype,   &
                  elset%lnods,elset%matno, elset%tmass,             &
                  elset%heat, elset%ngaus, elset%shape, elset%cartd, elset%dvolu)

   CASE ('WRTPOS')
     CALL mase10(nnode,nelem,elset%matno,elset%lnods,elset%sname,   &
                 elset%ntype, elset%nmast )

   CASE ('DELETE')
     IF( flag2 )EXIT
     flag2 = TRIM(elsnam) == TRIM(sname)
     IF (flag2) THEN
       CALL del_ele10 (head, tail, anter, elset)
       nelms(10) = nelms(10) - 1
       EXIT
     END IF
     IF ( ASSOCIATED (elset%next) ) anter => elset

   CASE ('SURFAC')
     IF( flag2 )EXIT
     flag2 = TRIM(elsnam) == TRIM(sname)
     IF (flag2) THEN
       ! get surface definition from the element set
       CALL surf10 (elset%lnods,elset%nelem,nnode,ndime,ntype)
       EXIT
     END IF

   CASE ('TLMASS')      !compute (lumped) capacity vector
     IF( elset%heat )CALL tlma10 (nelem,ntype,ndime,nnode,elset%ngaus,elset%lnods,elset%dvolu,     &
                  elset%shape,elset%matno,tmass,iftmp)

   CASE ('TRESID')      !compute internal residual heats
     IF( elset%heat )CALL tres10(elset%ngaus,nnode,nelem,ndime,ntype,elset%matno,elset%lnods,   &
                 elset%shape,elset%cartd,elset%dvolu,tresi)

   END SELECT

   IF ( ASSOCIATED (elset%next) ) THEN
     elset => elset%next
   ELSE
     EXIT
   ENDIF

 END DO
 !     delete null elements sets after Nodal Updating
 IF( TRIM(task) == 'UPDLON')THEN
   NULLIFY (anter)
   elset => head
   DO
     IF(.NOT.ASSOCIATED(elset))EXIT
     IF( elset%nelem == 0 )THEN
       CALL del_ele10 (head, tail, anter, elset)
       nelms(10) = nelms(10) - 1
     END IF
     IF (ASSOCIATED(elset))THEN
       anter => elset
       elset => elset%next
     ELSE IF(ASSOCIATED(anter))THEN
       elset => anter%next
     ELSE IF (ASSOCIATED(head)) THEN
       elset => head
     ELSE
       EXIT
     END IF
   END DO
 END IF

 RETURN
 END SUBROUTINE elem10

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
      MODULE vmec_input
      USE vparams, ONLY: rprec, dp, mpol1d, ntord, ndatafmax
      USE vsvd0
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!      dsiobt:   measured flux loop signals corresponding to the
!                combination of signals in iconnect array
!     indxflx:   array giving index of flux measurement in iconnect array
!    indxbfld:   array giving index of bfield measurement used in matching
!    bbc(m,n):   measured magnetic field at rbcoil(m,n),zbcoil(m,n) at
!                the orientation br*cos(abcoil) + bz*sin(abcoil)
!                of the m-th coil in the n-th set from mgrid file.
!       nflxs:   number of flux loop measurements used in matching
!    nbfld(n):   number of selected external bfield measurements in set n from nml file
!    pres_scale: factor used to scale pressure profile (default value = 1)
!                useful so user can fix profile and change beta without having to change
!                all AM coefficients separately
      INTEGER :: nfp, ncurr, nsin, niter, nstep, nvacskip, mpol, ntor,
     1           ntheta, nzeta, ipmass, ipiota, ipcurr
      INTEGER, DIMENSION(100) :: ns_array
      INTEGER :: imse, isnodes, itse, ipnodes, iopt_raxis,
     1   imatch_phiedge, nflxs
      INTEGER, DIMENSION(nbsetsp) :: nbfld
      INTEGER, DIMENSION(nfloops) :: indxflx
      INTEGER, DIMENSION(nbcoilsp,nbsetsp) :: indxbfld
      REAL(rprec), DIMENSION(-ntord:ntord,0:mpol1d) ::
     1   rbs, zbc, rbc, zbs
      REAL(rprec) :: time_slice, curtor, delt, ftol, tcon0,
     1   gamma, phiedge, phidiam, sigma_current, sigma_delphid, tensi,
     2   tensp, tensi2, fpolyi, presfac, mseangle_offset, pres_offset,
     3   mseangle_offsetm, spres_ped, bloat, pres_scale
!     REAL(rprec), DIMENSION(0:10) :: am, ai, ac, aphi
      REAL(rprec), DIMENSION(0:12) :: am, ai, ac, aphi  !change by J.Geiger
      REAL(rprec), DIMENSION(0:ntord) :: raxis, zaxis                !!Backwards compatibility: Obsolete
      REAL(rprec), DIMENSION(0:ntord) :: raxis_cc, raxis_cs,
     1                                   zaxis_cc, zaxis_cs
      REAL(rprec), DIMENSION(100) :: ftol_array
      REAL(rprec), DIMENSION(nigroup) :: extcur
      REAL(rprec), DIMENSION(nmse) :: mseprof
      REAL(rprec), DIMENSION(ntse) :: rthom, datathom, sigma_thom
      REAL(rprec), DIMENSION(nmse) :: rstark, datastark,
     1    sigma_stark
      REAL(rprec), DIMENSION(nfloops) :: dsiobt, sigma_flux
      REAL(rprec), DIMENSION(nbcoilsp,nbsetsp) :: bbc, sigma_b
      REAL(rprec), DIMENSION(ndatafmax) :: psa, pfa, isa, ifa
      LOGICAL :: lpofr, lmac, lfreeb, lrecon, loldout, ledge_dump,
     1           lasym, laddout, ldiagno, lmoreiter
     2         , lspectrum_dump, loptim           !!Obsolete
      CHARACTER*(100) :: mgrid_file
      CHARACTER*(120) :: arg1
      CHARACTER*(100) :: input_extension

      NAMELIST /indata/ mgrid_file, time_slice, nfp, ncurr, nsin,
     1   niter, nstep, nvacskip, delt, ftol, gamma, am, ai, ac, aphi,
     2   rbc, zbs, rbs, zbc, spres_ped, pres_scale, raxis_cc, zaxis_cs,
     3   raxis_cs, zaxis_cc, mpol, ntor, ntheta, nzeta,
     4   ns_array, ftol_array, tcon0, curtor, sigma_current, extcur,
     5   phiedge, psa, pfa, isa, ifa, imatch_phiedge, iopt_raxis,
     6   tensi, tensp, mseangle_offset, mseangle_offsetm, imse,
     7   isnodes, rstark, datastark, sigma_stark, itse, ipnodes,
     8   presfac, pres_offset, rthom, datathom, sigma_thom, phidiam,
     9   sigma_delphid, tensi2, fpolyi, nflxs, indxflx, dsiobt,
     A   sigma_flux, nbfld, indxbfld, laddout, ldiagno, lmoreiter,
     A   bbc, sigma_b, lpofr, lfreeb, lrecon, lmac, loldout, lasym,
     B   ipmass, ipiota, ipcurr,
     B   ledge_dump, lspectrum_dump, loptim, bloat, raxis, zaxis

      NAMELIST /mseprofile/ mseprof

      CONTAINS

      SUBROUTINE read_indata_namelist (iunit, istat)
      REAL(rprec), PARAMETER :: smallno = -HUGE(1._dp)
      INTEGER :: iunit, istat

!
!     BACKWARDS COMPATIBILITY
!
      raxis = smallno
      zaxis = smallno

      READ (iunit, nml=indata, iostat=istat)

      WHERE (raxis .gt. smallno) raxis_cc = raxis
      WHERE (zaxis .gt. smallno) zaxis_cs = zaxis

      END SUBROUTINE read_indata_namelist

      SUBROUTINE read_mse_namelist (iunit, istat)
      INTEGER :: iunit, istat

      READ (iunit, nml=mseprofile, iostat=istat)

      END SUBROUTINE read_mse_namelist

      END MODULE vmec_input


text:
        SUBROUTINE DG_NTIM  ( chngnv, coladd, time, nxtok, iret )
C************************************************************************
C* DG_NTIM                                                              *
C*                                                                      *
C* This subroutine initialize processing for the next time.             *
C*                                                                      *
C* DG_NTIM  ( CHNGNV, COLADD, TIME, NXTOK, IRET )			*
C*                                                                      *
C* Input parameters:                                                    *
C*      CHNGNV          LOGICAL         Flag to change navigation       *
C*	COLADD		LOGICAL		Flag to add a column of data	*
C*                                                                      *
C* Output parameters:                                                   *
C*	TIME (2)	CHAR*		Next time from dtmlst		*
C*	NXTOK		LOGICAL		Flag for existance of next time	*
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C*					-50 = navigation not same	*
C*					-57 = grid file open failed	*
C**                                                                     *
C* Log:                                                                 *
C* R. Tian/SAIC		 3/06	Fortran wrapper of DGC_NTIM		*
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
	LOGICAL		chngnv, coladd, nxtok
        CHARACTER*(*)  	time (2)
C------------------------------------------------------------------------
	CALL DGC_NTIM ( chngnv, coladd, time (1), time (2), nxtok, iret )
	CALL ST_RNUL ( time (1), time (1), lens, ier )
	CALL ST_RNUL ( time (2), time (2), lens, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module convert
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
! <doc>Convert from complex variable a(d1,d2,d3, ...) to a
! real variable ar(2,d1,d2,d3,...) and back.
! This is necessary for saving complex variables in NetCDF format</doc>
!
!     (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland,
!     P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED.
!
   implicit none
   private
   public :: c2r, r2c

   interface c2r
      module procedure x1c2r
      module procedure x2c2r
      module procedure x3c2r
      module procedure x4c2r
      module procedure x5c2r
      module procedure x6c2r
   end interface

   interface r2c
      module procedure x1r2c
      module procedure x2r2c
      module procedure x3r2c
      module procedure x4r2c
      module procedure x5r2c
   end interface

contains
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x5c2r(a, a_ri)

      complex, dimension(:, :, :, :, :), intent(in) :: a
      real, dimension(:, :, :, :, :, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x5c2r: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x5c2r: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x5c2r: size(a, 3) does not match size(a_ri, 4)')
      if (size(a, 4) /= size(a_ri, 5)) call aborter(6, 'x5c2r: size(a, 4) does not match size(a_ri, 5)')
      if (size(a, 5) /= size(a_ri, 6)) call aborter(6, 'x5c2r: size(a, 5) does not match size(a_ri, 6)')
      a_ri(1, :, :, :, :, :) = real(a(:, :, :, :, :))
      a_ri(2, :, :, :, :, :) = aimag(a(:, :, :, :, :))

   end subroutine x5c2r

   subroutine x6c2r(a, a_ri)

      complex, dimension(:, :, :, :, :, :), intent(in) :: a
      real, dimension(:, :, :, :, :, :, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x6c2r: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x6c2r: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x6c2r: size(a, 3) does not match size(a_ri, 4)')
      if (size(a, 4) /= size(a_ri, 5)) call aborter(6, 'x6c2r: size(a, 4) does not match size(a_ri, 5)')
      if (size(a, 5) /= size(a_ri, 6)) call aborter(6, 'x6c2r: size(a, 5) does not match size(a_ri, 6)')
      if (size(a, 6) /= size(a_ri, 7)) call aborter(6, 'x6c2r: size(a, 6) does not match size(a_ri, 7)')
      a_ri(1, :, :, :, :, :, :) = real(a(:, :, :, :, :, :))
      a_ri(2, :, :, :, :, :, :) = aimag(a(:, :, :, :, :, :))

   end subroutine x6c2r
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x5r2c(a, a_ri)

      real, dimension(:, :, :, :, :, :), intent(in) :: a_ri
      complex, dimension(:, :, :, :, :), intent(out) :: a

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x5r2c: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x5r2c: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x5r2c: size(a, 3) does not match size(a_ri, 4)')
      if (size(a, 4) /= size(a_ri, 5)) call aborter(6, 'x5r2c: size(a, 4) does not match size(a_ri, 5)')
      if (size(a, 5) /= size(a_ri, 6)) call aborter(6, 'x5r2c: size(a, 5) does not match size(a_ri, 6)')
      a(:, :, :, :, :) = cmplx(a_ri(1, :, :, :, :, :), a_ri(2, :, :, :, :, :))

   end subroutine x5r2c
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x4c2r(a, a_ri)

      complex, dimension(:, :, :, :), intent(in) :: a
      real, dimension(:, :, :, :, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x4c2r: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x4c2r: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x4c2r: size(a, 3) does not match size(a_ri, 4)')
      if (size(a, 4) /= size(a_ri, 5)) call aborter(6, 'x4c2r: size(a, 4) does not match size(a_ri, 5)')
      a_ri(1, :, :, :, :) = real(a(:, :, :, :))
      a_ri(2, :, :, :, :) = aimag(a(:, :, :, :))

   end subroutine x4c2r
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x4r2c(a, a_ri)

      real, dimension(:, :, :, :, :), intent(in) :: a_ri
      complex, dimension(:, :, :, :), intent(out) :: a

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x4r2c: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x4r2c: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x4r2c: size(a, 3) does not match size(a_ri, 4)')
      if (size(a, 4) /= size(a_ri, 5)) call aborter(6, 'x4r2c: size(a, 4) does not match size(a_ri, 5)')
      a(:, :, :, :) = cmplx(a_ri(1, :, :, :, :), a_ri(2, :, :, :, :))

   end subroutine x4r2c
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x3c2r(a, a_ri)

      complex, dimension(:, :, :), intent(in) :: a
      real, dimension(:, :, :, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x3c2r: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x3c2r: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x3c2r: size(a, 3) does not match size(a_ri, 4)')
      a_ri(1, :, :, :) = real(a(:, :, :))
      a_ri(2, :, :, :) = aimag(a(:, :, :))

   end subroutine x3c2r
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x3r2c(a, a_ri)

      real, dimension(:, :, :, :), intent(in) :: a_ri
      complex, dimension(:, :, :), intent(out) :: a

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x3r2c: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x3r2c: size(a, 2) does not match size(a_ri, 3)')
      if (size(a, 3) /= size(a_ri, 4)) call aborter(6, 'x3r2c: size(a, 3) does not match size(a_ri, 4)')
      a(:, :, :) = cmplx(a_ri(1, :, :, :), a_ri(2, :, :, :))

   end subroutine x3r2c
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x2c2r(a, a_ri)

      complex, dimension(:, :), intent(in) :: a
      real, dimension(:, :, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x2c2r: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x2c2r: size(a, 2) does not match size(a_ri, 3)')
      a_ri(1, :, :) = real(a(:, :))
      a_ri(2, :, :) = aimag(a(:, :))

   end subroutine x2c2r
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x2r2c(a, a_ri)

      real, dimension(:, :, :), intent(in) :: a_ri
      complex, dimension(:, :), intent(out) :: a

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x2r2c: size(a, 1) does not match size(a_ri, 2)')
      if (size(a, 2) /= size(a_ri, 3)) call aborter(6, 'x2r2c: size(a, 2) does not match size(a_ri, 3)')
      a(:, :) = cmplx(a_ri(1, :, :), a_ri(2, :, :))

   end subroutine x2r2c
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x1c2r(a, a_ri)

      complex, dimension(:), intent(in) :: a
      real, dimension(:, :), intent(out) :: a_ri

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x2c2r: size(a, 1) does not match size(a_ri, 2)')
      a_ri(1, :) = real(a(:))
      a_ri(2, :) = aimag(a(:))

   end subroutine x1c2r
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   subroutine x1r2c(a, a_ri)

      real, dimension(:, :), intent(in) :: a_ri
      complex, dimension(:), intent(out) :: a

      if (size(a, 1) /= size(a_ri, 2)) call aborter(6, 'x2r2c: size(a, 1) does not match size(a_ri, 2)')
      a(:) = cmplx(a_ri(1, :), a_ri(2, :))

   end subroutine x1r2c
!------------------------------------------------------------------------------
!                              AstroGK, 2009
!------------------------------------------------------------------------------
!
   Subroutine Aborter(iunit, ierrmsg)
!----------------------------------------------------------------------
!  ABORT A PROGRAM AFTER A FATAL ERROR CONDITION IS DETECTED.
!
!
! input: iunit  Unit Number of the file to write error messages to.
!        ierrmsg  An error message to write ilunerr
!
!       The advantage of using this subroutine is that it will
!       generate a traceback showing the chain of subroutines which
!       eventually bombed, and it forces an arithmetic error which
!       the job control system can detect as a fatal error.
!
      character ierrmsg * (*)
      real :: zz0, zz1
      integer :: iunit, ilunerr
      common / abortcmn / zz0, zz1
!
! zz0 is in a common block to prevent an optimizing compiler
! from evaluating 1.0/zz0 during compilation rather than during
! execution.
!
      write (iunit, 1001)
1001  format(//' %ABORTER:  ** FATAL ERROR.  ABORT SUBROUTINE CALLED **'//)

      write (iunit, 1002) ierrmsg
1002  format(1x, a, //)

! on CRAY's and VAXes:
! generate a divide-by-zero error:
      zz1 = 1.0 / zz0         !^
      ilunerr = 5           !^

! on the DecStation:
!^      call exit(1)

      write (ilunerr, 1010) zz0
1010  format(' ?ABORTER-- ZZ0= ', 1PE11.4, ' AND STILL EXECUTING...')

   end subroutine aborter
!------------------------------------------------------------------------------
end module convert

max_stars_repo_name: AntonioG-Jerez/stella

max_issues_repo_name: AntonioG-Jerez/stella

max_forks_repo_name: AntonioG-Jerez/stella


text:
module local_routines
    !!
    !! Setup the Merewether problem.
    !!

    use global_mod, only: dp, ip, force_double, charlen, wall_elevation, pi
    use domain_mod, only: domain_type, STG, UH, VH, ELV
    use read_raster_mod, only: read_gdal_raster
    use which_mod, only: which
    use ragged_array_mod, only: ragged_array_2d_ip_type
    use file_io_mod, only: read_character_file, read_csv_into_array
    use points_in_poly_mod, only: points_in_poly
    use iso_c_binding, only: c_f_pointer, c_loc
    implicit none

    ! Add rain
    ! Discharge of 19.7 m^3/s occurs over a cirle with radius 15
    real(dp), parameter :: rain_centre(2) = [382300.0_dp, 6354290.0_dp], rain_radius = 15.0_dp
    real(force_double) :: rain_rate =  real(19.7_dp, force_double)/(pi * rain_radius**2)
    integer(ip) :: num_input_discharge_cells

    contains

    !
    ! Main setup routine
    !
    subroutine set_initial_conditions_merewether(domain, reflective_boundaries)

        class(domain_type), target, intent(inout):: domain
        logical, intent(in):: reflective_boundaries

        integer(ip):: i, j, k
        real(dp), allocatable:: x(:,:), y(:,:)

        character(len=charlen):: input_elevation_file, polygon_filename

        ! Add this elevation to all points in houses/ csv polygons
        real(dp), parameter:: house_height = 3.0_dp
        ! friction parameters
        real(dp), parameter:: friction_road = 0.02_dp, friction_other = 0.04_dp

        ! things to help read polygons
        character(len=charlen), allocatable:: house_filenames(:)
        integer(ip):: house_file_unit, inside_point_counter, house_cell_count
        real(dp), allocatable:: polygon_coords(:,:)
        logical, allocatable:: is_inside_poly(:)
        type(ragged_array_2d_ip_type), pointer :: rainfall_region_indices
        integer(ip), allocatable :: i_inside(:)

        allocate(x(domain%nx(1),domain%nx(2)), y(domain%nx(1),domain%nx(2)))

        !
        ! Dry flow to start with. Later we clip stage to elevation.
        !
        domain%U(:,:,[STG, UH, VH]) = 0.0_dp

        !
        ! Set elevation with the raster
        !
        input_elevation_file = "./topography/topography1.tif"

        do j = 1, domain%nx(2)
          do i = 1, domain%nx(1)
            x(i,j) = domain%lower_left(1) + (i-0.5_dp)*domain%dx(1)
            y(i,j) = domain%lower_left(2) + (j-0.5_dp)*domain%dx(2)
          end dO
        end do
        call read_gdal_raster(input_elevation_file, x, y, domain%U(:,:,ELV), &
            domain%nx(1)*domain%nx(2), verbose=1_ip, bilinear=0_ip)

        print*, 'Elevation range: ', minval(domain%U(:,:,ELV)), maxval(domain%U(:,:,ELV))

        ! Get filenames for the houses
        open(newunit = house_file_unit, file='houses_filenames.txt')
        call read_character_file(house_file_unit, house_filenames, "(A)")
        close(house_file_unit)

        ! Read the houses and add house_height to the elevation for all points inside
        allocate(is_inside_poly(domain%nx(1)))
        house_cell_count = 0
        do k = 1, size(house_filenames)

            inside_point_counter = 0
            polygon_filename = './' // TRIM(house_filenames(k))
            call read_csv_into_array(polygon_coords, polygon_filename)

            do j = 1, domain%nx(2)
                ! Find points in the polygon in column j
                call points_in_poly(polygon_coords(1,:), polygon_coords(2,:), x(:,j), y(:,j), is_inside_poly)
                inside_point_counter = inside_point_counter + count(is_inside_poly)

                do i = 1, domain%nx(1)
                    if(is_inside_poly(i)) then
                        domain%U(i,j,ELV) = domain%U(i,j,ELV) + house_height
                        house_cell_count = house_cell_count + 1_ip
                    end if
                end do

            end do

            print*, trim(house_filenames(k)), ': ', inside_point_counter

        end do

        print*, '# cells in houses: ', house_cell_count

        if(reflective_boundaries) then
            ! Walls all sides
            domain%U(1,:,ELV) = wall_elevation
            domain%U(domain%nx(1), :, ELV) = wall_elevation
            domain%U(:, 1 ,ELV) = wall_elevation
            domain%U(:, domain%nx(2), ELV) = wall_elevation
        else
            ! Transmissive outflow -- but we prevent mass leaking out of the back of the domain.
            ! Also block off other sides -- only need outflow on east edge of domain.
            domain%U(:, 1:2, ELV) = domain%U(:,1:2, ELV) + 10.0_dp
            domain%U(1:2, :, ELV) = domain%U(1:2, :, ELV) + 10.0_dp
            domain%U(:, (domain%nx(2)-1):domain%nx(2), ELV) = 10.0_dp + &
                domain%U(:, (domain%nx(2)-1):domain%nx(2), ELV)
        end if

        !
        ! Set friction
        !

        ! Initial value -- later updated based on roads
        domain%manning_squared = friction_other * friction_other

        ! Find points in road polygon and set friction there
        polygon_filename = 'Road/RoadPolygon.csv'
        call read_csv_into_array(polygon_coords, polygon_filename)
        inside_point_counter = 0
        do j = 1, domain%nx(2)

            call points_in_poly(polygon_coords(1,:), polygon_coords(2,:), x(:,j), y(:,j), is_inside_poly)
            inside_point_counter = inside_point_counter + count(is_inside_poly)

            do i = 1, domain%nx(1)
                if(is_inside_poly(i)) domain%manning_squared(i,j) = friction_road * friction_road
            end do

        end do

        print*, ''
        print*, '# Points in road polygon :', inside_point_counter
        print*, ''

        ! Ensure stage >= elevation
        domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV) + 1.0e-08_dp)

        print*, 'Elevation range: ', minval(domain%U(:,:,ELV)), maxval(domain%U(:,:,ELV))
        print*, 'Stage range: ', minval(domain%U(:,:,STG)), maxval(domain%U(:,:,STG))

        ! Figure out indices that are inside the rainfall forcing region.
        allocate(rainfall_region_indices)
        allocate(rainfall_region_indices%i2(domain%nx(2)))
        num_input_discharge_cells = 0
        do j = 1, domain%nx(2)
            ! Find cells in this row that are within the rain circle
            call which( (domain%x - rain_centre(1))**2 < rain_radius**2 - (domain%y(j) - rain_centre(2))**2, &
                rainfall_region_indices%i2(j)%i1)

            ! For this routine we cross-check conservation by counting the number of input discharge cells and
            ! doing a separate mass balance. This assumes 1 domain (only)
            num_input_discharge_cells = num_input_discharge_cells + size(rainfall_region_indices%i2(j)%i1)
        end do
        domain%forcing_context_cptr = c_loc(rainfall_region_indices)

    end subroutine

    !
    ! This is the discharge source term, called every time-step. It's like rainfall
    ! in a "circle"
    !
    subroutine apply_rainfall_forcing(domain, dt)
        type(domain_type), intent(inout) :: domain
        real(dp), intent(in) :: dt

        integer(ip) :: j
        type(ragged_array_2d_ip_type), pointer :: rainfall_region_indices

        ! Unpack the forcing context pointer. In this case it is a ragged array giving the indices
        ! where we should apply the rainfall in this domain
        call c_f_pointer(domain%forcing_context_cptr, rainfall_region_indices)

        !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(domain, rain_rate, dt, rainfall_region_indices)
        !$OMP DO SCHEDULE(STATIC)
        do j = 1, domain%nx(2)
            if(size(rainfall_region_indices%i2(j)%i1) > 0) then
                ! Rainfall at cells within a circle of radius "rain_radius" about "rain_centre"
                domain%U(    rainfall_region_indices%i2(j)%i1, j,STG) = rain_rate * dt + &
                    domain%U(rainfall_region_indices%i2(j)%i1, j,STG)
            end if
        end do
        !$OMP END DO
        !$OMP END PARALLEL

    end subroutine

end module

!@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program merewether
    !! Urban flooding test case in Merewether from Australian Rainfall and Runoff.
    !! Smith, G. & Wasko, C. Revision Project 15: Two Dimensional Simulations In
    !! Urban Areas - Representation of Buildings in 2D Numerical Flood Models Australian
    !! Rainfall and Runoff, Engineers Australia, 2012

    use global_mod, only: ip, dp, minimum_allowed_depth
    use domain_mod, only: domain_type
    use boundary_mod, only: transmissive_boundary
    use local_routines

    implicit none

    integer(ip):: j
    real(dp):: last_write_time
    type(domain_type):: domain

    ! Approx timestep between outputs
    real(dp), parameter :: approximate_writeout_frequency = 60.0_dp
    real(dp), parameter :: final_time = 900.0_dp

    ! Length/width
    real(dp), parameter, dimension(2):: global_lw = [320.0_dp, 415.0_dp]
    ! Lower-left corner coordinate
    real(dp), parameter, dimension(2):: global_ll = [382251.0_dp, 6354266.5_dp]
    ! grid size (number of x/y cells)
    integer(ip), parameter, dimension(2):: global_nx = [320_ip, 415_ip] ![160_ip, 208_ip] ![321_ip, 416_ip]
    ! Use reflective or transmissive boundaries
    logical, parameter:: reflective_boundaries = .FALSE.
    ! Track mass
    real(force_double) :: volume_added, volume_initial
    ! Use a small time step at the start
    real(dp) :: startup_timestep = 0.05_dp
    real(dp) :: startup_time = 10.0_dp

    domain%timestepping_method = 'rk2n' !'euler' !'rk2n'

    ! Use a very small maximum timestep initially, because the domain is dry, but water is quickly
    ! added, which can cause 'first-step' instability otherwise. Re-set it after an evolve
    call domain%allocate_quantities(global_lw, global_nx, global_ll)

    if((.not. reflective_boundaries) .and. (.not. all(domain%is_nesting_boundary))) then
        ! Allow waves to propagate outside all edges
        domain%boundary_subroutine => transmissive_boundary
    end if
    ! Add rainfall to the domain
    domain%forcing_subroutine => apply_rainfall_forcing

    ! Call local routine to set initial conditions
    call set_initial_conditions_merewether(domain, reflective_boundaries)
    call domain%update_boundary()
    volume_initial = domain%mass_balance_interior()

    ! Trick to get the code to write out just after the first timestep
    last_write_time = -approximate_writeout_frequency

    ! Evolve the code
    do while (.TRUE.)

        if(domain%time - last_write_time >= approximate_writeout_frequency) then

            last_write_time = last_write_time + approximate_writeout_frequency

            call domain%print()
            call domain%write_to_output_files()

            print*, 'mass_balance: ', domain%mass_balance_interior()
            print*, 'volume_added + initial: ', volume_added + volume_initial
            print*, '.... difference: ', domain%mass_balance_interior() - (volume_added + volume_initial)
        end if

        if (domain%time > final_time) then
            exit
        end if

        ! At the start, evolve with a specified (small) timestep.
        ! Later allow an adaptive step
        if(domain%time < startup_time) then
            call domain%evolve_one_step(startup_timestep)
        else
            call domain%evolve_one_step()
        end if

        ! Track volume added 'in theory' based on the known rainfall rate that is going into
        ! a known number of cells
        volume_added = domain%time * rain_rate * num_input_discharge_cells * product(domain%dx)

        !print*, 'negative_depth_fix_counter: ', domain%negative_depth_fix_counter

    end do

    call domain%write_max_quantities()

    ! Print timing info
    call domain%timer%print()

    ! Simple test (mass conservation only)
    print*, 'MASS CONSERVATION TEST'
    if(abs(domain%mass_balance_interior() - (volume_added + volume_initial) ) < 1.0e-06_dp) then
        print*, 'PASS'
    else
        print*, 'FAIL -- mass conservation not good enough. ', &
            'This is expected with single precision, which will affect the computed direct rainfall volumes'
    end if

    call domain%finalise()
END PROGRAM

max_stars_repo_name: KayaGeoAus/ptha

max_issues_repo_name: KayaGeoAus/ptha

max_forks_repo_name: KayaGeoAus/ptha


text:
	SUBROUTINE getgrid(Nx,Ny,Nz,Lx,Ly,Lz,pi,name_config,x,y,z,kx,ky,kz)
	!--------------------------------------------------------------------
	!
	!
	! PURPOSE
	!
	! This subroutine gets grid points and fourier frequencies for a
	! pseudospectral simulation of the 2D nonlinear Klein-Gordon equation
	!
	! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*u^3
	!
	! The boundary conditions are u(x=-Lx*pi,y,z)=u(x=Lx*\pi,y,z),
	!	u(x,y=-Ly*pi,z)=u(x,y=Ly*pi,z),u(x,y,z=-Ly*pi)=u(x,y,z=Ly*pi),
	!
	! INPUT
	!
	! .. Scalars ..
	!  Nx				= number of modes in x - power of 2 for FFT
	!  Ny				= number of modes in y - power of 2 for FFT
	!  Ny				= number of modes in z - power of 2 for FFT
	!  pi				= 3.142....
	!  Lx				= width of box in x direction
	!  Ly				= width of box in y direction
	!  Lz				= width of box in z direction
	! OUTPUT
	!
	! .. Vectors ..
	!  kx				= fourier frequencies in x direction
	!  ky				= fourier frequencies in y direction
	!  kz				= fourier frequencies in z direction
	!  x				= x locations
	!  y				= y locations
	!  z				= z locations
	!
	! LOCAL VARIABLES
	!
	! .. Scalars ..
	!  i				= loop counter in x direction
	!  j				= loop counter in y direction
	!  k				= loop counter in z direction
	!
	! REFERENCES
	!
	! ACKNOWLEDGEMENTS
	!
	! ACCURACY
	!
	! ERROR INDICATORS AND WARNINGS
	!
	! FURTHER COMMENTS
	! Check that the initial iterate is consistent with the
	! boundary conditions for the domain specified
	!--------------------------------------------------------------------
	! External routines required
	!
	! External libraries required
	! OpenMP library
	USE omp_lib
	IMPLICIT NONE
	! Declare variables
	INTEGER(KIND=4), INTENT(IN)							:: Nx,Ny,Nz
	REAL(kind=8), INTENT(IN)							:: Lx,Ly,Lz,pi
	REAL(KIND=8), DIMENSION(1:NX), INTENT(OUT) 			:: x
	REAL(KIND=8), DIMENSION(1:NY), INTENT(OUT) 			:: y
	REAL(KIND=8), DIMENSION(1:NZ), INTENT(OUT) 			:: z
	COMPLEX(KIND=8), DIMENSION(1:NX), INTENT(OUT)		:: kx
	COMPLEX(KIND=8), DIMENSION(1:NY), INTENT(OUT)		:: ky
	COMPLEX(KIND=8), DIMENSION(1:NZ), INTENT(OUT)		:: kz
	CHARACTER*100, INTENT(OUT)							:: name_config
	INTEGER(kind=4)										:: i,j,k
	!$OMP PARALLEL
	!$OMP DO SCHEDULE(static)
	DO i=1,1+Nx/2
		kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx
	END DO
	!$OMP END DO
	kx(1+Nx/2)=0.0d0
	!$OMP DO SCHEDULE(static)
	DO i = 1,Nx/2 -1
		kx(i+1+Nx/2)=-kx(1-i+Nx/2)
	END DO
	!$OMP END DO
	!$OMP DO SCHEDULE(static)
	DO i=1,Nx
		x(i)=(-1.0d0 + 2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)))*pi*Lx
	END DO
	!$OMP END DO

	!$OMP DO SCHEDULE(static)
	DO j=1,1+Ny/2
		ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly
	END DO
	!$OMP END DO
	ky(1+Ny/2)=0.0d0
	!$OMP DO SCHEDULE(static)
	DO j = 1,Ny/2 -1
		ky(j+1+Ny/2)=-ky(1-j+Ny/2)
	END DO
	!$OMP END DO
	!$OMP DO SCHEDULE(static)
	DO j=1,Ny
		y(j)=(-1.0d0 + 2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)))*pi*Ly
	END DO
	!$OMP END DO

	!$OMP DO SCHEDULE(static)
	DO k=1,1+Nz/2
		kz(k)= cmplx(0.0d0,1.0d0)*REAL(k-1,kind(0d0))/Lz
	END DO
	!$OMP END DO
	kz(1+Nz/2)=0.0d0
	!$OMP DO SCHEDULE(static)
	DO k = 1,Nz/2 -1
		kz(k+1+Nz/2)=-kz(1-k+Nz/2)
	END DO
	!$OMP END DO
	!$OMP DO SCHEDULE(static)
	DO k=1,Nz
		z(k)=(-1.0d0 + 2.0d0*REAL(k-1,kind(0d0))/REAL(Nz,kind(0d0)))*pi*Lz
	END DO
	!$OMP END DO
	!$OMP END PARALLEL
	! Save x grid points in text format
	name_config = 'xcoord.dat'
	OPEN(unit=11,FILE=name_config,status="UNKNOWN")
	REWIND(11)
	DO i=1,Nx
		WRITE(11,*) x(i)
	END DO
	CLOSE(11)
	! Save y grid points in text format
	name_config = 'ycoord.dat'
	OPEN(unit=11,FILE=name_config,status="UNKNOWN")
	REWIND(11)
	DO j=1,Ny
		WRITE(11,*) y(j)
	END DO
	CLOSE(11)
	! Save z grid points in text format
	name_config = 'zcoord.dat'
	OPEN(unit=11,FILE=name_config,status="UNKNOWN")
	REWIND(11)
	DO k=1,Nz
		WRITE(11,*) z(k)
	END DO
	CLOSE(11)

	END SUBROUTINE getgrid

max_stars_repo_name: bcloutier/PSNM

max_issues_repo_name: bcloutier/PSNM

max_forks_repo_name: bcloutier/PSNM


text:
program test

  real :: a(3)
  real :: sq = 9

  a(1) = 1.0
  a(2) = 2.0
  a(3) = 4.0

  a = sqrt(sq) * a
  print *, a

end program

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
      logical function interior_edge(itin,iein,
     *  ipitetoff,ipjtetoff,ipitet,ipjtet,ipitettyp,ipiparent,nef,
     *  mbndry)
C
C #####################################################################
C
C     PURPOSE -
C
C        This routine finds returns true if the edge is an interior
C        edge false if it is a boundary edge
C
C     INPUT ARGUMENTS -
C
C       it   -- seed element
C       ie   -- local edge # of it that is edge to be tested
C       (ip)itetoff (ip)jtetoff (ip)itet (ip)jtet (ip)itettyp
c            standard mesh object pointers
C       (ip)iparent - Pointer to array of parent points of all nodes
C       nef - number of element faces for this cmo
C
C     OUTPUT ARGUMENTS -
C        nelts - number of elements containing the edge
C        (ip)elts - elements containing edge
C        (ip)edges _ local edge numbers
C
C     CHANGE HISTORY -
C
C        $Log: interior_edge.f,v $
C        Revision 2.00  2007/11/05 19:45:59  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.21   02 Oct 2007 12:40:28   spchu
CPVCS    original version
C
C ######################################################################
C
      implicit none
C
      character*132 logmess
C
      include "local_element.h"
C
C ######################################################################
C
C
      integer nelts,itin,iein
      pointer (ipitetoff, itetoff)
      pointer (ipjtetoff, jtetoff)
      pointer (ipitettyp, itettyp)
      integer  itetoff(*), jtetoff(*),itettyp(*)
C
      pointer (ipitet, itet1)
      pointer (ipjtet, jtet1)
      integer itet1(*), jtet1(*)
      pointer (ipiparent, iparent)
      integer iparent(*)
C
      integer icscode,iip1,iip2,
     *        ie,ip1,ip2,nef,
     *        i1,i2,itstart,itt,itold,itnew,
     *        iff,ifold,ifnew,j,mbndry,
     *        ittyp,nf,
     *        jt,idir,nfstart,nf1
      character*32 isubname
C
C#######################################################################
C
C
      isubname='interior_edge'
C
C  nodes on edge (itin,iein)
C
      iip1=itet1(itetoff(itin)+ielmedge1(1,iein,itettyp(itin)))
      iip2=itet1(itetoff(itin)+ielmedge1(2,iein,itettyp(itin)))
      ip1=iparent(iip1)
      ip2=iparent(iip2)
C
C  Loop through all neigbors of element to find one that shares
C  this edge
C
      interior_edge=.true.
      nelts = 1
      itstart=itin
      nfstart=1
      itold=0
      ifold=0
      itt=itin
      idir=0
      ittyp=itettyp(itt)
      nf1=1
 10   do nf =nf1,nelmnef(ittyp)
           do j=1,ielmface0(nf,ittyp)
             ie=ielmface2(j,nf,ittyp)
             i1=itet1(itetoff(itt)+ielmedge1(1,ie,ittyp))
             i2=itet1(itetoff(itt)+ielmedge1(2,ie,ittyp))
             i1=iparent(i1)
             i2=iparent(i2)
             if(i1.eq.ip1.and.i2.eq.ip2.or.
     *         i2.eq.ip1.and.i1.eq.ip2) then
               iff=nf
               if(itt.eq.itstart) nfstart=nf
C
C  check if external boundary face
C  make sure to check both directions for material
C  interfaces (use idir)
C
               jt=jtet1(jtetoff(itt)+iff)
               if(jt.eq.mbndry) then
                  interior_edge=.false.
                  go to 9999
               endif
               nf1=1
C
C  check if interface face - if so both edges
C  and node to be added are on an internal interface
C
               if(jt.gt.mbndry)then
                  jt=jt-mbndry
               endif
C
C  check that we didn't go backwards
C
               itnew=1+(jt-1)/nef
               ifnew=jt-nef*(itnew-1)
               if (itnew.ne.itold.or.ifnew.ne.ifold) then
C
C  see if back to starting element
C
                  if( itnew.eq.itstart) then
                      go to 9999
                  endif
C
C  found next element in chain
C
                  ittyp=itettyp(itnew)
                  ittyp=itettyp(itnew)
                  itold=itt
                  ifold=iff
C
C  still looking for more tets surrounding the edge
C
               itt=itnew
               go to 10
               endif
             endif
C end loop on edges
           enddo
C end loop on faces
      enddo
C
C  shouldn't ever get here
C
  19  write(logmess, 20)itin,iein
  20  format (' ATTENTION: error in interior_edge',
     *         i10,' edge ',i10 )
      call writloga('default',0,logmess,0,icscode)
C
9999  continue
      return
      end

max_forks_repo_name: daniellivingston/LaGriT