text:
      FUNCTION F_MAIN()

      IMPLICIT NONE
      INCLUDE 'IPhreeqc.f.inc'
      INTEGER(KIND=4) id

      INTEGER(KIND=4)    r
      INTEGER(KIND=4)    c
      INTEGER(KIND=4)    t
      REAL(KIND=8)       d
      CHARACTER(LEN=80)  s

      INTEGER(KIND=4) F_MAIN
      INTEGER(KIND=4) TestGetSet
      INTEGER(KIND=4) TestGetSetName

      INTEGER(KIND=4) EXIT_SUCCESS
      PARAMETER (EXIT_SUCCESS=0)

      INTEGER(KIND=4) EXIT_FAILURE
      PARAMETER (EXIT_FAILURE=1)

      EXTERNAL GetDumpFileOn
      EXTERNAL SetDumpFileOn

      EXTERNAL GetDumpStringOn
      EXTERNAL SetDumpStringOn

      EXTERNAL GetErrorFileOn
      EXTERNAL SetErrorFileOn

      EXTERNAL GetLogFileOn
      EXTERNAL SetLogFileOn

      EXTERNAL GetOutputFileOn
      EXTERNAL SetOutputFileOn

      EXTERNAL GetSelectedOutputFileOn
      EXTERNAL SetSelectedOutputFileOn

      EXTERNAL GetDumpFileName
      EXTERNAL SetDumpFileName

      EXTERNAL GetOutputFileName
      EXTERNAL SetOutputFileName

      id = CreateIPhreeqc()
      IF (id.LT.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Dump
      IF (TestGetSet(id,GetDumpFileOn,SetDumpFileOn).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Dump string
      IF (TestGetSet(id,GetDumpStringOn,SetDumpStringOn).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Dump filename
      IF (TestGetSetName(id,GetDumpFileName,SetDumpFileName).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Error
      IF (TestGetSet(id,GetErrorFileOn,SetErrorFileOn).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Log
      IF (TestGetSet(id,GetLogFileOn,SetLogFileOn).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Output
      IF (TestGetSet(id,GetOutputFileOn,SetOutputFileOn).NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Output filename
      IF (TestGetSetName(id,GetOutputFileName,SetOutputFileName)
     &     .NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

C     Selected output
      IF (TestGetSet(id,GetSelectedOutputFileOn,SetSelectedOutputFileOn)
     &     .NE.0) THEN
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

      IF (LoadDatabase(id, "phreeqc.dat").NE.0) THEN
         CALL OutputErrorString(id)
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

      IF (RunFile(id, "ex2").NE.0) THEN
         CALL OutputErrorString(id)
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

      DO 20 r=0,GetSelectedOutputRowCount(id)
         DO 10 c=1,GetSelectedOutputColumnCount(id)
            IF (GetSelectedOutputValue(id,r,c,t,d,s).NE.IPQ_OK) THEN
               CALL OutputErrorString(id)
               F_MAIN = EXIT_FAILURE
               RETURN
            END IF
 10      CONTINUE
 20   CONTINUE

      IF (DestroyIPhreeqc(id).NE.0) THEN
         CALL OutputErrorString(id)
         F_MAIN = EXIT_FAILURE
         RETURN
      END IF

      F_MAIN = EXIT_SUCCESS
      RETURN

      END FUNCTION F_MAIN

      FUNCTION TestGetSet(id,getFunc,setFunc)

      IMPLICIT NONE
      INCLUDE 'IPhreeqc.f.inc'
      INTEGER(KIND=4) id
      INTEGER(KIND=4) TestGetSet

      LOGICAL(KIND=4) getFunc
      INTEGER(KIND=4) setFunc

      INTEGER(KIND=4) EXIT_SUCCESS
      PARAMETER (EXIT_SUCCESS=0)

      INTEGER(KIND=4) EXIT_FAILURE
      PARAMETER (EXIT_FAILURE=1)

      IF (getFunc(id)) THEN
         TestGetSet = EXIT_FAILURE
         RETURN
      END IF

      IF (setFunc(id,.TRUE.).NE.IPQ_OK) THEN
         TestGetSet = EXIT_FAILURE
         RETURN
      END IF

      IF (.NOT.getFunc(id)) THEN
         TestGetSet = EXIT_FAILURE
         RETURN
      END IF

      IF (setFunc(id,.FALSE.).NE.IPQ_OK) THEN
         TestGetSet = EXIT_FAILURE
         RETURN
      END IF

      TestGetSet = EXIT_SUCCESS
      RETURN

      END FUNCTION TestGetSet

      FUNCTION TestGetSetName(id,getFuncName,setFuncName)

      IMPLICIT NONE
      INCLUDE 'IPhreeqc.f.inc'
      INTEGER(KIND=4) id
      INTEGER(KIND=4) TestGetSetName

      EXTERNAL getFuncName
      INTEGER(KIND=4) setFuncName

      INTEGER(KIND=4) EXIT_SUCCESS
      PARAMETER (EXIT_SUCCESS=0)

      INTEGER(KIND=4) EXIT_FAILURE
      PARAMETER (EXIT_FAILURE=1)

      CHARACTER(LEN=80) FILEN

      CALL getFuncName(id,FILEN)

      IF (setFuncName(id,'ABCDEFG').NE.IPQ_OK) THEN
         TestGetSetName = EXIT_FAILURE
         WRITE(*,*) "FAILURE"
         RETURN
      END IF

      CALL getFuncName(id,FILEN)
      IF (.NOT.LLE('ABCDEFG', FILEN)) THEN
         TestGetSetName = EXIT_FAILURE
         WRITE(*,*) "FAILURE"
         RETURN
      END IF

      IF (setFuncName(id,'XYZ').NE.IPQ_OK) THEN
         TestGetSetName = EXIT_FAILURE
         WRITE(*,*) "FAILURE"
         RETURN
      END IF

      CALL getFuncName(id,FILEN)
      IF (.NOT.LLE('XYZ', FILEN)) THEN
         TestGetSetName = EXIT_FAILURE
         WRITE(*,*) "FAILURE"
         RETURN
      END IF

      TestGetSetName = EXIT_SUCCESS
      RETURN

      END FUNCTION TestGetSetName

ext: f

max_stars_repo_name: usgs-coupled/webmod

max_issues_repo_name: usgs-coupled/webmod

max_forks_repo_name: usgs-coupled/webmod


text:
! Program to test the PRODUCT intrinsic
program testproduct
   implicit none
   integer, dimension (3, 3) :: a
   integer, dimension (3) :: b
   logical, dimension (3, 3) :: m, tr
   character(len=12) line

   a = reshape ((/1, 2, 3, 4, 5, 6, 7, 8, 9/), (/3, 3/));

   b = product (a, 1)

   tr = .true.

   if (any(b .ne. (/6, 120, 504/))) call abort

   write (line, 9000) product(a,1)
   if (line .ne. '   6 120 504') call abort

   if (product (a) .ne. 362880) call abort

   write (line, 9010) product(a)
   if (line .ne. '362880') call abort

   m = .true.
   m(1, 1) = .false.
   m(2, 1) = .false.

   b = product (a, 2, m)
   if (any(b .ne. (/28, 40, 162/))) call abort

   b = product (a, 2, m .and. tr)
   if (any(b .ne. (/28, 40, 162/))) call abort

   write (line, 9000) product(a, 2, m)
   if (line .ne. '  28  40 162') call abort

   if (product (a, mask=m) .ne. 181440) call abort

   if (product (a, mask=m .and. tr) .ne. 181440) call abort

   write (line, 9010) product(a, mask=m)
   if (line .ne. '181440') call abort

9000 format (3I4)
9010 format (I6)
end program

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
*
*     ------------------------------------------------------------------
*       P A C K
*     ------------------------------------------------------------------
*  Subroutine written by Bin LIU
*
*  Rules for encoding
*  1. All blanks deleted
*  2. If Qi=1, omit Qi
*  3. If Qi=1 or Qi>=4l+1, omit ALFAi
*  4. If i=1 or (Qi=4l+2 and i<>m), insert '.'; else _BETAi.
*
      SUBROUTINE PACK (M, EL, Q, COUPLE, STR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

          INTEGER FULL,Q(*),CONST
          CHARACTER*3 EL(*),COUPLE(*)
          CHARACTER   STR*64,CH1*1,CH3*3
*
*  FULL   :  4l+2
*  CONST  :   constant for converting lowercase to uppercase
*  CH*    :  temporary variables
*
          CONST = ICHAR('a') - ICHAR('A')
          STR=' '
          J = 1
*
*   -----  begin to encode  -----
*
          DO 100 I=1,M
            IF (Q(I) .EQ. 0) GO TO 100
            IF (EL(I)(1:1) .EQ. ' ') THEN
                STR(J:J+1) = EL(I)(2:3)
                K = 2
            ELSE
                STR(J:J+2) = EL(I)
                K = 3
                IF (EL(I)(3:3) .EQ. ' ') K = 2
            END IF
            CH1 = STR(J+1:J+1)
            IF ((CH1.GE.'A') .AND. (CH1.LE.'Z'))
     :            STR(J+1:J+1)=CHAR(ICHAR(CH1)+CONST)
            FULL=4*LVAL(CH1)+2
*
*  -----  convert Qi into character  -----
*
            J = J+K
            N=Q(I)
            IF (N .GT. 14) STOP 'Too many electrons'
*
*  -----  If Qi<>1, add Qi
*           If Qi<4l+1, add TERMi for the shell -----
*
            IF (N .NE. 1) THEN
               STR(J:J) = '('
               J = J + 1
               IF (N .GT. 9) THEN
                    STR(J:J) = '1'
                    J = J + 1
                    N = N - 10
               END IF
               STR(J:J+1) = CHAR(ICHAR('0')+N)//')'
               J = J + 2
               IF (N .LT. FULL-1 .AND. M .NE. 1) THEN
                   CH3=COUPLE(I)
                   CH1= CH3(2:2)
                   IF (CH1.GE.'a' .AND. CH1.LE.'z')
     :                 CH3(2:2) = CHAR(ICHAR(CH1)-CONST)
                   STR(J:J+2) = CH3
                   J = J+3
               ENDIF
            ENDIF
*
*  -----  If i=1 or Qi=4l+2 and i<>m,
*           insert '.'; else _RESULTANTi.  -----
*
            IF ((I .NE. 1 .AND. N .NE. FULL) .OR. M .EQ. I) THEN
                CH3=COUPLE(M+I-1)
                CH1 = CH3(2:2)
                IF (CH1.GE.'a' .AND. CH1.LE.'z')
     :                  CH3(2:2) = CHAR(ICHAR(CH1)-CONST)
                K = 2
                IF (M .EQ. 1) K = 3
                STR(J:J+K) = '_'//CH3(1:K)
                J = J + K + 1
            ENDIF
            IF (I .NE. M) THEN
                STR(J:J)='.'
                J = J + 1
            ENDIF
100       CONTINUE
*
*>>>>>    Because of a compiler error on the SUN, the following is
*         needed to have the string printed correctly.
*         STR = STR(1:J)
          RETURN
        END


!---------------------------------------------------------------------!
! OWNER: Ithaca Combustion Enterprise, LLC                            !
! COPYRIGHT: © 2012, Ithaca Combustion Enterprise, LLC                !
! LICENSE: BSD 3-Clause License (The complete text of the license can !
!  be found in the `LICENSE-ICE.txt' file included in the ISAT-CK7    !
!  source directory.)                                                 !
!---------------------------------------------------------------------!

subroutine ell_aff_pr( n, c, gg, m, d, T, cp, ggp )

!  Given an ellipsoid E (in R^n), and an m-dimensional affine space A,
!  determine the ellipsoid P(E) which is the orthogonal projection of E onto A,
!  for  1 <= m <= n.

!  E is defined by:  || G^T * (x-c) || <= 1, where the array gg
!  contains the lower Cholesky triangle G in packed format.

!  A is defined by { x | x = d + T * s }, where T is an n x m orthogonal
!  matrix, and s is any m-vector.

!  The orthogonal projection of the general point x is given by:
!  P(x) = { x | d + T * s,  s = T^T * ( x - d ) }.

!  The orthogonal projection of E onto A is given by:
!  P(E) = { x | d + T * s,  || GP^T * ( s - cp ) || <= 1 },  where the
!  array  ggp contains GP in packed format.

!  Method:
!    cp = T^T * ( c - d )
!    F  = G^{-1} * T = U * [ S ; 0 ] * V^T
!    B  = V * S^{-1}
!    GP * GP^T = B * B^T

implicit none

integer, parameter      :: k_dp = kind(1.d0)
integer, intent(in)     :: n, m
real(k_dp), intent(in)  :: c(n), gg((n*(n+1))/2), d(n), T(n,m)
real(k_dp), intent(out) :: cp(m),  ggp((m*(m+1))/2)

integer, parameter ::  ludiag = 0
integer    :: i, j, k, info, lwork
real(k_dp) :: G(n,n), F(n,m), S(m), B(m,m), VT(m,m), work(5*n*n+20*n)

!  check dimensions
if( n < 1 ) then
   write(ludiag,*)'ell_aff_pr: n < 1, n = ', n
   stop
endif

lwork = 5*n*n+20*n  ! amount of work space for SVD

if( m > n ) then
   write(ludiag,*)'ell_aff_pr: m > n, [m n] = ', m, n
   stop
elseif( m < 1 ) then
   write(ludiag,*)'ell_aff_pr: m < 1, m = ', m
   stop
endif

!  form cp = T^T * ( c - d )

cp = matmul( (c-d) , T )

!  unpack G

G = 0.d0

k  = 0
do j = 1, n
   do i = j, n
      k = k + 1
	  G(i,j) = gg(k)
   end do
end do

!  form F = G^{-1} * T
F = T
call dtrsm('L','L','N','N',n,m,1.d0,G,n,F,n)

!  form S and V^T in:  F = U * [ S ; 0 ] * V^T

call dgesvd('N','A',n,m,F,n,S,B,m,VT,m,work,lwork,info)

if( info /= 0 ) then
   write(ludiag,*)'ell_aff_pr: SVD failed, info = ', info
   stop
elseif( S(m) <= 0.d0 ) then
   write(ludiag,*)'ell_aff_pr: zero singular value, S = ', S
   stop
endif

!  form B = V * S^{-1}

do j = 1, m
   B(1:m,j) = VT(j,1:m) / S(j)
end do

!  form GP:  GP * GP^T = B * B^T

call ell_bbt2chol( m, B, ggp )

return
end subroutine ell_aff_pr


text:
!!# FUNCTION MODULE: <STRTIME>
MODULE FUN_STRTIME

!!## PURPOSE
!! Returns a runtime in a human-readable, compact format.

!!## USAGE
!
! WRITE(*,*)STRTIME(dt)
!
!! Note: <STRTIME(dt)> is 7 characters---always.

!!## ACCESS AND IMPLICIT
IMPLICIT NONE
PRIVATE

!!## DETAILS
!! The passed time <dt> (or $\Delta t$) is handled
!! in the following way:
!!
!! @ if $\Delta t < t_{us} $ or $\Delta t >t_{max}$ then use
!!   exponential format with seconds, otherwise
!!
!! @ if $\Delta t >= t_{s}$, then use seconds, otherwise
!!
!! @ if $\Delta t >= t_{ms}$, then use milliseconds, otherwise
!!
!! @ if $\Delta t >= t_{us}$, then use microseconds.
!!

!!## LOCAL VARIABLES
REAL :: tmax=60.
REAL :: ts=1.0
REAL :: tms=1.e-3
REAL :: tus=1.e-6

!!## PUBLIC ACCESS LIST
PUBLIC :: TEST_STRTIME,STRTIME

!!## PROCEDURES
CONTAINS

!!### FUNCTION <STRTIME>
FUNCTION STRTIME(dt)

!!#### REQUIRED INPUT
REAL,INTENT(IN) :: dt

!!#### REQUIRED OUTPUT
CHARACTER(7) :: STRTIME

!!#### LOCAL VARIABLES
CHARACTER(10) :: s

!!--begin--

!use seconds in exponential format
IF( dt>tmax .OR. dt<tus )THEN

 WRITE(s,"(Es10.2)")ABS(dt)
 s = ADJUSTL(s)

 !need both exponent digits
 IF( ABS(dt)<1.e-9 .OR. ABS(dt)>1.e+9 )THEN
  STRTIME = s(1:2)//"E"//s(6:6)//s(7:8)//"s"
 !just need one
 ELSE
  STRTIME = s(1:3)//"E"//s(6:6)//s(8:8)//"s"
 END IF
!use seconds
ELSE IF( dt>=ts )THEN

 WRITE(s,"(f10.4)")ABS(dt)
 s = ADJUSTL(s)
 STRTIME = s(1:6)//"s"

!use millseconds
ELSE IF( dt>=tms )THEN

 WRITE(s,"(f9.3)")ABS(dt*1.E3)
 s = ADJUSTL(s)
 STRTIME = s(1:5)//"ms"

!use micrseconds
ELSE IF( dt>=tus )THEN

 WRITE(s,"(f9.3)")ABS(dt*1.E6)
 s = ADJUSTL(s)
 STRTIME = s(1:5)//"us"

END IF

!!--end--
END FUNCTION

SUBROUTINE TEST_STRTIME()
REAL :: r
INTEGER :: i,j

!!--begin--
WRITE(*,"(a12,2x,a12)")"r (seconds)","STRTIME(r)"
DO i=1,10
 CALL RANDOM_NUMBER(r)
 r=r/1.E7
 DO j=1,10
  r = r*10.
 WRITE(*,"(e12.5,2x,a)")r,STRTIME(r)
END DO

END DO

!!--end--
END SUBROUTINE

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
module default_btc_factory_module
    use :: util_api, only : &
            string, &
            dictionary, &
            add_prefix_to_priorities

    use :: tensor_contraction_api, only : &
            tensor_contraction_factory

    use :: batched_tensor_contraction_dev, only : &
            btc_factory, &
            batched_tensor_contraction, &
            simple_loops_btc_driver

    implicit none
    private

    public :: default_btc_factory

    type, extends(btc_factory) :: default_btc_factory
    contains
        procedure :: create_from_key => create_from_key
        procedure :: get_available_btc_drivers => get_available_btc_drivers
        procedure :: build => build
        procedure, private :: count_available_drivers => count_available_drivers
    end type default_btc_factory

    character(len=*), parameter :: btc_driver_key = "btc_driver"
    character(len=*), parameter :: default_btc_driver = "simple_loops"
    character(len=*), parameter :: prefix = "btc-"
contains
    subroutine build(this, btc, contraction, options, priorities)
        class(default_btc_factory), intent(in) :: this
        class(batched_tensor_contraction), intent(inout) :: btc
        type(string), intent(in) :: contraction
        type(dictionary), intent(in), optional :: options
        type(string), dimension(:), intent(in), optional :: priorities

        type(string), dimension(:), allocatable :: local_priorities

        local_priorities = add_prefix_to_priorities(prefix, priorities)

        select type (btc)
        type is (simple_loops_btc_driver)
            call tensor_contraction_factory%create(btc%driver, contraction%char_array, &
                    options=options, priorities=local_priorities)
        class default
            error stop "default_btc_factory::build:Unknown btc implementation."
        end select
    end subroutine build

    subroutine create_from_key(this, btc, key)
        class(default_btc_factory), intent(in) :: this
        class(batched_tensor_contraction), allocatable, intent(inout) :: btc
        type(string), intent(in) :: key

        select case (key%char_array)
        case ("simple_loops")
            btc = simple_loops_btc_driver()
        case default
            error stop "default_btc_factory::create_from_key:Not a valid btc driver: "//key%char_array
        end select
    end subroutine create_from_key

    function get_available_btc_drivers(this) result(drivers)
        class(default_btc_factory), intent(in) :: this
        type(string), dimension(:), allocatable :: drivers

        allocate(drivers(this%count_available_drivers()))

        drivers(1) = "simple_loops"
    end function get_available_btc_drivers

    integer function count_available_drivers(this)
        class(default_btc_factory), intent(in) :: this

        count_available_drivers = 1
    end function count_available_drivers
end module default_btc_factory_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


! { dg-do run }
! { dg-options "-fdec-static -finit-local-zero" }
!
! Test AUTOMATIC and STATIC attributes.
!
subroutine assert(s, i1, i2)
  implicit none
  integer, intent(in)      :: i1, i2
  character(*), intent(in) :: s
  if (i1 .ne. i2) then
    print *, s, ": expected ", i2, " but was ", i1
    STOP 1
  endif
endsubroutine assert

function f (x, y)
  implicit none
  integer f
  integer, intent(in)  :: x, y
  integer              :: a    ! only a can actually be saved
  integer, automatic   :: c    ! should actually be automatic
  save

  ! a should be incremented by x every time and saved
  a = a + x
  f = a

  ! c should be zeroed every time, therefore equal y
  c = c + y
  call assert ("f%c", c, y)
  return
endfunction

implicit none
integer :: f

! Should return static value of a; accumulates x
call assert ("f()", f(1,3), 1)
call assert ("f()", f(1,4), 2)
call assert ("f()", f(1,2), 3)

end


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
!  Sorting routines
!  References:
!  D. Knuth
!  The Art of Computer Programming, vol. 3
!  Addison-Wesley
!
!  Aho, Hopcroft, Ullman
!  Data Structures and Algorithms
!  Addison-Wesley
!
module psb_i_hsort_x_mod
  use psb_const_mod
  use psb_e_hsort_mod
  use psb_m_hsort_mod

  type psb_i_heap
    integer(psb_ipk_) :: last, dir
    integer(psb_ipk_), allocatable    :: keys(:)
  contains
    procedure, pass(heap) :: init       => psb_i_init_heap
    procedure, pass(heap) :: howmany    => psb_i_howmany
    procedure, pass(heap) :: insert     => psb_i_insert_heap
    procedure, pass(heap) :: get_first  => psb_i_heap_get_first
    procedure, pass(heap) :: dump       => psb_i_dump_heap
    procedure, pass(heap) :: free       => psb_i_free_heap
  end type psb_i_heap

  type psb_i_idx_heap
    integer(psb_ipk_) :: last, dir
    integer(psb_ipk_), allocatable    :: keys(:)
    integer(psb_ipk_), allocatable :: idxs(:)
  contains
    procedure, pass(heap) :: init       => psb_i_idx_init_heap
    procedure, pass(heap) :: howmany    => psb_i_idx_howmany
    procedure, pass(heap) :: insert     => psb_i_idx_insert_heap
    procedure, pass(heap) :: get_first  => psb_i_idx_heap_get_first
    procedure, pass(heap) :: dump       => psb_i_idx_dump_heap
    procedure, pass(heap) :: free       => psb_i_idx_free_heap
  end type psb_i_idx_heap

contains

  subroutine psb_i_init_heap(heap,info,dir)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none
    class(psb_i_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)            :: info
    integer(psb_ipk_), intent(in), optional   :: dir

    info = psb_success_
    heap%last=0
    if (present(dir)) then
      heap%dir = dir
    else
      heap%dir = psb_sort_up_
    endif
    select case(heap%dir)
    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
      ! ok, do nothing
    case default
      write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_'
      heap%dir = psb_sort_up_
    end select
    call psb_ensure_size(psb_heap_resize,heap%keys,info)

    return
  end subroutine psb_i_init_heap

  function psb_i_howmany(heap) result(res)
    implicit none
    class(psb_i_heap), intent(in) :: heap
    integer(psb_ipk_) :: res
    res  = heap%last
  end function psb_i_howmany

  subroutine psb_i_insert_heap(key,heap,info)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none

    integer(psb_ipk_), intent(in)              :: key
    class(psb_i_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)                       :: info

    info = psb_success_
    if (heap%last < 0) then
      write(psb_err_unit,*) 'Invalid last in heap ',heap%last
      info = heap%last
      return
    endif

    call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize)
    if (info /= psb_success_) then
      write(psb_err_unit,*) 'Memory allocation failure in heap_insert'
      info = -5
      return
    end if
    call psi_insert_heap(key,&
         & heap%last,heap%keys,heap%dir,info)

    return
  end subroutine psb_i_insert_heap

  subroutine psb_i_heap_get_first(key,heap,info)
    implicit none

    class(psb_i_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)     :: info
    integer(psb_ipk_), intent(out)       :: key

    info = psb_success_

    call psi_heap_get_first(key,&
         & heap%last,heap%keys,heap%dir,info)

    return
  end subroutine psb_i_heap_get_first

  subroutine psb_i_dump_heap(iout,heap,info)

    implicit none
    class(psb_i_heap), intent(in) :: heap
    integer(psb_ipk_), intent(out)    :: info
    integer(psb_ipk_), intent(in)     :: iout

    info = psb_success_
    if (iout < 0) then
      write(psb_err_unit,*) 'Invalid file '
      info =-1
      return
    end if

    write(iout,*) 'Heap direction ',heap%dir
    write(iout,*) 'Heap size      ',heap%last
    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
         & (size(heap%keys)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else
      write(iout,*) heap%keys(1:heap%last)
    end if
  end subroutine psb_i_dump_heap

  subroutine psb_i_free_heap(heap,info)
    implicit none
    class(psb_i_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)           :: info

    info=psb_success_
    if (allocated(heap%keys)) deallocate(heap%keys,stat=info)

  end subroutine psb_i_free_heap

  subroutine psb_i_idx_init_heap(heap,info,dir)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none
    class(psb_i_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)            :: info
    integer(psb_ipk_), intent(in), optional   :: dir

    info = psb_success_
    heap%last=0
    if (present(dir)) then
      heap%dir = dir
    else
      heap%dir = psb_sort_up_
    endif
    select case(heap%dir)
    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
      ! ok, do nothing
    case default
      write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_'
      heap%dir = psb_sort_up_
    end select

    call psb_ensure_size(psb_heap_resize,heap%keys,info)
    call psb_ensure_size(psb_heap_resize,heap%idxs,info)
    return
  end subroutine psb_i_idx_init_heap

  function psb_i_idx_howmany(heap) result(res)
    implicit none
    class(psb_i_idx_heap), intent(in) :: heap
    integer(psb_ipk_) :: res
    res  = heap%last
  end function psb_i_idx_howmany

  subroutine psb_i_idx_insert_heap(key,index,heap,info)
    use psb_realloc_mod, only : psb_ensure_size
    implicit none

    integer(psb_ipk_), intent(in)              :: key
    integer(psb_ipk_), intent(in)                        :: index
    class(psb_i_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)                       :: info

    info = psb_success_
    if (heap%last < 0) then
      write(psb_err_unit,*) 'Invalid last in heap ',heap%last
      info = heap%last
      return
    endif

    call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize)
    if (info == psb_success_) &
         & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize)
    if (info /= psb_success_) then
      write(psb_err_unit,*) 'Memory allocation failure in heap_insert'
      info = -5
      return
    end if
    call psi_idx_insert_heap(key,index,&
         & heap%last,heap%keys,heap%idxs,heap%dir,info)

    return
  end subroutine psb_i_idx_insert_heap

  subroutine psb_i_idx_heap_get_first(key,index,heap,info)
    implicit none

    class(psb_i_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)       :: index
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_), intent(out)           :: key

    info = psb_success_

    call psi_idx_heap_get_first(key,index,&
         & heap%last,heap%keys,heap%idxs,heap%dir,info)

    return
  end subroutine psb_i_idx_heap_get_first

  subroutine psb_i_idx_dump_heap(iout,heap,info)

    implicit none
    class(psb_i_idx_heap), intent(in) :: heap
    integer(psb_ipk_), intent(out)    :: info
    integer(psb_ipk_), intent(in)     :: iout

    info = psb_success_
    if (iout < 0) then
      write(psb_err_unit,*) 'Invalid file '
      info =-1
      return
    end if

    write(iout,*) 'Heap direction ',heap%dir
    write(iout,*) 'Heap size      ',heap%last
    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
         & (size(heap%keys)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else    if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
         & (size(heap%idxs)<heap%last))) then
      write(iout,*) 'Inconsistent size/allocation status!!'
    else
      write(iout,*) heap%keys(1:heap%last)
      write(iout,*) heap%idxs(1:heap%last)
    end if
  end subroutine psb_i_idx_dump_heap

  subroutine psb_i_idx_free_heap(heap,info)
    implicit none
    class(psb_i_idx_heap), intent(inout) :: heap
    integer(psb_ipk_), intent(out)           :: info

    info=psb_success_
    if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
    if ((info == psb_success_).and.(allocated(heap%idxs))) &
         & deallocate(heap%idxs,stat=info)

  end subroutine psb_i_idx_free_heap

end module psb_i_hsort_x_mod


text:
BloodSource is a nonprofit blood center that collects donated blood from Volunteer volunteers at their many locations as well as from regular blood drives conducted on campus and in the nearby area.  The Davis BloodSource center is also somehow extremely well connected, being one of only two blood providers that was able to get usable blood to New York the very day of September 11th.

BloodSource is clean and bright and frequently decorated with seasonal themes. The staff are all very pleasant and accommodating, clients donating plasma have a nice selection of movies to watch while they donate. The snacks after donation are plentiful.

Appointments are easy to make. Simply call them up on the phone and arrive on time. If theyre not busy, they will accept walkins, but making an appointment is more efficient. Many potential donors wrongly feel discouraged from donating because they erroneously assume they dont qualify (due to Drugs illicit drug use, for example). If youre one of these potential donors, dont be afraid: simply walk in and ask if you can donate.

There is currently a shortage of blood in the country. This is due in part to one of the most reliable sources of public blood, the military, having to keep its donors to replenish their own reserves. Aside from the military, college and high school students are one of the most important demographics for blood donations, since they tend to be the most eligible and healthy. Save a life. Give blood!

http://www.bloodsource.org/HowYouCanHelp/DonateBlood/What_to_expect.cfm There are a number of other things you can do, or donate

They also take platelet and plasma donations by appointment. Read the information on their Web page, give them a call, or ask one of the staff the next time you donate whole blood. and after some tests.

Causeway Classic Blood Drive

The Causeway Classic is an annual, friendly competition started in 2008 between UC Davis and Sacramento State to register the most donors and collect the most pints over two days. The Causeway Classic blood drive is held in late Novemeber and coincides with other competitions between the two schools including the football and hockey games. The winning school enjoys their name on a trophy and bragging rights for a year.

Year Winner Footnote(Source: http://www.bloodsource.org/News/NewsReleases/131007CausewayClassic.aspx)
2011 Sac State
2012 Sac State
2013 Sac State
 list incomplete

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
 Beta Site

About
This picture is no longer representative because it is a small pond by now.  Its really too bad.
Unfortunately, Beta Site is a long since dead project, but the remnants of it still possibly exist where it once was.
Because of its questionable legality, the creators are not able to reveal the location.

20051022 16:58:11 nbsp Look into a http://cmms.cat.com/cmms/specgraphics/C055701.pdf 446D backhoe by www.cat.com Caterpillar. Then find a civil engineer who knows about gunite and rebar. Users/JasonAller

20051111 08:51:05 nbsp Seasonally, digging a hole in the fall around here results in a hole with rainwater in it. Good luck. Users/JasonAller

20051111 09:13:37 nbsp well, we hope to have the door + weather lining installed before the first big rain.  Thanks for the tip Though! Users/AlexNorris

20051111 13:20:08 nbsp It sounds like a cool idea, but is this legal? Users/GeoffJohnson

20051111 14:58:50 nbsp In response to GeoffJohnson, yes. Users/AlexNorris

20051113 20:49:13 nbsp Funny, I know of many male youth which have dug big holes. Must be some primal thing. Question is, what do you do with it once its finished? Users/SteveDavison

20051113 21:48:25 nbsp Beta Site is pretty much going to be a room or two undergroud. A place for sleepovers, hanging out, music writing, and acoustic shows. Its well underway now, as progress this weekend was very great. Users/JohnDudek

20051116 08:34:43 nbsp Daang N. Mouth u guys should have told me bout this. Let me know how I can help and Ill see what can be done. Users/JoshPelz

20051117 22:08:47 nbsp Mmm... If I had time and supplies Id help. If you guys need like publicizing or anything you think I could help with, give me a ring... or a comment Users/SamanthaSpada

20051118 17:47:05 nbsp publicizing it probably isnt the best plan right now,  davis wiki people can know about it sure, but if every one knew, and someone else figured out where it was, then they would do drugs there till the police found out.  so lets keep it on a downlowish basis right now. Users/DanielGonzales

20051118 22:06:07 nbsp Gah, not what i meant by publicizing. I meant like... try to get you guys help from people who you know... dont spend all their time on the wiki like we do Users/SamanthaSpada

20051118 23:03:33 nbsp Oh, then thats cool have fun with that, you should come look at our hole some time.  and we need to take a picture or us standing in it Users/DanielGonzales

20051119 18:05:15 nbsp Im sure the citys planning commission would have something to say about this. Users/ApolloStumpy

20051120 10:38:39 nbsp I dont think its OK to delete this page even ifyou created it. Is this hole legal? Does it fit the zoning ordinances? What about the environmental impact? I think the Davis community deserves to know. Users/ApolloStumpy

20051120 10:41:38 nbsp weve checked the environmental impact, Apollo, no problem there. As far as legality, where its at, Im sure its not illegal. Im just trying to keep this lowdown so nobody we dont like shows up. Thanks for your concern. Users/JohnDudek

20051120 12:03:36 nbsp What were worried about is that other holes in davis get taken over by people who use them for smoken the reefer and drinkin the bud, and then the police have to close it down.  We dont want that to happen to us, because this hole is a totally Strait Edge hole.  Thanks for your concern. Users/DanielGonzales

20051120 12:42:01 nbsp cmon, seriously...do you think some methhead is going to get on DavisWiki, see this, and automatically know where the hole is? If some druggie comes across your hole and uses it as his drug lair, it wont be because he heard of it here. Users/ApolloStumpy

20051120 13:04:41 nbsp well, we know some people who, if they get wind of this, WILL use it to smoke and stuff, so were just watching out. Users/JohnDudek

20051121 14:31:30 nbsp Whats this? Holes are now illegal because they could be used for someone to sit in and smoke weed? Why not make air illegal too? (whoops) Users/SteveDavison

20051121 22:02:11 nbsp Maybe air should be illegal... I hear it can give you cancer. Users/SamanthaSpada

20051121 22:14:40 nbsp If I had an underground lair, Id smoke weed in it.  What a waste of a good hole. Users/SummerSong

20051121 22:23:36 nbsp Sweet, this sounds like one of my childhood dreams come true. Users/RoyWright

20051122 00:45:04 nbsp so is this a hole or a tunnel? Seems to me a tunnel/subterranean lair would not be looked upon kindly by the city, whether you were smoking pot in it or having bible studies. Users/ApolloStumpy

20051122 00:58:57 nbsp HaHa. The comments conversation alone is awesome. Users/MichaelGiardina

20051122 09:25:17 nbsp it is a hole with a short tunnel. Clarified! Users/JohnDudek

20051122 11:22:45 nbsp photo! Users/ArlenAbraham

20051127 20:25:33 nbsp one word: wiki:WikiPedia:Buttress buttresses Users/ArlenAbraham

20051127 20:40:35 nbsp i say: fun in theory, dumb in reality. and where the heck is it, seriously.. in the pic it looks like its made out of ice. Users/AubreyJohnson

20051128 07:34:17 nbsp maybe it is... Users/AlexNorris

20051128 17:39:13 nbsp Guess what probably flooded like all hell today? Yeah... Beta Site. Probably. I havent checked yet, cant wait to slosh in that eight foot mud pit this week. Users/JohnDudek

20051128 21:36:55 nbsp I hope someone goes and checks on it toay, I want to see pictures of flooded/ not flooded beta site. Users/DanielGonzales

20051208 18:09:13 nbsp I have a dozen or so books on how to build underground houses.  You might look into the $50 and Up Underground House Book.  The website is http://www.undergroundhousing.com/ here, but you really want to get the book. Users/JabberWokky

20051208 21:37:41 nbsp Could you add the photo back?  I think it was pretty cool. Users/PhilipNeustrom

20051210 17:57:47 nbsp haha, knew it! sorry though guys, ah well eh. thats the way the pumpkin rolls. Users/AubreyJohnson

20060109 15:57:30 nbsp Looks like at least a couple of cal/osha violations, what with no shoring and no confined entry rescue plan. Users/GrumpyoldGeek

20060109 16:02:16 nbsp Since this isnt exactly an employer, how does cal/osha apply? Users/JaimeRaba

20060109 16:51:40 nbsp Is is your legal opinion that cal/osha does not apply? Users/GrumpyoldGeek

20080706 13:05:56 nbsp Did anyone ever find out where this was? Users/JoePomidor

20090506 18:37:43 nbsp Sorry, Joe, nobody will ever know where Beta Site was, Im afraid. Users/JohnDudek


text:
      SUBROUTINE GWF1ETS1ALP(ISUM,ISUMI,LCIETS,LCETSR,LCETSX,LCETSS,
     &                       NCOL,NROW,NETSOP,IN,IOUT,IETSCB,IFREFM,
     &                       NPETS,IETSPF,NETSEG,LCPXDP,LCPETM,NSEGAR)
C
C-----VERSION 20100315 ERB
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR EVAPOTRANSPIRATION SEGMENTS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
  500 FORMAT(1X,/
     &1X,'ETS1 -- EVAPOTRANSPIRATION SEGMENTS PACKAGE, VERSION 1,',
     &     ' 5/2/2000',/,9X,'INPUT READ FROM UNIT ',I4)
  510 FORMAT(
     &1X,I5,' SEGMENTS DEFINE EVAPOTRANSPIRATION RATE FUNCTION')
  520 FORMAT(' EVAPOTRANSPIRATION RATE FUNCTION IS LINEAR')
  530 FORMAT(
     &' ERROR: EVAPOTRANSPIRATION RATE FUNCTION MUST CONTAIN AT',/,
     &' LEAST ONE SEGMENT -- STOP EXECUTION (GWF1ETS1ALP)')
  540 FORMAT(1X,'ILLEGAL ET OPTION CODE. SIMULATION ABORTING')
  550 FORMAT(1X,'OPTION 1 -- EVAPOTRANSPIRATION FROM TOP LAYER')
  560 FORMAT(1X,'OPTION 2 -- EVAPOTRANSPIRATION FROM ONE SPECIFIED',
     &   ' NODE IN EACH VERTICAL COLUMN')
  564 FORMAT(1X,'OPTION 3 -- EVAPOTRANSPIRATION FROM UPPERMOST ACTIVE ',
     &   'CELL')
!  566 FORMAT(1X,'OPTION 4 -- EVAPOTRANSPIRATION FROM UPPERMOST ACTIVE ',
!     &   'CELL IN STACK OF',/,1X,
!     &   'CONTIGUOUS ACTIVE CELLS THAT INCLUDES BOTTOMMOST ACTIVE ',
!     &   'CELL')
  570 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
  580 FORMAT(1X,I10,' ELEMENTS IN RX ARRAY ARE USED BY ETS')
  590 FORMAT(1X,I10,' ELEMENTS IN IR ARRAY ARE USED BY ETS')
C
C1------IDENTIFY PACKAGE.
      IETSPF=20
      WRITE(IOUT,500)IN
C
C     READ COMMENT LINE(S) (ITEM 0)
      CALL URDCOM(IN,IOUT,LINE)
C
C2------READ ET OPTION (NETSOP), UNIT OR FLAG FOR CELL-BY-CELL FLOW
C       TERMS (IETSCB), NUMBER OF PARAMETERS (NPETS), AND NUMBER OF
C       SEGMENTS (NETSEG) (ITEM 1)
      IF (IFREFM.EQ.0) THEN
        READ(LINE,'(4I10)') NETSOP,IETSCB,NPETS,NETSEG
      ELSE
        LLOC=1
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NETSOP,R,IOUT,IN)
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IETSCB,R,IOUT,IN)
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NPETS,R,IOUT,IN)
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NETSEG,R,IOUT,IN)
      ENDIF
C
C3------CHECK TO SEE THAT ET OPTION IS LEGAL.
!      IF (NETSOP.GE.1 .AND. NETSOP.LE.4) GO TO 10
      IF (NETSOP.GE.1 .AND. NETSOP.LE.3) GO TO 10
C
C3A-----OPTION IS ILLEGAL -- PRINT A MESSAGE & ABORT SIMULATION.
      WRITE(IOUT,540)
      CALL USTOP(' ')
C
C4------OPTION IS LEGAL -- PRINT THE OPTION CODE.
   10 CONTINUE
      IF (NETSOP.EQ.1) WRITE(IOUT,550)
      IF (NETSOP.EQ.2) WRITE(IOUT,560)
      IF (NETSOP.EQ.3) WRITE(IOUT,564) ! Add option 3 ERB 9/7/2006
!      IF (NETSOP.EQ.4) WRITE(IOUT,566) ! Add option 4 ERB 9/7/2006
C
C5------IF CELL-BY-CELL FLOWS ARE TO BE SAVED, THEN PRINT UNIT NUMBER.
      IF (IETSCB.GT.0) WRITE(IOUT,570) IETSCB
C
C-----PRINT NUMBER OF PARAMETERS TO BE USED
      CALL UPARARRAL(-1,IOUT,LINE,NPETS)
C
C     PRINT MESSAGE IDENTIFYING NUMBER OF SEGMENTS IN ET VS. HEAD CURVE
      IF(NETSEG.GT.1) THEN
        WRITE(IOUT,510) NETSEG
        NSEGAR = NETSEG - 1
      ELSEIF (NETSEG.EQ.1) THEN
        WRITE(IOUT,520)
        NSEGAR = 1
      ELSE
        WRITE(IOUT,530)
        CALL USTOP(' ')
      ENDIF
C
C6------ALLOCATE SPACE FOR THE ARRAYS ETSR, ETSX, ETSS, PXDP, AND PETM.
      IRK=ISUM
      LCETSR=ISUM
      ISUM=ISUM+NCOL*NROW
      LCETSX=ISUM
      ISUM=ISUM+NCOL*NROW
      LCETSS=ISUM
      ISUM=ISUM+NCOL*NROW
      IF (NETSEG.GT.1) THEN
        LCPXDP=ISUM
        ISUM=ISUM+NCOL*NROW*(NETSEG-1)
        LCPETM=ISUM
        ISUM=ISUM+NCOL*NROW*(NETSEG-1)
      ELSE
        LCPXDP=1
        LCPETM=1
      ENDIF
C
C7------ALLOCATE SPACE FOR LAYER INDICATOR ARRAY (IETS) EVEN IF ET
C7------OPTION IS NOT 2, TO AVOID ERROR OF ARRAY (IR) NOT LARGE ENOUGH
      LCIETS=ISUMI
      ISUMI=ISUMI+NCOL*NROW
C
C8------CALCULATE & PRINT AMOUNT OF SPACE USED BY ET PACKAGE.
      IRK=ISUM-IRK
      WRITE(IOUT,580)IRK
      IRK=NCOL*NROW
      WRITE(IOUT,590)IRK
C
C9------RETURN.
      RETURN
      END
C=======================================================================
      SUBROUTINE GWF1ETS1RPPD(IN,IOUT,NPETS,ITERP,INAMLOC)
C
C-----VERSION 20011120 ERB
C     ******************************************************************
C     READ EVAPOTRANSPIRATION SEGMENTS PARAMETER DEFINITIONS
C     ******************************************************************
C     Modified 11/21/2001 to support parameter instances - ERB
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*4 PTYP
C     ------------------------------------------------------------------
C
C-------READ NAMED PARAMETERS
      IF (ITERP.EQ.1) WRITE(IOUT,5) NPETS
    5 FORMAT(1X,//1X,I5,' Evapotranspiration segments parameters')
      IF (NPETS.GT.0) THEN
        DO 20 K=1,NPETS
C         UPARARRRP READS PARAMETER NAME AND DEFINITION (ITEMS 2 AND 3)
          CALL UPARARRRP(IN,IOUT,N,0,PTYP,ITERP,1,INAMLOC)
          IF(PTYP.NE.'ETS') THEN
            WRITE(IOUT,7)
    7       FORMAT(1X,'Parameter type must be ETS')
            CALL USTOP(' ')
          ENDIF
   20   CONTINUE
      ENDIF
C
C8------RETURN
   60 RETURN
      END
C=======================================================================
      SUBROUTINE GWF1ETS1RPSS(NETSOP,IETS,ETSR,ETSX,ETSS,DELR,DELC,NCOL,
     &                        NROW,IN,IOUT,IFREFM,NPETS,RMLT,IZON,
     &                        NMLTAR,NZONAR,IETSPF,NETSEG,PXDP,PETM,
     &                        NSEGAR)
C
C     VERSION 20000620 ERB
C     ******************************************************************
C     READ EVAPOTRANSPIRATION DATA, AND PERFORM SUBSTITUTION USING
C     PARAMETER VALUES IF ETS PARAMETERS ARE DEFINED
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*24 ANAME(6)
      DIMENSION IETS(NCOL,NROW),ETSR(NCOL,NROW),ETSX(NCOL,NROW),
     &          ETSS(NCOL,NROW),DELR(NCOL),DELC(NROW),
     &          RMLT(NCOL,NROW,NMLTAR),IZON(NCOL,NROW,NZONAR),
     &          PXDP(NCOL,NROW,NSEGAR),PETM(NCOL,NROW,NSEGAR)
C
      DATA ANAME(1) /'   ET LAYER INDEX (IETS)'/
      DATA ANAME(2) /'       ET SURFACE (ETSS)'/
      DATA ANAME(3) /' EVAPOTRANS. RATE (ETSR)'/
      DATA ANAME(4) /' EXTINCTION DEPTH (ETSX)'/
      DATA ANAME(5) /'EXTINCT. DEP. PROPORTION'/
      DATA ANAME(6) /'      ET RATE PROPORTION'/
C     ------------------------------------------------------------------
C
C1------READ FLAGS SHOWING WHETHER DATA FROM PREVIOUS STRESS PERIOS ARE
C       TO BE REUSED.
      IF (NETSEG.GT.1) THEN
        IF(IFREFM.EQ.0) THEN
          READ(IN,'(5I10)') INETSS,INETSR,INETSX,INIETS,INSGDF
        ELSE
          READ(IN,*) INETSS,INETSR,INETSX,INIETS,INSGDF
        ENDIF
      ELSE
        IF(NETSOP.EQ.2) THEN
          IF(IFREFM.EQ.0) THEN
            READ(IN,'(4I10)') INETSS,INETSR,INETSX,INIETS
          ELSE
            READ(IN,*) INETSS,INETSR,INETSX,INIETS
          ENDIF
        ELSE
          IF(IFREFM.EQ.0) THEN
            READ(IN,'(3I10)') INETSS,INETSR,INETSX
          ELSE
            READ(IN,*) INETSS,INETSR,INETSX
          ENDIF
        ENDIF
      ENDIF
C
C2------TEST INETSS TO SEE WHERE SURFACE ELEVATION COMES FROM.
      IF (INETSS.LT.0) THEN
C2A------IF INETSS<0 THEN REUSE SURFACE ARRAY FROM LAST STRESS PERIOD
        WRITE(IOUT,10)
   10   FORMAT(1X,/1X,'REUSING ETSS FROM LAST STRESS PERIOD')
      ELSE
C3-------IF INETSS=>0 THEN CALL MODULE U2DREL TO READ SURFACE.
        CALL U2DREL(ETSS,ANAME(2),NROW,NCOL,0,IN,IOUT)
      ENDIF
C
C4------TEST INETSR TO SEE WHERE MAX ET RATE COMES FROM.
      IF (INETSR.LT.0) THEN
C4A-----IF INETSR<0 THEN REUSE MAX ET RATE.
        WRITE(IOUT,20)
   20   FORMAT(1X,/1X,'REUSING ETSR FROM LAST STRESS PERIOD')
      ELSE
C5------IF INETSR=>0 CALL MODULE U2DREL TO READ MAX ET RATE.
        IF(NPETS.EQ.0) THEN
          CALL U2DREL(ETSR,ANAME(3),NROW,NCOL,0,IN,IOUT)
        ELSE
C    INETSR is the number of parameters to use this stress period
          CALL PRESET('ETS')
          WRITE(IOUT,30)
   30     FORMAT(1X,///1X,
     &        'ETSR array defined by the following parameters:')
          IF (INETSR.EQ.0) THEN
            WRITE(IOUT,35)
   35       FORMAT(' ERROR: When parameters are defined for the ETS',
     &      ' Package, at least one parameter',/,' must be specified',
     &      ' each stress period -- STOP EXECUTION (GWF1ETS1RPSS)')
            CALL USTOP(' ')
          ENDIF
          CALL UPARARRSUB2(ETSR,NCOL,NROW,0,INETSR,IN,IOUT,'ETS',
     &                     ANAME(3),'ETS',IETSPF,RMLT,IZON,NMLTAR,
     &                     NZONAR)
        ENDIF
C
C6------MULTIPLY MAX ET RATE BY CELL AREA TO GET VOLUMETRIC RATE
        DO 50 IR=1,NROW
          DO 40 IC=1,NCOL
            ETSR(IC,IR)=ETSR(IC,IR)*DELR(IC)*DELC(IR)
   40     CONTINUE
   50   CONTINUE
      ENDIF
C
C7------TEST INETSX TO SEE WHERE EXTINCTION DEPTH COMES FROM
      IF (INETSX.LT.0) THEN
C7A------IF INETSX<0 REUSE EXTINCTION DEPTH FROM LAST STRESS PERIOD
        WRITE(IOUT,60)
   60   FORMAT(1X,/1X,'REUSING ETSX FROM LAST STRESS PERIOD')
      ELSE
C8-------IF INETSX=>0 CALL MODULE U2DREL TO READ EXTINCTION DEPTH
        CALL U2DREL(ETSX,ANAME(4),NROW,NCOL,0,IN,IOUT)
      ENDIF
C
C9------IF OPTION(NETSOP) IS 2 THEN WE NEED AN INDICATOR ARRAY.
      IF (NETSOP.EQ.2) THEN
C10------IF INIETS<0 THEN REUSE LAYER INDICATOR ARRAY.
        IF (INIETS.LT.0) THEN
          WRITE(IOUT,70)
   70     FORMAT(1X,/1X,'REUSING IETS FROM LAST STRESS PERIOD')
        ELSE
C11------IF INIETS=>0 THEN CALL MODULE U2DINT TO READ INDICATOR ARRAY.
          CALL U2DINT(IETS,ANAME(1),NROW,NCOL,0,IN,IOUT)
        ENDIF
      ENDIF
C
C12------IF ET FUNCTION IS SEGMENTED PXDP AND PETM ARRAYS ARE NEEDED.
      IF (NETSEG.GT.1) THEN
C13------IF INSGDF<0 THEN REUSE PXDP AND PETM ARRAYS.
        IF (INSGDF.LT.0) THEN
          WRITE(IOUT,80)
   80     FORMAT(1X,/1X,
     &           'REUSING PXDP AND PETM FROM LAST STRESS PERIOD')
C14------IF INSGDF=>0 THEN CALL MODULE U2DREL TO READ PXDP AND PETM
C        ARRAYS.
        ELSE
          DO 90 I = 1,NETSEG-1
            WRITE(IOUT,100) I
            CALL U2DREL(PXDP(1,1,I),ANAME(5),NROW,NCOL,0,IN,IOUT)
            CALL U2DREL(PETM(1,1,I),ANAME(6),NROW,NCOL,0,IN,IOUT)
   90     CONTINUE
        ENDIF
      ENDIF
  100 FORMAT(/,' PXDP AND PETM ARRAYS FOR INTERSECTION ',I4,
     &' OF HEAD/ET RELATION:')
C
C15-----RETURN
      RETURN
      END
C=======================================================================
      SUBROUTINE GWF1ETS1FM(NETSOP,IETS,ETSR,ETSX,ETSS,RHS,HCOF,IBOUND,
     &                      HNEW,NCOL,NROW,NLAY,NETSEG,PXDP,PETM,
     &                      NSEGAR)
C
C-----VERSION 20100315 ERB
C     ******************************************************************
C        ADD EVAPOTRANSPIRATION TO RHS AND HCOF
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW, HH, SS, XX, DD, PXDP1, PXDP2
      DIMENSION IETS(NCOL,NROW), ETSR(NCOL,NROW), ETSX(NCOL,NROW),
     &          ETSS(NCOL,NROW), RHS(NCOL,NROW,NLAY),
     &          HCOF(NCOL,NROW,NLAY), IBOUND(NCOL,NROW,NLAY),
     &          HNEW(NCOL,NROW,NLAY), PXDP(NCOL,NROW,NSEGAR),
     &          PETM(NCOL,NROW,NSEGAR)
C     ------------------------------------------------------------------
C
C1------PROCESS EACH HORIZONTAL CELL LOCATION
      DO 30 IR=1,NROW
        DO 20 IC=1,NCOL
C
          IF (NETSOP.EQ.1) THEN
C2----------SET THE LAYER INDEX EQUAL TO 1      .
            IL=1
          ELSEIF (NETSOP.EQ.2) THEN
C3----------GET LAYER INDEX FROM IETS ARRAY
            IL=IETS(IC,IR)
            IF (IL.EQ.0) GO TO 20  ! ERB 1/11/07
          ELSEIF (NETSOP.EQ.3) THEN
C3A---------FIND UPPERMOST ACTIVE CELL
            IL = 1 ! If stack is inactive, this is the default
            FINDFIRST: DO ILQ=1,NLAY
              IF (IBOUND(IC,IR,ILQ).NE.0) THEN
                IL = ILQ
                EXIT FINDFIRST
              ENDIF
            ENDDO FINDFIRST
!          ELSEIF (NETSOP.EQ.4) THEN
!C3B---------FIND UPPERMOST ACTIVE CELL IN CONTIGUOUS STACK OF ACTIVE CELLS
!C           THAT INCLUDES LOWERMOST ACTIVE CELL
!            IL = 1 ! If stack is inactive, this is the default
!            FINDLAY: DO KL=NLAY,1,-1
!              IF (IBOUND(IC,IR,KL)==0) THEN
!                ! Cell is inactive, go up a layer
!                CYCLE FINDLAY
!              ELSE
!                ! Cell is active (variable-head or constant-head)
!                IF (KL>1) THEN
!                  IF (IBOUND(IC,IR,KL-1) .NE. 0) THEN
!                    ! Cell above is active; go up a layer
!                    CYCLE FINDLAY
!                  ENDIF
!                ENDIF
!              ENDIF
!              !   Cell is active and cell above is either (a) inactive, or
!              !   (b) nonexistent because this cell is in layer 1:
!              !   Evapotranspiration is from this cell.
!              IL = KL
!              !   Check IL and higher cells for constant head
!              DO KLC=1,IL
!C5B-------------IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION.
!                IF(IBOUND(IC,IR,KLC).LT.0) GOTO 20
!              ENDDO
!              EXIT FINDLAY
!            ENDDO FINDLAY
          ENDIF
C
C4------IF THE CELL IS EXTERNAL IGNORE IT.
          IF (IBOUND(IC,IR,IL).GT.0) THEN
            C=ETSR(IC,IR)
            S=ETSS(IC,IR)
            SS=S
            HH=HNEW(IC,IR,IL)
C
C5------IF HEAD IN CELL IS GREATER THAN OR EQUAL TO ETSS, ET IS CONSTANT
            IF(HH.GE.SS) THEN
C
C5A-----SUBTRACT -ETSR FROM RHS
              RHS(IC,IR,IL)=RHS(IC,IR,IL) + C
            ELSE
C
C6------IF DEPTH TO WATER>=EXTINCTION DEPTH THEN ET IS 0
              DD=SS-HH
              X=ETSX(IC,IR)
              XX=X
              IF (DD.LT.XX) THEN
C7------VARIABLE RANGE. ADD ET TERMS TO BOTH RHS AND HCOF.
C
                IF (NETSEG.GT.1) THEN
C                 DETERMINE WHICH SEGMENT APPLIES BASED ON HEAD, AND
C                 CALCULATE TERMS TO ADD TO RHS AND HCOF
C
C                 SET PROPORTIONS CORRESPONDING TO ETSS ELEVATION
                  PXDP1 = 0.0
                  PETM1 = 1.0
                  DO 10 ISEG = 1,NETSEG
C                   SET PROPORTIONS CORRESPONDING TO LOWER END OF
C                   SEGMENT
                    IF (ISEG.LT.NETSEG) THEN
                      PXDP2 = PXDP(IC,IR,ISEG)
                      PETM2 = PETM(IC,IR,ISEG)
                    ELSE
                      PXDP2 = 1.0
                      PETM2 = 0.0
                    ENDIF
                    IF (DD.LE.PXDP2*XX) THEN
C                     HEAD IS IN DOMAIN OF THIS SEGMENT
                      GOTO 15
                    ENDIF
C                   PROPORTIONS AT LOWER END OF SEGMENT WILL BE FOR
C                   UPPER END OF SEGMENT NEXT TIME THROUGH LOOP
                    PXDP1 = PXDP2
                    PETM1 = PETM2
   10             CONTINUE
   15             CONTINUE
C                 CALCULATE TERMS TO ADD TO RHS AND HCOF BASED ON
C                 SEGMENT THAT APPLIES AT HEAD ELEVATION
                  THCOF = -(PETM1-PETM2)*C/((PXDP2-PXDP1)*X)
                  TRHS = THCOF*(S-PXDP1*X) + PETM1*C
                ELSE
C                 CALCULATE TERMS TO ADD TO RHS AND HCOF BASED ON SIMPLE
C                 LINEAR RELATION OF ET VS. HEAD
                  TRHS = C-C*S/X
                  THCOF = -C/X
                ENDIF
                RHS(IC,IR,IL)=RHS(IC,IR,IL)+TRHS
                HCOF(IC,IR,IL)=HCOF(IC,IR,IL)+THCOF
              ENDIF
            ENDIF
          ENDIF
   20   CONTINUE
   30 CONTINUE
C
C8------RETURN
      RETURN
      END
C=======================================================================
      SUBROUTINE GWF1ETS1BD(NETSOP,IETS,ETSR,ETSX,ETSS,IBOUND,HNEW,NCOL,
     &                      NROW,NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,
     &                      IETSCB,ICBCFL,BUFF,IOUT,PERTIM,TOTIM,NETSEG,
     &                      PXDP,PETM,NSEGAR)
C-----VERSION 201003215 ERB
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR EVAPOTRANSPIRATION SEGMENTS
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*16 VBNM(MSUM), TEXT
      DOUBLE PRECISION HNEW, RATOUT, QQ, HH, SS, DD, XX, HHCOF, RRHS,
     &                 PXDP1, PXDP2
      DIMENSION IETS(NCOL,NROW), ETSR(NCOL,NROW), ETSX(NCOL,NROW),
     &          ETSS(NCOL,NROW), IBOUND(NCOL,NROW,NLAY),
     &          VBVL(4,MSUM), HNEW(NCOL,NROW,NLAY),
     &          BUFF(NCOL,NROW,NLAY), PXDP(NCOL,NROW,NSEGAR),
     &          PETM(NCOL,NROW,NSEGAR)
C
      DATA TEXT /'     ET SEGMENTS'/
C     ------------------------------------------------------------------
C
C1------CLEAR THE RATE ACCUMULATOR.
      ZERO=0.
      RATOUT=ZERO
C
C2------SET CELL-BY-CELL BUDGET SAVE FLAG (IBD) AND CLEAR THE BUFFER.
      IBD=0
      IF(IETSCB.GT.0) IBD=ICBCFL
      DO 30 IL=1,NLAY
        DO 20 IR=1,NROW
          DO 10 IC=1,NCOL
            BUFF(IC,IR,IL)=ZERO
   10     CONTINUE
   20   CONTINUE
   30 CONTINUE
C
C3------PROCESS EACH HORIZONTAL CELL LOCATION.
      DO 70 IR=1,NROW
        DO 60 IC=1,NCOL
C
          IF (NETSOP.EQ.1) THEN
C2----------SET THE LAYER INDEX EQUAL TO 1      .
            IL=1
          ELSEIF (NETSOP.EQ.2) THEN
C3----------GET LAYER INDEX FROM IETS ARRAY
            IL=IETS(IC,IR)
            IF (IL.EQ.0) GO TO 60  ! ERB 1/11/07
          ELSEIF (NETSOP.EQ.3) THEN
C3A---------FIND UPPERMOST ACTIVE CELL
            IL = 1 ! If stack is inactive, this is the default
            FINDFIRST: DO ILQ=1,NLAY
              IF (IBOUND(IC,IR,ILQ).NE.0) THEN
                IL = ILQ
                EXIT FINDFIRST
              ENDIF
            ENDDO FINDFIRST
            IETS(IC,IR) = IL
!          ELSEIF (NETSOP.EQ.4) THEN
!C3B---------FIND UPPERMOST ACTIVE CELL IN CONTIGUOUS STACK OF ACTIVE CELLS
!C           THAT INCLUDES LOWERMOST ACTIVE CELL
!            IL = 1 ! If stack is inactive, this is the default
!            FINDLAY: DO KL=NLAY,1,-1
!              IF (IBOUND(IC,IR,KL)==0) THEN
!                ! Cell is inactive, go up a layer
!                CYCLE FINDLAY
!              ELSE
!                ! Cell is active (variable-head or constant-head)
!                IF (KL>1) THEN
!                  IF (IBOUND(IC,IR,KL-1) .NE. 0) THEN
!                    ! Cell above is active; go up a layer
!                    CYCLE FINDLAY
!                  ENDIF
!                ENDIF
!              ENDIF
!              !   Cell is active and cell above is either (a) inactive, or
!              !   (b) nonexistent because this cell is in layer 1:
!              !   Evapotranspiration is from this cell.
!              IL = KL
!              !   Check IL and higher cells for constant head
!              DO KLC=1,IL
!C5B-------------IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION.
!                IF(IBOUND(IC,IR,KLC).LT.0) GOTO 60
!              ENDDO
!              EXIT FINDLAY
!            ENDDO FINDLAY
!            IETS(IC,IR) = IL
          ENDIF
C
C6------IF CELL IS EXTERNAL THEN IGNORE IT.
          IF (IBOUND(IC,IR,IL).GT.0) THEN
            C=ETSR(IC,IR)
            S=ETSS(IC,IR)
            SS=S
            HH=HNEW(IC,IR,IL)
C
C7------IF HEAD IN CELL => ETSS,SET Q=MAX ET RATE.
            IF (HH.GE.SS) THEN
             QQ=-C
            ELSE
C
C8------IF DEPTH=>EXTINCTION DEPTH, ET IS 0.
              X=ETSX(IC,IR)
              XX=X
              DD=SS-HH
              IF (DD.LT.XX) THEN
C9------VARIABLE RANGE.  CALCULATE Q DEPENDING ON NUMBER OF SEGMENTS
C
                IF (NETSEG.GT.1) THEN
C                 DETERMINE WHICH SEGMENT APPLIES BASED ON HEAD, AND
C                 CALCULATE TERMS TO ADD TO RHS AND HCOF
C
C                 SET PROPORTIONS CORRESPONDING TO ETSS ELEVATION
                  PXDP1 = 0.0
                  PETM1 = 1.0
                  DO 40 ISEG = 1,NETSEG
C                   SET PROPORTIONS CORRESPONDING TO LOWER END OF
C                   SEGMENT
                    IF (ISEG.LT.NETSEG) THEN
                      PXDP2 = PXDP(IC,IR,ISEG)
                      PETM2 = PETM(IC,IR,ISEG)
                    ELSE
                      PXDP2 = 1.0
                      PETM2 = 0.0
                    ENDIF
                    IF (DD.LE.PXDP2*XX) THEN
C                     HEAD IS IN DOMAIN OF THIS SEGMENT
                      GOTO 50
                    ENDIF
C                   PROPORTIONS AT LOWER END OF SEGMENT WILL BE FOR
C                   UPPER END OF SEGMENT NEXT TIME THROUGH LOOP
                    PXDP1 = PXDP2
                    PETM1 = PETM2
   40             CONTINUE
   50             CONTINUE
C9------CALCULATE ET RATE BASED ON SEGMENT THAT APPLIES AT HEAD
C9------ELEVATION
                  HHCOF = -(PETM1-PETM2)*C/((PXDP2-PXDP1)*X)
                  RRHS = -HHCOF*(S-PXDP1*X) - PETM1*C
                ELSE
C10-----SIMPLE LINEAR RELATION.  Q=-ETSR*(HNEW-(ETSS-ETSX))/ETSX, WHICH
C10-----IS FORMULATED AS Q= -HNEW*ETSR/ETSX + (ETSR*ETSS/ETSX -ETSR).
                  HHCOF = -C/X
                  RRHS = (C*S/X) - C
                ENDIF
                QQ = HH*HHCOF + RRHS
              ELSE
                QQ = 0.0
              ENDIF
            ENDIF
C
C10-----ACCUMULATE TOTAL FLOW RATE.
            Q=QQ
            RATOUT=RATOUT-QQ
C
C11-----ADD Q TO BUFFER.
            BUFF(IC,IR,IL)=Q
          ENDIF
   60   CONTINUE
   70 CONTINUE
C
C12-----IF CELL-BY-CELL FLOW TO BE SAVED, CALL APPROPRIATE UTILITY
C12-----MODULE SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IETSCB,BUFF,NCOL,NROW,
     &                         NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV3(KSTP,KPER,TEXT,IETSCB,BUFF,IETS,NETSOP,
     &                   NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
C
C13-----MOVE TOTAL ET RATE INTO VBVL FOR PRINTING BY BAS1OT.
      ROUT=RATOUT
      VBVL(3,MSUM)=ZERO
      VBVL(4,MSUM)=ROUT
C
C14-----ADD ET(ET_RATE TIMES STEP LENGTH) TO VBVL.
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
C
C15-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS1OT.
      VBNM(MSUM)=TEXT
C
C16-----INCREMENT BUDGET TERM COUNTER.
      MSUM=MSUM+1
C
C17-----RETURN.
      RETURN
      END
C=======================================================================

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
c  tcspline -- dmc 20 Jan 1999
c
c  set up coefficients for bicubic spline with following BC's:
c  * LHS and RHS handled as in cubspl.for for 1st coordinate
c  * derivatives periodic in second coordinate (use pspline.for)
c
c workspace:
c  if phi bdy cond. is periodic, not-a-knot, df/dphi = 0 everywhere,
c  or d2f/dphi2 = 0 everywhere, then the phi boundary condition is
c  "linear" and a workspace of size at least:
c
c     nwk = 20*inx*inth + 10*max(inx,inth,inph)
c
c  will suffice.
c
c  if the phi bdy cond. involves specification of df/dphi .ne. 0 or
c  d2f/dphi .ne. 0 at any (x,theta) grid point, then, the phi boundary
c  condition is "non-linear", a correction step is needed, and a workspace
c  of size at least:
c
c     nwk = 16*inx*inth*inph
c
c  is required.
c
      subroutine tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
     >                    ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
     >                    ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th,
     >                    ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph,
     >                    wk,nwk,ilinx,ilinth,ilinph,ier)
c
      real x(inx),th(inth),ph(inph)
      real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
      real bcxmin(inb1x,*),bcxmax(inb1x,*) ! inth x inph defined (if used)
      real bcthmin(inb1th,*),bcthmax(inb1th,*) ! inx x inph defined (if used)
      real bcphmin(inb1ph,*),bcphmax(inb1ph,*) ! inx x inth defined (if used)
c
c  input:
c    x(1...inx) -- abscissae, first dimension of data
c   th(1...inth) -- abscissae, second (periodic) dimension of data
c   ph(1...inph) -- abscissae, third (periodic) dimension of data
c   fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
c   inf4 -- fspl dimensioning, inf4.ge.inx required.
c   inf5 -- fspl dimensioning, inf5.ge.inth required.
c
c  boundary conditions input:
c
c   bc data at xmin, xmax  vs.  theta,phi
c   bc data at thmin, thmax  vs.  x,phi
c   bc data at phmin, phmax  vs.  x,theta
c
c   ibcxmin -- indicator for boundary condition at x(1):
c    bcxmin(...) -- boundary condition data
c     =-1 -- use periodic boundary condition
c     =0 -- use "not a knot", bcxmin(...) ignored
c     =1 -- match slope, specified at x(1),th(ith),ph(iph) by bcxmin(ith,iph)
c     =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
c           by bcxmin(ith,iph
c     =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c     =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c     =5 -- match 1st derivative to 1st divided difference
c     =6 -- match 2nd derivative to 2nd divided difference
c     =7 -- match 3rd derivative to 3rd divided difference
c           (for more detailed definition of BCs 5-7, see the
c           comments of subroutine mkspline)
c   NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c   ibcxmax -- indicator for boundary condition at x(nx):
c    bcxmax(...) -- boundary condition data
c     (interpretation as with ibcxmin, bcxmin)
c     NOTE:  if ibcxmin=-1 then the periodic BC applies on both sides
c            and ibcxmax is ignored.
c   inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax .gt. 0
c            this must be .ge. inth:
c
c   interpretation of ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th
c     is same as with ibcxmin,...
c
c   interpretation of ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph
c     is same as with ibcxmin,...
c
c   the explicit bdy condition arrays are referenced only if the
c     corresponding "ibc" flag values are set to 1 or 2.
c
c  output:
c   fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
c   ...fspl(1,1,1,*,*,*) is not replaced.
c
c   ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
c   ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c   ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
c
c   ier -- completion code, 0 for normal
c
c  workspace:
c   wk -- must be at least 5*max(inx,inth,inph) large -- or more, see
c         comments, above.
c  nwk -- size of workspace
c
c---------------------------------
c  ** in what follows, f is an abbreviation for fspl **
c
c  compute tricubic spline of 3d function, given values at the
c  grid crossing points, f(1,1,1,i,j,k)=f(x(i),th(j),ph(k)).
c
c  on evaluation:  for point x btw x(i) and x(i+1), dx=x-x(i)
c                       and th btw th(j) and th(j+1), dt=th-th(j),
c                       and ph btw ph(k) and ph(k+1), dp=ph-ph(k),
c
c      spline =
c        f(1,1,1,i,j,k)+dx*f(2,1,1,i,j,k)+dx2*f(3,1,1,i,j,k)+dx3*f(4,1,1,i,j,k)
c  +dt*(f(1,2,1,i,j,k)+dx*f(2,2,1,i,j,k)+dx2*f(3,2,1,i,j,k)+dx3*f(4,2,1,i,j,k))
c +dt2*(f(1,3,1,i,j,k)+dx*f(2,3,1,i,j,k)+dx2*f(3,3,1,i,j,k)+dx3*f(4,3,1,i,j,k))
c +dt3*(f(1,4,1,i,j,k)+dx*f(2,4,1,i,j,k)+dx2*f(3,4,1,i,j,k)+dx3*f(4,4,1,i,j,k))
c        +dp*(
c        f(1,1,2,i,j,k)+dx*f(2,1,2,i,j,k)+dx2*f(3,1,2,i,j,k)+dx3*f(4,1,2,i,j,k)
c  +dt*(f(1,2,2,i,j,k)+dx*f(2,2,2,i,j,k)+dx2*f(3,2,2,i,j,k)+dx3*f(4,2,2,i,j,k))
c +dt2*(f(1,3,2,i,j,k)+dx*f(2,3,2,i,j,k)+dx2*f(3,3,2,i,j,k)+dx3*f(4,3,2,i,j,k))
c +dt3*(f(1,4,2,i,j,k)+dx*f(2,4,2,i,j,k)+dx2*f(3,4,2,i,j,k)+dx3*f(4,4,2,i,j,k)))
c        +dp2*(
c        f(1,1,3,i,j,k)+dx*f(2,1,3,i,j,k)+dx2*f(3,1,3,i,j,k)+dx3*f(4,1,3,i,j,k)
c  +dt*(f(1,2,3,i,j,k)+dx*f(2,2,3,i,j,k)+dx2*f(3,2,3,i,j,k)+dx3*f(4,2,3,i,j,k))
c +dt2*(f(1,3,3,i,j,k)+dx*f(2,3,3,i,j,k)+dx2*f(3,3,3,i,j,k)+dx3*f(4,3,3,i,j,k))
c +dt3*(f(1,4,3,i,j,k)+dx*f(2,4,3,i,j,k)+dx2*f(3,4,3,i,j,k)+dx3*f(4,4,3,i,j,k)))
c        +dp3*(
c        f(1,1,4,i,j,k)+dx*f(2,1,4,i,j,k)+dx2*f(3,1,4,i,j,k)+dx3*f(4,1,4,i,j,k)
c  +dt*(f(1,2,4,i,j,k)+dx*f(2,2,4,i,j,k)+dx2*f(3,2,4,i,j,k)+dx3*f(4,2,4,i,j,k))
c +dt2*(f(1,3,4,i,j,k)+dx*f(2,3,4,i,j,k)+dx2*f(3,3,4,i,j,k)+dx3*f(4,3,4,i,j,k))
c +dt3*(f(1,4,4,i,j,k)+dx*f(2,4,4,i,j,k)+dx2*f(3,4,4,i,j,k)+dx3*f(4,4,4,i,j,k)))
c
c      where dx2=dx**2 and dx3=dx**3.
c      where dt2=dt**2 and dt3=dt**3.
c      where dp2=dp**2 and dp3=dp**3.
c
c---------------------------------
      integer iselect1(10)
      integer iselect2(10)
c
      real z0,z1,ztol
      real zcur(1)
c
      data z0/0.0e0/
      data z1/1.0e0/
      data ztol/1.0e-3/
c
c---------------------------------
c
      ier=0
c
      iflg=0
c
c  check phi bdy condition "linearity"
c
      if(ibcphmin.ne.-1) then
         if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) then
            do ith=1,inth
               do ix=1,inx
                  if(bcphmin(ix,ith).ne.z0) iflg=1
               enddo
            enddo
         endif
         if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) then
            do ith=1,inth
               do ix=1,inx
                  if(bcphmax(ix,ith).ne.z0) iflg=1
               enddo
            enddo
         endif
      endif
c
      itest=10*max(inx,inth,inph) + 20*inx*inth
      if(iflg.eq.1) then
         itest=16*inx*inth*inph
      endif
C
      if(nwk.lt.itest) then
         write(6,9901) nwk,itest
         ier=1
 9901    format(' ?tcspline:  workspace too small.'/
     >      '  user supplied nwk=',i7,'; need at least: ',i7/
     >      '  If no explicit df/dph boundary condition is set,'/
     >      '  nwk = 20*inx*inth + 10*max(inx,inth,inph) can be used.'/
     >      '  If an explicit df/dph or d2f/dph2 boundary condition'/
     >      '  is set, nwk=16*inx*inth*inph is required.')
      endif
      if(inx.lt.2) then
         write(6,'('' ?tcspline:  at least 2 x points required.'')')
         ier=1
      endif
      if(inth.lt.2) then
         write(6,'('' ?tcspline:  need at least 2 theta points.'')')
         ier=1
      endif
      if(inph.lt.2) then
         write(6,'('' ?tcspline:  need at least 2 phi points.'')')
         ier=1
      endif
c
      if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
     >   (ibcxmax.eq.2)) then
         if(inb1x.lt.inth) then
            ier=1
            write(6,
     >'('' ?tcspline:  1st dim of bcxmin/max arrays .lt. inth'')')
         endif
      endif
c
      if((ibcthmin.eq.1).or.(ibcthmax.eq.1).or.(ibcthmin.eq.2).or.
     >   (ibcthmax.eq.2)) then
         if(inb1th.lt.inx) then
            ier=1
            write(6,
     >'('' ?tcspline:  1st dim of bcthmin/max arrays .lt. inx'')')
         endif
      endif
c
      if((ibcphmin.eq.1).or.(ibcphmax.eq.1).or.(ibcphmin.eq.2).or.
     >   (ibcphmax.eq.2)) then
         if(inb1ph.lt.inx) then
            ier=1
            write(6,
     >'('' ?tcspline:  1st dim of bphmin/max arrays .lt. inx'')')
         endif
      endif
c
      call ibc_ck(ibcxmin,'tcspline','xmin',-1,7,ier)
      if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'tcspline','xmax',0,7,ier)
c
      call ibc_ck(ibcthmin,'tcspline','thmin',-1,7,ier)
      if(ibcthmin.ge.0) call ibc_ck(ibcthmax,'tcspline','thmax',0,7,ier)
c
      call ibc_ck(ibcphmin,'tcspline','phmin',-1,7,ier)
      if(ibcphmax.ge.0) call ibc_ck(ibcphmax,'tcspline','phmax',0,7,ier)
c
c  check ilinx & x vector
c
      call splinck(x,inx,ilinx,ztol,ierx)
      if(ierx.ne.0) ier=2
c
      if(ier.eq.2) then
         write(6,'('' ?tcspline:  x axis not strict ascending'')')
      endif
c
c  check ilinth & th vector
c
      call splinck(th,inth,ilinth,ztol,ierth)
      if(ierth.ne.0) ier=3
c
      if(ier.eq.3) then
         write(6,'('' ?tcspline:  theta axis not strict ascending'')')
      endif
c
c  check ilinth & th vector
c
      call splinck(ph,inph,ilinph,ztol,ierph)
      if(ierph.ne.0) ier=4
c
      if(ier.eq.4) then
         write(6,'('' ?tcspline:  phi axis not strict ascending'')')
      endif
c
      if(ier.ne.0) return
c
c------------------------------------
c
c  part 1.  compute (x,theta) spline coeffs via an intermediate
c  routine that call bcspline
c
c  workspace addresses
c
      iaspl2=1
      iabcx1=iaspl2+16*inx*inth
      iabcx2=iabcx1+inth
      iabcth1=iabcx2+inth
      iabcth2=iabcth1+inx
      iawk=iabcth2+inx
      inwk=nwk-iawk+1
c
      do iph=1,inph
c
c  copy bc data
c
         do ix=1,inx
            wk(iabcth1+ix-1)=0.0
            wk(iabcth2+ix-1)=0.0
            if((ibcthmin.eq.1).or.(ibcthmin.eq.2)) then
               wk(iabcth1+ix-1)=bcthmin(ix,iph)
            endif
            if((ibcthmin.ne.-1).and.
     >         ((ibcthmax.eq.1).or.(ibcthmax.eq.2))) then
               wk(iabcth2+ix-1)=bcthmax(ix,iph)
            endif
         enddo
         do ith=1,inth
            wk(iabcx1+ith-1)=0.0
            wk(iabcx2+ith-1)=0.0
            if((ibcxmin.eq.1).or.(ibcxmin.eq.2)) then
               wk(iabcx1+ith-1)=bcxmin(ith,iph)
            endif
            if((ibcxmin.ne.-1).and.
     >         ((ibcxmax.eq.1).or.(ibcxmax.eq.2))) then
               wk(iabcx2+ith-1)=bcxmax(ith,iph)
            endif
         enddo
c
c  call 2d spline intermediary routine
c
         call tcsp23(x,inx,th,inth,fspl(1,1,1,1,1,iph),inf4,
     >      ibcxmin,wk(iabcx1),ibcxmax,wk(iabcx2),
     >      ibcthmin,wk(iabcth1),ibcthmax,wk(iabcth2),
     >      wk(iaspl2),wk(iawk),inwk,ilinx,ilinth,ier)
c
         if(ier.ne.0) then
            write(6,*) ' ?tcspline:  error in 2d spline, exiting.'
            return
         endif
c
      enddo
c
c  ok now fspl(*,*,1,*,*,*) have been evaluated and C2 in (x,theta)
c  now need to extend to coeffs in phi direction.
c
      xo2=0.5e0
      xo6=1.0e0/6.0e0
c
c  spline each (x,th) coeff in the phi direction
c
      inpho=4*(inph-1)
      do ith=1,inth-1
         do ix=1,inx-1
c
            do ic1=1,4
               do ic2=1,4
c
c  copy coeff. ordinates in
c
                  do iph=1,inph
                     wk(4*(iph-1)+1)=fspl(ic1,ic2,1,ix,ith,iph)
                  enddo
c
c  use linear BC on this first pass; will correct later if
c  necessary
c
                  wk(2)=0.0
                  wk(3)=0.0
                  wk(inpho+2)=0.0
                  wk(inpho+3)=0.0
c
                  ibcphmina=ibcphmin
                  ibcphmaxa=ibcphmax
                  if(iflg.eq.1) then
                     if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) ibcphmina=0
                     if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) ibcphmaxa=0
                  endif
c
                  call v_spline(ibcphmina,ibcphmaxa,inph,ph,wk,
     >               wk(4*inph+1))
c
c  copy coeffs out
c
                  do iph=1,inph-1
                     fspl(ic1,ic2,2,ix,ith,iph)=wk(4*(iph-1)+2)
                     fspl(ic1,ic2,3,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
                     fspl(ic1,ic2,4,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
                  enddo
c
               enddo                    ! ic2
            enddo                       ! ic1
c
         enddo                          ! ix
      enddo                             ! ith
c
c  if there are "non-linear" BCs requiring correction...
c
c  at each (x(ix),th(ith)) get the d/dph BC's right while preserving C2
c  everywhere...
c
      if(iflg.eq.1) then
c
c  first get BC correction numbers
c
         iabcph1=1
         iabcph2=iabcph1+inx*inth
         iaccoef=iabcph2+inx*inth
         iawk=iaccoef+12*inx*inth*inph
         inwk=nwk-iawk+1
c
         do i=1,10
            iselect1(i)=0
            iselect2(i)=0
         enddo
c
c  note because iflg=1, we know at least one of ibcphmin/max = 1 or 2
c
         iskip1=0
         if(ibcphmin.eq.1) then
            iselect1(4)=1               ! df/dph
         else if(ibcphmin.eq.2) then
            iselect1(7)=1               ! d2f/dph2
         else
            iskip1=1
         endif
c
         iskip2=0
         if(ibcphmax.eq.1) then
            iselect2(4)=1               ! df/dph
         else if(ibcphmax.eq.2) then
            iselect2(7)=1               ! d2f/dph2
         else
            iskip2=1
         endif
c
         ia1=iabcph1-1
         ia2=iabcph2-1
         do ith=1,inth
            do ix=1,inx
               ia1=ia1+1
               ia2=ia2+1
c
               if(iskip1.eq.0) then
                  call tcspeval(x(ix),th(ith),ph(1),iselect1, zcur,
     >               x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
     >               fspl,inf4,inf5,ier)
                  if(ier.ne.0) then
                     write(6,*) ' ?? tcspline:  error in tcspeval call'
                     return
                  endif
                  wk(ia1)=bcphmin(ix,ith)-zcur(1) ! correction needed
               else
                  wk(ia1)=z0
               endif
c
               if(iskip2.eq.0) then
                  call tcspeval(x(ix),th(ith),ph(inph),iselect2, zcur,
     >               x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
     >               fspl,inf4,inf5,ier)
                  if(ier.ne.0) then
                     write(6,*) ' ?? tcspline:  error in tcspeval call'
                     return
                  endif
                  wk(ia2)=bcphmax(ix,ith)-zcur(1) ! correction needed
               else
                  wk(ia2)=z0
               endif
            enddo
         enddo
c
         call tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
     >      ibcxmin,ibcxmax,ibcthmin,ibcthmax,
     >      ibcphmin,wk(iabcph1),ibcphmax,wk(iabcph2),
     >      wk(iaccoef),wk(iawk),inwk)
c
      endif
c
      return
      end
c-----------------------------------------
      subroutine tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
     >   ibcxmin,ibcxmax,ibcthmin,ibcthmax,
     >   ibcphmin,bcph1,ibcphmax,bcph2,ccorr,wk,nwk)
c
c  intermediary routine for tcspline:
c  do correction needed to get C2 3d spline with phi bdy conditions
c  matched.
c
c  all input unless noted:
c
      real x(inx)                       ! x axis
      real th(inth)                     ! th axis
      real ph(inph)                     ! ph axis
c
      real fspl(4,4,4,inf4,inf5,inph)   ! spline coeffs -- adjusted
c
      integer ibcxmin,ibcxmax           ! x BC flags
      integer ibcthmin,ibcthmax         ! th BC flags
      integer ibcphmin,ibcphmax         ! ph BC flags
c
      real bcph1(inx,inth)              ! ph BC correction array @ phi(1)
      real bcph2(inx,inth)              ! ph BC correction array @ phi(inph)
c
c  workspaces:
c
      real ccorr(3,4,inx,inth,inph)     ! correction coefficients (partial)
c
      real wk(nwk)                      ! workspace
c
c---------------------------
c
      xo2=0.5e0
      xo6=1.0e0/6.0e0
c
c  1.  splines in phi -- fcns are zero everywhere but have non-zero BCs
c
      if(nwk.lt.10*max(inx,inth,inph)) then
         write(6,*) ' ?? programming error in tcspcorr (tcspline)'
         return
      endif
c
      z0=0.0
      iawk2=4*inph+1
c
      inpho=4*(inph-1)
      do ith=1,inth
         do ix=1,inx
c
            do iph=1,inph
               wk(4*(iph-1)+1)=z0
            enddo
c
c  set BC for this 1d spline correction
c
            if(ibcphmin.eq.1) then
               wk(2)=bcph1(ix,ith)
            else if(ibcphmin.eq.2) then
               wk(3)=bcph1(ix,ith)
            endif
c
            if(ibcphmax.eq.1) then
               wk(inpho+2)=bcph2(ix,ith)
            else if(ibcphmax.eq.2) then
               wk(inpho+3)=bcph2(ix,ith)
            endif
c
            call v_spline(ibcphmin,ibcphmax,inph,ph,wk,wk(iawk2))
c
c  copy non-zero coeffs out to ccorr
c
            do iph=1,inph-1
               ccorr(1,1,ix,ith,iph)=wk(4*(iph-1)+2)
               ccorr(2,1,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
               ccorr(3,1,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
            enddo
c
         enddo
      enddo
c
c  2. spline the coeffs in x -- use ibcx flags & zero for derivative
c  bc if necessary
c
      iawk2=4*inx+1
c
      inxo=4*(inx-1)
      do iph=1,inph-1
         do ith=1,inth
c
            do icph=1,3
c
               do ix=1,inx
                  wk(4*(ix-1)+1)=ccorr(icph,1,ix,ith,iph)
               enddo
c
c  zero BC:  correction spline
c
               wk(2)=0.0
               wk(3)=0.0
               wk(inxo+2)=0.0
               wk(inxo+3)=0.0
c
               call v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(iawk2))
c
               do ix=1,inx-1
                  ccorr(icph,2,ix,ith,iph)=wk(4*(ix-1)+2)
                  ccorr(icph,3,ix,ith,iph)=wk(4*(ix-1)+3)*xo2
                  ccorr(icph,4,ix,ith,iph)=wk(4*(ix-1)+4)*xo6
               enddo
c
            enddo
c
         enddo
      enddo
c
c  3.  spline all the ccorr coefs in th -- use ibcth flags & zero for
c      derivative correction BC if necessary
c
c      add the results into fspl
c
      iawk2=4*inth+1
c
      intho=4*(inth-1)
      do iph=1,inph-1
         do ix=1,inx-1
c
            do icx=1,4
               do icph=1,3
c
                  do ith=1,inth
                     wk(4*(ith-1)+1)=ccorr(icph,icx,ix,ith,iph)
                  enddo
c
c  zero BC:  correction spline
c
                  wk(2)=0.0
                  wk(3)=0.0
                  wk(intho+2)=0.0
                  wk(intho+3)=0.0
c
                  call v_spline(ibcthmin,ibcthmax,inth,th,wk,
     >               wk(iawk2))
c
                  do ith=1,inth-1
                     do icth=1,4
                        zfac=1.0
                        if(icth.eq.3) zfac=xo2
                        if(icth.eq.4) zfac=xo6
                        fspl(icx,icth,icph+1,ix,ith,iph)=
     >                     fspl(icx,icth,icph+1,ix,ith,iph)+
     >                     wk(4*(ith-1)+icth)*zfac
                     enddo
                  enddo
c
               enddo
            enddo
c
         enddo
      enddo
c
      return
      end
c-----------------------------------------
      subroutine tcsp23(x,inx,th,inth,fspl,inf4,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,
     >   ibcthmin,bcthmin,ibcthmax,bcthmax,
     >   fspl2,wk,nwk,ilinx,ilinth,ier)
c
c  intermediary routines
c  call bcspline from tcspline loop
c  to set up 2d splines in each phi plane
c
      real x(inx)                       ! x axis
      real th(inth)                     ! th axis
c
      real fspl(4,4,4,inf4,inth)        ! fspl array at one phi pt.
      real fspl2(4,4,inx,inth)          ! temp fspl array for bcspline
c
      real bcxmin(inth),bcxmax(inth)    ! d/dx BC's @ x(1),x(inx), th(*)
      real bcthmin(inx),bcthmax(inx)    ! d/dth BC's @ th(1),th(inth), x(*)
c
      real wk(nwk)
c
c--------------------
c
c  1.  copy spline data in
c
      do ith=1,inth
         do ix=1,inx
            fspl2(1,1,ix,ith)=fspl(1,1,1,ix,ith)
         enddo
      enddo
c
c  2.  compute the 2d spline
c
      call bcspline(x,inx,th,inth,fspl2,inx,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,
     >   ibcthmin,bcthmin,ibcthmax,bcthmax,
     >   wk,nwk,ilinx,ilinth,ier)
      if(ier.ne.0) return
c
c  3.  copy spline coeff results out
c
      do ith=1,inth-1
         do ix=1,inx-1
            do j=1,4
               do i=1,4
                  fspl(i,j,1,ix,ith)=fspl2(i,j,ix,ith)
               enddo
            enddo
         enddo
      enddo
c
      return
      end

max_stars_repo_name: ORNL-Fusion/RFSciDAC-testing

max_issues_repo_name: ORNL-Fusion/RFSciDAC-testing

max_forks_repo_name: ORNL-Fusion/RFSciDAC-testing


text: The Department of French and Italian, part of the College of Letters and Science, offers both a majors major and minors minor program in Italian. Courses include Italian language, literature, and culture.


text:
      subroutine urqmd(n)
      include 'epos.inc'

      print *,n
      if(iurqmd.eq.1)stop'compile with eposu.f instead of eposu_no.f '

      end

max_stars_repo_name: ckamtsikis/cmssw

max_issues_repo_name: ckamtsikis/cmssw

max_forks_repo_name: ckamtsikis/cmssw


! { dg-do run }
! { dg-require-effective-target lto }
! { dg-options "-flto" }
!
! Checks that the results of module procedures have the correct characteristics
! and that submodules use the module version of vtables (PR66762). This latter
! requires the -flto compile option.
!
! Contributed by Reinhold Bader  <reinhold.bader@lrz.de>
!
module mod_a
  implicit none
  type, abstract :: t_a
  end type t_a
  interface
    module subroutine p_a(this, q)
      class(t_a), intent(inout) :: this
      class(*), intent(in) :: q
    end subroutine
    module function create_a() result(r)
      class(t_a), allocatable :: r
    end function
    module subroutine print(this)
      class(t_a), intent(in) :: this
    end subroutine
  end interface
end module mod_a

module mod_b
  implicit none
  type t_b
    integer, allocatable :: I(:)
  end type t_b
  interface
    module function create_b(i) result(r)
      type(t_b) :: r
      integer :: i(:)
    end function
  end interface
end module mod_b

submodule(mod_b) imp_create
contains
  module procedure create_b
    if (allocated(r%i)) deallocate(r%i)
    allocate(r%i, source=i)
  end procedure
end submodule imp_create

submodule(mod_a) imp_p_a
  use mod_b
  type, extends(t_a) :: t_imp
    type(t_b) :: b
  end type t_imp
  integer, parameter :: ii(2) = [1,2]
contains
  module procedure create_a
    type(t_b) :: b
    b = create_b(ii)
    allocate(r, source=t_imp(b))
  end procedure

  module procedure  p_a
    select type (this)
      type is (t_imp)
        select type (q)
          type is (t_b)
            this%b = q
          class default
            call abort
         end select
      class default
        call abort
      end select
  end procedure p_a
  module procedure print
    select type (this)
      type is (t_imp)
        if (any (this%b%i .ne. [3,4,5])) call abort
      class default
        call abort
    end select
  end procedure
end submodule imp_p_a

program p
  use mod_a
  use mod_b
  implicit none
  class(t_a), allocatable :: a
  allocate(a, source=create_a())
  call p_a(a, create_b([3,4,5]))
  call print(a)
end program p
! { dg-final { cleanup-submodules "mod_a@imp_p_a" } }
! { dg-final { cleanup-submodules "mod_b@imp_create" } }


!> A module for for writing help, metadata and the input file
!! to the netcdf file.
!!
!! To extract the input file use the bash utility scripts/extract_input_file
module diagnostics_metadata
  use diagnostics_config, only: diagnostics_type
  implicit none
  character, dimension(100000) :: inputfile_array
  integer :: inputfile_length

contains
  subroutine write_metadata(gnostics)
    use simpledataio, only: add_metadata
    use simpledataio, only: add_standard_metadata
    use runtime_tests, only: build_identifier
    use runtime_tests, only: get_svn_rev
    use runtime_tests, only: is_release, release
    use run_parameters, only: user_comments
    type(diagnostics_type), intent(in) :: gnostics
    character (20) :: datestamp, timestamp, timezone
    character (31) :: strtime
    !integer :: inttime
    integer, dimension(8) :: values

    call add_metadata(gnostics%sfile, "title", &
      jline("")//&
      jline("=================================================")//&
      jline("                GS2 Output File")//&
      jline("=================================================")//&
      jline(""))
    call add_metadata(gnostics%sfile, "file_description", &
      jline("")//&
      jline(" This file contains output from the gyrokinetic")//&
      jline(" flux tube code GS2. GS2 is a code")//&
      jline(" which evolves the gyrokinetic equation, ")//&
      jline(" which describes the behaviour of a strongly")//&
      jline(" magnetized plasma. ")// &
      jline(" Among other information this file may contain: ")//&
      jline(" (*) Values of geometric coefficients. ")//&
      jline(" (*) Values of the electromagnetic fields. ")//&
      jline(" (*) Values of linear growth rates. ")//&
      jline(" (*) Values of turbulent fluxes. ")//&
      jline(" (*) Values of moments such as density, ")//&
      jline("         parallel flow and temperature. ")//&
      jline(" The details of what is contained are controlled")//&
      jline(" by the namelist diagnostics_config. ")//&
      jline(""))
    call add_metadata(gnostics%sfile, "data_description", &
      jline("")//&
      jline(" Most data consists of physical quantities given")//&
      jline(" as a function of one more of the five dimensions")//&
      jline(" in the GK eqn: theta, kx, ky, energy and lamda. ")// &
      jline(" Data is typically double precision real numbers.  ")//&
      jline(" Complex numbers are handled by having an extra")//&
      jline(" dimension ri. ")//&
      jline(""))
    call add_metadata(gnostics%sfile, "normalization_description", &
      jline("")//&
      jline(" Quantities in this file are given in ")//&
      jline(" dimensionless form.  Dimensional quantities can")//&
      jline(" be reconstructed from these dimensionless ")//&
      jline(" quantities by multiplying by the expression ")// &
      jline(" given in the units attribute. This expression")//&
      jline(" is constructed from one or more of the following")//&
      jline(" normalising quantities: ")//&
      jline(" (*) a_ref: the normalising length ")//&
      jline("       (this is half-diameter of ")//&
      jline("       the LCFS for numerical equilibria ")//&
      jline("       but has no physical meaning  ")//&
      jline("       in Miller, circular and slab geometries) ")//&
      jline(" (*) B_ref: the normalising field ")//&
      jline("       (this is the field on the magnetic ")//&
      jline("       axis for numerical equilibria ")//&
      jline("       but has no physical meaning in ")//&
      jline("       Miller, circular and slab geometries) ")//&
      jline(" (*) Properties of the reference species (a ")//&
      jline("       hypothetical species whose properties  ")//&
      jline("       may be equal to one or none of the species ")//&
      jline("       in the simulation): ")//&
      jline("    (*) T_ref: the reference temperature ")//&
      jline("    (*) n_ref: the reference density ")//&
      jline("    (*) n_ref: the reference mass ")//&
      jline("    (*) Z_ref: the reference charge, which ")//&
      jline("          is always equal to the proton charge")//&
      jline(" (*) vth_ref = sqrt(2 T_ref/m_ref). ")//&
      jline(" (*) rho_ref = vth_ref / (Z_ref B_ref / m_ref c)")//&
      jline(""))

    call add_metadata(gnostics%sfile, "gs2_help", &
      jline("")//&
      jline(" At the time of writing you can obtain ")//&
      jline(" help for using GS2 in the following")//&
      jline(" places: ")//&
      jline(" (*) http://gyrokinetics.sourceforge.net/wiki")//&
      jline("       (User help: installing, running) ")//&
            new_line("")//&
            " (*) http://gyrokinetics.sourceforge.net/gs2_documentation/ "//&
      jline("       (Doxygen documentation) ")//&
      jline(" (*) http://sourceforge.net/projects/gyrokinetics ")//&
      jline("       (Downloads, bug tracker) ")//&
      jline(""))
    call add_metadata(gnostics%sfile, "input_file_extraction", &
      jline("")//&
      jline(" A bash utility for extracting the input file ")//&
      jline(" from this file is included in the scripts ")//&
      jline(" folder in the GS2 source. To invoke it:")//&
      jline(" $ ./extract_input_file <netcdf_file>")//&
      jline(""))
    call add_metadata(gnostics%sfile, "user_comments", &
      trim(user_comments))
    call add_metadata(gnostics%sfile, "svn_revision", &
      trim(get_svn_rev()))
    call add_metadata(gnostics%sfile, "build_identifier", &
      trim(build_identifier()))
    if (is_release()) then
      call add_metadata(gnostics%sfile, "release_info", &
        "This is official release number "//trim(build_identifier()))
    else
      call add_metadata(gnostics%sfile, "release_info", &
        jline("This is not an official release (presumably")//&
        jline("a development version)"))
    end if

    call add_standard_metadata(gnostics%sfile)
  end subroutine write_metadata

  function jline(inputline)
    use mp, only: mp_abort
    character(*), intent(in) :: inputline
    character(len=52) :: jline
    integer :: lenin
    lenin = len(trim(inputline))
    if (lenin .gt. 50) then
      call mp_abort("Long line: "//inputline, .true.)
    end if

    jline = ''

    jline(2:lenin+1) = inputline(1:lenin)
    jline(1:1) = new_line("")

  end function jline

  subroutine read_input_file_and_get_size(gnostics)
    use file_utils, only: get_input_unit
    use mp, only: proc0, broadcast
    type(diagnostics_type), intent(in) :: gnostics
    character(len=1000) :: line
    integer :: inputunit
    integer :: ios, i
    !Note : The input file used here is not directly the input file passed to GS2
    !as during init_input_unit (file_utils) the passed input file is processed to
    !a) Strip out any comments
    !b) Deal with any included files.
    !It is this processed file (written to .<run_name>.in) which is stored.
    !We may wish to add the option to retain comments as these may be useful
    !when revisting old files to remind the user why the file is setup as it is.

    ! Find out the size of the input file and read it in
    ! This is unnecessarily done twice, as this function
    ! is called once to create the variable and once to
    ! write, but the overhead is trivial
    if (proc0) then
      call get_input_unit (inputunit)
      rewind (unit=inputunit)
      ios = 0
      inputfile_length = 0
      do while (ios .eq. 0)
        line = ''
        read(unit=inputunit, iostat=ios, fmt=('(A,$)')) &
          line

        !if (len(trim(line)) .gt. 1) then
          do i = 1, len(trim(line))
            inputfile_length = inputfile_length + 1
            inputfile_array(inputfile_length) = line(i:i)
          end do
          inputfile_length = inputfile_length + 1
          inputfile_array(inputfile_length) = " "
          inputfile_length = inputfile_length + 1
          inputfile_array(inputfile_length) = "\" !" !These comment characters are for editors which think \ is escaping "
          inputfile_length = inputfile_length + 1
          inputfile_array(inputfile_length) = "n"
        !end if
      end do
    end if
    call broadcast(inputfile_length)
    call broadcast(inputfile_array)
  end subroutine read_input_file_and_get_size

  !> Write the input file.

  subroutine write_input_file(gnostics)
    use file_utils, only: get_input_unit
    use diagnostics_create_and_write, only: create_and_write_variable
    use simpledataio, only: SDATIO_CHAR
    use diagnostics_dimensions, only: dim_string
    type(diagnostics_type), intent(in) :: gnostics

    call create_and_write_variable(gnostics, SDATIO_CHAR, "input_file", &
      trim(dim_string(gnostics%dims%input_file_dim)), &
       "Full text of the input file", "text", inputfile_array)

  end subroutine write_input_file
end module diagnostics_metadata


text:
C***********************************************************************
C*
C*	Copyright 1996, University Corporation for Atmospheric Research.
C*
C*	wprofbor.f
C*
C*	Used by wind profiler function to draw background. Derived from
C*	the GEMPAK program SNSBOR.
C*
C*	History:
C*
C*	11/96	COMET	Original copy
C*
C************************************************************************

	SUBROUTINE wprofbor ( border, ybot, ytop, yaxval, nyval, iylbfr,
     +			      iyglfr, iytmfr, xmin, xmax, nstn, stns,
     +			      sloc, x, pontha, toptmp, topwnd, xtlbl,
     +			      ctlbl, nxlbl, ixlbfr, ixglfr, ixtmfr,
     +			      timflg, iret )
C************************************************************************
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	border, stns (*), ctlbl (*)
	REAL		sloc(*), x(*), pontha(40,*), toptmp(*),
     +			topwnd(*), yaxval (*), xtlbl (*)
	LOGICAL		timflg, scflag
C*
	REAL		xx (3), yy (3)
	REAL		smooth, filter
csd	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
	values = 0.
C
C*	Get input values; the "surface" line is set to the background
C*	color.
C
	CALL IN_LINE  ( border, values, 1, icbor, ibtyp, ibwid, iblab,
     +                  smooth, filter, scflag, ier )
C
C*	Draw background.
C
	IF  ( icbor .ne. 0 ) THEN
C
C*	    Set background color and line attributes.
C
	    CALL GSCOLR  ( icbor, ier )
	    CALL GQLINE  ( i1, i2, i3, i4, ier )
	    CALL GSLINE  ( ibtyp, 0, ibwid, 0, ier )
C
C*	    Draw y axis.
C
	    CALL GDAXIS ( 2, xmin, .false., iylbfr, iytmfr, iyglfr, -1,
     +			   nyval, yaxval, ier )
C
C*	    Draw x axis.
C
	    IF  ( timflg )  THEN
		CALL GAAXIS  ( 1, ybot, .true., ixlbfr, ixtmfr, ixglfr,
     +			       nxlbl, xtlbl, ctlbl, ier )
	      ELSE
		CALL GAAXIS ( 1, ybot, .true., 101, 0, 0, nstn, sloc,
     +	 		      stns, ier )
	    END IF
C
C*	    Draw station line.
C
       if( .not. timflg ) then
	    DO  i = 1, nstn
		xx (1) = sloc (i)
		xx (2) = sloc (i)
		yy (1) = ybot
csd		IF  ( .not. ERMISS ( topwnd (i) ) )  THEN
		    yy (2) = topwnd (i)
		    CALL GLINE  ( 'M', 2, xx, yy, ier )
		    yy (1) = yy (2)
csd		END IF
csd		IF  ( ( .not. ERMISS (toptmp (i) ) ) .and.
csd     +		      ( ( toptmp (i) .lt. topwnd (i) ) .or.
csd     +			( ERMISS (topwnd (i) ) ) ) )  THEN
		    yy (2) = toptmp (i)
		    CALL GQLINE  ( i1, i2, i3, i4, ier )
		    CALL GSLINE  (  3,  0, i3,  0, ier )
		    CALL GLINE   ( 'M', 2, xx, yy, ier )
		    CALL GSLINE  ( i1,  0, i3,  0, ier )
csd		END IF
	    END DO
          endif
C
C*	    Reset solid line.
C
	    CALL GSLINE  (  i1, 0, i3, 0, ier )
C
C*	    Draw surface.
C
csd	    i=0
csd	    yy (1) = ybot
csd	    DO WHILE ( i .lt. 40 )
csd		i = i + 1
csd		IF ( ERMISS ( pontha(i,1) ) ) THEN
C
C*		    Skip this point.
C
csd		ELSE
csd		    xx (1) = x (i)
csd		    xx (2) = x (i)
csd		    yy (2) = pontha (i,1)
csd		    IF ( ERMISS ( pontha (i+1,1) ) ) THEN
csd			npts = 2
csd		    ELSE
csd			xx (3) = x (i+1)
csd			yy (3) = pontha (i+1,1)
csd			npts = 3
csd		    END IF
csd	            CALL GLINE  ( 'M', npts, xx, yy, ier )
csd		END IF
csd	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: Unidata/Garp

max_issues_repo_name: Unidata/Garp

max_forks_repo_name: Unidata/Garp


text:
! from PR 15962, we used to require constant expressions instead of
! initialization expressions in case-statements
function j(k)
integer :: k
integer :: j
integer, parameter :: i(2) = (/1,2/)

select case(k)
case (1:size(i))
    j = i(k)
case default
    j = 0
end select
end function

if (j(2).NE.2 .OR. j(11).NE.0) call abort()
end