text:
MODULE mympi
  !
  ! Description:
  !   Module for mpi to help associating data to multiple processors
  !
  ! Current Code Owner: [Name of person responsible for this code]
  !
  ! History:
  !
  ! Version   Date     Comment
  ! -------   ----     -------
  ! 1.0       11-11-08 Original code. Thomas Hauser
  !
  ! Code Description:
  !   Language:           Fortran 90.
  !   Software Standards: "European Standards for Writing and
  !     Documenting Exchangeable Fortran 90 Code".
  !

  IMPLICIT NONE

CONTAINS

  INTEGER FUNCTION block_owner(j,p,n)
    INTEGER, INTENT(in) :: j,p,n
    block_owner = ((p*(j+1)-1)/n)
    RETURN
  END FUNCTION block_owner

  INTEGER FUNCTION block_low(id,p,n)
    INTEGER, INTENT(in) :: id,p,n
    block_low =  (id*n)/p
    RETURN
  END FUNCTION block_low

  INTEGER FUNCTION block_high(id,p,n)
    INTEGER, INTENT(in) :: id,p,n
    block_high = ((id+1)*n)/p-1
    RETURN
  END FUNCTION block_high

  INTEGER FUNCTION block_size(id,p,n)
    INTEGER, INTENT(in) :: id,p,n
    block_size = (block_low(id+1,p,n)-block_low(id,p,n))
    RETURN
  END FUNCTION block_size

END MODULE mympi

max_stars_repo_name: rctraining/USGS-Advanced-HPC-Workshop

max_issues_repo_name: rctraining/USGS-Advanced-HPC-Workshop

max_forks_repo_name: rctraining/USGS-Advanced-HPC-Workshop


text:
Gary Sandy was Mayor of wiki:Woodland. He currently serves as the Director of the UC Davis Office of Local Government Relations. In that capacity he serves as the Chancellors representative to the City  UCD Student Liaison Commission.

In 2009, Gary was elected a Trustee of the Yuba Community College District, representing Area 3.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE ESTIM(J)
!                                                                      *
!   This  subprogram implements  Part 1 of Algorithm 7.1 of C Froese   *
!   Fischer, Comput Phys Rep, 3 (1986) 320-321.                        *
!                                                                      *
!   Written by Farid A Parpia, at Oxford    Last update: 26 Sep 1993   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:06:32   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE def_C
      USE grid_C
      USE orb_C
      USE scf_C
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: J
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: NPJ, NAKABS
      REAL(DOUBLE) :: ALPHA, CSQ, FNREL, FKABS, FKAP2, ZALPHA, GAMMA, EBYM
!-----------------------------------------------
!
!   Initializations
!
      ALPHA = 1.0D00/C
      CSQ = C*C
      NPJ = NP(J)
      NAKABS = ABS(NAK(J))
      FNREL = DBLE(NPJ - NAKABS)
      FKABS = DBLE(NAKABS)
      FKAP2 = FKABS*FKABS
!
!   Set ZINF, the asymptotic charge seen by the electron
!
!     ZINF = DBLE (IONCTY+1)
!
!   Changed on 07/06/93 by WPW
!
      ZINF = Z + DBLE((-NELEC) + 1)
!
!   Set the lower bound
!
      ZALPHA = ZINF*ALPHA
      IF (ZALPHA < FKABS) THEN
         GAMMA = SQRT(FKAP2 - ZALPHA*ZALPHA)
         EBYM = 1.0D00/SQRT(1.0D00 + (ZALPHA/(GAMMA + FNREL + 0.5D00))**2)
         EPSMIN = (1.0D00 - EBYM)*CSQ
      ELSE
         EPSMIN = 0.25D00*CSQ/DBLE(NPJ*NPJ)
      ENDIF
      EMIN = EPSMIN
!
!   Set the upper bound
!
      ZALPHA = Z*ALPHA
      IF (ZALPHA < FKABS) THEN
         GAMMA = SQRT(FKAP2 - ZALPHA*ZALPHA)
         EBYM = 1.0D00/SQRT(1.0D00 + (ZALPHA/(GAMMA + FNREL - 0.5D00))**2)
         EPSMAX = (1.0D00 - EBYM)*CSQ
      ELSE
         EPSMAX = CSQ + CSQ
      ENDIF
      EMAX = EPSMAX
!
      RETURN
      END SUBROUTINE ESTIM

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
Describe RobertFosselman here.

This is my second year living in Davis and I really enjoy it.  My main hobby is running, ultra running to be more specific.  I am not fast so I like to go the distance.  My big race I am currently training for is Western States, it is a 100 mile endurance race in June.  I recently got a car which is kind of cool so now I can go up to Auburn to train.  Although around town I am usually on my old Nishiki Superb, although I did just recently become the owner of a very cool early 90s Specialized Epic so that may become my new ride for awhile.

I am not an environmentalist but I am not a wasteful person.  I get very irritated when people and the city water the streets.  That realllllllllly gets me mad.  I also would appreciate it if people wouldnt put their Green Waste clippings on the street.  I get really irritated when I have to move from the bike lane to the street just so I am not riding my bike over a bunch of clippings or tree branches having to worry about popping my tubes.


text:
! =======================================================================
! WEIGHTS   Compute coefficients for finite difference approximation for
!           the derivatives 1 to m at point z assuming data is known at
!           points in array x. Based on the program "weights" in
!           B. Fornberg, "Calculation of weights in finite difference
!           formulas", SIAM Review 40 (1998), pp. 685-691.
! Arguments
!
!   z   (input) DOUBLE
!       location where approximations are to be accurate
!
!   x   (input) DOUBLE Array
!       array containing interpolation points
!
!   m   (input) INTEGER
!       highest derivative for which weights are sought
!
!   W   (output) DOUBLE array, dimension(size(x),m+1)
!       matrix that gives weights at grid locations x for
!       derivative of order j<=m are found in c(:,j)
! =======================================================================

subroutine weights_z(z,x,m,W)

    ! Arguments
    complex(dp), intent(in)    :: z
    complex(dp), intent(in)    :: x(:)
    integer,  intent(in)       :: m
    complex(dp), intent(out)   :: W(size(x),m+1)

    ! Variable Declarations
    complex(dp) :: c1, c2, c3, c4, c5
    integer  :: i,j,k,n,mn

    c1 = 1.0_dp
    c4 = x(1) - z
    W  = 0.0_dp
    W(1,1) = 1.0_dp

    n = size(x)
    do i=2,n
        mn = min(i,m+1)
        c2 = 1.0_dp
        c5 = c4
        c4 = x(i) - z
        do j=1,i-1
            c3 = x(i) - x(j)
            c2 = c2*c3;
            if(j == i-1) then
                do k=mn,2,-1
                    W(i,k) = c1*(real(k-1,dp)*W(i-1,k-1) - c5*W(i-1,k))/c2;
                enddo
                W(i,1) = -c1*c5*W(i-1,1)/c2;
            endif
            do k=mn,2,-1
                W(j,k) = (c4*W(j,k) - real(k-1,dp)*W(j,k-1))/c3;
            enddo
            W(j,1) = c4*W(j,1)/c3;
        enddo
        c1 = c2;
    enddo

end subroutine weights_z

max_stars_repo_name: buvoli/epbm

max_issues_repo_name: buvoli/epbm

max_forks_repo_name: buvoli/epbm


text:
High Density Networks (HDN) is a High Speed Wireless Internet provider that is offered at various housing communities in the Davis area.  It also supplies DIRECTV Packages starting at $29.99.

They charge $24.95/month for wifi access in your housing complex.
Speeds are claimed as being up to 1.5 Mbps.

The service is plagued by outages, lack of connectivity, stability issues, and faulty/antiquated equipment.

What they dont tell you is that your money actually goes to your apartment complex,  as the complex pays them a flat fee for the service.  So they really dont care how many people use it or the level of customer satisfaction.
Calling customer service is generally met with a voicemail recording which seemingly is never checked, and their response rates are only slightly better via email.  Should you ever reach a live person, they refuse to accept that the problem lies in their network, and at best will promise to send a technician out sometime.  Dont expect any promptness, nor improvement.

Heres a preview of the various error messages you will regularly be subjected to:

 ||Image(TimeoutError.JPG, thumbnail)||Image(AcctError.JPG, thumbnail)||  Image(AddError.JPG, thumbnail)||


module view_m
#ifdef dnad
    use dnadmod
#define real type(dual)
#endif
    use plane_m
    implicit none

contains

!-----------------------------------------------------------------------------------------------------------
subroutine view_plotmtv(t)
    type(plane_t) :: t
    type(section_t),pointer :: si
    character(100) :: filename
    character(140) :: command
    integer :: i,iwing,isec,gnum,iforce,ipoint,min_isec
    integer :: ierror = 0
    real :: P0(3),P1(3),P2(3),gpsize,delta,force_pos(3),force_dir(3),force_mag,dist,min_dist

    filename = trim(adjustl(t%master_filename))//'_PlotData.txt'
    open(unit = 10, File = filename, action = 'write', iostat = ierror)

    write(10,*) ' $ DATA=CURVE3D'
    write(10,*) ' % toplabel   = "Planform"'
    write(10,*) ' % xlabel     = "-x"'
    write(10,*) ' % ylabel     = "-y"'
    write(10,*) ' % zlabel     = "-z"'
    write(10,*) ' % grid       = False'
!    write(10,*) ' % axislabel  = False'
    write(10,*) ' % equalscale = True'
    write(10,*) ' % axisscale  = False'
!    write(10,*) ' % axisguides = False'
    write(10,*) ' % eyepos.x   = 0.50'
    write(10,*) ' % eyepos.y   = 0.75'
    write(10,*) ' % eyepos.z   = 0.25'
!    write(10,*) ' % xmin   = -0.60'
!    write(10,*) ' % xmax   =  0.60'
!    write(10,*) ' % ymin   = -0.70'
!    write(10,*) ' % ymax   =  0.2'
    write(10,*) ' % dlinecolor = 6'
    write(10,*) ' '

    do iwing=1,t%nrealwings !real wings
        do isec=1,t%wings(iwing)%nSec
            si => t%wings(iwing)%sec(isec)
            write(10,*) ' % linecolor = 5' !Quarter Chord
            write(10,*) -si%P1(:)
            write(10,*) -si%P2(:)
            write(10,*) ' '

            write(10,*) ' % linecolor = 5' !leading edge
            write(10,*) -si%P1(:)
            write(10,*) -(si%P1(:) - 0.25*si%chord_1*si%ua(:))
            write(10,*) -(si%P2(:) - 0.25*si%chord_2*si%ua(:))
            write(10,*) -si%P2(:)
            write(10,*) ' '

            write(10,*) ' % linecolor = 5' !trailing edge
            write(10,*) -si%P1(:)
            write(10,*) -(si%P1(:) + (0.75-si%cf_c)*si%chord_1*si%ua(:))
            write(10,*) -(si%P2(:) + (0.75-si%cf_c)*si%chord_2*si%ua(:))
            write(10,*) -si%P2(:)
            write(10,*) ' '

            if(si%cf_c > 0.0) then
                write(10,*) ' % linecolor = 7' !flaps
                P1 = -(si%P1(:) + (0.75-si%cf_c)*si%chord_1*si%ua(:))
                P2 = -(si%P2(:) + (0.75-si%cf_c)*si%chord_2*si%ua(:))
                write(10,*) P1
                write(10,*) P1 - si%cf_c*si%chord_1*si%uf(:)
                write(10,*) P2 - si%cf_c*si%chord_2*si%uf(:)
                write(10,*) P2
                write(10,*) ' '
            end if
        end do
    end do

    do iwing=t%nrealwings+1,t%nwings !reflected wings
        do isec=1,t%wings(iwing)%nSec
            si => t%wings(iwing)%sec(isec)
            write(10,*) ' % linecolor = 3' !Quarter Chord
            write(10,*) -si%P1(:)
            write(10,*) -si%P2(:)
            write(10,*) ' '

            write(10,*) ' % linecolor = 3' !leading edge
            write(10,*) -si%P1(:)
            write(10,*) -(si%P1(:) - 0.25*si%chord_1*si%ua(:))
            write(10,*) -(si%P2(:) - 0.25*si%chord_2*si%ua(:))
            write(10,*) -si%P2(:)
            write(10,*) ' '

            write(10,*) ' % linecolor = 3' !trailing edge
            write(10,*) -si%P1(:)
            write(10,*) -(si%P1(:) + (0.75-si%cf_c)*si%chord_1*si%ua(:))
            write(10,*) -(si%P2(:) + (0.75-si%cf_c)*si%chord_2*si%ua(:))
            write(10,*) -si%P2(:)
            write(10,*) ' '

            if(si%cf_c > 0.0) then
                write(10,*) ' % linecolor = 9' !flaps
                P1 = -(si%P1(:) + (0.75-si%cf_c)*si%chord_1*si%ua(:))
                P2 = -(si%P2(:) + (0.75-si%cf_c)*si%chord_2*si%ua(:))
                write(10,*) P1
                write(10,*) P1 - si%cf_c*si%chord_1*si%uf(:)
                write(10,*) P2 - si%cf_c*si%chord_2*si%uf(:)
                write(10,*) P2
                write(10,*) ' '
            end if
        end do
    end do

    if(t%Gammas(1) .ne. 0.0) then
        do i=1,t%nSize ! %normal
            si => t%sec(i)%myp
            write(10,*) ' % linecolor = 6'
            write(10,*) -si%PC(:)
!            write(10,*) -(si%PC(:) + 10.0*sec_CL(si)*si%un(:))
            write(10,*) -(si%PC(:) + 10.0*t%Gammas(i)*si%un(:))
            write(10,*) ' '
        end do
    end if

    if(t%groundplane .eq. 1) then !draw a ground plane
        gpsize = 10.0 !length of ground plane
        gnum = 5 !number of squares on plane
        delta = gpsize/REAL(gnum)

        P0(1) = t%CG(1) - t%hag*sin(t%alpha) !offset from CG !this code copied from wing.f90
        P0(2) = t%CG(2)
        P0(3) = t%CG(3) + t%hag*cos(t%alpha)

        do i=-gnum,gnum
            write(10,*) ' % linecolor = 8'
            P1(1) = P0(1) - gpsize*cos(t%alpha); P1(2) = P0(2) + REAL(i)*delta; P1(3) = P0(3) - gpsize*sin(t%alpha)
            P2(1) = P0(1) + gpsize*cos(t%alpha); P2(2) = P0(2) + REAL(i)*delta; P2(3) = P0(3) + gpsize*sin(t%alpha)
            write(10,*) -P1(:)
            write(10,*) -P2(:)
            write(10,*) ' '
        end do
        do i=-gnum,gnum
            write(10,*) ' % linecolor = 8'
            P1(1) = P0(1) - REAL(i)*delta*cos(t%alpha); P1(2) = P0(2) + gpsize; P1(3) = P0(3) - REAL(i)*delta*sin(t%alpha)
            P2(1) = P0(1) - REAL(i)*delta*cos(t%alpha); P2(2) = P0(2) - gpsize; P2(3) = P0(3) - REAL(i)*delta*sin(t%alpha)
            write(10,*) -P1(:)
            write(10,*) -P2(:)
            write(10,*) ' '
        end do
    end if

    !show external forces tied to closest structure
    do iforce = 1,size(t%external_forces)
        do ipoint=1,t%external_forces(iforce)%datasize
            force_pos(:) = t%external_forces(iforce)%RawData(ipoint,1:3)
            force_dir(:) = t%external_forces(iforce)%RawData(ipoint,4:6)
            force_mag    = t%external_forces(iforce)%RawData(ipoint,7)
            min_dist = 1.0e16
            do isec=1,t%nSize
                si => t%sec(isec)%myp
                dist = math_length(3,force_pos(:),si%PC(:))
                if(dist < min_dist) then
                    min_dist = dist
                    min_isec = isec
                end if
            end do
            si => t%sec(min_isec)%myp
            write(10,*) ' % linecolor = 9'
            write(10,*) -si%PC(:)
            write(10,*) -force_pos(:)
            write(10,*)
            write(10,*) ' % linecolor = 10'
            write(10,*) -force_pos(:)
            write(10,*) -(force_pos(:) + force_dir(:))
            write(10,*)
        end do
    end do

    close(10)
    command = '/sw/bin/Plotmtv -fg white -bg black '//filename
    call system(command)
    command = 'Plotmtv -fg white -bg black '//filename
    call system(command)
end subroutine view_plotmtv

!-----------------------------------------------------------------------------------------------------------

subroutine view_vtk(t)
    type(plane_t) :: t
    type(section_t),pointer :: si
    character(100) :: filename,str1
    integer :: iwing,isec
    integer :: ierror = 0
    120 Format(3ES25.13)
    130 Format(3I5)

    filename = trim(adjustl(t%master_filename))//'_view.vtk'
    open(unit = 10, File = filename, action = 'write', iostat = ierror)

    write(10,'(A)') '# vtk DataFile Version 3.0'
    write(10,'(A)') 'vtk output'
    write(10,'(A)') 'ASCII'
    write(10,'(A)') 'DATASET POLYDATA'
    write(10,'(A)')
    write(10,'(A,I5,A)') 'POINTS ',t%nSize*2,' float'

    do iwing=1,t%nrealwings !real wings
        do isec=1,t%wings(iwing)%nSec
            si => t%wings(iwing)%sec(isec)
            write(str1,120) si%PC(:); write(10,'(A)') trim(adjustl(str1))
            write(str1,120) si%PC(:) + (0.75-si%cf_c)*si%chord_c*si%ua(:); write(10,'(A)') trim(adjustl(str1))
        end do
    end do

    write(10,*)
    write(10,'(A,I5,I5)') 'LINES ',t%nSize,t%nSize*3
    do iwing=0,t%nSize-1
        write(str1,130) 2,2*iwing,2*iwing+1; write(10,'(A)') trim(adjustl(str1))
    end do

    close(10)
end subroutine view_vtk

!-----------------------------------------------------------------------------------------------------------

subroutine view_stl(t)
    type(plane_t) :: t
    type(section_t),pointer :: si
    type(airfoil_t),pointer :: af1
    type(airfoil_t),pointer :: af2
    real,allocatable,dimension(:,:) :: af_points1,af_points2
    character(100) :: filename
    integer :: i,iwing,isec,af_datasize
    integer :: ierror = 0

    do i=1,size(airfoils)
        call af_create_geom_from_file(airfoils(i),DB_Airfoil)
    end do

    filename = trim(adjustl(t%master_filename))//'_view.stl'
    open(unit = 10, File = filename, action = 'write', iostat = ierror)
!    filename = trim(adjustl(t%master_filename))//'_view.web'
!    open(unit = 20, File = filename, action = 'write', iostat = ierror)

    write(10,'(A)') 'solid geom'

    do iwing=1,t%nrealwings !real wings

        af_datasize = t%wings(iwing)%airfoils(1)%p%geom%datasize
        do i=1,t%wings(iwing)%nairfoils
            if(t%wings(iwing)%airfoils(i)%p%geom%datasize .ne. af_datasize) then
                write(*,*) 'All airfoils for wing ',t%wings(iwing)%name,' must have same number of nodes.'
                stop
            end if
        end do

        allocate(af_points1(af_datasize,3))
        allocate(af_points2(af_datasize,3))

        if(t%wings(iwing)%is_linear.ne.1) then
            do isec=1,t%wings(iwing)%nSec
                si => t%wings(iwing)%sec(isec)
                call view_create_local_airfoil(si%af1_a,si%af1_b,t%wings(iwing)%side,si%af_weight_1,si%chord_1,&
                                            & si%twist1,si%dihedral1,af_datasize,si%P1(:),af_points1)
                call view_create_local_airfoil(si%af2_a,si%af2_b,t%wings(iwing)%side,si%af_weight_2,si%chord_2,&
                                            & si%twist2,si%dihedral2,af_datasize,si%P2(:),af_points2)
                call view_create_stl_shell(af_datasize,af_points1,af_points2)
            end do
        end if

        if(t%wings(iwing)%side == 'right') then
            si => t%wings(iwing)%sec(1)
            call view_create_local_airfoil(si%af1_a,si%af1_b,t%wings(iwing)%side,si%af_weight_1,si%chord_1,&
            & si%twist1,si%dihedral1,af_datasize,si%P1(:),af_points1)
            si => t%wings(iwing)%sec(t%wings(iwing)%nSec)
            call view_create_local_airfoil(si%af2_a,si%af2_b,t%wings(iwing)%side,si%af_weight_2,si%chord_2,&
            & si%twist2,si%dihedral2,af_datasize,si%P2(:),af_points2)
        else
            si => t%wings(iwing)%sec(1)
            call view_create_local_airfoil(si%af2_a,si%af2_b,t%wings(iwing)%side,si%af_weight_2,si%chord_2,&
                                        & si%twist2,si%dihedral2,af_datasize,si%P2(:),af_points2)
            si => t%wings(iwing)%sec(t%wings(iwing)%nSec)
            call view_create_local_airfoil(si%af1_a,si%af1_b,t%wings(iwing)%side,si%af_weight_1,si%chord_1,&
                                        & si%twist1,si%dihedral1,af_datasize,si%P1(:),af_points1)
        end if

        if(t%wings(iwing)%is_linear.eq.1) call view_create_stl_shell(af_datasize,af_points1,af_points2)
!        call view_create_stl_rib(af_datasize,af_points1)
!        call view_create_stl_rib(af_datasize,af_points2)
        deallocate(af_points1)
        deallocate(af_points2)
    end do

    write(10,*) 'endsolid geom'
    close(10)
!    close(20)
end subroutine view_stl

!-----------------------------------------------------------------------------------------------------------
subroutine view_create_stl_shell(datasize,af_points1,af_points2)
    integer :: datasize,i
    real :: af_points1(datasize,3),af_points2(datasize,3)
    real :: P1(3),P2(3),P3(3)

    !Outer surface
    do i=1,datasize-1
        P1(:) = af_points1(i,:)
        P2(:) = af_points2(i,:)
        P3(:) = af_points1(i+1,:)

        call view_add_stl_triangle(P1,P3,P2)

        P1 = P2
        P2(:) = af_points2(i+1,:)
        call view_add_stl_triangle(P1,P3,P2)
    end do
    !Close Surface
    i = datasize
    P1(:) = af_points1(i,:)
    P2(:) = af_points1(1,:)
    P3(:) = af_points2(i,:)

    call view_add_stl_triangle(P1,P3,P2)

    P1 = P2
    P2(:) = af_points2(1,:)
    call view_add_stl_triangle(P1,P3,P2)
end subroutine view_create_stl_shell

!-----------------------------------------------------------------------------------------------------------
subroutine view_create_stl_rib(datasize,af_points)
    integer :: datasize,i
    real :: af_points(datasize,3)
    real :: P1(3),P2(3),P3(3)

    do i=1,(datasize-2)/2
        P1(:) = af_points(i,:)
        P2(:) = af_points(datasize-i+1,:)
        P3(:) = af_points(i+1,:)
        call view_add_stl_triangle(P1,P3,P2)

        P1 = P2
        P2(:) = af_points(datasize-i,:)
        call view_add_stl_triangle(P1,P3,P2)
    end do
end subroutine view_create_stl_rib

!-----------------------------------------------------------------------------------------------------------
subroutine view_add_stl_triangle(P1,P2,P3)
    real :: P1(3),P2(3),P3(3),norm(3)
    110 Format(A15, 3ES25.13)

    call math_plane_normal(P1,P2,P3,norm)

#ifndef dnad
    write(10,110) 'facet normal ',norm(:)
    write(10,*) 'outer loop'
    write(10,110) 'vertex ',P1(:)
    write(10,110) 'vertex ',P2(:)
    write(10,110) 'vertex ',P3(:)
    write(10,*) 'endloop'
    write(10,*) 'endfacet'
#else
    write(10,110) 'facet normal ',norm(:)%x
    write(10,*) 'outer loop'
    write(10,110) 'vertex ',P1(:)%x
    write(10,110) 'vertex ',P2(:)%x
    write(10,110) 'vertex ',P3(:)%x
    write(10,*) 'endloop'
    write(10,*) 'endfacet'
#endif

!    write(20,120) P1(:),P2(:),P3(:)
end subroutine view_add_stl_triangle

!-----------------------------------------------------------------------------------------------------------
subroutine view_create_local_airfoil(af1,af2,side,percent,chord,twist,dihedral,datasize,point,output)
    type(airfoil_t),pointer :: af1
    type(airfoil_t),pointer :: af2
    character(5) :: side
    real :: percent,chord,twist,dihedral
    integer :: datasize,i
    real :: point(3),output(datasize,3)

    output(:,1:2) = af1%geom%RawData(:,:) + percent*(af2%geom%RawData(:,:) - af1%geom%RawData(:,:))
    output(:,1) = output(:,1) - 0.25
    output = chord*output
    output(:,1) = -output(:,1) !flip x
    output(:,3) = -output(:,2) !assign y to z
    output(:,2) = 0.0 !set y to zero

    if(side.eq.'left') then
!        twist = -twist
        dihedral = -dihedral
    end if

    do i=1,datasize
        call math_rot_y(output(i,:),twist)
        call math_rot_x(output(i,:),-dihedral)
        output(i,:) = output(i,:) + point(:)
    end do

end subroutine view_create_local_airfoil

end module view_m


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE BRINT1 (IA,IB,IC,ID,K,TEGRAL)
!                                                                      *
!   Returns integrals for the transverse photon interaction.           *
!                                                                      *
!   Written by Per Jonsson                   Octaober 2014             *
!                                                                      *
!***********************************************************************
!...Translated by Charlotte Froese Fischer
!                       Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE bilst_C
      USE orb_C,           ONLY: NW
      USE kkstartbreit_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: ia, ib, ic, id, k
      REAL(DOUBLE), INTENT(out) :: tegral
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER ::  index, loc, ju, jl, jm, key
!-----------------------------------------------
!
      KEY = NW + 1
      INDEX = ((IA*KEY+IB)*KEY+IC)*KEY+ID
!
      JL = KSTARTBREIT1(K)
      JU = KSTARTBREIT1(K+1) - 1
!
      IF (INDEX < INDTP1(JL).OR.INDEX > INDTP1(JU)) THEN
        WRITE(*,*) 'Something wrong in brint1'
        STOP
      ENDIF
!
!   The index is within the range of the indices stored; search
!   for it in the list of indices
!
    1 IF (JU-JL .GT. 1) THEN
         JM = (JU+JL)/2
         IF (INDTP1(JM) > INDEX) THEN
            JU = JM
         ELSE
            JL = JM
         ENDIF
         GOTO 1
      ENDIF
!
!   The range is bracketed to the extent possible
!
      IF (INDEX .EQ. INDTP1(JU)) THEN
         LOC = JU
      ELSEIF (INDEX .EQ. INDTP1(JL)) THEN
         LOC = JL
      ELSE
         WRITE(*,*) K,IA,IB,IC,ID,INDEX
         WRITE(*,*) 'Brint1 Integral not found'
         STOP
      ENDIF
!
!   Return the value of the integral
!   from storage

      TEGRAL = VALTP1(LOC)
!
      RETURN
      END SUBROUTINE BRINT1

max_issues_repo_name: compas/grasp-rci-qed

max_forks_repo_name: compas/grasp-rci-qed


text:
      SUBROUTINE QUIKVIS3
      IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE IS PART OF THE QUIKVIS PROGRAM.  IT INITIALIZES ALL
C QUIKVIS-UNIQUE ROUTINES NOT USED FOR PROGRAM INITIALIZATION OR
C USER INPUT.
C
C ROUTINES ARE INITIALIZED WHEN QVINIT'S LOAD THE GLOBAL VARIABLE VALUES
C INTO THEIR OWN LOCAL VARIABLES.
C
C QVINIT4, FOR EXAMPLE, IS AN ENTRY POINT IN QUIKVIS4.
C
C***********************************************************************
C
C BY C PETRUZZO/GFSC/742.   2/86.
C       MODIFIED.... 9/86. CJP. ADDED CALLS TO FOR NEWLY WRITTEN ENTRIES
C                               QVINIT5A2C3, ..5A3, ..5U, ..5V, AND ..5W
C
C
C***********************************************************************
C
      CALL QVINIT4
C
      CALL QVINIT5
        CALL QVINIT5A
          CALL QVINIT5A1
          CALL QVINIT5A2
            CALL QVINIT5A2A
            CALL QVINIT5A2B
            CALL QVINIT5A2C
              CALL QVINIT5A2C1
              CALL QVINIT5A2C2
              CALL QVINIT5A2C3
            CALL QVINIT5A2D
              CALL QVINIT5A2D1
          CALL QVINIT5A3
        CALL QVINIT5B
          CALL QVINIT5B1
            CALL QVINIT5B1A
            CALL QVINIT5B1B
            CALL QVINIT5B1X
        CALL QVINIT5U
        CALL QVINIT5V
        CALL QVINIT5W
        CALL QVINIT5X
        CALL QVINIT5Y
        CALL QVINIT5Z
C
      RETURN
      END

ext: for


text:
      MODULE m_spratm
c------------------------------------------------------------------
c
c    This is the driver subroutine for a full-relativistic spin-polarized
c    core charge and spin density calculation using the 4 fold coupled
c    dirac equation : collection of references on this story see in
c    H. Ebert, J.Phys.: Condens. Matter 1 (1989) 9111.
c
c    Attention : The algorithm uses Ry-units. Therefore , potential
c                is multiplied by factor of 2 before it is used in the
c                dirac equation
c
c---> input:
c    Vr      =   spherical potential
c    Br      =   spherical magn. field
c    z       =   atomic charge
c    rnot    =   radial mesh starting point
c    dx      =   radial mesh logariphmic increment
c    jtop    =   upper bond for core radial mesh
c---> i/o
c    ectab   =   atomic energy levels for (\kappa,\mu) (in Hr)
c---> output
c    sume    =   sum of atomic eigenvalues (in Hr)
c    rhochr  =   core charge density
c    rhospn  =   core spin density
c
c........................................................ spratm
      CONTAINS
      SUBROUTINE spratm(
     >                  msh,vr,br,z,rnot,dx,jtop,ectab,ntab,ltab,
     <                  sume,rhochr,rhospn)
c
      USE m_core
      IMPLICIT NONE
c
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: msh,jtop
      REAL,    INTENT (IN) :: dx,rnot,z
      REAL,    INTENT (OUT):: sume
C     ..
C     .. Array Arguments ..
      INTEGER, INTENT (IN) :: ntab(100),ltab(100)
      REAL,    INTENT (IN) :: br(msh),vr(msh)
      REAL,    INTENT (OUT):: rhochr(msh),rhospn(msh)
      REAL,    INTENT (INOUT):: ectab(100)
C     ..
C     .. Local Scalars ..
      REAL rr,stval
      INTEGER ic,ir,nshell,n_old,l_old
C     ..
C     .. Local Arrays ..
      REAL bt(msh),vt(msh)
      INTEGER nqntab(15),lqntab(15)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC exp,log
C     ..
      nshell = 0
      ic = 0 ; n_old = -1 ; l_old = -1
      DO WHILE (ntab(ic+1).GT.0)
        ic = ic + 1
        IF  (ntab(ic).NE.n_old) THEN
           nshell = nshell + 1
           nqntab(nshell) = ntab(ic)
           lqntab(nshell) = ltab(ic)
           n_old = ntab(ic)
           l_old = ltab(ic)
        ELSEIF (ltab(ic).NE.l_old) THEN
           nshell = nshell + 1
           nqntab(nshell) = ntab(ic)
           lqntab(nshell) = ltab(ic)
           n_old = ntab(ic)
           l_old = ltab(ic)
        ENDIF
      ENDDO
c Hr -> Ry
      ic = 0
      DO ic = 1,100
         ectab(ic) = 2.*ectab(ic)
      END DO
c potential and field redefinition
      rr = rnot
      DO ir = 1,msh
         vt(ir) = 2.*vr(ir)/rr
         bt(ir) = 2.*br(ir)/rr
         rr = rr*exp(dx)
      END DO
      stval = log(rnot)
c
      CALL core(
     >          msh,vt,bt,z,stval,dx,nshell,nqntab,lqntab,jtop,
     X          ectab,
     <          rhochr,rhospn)

c Ry -> Hr
      sume = 0.0
      DO ic = 1,100
         ectab(ic) = ectab(ic)/2.
         sume = sume + ectab(ic)
      END DO

      END SUBROUTINE spratm
      END MODULE m_spratm

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
! Depending on the OpenCL Platform selected at compile time (see Makefile)
! a FILENAME with the appropriate string value is chosen.
! This can be an OpenCL source file or a compiled binary file.
#ifdef JIT
#define FILENAME "vadddevice.cl"
#elif INTEL_FPGA
#define FILENAME "vadddevice.aocx"
#elif XILINX_FPGA
#define FILENAME "vadddevice.xclbin"
#endif

program vadd
  use fortcl
  use clfortran
  use ocl_utils_mod, only: check_status
  use mpi
  use iso_c_binding
  implicit none

  integer, parameter :: N = 100
  integer(c_intptr_t), target :: deviceA, deviceB, deviceC
  integer(c_intptr_t), target :: vaddkernel
  integer(c_intptr_t), pointer :: cmd_queues(:)
  integer(c_size_t), target :: globalsize(1) = (/N/)
  real, dimension(N), target :: hostA, hostB, hostC, AfromDevice
  integer(c_size_t) :: arraysize = c_sizeof(hostA)
  integer(c_size_t), parameter :: zero = 0
  integer :: i, ierr, myrank, num_of_ranks

  ! Initialize MPI
  call MPI_INIT(ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD, num_of_ranks, ierr)

  ! Initialize B and C with random numbers
  call RANDOM_NUMBER(hostB)
  call RANDOM_NUMBER(hostC)

  ! Perform host vector add
  hostA = hostB + hostC

  ! Initialise the OpenCL device and environment (serialized for nicer output)
  do i=0, num_of_ranks
    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
    if (myrank == i) then
      write(*,*) "Initializing OpenCL context for MPI rank ", myrank
      call ocl_env_init(1, myrank + 1)
      cmd_queues => get_cmd_queues()
    endif
  enddo

  ! Allocate buffers on the device
  deviceA = create_wonly_buffer(arraysize)
  deviceB = create_ronly_buffer(arraysize)
  ! deviceC is only read, but we use rw_buffer to test FortCL API
  deviceC = create_rw_buffer(arraysize)

  ! Load the kernel and set up the arguments
  call add_kernels(1, "vadd", FILENAME)
  vaddkernel = get_kernel_by_name("vadd")
  ierr = clSetKernelArg(vaddkernel, 0, C_SIZEOF(deviceB), C_LOC(deviceB))
  call check_status('clSetKernelArg: arg 1 of vadd', ierr)
  ierr = clSetKernelArg(vaddkernel, 1, C_SIZEOF(deviceC), C_LOC(deviceC))
  call check_status('clSetKernelArg: arg 2 of vadd', ierr)
  ierr = clSetKernelArg(vaddkernel, 2, C_SIZEOF(deviceA), C_LOC(deviceA))
  call check_status('clSetKernelArg: arg 3 of vadd', ierr)

  ! Copy B and C to the device
  ierr = clEnqueueWriteBuffer(cmd_queues(1), deviceB, CL_TRUE, zero, &
        arraysize, C_LOC(hostB), 0, C_NULL_PTR, C_NULL_PTR)
  call check_status('enqueue write buffer B', ierr)
  ierr = clEnqueueWriteBuffer(cmd_queues(1), deviceC, CL_TRUE, zero, &
       arraysize, C_LOC(hostC), 0, C_NULL_PTR, C_NULL_PTR)
  call check_status('enqueue write buffer C', ierr)
  ierr = clFinish(cmd_queues(1))
  call check_status('clFinish', ierr)

  ! Execute kernel
   ierr = clEnqueueNDRangeKernel(cmd_queues(1), vaddkernel, 1, &
      C_NULL_PTR, C_LOC(globalsize), C_NULL_PTR, 0, C_NULL_PTR, C_NULL_PTR)
  call check_status('NDRangeKernel', ierr)
  ierr = clFinish(cmd_queues(1))
  call check_status('clFinish', ierr)

  ! Copy A from the device
  ierr = clEnqueueReadBuffer(cmd_queues(1), deviceA, CL_TRUE, zero, &
       arraysize, C_LOC(AfromDevice), 0, C_NULL_PTR, C_NULL_PTR)
  call check_status('Read buffer', ierr)
  ierr = clFinish(cmd_queues(1))
  call check_status('clFinish', ierr)

  ! Clean up OpenCL device
  call ocl_release()

  ! Finalise MPI
  call MPI_FINALIZE(ierr)

  write(*,*) "The MPI+OpenCL example has finalised successfully!"
end program vadd

max_stars_repo_name: stfc/FortCL

max_issues_repo_name: stfc/FortCL

max_forks_repo_name: stfc/FortCL


text:
      SUBROUTINE  DSWAP (N,DX,INCX,DY,INCY)
C
C     INTERCHANGES TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DY(1),DTEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        DTEMP = DX(IX)
        DX(IX) = DY(IY)
        DY(IY) = DTEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C       CLEAN-UP LOOP
C
   20 M = MOD(N,3)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
   30 CONTINUE
      IF( N .LT. 3 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,3
        DTEMP = DX(I)
        DX(I) = DY(I)
        DY(I) = DTEMP
        DTEMP = DX(I + 1)
        DX(I + 1) = DY(I + 1)
        DY(I + 1) = DTEMP
        DTEMP = DX(I + 2)
        DX(I + 2) = DY(I + 2)
        DY(I + 2) = DTEMP
   50 CONTINUE
      RETURN
      END


text:
      SUBROUTINE iau_TPSTV ( XI, ETA, V0, V )
*
*  - - - - - - - - - -
*   i a u _ T P S T V
*  - - - - - - - - - -
*
*  In the tangent plane projection, given the star's rectangular
*  coordinates and the direction cosines of the tangent point, solve
*  for the direction cosines of the star.
*
*  This routine is part of the International Astronomical Union's
*  SOFA (Standards of Fundamental Astronomy) software collection.
*
*  Status:  support routine.
*
*  Given:
*     XI,ETA    d       rectangular coordinates of star image (Note 2)
*     V0        d(3)    tangent point's direction cosines (Note 4)
*
*  Returned:
*     V         d(3)    star's direction cosines
*
*  1) The tangent plane projection is also called the "gnomonic
*     projection" and the "central projection".
*
*  2) The eta axis points due north in the adopted coordinate system.
*     If the direction cosines represent observed (RA,Dec), the tangent
*     plane coordinates (xi,eta) are conventionally called the "standard
*     coordinates".  If the direction cosines are with respect to a
*     right-handed triad, (xi,eta) are also right-handed.  The units of
*     (xi,eta) are, effectively, radians at the tangent point.
*
*  3) The method used is to complete the star vector in the (xi,eta)
*     based triad and normalize it, then rotate the triad to put the
*     tangent point at the pole with the x-axis aligned to zero
*     longitude.  Writing (a0,b0) for the celestial spherical
*     coordinates of the tangent point, the sequence of rotations is
*     (b0-pi/2) around the x-axis followed by (-a0-pi/2) around the
*     z-axis.
*
*  4) If vector V0 is not of unit length, the returned vector V will
*     be wrong.
*
*  5) If vector V0 points at a pole, the returned vector V will be
*     based on the arbitrary assumption that the longitude coordinate
*     of the tangent point is zero.
*
*  6) This routine is a member of the following set:
*
*         spherical       vector       solve for
*
*         iau_TPXES      iau_TPXEV      xi,eta
*         iau_TPSTS    > iau_TPSTV <     star
*         iau_TPORS      iau_TPORV      origin
*
*  References:
*
*     Calabretta M.R. & Greisen, E.W., 2002, "Representations of
*     celestial coordinates in FITS", Astron.Astrophys. 395, 1077
*
*     Green, R.M., "Spherical Astronomy", Cambridge University Press,
*     1987, Chapter 13.
*
*  This revision:   2018 January 2
*
*  SOFA release 2020-07-21
*
*  Copyright (C) 2020 IAU SOFA Board.  See notes at end.
*
*-----------------------------------------------------------------------

      IMPLICIT NONE

      DOUBLE PRECISION XI, ETA, V0(3), V(3)

      DOUBLE PRECISION X, Y, Z, R, F

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  Tangent point.
      X = V0(1)
      Y = V0(2)
      Z = V0(3)

*  Deal with polar case.
      R = SQRT(X*X+Y*Y)
      IF ( R .EQ. 0D0 ) THEN
         R = 1D-20
         X = R
      END IF

*  Star vector length to tangent plane.
      F = SQRT(1D0+XI*XI+ETA*ETA)

*  Apply the transformation and normalize.
      V(1) = ( X - (XI*Y+ETA*X*Z) / R ) / F
      V(2) = ( Y + (XI*X-ETA*Y*Z) / R ) / F
      V(3) = ( Z + ETA*R ) / F

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2020
*  Standards Of Fundamental Astronomy Board
*  of the International Astronomical Union.
*
*  =====================
*  SOFA Software License
*  =====================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
*  CONDITIONS WHICH APPLY TO ITS USE.
*
*  1. The Software is owned by the IAU SOFA Board ("SOFA").
*
*  2. Permission is granted to anyone to use the SOFA software for any
*     purpose, including commercial applications, free of charge and
*     without payment of royalties, subject to the conditions and
*     restrictions listed below.
*
*  3. You (the user) may copy and distribute SOFA source code to others,
*     and use and adapt its code and algorithms in your own software,
*     on a world-wide, royalty-free basis.  That portion of your
*     distribution that does not consist of intact and unchanged copies
*     of SOFA source code files is a "derived work" that must comply
*     with the following requirements:
*
*     a) Your work shall be marked or carry a statement that it
*        (i) uses routines and computations derived by you from
*        software provided by SOFA under license to you; and
*        (ii) does not itself constitute software provided by and/or
*        endorsed by SOFA.
*
*     b) The source code of your derived work must contain descriptions
*        of how the derived work is based upon, contains and/or differs
*        from the original SOFA software.
*
*     c) The names of all routines in your derived work shall not
*        include the prefix "iau" or "sofa" or trivial modifications
*        thereof such as changes of case.
*
*     d) The origin of the SOFA components of your derived work must
*        not be misrepresented;  you must not claim that you wrote the
*        original software, nor file a patent application for SOFA
*        software or algorithms embedded in the SOFA software.
*
*     e) These requirements must be reproduced intact in any source
*        distribution and shall apply to anyone to whom you have
*        granted a further right to modify the source code of your
*        derived work.
*
*     Note that, as originally distributed, the SOFA software is
*     intended to be a definitive implementation of the IAU standards,
*     and consequently third-party modifications are discouraged.  All
*     variations, no matter how minor, must be explicitly marked as
*     such, as explained above.
*
*  4. You shall not cause the SOFA software to be brought into
*     disrepute, either by misuse, or use for inappropriate tasks, or
*     by inappropriate modification.
*
*  5. The SOFA software is provided "as is" and SOFA makes no warranty
*     as to its use or performance.   SOFA does not and cannot warrant
*     the performance or results which the user may obtain by using the
*     SOFA software.  SOFA makes no warranties, express or implied, as
*     to non-infringement of third party rights, merchantability, or
*     fitness for any particular purpose.  In no event will SOFA be
*     liable to the user for any consequential, incidental, or special
*     damages, including any lost profits or lost savings, even if a
*     SOFA representative has been advised of such damages, or for any
*     claim by any third party.
*
*  6. The provision of any version of the SOFA software under the terms
*     and conditions specified herein does not imply that future
*     versions will also be made available under the same terms and
*     conditions.
*
*  In any published work or commercial product which uses the SOFA
*  software directly, acknowledgement (see www.iausofa.org) is
*  appreciated.
*
*  Correspondence concerning SOFA software should be addressed as
*  follows:
*
*      By email:  sofa@ukho.gov.uk
*      By post:   IAU SOFA Center
*                 HM Nautical Almanac Office
*                 UK Hydrographic Office
*                 Admiralty Way, Taunton
*                 Somerset, TA1 2DN
*                 United Kingdom
*
*-----------------------------------------------------------------------

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
c*****************************************************************************
c
c      Manifold pointer documentation and setup routines
c
c******************************************************************************
c
       SUBROUTINE  SPTRXX
c
c      Sets pointers to variables in the manifold work arrays
c============================================================================
c      Explanation of pointers
c
c      In order to enable the reference to a large block of data with
c      a single subroutine argument, work arrays are used. Each manifold
c      has a real work array RW and an integer work array IW.
c
c      Quantities in the work arrays are indexed by pointers in the
c      common block /MPTRXX/. Pointers have names IoWXYX, where Io stands
c      for "index of" and WXYX is the four-character (max) variable name.
c      Thus, WXYZ is referred to as RW(IoWXYZ). This makes the code easy
c      to read.
c
c      Scalar pointers are direct. The pointer for an array points to the
c      position just before the first array element. Then the Ith element
c      of a vector WXYZ is referred to as RW(IoWXYZ+I). The (I,J) element
c      of an KxN matrix WXYZ is referred to as RW(IoWXYZ+I+K*J).  Hence,
c      for a scalar the pointer is to its location in the array.  For a
c      vector the pointer is to one place before the first element of the
c      vector. For an KxN matrix the pointer is to a place K+1 before
c      the first element.
c
c      The intake and exhaust manifolds have the same number of pointers,
c      so those in MPTRXX apply either.
c
c      The intake manfold is MANifold 1
c          Point 1: start of feeder downstream of ambient plenum
c          Point 2: discharge of feeder into junction
c          Point 3: entrance to runner from junction
c          Point 4: discharge from runner to cylinder (before EGR)
c
c      The intake exhaust manifold is MANifold 2
c          Point 1: intake from cylinder (after EGR)
c          Point 2: discharge of runner into junction
c          Point 3: entrance to collector from junctiom
c          Point 4: discharge from collector to ambient
c============================================================================
c      The layout of the data in the work arrays is described below.
c
c  ** Real work array RW  (RWI or RWE)
c    location  symbol          quantity
c------------------------------------------------------------------------------
c  k-dependent constants
c      1       K               k, specific heat ratio (1.4 for air)
c      2       KM1             k - 1
c      3       KP1             k + 1
c      4       KMKP            (k - 1)/(k + 1)
c      5       KPKM            (k + 1)/(k - 1)
c      6       KM12            (k - 1)/2
c      7       KM3M            (k - 1)/(3 - k)
c      8       KKM1            k/(k - 1)
c      9       KM1K            (k - 1)/k
c     10       KPO2            (k+1)/2k
c     11       RK              1/k
c     12       R2K             1/2k
c     13       RKM1            1/(k - 1)
c     14       RKKM            k^(-1/(k-1))
c     15       TK              2k
c     16       TRKM            2/(k - 1)
c     17       TRKP            2/(k + 1)
c     18       TOKS            [2/(k - 1)]^2
c     19       TKOM            2*k/(k-1)
c     20       FRKM            4/(k - 1)
c     21       SKMS            16/(k - 1)^2
c     22       HKPM            (k+1)/(2(k-1))
c     23       HKMK            (k-1)/(2k)
c     24       C0CS            sqrt{(k + 1)/2}
c     25       XMAX            sqrt{2/(k + 1)}[(k + 1)/(k - 1)]
c     26       XMIN            sqrt{2/(k + 1)}[(3 - k)/(k - 1)]
c     27       PRPC            1/[(k + 1)/2]^{(k/(k-1)} = P_*/P_0 (s)
c     28       CRCC            1/sqrt((k+1)/2) =  C*/C0
c     29       RCHM            1/( (k+1)/2)**((k+1)/2(k-1)) )
c     30       SQTK            sqrt(k)
c     31       PHMX            1/((k+1)/(k-1))**((k+1)/(k-1))
c     32       SCKF            1/(k(k-1)) (for entropy friction correction)
c     33       TKOP            2k/(k+1)
c     34       K2M1            k^2-1
c     35       HKP1            (k+1)/2
c     36       SFMC            sqrt{2/(3-k)}
c     37       CMD(1)          coefficient in expansion for valve outflow MD
c     38       CMD(2)           "
c     39       CMD(3)           "
c     40       CMD(4)           "
c     41       CMD(5)           "
c     42       CMD(6)           "
c------------------------------------------------------------------------------
c     43-49    reserved
c------------------------------------------------------------------------------
c  constants involving the manifold blockage at N=1,4  (4 not used)
c     50       reserved
c     51-54    MBF(N)           blockage fraction
c     55-58    AMR(N)           restriction flow area
c------------------------------------------------------------------------------
c   averaging weights for the acoustic time delays
c     59       FTDO             delay-time averaging weight for old D cycles
c     60       FTDN             delay-time averaging weight for new D cycles
c     61       FTEO             delay-time averaging weight for old E cycles
c     62       FTEN             delay-time averaging weight for new E cycles
c------------------------------------------------------------------------------
c   63-70      reserved
c------------------------------------------------------------------------------
c  array data at angle I at point N, i=1,720   N=1,4
c    71-2950    C(I,N)         sound speed
c  2951-5830    V(I,N)         velocity in the nominal direction
c  5831-8710    Y(I,N)         entropy parameter P/rho^k
c  8711-11590   P(I,N)         pressure
c 11591-14470   W(I,N)         mass/area that has passed point J since t=0
c 14471-17350   Z(I,N)         phase-averaged arriving acoustic parameter
c------------------------------------------------------------------------------
c  array data at angle I=1,720
c 17351-18070   P0J(I)         junction stagnation pressure
c 18071-18790   YJ(I)          junction stagnation entropy parameter P/rho^k
c------------------------------------------------------------------------------
c  array data at angle I=1,720
c 18791-19510    YPAV(I)       phase-averaged Y, used in junction analysis
c 19511-20230    ZPAC(I)       phase-averaged Z, used in junction analysis
c 20231-20950    reserved
c------------------------------------------------------------------------------
c  data for Z correction at point N
c 20951-20954    DZF(N)        friction correction to Z
c------------------------------------------------------------------------------
c  data for acoustic characteristic arriving at point N=1,4
c 20955-20958    CL(N)         speed of sound at launch
c 20959-20962    VL(N)         velocity at launch
c 20963-20966    YL(N)         entropy measure at launch
c------------------------------------------------------------------------------
c  data for point N=1,4 at the current time
c 20967-20970    G(N)          mass velocity rho*V
c 20971-20974    WCYC(N)       mass/area past point N for last 720 degrees
c 20975-20978    F(N)          total mass flow past N for current cycle (~/DT)
c 20979-20982    RSUM(N)       accumulating cycle density sum
c 20983-20986    DZC(N)        Z correction for mass unbalance
c------------------------------------------------------------------------------
c  data for for duct J=1,2   variables at the current time
c 20987-20988    LD(J)         duct length
c 20989-20990    AD(J)         flow area
c 20991-20992    TDP(J)        time delay for Z+ characteristics
c 20993-20994    TDM(J)        time delay for Z- characteristics
c 20995-20996    TDPA(J)       cycle average TDP at start of cycle
c 20997-20998    TDMA(J)       cycle average TDM at start of cycle
c 20999-21000    DWC(J)        W correction for unbalanced flows
c------------------------------------------------------------------------------
c  constants involving the friction factor and k in duct J=1,2
c 21001-21002    CF(J)         friction factor
c 21003-21004    CFCB(J)       convection backflow coefficient
c 21005-21006    CFCT(J)       convection throughtflow coefficient
c 21007-21008    CFA(J)        acoustic chracteristic friction coefficient
c------------------------------------------------------------------------------
c  other scalars
c    21009       FAPO          weight to old data in phase averages
c    21010       FAPN          weight to new data in phase averages
c    21011       DZCO          Z  correction for non-cumputed runners
c    21012       FOUT          total junct. flow to other cyls. current cycle
c    21013       EMB           manifold mass unbalance, current cycle
c    21014       R             gas constant (SI units)
c------------------------------------------------------------------------------
c 21015-21024   reserved
c==============================================================================
c  ** Integer work array IW  (IWI or IWE)
c    location  symbol          quantity
c------------------------------------------------------------------------------
c       1-2880   JUMP(I,N)     jump at angle I=1,720 for point N=1,4
c    2881-2884   LAG(N)        lag at point N=1,4
c    2885-2888                 reserved
c         2889   JTAU          degrees per cylinder cycle (720/Ncyl)
c         2890   JCAD          degrees into current feeder/collector D cycle
c         2891   NCYE          engine cycle number
c         2892   NCYD          feeder/collector duct cycle number
c         2893   NCYL          number of cylinders on the manifold
c         2894   IBEG          index at start of cycle calculation
c         2895   MAN           manfold type; 1=intake, 2=exhaust
c         2896   I             current angle index (1-720) for Y and Z
c    2897-2900   IV(N)         index of min Win point
c    2901-2904   IP(N)         index of max Win point
c============================================================================
c   code for memory reservations and pointer locations
c
       PARAMETER           (NMPTR = 96)
c------------------------------------------------------------------------------
       DIMENSION           MPTR(NMPTR)
       COMMON  /MPTRXX/    MPTR
c------------------------------------------------------------------------------
c    Pointers for the work arrays
       EQUIVALENCE (IoCMD2,MPTR(37))
C      EQUIVALENCE (,MPTR(38))
       EQUIVALENCE (IoMBF ,MPTR(39))
       EQUIVALENCE (IoAMR ,MPTR(40))
       EQUIVALENCE (IoFTDO,MPTR(41))
       EQUIVALENCE (IoFTDN,MPTR(42))
C      EQUIVALENCE (IoZPAV,MPTR(43))
       EQUIVALENCE (IoC   ,MPTR(44))
       EQUIVALENCE (IoV   ,MPTR(45))
       EQUIVALENCE (IoY   ,MPTR(46))
       EQUIVALENCE (IoP   ,MPTR(47))
       EQUIVALENCE (IoW   ,MPTR(48))
       EQUIVALENCE (IoZ   ,MPTR(49))
       EQUIVALENCE (IoP0J ,MPTR(50))
       EQUIVALENCE (IoYJ  ,MPTR(51))
       EQUIVALENCE (IoFI  ,MPTR(52))
       EQUIVALENCE (IoQI  ,MPTR(53))
       EQUIVALENCE (IoFO  ,MPTR(54))
       EQUIVALENCE (IoDZF ,MPTR(55))
       EQUIVALENCE (IoCL  ,MPTR(56))
       EQUIVALENCE (IoVL  ,MPTR(57))
       EQUIVALENCE (IoYL  ,MPTR(58))
       EQUIVALENCE (IoG   ,MPTR(59))
       EQUIVALENCE (IoWCYC,MPTR(60))
       EQUIVALENCE (IoF   ,MPTR(61))
       EQUIVALENCE (IoRSUM,MPTR(62))
       EQUIVALENCE (IoDZC ,MPTR(63))
       EQUIVALENCE (IoLD  ,MPTR(64))
       EQUIVALENCE (IoAD  ,MPTR(65))
       EQUIVALENCE (IoTDP ,MPTR(66))
       EQUIVALENCE (IoTDM ,MPTR(67))
       EQUIVALENCE (IoTDPA,MPTR(68))
       EQUIVALENCE (IoTDMA,MPTR(69))
       EQUIVALENCE (IoDWC ,MPTR(70))
       EQUIVALENCE (IoCF  ,MPTR(71))
       EQUIVALENCE (IoCFCB,MPTR(72))
       EQUIVALENCE (IoCFCT,MPTR(73))
       EQUIVALENCE (IoCFA ,MPTR(74))
       EQUIVALENCE (IoFTEO,MPTR(75))
       EQUIVALENCE (IoFTEN,MPTR(76))
       EQUIVALENCE (IoYPAV,MPTR(77))
       EQUIVALENCE (IoZPAV,MPTR(78))
       EQUIVALENCE (IoEMB ,MPTR(79))
       EQUIVALENCE (IoR   ,MPTR(80))
       EQUIVALENCE (IoJUMP,MPTR(81))
       EQUIVALENCE (IoLAG ,MPTR(82))
       EQUIVALENCE (IoFOUT,MPTR(83))
       EQUIVALENCE (IoJTAU,MPTR(84))
       EQUIVALENCE (IoJCAD,MPTR(85))
       EQUIVALENCE (IoNCYE,MPTR(86))
       EQUIVALENCE (IoNCYD,MPTR(87))
       EQUIVALENCE (IoNCYL,MPTR(88))
       EQUIVALENCE (IoIBEG,MPTR(89))
       EQUIVALENCE (IoMAN ,MPTR(90))
       EQUIVALENCE (IoI   ,MPTR(91))
       EQUIVALENCE (IoIV  ,MPTR(92))
       EQUIVALENCE (IoIP  ,MPTR(93))
       EQUIVALENCE (IoDZCO,MPTR(94))
       EQUIVALENCE (IoFAPO,MPTR(95))
       EQUIVALENCE (IoFAPN,MPTR(96))
c-----------------------------------------------------------------------------
c ** real work array
c
c    set all pointers null
       DO 3 I=1,NMPTR
           MPTR(I) = 0
3          CONTINUE
c    k-dependent scalar pointers
       DO 9 I=1,36
           MPTR(I) = I
9          CONTINUE
c    MD coefficients
       IoCMD2 = 36
c    vector pointers
       IoMBF  = 50
       IoAMR  = 54
c    more scalar pointers
       IoFTDO = 59
       IoFTDN = 60
       IoFTEO = 61
       IoFTEN = 62
c    matrix pointers 720x4
       IoC    = 71 - 721
       IoV    = 2951 - 721
       IoY    = 5831 - 721
       IoP    = 8711 - 721
       IoW    = 11591 - 721
       IoZ    = 14471 - 721
c    vector pointers
       IoP0J  = 17350
       IoYJ   = 18070
       IoYPAV = 18790
       IoZPAV = 19510
       IoDZF  = 20950
       IoCL   = 20954
       IoVL   = 20958
       IoYL   = 20962
       IoG    = 20966
       IoWCYC = 20970
       IoF    = 20974
       IoRSUM = 20978
       IoDZC  = 20982
       IoLD   = 20986
       IoAD   = 20988
       IoTDP  = 20990
       IoTDM  = 20992
       IoTDPA = 20994
       IoTDMA = 20996
       IoDWC  = 20998
c    k-dependent vector pointers 2
       IoCF   = 21000
       IoCFCB = 21002
       IoCFCT = 21004
       IoCFA  = 21006
c    scalars
       IoFAPO = 21009
       IoFAPN = 21010
       IoDZCO = 21011
       IoFOUT = 21012
       IoEMB  = 21013
       IoR    = 21014
c
c ** integer work array
       IoJUMP = 1 - 721
       IoLAG  = 2880
       IoIP   = 2884
       IoIV   = 2886
       IoJTAU = 2889
       IoJCAD = 2890
       IoNCYE = 2891
       IoNCYD = 2892
       IoNCYL = 2893
       IoIBEG = 2894
       IoMAN  = 2895
       IoI    = 2896
       IoIV   = 2896
       IoIP   = 2900
c
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE  SETKXX(RW)
c
c      Loads work array RW with k-dependent constants used repeatedly.
c      K must be loaded in RW(IoK) at call.
c------------------------------------------------------------------------------
c      K               k, specific heat ratio
c      KM1             k - 1
c      KP1             k + 1
c      KMKP            (k - 1)/(k + 1)
c      KPKM            (k + 1)/(k - 1)
c      KM12            (k - 1)/2
c      KM3M            (k - 1)/(3 - k)
c      KKM1            k/(k - 1)
c      KM1K            (k - 1)/k
c      KPO2            (k+1)/2k
c      RK              1/k
c      R2K             1/2k
c      RKM1            1/(k - 1)
c      RKKM            k^(-1/(k-1))
c      TK              2k
c      TRKM            2/(k - 1)
c      TRKP            2/(k + 1)
c      TOKS            [2/(k - 1)]^2
c      TKOM            2*k/(k-1)
c      FRKM            4/(k - 1)
c      SKMS            16/(k - 1)^2
c      HKPM            (k+1)/(2(k-1))
c      HKMK            (k-1)/(2k)
c      C0CS            sqrt{(k + 1)/2}
c      XMAX            sqrt{2/(k + 1)}[(k + 1)/(k - 1)]
c      XMIN            sqrt{2/(k + 1)}[(3 - k)/(k - 1)]
c      PRPC            1/[(k + 1)/2]^{(k/(k-1)} = P_*/P_0 (s)
c      CRCC            1/sqrt((k+1)/2) =  C*/C0
c      RCHM            1/( (k+1)/2)**((k+1)/2(k-1)) )
c      SQTK            sqrt(k)
c      PHMX            1/((k+1)/(k-1))**((k+1)/(k-1))
c      SCKF            1/(k(k-1))  (entropy correction K factor)
c      TKOP            2k/(k+1)
c      K2M1            k^2-1
c      HKP1            (k+1)/2
c      SFMC            sqrt{2/(3-k)}
c      CMD2(N)         coefficients in expansion for MD2
c------------------------------------------------------------------------------
       PARAMETER           (NMPTR = 96)
       PARAMETER           (NDRW = 21024)
       DIMENSION           MPTR(NMPTR)
       DIMENSION           RW(NDRW)
       COMMON  /MPTRXX/    MPTR
c------------------------------------------------------------------------------
       EQUIVALENCE (IoK   ,MPTR(1))
       EQUIVALENCE (IoKM1 ,MPTR(2))
       EQUIVALENCE (IoKP1 ,MPTR(3))
       EQUIVALENCE (IoKMKP,MPTR(4))
       EQUIVALENCE (IoKPKM,MPTR(5))
       EQUIVALENCE (IoKM12,MPTR(6))
       EQUIVALENCE (IoKM3M,MPTR(7))
       EQUIVALENCE (IoKKM1,MPTR(8))
       EQUIVALENCE (IoKM1K,MPTR(9))
       EQUIVALENCE (IoKPO2,MPTR(10))
       EQUIVALENCE (IoRK  ,MPTR(11))
       EQUIVALENCE (IoR2K ,MPTR(12))
       EQUIVALENCE (IoRKM1,MPTR(13))
       EQUIVALENCE (IoRKKM,MPTR(14))
       EQUIVALENCE (IoTK  ,MPTR(15))
       EQUIVALENCE (IoTRKM,MPTR(16))
       EQUIVALENCE (IoTRKP,MPTR(17))
       EQUIVALENCE (IoTOKS,MPTR(18))
       EQUIVALENCE (IoTKOM,MPTR(19))
       EQUIVALENCE (IoFRKM,MPTR(20))
       EQUIVALENCE (IoSKMS,MPTR(21))
       EQUIVALENCE (IoHKPM,MPTR(22))
       EQUIVALENCE (IoHKMK,MPTR(23))
       EQUIVALENCE (IoC0CS,MPTR(24))
       EQUIVALENCE (IoXMAX,MPTR(25))
       EQUIVALENCE (IoXMIN,MPTR(26))
       EQUIVALENCE (IoPRPC,MPTR(27))
       EQUIVALENCE (IoCRCC,MPTR(28))
       EQUIVALENCE (IoRCHM,MPTR(29))
       EQUIVALENCE (IoSQTK,MPTR(30))
       EQUIVALENCE (IoPHMX,MPTR(31))
       EQUIVALENCE (IoSCKF,MPTR(32))
       EQUIVALENCE (IoTKOP,MPTR(33))
       EQUIVALENCE (IoK2M1,MPTR(34))
       EQUIVALENCE (IoHKP1,MPTR(35))
       EQUIVALENCE (IoSFMC,MPTR(36))
       EQUIVALENCE (IoCMD2,MPTR(37))
c------------------------------------------------------------------------------
       REAL    K,K2
c------------------------------------------------------------------------------
c    diagnostic monitor
       LOGICAL OUT
       COMMON /MONMXX/IUMONM,MON,OUT
c------------------------------------------------------------------------------
c    compute the constants
       K = RW(IoK)
       RW(IoTK) = 2*K
       RW(IoKM1) = K - 1
       RW(IoKP1) = K + 1
       RW(IoKMKP) = RW(IoKM1)/RW(IoKP1)
       RW(IoKPKM) = RW(IoKP1)/RW(IoKM1)
       RW(IoKM12) = RW(IoKM1)/2
       RW(IoKM3M) = RW(IoKM1)/(3 - K)
       RW(IoKKM1) = K/RW(IoKM1)
       RW(IoKM1K) = 1/RW(IoKKM1)
       RW(IoKPO2) = RW(IoKP1)/RW(IoTK)
       RW(IoRK) = 1/K
       RW(IoR2K) = 1/RW(IoTK)
       RW(IoRKM1) = 1/RW(IoKM1)
       RW(IoRKKM) = 1/(K**RW(IoRKM1))
       RW(IoTRKM) = 2*RW(IoRKM1)
       RW(IoTRKP) = 2/RW(IoKP1)
       RW(IoTOKS) = RW(IoTRKM)*RW(IoTRKM)
       RW(IoTKOM) = 2*RW(IoKKM1)
       RW(IoSKMS) = 16*RW(IoRKM1)*RW(IoRKM1)
       RW(IoFRKM) = 4*RW(IoRKM1)
       RW(IoHKPM) = RW(IoKPKM)/2
       RW(IoHKMK) = RW(IoKM1K)/2
       RW(IoC0CS) = SQRT(RW(IoKP1)/2)
       RW(IoCRCC) = 1/RW(IoC0CS)
       RW(IoXMAX) = RW(IoCRCC)*RW(IoKPKM)
       RW(IoXMIN) = RW(IoCRCC)/RW(IoKM3M)
       RW(IoPRPC) = RW(IoTRKP)**RW(IoKKM1)
       RW(IoRCHM) = RW(IoTRKP)**RW(IoHKPM)
       RW(IoSQTK) = SQRT(K)
       RW(IoPHMX) = 1/RW(IoKPKM)**RW(IoKPKM)
       RW(IoSCKF) = 1/(K*(K-1))
       RW(IoTKOP) = 1/RW(IoKPO2)
       RW(IoK2M1) = RW(IoK)**2 - 1
       RW(IoHKP1) = RW(IoKP1)/2
       RW(IoSFMC) = SQRT(2/(3 - RW(IoK)))
c
c    coefficients in expansion for MD2 about the quadratic singular point
       COEFR = RW(IoRKM1)*RW(IoRK)
       K2 = K*K
       OPKS = 1 + K2
       OPKS2 = OPKS*OPKS
       TERM = K2 - OPKS2/4
c    series developed by Mathematica and checked numerically for k=1.4
       RW(IoCMD2+1) = RW(IoHKMK)
       RW(IoCMD2+2) = COEFR*TERM/2
       RW(IoCMD2+3) = COEFR*(- OPKS*TERM)/4.
       RW(IoCMD2+4) = COEFR*(-0.5*K2*TERM + 0.625*OPKS2*TERM)/4.
       RW(IoCMD2+5) = COEFR*(0.5*K2*OPKS*TERM
     ;         -0.875*OPKS*(-0.5*K2*TERM
     ;         + 0.625*OPKS2*TERM))/5.
       RW(IoCMD2+6) = COEFR*(-0.75*K2*(-0.5*K2*TERM
     ;                                +0.625*OPKS2*TERM)
     ;             -0.9*OPKS*(0.5*K2*OPKS*TERM
     ;             -0.875*OPKS*(-0.5*K2*TERM
     ;             +0.625*OPKS2*TERM)))/6.
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE SCFLXX(RW)
c
c      Loads work array RW with flow constants used repeatedly.
c      CF(J) is the friction factor list and must be loaded in
c      RW(IoCF+JN) J=1,24 at call (after SETKXX).
c------------------------------------------------------------------------------
       PARAMETER           (NMPTR = 96)
       PARAMETER           (NDRW = 21024)
       DIMENSION           MPTR(NMPTR)
       DIMENSION           RW(NDRW)
       COMMON  /MPTRXX/    MPTR
c------------------------------------------------------------------------------
       EQUIVALENCE (IoSCKF,MPTR(32))
       EQUIVALENCE (IoMBF ,MPTR(39))
       EQUIVALENCE (IoAMR ,MPTR(40))
       EQUIVALENCE (IoLD  ,MPTR(64))
       EQUIVALENCE (IoAD  ,MPTR(65))
       EQUIVALENCE (IoCF  ,MPTR(71))
       EQUIVALENCE (IoCFCB,MPTR(72))
       EQUIVALENCE (IoCFCT,MPTR(73))
       EQUIVALENCE (IoCFA ,MPTR(74))
c------------------------------------------------------------------------------
       DATA PI/ 3.14159/
c------------------------------------------------------------------------------
c    friction
       DO 9 J=1,2
           BOA = PI*SQRT(4*RW(IoAD+J)/PI)/RW(IoAD+J)
           TERM = 0.5*RW(IoCF+J)*BOA*RW(IoLD+J)/2
           RW(IoCFCB+J) = RW(IoSCKF)*0.5*RW(IoCF+J)*BOA/32
           RW(IoCFCT+J) = RW(IoSCKF)*TERM
           RW(IoCFA+J) = TERM
9          CONTINUE
       RETURN
       END
c******************************************************************************

max_stars_repo_name: kjfergu/esp

max_issues_repo_name: fergu/esp

max_forks_repo_name: fergu/esp


text:
C    Copyright(C) 2014-2017 National Technology & Engineering Solutions of
C    Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: wrgens.f,v 1.9 1998/03/23 05:17:54 gdsjaar Exp $
C $Log: wrgens.f,v $
C Revision 1.9  1998/03/23 05:17:54  gdsjaar
C Fixed data statement ordering
C
C Revision 1.8  1997/09/02 18:51:45  caforsy
C Increase inital value of NIN in rdfsq.f and IGUESS in wrgens.f
C
C Revision 1.7  1992/11/17 15:06:39  gdsjaar
C Fixed problem with version number output.
C
c Revision 1.6  1992/11/16  23:32:47  gdsjaar
c Fixed problem with outputting too many attributes for some elements,
c upped version number to 2.2X
c
c Revision 1.5  1992/06/09  22:31:18  gdsjaar
c Fixed problem with sideset node and element numbering for 8 and 9 node elements
c
c Revision 1.4  1992/04/15  22:56:38  gdsjaar
c Fixed EXODUS output of 3-node beams--connectivity
c conforms to EXODUS document
c Upped version number to 2.1X to reflect change
c
c Revision 1.3  1991/12/18  21:12:03  gdsjaar
c Removed the info field
c
c Revision 1.2  1990/11/30  11:30:08  gdsjaar
c Rewrote indexing for reads and writes
c
c Revision 1.1.1.1  90/11/30  11:17:52  gdsjaar
c FASTQ Version 2.0X
c
c Revision 1.1  90/11/30  11:17:51  gdsjaar
c Initial revision
c
C
CC* FILE: [.MAIN]WRGENS.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE WRGENS (MS, MR, NPNODE, NPELEM, MXNFLG, MXSFLG, NPREGN,
     &   NPNBC, NPSBC, IUNIT, NNN, KKK, NNXK, NODES, NELEMS, NNFLG,
     &   NNPTR, NNLEN, NSFLG, NSPTR, NSLEN, NVPTR, NVLEN, NSIDEN,
     &   MAPDXG, XN, YN, NXK, MAT, MAPGXD, MATMAP, WTNODE, WTSIDE,
     &   NBCNOD, NNLIST, NBCSID, NSLIST, NVLIST, NUMMAT, LINKM, TITLE,
     &   ERR, EIGHT, NINE, VERSN)
C************************************************************************
C
C  SUBROUTINE WRGENS = WRITES GENESIS DATABASE MESH OUTPUT
C
C***********************************************************************
C
      PARAMETER (IGUESS = 1000)
C
C  IGUESS IS THE NUMBER OF ELEMENT BLOCKS,  FOR USE WITH THE ENAME
C  VARIABLE.  IF THIS VARIABLE IS NOT AS LARGE AS NUMMAT,  IT WILL NOT
C  RESULT IN A FATAL ERROR,  BUT SIMPLY A WARNING,  AND NO ELEMENT
C  NAMES WILL BE WRITTEN.
C
      DIMENSION XN (NPNODE), YN (NPNODE), NXK (NNXK, NPELEM)
      DIMENSION MAT (NPELEM)
      DIMENSION NODES (NPNBC), NELEMS (NPSBC), NSIDEN (NPSBC)
      DIMENSION NNFLG (MXNFLG), NNLEN (MXNFLG)
      DIMENSION NNPTR (MXNFLG), WTNODE (NPNBC)
      DIMENSION NSFLG (MXSFLG), NSLEN (MXSFLG)
      DIMENSION NSPTR (MXSFLG), WTSIDE (NPSBC)
      DIMENSION NVLEN (MXSFLG), NVPTR (MXSFLG), LINKM (2,  (MS+MR))
      DIMENSION MAPDXG (NPNODE), MAPGXD (NPNODE), MATMAP (3, NPREGN)
C
      CHARACTER*72 TITLE, HOLD*80
      CHARACTER*8 DATE, TIME, VERSN1, VERSN2, XNAME, YNAME
      CHARACTER*8 ENAME (IGUESS)
      CHARACTER*10 VERSN
C
      LOGICAL ERR, EIGHT, NINE
C
      integer lcon(9)
      integer lbar(3)

      data lcon /1,3,5,7,2,4,6,8,9/
      data lbar /1,3,2/

      CALL EXDATE (DATE)
      CALL EXTIME (TIME)
      VERSN1 = '        '
      VERSN2 = '        '
      VERSN1 = VERSN (1:5)
      VERSN2 = VERSN (6:10)
      ERR = .TRUE.
      HOLD = TITLE
      XNAME = 'X'
      YNAME = 'Y'
C
C  CHECK TO MAKE SURE THAT THERE IS ENOUGH ROOM FOR ELEMENT NAMES
C
      IF (NUMMAT.GT.IGUESS) THEN
         CALL MESAGE ('WARNING:  THE NUMBER OF ELEMENT BLOCKS EXCEEDS')
         CALL MESAGE ('          THE CAPACITY TO NAME EACH BLOCK.')
         CALL MESAGE ('          NO ELEMENT NAMES WILL BE WRITTEN.')
      ENDIF
C
C  WRITE OUT HEADER INFORMATION
C
      WRITE (IUNIT, ERR = 110)HOLD
      IJK = 2
      IVERS = 1
C ... Fix up side set nodes for 8 and 9 node elements.
C ... At this point, they are treated as two linear segments,
C ... They should be ends followed by middle

C    1-----3-----2 Now: 1 3 3 2 Correct: 1 2 3

      if (eight .or. nine) then
         nvlst = nvlist / 4 * 3
         if (nslist .gt. 1) then
            nslst = nslist / 2
         else
            nslst = nslist
         end if
      else
         nvlst = nvlist
         nslst = nslist
      end if

      WRITE (IUNIT, ERR = 110)NNN, IJK, KKK, NUMMAT, NBCNOD, NNLIST,
     &   NBCSID, NSLST, NVLST, IVERS
C
C  WRITE OUT NODE BLOCK
C
      WRITE (IUNIT, ERR = 110) (XN (I), I = 1, NNN),
     &   (YN (I), I = 1, NNN)
      WRITE (IUNIT, ERR = 110) (MAPDXG (I), I = 1, KKK)
C
C  WRITE OUT ELEMENT BLOCKS
C
      DO 100 I = 1, NUMMAT
         IF (NXK (3, MATMAP (2, I)) .EQ. 0) THEN
            INODE = 2
            NATTR = 1
            ATTR = 1.
            IF (I.LE.IGUESS)ENAME (I) = 'BEAM'
         ELSEIF (NXK (4, MATMAP (2, I)) .EQ. 0) THEN
            INODE = 3
            NATTR = 1
            ATTR = 1.
            IF (I.LE.IGUESS)ENAME (I) = 'BEAM'
            CALL MESAGE ('NOTE:  The connectivity numbering for 3-node')
            CALL MESAGE ('       beams/trusses has been fixed to')
            CALL MESAGE ('       conform to EXODUS convention (1-3-2)')
         ELSEIF (EIGHT) THEN
            INODE = 8
            NATTR = 0
            IF (I.LE.IGUESS)ENAME (I) = 'QUAD'
         ELSEIF (NINE) THEN
            INODE = 9
            NATTR = 0
            IF (I.LE.IGUESS)ENAME (I) = 'QUAD'
         ELSE
            INODE = 4
            NATTR = 0
            IF (I.LE.IGUESS)ENAME (I) = 'QUAD'
         ENDIF
C
C  NLOOP IS NEEDED TO WRITE SOMETHING OUT THE CURRENT COUNTER IS ZERO.
C  THIS IS DONE TO SOLVE A CRAY OPERATING SYSTEM [CTSS] PROBLEM
C  WHERE NULL RECORD WRITES ARE NOT DONE APPROPRIATELY
C
         WRITE (IUNIT, ERR = 110) MATMAP (1, I),
     &      MATMAP (3, I) - MATMAP (2, I)+1, INODE, NATTR
C... 8 or 9 node quads
         IF (INODE .EQ. 8 .or. inode .eq. 9) THEN
            write (iunit, err = 110) ((nxk(lcon(ii), k), ii=1, nnxk),
     $           k = matmap(2,i), matmap(3,i))
C... 3 node beam/truss
         ELSEIF (INODE .EQ. 3) THEN
            write (iunit, err = 110) ((nxk(lbar(ii), k), ii=1, inode),
     $           k = matmap(2,i), matmap(3,i))
C... 4 node quad or 2 node beam/truss
         ELSE
            WRITE (IUNIT, ERR = 110) ((NXK (J, K), J = 1, INODE),
     &         K = MATMAP (2, I), MATMAP (3, I))
         ENDIF
         NLOOP = MAX0 (1, NATTR*KKK)
         WRITE (IUNIT, ERR = 110) (ATTR, J = 1, NLOOP)
  100 CONTINUE
C
C  WRITE OUT NODAL BOUNDARY FLAGS
C
      NLOOP = MAX0 (1, NBCNOD)
      WRITE (IUNIT, ERR = 110) (NNFLG (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NNLEN (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NNPTR (I), I = 1, NLOOP)
      NLOOP = MAX0 (1, NNLIST)
      WRITE (IUNIT, ERR = 110) (NODES (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (WTNODE (I), I = 1, NLOOP)
C
C  WRITE OUT SIDE BOUNDARY FLAGS
C

C ... Fix up side set nodes and elements for 8 and 9 node elements.
C ... At this point, they are treated as two linear segments,
C ... They should be ends followed by middle

C    1-----3-----2 Now: 1 3 3 2 Correct: 1 2 3

      if (eight .or. nine) then
         ipn = 1
         ipe = 1
         do 130 ibc = 1, nbcsid

            ibee = nsptr(ibc)
            iene = nsptr(ibc) + nslen(ibc) - 1
            nsptr(ibc) = ipe
            nslen(ibc) = max(1, nslen(ibc)/2)
            do 115 iel = ibee, iene, 2
               nelems(ipe)   = nelems(iel)
               ipe = ipe + 1
 115        continue

            ibeg = nvptr(ibc)
            iend = nvptr(ibc) + nvlen(ibc) - 1
            nvptr(ibc) = ipn
            nvlen(ibc) = nvlen(ibc) / 4 * 3

            do 120 inod = ibeg, iend, 4
               nsiden(ipn)   = nsiden(inod)
               nsiden(ipn+2) = nsiden(inod+1)
               nsiden(ipn+1) = nsiden(inod+3)
               wtside(ipn)   = wtside(inod)
               wtside(ipn+2) = wtside(inod+1)
               wtside(ipn+1) = wtside(inod+3)
               ipn = ipn + 3
 120        continue
 130     continue
         nvlist = ipn - 1
         nslist = ipe - 1
      end if

      NLOOP = MAX0 (1, NBCSID)
      WRITE (IUNIT, ERR = 110) (NSFLG (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NSLEN (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NVLEN (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NSPTR (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (NVPTR (I), I = 1, NLOOP)
      NLOOP = MAX0 (1, NSLIST)
      WRITE (IUNIT, ERR = 110) (NELEMS (I), I = 1, NLOOP)
      NLOOP = MAX0 (1, NVLIST)
      WRITE (IUNIT, ERR = 110) (NSIDEN (I), I = 1, NLOOP)
      WRITE (IUNIT, ERR = 110) (WTSIDE (I), I = 1, NLOOP)
C
C  WRITE OUT THE QA INFORMATION
C
      IHOLD = 1
      WRITE (IUNIT, ERR = 110)IHOLD
      WRITE (IUNIT, ERR = 110)VERSN1, VERSN2, DATE, TIME
C
C  WRITE THE HEADER INFORMATION
C
      IHOLD = 0
      WRITE (IUNIT, ERR = 110)IHOLD
C
C  WRITE THE COORDINATE NAMES AND ELEMENT NAMES
C
      WRITE (IUNIT)XNAME, YNAME
      IF (NUMMAT.LE.IGUESS)WRITE (IUNIT) (ENAME (I), I = 1, NUMMAT)
C
C  SUCCESSFULL WRITE COMPLETED
C
      CALL MESAGE (' ')
      CALL MESAGE (' ')
      CALL MESAGE ('GENESIS OUTPUT FILE SUCCESSFULLY WRITTEN')
      CALL MESAGE (' ')
      ERR = .FALSE.
      RETURN
C
C  ERR DURING WRITE PROBLEMS
C
  110 CONTINUE
      CALL MESAGE ('ERR DURING WRITE TO OUTPUT FILE')
      CALL MESAGE ('      - NO FILE SAVED -')
      RETURN
C
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
      DOUBLE PRECISION FUNCTION erfc1(ind,x)
C-----------------------------------------------------------------------
C         EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
C
C          ERFC1(IND,X) = ERFC(X)            IF IND = 0
C          ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
C-----------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
      INTEGER ind
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ax,bot,c,e,t,top,w
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION exparg
      EXTERNAL exparg
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,dble,exp
C     ..
C     .. Data statements ..
C-------------------------
C-------------------------
C-------------------------
C-------------------------
      DATA c/.564189583547756D0/
      DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/,
     +     a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/,
     +     a(5)/.128379167095513D+00/
      DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/,
     +     b(3)/.375795757275549D+00/
      DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/,
     +     p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/,
     +     p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/,
     +     p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/
      DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/,
     +     q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/,
     +     q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/,
     +     q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/
      DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/,
     +     r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/,
     +     r(5)/2.82094791773523D-01/
      DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/,
     +     s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/
C     ..
C     .. Executable Statements ..
C-------------------------
C
C                     ABS(X) .LE. 0.5
C
      ax = abs(x)
      IF (ax.GT.0.5D0) GO TO 10
      t = x*x
      top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0
      bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0
      erfc1 = 0.5D0 + (0.5D0-x* (top/bot))
      IF (ind.NE.0) erfc1 = exp(t)*erfc1
      RETURN
C
C                  0.5 .LT. ABS(X) .LE. 4
C
   10 IF (ax.GT.4.0D0) GO TO 20
      top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+
     +      p(7))*ax + p(8)
      bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+
     +      q(7))*ax + q(8)
      erfc1 = top/bot
      GO TO 40
C
C                      ABS(X) .GT. 4
C
   20 IF (x.LE.-5.6D0) GO TO 60
      IF (ind.NE.0) GO TO 30
      IF (x.GT.100.0D0) GO TO 70
      IF (x*x.GT.-exparg(1)) GO TO 70
C
   30 t = (1.0D0/x)**2
      top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5)
      bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0
      erfc1 = (c-t*top/bot)/ax
C
C                      FINAL ASSEMBLY
C
   40 IF (ind.EQ.0) GO TO 50
      IF (x.LT.0.0D0) erfc1 = 2.0D0*exp(x*x) - erfc1
      RETURN

   50 w = dble(x)*dble(x)
      t = w
      e = w - dble(t)
      erfc1 = ((0.5D0+ (0.5D0-e))*exp(-t))*erfc1
      IF (x.LT.0.0D0) erfc1 = 2.0D0 - erfc1
      RETURN
C
C             LIMIT VALUE FOR LARGE NEGATIVE X
C
   60 erfc1 = 2.0D0
      IF (ind.NE.0) erfc1 = 2.0D0*exp(x*x)
      RETURN
C
C             LIMIT VALUE FOR LARGE POSITIVE X
C                       WHEN IND = 0
C
   70 erfc1 = 0.0D0
      RETURN

      END

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


text:
PROGRAM problematic_if

! IF (printlev_loc >=2)  WRITE(numout,*),'tau_outflow, coeff_rel = ', tau_outflow, coeff_rel
  IF (printlev_loc >=2)  WRITE(numout,*) 'tau_outflow, coeff_rel = ', tau_outflow, coeff_rel

END

max_stars_repo_name: mbdevpl/open-fortran-parser-xml

max_issues_repo_name: mbdevpl/open-fortran-parser-xml

max_forks_repo_name: mbdevpl/open-fortran-parser-xml