!=============================================================================!
        subroutine csolve2bc( n, nsys, nblk, mat, rhs, bc, mult, fact )
!=============================================================================!
!
!  Solves a block tridiagonal system of equations without pivoting.
!  Vectorized over the second index.
!
!  This routine has been modified to allow for high-order wall boundary
!  conditions.
!
!  WARNING: Since pivoting is not done, this algorithm could fail
!
!  WARNING: This routine works for nblk = 5 only!
!
!  input:
!    nsys       : number of systems to vectorize over
!    n          : number of rows in the systems
!    nblk       : dimension of the blocks in each row
!
!  inout:
!    mat        : the block tridiagonal matrix to be solved
!    rhs        : the right hand side
!    bc         : the bc array
!=============================================================================!
        implicit none

        integer, intent(in) :: n, nblk, nsys
        complex, intent(inout) :: mat(n,nsys,nblk,nblk*3), rhs(n,nsys,nblk)
        complex, intent(inout) :: bc(nsys,nblk,nblk,14)
        complex, intent(inout) :: mult(nsys), fact(nsys)

!.... useful constants

        real, parameter :: zero = 0.0, one = 1.0

!.... local variables

        integer :: nc, nr, i, j, l, m, p, q, s, r, t, iv
!=============================================================================!
!                     F O R W A R D   E L I M I N A T I O N
!=============================================================================!

!.... first 3 rows are special due to boundary conditions

        i = 1

        do nc = 1, nblk - 1

          j = nblk + nc

          do m = nc, nblk - 1
            l = m + 1
            mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
            do q = nc + nblk + 1, nblk * 2
              mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
            end do

            do iv = 1, nsys
              mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
              mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
              mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
              mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
              mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

              bc(iv,l,1,1) = bc(iv,l,1,1) + mult(iv) * bc(iv,nc,1,1)
              bc(iv,l,2,1) = bc(iv,l,2,1) + mult(iv) * bc(iv,nc,2,1)
              bc(iv,l,3,1) = bc(iv,l,3,1) + mult(iv) * bc(iv,nc,3,1)
              bc(iv,l,4,1) = bc(iv,l,4,1) + mult(iv) * bc(iv,nc,4,1)
              bc(iv,l,5,1) = bc(iv,l,5,1) + mult(iv) * bc(iv,nc,5,1)

              bc(iv,l,1,2) = bc(iv,l,1,2) + mult(iv) * bc(iv,nc,1,2)
              bc(iv,l,2,2) = bc(iv,l,2,2) + mult(iv) * bc(iv,nc,2,2)
              bc(iv,l,3,2) = bc(iv,l,3,2) + mult(iv) * bc(iv,nc,3,2)
              bc(iv,l,4,2) = bc(iv,l,4,2) + mult(iv) * bc(iv,nc,4,2)
              bc(iv,l,5,2) = bc(iv,l,5,2) + mult(iv) * bc(iv,nc,5,2)

              bc(iv,l,1,3) = bc(iv,l,1,3) + mult(iv) * bc(iv,nc,1,3)
              bc(iv,l,2,3) = bc(iv,l,2,3) + mult(iv) * bc(iv,nc,2,3)
              bc(iv,l,3,3) = bc(iv,l,3,3) + mult(iv) * bc(iv,nc,3,3)
              bc(iv,l,4,3) = bc(iv,l,4,3) + mult(iv) * bc(iv,nc,4,3)
              bc(iv,l,5,3) = bc(iv,l,5,3) + mult(iv) * bc(iv,nc,5,3)

              rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
            end do

          end do

        end do  ! end of loop on nc

!.... for all columns

        p = i + 1
        r = i + 2

        do nc = 1, nblk

          j = nc + nblk
          fact(:) = one / mat(i,:,nc,j)

          do m = 1, nblk

!.... second row

            mult(:) = -mat(p,:,m,nc) * fact(:)

            do q = nc + 1, nblk
              l = q + nblk
              mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
            end do

            do iv = 1, nsys
              mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
              mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
              mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
              mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
              mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

              mat(p,iv,m,11) = mat(p,iv,m,11) + mult(iv) * bc(iv,nc,1,1)
              mat(p,iv,m,12) = mat(p,iv,m,12) + mult(iv) * bc(iv,nc,2,1)
              mat(p,iv,m,13) = mat(p,iv,m,13) + mult(iv) * bc(iv,nc,3,1)
              mat(p,iv,m,14) = mat(p,iv,m,14) + mult(iv) * bc(iv,nc,4,1)
              mat(p,iv,m,15) = mat(p,iv,m,15) + mult(iv) * bc(iv,nc,5,1)

              bc(iv,m,1,4)  = bc(iv,m,1,4)  + mult(iv) * bc(iv,nc,1,2)
              bc(iv,m,2,4)  = bc(iv,m,2,4)  + mult(iv) * bc(iv,nc,2,2)
              bc(iv,m,3,4)  = bc(iv,m,3,4)  + mult(iv) * bc(iv,nc,3,2)
              bc(iv,m,4,4)  = bc(iv,m,4,4)  + mult(iv) * bc(iv,nc,4,2)
              bc(iv,m,5,4)  = bc(iv,m,5,4)  + mult(iv) * bc(iv,nc,5,2)

              bc(iv,m,1,5)  = bc(iv,m,1,5)  + mult(iv) * bc(iv,nc,1,3)
              bc(iv,m,2,5)  = bc(iv,m,2,5)  + mult(iv) * bc(iv,nc,2,3)
              bc(iv,m,3,5)  = bc(iv,m,3,5)  + mult(iv) * bc(iv,nc,3,3)
              bc(iv,m,4,5)  = bc(iv,m,4,5)  + mult(iv) * bc(iv,nc,4,3)
              bc(iv,m,5,5)  = bc(iv,m,5,5)  + mult(iv) * bc(iv,nc,5,3)

              rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
            end do

!.... third row

            mult(:) = -bc(:,m,nc,7) * fact(:)

            do q = nc + 1, nblk
              l = q + nblk
              bc(:,m,q,7) = bc(:,m,q,7) + mult(:) * mat(i,:,nc,l)
            end do

            do iv = 1, nsys
              mat(r,iv,m, 1) = mat(r,iv,m, 1) + mult(iv) * mat(i,iv,nc,11)
              mat(r,iv,m, 2) = mat(r,iv,m, 2) + mult(iv) * mat(i,iv,nc,12)
              mat(r,iv,m, 3) = mat(r,iv,m, 3) + mult(iv) * mat(i,iv,nc,13)
              mat(r,iv,m, 4) = mat(r,iv,m, 4) + mult(iv) * mat(i,iv,nc,14)
              mat(r,iv,m, 5) = mat(r,iv,m, 5) + mult(iv) * mat(i,iv,nc,15)

              mat(r,iv,m, 6) = mat(r,iv,m, 6) + mult(iv) * bc(iv,nc,1,1)
              mat(r,iv,m, 7) = mat(r,iv,m, 7) + mult(iv) * bc(iv,nc,2,1)
              mat(r,iv,m, 8) = mat(r,iv,m, 8) + mult(iv) * bc(iv,nc,3,1)
              mat(r,iv,m, 9) = mat(r,iv,m, 9) + mult(iv) * bc(iv,nc,4,1)
              mat(r,iv,m,10) = mat(r,iv,m,10) + mult(iv) * bc(iv,nc,5,1)

              mat(r,iv,m,11) = mat(r,iv,m,11) + mult(iv) * bc(iv,nc,1,2)
              mat(r,iv,m,12) = mat(r,iv,m,12) + mult(iv) * bc(iv,nc,2,2)
              mat(r,iv,m,13) = mat(r,iv,m,13) + mult(iv) * bc(iv,nc,3,2)
              mat(r,iv,m,14) = mat(r,iv,m,14) + mult(iv) * bc(iv,nc,4,2)
              mat(r,iv,m,15) = mat(r,iv,m,15) + mult(iv) * bc(iv,nc,5,2)

              bc(iv,m,1,6)  = bc(iv,m,1,6)  + mult(iv) * bc(iv,nc,1,3)
              bc(iv,m,2,6)  = bc(iv,m,2,6)  + mult(iv) * bc(iv,nc,2,3)
              bc(iv,m,3,6)  = bc(iv,m,3,6)  + mult(iv) * bc(iv,nc,3,3)
              bc(iv,m,4,6)  = bc(iv,m,4,6)  + mult(iv) * bc(iv,nc,4,3)
              bc(iv,m,5,6)  = bc(iv,m,5,6)  + mult(iv) * bc(iv,nc,5,3)

              rhs(r,iv,m) = rhs(r,iv,m) + mult(iv) * rhs(i,iv,nc)
            end do
          end do

        end do  ! end of loop on nc

!.... second special row

        i = 2

        do nc = 1, nblk - 1

          j = nblk + nc

          do m = nc, nblk - 1
            l = m + 1
            mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
            do q = nc + nblk + 1, nblk * 2
              mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
            end do

            do iv = 1, nsys
              mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
              mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
              mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
              mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
              mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

              bc(iv,l,1,4) = bc(iv,l,1,4) + mult(iv) * bc(iv,nc,1,4)
              bc(iv,l,2,4) = bc(iv,l,2,4) + mult(iv) * bc(iv,nc,2,4)
              bc(iv,l,3,4) = bc(iv,l,3,4) + mult(iv) * bc(iv,nc,3,4)
              bc(iv,l,4,4) = bc(iv,l,4,4) + mult(iv) * bc(iv,nc,4,4)
              bc(iv,l,5,4) = bc(iv,l,5,4) + mult(iv) * bc(iv,nc,5,4)

              bc(iv,l,1,5) = bc(iv,l,1,5) + mult(iv) * bc(iv,nc,1,5)
              bc(iv,l,2,5) = bc(iv,l,2,5) + mult(iv) * bc(iv,nc,2,5)
              bc(iv,l,3,5) = bc(iv,l,3,5) + mult(iv) * bc(iv,nc,3,5)
              bc(iv,l,4,5) = bc(iv,l,4,5) + mult(iv) * bc(iv,nc,4,5)
              bc(iv,l,5,5) = bc(iv,l,5,5) + mult(iv) * bc(iv,nc,5,5)

              rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
            end do
          end do

        end do  ! end of loop on nc

!.... for all columns

        do nc = 1, nblk

          p  = i + 1

          j = nc + nblk
          fact(:) = one / mat(i,:,nc,j)

          do m = 1, nblk
            mult(:) = -mat(p,:,m,nc) * fact(:)

            do q = nc + 1, nblk
              l = q + nblk
              mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
            end do

            do iv = 1, nsys
              mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
              mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
              mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
              mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
              mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

              mat(p,iv,m,11) = mat(p,iv,m,11) + mult(iv) * bc(iv,nc,1,4)
              mat(p,iv,m,12) = mat(p,iv,m,12) + mult(iv) * bc(iv,nc,2,4)
              mat(p,iv,m,13) = mat(p,iv,m,13) + mult(iv) * bc(iv,nc,3,4)
              mat(p,iv,m,14) = mat(p,iv,m,14) + mult(iv) * bc(iv,nc,4,4)
              mat(p,iv,m,15) = mat(p,iv,m,15) + mult(iv) * bc(iv,nc,5,4)

              bc(iv,m,1,6)  = bc(iv,m,1,6)  + mult(iv) * bc(iv,nc,1,5)
              bc(iv,m,2,6)  = bc(iv,m,2,6)  + mult(iv) * bc(iv,nc,2,5)
              bc(iv,m,3,6)  = bc(iv,m,3,6)  + mult(iv) * bc(iv,nc,3,5)
              bc(iv,m,4,6)  = bc(iv,m,4,6)  + mult(iv) * bc(iv,nc,4,5)
              bc(iv,m,5,6)  = bc(iv,m,5,6)  + mult(iv) * bc(iv,nc,5,5)

              rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
            end do
          end do

        end do  ! end of loop on nc

!.... the third and final special row

        i = 3

        do nc = 1, nblk - 1

          j = nblk + nc

          do m = nc, nblk - 1
            l = m + 1
            mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
            do q = nc + nblk + 1, nblk * 2
              mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
            end do

            do iv = 1, nsys
              mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
              mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
              mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
              mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
              mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

              bc(iv,l,1,6) = bc(iv,l,1,6) + mult(iv) * bc(iv,nc,1,6)
              bc(iv,l,2,6) = bc(iv,l,2,6) + mult(iv) * bc(iv,nc,2,6)
              bc(iv,l,3,6) = bc(iv,l,3,6) + mult(iv) * bc(iv,nc,3,6)
              bc(iv,l,4,6) = bc(iv,l,4,6) + mult(iv) * bc(iv,nc,4,6)
              bc(iv,l,5,6) = bc(iv,l,5,6) + mult(iv) * bc(iv,nc,5,6)

              rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
            end do
          end do

        end do  ! end of loop on nc

!.... for all columns

        do nc = 1, nblk

          p  = i + 1

          j = nc + nblk
          fact(:) = one / mat(i,:,nc,j)

          do m = 1, nblk
            mult(:) = -mat(p,:,m,nc) * fact(:)

            do q = nc + 1, nblk
              l = q + nblk
              mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
            end do

            do iv = 1, nsys
              mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
              mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
              mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
              mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
              mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

              mat(p,iv,m,11) = mat(p,iv,m,11) + mult(iv) * bc(iv,nc,1,6)
              mat(p,iv,m,12) = mat(p,iv,m,12) + mult(iv) * bc(iv,nc,2,6)
              mat(p,iv,m,13) = mat(p,iv,m,13) + mult(iv) * bc(iv,nc,3,6)
              mat(p,iv,m,14) = mat(p,iv,m,14) + mult(iv) * bc(iv,nc,4,6)
              mat(p,iv,m,15) = mat(p,iv,m,15) + mult(iv) * bc(iv,nc,5,6)

              rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
            end do
          end do

        end do  ! end of loop on nc

!.... the rest of the matrix is a standard tridiagonal up to the other boundary

        do i = 4, n - 5

!.... do the first four columns

          do nc = 1, nblk - 1

            j = nblk + nc

            do m = nc, nblk - 1
              l = m + 1
              mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
              do q = nc + nblk + 1, nblk * 2
                mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
              end do

              do iv = 1, nsys
                mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
                mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
                mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
                mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
                mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

                rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p  = i + 1

            j = nc + nblk
            fact(:) = one / mat(i,:,nc,j)

            do m = 1, nblk
              mult(:) = -mat(p,:,m,nc) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
                mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
                mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
                mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
                mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

                rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

        end do          ! end of loop on i

!.... do the FIFTH to the last row

        i = n - 4

!.... first four columns

          do nc = 1, nblk - 1

            j = nblk + nc

            do m = nc, nblk - 1
              l = m + 1
              mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
              do q = nc + nblk + 1, nblk * 2
                mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
              end do

              do iv = 1, nsys
                mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
                mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
                mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
                mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
                mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

                rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          p = i + 1
          r = i + 2
          s = i + 3
          t = i + 4

          do nc = 1, nblk

            j = nc + nblk
            fact(:) = one / mat(i,:,nc,j)

            do m = 1, nblk
              mult(:) = -mat(p,:,m,nc) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
                mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
                mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
                mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
                mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

                rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary conditions: bc(8)

              mult(:) = -bc(:,m,nc,8) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,8) = bc(:,m,q,8) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(r,iv,m, 1) = mat(r,iv,m, 1) + mult(iv) * mat(i,iv,nc,11)
                mat(r,iv,m, 2) = mat(r,iv,m, 2) + mult(iv) * mat(i,iv,nc,12)
                mat(r,iv,m, 3) = mat(r,iv,m, 3) + mult(iv) * mat(i,iv,nc,13)
                mat(r,iv,m, 4) = mat(r,iv,m, 4) + mult(iv) * mat(i,iv,nc,14)
                mat(r,iv,m, 5) = mat(r,iv,m, 5) + mult(iv) * mat(i,iv,nc,15)

                rhs(r,iv,m) = rhs(r,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary conditions: bc(10)

              mult(:) = -bc(:,m,nc,10) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,10) = bc(:,m,q,10) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                bc(iv,m, 1,11) = bc(iv,m, 1,11) + mult(iv) * mat(i,iv,nc,11)
                bc(iv,m, 2,11) = bc(iv,m, 2,11) + mult(iv) * mat(i,iv,nc,12)
                bc(iv,m, 3,11) = bc(iv,m, 3,11) + mult(iv) * mat(i,iv,nc,13)
                bc(iv,m, 4,11) = bc(iv,m, 4,11) + mult(iv) * mat(i,iv,nc,14)
                bc(iv,m, 5,11) = bc(iv,m, 5,11) + mult(iv) * mat(i,iv,nc,15)

                rhs(s,iv,m) = rhs(s,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary conditions: bc(12)

              mult(:) = -bc(:,m,nc,12) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,12) = bc(:,m,q,12) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                bc(iv,m, 1,13) = bc(iv,m, 1,13) + mult(iv) * mat(i,iv,nc,11)
                bc(iv,m, 2,13) = bc(iv,m, 2,13) + mult(iv) * mat(i,iv,nc,12)
                bc(iv,m, 3,13) = bc(iv,m, 3,13) + mult(iv) * mat(i,iv,nc,13)
                bc(iv,m, 4,13) = bc(iv,m, 4,13) + mult(iv) * mat(i,iv,nc,14)
                bc(iv,m, 5,13) = bc(iv,m, 5,13) + mult(iv) * mat(i,iv,nc,15)

                rhs(t,iv,m) = rhs(t,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

            end do

          end do        ! end of loop on nc

!.... do the FOURTH to the last row

        i = n - 3

!.... first four columns

          do nc = 1, nblk - 1

            j = nblk + nc

            do m = nc, nblk - 1
              l = m + 1
              mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
              do q = nc + nblk + 1, nblk * 2
                mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
              end do

              do iv = 1, nsys
                mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
                mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
                mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
                mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
                mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

                rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          p = i + 1
          r = i + 2
          s = i + 3

          do nc = 1, nblk

            j = nc + nblk
            fact(:) = one / mat(i,:,nc,j)

            do m = 1, nblk
              mult(:) = -mat(p,:,m,nc) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
                mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
                mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
                mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
                mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

                rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary conditions: bc(11)

              mult(:) = -bc(:,m,nc,11) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,11) = bc(:,m,q,11) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(r,iv,m, 1) = mat(r,iv,m, 1) + mult(iv) * mat(i,iv,nc,11)
                mat(r,iv,m, 2) = mat(r,iv,m, 2) + mult(iv) * mat(i,iv,nc,12)
                mat(r,iv,m, 3) = mat(r,iv,m, 3) + mult(iv) * mat(i,iv,nc,13)
                mat(r,iv,m, 4) = mat(r,iv,m, 4) + mult(iv) * mat(i,iv,nc,14)
                mat(r,iv,m, 5) = mat(r,iv,m, 5) + mult(iv) * mat(i,iv,nc,15)

                rhs(r,iv,m) = rhs(r,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary conditions: bc(13)

              mult(:) = -bc(:,m,nc,13) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,13) = bc(:,m,q,13) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                bc(iv,m, 1,14) = bc(iv,m, 1,14) + mult(iv) * mat(i,iv,nc,11)
                bc(iv,m, 2,14) = bc(iv,m, 2,14) + mult(iv) * mat(i,iv,nc,12)
                bc(iv,m, 3,14) = bc(iv,m, 3,14) + mult(iv) * mat(i,iv,nc,13)
                bc(iv,m, 4,14) = bc(iv,m, 4,14) + mult(iv) * mat(i,iv,nc,14)
                bc(iv,m, 5,14) = bc(iv,m, 5,14) + mult(iv) * mat(i,iv,nc,15)

                rhs(s,iv,m) = rhs(s,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... do the THIRD to the last row

        i = n - 2

!.... first four columns

          do nc = 1, nblk - 1

            j = nblk + nc

            do m = nc, nblk - 1
              l = m + 1
              mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
              do q = nc + nblk + 1, nblk * 2
                mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
              end do

              do iv = 1, nsys
                mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
                mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
                mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
                mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
                mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

                bc(iv,l, 1,9) = bc(iv,l, 1,9) + mult(iv) * bc(iv,nc,1,9)
                bc(iv,l, 2,9) = bc(iv,l, 2,9) + mult(iv) * bc(iv,nc,2,9)
                bc(iv,l, 3,9) = bc(iv,l, 3,9) + mult(iv) * bc(iv,nc,3,9)
                bc(iv,l, 4,9) = bc(iv,l, 4,9) + mult(iv) * bc(iv,nc,4,9)
                bc(iv,l, 5,9) = bc(iv,l, 5,9) + mult(iv) * bc(iv,nc,5,9)

                rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          p = i + 1
          r = i + 2

          do nc = 1, nblk

            j = nc + nblk
            fact(:) = one / mat(i,:,nc,j)

            do m = 1, nblk
              mult(:) = -mat(p,:,m,nc) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
                mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
                mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
                mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
                mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

                mat(p,iv,m,11) = mat(p,iv,m,11) + mult(iv) * bc(iv,nc,1,9)
                mat(p,iv,m,12) = mat(p,iv,m,12) + mult(iv) * bc(iv,nc,2,9)
                mat(p,iv,m,13) = mat(p,iv,m,13) + mult(iv) * bc(iv,nc,3,9)
                mat(p,iv,m,14) = mat(p,iv,m,14) + mult(iv) * bc(iv,nc,4,9)
                mat(p,iv,m,15) = mat(p,iv,m,15) + mult(iv) * bc(iv,nc,5,9)

                rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do

!.... eliminate the boundary condition:  bc(14)

              mult(:) = -bc(:,m,nc,14) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                bc(:,m,q,14) = bc(:,m,q,14) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(r,iv,m, 1) = mat(r,iv,m, 1) + mult(iv) * mat(i,iv,nc,11)
                mat(r,iv,m, 2) = mat(r,iv,m, 2) + mult(iv) * mat(i,iv,nc,12)
                mat(r,iv,m, 3) = mat(r,iv,m, 3) + mult(iv) * mat(i,iv,nc,13)
                mat(r,iv,m, 4) = mat(r,iv,m, 4) + mult(iv) * mat(i,iv,nc,14)
                mat(r,iv,m, 5) = mat(r,iv,m, 5) + mult(iv) * mat(i,iv,nc,15)

                mat(r,iv,m, 6) = mat(r,iv,m, 6) + mult(iv) * bc(iv,nc,1,9)
                mat(r,iv,m, 7) = mat(r,iv,m, 7) + mult(iv) * bc(iv,nc,2,9)
                mat(r,iv,m, 8) = mat(r,iv,m, 8) + mult(iv) * bc(iv,nc,3,9)
                mat(r,iv,m, 9) = mat(r,iv,m, 9) + mult(iv) * bc(iv,nc,4,9)
                mat(r,iv,m,10) = mat(r,iv,m,10) + mult(iv) * bc(iv,nc,5,9)

                rhs(r,iv,m) = rhs(r,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... do the SECOND to the last row (which is normal tridiagonal)

        i = n - 1

!.... first four columns

          do nc = 1, nblk - 1

            j = nblk + nc

            do m = nc, nblk - 1
              l = m + 1
              mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
              do q = nc + nblk + 1, nblk * 2
                mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
              end do

              do iv = 1, nsys
                mat(i,iv,l,11) = mat(i,iv,l,11) + mult(iv) * mat(i,iv,nc,11)
                mat(i,iv,l,12) = mat(i,iv,l,12) + mult(iv) * mat(i,iv,nc,12)
                mat(i,iv,l,13) = mat(i,iv,l,13) + mult(iv) * mat(i,iv,nc,13)
                mat(i,iv,l,14) = mat(i,iv,l,14) + mult(iv) * mat(i,iv,nc,14)
                mat(i,iv,l,15) = mat(i,iv,l,15) + mult(iv) * mat(i,iv,nc,15)

                rhs(i,iv,l) = rhs(i,iv,l) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... for all columns

          p = i + 1

          do nc = 1, nblk

            j = nc + nblk
            fact(:) = one / mat(i,:,nc,j)

            do m = 1, nblk
              mult(:) = -mat(p,:,m,nc) * fact(:)
              do q = nc + 1, nblk
                l = q + nblk
                mat(p,:,m,q) = mat(p,:,m,q) + mult(:) * mat(i,:,nc,l)
              end do

              do iv = 1, nsys
                mat(p,iv,m, 6) = mat(p,iv,m, 6) + mult(iv) * mat(i,iv,nc,11)
                mat(p,iv,m, 7) = mat(p,iv,m, 7) + mult(iv) * mat(i,iv,nc,12)
                mat(p,iv,m, 8) = mat(p,iv,m, 8) + mult(iv) * mat(i,iv,nc,13)
                mat(p,iv,m, 9) = mat(p,iv,m, 9) + mult(iv) * mat(i,iv,nc,14)
                mat(p,iv,m,10) = mat(p,iv,m,10) + mult(iv) * mat(i,iv,nc,15)

                rhs(p,iv,m) = rhs(p,iv,m) + mult(iv) * rhs(i,iv,nc)
              end do
            end do

          end do        ! end of loop on nc

!.... do the last row

        i = n

        do nc = 1, nblk - 1

          j = nblk + nc

          do m = nc, nblk - 1
            l = m + 1
            mult(:) = -mat(i,:,l,j) / mat(i,:,nc,j)
            do q = nc + nblk + 1, nblk * 2
              mat(i,:,l,q) = mat(i,:,l,q) + mult(:) * mat(i,:,nc,q)
            end do
            rhs(i,:,l) = rhs(i,:,l) + mult(:) * rhs(i,:,nc)
          end do

        end do  ! end of loop on nc

!=============================================================================!
!                   B A C K W A R D   S U B S T I T U T I O N
!=============================================================================!

!.... do the LAST row first

        i = n

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do
          rhs(i,:,m) = rhs(i,:,m) / mat(i,:,m,j)
        end do

!.... do the SECOND to the last row (normal row)

        i = n - 1
        p = i + 1

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do

          do iv = 1, nsys
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

            rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
          end do
        end do

!.... do the THIRD to the last row

        i = n - 2
        p = i + 1
        r = i + 2

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do

          do iv = 1, nsys
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,1) * bc(iv,m,1,9)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,2) * bc(iv,m,2,9)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,3) * bc(iv,m,3,9)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,4) * bc(iv,m,4,9)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,5) * bc(iv,m,5,9)

            rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
          end do

        end do

!.... now do the rest of the rows in reverse order upto the boundary

        do i = n - 3, 4, -1

          p = i + 1

          do m = nblk, 1, -1
            j = m + nblk
            do q = m+1, nblk
              l = q + nblk
              rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
            end do

            do iv = 1, nsys
              rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
              rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
              rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
              rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
              rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

              rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
            end do
          end do

        end do

!.... now do the boundary rows

        i = 3

        p = i + 1
        r = i + 2

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do

          do iv = 1, nsys
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,1) * bc(iv,m,1,6)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,2) * bc(iv,m,2,6)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,3) * bc(iv,m,3,6)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,4) * bc(iv,m,4,6)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,5) * bc(iv,m,5,6)

            rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
          end do
        end do

        i = 2

        p = i + 1
        r = i + 2
        s = i + 3

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do

          do iv = 1, nsys
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,1) * bc(iv,m,1,4)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,2) * bc(iv,m,2,4)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,3) * bc(iv,m,3,4)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,4) * bc(iv,m,4,4)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,5) * bc(iv,m,5,4)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,1) * bc(iv,m,1,5)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,2) * bc(iv,m,2,5)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,3) * bc(iv,m,3,5)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,4) * bc(iv,m,4,5)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,5) * bc(iv,m,5,5)

            rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
          end do
        end do

        i = 1

        p = i + 1
        r = i + 2
        s = i + 3
        t = i + 4

        do m = nblk, 1, -1
          j = m + nblk
          do q = m+1, nblk
            l = q + nblk
            rhs(i,:,m) = rhs(i,:,m) - rhs(i,:,q) * mat(i,:,m,l)
          end do

          do iv = 1, nsys
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,1) * mat(i,iv,m,11)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,2) * mat(i,iv,m,12)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,3) * mat(i,iv,m,13)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,4) * mat(i,iv,m,14)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(p,iv,5) * mat(i,iv,m,15)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,1) * bc(iv,m,1,1)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,2) * bc(iv,m,2,1)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,3) * bc(iv,m,3,1)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,4) * bc(iv,m,4,1)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(r,iv,5) * bc(iv,m,5,1)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,1) * bc(iv,m,1,2)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,2) * bc(iv,m,2,2)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,3) * bc(iv,m,3,2)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,4) * bc(iv,m,4,2)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(s,iv,5) * bc(iv,m,5,2)

            rhs(i,iv,m) = rhs(i,iv,m) - rhs(t,iv,1) * bc(iv,m,1,3)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(t,iv,2) * bc(iv,m,2,3)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(t,iv,3) * bc(iv,m,3,3)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(t,iv,4) * bc(iv,m,4,3)
            rhs(i,iv,m) = rhs(i,iv,m) - rhs(t,iv,5) * bc(iv,m,5,3)

            rhs(i,iv,m) = rhs(i,iv,m) / mat(i,iv,m,j)
          end do
        end do
!=============================================================================!
        return
        end


!***********************************/
!	Name:FDS法による1次元風上流束を求めるプログラム
!	Alias:UpwindFlux_Dim1
!	Description:FDS法を用いて1次元系における風上流束を求める
!	Type:CE
!	Input:Geometory,CellEdge,CellCenter
!	Output:CellEdge
!	Note:次元ごとにサブルーチンを分けているのは1つだけでも割と長いプログラムになってしまったため
!	Author:Akitaka Toyota
!	Date:2017.10.26
!	Update:2017.11.09
!	Other:1次元と2次元に対応
!***********************************/
    subroutine UpwindFlux_Dim1(Geom,CE,CC) !MUSCL経由の場合CEのみ，1次精度の場合CCを渡す
    use StructVar_Mod
    use LoopVar_Mod
    use ConstantVar_Mod, Gamma => SpecificOfHeatRatio, Epsilon => EntropyCorrection
    implicit none
    type(Geometry), intent(in) :: Geom
    type(CellEdge), intent(inout) :: CE
    type(CellCenter), intent(in), optional :: CC

    type(FDSMatrix) :: FM
    type(RoeAverage) :: RA
    double precision, allocatable :: DeltaQuantity(:)
    double precision, allocatable :: RoeFlux(:)
    integer iTestMode
    iTestMode = 0
    allocate(RA%GeneralFlux(3,1,2), DeltaQuantity(3),RoeFlux(3))
    allocate(FM%RightEigenMatrix(3,3),FM%LeftEigenMatrix(3,3),FM%AbsEigenValues(3,3))
    allocate(RA%RoeAverage(3),RA%SideQuantity(3,2))

    !Roeスキームの流束を求める
    do iFaceX=0,Geom%CellNumber(1) !すべての界面についての総和
        if( present(CC)) then !保存変数からの流束計算
            RA%SideQuantity(:,1) = CC%ConservedQuantity(:,iCX,1,1) !負側の保存量
            RA%SideQuantity(:,2) = CC%ConservedQuantity(:,iCX+1,1,1) !正側の保存量
            RA%TypeOfQuantity = 1
        else  !基礎変数からの流束計算
            RA%SideQuantity(1,1) = CE%RebuildQunatity(1,iCX,1,1,2)
            RA%SideQuantity(2,1) = CE%RebuildQunatity(2,iCX,1,1,2)*CE%RebuildQunatity(1,iCX,1,1,2)
            RA%SideQuantity(3,1) = InverseGmin1*CE%RebuildQunatity(3,iCX,1,1,2) &
                &   + 0.5d0*CE%RebuildQunatity(1,iCX,1,1,2)*CE%RebuildQunatity(2,iCX,1,1,2)**2

            RA%SideQuantity(1,2) = CE%RebuildQunatity(1,iCX+1,1,1,1)
            RA%SideQuantity(2,2) = CE%RebuildQunatity(2,iCX+1,1,1,1)*CE%RebuildQunatity(1,iCX+1,1,1,1)
            RA%SideQuantity(3,2) = InverseGmin1*CE%RebuildQunatity(3,iCX+1,1,1,1)  &
                &   + 0.5d0*CE%RebuildQunatity(1,iCX+1,1,1,1)*CE%RebuildQunatity(2,iCX+1,1,1,1)**2
            RA%TypeOfQuantity = 1
        end if

        DeltaQuantity(:) = RA%SideQuantity(:,2) - RA%SideQuantity(:,1)
        call GetGeneralFlux(Geom,RA)

        !call GetRoeAverage(Geom,RA)

!左右の固有ベクトルから風上化流束を求める
        !call GetEigenMatrix(Geom,RA,FM)

!エントロピー補正
        do iLoop =1,3
        if(FM%AbsEigenValues(iLoop,iLoop) < Epsilon) then
            FM%AbsEigenValues(iLoop,iLoop) = 0.5d0*(FM%AbsEigenValues(iLoop,iLoop)**2 + Epsilon**2)/Epsilon
        end if
        end do
        !TestCode
        iTestMode = 0
        if(iTestMode == 0) then
        RoeFlux = matmul(matmul(matmul(FM%RightEigenMatrix,FM%AbsEigenValues),FM%LeftEigenMatrix),DeltaQuantity)
        else
            call RoeMultiply
        end if
        CE%NormalFluxDiff(:,iCX,1,1,1) = &
            &   0.5d0 * ((RA%GeneralFlux(:,1,1)+RA%GeneralFlux(:,1,2)) - RoeFlux(:))

    end do
    return
contains
    subroutine RoeMultiply
    implicit none
    double precision :: CalcAssistB2,CalcAssistB1
        CalcAssistB2 = Gmin1/RA%RoeAverage(1)**2
        CalcAssistB1 = 0.5d0*RA%RoeAverage(2)**2 * CalcAssistB2

        RoeFlux(1) = &
        & (0.5d0*(CalcAssistB1 + RA%RoeAverage(2)/RA%RoeAverage(1)) &
        &+(1.0d0 - CalcAssistB1) &
        &+0.5d0*(CalcAssistB1 - RA%RoeAverage(2)/RA%RoeAverage(1)))*DeltaQuantity(1) &
        &*dabs(RA%RoeAverage(2) - RA%RoeAverage(1))

        RoeFlux(2) = &
        & ((RA%RoeAverage(2) - RA%RoeAverage(1))&
        &*(-0.5d0*(1.0d0/RA%RoeAverage(1) + CalcAssistB2*RA%RoeAverage(2))) &
        &+RA%RoeAverage(2)*CalcAssistB2 * RA%RoeAverage(2) &
        &+(RA%RoeAverage(2) + RA%RoeAverage(1))&
        &*0.5d0*(1.0d0/RA%RoeAverage(1) - CalcAssistB2*RA%RoeAverage(2))) &
        &*DeltaQuantity(2)*dabs(RA%RoeAverage(2))

        RoeFlux(3) = &
        & ((RA%RoeAverage(3) - RA%RoeAverage(2)*RA%RoeAverage(1))*0.5d0*CalcAssistB2 &
        &+0.5d0 * RA%RoeAverage(2)**2 *(-CalcAssistB2) &
        &+(RA%RoeAverage(3) + RA%RoeAverage(2)*RA%RoeAverage(1))*0.5d0*CalcAssistB2) &
        &*DeltaQuantity(3)*dabs(RA%RoeAverage(2) + RA%RoeAverage(1))
    return
    end subroutine RoeMultiply

    end subroutine UpwindFlux_Dim1


text:
  4   IF (ICHAR(STORE(D)).NE.0) GO TO 3
      IF (ICHAR(STORE(D)).NE.0) GO TO 3

ext: f

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text: 20071019 15:52:39 nbsp hello. thanks for adding info to the Tennis tennis pages. if youd like to hit sometime, let me know and we can play on campus. Users/JessicaRockwell

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine checkversion(inpunit,filename)
************************************************************************
*  Checks that the version of MCFM specified in the next line of unit  *
*  "inpunit" agrees with the version number of the code                *
************************************************************************
      implicit none
      include 'codeversion.f'
      integer inpunit,j,dat
      character*6 fileversion
      character*(*) filename

      read(inpunit,*) fileversion

      if (fileversion .ne. codeversion) then
        dat=18
        do j=1,20
          if (filename(j:j) .eq. 'D') dat=j
        enddo
        write(6,*)
        write(6,*) 'Sorry, the version of this input file does not'
        write(6,*) 'match with the code version number. Please refer'
        write(6,*) 'to the documentation and adjust accordingly.'
        write(6,*)
        write(6,*) '     Filename: ',filename(1:dat+2)
        write(6,*) ' File version: ',fileversion
        write(6,*) ' Code version: ',codeversion
        write(6,*)
        stop
      endif

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
    module CalcLike_Cosmology
    use CalcLike
    use DataLikelihoodList
    use CosmologyTypes
    use CosmoTheory
    use Calculator_Cosmology
    use ParamPointSet
    implicit none

    private

    type, extends(TTheoryLikeCalculator) :: TCosmoLikeCalculator
        integer :: semislow_changes=0, slow_changes=0
        class(TCosmologyCalculator), pointer :: CosmoCalc
    contains
    procedure :: CalculateRequiredTheoryChanges =>Cosmo_CalculateRequiredTheoryChanges
    procedure :: GetTheoryForLike=>Cosmo_GetTheoryForLike
    procedure :: WritePerformanceStats => Cosmo_WritePerformanceStats
    procedure :: InitConfig => TCosmoLikeCalculator_InitConfig
    procedure :: UpdateTheoryForLikelihoods => Cosmo_UpdateTheoryForLikelihoods
    end Type TCosmoLikeCalculator

    public TCosmoLikeCalculator
    contains

    subroutine TCosmoLikeCalculator_InitConfig(this, Config)
    class(TCosmoLikeCalculator) :: this
    class(TGeneralConfig), target :: config

    this%Config => config
    select type (p=>config%Calculator)
    class is (TCosmologyCalculator)
        this%CosmoCalc => p
    end select
    end subroutine TCosmoLikeCalculator_InitConfig

    subroutine Cosmo_GetTheoryForLike(this,Like)
    class(TCosmoLikeCalculator) :: this
    class(TDataLikelihood), pointer :: like
    logical, save :: backgroundSet

    if (.not. associated(like)) then
        !initialize
        backgroundSet = this%slowChanged
        return
    end if

    if (any(like%dependent_params(1:num_hard)) .and. .not. backgroundSet) then
        select type (CMB=>this%TheoryParams)
        class is (CMBParams)
            call this%CosmoCalc%SetParamsForBackground(CMB)
        end select
        backgroundSet = .true.
    end if

    end subroutine Cosmo_GetTheoryForLike

    logical function Cosmo_CalculateRequiredTheoryChanges(this)
    ! EFTCOSMOMC MOD START: use the required modules
    use EFTinitialization, only: EFTCAMB_initialization
    use CAMB, only: CAMBparams, CAMBParams_Set
    use ModelParams, only: feedbacklevel
    ! EFTCOSMOMC MOD END.
    class(TCosmoLikeCalculator) :: this
    integer error
    Type(TTimer) :: Timer

    ! EFTCOSMOMC MOD START: define
    logical :: EFTsuccess
    Type(CAMBParams) :: EFT_CP_Temp
    real(mcp) :: Haux, phi_aux
    ! EFTCOSMOMC MOD END.

    this%SlowChanged = any(this%changeMask(1:num_hard))
    this%SemiSlowchanged = any(this%changeMask(index_initpower:index_initpower+num_initpower-1))
    error=0
    if (this%timing) call Timer%Start

    select type (Theory=>this%Params%Theory)
    class is (TCosmoTheoryPredictions)
        select type (CMB=>this%TheoryParams)
        class is (CMBParams)
            if (CosmoSettings%Use_CMB .or. CosmoSettings%Use_LSS .or. CosmoSettings%get_sigma8) then
                if (this%SlowChanged .or. this%SemiSlowchanged .and. .not. BaseParams%block_semi_fast) then
                    ! EFTCOSMOMC MOD START: run the stability check for the choosen parameter.
                    call this%CosmoCalc%CMBToCAMB(CMB,EFT_CP_Temp)
                    call CAMBParams_Set(EFT_CP_Temp,Haux,phi_aux)
                    ! Call EFTCAMB initialization

                    if (EFT_CP_Temp%EFTflag/=0) then
                        if (Feedback>2) feedbacklevel = 3
                        call EFTCAMB_initialization(EFTsuccess)
                        !Modified for EFT start
                        if (CosmoSettings%EFTFlag .eq. 2 .and. (CosmoSettings%designerEFTmodel .eq. 4 .or. CosmoSettings%designerEFTmodel .eq. 5)) then
                             EFTsuccess = .true.
                        end if
                        !Modified for EFT end
                        if (.not.EFTsuccess) then
                            Cosmo_CalculateRequiredTheoryChanges = .false.
                            return
                        end if
                    end if
                    ! EFTCOSMOMC MOD END.
                    this%slow_changes = this%slow_changes + 1
                    this%Params%validInfo = .false.
                    call this%CosmoCalc%GetNewTransferData(CMB, this%Params%Info,Theory, error)
                end if
                if ((this%SlowChanged .or. this%SemiSlowchanged) .and. error==0) then
                    if (.not. this%SlowChanged) this%semislow_changes = this%semislow_changes  + 1
                    call this%CosmoCalc%GetNewPowerData(CMB, this%Params%Info, Theory,error)
                end if
            else
                if (this%SlowChanged) then

                    ! EFTCOSMOMC MOD START: run the stability check for the choosen parameter.
                    call this%CosmoCalc%CMBToCAMB(CMB,EFT_CP_Temp)
                    call CAMBParams_Set(EFT_CP_Temp,Haux,phi_aux)
                    ! Call EFTCAMB initialization
                    if (EFT_CP_Temp%EFTflag/=0) then
                        if (Feedback>2) feedbacklevel = 3
                        call EFTCAMB_initialization(EFTsuccess)
                        if (.not.EFTsuccess) then
                            Cosmo_CalculateRequiredTheoryChanges = .false.
                            return
                        end if
                    end if
                    ! EFTCOSMOMC MOD END.

                    call this%CosmoCalc%GetNewBackgroundData(CMB, Theory, error)
                end if
            end if
        end select
    end select

    if (this%timing) call Timer%WriteTime('Time for theory')

    this%Params%validInfo = .true.
    Cosmo_CalculateRequiredTheoryChanges = error==0

    end function Cosmo_CalculateRequiredTheoryChanges

    subroutine Cosmo_WritePerformanceStats(this, unit)
    class(TCosmoLikeCalculator) :: this
    integer, intent(in) :: unit

    write(unit,*) 'slow changes', this%slow_changes, ' power changes', this%semislow_changes

    end subroutine Cosmo_WritePerformanceStats

    subroutine Cosmo_UpdateTheoryForLikelihoods(this, Params)
    class(TCosmoLikeCalculator) :: this
    class(TCalculationAtParamPoint) :: Params
    integer :: index_error

    select type (Theory=>Params%Theory)
    class is (TCosmoTheoryPredictions)
        if (CosmoSettings%Use_LSS) then
            if(Theory%sigma_8==0) &
                call MpiStop('ERROR: Matter power/sigma_8 have not been computed. Use redo_theory and redo_pk')
            if (allocated(Theory%MPK)) then
                if((CosmoSettings%power_redshifts(CosmoSettings%num_power_redshifts)-Theory%MPK%y(Theory%MPK%ny))>1.d-3)then
                    write(*,*) 'ERROR: Thes elected datasets call for a higher redshift than has been calculated'
                    write(*,*) '       Use redo_theory and redo_pk'
                    call MpiStop()
                end if
                if(CosmoSettings%num_power_redshifts > Theory%MPK%ny)then
                    write(*,*) 'ERROR: The selected datasets call for more redshifts than are calculated'
                    write(*,*) '       Use redo_theory and redo_pk'
                    call MpiStop()
                end if
            end if
            index_error =0
            !AL for the moment only allow reading with matching mpk
            !call IndexExactRedshifts(Theory%MPK%y,index_error)
            !if(index_error>0)then
            !    write(*,*) 'ERROR: One of the datasets needs an exact redshift that is not present '
            !    write(*,*) '       Use redo_theory and redo_pk'
            !    call MpiStop()
            !end if
            if(CosmoSettings%use_nonlinear .and. .not. allocated(Theory%NL_MPK))then
                write(*,*) 'ERROR: One of the datasets wants a non-linear MPK which is not present '
                write(*,*) '       Use redo_theory and redo_pk or turn off non-linear'
                call MpiStop()
            end if
        end if
    end select

    end subroutine Cosmo_UpdateTheoryForLikelihoods

    end module CalcLike_Cosmology

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD


text:
!! NOTE: this file is autogenerated from adjac_fft.f95.in: do not edit manually
module adjac_fft
  use adjac

  private
  public fft, ifft

  interface fft
     module procedure fft_a, fft_d, fft_b, fft_z
  end interface fft

  interface ifft
     module procedure ifft_a, ifft_d, ifft_b, ifft_z
  end interface ifft

  type(adjac_double), dimension(:), allocatable :: ch_a
  double precision, dimension(:), allocatable :: wa, ch_d
  integer, dimension(15) :: ifac
contains
  subroutine fft_a(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    type(adjac_double), dimension(:), contiguous, intent(inout) :: q
    integer :: n

    n = size(q)/2
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_a) .or. size(ch_a) .ne. 2*n) then
       if (allocated(ch_a)) deallocate(ch_a)
       allocate(ch_a(2*n))
    end if
    call zfftf1a(n, q, ch_a, wa, ifac)
  end subroutine fft_a

  subroutine ifft_a(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    type(adjac_double), dimension(:), contiguous, intent(inout) :: q
    integer :: i, n

    n = size(q)/2
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_a) .or. size(ch_a) .ne. 2*n) then
       if (allocated(ch_a)) deallocate(ch_a)
       allocate(ch_a(2*n))
    end if
    call zfftb1a(n, q, ch_a, wa, ifac)

    do i = 1, size(q)
       q(i) = q(i) / n
    end do
  end subroutine ifft_a
  subroutine fft_d(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    double precision, dimension(:), contiguous, intent(inout) :: q
    integer :: n

    n = size(q)/2
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_d) .or. size(ch_d) .ne. 2*n) then
       if (allocated(ch_d)) deallocate(ch_d)
       allocate(ch_d(2*n))
    end if
    call zfftf1d(n, q, ch_d, wa, ifac)
  end subroutine fft_d

  subroutine ifft_d(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    double precision, dimension(:), contiguous, intent(inout) :: q
    integer :: i, n

    n = size(q)/2
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_d) .or. size(ch_d) .ne. 2*n) then
       if (allocated(ch_d)) deallocate(ch_d)
       allocate(ch_d(2*n))
    end if
    call zfftb1d(n, q, ch_d, wa, ifac)

    do i = 1, size(q)
       q(i) = q(i) / n
    end do
  end subroutine ifft_d
  subroutine fft_b(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    type(adjac_complex), dimension(:), contiguous, intent(inout) :: q
    integer :: n

    n = size(q)
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_a) .or. size(ch_a) .ne. 2*n) then
       if (allocated(ch_a)) deallocate(ch_a)
       allocate(ch_a(2*n))
    end if
    call zfftf1a(n, q, ch_a, wa, ifac)
  end subroutine fft_b

  subroutine ifft_b(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    type(adjac_complex), dimension(:), contiguous, intent(inout) :: q
    integer :: i, n

    n = size(q)
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_a) .or. size(ch_a) .ne. 2*n) then
       if (allocated(ch_a)) deallocate(ch_a)
       allocate(ch_a(2*n))
    end if
    call zfftb1a(n, q, ch_a, wa, ifac)

    do i = 1, size(q)
       q(i) = q(i) / n
    end do
  end subroutine ifft_b
  subroutine fft_z(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    complex(kind=kind(0d0)), dimension(:), contiguous, intent(inout) :: q
    integer :: n

    n = size(q)
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_d) .or. size(ch_d) .ne. 2*n) then
       if (allocated(ch_d)) deallocate(ch_d)
       allocate(ch_d(2*n))
    end if
    call zfftf1d(n, q, ch_d, wa, ifac)
  end subroutine fft_z

  subroutine ifft_z(q)
    use adjac
    implicit none
    external :: zffti1, zfftf1
    complex(kind=kind(0d0)), dimension(:), contiguous, intent(inout) :: q
    integer :: i, n

    n = size(q)
    if (n.le.1) return

    if (.not.allocated(wa) .or. size(wa) .ne. 2*n) then
       if (allocated(wa)) deallocate(wa)
       allocate(wa(2*n))
       call zffti1(n, wa, ifac)
    end if
    if (.not.allocated(ch_d) .or. size(ch_d) .ne. 2*n) then
       if (allocated(ch_d)) deallocate(ch_d)
       allocate(ch_d(2*n))
    end if
    call zfftb1d(n, q, ch_d, wa, ifac)

    do i = 1, size(q)
       q(i) = q(i) / n
    end do
  end subroutine ifft_z
end module adjac_fft

max_stars_repo_name: gyzhangqm/adjac

max_issues_repo_name: pv/adjac

max_forks_repo_name: pv/adjac


text:
      SUBROUTINE ZGLOBC(PASS,IDCCRC,ROW,CARPOS,TNAMES,TCS,JD,TEMPS,
     &     CAP,COEF,ISTAT)
*
*  Module number:
*
*  Module name: zglobc
*
*  Keyphrase:
*  ---------
*     Calculate dispersion coefficients from global coefficients.
*
*  Description:
*  -----------
*     This routine reads the global coefficients from the CCRC table,
*     calculates the dispersion coefficients for the equation:
*
*     s = a0 + a1*m*w + a2*m*m*w*w + a3*m + a4*w +
*         a5*m*m*w + a6*m*w*w + a7*m*m*m*w*w*w
*
*     and applies the following dispersion corrections:
*
*     COR(1) = TCS(1)*(TEMPS(1)-REFTMP(1)) + TCS(2)*(TEMPS(2)-REMTMP(2))
*     COR(2) = JD(3)*(obstime - JD(1))
*     COR(3) = TCS(3)*(TEMPS(3)-REFTMP(3))
*
*     F00 = COR(1) + COR(2)
*     F10 = COR(3)
*
*  FORTRAN name: ztherm.f
*
*  Keywords of accessed files and tables:
*  -------------------------------------
*  Name                          I/O    Description / Comments
*     CCRC                       I      Dispersion coef. table
*
*  Subroutines Called:
*
*  History:
*  --------
*  Version      Date        Author          Description
*
*    2.0        Feb 98     M.D. De La Pena  Fixed calls to UTCFND & UTRGTD.
*                                           An array was being passed for
*                                           the status when a scalar was
*                                           required. Put the calls inside the
*                                           loop so each call can be checked.
*
*---------------------------------------------------------------------------
*
*  INPUTS:
*     PASS   - How many times this routine is called.
*     IDCCRC - Table descriptor for CCRC.
*     ROW    - Row number from CCRC table to read values from.
*     CARPOS - Carrousel position
*     TNAMES - Names of the temperature monitors to use.
*     TCS    - Thermal scale factors.
*     JD     - Time scale factors.
*     CAP    - Values of CAP_A and CAP_C.
*
*  OUTPUTS:
*     COEF   - The dispersion coefficients.
*     ISTAT  - If 0 all is OK.
*
*---------------------------------------------------------------------------
      INTEGER PASS,IDCCRC,ROW,CARPOS,ISTAT
      CHARACTER*8 TNAMES(3)
      DOUBLE PRECISION TCS(3),JD(3),CAP(2),TEMPS(3)
      DOUBLE PRECISION COEF(8)
C
C     ZMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Number of columns to read from the table.
C
      INTEGER NCOLS
      PARAMETER (NCOLS=19)
C
      INTEGER FIRST
      PARAMETER (FIRST = 1)
      INTEGER LAST
      PARAMETER (LAST = -1)
C
C Local variables
C
      CHARACTER*15 COLNAM(NCOLS)
      INTEGER COLIDS(NCOLS),TSTAT,I
      CHARACTER*161 CONTXT
      DOUBLE PRECISION GLOBAL(NCOLS)
      DOUBLE PRECISION CPOS,CPOS2,COR(3)
      DOUBLE PRECISION F0,F1,F2,F3,F4,F5
      DOUBLE PRECISION MCEN,WCEN,K,K2,K3
      LOGICAL NULL(NCOLS)
C
      DATA COLNAM/'MCENTER','F00','F01','F02',
     &     'F10','F11','F12','F20','F21','F22','F30','F31','F32',
     &     'F40','F41','F42','F50','F51','F52'/
C---------------------------------------------------------------------------
C
C Initializations
C
      IF(PASS.EQ.FIRST)THEN
C
C Read in the global coefficients from the table.
C
         TSTAT = 0
         DO 1000 I = 1, NCOLS
            CALL UTCFND(IDCCRC,COLNAM(I),1,COLIDS(I),TSTAT)
            IF(TSTAT.NE.0)THEN
               WRITE(CONTXT,9)COLNAM(I)
 9             FORMAT('ERROR: Could not find column ',a15,
     &              ' in CCRC table')
               CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
               ISTAT = TSTAT
               GO TO 999
            ENDIF
 1000    CONTINUE
         DO 2000 I = 1, NCOLS
            CALL UTRGTD(IDCCRC,COLIDS(I),1,ROW,GLOBAL(I),NULL(I),TSTAT)
            IF(TSTAT.NE.0)THEN
               WRITE(CONTXT,1009)COLNAM(I)
 1009          FORMAT('ERROR: Could read value from column ',a15,
     &              ' in CCRC table')
               CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
               ISTAT = TSTAT
               GO TO 999
            ENDIF
 2000    CONTINUE
C
C Get the thermal correction.
C
         CALL ZTHERM(IDCCRC,ROW,TNAMES,TCS,TEMPS,COR(1),ISTAT)
         IF(ISTAT.NE.0)THEN
            WRITE(CONTXT,1019)
 1019       FORMAT('No thermal correction will be applied')
            CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
            COR(1)=0.D0
         ENDIF
C
C Get the time correction.
C
         CALL ZTIME(JD,COR(2),ISTAT)
         IF(ISTAT.NE.0)THEN
            WRITE(CONTXT,1028)
 1028       FORMAT('No time correction will be applied')
            CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
            COR(1)=0.D0
         ENDIF
C
C Get the thermal linear correction.
C
         CALL ZTEMP(IDCCRC,ROW,TNAMES(3),TCS(3),TEMPS(3),
     &        COR(3),ISTAT)
         IF(ISTAT.NE.0)THEN
            WRITE(CONTXT,1029)
 1029       FORMAT('No thermal linear correction will be applied')
            CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
            COR(1)=0.D0
         ENDIF
C
C And of first pass initializations.
C
      ENDIF
C
C Compute global coefficients based on carrousel position.
C
      CPOS = CARPOS
      CPOS2 = CPOS * CPOS
      F0 = COR(1) + COR(2) + GLOBAL(2) + GLOBAL(3)*CPOS +
     &     GLOBAL(4) * CPOS2
      F1 = COR(3) + GLOBAL(5) + GLOBAL(6)*CPOS +
     &     GLOBAL(7)*CPOS2
      F2 = GLOBAL(8) + GLOBAL(9)*CPOS +GLOBAL(10)*CPOS2
      F3 = GLOBAL(11) + GLOBAL(12)*CPOS +GLOBAL(13)*CPOS2
      F4 = GLOBAL(14) + GLOBAL(15)*CPOS +GLOBAL(16)*CPOS2
      F5 = GLOBAL(17) + GLOBAL(18)*CPOS +GLOBAL(19)*CPOS2
C
C Compute the dispersion coefficients.
C
      MCEN = GLOBAL(1)
      WCEN = (CAP(1)/MCEN)*SIN((CAP(2)-CPOS)/10430.378D0)
      K = MCEN * WCEN
      K2 = K * K
      K3 = K * K2
      COEF(1) = F0 - F1*K + F2*K2 - F3*K3 - F4*WCEN - F5*MCEN
      COEF(2) = F1 - 2*F2*K + 3*F3*K2
      COEF(3) = F2 - 3*F3*K
      COEF(4) = F5
      COEF(5) = F4
      COEF(6) = 0.0D0
      COEF(7) = 0.0D0
      COEF(8) = F3
C
      ISTAT = 0
 999  RETURN
      END

max_stars_repo_name: iraf-community/stsdas


text:
!-------------------------------------------------------------------
! Fast Spectral Method
!-------------------------------------------------------------------

subroutine boltzmann_fft(Q, f, kn, phi, psi, phipsi, unum, vnum, wnum, m)

    !include 'fftw3.f'
    include '/opt/intel/oneapi/mkl/latest/include/fftw/fftw3.f'

    integer, intent(in) :: unum, vnum, wnum, m
    real(kind=8), intent(inout) :: Q(unum, vnum, wnum)
    real(kind=8), intent(in) :: f(unum, vnum, wnum)
    real(kind=8), intent(in) :: kn
    real(kind=8), intent(in) :: phi(unum, vnum, wnum, m*(m-1)), psi(unum, vnum, wnum, m*(m-1))
    real(kind=8), intent(in) :: phipsi(unum, vnum, wnum)
    complex(kind=8) :: f_spec(unum, vnum, wnum), fc(unum, vnum, wnum), fc1(unum, vnum, wnum), fc2(unum, vnum, wnum), &
                       f_temp(unum, vnum, wnum), fc11(unum, vnum, wnum), fc22(unum, vnum, wnum)
    integer(kind=8) :: plan_backward, plan_backward2, plan_forward, plan_conv_f1, plan_conv_f2
    integer :: i

    ! initialize DFT
    ! forward(-1)/backward(+1) Fourier transformation
    call dfftw_plan_dft_3d(plan_backward, unum, vnum, wnum, f_spec, f_spec, FFTW_BACKWARD, FFTW_MEASURE)
    call dfftw_plan_dft_3d(plan_backward2, unum, vnum, wnum, f_temp, f_temp, FFTW_BACKWARD, FFTW_MEASURE)
    call dfftw_plan_dft_3d(plan_forward, unum, vnum, wnum, f_spec, fc, FFTW_FORWARD, FFTW_MEASURE)
    call dfftw_plan_dft_3d(plan_conv_f1, unum, vnum, wnum, fc1, fc11, FFTW_FORWARD, FFTW_MEASURE)
    call dfftw_plan_dft_3d(plan_conv_f2, unum, vnum, wnum, fc2, fc22, FFTW_FORWARD, FFTW_MEASURE)

    f_spec = cmplx(f, 0.0d0)
    call dfftw_execute_dft(plan_backward, f_spec, f_spec) ! backward
    f_spec = f_spec / (unum * vnum * wnum)

    ! shift zero-frequency component to center of spectrum
    call fftshift3d(f_spec)

    ! gain term
    f_temp = cmplx(0.0)
    do i=1,M*(M-1)
        fc1 = cmplx(f_spec * phi(:, :, :, i))
        fc2 = cmplx(f_spec * psi(:, :, :, i))
        call dfftw_execute_dft(plan_conv_f1, fc1, fc11)
        call dfftw_execute_dft(plan_conv_f2, fc2, fc22)
        !convolution
        f_temp = f_temp + fc11*fc22
    enddo

    ! loss term
    fc1 = cmplx(f_spec * phipsi)
    fc2 = f_spec
    call dfftw_execute_dft(plan_conv_f1, fc1, fc11)
    call dfftw_execute_dft(plan_conv_f2, fc2, fc22)

    ! gain-loss
    f_temp = f_temp - fc11 * fc22
    Q = 4.0d0 * PI**2 / kn / M**2 * real(f_temp)

    call dfftw_destroy_plan(plan_backward)
    call dfftw_destroy_plan(plan_backward2)
    call dfftw_destroy_plan(plan_forward)
    call dfftw_destroy_plan(plan_conv_f1)
    call dfftw_destroy_plan(plan_conv_f2)

end

max_stars_repo_name: vavrines/Kifort

max_issues_repo_name: vavrines/Kifort

max_forks_repo_name: vavrines/Kifort


text:
On April 1^st^, otherwise known as wiki:WikiPedia:April_Fools_Day April Fools Day, the Wiki Developers Wiki Administrators decided to play a fairly nasty trick on the wiki users.  When you entered the site via http://daviswiki.org you were presented with http://daviswiki.org/cool_files/intro.html an extravagant flash intro movie.  The movie showcased a few images from the site, had the Wiki Logo twirl around a bit, and ended with the http://msn.com MSN butterfly being stamped onto the Wiki Logo Davis Wiki Logo.

Once the main page was brought up, there were a few noticable differences: http://google.com/adsense Ads were pasted above and to the side of every page, the Wiki Logo had the http://msn.com MSN butterfly next to it, search results had wiki:WikiPedia:eBay eBay results laced in, and there was a news item noting:
  We are pleased to announce that in our continuing effort to bring quality online content to Davisites, Davis Wiki has partnered with MSN.  In doing so, we hope to offer a more pleasurable experience to all of our readers and editors.  Rest assured that aside from a few unobtrusive, highly targeted ads that improve your browsing experience, we remain true to our commitment to provide this unique medium as a totally free collaborative online community.  Davis Wiki by MSN Its better with the butterfly.

 April fools featured page

On april fools day, a page called Manholes was made featured.  A play on pages that catalog locations (such as Payphone Project), this page was entirely made up in a matter of minutes.  While Users/PhilipNeustrom I realize that a cool Manholes page could be made, I thought it would be funny to pretend to have a page listing all the manhole locations in Davis along with what its like to crawl into a sewer.  If someone actually wants to document all the manholes, the http://www.cityofdavis.org/gis/library/ Davis GIS data has the manholes listed.

A manhole (also known as a maintenance hole) is an entry way from the street into the great underground depths.  While almost always unnoticed by the casual pedestrian or driver, these manholes lead to an array of treasures such as sewers, storm drains, telephone switches, and the eversoelusive gas pipes.

Until now, citizens of Davis had to stumble helplessly through the streets in search of manholes.  The City and various corporations attempt to limit information on this important topic.  In the interest of the public good, we press on with a critical analysis of Davis vast and often chaotic manhole and manhole cover system.

After reading this page you should feel comfortable:
 1. Spotting manholes as you walk by them.
 2. Identifying different kinds of manholes by their cover design.
 3. Cataloging newly identified manhole covers for future scientific and public inquiry.
 4. Finding manholes with the use of GIS data, GPS receivers, and the Wiki Map functionality.
 5. Opening manhole covers and exploring sewers.

 Manhole identification

Manholes can usually be spotted by means of finding their manhole cover  a (usually round) region cut out of the normal sidewalk or street, typically with a few holes through it.  The cover usually reads Sewer, Water, or Telephone.  However, the only real way to know whats inside of the manhole is to open it yourself (covered below!)

 What to do if you find a manhole
In the rare event you encounter a manhole, what is the proper way to report it?  Here is the recommended approach, used by __Citizens for Open Manholes__:

 Step 1

 Step 2

 Step 3

 Step 4

 Step 5

 Opening manholes

While taking an interest in the manhole cover itself may often reveal important details essential to urban survival, there is often the need to adjust things like telephone service, sewer systems, and water pipes in a more prompt and effective manner (after all, when was the last time city services got this right on their own?)  Unfortunately, like so many other things, matters are best taken into your own, skilled hands.

The first thing youll need to do when opening a manhole is identify the holes in the cover.  These are typically circular and toward the edges, but they can be found in different shapes and locations.  The manhole at Third and B (pictured left) features 4 holes, each on the far edge of the disc.   In order to open the manhole it was required to use a pair of crowbars.  In this case, two people worked together to open the manhole, but in the case youre going at it alone youll need to either get both crowbars in independently then push on both at once (somewhat difficult to do but requires no prep work) or weld two crowbars together at the base and press a single bar down (easy to do but requires metal shop prep work).

 Entering the manhole

An entire world opens up before your eyes once the manholes cover has been lifted  but beware!  Many sewers are infested with rats (which carry the plague) and possibly alligators (please update if you can confirm alligators in Davis sewers).  There may or may not be a ladder to ease entry into the sewer system.  You should be prepared with a headmounted light, as well.  After youre inside, make sure you dont let the manhole cover slip back over the hole!  It can be very, very difficult to remove a manhole cover from the inside (trust me on this one.)

 Finding manholes

Because manholes are so rare, we attempt to compile a map of local manholes.  Through the amazing collaborative powers of wiki, we can maintain a record of all manholes, information about each one, some designs (see What to do if I find a manhole), photographs, and personal stories/accounts of various manholes.

  {{{Downtown Manhole 1}}}
  {{{Downtown Manhole 2}}}
  {{{Downtown Manhole 3}}} (a favorite among many)
  {{{Downtown Manhole 4}}}
  {{{Downtown Manhole 5}}}
  {{{Downtown Manhole 6}}}

 Some reactions
20050401 08:38:02 nbsp Hmmm.... I hate the ads. I hate MSN probably about as much as I hate ads. I realize the need for funding and that some of us yes, I mean myself have been extremely lame in that department. But of all the things Ive chosen to live with, inviting corporate into the community is not something Im sure I can handle. Im actually kinda surprised at how much this disturbs me. Maybe it would be helpful to know what benefits msn, its ads, and data collection provide... shouldnt we at least get some sort of privacy notification now that weve been sold to corporate? Users/AlphaDog

20050401 08:41:23 nbsp Users/AlphaDog, look at your calendar. Users/KenBloom

20050401 12:40:23 nbsp Personally, I have the https://addons.mozilla.org/extensions/moreinfo.php?applicationfirefox&version1.0&categoryHumor&numpg10&id217 Adbar on my version of Firefox.  I just dont feel comfortable browsing places like daviswiki without the familiar wiki:WikiPedia:Fnord fnords. Users/BrentLaabs

20050401 13:03:03 nbsp I dont like the new intro at all. It took too long to load on high speed and its excessive. Id rather have it the way it was. SS

20050401 13:37:38 nbsp WTF. Better with the butterfly? It makes it seem like the wiki is being coded by people from microsoft? Can someone please explain a little bit better? Users/MattJojola

20050401 14:09:11 nbsp Chill out. Please everyone check your calendar, and wait a couple days before complaining about all of the new features. Users/KenBloom
  When I saw your first response, I just figured you were being, uh, smart. Now I realize that Im not being so smart... whats this calendar that it sounds like I should know something about?? Maybe we could have some sort of discussion about financing the wiki. I really enjoy it to the point that I spend way too much time here, and definitely appreciate all the time and money that you guys have poured in to create this community. Im curious to know if this might be something that would be eligible for community funding from the City of Davis  that might at least offset the direct costs of hosting the wiki. BTW: Shouldnt that msn butterflies & moths flutterby be linked to one of our internal pages?? Users/AlphaDog

20050401 15:33:49 nbsp http://daviswiki.cheshirehall.org/  the new AD FREE WIKI!!! Users/JabberWokky

20050401 16:04:22 nbsp 1) look what day it is 2) i cant see the ads, get the adblock extension for firefox Users/ArlenAbraham

20050401 21:23:36 nbsp This is some funny shit, man! :) Users/SummerSong

20050401 21:58:26 nbsp Crap... its payback time. Users/AlphaDog


text:
      program t
      print '(f9.2)',1234e+2
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


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

! Check for semantic errors in ALLOCATE statements

subroutine C935(l, ac1, ac2, ac3, dc1, dc2, ec1, ec2, aa, ab, ea, eb, da, db, whatever, something)
! A type-param-value in a type-spec shall be an asterisk if and only if each
! allocate-object is a dummy argument for which the corresponding type parameter
! is assumed.

  type A(la)
    integer, len :: la
    integer vector(la)
  end type

  type, extends(A) :: B(lb)
    integer, len :: lb
    integer matrix(lb, lb)
  end type

  type, extends(B) :: C(lc1, lc2, lc3)
    integer, len :: lc1, lc2, lc3
    integer array(lc1, lc2, lc3)
  end type

  integer l
  character(len=*), pointer :: ac1, ac2(:)
  character*(*), allocatable :: ac3(:)
  character*(:), allocatable :: dc1
  character(len=:), pointer :: dc2(:)
  character(len=l), pointer :: ec1
  character*5, allocatable :: ec2(:)

  class(A(*)), pointer :: aa
  type(B(* , 5)), allocatable :: ab(:)
  type(B(* , *)), pointer :: ab2(:)
  class(A(l)), allocatable :: ea
  type(B(5 , 5)), pointer :: eb(:)
  class(A(:)), allocatable :: da
  type(B(: , 5)), pointer :: db(:)
  class(*), allocatable :: whatever
  type(C(la=*, lb=:, lc1=*, lc2=5, lc3=*)), pointer :: something(:)
  type(C(la=*, lb=:, lc1=5, lc2=5, lc3=*)), pointer :: something_else(:)

  ! OK
  allocate(character(len=*):: ac1, ac3(3))
  allocate(character*(*):: ac2(5))
  allocate(B(*, 5):: aa, ab(2)) !OK but segfault GCC
  allocate(B(*, *):: ab2(2))
  allocate(C(la=*, lb=10, lc1=*, lc2=5, lc3=*):: something(5))
  allocate(C(la=*, lb=10, lc1=2, lc2=5, lc3=3):: aa)
  allocate(character(5):: whatever)

  ! Not OK

  ! Should be * or no type-spec
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=5):: ac1)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=5):: ac2(3), ac3)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=l):: ac1)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=l):: ac2(3), ac3)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(A(5):: aa)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(B(5, 5):: ab(5))
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(B(l, 5):: aa, ab(5))

  ! Must not be *
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=*):: ac1, dc1, ac3(2))
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character*(*):: dc2(5))
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character*(*):: ec1)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(*):: whatever)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(character(len=*):: ac2(5), ec2(5))
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(A(*):: ea) !segfault gfortran
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(B(*, 5):: eb(2)) !segfault gfortran
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(A(*):: da) !segfault gfortran
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(B(*, 5):: db(2)) !segfault gfortran
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(A(*):: aa, whatever)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(B(*, *):: aa)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(C(la=*, lb=10, lc1=*, lc2=5, lc3=*):: something_else(5))
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(C(la=5, lb=10, lc1=4, lc2=5, lc3=3):: aa)
  !ERROR: Type parameters in type-spec must be assumed if and only if they are assumed for allocatable object in ALLOCATE
  allocate(C(la=*, lb=10, lc1=*, lc2=5, lc3=*):: aa)
end subroutine


text:
Davis is adorned with numerous parks, most of which are networked by The Greenbelt and thus accessible by Bicycles bicycle. In Davis parks, youll find Frisbee golf frisbee golf, Horseshoes horseshoe pits, archery, a fishing pond, Soccer soccer/Baseball baseball/Softball softball fields, Tennis Courts tennis courts, an area for R/C gliders, picnic BBQ grills and picnic tables, Basketball basketball courts, batting cages, Dogs dog areas, Town Art town art, volleyball courts, a rock climbing wall, exercise courses, a fire pit, swimming pools, Roller Hockey roller and Street Hockey street hockey rinks, amphitheatres, a Public Skate Park skate park, toddler play areas, a The Flying Carousel of the Delta Breeze carousel, fountains, Drinking Fountains drinking fountains, wildlife, and all of the usual amenities.

Map all the parks by going to Davis Map (select only Parks and Recreation) or visit the official http://cityofdavis.org/pgs/facility/ parks page provided by the City of Davis.
An http://maps.google.com/maps/ms?msid201144637277324396936.0004b60bf93a20421320a&msa0&ll38.553938,121.762934&spn0.080143,0.06094 editable google map of parks and playgrounds in davis is under construction.  Eventually all locations should link to park pages on daviswiki.

North Davis Parks
  Covell Park  tennis courts, amphitheatre
  Flicker Pocket Park
  Hacienda Park
  Northstar Park  fishing pond, soccer field
  Northstar Pocket Park Northstar Pocket Park aka. Teacup Park  small play area, bbq

South Davis Parks
  East Putah Creek Park
  Pioneer Park  tennis/basketball courts, soccer fields, toddler play area, volleyball posts, offleash dog area
  Playfields Park  baseball/softball/soccer fields, basketball/volleyball courts
  Putah Creek Park
  Putah Creek Riparian Reserve
  Putah Creek South Fork Preserve
  Rosecreek Greenbelt Play Area
  Southfield Greenbelt Play Area
  Village Park  basketball court
  Walnut Park  offleash dog area, volleyball posts, toddler play area
  Willowcreek Park  basketball court, backstop
  Woodbridge Park  picnic area, toddler play area
  Yolo County Grasslands Regional Park  wildlife, RC glider area, archery, horseshoes, offleash dog area

East Davis Parks
  Cedar Park  Horseshoes horseshoe pit
  Chestnut Park  tennis/basketball courts, softball field
  Frank Lloyd Wright Park   field
  John Barovetto Park  Baseball, Soccer Field, Basketball, swings and play area
  La Playa Park  basketball courts, backstop
  Mace Ranch Park  basketball courts
  Manor Park or Slide Hill Park  tennis/basketball courts, toddler play area, backstop, pool, huge slide
  Oak Grove Park  basketball courts
  Pirate Park
  Robert Arneson Park
  Sandy Motley Park  soccer fields, backstop
  Toad Hollow Dog Park  for dogs
  Wildhorse Park

West Davis Parks
  Arroyo Park  rock climbing wall, fire pit, pool
  Aspen Greenbelt  wildlife, offleash dog area
  Dinosaur Park  toddler play area, swings, picnic benches, sand
  Stonegate Dog Park  for dogs
  West Manor Park  softball field, tennis/basketball courts, roller/street hockey rink
  Westwood Park  soccer/softball fields
  Whaleback Park

Central Davis Parks
  Civic Center Park
  College Park
  Community Park  offleash dog area, skate park, pool
  N Street Park  volleyball posts
  Oxford Circle Park (Little Big Peoples Park)  amphitheatre, frisbee golf course, Horseshoes horseshoe pit
  Redwood Park
  Sycamore Park  tennis courts, offleash dog area

Downtown Parks
  Central Park  basketball court, carousel, horseshoe pits, fountain

Quick Features/Shade Guide
Many people wonder which parks are shaded or have certain features. Use this table to help compare shade and basic features of the parks.

Legend
Morning shade checked around 10:00am; Afternoon shade checked around 3:30pm

#Playgrounds/Play_Areas Playgrounds/Play Areas
Many parks include play structures, slides, swings, and sandboxes.  There is a http://cityofdavis.org/media/default/documents/pdf/gis/maps/park%20play%20areas.pdf pdf map on the city website that categorizes play areas based on their suitability for tots (ages 25) and/or older kids (612).


!***********************************/
!	Name:影響域分離法のための配列確保
!	Alias:InitializeInfluenceSeparation
!	Description:
!	Type:
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.01.31
!	Update:
!	Other:
!***********************************/
    subroutine InitializeInfluenceSeparation(UG,MG,CC,CC4MB)
    use StructVar_Mod
    implicit none
    type(UnstructuredGrid), intent(in) :: UG
    type(MoveGrid), intent(inout) :: MG
    type(CellCenter), intent(in) :: CC
    type(CellCenter), intent(inout) :: CC4MB

    integer :: iCell
    CC4MB = CC
    CC4MB%InterpolatedQuantity = CC%ConservedQuantity

    do iCell = 1, UG%GI%RealCells
        MG%IS(iCell)%InfluenceDepth = 0
        MG%IS(iCell)%BoundaryAdjacent = 0
    end do

    MG%MaxInfluenceRange = 0.0d0

    return
    end subroutine InitializeInfluenceSeparation


text:
      SUBROUTINE PLA31
C
C     THIS ROUTINE READS THE INCREMENTAL DISPLACEMENT VECTOR INTO CORE
C     AND APPENDS THE PROPER DISPLACEMENT VALUES TO THE ESTNL ENTRY FOR
C     EACH ELEMENT, THEREBY CREATING THE ESTNLS, THE ESTNL SCRATCH FILE,
C     WHICH IS PROCESSED BY SUBROUTINE PLA32.
C
      INTEGER         BUFSZ,BUFR1,BUFR2,DELUGV,ESTNL,ESTNLS,FILE,EOR,
     1                CLSRW,IZ(1),IESTBK(100),ESTWDS(40),ELTYPE
      DIMENSION       NAME(2),NGPTS(40),MCBUGV(7),ESTBK(100)
      COMMON /BLANK / ICOM
      COMMON /SYSTEM/ BUFSZ
      COMMON /ZZZZZZ/ Z(1)
      COMMON /UNPAKX/ ITYPEB,IUNPK,JUNPK,INCUPK
      EQUIVALENCE     (Z(1),IZ(1)),(ESTBK(1),IESTBK(1))
      DATA    NAME  / 4HPLA3,4H1   /
      DATA    DELUGV, ESTNL,ESTNLS / 104,105,301/
      DATA    EOR   , NEOR,CLSRW   / 1,0,1      /
C
C    1        ROD       BEAM      TUBE      SHEAR     TWIST
C    2        TRIA1     TRBSC     TRPLT     TRMEM     CONROD
C    3        ELAS1     ELAS2     ELAS3     ELAS4     QDPLT
C    4        QDMEM     TRIA2     QUAD2     QUAD1     DAMP1
C    5        DAMP2     DAMP3     DAMP4     VISC      MASS1
C    6        MASS2     MASS3     MASS4     CONM1     CONM2
C    7        PLOTEL    REACT     QUAD3     BAR       CONE
C    8        TRIARG    TRAPRG    TORDRG    CORE      CAP
C
      DATA    ESTWDS/
     1           21,        0,       20,        0,        0,
     2           38,        0,        0,       27,       21,
     3            0,        0,        0,        0,        0,
     4           32,       32,       37,       43,        0,
     5            0,        0,        0,        0,        0,
     6            0,        0,        0,        0,        0,
     7            0,        0,        0,       50,        0,
     8            0,        0,        0,        0,        0 /
      DATA    NGPTS /
     1            2,        2,        2,        4,        4,
     2            3,        3,        3,        3,        2,
     3            2,        2,        2,        2,        4,
     4            4,        3,        4,        4,        2,
     5            2,        2,        2,        2,        2,
     6            2,        2,        2,        2,        2,
     7            2,        0,        0,        2,        2,
     8            3,        4,        2,        4,        2 /
C
C     DETERMINE SIZE OF CORE, DEFINE BUFFERS AND INITIALIZE CORE
C     POINTERS AND COUNTERS
C
      IZMAX = KORSZ (Z)
      BUFR1 = IZMAX - BUFSZ
      BUFR2 = BUFR1 - BUFSZ
      LEFT  = BUFR2 - 1
      IDISP = 0
C
C     OPEN THE DISPLACEMENT VECTOR FILE AND READ THE DISPLACEMENT VECTOR
C     INTO OPEN CORE.
C
      FILE = DELUGV
      CALL GOPEN (DELUGV,Z(BUFR1),0)
      MCBUGV(1) = DELUGV
      CALL RDTRL (MCBUGV)
      IF (LEFT .LT. MCBUGV(3)) CALL MESAGE (-8,0,NAME)
      ITYPEB = 1
      IUNPK  = 1
      JUNPK  = MCBUGV(3)
      INCUPK = 1
      CALL UNPACK (*9040,DELUGV,Z(IDISP+1))
      CALL CLOSE (DELUGV,CLSRW)
C
C     BUILD THE SCRATCH FILE ESTNLS
C
      CALL GOPEN (ESTNL,Z(BUFR1),0)
      CALL GOPEN (ESTNLS,Z(BUFR2),1)
C
C     READ AN ELEMENT TYPE FROM ESTNL AND WRITE IT ON ESTNLS.
C
   10 CALL READ (*60,*9030,ESTNL,ELTYPE,1,NEOR,IFLAG)
      NWDSRD = ESTWDS(ELTYPE)
      IF (NWDSRD .LE. 0) CALL MESAGE (-30,91,ELTYPE)
      CALL WRITE (ESTNLS,ELTYPE,1,NEOR)
C
C     READ AN ESTNL ENTRY
C
   20 J = NWDSRD
      CALL READ (*9020,*50,ESTNL,ESTBK,J,NEOR,IFLAG)
      NOGPTS = NGPTS(ELTYPE)
      IF (NOGPTS .LE. 0) CALL MESAGE (-30,92,ELTYPE)
C
C     APPEND THE DISPLACEMENT VECTORS ONTO THE ESTBK.
C
      NWDS = 3
      J = J + 1
      IF (ELTYPE.EQ. 1 .OR. ELTYPE.EQ. 3 .OR. ELTYPE.EQ.10 .OR.
     1    ELTYPE.EQ. 6 .OR. ELTYPE.EQ.17 .OR. ELTYPE.EQ.18 .OR.
     2    ELTYPE.EQ.19 .OR. ELTYPE.EQ.34) NWDS = 6
      DO 40 I = 1,NOGPTS
      INDEX = IDISP + IESTBK(I+1)
      DO 30 K = 1,NWDS
      ESTBK(J) = Z(INDEX)
      INDEX = INDEX + 1
   30 J = J + 1
   40 CONTINUE
C
C     THE APPENDED ESTNL ENTRY, WHICH IS AT ESTBK IS NOW COMPLETE.
C
      CALL WRITE (ESTNLS,ESTBK,J-1,NEOR)
      GO TO 20
C
C     WRITE AN EOR ON THE ESTNLS FILE.
C
   50 CALL WRITE (ESTNLS,0,0,EOR)
      GO TO 10
C
C     PROCESSING IS NOW COMPLETE
C
   60 CALL CLOSE (ESTNL,CLSRW)
      CALL CLOSE (ESTNLS,CLSRW)
      RETURN
C
C     FATAL ERRORS
C
 9020 CALL MESAGE (-2,FILE,NAME)
 9030 CALL MESAGE (-3,FILE,NAME)
 9040 CALL MESAGE (-5,DELUGV,NAME)
      RETURN
      END


text:
Abortion is the termination of a pregnancy, whether voluntary or involuntary. Some view it as terminating a life while others do not. It is a choice that many may face.  If you believe yourself to be pregnant and would like further information or to know the resources available, you can go to the WRRC, Student Health and Wellness Center, Sutter Davis Hospital or Planned Parenthood. Ways of prevention include abstinence and Contraception (before, during, and after Sex).  A full explanation and other alternatives are listed on the Pregnancy page.

In 1989, the Davis City Council approved a resolution to declare Davis a prochoice city.  The LA Times http://articles.latimes.com/19891118/news/mn1577_1_citycouncil ran a short piece about it.

Oftentimes the cost of having an abortion increases with the amount of weeks one is pregnant.

While abortion has stirred up a few protests in Davis, the Davis Wiki is not a good place to host debates.

 Getting an Abortion in Davis

If you have decided to seek an abortion near Davis, you will have the following choices to make:
 If you are in the first trimester, do you want to have a medical (by pills) or surgical (vacuum aspiration) abortion? Both types of services are available locally, but some providers only offer one or the other.
 If you opt for or must have a surgical abortion, is it important to you to be able to have someone in the exam room with you to provide support? Some allow this (e.g. Sutter Hospitals Womens Clinic and The Womens Health Specialists of Sacramento); Planned Parenthood does not.
 Will you pay for this procedure with insurance or not? Even if you have health insurance, paying out of pocket can sometimes be cheaper (e.g. at Planned Parenthood and The Womens Health Specialists of Sacramento) depending on your insurance policy and income, due to slidingscale cost structures.

Information you will want to have on hand:
 The date of the first day of your last period
 Blood type (women with some blood types will be given antibody injections; if you dont know it they will give you a blood test)
 Your income (some providers charge on a sliding scale)

Types of abortion in the first trimester: medical (by pills) or surgical (vacuum aspiration)
If you are up to seven weeks pregnant (nine weeks from the start of your last period) you may opt for a medical abortion, in which case, the doctor will give you a pill at the office following an exam, and you will take the second pill (actually four pills of one medication) two days later at home. In the case of a medical abortion, the abortion process will last about seven days with the heaviest cramping and bleeding during the six hours following pill ingestion. You may be prescribed pain medication and/or be advised to take over the counter pain medication. A week later, you will follow up at the clinic to verify that the abortion was successful. 5–8% of women will discover the abortion failed, and will need to then have a surgical abortion.

Surgical abortions cannot be performed as early on as can medical abortions but they can be performed later on than medical abortions. They are generally are not performed on women who are less than 6 weeks from the first day of their last period. Surgical abortions are quickereverything from your initial consultation to the abortion itself can be completed in one day. Two procedures are usually performed in this day: in the morning the patient receives counseling, an ultrasound is performed to verify the stage and location of the pregnancy, and the patient receives treatment to dilate and soften her cervix; in the afternoon the patient undergoes the abortion procedure. A variety of pain medications may be prescribed, and can be accompanied by local anesthesia (a cervical block). While the procedure itself is painful, subsequent cramping is less prolonged than in a medical abortion and bleeding is lighter.

After the first trimester, only surgical abortions are performed, and procedures other than vacuum aspiration are employed.

Planned Parenthood
You can http://www.plannedparenthood.org/findCenterProcess.asp search for planned parenthood locations near Davis. You can also call the national Planned Parenthood number for more information: 1800230PLAN. The various local Planned Parenthood clinics offer abortion services on different days of the week. As of Aug 2008, the nearest Planned Parenthood that provides abortion services is in Sacramento.  Its location is:

Some Planned Parenthood clinics offer both medical (by pills) and surgical (vacuum aspiration) abortions; others provide only one of these types of services. For example, Fairfield Planned Parenthood only provides surgical abortions. In the case of a surgical abortion at a Planned Parenthood clinic, you can expect to arrive in the morning (between 7am and 10am depending on your appointment time), receive initial counseling and medical preparation (cervical dilation), and wait onsite for 15 hours until your cervix is prepared and the practitioner is available (the wait may be hours longer if it is a busy day) for the abortion procedure. The abortion procedure itself will only take a few minutes. Planned Parenthoods policy is not to allow anyone (e.g. a friend or partner) to accompany the patient in the exam room during the procedure, citing patient safety, although they can accompany you in the waiting room. Other providers do allow the patient to be accompanied.

 Sutter Davis Hospital

The Womens Clinic at Sutter Davis Hospital is the most local option in Davis aside from private practices. Both medical (by pills) and surgical (vacuum aspiration) abortions are available here. Patients undergoing a surgical abortion are not expected to remain at the clinic between their morning and afternoon appointments (like at Planned Parenthood), and can go home during this interval. The base price of either procedure for someone without insurance is $450, but the negotiated rate charged to insurance companies, of which you will pay a percentage, is higher (likely about $600 for someone with the Aetna Student health Plan, but it will vary). Your insurance may cover partial or full cost.

The clinic is in the building immediately to the west of the hospital which is about five steps away the Unitrans P/Q Bus Line Pline stop (location #44 on the route map; stop #199). The Qline stop is across the street.

 The Feminist Womens Health Center / Womens Health Specialists of Sacramento
http://www.womenshealthspecialists.org/
Like Planned Parenthood, this organization provides abortion services on an incomedependent sliding scale to those without insurance. Cost appears to be higher than Planned Parenthood, but unlike at Planned Parenthood patients opting for a surgical abortion may be accompanied by a friend in the exam room during the abortion procedure. As of January, 2014, surgical abortions are performed on Wednesdays.

 Information for UC Davis Students

Availability: UC Davis Student Health and Counseling Services (SHCS) does not provide abortion services. Instead, a primary care provider will provide you with a referral to seek these services nearby. If you are not going to use your student health plan insurance, you do not need to get a referral from SHCS before contacting abortion providers to schedule your appointment.

Cost: Most students are on the Aetna Student Health Plan. Depending on what provider of abortion services you choose, it might be cheaper for your to use your Aetna Student Health Plan insurance or not. Under this insurance, you have a $300 deductible, and then your insurance will cover either 80% or 60% of the remaining negotiated cost of the procedure, depending on whether you choose a provider that is innetwork or outofnetwork. As of January, 2014, the only local innetwork provider of abortion services is the Sutter Womens Clinic. While Planned Parenthood is not an innetwork provider, they charge for abortion services on a sliding scale to those not using insurance, so depending on your income the cost to you might be less if you do not use your insurance. The Insurance Services desk can give you more information.

 Groups opposing Abortion Rights

  Students for Life at UCD
  Christian Democratic Party
  http://www.saclife.org/ The Sacramento Life Center provides free pregnancy tests, ultrasounds, counseling
  http://www.bgmh.org/ Bishop Gallegos Maternity Home  gives a loving, safe place for women to live and have their child.
  http://www.rachelsvineyard.org/ Rachels Vineyard offers help to those healing from abortion.

 Groups supporting Abortion Rights

  National Organization for Women
  Planned Parenthood provides abortions, as well as other prenatal options. To find the nearest location: http://www.plannedparenthood.org/findCenterProcess.asp Search
  Womens Health Specialists/ The Feminist Womens Health Centers of California http://www.womenshealthspecialists.org/