text:
      SUBROUTINE TRIO (J)
C     COMPUTE THE COMPONENTS OF ALL BASIS FUNCTIONS ON SEGMENT J
      COMMON /DATA/ LD,N1,N2,N,NP,M1,M2,M,MP,X(300),Y(300),Z(300),SI(300
     1),BI(300),ALP(300),BET(300),ICON1(300),ICON2(300),ITAG(300),ICONX(
     2300),WLAM,IPSYM
      COMMON /SEGJ/ AX(30),BX(30),CX(30),JCO(30),JSNO,ISCON(50),NSCON,IP
     1CON(10),NPCON
      DATA JMAX/30/
      JSNO=0
      JCOX=ICON1(J)
      IF (JCOX.GT.10000) GO TO 7
      JEND=-1
      IEND=-1
      IF (JCOX) 1,7,2
1     JCOX=-JCOX
      GO TO 3
2     JEND=-JEND
3     IF (JCOX.EQ.J) GO TO 6
      JSNO=JSNO+1
      IF (JSNO.GE.JMAX) GO TO 9
      CALL SBF (JCOX,J,AX(JSNO),BX(JSNO),CX(JSNO))
      JCO(JSNO)=JCOX
      IF (JEND.EQ.1) GO TO 4
      JCOX=ICON1(JCOX)
      GO TO 5
4     JCOX=ICON2(JCOX)
5     IF (JCOX) 1,9,2
6     IF (IEND.EQ.1) GO TO 8
7     JCOX=ICON2(J)
      IF (JCOX.GT.10000) GO TO 8
      JEND=1
      IEND=1
      IF (JCOX) 1,8,2
8     JSNO=JSNO+1
      CALL SBF (J,J,AX(JSNO),BX(JSNO),CX(JSNO))
      JCO(JSNO)=J
      RETURN
9     PRINT 10, J
      STOP
C
10    FORMAT (44H TRIO - SEGMENT CONNENTION ERROR FOR SEGMENT,I5)
      END

ext: f

max_stars_repo_name: ViperScience/ViperNEC

max_issues_repo_name: ViperScience/ViperNEC

max_forks_repo_name: ViperScience/ViperNEC


text:
Address(Monarch Lane) is a residential street in East Davis that runs from Loyola Drive up to East Covell Boulevard.

  Loyola Drive
  Two connections with Layton Drive
  Two connections with Grinnell Drive and one with Grinnell Place
  Temple Drive
  Campbell Place
  Bryant Place
  East Covell Boulevard

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module struct_fact
!! Computes the struture factor from the radial distribution functions.
use iso_fortran_env
use iso_c_binding
use sinft
implicit none

real(8), parameter, private :: PI=3.141592653589793

private :: eleFFS, do_dst

contains

function EleFFS(q, r)
  !! Proportional to the form factor of a sphere.
  real(8) :: EleFFS
  real(8), intent(in) :: q,r
  real(8) :: x
  x=q*r
  eleFFs = sin(x)-x*cos(x)
end function

subroutine do_dst(n, y, out)
  !! Do a sinus tranform of y-1.
  integer :: n !! dat length
  real(8), intent(in) :: y(n) !! input data
  real(8) :: out(n) !! result
  real(8) :: in(n)
  integer :: i
  do i=1,n
    in(i) = (y(i)-1.D0)*(real(i,8)-0.5D0)
  end do
  call sinft_calc(n, in, out)
end subroutine

subroutine calc_sdeq( dx, n,  q,  sf,  pf, nqext,qext,  sfext, ng,  gder, nf,  pop&
	      &,  rad,  phi)
  !! Calculates S(q) from RDFs.
  use spline
  integer, intent(in) :: n !! length of q array
  integer, intent(in) :: nqext !! length of qext
  integer, intent(in) :: ng  !! lenght of RDFs
  integer, intent(in) :: nf  !! number of families
  real(8), intent(in) :: dx  !! space bin of RDFs
  real(8), intent(in) :: phi !! volume fraction
  real(8) :: q(n) !! wave numbers
  real(8) :: sf(n) !! struture factor
  real(8) :: pf(n) !! form factor
  real(8) :: qext(nqext) !! new q values
  real(8) :: sfext(nqext) !! interpolated struture factor
  real(8) :: gder(n,ng)  !! RDFs
  real(8) :: pop(nf)     !! family populations
  real(8) :: rad(nf)     !! family radii

  integer :: i, j,k,ind
  real(8) :: vp, gsk, dump,dsq
  real(8), allocatable :: s(:),g(:),rho(:),sf2ip(:)
  allocate(s(n),g(n),rho(nf))

  if(ng /= (nf*(nf+1))/2) stop "calc_sqdeq: ng != (nf*(nf+1))/2"

  vp=0.0
  do i=1,nf
    vp = vp + pop(i)*4.D0/3.D0*PI*(rad(i)**3)
  end do
  do i=1,nf
    rho(i) = pop(i)*phi/vp
  end do

  do k=1,n
    q(k) = (PI*(k))/(dx*n);
  end do

  sf = 0.D0
  pf = 0.D0

  ind = 0
  do i=1,nf
    ind = ind + 1
    do k=1,n
      g(k) = gder(k,ind)
    end do
    call do_dst(n, g, s)
    do k=1,n
      gsk = 1.d0 + (2.d0*PI*(rho(i)*dx*dx)/q(k))*s(k)
      dump = EleFFS(q(k), rad(i))*EleFFS(q(k), rad(i))*rho(i)
      sf(k) = sf(k) + gsk*dump
      pf(k) = pf(k) + dump
    end do
  end do
  do i=1,nf
    do j=i+1,nf
      ind = ind + 1
      do k=1,n
        g(k) = gder(k,ind)
      end do
      call do_dst(n, g, s)
      do k=1,n
        dsq = sqrt(rho(i)*rho(j))
        gsk = 2.D0*PI*(dsq*dx*dx)/q(k)*s(k)
        dump = EleFFS(q(k), rad(i))*EleFFS(q(k), rad(j))*dsq;
        sf(k) = sf(k) + 2.D0*gsk*dump
      end do
    end do
  end do
  call sinft_cleanup()
  sf = sf / pf

  allocate(sf2ip(n))
  call cubspline_init(n, q, sf,sf2ip)
  call cubspline_equalsep_eval(n,q,sf,sf2ip,nqext,qext, sfext)

end subroutine

end module struct_fact


text:
http://www.certifiedfoods.com/jbflours.htm Organic, Kosher bulk flours and http://www.certifiedfoods.com/brewers.htm brewers malts.

20120731 14:51:51 nbsp According to their website (and my friends), they sell in pallet lots, truck loads or rail cars.  It was unusual, however, my friends were able to buy their flour by the individual (albeit large) sack.   Users/Angel.York


text:
      Subroutine gespd(a,rhs,sol,n,m)
!
!  Solves Ax = d with A a nxn SPD and d a nxm.
!
      implicit none
      real*8, dimension(n,n), intent(inout):: a
      real*8, dimension(n,m), intent(inout):: rhs
      real*8, dimension(n,n+m):: aa
      real*8, dimension(n,m) :: y
      real*8, dimension(n,m),intent(out)::sol
      integer ::k,i,j,l
      integer,intent(in)::n,m
      aa(1:n,1:n)= a
      aa(1:n,(n+1):n+m) = rhs
!
!  Factor A via column version and
!  write over the matrix.
!
      do k=1,n-1
        aa(k+1:n,k) = aa(k+1:n,k)/aa(k,k)
        do j=k+1,n
          do i=k+1,n
            aa(i,j) = aa(i,j) - aa(i,k)*aa(k,j)
          enddo
        enddo
      enddo
!
!  Solve Ly = d via column version and
!  multiple right sides.
!
      do j=1,n-1
        do l =1,m
          y(j,l)=aa(j,n+l)
        enddo
        do i = j+1,n
          do l=1,m
            aa(i,n+l) = aa(i,n+l) - aa(i,j)*y(j,l)
          enddo
        enddo
      enddo
!
!  Solve Ux = y via column version and
!  multiple right sides.
!
      do j=n,2,-1
        do l = 1,m
          sol(j,l) = aa(j,n+l)/aa(j,j)
        enddo
        do i = 1,j-1
          do l=1,m
            aa(i,n+l)=aa(i,n+l)-aa(i,j)*sol(j,l)
          enddo
        enddo
      enddo
      do l=1,m
        sol(1,l) = aa(1,n+l)/a(1,1)
      enddo
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       Written by: Anshuk Attri
C
C       Contact: contact@anshukattri.in
C       Website: www.anshukattri.in/research
C       GITHUB: github.com/aa-tree/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        SUBROUTINE TRANS_ORIGIN_FROMSSB(ET, TOFRAME, VIN, VOUT)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        IMPLICIT NONE
        INCLUDE 'common.dat'

        REAL*8 VIN(6), VOUT(6), ET
        INTEGER FROMF, TOFRAME

        REAL*8 R0(6), RIBC(6)
        REAL*8 LTIME
        INTEGER I

C
C       Intitialising output vector with zereos.
C
        DO I=1,6
          VOUT(I)=0.00D0
        END DO

C
C       Getting coordinates of the frame, with respect to Solar System
C       Barycentre.
C
        CALL SPKEZ(TOFRAME, ET, 'J2000', 'NONE', 0,
     &RIBC, LTIME)

        DO I=1,6
          VOUT(I)=VIN(I)-RIBC(I)
        END DO

        RETURN
        END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

max_stars_repo_name: aa-tree/SPICEReader

max_issues_repo_name: aa-tree/SPICEReader

max_forks_repo_name: aa-tree/SPICEReader


text: Yolo County Buyers Guide is a professional media publication created by Davis Life Magazine for Coldwell Banker Doug Arnold Real Estate, Inc. The magazine features stories about Davis, the real estate market, homes for sale, real estate agents, preferred vendors, etc. On the first of each month a new magazine is published. Visit the real estate sections to see listings, virtual tours, agent directories and open houses. The Yolo County Buyers Guide is another great media resource for Davis. You may sign up for a subscription at no cost by visiting http://www.yolocountybuyersguide.com their website.


text:
PROGRAM quad_surf
! Fit a quadratic surface to a set of data

USE lsq
IMPLICIT NONE
INTEGER              :: iostatus, n, ier, row, col
INTEGER, PARAMETER   :: nmax = 100
REAL (dp), PARAMETER :: wt = 1.0_dp
REAL (dp)            :: rh, temp, yval, x(nmax,5), y(nmax), xrow(0:5), b(0:5), &
                        yfit(3)

OPEN(UNIT=8, FILE='dry_bulb.dat', STATUS='OLD')
n = 0
CALL startup(5, .TRUE.)
DO
  READ(8, *, IOSTAT=iostatus) rh, temp, yval
  IF (iostatus < 0) EXIT
  IF (iostatus > 0) CYCLE
  n = n + 1
  rh = rh - 70._dp
  temp = temp - 95._dp
  y(n) = yval
  x(n,1) = rh
  x(n,2) = temp
  x(n,3) = rh**2
  x(n,4) = rh*temp
  x(n,5) = temp**2
  xrow(1:5) = x(n,:)
  xrow(0) = 1.0_dp
  CALL includ(wt, xrow, y(n))
END DO

CALL regcf(b, 6, ier)
WRITE(*, *) 'Fitted quadratic surface:'
WRITE(*, '(a, 2f8.3, a, f8.4, a, f8.4, a, f8.4, a, f8.4, a)')  &
         ' Y = ', b(0), b(1), '*X1 ', b(2), '*X2 ', b(3), '*X1^2 ',  &
         b(4), '*X1*X2 ', b(5), '*X2^2'
WRITE(*, *)
WRITE(*, *) 'Fitted values:'
WRITE(*, *) '   %RH          90       95      100'
WRITE(*, *) '-------------------------------------'
DO row = 1, 4
  xrow(0) = 1.0_dp
  xrow(1) = 20.0 - 10.0*row
  xrow(3) = xrow(1)**2
  DO col = 1, 3
    xrow(2) = -10.0 + 5.0*col
    xrow(4) = xrow(1)*xrow(2)
    xrow(5) = xrow(2)**2
    yfit(col) = DOT_PRODUCT( b, xrow )
  END DO
  WRITE(*, '(f8.0, "   ", 3f9.1)') xrow(1) + 70., yfit
END DO

STOP
END PROGRAM quad_surf

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


! =================================================================================================
! =================================================================================================
! ==========================                                           ============================
! ==========================     Legendre Polynomial Evaluation        ============================
! ==========================                                           ============================
! =================================================================================================
! =================================================================================================
!
! Author: Christopher A. Wong
! Date: 6 February 2016
! Format: Fortran 2008
!
! This module contains the procedures evaluate the Legendre polynomials and their derivatives up to
! a given degree on the d-dimensional unit cube.
!
! Dependencies: LAPACK
!
! Externally called procedures:
!   leg_buildpars   --- Initial subroutine called to form Legendre polynomial parameters
!   leg_eval        --- Main subroutine for evaluating Legendre polynomials.
!   leg_quad2       --- Computes Legendre matrix for 2D bilinear quadrature computations.
!
! Internal procedures:
!   choose          --- Computes binomial coefficients.

module legendre

    implicit none
    private
    integer,    parameter                       :: dp = kind(1.0d0)
    integer                                     :: maxdeg, d, numfuncs

    public  :: leg_eval, leg_quad2

contains

! =================================================================================================

! Binomial coefficient function

pure integer function choose (n, k)

    implicit none
    integer,    intent(in)  :: n
    integer,    intent(in)  :: k
    integer                 :: i,kk

    if (n-k < k) then
        kk = n-k
    else
        kk = k
    end if

    choose = 1

    if (kk >=1) then
        do i = 1, kk
            choose = (choose * (n-kk+i))/i
        end do
    end if

end function choose

! =================================================================================================

subroutine leg_buildpars(deg_in,d_in)

    implicit none

    integer,    intent(in)  :: deg_in, d_in

    maxdeg = deg_in
    d = d_in
    numfuncs = choose(maxdeg+d,d)

end subroutine leg_buildpars

! =================================================================================================
!
! Legendre polynomial evaluation
!
! Inputs:
!   x           --- input point as a dx1 array
!   need_deriv  --- optional, set to .TRUE. to compute derivatives, default is .FALSE.
!
! Outputs:
!   func_out    --- evaluations of multivariate orthogonal polynomials at x
!   deriv_out   --- values of the 1st order partial derivatives of the multivariate orthogonal
!       polynomials at x, given as a matrix with rows as functions
!   flag        --- set to .true. if an error occurs

subroutine leg_eval(x, need_deriv, func_out,deriv_out,flag)

implicit none

real(dp),   dimension(:),   intent(in)  :: x
logical,    optional,       intent(in)  :: need_deriv
real(dp),   dimension(:),   intent(out) :: func_out
real(dp),   dimension(:,:), intent(out) :: deriv_out
logical,                    intent(out) :: flag

logical                                 :: deriv_flag   ! Determines if derivatives are computed
integer                                 :: i, ii, j, jj, k, kk ! loop variables

integer,    dimension(1:choose(maxdeg+d,d),1:d) :: degidx ! stores degrees of multivar. polys
real(dp),   dimension(0:maxdeg,1:d)             :: leg1D, legder1D ! stores 1d Leg poly evaluations

! Initialize
func_out    = 1.0d0
deriv_out   = 1.0d0
flag        = .false.

if (present(need_deriv)) then
    deriv_flag = need_deriv
else
    deriv_flag = .false.
end if

! Check whether inputs have the correct sizes
if (size(func_out) /= numfuncs) then
    flag = .true.
else if ( (deriv_flag) .and. (size(deriv_out,1) /= numfuncs)) then
    flag = .true.
else if ( (deriv_flag) .and. (size(deriv_out,2) /=d) ) then
    flag = .true.
end if

! form the list of polynomial degrees in the tensor product of legendre polys
! this is equivalent to forming the set of all ordered partitions of positive
! integers up to maxdeg

if (d == 1) then
    degidx(1:numfuncs,1) = [ (j, j = 0, numfuncs-1) ]
else if (d == 2) then
    k = 1                               ! index counter
    do i = 0, maxdeg                    ! total degree
        do j = 0,i                      ! degree in first variable
            degidx(k,1:2) = [ j,i - j ]
            k = k + 1
        end do
    end do
else if (d == 3) then
    kk = 1                              ! index counter
    do i = 0, maxdeg                    ! total degree
        do j = 0, i                     ! total degree of the first two terms
            do k = 0, j                 ! degree in the first variable
                degidx(kk,1:3) = [ k, j-k, i-j ]
                kk = kk + 1
            end do
        end do
    end do

else
    flag = .true.
end if

! Evaluate 1D Legendre polynomials and derivatives at the d coordinates of x
! using the three-term recurrence formula for orthogonal polynomials

! set values for degrees 0 and 1
leg1D(0,1:d)    = 1.0d0
leg1D(1,1:d)    = x

if (deriv_flag) then
    legder1D(0,1:d) = 0.0d0
    legder1D(1,1:d) = 1.0d0
end if

do concurrent (i = 1:d) ! looping through d dimensions
    do j = 1, maxdeg -1 ! looping through degrees

        leg1D(j+1,i) =( (2*j+1) * x(i) * leg1D(j,i) - j*leg1D(j-1,i) ) / (j + 1.0d0)

        if (deriv_flag) then
            legder1D(j+1,i) = ( (2*j+1) * leg1D(j,i) + (2*j+1) * x(i) * legder1D(j,i) - &
                j* legder1D(j-1,i) ) / (j + 1.0d0)
        end if

    end do
end do

! renormalization
if (deriv_flag) then
    do concurrent (i = 0:maxdeg)
        leg1D(i,1:d)    = leg1D(i,1:d) * sqrt((2.0d0*i +1)/2)
        legder1D(i,1:d) = legder1D(i,1:d) * sqrt((2.0d0*i +1)/2)
    end do
else
    do concurrent (i = 0:maxdeg)
        leg1D(i,1:d)    = leg1D(i,1:d) * sqrt((2.0d0*i +1)/2)
    end do
end if

! Evaluate d-dimensional Legendre polynomials

if (d == 1) then
    func_out(1:numfuncs) = leg1D(0:numfuncs-1,1)

    if (deriv_flag) then
        deriv_out = legder1D
    end if
else
    do i = 1, numfuncs ! loop over all functions
        do j = 1,d ! loop over all dimensions
            func_out(i) = func_out(i) * leg1D(degidx(i,j),j)

            if (deriv_flag) then
                do k = 1,d ! loop over all partial derivatives

                    if (j /= k) then
                        deriv_out(i,k) = deriv_out(i,k) * leg1D(degidx(i,j),j)
                    else if (j == k) then
                        deriv_out(i,k) = deriv_out(i,k) * legder1D(degidx(i,j),j)
                    end if

                end do
            end if
        end do
    end do
end if

end subroutine leg_eval

! =================================================================================================
!
! This special subroutine is only used for computing bilinear quadrature that is exact on a space
! of polynomials and minimal on polynomials of the next degree. It is designed to be called with
! the NLopt optimization library.
!
subroutine leg_quad2(x,numpt,A_out)

    implicit none
    real(dp),   dimension(:),                   intent(in)  :: x        ! input array of points
    integer,                                    intent(in)  :: numpt    ! number of points
    real(dp),   dimension(:,:),                 intent(out) :: A_out    ! output matrix

    real(dp),   dimension(1:numpt,1:numfuncs)               :: output
    real(dp),   dimension(1:numpt,1:choose(maxdeg+1,2))     :: F
    real(dp),   dimension(1:numpt,1:(maxdeg+1))             :: G
    integer,    dimension(1:numpt)                          :: IPIV     ! LAPACK pivot info
    integer                                                 :: INFO     ! LAPACK error value
    logical                                                 :: flag

    ! ---------

    ! Check that a symmetric bilinear quadrature has been requested
    if (numpt /= (numfuncs - maxdeg - 1)) then
        print *, 'Error: Non-symmetric bilinear quadrature requested.'
        stop
    end if

    ! Evaluate matrix functions
    call leg_eval(reshape(x, (/ 2, numpt /) ), need_deriv = .false., output, 0.0d0, flag)
    F = output(:,1:(numfuncs - maxdeg -1))
    G = output(:,(numfuncs - maxdeg):numfuncs)

    ! Compute objective function A(x) = inv(F(x)) * G(x) using LAPACK routines
    call dgetrf(numpt,numpt,F,numpt,IPIV,INFO)
    call dgetrs('N',numpt,maxdeg+1,F,numpt,IPIV,G,numpt,INFO)

    A_out = G

end subroutine leg_quad2

! =================================================================================================

end module legendre


text:
      subroutine compute_jsxx(n_PHI2,n_PI2,n_U,n_a2,Nx,res)
      implicit none
      integer i
      integer Nx
      real*8 n_PHI2(Nx)
      real*8 n_PI2(Nx)
      real*8 n_U(Nx)
      real*8 n_a2(Nx)
      real*8 res(Nx)
      real*8 qb
      do i=1, Nx, 1
      qb = 0.5000000000000000D0 * n_PI2(i) + 0.5000000000000000D0 * n_PH
     #I2(i) - 0.5000000000000000D0 * n_a2(i) * n_U(i)
      res(i)=qb
      end do
      END


text:
Hello, Im Kurt Sassenrath, a Computer Science and Engineering Computer Science Engineering major at UC Davis since 2008.  Ill be the first to say that Im an academic at heart, with a passion for most of the sciences, though engineering ties me down.  I enjoy discussions on just about any topic provided that arguments do not deteriorate into debates, excluding the obvious facetious debates friends have with each other.  My typical activities include programming projects, gaming (video and board), bowling, and cooking.

Programming/scripting languages I have experience in include:
C/C++
REBOL
Java
Javascript
Python
Perl
PHP
LISP
Prolog

I have experience in web design (HTML/CSS) and Photoshop.  I enjoy creating art, and have created several brush set resources for Photoshop, which can be found http://dragonfire3k3.deviantart.com here.

Classes taken include, but are not limited to the following:
ECS  30, 40, 50, 60, 124A, 140A, 154AB, 165AB, 166
MAT  21AD, 22AB
PHY  9AD
CHE  2AB
ENG  17
NAS  5

A belated Welcome to the Wiki! Users/TomGarberson


text:
!>------------------------------------------------------------
!! Simple test program to verify that fftw libraries are properly linked
!!
!!  @author
!!  Ethan Gutmann
!!
!!------------------------------------------------------------
program fftw_test
    use, intrinsic :: iso_c_binding
    use fft

    implicit none

    integer*8 :: n_elements, nx, ny, nz

    type(C_PTR) :: fftw_aligned_data
    complex(C_DOUBLE_COMPLEX), pointer, dimension(:,:,:) :: fftw_fdata

    nx = 100
    ny = 100
    nz = 20

    n_elements = nx * ny * nz
    print*, "Testing fftw_alloc_complex"
    fftw_aligned_data = fftw_alloc_complex(n_elements)
    call c_f_pointer(fftw_aligned_data, fftw_fdata, [nx,ny,nz])

    fftw_fdata(1,1,1) = 10
    ! print*, fftw_fdata(1,1,1)
    if (real(fftw_fdata(1,1,1)) == 10) then
        print*, "PASSED"
    endif

end program fftw_test

max_stars_repo_name: scrasmussen/icar

max_issues_repo_name: scrasmussen/icar

max_forks_repo_name: scrasmussen/icar


text: Weirdo Begierdos Kids Korner is a radio show on KDVS running from 46 AM on thursdays featuring Sharmi Basu and Parijata Prabhakara.  Playing the finest in International, Eclectic, Electronic, Folk, Rock, Indie, HipHop, Noise, Experimental, and even some Bollywood classics, its the perfect blend for latenight studying, driving, or even Tweakers Tweaking.  Their stream can be foundhttp://www.kdvs.org/showme.cfm?show2024&titleWeirdo%20Begierdo%5Cs%20Kids%20Korner&back_grid here


text:
needs revagui/window
needs revagui/button
needs callbacks

with~ ~revagui

0 value w
0 value b
variable has-clicked

: mycallback callback ( object dataptr -- )
	has-clicked @ 0if
		b z" Don't click me again!" setlabel
		has-clicked on
	else
		." ARRGH!!!" cr bye
	then
	;

100 100 400 200 z" Hi there" window_xy to w
	w window_begin
	20 40 260 100 z" Click me!" button to b
	w window_end

b ' mycallback 1000 setcallback

w window_show run bye

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     This is for solving fluid transport equation solving Navier-Stokes
!     equations. Dirichlet boundary conditions are either treated
!     strongly, weakly or using penalty methods for immersed boundaries.
!
!--------------------------------------------------------------------

      SUBROUTINE FLUID3D (eNoN, w, N, Nx, al, yl, bfl, ksix, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN), al(tDof,eNoN),
     2   yl(tDof,eNoN), bfl(3,eNoN), ksix(3,3)
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      REAL(KIND=8), PARAMETER :: ct(2) = (/1D0,36D0/)
      INTEGER a, b
      REAL(KIND=8) tauM, tauC, tauB, kS, kU, nu, rho, T1, T2, T3, divU,
     2  amd, wl, wr, wrl, s, p, u(3), ud(3), px(3), f(3), up(3), ua(3),
     3  ux(3,3), udNx(eNoN), updNx(eNoN), uadNx(eNoN), rV(3), rM(3,3),
     4  NxdNx, nu_s, es(3,3), gam, nu_x, es_x(3,eNoN)

      rho  = eq(cEq)%dmn(cDmn)%prop(fluid_density)
      f(1) = eq(cEq)%dmn(cDmn)%prop(f_x)
      f(2) = eq(cEq)%dmn(cDmn)%prop(f_y)
      f(3) = eq(cEq)%dmn(cDmn)%prop(f_z)

      wr   = w*rho
      T1   = eq(cEq)%af*eq(cEq)%gam*dt
      amd  = eq(cEq)%am/T1
      wrl  = wr*T1
      wl   = w*T1

!     Indices are not selected based on the equation only
!     because fluid equation always come first
      p  = 0D0
      u  = 0D0
      ud = -f
      px = 0D0
      ux = 0D0
      DO a=1, eNoN
         p  = p + N(a)*yl(4,a)

         ud(1) = ud(1) + N(a)*(al(1,a)-bfl(1,a))
         ud(2) = ud(2) + N(a)*(al(2,a)-bfl(2,a))
         ud(3) = ud(3) + N(a)*(al(3,a)-bfl(3,a))

         px(1) = px(1) + Nx(1,a)*yl(4,a)
         px(2) = px(2) + Nx(2,a)*yl(4,a)
         px(3) = px(3) + Nx(3,a)*yl(4,a)

         u(1) = u(1) + N(a)*yl(1,a)
         u(2) = u(2) + N(a)*yl(2,a)
         u(3) = u(3) + N(a)*yl(3,a)

         ux(1,1) = ux(1,1) + Nx(1,a)*yl(1,a)
         ux(2,1) = ux(2,1) + Nx(2,a)*yl(1,a)
         ux(3,1) = ux(3,1) + Nx(3,a)*yl(1,a)
         ux(1,2) = ux(1,2) + Nx(1,a)*yl(2,a)
         ux(2,2) = ux(2,2) + Nx(2,a)*yl(2,a)
         ux(3,2) = ux(3,2) + Nx(3,a)*yl(2,a)
         ux(1,3) = ux(1,3) + Nx(1,a)*yl(3,a)
         ux(2,3) = ux(2,3) + Nx(2,a)*yl(3,a)
         ux(3,3) = ux(3,3) + Nx(3,a)*yl(3,a)
      END DO
      IF (mvMsh) THEN
         DO a=1, eNoN
            u(1) = u(1) - N(a)*yl(5,a)
            u(2) = u(2) - N(a)*yl(6,a)
            u(3) = u(3) - N(a)*yl(7,a)
         END DO
      END IF

!     Strain rate tensor 2*e_ij := (u_ij + u_ji)
      es(1,1) = ux(1,1) + ux(1,1)
      es(2,2) = ux(2,2) + ux(2,2)
      es(3,3) = ux(3,3) + ux(3,3)

      es(1,2) = ux(1,2) + ux(2,1)
      es(1,3) = ux(1,3) + ux(3,1)
      es(2,3) = ux(2,3) + ux(3,2)

      es(2,1) = es(1,2)
      es(3,1) = es(1,3)
      es(3,2) = es(2,3)

      DO a=1, eNoN
        es_x(1,a) = es(1,1)*Nx(1,a) + es(1,2)*Nx(2,a) + es(1,3)*Nx(3,a)
        es_x(2,a) = es(2,1)*Nx(1,a) + es(2,2)*Nx(2,a) + es(2,3)*Nx(3,a)
        es_x(3,a) = es(3,1)*Nx(1,a) + es(3,2)*Nx(2,a) + es(3,3)*Nx(3,a)
      END DO

!     Shear-rate := (2*e_ij*e_ij)^.5
      gam = es(1,1)*es(1,1) + es(1,2)*es(1,2) + es(1,3)*es(1,3) +
     2      es(2,1)*es(2,1) + es(2,2)*es(2,2) + es(2,3)*es(2,3) +
     3      es(3,1)*es(3,1) + es(3,2)*es(3,2) + es(3,3)*es(3,3)
      gam = SQRT(0.5D0*gam)

!     Compute viscosity based on shear-rate and chosen viscosity model
      CALL GETVISCOSITY(eq(cEq)%dmn(cDmn), gam, nu, nu_s, nu_x)
      nu   = nu/rho
      nu_s = nu_s/rho
      IF (ISZERO(gam)) THEN
         nu_x = 0D0
      ELSE
         nu_x = nu_x/rho/gam
      END IF
      s  = nu/eq(cEq)%dmn(cDmn)%prop(permeability)

      kU = u(1)*u(1)*ksix(1,1) + u(2)*u(1)*ksix(2,1)
     2   + u(3)*u(1)*ksix(3,1) + u(1)*u(2)*ksix(1,2)
     3   + u(2)*u(2)*ksix(2,2) + u(3)*u(2)*ksix(3,2)
     4   + u(1)*u(3)*ksix(1,3) + u(2)*u(3)*ksix(2,3)
     5   + u(3)*u(3)*ksix(3,3)

      kS = ksix(1,1)*ksix(1,1) + ksix(2,1)*ksix(2,1)
     2   + ksix(3,1)*ksix(3,1) + ksix(1,2)*ksix(1,2)
     3   + ksix(2,2)*ksix(2,2) + ksix(3,2)*ksix(3,2)
     4   + ksix(1,3)*ksix(1,3) + ksix(2,3)*ksix(2,3)
     5   + ksix(3,3)*ksix(3,3)

      tauM = 1D0/SQRT((2D0*ct(1)/dt)**2D0 + kU +
     2   ct(2)*nu_s*nu_s*kS + s*s)

      up(1) = -tauM*(ud(1) + px(1)/rho + u(1)*ux(1,1) + u(2)*ux(2,1)
     2      + u(3)*ux(3,1) + s*u(1))
      up(2) = -tauM*(ud(2) + px(2)/rho + u(1)*ux(1,2) + u(2)*ux(2,2)
     2      + u(3)*ux(3,2) + s*u(2))
      up(3) = -tauM*(ud(3) + px(3)/rho + u(1)*ux(1,3) + u(2)*ux(2,3)
     2      + u(3)*ux(3,3) + s*u(3))

      tauC = ksix(1,1) + ksix(2,2) + ksix(3,3)
      tauC = 1D0/tauM/tauC

      tauB = up(1)*up(1)*ksix(1,1) + up(2)*up(1)*ksix(2,1)
     2     + up(3)*up(1)*ksix(3,1) + up(1)*up(2)*ksix(1,2)
     3     + up(2)*up(2)*ksix(2,2) + up(3)*up(2)*ksix(3,2)
     4     + up(1)*up(3)*ksix(1,3) + up(2)*up(3)*ksix(2,3)
     5     + up(3)*up(3)*ksix(3,3)

      IF (ISZERO(tauB)) tauB = eps
      tauB = 1D0/SQRT(tauB)

      divU = ux(1,1) + ux(2,2) + ux(3,3)

      ua(1) = u(1) + up(1)
      ua(2) = u(2) + up(2)
      ua(3) = u(3) + up(3)

      rV(1) = tauB*(up(1)*ux(1,1) + up(2)*ux(2,1) + up(3)*ux(3,1))
      rV(2) = tauB*(up(1)*ux(1,2) + up(2)*ux(2,2) + up(3)*ux(3,2))
      rV(3) = tauB*(up(1)*ux(1,3) + up(2)*ux(2,3) + up(3)*ux(3,3))

      rM(1,1) = nu*es(1,1) - up(1)*ua(1) + rV(1)*up(1) + divU*tauC
     2        - p/rho
      rM(2,1) = nu*es(2,1) - up(2)*ua(1) + rV(2)*up(1)
      rM(3,1) = nu*es(3,1) - up(3)*ua(1) + rV(3)*up(1)

      rM(1,2) = nu*es(1,2) - up(1)*ua(2) + rV(1)*up(2)
      rM(2,2) = nu*es(2,2) - up(2)*ua(2) + rV(2)*up(2) + divU*tauC
     2        - p/rho
      rM(3,2) = nu*es(3,2) - up(3)*ua(2) + rV(3)*up(2)

      rM(1,3) = nu*es(1,3) - up(1)*ua(3) + rV(1)*up(3)
      rM(2,3) = nu*es(2,3) - up(2)*ua(3) + rV(2)*up(3)
      rM(3,3) = nu*es(3,3) - up(3)*ua(3) + rV(3)*up(3) + divU*tauC
     2        - p/rho

      rV(1) = ud(1) + ua(1)*(s+ux(1,1)) + ua(2)*ux(2,1) + ua(3)*ux(3,1)
      rV(2) = ud(2) + ua(1)*ux(1,2) + ua(2)*(s+ux(2,2)) + ua(3)*ux(3,2)
      rV(3) = ud(3) + ua(1)*ux(1,3) + ua(2)*ux(2,3) + ua(3)*(s+ux(3,3))

      DO a=1, eNoN
         udNx(a)  = u(1)*Nx(1,a)  + u(2)*Nx(2,a)  + u(3)*Nx(3,a)
         updNx(a) = up(1)*Nx(1,a) + up(2)*Nx(2,a) + up(3)*Nx(3,a)
         uadNx(a) = updNx(a) + udNx(a)

         lR(1,a) = lR(1,a) + wr*(rV(1)*N(a) + rM(1,1)*Nx(1,a)
     2      + rM(1,2)*Nx(2,a) + rM(1,3)*Nx(3,a))

         lR(2,a) = lR(2,a) + wr*(rV(2)*N(a) + rM(2,1)*Nx(1,a)
     2      + rM(2,2)*Nx(2,a) + rM(2,3)*Nx(3,a))

         lR(3,a) = lR(3,a) + wr*(rV(3)*N(a) + rM(3,1)*Nx(1,a)
     2      + rM(3,2)*Nx(2,a) + rM(3,3)*Nx(3,a))

         lR(4,a) = lR(4,a) + w*(N(a)*divU - updNx(a))
      END DO

      DO a=1, eNoN
         DO b=1, eNoN
            rM(1,1) = Nx(1,a)*Nx(1,b)
            rM(2,1) = Nx(2,a)*Nx(1,b)
            rM(3,1) = Nx(3,a)*Nx(1,b)
            rM(1,2) = Nx(1,a)*Nx(2,b)
            rM(2,2) = Nx(2,a)*Nx(2,b)
            rM(3,2) = Nx(3,a)*Nx(2,b)
            rM(1,3) = Nx(1,a)*Nx(3,b)
            rM(2,3) = Nx(2,a)*Nx(3,b)
            rM(3,3) = Nx(3,a)*Nx(3,b)

            NxdNx = Nx(1,a)*Nx(1,b) + Nx(2,a)*Nx(2,b) + Nx(3,a)*Nx(3,b)

            T1 = nu*NxdNx + tauB*updNx(a)*updNx(b)
     2         + N(a)*((s+amd)*N(b) + uadNx(b))
     3         + tauM*uadNx(a)*(udNx(b) + (s+amd)*N(b))

            T2 = tauM*udNx(a)
            T3 = tauM*(amd*N(b) + udNx(b))
!           dM/dU
            lK(1,a,b)  = lK(1,a,b)  + wrl*((nu + tauC)*rM(1,1) + T1 +
     2         nu_x*es_x(1,a)*es_x(1,b))
            lK(2,a,b)  = lK(2,a,b)  + wrl*(nu*rM(2,1) + tauC*rM(1,2) +
     2         nu_x*es_x(1,a)*es_x(2,b))
            lK(3,a,b)  = lK(3,a,b)  + wrl*(nu*rM(3,1) + tauC*rM(1,3) +
     2         nu_x*es_x(1,a)*es_x(3,b))

            lK(5,a,b)  = lK(5,a,b)  + wrl*(nu*rM(1,2) + tauC*rM(2,1) +
     2         nu_x*es_x(2,a)*es_x(1,b))
            lK(6,a,b)  = lK(6,a,b)  + wrl*((nu + tauC)*rM(2,2) + T1 +
     2         nu_x*es_x(2,a)*es_x(2,b))
            lK(7,a,b)  = lK(7,a,b)  + wrl*(nu*rM(3,2) + tauC*rM(2,3) +
     2         nu_x*es_x(2,a)*es_x(3,b))

            lK(9,a,b)  = lK(9,a,b)  + wrl*(nu*rM(1,3) + tauC*rM(3,1) +
     2         nu_x*es_x(3,a)*es_x(1,b))
            lK(10,a,b) = lK(10,a,b) + wrl*(nu*rM(2,3) + tauC*rM(3,2) +
     2         nu_x*es_x(3,a)*es_x(2,b))
            lK(11,a,b) = lK(11,a,b) + wrl*((nu + tauC)*rM(3,3) + T1 +
     2         nu_x*es_x(3,a)*es_x(3,b))

!           dM/dP
            lK(4,a,b)  = lK(4,a,b)  - wl*(Nx(1,a)*N(b) - Nx(1,b)*T2)
            lK(8,a,b)  = lK(8,a,b)  - wl*(Nx(2,a)*N(b) - Nx(2,b)*T2)
            lK(12,a,b) = lK(12,a,b) - wl*(Nx(3,a)*N(b) - Nx(3,b)*T2)
!           dC/dU
            lK(13,a,b) = lK(13,a,b) + wl*(Nx(1,b)*N(a) + Nx(1,a)*T3)
            lK(14,a,b) = lK(14,a,b) + wl*(Nx(2,b)*N(a) + Nx(2,a)*T3)
            lK(15,a,b) = lK(15,a,b) + wl*(Nx(3,b)*N(a) + Nx(3,a)*T3)
!           dC/dP
            lK(16,a,b) = lK(16,a,b) + wl*(tauM*NxdNx)/rho
         END DO
      END DO

      RETURN
      END SUBROUTINE FLUID3D
!--------------------------------------------------------------------
      SUBROUTINE FLUID2D (eNoN, w, N, Nx, al, yl, bfl, ksix, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(2,eNoN), al(tDof,eNoN),
     2   yl(tDof,eNoN), bfl(2,eNoN), ksix(2,2)
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN),
     2   lK(dof*dof,eNoN,eNoN)

      REAL(KIND=8), PARAMETER :: ct(2) = (/1D0,36D0/)
      INTEGER a, b
      REAL(KIND=8) tauM, tauC, tauB, kS, kU, nu, rho, T1, T2, T3, divU,
     2  amd, wl, wr, wrl, s, p, u(2), ud(2), px(2), f(2), up(2), ua(2),
     3  ux(2,2), udNx(eNoN), updNx(eNoN), uadNx(eNoN), rV(2), rM(2,2),
     4  NxdNx, nu_s, es(2,2), gam, nu_x, es_x(2,eNoN)

      rho  = eq(cEq)%dmn(cDmn)%prop(fluid_density)
      f(1) = eq(cEq)%dmn(cDmn)%prop(f_x)
      f(2) = eq(cEq)%dmn(cDmn)%prop(f_y)

      wr   = w*rho
      T1   = eq(cEq)%af*eq(cEq)%gam*dt
      amd  = eq(cEq)%am/T1
      wrl  = wr*T1
      wl   = w*T1

!     Indices are not selected based on the equation only
!     because fluid equation always come first
      p  = 0D0
      u  = 0D0
      ud = -f
      px = 0D0
      ux = 0D0
      DO a=1, eNoN
         p = p + N(a)*yl(3,a)

         ud(1) = ud(1) + N(a)*(al(1,a)-bfl(1,a))
         ud(2) = ud(2) + N(a)*(al(2,a)-bfl(2,a))

         px(1) = px(1) + Nx(1,a)*yl(3,a)
         px(2) = px(2) + Nx(2,a)*yl(3,a)

         u(1) = u(1) + N(a)*yl(1,a)
         u(2) = u(2) + N(a)*yl(2,a)

         ux(1,1) = ux(1,1) + Nx(1,a)*yl(1,a)
         ux(2,1) = ux(2,1) + Nx(2,a)*yl(1,a)
         ux(1,2) = ux(1,2) + Nx(1,a)*yl(2,a)
         ux(2,2) = ux(2,2) + Nx(2,a)*yl(2,a)
      END DO
      IF (mvMsh) THEN
         DO a=1, eNoN
            u(1) = u(1) - N(a)*yl(4,a)
            u(2) = u(2) - N(a)*yl(5,a)
         END DO
      END IF

!     Strain rate tensor 2*e_ij := (u_ij + u_ji)
      es(1,1) = ux(1,1) + ux(1,1)
      es(2,2) = ux(2,2) + ux(2,2)

      es(1,2) = ux(1,2) + ux(2,1)
      es(2,1) = es(1,2)

      DO a=1, eNoN
        es_x(1,a) = es(1,1)*Nx(1,a) + es(1,2)*Nx(2,a)
        es_x(2,a) = es(2,1)*Nx(1,a) + es(2,2)*Nx(2,a)
      END DO

!     Shear-rate := (2*e_ij*e_ij)^.5
      gam = es(1,1)*es(1,1) + es(1,2)*es(1,2) + es(2,1)*es(2,1)
     2    + es(2,2)*es(2,2)
      gam = SQRT(0.5D0*gam)

!     Compute viscosity based on shear-rate and chosen viscosity model
      CALL GETVISCOSITY(eq(cEq)%dmn(cDmn), gam, nu, nu_s, nu_x)
      nu   = nu/rho
      nu_s = nu_s/rho
      IF (ISZERO(gam)) THEN
         nu_x = 0D0
      ELSE
         nu_x = nu_x/rho/gam
      END IF
      s  = nu/eq(cEq)%dmn(cDmn)%prop(permeability)

      kU = u(1)*u(1)*ksix(1,1) + u(2)*u(1)*ksix(2,1)
     2   + u(1)*u(2)*ksix(1,2) + u(2)*u(2)*ksix(2,2)

      kS = ksix(1,1)*ksix(1,1) + ksix(2,1)*ksix(2,1)
     2   + ksix(1,2)*ksix(1,2) + ksix(2,2)*ksix(2,2)

      tauM = 1D0/SQRT((2D0*ct(1)/dt)**2D0 + kU +
     2   ct(2)*nu_s*nu_s*kS + s*s)

      up(1) = -tauM*(ud(1) + px(1)/rho + u(1)*ux(1,1) + u(2)*ux(2,1)
     2      + s*u(1))
      up(2) = -tauM*(ud(2) + px(2)/rho + u(1)*ux(1,2) + u(2)*ux(2,2)
     2      + s*u(2))

      tauC = ksix(1,1) + ksix(2,2)
      tauC = 1D0/tauM/tauC

      tauB = up(1)*up(1)*ksix(1,1) + up(2)*up(1)*ksix(2,1)
     2     + up(1)*up(2)*ksix(1,2) + up(2)*up(2)*ksix(2,2)

      IF (ISZERO(tauB)) tauB = eps
      tauB = 1D0/SQRT(tauB)

      divU = ux(1,1) + ux(2,2)

      ua(1) = u(1) + up(1)
      ua(2) = u(2) + up(2)

      rV(1) = tauB*(up(1)*ux(1,1) + up(2)*ux(2,1))
      rV(2) = tauB*(up(1)*ux(1,2) + up(2)*ux(2,2))

      rM(1,1) = nu*es(1,1) - up(1)*ua(1) + rV(1)*up(1) + divU*tauC
     2        - p/rho
      rM(2,1) = nu*es(2,1) - up(2)*ua(1) + rV(2)*up(1)

      rM(1,2) = nu*es(1,2) - up(1)*ua(2) + rV(1)*up(2)
      rM(2,2) = nu*es(2,2) - up(2)*ua(2) + rV(2)*up(2) + divU*tauC
     2        - p/rho

      rV(1) = ud(1) + ua(1)*(s + ux(1,1)) + ua(2)*ux(2,1)
      rV(2) = ud(2) + ua(1)*ux(1,2) + ua(2)*(s + ux(2,2))

      DO a=1, eNoN
         udNx(a)  = u(1)*Nx(1,a)  + u(2)*Nx(2,a)
         updNx(a) = up(1)*Nx(1,a) + up(2)*Nx(2,a)
         uadNx(a) = updNx(a) + udNx(a)

         lR(1,a) = lR(1,a) + wr*(rV(1)*N(a) + rM(1,1)*Nx(1,a)
     2      + rM(1,2)*Nx(2,a))

         lR(2,a) = lR(2,a) + wr*(rV(2)*N(a) + rM(2,1)*Nx(1,a)
     2      + rM(2,2)*Nx(2,a))

         lR(3,a) = lR(3,a) + w*(N(a)*divU - updNx(a))
      END DO

      DO a=1, eNoN
         DO b=1, eNoN
            rM(1,1) = Nx(1,a)*Nx(1,b)
            rM(2,1) = Nx(2,a)*Nx(1,b)
            rM(1,2) = Nx(1,a)*Nx(2,b)
            rM(2,2) = Nx(2,a)*Nx(2,b)

            NxdNx = Nx(1,a)*Nx(1,b) + Nx(2,a)*Nx(2,b)

            T1 = nu*NxdNx + tauB*updNx(a)*updNx(b)
     2         + N(a)*((s+amd)*N(b) + uadNx(b))
     3         + tauM*uadNx(a)*(udNx(b) + (s+amd)*N(b))

            T2 = tauM*udNx(a)
            T3 = tauM*(amd*N(b) + udNx(b))

!           dM/dU
            lK(1,a,b)  = lK(1,a,b) + wrl*((nu + tauC)*rM(1,1) + T1 +
     2         nu_x*es_x(1,a)*es_x(1,b))
            lK(2,a,b)  = lK(2,a,b) + wrl*(nu*rM(2,1) + tauC*rM(1,2) +
     2         nu_x*es_x(1,a)*es_x(2,b))

            lK(4,a,b)  = lK(4,a,b) + wrl*(nu*rM(1,2) + tauC*rM(2,1) +
     2         nu_x*es_x(2,a)*es_x(1,b))
            lK(5,a,b)  = lK(5,a,b) + wrl*((nu + tauC)*rM(2,2) + T1 +
     2         nu_x*es_x(2,a)*es_x(2,b))

!           dM/dP
            lK(3,a,b)  = lK(3,a,b) - wl*(Nx(1,a)*N(b) - Nx(1,b)*T2)
            lK(6,a,b)  = lK(6,a,b) - wl*(Nx(2,a)*N(b) - Nx(2,b)*T2)

!           dC/dU
            lK(7,a,b) = lK(7,a,b)  + wl*(Nx(1,b)*N(a) + Nx(1,a)*T3)
            lK(8,a,b) = lK(8,a,b)  + wl*(Nx(2,b)*N(a) + Nx(2,a)*T3)

!           dC/dP
            lK(9,a,b) = lK(9,a,b)  + wl*tauM*NxdNx/rho
         END DO
      END DO

      RETURN
      END SUBROUTINE FLUID2D
!####################################################################
      PURE SUBROUTINE BFLUID (eNoN, w, N, y, h, nV, lR, lK)
      USE COMMOD
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), y(tDof), h, nV(nsd)
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN), lK(dof*dof,eNoN,eNoN)

      INTEGER a, b, i, j
      REAL(KIND=8) T1, wl, hc(nsd), udn, u(nsd)

      wl  = w*eq(cEq)%af*eq(cEq)%gam*dt
      udn = 0D0
      IF (mvMsh) THEN
         DO i=1, nsd
            j    = i + nsd + 1
            u(i) = y(i) - y(j)
            udn  = udn + u(i)*nV(i)
         END DO
      ELSE
         DO i=1, nsd
            u(i) = y(i)
            udn  = udn + u(i)*nV(i)
         END DO
      END IF

      udn = eq(cEq)%dmn(cDmn)%prop(backflow_stab)/2D0*
     2   eq(cEq)%dmn(cDmn)%prop(fluid_density)*(udn - ABS(udn))
      hc  = h*nV + udn*u

!     Here the loop is started for constructing left and right hand side
      IF (nsd .EQ. 2) THEN
         DO a=1, eNoN
            lR(1,a) = lR(1,a) - w*N(a)*hc(1)
            lR(2,a) = lR(2,a) - w*N(a)*hc(2)
            DO b=1, eNoN
               T1        = wl*N(a)*N(b)*udn
               lK(1,a,b) = lK(1,a,b) - T1
               lK(5,a,b) = lK(5,a,b) - T1
            END DO
         END DO
      ELSE
         DO a=1, eNoN
            lR(1,a) = lR(1,a) - w*N(a)*hc(1)
            lR(2,a) = lR(2,a) - w*N(a)*hc(2)
            lR(3,a) = lR(3,a) - w*N(a)*hc(3)
            DO b=1, eNoN
               T1 = wl*N(a)*N(b)*udn
               lK(1,a,b)  = lK(1,a,b)  - T1
               lK(6,a,b)  = lK(6,a,b)  - T1
               lK(11,a,b) = lK(11,a,b) - T1
            END DO
         END DO
      END IF

      RETURN
      END SUBROUTINE BFLUID
!####################################################################
      SUBROUTINE BWFLUID3D (eNoN, w, N, Nx, yl, ub, nV, tauB, lR, lK)
      USE COMMOD
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN), yl(tDof,eNoN),
     2   ub(3), nV(3), tauB(2)
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN), lK(dof*dof,eNoN,eNoN)

      INTEGER :: a, b
      REAL(KIND=8) :: rho, nu, T1, wr, wl, wrl, tauT, tauN, p, uhn, ubn,
     2   u(3), uh(3), ux(3,3), sigman(3), Nxn(eNoN), rV(3), rM(3,3),
     3   nu_s, es(3,3), gam, nu_x

      rho  = eq(cEq)%dmn(cDmn)%prop(fluid_density)
      T1   = eq(cEq)%af * eq(cEq)%gam * dt
      wr   = w * rho
      wl   = w * T1
      wrl  = wr * T1
      tauT = tauB(1) / rho
      tauN = tauB(2) / rho

      p    = 0D0
      u    = 0D0
      ux   = 0D0
      DO a=1, eNoN
         p = p + N(a)*yl(4,a)

         u(1) = u(1) + N(a)*yl(1,a)
         u(2) = u(2) + N(a)*yl(2,a)
         u(3) = u(3) + N(a)*yl(3,a)

         ux(1,1) = ux(1,1) + Nx(1,a)*yl(1,a)
         ux(2,1) = ux(2,1) + Nx(2,a)*yl(1,a)
         ux(3,1) = ux(3,1) + Nx(3,a)*yl(1,a)
         ux(1,2) = ux(1,2) + Nx(1,a)*yl(2,a)
         ux(2,2) = ux(2,2) + Nx(2,a)*yl(2,a)
         ux(3,2) = ux(3,2) + Nx(3,a)*yl(2,a)
         ux(1,3) = ux(1,3) + Nx(1,a)*yl(3,a)
         ux(2,3) = ux(2,3) + Nx(2,a)*yl(3,a)
         ux(3,3) = ux(3,3) + Nx(3,a)*yl(3,a)

         Nxn(a)  = Nx(1,a)*nV(1) + Nx(2,a)*nV(2) + Nx(3,a)*nV(3)
      END DO

      uh = 0D0
      IF (mvMsh) THEN
         DO a=1, eNoN
            uh(1) = uh(1) + N(a)*yl(5,a)
            uh(2) = uh(2) + N(a)*yl(6,a)
            uh(3) = uh(3) + N(a)*yl(7,a)
         END DO
      END IF
      ubn = (u(1)-ub(1))*nV(1) + (u(2)-ub(2))*nV(2) + (u(3)-ub(3))*nV(3)
      uhn = (u(1)-uh(1))*nV(1) + (u(2)-uh(2))*nV(2) + (u(3)-uh(3))*nV(3)
      uhn = (ABS(uhn) - uhn) / 2D0

!     Strain rate tensor 2*e_ij := (u_ij + u_ji)
      es(1,1) = ux(1,1) + ux(1,1)
      es(2,2) = ux(2,2) + ux(2,2)
      es(3,3) = ux(3,3) + ux(3,3)

      es(1,2) = ux(1,2) + ux(2,1)
      es(1,3) = ux(1,3) + ux(3,1)
      es(2,3) = ux(2,3) + ux(3,2)

      es(2,1) = es(1,2)
      es(3,1) = es(1,3)
      es(3,2) = es(2,3)

!     Shear-rate := (2*e_ij*e_ij)^.5
      gam = es(1,1)*es(1,1) + es(1,2)*es(1,2) + es(1,3)*es(1,3) +
     2      es(2,1)*es(2,1) + es(2,2)*es(2,2) + es(2,3)*es(2,3) +
     3      es(3,1)*es(3,1) + es(3,2)*es(3,2) + es(3,3)*es(3,3)
      gam = SQRT(0.5D0*gam)

!     Compute viscosity based on shear-rate and chosen viscosity model
      CALL GETVISCOSITY(eq(cEq)%dmn(cDmn), gam, nu, nu_s, nu_x)
      nu = nu/rho

      sigman(1) = nu*(es(1,1)*nV(1) + es(1,2)*nV(2) + es(1,3)*nV(3))
     2   - (p/rho)*nV(1)
      sigman(2) = nu*(es(2,1)*nV(1) + es(2,2)*nV(2) + es(2,3)*nV(3))
     2   - (p/rho)*nV(2)
      sigman(3) = nu*(es(3,1)*nV(1) + es(3,2)*nV(2) + es(3,3)*nV(3))
     2   - (p/rho)*nV(3)

      rV(1) = -sigman(1) + (tauT + uhn)*(u(1)-ub(1))
      rV(2) = -sigman(2) + (tauT + uhn)*(u(2)-ub(2))
      rV(3) = -sigman(3) + (tauT + uhn)*(u(3)-ub(3))

      DO a=1, eNoN
c         lR(1,a) = lR(1,a) + wr*( N(a)*(rV(1) + (tauT-tauN)*ubn*nV(1)) -
c     2      nu*(Nxn(a)*(u(1)-ub(1)) + Nx(1,a)*ubn) )
         lR(1,a) = lR(1,a) + wr*( N(a)*(rV(1) + (tauT-tauN)*ubn*nV(1)) -
     2      nu*(u(1)-ub(1))*(Nxn(a) + Nx(1,a)*nV(1)) )

c         lR(2,a) = lR(2,a) + wr*( N(a)*(rV(2) + (tauT-tauN)*ubn*nV(2)) -
c     2      nu*(Nxn(a)*(u(2)-ub(2)) + Nx(2,a)*ubn) )
         lR(2,a) = lR(2,a) + wr*( N(a)*(rV(2) + (tauT-tauN)*ubn*nV(2)) -
     2      nu*(u(2)-ub(2))*(Nxn(a) + Nx(2,a)*nV(2)) )

c         lR(3,a) = lR(3,a) + wr*( N(a)*(rV(3) + (tauT-tauN)*ubn*nV(3)) -
c     2      nu*(Nxn(a)*(u(3)-ub(3)) + Nx(3,a)*ubn) )
         lR(3,a) = lR(3,a) + wr*( N(a)*(rV(3) + (tauT-tauN)*ubn*nV(3)) -
     2      nu*(u(3)-ub(3))*(Nxn(a) + Nx(3,a)*nV(3)) )

         lR(4,a) = lR(4,a) - w*N(a)*ubn
      END DO

      DO a=1, eNoN
         DO b=1, eNoN
            T1 = (tauT + uhn)*N(a)*N(b) - nu*(N(a)*Nxn(b) + N(b)*Nxn(a))

            rM(1,1) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(1) -
     2         nu*(N(a)*Nx(1,b)*nV(1) + Nx(1,a)*N(b)*nV(1))
c            rM(1,2) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(2) -
c     2         nu*(N(a)*Nx(1,b)*nV(2) + Nx(1,a)*N(b)*nV(2))
c            rM(1,3) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(3) -
c     2         nu*(N(a)*Nx(1,b)*nV(3) + Nx(1,a)*N(b)*nV(3))
            rM(1,2) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(2) -
     2         nu*N(a)*Nx(1,b)*nV(2)
            rM(1,3) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(3) -
     2         nu*N(a)*Nx(1,b)*nV(3)

c            rM(2,1) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(1) -
c     2         nu*(N(a)*Nx(2,b)*nV(1) + Nx(2,a)*N(b)*nV(1))
            rM(2,1) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(1) -
     2         nu*N(a)*Nx(2,b)*nV(1)
            rM(2,2) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(2) -
     2         nu*(N(a)*Nx(2,b)*nV(2) + Nx(2,a)*N(b)*nV(2))
c            rM(2,3) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(3) -
c     2         nu*(N(a)*Nx(2,b)*nV(3) + Nx(2,a)*N(b)*nV(3))
            rM(2,3) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(3) -
     2         nu*N(a)*Nx(2,b)*nV(3)

c            rM(3,1) = (tauT-tauN)*N(a)*N(b)*nV(3)*nV(1) -
c     2         nu*(N(a)*Nx(3,b)*nV(1) + Nx(3,a)*N(b)*nV(1))
c            rM(3,2) = (tauT-tauN)*N(a)*N(b)*nV(3)*nV(2) -
c     2         nu*(N(a)*Nx(3,b)*nV(2) + Nx(3,a)*N(b)*nV(2))
            rM(3,1) = (tauT-tauN)*N(a)*N(b)*nV(3)*nV(1) -
     2         nu*N(a)*Nx(3,b)*nV(1)
            rM(3,2) = (tauT-tauN)*N(a)*N(b)*nV(3)*nV(2) -
     2         nu*N(a)*Nx(3,b)*nV(2)
            rM(3,3) = (tauT-tauN)*N(a)*N(b)*nV(3)*nV(3) -
     2         nu*(N(a)*Nx(3,b)*nV(3) + Nx(3,a)*N(b)*nV(3))

!           dM/dU
            lK(1,a,b) = lK(1,a,b) + wrl*(T1 + rM(1,1))
            lK(2,a,b) = lK(2,a,b) + wrl*rM(1,2)
            lK(3,a,b) = lK(3,a,b) + wrl*rM(1,3)

            lK(5,a,b) = lK(5,a,b) + wrl*rM(2,1)
            lK(6,a,b) = lK(6,a,b) + wrl*(T1 + rM(2,2))
            lK(7,a,b) = lK(7,a,b) + wrl*rM(2,3)

            lK(9,a,b)  = lK(9,a,b)  + wrl*rM(3,1)
            lK(10,a,b) = lK(10,a,b) + wrl*rM(3,2)
            lK(11,a,b) = lK(11,a,b) + wrl*(T1 + rM(3,3))

!           dM/dP
            lK(4,a,b)  = lK(4,a,b)  + wl*N(a)*N(b)*nV(1)
            lK(8,a,b)  = lK(8,a,b)  + wl*N(a)*N(b)*nV(2)
            lK(12,a,b) = lK(12,a,b) + wl*N(a)*N(b)*nV(3)

!           dC/dU
            lK(13,a,b) = lK(13,a,b) - wl*N(a)*N(b)*nV(1)
            lK(14,a,b) = lK(14,a,b) - wl*N(a)*N(b)*nV(2)
            lK(15,a,b) = lK(15,a,b) - wl*N(a)*N(b)*nV(3)
         END DO
      END DO

      RETURN
      END SUBROUTINE BWFLUID3D
!--------------------------------------------------------------------
      SUBROUTINE BWFLUID2D (eNoN, w, N, Nx, yl, ub, nV, tauB, lR, lK)
      USE COMMOD
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(2,eNoN), yl(tDof,eNoN),
     2   ub(2), nV(2), tauB(2)
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN), lK(dof*dof,eNoN,eNoN)

      INTEGER :: a, b
      REAL(KIND=8) :: rho, nu, T1, wr, wl, wrl, tauT, tauN, p, uhn, ubn,
     2   u(2), uh(2), ux(2,2), sigman(2), Nxn(eNoN), rV(2), rM(2,2),
     3   nu_s, es(2,2), gam, nu_x

      rho  = eq(cEq)%dmn(cDmn)%prop(fluid_density)
      T1   = eq(cEq)%af * eq(cEq)%gam * dt
      wr   = w*rho
      wl   = w  * T1
      wrl  = wr * T1
      tauT = tauB(1) / rho
      tauN = tauB(2) / rho

      p    = 0D0
      u    = 0D0
      ux   = 0D0
      DO a=1, eNoN
         p = p + N(a)*yl(3,a)

         u(1) = u(1) + N(a)*yl(1,a)
         u(2) = u(2) + N(a)*yl(2,a)

         ux(1,1) = ux(1,1) + Nx(1,a)*yl(1,a)
         ux(2,1) = ux(2,1) + Nx(2,a)*yl(1,a)
         ux(1,2) = ux(1,2) + Nx(1,a)*yl(2,a)
         ux(2,2) = ux(2,2) + Nx(2,a)*yl(2,a)

         Nxn(a)  = Nx(1,a)*nV(1) + Nx(2,a)*nV(2)
      END DO

      uh = 0D0
      IF (mvMsh) THEN
         DO a=1, eNoN
            uh(1) = uh(1) + N(a)*yl(4,a)
            uh(2) = uh(2) + N(a)*yl(5,a)
         END DO
      END IF
      uhn = (u(1)-uh(1))*nV(1) + (u(2)-uh(2))*nV(2)
      ubn = (u(1)-ub(1))*nV(1) + (u(2)-ub(2))*nV(2)
      uhn = (ABS(uhn) - uhn) / 2D0

!     Strain rate tensor 2*e_ij := (u_ij + u_ji)
      es(1,1) = ux(1,1) + ux(1,1)
      es(2,2) = ux(2,2) + ux(2,2)

      es(1,2) = ux(1,2) + ux(2,1)
      es(2,1) = es(1,2)

!     Shear-rate := (2*e_ij*e_ij)^.5
      gam = es(1,1)*es(1,1) + es(1,2)*es(1,2) + es(2,1)*es(2,1)
     2    + es(2,2)*es(2,2)
      gam = SQRT(0.5D0*gam)

!     Compute viscosity based on shear-rate and chosen viscosity model
      CALL GETVISCOSITY(eq(cEq)%dmn(cDmn), gam, nu, nu_s, nu_x)
      nu = nu/rho

      sigman(1) = -(p/rho)*nV(1) + nu*(es(1,1)*nV(1) + es(1,2)*nV(2))
      sigman(2) = -(p/rho)*nV(2) + nu*(es(2,1)*nV(1) + es(2,2)*nV(2))

      rV(1) = -sigman(1) + (tauT + uhn)*(u(1)-ub(1))
      rV(2) = -sigman(2) + (tauT + uhn)*(u(2)-ub(2))

      DO a=1, eNoN
c         lR(1,a) = lR(1,a) + wr*( N(a)*(rV(1) + (tauT-tauN)*ubn*nV(1)) -
c     2      nu*(Nxn(a)*(u(1)-ub(1)) + Nx(1,a)*ubn) )
         lR(1,a) = lR(1,a) + wr*( N(a)*(rV(1) + (tauT-tauN)*ubn*nV(1)) -
     2      nu*(u(1)-ub(1))*(Nxn(a) + Nx(1,a)*nV(1)) )

c         lR(2,a) = lR(2,a) + wr*( N(a)*(rV(2) + (tauT-tauN)*ubn*nV(2)) -
c     2      nu*(Nxn(a)*(u(2)-ub(2)) + Nx(2,a)*ubn) )
         lR(2,a) = lR(2,a) + wr*( N(a)*(rV(2) + (tauT-tauN)*ubn*nV(2)) -
     2      nu*(u(2)-ub(2))*(Nxn(a) + Nx(2,a)*nV(2)) )

         lR(3,a) = lR(3,a) - w*N(a)*ubn
      END DO

      DO a=1, eNoN
         DO b=1, eNoN
            T1 = (tauT + uhn)*N(a)*N(b) - nu*(N(a)*Nxn(b) + N(b)*Nxn(a))

            rM(1,1) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(1) -
     2         nu*(N(a)*Nx(1,b)*nV(1) + Nx(1,a)*N(b)*nV(1))
c            rM(1,2) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(2) -
c     2         nu*(N(a)*Nx(1,b)*nV(2) + Nx(1,a)*N(b)*nV(2))
            rM(1,2) = (tauT-tauN)*N(a)*N(b)*nV(1)*nV(2) -
     2         nu*N(a)*Nx(1,b)*nV(2)

c            rM(2,1) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(1) -
c     2         nu*(N(a)*Nx(2,b)*nV(1) + Nx(2,a)*N(b)*nV(1))
            rM(2,1) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(1) -
     2         nu*N(a)*Nx(2,b)*nV(1)
            rM(2,2) = (tauT-tauN)*N(a)*N(b)*nV(2)*nV(2) -
     2         nu*(N(a)*Nx(2,b)*nV(2) + Nx(2,a)*N(b)*nV(2))

!           dM/dU
            lK(1,a,b) = lK(1,a,b) + wrl*(T1 + rM(1,1))
            lK(2,a,b) = lK(2,a,b) + wrl*rM(1,2)

            lK(4,a,b) = lK(4,a,b) + wrl*rM(2,1)
            lK(5,a,b) = lK(5,a,b) + wrl*(T1 + rM(2,2))

!           dM/dP
            lK(3,a,b) = lK(4,a,b) + wl*N(a)*N(b)*nV(1)
            lK(6,a,b) = lK(8,a,b) + wl*N(a)*N(b)*nV(2)

!           dC/dU
            lK(7,a,b) = lK(7,a,b) - wl*N(a)*N(b)*nV(1)
            lK(8,a,b) = lK(8,a,b) - wl*N(a)*N(b)*nV(2)
         END DO
      END DO

      RETURN
      END SUBROUTINE BWFLUID2D
!####################################################################
      SUBROUTINE GETVISCOSITY(lDmn, gamma, mu, mu_s, mu_x)
      USE COMMOD
      IMPLICIT NONE

      TYPE(dmnType), INTENT(IN) :: lDmn
      REAL(KIND=8), INTENT(INOUT)  :: gamma
      REAL(KIND=8), INTENT(OUT) :: mu, mu_s, mu_x

      REAL(KIND=8) :: mu_i, mu_o, lam, a, n, T1, T2

      SELECT CASE (lDmn%visc%viscType)
      CASE (viscType_Const)
         mu   = lDmn%visc%mu_i
         mu_s = mu
         mu_x = 0D0

      CASE (viscType_CY)
         mu_i = lDmn%visc%mu_i
         mu_o = lDmn%visc%mu_o
         lam  = lDmn%visc%lam
         a    = lDmn%visc%a
         n    = lDmn%visc%n

         T1   = 1.0D0 + (lam*gamma)**a
         T2   = T1**((n-1.0D0)/a)
         mu   = mu_i + (mu_o-mu_i)*T2
         mu_s = mu_i

         T1   = T2/T1
         T2   = lam**a * gamma**(a-1.0D0) * T1
         mu_x = (mu_o-mu_i)*(n-1.0D0)*T2

      CASE (viscType_Cass)
         mu_i = lDmn%visc%mu_i
         mu_o = lDmn%visc%mu_o
         lam  = lDmn%visc%lam

         IF (gamma .LT. lam) THEN
            mu_o  = mu_o/SQRT(lam)
            gamma = lam
         ELSE
            mu_o = mu_o/SQRT(gamma)
         END IF
         mu   = (mu_i + mu_o) * (mu_i + mu_o)
         mu_s = mu_i*mu_i
         mu_x = 2.0D0*mu_o*(mu_o + mu_i)/gamma

      END SELECT

      RETURN
      END SUBROUTINE GETVISCOSITY
!####################################################################


text:
program test_cmake_export
  use ${project_name_lower}$, only : get_version
  implicit none

  integer :: major, minor, patch

  call get_version(major, minor, patch)
  print "('Library version ',I0,'.',I0,'.',I0,' succesfully linked')", major, minor, patch

end program test_cmake_export

max_stars_repo_name: awvwgk/fypro

max_issues_repo_name: awvwgk/fypro

max_forks_repo_name: awvwgk/fypro


text:
      subroutine gendenspsibulge

      include 'commonblocks'

      real psi

      open(file='denspsibulge.dat',unit=50,status='replace')

      fcut_bulge = dfbulge(psic)
      fcut_bulge = 0.

      do i=1,npsi
         psi = tableE(i)
         call getbulgedens(psi,denspsibulge(i))
         if(i.eq.npsi) denspsibulge(i) = 0.
         write(50,*) tableE(i),denspsibulge(i)
      enddo

      close(50)
      return

      end

      subroutine getbulgedens(psi,rho)

      include 'commonblocks'
      real psi

      m = 100
      v02 = v0*v0
      vmin = -10.
      vmax = log(sqrt(2.*(psi-psic)))
      dlogv = (vmax - vmin)/float(m-1)
      rho = 0.
      do j=1,4
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfbulge(e)-fcut_bulge
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*coef(j)*df_lowered
      enddo
      do j=5,m-4
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfbulge(e)-fcut_bulge
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*df_lowered
      enddo
      do jj=4,2,-1
         j = m-jj+1
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfbulge(e)-fcut_bulge
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*coef(jj)*df_lowered
      enddo
      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
Renters Insurance is like homeowners insurance for renters.  Some apartments apartment complexes in town require tenants to carry a renters insurance policy, but dont despairits pretty cheap.  Many (most?  all?) auto insurance companies give hefty discounts on your auto policy if you pick up a second policyhomeowners, renters, whatever.  You may even end up with net savings, actually saving more on your auto policy than you pay for your renters policy (State Farm Insurance example: with a twocar policy, a baseline $80/yr. renters policy knocked more than $100/yr. off one of the car policies).

Renters insurance is different from homeowners insurance in a few ways. Policies usually cover damage to your personal belongings and liability for injury in your rental property. The policy doesnt cover damage to the building itself—your landlord should have a policy for that.


text:
	subroutine offscheck(ch,iraw,vdata,dbspec)
c
c	include 'efflen.for'                     ! file giving effective lengths
c
c	this routine is to be called when an FFT (FFTH,FFTM or FFTL) event
c	has been found.
c	ch is the channel number, from 'w_channel_open'in the calling prog
c	It returns iraw, vdata, dbspec, and tmspec
c	iraw = 0 for an fft spectrum, iraw =
c		1 for "raw" = time series data,
c	if iraw = 1, vdata contains the time series in volts/m or nT,
c		corrected for the frequency response of the system.
c	dbspec is spectrum, corrected for the frequency response, in dB
c		with respect to  1. (volt/m)**2/hz rms,
c		or 1. nT**2/Hz in the case of the search coils.
c	dbspec(1) is the DC term, (not logarithmic) in volt/m or nT,
c	dbspec(2) for is the power at the fundamental frequency
c	in dbspec(n), n ranges from 1 to 513
c	tmspec is equivalent to the item returned when on-board fft's
c	are returned, i.e. 1/2 dB steps above an arbitrary reference.
c	not corrected for frequency response
c
c	for my own use, two data sets are returned in common blocks:
c		in common /rawvolts/ are returned the set of voltage
c		samples (at preamp input) not corrected for freq. response
c		in common /telem/ tmspec is the spectrum equivalent to
c		the uncorrected telemetry spectrum (an integer)
c
	common /rawvolts/ DATA(1026)
	common /telem/ tmspec(1024)
	integer*4 ch,iraw,ifftch,ipa,ok,okl,ichgp,size,return_size
	integer*4 mantissa,gain,tmspec,pwrw
	integer*4 len_index
	integer*4 w_item_i4,w_item_r4
	complex preamp,fftamp,fft_filter
	real ffund(3),dbspec(513),vdata(1024)
	character*32 item
	data size /1024/
	data ffund /21.34766,5.336914,.3335571/
c
c	This is  (3/1.05) 10**-3 volts /step.
c	This is at gain step 3, so divide by 16.3**3
c	DATA VPSTEP/ .24E-6/              ! used to 17 FEb 1996
c	The above was in error.
C	REDONE 30 JAN 1997.  STEVES COUNTS ARE ALSO P-P, AND FOR
C	FLIGHT PREAMPS, IT IS AVERAGE OF 1.077 COUNTS PER 3 mV
C	SEE N.B. VOL II P.80.  VOLTS/COUNT IS THEN (3 10^-3/1.077)
C	DIVIDED BY (16.3)^3 = .643E-6

C	data Volts_per_step /.24e-6/		!old value
C	data Volts_per_step /.66e-6/		!latest value from PJK 27.2.96
	data Vpstep /.643e-6/		!latest value from PJK 30.1.97
c
c	   ! now get the items
c
	   item = 'channel_number'
	   ok = w_item_i4(ch, item, ifftch, 1, return_size)
	   item = 'source'
	   ok = w_item_i4(ch, item, ipa, 1, return_size)
c
	   ichgp = (ifftch+1)/4 + 1               ! 1 = hi, 2 = mid, 3 = lo
	   fundfr = ffund(ichgp)
c
c
	IF(ipa.ge.7) THEN			! SEARCH COILS
	   item = 'BZ_LENGTH_EFF'
	   okl = w_item_R4(ch, item, efflen, 1, return_size)
	ELSEIF(ipa.eq.1.or.ipa.eq.4) then
	   item = 'EX_LENGTH_EFF'
	   okl = w_item_R4(ch, item, efflen, 1, return_size)
	ELSEIF(ipa.eq.2.or.ipa.eq.5) then
	   item = 'EY_LENGTH_EFF'
	   okl = w_item_R4(ch, item, efflen, 1, return_size)
	ELSE
	   item = 'EZ_LENGTH_EFF'
	   okl = w_item_R4(ch, item, efflen, 1, return_size)
	ENDIF
C	print*, 'ipa,efflen',ipa,efflen
C	write(37,*) 'ipa,efflen',ipa,efflen
C
	   ok = w_item_i4(ch, 'PACKET_SUBTYPE', iraw, 1, return_size)
C	 print*,'packet subtype, iraw =', iraw
	   if(iraw.ne.1) iraw = 0

c
	   if(iraw.eq.0) then
	     item = 'DATA'
	     ok = w_item_i4(ch, item, tmspec, size, return_size)
c
	     dbspec(1) = .5*tmspec(1)
c
	     const   =    40.*(alog10(512./vpstep))
     1			+ 40.*alog10(efflen)
     1			+ 20.*alog10(fundfr) - 124.	    ! half dB steps
c
	     do n = 2,513
	  	f = (n-1)*fundfr
	  	pgain = CABS(preamp(ipa,f))
	 	pgain = PGAIN*CABS(fftamp(ifftch,f))*cabs(fft_filter(ifftch,f))
	  	dbspec(n) = .5*(tmspec(n) - const)-20.*ALOG10(PGAIN)
	     enddo
c
	   endif
C
	if (iraw.eq.1) then
C
C	  WINDFFT DOES FFT AS IN FLIGHT SOFTWARE
C
	      CALL CHECKFFT(CH,IFFTCH,IPA,VDATA,TMSPEC,DBSPEC)
C
c	  print*,'ch no.,etc',ifftch,fmax,dbmax,dbspec(maxn-1),dbspec(maxn+1)

c     1    ,PREAMP(IPA,FMAX),FFTAMP(IFFTCH,FMAX),SPMAX
	endif
c
	return
	end
	SUBROUTINE CHECKFFT(TMCH,CHAN,IRX,VDATA,TMSPEC,SPECT)
C
C	THIS WINDFFT DOES FFT AS IN FLIGHT SOFTWARE
C	ON RETURN, DATA IS A TIME SERIES, IN TELEMETRY COUNTS BUT CORRECTED
C	FOR GAINS, TMSPEC IS IN HALF DB STEPS, IDENTICAL TO THE RETURN
C	FROM ON-BOARD FFT, AND SPECT IS THE SPECTRUM IN DB, CORRECTED FOR
C	GAINS BUT STILL RELATIVE TO AN ARBITRARY LEVEL BASED ON TELEMETRY NUMBERS
C	THE MAIN PROGRAM CONVERTS THESE TO VOLTS/M OR nT
C
	REAL*4	 WINDOW(1024)
	COMMON /SPBLK/ PWRW(1024),PWRNW(1024),FDATA(1024),CHDATA(1024)
	common /rawvolts/ DATA(1026)
	integer*4	gain(1024),mantissa(1024),TMSPEC(1024)
c
	COMPLEX 	PREAMP
	COMPLEX 	FFTAMP
	COMPLEX 	FFT_FILTER
	COMPLEX		CCORR,FCOEF,FCOEFT
	character*32	item
	REAL*4 		SPECT(513),DATA,VOLT(1030),OFFS(4,10)
	REAL*4		SPHASE(1024),VDATA(1024)
	REAL*4		FREQMAX,FFUND(3)
	REAL*4		HIPASS(10)
	INTEGER*4 	TMCH,size,return_size,w_item_i4
	integer*4	IRX
	integer*4	CHAN
	data size /1024/
	data ffund /21.34766,5.336914,.3335571/
	DATA HIPASS /2*85.,4*30.,4*.45/
C	DATA HIPASS /2*150.,4*8.,4*.45/     ! fftm 3 dB pt.
C
	DATA TWOPI /6.2831853/
	DATA WINDOW /1024*1./		! PLACEHOLDER, NOT USED
C
C	OFFSETS IN FLIGHT PROMS, ACCORDING TO RED NOTEBOOK
C
	DATA OFFS /2048. ,2050. ,2046. ,2048.,   ! CH 1 GN STP 0-3
     2             2048. ,2050. ,2046. ,2051.,   ! CH 2 GN STP 0-3
     3             2048. ,2082. ,2041. ,2043.,   ! etc.
     4             2048. ,2078. ,2042. ,2046.,
     5             2048. ,2048. ,2048. ,2047.,
     6             2048. ,2075. ,2043. ,2045.,
     7             2048. ,2048. ,2048. ,2037.,
     8             2048. ,2048. ,2048. ,2040.,
     9             2048. ,2048. ,2048. ,2036.,
     X             2048. ,2048. ,2287. ,2048./
C
	NPROS = 1024
	NFRQ = NPROS/2
C******
	OFFS(1,3) = 1500.
	      item = 'MANTISSA'
	      ok = w_item_i4(tmch,item,mantissa, size, return_size)

	      item = 'EXPONENT'
	      ok = w_item_i4(tmch, item, gain, size, return_size)

C
	datamax = 0.
	DO IK = 1,1024
	  NTEMP = MANTISSA(IK).AND.4095
	  NEX = GAIN(IK).AND.3
	  TEMP = -(NTEMP-OFFS(NEX+1,CHAN))
	  CHDATA(IK) = TEMP*(-16.3)**NEX
	  datamax = amax1(abs(chdata(ik)),datamax)
	ENDDO
C	print*,'max chdata',datamax
C
	CALL TDPLOT(TMCH)
C
	RETURN
	END

ext: for


module ModpathBasicDataModule
  use UtilMiscModule,only : utrimall
  use ModflowRectangularGridModule,only : ModflowRectangularGridType
  implicit none

! Set default access status to private
  private

! Public derived data type definitions
!--------------------------------------
! type:
!--------------------------------------
  type,public :: ModpathBasicDataType
      integer :: DefaultIfaceCount
      doubleprecision :: HNoFlow
      doubleprecision :: HDry
      integer,dimension(:),allocatable :: IBound
      integer,dimension(:),allocatable :: DefaultIfaceValues
      character(len=16),dimension(:),allocatable :: DefaultIFaceLabels
      doubleprecision,dimension(:),allocatable :: Porosity
  contains
    procedure :: ReadData=>pr_ReadData
  end type

contains

  subroutine pr_ReadData(this, inUnit, outUnit, grid)
  use utl7module,only : urdcom, upcase
  use UTL8MODULE,only : urword, ustop, u1dint, u1drel, u1ddbl, u8rdcom,         &
    u3dintmp, u3dintmpusg, u3ddblmp, u3ddblmpusg
  implicit none
  class(ModpathBasicDataType) :: this
  integer,intent(in) :: inUnit, outUnit
  integer,dimension(:),allocatable :: cellsPerLayer
  integer,dimension(:),allocatable :: layerTypes
  class(ModflowRectangularGridType),intent(inout) :: grid
  character(len=200) :: line
  character(len=16) :: txt
  integer :: n, m, nn, length, iface, errorCode, layer
  character(len=24),dimension(2) :: aname
  data aname(1) /'          BOUNDARY ARRAY'/
  data aname(2) /'                POROSITY'/

  if(allocated(this%IBound)) deallocate(this%IBound)
  if(allocated(this%Porosity)) deallocate(this%Porosity)
  allocate(this%IBound(grid%CellCount))
  allocate(this%Porosity(grid%CellCount))
  allocate(cellsPerLayer(grid%LayerCount))
  do n = 1, grid%LayerCount
      cellsPerLayer(n) = grid%GetLayerCellCount(n)
  end do

  ! Write header to the listing file
  write(outUnit, *)
  write(outUnit, '(1x,a)') 'MODPATH basic data file data'
  write(outUnit, '(1x,a)') '----------------------------'

  ! READ MPBAS DATA
  ! READ AND PRINT COMMENTS.
  call u8rdcom(inUnit,outUnit,line,errorCode)

  ! No flow and dry cell head flags
  if(grid%GridType .gt. 2) then
      this%HNoFlow = 1.0E+30
      this%HDry = -1.0E+30
      read(line, *) this%DefaultIfaceCount
      write(outUnit,'(1X,A,1PG12.5,A)') 'Aquifer head is set to ', this%HNoFlow,    &
        ' at cells with IDOMAIN = 0.'
      write(outUnit,'(1X,A,1PG12.5,A)') 'Aquifer head is set to ', this%Hdry,       &
        ' at all dry cells when Newton-Raphson rewetting is not active.'
      write(outUnit, *)
  else
      read(line,*) this%HNoFlow, this%HDry
      read(inUnit, *) this%DefaultIfaceCount
      write(outUnit,'(1X,A,1PG12.5,A)') 'Aquifer head is set to ', this%HNoFlow,    &
        ' at all no-flow cells (IBOUND=0).'
      write(outUnit,'(1X,A,1PG12.5,A)') 'Aquifer head is set to ', this%Hdry,       &
        ' at all dry cells.'
      write(outUnit, *)
  end if

      ! READ NUMBER OF STRESS PACKAGES THAT DO NOT HAVE AUXILIARY IFACE SUPPORT
      allocate (this%DefaultIfaceValues(this%DefaultIfaceCount))
      allocate (this%DefaultIfaceLabels(this%DefaultIfaceCount))

      if(this%DefaultIfaceCount .gt. 0) then
        ! READ BUDGET LABELS AND DEFAULT IFACE
        do n = 1, this%DefaultIfaceCount
          read(inUnit,'(A)') line
          call utrimall(line)
          length = len_trim(line)
          if(length .eq. 0) THEN
            call ustop('Stress package labels cannot be blank strings.')
          else
            if(length .gt. 16) length = 16
            this%DefaultIfaceLabels(n) = line
            call upcase(this%DefaultIfaceLabels(n))
          end if

          ! READ DEFAULT IFACE
          read(inUnit,*) this%DefaultIfaceValues(n)

        end do

        ! Check for duplicate labels
        do n = 1, this%DefaultIfaceCount
          txt = this%DefaultIfaceLabels(n)
          do m = 1, this%DefaultIfaceCount
            if(m .NE. n) then
              if(txt .eq. this%DefaultIfaceLabels(m)) THEN
               call ustop('Duplicate stress package labels are not allowed.')
              end if
            end if
          end do
        end do

        write(outUnit, *)
        write(outUnit,'(1X,A,A)') 'Default stress package boundary face options (IFACE):'
        if(this%DefaultIfaceCount .eq. 0) THEN
          write(outUnit,'(3X,A)') 'None were specified'
        else
          do n = 1, this%DefaultIfaceCount
            iface = this%DefaultIfaceValues(n)
            if(iface .eq. 0) THEN
              write(outUnit,'(3X,A,A)') trim(this%DefaultIfaceLabels(n)),       &
                ' will be treated as internal stresses (IFACE = 0)'
            else if((iface .lt. 0) .or. (iface .gt. 6)) THEN
              call ustop(' IFACE must have a value between 0 and 6.')
            else
              write(outUnit,'(3X,A,A,I2)') trim(this%DefaultIfaceLabels(n)),    &
                ' will be assigned to face ',IFACE
            end if
          end do
        end if

      end if

      ! LAYTYP
      ! Read LAYTYP only for MODFLOW-2005 and MODFLOW-USG datasets.
      ! Loop through layers and set the cell type for all cells in a
      ! layer equal to the LAYTYP value for the layer.
      ! LAYTYP is not read for MODFLOW-6 datasets. Instead, cell types
      ! for all cells are read directly for from the binary grid file.
      if((grid%GridType .eq. 1) .or. (grid%GridType .eq. 2)) then
          write(outUnit, *)
          write(outUnit,'(1X,A)') 'Layer type (LAYTYP)'
          allocate(layerTypes(grid%LayerCount))
          read(inUnit,*) (layerTypes(layer), layer = 1, grid%LayerCount)
          write(outUnit,'(1X,80I2)') (layerTypes(layer), layer = 1, grid%LayerCount)
          n = 0
          do layer = 1, grid%LayerCount
            if(layerTypes(layer) .lt. 0) layerTypes(layer) = 1
            do m = 1, cellsPerLayer(layer)
                n = n + 1
                grid%CellType(n) = layerTypes(layer)
            end do
          end do
      end if

      ! IBOUND
      if(grid%GridType .eq. 1) then
            call u3dintmp(inUnit, outUnit, grid%LayerCount, grid%RowCount,      &
              grid%ColumnCount, grid%CellCount, this%IBound, ANAME(1))
      else if(grid%GridType .eq. 2) then
          call u3dintmpusg(inUnit, outUnit, grid%CellCount, grid%LayerCount, this%IBound, &
            aname(1), cellsPerLayer)
      else if(grid%GridType .ge. 3) then
          do n = 1, grid%CellCount
              if(grid%GetIDomain(n) .lt. 1) then
                  this%IBound(n) = 0
              else
                  this%IBound(n) = 1
              end if
          end do
      else
            write(outUnit,*) 'Invalid grid type specified when reading IBOUND array data.'
            write(outUnit,*) 'Stopping.'
            call ustop(' ')
      end if

      ! POROSITY
      if((grid%GridType .eq. 1) .or. (grid%GridType .eq. 3)) then
            call u3ddblmp(inUnit, outUnit, grid%LayerCount, grid%RowCount,      &
              grid%ColumnCount, grid%CellCount, this%Porosity, ANAME(2))
      else if((grid%GridType .eq. 2) .or. (grid%GridType .eq. 4)) then
          call u3ddblmpusg(inUnit, outUnit, grid%CellCount, grid%LayerCount,              &
            this%Porosity, aname(2), cellsPerLayer)
      else
            write(outUnit,*) 'Invalid grid type specified when reading POROSITY array data.'
            write(outUnit,*) 'Stopping.'
            call ustop(' ')
      end if

  end subroutine pr_ReadData

end module  ModpathBasicDataModule


text:
! { dg-do compile }
! This tests a patch for a regression caused by the second part of
! the fix for PR30554.  The linked derived types dummy_atom and
! dummy_atom_list caused a segment fault because they do not have
! a namespace.
!
! Contributed by Daniel Franke <franke.daniel@gmail.com>
!
MODULE types
TYPE :: dummy_atom_list
  TYPE(dummy_atom), DIMENSION(:), POINTER :: table => null()
END TYPE

TYPE :: dummy_atom
  TYPE(dummy_atom_private), POINTER :: p => null()
END TYPE

TYPE :: dummy_atom_private
  INTEGER                     :: id
END TYPE
END MODULE

MODULE atom
USE types, ONLY: dummy_atom
INTERFACE
  SUBROUTINE dummy_atom_insert_symmetry_mate(this, other)
    USE types, ONLY: dummy_atom
    TYPE(dummy_atom), INTENT(inout) :: this
    TYPE(dummy_atom), INTENT(in)    :: other
  END SUBROUTINE
END INTERFACE
END MODULE

MODULE list
INTERFACE
  SUBROUTINE dummy_atom_list_insert(this, atom2)
    USE types, ONLY: dummy_atom_list
    USE atom, ONLY: dummy_atom

    TYPE(dummy_atom_list), INTENT(inout) :: this
    TYPE(dummy_atom), INTENT(in)         :: atom2
  END SUBROUTINE
END INTERFACE
END MODULE
! { dg-final { cleanup-modules "atom types list" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program mag2geo

use, intrinsic :: iso_fortran_env, only: real64
use apex, only : apex_setup,gm2gc

implicit none (type, external)

integer :: i, u, N
character(256) :: buf, nmlfn
real :: yeardec, r
real, allocatable, dimension(:) :: glat, glon, gmlat, gmlon
real(real64) :: tic, toc

namelist / input_shape / N
namelist / output / glat, glon
namelist / input / gmlat, gmlon

call get_command_argument(1, buf, status=i)
if(i/=0) error stop 'missing decimal year e.g. 2020.35'
read(buf, '(f8.4)', iostat=i) yeardec
if(i/=0) error stop 'command line should be like "./geo2mag yeardec in.nml out.nml"'

call get_command_argument(2, nmlfn, status=i)
if(i/=0) error stop 'specify input namelist file'

! -------------------------
!! read first, in case there is an IO issue, we don't want to waste time
!! initializing Apex if it will fail anyway
open(newunit=u, file=nmlfn, action='readwrite', status='old')
read(u, nml=input_shape)

allocate(glat(N), glon(N), gmlat(N), gmlon(N))

read(u, nml=input)

call cpu_time(tic)
call apex_setup(yeardec,altmax=1000.)
call cpu_time(toc)
print '(F5.1,A)', toc-tic, ' seconds to setup Apex'

call cpu_time(tic)
do i = 1,size(gmlat)
  call gm2gc(gmlat(i), gmlon(i), glat(i), glon(i))
end do
call cpu_time(toc)
print '(ES12.3,A,I8,A)', toc-tic, ' seconds to convert ', N, ' points.'

write(u, nml=output)
close(u)

end program

max_stars_repo_name: space-physics/apex-python

max_issues_repo_name: space-physics/apex-python

max_forks_repo_name: space-physics/apex-python


! { dg-do link }
!
! PR fortran/40568
!
! Module checks for C_SIZEOF (part of ISO_C_BINDING)
!

implicit none
intrinsic c_sizeof ! { dg-error "does not exist" }
end


text:
Paintings by Dave is owned by Dave Reynolds, a local Davis artist.  He sells greeting cards and paintings.

His greeting cards are on sale at:

  French Cuff Consignment
  Nestware
  Hallmark Inn Hotel
  All Things Right and Relevant


text:
      program riemann
      implicit none

c..solves the exact riemann problem
c..original code from bruce fryxell

c..declare
      integer    n,npts,itmax,iter,i
      parameter (n = 1000)
      real       x(n),rho(n),u(n),p(n),
     1           rhol,pl,ul,rhor,pr,ur,gamma,xi,t,xl,xr,
     2           rho1,p1,u1,rho5,p5,u5,p40,p41,f0,f,eps,
     3           f1,p4,error,z,c5,gm1,gp1,gmfac1,gmfac2,fact,
     4           u4,rho4,w,p3,u3,rho3,c1,c3,xsh,xcd,xft,xhd,dx

c..define initial conditions
c..state at left of discontinuity
         rhol = 10.0
         pl   = 100.0
         ul   = 0.

c..state at right of discontinuity
         rhor = 1.0
         pr   = 1.0
         ur   = 0.

         if (ul .ne. 0. .or. ur .ne. 0.) then
          write (6,*) 'must have ul = ur = 0.'
          stop
         endif

c.. equation of state
         gamma = 1.4

c..location of discontinuity at t = 0
         xi = 2.0

c..time at which solution is desired
         t = 0.40

c..number of points in solution
         npts = 500
         if (npts .gt. n) then
          write (6,*) 'number of points exceeds array size'
          stop
         endif

c..spatial interval over which to compute solution
         xl = 0.0
         xr = 5.0
         if (xr .lt. xl) then
          write (6,*) 'xr must be greater than xl'
          stop
         endif

c..begin solution
      if (pl .gt. pr) then
         rho1 = rhol
         p1   = pl
         u1   = ul
         rho5 = rhor
         p5   = pr
         u5   = ur
      else
         rho1 = rhor
         p1   = pr
         u1   = ur
         rho5 = rhol
         p5   = pl
         u5   = ul
      endif

c..solve for post-shock pressure by secant method
c..initial guesses

         p40 = p1
         p41 = p5
         f0  = f(p40, p1, p5, rho1, rho5, gamma)

c..maximum number of iterations and maxium allowable relative error
         itmax = 20
         eps   = 1.e-5

         do iter = 1, itmax
          f1 = f(p41, p1, p5, rho1, rho5, gamma)
          if (f1 .eq. f0) go to 10

          p4 = p41 - (p41 - p40) * f1 / (f1 - f0)

          error = abs (p4 - p41) / p41
          if (error .lt. eps) go to 10

          p40 = p41
          p41 = p4
          f0  = f1
         enddo
         write (6,*) 'iteration failed to converge'
         stop 'abnormal termination'

10       continue

c..compute post-shock density and velocity
         z  = (p4 / p5 - 1.)
         c5 = sqrt (gamma * p5 / rho5)

         gm1 = gamma - 1.
         gp1 = gamma + 1.
         gmfac1 = 0.5 * gm1 / gamma
         gmfac2 = 0.5 * gp1 / gamma

         fact = sqrt (1. + gmfac2 * z)

         u4 = c5 * z / (gamma * fact)
         rho4 = rho5 * (1. + gmfac2 * z)
     1               / (1. + gmfac1 * z)

c..shock speed
         w = c5 * fact

c..compute values at foot of rarefaction
         p3 = p4
         u3 = u4
         rho3 = rho1 * (p3 / p1)**(1. /gamma)

c..compute positions of waves
      if (pl .gt. pr) then
       c1 = sqrt (gamma * p1 / rho1)
       c3 = sqrt (gamma * p3 / rho3)

       xsh = xi + w * t
       xcd = xi + u3 * t
       xft = xi + (u3 - c3) * t
       xhd = xi - c1 * t

c..and do say what we found
       write (6, 500)
       write (6, 501) rho1, p1, u1
       write (6, 502)
       write (6, 503) rho3, p3, u3
       write (6, 504) rho4, p4, u4
       write (6, 505) rho5, p5, u5

       write (6, 506) xhd
       write (6, 507) xft
       write (6, 508) xcd
       write (6, 509) xsh

500    format (// 2x, 'Region', 4x, 'Density', 8x, 'Pressure',
     1                          8x, 'Velocity')
501    format (5x, '1', 3(2x,1pe14.7))
502    format (5x, '2' ,20x, 'RAREFACTION')
503    format (5x, '3', 3(2x,1pe14.7))
504    format (5x, '4', 3(2x,1pe14.7))
505    format (5x, '5', 3(2x,1pe14.7)//)

506    format (2x, 'Head Of Rarefaction    x = ', 1pe14.7)
507    format (2x, 'Foot Of Rarefaction    x = ', 1pe14.7)
508    format (2x, 'Contact Discontinuity  x = ', 1pe14.7)
509    format (2x, 'Shock                  x = ', 1pe14.7//)

c..compute solution as a function of position
       dx = (xr - xl) / (npts - 1)
       do i = 1, npts
        x(i) = xl + dx * (i - 1)
       enddo

       do i = 1, npts
        if (x(i) .lt. xhd) then
         rho(i) = rho1
         p(i)   = p1
         u(i)   = u1
        else if (x(i) .lt. xft) then
         u(i)   = 2. / gp1 * (c1 + (x(i) - xi) / t)
         fact   = 1. - 0.5 * gm1 * u(i) / c1
         rho(i) = rho1 * fact ** (2. / gm1)
         p(i)   = p1 * fact ** (2. * gamma / gm1)
        else if (x(i) .lt. xcd) then
         rho(i) = rho3
         p(i)   = p3
         u(i)   = u3
        else if (x(i) .lt. xsh) then
         rho(i) = rho4
         p(i)   = p4
         u(i)   = u4
        else
         rho(i) = rho5
         p(i)   = p5
         u(i)   = u5
        endif
       enddo
      endif

c..if pr > pl, reverse solution
      if (pr .gt. pl) then
       c1 = sqrt (gamma * p1 / rho1)
       c3 = sqrt (gamma * p3 / rho3)

       xsh = xi - w * t
       xcd = xi - u3 * t
       xft = xi - (u3 - c3) * t
       xhd = xi + c1 * t

c..and do say what we found
       write (6, 500)
       write (6, 501) rho5, p5, u5
       write (6, 602) rho4, p4, u4
       write (6, 503) rho3, p3, u3
       write (6, 604)
       write (6, 505) rho1, p1, u1

       write (6, 609) xsh
       write (6, 508) xcd
       write (6, 507) xft
       write (6, 606) xhd

602    format (5x, '2', 3(2x,1pe14.7))
604    format (5x, '4' ,20x, 'RAREFACTION')
606    format (2x, 'Head Of Rarefaction    x = ', 1pe14.7//)
609    format (2x, 'Shock                  x = ', 1pe14.7)

       dx = (xr - xl) / (npts - 1)
       do i = 1, npts
         x(i) = xl + dx * (i - 1)
       enddo

       do i = 1, npts
        if (x(i) .lt. xsh) then
           rho(i) = rho5
           p(i)   = p5
           u(i)   = -u5
        else if (x(i) .lt. xcd) then
           rho(i) = rho4
           p(i)   = p4
           u(i)   = -u4
        else if (x(i) .lt. xft) then
           rho(i) = rho3
           p(i)   = p3
           u(i)   = -u3
        else if (x(i) .lt. xhd) then
           u(i)   = -2. / gp1 * (c1 + (xi - x(i)) / t)
           fact   = 1. + 0.5 * gm1 * u(i) / c1
           rho(i) = rho1 * fact ** (2. / gm1)
           p(i)   = p1 * fact ** (2. * gamma / gm1)
        else
           rho(i) = rho1
           p(i)   = p1
           u(i)   = -u1
        endif
       enddo
      endif

      open (unit=1, file = 'output',status='unknown')
      write (1, 1000)
1000  format (2x, 'i', 10x, 'x', 12x, 'density', 8x, 'pressure',
     1             8x, 'velocity'/)
1001  format (i4, 4(2x, 1pe14.7))
      do i = 1, npts
       write (1,1001) i, x(i), rho(i), p(i), u(i)
      enddo
      close (1)
      stop 'normal termination'
      end

      real function f(p4, p1, p5, rho1, rho5, gamma)
      implicit none

c..shock tube equation

c..declare the pass
      real p4,p1,p5,rho1,rho5,gamma

c..local variables
      real z,c1,c5,gm1,gp1,g2,fact

      z = (p4 / p5 - 1.)
      c1 = sqrt (gamma * p1 / rho1)
      c5 = sqrt (gamma * p5 / rho5)

      gm1 = gamma - 1.
      gp1 = gamma + 1.
      g2  = 2. * gamma

      fact = gm1 / g2 * (c5 / c1) * z
     1     / sqrt (1. + gp1 / g2 * z)
      fact = (1. - fact) ** (g2 / gm1)

      f = p1 * fact - p4

      return
      end

max_stars_repo_name: ADACS-Australia/quokka

max_issues_repo_name: ADACS-Australia/quokka

max_forks_repo_name: ADACS-Australia/quokka


text:
Specialized in dTMS treatment in association with Dr. Guohua Xia, offering comprehensive mental health care from medication treatment, psychotherapy, to brain stimulation for all major mental disorders.
A research site for international randomized controlled clinical trials. Ongoing phase III clinical trials currently recruiting volunteer participants include: A study on Safety and Efficacy of deep TMS in treatment of Bipolar Depression,
A study for Evaluation of Deep Transcranial Magnetic Stimulation (DTMS) With the HADD Coil as an Aid to Smoking Cessation.


text:
c=======================================================================
c
c     subroutine ALLOCSRRFR
c
c     Sharpe Ratio (SR) maximisation strategy with a risk-free asset
c
c     Max[ Sharpe Ratio ]
c      s.t.
c     C*w <= b                       (linear constraints)
c     Cinf    <= w    <= Csup        (lower/upper bounds risky assets)
c     Cinfrfr <= wrfr <= Csuprfr     (lower/upper bounds risk-free asset)
c
c     where, Sharpe Ratio := (rho'*w - rfr)/sqrt(w'Qw)
c
c        w    : risky assets weights
c        wrfr : risk-free asset weights
c        Q    : covariance matrix
c        rho  : assets performance
c        rfr  : risk-free rate
c
c-----------------------------------------------------------------------
      SUBROUTINE allocsrrfr (n, cov, rho, rfr,
     &                       neq, nin, ccst, bcst,
     &                       cinf, csup, cinfrfr, csuprfr,
     &                       iwork, dwork, wopt, wrf, info )
c-----------------------------------------------------------------------
c
c     INPUT
c        n       : number of assets                              integer
c        cov     : covariance matrix (n*n)                        double
c        rho     : expected returns (n)                           double
c        rfr     : risk-free rate                                 double
c        neq     : number equality constraints                   integer
c        nin     : number inequality constraints                 integer
c        ccst    : matrix of constraints (n*(neq+nin))            double
c        bcst    : vector initial of constraints (neq+nin)        double
c        cinf    : lower bound (n)                                double
c        csup    : upper bound (n)                                double
c        cinfrfr : risk-free asset lower bound                     double
c        csuprfr : risk-free asset upper bound                     double
c
c     WORKSPACE
c       iwork    : 7*n + 2*nin + neq + 7                         integer
c       dwork    : n*(5*n+2*neq+3*nin+35)+4*neq+9*nin+27          double
c
c     OUTPUT
c       wopt     : optimal portfolio (n)                          double
c       wrf      : optimal risk-free asset                        double
c       info     : diagnostic arguments                          integer
c
c     CALL
c       TESTSDP : matrix SDP test
c       OPSR    : computing optimization for ALLOCSR
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin
      DOUBLE PRECISION rfr, cinfrfr, csuprfr, wrf
      DOUBLE PRECISION cov(*), rho(*), cinf(*), csup(*),
     &                 ccst(n,*), bcst(*), wopt(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER p, i, j, neqtot, nintot, ncst, ntot
      INTEGER piwo, pdwo, pdcov, pdrho, pdlagr, pdc, pdb, pdx,
     &        pdcinf, pdcsup
      DOUBLE PRECISION coef, mu, t, sum, rhomax, ZERO, EPS, NEPS, INF,
     &                 NINF
      PARAMETER (ZERO = 0.0, EPS = 1.E-12, NEPS = -1.E-12)
      PARAMETER (INF = 1.E+12, NINF = -1.E+12)
c
c     external subroutines
      EXTERNAL testsdp, opmv, qp, YMCPIR, CTSR, IMX, IVX
c
c-----------------------------------------------------------------------
c
c     initializations
      info = 0      ! diagnostic argument
      p    = n + 1  ! problem size
      CALL IVX ( n, wopt, ZERO )
      wrf = 0.0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piwo  = 1
c     piwo  : pointer for workspaces of TESTSDP, OPMV, QP
c
c     Total size of iwork array = 7*n + 2*nin + neq + 7
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
c      pdwo  = 1
c     pdwo  : pointer for workspaces of TESTSDP, OPMV, QP
c             n*( 2*n + neq + nin + 12 ) + 2*neq + 4*nin + 7
c     pdwo + n*( 2*n + neq + nin + 12 ) + 2*neq + 4*nin + 7
      pdcov = 1
c     pdcov : pointer for temporary cov. matrix (n+1)*(n+1)
      pdrho = pdcov + ( p*p )
c     pdrho : pointer for temporary rho. verctor (p)
      pdc   = pdrho + ( p )
c     pdc   : pointer for constraint matrix p*(2*n + nin + 2)
      pdb   = pdc + (  p*(2*n + nin + 2) )
c     pdb   : pointer for constraint vector (2*n + nin + 2)
      pdcinf= pdb + ( 2*n + nin + 2 )
c     pdcinf: pointer for lower bounds (p)
      pdcsup= pdcinf + ( p )
c     pdcsup: pointer for upper bounds (p)
      pdlagr= pdcsup + ( p )
c     pdlagr : pointer for Lagrange multipliers (p + 2 + 2*n + nin)
      pdx   = pdlagr + ( p + 2 + 2*n + nin )
c     pdx    : pointer for quadratic optimal points ( p )
      pdwo  = pdx + p
c     needs
c           Testsdp needs n*(2*n+7)
c           Quapro needs p*p + 6*p + 2*neq
c           OPMV needs 2*n*(neq+nin+n+9)+4*neq+6*nin+2*n+15
c           so  2*n*(neq+nin+n+9)+4*neq+6*nin+2*n+15

c     Total size of dwork array = n*(2*neq+2*nin+2*n+18)+4*neq+6*nin+2*n+15 + n*(3*n+15+nin) + 12 +3*nin
c                               = n*( 6*n + 2*neq + 3*nin + 35 ) + 4*neq + 9*nin + 27
c
c-----------------------------------------------------------------------
c
c     case cinfrfr = 100%
      IF ((cinfrfr .GT. 1-EPS).AND.(csuprfr .GT. 1-EPS)) THEN
        wrf = 1.0
        RETURN
      ENDIF
      wrf  = cinfrfr
      coef = 1.0/(1.0 - cinfrfr)
c
c     covariance matrix SDP test
      CALL testsdp (n, cov, iwork(piwo), dwork(pdwo), info)
      IF (info .NE. 0) THEN
         info = -108
         RETURN
      ENDIF
c
c     general linear constraints
      ncst   = neq + nin
      neqtot = neq + 1
      ntot   = ncst + 1
      DO i = 1,n
        dwork(pdc + i - 1) =  1.0
      ENDDO
      dwork(pdb) = coef*1.0
c     case ncst = 0 (neq=0, nin=0)
      IF (ncst .NE. 0) THEN
        DO i = 1,n
            DO j = 1,ncst
                dwork(pdc + n + n*(j-1) + i - 1) = ccst(i,j)
            ENDDO
        ENDDO
        DO j = 1,ncst
            dwork(pdb + j) = coef*bcst(j)
        ENDDO
      ENDIF
c
c     tests constraints compatibility
c      mu = -1.E+8
      DO i = 1,n
        dwork(pdcinf + i - 1) = coef*cinf(i)
        dwork(pdcsup + i - 1) = coef*csup(i)
      ENDDO
c
c     test if max. returns < rfr
      CALL EVMAX ( n, rho, rhomax )
      IF (rhomax .LT. rfr) THEN
        CALL IVX ( n, wopt, ZERO )
        wrf = 1.0
        info = 0
        IF (wrf .GT. (csuprfr+1.E-8)) THEN
            info = -113
        ENDIF
        RETURN
      ENDIF
c
c     test constraints
      mu = rfr + 1.E-10
      CALL opmv ( n, cov, rho, mu,
     &            neqtot, nin, dwork(pdc), dwork(pdb),
     &            dwork(pdcinf), dwork(pdcsup),
     &            iwork(piwo), dwork(pdwo), wopt, info )
      IF (info .LT. 0) THEN
        info = -100
        RETURN
      ENDIF
c
c     construction of the quadratic part
c     Q = | cov   0 |
c         | 0 ... 0 |
      CALL IMX ( p, p, dwork(pdcov), ZERO )
      DO i = 1,n
        iwork(piwo + i - 1) = i
      ENDDO
      CALL YMCPIR ( n, cov, p, iwork(piwo), dwork(pdcov), info )
      dwork(pdcov + p*p - 1) = 1.E-12
c
c     initialize linear part to 0
      CALL IVX ( p, dwork(pdrho), ZERO )
c
c     construction of the constraints matrix and vector
      IF ((neq+nin) .GT. 0) THEN
        DO i = 1, (neq+nin)
            dwork(pdwo + i - 1) = coef*bcst(i)
        ENDDO
      ENDIF
      CALL CTSR ( n, rho, rfr, neq, nin, ccst, dwork(pdwo),
     &            dwork(pdcinf), dwork(pdcsup),
     &            neqtot, nintot, dwork(pdc), dwork(pdb))
c
c     initialize lower/upper bounds
      CALL IVX ( p, dwork(pdcinf), ZERO )
      CALL IVX ( p, dwork(pdcsup), INF  )
      dwork(pdcsup + p - 1) = n*INF
c
c     quadratic solver
      CALL qp ( p, dwork(pdcov), dwork(pdrho),
     &          neqtot, nintot, dwork(pdc), dwork(pdb),
     &          dwork(pdcinf), dwork(pdcsup),
     &          iwork(piwo), dwork(pdwo),
     &          dwork(pdlagr), dwork(pdx), info)
      IF (info .NE. 0) THEN
        CALL IVX ( n, wopt, ZERO )
        RETURN
      ENDIF
c
c     re-scale weights
      t = dwork(pdx + p - 1)
c
c     if |t| < EPS then (no solution)
c     return min. vol. solution
      IF (ABS(t) .LT. 100*EPS) THEN
c        info = -113
         DO i = 1,n
            wopt(i) = (1.0 - wrf)*wopt(i)
        ENDDO
        RETURN
      ENDIF
      sum = 0.0
      DO i = 1,n
        wopt(i) = dwork(pdx + i - 1)/t
        sum = sum + wopt(i)
      ENDDO
      DO i = 1,n
        wopt(i) = wopt(i)/sum
      ENDDO
c
c     re-scale result
      DO i = 1,n
        wopt(i) = (1.0 - wrf)*wopt(i)
      ENDDO
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
! Time of File Save by ERB: 6/6/2006 4:22PM
C     Last change:  ERB   8 Nov 2001    3:06 pm
      SUBROUTINE SEN1DRT1FM(MXDRT,DRTF,HNEW,NCOL,NROW,NLAY,IBOUND,RHS,
     &                      IP,NDRTVL,IDRTFL)
C-----VERSION 20060606 ERB
C     ******************************************************************
C     FOR DRAIN-RETURN AND RECIPIENT CELLS: CALCULATE MATRIX AND VECTOR
C     DERIVATIVES, MULTIPLY BY HEADS, AND ADD COMPONENTS TO RHS.
C     ******************************************************************
C     Modified 11/8/2001 to support parameter instances - ERB
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL DRTF, RHS
      INTEGER I, IBOUND, II, J, K, MXDRT, NCOL, NLAY, NROW
      DOUBLE PRECISION DF, HH, HNEW(NCOL,NROW,NLAY)
      DIMENSION DRTF(NDRTVL,MXDRT), IBOUND(NCOL,NROW,NLAY),
     &          RHS(NCOL,NROW,NLAY)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
C
C     FIND INSTANCE NUMBER
      NI = IACTIVE(IP)
      IF (NI.EQ.0) RETURN
C
C     DETERMINE RANGE OF ELEMENTS CURRENTLY ACTIVE FOR THIS PARAMETER
      IPOS1 = IPLOC(1,IP)
      NPC = IPLOC(2,IP)-IPOS1+1
      NUMINST = IPLOC(3,IP)
      IF (NUMINST.GT.1) NPC = NPC/NUMINST
C-----LOOP THROUGH PARAMETER CELLS
      DO 20 II = 1, NPC
        ICP = IPOS1-1+II+(NI-1)*NPC
        K = DRTF(1,ICP)
        I = DRTF(2,ICP)
        J = DRTF(3,ICP)
        IF (IBOUND(J,I,K).GT.0) THEN
C---------CALCULATE CONTRIBUTION TO SENSITIVITY
          ELEVD = DRTF(4,ICP)
          HH = HNEW(J,I,K)
          IF (HH.GT.ELEVD) THEN
            FACTOR = DRTF(5,ICP)
            DF = FACTOR*(ELEVD-HH)
            RHS(J,I,K) = RHS(J,I,K) - DF
C  CALCULATE CONTRIBUTION TO SENSITIVITY OF RETURN FLOW
            IF (IDRTFL.GT.0) THEN
              KR = DRTF(6,ICP)
              IF (KR.NE.0) THEN
                IR = DRTF(7,ICP)
                JR = DRTF(8,ICP)
                IF (IBOUND(JR,IR,KR) .GT. 0) THEN
                  RFPROP = DRTF(9,ICP)
                  DFR = RFPROP*FACTOR*(HH-ELEVD)
                  RHS(JR,IR,KR) = RHS(JR,IR,KR) - DFR
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDIF
   20 CONTINUE
C
      RETURN
      END

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
The UCDMC or UCD Medical Center Bus Line stops on campus behind Haring Hall across from the Silo (Lot 29), out by the Health Sciences Complex (Lot 56), and in front of Mondavi Center   then its off to the UCDMC!

Riding this transportation option is limited to university affiliates (faculty, students, staff, interns, volunteers, patients, etc...)

An unlimited pass is issued to you if you volunteer at the medical center. Otherwise, you will need to buy a ticket (drivers do not take cash). Tickets are $1.50 each way and can be purchased with a campus I.D. at the following locations:

  UCDMC Parking and Transportation Services located in the Facilities Services Building; 7:30 a.m.  5:00 p.m.
  UCDMC Parking lot 4 kiosk, south of the Davis Tower; 7:30 a.m.  5:30 p.m.
  UC Davis Transportation and Parking Services Office on Extension Center Drive; 7:00 a.m.  5:00 p.m.
  UC Davis Cashiers Office in 1200 Dutton Hall (first floor); 9:00 a.m.  4:00 p.m.

Commuters can also purchase a 10ride ticket for $15, or an unlimited monthly pass for $45.

It may be an alternative to Yolobus for some  its 2030 minutes into Sacramento vs. Yolobuss 60 minutes.

During rush hour in the morning and evening, expect to arrive 1530 minutes later going either direction than what the schedule below says!  Also, if the bus fills up with riders you may have to remain standing for the ride  although they have been known to provide another shuttle if there are many extra passengers.

The buses are equipped with bike racks for bringing your bicycle with you across the Yolo Causeway Causeway.  Each bus only has space for two bicycles  first come, first served.

 UCD to UCDMC

 UCDMC to UCD

 Service for Mondavi Center only.