text:
      subroutine rdd2l(nfldmx,nfldt,nfldtx,ninpts,nlchmx,
     $ nttyo,qrderr,ufield,uline1,uline2,ulscr)
c
c     This subroutine reads a "2 line" block containing an unknown
c     header or tag string from an EQ3/6 input in menu-style ("D")
c     format. The first line contains the header and data. The second
c     line is a separator line containing a single field of dashes.
c
c     EQLIBU/rdd2lh.f is very similar to the present subroutine. It
c     differs in that the header is required to match a defined string.
c
c     This subroutine is called by:
c
c       XCON3/rd3d7.f
c       XCON6/rd6d7.f
c       XCON6/rd6d8.f
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       nfldmx = the maximum number of fields on line uline1 or
c                  line uline2
c       nfldtx = the expected number of fields on a line
c       nlchmx = the maximum number of characters on a line
c       ninpts = the unit number of the stripped input file
c       nttyo  = the unit number of the screen file
c       ulscr  = scratch character variable
c
c     Output:
c
c       nfldt  = the number of fields on line uline1 or line uline2
c       qrderr = logical flag, true if an error was encountered
c       ufield = the array of fields on a line
c       uline1 = the first line to be read
c       uline2 = the second line to be read
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nttyo
c
      integer nfldmx,nfldt,nfldtx,ninpts,nlchmx
c
      logical qrderr
c
      character*(*) ufield(nfldmx)
      character*(*) uline1,uline2,ulscr
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer j2
c
      integer ilnobl
c
c-----------------------------------------------------------------------
c
      qrderr = .false.
c
c     Read and parse the first line.
c
      read (ninpts,1000,err=990) uline1
 1000 format(a80)
c
      call parslj(nfldmx,nfldt,nlchmx,ufield,uline1,ulscr)
c
c     Compare the number of fields found with that expected.
c
      if (nfldtx .gt. 0) then
        if (nfldt .ne. nfldtx) then
          j2 = ilnobl(uline1(1:70))
          write (nttyo,1010) nfldt,nfldtx,uline1(1:j2)
 1010     format(/' * Warning - (EQLIBU/rdd2l) Found ',i2,' fields',
     $    /7x,'where ',i2,' were expected on the line beginning with',
     $    /7x,'"',a,'".')
        endif
      endif
c
c     Read the second line.
c
      read (ninpts,1000,err=990) uline2
      if (uline2(2:9) .eq. '--------') go to 999
c
      qrderr = .true.
      j2 = ilnobl(uline2(1:70))
      write (nttyo,1030) uline2(1:j2)
 1030 format(/' * Error - (EQLIBU/rdd2l) Found the line beginning',
     $ ' with',/7x,'"',a,'"',
     $ /7x,'where a dashed separator line was expected.')
      go to 999
c
  990 qrderr = .true.
c
  999 continue
      end


text:
!!! FROM sort.f90 (https://www.mjr19.org.uk/IT/sorts/sorts.f90)

MODULE MODULE_SORT

CONTAINS

RECURSIVE SUBROUTINE QUICKSORT(array)

    IMPLICIT NONE

    DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT)::array
    DOUBLE PRECISION :: temp,pivot
    INTEGER :: ii, jj, last, left, right

    last=SIZE(array)

    IF (last.lt.50) THEN
        DO ii=2,last
            temp=array(ii)
            DO jj=ii-1,1,-1
                IF (array(jj).le.temp) EXIT
                array(jj+1)=array(jj)
            END DO
            array(jj+1)=temp
        END DO
        RETURN
    END IF

    temp=array(last/2)
    array(last/2)=array(2)
    IF (temp.gt.array(last)) THEN
        array(2)=array(last)
        array(last)=temp
    ELSE
        array(2)=temp
    END IF
    IF (array(1).gt.array(last)) THEN
        temp=array(1)
        array(1)=array(last)
        array(last)=temp
    END IF
    IF (array(1).gt.array(2)) THEN
        temp=array(1)
        array(1)=array(2)
        array(2)=temp
    END IF
    pivot=array(2)

    left=3
    right=last-1
    DO
        DO WHILE(array(left).lt.pivot)
            left=left+1
        END DO
        DO WHILE(array(right).gt.pivot)
          right=right-1
        END DO
        IF (left.ge.right) EXIT
        temp=array(left)
        array(left)=array(right)
        array(right)=temp
        left=left+1
        right=right-1
    END DO
    IF (left.eq.right) left=left+1
    CALL QUICKSORT(array(1:left-1))
    CALL QUICKSORT(array(left:))

END SUBROUTINE QUICKSORT

RECURSIVE SUBROUTINE QUICKSORT_ARG(array, arg)

    IMPLICIT NONE

    DOUBLE PRECISION, DIMENSION(:), INTENT(INOUT)::array
    INTEGER, DIMENSION(:), INTENT(INOUT)::arg
    DOUBLE PRECISION :: temp,pivot
    INTEGER :: temparg
    INTEGER :: ii, jj, last, left, right

    last=SIZE(array)

    IF (last.lt.50) THEN
        DO ii=2,last
            temp=array(ii)
            temparg=arg(ii)
            DO jj=ii-1,1,-1
                IF (array(jj).le.temp) EXIT
                array(jj+1)=array(jj)
                arg(jj+1)=arg(jj)
            END DO
            array(jj+1)=temp
            arg(jj+1)=temparg
        END DO
        RETURN
    END IF

    temp=array(last/2)
    temparg=arg(last/2)
    array(last/2)=array(2)
    arg(last/2)=arg(2)
    IF (temp.gt.array(last)) THEN
        array(2)=array(last)
        arg(2)=arg(last)
        array(last)=temp
        arg(last)=temparg
    ELSE
        array(2)=temp
        arg(2)=temparg
    END IF
    IF (array(1).gt.array(last)) THEN
        temp=array(1)
        temparg=arg(1)
        array(1)=array(last)
        arg(1)=arg(last)
        array(last)=temp
        arg(last)=temparg
    END IF
    IF (array(1).gt.array(2)) THEN
        temp=array(1)
        temparg=arg(1)
        array(1)=array(2)
        arg(1)=arg(2)
        array(2)=temp
        arg(2)=temparg
    END IF
    pivot=array(2)

    left=3
    right=last-1
    DO
        DO WHILE(array(left).lt.pivot)
            left=left+1
        END DO
        DO WHILE(array(right).gt.pivot)
          right=right-1
        END DO
        IF (left.ge.right) EXIT
        temp=array(left)
        temparg=arg(left)
        array(left)=array(right)
        arg(left)=arg(right)
        array(right)=temp
        arg(right)=temparg
        left=left+1
        right=right-1
    END DO
    IF (left.eq.right) left=left+1
    CALL QUICKSORT_ARG(array(1:left-1), arg(1:left-1))
    CALL QUICKSORT_ARG(array(left:), arg(left:))

END SUBROUTINE QUICKSORT_ARG

SUBROUTINE ARGSORT(array, arg)

    DOUBLE PRECISION, DIMENSION(:), INTENT(IN)::array
    DOUBLE PRECISION, DIMENSION(SIZE(array))::aux_array
    INTEGER, DIMENSION(:), INTENT(OUT)::arg
    INTEGER:: n, ii

    aux_array(:)=array(:)

    n=SIZE(array)

    arg = (/ (ii, ii=1, n) /)

    CALL QUICKSORT_ARG(aux_array, arg)

END SUBROUTINE ARGSORT

SUBROUTINE ARGSORT_REVERSE(array,arg)

    IMPLICIT NONE

    DOUBLE PRECISION, DIMENSION(:), INTENT(IN)::array
    DOUBLE PRECISION, DIMENSION(SIZE(array))::aux_array
    INTEGER, DIMENSION(:), INTENT(OUT):: arg
    INTEGER:: n, ii
    INTEGER::head, tail, temp

    aux_array(:)=array(:)

    n=SIZE(array)

    arg = (/ (ii, ii=1,n) /)

    CALL QUICKSORT_ARG(aux_array, arg)

    head = 1
    tail = n
    DO
       IF (head >= tail)  EXIT
       temp    = arg(head)
       arg(head) = arg(tail)
       arg(tail) = temp
       head    = head + 1
       tail    = tail - 1
    END DO

END SUBROUTINE ARGSORT_REVERSE

END MODULE MODULE_SORT

max_stars_repo_name: uibcdf/OpenKinNet

max_issues_repo_name: uibcdf/OpenKinNet

max_forks_repo_name: uibcdf/OpenKinNet


!-----------------------------------------------------------------------
!
! RAdS May2014.
! Copyright (C) 2011-2014 by Rodrigo Alvares de Souza.
! Mail: <rsouza01@gmail.com>. Web: "http://www.astro.iag.usp.br/~rsouza/".
! This program may be copied and/or distributed freely. See the
! _ terms and conditions in the files in the doc/ subdirectory.
!
! MODULE: linked_list_class
!
!> @author
!> Rodrigo Souza
!
! DESCRIPTION:
!> Object oriented linked list.
!
!-----------------------------------------------------------------------

module linked_list_class

	implicit none

        integer, parameter :: INSERT_ASC = 1
        integer, parameter :: INSERT_DESC = 2

        !---------------------------------------------------------------
        ! The linked list
        !---------------------------------------------------------------
        type, public :: linked_list

        private

                integer :: insert_type = 0

                !First element
                class(linked_list_node), pointer :: head => null()

                !Last element
                class(linked_list_node), pointer :: tail => null()

        contains

                procedure, public :: set_insert_type =>    &
                        set_insert_type_sub

                procedure, public :: get_insert_type => &
                        get_insert_type_sub

                procedure, public :: add_linked_list_node =>    &
                        add_linked_list_node_sub

                procedure, public :: deallocate_linked_list =>  &
                        deallocate_linked_list_sub

                procedure, public :: to_string =>       &
                        linked_list_to_string_sub

        end type linked_list

        !---------------------------------------------------------------
        !Linked list node
        !---------------------------------------------------------------
        type, public :: linked_list_node
                private
                !Previous element
                class(linked_list_node), pointer :: previous_element => null()
                !Next element
                class(linked_list_node), pointer :: next_element => null()

        contains
                procedure, public :: to_string => &
                        linked_list_node_to_string_sub

        end type linked_list_node

        !---------------------------------------------------------------
        !Restrict access from the actual procedure names
        !---------------------------------------------------------------
        private :: set_insert_type_sub, add_linked_list_node_sub, &
                linked_list_to_string_sub, get_insert_type_sub
        private :: deallocate_linked_list_sub

	contains

        !---------------------------------------------------------------
        !Sets the insertion type
        !---------------------------------------------------------------
        subroutine set_insert_type_sub(this, insert_type)
                !Parameters
                class(linked_list) :: this
                integer :: insert_type

                if(this%insert_type == 0) then
                        this%insert_type = insert_type
                end if

        end subroutine set_insert_type_sub

        !---------------------------------------------------------------
        !Gets the insertion type
        !---------------------------------------------------------------
        integer function get_insert_type_sub(this)
                !Parameters
                class(linked_list) :: this

                get_insert_type_sub = this%insert_type

        end function get_insert_type_sub

        !--------------------------------
        !---------------------------------------------------------------
        !Allocates the linked list node
        !---------------------------------------------------------------
        subroutine add_linked_list_node_sub(this, node)
                !Parameters
                class(linked_list) :: this
                class(linked_list_node), pointer :: node
                !Variables
                class(linked_list_node), pointer :: temp_node => null()

                if(.not. associated(this%head)) then

                        this%head => node
                        this%tail => node

                        nullify(this%head%previous_element)
                else
                        this%tail%next_element => node
                        this%tail => node

                        nullify(this%tail%next_element)
                end if

        end subroutine add_linked_list_node_sub

        !---------------------------------------------------------------
        !Dellocates the linked list
        !---------------------------------------------------------------
        subroutine deallocate_linked_list_sub(this)
                class(linked_list) :: this

                !Variables
                class(linked_list_node), pointer :: current_node
                class(linked_list_node), pointer :: temp_node

                current_node => this%head

                do while (associated(current_node))

                        temp_node => current_node
                        current_node => current_node%next_element

                        deallocate(temp_node)
                end do

                nullify(this%head)
                nullify(this%tail)

        end subroutine deallocate_linked_list_sub

        !---------------------------------------------------------------
        !To string
        !---------------------------------------------------------------
        subroutine linked_list_to_string_sub(this)
                !Parameters
                class(linked_list), intent(in) :: this
                !Variables
                class(linked_list_node), pointer :: current_node

                current_node => this%head

                do while (associated(current_node))
                        call current_node%to_string
                        current_node => current_node%next_element
                end do

        end subroutine linked_list_to_string_sub

        !---------------------------------------------------------------
        !To string
        !---------------------------------------------------------------
        subroutine linked_list_node_to_string_sub(this)
                class(linked_list_node), intent(in) :: this

                write (*,*) "linked_list_node%to_string"

        end subroutine linked_list_node_to_string_sub

end module linked_list_class


text:
There are two active landfills in the Davis outskirts. One is the Yolo County Central Landfill and the other is the UC Davis Landfill.

The hill between the Blue Max Kart Club and PVUSA was the active landfill for Davis through at least the early 1970s.

There are also numerous retired landfills that should be documented here.

 Message from Paul Palmer  March 26, 2007

The word landfill was invented some time in the 1970s by the garbage industry to replace the less savory but very descriptive word dump. It makes the concept of dump seem less repulsive and is therefore part of a promotional effort by the garbage industry to sell the dump concept. As such, no one in the recycling or zero waste or environmental movement should ever use the term landfill but should insist on the more raw term dump. Those who oppose the idea of throwing garbage into holes in the ground should not be swayed by the governments use of landfill, such as in Yolo County Landfill but should insist on replacing the name with Yolo County Dump in all communications and public testimony.

The pictures in this Wiki show UC Davis Landfill seagulls alighting on the UCD dump. In California, every dump is required to exclude sea gulls (known along with rodents in garbage terminology as vermin) from all dumps. To accomplish this, the dump operators traditionally covered the working face i.e. the days garbage, with dirt every day. In 1995, the dumpers realized that they could use this habit to create bogus recycling numbers that would meet the requirements of AB 939 (1989) that mandate constantly rising recycling rates. They prevailed upon Bustamonte to carry a bill, known as AB 1647, which gave them permission to take any element of garbage (except organic which might attract vermin but including compost), spread it on the dump (at the working face) and call it Alternative Daily Cover (ADC). This in lieu of covering with dirt. The payoff is that when they do this, the ADC is allowed by law to be called RECYCLING. In other words, they can take a large tonnage out of the garbage column of their reports and arbitrarily move it into the Recycling column. As a result of this clear scam, the recycling rates have skyrocketed, even though nothing changed in the dump. The recyclers need to familiarize themselves with this situation as governments continually use these bogus recycling numbers to sit on their hands and preen themselves for their high recycling rates. In Marin County for example, the recycling rate in 2006 is reported to be 77% which is outrageously inflated ( a more realistic rate might be 10 to 20 percent) even though every obstacle to recycling is in place there, such as single, mixed garbage pickup and domination by a garbage company. There is not even any definition of how long a day is so the dumper can transfer garbage to ADC every hour (for example) and call it daily cover (I dont know if this is actually done but the law allows it).

The worst aspect of this garbage and dump friendly legislation is that it is completely unnecessary. Composting is a well known preferred alternative to dumping. If every single piece of food or organic matter were to be rigorously excluded from every dump everywhere by law, the socalled vermin would have nothing to attract them and there would be no bogus rationalization for ADC. In addition, the nutrient and resource values inherent in organic matter would be productively returned to the soil where they belong, not be languishing in dumps uselessly and harmfully.

Paul Palmer, author of http://www.gettingtozerowaste.com/ Getting To Zero Waste.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c**********************************************************************
c			IO_GETNS.F
c**********************************************************************
c Determines that number of istat variables there are in a tp.in file
c
c             Input:
c                 iu            ==> unit number (integer)
c
c             Output:
c                 ns            ==>  number of istat variable (int scalar)
c
c
c Remarks:
c Authors:  Hal Levison
c Date:    10/1/96
c Last revision:  3/18/97

      subroutine io_getns(iu,ns)

      include '../swift.inc'
      include 'io.inc'

c...  Input
      integer iu

c...  Output
      integer ns

c...  Internal
      character*1024 line
      integer i,i1,ib
      real*8 xht,yht,zht
      real*8 vxht,vyht,vzht

c-----
c...  Executable code

c...  get the irrelavant stuff
      read(7,*) xht,yht,zht
      read(7,*) vxht,vyht,vzht

      ns = 0
      do while(.true.)
         read(7,fmt='(a)') line

c...     if there are `.' then it is not a istat line. Therefore leave
         do i = 1,1024
            if(line(i:i).eq.'.') goto 99
         enddo

c...     Find the first non-blank character
         i1 = 0
         do i=1,1024
            if( (i1.eq.0) .and. (line(i:i).ne.' ') ) then
               i1 = i
            endif
         enddo

         ib = 1
         do i=i1+1,1024
            if( (ib.eq.1) .and. (line(i:i).eq.' ') ) then
               ns = ns + 1
            endif
            if (line(i:i).eq.' ') then
               ib = 0
            else
               ib = 1
            endif
         enddo

      enddo
 99   continue

      rewind(7)

      return
      end     ! io_getns
c----------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
! { dg-do compile }
!
! NULL(...) pointer is not allowed as operand
! PR fortran/20888
!
! Contributed by Joost VandeVondele
!
PROGRAM main
  IMPLICIT NONE
  REAL, POINTER :: TEST
  NULLIFY(TEST)
  TEST => -NULL(TEST) ! { dg-error "Invalid context for NULL" }
  IF (TEST .EQ. NULL(TEST)) TEST=>NULL() ! { dg-error "Invalid context for NULL" }
  IF (NULL(TEST) .EQ. TEST) TEST=>NULL() ! { dg-error "Invalid context for NULL" }
END PROGRAM main

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Relaxation is the goal of the soul.

I believe your session should be all about you! I offer Massage, Reiki, Guided Visualization, and Breathing Exercises to help you reduce stress and increase relaxation. Together we choose the most effective approach for nurturing your Body, Mind & Spirit.

Massage: Ahhh... the benefits of massage include increased circulation, relaxed muscles, reduced anxiety and depression, improved posture, less stress, relaxation & much more....

Reiki: Reiki is a type of energy work that is applied through light touch. It balances the energy centers in the body and triggers your natural healing abilities.

Guided Visualization: Let your mind and body relax as you are guided in a visualization exercise that will be a calming experience and one you will learn to access whenever you desire.

Breathing Exercises: When we become aware of our breathing and practice breathing exercises we can either calm or energize ourselves.

RATES

60 Minutes (23 Modalities)
$70.00

90 Minutes (34 modalities)
$95.00

CONTACT INFO

(530) 3048934
2657 Portage Bay, East, Davis, CA 95618

Men by referral only.


text:
Helicopters can frequently be seen hovering over Davis, especially at Night Time.

A search in the DPD website for the word helicopter yields no results.

Whos helicopters are these?  Are they simply consuming fuel to justify the existence of helicopters in their (whoever they are) budgets?  Or is this a plot even more sinister? Users/RishiTrivedi
Davis can request assistants from Sac (city and county) who have a fleet, its part of interagency cooperation Daubert

Ive seen them practicing approaches and landings (at night) at the UC Davis airport. I suspect that someone may be doing some flight training recently, and has been using the Davis area as part of their training. Im not sure, but I think there is a helicopter training school at one of the Sacramento airports. I know for sure there is one in Concord. Users/EricKlein

About a year ago there was a helicopter that checked the canal near Fst for a missing senior Daubert

What was with the circling helicopter for over an hour last night (72110) from 1011pm? It seemed to be over the old cannery area.  Users/NoelBruening
 The closest fit in the police log...
  22:38 SUSPICIOUS CIRCUMSTANCE 100721137
  Occurred on Cranbrook Ct. UNKNOWN SUBJECT WAS TRYING TO OPEN LOCKED DOOR  Disposition: Service
 Users/WilliamLewis

I saw him as well, it looked as if they were checking H street and the train tracks, after circling for a while it peeled of and headed back to sac Daubert
Helicopter over head tonight broadcasting with loudspeaker ~SD
Also last week a squad of attack helicopters & some hueys were seen heading towards Travis


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

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

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

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

! End MIT license text.

      SUBROUTINE FBS_SUPRLU ( CALLING_SUBR, MATIN_NAME, NROWS, NTERMS, I_MATIN, J_MATIN, MATIN, ICOL, RHS_COL, INFO )

! FBS_SUPRLU performs the forward-backward substitution to get displacements after the stiffness matrix has been decomposed using
! subr SYM_MAT_DECOMP_SUPRLU

! (1) Scales the INOUT (e.g. load) vector if requested
! (2) Does the forward-backward solution to solve for a left-hand side vector given 1 right-hand side input vector (INOUT_COL)
!     The reult is returned in INOUT_COL
! (3) Scales the solution vector if requested

      USE PENTIUM_II_KIND, ONLY       :  LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_LOG, ERR, F04, F06, SC1
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE CONSTANTS_1, ONLY           :  ZERO
      USE PARAMS, ONLY                :  CRS_CCS
      USE SCRATCH_MATRICES, ONLY      :  I_CCS1, J_CCS1, CCS1
      USE SUBR_BEGEND_LEVELS, ONLY    :  FBS_SUPRLU_BEGEND
      USE SuperLU_STUF, ONLY          :  SLU_FACTORS
      USE FBS_SUPRLU_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'FBS_SUPRLU'
      CHARACTER(LEN=*), INTENT(IN)    :: CALLING_SUBR      ! The subr that called this subr (used for output error purposes)
      CHARACTER(LEN=*), INTENT(IN)    :: MATIN_NAME        ! Name of matrix to be decomposed

      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG), INTENT(IN)       :: ICOL              ! Internal subcase or row number for which the FBS is being performed
      INTEGER(LONG), INTENT(IN)       :: NROWS             ! Number of rows in sparse matrix MATIN
      INTEGER(LONG), INTENT(IN)       :: NTERMS            ! Number of nonzeros in sparse matrix MATIN
      INTEGER(LONG), INTENT(IN)       :: I_MATIN(NROWS+1)  ! Indicators of number of nonzero terms in rows of matrix MATIN
      INTEGER(LONG), INTENT(IN)       :: J_MATIN(NTERMS)   ! Col numberts of nonzero terms in matrix MATIN

      INTEGER(LONG), INTENT(INOUT)    :: INFO              ! Output from SuperLU routine

      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = FBS_SUPRLU_BEGEND

      REAL(DOUBLE) , INTENT(IN)       :: MATIN(NTERMS)     ! A small number to compare real zero
      REAL(DOUBLE) , INTENT(IN)       :: RHS_COL(NROWS)    ! RHS column for which the FBS is solving

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

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

      IF      (CRS_CCS == 'CRS') THEN                      ! Use KLL stored in Compressed Row Storage (CRS) format

          CALL C_FORTRAN_DGSSV( 2, NROWS, NTERMS, 1, MATIN, J_MATIN, I_MATIN, RHS_COL, NROWS, SLU_FACTORS, INFO )

      ELSE IF (CRS_CCS == 'CCS') THEN                      ! Use KLL stored in Compressed Col Storage (CCS) format

         CALL ALLOCATE_SCR_CCS_MAT ( 'CCS1', NROWS, NTERMS, SUBR_NAME )
         CALL SPARSE_CRS_SPARSE_CCS ( NROWS, NROWS, NTERMS, MATIN_NAME, I_MATIN, J_MATIN, MATIN, 'CCS1', J_CCS1, I_CCS1, CCS1, 'Y' )
         CALL C_FORTRAN_DGSSV( 2, NROWS, NTERMS, 1, CCS1, I_CCS1, J_CCS1, RHS_COL, NROWS, SLU_FACTORS, INFO )

      ELSE

         WRITE(ERR,933) SUBR_NAME, 'CRS_CCS'
         WRITE(F06,933) SUBR_NAME, 'CRS_CCS'
         FATAL_ERR = FATAL_ERR + 1
         CALL OUTA_HERE ( 'Y' )

      ENDIF

      IF (INFO .EQ. 0) THEN
         WRITE (SC1,9904) ICOL, TRIM(SUBR_NAME)
         WRITE (F06,9904) ICOL, TRIM(SUBR_NAME)
      ELSE
         WRITE( * ,*) 'INFO FROM SUPERLU TRIANGULAR SOLVE = ', INFO
         WRITE(ERR,9903) INFO, TRIM(SUBR_NAME), TRIM(CALLING_SUBR)
         WRITE(F06,9903) INFO, TRIM(SUBR_NAME), TRIM(CALLING_SUBR)
         FATAL_ERR = FATAL_ERR + 1
         CALL OUTA_HERE ( 'Y' )
      ENDIF

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

      RETURN

!***********************************************************************************************************************************
  932 FORMAT(' *ERROR   932: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' PARAMETER ', A, ' MUST BE EITHER "SYM" OR "NONSYM" BUT VALUE IS ',A)

  933 FORMAT(' *ERROR   933: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' PARAMETER ', A, ' MUST BE EITHER "CRS" OR "CCS" BUT VALUE IS ',A)

 9903 FORMAT(' *ERROR  9902: SUPERLU SPARSE SOLVER HAS FAILED WITH INFO = ', I8,' IN SUBR ', A, ' CALLED BY SUBR ', A)

 9904 FORMAT(' SUPERLU SPARSE SOLVER SUCCESSFUL FOR CASE ', I8,' IN SUBR ', A)

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

      END SUBROUTINE FBS_SUPRLU

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


text:
!!!    This program sloves Lid Driven Cavity Flow problem using Lattice Boltzmann Method
!!!    Lattice Boltzmann Equation with MRT-LBE model

module commondata
    implicit none

    integer, parameter :: nx=201,ny=201
    integer, parameter :: nxHalf=(nx-1)/2+1, nyHalf=(ny-1)/2+1
    real(8), parameter :: Reynolds=1000.0d0
    real(8), parameter :: rho0=1.0d0
    real(8), parameter :: U0=0.1d0
    real(8), parameter :: tauf=U0*dble(nx)/Reynolds*3.0d0+0.5d0

    integer :: itc
    integer, parameter :: itc_max=INT(50000000)

    real(8) :: errorU
    real(8), parameter :: eps=1e-6

    real(8) :: xp(0:nx+1), yp(0:ny+1)
    real(8), allocatable :: u(:,:), v(:,:)
    real(8), allocatable :: rho(:,:)
    real(8), allocatable :: up(:,:), vp(:,:)

    real(8), allocatable :: f(:,:,:), f_post(:,:,:)

    real(8) :: omega(0:8)
    integer :: ex(0:8), ey(0:8)
    data ex/0, 1, 0, -1,  0, 1, -1, -1,  1/
    data ey/0, 0, 1,  0, -1, 1,  1, -1, -1/
    real(8) :: Uwall(1:2)
    data Uwall/U0,  0/

    real(8), parameter :: Snu=1.0d0/tauf, Sq=8.0d0*(2.0d0*tauf-1.0d0)/(8.0d0*tauf-1.0d0)
end module commondata

program main
    use commondata
    implicit none
    real(8) :: start, finish
    real(8) :: start2, finish2
    integer :: myMaxThreads

    write(*,*) "--------------------------"
    write(*,*) "I am MRT collision operator"
    write(*,*) "--------------------------"

    call initial()

    call CPU_TIME(start)

    do while((errorU.GT.eps).AND.(itc.LE.itc_max))

        itc = itc+1

        call collision()

        call streaming()

        call bounceback()

        call macro()

        if(MOD(itc,2000).EQ.0) then
            call check()
            ! call output_ASCII()
        endif
        !if(MOD(itc,100000).EQ.0) then
        !    call output_Tecplot()
        !endif

    enddo

    call CPU_TIME(finish)

    write(*,*) "Time (CPU) = ", real(finish-start), "s"

    itc = itc+1
    call output_Tecplot()
    call output_ASCII()
    ! call output_binary()
    call getVelocity()

    write(*,*) "Deallocate Array..."
    deallocate(u)
    deallocate(v)
    deallocate(rho)
    deallocate(up)
    deallocate(vp)
    deallocate(f)
    deallocate(f_post)
    write(*,*) "    "

    write(*,*) "Successfully: DNS completed!"

    stop
end program main

subroutine initial()
    use commondata
    implicit none
    integer :: i, j
    integer :: alpha
    real(8) :: un(0:8)
    real(8) :: us2

    itc = 0
    errorU = 100.0d0

    write(*,*) "nx=",nx,", ny=",ny
    write(*,*) "Reynolds=",real(Reynolds)
    write(*,*) "U0=",real(U0),",    tauf=",real(tauf)

    xp(0) = 0.0d0
    xp(nx+1) = dble(nx)
    do i=1,nx
        xp(i) = dble(i)-0.5d0
    enddo
    yp(0) = 0.0d0
    yp(ny+1) = dble(ny)
    do j=1,ny
        yp(j) = dble(j)-0.5d0
    enddo

    allocate (u(nx,ny))
    allocate (v(nx,ny))
    allocate (rho(nx,ny))
    allocate (up(nx,ny))
    allocate (vp(nx,ny))

    allocate (f(0:8,nx,ny))
    allocate (f_post(0:8,0:nx+1,0:ny+1))

    rho = rho0
    u = 0.0d0
    v = 0.0d0
    up = 0.0d0
    vp = 0.0d0

    do i=1,nx
        u(i,ny) = U0
    enddo

    omega(0) = 4.0d0/9.0d0
    do alpha=1,4
        omega(alpha) = 1.0d0/9.0d0
    enddo
    do alpha=5,8
        omega(alpha) = 1.0d0/36.0d0
    enddo

    do j=1,ny
        do i=1,nx
            us2 = u(i,j)*u(i,j)+v(i,j)*v(i,j)
            do alpha=0,8
                un(alpha) = u(i,j)*ex(alpha)+v(i,j)*ey(alpha)
                f(alpha,i,j) = rho(i,j)*omega(alpha)*(1.0d0+3.0d0*un(alpha)+4.5d0*un(alpha)*un(alpha)-1.5d0*us2)
            enddo
        enddo
    enddo

    return
end subroutine initial

subroutine collision()
    use commondata
    implicit none
    integer :: i, j
    integer :: alpha
    !---------------------------
    real(8) :: s(0:8)
    real(8) :: m(0:8)
    real(8) :: m_post(0:8)
    real(8) :: meq(0:8)
    !---------------------------

    do j=1,ny
        do i=1,nx

    m(0) = f(0,i,j)+f(1,i,j)+f(2,i,j)+f(3,i,j)+f(4,i,j)+f(5,i,j)+f(6,i,j)+f(7,i,j)+f(8,i,j)
    m(1) = -4.0d0*f(0,i,j)-f(1,i,j)-f(2,i,j)-f(3,i,j)-f(4,i,j)+2.0d0*(f(5,i,j)+f(6,i,j)+f(7,i,j)+f(8,i,j))
    m(2) = 4.0d0*f(0,i,j)-2.0d0*(f(1,i,j)+f(2,i,j)+f(3,i,j)+f(4,i,j))+f(5,i,j)+f(6,i,j)+f(7,i,j)+f(8,i,j)
    m(3) = f(1,i,j)-f(3,i,j)+f(5,i,j)-f(6,i,j)-f(7,i,j)+f(8,i,j)
    m(4) = -2.0d0*f(1,i,j)+2.0d0*f(3,i,j)+f(5,i,j)-f(6,i,j)-f(7,i,j)+f(8,i,j)
    m(5) = f(2,i,j)-f(4,i,j)+f(5,i,j)+f(6,i,j)-f(7,i,j)-f(8,i,j)
    m(6) = -2.0d0*f(2,i,j)+2.0d0*f(4,i,j)+f(5,i,j)+f(6,i,j)-f(7,i,j)-f(8,i,j)
    m(7) = f(1,i,j)-f(2,i,j)+f(3,i,j)-f(4,i,j)
    m(8) = f(5,i,j)-f(6,i,j)+f(7,i,j)-f(8,i,j)

            meq(0) = rho(i,j)
            meq(1) = rho(i,j)*( -2.0d0+3.0d0*(u(i,j)*u(i,j)+v(i,j)*v(i,j)) )
            meq(2) = rho(i,j)*( 1.0d0-3.0d0*(u(i,j)*u(i,j)+v(i,j)*v(i,j)) )
            meq(3) = rho(i,j)*u(i,j)
            meq(4) = -rho(i,j)*u(i,j)
            meq(5) = rho(i,j)*v(i,j)
            meq(6) = -rho(i,j)*v(i,j)
            meq(7) = rho(i,j)*( u(i,j)*u(i,j)-v(i,j)*v(i,j) )
            meq(8) = rho(i,j)*( u(i,j)*v(i,j) )

            s(0) = 0.0d0      !!s_{\rho}
            s(1) = Snu !!s_{e}
            s(2) = Snu !!s_{\epsilon}
            s(3) = 0.0d0      !!s_{j}
            s(4) = Sq !!s_{q}
            s(5) = 0.0d0      !!s_{j}
            s(6) = Sq       !!s_{q}
            s(7) = Snu !!s_{\nu}
            s(8) = Snu       !!s_{\nu}

            do alpha=0,8
                m_post(alpha) = m(alpha)-s(alpha)*(m(alpha)-meq(alpha))
            enddo

    f_post(0,i,j) = ( m_post(0)-m_post(1)+m_post(2) )/9.0d0
    f_post(1,i,j) = m_post(0)/9.0d0-m_post(1)/36.0d0-m_post(2)/18.0d0+m_post(3)/6.0d0-m_post(4)/6.0d0 &
                    +m_post(7)*0.25d0
    f_post(2,i,j) = m_post(0)/9.0d0-m_post(1)/36.0d0-m_post(2)/18.0d0 &
                    +m_post(5)/6.0d0-m_post(6)/6.0d0-m_post(7)*0.25d0
    f_post(3,i,j) = m_post(0)/9.0d0-m_post(1)/36.0d0-m_post(2)/18.0d0-m_post(3)/6.0d0+m_post(4)/6.0d0 &
                    +m_post(7)*0.25d0
    f_post(4,i,j) = m_post(0)/9.0d0-m_post(1)/36.0d0-m_post(2)/18.0d0 &
                    -m_post(5)/6.0d0+m_post(6)/6.0d0-m_post(7)*0.25d0
    f_post(5,i,j) = m_post(0)/9.0d0+m_post(1)/18.0d0+m_post(2)/36.0d0+m_post(3)/6.0d0+m_post(4)/12.0d0 &
                    +m_post(5)/6.0d0+m_post(6)/12.0d0+m_post(8)*0.25d0
    f_post(6,i,j) = m_post(0)/9.0d0+m_post(1)/18.0d0+m_post(2)/36.0d0-m_post(3)/6.0d0-m_post(4)/12.0d0 &
                    +m_post(5)/6.0d0+m_post(6)/12.0d0-m_post(8)*0.25d0
    f_post(7,i,j) = m_post(0)/9.0d0+m_post(1)/18.0d0+m_post(2)/36.0d0-m_post(3)/6.0d0-m_post(4)/12.0d0 &
                    -m_post(5)/6.0d0-m_post(6)/12.0d0+m_post(8)*0.25d0
    f_post(8,i,j) = m_post(0)/9.0d0+m_post(1)/18.0d0+m_post(2)/36.0d0+m_post(3)/6.0d0+m_post(4)/12.0d0 &
                    -m_post(5)/6.0d0-m_post(6)/12.0d0-m_post(8)*0.25d0
        enddo
    enddo

    return
end subroutine collision

subroutine streaming()
    use commondata
    implicit none
    integer :: i, j
    integer :: ip, jp
    integer :: alpha

    do j=1,ny
        do i=1,nx
            do alpha=0,8
                ip = i-ex(alpha)
                jp = j-ey(alpha)

                f(alpha,i,j) = f_post(alpha,ip,jp)

            enddo
        enddo
    enddo

    return
end subroutine streaming

subroutine bounceback()
    use commondata
    implicit none
    integer :: i, j

    do j=1,ny
        !Left side (i=1)
        f(1,1,j) = f_post(3,1,j)
        f(5,1,j) = f_post(7,1,j)
        f(8,1,j) = f_post(6,1,j)

        !Right side (i=nx)
        f(3,nx,j) = f_post(1,nx,j)
        f(6,nx,j) = f_post(8,nx,j)
        f(7,nx,j) = f_post(5,nx,j)
    enddo

    do i=1,nx
        !Bottom side (j=1)
        f(2,i,1) = f_post(4,i,1)
        f(5,i,1) = f_post(7,i,1)
        f(6,i,1) = f_post(8,i,1)

        !Top side (j=ny)
        f(4,i,ny) = f_post(2,i,ny)
        f(7,i,ny) = f_post(5,i,ny)-rho(i,ny)*(Uwall(1))/6.0d0
        f(8,i,ny) = f_post(6,i,ny)-rho(i,ny)*(-Uwall(1))/6.0d0
    enddo

    return
end subroutine bounceback

subroutine macro()
    use commondata
    implicit none
    integer :: i, j

    do j=1,ny
        do i=1,nx
            rho(i,j) = f(0,i,j)+f(1,i,j)+f(2,i,j)+f(3,i,j)+f(4,i,j)+f(5,i,j)+f(6,i,j)+f(7,i,j)+f(8,i,j)
            u(i,j) = ( f(1,i,j)-f(3,i,j)+f(5,i,j)-f(6,i,j)-f(7,i,j)+f(8,i,j) )/rho(i,j)
            v(i,j) = ( f(2,i,j)-f(4,i,j)+f(5,i,j)+f(6,i,j)-f(7,i,j)-f(8,i,j) )/rho(i,j)
        enddo
    enddo

    return
end subroutine macro

subroutine check()
    use commondata
    implicit none
    integer :: i, j
    real(8) :: error1, error2

    error1 = 0.0d0
    error2 = 0.0d0

    do j=1,ny
        do i=1,nx
            error1 = error1+(u(i,j)-up(i,j))*(u(i,j)-up(i,j))+(v(i,j)-vp(i,j))*(v(i,j)-vp(i,j))
            error2 = error2+u(i,j)*u(i,j)+v(i,j)*v(i,j)

            up(i,j) = u(i,j)
            vp(i,j) = v(i,j)
        enddo
    enddo

    errorU = sqrt(error1)/sqrt(error2)

    write(*,*) itc,' ',errorU

    return
end subroutine check

subroutine output_ASCII()
    use commondata
    implicit none
    integer :: i, j
    character(len=100) :: filename

    write(filename,*) itc
    filename = adjustl(filename)

    open(unit=02,file='MRTcavity-'//trim(filename)//'.dat',status='unknown')
        write(02,*) 'TITLE="Lid Driven Cavity"'
        write(02,*) 'VARIABLES="X" "Y" "U" "V" "Pressure" '
        write(02,101) nx, ny
        do j=1,ny
            do i=1,nx
                write(02,100) xp(i), yp(j), u(i,j), v(i,j), rho(i,j)/3.0d0
            enddo
        enddo
100     format(1x,2(e11.4,' '),10(e13.6,' '))
101     format('ZONE',1x,'I=',1x,i5,2x,'J=',1x,i5,1x,'F=POINT')
    close(02)

    return
end subroutine output_ASCII

subroutine output_binary()
    use commondata
    implicit none
    integer :: i, j
    character(len=100) :: filename

    write(filename,*) itc
    filename = adjustl(filename)

    open(unit=01,file='MRTcavity-'//trim(filename)//'.bin',form="unformatted",access="sequential")
    write(01) ((u(i,j),i=1,nx),j=1,ny)
    write(01) ((v(i,j),i=1,nx),j=1,ny)
    write(01) ((rho(i,j),i=1,nx),j=1,ny)
    close(01)

    return
    end subroutine output_binary

!!!c--------------------------------
!!!c--------------------------------
    subroutine output_Tecplot()
    use commondata
    implicit none
    integer :: i, j, k
    character(len=9) :: B2
    REAL(4) :: zoneMarker, eohMarker
    character(len=40) :: title
    character(len=40) :: V1,V2,V3,V4,V5
    integer, parameter :: kmax=1
    character(len=40) :: zoneName

    write(B2,'(i9.9)') itc

    open(41,file='MRTcavity-'//B2//'.plt', access='stream', form='unformatted')
    !---------------------------------------------
    zoneMarker= 299.0
    eohMarker = 357.0

    !I. HEAD SECTION--------------------------------------
    !c--Magic number, Version number
    write(41) "#!TDV101"

    !c--Integer value of 1
    write(41) 1

    Title="MyFirst"
    call dumpstring(title)

    !c-- Number of variables in this data file (here 5 variables)
    write(41) 5

    !c-- Variable names.
    V1='X'
    call dumpstring(V1)
    V2='Y'
    call dumpstring(V2)
    V3='U'
    call dumpstring(V3)
    V4='V'
    call dumpstring(V4)
    V5='Pressure'
    call dumpstring(V5)

    !c-----Zones-----------------------------

    !c--------Zone marker. Value = 299.0
    write(41) zoneMarker

    !--------Zone name.
    zoneName="ZONE 001"
    call dumpstring(zoneName)

    !---------Zone Color
    write(41) -1

    !---------ZoneType
    write(41) 0

    !---------DataPacking 0=Block, 1=Point
    write(41) 1

    !---------Specify Var Location. 0 = Do not specify, all data
    !---------is located at the nodes. 1 = Specify
    write(41) 0

    !---------Number of user defined face neighbor connections
    ! (value >= 0)
    write(41) 0

    !---------IMax,JMax,KMax
    write(41) nx
    write(41) ny
    write(41) kmax

    !-----------1=Auxiliary name/value pair to follow
    !-----------0=No more Auxiliar name/value pairs.
    write(41) 0
    write(41) eohMarker

    !----zone ------------------------------------------------------------
    write(41) zoneMarker

    !--------variable data format, 1=Float, 2=Double, 3=LongInt,4=ShortInt, 5=Byte, 6=Bit
    write(41) 1
    write(41) 1
    write(41) 1
    write(41) 1
    write(41) 1

    !--------Has variable sharing 0 = no, 1 = yes.
    write(41) 0

    !----------Zone number to share connectivity list with (-1 = no
    ! sharing).
    write(41) -1

    !---------------------------------------------------------------------
    do k=1,kmax
        do j=1,ny
            do i=1,nx
                write(41) real(xp(i))
                write(41) real(yp(j))
                write(41) real(u(i,j))
                write(41) real(v(i,j))
                write(41) real(rho(i,j)/3.0d0)
            end do
        end do
    enddo
    close(41)
    !---------------------------------------------------------------------

    return
    end subroutine output_Tecplot
!!!c--------------------------------
    subroutine dumpstring(instring)
    implicit none
    character(len=40) instring
    integer :: stringLength
    integer :: ii
    integer :: I

    stringLength=LEN_TRIM(instring)
    do ii=1,stringLength
        I=ICHAR(instring(ii:ii))
        write(41) I
    end do
    write(41) 0

    return
    end subroutine dumpstring
!!!c--------------------------------
!!!c--------------------------------

subroutine getVelocity()
    use commondata
    implicit none
    integer :: i, j

    open(unit=02,file='./u-y.dat',status='unknown')
    do j=1,ny
        write(02,*) u(nxHalf,j)/U0, yp(j)/dble(nx)
    enddo
    close(02)

    open(unit=03,file='./x-v.dat',status='unknown')
    do i=1,nx
        write(03,*) xp(i)/dble(nx), v(i,nyHalf)/U0
    enddo
    close(03)

    return
end subroutine getVelocity

max_stars_repo_name: cheryli/parallel

max_issues_repo_name: cheryli/parallel

max_forks_repo_name: cheryli/MGLC


text:
!==============================================================================!
  real function Turb_Mod_Ebf_Scalar(turb, c, coef)
!------------------------------------------------------------------------------!
!   Calculates elliptic blending function for scalars.                         !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Turb_Type) :: turb
  integer         :: c
  real            :: coef    ! Prandtl number (maybe even Shmidt number)
!==============================================================================!

  Turb_Mod_Ebf_Scalar = 0.01 * (coef * turb % y_plus(c) ** 4  &
                        / ((1.0 + 5.0 * coef**3 * turb % y_plus(c)) + TINY))

  end function

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
program main
use utils
use g_elimination
use num_integration
use bspline
use finite_elements
use ogpf

implicit none

!---------------
! Variables
!---------------

integer :: i,j,romberg_iterations
integer, allocatable :: n(:)
real(wp) :: x,error
real(wp), allocatable :: error_array_linear(:),error_array_cubic(:)
type(gpf):: gplot

!---------------
! Logic
!---------------

write(*,*) "Program started!"

! Array of number of internal nodes for the integration
allocate(n(7))

! Array for plotting errors
allocate(error_array_linear(7))
allocate(error_array_cubic(7))

! Number of internal nodes for the integration
n(1:) = [7, 15, 31, 63, 127, 255, 511]

romberg_iterations = 2
! ---------------------------
! Solving with linear splines
! ---------------------------
do i = 1, 1
  call solve_finite_elements(n(i),linear_spline_basis,&
    linear_spline_basis,.false.,romberg_iterations,error)
  error_array_linear(i) = error
end do

!call gplot%title('Linear Spline Finite Elements Error')
!call gplot%ylabel('||u_n(x) - u(x)||_{\infty}')
!call gplot%xlabel('Number of nodes')
!call gplot%plot(linspace(7.0_wp,63.0_wp,4),error_array_linear(1:4))

! ---------------------------
! Solving with cubic splines
! ---------------------------
!do i = 1, 4
!  call solve_finite_elements(n(i),cubic_spline_basis,&
!    cubic_spline_basis_deriv,.true.,romberg_iterations,error)
!  error_array_cubic(i) = error
!end do

!call gplot%title('Cubic Spline Finite Elements Error')
!call gplot%ylabel('||u_n(x) - u(x)||_{\infty}')
!call gplot%xlabel('Number of nodes')
!call gplot%plot(linspace(7.0_wp,63.0_wp,4),error_array_cubic(1:4))

contains

end program main

max_stars_repo_name: AndreTGMello/numerical-analysis-course

max_issues_repo_name: AndreTGMello/numerical-analysis-course

max_forks_repo_name: AndreTGMello/numerical-analysis-course


text:
      SUBROUTINE mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,nalp,chisq,
     *funcs)
      INTEGER ma,nalp,ndata,ia(ma),MMAX
      REAL*8 chisq,a(ma),alpha(nalp,nalp),beta(ma),sig(ndata),x(ndata),
     *y(ndata)
      EXTERNAL funcs
      PARAMETER (MMAX=6)
      INTEGER mfit,i,j,k,l,m
      REAL*8 dy,sig2i,wt,ymod,dyda(MMAX)
      mfit=0
      do 11 j=1,ma
        if (ia(j).ne.0) mfit=mfit+1
11    continue
      do 13 j=1,mfit
        do 12 k=1,j
          alpha(j,k)=0.
12      continue
        beta(j)=0.
13    continue
      chisq=0.
      do 16 i=1,ndata
        call funcs(x(i),a,ymod,dyda,ma)
        sig2i=1./(sig(i)*sig(i))
        dy=y(i)-ymod
        j=0
        do 15 l=1,ma
          if(ia(l).ne.0) then
            j=j+1
            wt=dyda(l)*sig2i
            k=0
            do 14 m=1,l
              if(ia(m).ne.0) then
                k=k+1
                alpha(j,k)=alpha(j,k)+wt*dyda(m)
              endif
14          continue
            beta(j)=beta(j)+dy*wt
          endif
15      continue
        chisq=chisq+dy*dy*sig2i
16    continue
      do 18 j=2,mfit
        do 17 k=1,j-1
          alpha(k,j)=alpha(j,k)
17      continue
18    continue
      return
      END

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
!
!The module file contains two modules: reservoir_h and lake_h
!

module lake_h
save

! Logical for reading or calculation ofthe number of reservoir per size class per year
LOGICAL :: doacudyear
! number of small reservoirs per volume class (5 volume classes)
! rounded to nearest integer
!Allocatable      real acud(nmun,5)
real, allocatable :: acud(:,:)
! number of acudes per volume class (5 volume classes)
! not rounded, for calculation of growth per year
!Allocatable       real acudfloat(nmun,5)
real, allocatable :: acudfloat(:,:)
! Andreas
! number of acudes per volume class (5 volume classes) and year
! specific values for each year, read from small_rservoirs_year.dat
! calculation of growth per year is not used if data are in the file
!Allocatable       real acudfloatyear(nmun,5,200)
real, allocatable :: acudfloatyear(:,:,:)
! initial maximum storage capacity in one acude of each volume classes (m3)
real :: maxlake0(5)
! initial maximum storage capacity in one acude of each volume classes (m3)
real, allocatable :: maxlakesub0(:,:)
! trend in number of acudes for each volume class (+-nbr/year)
!Allocatable       real laketrend(nmun,5)
real, allocatable :: laketrend(:,:)
! maximum capacity of water storage (m3)
!Allocatable       real maxlakewater(nmun,5)
real, allocatable ::  maxlakewater(:,:)
! fraction of generated river runoff being routed through small acudes
real :: intercepted
! water stored in small reservoirs (m3)
!Allocatable      real lakewater(366*nt,nmun,5)
real, allocatable :: lakewater(:,:,:)
! water stored (m3)
!Allocatable      real laketot(366*nt,nmun)
real, allocatable :: laketot(:,:)
! relative water storage (-)
!Allocatable      real laketotfrac(366*nt,nmun)
real, allocatable :: laketotfrac(:,:)
! evaporation from lakes, all classes (m3)
!Allocatable      real lakeevap(366*nt,nmun,5)
real, allocatable :: lakeevap(:,:,:)
! precipitation on lakes, all classes (m3)
!Allocatable      real lakeprec(366*nt,nmun,5)
real, allocatable :: lakeprec(:,:,:)
! consumptive water use (effective extraction) from lakes, all classes (m3)
!Allocatable      real lakeex(366*nt,nmun,5)
real, allocatable :: lakeex(:,:,:)
! area of water storages (km2)
!Allocatable      real lakearea(nmun,5)
real, allocatable :: lakearea(:,:)
! Percentage  of upper limit of reservoir size class that represents the maximum storage capacity
! of the hypothetical representative reservoir of each class in subbasin(-)
real, allocatable :: maxlake_factor(:,:)
! Percentage  of maximum volume per reservoir size class (-)
real :: lake_vol0_factor(5)
! Change on the number of reservoir per size class (-)
real :: lake_change(5)
! river inflow (m3)
!Allocatable       real lakeinflow(366*nt,subasin)
real, allocatable :: lakeinflow(:,:)
! river outflow (m3)
!Allocatable       real lakeoutflow(366*nt,subasin)
real, allocatable ::  lakeoutflow(:,:)
!Allocatable       real muniret(366*nt,subasin)
real, allocatable :: muniret(:,:)
!alpha factor for area-volume relationship (Molle,1989a)[-]
real :: alpha_Molle(5)
!K factor for area-volume relationship (Molle,1989a)[-]
real :: damk_Molle(5)
!Qout=c*Hv**d rating curve of the spillway, where Hv is the water height above the spillway of hypothetical representative reservoir from each class
real :: damc_hrr(5)
real :: damd_hrr(5)
! maximum storage capacity of the hypothetical representative reservoir from each class (m**3)
real, allocatable :: maxlake(:,:)
! maximum storage capacity of the hypothetical representative reservoir from each class (m**3) time-dependent variable
real, allocatable :: maxstorcap_hrr(:,:,:)
! inflow discharge into the hypothetical representative reservoir from each class (m**3)
real, allocatable :: lakeinflow_hrr(:,:,:)
! outflow discharge from the hypothetical representative reservoir from each class (m**3)
real, allocatable :: lakeoutflow_hrr(:,:,:)
! water retention in the hypothetical representative reservoir from each class (m**3)
real, allocatable :: lakeretention_hrr(:,:,:)
! storage volume in the hypothetical representative reservoir from each class (m**3)
real, allocatable :: lakewater_hrr(:,:,:)
!Spillway overflow on the day before in the hypothetical representative reservoir from each class (m3/timestep)
real, allocatable :: outflow_last_hrr(:,:)
!Maximum water depth of the hypothetical representative reservoir from each class used for the overflow calculation (m)
real, allocatable :: hmax_hrr(:,:)
!Water volume above the spillway elevation on the day before in the hypothetical representative reservoir from each class (m**3/timestep)
real, allocatable :: volume_last_hrr(:,:)
!fraction of contributing area of a reservoir class within the sub-basin  (-)
real, allocatable :: lakefrarea(:,:)
!fraction of sub-basin area not controlled by reservoirs (-)
real, allocatable :: subfrarea(:)
!fraction of sub-basin area not controlled by reservoirs receiving outflow discharge from reservoir classes of smaller storage capacity (-)
real, allocatable :: subfrout(:)
!fraction of contributing area of a reservoir class receiving outflow discharge from reservoir classes of smaller storage capacity (-)
real, allocatable :: lakefrout(:,:)
!fraction of sub-basin area receiving outflow discharge from reservoir classes of smaller storage capacity (-)
real, allocatable :: lakecumfrout(:,:)
!generated runoff within the subbasin (m**3)
real, allocatable :: lakerunoff(:,:)

!sediment balance of small reservoirs
! sediment inflow into the small reservoirs (ton)
real, allocatable :: lakesedin(:,:)
! sediment outflow without retention in the small reservoirs (ton)
real, allocatable :: lakesedout(:,:)
! percentage of grain size g that flows into the small reservoirs (-)
real, allocatable :: lakefracsedin(:,:,:)
! percentage of grain size g that flows out of the sub-basin(-)
real, allocatable :: lakefracsedout(:,:,:)
! sediment inflow into the hypothetical representative reservoir from each class (ton)
real, allocatable :: lakesedin_hrr(:,:,:)
! sediment outflow from the hypothetical representative reservoir from each class (ton)
real, allocatable :: lakesedout_hrr(:,:,:)
! percentage of grain size g that flows into hypothetical representative reservoir from each class (-)
real, allocatable :: lakefracsedin_hrr(:,:,:,:)
! percentage of grain size g that flows out of hypothetical representative reservoir from each class(-)
real, allocatable :: lakefracsedout_hrr(:,:,:,:)
! sediment deposition at the hypothetical representative reservoir from each class(-)
real, allocatable :: lakedep_hrr(:,:,:)
! cumulative sediment deposition at the hypothetical representative reservoir from each class(-)
real, allocatable :: lakecumdep_hrr(:,:)
! storage capacity reduction of the hypothetical representative reservoir from each class(-)
real, allocatable :: lake_vollost(:,:,:)
! cumulative sediment deposition at the hypothetical representative reservoir from each class(-)
real, allocatable :: cumseddep(:)
! cumulative sediment deposition at the reservoir class(-)
real, allocatable :: lake_cumseddep(:,:,:)

end module lake_h

max_stars_repo_name: plvoit/WASA-SED

max_issues_repo_name: TillF/WASA-SED

max_forks_repo_name: ERottler/WASA-SED


text:
subroutine example
     real ar(10)
     pointer(ipt, arpte(10))
     real arpte
     ipt = loc(ar)  ! Makes arpte is an alias for ar
     arpte(1) = 1.0 ! Sets ar(1) to 1.0

end subroutine example

max_stars_repo_name: maurizioabba/rose


text:
C*********************************************************************
C
C Program  to evolve the initial conditions created with
C the cluster-filament SETUP program.
C
C From	Peter Quinn Sept. '87
C*********************************************************************
	PROGRAM EVOLVE_CLUSTER
	PARAMETER (max_bodies=10000)
	PARAMETER (G=1.,unit_vel=976.6,unit_mass=2.217e+11)

	REAL xp(max_bodies),yp(max_bodies),zp(max_bodies)
	REAL vxp(max_bodies),vyp(max_bodies),vzp(max_bodies)
	REAL accx(max_bodies),accy(max_bodies),accz(max_bodies)

	LOGICAL is_member

	REAL lum_mass_filgal,lum_mass_member

	CHARACTER*2 orbit_type
	CHARACTER*40 params_file,model_file,orbit_file
	CHARACTER*40 audit_file

	params_file='params'
	model_file='model'
	orbit_file='orbit'
	audit_file='audit'

	OPEN(unit=11,file=model_file,status='unknown')
	READ(11,*)nbodies
	READ(11,*)it,time
	DO i=1,nbodies
	   READ(11,*)j,xp(i),yp(i),work,vxp(i),vyp(i),work
	   zp(i)=0.
	   vzp(i)=0.
	   ENDDO

	PRINT *,'evolve>> Model file read in <<'

C Read PARAMS.DAT
 	OPEN (unit=10,file=params_file,status='unknown')
	READ(10,*)sigma_cluster,core_radius_cluster,
     &            outer_radius_cluster
	READ(10,100)is_member
 100	FORMAT(l1)
	READ(10,*)x_member,y_member,z_member
	READ(10,*)lum_mass_member,re_member
	READ(10,*)nbodies,lum_mass_filgal
	READ(10,*)re_filgal,gamma_halo,extent
        READ(10,*)r_filgal
	READ(10,101)orbit_type
 101	FORMAT(a2)
	READ(10,*)e_fraction,peri_distance
	CLOSE(unit=10)
	PRINT *,'evolve>> Parameters file read in <<'

C Normalization :
	sigma_cluster=sigma_cluster/unit_vel
	lum_mass_member=lum_mass_member/unit_mass
	lum_mass_filgal=lum_mass_filgal/unit_mass

C Parameters of the simulation
	PRINT *,'evolve>> dt (in E+6 years), nsteps',
     1	' and nsteps between output ?'
	READ(5,*)dt,nsteps,nout

	PRINT *,' '
	PRINT *,'evolve>---> dt = ',dt
	PRINT *,'evolve>---> nsteps = ',nsteps
	PRINT *,'evolve>---> nout = ',nout

C Translation of the model :
	x_filgal=r_filgal/sqrt(2.)
	y_filgal=r_filgal/sqrt(2.)
	z_filgal=0.
	DO i=1,nbodies
	   xp(i)=xp(i)+x_filgal
	   yp(i)=yp(i)+y_filgal
	   ENDDO

	PRINT *,'evolve>> Model translated <<'

C Setting up the orbit:
	IF(orbit_type.eq.' c')THEN
	   PRINT *,'evolve>> Setting up a circular orbit <<'
	   CALL cluster_acc(x_filgal,y_filgal,z_filgal,G,
     &                      core_radius_cluster,
     &                      outer_radius_cluster,sigma_cluster,
     &                      accx_filgal,
     &                      accy_filgal,accz_filgal)
	   radial_acc=(x_filgal/r_filgal)*accx_filgal +
     &                (y_filgal/r_filgal)*accy_filgal +
     &                (z_filgal/r_filgal)*accz_filgal
	   v_circ=sqrt(abs(radial_acc)*r_filgal)
	   vx_filgal=-(y_filgal/r_filgal)*v_circ
	   vy_filgal=(x_filgal/r_filgal)*v_circ
	   vz_filgal=0.
	   v_circ_scaled=v_circ*unit_vel
	   PRINT *,'evolve>> Circular velocity = ',v_circ_scaled
	   ELSEIF(orbit_type.eq.' p')THEN
	 	  PRINT *,'evolve>> Setting up a parabolic orbit <<'
	          CALL elements2(r_filgal,peri_distance,core_radius_cluster,
     &                           sigma_cluster,x_filgal,y_filgal,
     &                           vx_filgal,vy_filgal)
		  ELSEIF(orbit_type.eq.' e')THEN
	                 PRINT *,'evolve>> Setting up an elliptical orbit <<'
	                 PRINT *,'evolve>> Energy fraction = ',e_fraction
                         CALL cluster_acc(x_filgal,y_filgal,z_filgal,G,
     &                      core_radius_cluster,
     &                      outer_radius_cluster,sigma_cluster,
     &                      accx_filgal,
     &                      accy_filgal,accz_filgal)
	                  radial_acc=(x_filgal/r_filgal)*accx_filgal +
     &                               (y_filgal/r_filgal)*accy_filgal +
     &                               (z_filgal/r_filgal)*accz_filgal
	                  v_circ=sqrt(e_fraction*abs(radial_acc)*r_filgal)
                          vx_filgal=-(y_filgal/r_filgal)*v_circ
	                  vy_filgal=(x_filgal/r_filgal)*v_circ
	                  vz_filgal=0.
	                  v_circ_scaled=v_circ*unit_vel
			  PRINT *,'evolve>> Initial tangential velocity = ',
     &                                     v_circ_scaled
                          ENDIF
	PRINT *,'evolve>> Orbit established <<'

C Initial velocities:
	DO i=1,nbodies
	   vxp(i)=vxp(i)+vx_filgal
	   vyp(i)=vyp(i)+vy_filgal
	   ENDDO

C Initial accelerations:
	CALL cluster_acc(x_filgal,y_filgal,z_filgal,G,
     &                   core_radius_cluster,
     &                   outer_radius_cluster,sigma_cluster,accx_filgal,
     &                   accy_filgal,accz_filgal)

	DO i=1,nbodies
           CALL filgal_acc(xp(i),yp(i),zp(i),x_filgal,y_filgal,z_filgal,
     &                     G,lum_mass_filgal,
     &                     re_filgal,extent,
     &                     gamma_halo,accx(i),accy(i),accz(i))
           CALL cluster_acc(xp(i),yp(i),zp(i),G,
     &                      core_radius_cluster,
     &                      outer_radius_cluster,sigma_cluster,
     &                      accx_cluster,accy_cluster,accz_cluster)

	   accx(i)=accx(i)+accx_cluster
	   accy(i)=accy(i)+accy_cluster
C	   accz(i)=accz(i)+accz_cluster
	   accz(i)=0.

	   ENDDO

	   PRINT *,'evolve>> Initial accelerations computed <<'

C First half step :
	CALL forward_half_step_v(nbodies,vx_filgal,vy_filgal,vz_filgal,
     &                           vxp,vyp,vzp,
     &                           accx,accy,accz,accx_filgal,accy_filgal,
     &                           accz_filgal,dt)

	PRINT *,'evolve>> Half step complete <<'
	PRINT *,'evolve>> Initializing orbit and audit file <<'
	ncount=1
	OPEN(unit=13,file=orbit_file,status='unknown')
	WRITE(13,*)it,x_filgal,y_filgal,z_filgal
	OPEN(unit=14,file=audit_file,status='unknown')
	CALL audit(G,x_filgal,y_filgal,z_filgal,vx_filgal,vy_filgal,
     &             vz_filgal,sigma_cluster,core_radius_cluster,
     &             outer_radius_cluster,e_total,am_total)
	WRITE(14,*)it,time,am_total,e_total
	PRINT *,' '
	PRINT *,'evolve>> Beginning evolution <<'
	PRINT *,' '

C Main loop :
	DO it=1,nsteps
	   PRINT *,'evolve>> Step : ',it,' Time : ',time
	   IF(it.ne.1)THEN
             CALL cluster_acc(x_filgal,y_filgal,z_filgal,G,
     &                         core_radius_cluster,
     &                         outer_radius_cluster,sigma_cluster,
     &                         accx_filgal,
     &                         accy_filgal,accz_filgal)
	      DO i=1,nbodies
                 CALL filgal_acc(xp(i),yp(i),zp(i),x_filgal,
     &                           y_filgal,z_filgal,
     &                           G,lum_mass_filgal,
     &                           re_filgal,extent,
     &                           gamma_halo,accx(i),accy(i),accz(i))
                 CALL cluster_acc(xp(i),yp(i),zp(i),G,
     &                         core_radius_cluster,
     &                         outer_radius_cluster,sigma_cluster,
     &                         accx_cluster,
     &                         accy_cluster,accz_cluster)
	         accx(i)=accx(i)+accx_cluster
	         accy(i)=accy(i)+accy_cluster
C	         accz(i)=accz(i)+accz_cluster
	         accz(i)=0.
	         ENDDO
	      ENDIF

C Leap frog integration :
           CALL leap_frog(nbodies,x_filgal,y_filgal,z_filgal,
     &                    vx_filgal,
     &                    vy_filgal,vz_filgal,accx_filgal,
     &                    accy_filgal,accz_filgal,dt,time,
     &                    xp,yp,zp,vxp,vyp,vzp,
     &                    accx,accy,accz)

	   WRITE(13,*)it,x_filgal,y_filgal,z_filgal

           IF(ncount.eq.nout)THEN
	      CALL back_half_step_v(nbodies,vx_filgal,vy_filgal,vz_filgal,
     &                              vxp,vyp,vzp,
     &                              accx,accy,accz,accx_filgal,
     &                              accy_filgal,accz_filgal,dt)

	      PRINT *,'evolve>> Doing Audit <<'
	      CALL audit(G,x_filgal,y_filgal,z_filgal,vx_filgal,vy_filgal,
     &                   vz_filgal,sigma_cluster,core_radius_cluster,
     &                   outer_radius_cluster,e_total,am_total)
	      PRINT *,'evolve>---> Energy : ',e_total,' AM : ',am_total
	      WRITE(14,*)it,time,am_total,e_total

              WRITE(11,*)it,time
              DO i=1,nbodies
	         WRITE(11,*)i,xp(i),yp(i),zp(i),vxp(i),vyp(i),vzp(i)
	         ENDDO

              CALL forward_half_step_v(nbodies,vx_filgal,
     &                                 vy_filgal,vz_filgal,
     &                                 vxp,vyp,vzp,
     &                                 accx,accy,accz,
     &                                 accx_filgal,accy_filgal,
     &                                 accz_filgal,dt)
	      ncount=1
	      ELSE
	      ncount=ncount+1
	      ENDIF
	   ENDDO

	CLOSE(unit=11)
	CLOSE(unit=13)
	CLOSE(unit=14)
	PRINT *,'evolve>> Evolution complete <<'
	END

C*****************************************************************

C	Subroutine to compute orbital elements for
C	parabolic orbits in a modified halo field.

C*****************************************************************

	SUBROUTINE elements2(radius,peri_distance,gamma,sigma_v,
     &                       x1,y1,vx1,vy1)

	REAL m1,m2

	va=sqrt(2.)*sigma_v
	m1=0.
	p1=peri_distance

	m2=va*va*(radius-gamma*atan(radius/gamma))
	r1=radius/(1.+(m1/m2))
	x1=r1-p1
	y1=sqrt(r1*r1-x1*x1)
	vy1=-(m2/(m1+m2))*sqrt((2./r1)*m2/((y1*y1/(p1*p1))+1.))
	vx1=vy1*y1/p1

	END

C****************************************************************************

C	Subroutine to compute the acceleration due to the cluster

C****************************************************************************

	SUBROUTINE cluster_acc(x,y,z,G,
     &                         core_radius_cluster,
     &                         outer_radius_cluster,sigma_cluster,
     &                         accx,accy,accz)

	PARAMETER (epsilon=0.01)

	radius=sqrt(x*x+y*y+z*z)
	scale=radius/core_radius_cluster
	scale2=outer_radius_cluster/core_radius_cluster
	total_mass_cluster=2.*sigma_cluster*sigma_cluster*
     &                     core_radius_cluster*
     &                     (scale2-atan(scale2))
	contained_mass=2.*sigma_cluster*sigma_cluster*
     &                 core_radius_cluster*
     &                 (scale-atan(scale))
	IF(radius.gt.outer_radius_cluster)THEN
	   the_mass=total_mass_cluster
	   ELSE
	   the_mass=contained_mass
	   ENDIF
	radial_acc=-G*the_mass/(radius*radius+epsilon*epsilon)
	accx=radial_acc*x/radius
	accy=radial_acc*y/radius
	accz=radial_acc*z/radius
	RETURN
	END

C*****************************************************************************

C	Subroutine to compute the acceleration due to the filament galaxy

C*****************************************************************************

	SUBROUTINE filgal_acc(x,y,z,x_filgal,y_filgal,z_filgal,
     &                        G,lum_mass_filgal,
     &                        re_filgal,extent,
     &                        gamma_halo,accx,accy,accz)

	PARAMETER (epsilon=0.01)
	REAL lum_mass_filgal

	pi=2.*asin(1.)

	xx=x-x_filgal
	yy=y-y_filgal
	zz=z-z_filgal

	radius=sqrt(xx*xx+yy*yy+zz*zz)

	radial_acc_plummer=-G*lum_mass_filgal*radius/(
     &                      radius*radius+re_filgal*re_filgal)**1.5

	scale=radius/gamma_halo
	scale2=extent
	outer_radius=extent*gamma_halo
	sigma_halo=pi*G*lum_mass_filgal/(32.*re_filgal)
	contained_mass_halo=2.*sigma_halo*sigma_halo*gamma_halo*
     &                 (scale-atan(scale))
	total_mass=2.*sigma_halo*sigma_halo*gamma_halo*
     &                 (scale2-atan(scale2))
	IF(radius.gt.outer_radius)THEN
	   radial_acc_halo=-G*total_mass/(radius*radius)
	   ELSE
           radial_acc_halo=-G*contained_mass_halo/
     &                      (radius*radius+epsilon*epsilon)
	   ENDIF
C	radial_acc_halo=0.
	total_radial_acc=radial_acc_plummer+radial_acc_halo

	accx=total_radial_acc*xx/radius
	accy=total_radial_acc*yy/radius
	accz=total_radial_acc*zz/radius

	RETURN

	END

C******************************************************************************

C	Subroutine to advance the velocities by half a time step

C******************************************************************************

	SUBROUTINE forward_half_step_v(nbodies,vx_filgal,vy_filgal,vz_filgal,
     &                                 vxp,vyp,vzp,
     &                                 accx,accy,accz,accx_filgal,accy_filgal,
     &                                 accz_filgal,dt)

	DIMENSION vxp(1),vyp(1),vzp(1)
	DIMENSION accx(1),accy(1),accz(1)

	vx_filgal=vx_filgal+accx_filgal*dt/2.
	vy_filgal=vy_filgal+accy_filgal*dt/2.
	vz_filgal=vz_filgal+accz_filgal*dt/2.

	DO i=1,nbodies
	   vxp(i)=vxp(i)+accx(i)*dt/2.
	   vyp(i)=vyp(i)+accy(i)*dt/2.
	   vzp(i)=vzp(i)+accz(i)*dt/2.
	   ENDDO

	RETURN

	END

C****************************************************************************

C	Subroutine to advance the positions and velocities one time step

C***************************************************************************

	SUBROUTINE leap_frog(nbodies,x_filgal,y_filgal,z_filgal,
     &                       vx_filgal,
     &                       vy_filgal,vz_filgal,accx_filgal,
     &                       accy_filgal,accz_filgal,dt,time,
     &                       xp,yp,zp,vxp,vyp,vzp,
     &                       accx,accy,accz)

	DIMENSION xp(1),yp(1),zp(1)
	DIMENSION vxp(1),vyp(1),vzp(1)
	DIMENSION accx(1),accy(1),accz(1)

	vx_filgal=vx_filgal+dt*accx_filgal
	vy_filgal=vy_filgal+dt*accy_filgal
	vz_filgal=vz_filgal+dt*accz_filgal

	x_filgal=x_filgal+dt*vx_filgal
	y_filgal=y_filgal+dt*vy_filgal
	z_filgal=z_filgal+dt*vz_filgal

C	print *,x_filgal,y_filgal,v_filgal
C	print *,vx_filgal,vy_filgal,vz_filgal
C	print *,accx_filgal,accy_filgal,accz_filgal
	DO i=1,nbodies
	   IF(i.eq.10)THEN
	       xx=xp(i)-x_filgal
	       yy=yp(i)-y_filgal
	       zz=zp(i)-z_filgal
	       vx=vxp(i)-vx_filgal
	       vy=vyp(i)-vy_filgal
	       vz=vzp(i)-vz_filgal
	       rad=sqrt(xx*xx+yy*yy+zz*zz)
	       radv=(xx*vx/rad)+(yy*vy/rad)+(zz*vz/rad)
	       radacc=(xx*accx(i)/rad)+(yy*accy(i)/rad)+
     &	              (zz*accz(i)/rad)
C	       print *,i,rad,radv,radacc
C	       print *,i,accx(i),accy(i),accz(i)
C	       print *,i,vx,vy,vz
C	       print *,i,xx,yy,zz
C	   print *,i,vxp(i),accx(i),dt,vx_filgal
	       ENDIF
	   vxp(i)=vxp(i)+dt*accx(i)
	   vyp(i)=vyp(i)+dt*accy(i)
	   vzp(i)=vzp(i)+dt*accz(i)
	   xp(i)=xp(i)+dt*vxp(i)
	   yp(i)=yp(i)+dt*vyp(i)
	   zp(i)=zp(i)+dt*vzp(i)
	   ENDDO

	time=time+dt
	RETURN
	END

C******************************************************************************

C	Subroutine to backup the velocities by half a time step

C******************************************************************************

	SUBROUTINE back_half_step_v(nbodies,vx_filgal,vy_filgal,vz_filgal,
     &                                 vxp,vyp,vzp,
     &                                 accx,accy,accz,accx_filgal,accy_filgal,
     &                                 accz_filgal,dt)

	DIMENSION vxp(1),vyp(1),vzp(1)
	DIMENSION accx(1),accy(1),accz(1)

	vx_filgal=vx_filgal-accx_filgal*dt/2.
	vy_filgal=vy_filgal-accy_filgal*dt/2.
	vz_filgal=vz_filgal-accz_filgal*dt/2.

	DO i=1,nbodies
	   vxp(i)=vxp(i)-accx(i)*dt/2.
	   vyp(i)=vyp(i)-accy(i)*dt/2.
	   vzp(i)=vzp(i)-accz(i)*dt/2.
	   ENDDO

	RETURN
	END

C**************************************************************************

C	Subroutine to check the energy and angular momentum of the
C	filament galaxies orbit in the cluster field.

C**************************************************************************

	SUBROUTINE audit(G,x_filgal,y_filgal,z_filgal,vx_filgal,vy_filgal,
     &                   vz_filgal,sigma_cluster,core_radius_cluster,
     &                   outer_radius_cluster,e_total,am_total)

	amx=y_filgal*vz_filgal-z_filgal*vy_filgal
	amy=z_filgal*vx_filgal-x_filgal*vz_filgal
	amz=x_filgal*vy_filgal-y_filgal*vx_filgal

	am_total=sqrt(amx*amx+amy*amy+amz*amz)

	radius=sqrt(x_filgal*x_filgal+y_filgal*y_filgal+z_filgal*z_filgal)
	scale=radius/core_radius_cluster
	scale2=outer_radius_cluster/core_radius_cluster

	t1=(1./scale)*atan(scale)-1.
	t2=0.5*alog((1.+scale*scale)/(1.+scale2*scale2))
	cluster_mass=2.*sigma_cluster*sigma_cluster*core_radius_cluster*(
     &                  scale2-atan(scale2))
        t3=cluster_mass/(core_radius_cluster*(scale2-atan(scale2)))
	IF(radius.lt.outer_radius_cluster)THEN
           cluster_potential=t3*(t1+t2)
	   ELSE
	   cluster_potential=-G*cluster_mass/radius
	   ENDIF
	energy_kinetic=0.5*(vx_filgal*vx_filgal+vy_filgal*vy_filgal+
     &                      vz_filgal*vz_filgal)

	e_total=cluster_potential+energy_kinetic

	RETURN
	END

ext: for


! { dg-do run }
!
! PR 47024: [OOP] STORAGE_SIZE (for polymorphic types): Segfault at run time
! PR 47189: [OOP] calling STORAGE_SIZE on a NULL-initialized class pointer
! PR 47194: [OOP] EXTENDS_TYPE_OF still returns the wrong result if the polymorphic variable is unallocated
!
! Contributed by Tobias Burnus <burnus@gcc.gnu.org>

type t
  integer(kind=4) :: a
end type

class(t), pointer :: x => null()
class(t), allocatable :: y

if (storage_size(x)/=32) STOP 1
if (storage_size(y)/=32) STOP 2

allocate(y)

if (storage_size(y)/=32) STOP 3

deallocate(y)

if (storage_size(y)/=32) STOP 4

end


text:
! RUN: %S/test_errors.sh %s %t %f18 -fopenmp
! OpenMP Version 4.5
! 2.15.4.1 copyin Clause
! A list item that appears in a copyin clause must be threadprivate

program omp_copyin

  integer :: i
  integer, save :: j, k
  integer :: a(10), b(10)

  !$omp threadprivate(j, k)

  j = 20
  k = 10

  !$omp parallel do copyin(j, k)
  do i = 1, 10
    a(i) = k + i
    b(i) = j + i
  end do
  !$omp end parallel do

  print *, a, b

end program omp_copyin

max_stars_repo_name: Machiry/checkedc-clang


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: How to create a new BUFR message by cloning
!               an existing message.
!
!
program bufr_clone
   use eccodes
   implicit none
   integer       :: i, iret
   integer       :: infile, outfile
   integer       :: ibufr_in
   integer       :: ibufr_out

   ! Open source file
   call codes_open_file(infile, '../../data/bufr/syno_multi.bufr', 'r')

   ! Open target file
   call codes_open_file(outfile, 'bufr_clone_test_f.clone.bufr', 'w')

   ! The first bufr message is loaded from file,
   ! ibufr is the bufr id to be used in subsequent calls
   call codes_bufr_new_from_file(infile, ibufr_in, iret)

   ! Create several clones of this message and alter them
   ! in different ways
   do i = 1, 3

      ! Clone the current handle
      call codes_clone(ibufr_in, ibufr_out)

      ! This is the place where you may wish to modify the clone
      ! E.g. we change the bufrHeaderCentre
      call codes_set(ibufr_out, 'bufrHeaderCentre', 222)

      ! Write cloned messages to a file
      call codes_write(ibufr_out, outfile)

      ! Release the clone's handle
      call codes_release(ibufr_out)
   end do

   ! Release the original handle
   call codes_release(ibufr_in)

   call codes_close_file(infile)
   call codes_close_file(outfile)

end program bufr_clone

max_stars_repo_name: mtasaka/eccodes

max_issues_repo_name: mtasaka/eccodes

max_forks_repo_name: mtasaka/eccodes


text:
      subroutine rteinit

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    This subroutine reads in the areas associated with files processed with
!!    the recday, recepic, recmon and recyear commands, calculates subbasin
!!    areas, calculates reach and hydrograph node drainage areas.

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    da_km       |km2           |area of the watershed in square kilometers
!!    icodes(:)   |none          |routing command code:
!!                               |0 = finish       9 = save
!!                               |1 = subbasin    10 = recday
!!                               |2 = route       11 = reccnst
!!                               |3 = routres     12 = structure
!!                               |4 = transfer    13 =
!!                               |5 = add         14 = saveconc
!!                               |6 = rechour     15 =
!!                               |7 = recmon
!!                               |8 = recyear
!!    ihouts(:)   |none          |For ICODES equal to
!!                               |0: not used
!!                               |1,2,3,5,7,8,10,11: hydrograph storage
!!                               |                     location number
!!                               |4: departure type
!!                               |   (1=reach)
!!                               |   (2=reservoir)
!!                               |9: hydrograph storage location of data to
!!                               |   be printed to event file
!!                               |14:hydrograph storage location of data to
!!                               |   be printed to saveconc file
!!    inum1s(:)   |none          |For ICODES equal to
!!                               |0: not used
!!                               |1: subbasin number
!!                               |2: reach number
!!                               |3: reservoir number
!!                               |4: reach or res # flow is diverted from
!!                               |5: hydrograph storage location of 1st
!!                               |   dataset to be added
!!                               |7,8,9,10,11,14: file number
!!    inum2s(:)   |none          |For ICODES equal to
!!                               |0,1,7,8,10,11: not used
!!                               |2,3: inflow hydrograph storage location
!!                               |4: destination type
!!                               |   (1=reach)
!!                               |   (2=reservoir)
!!                               |5: hydrograph storage location of 2nd
!!                               |   dataset to be added
!!                               |9,14:print frequency
!!                               |   (0=daily)
!!                               |   (1=hourly)
!!    inum3s(:)   |none          |For ICODES equal to
!!                               |0,1,2,3,5,7,8,10,11: not used
!!                               |4: destination number. Subbasin, reach or
!!                               |   reservoir receiving water
!!                               |9: print format
!!                               |   (0=normal, fixed format)
!!                               |   (1=txt format for AV interface,recday)
!!    inum4s(:)   |none          |For ICODES equal to
!!                               |0,2,3,5,7,8,9,10,11: not used
!!                               |1: GIS code printed to output file
!!                               |   (optional)
!!                               |4: rule code governing transfer of water
!!                               |   (1=fraction transferred out)
!!                               |   (2=min flow left)
!!                               |   (3=exact amount transferred)
!!    mhyd        |none          |maximum number of hydrograph nodes
!!    rnum1s(:)   |none          |For ICODES equal to
!!                               |0,1,3,5,9: not used
!!                               |2: Fraction of flow in channel
!!                               |4: amount of water transferred (as
!!                               |   defined by INUM4S)
!!                               |7,8,10,11: drainage area in square kilometers
!!                               |   associated with the record file
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    hyd_dakm(:) |km**2         |total drainage area of hydrograph in square
!!                               |kilometers
!!    rch_dakm(:) |km**2         |total drainage area contributing to flow at
!!                               |the outlet (pour point) of the reach in
!!                               |square kilometers
!!    sub_fr(:)   |km2/km2       |fraction of total watershed area contained
!!                               |in subbasin
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    icd         |none          |variable to hold value for icodes(:)
!!    idum        |none          |hydrograph node counter
!!    iht         |none          |variable to hold value for ihouts(:)
!!    inm1        |none          |variable to hold value for inum1s(:)
!!    inm2        |none          |variable to hold value for inum2s(:)
!!    inm3        |none          |variable to hold value for inum3s(:)
!!    inm4        |none          |variable to hold value for inum4s(:)
!!    rnm1        |none          |variable to hold value for rnum1s(:)
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

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

      use parm

      integer :: icd, inm1, inm2, inm3, inm4, iht
      real :: rnm1

!!    calculate watershed area in hectares
      da_ha = da_km * 100.
!!    calculate subbasin fraction of watershed area
      do idum = 1, subtot
        sub_fr(idum) = sub_km(idum) / da_km
      end do
!!    calculate HRU fraction of watershed area
      do idum = 1, nhru
        hru_dafr(idum) = hru_km(idum) / da_km
        isb = hru_sub(idum)
	  if (idplt(idum) > 0) then
         if (bio_e(idplt(idum)) > 1.e-6) then
          subfr_nowtr(isb) = subfr_nowtr(isb) + hru_fr(idum)
	   end if
	  end if
      end do

!!    read in areas associated with .fig record files
      do idum = 1, mhyd
        icd = 0
        iht = 0
        rnm1 = 0.

        icd = icodes(idum)
        rnm1 = rnum1s(idum)
        iht = ihouts(idum)

        select case (icd)
        case (0)   !! FINISH command
          exit

        case (6)   !! RECHOUR command
          hyd_dakm(iht) = rnm1

        case (7)   !! RECMON command
          hyd_dakm(iht) = rnm1

        case (8)   !! RECYEAR command
          hyd_dakm(iht) = rnm1

        case (10)  !! RECDAY command
          hyd_dakm(iht) = rnm1

        case (11)  !! RECCNST command
          hyd_dakm(iht) = rnm1
        end select
      end do

!!    calculate subbasin areas
!!    calculate reach and hydrograph drainage areas
      do idum = 1, mhyd
        icd = 0
        iht = 0
        inm1 = 0
        inm2 = 0
        inm3 = 0
        rnm1 = 0.
        inm4 = 0

        icd = icodes(idum)
        iht = ihouts(idum)
        inm1 = inum1s(idum)
        inm2 = inum2s(idum)
        inm3 = inum3s(idum)
        rnm1 = rnum1s(idum)
        inm4 = inum4s(idum)

        select case (icd)
        case (0)   !! FINISH command
          exit

        case (1)   !! SUBBASIN command
          hyd_dakm(iht) = da_km * sub_fr(inm1)
          rch_dakm(inm1) = hyd_dakm(iht)
          if (iroutunit == 1) then
            j1 = hru1(inm1) - 1
          do j = 1, hrutot(inm1)
            jseq = j1 + j
            hru_fr(jseq) = 0.
            hru_dafr(jseq) = 0.
          end do
          end if

        case (2)   !! ROUTE command
          hyd_dakm(iht) = hyd_dakm(inm2)
          rch_dakm(inm1) = hyd_dakm(inm2)

        case (3)   !! ROUTRES command
          hyd_dakm(iht) = hyd_dakm(inm2)

        case (5)   !! ADD command
          hyd_dakm(iht) = hyd_dakm(inm1) + hyd_dakm(inm2)

       case (17)   !! ROUTEUNIT command
          hyd_dakm(iht) = rnm1
          j1 = hru1(inm2) - 1
          do j = 1, hrutot(inm2)
            jseq = j1 + j
            hru_fr(jseq) = hru_fr(jseq) + hru_rufr(inm1,j) *
     &                                daru_km(inm2,inm1) / sub_km(inm2)
            hru_dafr(jseq) = hru_dafr(jseq) + hru_rufr(inm1,j) *
     &                                      daru_km(inm2,inm1) / da_km
          end do

        case (18)   !! ROUTEOVER command
          hyd_dakm(iht) = hyd_dakm(inm2)

        end select
      end do

      return
      end


text:
program sqrt_2
    use, intrinsic :: iso_fortran_env, only : SP => REAL32, DP => REAL64
    implicit none
    real(kind=SP) :: sqr_root
    real(kind=DP) :: sqr_root_dp

    sqr_root = sqrt(2.0_SP)
    sqr_root_dp = sqrt(2.0_DP)
    print '(A, L)', 'sqrt(2)**2 == 2 (SP): ', sqr_root**2 == 2.0_SP
    print '(A, L)', 'sqrt(2)**2 =~ 2 (SP): ', abs(sqr_root**2 - 2.0_SP) <= 1.0e-6_SP
    print '(A, L)', 'sqrt(2)**2 == 2 (DP): ', sqr_root_dp**2 == 2.0_DP
    print '(A, L)', 'sqrt(2)**2 =~ 2 (DP): ', abs(sqr_root_dp**2 - 2.0_DP) <= 1.0e-14_DP
end program sqrt_2

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
c=======================================================================
c
c     subroutine ALLOCRBRFR
c
c     Risk Budgeting allocation strategy with a risk-free asset
c
c     w(i), i=1,...,n the riskly assets (portfolio weights)
c     w(0), the risk-free asset
c
c     Max[ w'*rho - k*w'*Q*w]
c      s.t.
c     sqrt[w'*cov_i*w] <= sigma(i)  (i-th risk budgeting constraint)
c     C*w <= b                      (general linear constraints)
c     Cinf(i) <= w(i) <= Csup(i)    (riskly assets lower/upper bounds)
c     Cinf(0) <= w(0) <= Csup(0)    (risk-free asset lower/upper bounds)
c
c        Q   : covariance matrix
c        rho : expected returns
c        rfr : risk-free rate
c
c-----------------------------------------------------------------------
      SUBROUTINE allocrbrfr ( n, cov, Q, kappa, rho, rfr,
     &                        neq, nin, ccst, bcst,
     &                        cinf, csup, cinfrfr, csuprfr,
     &                        nbc, class, sigma,
     &                        iwork, dwork, wopt, wrf, lambda, info )
c-----------------------------------------------------------------------
c
c     INPUT
c       n      : portfolio size                                  integer
c       cov    : covariance matrix (n*n)                          double
c       Q      : matrix (n*n)                                     double
c       kappa  : risk aversion coefficient                        double
c       rho    : expected returns vector (n)                      double
c       rfr    : expected risk-free rate                          double
c       neq    : number equality constraints                     integer
c       nin    : number inequality constraints                   integer
c       ccst   : matrix of constraints (nasset*(neq+nin))         double
c       bcst   : vector initial of constraints (neq+nin)          double
c       cinf   : lower bound (n)                                  double
c       csup   : upper bound (n)                                  double
c       cinfrfr: risk-free lower bound                            double
c       csuprfr: risk-free upper bound                            double
c       nbc    : number of risk budgeting constraints (>=1)      integer
c       class  : block definition (n)                            integer
c       sigma  : volatilities budgets constraints (nbc)           double
c
c     WORKSPACE
c       iwork  : 2*nbc + 5*n + neq + 2*tnin + 7                  integer
c              = 2*nbc + 5*n + neq + 2*nin + 11
c       dwork  : nbc*(nbc+21)/2
c              + n*(8*n+2*neq+2*nin+26)+5*neq+7*nin+30            double
c
c     OUTPUT
c       wopt   : optimal portfolio (n)                            double
c       wrf    : optimal risk-free asset                          double
c       lambda : optimal lambda (nbc)                             double
c       info   : diagnostic argument                             integer
c
c     CALL
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin, nbc
      INTEGER class(*)
      DOUBLE PRECISION wrf, rfr, kappa, cinfrfr, csuprfr
      DOUBLE PRECISION cov(*), Q(*), rho(*), cinf(*), csup(*),
     &                 ccst(n,*), bcst(*), wopt(*), sigma(*), lambda(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER i, infotmp, iter, piw, pdw
      DOUBLE PRECISION ncinfrfr, ncsuprfr, coef, sum, vol, EPS, ZERO
      PARAMETER (ZERO=0.E0, EPS=1.E-15)
c
c     external subroutines
      EXTERNAL multivolrfr, PVX2, SEV, testlinearcst
c
c-----------------------------------------------------------------------
c
c     initializations
      info = 0          ! diagnostic argument
      iter = 0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piw = 1
c
c     Total size of dwork array = 2*nbc + 5*n + neq + 2*nin + 11
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pdw = 1
c
c     Total size of dwork array = nbc*(nbc+21)/2
c                               + n*(8*n+3*neq+3*nin+28)
c                               + 5*neq + 7*nin + 31
c
c-----------------------------------------------------------------------
c
c     case sigma(1)=0
      IF (nbc .EQ. 1) THEN
        IF (ABS(sigma(1)) .LT. EPS) THEN
            CALL IVX ( n, wopt, ZERO )
            IF (csuprfr .GT. (1.0-1.E-8)) THEN
                wrf = 1.0
            ELSE
                wrf  = cinfrfr
                info = 110
            ENDIF
            RETURN
        ENDIF
      ENDIF
c
      CALL multivolrfr ( n, cov, Q, kappa, rho, rfr,
     &                   neq, nin, ccst, bcst,
     &                   cinf, csup, cinfrfr, csuprfr,
     &                   nbc, class, sigma,
     &                   iwork(piw), dwork(pdw),
     &                   wopt, wrf, lambda, info )
      infotmp = info
      GOTO 10
c
c     case info = -120 (risk budget cst. not repected)
 10   iter = iter + 1
      IF (iter .GT. 2) GOTO 200
      IF (info .EQ. -120) THEN
        IF (nbc .EQ. 1) THEN
c           portfolio volatility
            info = 0
            CALL OVTMCV(n, cov, wopt, vol)
            IF (vol .GT. 1.E-15) THEN
                vol = SQRT(vol)
                coef = sigma(1)/vol
            ELSE
                coef = 0.0
            ENDIF
            CALL PVX2 ( n, wopt, coef )
            CALL SEV ( n, wopt, sum )
            wrf = 1.0 - sum
            IF ((wrf .LT. (cinfrfr - 1.E-8)) .OR.
     &          (wrf .GT. (csuprfr + 1.E-8))) THEN
                info = 110
                GOTO 200
            ENDIF
            CALL testlinearcst ( n, wopt, neq, nin, ccst, bcst, cinf,
     &                           csup, dwork(pdw), info )
            IF (info .LT. 0) THEN
                info = -115
                GOTO 50
            ENDIF
            GOTO 200
        ELSE
            GOTO 60
        ENDIF

 50     coef = 0.0 ! risk-free rate = 0%
        CALL multivolrfr ( n, cov, Q, kappa, rho, coef,
     &                     neq, nin, ccst, bcst,
     &                     cinf, csup, cinfrfr, csuprfr,
     &                     nbc, class, sigma,
     &                     iwork(piw), dwork(pdw),
     &                     wopt, wrf, lambda, info )
        GOTO 10

 60     IF (cinfrfr .EQ. 0) THEN
            ncinfrfr = 0.0
            ncsuprfr = 0.0
            GOTO 100
        ENDIF
        IF (csuprfr .LT. 1) THEN
            ncinfrfr = csuprfr
            ncsuprfr = csuprfr
            GOTO 100
        ENDIF
100     CALL multivolrfr ( n, cov, Q, kappa, rho, rfr,
     &                     neq, nin, ccst, bcst,
     &                     cinf, csup, ncinfrfr, ncsuprfr,
     &                     nbc, class, sigma,
     &                     iwork(piw), dwork(pdw),
     &                     wopt, wrf, lambda, info )
        IF (info .EQ. -120) THEN
            info = 110
        ENDIF
      ENDIF
      GOTO 200
200   IF ((info.NE.1001).AND.(info.NE.-108).AND.(info.NE.-109)) THEN
        DO i = 1,nbc
            IF (ABS(lambda(i)) .LT. 1.E-14) THEN
                info = 111
            ENDIF
            IF (ABS(lambda(i)) .GT. 1.E+14) THEN
                info = 110
            ENDIF
        ENDDO
      ENDIF
      IF ((info .EQ. 0).AND.(infotmp.NE.-120).AND.(infotmp.NE.0)) THEN
        info = infotmp
      ENDIF

      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
Subroutine graduk(ix, iy, iz, gradxk, gradyk, gradzk)
  Parameter (maxx=20, maxz=24)
  Parameter (rho0=0.168)
  Common /dd/rho(-maxx:maxx, -maxx:maxx, -maxz:maxz), rhop(-maxx:maxx, -maxx:maxx, -maxz:maxz), rhon(-maxx:maxx, -maxx:maxx, -maxz:maxz)
  Common /ss/inout(20)
  Save
  rxplus = rho(ix+1, iy, iz)
  rxmins = rho(ix-1, iy, iz)
  ryplus = rho(ix, iy+1, iz)
  rymins = rho(ix, iy-1, iz)
  rzplus = rho(ix, iy, iz+1)
  rzmins = rho(ix, iy, iz-1)
  gradxk = (rxplus-rxmins)/2.
  gradyk = (ryplus-rymins)/2.
  gradzk = (rzplus-rzmins)/2.
  Return
End Subroutine graduk

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
      Function KRandomC(J,M,K,L,R)        ! Random, Row Rule
      KRandomC = 0
      if (R .LE. .5) KRandomC = 1
      Return
      End     ! Random Col Rule

max_stars_repo_name: Axelrod-Python/TourExec

max_issues_repo_name: Axelrod-Python/TourExec

max_forks_repo_name: Axelrod-Python/TourExec


text:
One of several Banks in Davis. A branch used to be located in Anderson Plaza where Chamonix Coffee used to reside. It may have moved after being Bank Robberies robbed 3 times.  Users/StevenDaubert says the location that would constantly get robbed had an ATMs ATM that would give out 10s as well as 20s!

Several years ago UC Davis announced a joint project with US Bank for a full US Bank branch inside the MU, along with six new ATMs. http://www.news.ucdavis.edu/search/news_detail.lasso?id9248 US Bank made a deal with UC Davis to place the US Bank logo on all student ID cards in exchange for the exclusive right for the next generation of student ID cards to be used as US Bank ATM cards. Reg Cards were redesigned and renamed to the new Aggie Card. The new card has the ability to be used as a US Bank ATM card if you sign up for a US Bank account. It can be used for PINbased purchases, but it cannot be used for signature purchases like a credit card.

The new US Bank branch opened up in the MU, however the branch was the scene of repeated protests and sitins. US Bank decided to abandon the branch and leave campus on February 28, 2012. Some people felt that corporate advertising had gone too far with US Bank logos being plastered onto every single student ID card. Other people felt that US Banks goal was to crowd out all other bank and credit union business and get as much student business as they could. However, the apparent goal of being the campus bank monopoly was thwarted by repeated protests. Students didnt appreciate US Bank and US Bank decided to abandon their campus bank branch since it was not profitable for them. Chancellor Katehi has mentioned that several other banks immediately started asking for the spot when US Bank left.

Misdemeanor charges have been filed against 12 people connected to the oncampus U.S. Bank protests. They are known as The Davis Dozen.

http://davisvanguard.org/index.php?optioncom_content&viewarticle&id5347:bankprotesterspleadnotguiltyrejectpleaagreement&catid74:courtwatch&Itemid100 Davis Dozen Student Protesters Plead Not Guilty  To be Defended by Legendary Defense Lawyer Tony Serra, aka the Semantic Samurai, May 10, 2012

The only place to find a US Bank ATM in Davis is at the US Bank branch in downtown.

US Bank used to have ATMs located at the MU ATM Wall, the ARC, the Health Sciences Complex, the MU Games Area, Silo, Mondavi Center, and at the former US Bank branch at the MU. However, in 2012 US Bank got rid of all of their ATMs on campus. This now forces US Bank customers to get cash at other ATMs which amounts to a total fee of about $4.50.

All US Bank ATMs are part of the MoneyPass network. Citibank customers may use MoneyPass and US Banks ATMs without a fee. Yolo Federal Credit Union also is part of the MoneyPass network and allows free use of US Bank ATMs.

Warning: US Bank conducts a hard credit inquiry when you open accounts with them. This lowers your credit score and can damage your credit rating. So be careful about opening accounts with them.

If you open up a checking account, your credit score drops. If you open up a savings account too, your score drops even lower. Opening up a credit card account also makes your score drop lower. These credit inquiries add up and your credit score just keeps going lower. A lower credit score is bad news.

I have always found it annoying that US Bank requires parent signatures for minors to withdraw money or cash checks larger than $20. Users/StephenHudson
Since when did they start this?  Users/StevenDaubert
Since my parents purposefully set up a limited savings subaccount for me that I was only a cosigner of. My mistake. Users/StephenHudson

Kerrian(sp?) is simply one of the coolest women around.  Gave me great pursecleaning tips.  Thanks, girl!

20080922 14:17:30 nbsp Hassan is the man when it comes to customer service! Users/StevenDaubert

20100929 08:53:34 nbsp By far the ugliest bank in Davis, both outside and inside. Their hours (closing at 5 p.m.) are pretty inconvenient. And from what I understand (if I overheard the teller correctly) the Aggie Card is NOT a debit card, its an ATM card. That means you can use it to get money out of a US Bank ATM, but you cant use it to buy food at the CoHo or books at the bookstore. Useless. Users/jsogul

20120125 01:19:22 nbsp I went to the MU branch today, and theyve hired a security guard, which was a surprise to me. Apparently they were shut down by sitin protests so many days out of the last two weeks that (at least for now) they are keeping out anyone who doesnt have a US Bank account and card. Users/StephenHudson
  Since the university decided to emblazon the US Bank logo across my student ID, I wonder how they are going to weed protesters out. M

20120314 14:30:43 nbsp Why are they protesting? Users/halfhiphop
  Why dont you ask them?
         i did, didnt really hear them.

20120601 01:23:32 nbsp I think its important to acknowledge that much of this information is not accurate. The university told US Bank to leave because UCD is trying to stay out of headlines after the pepperspraying incident. The students brought a lot of attention to US Bank and the university didnt want to take any action to allow US Bank customers access into the branch (as protestors made it impossible for US Bank to conduct business). That being said, I do not support corporations on the university campus. US Bank is just another corporation. Many people look at banks as being something evil, but just like all businesses, they offer products our culture demands. For those interested in really making change, work on banishing all corporations from conducting business on public university property and really make the university represent public education, not privatization. The blame should really be placed on the university, not US Bank. Users/Smitty

US Bank also recently got rid of all their ATMs on campus  there was one in the ARC and one recently added at the MU  now i have to pay the $4.50 fee to get out some cash.  WTF?

20130326 18:58:20 nbsp The staff is very rude and unprofessional Users/NobleMason