text:
C***********************************************************************
C Copyright (C) 1994-2016 Lawrence Livermore National Security, LLC.
C LLNL-CODE-425250.
C All rights reserved.
C
C This file is part of Silo. For details, see silo.llnl.gov.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions
C are met:
C
C    * Redistributions of source code must retain the above copyright
C      notice, this list of conditions and the disclaimer below.
C    * Redistributions in binary form must reproduce the above copyright
C      notice, this list of conditions and the disclaimer (as noted
C      below) in the documentation and/or other materials provided with
C      the distribution.
C    * Neither the name of the LLNS/LLNL 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  LAWRENCE
C LIVERMORE  NATIONAL SECURITY, LLC,  THE U.S.  DEPARTMENT OF  ENERGY OR
C CONTRIBUTORS BE LIABLE FOR  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
C EXEMPLARY, OR  CONSEQUENTIAL DAMAGES  (INCLUDING, BUT NOT  LIMITED TO,
C PROCUREMENT OF  SUBSTITUTE GOODS  OR SERVICES; LOSS  OF USE,  DATA, OR
C PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
C LIABILITY, WHETHER  IN CONTRACT, STRICT LIABILITY,  OR TORT (INCLUDING
C NEGLIGENCE OR  OTHERWISE) ARISING IN  ANY WAY OUT  OF THE USE  OF THIS
C SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
C This work was produced at Lawrence Livermore National Laboratory under
C Contract No.  DE-AC52-07NA27344 with the DOE.
C
C Neither the  United States Government nor  Lawrence Livermore National
C Security, LLC nor any of  their employees, makes any warranty, express
C or  implied,  or  assumes  any  liability or  responsibility  for  the
C accuracy, completeness,  or usefulness of  any information, apparatus,
C product, or  process disclosed, or  represents that its use  would not
C infringe privately-owned rights.
C
C Any reference herein to  any specific commercial products, process, or
C services by trade name,  trademark, manufacturer or otherwise does not
C necessarily  constitute or imply  its endorsement,  recommendation, or
C favoring  by  the  United  States  Government  or  Lawrence  Livermore
C National Security,  LLC. The views  and opinions of  authors expressed
C herein do not necessarily state  or reflect those of the United States
C Government or Lawrence Livermore National Security, LLC, and shall not
C be used for advertising or product endorsement purposes.
C***********************************************************************

      program testmat

C...This program tests the use of the UCD Fortran jacket routines, plus
C...the underlying UCD C routines.
      include 'silo.inc'
      integer driver, nargs
      character*256 cloption

      driver = DB_PDB
      nargs = iargc()
      call getarg(1, cloption)
      if (cloption .eq. "DB_HDF5") then
          driver = DB_HDF5
      end if

      call writeit ("quad_mat.silo", driver)
      print *,'Created quad_mat.silo'

      stop
      end

      subroutine writeit (fname, drvr)
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C Writes the following 2x3 quadmesh (note colmajor ordering)
C
C    ^Y axis
C    |
C   3-     1--------4--------7-------10
C    |     |        |        |       |
C    |     |   1    |   3    |   5   |
C   2-     2--------5--------8-------11
C    |     |        |        |       |
C    |     |   2    |   4    |   6   |
C   1-     3--------6--------9-------12
C    |
C   0------1--------2--------3-------4-->X axis
C
C with 2-material composition. . .
C
C          +--------+-------/+-------+
C          |   1    |  1   / |       |
C          |        |     / 2|   2   |
C          +--------+----/---+-------+
C          |   1    | 1 /    |       |
C          |        |  /  2  |   2   |
C          +--------+-/------+-------+
C
C using 2D Fortran arrays A(#rows,#cols) stored column major (e.g.
C row index varies fastest as you march through A's memory). Note
C that #rows is size in Y direction while #cols is size in X
C direction. So, A is dimensioned (nY,nX). In Fortran, the left-most
C dimension in an array declaration varies fastest when marching
C through array elements in memory. However, the 'dims' array in
C Silo's DBPut calls is arranged such that dims[0] is the size in
C X, dims[1] the size in Y and dims[2], if applicable, the size in Z.
C So, although we have declarations for A(3,4), the dims array is
C passed in as dims[0] = 4, dims[1] = 3
C----------------------------------------------------------------------
Cc      implicit double precision (a-h,o-z)

C...Include SILO definitions.

      include 'silo.inc'

      character*(*) fname
      integer       drvr

      integer  dbid, dbstat
      integer  tcycle
      real     ttime
      real     x(3,4), y(3,4)
C...For co-linear case
      real xcl(4), ycl(3)

      integer  matlist(2,3), matlistcl(3,2)
      integer  dims(2), dimscl(2)
      integer  zdims(2), zdimscl(2)

      integer  matnos(2)
      integer  mix_next(4)
      integer  mix_zone(4)
      integer  mix_mat(4)
      real     mix_vf(4)

      data matlist/1, 1, -1, -3, 2, 2/
C      data matlistcl/1, -1, 2, 1, -3, 2/
      data matlistcl/1, -3, 2, 1, -1, 2/
C      data dims/4,3/
      data dims/3,4/
      data dimscl/4,3/
C      data zdims/3,2/
      data zdims/2,3/
      data zdims/2,3/
      data zdimscl/3,2/
      data x/1.0,1.0,1.0,2.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,4.0/
      data y/3.0,2.0,1.0,3.0,2.0,1.0,3.0,2.0,1.0,3.0,2.0,1.0/
      data xcl/1.0,2.0,3.0,4.0/
      data ycl/1.0,2.0,3.0/
      data mix_next/2,0,4,0/
      data mix_zone/2,2,4,4/
      data mix_mat/1,2,1,2/
      data mix_vf/0.7,0.3,0.4,0.6/
      data matnos  /1, 2/
      data ttime   /1.23/
      data tcycle  /200/

      ttime = 2.345

      err = dbshowerrors(DB_ABORT)

      err = dbcreate(fname, len(fname), 0, DB_LOCAL,
     .               "foo-info", 8, drvr, dbid)
      if (err .ne. 0) then
         print *,'Error from dbcreate'
         stop
      endif

      err = dbputqm(dbid, "qmesh", 5, "x", 1, "y", 1,
     . DB_F77NULLSTRING, 0,
     . x, y, DB_F77NULL, dims, 2, DB_FLOAT, DB_NONCOLLINEAR,
     . DB_F77NULL, dbstat)

      err = dbputqm(dbid, "qmeshcl", 7, "xcl", 3, "ycl", 3,
     . DB_F77NULLSTRING, 0,
     . xcl, ycl, DB_F77NULL, dimscl, 2, DB_FLOAT, DB_COLLINEAR,
     . DB_F77NULL, dbstat)

      err = dbputmat (dbid, "mat", 3, "qmesh", 5,
     . 2, matnos, matlist, zdims, 2,
     . mix_next, mix_mat, mix_zone, mix_vf, 4,
     . DB_FLOAT, DB_F77NULL, dbstat)

      err = dbputmat (dbid, "matcl", 5, "qmeshcl", 7,
     . 2, matnos, matlistcl, zdimscl, 2,
     . mix_next, mix_mat, mix_zone, mix_vf, 4,
     . DB_FLOAT, DB_F77NULL, dbstat)

C...Close file before quitting.

      err = dbclose(dbid)

      return
      end

max_stars_repo_name: EMinsight/Silo

max_issues_repo_name: EMinsight/Silo

max_forks_repo_name: EMinsight/Silo


text:
      SUBROUTINE YGRATL
C
C          ------------------------------------------------
C          ROUTINE NO. ( 230)   VERSION (A8.6)    17:JUN:88
C          ------------------------------------------------
C
C          THIS DRAWS A Y-GRATICULE WITH NO ANNOTATION,
C          SETTING THE AXIS INTERVALS AUTOMATICALLY.
C
C
C          [Y-LOGARITHMIC MARKING].
C
C
      REAL    RDATA(1)
      INTEGER IDATA(1)
C
      COMMON /T0ADIY/ DIVLY,NTIKLY,NTIKHY
      COMMON /T0APOS/ AXPOSX,AXPOSY
      COMMON /T0ASKY/ NSKIPY,NDECSY
      COMMON /T0NOTA/ NOTATA
      COMMON /T0PPOS/ XPLOT0,YPLOT0
      COMMON /T0TRAC/ IPRINT
      COMMON /T0WNDO/ X1WND0,X2WND0,Y1WND0,Y2WND0
      COMMON /T3LIMS/ IMAXI,RMAXI,RMINI
C
      DATA RDATA /0.0/
C
C
      IF (IPRINT.EQ.1) CALL G0MESG(108,0)
C
C          IF THE WINDOW AREA IS WRONG, NOTHING MORE IS DONE.
C          THE ANNOTATION MARKER IS SAVED AND THEN SET = 0
C          SO THAT ANNOTATION IS NOT DONE WITH A GRATICULE.
C          THE CURRENT PLOTTING POSITION IS THEN SAVED.
C
      IF (ABS(Y2WND0-Y1WND0).LT.RMINI) RETURN
C
      NTATSV= NOTATA
      NOTATA= 0
      XHERE= XPLOT0
      YHERE= YPLOT0
C
C          THE Y-AXIS GRATICULE LINES ARE DRAWN.
C          IF THERE ARE MORE THAN 37 SUB-INTERVALS IN THE Y-AXIS,
C          THESE ARE DRAWN AS EDGE TICKS; OTHERWISE THEY ARE DONE
C          AS FULL GRATICULE LINES (JUST AS THE DECADE INTERVALS).
C
      IDATA(1)= 1
      CALL G3LINK(0,14,-1,IDATA,RDATA)
      CALL G0DIVL(2)
      IPOS= NTIKLY
      IDEC= NSKIPY
      NDECS= NDECSY
      NTICKS= NTIKHY
      IF (NTICKS.LE.37) GO TO 1
      IF (IPOS.EQ.1) IDEC= IDEC-1
C
      DO 100 ITICK= 1,NDECS
        YPOS= DIVLY*(10.0**(ITICK+IDEC))
        CALL POSITN(X1WND0,YPOS)
        CALL JOIN(X2WND0,YPOS)
  100 CONTINUE
C
      AXYSAV= AXPOSY
      AXPOSY= X1WND0
      CALL G0TICK
      CALL G0PLYL(0)
      AXPOSY= X2WND0
      CALL G0TICK
      CALL G0PLYL(0)
      AXPOSY= AXYSAV
      GO TO 2
C
    1 IF (NTICKS.LE.0) GO TO 2
C
      DO 200 ITICK= 1,NTICKS
        YPOS= DIVLY*IPOS*(10.0**IDEC)
        CALL POSITN(X1WND0,YPOS)
        CALL JOIN(X2WND0,YPOS)
        IPOS= IPOS+1
        IF (IPOS.LT.10) GO TO 200
C
        IPOS= 1
        IDEC= IDEC+1
  200 CONTINUE
C
    2 IDATA(1)= 0
      CALL G3LINK(0,14,-1,IDATA,RDATA)
C
C          THE ANNOTATION MARKER AND THE CURRENT
C          PLOTTING POSITION ARE THEN RESTORED.
C
      NOTATA= NTATSV
      CALL POSITN(XHERE,YHERE)
C
      RETURN
      END

ext: f

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
subroutine SHCilmToVector(cilm, vector, lmax)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This subroutine will convert an array of spherical harmonic
!	coefficients Cilm to an ordered 1D vector.
!
!	Calling Parameters
!		IN
!			cilm	Input spherical harmonic coefficients with dimension
!				cilm(2, lmax+1, lmax+1).
!			lmax	Maximum spherical harmonic degree of the input coefficients
!		OUT
!			vector	1D vector of ordered spherical harmonic coefficients with
!				dimension (lmax+1)**2. The ordering is described in YilmIndex.
!
!
!	Written by Mark Wieczorek (August 2009)
!
!	Copyright (c) 2009, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	implicit none
	real*8, intent(in) ::	cilm(:,:,:)
	real*8, intent(out) ::	vector(:)
	integer, intent(in) ::	lmax
	integer ::		i, l, m

	if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax+1 .or. size(cilm(1,1,:)) < lmax+1) then
		print*, "Error --- SHCilmToVector"
		print*, "CILM must be dimensioned as (2, LMAX+1, LMAX + 1)."
		print*, "LMAX = ", lmax
		print*, "Dimension of CILM = ", size(cilm(:,1,1)),  size(cilm(1,:,1)), size(cilm(1,1,:))
		stop
	endif

	if (size(vector(:)) < (lmax+1)**2) then
		print*, "Error --- SHCilmToVector"
		print*, "VECTOR must be have dimension (LMAX+1)**2."
		print*, "LMAX = ", lmax
		print*, "Dimension of VECTOR = ", size(vector(:))
		stop
	endif

	if (lmax < 0) then
		print*, "Error --- SHCilmToVector"
		print*, "LMAX must be positive."
		print*, "LMAX = ", lmax
		stop
	endif

	i=0
	do l=0, lmax
		do m=0, l
			i = i + 1
			vector(i) = cilm(1, l+1, m+1)
		enddo

		do m=1, l
			i = i + 1
			vector(i) = cilm(2, l+1, m+1)
		enddo
	enddo

end subroutine SHCilmToVector

subroutine SHVectorToCilm(vector, cilm, lmax)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This subroutine will convert a 1D ordered vector of spherical harmonic
!	coefficients to a 3D array Cilm(i, l, m).
!
!	Calling Parameters
!		IN
!			vector	1D vector of ordered spherical harmonic coefficients with
!				dimension (lmax+1)**2. The ordering is described in YilmIndex.
!			lmax	Maximum spherical harmonic degree of the input coefficients
!		OUT
!			cilm	Output spherical harmonic coefficients with dimension
!				cilm(2, lmax+1, lmax+1).
!
!
!	Written by Mark Wieczorek (August 2009)
!
!	Copyright (c) 2009, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	implicit none
	real*8, intent(out) ::	cilm(:,:,:)
	real*8, intent(in) ::	vector(:)
	integer, intent(in) ::	lmax
	integer ::		k, i, l, m

	if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax+1 .or. size(cilm(1,1,:)) < lmax+1) then
		print*, "Error --- SHVectorToCilm"
		print*, "CILM must be dimensioned as (2, LMAX+1, LMAX + 1)."
		print*, "LMAX = ", lmax
		print*, "Dimension of CILM = ", size(cilm(:,1,1)),  size(cilm(1,:,1)), size(cilm(1,1,:))
		stop
	endif

	if (size(vector(:)) < (lmax+1)**2) then
		print*, "Error --- SHVectorToCilm"
		print*, "VECTOR must be have dimension (LMAX+1)**2."
		print*, "LMAX = ", lmax
		print*, "Dimension of VECTOR = ", size(vector(:))
		stop
	endif

	if (lmax < 0) then
		print*, "Error --- SHVectorToCilm"
		print*, "LMAX must be positive."
		print*, "LMAX = ", lmax
		stop
	endif

	l = 0
	m = 0
	i = 1

	cilm(1,1,1) = vector(1)

	i = 2

	do k=2, (lmax+1)**2
		m = m + 1
		if (m > l .and. i == 1) then
			i = 2
			m = 1
		elseif (m > l .and. i == 2) then
			l = l + 1
			m = 0
			i = 1
		endif

		cilm(i,l+1,m+1) = vector(k)

	enddo

end subroutine SHVectorToCilm

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.

C ============================================================================
C
      SUBROUTINE APPSLN(T, NZ, NY, NU, NELE, NCOL, TI, X,
     1     ZVAL, DZVAL, YVAL, UVAL)
C !DEC$ ATTRIBUTES DLLEXPORT :: APPSLN
C
C     $Id: appsln.f,v 1.1 2002/05/02 18:52:15 andreasw Exp $
C
C     Subroutine for evaluating the DAE variables at given time T
C
C     Author:  Andreas Waechter
C              c/o Group of Larry Biegler
C              Department of Chemical Engineering
C              Carnegie Mellon University
C              Pittsburgh, PA
C
C     Version 1.0    ( 10-30-00 )
C
C     Description:
C     ============
C
C       evaluate Z, Y and U on time T based interpolation;
C          results in ZVAL, YVAL, UVAL
C
C     Parameters:
C
C     T        Input    time on which variable are to evaluated
C     NZ       Input    number of differential variables
C     NY       Input    number of algebraic variables
C     NU       Input    number of control variables
C     NELE     Input    number of elements
C     NCOL     Input    number of collocation points
C     TI       Input    boundaries of finite elements
C     X        Input    coefficients of collocation polynomials
C                         (for ordering see dae2nlp.f)
C     ZVAL     Output   value of differential variables at T
C     DZVAL    Output   value of derivative of differential variables at T
C     YVAL     Output   value of algebraic variables at T
C     UVAL     Output   value of control variables at T
C
      implicit none
      double precision T
      integer NZ, NY, NU, NELE, NCOL
      double precision TI(NELE+1)
      double precision X(NZ + NELE*(NCOL*(NZ+NY+NU)+NZ))
      double precision ZVAL(NZ)
      double precision DZVAL(NZ)
      double precision YVAL(NY)
      double precision UVAL(NU)

      include 'DAE2NLP.INC'
C !DEC$ ATTRIBUTES DLLEXPORT :: /DAENLP/

      integer iele, il, iu, lw, lu
      double precision theta(7), omega(7), dummy
C
C     Find element in which T is located (binary search)
C
      if( TI(2).ge.T ) then
         iele = 1
      elseif( TI(NELE).lt.T ) then
         iele = NELE
      else
         il = 2
         iu = NELE
 100     continue
         if( iu-il.gt.1 ) then
            iele = il+(iu-il)/2
            if( TI(iele).lt.T ) then
               il = IELE
            else
               iu = iele
            endif
            goto 100
         endif
         iele = il
      endif
C
C     Compute values
C
      lw = 1 + (iele-1)*(NCOL*(NZ+NY)+NZ)
      lu = 1 + NELE*(NCOL*(NZ+NY)+NZ) + NZ + (iele-1)*NCOL*NU
      call APPROX (3, NZ, NY, NU, NELE, NCOL, COEF, X(lw), X(lu), TI,
     1     iele, T, omega, theta, ZVAL, YVAL, UVAL)
      call APPROX (4, NZ, NY, NU, NELE, NCOL, COEF, X(lw), X(lu), TI,
     1     iele, T, omega, theta, DZVAL, dummy ,dummy)
C
C     That's it
C
      return
      end

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
c fishpk27 from portlib                                  12/30/83
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                        f i s h p a k                          *
c     *                                                               *
c     *                                                               *
c     *     a package of fortran subprograms for the solution of      *
c     *                                                               *
c     *      separable elliptic partial differential equations        *
c     *                                                               *
c     *                  (version 3.1 , october 1980)                  *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *        john adams, paul swarztrauber and roland sweet         *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the national center for atmospheric research          *
c     *                                                               *
c     *                boulder, colorado  (80307)  u.s.a.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the national science foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c
c          program to illustrate the use of subroutine hwscyl to solve
c     the equation
c
c     (1/r)(d/dr)(r*(du/dr)) + (d/dz)(du/dz)
c
c     = (2*r*z)**2*(4*z**2 + 3*r**2)
c
c     on the rectangle 0 .lt. r .lt. 1, 0 .lt. z .lt. 1 with the
c     boundary conditions
c
c     u(0,z) unspecified
c                                            0 .le. z .le. 1
c     (du/dr)(1,z) = 4*z**4
c
c     and
c
c     (du/dz)(r,0) = 0
c                                            0 .le. r .le. 1
c     (du/dz)(r,1) = 4*r**4 .
c
c          the r-interval will be divided into 50 panels and the
c     z-interval will be divided into 100 panels.
c
      program thwscyl
      dimension       f(75,105)  ,bda(101)   ,bdb(101)   ,bdc(51)    ,
     1                bdd(51)    ,w(1373)    ,r(51)      ,z(101)
c
c     from dimension statement we get value of idimf.  also note that w
c     is dimensioned 4*(n+1) + (13 + int(log2(n+1)))*(m+1) .
c
      idimf = 75
      a = 0.
      b = 1.
      m = 50
      mbdcnd = 6
      c = 0.
      d = 1.
      n = 100
      nbdcnd = 3
      elmbda = 0.
c
c     auxiliary quantities.
c
      mp1 = m+1
      np1 = n+1
c
c     generate and store grid points for the purpose of computing
c     boundary data and the right side of the poisson equation.
c
      do 101 i=1,mp1
         r(i) = float(i-1)/50.
  101 continue
      do 102 j=1,np1
         z(j) = float(j-1)/100.
  102 continue
c
c     generate boundary data.
c
      do 103 j=1,np1
         bdb(j) = 4.*z(j)**4
  103 continue
      do 104 i=1,mp1
         bdc(i) = 0.
         bdd(i) = 4.*r(i)**4
  104 continue
c
c     bda is a dummy variable.
c
c
c     generate right side of equation.
c
      do 106 i=1,mp1
         do 105 j=1,np1
            f(i,j) = 4.*r(i)**2*z(j)**2*(4.*z(j)**2+3.*r(i)**2)
  105    continue
  106 continue
      call hwscyl (a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,
     1             idimf,pertrb,ierror,w)
c
c     compute discretization error by minimizing over all a the function
c     norm(f(i,j) - a*1 - u(r(i),z(j))).  the exact solution is
c                u(r,z) = (r*z)**4 + arbitrary constant.
c
      x = 0.
      do 108 i=1,mp1
         do 107 j=1,np1
            x = x+f(i,j)-(r(i)*z(j))**4
  107    continue
  108 continue
      x = x/float(np1*mp1)
      do 110 i=1,mp1
         do 109 j=1,np1
            f(i,j) = f(i,j)-x
  109    continue
  110 continue
      err = 0.
      do 112 i=1,mp1
         do 111 j=1,np1
            x = abs(f(i,j)-(r(i)*z(j))**4)
            if (x .gt. err) err = x
  111    continue
  112 continue
      print 1001 , ierror,pertrb,err,w(1)
      stop
c
 1001 format (1h1,20x,25hsubroutine hwscyl example///
     1        10x,46hthe output from the ncar control data 7600 was//
     2        32x,10hierror = 0/32x,20hpertrb = 2.26734e-04/
     3        18x,34hdiscretization error = 3.73672e-04/
     4        12x,33hrequired length of w array = 1118//
     5        10x,32hthe output from your computer is//
     6        32x,8hierror =,i2/32x,8hpertrb =,e12.5/
     7        18x,22hdiscretization error =,e12.5/
     8        12x,28hrequired length of w array =,f5.0)
c
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text: Delaine Eastin is a former California Assembly representative from Alameda County, and former State Superintendent of Public Instruction (19952003) who currently lives in Davis.  More information about Ms. Eastin can be found at her http://en.wikipedia.org/wiki/Delaine_Eastin wikipedia entry.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
great_ne(q1,m1).
great_ne(m1,a1).
great_ne(b1,c1).
great_ne(v1,ii1).
great_ne(z1,ff1).
great_ne(y1,w1).
great_ne(z1,p1).
great_ne(p1,m1).
great_ne(ii1,w1).
great_ne(v1,c1).
great_ne(v1,b1).
great_ne(cc1,ee1).
great_ne(bb1,w1).
great_ne(z1,ii1).
great_ne(cc1,dd1).
great_ne(p1,x1).
great_ne(m1,ii1).
great_ne(p1,b1).
great_ne(dd1,v1).
great_ne(aa1,y1).
great_ne(aa1,i1).
great_ne(kk1,hh1).
great_ne(k1,x1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
MODULE MF_InitializationModule_NonRelativistic_TABLE

  ! --- AMReX Modules ---

  USE amrex_fort_module,       ONLY: &
    AR => amrex_real
  USE amrex_box_module,        ONLY: &
    amrex_box
  USE amrex_geometry_module,   ONLY: &
    amrex_geometry
  USE amrex_multifab_module,   ONLY: &
    amrex_multifab,     &
    amrex_mfiter,       &
    amrex_mfiter_build, &
    amrex_mfiter_destroy
  USE amrex_parallel_module,   ONLY: &
    amrex_parallel_ioprocessor
  USE amrex_parmparse_module,  ONLY: &
    amrex_parmparse,       &
    amrex_parmparse_build, &
    amrex_parmparse_destroy

  ! --- thornado Modules ---

  USE UnitsModule,             ONLY: &
    Gram,       &
    Kilogram,   &
    Centimeter, &
    Meter,      &
    Kilometer,  &
    Erg,        &
    Joule,      &
    Second,     &
    Kelvin,     &
    SpeedOfLight
  USE ProgramHeaderModule,     ONLY: &
    nDOFX,   &
    nX,      &
    nNodesX, &
    swX,     &
    nDimsX
  USE ReferenceElementModuleX, ONLY: &
    NodeNumberTableX
  USE MeshModule,              ONLY: &
    MeshType,    &
    CreateMesh,  &
    DestroyMesh, &
    NodeCoordinate
  USE GeometryFieldsModule,    ONLY: &
    nGF,          &
    iGF_Gm_dd_11, &
    iGF_Gm_dd_22, &
    iGF_Gm_dd_33
  USE FluidFieldsModule,       ONLY: &
    nCF,    &
    iCF_D,  &
    iCF_S1, &
    iCF_S2, &
    iCF_S3, &
    iCF_E,  &
    iPF_Ne, &
    nPF,    &
    iPF_D,  &
    iPF_V1, &
    iPF_V2, &
    iPF_V3, &
    iPF_E,  &
    iCF_Ne, &
    nAF,    &
    iAF_P,  &
    iAF_E,  &
    iAF_T,  &
    iAF_Ye
  USE Euler_UtilitiesModule,   ONLY: &
    ComputeConserved_Euler
  USE EquationOfStateModule,   ONLY: &
    ComputePressureFromPrimitive,   &
    ComputeTemperatureFromPressure, &
    ComputeThermodynamicStates_Primitive
  USE UtilitiesModule,         ONLY: &
    Locate
  USE Euler_ErrorModule,       ONLY: &
    DescribeError_Euler

  ! --- Local Modules ---

  USE InputParsingModule,      ONLY: &
    nLevels, &
    xL,      &
    xR,      &
    Gamma_IDEAL, &
    UseTiling

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: MF_InitializeFields_NonRelativistic_TABLE

  REAL(AR), PARAMETER :: Zero   = 0.0_AR
  REAL(AR), PARAMETER :: Half   = 0.5_AR
  REAL(AR), PARAMETER :: One    = 1.0_AR
  REAL(AR), PARAMETER :: Pi     = ACOS( -1.0_AR )
  REAL(AR), PARAMETER :: TwoPi  = 2.0_AR * Pi
  REAL(AR), PARAMETER :: FourPi = 4.0_AR * Pi

CONTAINS

  SUBROUTINE MF_InitializeFields_NonRelativistic_TABLE &
    ( ProgramName, MF_uGF, MF_uCF, GEOM )

    CHARACTER(LEN=*),     INTENT(in)    :: ProgramName
    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)
    TYPE(amrex_geometry), INTENT(in)    :: GEOM  (0:nLevels-1)

    IF( amrex_parallel_ioprocessor() )THEN

      WRITE(*,*)
      WRITE(*,'(A4,A,A)') '', 'Initializing: ', TRIM( ProgramName )

    END IF

    SELECT CASE ( TRIM( ProgramName ) )

      CASE( 'Advection_TABLE' )

        CALL InitializeFields_Advection_TABLE( MF_uGF, MF_uCF )

      CASE( 'Sod_TABLE' )

        CALL InitializeFields_Sod_TABLE( MF_uGF, MF_uCF )

      CASE( 'SphericalSod_TABLE' )

        CALL InitializeFields_SphericalSod_TABLE( MF_uGF, MF_uCF )

      CASE( 'CylindricalSod_TABLE' )

        CALL InitializeFields_CylindricalSod_TABLE( MF_uGF, MF_uCF )

      CASE( 'Jet_TABLE' )

        CALL InitializeFields_Jet_TABLE( MF_uGF, MF_uCF )

      CASE( 'Implosion_TABLE' )

        CALL InitializeFields_Implosion_TABLE( MF_uGF, MF_uCF )

      CASE DEFAULT

        IF( amrex_parallel_ioprocessor() )THEN

          WRITE(*,*)
          WRITE(*,'(4x,A,A)') 'Unknown Program: ', TRIM( ProgramName )
          WRITE(*,'(4x,A)')   'Valid Options:'
          WRITE(*,'(6x,A)')     'Advection_TABLE'
          WRITE(*,'(6x,A)')     'Sod_TABLE'
          WRITE(*,'(6x,A)')     'SphericalSod_TABLE'
          WRITE(*,'(6x,A)')     'CylindricalSod_TABLE'
          WRITE(*,'(6x,A)')     'Jet_TABLE'
          WRITE(*,'(6x,A)')     'Implosion_TABLE'

        END IF

        CALL DescribeError_Euler( 99 )

    END SELECT

  END SUBROUTINE MF_InitializeFields_NonRelativistic_TABLE

  SUBROUTINE InitializeFields_Advection_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1
    REAL(AR)       :: X1
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    ! --- Problem-dependent Parameters ---
    REAL(AR), PARAMETER :: D_0 = 1.0e12_AR * Gram / Centimeter**3
    REAL(AR), PARAMETER :: Amp = 1.0e11_AR * Gram / Centimeter**3
    REAL(AR), PARAMETER :: L   = 1.0e02_AR * Kilometer

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

            uPF_K(iNodeX,iPF_D ) &
              = D_0 + Amp * SIN( TwoPi * X1 / L )
            uPF_K(iNodeX,iPF_V1) &
              = 0.1_AR * SpeedOfLight
            uPF_K(iNodeX,iPF_V2) &
              = 0.0_AR * Kilometer / Second
            uPF_K(iNodeX,iPF_V3) &
              = 0.0_AR * Kilometer / Second
            uAF_K(iNodeX,iAF_P ) &
              = 1.0e-2_AR * D_0 * SpeedOfLight**2
            uAF_K(iNodeX,iAF_Ye) &
              = 0.3_AR

          END DO

          CALL ComputeTemperatureFromPressure &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                   uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

          CALL ComputeThermodynamicStates_Primitive &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                   uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                   uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

          CALL ComputeConserved_Euler &
                 ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                   uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                   uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                   uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                   uGF_K(:,iGF_Gm_dd_11), &
                   uGF_K(:,iGF_Gm_dd_22), &
                   uGF_K(:,iGF_Gm_dd_33), &
                   uAF_K(:,iAF_P) )

          uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
            = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_Advection_TABLE

  SUBROUTINE InitializeFields_Sod_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1
    REAL(AR)       :: X1
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

            IF( X1 .LE. Zero ) THEN

              uPF_K(iNodeX,iPF_D ) = 1.0e12_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero      * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.4_AR

            ELSE

              uPF_K(iNodeX,iPF_D ) = 1.25e11_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e31_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.3_AR

            END IF

          END DO

          CALL ComputeTemperatureFromPressure &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                   uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

          CALL ComputeThermodynamicStates_Primitive &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                   uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                   uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

          CALL ComputeConserved_Euler &
                 ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                   uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                   uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                   uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                   uGF_K(:,iGF_Gm_dd_11), &
                   uGF_K(:,iGF_Gm_dd_22), &
                   uGF_K(:,iGF_Gm_dd_33), &
                   uAF_K(:,iAF_P) )

          uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
            = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_Sod_TABLE

  SUBROUTINE InitializeFields_SphericalSod_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1
    REAL(AR)       :: X1
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

            IF( X1 .LE. 5.0_AR * Kilometer ) THEN

              uPF_K(iNodeX,iPF_D ) = 1.0e12_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero      * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.4_AR

            ELSE

              uPF_K(iNodeX,iPF_D ) = 1.25e11_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e31_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.4_AR

            END IF

          END DO

          CALL ComputeTemperatureFromPressure &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                   uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

          CALL ComputeThermodynamicStates_Primitive &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                   uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                   uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

          CALL ComputeConserved_Euler &
                 ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                   uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                   uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                   uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                   uGF_K(:,iGF_Gm_dd_11), &
                   uGF_K(:,iGF_Gm_dd_22), &
                   uGF_K(:,iGF_Gm_dd_33), &
                   uAF_K(:,iAF_P) )

          uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
            = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_SphericalSod_TABLE

  SUBROUTINE InitializeFields_CylindricalSod_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1, iNodeX2
    REAL(AR)       :: X1, X2
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
            X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

            IF( SQRT( X1**2 + X2**2 ) <= 5.0_AR * Kilometer ) THEN

              uPF_K(iNodeX,iPF_D ) = 1.0e12_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero      * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero      * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.4_AR

            ELSE

              uPF_K(iNodeX,iPF_D ) = 1.25e11_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e31_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.4_AR

            END IF

          END DO

          CALL ComputeTemperatureFromPressure &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                   uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

          CALL ComputeThermodynamicStates_Primitive &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                   uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                   uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

          CALL ComputeConserved_Euler &
                 ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                   uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                   uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                   uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                   uGF_K(:,iGF_Gm_dd_11), &
                   uGF_K(:,iGF_Gm_dd_22), &
                   uGF_K(:,iGF_Gm_dd_33), &
                   uAF_K(:,iAF_P) )

          uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
            = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_CylindricalSod_TABLE

  SUBROUTINE InitializeFields_Jet_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1, iNodeX2
    REAL(AR)       :: X1, X2
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
            X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

            IF( X1 .LE. Half * Kilometer .AND. X2 .LE. Half * Kilometer )THEN

              !SW
              uPF_K(iNodeX,iPF_D ) = 0.80e12_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.3e0_AR

            ELSE IF( X1 .LE. Half * Kilometer &
                       .AND. X2 .GT. Half * Kilometer )THEN
              !NW
              uPF_K(iNodeX,iPF_D ) = 1.0e12_AR  * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = 7.275e4_AR * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.3e0_AR

            ELSE IF( X1 .GT. Half * Kilometer &
                       .AND. X2 .GT. Half * Kilometer )THEN

              !NE
              uPF_K(iNodeX,iPF_D ) = 0.5313e12_AR * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero         * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = Zero         * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero         * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 0.4e32_AR    * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.3e0_AR

            ELSE

              !SE
              uPF_K(iNodeX,iPF_D ) = 1.0e12_AR  * Gram / Centimeter**3
              uPF_K(iNodeX,iPF_V1) = Zero       * Kilometer / Second
              uPF_K(iNodeX,iPF_V2) = 7.275e4_AR * Kilometer / Second
              uPF_K(iNodeX,iPF_V3) = Zero       * Kilometer / Second
              uAF_K(iNodeX,iAF_P ) = 1.0e32_AR  * Erg / Centimeter**3
              uAF_K(iNodeX,iAF_Ye) = 0.3e0_AR

            END IF

          END DO

        CALL ComputeTemperatureFromPressure &
               ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                 uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

        CALL ComputeThermodynamicStates_Primitive &
               ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                 uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                 uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

        CALL ComputeConserved_Euler &
               ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                 uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                 uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                 uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                 uGF_K(:,iGF_Gm_dd_11), &
                 uGF_K(:,iGF_Gm_dd_22), &
                 uGF_K(:,iGF_Gm_dd_33), &
                 uAF_K(:,iAF_P) )

        uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
          = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_Jet_TABLE

  SUBROUTINE InitializeFields_Implosion_TABLE( MF_uGF, MF_uCF )

    TYPE(amrex_multifab), INTENT(in)    :: MF_uGF(0:nLevels-1)
    TYPE(amrex_multifab), INTENT(inout) :: MF_uCF(0:nLevels-1)

    ! --- thornado ---
    INTEGER        :: iDim
    INTEGER        :: iX1, iX2, iX3
    INTEGER        :: iNodeX, iNodeX1, iNodeX2
    REAL(AR)       :: X1, X2
    REAL(AR)       :: uGF_K(nDOFX,nGF)
    REAL(AR)       :: uCF_K(nDOFX,nCF)
    REAL(AR)       :: uPF_K(nDOFX,nPF)
    REAL(AR)       :: uAF_K(nDOFX,nAF)
    TYPE(MeshType) :: MeshX(3)

    ! --- AMReX ---
    INTEGER                       :: iLevel
    INTEGER                       :: lo_G(4), hi_G(4)
    INTEGER                       :: lo_F(4), hi_F(4)
    TYPE(amrex_box)               :: BX
    TYPE(amrex_mfiter)            :: MFI
    REAL(AR), CONTIGUOUS, POINTER :: uGF(:,:,:,:)
    REAL(AR), CONTIGUOUS, POINTER :: uCF(:,:,:,:)

    ! --- Problem-dependent Parameters ---
    REAL(AR), PARAMETER :: D_0  = 1.25e13_AR * ( Gram / Centimeter**3 )
    REAL(AR), PARAMETER :: P_0  = 1.0e32_AR  * ( Erg / Centimeter**3 )
    REAL(AR), PARAMETER :: Ye_0 = 1.35e-1_AR
    REAL(AR), PARAMETER :: D_1  = 1.0e14_AR  * ( Gram / Centimeter**3 )
    REAL(AR), PARAMETER :: P_1  = 1.0e33_AR  * ( Erg / Centimeter**3 )
    REAL(AR), PARAMETER :: Ye_1 = 1.5e-1_AR

    uGF_K = Zero
    uCF_K = Zero
    uPF_K = Zero
    uAF_K = Zero

    DO iDim = 1, 3

      CALL CreateMesh &
             ( MeshX(iDim), nX(iDim), nNodesX(iDim), 0, &
               xL(iDim), xR(iDim) )

    END DO

    DO iLevel = 0, nLevels-1

      CALL amrex_mfiter_build( MFI, MF_uGF(iLevel), tiling = UseTiling )

      DO WHILE( MFI % next() )

        uGF => MF_uGF(iLevel) % DataPtr( MFI )
        uCF => MF_uCF(iLevel) % DataPtr( MFI )

        BX = MFI % tilebox()

        lo_G = LBOUND( uGF )
        hi_G = UBOUND( uGF )

        lo_F = LBOUND( uCF )
        hi_F = UBOUND( uCF )

        DO iX3 = BX % lo(3), BX % hi(3)
        DO iX2 = BX % lo(2), BX % hi(2)
        DO iX1 = BX % lo(1), BX % hi(1)

          uGF_K &
            = RESHAPE( uGF(iX1,iX2,iX3,lo_G(4):hi_G(4)), [ nDOFX, nGF ] )

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
            X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

            IF( X1 + X2 .LT. 0.15_AR * Kilometer )THEN

               uPF_K(iNodeX,iPF_D) &
                 = D_0
               uAF_K(iNodeX,iAF_P) &
                 = P_0
               uAF_K(iNodeX,iAF_Ye) &
                 = Ye_0

             ELSE

               uPF_K(iNodeX,iPF_D) &
                 = D_1
               uAF_K(iNodeX,iAF_P) &
                 = P_1
               uAF_K(iNodeX,iAF_Ye) &
                 = Ye_1

             ENDIF

             uPF_K(iNodeX,iPF_V1) &
               = Zero
             uPF_K(iNodeX,iPF_V2) &
               = Zero
             uPF_K(iNodeX,iPF_V3) &
               = Zero

          END DO

          CALL ComputeTemperatureFromPressure &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_P), &
                   uAF_K(:,iAF_Ye), uAF_K(:,iAF_T) )

          CALL ComputeThermodynamicStates_Primitive &
                 ( uPF_K(:,iPF_D ), uAF_K(:,iAF_T ), &
                   uAF_K(:,iAF_Ye), uPF_K(:,iPF_E ), &
                   uAF_K(:,iAF_E ), uPF_K(:,iPF_Ne) )

          CALL ComputeConserved_Euler &
                 ( uPF_K(:,iPF_D ), uPF_K(:,iPF_V1), uPF_K(:,iPF_V2), &
                   uPF_K(:,iPF_V3), uPF_K(:,iPF_E ), uPF_K(:,iPF_Ne), &
                   uCF_K(:,iCF_D ), uCF_K(:,iCF_S1), uCF_K(:,iCF_S2), &
                   uCF_K(:,iCF_S3), uCF_K(:,iCF_E ), uCF_K(:,iCF_Ne), &
                   uGF_K(:,iGF_Gm_dd_11), &
                   uGF_K(:,iGF_Gm_dd_22), &
                   uGF_K(:,iGF_Gm_dd_33), &
                   uAF_K(:,iAF_P) )

          uCF(iX1,iX2,iX3,lo_F(4):hi_F(4)) &
            = RESHAPE( uCF_K, [ hi_F(4) - lo_F(4) + 1 ] )

        END DO
        END DO
        END DO

      END DO

      CALL amrex_mfiter_destroy( MFI )

    END DO

    DO iDim = 1, 3

      CALL DestroyMesh( MeshX(iDim) )

    END DO

  END SUBROUTINE InitializeFields_Implosion_TABLE

END MODULE MF_InitializationModule_NonRelativistic_TABLE

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


text:
 The Basics
Foosball is just Table games table soccer. It it is made up of 4 poles per each player. One for the goalie (13 men), another for the defense (2 men), one more for the midfield (5 men), and the final one for the forwards (3 men). Spinning is usually not allowed for the more serious players.

 Foosball table locations

  The California Aggie
  AmyLee Amy Lees house
  Any Fraternities Fraternity worth its salt.
  Tercero: inside Trudys at the back (prior renovation; dont know if its still there (btw I was there Spring of 04))
  Sylvan Learning Center at the U Mall
  Davis Academy Driving School
  Organon Professional Services
  Tako Korean BBQ Tacos
  Please God tell me there are more.
  Sudwerk Brewery Dock Store


submodule (modsparse) modsparse_ll
 !$ use omp_lib
 implicit none

contains

!**CONSTRUCTOR
module subroutine constructor_sub_ll(sparse,m,n,lupper,unlog)
 class(llsparse),intent(out)::sparse
 integer(kind=int32),intent(in)::m
 integer(kind=int32),intent(in),optional::n,unlog
 logical,intent(in),optional::lupper

 call sparse%initialize('LINKED LIST',m,m)

 if(present(n))sparse%dim2=n
 if(present(lupper))sparse%lupperstorage=lupper
 if(present(unlog))sparse%unlog=unlog

 sparse%lsymmetric=.false.

 allocate(sparse%heads(sparse%dim1))

end subroutine

!**DESTROY
module subroutine destroy_scal_ptrnode(pnode)
 type(ptrnode)::pnode

 type(ptrnode)::cursor

 do while(associated(pnode%p))
  cursor=pnode
  pnode=pnode%p%next
  deallocate(cursor%p)
  nullify(cursor%p)
 enddo

end subroutine

module subroutine destroy_ll(sparse)
 class(llsparse),intent(inout)::sparse
 integer(kind=int32)::i

 call sparse%destroy_gen_gen()

 if(allocated(sparse%heads))then
  do i=1,size(sparse%heads)
   call destroy_scal_ptrnode(sparse%heads(i))
  enddo
  deallocate(sparse%heads)
 endif

end subroutine

!**DIAGONAL ELEMENTS

!EQUALITIES
module subroutine equal_node(nodeout,nodein)
 class(node),intent(out)::nodeout
 class(node),intent(in)::nodein

 nodeout%col=nodein%col
 nodeout%val=nodein%val

end subroutine

!**ADD ELEMENTS
module subroutine addtohead_ptrnode(pnode,col,val)
 type(ptrnode),intent(inout),pointer::pnode
 integer(kind=int32),intent(in)::col
 real(kind=wp),intent(in)::val

 type(ptrnode)::cursor

 allocate(cursor%p)
 cursor%p%next=pnode
 cursor%p%col=col
 cursor%p%val=val
 pnode=cursor

end subroutine

module subroutine addtohead_ll(sparse,row,col,val)
 class(llsparse),intent(inout),target::sparse
 integer(kind=int32),intent(in)::row,col
 real(kind=wp),intent(in)::val

 type(ptrnode)::cursor

 if(.not.validvalue_gen(sparse,row,col))return
 if(.not.validnonzero_gen(sparse,val))return
 if(sparse%lupperstorage.and..not.uppervalue_gen(row,col))return

 allocate(cursor%p)
 cursor%p%next=sparse%heads(row)
 cursor%p%col=col
 cursor%p%val=val
 sparse%heads(row)=cursor

end subroutine

module subroutine addinorder_ll(sparse,row,col,val)
 class(llsparse),intent(inout),target::sparse
 integer(kind=int32),intent(in)::row,col
 real(kind=wp),intent(in)::val

 type(ptrnode),pointer::cursor

 if(.not.validvalue_gen(sparse,row,col))return
 if(.not.validnonzero_gen(sparse,val))return
 if(sparse%lupperstorage.and..not.uppervalue_gen(row,col))return

 cursor=>sparse%heads(row)
 do while(associated(cursor%p))
  if(cursor%p%col.ge.col)then
   if(.not.col.ge.cursor%p%col)then
    call addtohead_ptrnode(cursor,col,val)
   else
    cursor%p%val=cursor%p%val+val
   endif
   return
  endif
  cursor=>cursor%p%next
 enddo
 allocate(cursor%p)
 cursor%p%col=col
 cursor%p%val=val

end subroutine

module subroutine addtotail_ll(sparse,row,col,val)
 class(llsparse),intent(inout),target::sparse
 integer(kind=int32),intent(in)::row,col
 real(kind=wp),intent(in)::val

 type(ptrnode),pointer::cursor

 if(.not.validvalue_gen(sparse,row,col))return
 if(.not.validnonzero_gen(sparse,val))return
 if(sparse%lupperstorage.and..not.uppervalue_gen(row,col))return

 cursor=>sparse%heads(row)
 do while(associated(cursor%p))
  cursor=>cursor%p%next
 enddo
 allocate(cursor%p)
 cursor%p%col=col
 cursor%p%val=val

end subroutine

!**GET ELEMENTS
module function get_ll(sparse,row,col) result(val)
 class(llsparse),intent(in)::sparse
 integer(kind=int32),intent(in)::row,col
 real(kind=wp)::val

 integer(kind=int32)::trow,tcol
 type(ptrnode),pointer::cursor
 type(ptrnode),target::replacecursor

 val=0.0_wp

 trow=row
 tcol=col
 if(sparse%lupperstorage.and.row.gt.col)then
  !swap row-col
  trow=col
  tcol=row
 endif

 !cursor=>sparse%heads(trow)
 replacecursor=sparse%heads(trow)
 cursor=>replacecursor

 do while(associated(cursor%p))
  if(cursor%p%col.eq.tcol)then
   val=cursor%p%val
   exit
  endif
  cursor=>cursor%p%next
 enddo

end function

!**EXTERNAL

!**LOAD

!**MULTIPLICATIONS
module subroutine multgenv_ll(sparse,alpha,trans,x,val,y)
 !Computes y=val*y+alpha*sparse(tranposition)*x
 class(llsparse),intent(in)::sparse
 real(kind=wp),intent(in)::val,alpha
 real(kind=wp),intent(in)::x(:)
 real(kind=wp),intent(out)::y(:)
 character(len=1),intent(in)::trans

 y=0._wp
 write(sparse%unlog,'(a)')' ERROR: Multiplication mult not implemented for llsparse'
 error stop

end subroutine

module subroutine multgenm_ll(sparse,alpha,trans,x,val,y)
 !Computes y=val*y+alpha*sparse(tranposition)*x
 class(llsparse),intent(in)::sparse
 real(kind=wp),intent(in)::val,alpha
 real(kind=wp),intent(in)::x(:,:)
 real(kind=wp),intent(out)::y(:,:)
 character(len=1),intent(in)::trans

 y=0._wp
 write(sparse%unlog,'(a)')' ERROR: Multiplication mult not implemented for llsparse'
 error stop

end subroutine

!**NUMBER OF ELEMENTS
module function totalnumberofelements_ptrnode(pnode) result(nel)
 class(ptrnode),intent(in),target::pnode
 integer(kind=int64)::nel

 type(ptrnode),pointer::cursor

 nel=0
 cursor=>pnode
 do while(associated(cursor%p))
  nel=nel+1
  cursor=>cursor%p%next
 enddo

end function

module function totalnumberofelements_ll(sparse) result(nel)
 class(llsparse),intent(in)::sparse
 integer(kind=int64)::nel

 integer(kind=int32)::i

 nel=0
 do i=1,sparse%dim1
  nel=nel+sparse%heads(i)%size()
 enddo

end function

!**SAVE

!**SCALE ALL ENTRIES
module subroutine scale_ll(sparse,val)
 class(llsparse),intent(inout)::sparse
 real(kind=wp),intent(in)::val

 write(sparse%unlog,'(a)')' ERROR: scale not implemented for llsparse'
 error stop

end subroutine

!**SET ELEMENTS

!**SOLVE

!**SORT ARRAY

!**SUBMATRIX

!**PRINT
module subroutine print_ll(sparse,lint,output)
 class(llsparse),intent(in)::sparse
 integer(kind=int32),intent(in),optional::output
 logical,intent(in),optional::lint

 integer(kind=int32)::i,un
 character(len=20)::frm='(2(i0,1x),g0)'
 logical::linternal
 type(ptrnode),pointer::cursor
 type(ptrnode),target::replacecursor

 linternal=.true.
 if(present(lint))linternal=lint

 un=sparse%unlog
 if(present(output))un=output

 do i=1,sparse%dim1
  !cursor=>sparse%heads(i)
  replacecursor=sparse%heads(i)
  cursor=>replacecursor
  do while(associated(cursor%p))
   write(un,frm)i,cursor%p%col,cursor%p%val
   if(.not.linternal.and.sparse%lupperstorage.and.sparse%lsymmetric.and.cursor%p%col.ne.i)then
    write(un,frm)cursor%p%col,i,cursor%p%val
   endif
   cursor=>cursor%p%next
  enddo
 enddo

end subroutine

module subroutine printsquare_ll(sparse,output)
 class(llsparse),intent(inout)::sparse
 integer(kind=int32),intent(in),optional::output

 integer(kind=int32)::un
 real(kind=wp),allocatable::tmp(:)

 un=sparse%unlog
 if(present(output))un=output

 write(un,'(a)')' Warning: this subroutine is not implemented yet!'

 allocate(tmp(sparse%dim2))

end subroutine

end submodule


subroutine map_task(ilevel)
USE Multigrid_Levels
USE Parallel_Variables
use kind_parameters
implicit none
include "mpif.h"

! Local variables
integer :: imax, jmax
integer, intent(in) :: ilevel
real(kind=dp):: res

!write(*,*) 'num_levels =', num_level

 imax = mg_level(ilevel)%nx
 jmax = mg_level(ilevel)%ny

block_size = (jmax-2)/numtasks
res = mod((jmax-2), numtasks)

if (taskid < res) then

 block_size = block_size + 1
 start_x = 2 + taskid*block_size

else

 start_x = 2 + res + taskid*block_size

end if

end_x = min(start_x + block_size - 1, jmax-1)

mg_level(ilevel)%xi = start_x
mg_level(ilevel)%xf = end_x

end subroutine


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE DLWQ11 ( IOUT   , CONC   , BOUND  , ITIME  , IDT    ,
     &                    IDSTRT , IDSTOP , IDSTEP , SNAME  , MNAME  ,
     &                    CGRID  , LGRID  , NX     , NY     , NOTOT  ,
     &                             NOSYS  , IP     , ISFLAG , IDFLAG )
!
!
!     Deltares        SECTOR WATERRESOURCES AND ENVIRONMENT
!
!     CREATED            : june 1988  BY M.E.Sileon
!
!     FUNCTION           : Writes concentrations in grid-layout and
!                          writes the result in IOUT .
!
!     LOGICAL UNITS      : IOUT = number of dump file
!
!     SUBROUTINES CALLED : none
!
!     PARAMETERS         :
!
!     NAME    KIND     LENGTH      FUNCT.  DESCRIPTION
!     ---------------------------------------------------------
!     IOUT    INTEGER  1           INPUT   unit number output file
!     CONC    REAL     NOTOT*?     INPUT   concentration values
!     BOUND   REAL     NOTOT*?     INPUT   boundary      values
!     ITIME   INTEGER  1           INPUT   present time in clock units
!     IDT     INTEGER  1           INPUT   time step of simulation
!     IDSTRT  INTEGER  1           INPUT   start time of dump output
!     IDSTOP  INTEGER  1           INPUT   stop time of dump output
!     IDSTEP  INTEGER  1           INPUT   time step of dump output
!     SNAME   CHAR*20  NOTOT       INPUT   names of substances
!     MNAME   CHAR*40  4           INPUT   model identhification
!     CGRID   CHAR*6   20*NY       IN/OUT  concentrations in grid layout
!     LGRID   INTEGER  NX*NY       INPUT   grid layout
!     NX      INTEGER  1           INPUT   number of columns in grid
!     NY      INTEGER  1           INPUT   number of rows in grid
!     NOTOT   INTEGER  1           INPUT   total number of systems
!     NOSYS   INTEGER  1           INPUT   number of active systems
!     IP      INTEGER  6           IN/OUT  paging structure
!     ISFLAG  INTEGER  1           INPUT   if 1 then dd-hh:mm'ss'
!     IDFLAG  INTEGER  1           OUTPUT  TRUE if dump took place
!
!
      use timers
      DIMENSION     CONC  ( NOTOT, * ) , BOUND( NOSYS, * ), IP (*)
      CHARACTER*6   CGRID ( 20, * )    , POINT , PADDER
      CHARACTER*20  SNAME ( * )
      CHARACTER*40  MNAME ( * )
      DIMENSION     LGRID ( * )
      LOGICAL       IDFLAG
      DATA          POINT / '  .   ' /
      integer(4) ithandl /0/
      if ( timon ) call timstrt ( "dlwq11", ithandl )
!
      IDFLAG = .FALSE.
      IF ( IDSTEP                     .LE. 0      ) goto 9999
      IF ( ITIME                      .LT. IDSTRT ) goto 9999
      IF ( ITIME-IDT                  .GE. IDSTOP ) goto 9999
      IF ( MOD( ITIME-IDSTRT, IDSTEP ).GE. IDT    ) goto 9999
      IDFLAG = .TRUE.
      IF ( NX*NY                      .EQ. 0      ) goto 9999
!
!         initialise the paging
!
      IF ( IP(3) .EQ. 0 ) THEN
           IP(3) = MAX(1,IP(1)/(7+(NY+5)*((NX+IP(2)-1)/IP(2))))
           IP(4) = 0
      ENDIF
!
!      repeat output for every substance
!
      DO 50 ITOT =  1 , NOTOT
!
!      Calculate maximum concentration of displayed segments
!
      CMAX = 0.0
      DO 10 I1 = 1, NY
      DO 10 I2 = 1, NX
      I3 = LGRID ( (I1-1)*NX+I2 )
      IF ( I3 .GT. 0 ) CMAX = AMAX1 ( CMAX, CONC (ITOT, I3) )
      IF ( I3 .LT. 0 .AND. ITOT .LE. NOSYS )
     *                 CMAX = AMAX1 ( CMAX, BOUND(ITOT,-I3) )
   10 CONTINUE
!
!      Calculate scale factor
!
      IF ( CMAX .LE. 0.0 ) THEN
           ISCALE = 0
      ELSE IF ( CMAX .GE. 0.5 ) THEN
              ISCALE =  AINT( ALOG10(CMAX) + 2.2E-5 )
           ELSE
              ISCALE = -AINT(-ALOG10(CMAX) - 2.2E-5 ) - 1
      ENDIF
!
!         start printing
!
      IF ( MOD(IP(4),IP(3)) .EQ. 0 ) THEN
           WRITE (IOUT,'(''1'')')
           WRITE (IOUT,2040 ) ( MNAME(K),K=1,4)
      ENDIF
      IP(4) = IP(4) + 1
      CALL REPTIM ( IOUT  , ITIME , ISFLAG, -999.)
      WRITE(IOUT,2000) SNAME (ITOT) , ISCALE
!
!      Put concentration values in grid layout
!
      FACTOR = 10.0**ISCALE
      DO 40 I = 1, NX, IP(2)
         DO 30 J = 1, NY
            NEND = MIN ( NX, I+IP(2)-1 )
            DO 20 K = I, NEND
            CGRID ( K-I+1,J ) = POINT
            I3 = LGRID ( (J-1)*NX+K )
            IF ( I3 .GT. 0 ) THEN
                 WRITE ( PADDER, '(F6.3)')   CONC (ITOT, I3 )/FACTOR
                 CGRID(K-I+1,J) = PADDER
            ENDIF
            IF ( I3 .LT. 0 .AND. ITOT .LE. NOSYS ) THEN
                 WRITE ( PADDER, '(F6.3)')   BOUND(ITOT,-I3 )/FACTOR
                 CGRID(K-I+1,J) = PADDER
            ENDIF
   20       CONTINUE
            WRITE ( IOUT, 2030 ) ( CGRID ( K-I+1, J ) , K=I, NEND )
   30   CONTINUE
        WRITE ( IOUT, '('' '')' )
   40 CONTINUE
!
   50 CONTINUE
!
 9999 if ( timon ) call timstop ( ithandl )
      RETURN
!
 2000 FORMAT (/' DUMP OF RESULTS OF ',A20,
     *          ' SCALE FACTOR = 10.0 EXP (',I3,' ).'//)
 2030 FORMAT (   6X , 20A6  )
 2040 FORMAT (  45X ,   A40 )
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Bad coding example 1
!
! Shamefully written by Ross Walker (SDSC, 2006)
!
! This code reads a series of coordinates and charges from the file
! specified as argument $1 on the command line.
!
! This file should have the format:
!  I9
! 4F10.4   (repeated I9 times representing x,y,z,q)
!
! It then calculates the following fictional function:
!
!             exp(rij*qi)*exp(rij*qj)   1
!    E = Sum( ----------------------- - - )  (rij <= cut)
!        j<i           r(ij)            a
!
! where cut is a cut off value specified on the command line ($2),
! r(ij) is a function of the coordinates read in for each atom and
! a is a constant.
!
! The code prints out the number of atoms, the cut off, total number of
! atom pairs which were less than or equal to the distance cutoff, the
! value of E, the time take to generate the coordinates and the time
! taken to perform the calculation of E.
!
! All calculations are done in double precision.
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

program calc_energy

     implicit none

     integer*4 :: iarg  !Intrinsic getarg require 4 byte integers so this
                  ! provides a guard against 64 bit integer compilation.
     integer :: natom, ier, ios, i, j, cut_count

     !Timer variables
     integer :: int_time, int_count_rate
     REAL*8 :: time0, time1, time2, count_rate

     character(len=80) :: arg, filename
     REAL*8 :: cut     !Cut off for Rij in distance units
     REAL*8, allocatable, dimension(:,:) :: coords
     REAL*8, allocatable, dimension(:) :: q
     REAL*8 :: total_e, current_e, vec2, rij
     REAL*8 :: a

     a = 3.2d0

     call system_clock(int_time, int_count_rate)  !Start time
     time0 = real(int_time)
     write(6,'(''Value of system clock at start = '',F14.4,/)') time0

     ! Step 1 - obtain the filename of the coord file and the value of
     ! cut from the command line.
     !         Argument 1 should be the filename of the coord file (char).
     !         Argument 2 should be the cut off (float).
     iarg = 1
     call getarg(iarg, arg)
     filename = arg
     write(6,'(''Coordinates will be read from file: '',A,/)') filename
     iarg=2
     call getarg(iarg, arg)
     read(arg,'(f10.5)') cut

     ! Step 2 - Open the coordinate file and read the first line to
     ! obtain the number of atoms
     open(unit=10,file=filename,form='FORMATTED',iostat=ios)
     read(10,'(i9)') natom

     write(6,*) 'Natom = ', natom
     write(6,*) '  cut = ', cut

     ! Step 3 - Allocate the arrays to store the coordinate and charge data
     ier = 0
     allocate ( coords(natom,3), stat=ier )
     if ( ier /= 0 ) then
       write(6,*) 'Allocation error coords'
       stop
     end if
     allocate ( q(natom), stat=ier )
     if ( ier /= 0 ) then
       write(6,*) 'Allocation error q'
       stop
     end if

     ! Step 4 - read the coordinates and charges.
     do i = 1, natom
        read(10,'(4f10.4)') coords(i,1:3),q(i)
     end do

     call system_clock(int_time, int_count_rate)  !time after reading of file
     time1 = real(int_time)
     write(6,'(''Value of system clock after coord read = '',F14.4,/)') time1

     ! Step 5 - calculate the number of pairs and E. - this is the
     ! majority of the work.
     total_e = 0.0d0
     cut_count = 0
     do i = 1, natom
       do j = 1, natom
         if ( j < i ) then   !Avoid double counting.
           vec2 = (coords(i,1)-coords(j,1))**2 + (coords(i,2)-coords(j,2))**2 &
                + (coords(i,3)-coords(j,3))**2
                  !X^2 + Y^2 + Z^2
           rij = sqrt(vec2)
           !Check if this is below the cut off
           if ( rij <= cut ) then
              !Increment the counter of pairs below cutoff
             cut_count = cut_count+1
             current_e = (exp(rij*q(i))*exp(rij*q(j)))/rij
             total_e = total_e + current_e - 1.0d0/a
           end if
         end if
       end do
     end do

     !time after reading of file and calculation
     call system_clock(int_time, int_count_rate)
     time2 = real(int_time)
     count_rate = real(int_count_rate)
     write(6,'(''Value of system clock after coord read and E calc = '',F14.4,/)') &
          time2

     ! Step 6 - write out the results
     write(6,*) '                         Final Results'
     write(6,*) '                         -------------'
     write(6,*) '                   Num Pairs = ',cut_count
     write(6,*) '                     Total E = ',total_e
     write(6,'(''     Time to read coord file = '',F14.4, '' Seconds'')') &
          (time1-time0)/count_rate
     write(6,'(''         Time to calculate E = '',F14.4, '' Seconds'')') &
          (time2-time1)/count_rate
     write(6,'(''        Total Execution Time = '',F14.4, '' Seconds'')') &
          (time2-time0)/count_rate

     ! Step 7 - Deallocate the arrays - we should strictly check the
     ! return values here but for the purposes of this tutorial we
     ! can ignore this.
     deallocate(q, stat=ier)
     deallocate(coords, stat=ier)

     close (unit=10)

     stop

end program calc_energy

max_stars_repo_name: stevenstetzler/parallel-computing

max_issues_repo_name: stevenstetzler/parallel-computing

max_forks_repo_name: stevenstetzler/parallel-computing


text:
!
!!  Copyright (C) 1996 Ernst A. Meese
!!    Refer to the file copyright.doc for details and important disclaimer.
!
!!  Created: January 1996
!!  Development version
!
MODULE USR_SMlib_Precision
  !
  !!  This module is part of SMLIB v. 1.1.  It defines the real precision USEd
  !!  throughout the library.  Redefine prec in this file to change working
  !!  precision.
  !
  !! ------------------------------------------------------------------------
  !
  !!  Copyright (C) 1996 Ernst A. Meese
  !!    Refer to the file copyright.doc for details and important disclaimer.
  !
  !!  Created: January 1996
  !!  Development version
  !
  !! ------------------------------------------------------------------------
  !!  Module contents
  !
  !!   Parameters:
  !!      single_prec - Kind value for single precision.
  !!      double_prec - Kind value for double precision.
  !!      prec        - Working precision of SMLIB.
  !
  !! -----------------------------------------------------------------------
  !
  IMPLICIT NONE

  INTEGER, PARAMETER :: single_prec = KIND(1.0E0)
  INTEGER, PARAMETER :: double_prec = KIND(1.0D0)
  INTEGER, PARAMETER ::        prec = double_prec

END MODULE USR_SMlib_Precision

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
C
C***********************************************************************
C                                                                      *
C            NACHOS II - A Finite Element Computer Program             *
C                        for Incompressible Flow Problems              *
C                                                                      *
C     Copyright (c) 1986,2019   National Technology & Engineering      *
C                               Solutions of Sandia, LLC (NTESS)       *
C                                                                      *
C                            All rights reserved.                      *
C                                                                      *
C     This software is distributed under the BSD 3-Clause License.     *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE FLUX8 (NSTRSQ,SQ,TQ)
C
C     ******************************************************************
C
C     SUBROUTINE TO CALCULATE THE VALUES OF THE EIGHT NODE 'SERENDIPITY'
C     SHAPE FUNCTIONS AND THEIR DERIVATIVES AT THE FLUX EVALUATION
C     POINTS.
C     THE POINTS MAY BE ON THE ELEMENT BOUNDARY OR AT INTERIOR POINTS
C
C     ******************************************************************
C
      COMMON /QUAD8/ F8Q(8,9),F8L(4,9),DF8QDS(8,9),DF8QDT(8,9)
C
      DIMENSION SQ(*), TQ(*)
      DIMENSION SI(8), TI(8)
C
      DATA (SI(I),I=1,8)/-1.,1.,1.,-1.,0.,1.,0.,-1./
      DATA (TI(I),I=1,8)/-1.,-1.,1.,1.,-1.,0.,1.,0./
C
C     STATEMENT FUNCTIONS FOR THE SHAPE FUNCTIONS AND THEIR DERIVATIVES
C
      FC(I,S,T)=.25*(1.+S*SI(I))*(1.+T*TI(I))*(S*SI(I)+T*TI(I)-1.)
      FM(I,S,T)=.50*((1.-S**2)*(1.+T*TI(I))*TI(I)**2+(1.-T**2)*(1.+S*SI
     1(I))*SI(I)**2)
      G(I,S,T)=.25*(1.+S*SI(I))*(1.+T*TI(I))
      DFCDS(I,S,T)=.25*(1.+T*TI(I))*(2.*S*SI(I)+T*TI(I))*SI(I)
      DFMDS(I,S,T)=.50*((1.+T*TI(I))*(-2.*S)*TI(I)**2+(1.-T**2)*SI(I))
      DFCDT(I,S,T)=.25*(1.+S*SI(I))*(2.*T*TI(I)+S*SI(I))*TI(I)
      DFMDT(I,S,T)=.50*((1.+S*SI(I))*(-2.*T)*SI(I)**2+(1.-S**2)*TI(I))
C
C     ******************************************************************
C
C     EVALUATE CORNER FUNCTIONS
C
      DO 10 I=1,4
C
      DO 10 J=1,NSTRSQ
      DF8QDS(I,J)=DFCDS(I,SQ(J),TQ(J))
      DF8QDT(I,J)=DFCDT(I,SQ(J),TQ(J))
      F8Q(I,J)=FC(I,SQ(J),TQ(J))
      F8L(I,J)=G(I,SQ(J),TQ(J))
   10 CONTINUE
C
C     EVALUATE MID-SIDE FUNCTIONS
C
      DO 20 I=5,8
C
      DO 20 J=1,NSTRSQ
      DF8QDS(I,J)=DFMDS(I,SQ(J),TQ(J))
      DF8QDT(I,J)=DFMDT(I,SQ(J),TQ(J))
      F8Q(I,J)=FM(I,SQ(J),TQ(J))
   20 CONTINUE
C
      RETURN
      END

max_stars_repo_name: nachos-ii/nachos-ii

max_issues_repo_name: nachos-ii/nachos-ii

max_forks_repo_name: nachos-ii/nachos-ii


text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                       .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c ... visgeo.f
c
c     contains documentation and code for subroutine visgeo
c
      SUBROUTINE VISGEO (M,IDP,JDP,X,Y,Z,H,EYER,EYELAT,EYELON,
     1                          WORK,LWORK,IWORK,LIWORK,IERROR)
c
c     subroutine visgeo will display a function on the sphere
c     as a solid. ie. as a "lumpy" sphere. visgeo calls subroutine
c     vsurf to produce the visible surface rendering. X, Y, and Z
c     are the points on an icosahedral geodesic computed by
c     subroutine geopts available in spherepack.
c
c     requires routines visgeo1 ctos stoc vsurf vsurf1
c                       prjct box
c
c     visgeo uses the ncar graphics package.
c     compile with: ncargf77 (all programs above)
c
c     execute with:  a.out
c
c     on screen display with:  ctrans -d x11 gmeta
c
c     print with:  ctrans -d ps.color gmeta > gmeta.ps
c                  lpr -P(your printer) gmeta.ps
c
c
c     input parameters
c
c     m        the number of points on one edge of the icosahedron
c
c     idp,jdp  the first and second dimensions of the three
c              dimensional arrays x, y, z, and h.
c
c     x,y,z    the coordinates of the geodesic points on
c              the unit sphere computed by subroutine geopts.
c              the indices are defined on the unfolded
c              icosahedron as follows for the case m=3
c
c                north pole
c
c                 (5,1)          0      l
c        i     (4,1) (5,2)              a    (repeated for
c           (3,1) (4,2) (5,3)  theta1   t    k=2,3,4,5 in
c        (2,1) (3,2) (4,3)              i        -->
c     (1,1) (2,2) (3,3)        theta2   t    the longitudinal
c        (1,2) (2,3)                    u    direction)
c           (1,3)                pi     d
c      j                                e
c         south pole
c
c            total number of vertices is  10*(m-1)**2+2
c            total number of triangles is 20*(m-1)**2
c
c     h      a three dimensional array that contains the discrete
c            function to be displayed. h(i,j,k) is the distance from
c            the center of the sphere to the "lumpy" surface at the
c             point [x(i,j,k),y(i,j,k),z(i,j,k)] on the unit sphere.
c
c     eyer   the distance from the center of the sphere to the eye.
c
c     eyelat the colatitudinal coordinate of the eye (in degrees).
c
c     eyelon the longitudinal  coordinate of the eye (in degrees).
c
c     idp    the first dimension of the array h as it appears in
c            the program that calls visgeo
c
c     jdp    the second dimension of the array h as it appears in
c            the program that calls visgeo
c
c     work   a real work array
c
c     lwork  the dimension of the array work as it appears in the
c            program that calls visgeo. lwork must be at least
c                       480*(m-1)**2.
c
c     iwork  an integer work array
c
c     liwork the dimension of the array iwork as it appears in the
c            program that calls visgeo. liwork must be at least
c                       140*(m-1)**2.
c
c     input parameter
c
c     ierror = 0    no error
c            = 1    h(i,j,k) is less than zero for some i,j,k.
c            = 2    eyer is less than h(i,j,k) for some i,k,k.
c            = 3    lwork  is less than 480*(m-1)**2
c            = 4    liwork is less than 140*(m-1)**2
c
      dimension h(idp,jdp,5),x(idp,jdp,5),y(idp,jdp,5),z(idp,jdp,5),
     1                                                        work(*)
      INTEGER IWORK(*)
      mmsq = (m-1)**2
      ierror = 3
      if(lwork .lt. 480*mmsq) return
      ierror = 4
      if(liwork .lt. 140*mmsq) return
      do 10 k=1,5
      do 10 j=1,m
      do 10 i=1,m+m-1
      if(h(i,j,k) .ge. 0.) go to 15
      ierror = 1
      return
   15 if(eyer .gt. h(i,j,k)) go to 10
      ierror = 2
      return
   10 continue
      ierror = 0
      lt = 20*(m-1)**2
      lg = 5*m*(m+m-1)
      i1 = 1
      i2 = i1+lt
      i3 = i2+lt
      i4 = i3+lt
      i5 = i4+lt
      i6 = i5+lt
      i7 = i6+lt
      i8 = i7+lt
      i9 = i8+lt
      i10 = i9+lt
      i11 = i10+lt
      i12 = i11
      i13 = i12+lg
      i14 = i13+lg
      call visgeo1 (m,idp,jdp,h,eyer,eyelat,eyelon,x,y,z,
     1work(i1),work(i2),work(i3),work(i4),work(i5),work(i6),
     2work(i7),work(i8),work(i9),IWORK(1),work(i11),
     3work(i12),work(i13),work(i14),IWORK(lt+1))
      return
      end
      SUBROUTINE VISGEO1(M,IDP,JDP,H,EYER,EYELAT,EYELON,
     1    XI,YI,ZI,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,ITYPE,WORK,X,Y,Z,IWORK)
      dimension h(idp,jdp,5),xi(idp,jdp,5),yi(idp,jdp,5),zi(idp,jdp,5),
     1x1(*),y1(*),z1(*),x2(*),y2(*),z2(*),x3(*),y3(*),z3(*),itype(*),
     2work(*),x(m+m-1,m,5),y(m+m-1,m,5),z(m+m-1,m,5)
      INTEGER IWORK(*)
c
c     the * above refers to 20*(m-1)**2 locations which is the
c     number of triangles
c
      do 10 k=1,5
      do 10 j=1,m
      do 10 i=1,m+m-1
      call ctos(xi(i,j,k),yi(i,j,k),zi(i,j,k),rad,theta,elambda)
      call stoc(h(i,j,k),theta,elambda,x(i,j,k),y(i,j,k),z(i,j,k))
   10 continue
      ntri = 0
      do 20 k=1,5
      do 20 j=1,m-1
      do 20 i=1,m+m-2
      ntri = ntri+1
      x1(ntri) = x(i,j,k)
      y1(ntri) = y(i,j,k)
      z1(ntri) = z(i,j,k)
      x2(ntri) = x(i+1,j+1,k)
      y2(ntri) = y(i+1,j+1,k)
      z2(ntri) = z(i+1,j+1,k)
      x3(ntri) = x(i+1,j,k)
      y3(ntri) = y(i+1,j,k)
      z3(ntri) = z(i+1,j,k)
      itype(ntri) = 13
      ntri = ntri+1
      x1(ntri) = x(i,j,k)
      y1(ntri) = y(i,j,k)
      z1(ntri) = z(i,j,k)
      x2(ntri) = x(i+1,j+1,k)
      y2(ntri) = y(i+1,j+1,k)
      z2(ntri) = z(i+1,j+1,k)
      x3(ntri) = x(i,j+1,k)
      y3(ntri) = y(i,j+1,k)
      z3(ntri) = z(i,j+1,k)
      itype(ntri) = 3
 20   continue
c      write(6,22) ntri
 22   format(i10)
c      write(6,23) (x1(l2),y1(l2),z1(l2),x2(l2),y2(l2),z2(l2),
c     1             x3(l2),y3(l2),z3(l2),l2=1,ntri)
c 23   format(9f10.7)
c
      pi = 4.*atan(1.)
      dtr = pi/180.
      xeye=eyer*sin(dtr*eyelat)
      yeye=xeye*sin(dtr*eyelon)
      xeye=xeye*cos(dtr*eyelon)
      zeye=eyer*cos(dtr*eyelat)
      CALL VSURF(XEYE,YEYE,ZEYE,NTRI,X1,Y1,Z1,X2,Y2,Z2,
     1                 X3,Y3,Z3,ITYPE,WORK,IWORK)
      return
      end
      subroutine ctos(x,y,z,r,theta,phi)
      r1 = x*x+y*y
      if(r1 .ne. 0.) go to 10
      phi = 0.
      theta = 0.
      if(z .lt. 0.) theta = 4.*atan(1.)
      return
   10 r = sqrt(r1+z*z)
      r1 = sqrt(r1)
      phi = atan2(y,x)
      theta = atan2(r1,z)
      return
      end
      subroutine stoc(r,theta,phi,x,y,z)
      st = sin(theta)
      x = r*st*cos(phi)
      y = r*st*sin(phi)
      z = r*cos(theta)
      return
      end
      SUBROUTINE VSURF(XEYE,YEYE,ZEYE,NTRI,X1,Y1,Z1,X2,Y2,Z2,
     1                 X3,Y3,Z3,ITYPE,WORK,IWORK)
c
c    subroutine vsurf is like subroutine hidel except the triangles
c    are categorized. vsurf is also like solid except triangles rather
c    than lines are covered.
c
c     written by paul n. swarztrauber, national center for atmospheric
c     research, p.o. box 3000, boulder, colorado, 80307
c
c    this program plots visible lines for the surface defined
c    by the input 3-d triangles with corners at (x1,y1,z1), (x2,y2,z2)
c    and (x3,y3,z3). the sides of these these triangles may or
c    may not be plotted depending on itype. if itype is 1 then the
c    side between points (x1,y1,z1) and (x2,y2,z2) is plotted if it
c    is visible. if itype is 2 then the side between (x2,y2,z2)
c    and (x3,y3,z3) is plotted. if itype is 3 then the visible portion
c    of the side between (x3,y3,z3) and (x1,y1,z1) is plotted.
c    any combination is possible by specifying itype to be one
c    of the following values: 0,1,2,3,12,13,23,123.
c
c    the length of real    array  work must be at least 14*ntri
c
c    the length of integer array iwork must be at least  6*ntri
c
c
c    the vertices of the triangles are renumbered by vsurf so that
c    their projections are orientated counterclockwise. the user need
c    only be aware that the vertices may be renumbered by vsurf.
c
      dimension x1(ntri),y1(ntri),z1(ntri),x2(ntri),y2(ntri),z2(ntri),
     1          x3(ntri),y3(ntri),z3(ntri),itype(ntri),work(14*ntri)
      INTEGER IWORK(6*NTRI)
c
      call vsurf1(xeye,yeye,zeye,ntri,x1,y1,z1,x2,y2,z2,x3,y3,z3,
     1 itype,work,work(ntri+1),work(2*ntri+1),work(3*ntri+1),
     2 work(4*ntri+1),work(5*ntri+1),work(6*ntri+1),work(7*ntri+1),
     3 work(8*ntri+1),work(9*ntri+1),work(10*ntri+1),work(11*ntri+1),
     3 work(12*ntri+1),work(13*ntri+1),IWORK,IWORK(NTRI+1),
     4 IWORK(2*NTRI+1),IWORK(4*NTRI+1))
      return
      end
      subroutine vsurf1(xeye,yeye,zeye,ntri,x1,y1,z1,x2,y2,z2,x3,y3,z3,
     1 itype,px1,py1,px2,py2,px3,py3,vx1,vy1,vx2,vy2,vx3,vy3,tl,tr,kh,
     2 next,istart,ifinal)
c
      dimension x1(ntri),y1(ntri),z1(ntri),x2(ntri),y2(ntri),z2(ntri),
     1          x3(ntri),y3(ntri),z3(ntri),itype(ntri),
     2          px1(ntri),py1(ntri),px2(ntri),py2(ntri),
     3          px3(ntri),py3(ntri),vx1(ntri),vy1(ntri),
     4          vx2(ntri),vy2(ntri),vx3(ntri),vy3(ntri),
     5          tl(ntri),tr(ntri),next(ntri),kh(ntri),
     6          istart(2*ntri),ifinal(2*ntri),ltp(3),
     7          ird(11),ip2(11),nct(11),ncv(11),last(11)
c
      real l2e
      double precision le2
c
c     compute projections of 3-d points
c
      le2 = .6931471805599453094172321d0
      l2e = 1.d0/le2
      fntri = ntri
      irmax = .5*l2e*log(fntri)
      irmax = min(irmax,10)
      irmp1 = irmax+1
      do 4 icv=1,11
      ncv(icv) = 0
    4 continue
      nct(1) = 0
      ip2(1) = 1
      ird(1) = 0
      isize = 4
      do 7 irp1=2,irmp1
      ir = irp1-1
      nct(irp1) = 0
      ip2(irp1) = 2**ir
      ird(irp1) = ird(ir)+isize
      isize = (ip2(irp1)+1)**2
    7 continue
      isxm = ird(irmp1)+isize+1
      do 8 isx=1,isxm
      istart(isx) = 0
      ifinal(isx) = 0
    8 continue
      do 6 i=1,ntri
      next(i) = 0
    6 continue
      call prjct(0,xeye,yeye,zeye,x,y,z,dum1,dum2)
c      write(6,127) ntri
  127 format(' ntri in hidel', i5)
      do 86 k=1,ntri
      call prjct(1,xeye,yeye,zeye,x1(k),y1(k),z1(k),px1(k),py1(k))
      call prjct(1,xeye,yeye,zeye,x2(k),y2(k),z2(k),px2(k),py2(k))
      call prjct(1,xeye,yeye,zeye,x3(k),y3(k),z3(k),px3(k),py3(k))
        if (k .lt. 3) then
c          write(6,333) xeye,yeye,zeye,x1(k),y1(k),z1(k),px1(k),py1(k)
  333     format(' xeye, etc.',8e8.1)
        endif
   86 continue
c
c     orientate triangles counter clockwise
c
      do 70 k=1,ntri
      cprod = (px2(k)-px1(k))*(py3(k)-py1(k))-(py2(k)-py1(k))
     1       *(px3(k)-px1(k))
c      if(cprod.eq.0.) write(6,79) k,px1(k),px2(k),px3(k),
c     -                              py1(k),py2(k),py3(k)
   79 format('  cprod=0 at k=', i5,6e9.2)
      if(cprod.ge.0.) go to 70
      px1h = px1(k)
      py1h = py1(k)
      px1(k) = px2(k)
      py1(k) = py2(k)
      px2(k) = px1h
      py2(k) = py1h
      x1hold = x1(k)
      y1hold = y1(k)
      z1hold = z1(k)
      x1(k) = x2(k)
      y1(k) = y2(k)
      z1(k) = z2(k)
      x2(k) = x1hold
      y2(k) = y1hold
      z2(k) = z1hold
      ityp = itype(k)
      if(ityp.eq.2) itype(k) = 3
      if(ityp.eq.3) itype(k) = 2
      if(ityp.eq.12) itype(k) = 13
      if(ityp.eq.13) itype(k) = 12
   70 continue
c
c     set screen limits
c
      pmax = px1(1)
      pmin = px1(1)
      do 87 k=1,ntri
      pmin = amin1(pmin,px1(k),py1(k),px2(k),py2(k),px3(k),py3(k))
      pmax = amax1(pmax,px1(k),py1(k),px2(k),py2(k),px3(k),py3(k))
   87 continue
      pmin = 1.1*pmin
      pmax = 1.1*pmax
      call set(0.,1.,0.,1.,pmin,pmax,pmin,pmax,1)
      xmin = amin1(px1(1),px2(1),px3(1))
      xmax = amax1(px1(1),px2(1),px3(1))
      ymin = amin1(py1(1),py2(1),py3(1))
      ymax = amax1(py1(1),py2(1),py3(1))
      do 1 i=2,ntri
      xmin = amin1(xmin,px1(i),px2(i),px3(i))
      xmax = amax1(xmax,px1(i),px2(i),px3(i))
      ymin = amin1(ymin,py1(i),py2(i),py3(i))
      ymax = amax1(ymax,py1(i),py2(i),py3(i))
    1 continue
      dmx = xmax-xmin
      dmy = ymax-ymin
      if(dmx .gt. dmy) go to 2
      c = ymin
      d = ymax
      xmid = .5*(xmin+xmax)
      hdy = .5*dmy
      a = xmid-hdy
      b = xmid+hdy
      go to 3
    2 a = xmin
      b = xmax
      ymid = .5*(ymin+ymax)
      hdx = .5*dmx
      c = ymid-hdx
      d = ymid+hdx
    3 hgr = b-a
c
c     categorize triangles
c
      do 100 i=1,ntri
      xmin = amin1(px1(i),px2(i),px3(i))
      xmax = amax1(px1(i),px2(i),px3(i))
      ymin = amin1(py1(i),py2(i),py3(i))
      ymax = amax1(py1(i),py2(i),py3(i))
      dxt = amax1(xmax-xmin,ymax-ymin)
      if(dxt .gt. 0.) go to 10
      ir = irmax
      go to 20
   10 ir = l2e*log(hgr/dxt)
      ir = min(ir,irmax)
   20 irp1 = ir+1
      nct(irp1) = nct(irp1)+1
      hr = hgr/ip2(irp1)
      xmid = .5*(xmin+xmax)
      id = (xmid-a)/hr+1.5
      ymid = .5*(ymin+ymax)
      jd = (ymid-c)/hr+1.5
      ijd = ip2(irp1)+1
      isx = id+(jd-1)*ijd+ird(irp1)
      ifx = ifinal(isx)
      if(ifx .gt. 0) go to 50
      istart(isx) = i
      go to 60
   50 next(ifx) = i
   60 ifinal(isx) = i
  100 continue
c      write(6,106) tcat,(irp1,nct(irp1),irp1=1,irmp1)
  106 format(' time to categorize   ', e15.6/(' ir+1',i3,' ntri',i7))
c
c     sort triangles into boxes
c
      l = 0
      do 30 irp1=1,irmp1
      if(nct(irp1) .eq. 0) go to 30
      ist = ird(irp1)+1
      isd = ip2(irp1)+1
      call box(isd,istart(ist),next,l,ifinal)
      last(irp1) = l+1
   30 continue
      do 35 irp1=1,irmp1
      il = ird(irp1)+(ip2(irp1)+1)**2+1
      if(istart(il) .eq. 0) istart(il) = last(irp1)
   35 continue
c      write(6,31) tsort,l,ntri
   31 format(' time to sort  ', e15.6,'   l', i8,'   ntri',i8)
      do 90 k=1,ntri
      vx1(k) = px2(k)-px1(k)
      vy1(k) = py2(k)-py1(k)
      vx2(k) = px3(k)-px2(k)
      vy2(k) = py3(k)-py2(k)
      vx3(k) = px1(k)-px3(k)
      vy3(k) = py1(k)-py3(k)
   90 continue
      tl1 = 0.
      tl2 = 0.
      maxs = 0
      do 500 ir2=1,irmp1
      if(nct(ir2) .eq. 0) go to 500
      ist = ird(ir2)
      isd = ip2(ir2)+1
      do 490 j2=1,isd
      do 480 i2=1,isd
      ist = ist+1
      ls = istart(ist)
      lf = istart(ist+1)-1
      if(lf .lt. ls) go to 480
c
c     define coverings
c
      kcv = 0
      i2m = i2-1
      j2m = j2-1
      do 300 ir1=1,irmp1
      if(nct(ir1) .eq. 0) go to 300
      if(ir1 .ge. ir2) go to 260
      irdp = 2**(ir2-ir1)
      i1s = (i2m-1)/irdp
      i1f = (i2m+1)/irdp
      if = i2m+1-i1f*irdp
      if(if .gt. 0) i1f = i1f+1
      j1s = (j2m-1)/irdp
      j1f = (j2m+1)/irdp
      jf = j2m+1-j1f*irdp
      if(jf .gt. 0) j1f = j1f+1
      go to 270
  260 irdp = 2**(ir1-ir2)
      i1s = irdp*(i2m-1)
      i1f = irdp*(i2m+1)
      j1s = irdp*(j2m-1)
      j1f = irdp*(j2m+1)
  270 ijd = ip2(ir1)+1
      i1s = max(i1s+1,1)
      i1f = min(i1f+1,ijd)
      j1s = max(j1s+1,1)
      j1f = min(j1f+1,ijd)
      ixh = (j1s-2)*ijd+ird(ir1)
      ixs = i1s+ixh
      ixf = i1f+ixh
      do 290 j1=j1s,j1f
      ixs = ixs+ijd
      kds = istart(ixs)
      ixf = ixf+ijd
      kdf = istart(ixf+1)-1
      if(kdf .lt. kds) go to 290
      do 280 kd=kds,kdf
      kcv = kcv+1
      kh(kcv) = ifinal(kd)
  280 continue
  290 continue
  300 continue
      do 310 icv=1,10
      if(kcv .le. ncv(icv)) go to 310
      ncv(icv) = kcv
      go to 320
  310 continue
c
c
  320 do 470 ldo=ls,lf
      l = ifinal(ldo)
      ith = itype(l)
      if(ith .eq. 0) go to 470
      ltp(1) = 0
      ltp(2) = 0
      ltp(3) = 0
      id1 = ith/100
      ith = ith-100*id1
      id2 = ith/10
      id3 = ith-10*id2
      if(id1 .ne. 0) ltp(id1) = 1
      if(id2 .ne. 0) ltp(id2) = 1
      if(id3 .ne. 0) ltp(id3) = 1
c     if((ith.eq.123) .or. (ith.eq.12) .or.(ith.eq.13)) ltp(1) = 1
c     if((ith.eq.123) .or. (ith.eq.23) .or.(ith.eq.12)) ltp(2) = 1
c     if((ith.eq.123) .or. (ith.eq.13) .or.(ith.eq.23)) ltp(3) = 1
      do 460 ns=1,3
      go to (101,102,103),ns
  101 if(ltp(ns) .eq. 0) go to 460
      px4 = px1(l)
      py4 = py1(l)
      px5 = px2(l)
      py5 = py2(l)
      x4 = x1(l)
      y4 = y1(l)
      z4 = z1(l)
      x5 = x2(l)
      y5 = y2(l)
      z5 = z2(l)
      go to 105
  102 if(ltp(ns) .eq. 0) go to 460
      px4 = px2(l)
      py4 = py2(l)
      px5 = px3(l)
      py5 = py3(l)
      x4 = x2(l)
      y4 = y2(l)
      z4 = z2(l)
      x5 = x3(l)
      y5 = y3(l)
      z5 = z3(l)
      go to 105
  103 if(ltp(ns) .eq. 0) go to 460
      px4 = px1(l)
      py4 = py1(l)
      px5 = px3(l)
      py5 = py3(l)
      x4 = x1(l)
      y4 = y1(l)
      z4 = z1(l)
      x5 = x3(l)
      y5 = y3(l)
      z5 = z3(l)
  105 x54 = px5-px4
      y54 = py5-py4
      nseg = 0
      do 440 kd=1,kcv
      k = kh(kd)
      c17 = vx1(k)*y54-vy1(k)*x54
      c27 = vx2(k)*y54-vy2(k)*x54
      c37 = vx3(k)*y54-vy3(k)*x54
      c14 = vy1(k)*(px4-px1(k))-vx1(k)*(py4-py1(k))
      c25 = vy2(k)*(px4-px2(k))-vx2(k)*(py4-py2(k))
      c36 = vy3(k)*(px4-px3(k))-vx3(k)*(py4-py3(k))
      tmin = 0.
      tmax = 1.
      if(c17) 151,152,153
  151 tmax = amin1(c14/c17,tmax)
      go to 154
  152 if(c14) 154,440,440
  153 tmin = amax1(c14/c17,tmin)
  154 if(c27) 155,156,157
  155 tmax = amin1(c25/c27,tmax)
      go to 158
  156 if(c25) 158,440,440
  157 tmin = amax1(c25/c27,tmin)
  158 if(c37) 159,160,161
  159 tmax = amin1(c36/c37,tmax)
      go to 162
  160 if(c36) 162,440,440
  161 tmin = amax1(c36/c37,tmin)
  162 if(tmax-tmin .lt. .00001) go to 440
      xpl = x4+tmin*(x5-x4)
      ypl = y4+tmin*(y5-y4)
      zpl = z4+tmin*(z5-z4)
      xpr = x4+tmax*(x5-x4)
      ypr = y4+tmax*(y5-y4)
      zpr = z4+tmax*(z5-z4)
c
c     the projections of line and plane intersect
c     now determine if plane covers line
c
      vx1t = x2(k)-x1(k)
      vy1t = y2(k)-y1(k)
      vz1t = z2(k)-z1(k)
      vx2t = x3(k)-x1(k)
      vy2t = y3(k)-y1(k)
      vz2t = z3(k)-z1(k)
      apl = vy1t*vz2t-vy2t*vz1t
      bpl = vx2t*vz1t-vx1t*vz2t
      cpl = vx1t*vy2t-vx2t*vy1t
      dpl = apl*x1(k)+bpl*y1(k)+cpl*z1(k)
      vx3t = xpl-xeye
      vy3t = ypl-yeye
      vz3t = zpl-zeye
      den = apl*vx3t+bpl*vy3t+cpl*vz3t
      til = 0.
      if(den .eq. 0.) go to 410
      til = (dpl-apl*xeye-bpl*yeye-cpl*zeye)/den
  410 vx3t = xpr-xeye
      vy3t = ypr-yeye
      vz3t = zpr-zeye
      den = apl*vx3t+bpl*vy3t+cpl*vz3t
      tir = 0.
      if(den .eq. 0.) go to 412
      tir = (dpl-apl*xeye-bpl*yeye-cpl*zeye)/den
  412 if(til.ge..99999.and.tir.ge..99999) go to 440
      if(til.lt.1..and.tir.lt.1.) go to 164
      vx3t = xpr-xpl
      vy3t = ypr-ypl
      vz3t = zpr-zpl
      den = apl*vx3t+bpl*vy3t+cpl*vz3t
      tim = 0.
      if(den .eq. 0.) go to 414
      tim = (dpl-apl*xpl-bpl*ypl-cpl*zpl)/den
  414 thold = tmin+tim*(tmax-tmin)
      if(til.ge.1.) go to 163
      tmax = thold
      go to 164
  163 tmin = thold
  164 nseg = nseg+1
      tl(nseg) = tmin
      tr(nseg) = tmax
  440 continue
      maxs = max0(maxs,nseg)
      if(nseg-1)171,180,172
  171 call line(px4,py4,px5,py5)
      go to 460
c
c     order the segments according to left end point tl(k)
c
  172 do 173 k=2,nseg
      do 173 i=k,nseg
      if(tl(k-1).le.tl(i)) go to 173
      tlh = tl(k-1)
      trh = tr(k-1)
      tl(k-1) = tl(i)
      tr(k-1) = tr(i)
      tl(i) = tlh
      tr(i) = trh
  173 continue
c
c     eliminate segment overlap
c
      k1 = 1
      k2 = 1
  174 k2 = k2+1
      if(k2.gt.nseg) go to 176
      if(tr(k1).lt.tl(k2)) go to 175
      tr(k1) = amax1(tr(k1),tr(k2))
      go to 174
  175 k1 = k1+1
      tl(k1) = tl(k2)
      tr(k1) = tr(k2)
      go to 174
  176 nseg = k1
c
c     plot all segments of the line
c
  180 do 181 ks =1,nseg
      kb = nseg-ks+1
      tl(kb+1) = tr(kb)
      tr(kb) = tl(kb)
  181 continue
      tl(1) = 0.
      tr(nseg+1) = 1.
      nsegp = nseg+1
      do 450 k=1,nsegp
      if(abs(tr(k)-tl(k)).lt..000001) go to 450
      xa = px4+tl(k)*(px5-px4)
      ya = py4+tl(k)*(py5-py4)
      xb = px4+tr(k)*(px5-px4)
      yb = py4+tr(k)*(py5-py4)
      call line(xa,ya,xb,yb)
  450 continue
  460 continue
  470 continue
  480 continue
  490 continue
  500 continue
c      write(6,903) tl1,tl2
  903 format(' time to cover', e15.6/
     1       ' time to test ', e15.6)
c      write(6,904) maxs
  904 format(' maximum number of segments', i5)
c      write(6,250) (ncv(icv),icv=1,10)
  250 format('  the ten largest coverings'/(10i5))
      call frame
      end
      subroutine prjct(init,xeye,yeye,zeye,x,y,z,px,py)
c
c     subroutine prjct projects the point x,y,z onto a plane through
c     the origin that is perpendicular to a line between the origin
c     and the eye. the projection is along the line between the eye
c     and the point x,y,z. px and py are the coordinates of the
c     projection in the plane.
c     (version 2 , 12-10-82)
c
      save
      if(init.ne.0) go to 1
      rads1 = xeye**2+yeye**2
      rads2 = rads1+zeye**2
      d1 = sqrt(rads1)
      d2 = sqrt(rads2)
      if(d1.ne.0.) go to 2
      cx1 = 1.
      cy1 = 0.
      cx2 = 0.
      cy2 = 1.
      cz2 = 0.
      cx3 = 0.
      cy3 = 0.
      cz3 = 1.
      return
    2 cx1 = -yeye/d1
      cy1 = xeye/d1
      cx2 = -xeye*zeye/(d1*d2)
      cy2 = -yeye*zeye/(d1*d2)
      cz2 = d1/d2
      cx3 = xeye/d2
      cy3 = yeye/d2
      cz3 = zeye/d2
      return
    1 x1 = cx1*x+cy1*y
      y1 = cx2*x+cy2*y+cz2*z
      z1 = cx3*x+cy3*y+cz3*z
      ratio = d2/(d2-z1)
      px = ratio*x1
      py = ratio*y1
      return
      end
      subroutine box(isd,istart,next,l,list)
      dimension istart(isd,isd),next(1),list(1)
      do 30 jd=1,isd
      do 10 id=1,isd
      idx = istart(id,jd)
      istart(id,jd) = l+1
      if(idx .eq. 0) go to 10
   20 l = l+1
      list(l) = idx
      if(next(idx) .eq. 0) go to 10
      idx = next(idx)
      go to 20
   10 continue
   30 continue
      return
      end


text: Im ketan. Im Economics economics and Statistics statistics at davis 10. I make music. http://www.soundcloud.com/caden. lolz.