text:
include '../functionsS.f90'
!
! Diagonal tension criteria for curved stiffened panels
! Written by Saullo Castro, July 2015
!
subroutine R3SVALS(GRPID, TYPNAM, NITEMS, ARGLIS, NSIZE, ARGVAL, &
                   NWRDA8, ARGCHR, forg, nresp, narg, DR3VAL, senval, error)
    ! ----------------------------------------------------------------------
    !
    !     PURPOSE: COMPUTE THE EXTERNAL RESPONSE
    !
    !     GRPID   INPUT INTEGER        - GROUP ID
    !     TYPNAM  INPUT CHARACTER*8    - NAME OF EXTERNAL RESPONSE TYPE
    !     NITEMS  INPUT INTEGER        - DIMENSION OF ARRAY ARGLIS
    !     NSIZE   INPUT INTEGER        - DIMENSION OF ARRAY ARGVAL
    !     NWRDA8  INPUT INTEGER        - DIMENSION OF CHARACTER ARRAY ARGCHR
    !     ARGLIS  INPUT INTEGER        - ARRAY OF NO. OF ITEMS FOR EACH
    !                                     ARGUMENT type
    !     ARGVAL  INPUT DOUBLE         - ARRAY OF ARGUMENT VALUES (except
    !                                            characters)
    !     ARGCHR  INPUT CHARACTER*8    - ARRAY OF CHARACTER VALUES
    !     nresp   input integer        - number of responses
    !     forg    input integer        - type of call
    !                                     = 0 function evaluation
    !                                     = 1 sensitivity evaluation
    !     narg    input integer        - number of arguments needing gradients
    !     DR3VAL  OUTPUT DOUBLE        - VALUE OF THE EXTERNAL RESPONSEs
    !     senval  output double        - matrix of the sensitivity of the IRth
    !                                    response to the IARGth argument
    !     ERROR   INPUT/OUTPUT INTEGER -ERROR CODE FOR THE CALL.
    !                                   0 = PRINT ERROR MESSAGES
    !                                   1 = DO NOT PRINT ERROR MESSAGES.
    !
    !   METHOD
    !       A)SET UP VARIOUS PARAMETERS FROM THE ARGUMENT LIST
    !       B)if forg = 0 EVALUATE THE EXTERNAL RESPONSE BASED ON THE
    !                     GIVEN TYPNAM
    !       C)else if forg = 1 EVALUATE THE sensitivities of the external
    !                          responses to the arguments that can vary for
    !                          the given typnam
    !       D)RETURN BADTYP ERROR IF TYPNAM IS NOT MATCHED HERE.
    !
    !     nsize - the number of arguments or values in a dresp3 entry
    !
    !     nsize=nv+nc+nr+nnc+ndvp1+ndvp2+ndvc1+ndvc2+ndvm1+ndvm2+nrr2
    !        where:
    !        nv     = number of desvars    nc    = number of dtables
    !        nr     = number of dresp1s    nnc   = number of dnode pairs
    !        ndvp1  = number of dvprel1s   ndvp2 = number dvprel2s
    !        ndvc1  = number of dvcrel1s   ndvc2 = number dvcrel2s
    !        ndvm1  = number of dvmrel1s   ndvm2 = number dvmrel2s
    !        nrr2   = number of dresp2s
    !     narg = nsize - nc
    !
    !   CALLED BY
    !       SENDR3SVALD
    ! ----------------------------------------------------------------------
    !
    ! VARIABLES PASSED IN
    !
    implicit none

    character*8 TYPNAM, ARGCHR(NWRDA8)
    integer forg, nresp, narg
    integer GRPID, NITEMS, NSIZE, ARGLIS(NITEMS), ERROR, NWRDA8
    real*4 ARGVAL(NSIZE), DR3VAL(*), senval(nresp,*)
    integer BADTYP
    data BADTYP /7554/

    real*4, external :: cot, fSScr_skin, fSCcr_skin

    real*4 Nxx, Nxy, FC, FS, nu, Ec
    real*4 FScr, FCcr
    real*4 A_ratio, B_ratio, Rc, Rs

    real*4 r, t, a, b
    real*4 ms

    real*4 pi
    pi = atan(1.)*4

    ERROR = 0

    if (TYPNAM .eq. 'METHOD1') then
        ! r --> radius
        ! b --> circumferential distance between stringers
        ! a --> distance between frame rings

        ! Variables (DESVAR) from panel
        t = ARGVAL(1) ! Thickness of Panel 1

        ! Constants (DTABLE)
        r = ARGVAL(2) ! Curved panel radius
        a = ARGVAL(3) ! Panel length
        b = ARGVAL(4) ! Panel circumferential width

        Ec = ARGVAL(5) ! Compressive Young's module
        nu = ARGVAL(6) ! Poisson's ratio

        ! Mid-surface stress data (DRESP2)
        Nxx = ARGVAL(7) ! Nxx membrane force
        Nxy = ARGVAL(8) ! Nxy membrane force

        FC = Nxx/t ! FC compressive stress
        FS = Nxy/t ! FS shear stress

        FC = FC
        FS = ABS(FS)

        FCcr = fSCcr_skin(a, b, t, r, Ec, nu)
        FScr = fSScr_skin(a, b, t, r, Ec, nu)

        Rc = FC/FCcr
        Rs = FS/FScr

        ms = 2./(Rc + SQRT(Rc**2 + 4*Rs**2)) - 1.
    end if

    DR3VAL(1) = ms

    return
end subroutine R3SVALS

max_stars_repo_name: saullocastro/structmanager

max_issues_repo_name: saullocastro/structmanager

max_forks_repo_name: saullocastro/structmanager


text:
*> \brief \b DGET39
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
*
*       .. Scalar Arguments ..
*       INTEGER            KNT, LMAX, NINFO
*       DOUBLE PRECISION   RMAX
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGET39 tests DLAQTR, a routine for solving the real or
*> special complex quasi upper triangular system
*>
*>      op(T)*p = scale*c,
*> or
*>      op(T + iB)*(p+iq) = scale*(c+id),
*>
*> in real arithmetic. T is upper quasi-triangular.
*> If it is complex, then the first diagonal block of T must be
*> 1 by 1, B has the special structure
*>
*>                B = [ b(1) b(2) ... b(n) ]
*>                    [       w            ]
*>                    [           w        ]
*>                    [              .     ]
*>                    [                 w  ]
*>
*> op(A) = A or A', where A' denotes the conjugate transpose of
*> the matrix A.
*>
*> On input, X = [ c ].  On output, X = [ p ].
*>               [ d ]                  [ q ]
*>
*> Scale is an output less than or equal to 1, chosen to avoid
*> overflow in X.
*> This subroutine is specially designed for the condition number
*> estimation in the eigenproblem routine DTRSNA.
*>
*> The test code verifies that the following residual is order 1:
*>
*>      ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
*>    -----------------------------------------
*>        max(ulp*(||T||+||B||)*(||x1||+||x2||),
*>            (||T||+||B||)*smlnum/ulp,
*>            smlnum)
*>
*> (The (||T||+||B||)*smlnum/ulp term accounts for possible
*>  (gradual or nongradual) underflow in x1 and x2.)
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[out] RMAX
*> \verbatim
*>          RMAX is DOUBLE PRECISION
*>          Value of the largest test ratio.
*> \endverbatim
*>
*> \param[out] LMAX
*> \verbatim
*>          LMAX is INTEGER
*>          Example number where largest test ratio achieved.
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER
*>          Number of examples where INFO is nonzero.
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*>          KNT is INTEGER
*>          Total number of examples tested.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date December 2016
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
*
*  -- LAPACK test routine (version 3.7.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     December 2016
*
*     .. Scalar Arguments ..
      INTEGER            KNT, LMAX, NINFO
      DOUBLE PRECISION   RMAX
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            LDT, LDT2
      PARAMETER          ( LDT = 10, LDT2 = 2*LDT )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
     $                   NDIM
      DOUBLE PRECISION   BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
     $                   SCALE, SMLNUM, W, XNORM
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DLANGE
      EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           DCOPY, DGEMV, DLABAD, DLAQTR
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, COS, DBLE, MAX, SIN, SQRT
*     ..
*     .. Local Arrays ..
      INTEGER            IDIM( 6 ), IVAL( 5, 5, 6 )
      DOUBLE PRECISION   B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
     $                   VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
     $                   VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
*     ..
*     .. Data statements ..
      DATA               IDIM / 4, 5*5 /
      DATA               IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
     $                   4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
     $                   0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
     $                   0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
     $                   4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
     $                   -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
     $                   5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
     $                   4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
*     ..
*     .. Executable Statements ..
*
*     Get machine parameters
*
      EPS = DLAMCH( 'P' )
      SMLNUM = DLAMCH( 'S' )
      BIGNUM = ONE / SMLNUM
      CALL DLABAD( SMLNUM, BIGNUM )
*
*     Set up test case parameters
*
      VM1( 1 ) = ONE
      VM1( 2 ) = SQRT( SMLNUM )
      VM1( 3 ) = SQRT( VM1( 2 ) )
      VM1( 4 ) = SQRT( BIGNUM )
      VM1( 5 ) = SQRT( VM1( 4 ) )
*
      VM2( 1 ) = ONE
      VM2( 2 ) = SQRT( SMLNUM )
      VM2( 3 ) = SQRT( VM2( 2 ) )
      VM2( 4 ) = SQRT( BIGNUM )
      VM2( 5 ) = SQRT( VM2( 4 ) )
*
      VM3( 1 ) = ONE
      VM3( 2 ) = SQRT( SMLNUM )
      VM3( 3 ) = SQRT( VM3( 2 ) )
      VM3( 4 ) = SQRT( BIGNUM )
      VM3( 5 ) = SQRT( VM3( 4 ) )
*
      VM4( 1 ) = ONE
      VM4( 2 ) = SQRT( SMLNUM )
      VM4( 3 ) = SQRT( VM4( 2 ) )
      VM4( 4 ) = SQRT( BIGNUM )
      VM4( 5 ) = SQRT( VM4( 4 ) )
*
      VM5( 1 ) = ONE
      VM5( 2 ) = EPS
      VM5( 3 ) = SQRT( SMLNUM )
*
*     Initialization
*
      KNT = 0
      RMAX = ZERO
      NINFO = 0
      SMLNUM = SMLNUM / EPS
*
*     Begin test loop
*
      DO 140 IVM5 = 1, 3
         DO 130 IVM4 = 1, 5
            DO 120 IVM3 = 1, 5
               DO 110 IVM2 = 1, 5
                  DO 100 IVM1 = 1, 5
                     DO 90 NDIM = 1, 6
*
                        N = IDIM( NDIM )
                        DO 20 I = 1, N
                           DO 10 J = 1, N
                              T( I, J ) = DBLE( IVAL( I, J, NDIM ) )*
     $                                    VM1( IVM1 )
                              IF( I.GE.J )
     $                           T( I, J ) = T( I, J )*VM5( IVM5 )
   10                      CONTINUE
   20                   CONTINUE
*
                        W = ONE*VM2( IVM2 )
*
                        DO 30 I = 1, N
                           B( I ) = COS( DBLE( I ) )*VM3( IVM3 )
   30                   CONTINUE
*
                        DO 40 I = 1, 2*N
                           D( I ) = SIN( DBLE( I ) )*VM4( IVM4 )
   40                   CONTINUE
*
                        NORM = DLANGE( '1', N, N, T, LDT, WORK )
                        K = IDAMAX( N, B, 1 )
                        NORMTB = NORM + ABS( B( K ) ) + ABS( W )
*
                        CALL DCOPY( N, D, 1, X, 1 )
                        KNT = KNT + 1
                        CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
     $                               DUMM, SCALE, X, WORK, INFO )
                        IF( INFO.NE.0 )
     $                     NINFO = NINFO + 1
*
*                       || T*x - scale*d || /
*                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
*
                        CALL DCOPY( N, D, 1, Y, 1 )
                        CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
     $                              X, 1, -SCALE, Y, 1 )
                        XNORM = DASUM( N, X, 1 )
                        RESID = DASUM( N, Y, 1 )
                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
     $                          ( NORM*EPS )*XNORM )
                        RESID = RESID / DOMIN
                        IF( RESID.GT.RMAX ) THEN
                           RMAX = RESID
                           LMAX = KNT
                        END IF
*
                        CALL DCOPY( N, D, 1, X, 1 )
                        KNT = KNT + 1
                        CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
     $                               DUMM, SCALE, X, WORK, INFO )
                        IF( INFO.NE.0 )
     $                     NINFO = NINFO + 1
*
*                       || T*x - scale*d || /
*                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
*
                        CALL DCOPY( N, D, 1, Y, 1 )
                        CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
     $                              1, -SCALE, Y, 1 )
                        XNORM = DASUM( N, X, 1 )
                        RESID = DASUM( N, Y, 1 )
                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
     $                          ( NORM*EPS )*XNORM )
                        RESID = RESID / DOMIN
                        IF( RESID.GT.RMAX ) THEN
                           RMAX = RESID
                           LMAX = KNT
                        END IF
*
                        CALL DCOPY( 2*N, D, 1, X, 1 )
                        KNT = KNT + 1
                        CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
     $                               SCALE, X, WORK, INFO )
                        IF( INFO.NE.0 )
     $                     NINFO = NINFO + 1
*
*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
*                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
*                                  smlnum/ulp * (||T||+||B||), smlnum )
*
*
                        CALL DCOPY( 2*N, D, 1, Y, 1 )
                        Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) +
     $                           SCALE*Y( 1 )
                        DO 50 I = 2, N
                           Y( I ) = W*X( I+N ) + SCALE*Y( I )
   50                   CONTINUE
                        CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
     $                              X, 1, -ONE, Y, 1 )
*
                        Y( 1+N ) = DDOT( N, B, 1, X, 1 ) -
     $                             SCALE*Y( 1+N )
                        DO 60 I = 2, N
                           Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
   60                   CONTINUE
                        CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
     $                              X( 1+N ), 1, ONE, Y( 1+N ), 1 )
*
                        RESID = DASUM( 2*N, Y, 1 )
                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
     $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
                        RESID = RESID / DOMIN
                        IF( RESID.GT.RMAX ) THEN
                           RMAX = RESID
                           LMAX = KNT
                        END IF
*
                        CALL DCOPY( 2*N, D, 1, X, 1 )
                        KNT = KNT + 1
                        CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
     $                               SCALE, X, WORK, INFO )
                        IF( INFO.NE.0 )
     $                     NINFO = NINFO + 1
*
*                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
*                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
*                                  smlnum/ulp * (||T||+||B||), smlnum )
*
                        CALL DCOPY( 2*N, D, 1, Y, 1 )
                        Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
                        DO 70 I = 2, N
                           Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
     $                              SCALE*Y( I )
   70                   CONTINUE
                        CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
     $                              1, ONE, Y, 1 )
*
                        Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
                        DO 80 I = 2, N
                           Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
     $                                SCALE*Y( I+N )
   80                   CONTINUE
                        CALL DGEMV( 'Transpose', N, N, ONE, T, LDT,
     $                              X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
*
                        RESID = DASUM( 2*N, Y, 1 )
                        DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
     $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
                        RESID = RESID / DOMIN
                        IF( RESID.GT.RMAX ) THEN
                           RMAX = RESID
                           LMAX = KNT
                        END IF
*
   90                CONTINUE
  100             CONTINUE
  110          CONTINUE
  120       CONTINUE
  130    CONTINUE
  140 CONTINUE
*
      RETURN
*
*     End of DGET39
*
      END

ext: f

max_stars_repo_name: dnoan/OpenBLAS

max_issues_repo_name: dnoan/OpenBLAS

max_forks_repo_name: dnoan/OpenBLAS


text:
! gro2cfg
!
! Created on: Mar 9, 2012
! Author: Pedro Duarte and Pedro Morgado
! License: MIT License
! Version: 0.2.1
! Description: Convert a GROMACS .gro file into a DL Poly CONFIG file.

program gro2cfg

	implicit none
	! Variable declaration
	character :: arg*20, filename*80, outfilename*80, name*20, molecule*8, fmt*20, fmt_box*20, skip, atom*8
	integer :: narg, i, line, totlines, indexno
	real :: xbox, ybox, zbox, ax, ay, az, nullval=0.

	! Handle command line arguments
	! Check if any arguments are found
	narg = command_argument_count()
	! Loop over the arguments
	if(narg>0)then
		! loop across options
		do i=1,narg
			call get_command_argument(i,arg)
			select case(adjustl(arg))
				case("--help", "-h")
					write(*,'(A33,(/))') "gro2cfg"
					write(*,*) "Transform .gro (GROMACS) into a CONFIG file (DL_POLY)."
					stop
				case("--version", "-v")
					write(*,*) "v0.2.1"
					stop
				case default
					write(*,*)"Option unknown: ", adjustl(arg)
					stop
			end select
		end do
	end if

	!Ask user for the path to the file to convert and output filename
	!print *, "Enter file path without the file (can be relative to this file location):"
	!read *, filepath
	filename = ""
    do while (filename == "")
        print *, "Enter filename (mustn't be null):"
        read *, filename
    end do
    print *, "Enter output filename:"
    read *, outfilename
	! Set filename to CONFIG if null input
	if (outfilename == "") filename="CONFIG"

	open(7, file=filename, status="old", action="read") ! To determine number of lines
	! Find number of lines of file and get box dimensions
	do line=1, 3
	    if (line == 1) then
	        read(7,*) name
	    else if (line == 2) then
	        read(7,*) totlines
	        totlines = totlines + 2 + 1 !Total lines is equal to all atoms plus 2 on top and 1 at the bottom
	    end if
	end do
	do line=3, totlines
	    if (line == totlines) then
            read(7,*) xbox, ybox, zbox
            xbox = xbox*10
            ybox = ybox*10
            zbox = zbox*10
        else
            read(7,*) skip
        end if
    end do

	close(7)
	open(8, file=filename, status="old", action="read") ! To do the actual conversion
	open(9, file=outfilename, status="replace", action="write") ! File to write output

	! Read positions and make CONFIG file
	fmt="(3(f20.12))"
    fmt_box="(3(f20.12))"
	!File HEADER
	write(9,'(A80)') name
	write(9,*) "        0","         1"
	write(9,fmt_box) xbox, nullval, nullval
    write(9,fmt_box) nullval, ybox, nullval
    write(9,fmt_box) nullval, nullval, zbox

	!File BODY
	do line=1, totlines
        if (line == 1 .or. line == 2 .or. line == totlines) then
            read(8,*)
        else
            read(8,*) molecule, atom, indexno, ax, ay, az
            ! Convert to Angstroms
            ax = ax*10
            ay = ay*10
            az = az*10

            if (ax<0) then
                ax = ax+xbox
            else if (ax>xbox) then
                ax = ax-xbox
            end if

	        if (ay<0) then
                ay = ay+ybox
            else if (ay>ybox) then
                ay = ay-ybox
            end if

            if (az<0) then
                az = az+zbox
            else if (az>zbox) then
                az = az-zbox
            end if

            ax = ax-(xbox/2)
            ay = ay-(ybox/2)
            az = az-(zbox/2)
            write(9,'(A8,I10)') atom, indexno
            write(9,fmt) ax, ay, az
        end if
	end do

	close(8)
	close(9)

end program gro2cfg

max_stars_repo_name: xipasduarte/simutil

max_issues_repo_name: xipasduarte/simutil

max_forks_repo_name: xipasduarte/simutil


text:
Trisha has lived in Davis for a long time.

Trisha has a son named Matt, who works at Hibbert Lumber.

Trishas husband works at Ganesh Fixit or otherwise known as Ganesh Works, he has worked there a long time.

Matt no longer works at Hibbert he now lives in San Francisco and is currently a bum.

Matt is not such bumish anymore he is now a sandwich artiste.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,SLAB,TISO,IFLAG)
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
C
C
C DIMENSION OF           T(LU,LV,MW),EYE(3),SLAB(MUVWP2,MUVWP2)
C ARGUMENTS
C
C LATEST REVISION        DECEMBER 1984
C
C PURPOSE                ISOSRF DRAWS AN APPROXIMATION OF AN ISO-VALUED
C                        SURFACE FROM A THREE-DIMENSIONAL ARRAY WITH
C                        HIDDEN LINES REMOVED.
C
C USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
C
C                            CALL EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
C
C                               ASSUMPTIONS:
C                                 -- ALL OF THE T ARRAY IS TO BE USED.
C                                 -- IFLAG IS CHOSEN INTERNALLY.
C                                 -- FRAME IS CALLED BY EZISOS.
C
C                        IF THE ASSUMPTIONS ARE NOT MET, USE
C
C                            CALL ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,
C                                         SLAB,TISO,IFLAG)
C
C ARGUMENTS
C
C ON INPUT               T
C                          THREE DIMENSIONAL ARRAY OF DATA THAT DEFINES
C                          THE ISO-VALUED SURFACE.
C
C                        LU
C                          FIRST DIMENSION OF T IN THE CALLING PROGRAM.
C
C                        MU
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE U DIRECTION (THE FIRST
C                          SUBSCRIPT DIRECTION).  WHEN PROCESSING THE
C                          ENTIRE ARRAY, LU = MU (AND LV = MV).
C
C                        LV
C                          SECOND DIMENSION OF T IN THE CALLING PROGRAM.
C
C                        MV
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE V DIRECTION (THE SECOND
C                          SUBSCRIPT DIRECTION).
C
C                        MV
C                          THE NUMBER OF DATA VALUES OF T TO BE
C                          PROCESSED IN THE W DIRECTION (THE THIRD
C                          SUBSCRIPT DIRECTION).
C
C                        EYE
C                          THE POSITION OF THE EYE IN THREE-SPACE.  T IS
C                          CONSIDERED TO BE IN A BOX WITH OPPOSITE
C                          CORNERS (1,1,1) AND (MU,MV,MW).  THE EYE IS
C                          AT (EYE(1),EYE(2),EYE(3)), WHICH MUST BE
C                          OUTSIDE THE BOX THAT CONTAINS T.  WHILE GAINING
C                          EXPERIENCE WITH THE ROUTINE, A GOOD CHOICE
C                          FOR EYE MIGHT BE (5.0*MU,3.5*MV,2.0*MW).
C
C                        MUVWP2
C                          THE MAXIMUM OF (MU,MV,MW)+2; THAT IS,
C                          MUVWP2 = MAX(MU,MV,MW)+2).
C
C                        SLAB
C                          A WORK SPACE USED FOR INTERNAL STORAGE.  SLAB
C                          MUST BE AT LEAST MUVWP2*MUVWP2 WORDS LONG.
C
C                        TISO
C                          THE ISO-VALUE USED TO DEFINE THE SURFACE.  THE
C                          SURFACE DRAWN WILL SEPARATE VOLUMES OF T THAT
C                          HAVE VALUES GREATER THAN OR EQUAL TO TISO FROM
C                          VOLUMES OF T THAT HAVE VALUES LESS THAN TISO.
C
C                        IFLAG
C                          THIS FLAG SERVES TWO PURPOSES.
C                          .  FIRST, THE ABSOLUTE VALUE OF IFLAG
C                             DETERMINES WHICH TYPES OF LINES ARE DRAWN
C                             TO APPROXIMATE THE SURFACE.  THREE TYPES
C                             OF LINES ARE CONSIDERED:  LINES OF
C                             CONSTANT U, LINES OF CONSTANT V AND LINES
C                             OF CONSTANT W.  THE FOLLOWING TABLE LISTS
C                             THE TYPES OF LINES DRAWN.
C
C                                                LINES OF CONSTANT
C                                                -----------------
C                             IABS(IFLAG)           U    V    W
C                             -----------          ---  ---  ---
C                                  1               NO   NO   YES
C                                  2               NO   YES  NO
C                                  3               NO   YES  YES
C                                  4               YES  NO   NO
C                                  5               YES  NO   YES
C                                  6               YES  YES  NO
C                             0, 7 OR MORE         YES  YES  YES
C
C                          .  SECOND, THE SIGN OF IFLAG DETERMINES WHAT
C                             IS INSIDE AND WHAT IS OUTSIDE, HENCE,
C                             WHICH LINES ARE VISIBLE AND WHAT IS DONE
C                             AT THE BOUNDARY OF T.  FOR IFLAG:
C
C                             POSITIVE   T VALUES GREATER THAN TISO ARE
C                                        ASSUMED TO BE INSIDE THE SOLID
C                                        FORMED BY THE DRAWN SURFACE.
C                             NEGATIVE   T VALUES LESS THAN TISO ARE
C                                        ASSUMED TO BE INSIDE THE SOLID
C                                        FORMED BY THE DRAWN SURFACE.
C                             IF THE ALGORITHM DRAWS A CUBE, REVERSE THE
C                             SIGN OF IFLAG.
C
C ON OUTPUT              T,LU,MU,LV,MV,MW,EYE,MUVWP2,TISO AND IFLAG ARE
C                        UNCHANGED.  SLAB HAS BEEN WRITTEN IN.
C
C NOTE                   .   THIS ROUTINE IS FOR LOWER RESOLUTION ARRAYS
C                            THAN ISOSRFHR.  40 BY 40 BY 40 IS A
C                            PRACTICAL MAXIMUM.
C                        .   TRANSFORMATIONS CAN BE ACHIEVED BY
C                            ADJUSTING SCALING STATEMENT FUNCTIONS IN
C                            ISOSRF, SET3D AND TR32.
C                        .   THE HIDDEN-LINE ALGORITHM IS NOT EXACT, SO
C                            VISIBILITY ERRORS CAN OCCUR.
C                        .   THREE-DIMENSIONAL PERSPECTIVE CHARACTER
C                            LABELING OF ISOSRF IS POSSIBLE BY USING
C                            THE UTILITY PWRZI.  FOR A DESCRIPTION OF
C                            THE USAGE, SEE THE PWRZI DOCUMENTATION.
C
C ENTRY POINTS           ISOSRF, EZISOS, SET3D, TRN32I, ZEROSC,
C                        STCNTR, DRCNTR, TR32, FRSTS, KURV1S, KURV2S,
C                        FRSTC, FILLIN, DRAWI, ISOSRB, MMASK
C
C COMMON BLOCKS          ISOSR1, ISOSR2, ISOSR3, ISOSR4, ISOSR5,
C                        ISOSR6, ISOSR7, ISOSR8, ISOSR9, TEMPR,
C                        PWRZ1I
C
C REQUIRED LIBRARY       THE ERPRT77 PACKAGE AND THE SPPS.
C ROUTINES
C
C I/O                    PLOTS SURFACE
C
C PRECISION              SINGLE
C
C LANGUAGE               FORTRAN 77
C
C HISTORY                DEVELOPED FOR USERS OF ISOSRFHR WITH SMALLER
C                        ARRAYS.
C
C ALGORITHM              CUTS THROUGH THE THREE-DIMENSIONAL ARRAY ARE
C                        CONTOURED WITH A SMOOTHING CONTOURER WHICH ALSO
C                        MARKS A MODEL OF THE PLOTTING PLANE.  INTERIORS
C                        OF BOUNDARIES ARE FILLED IN AND THE RESULT IS
C                        .OR.ED INTO ANOTHER MODEL OF THE PLOTTING PLANE
C                        WHICH IS USED TO TEST SUBSEQUENT CONTOUR LINES
C                        FOR VISIBILITY.
C
C TIMING                 VARIES WIDELY WITH SIZE OF T AND THE VOLUME OF
C                        THE SPACE ENCLOSED BY THE SURFACE DRAWN.
C
C  **NOTE**                  SPACE REQUIREMENTS CAN BE REDUCED BY
C                            CHANGING THE SIZE OF THE ARRAYS ISCR, ISCA
C                            (FOUND IN COMMON ISOSR2), MASK(FOUND IN
C                            COMMON ISOSR5) AND THE VARIABLE NBPW
C                            (COMMON ISOSR5).
C                            ISCR AND ISCA NEED 128X128 BITS.  SO ON A
C                            64 BIT MACHINE ISCR, ISCA CAN BE
C                            DIMENSIONED TO (2,128). NBPW SET IN
C                            SUBROUTINE MMASK SHOULD CONTAIN THE
C                            NUMBER OF BITS PER WORD YOU WISH TO
C                            UTILIZE.
C                            THE DIMENSION OF MASK AND NMASK SHOULD
C                            EQUAL THE VALUE OF NBPW.
C                            LS SHOULD BE SET TO THE FIRST DIMENSION
C                            OF ISCA AND ISCR.
C
C               EXAMPLES:
C                            ON A 60 BIT MACHINE:
C                                 DIMENSION ISCA(4,128), ISCR(4,128)
C                                 DIMENSION MASK(32)
C                                 NBPW = 32
C                            ON A 64 BIT MACHINE:
C                                 DIMENSION ISCA(2,128), ISCR(2,128)
C                                 DIMENSION MASK(64)
C                                 NBPW = 64
C
C INTERNAL PARAMETERS    NAME  DEFAULT           FUNCTION
C                        ----  -------           --------
C                        IREF     1     FLAG TO CONTROL DRAWING OF AXES.
C                                       .IREF=NONZERO   DRAW AXES.
C                                       .IREF=ZERO      DO NOT DRAW AXES.
C
C
      SAVE
      DIMENSION       T(LU,LV,MW),EYE(3)     ,SLAB(MUVWP2,MUVWP2)
C
      COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
      COMMON /ISOSR4/ RX         ,RY
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR9/ BIG        ,IXBIT
C
      LOGICAL         GENDON
      DATA IREF/1/
C
      AVE(A,B) = (A+B)*.5
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
C +NOAO - Blockdata ISOSRB rewritten as run time initialization
C     EXTERNAL ISOSRB
      call isosrb
C -NOAO
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('NSSL','ISOSRF','ISOSRF','VERSION 12')
      NERR = 0
C
C 3-SPACE  U,V,W,IU,IV,IW,ETC
C 2-SPACE  X,Y,IX,IY,ETC
C
C  INITIALIZE MASKS
C
      IF (.NOT.GENDON) CALL MMASK
C
C  SET SHIFT VALUE FOR X,Y PACKING
C
C  IF YOUR MACHINE HAS MORE THAN 16 BITS PER WORD THIS CHECK MAY BE
C  MODIFIED
C
      IF (LU .LE. 256) GO TO  10
      NERR = NERR + 1
      CALL SETER('DIMENSION OF CUBE EXCEEDS 256',NERR,2)
      RETURN
   10 DO  20 J=1,30
         IF (LU .LE. 2**(J-1)) GO TO  30
   20 CONTINUE
   30 IXBIT = J
      NU = MU
      NUP2 = NU+2
      NV = MV
      NVP2 = NV+2
      NW = MW
      NWP2 = NW+2
      FNU = NU
      FNV = NV
      FNW = NW
      SU1 = SU(1.)
      SV1 = SV(1.)
      SW1 = SW(1.)
      SUNU = SU(FNU)
      SVNV = SV(FNV)
      SWNW = SW(FNW)
      AVEU = AVE(SU1,SUNU)
      AVEV = AVE(SV1,SVNV)
      AVEW = AVE(SW1,SWNW)
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
      NUVWP2 = MUVWP2
      TVAL = TISO
      NFLAG = IABS(IFLAG)
      IF (NFLAG.EQ.0 .OR. NFLAG.GE.8) NFLAG = 7
C
C SET UP SCALING
C
      FACT = -ISIGN(1,IFLAG)
      CALL SET3D (EYE,1.,FNU,1.,FNV,1.,FNW)
C
C BOUND LOWER AND LEFT EDGE OF SLAB
C
      EDGE = SIGN(BIG,FACT)
      DO  40 IUVW=1,NUVWP2
         SLAB(IUVW,1) = EDGE
         SLAB(1,IUVW) = EDGE
   40 CONTINUE
C
C SLICES PERPENDICULAR TO U. THAT IS, V W SLICES. T OF CONSTANT U.
C
      IF (NFLAG .LT. 4) GO TO 100
      CALL ZEROSC
      ISLBT = -1
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO  50 IV=2,NVP2
         SLAB(IV,NWP2) = EDGE
   50 CONTINUE
      DO  60 IW=2,NWP2
         SLAB(NVP2,IW) = EDGE
   60 CONTINUE
C
C GO THRU 3-D ARRAY IN U DIRECTION.  IUEW=IU EITHER WAY.
C PICK IU BASED ON EYEU.
C
      DO  90 IUEW=1,NU
         IU = IUEW
         IF (EYEU .GT. AVEU) IU = NU+1-IUEW
         U = IU
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO  80 IV=1,NV
            DO  70 IW=1,NW
               SLAB(IV+1,IW+1) = T(IU,IV,IW)
   70       CONTINUE
   80    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NVP2,NWP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
   90 CONTINUE
C
C SLICES PERPENDICULAR TO V. U W SLICES. T OF CONSTANT V.
C
  100 IF (MOD(NFLAG/2,2) .EQ. 0) GO TO 160
      CALL ZEROSC
      ISLBT = 0
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO 110 IU=2,NUP2
         SLAB(IU,NWP2) = EDGE
  110 CONTINUE
      DO 120 IW=2,NWP2
         SLAB(NUP2,IW) = EDGE
  120 CONTINUE
C
C GO THRU T IN V DIRECTION.  IVEW=IV EITHER WAY.
C
      DO 150 IVEW=1,NV
         IV = IVEW
         IF (EYEV .GT. AVEV) IV = NV+1-IVEW
         V = IV
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO 140 IU=1,NU
            DO 130 IW=1,NW
               SLAB(IU+1,IW+1) = T(IU,IV,IW)
  130       CONTINUE
  140    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NUP2,NWP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
  150 CONTINUE
C
C SLICES PERPENDICULAR TO W. U V SLICES. T OF CONSTANT W.
C
  160 IF (MOD(NFLAG,2) .EQ. 0) GO TO 220
      CALL ZEROSC
C
      ISLBT = 1
C
C BOUND UPPER AND RIGHT EDGE OF SLAB.
C
      DO 170 IU=2,NUP2
         SLAB(IU,NVP2) = EDGE
  170 CONTINUE
      DO 180 IV=2,NVP2
         SLAB(NUP2,IV) = EDGE
  180 CONTINUE
C
C GO THRU T IN W DIRECTION.
C
      DO 210 IWEW=1,NW
         IW = IWEW
         IF (EYEW .GT. AVEW) IW = NW+1-IWEW
         W = IW
C
C LOAD THIS SLICE OF T INTO SLAB.
C
         DO 200 IU=1,NU
            DO 190 IV=1,NV
               SLAB(IU+1,IV+1) = T(IU,IV,IW)
  190       CONTINUE
  200    CONTINUE
C
C CONTOUR THIS SLAB.
C
         CALL STCNTR (SLAB,NUVWP2,NUP2,NVP2,TVAL)
C
C CONSTRUCT VISIBILITY ARRAY.
C
         CALL FILLIN
  210 CONTINUE
C
C DRAW REFERENCE PLANE EDGES AND W AXIS.
C
  220 IF (IREF .EQ. 0) RETURN
      CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
      IF (EYEV .LT. SV1) GO TO 240
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 230 IU=2,NU
         CALL TRN32I (SU(FLOAT(IU)),SV1,SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  230 CONTINUE
      GO TO 250
  240 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SUNU,SV1,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  250 IF (EYEU .GT. SUNU) GO TO 270
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 260 IV=2,NV
         CALL TRN32I (SUNU,SV(FLOAT(IV)),SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  260 CONTINUE
      GO TO 280
  270 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SUNU,SVNV,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  280 IF (EYEV .GT. SVNV) GO TO 300
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 290 IUOW=2,NU
         CALL TRN32I (SU(FLOAT(NU-IUOW+1)),SVNV,SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  290 CONTINUE
      GO TO 310
  300 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SVNV,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  310 IF (EYEU .LT. SU1) GO TO 330
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 320 IVOW=2,NV
         CALL TRN32I (SU1,SV(FLOAT(NV-IVOW+1)),SW1,XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  320 CONTINUE
      GO TO 340
  330 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  340 IF (EYEU.LE.SU1 .OR. EYEV.LE.SV1) GO TO 360
      CALL FRSTC (IFIX(XT),IFIX(YT),1)
      DO 350 IW=2,NW
         CALL TRN32I (SU1,SV1,SW(FLOAT(IW)),XT,YT,DUM,2)
         CALL FRSTC (IFIX(XT),IFIX(YT),2)
  350 CONTINUE
C +NOAO - Plotit buffer needs to be flushed before returning.
      call plotit (0, 0, 2)
C -NOAO
      RETURN
  360 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
      CALL TRN32I (SU1,SV1,SWNW,XT,YT,DUM,2)
      CALL PLOTIT (IFIX(XT),IFIX(YT),1)
C +NOAO - Plotit buffer needs to be flushed before returning.
      call plotit (0, 0, 2)
C -NOAO
      RETURN
      END
      SUBROUTINE EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
C
      SAVE
      DIMENSION       T(MU,MV,MW),EYE(3)
C
        DATA ANG,PI/.35,3.141592/
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('NSSL','ISOSRF','EZISOS','VERSION 12')
C
C ARGUMENTS DESCRIBED IN ISOSRF
C
C PICK TYPES OF LINES TO DRAW
C
      NU = MU
      NV = MV
      NW = MW
      TVAL = TISO
      MAX = MAX0(NU,NV,NW)+2
      ATU = NU/2
      ATV = NV/2
      ATW = NW/2
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
      RU = EYEU-ATU
      RV = EYEV-ATV
      RW = EYEW-ATW
      RU2 = RU*RU
      RV2 = RV*RV
      RW2 = RW*RW
      DU = SQRT(RV2+RW2)
      DV = SQRT(RU2+RW2)
      DW = SQRT(RU2+RV2)
      DR = 1./SQRT(RU2+RV2+RW2)
C
C  COMPUTE THE ARCCOSINE
C
      TU = DU*DR
      ANGU = ATAN(ABS(SQRT(1.-TU*TU)/TU))
      IF (TU .LE. 0.) ANGU = PI-ANGU
      TV = DV*DR
      ANGV = ATAN(ABS(SQRT(1.-TV*TV)/TV))
      IF (TV .LE. 0.) ANGV = PI-ANGV
      TW = DW*DR
      ANGW = ATAN(ABS(SQRT(1.-TW*TW)/TW))
      IF (TW .LE. 0.) ANGW = PI-ANGW
C
C BREAK POINT IS ABOUT 20 DEGREES OR ABOUT .35 RADIANS
C
      IFLAG = 0
      IF (ANGU .GT. ANG) IFLAG = IFLAG+4
      IF (ANGV .GT. ANG) IFLAG = IFLAG+2
      IF (ANGW .GT. ANG) IFLAG = IFLAG+1
C
C FIND SIGN OF IFLAG
C
      ICNT = 0
      IF (ABS(RU) .LE. ATU) GO TO  30
      IU = 1
      IF (EYEU .GT. ATU) IU = NU
      DO  20 IW=1,NW
         DO  10 IV=1,NV
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   10    CONTINUE
   20 CONTINUE
   30 IF (ABS(RV) .LE. ATV) GO TO  60
      IV = 1
      IF (EYEV .GT. ATV) IV = NV
      DO  50 IW=1,NW
         DO  40 IU=1,NU
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   40    CONTINUE
   50 CONTINUE
   60 IF (ABS(RW) .LE. ATW) GO TO  90
      IW = 1
      IF (EYEW .GT. ATW) IW = NW
      DO  80 IV=1,NV
         DO  70 IU=1,NU
            IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
            ICNT = ICNT+1
   70    CONTINUE
   80 CONTINUE
   90 IFLAG = ISIGN(IFLAG,ICNT)
      CALL ISOSRF (T,NU,NU,NV,NV,NW,EYE,MAX,SLAB,TVAL,IFLAG)
C +NOAO - Call to frame is suppressed.
C     CALL FRAME
C -NOAO
      RETURN
      END
      SUBROUTINE SET3D (EYE,ULO,UHI,VLO,VHI,WLO,WHI)
      SAVE
      COMMON /TEMPR/  RZERO
C
      DIMENSION       EYE(3)
C
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
      COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
     1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
     2                EYEV       ,EYEW
C
C
      AVE(A,B) = (A+B)*.5
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
C CONSTANTS FOR PWRZ
C
      UUMIN = ULO
      UUMAX = UHI
      VVMIN = VLO
      VVMAX = VHI
      WWMIN = WLO
      WWMAX = WHI
      EYEU = EYE(1)
      EYEV = EYE(2)
      EYEW = EYE(3)
C
C FIND CORNERS IN 2-SPACE FOR 3-SPACE BOX CONTAINING OBJECT
C
      ISCALE = 0
      ATU = AVE(SU(UUMIN),SU(UUMAX))
      ATV = AVE(SV(VVMIN),SV(VVMAX))
      ATW = AVE(SW(WWMIN),SW(WWMAX))
      BIGD = 0.
      IF (RZERO .LE. 0.) GO TO  10
C
C RELETIVE SIZE FEATURE IN USE.
C GENERATE EYE POSITION THAT MAKES BOX HAVE MAXIMUM PROJECTED SIZE.
C
      ALPHA = -(VVMIN-ATV)/(UUMIN-ATU)
      VVEYE = -RZERO/SQRT(1.+ALPHA*ALPHA)
      UUEYE = VVEYE*ALPHA
      VVEYE = VVEYE+ATV
      UUEYE = UUEYE+ATU
      WWEYE = ATW
      CALL TRN32I (ATU,ATV,ATW,UUEYE,VVEYE,WWEYE,1)
      CALL TRN32I (UUMIN,VVMIN,ATW,XMIN,DUMM,DUMM,2)
      CALL TRN32I (UUMAX,VVMIN,WWMIN,DUMM,YMIN,DUMM,2)
      CALL TRN32I (UUMAX,VVMAX,ATW,XMAX,DUMM,DUMM,2)
      CALL TRN32I (UUMAX,VVMIN,WWMAX,DUMM,YMAX,DUMM,2)
      BIGD = SQRT((UUMAX-UUMIN)**2+(VVMAX-VVMIN)**2+(WWMAX-WWMIN)**2)*.5
      R0 = RZERO
      GO TO  20
   10 CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
      CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMIN),X1,Y1,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMAX),X2,Y2,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMIN),X3,Y3,DUM,2)
      CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMAX),X4,Y4,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMIN),X5,Y5,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMAX),X6,Y6,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMIN),X7,Y7,DUM,2)
      CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMAX),X8,Y8,DUM,2)
      XMIN = AMIN1(X1,X2,X3,X4,X5,X6,X7,X8)
      XMAX = AMAX1(X1,X2,X3,X4,X5,X6,X7,X8)
      YMIN = AMIN1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
      YMAX = AMAX1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
C
C ADD RIGHT AMOUNT TO KEEP PICTURE SQUARE
C
   20 WIDTH = XMAX-XMIN
      HIGHT = YMAX-YMIN
      DIF = .5*(WIDTH-HIGHT)
      IF (DIF)  30, 50, 40
   30 XMIN = XMIN+DIF
      XMAX = XMAX-DIF
      GO TO  50
   40 YMIN = YMIN-DIF
      YMAX = YMAX+DIF
   50 ISCALE = 1
      CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
      RETURN
      END
      SUBROUTINE TRN32I (U,V,W,XT,YT,ZT,IENT)
C
C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR-
C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE
C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15,
C 2, 193-204,1968.
C ARGUMENTS FOR SET
C U,V,W    ARE THE 3-SPACE COORDINATES OF THE INTERSECTION
C          OF THE LINE OF SIGHT AND THE IMAGE PLANE.  THIS
C          POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT.
C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION.
C
C TRN32 ARGUMENTS
C U,V,W    ARE THE 3-SPACE COORDINATES OF A POINT TO BE
C          TRANSFORMED.
C XT,YT    THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR-
C          MATION.  WHEN ISCALE=0, XT AND YT ANR IN THE SAME
C          UNITS AS U,V, AND W.  WHEN ISCALE'0, XT AND YT
C          ARE IN PLOTTER COORDINATES.
C ZT       NOT USED.
C
      SAVE
      COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
     1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
     2                EYEV       ,EYEW
      COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
     1                YMAX       ,BIGD       ,R0
C
C RANGE OF PLOTTER COORDINATES
C
C
C  WARNING
C  IF PLOTTER MAXIMUM VALUE RANGES (IN X OR Y DIRECTION) FALL BELOW
C  101, THEN CHANGES MUST BE MADE IN SUBROUTINE FRSTC. THE REQUIRED
C  CHANGES ARE MARKED BY WARNING COMMENTS IN FRSTC.
      DATA NLX,NBY,NRX,NTY/10,10,32760,32760/
        DATA PI/3.1411592/
C
C STORE THE PARAMETERS OF THE SET CALL FOR USE
C WITH THE TRANSLATE CALL
C
C DECIDE IF SET OR TRANSLATE CALL
C
      IF (IENT .NE. 1) GO TO  50
      AU = U
      AV = V
      AW = W
      EU = XT
      EV = YT
      EW = ZT
C
C
C
C
C
      DU = AU-EU
      DV = AV-EV
      DW = AW-EW
      D = SQRT(DU*DU+DV*DV+DW*DW)
      COSAL = DU/D
      COSBE = DV/D
      COSGA = DW/D
C
C  COMPUTE THE ARCCOSINE
C
      AL = ATAN(ABS(SQRT(1.-COSAL*COSAL)/COSAL))
      IF (COSAL .LE. 0.) AL = PI-AL
      BE = ATAN(ABS(SQRT(1.-COSBE*COSBE)/COSBE))
      IF (COSBE .LE. 0.) BE = PI-BE
      GA = ATAN(ABS(SQRT(1.-COSGA*COSGA)/COSGA))
      IF (COSGA .LE. 0.) GA = PI-GA
      SINGA = SIN(GA)
C
C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF
C THE 2-SPACE.  THE 3-SPACE W AXIS IS TRANSFORMED INTO THE
C 2-SPACE Y AXIS.  IF THE LINE OF SIGHT IS CLOSE TO PARALLEL
C TO THE 3-SPACE W AXIS, THE 3-SPACE V AXIS IS CHOSEN (IN-
C STEAD OF THE 3-SPACE W AXIS) TO BE TRANSFORMED INTO THE
C 2-SPACE Y AXIS.
C
      ASSIGN  90 TO JDONE
      IF (ISCALE)  10, 30, 10
   10 X0 = XMIN
      Y0 = YMIN
      X1 = NLX
      Y1 = NBY
      X2 = NRX-NLX
      Y2 = NTY-NBY
      X3 = X2/(XMAX-XMIN)
      Y3 = Y2/(YMAX-YMIN)
      X4 = NRX
      Y4 = NTY
      FACT = 1.
      IF (BIGD .LE. 0.) GO TO  20
      X0 = -BIGD
      Y0 = -BIGD
      X3 = X2/(2.*BIGD)
      Y3 = Y2/(2.*BIGD)
      FACT = R0/D
   20 DELCRT = X2
      ASSIGN  80 TO JDONE
   30 IF (SINGA .LT. 0.0001) GO TO  40
      R = 1./SINGA
      ASSIGN  70 TO JUMP
      RETURN
   40 SINBE = SIN(BE)
      R = 1./SINBE
      ASSIGN  60 TO JUMP
      RETURN
C
C********************  ENTRY TRN32  ************************
C     ENTRY TRN32 (U,V,W,XT,YT,ZT)
C
   50 UU = U
      VV = V
      WW = W
      Q = D/((UU-EU)*COSAL+(VV-EV)*COSBE+(WW-EW)*COSGA)
      GO TO JUMP,( 60, 70)
   60 UU = ((EW+Q*(WW-EW)-AW)*COSAL-(EU+Q*(UU-EU)-AU)*COSGA)*R
      VV = (EV+Q*(VV-EV)-AV)*R
      GO TO JDONE,( 80, 90)
   70 UU = ((EU+Q*(UU-EU)-AU)*COSBE-(EV+Q*(VV-EV)-AV)*COSAL)*R
      VV = (EW+Q*(WW-EW)-AW)*R
      GO TO JDONE,( 80, 90)
   80 XT = AMIN1(X4,AMAX1(X1,X1+X3*(FACT*UU-X0)))
      YT = AMIN1(Y4,AMAX1(Y1,Y1+Y3*(FACT*VV-Y0)))
      RETURN
   90 XT = UU
      YT = VV
      RETURN
      END
      SUBROUTINE ZEROSC
      SAVE
C
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
C
C ZERO BOTH SCRENE MODELS.
C
      DO  20 I=1,LX
         DO  10 J=1,NY
            ISCR(I,J) = 0
            ISCA(I,J) = 0
   10    CONTINUE
   20 CONTINUE
      RETURN
      END
      SUBROUTINE STCNTR (Z,L,M,N,CONV)
C
      SAVE
      DIMENSION       Z(L,N)
C
C THIS ROUTINE FINDS THE BEGINNINGS OF ALL CONTOUR LINES AT LEVEL CONV.
C FIRST THE EDGES ARE SEARCHED FOR LINES INTERSECTING THE EDGE (OPEN
C LINES) THEN THE INTERIOR IS SEARCHED FOR LINES WHICH DO NOT INTERSECT
C THE EDGE (CLOSED LINES).  BEGINNINGS ARE STORED IN IR TO PREVENT RE-
C TRACING OF LINES.  IF IR IS FILLED, THE SEARCH IS STOPPED FOR THIS
C CONV.
C
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR9/ BIG        ,IXBIT
C
C  PACK X AND Y
C
      IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
C
      IENTRY = 0
      NP = 0
      CV = CONV
C
C THE FOLLOWING CODE SHOULD BE RE-ENABLED IF THIS ROUTINE IS USED FOR
C GENERAL CONTOURING
C
C     ISS=0
C     DO 2 IP1=2,M
C     I=IP1-1
C     IF(Z(I,1).GE.CV.OR.Z(IP1,1).LT.CV) GO TO 1
C     IX=IP1
C     IY=1
C     IDX=-1
C     IDY=0
C     IS=1
C     CALL DRLINE(Z,L,M,N)
C   1 IF(Z(IP1,N).GE.CV.OR.Z(I,N).LT.CV) GO TO 2
C     IX=I
C     IY=N
C     IDX=1
C     IDY=0
C     IS=5
C     CALL DRLINE(Z,L,M,N)
C   2 CONTINUE
C     DO 4 JP1=2,N
C     J=JP1-1
C     IF(Z(M,J).GE.CV.OR.Z(M,JP1).LT.CV) GO TO 3
C     IX=M
C     IY=JP1
C     IDX=0
C     IDY=-1
C     IS=7
C     CALL DRLINE(Z,L,M,N)
C   3 IF(Z(1,JP1).GE.CV.OR.Z(1,J).LT.CV) GO TO 4
C     IX=1
C     IY=J
C     IDX=0
C     IDY=1
C     IS=3
C     CALL DRLINE(Z,L,M,N)
C   4 CONTINUE
C
      ISS = 1
      DO  40 JP1=3,N
         J = JP1-1
         DO  30 IP1=2,M
            I = IP1-1
            IF (Z(I,J).GE.CV .OR. Z(IP1,J).LT.CV) GO TO  30
            IXY = IPXY(IP1,J)
            IF (NP .EQ. 0) GO TO  20
            DO  10 K=1,NP
               IF (IR(K) .EQ. IXY) GO TO  30
   10       CONTINUE
   20       NP = NP+1
            IF (NP .GT. NR) RETURN
            IR(NP) = IXY
            IX = IP1
            IY = J
            IDX = -1
            IDY = 0
            IS = 1
            CALL DRCNTR (Z,L,M,N)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE DRCNTR (Z,L,MM,NN)
      SAVE
C
      DIMENSION       Z(L,NN)
C
C THIS ROUTINE TRACES A CONTOUR LINE WHEN GIVEN THE BEGINNING BY STLINE.
C TRANSFORMATIONS CAN BE ADDED BY DELETING THE STATEMENT FUNCTIONS FOR
C FX AND FY IN DRLINE AND MINMAX AND ADDING EXTERNAL FUNCTIONS.
C X=1. AT Z(1,J), X=FLOAT(M) AT Z(M,J). X TAKES ON NON-INTEGER VALUES.
C Y=1. AT Z(I,1), Y=FLOAT(N) AT Z(I,N). Y TAKES ON NON-INTEGER VALUES.
C
      COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
     1                IS         ,ISS        ,NP         ,CV         ,
     2                INX(8)     ,INY(8)     ,IR(500)    ,NR
      COMMON /ISOSR9/ BIG        ,IXBIT
C
      LOGICAL         IPEN       ,IPENO
C
      DATA IOFFP,SPVAL/0,0./
      DATA IPEN,IPENO/.TRUE.,.TRUE./
C
C  PACK X AND Y
C
      IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
      FX(X1,Y1) = X1
      FY(X1,Y1) = Y1
      C(P11,P21) = (P11-CV)/(P11-P21)
C
      M = MM
      N = NN
      IF (IOFFP .EQ. 0) GO TO  10
      ASSIGN 100 TO JUMP1
      ASSIGN 150 TO JUMP2
      GO TO  20
   10 ASSIGN 120 TO JUMP1
      ASSIGN 160 TO JUMP2
   20 IX0 = IX
      IY0 = IY
      IS0 = IS
      IF (IOFFP .EQ. 0) GO TO  30
      IX2 = IX+INX(IS)
      IY2 = IY+INY(IS)
      IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL
      IPENO = IPEN
   30 IF (IDX .EQ. 0) GO TO  40
      Y = IY
      ISUB = IX+IDX
      X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
      GO TO  50
   40 X = IX
      ISUB = IY+IDY
      Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
   50 IF (IPEN) CALL FRSTS (FX(X,Y),FY(X,Y),1)
   60 IS = IS+1
      IF (IS .GT. 8) IS = IS-8
      IDX = INX(IS)
      IDY = INY(IS)
      IX2 = IX+IDX
      IY2 = IY+IDY
      IF (ISS .NE. 0) GO TO  70
      IF (IX2.GT.M .OR. IY2.GT.N .OR. IX2.LT.1 .OR. IY2.LT.1) GO TO 190
   70 IF (CV-Z(IX2,IY2))  80, 80, 90
   80 IS = IS+4
      IX = IX2
      IY = IY2
      GO TO  60
   90 IF (IS/2*2 .EQ. IS) GO TO  60
      GO TO JUMP1,(100,120)
  100 ISBIG = IS+(8-IS)/6*8
      IX3 = IX+INX(ISBIG-1)
      IY3 = IY+INY(ISBIG-1)
      IX4 = IX+INX(ISBIG-2)
      IY4 = IY+INY(ISBIG-2)
      IPENO = IPEN
      IF (ISS .NE. 0) GO TO 110
      IF (IX3.GT.M .OR. IY3.GT.N .OR. IX3.LT.1 .OR. IY3.LT.1) GO TO 190
      IF (IX4.GT.M .OR. IY4.GT.N .OR. IX4.LT.1 .OR. IY4.LT.1) GO TO 190
  110 IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL .AND.
     1       Z(IX3,IY3).NE.SPVAL .AND. Z(IX4,IY4).NE.SPVAL
  120 IF (IDX .EQ. 0) GO TO 130
      Y = IY
      ISUB = IX+IDX
      X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
      GO TO 140
  130 X = IX
      ISUB = IY+IDY
      Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
  140 GO TO JUMP2,(150,160)
  150 IF (.NOT.IPEN) GO TO 170
      IF (IPENO) GO TO 160
C
C END OF LINE SEGMENT
C
      CALL FRSTS (D1,D2,3)
      CALL FRSTS (FX(XOLD,YOLD),FY(XOLD,YOLD),1)
C
C CONTINUE LINE SEGMENT
C
  160 CALL FRSTS (FX(X,Y),FY(X,Y),2)
  170 XOLD = X
      YOLD = Y
      IF (IS .NE. 1) GO TO 180
      NP = NP+1
      IF (NP .GT. NR) GO TO 190
      IR(NP) = IPXY(IX,IY)
  180 IF (ISS .EQ. 0) GO TO  60
      IF (IX.NE.IX0 .OR. IY.NE.IY0 .OR. IS.NE.IS0) GO TO  60
C
C END OF LINE
C
  190 CALL FRSTS (D1,D2,3)
      RETURN
      END
      SUBROUTINE TR32 (X,Y,MX,MY)
      SAVE
C
      COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
C
C A.S.F. FOR SCALING
C
      SU(UTEMP) = UTEMP
      SV(VTEMP) = VTEMP
      SW(WTEMP) = WTEMP
C
      XX = X
      YY = Y
      IF (ISLBT)  10, 20, 30
   10 CALL TRN32I (SU(U),SV(XX-1.),SW(YY-1.),XT,YT,DUM,2)
      GO TO  40
   20 CALL TRN32I (SU(XX-1.),SV(V),SW(YY-1.),XT,YT,DUM,2)
      GO TO  40
   30 CALL TRN32I (SU(XX-1.),SV(YY-1.),SW(W),XT,YT,DUM,2)
   40 MX = XT
      MY = YT
      RETURN
      END
      SUBROUTINE FRSTS (XX,YY,IENT)
C
C THIS IS A SPECIAL VERSION OF THE SMOOTHING DASHED LINE PACKAGE.  LINES
C ARE SMOOTHED IN THE SAME WAY, BUT NO SOFTFARE DASHED LINES ARE USED.
C CONDITIONAL PLOTTING ROUTINES ARE CALL WHICH DETERMINE THE VISIBILITY
C OF A LINE SEGMENT BEFORE PLOTTING.
C
      SAVE
      DIMENSION       XSAVE(70)  ,YSAVE(70)  ,XP(70)     ,YP(70)     ,
     1                TEMP(70)
C
      COMMON /ISOSR7/ IENTRY     ,IONES
C
      DATA NP/150/
      DATA L1/70/
      DATA TENSN/2.5/
      DATA PI/3.14159265358/
      DATA SMALL/128./
C
      AVE(A,B) = .5*(A+B)
C
C  DECIDE IF FRSTS,VECTS,LASTS CALL
C
      GO TO ( 10, 20, 40),IENT
   10 DEG = 180./PI
      X = XX
      Y = YY
      LASTFL = 0
      SSLP1 = 0.0
      SSLPN = 0.0
      XSVN = 0.0
      YSVN = 0.0
C
C INITIALIZE THE POINT AND SEGMENT COUNTER
C N COUNTS THE NUMBER OF POINTS/SEGMENT
C
      N = 0
C
C NSEG = 0       FIRST SEGMENT
C NSEG = 1       MORE THAN ONE SEGMENT
C
      NSEG = 0
      CALL TR32 (X,Y,MX,MY)
C
C SAVE THE X,Y COORDINATES OF THE FIRST POINT
C XSV1           CONTAINS THE X COORDINATE OF THE FIRST POINT
C                OF A LINE
C YSV1           CONTAINS THE Y COORDINATE OF THE FIRST POINT
C                OF A LINE
C
      XSV1 = MX
      YSV1 = MY
      GO TO  30
C
C ************************* ENTRY VECTS *************************
C     ENTRY VECTS (XX,YY)
C
   20 X = XX
      Y = YY
C
C VECTS          SAVES THE X,Y COORDINATES OF THE ACCEPTED
C                POINTS ON A LINE SEGMENT
C
      CALL TR32 (X,Y,MX,MY)
C
CIF THE NEW POINT IS TOO CLOSE TO THE PREVIOUS POINT, IGNORE IT
C
      IF (ABS(FLOAT(IFIX(XSVN)-MX))+ABS(FLOAT(IFIX(YSVN)-MY)) .LT.
     1    SMALL) RETURN
      IFLAG = 0
   30 N = N+1
C
C SAVE THE X,Y COORDINATES OF EACH POINT OF THE SEGMENT
C XSAVE          THE ARRAY OF X COORDINATES OF LINE SEGMENT
C YSAVE          THE ARRAY OF Y COORDINATES OF LINE SEGMENT
C
      XSAVE(N) = MX
      YSAVE(N) = MY
      XSVN = XSAVE(N)
      YSVN = YSAVE(N)
      IF (N .GE. L1-1) GO TO  50
      RETURN
C
C ************************* ENTRY LASTS *************************
C     ENTRY LASTS
C
   40 LASTFL = 1
C
C LASTS          CHECKS FOR PERIODIC LINES AND SETS UP
C                  THE CALLS TO KURV1S AND KURV2S
C
C IFLAG = 0      OK TO CALL LASTS DIRECTLY
C IFLAG = 1      LASTS WAS JUST CALLED FROM BY VECTS
C                IGNORE CALL TO LASTS
C
      IF (IFLAG .EQ. 1) RETURN
C
C COMPARE THE LAST POINT OF SEGMENT WITH FIRST POINT OF LINE
C
   50 IFLAG = 1
C
C IPRD = 0       PERIODIC LINE
C IPRD = 1       NON-PERIODIC LINE
C
      IPRD = 1
      IF (ABS(XSV1-XSVN)+ABS(YSV1-YSVN) .LT. SMALL) IPRD = 0
C
C TAKE CARE OF THE CASE OF ONLY TWO DISTINCT P0INTS ON A LINE
C
      IF (NSEG .GE. 1) GO TO  70
      IF (N-2) 160,150, 60
   60 IF (N .GE. 4) GO TO  70
      DX = XSAVE(2)-XSAVE(1)
      DY = YSAVE(2)-YSAVE(1)
      SLOPE = ATAN2(DY,DX)*DEG+90.
      IF (SLOPE .GE. 360.) SLOPE = SLOPE-360.
      IF (SLOPE .LE. 0.) SLOPE = SLOPE+360.
      SLP1 = SLOPE
      SLPN = SLOPE
      ISLPSW = 0
      SIGMA = TENSN
      GO TO 110
   70 SIGMA = TENSN
      IF (IPRD .GE. 1) GO TO  90
      IF (NSEG .GE. 1) GO TO  80
C
C SET UP FLAGS FOR A  1  SEGMENT, PERIODIC LINE
C
      ISLPSW = 4
      XSAVE(N) = XSV1
      YSAVE(N) = YSV1
      GO TO 110
C
C SET UP FLAGS FOR AN N-SEGMENT, PERIODIC LINE
C
   80 SLP1 = SSLPN
      SLPN = SSLP1
      ISLPSW = 0
      GO TO 110
   90 IF (NSEG .GE. 1) GO TO 100
C
C SET UP FLAGS FOR THE 1ST SEGMENT OF A NON-PERIODIC LINE
C
      ISLPSW = 3
      GO TO 110
C
C SET UP FLAGS FOR THE NTH SEGMENT OF A NON-PERIODIC LINE
C
  100 SLP1 = SSLPN
      ISLPSW = 1
C
C CALL THE SMOOTHING ROUTINES
C
  110 CALL KURV1S (N,XSAVE,YSAVE,SLP1,SLPN,XP,YP,TEMP,S,SIGMA,ISLPSW)
      IF (IPRD.EQ.0 .AND. NSEG.EQ.0 .AND. S.LT.70.) GO TO 170
      IENTRY = 1
C
C DETERMINE THE NUMBER OF POINTS TO INTERPOLATE FOR EACH SEGMENT
C
      IF (NSEG.GE.1 .AND. N.LT.L1-1) GO TO 120
      NPRIME = FLOAT(NP)-(S*FLOAT(NP))/(2.*32768.)
      IF (S .GE. 32768.) NPRIME = .5*FLOAT(NP)
      NPL = FLOAT(NPRIME)*S/32768.
      IF (NPL .LT. 2) NPL = 2
  120 DT = 1./FLOAT(NPL)
      IF (NSEG .LE. 0) CALL FRSTC (IFIX(XSAVE(1)),IFIX(YSAVE(1)),1)
      T = 0.0
      NSLPSW = 1
      IF (NSEG .GE. 1) NSLPSW = 0
      NSEG = 1
      CALL KURV2S (T,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
C
C SAVE SLOPE AT THE FIRST POINT OF THE LINE
C
      IF (NSLPSW .GE. 1) SSLP1 = SLP
      NSLPSW = 0
      XSOLD = XSAVE(1)
      YSOLD = YSAVE(1)
      DO 130 I=1,NPL
         T = T+DT
         TT = -T
         IF (I .EQ. NPL) NSLPSW = 1
         CALL KURV2S (TT,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
C
C SAVE THE LAST SLOPE OF THIS LINE SEGMENT
C
         IF (NSLPSW .GE. 1) SSLPN = SLP
C
C DRAW EACH PART OF THE LINE SEGMENT
C
         CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
         CALL FRSTC (IFIX(XS),IFIX(YS),2)
         XSOLD = XS
         YSOLD = YS
  130 CONTINUE
      IF (IPRD .NE. 0) GO TO 140
C
C CONNECT THE LAST POINT WITH THE FIRST POINT OF A PERIODIC LINE
C
      CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
      CALL FRSTC (IFIX(XSV1),IFIX(YSV1),2)
C
C BEGIN THE NEXT LINE SEGMENT WITH THE LAST POINT OF THIS SEGMENT
C
  140 XSAVE(1) = XS
      YSAVE(1) = YS
      N = 1
  150 CONTINUE
  160 RETURN
  170 N = 0
      RETURN
      END
      SUBROUTINE FRSTC (MX,MY,IENT)
      SAVE
C
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR4/ RX         ,RY
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
     1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
     2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
     3                IHB        ,IHS        ,IHV        ,IVOLD      ,
     4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
     5                IHF
      LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
C
C
C  DRAW LINE TO THE POINT MX,MY
C
C  ENTER THE POINT INTO THE CURRENT SCREEN, ISCR, IF THE POINT CONFORMS
C  TO THE SHADING ALGORITHM.
C  THE POINT IS NOT ENTERED WHEN;
C  1. IT IS THE SAME POINT USED IN THE LAST CALL, RESOLUTION PROBLEM
C  2. IT IS PART OF A HORIZONTAL LINE BUT NOT AN END POINT
C  3.  THE ENTIRE CONTOUR RESTS ON A HORIZONTAL PLANE
C
C  WHEN DRAWING A HORIZONTAL LINE THREE CONDITIONS EXIST;
C  1. WHEN THE LINE IS A HORIZONTAL STEP ENTER ONLY THE OUTSIDE POINT.
C     A HORIZONTAL STEP IS DEFINED BY THE ENTERING AND EXITING Y
C     DIRECTION THAT IS THE SAME.
C  2. ENTER BOTH END POINTS OF A HORIZONTAL TURNING POINT. A HORIZONTAL
C     TURNING POINT IS A LINE WITH GREATER THAN 1 HORIZONTAL BITS
C     AND THE ENTERING AND EXITING Y DIRECTION IS DIFFIRENT.
C  3.  WHEN THE ENTIRE CONTOUR IS A HORIZONTAL LINE NO POINTS ARE
C      ENTERED. THIS CONDITION IS DETECTED BY THE STATUS OF YCHANG.
C      IF IT IS TRUE THEN THE CONTOUR IS NOT A SINGLE HORIZONTAL LINE.
C
C  THE PREVIOUS POINT IS ERASED IF IT IS A VERTICAL TURNING POINT.
C  A VERTICAL TURNING POINT IS A HORIZONTAL LINE WITH ONLY 1 POINT
C  AND THE ENTERING AND EXITING Y DIRECTION DIFFERS.THIS DATA IS
C  IN THE VARIABLES IOSLSN-OLD SLOPE AND ISLSGN-NEW SLOPE.
C  THE CHANGE IN SLOPE MUST BE -1 TO 1 OR 1 TO -1.
C
C  OTHERWISE THE POINT IS ENTERED INTO ISCR.
C
C  THE TWO ENTRY POINTS ARE REQUIRED BY THE HARDWARE DRAWING ROUTINES.
C  FIRSTC IS USED FOR THE FIRST POINT ON THE CONTOUR. THE REMAINING
C  POINTS ON THE SAME CONTOUR ARE ENTERED VIA VECTC.
C
      DATA IONE/1/
      AVE(A,B) = (A+B)*.5
C
C  COMPUTE VISIBILITY OF THIS POINT
C
C  WARNING
C  IF X OR Y PLOTTER MAXIMUM VALUE RANGES FALL BELOW 101 THEN THE
C  FOLLOWING TWO STATEMENTS WHICH SET IX AND IY MUST BE CHANGED.
C  REPLACE THE CONSTANT 1.0 BY 0.5 IN THE STATEMENTS WHERE THE
C  MAXIMUM PLOTTER VALUE IS LESS THAN 101 FOR THAT DIRECTION.  THE
C  PLOTTER CORDINATE RANGES ARE SET IN SET32.
C
      IX = FLOAT(MX-1)*RX+1.0
      NRLX = IX
      IY = FLOAT(MY-1)*RY+1.0
      IBIT = NBPW-MOD(IX,NBPW)
      IX = IX/NBPW+1
      IVNOW = IAND(ISHIFT(ISCA(IX,IY),1-IBIT),IONE)
C
C  DECIDE IF FRSTC OR VECTC CALL
C
      IF (IENT .NE. 1) GO TO  10
C
      XOLD = MX
      YOLD = MY
C
C
C  SET INITIAL VALUES
C
      IHF = .FALSE.
      IYDIR = 0
      ITPD = 0
      IVAL = 0
      IOSLSN = 0
      IFSX = NRLX
      IFSY = IY
      LASTV = IVNOW
      HBFLAG = .FALSE.
      YCHANG = .FALSE.
      CALL PLOTIT (IFIX(XOLD),IFIX(YOLD),0)
      GO TO 180
C
C****************************  ENTRY VECTC  ****************************
C     ENTRY VECTC (MX,MY)
C
   10 XNOW = MX
      YNOW = MY
      JUMP = IVNOW*2+LASTV+1
      GO TO ( 20, 30, 40, 50),JUMP
C
C BOTH VISIBLE
C
   20 CALL PLOTIT (IFIX(XNOW),IFIX(YNOW),1)
      GO TO  50
C
C JUST TURNED VISIBLE
C
   30 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),0)
      GO TO  50
C
C JUST TURNED INVISIBLE
C
   40 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),1)
C
C BOTH INVISIBLE
C
   50 XOLD = XNOW
      YOLD = YNOW
      LASTV = IVNOW
C
C  TEST FOR RESOLUTION PROBLEM
C
      IF (NRLX.EQ.LRLX .AND. IY.EQ.IYOLD) RETURN
C
C  TEST FOR HORIZONTAL BITS
C
      IF (IYOLD .NE. IY) GO TO  70
C
C  HORIZONTAL BITS DETECTED. SET FLAG AND EXIT.
C  THIS AND THE NEXT HORIZONTAL BIT TEST IS NECESSARY FOR ISCR TO
C  CONFORM TO THE SHADING ALGORITHM IN SUBROUTINE FILLIN
C
C
C  IF HORIZONTAL LINE PREVIOUSLY DETECTED EXIT
C
      IF (.NOT.HBFLAG) GO TO  60
C
C  IF END OF CONTOUR ON A HORIZONTAL LINE BRANCH FOR SPECIAL PROCESSING.
C
      IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
      GO TO 200
C
C  SAVE SLOPE PRIOR TO HORIZONTAL LINE
C
   60 IHX = IXOLD
      IHB = IBTOLD
      IHS = IOSLSN
      IOSLSN = 0
      HBFLAG = .TRUE.
      IHRX = LRLX
      IHV = IVOLD
      IF (LRLX.EQ.IFSX .AND. IYOLD.EQ.IFSY) IHF = .TRUE.
C
C  THIS IS THE SECOND TRAP FOR END OF CONTOUR ON A HORIZONTAL LINE.
C
      IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
      GO TO 200
C
C  COMPUTE THE SLOPE TO THIS POINT
C
   70 IF (IY-IYOLD)  80, 90,100
   80 ISLSGN = 1
      GO TO 110
   90 ISLSGN = 0
      GO TO 120
  100 ISLSGN = -1
  110 IF (IYDIR .EQ. 0) IYDIR = ISLSGN
  120 CONTINUE
C
C  IF PROCESS REACHES THIS CODE THE CONTOUR IS NOT CONTAINED ON A SINGLE
C  HORIZONTAL  PLANE, SO RECORD THIS FACT BY SETTING Y CHANGE FLAG.
C
      YCHANG = .TRUE.
C
C  TEST FOR END OF HORIZONTAL LINE
C
      IF (.NOT.HBFLAG) GO TO 160
      HBFLAG = .FALSE.
C
C  HORIZONTAL LINE JUST ENDED
C
C  TEST FOR REDRAW
C
      ITEMP = IAND(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
      IF ((IHV .EQ. 0) .AND. (ITEMP .EQ. 0)) GO TO 130
C
C  REDRAWING ERASE THIS POINT
C
      ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
      ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
      GO TO 170
C
C  TEST FOR STEP PROBLEM
C
  130 IF (IHS .NE. ISLSGN) GO TO 140
C
C  STEP PROBLEM
C
      GO TO 170
C
C  TURNING PROBLEM HORIZONTAL LINE IS A TURNING POINT
C
  140 CONTINUE
C
C  ENTER THE TURNING POINT ONLY IF IT IS NOT THE SECOND SUCCEEDING
C  EVENT IN A ROW
C
      ICTPD = 1
      IF (IHRX .GT. NRLX) ICTPD = -1
      IF (ICTPD .NE. ITPD) GO TO 150
      ITPD = 0
C
C  ERASE THE FIRST POINT
C
      ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
      GO TO 170
C
C  ENTER THE TURNING POINT
C
  150 CONTINUE
      ITPD = ICTPD
C
C  ENTER THE SECOND POINT
C
      ISCR(IXOLD,IYOLD) = IOR(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
      GO TO 170
C
C  CHECK IF PREVIOUS ENTRY WAS A VERTICAL TURNING POINT.
C  IF SO ERASE IT.
C
  160 IF (ISLSGN.EQ.IOSLSN .OR. (IOSLSN.EQ.0 .OR. ISLSGN.EQ.0))
     1    GO TO 170
      ITPD = 0
      ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
C
  170 IOSLSN = ISLSGN
C
C  CHECK IF THIS GRID POINT PREVIOUSLY ACTIVATED
C
      IVAL = IAND(ISCR(IX,IY),MASK(IBIT))
C
C  IF GRID POINTS ACTIVATED BRANCH
C
      IF (IVAL .NE. 0) GO TO 190
C
C  GRID POINT NOT ACTIVATED   SET AND EXIT
C
  180 CONTINUE
      ISCR(IX,IY) = IOR(ISCR(IX,IY),MASK(IBIT))
      GO TO 200
C
C  THIS POINT IS BEING REDRAWN SO ERASE IT.
C  (THIS IS TO CONFORM WITH THE SHADING ALGORITHM, FILLIN.
C  HOWEVER IF BACK TO STARTING POINT DO NOT ERASE
C
  190 IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) RETURN
      ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
C
C
  200 IXOLD = IX
      LRLX = NRLX
      IYOLD = IY
      IBTOLD = IBIT
      IVOLD = IVAL
      RETURN
C
C  PERFORM THIS OPERATION WHEN A CONTOUR STARTS OR ENDS ON A HORIZONTAL
C  LINE.
C
  210 CONTINUE
C
C  ERASE THE FIRST POINT OF A CONTOUR WHEN IT IS PART OF A HORIZONTAL
C  LINE SEGMENT AND IS NOT THE ENDPOINT OF THE SEGMENT
C
      IF (.NOT.IHF) GO TO 220
      ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
  220 CONTINUE
C
C  ERASE THE FIRST POINT OF A HORIZONTAL LINE SEGMENT WHEN IT ENDS
C  THE CONTOUR AND IS NOT THE HIGHEST LINE SEG ON THS SIDE.
C
      IF (.NOT.YCHANG) GO TO 230
      IF (IYDIR .NE. IHS) GO TO 200
  230 ISCR(IHX,IY) = IAND(ISCR(IHX,IY),NMASK(IHB))
      GO TO 200
      END
      SUBROUTINE FILLIN
C
      SAVE
      COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
     1                ISCA(8,128)
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR7/ IENTRY     ,IONES
C
      IF (IENTRY .EQ. 0) RETURN
C
C  THIS IS A SHADING ALGORITHM IT IS USED TO DETERMINE CONTOUR LINES
C  THAT ARE HIDDEN BY THE PRESENT LINE. THE ALGORITHM PROCESSES
C  HORIZONTAL ROWS. IT ASSUMES THAT THE BIT PATTERN PASSED TO IT
C  HAS ONLY BITS SET TO MARK THE START AND END OF SHADING. THE
C  ALGORITHM ALSO ASSUMES THAT WHEN AN ON BIT IS ENCOUNTERED THAT A
C  CORRESPONDING OFF BIT IS INCLUDED IN THE SAME ROW.
C
C
C  PULL OUT ROWS OF THE CONTOUR PATTERN
C
      IBVAL = 0
      DO  80 IYNOW=1,NY
         DO  40 IXNOW=1,LX
C
C  IF NO ACTIVATED BITS BRANCH
C
            ICRWD = ISCR(IXNOW,IYNOW)
            IF (ICRWD .EQ. 0) GO TO  30
C
C  ACTIVATED BITS IN WORD SET SHADING FLAG
C
C  CHECK BIT BY BIT FOR ON/OFF FLAGS
C
            DO  20 IB=1,NBPW
               IBIT = (NBPW+1)-IB
C
C
C  PULL OUT THE CURRENT GRID POINT VALUE
C
               IVAL = IAND(ICRWD,MASK(IBIT))
C
C  IF IVAL SET, THIS IS AN ON/OFF FLAG
C
               IF (IVAL .EQ. 0) GO TO  10
C
C  FLAG BIT, ALWAYS SET
C
               IBVAL = MOD(IBVAL+1,2)
               GO TO  20
C
C  SHADE THE SCREEN ACCORDING TO THE STATUS OF IBVAL
C
   10          IF (IBVAL .NE. 0) ICRWD = IOR(ICRWD,MASK(IBIT))
C
   20       CONTINUE
C
C  ZERO OUT THE SCREEN
C
            ISCR(IXNOW,IYNOW) = 0
            ISCA(IXNOW,IYNOW) = IOR(ICRWD,ISCA(IXNOW,IYNOW))
            GO TO  40
C
   30       IF (IBVAL .NE. 0) ISCA(IXNOW,IYNOW) = IONES
   40    CONTINUE
C
C  FIX FOR NONCORRECTABLE RUNAWAYS
C
         IF (IBVAL .EQ. 0) GO TO  80
         IBVAL = 0
         DO  70 K=1,LX
            ITEST = 0
            IF (IYNOW .EQ. 1) GO TO  50
            ITEST = ISCA(K,IYNOW-1)
            IF (IYNOW .EQ. NY) GO TO  60
   50       ITEST = IOR(ITEST,ISCA(K,IYNOW+1))
   60       ISCA(K,IYNOW) = ITEST
   70    CONTINUE
C
   80 CONTINUE
      RETURN
      END
      SUBROUTINE DRAWI (IXA,IYA,IXB,IYB)
C
C INCLUDED FOR USE BY PWRZ
C
      SAVE
      CALL FRSTC (IXA,IYA,1)
      CALL FRSTC (IXB,IYB,2)
      RETURN
      END
      SUBROUTINE MMASK
C
C  MAKE THE MACHINE DEPENDENT MASKS USED IN THE CONTOUR DRAWING
C  AND SHADING ALGORITHMS
C
      SAVE
      COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
      LOGICAL         GENDON
      COMMON /ISOSR7/ IENTRY     ,IONES
      COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
     1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
     2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
     3                IHB        ,IHS        ,IHV        ,IVOLD      ,
     4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
     5                IHF
      COMMON /ISOSR9/ BIG        ,IXBIT
      LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
      GENDON = .TRUE.
      NBPW = 16
C
C  GET BIGGEST REAL NUMBER
C
      BIG = R1MACH(2)
C
C  MASKS TO SELECT A SPECIFIC BIT
C
      DO  10 K=1,NBPW
         MASK(K) = ISHIFT(1,K-1)
   10 CONTINUE
C
C  GENERATE  THE BIT PATTERN 177777 OCTAL
C
      ITEMP1 = 0
      ITEMP = MASK(NBPW)
      IST = NBPW-1
      DO  20 K=1,IST
         ITEMP1 = IOR(ITEMP,ISHIFT(ITEMP1,-1))
   20 CONTINUE
      MFIX = IOR(ITEMP1,1)
C
C  MASKS TO CLEAR A SPECIFIC BIT
C
      DO  30 K=1,NBPW
         NMASK(K) = IAND(ITEMP1,MFIX)
         ITEMP1 = IOR(ISHIFT(ITEMP1,1),1)
   30 CONTINUE
      IONES = MFIX
      RETURN
C
C REVISION HISTORY---
C
C JANUARY 1978     DELETED REFERENCES TO THE  *COSY  CARDS AND
C                  ADDED REVISION HISTORY
C JANUARY 1979     NEW SHADING ALGORITHM
C MARCH 1979       MADE CODE MACHINE INDEPENDENT AND CONFORM
C                  TO 66 FORTRAN STANDARD
C JUNE 1979        THIS VERSION PLACED ON ULIB.
C SEPTEMBER 1979   FIXED PROBLEM IN EZISOS DEALING WITH
C                  DETERMINATION OF VISIBILITY OF W PLANE.
C DECEMBER 1979    FIXED PROBLEM WITH PEN DOWN ON CONTOUR
C                  INITIALIZATION IN SUBROUTINE FRSTC
C MARCH            CHANGED ROUTINE NAMES  TRN32I  AND  DRAW  TO
C                  TRN32I  AND  DRAWI  TO BE CONSISTENT WITH THE
C                  USAGE OF THE NEW ROUTINE  PWRZI.
C JUNE  1980       FIXED PROBLEM WITH ZERO INDEX COMPUTATION IN
C                  SUBROUTINE FRSTC.  ADDED INPUT PARAMETER
C                  DIMENSION STATEMENT MISSING IN EZISOS.
C                  FIXED ERROR IN COMPUTATION OF ARCCOSINE
C                  IN EZISOS AND TRN32I.
C DECEMBER 1984    CONVERTED TO GKS LEVEL 0A AND STANDARD FORTRAN 77
C-----------------------------------------------------------------------
C
      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
c**********************************************************************
      SUBROUTINE FFT2(C,MM,INV)
c**********************************************************************
!...............................................................................
! Standard Fast Fourier Transform routine from Cooley et al 1969 (reproduced in
! Aki & Richards Quantitative Seismology p596)
!      the inverse part of this from Bob Parker
! MM  - There are 2**MM points on the trace being transformed. This trace is
!       passed to the routine as complex array C
!...............................................................................
!    used by tx2tp and tp2tx
!    this is the same as fftfwd from Oppenheim and Shaefer
!
      COMPLEX C(*),U,W,T,CI
      SAVE PI
      DATA PI/3.1415927/
!
      N=2**MM
      NBY2=N/2
      NM1=N-1
      XN = FLOAT(N)
      J=1
!
!      DO 7 I=1,NM1
      DO I=1,NM1
         IF(I.GE.J) GO TO 5
         T=C(J)
         C(J)=C(I)
         C(I)=T
   5     K=NBY2
   6     IF(K.GE.J) GO TO 7
         J=J-K
         K=K/2
         GO TO 6
   7     J=J+K
      ENDDO
!
!      DO 20 L=1,MM
      DO L=1,MM
         LE=2**L
         LE1=LE/2
         U=(1.0,0.0)
         W=CMPLX(COS(PI/FLOAT(LE1)),SIN(PI/FLOAT(LE1)))
         DO J=1,LE1   ! was do 20
!           DO 10 I=J,N,LE
            DO I=J,N,LE
               IP=I+LE1
               T=C(IP)*U
               C(IP)=C(I)-T
  10           C(I)=C(I)+T
            ENDDO
  20        U=U*W
         ENDDO
         U=U*W
      ENDDO
!
      IF(INV.NE.0) THEN
        N1=1 + N/2
!        DO 30 I = 2,N1
        DO I = 2,N1
           CI=C(I)
           C(I)=C(N-I+2)/XN
  30       C(N-I+2) = CI/XN
        ENDDO
        C(1)=C(1)/XN
      ENDIF
      RETURN
      END
c*********************************************************************
      subroutine realtr(a,b,n,isn)
c*********************************************************************
c              rearranges array for fft of a real series in place
c              called after a forward fft, before an inverse fft
c              example: forward transform of real series, length n
c              call fft (x,n/2)
c              call realtr(x,x(2),n/2,2)
c                   here fft expects complex array with real and imag
c                   parts alternating. the result appears as cos,sin,
c                   cos,sin,... in x(1),x(2),...,x(n+1),x(n+2)
c                   thus array x must be dimensioned at least n+2 in
c                   order to hold the transform of the Nyquist
c             from Bob Parker
c             now with proper scaling!
c           **CALLS NO OTHER SUBROUTINES**
c
      dimension a(1),b(1)
      real im
      inc = iabs(isn)
      nk = n * inc + 2
      nh = nk/2
      sd = 2. * atan(1.)/float(n)
      cd = 2. * sin(sd)**2
      sd = sin(sd+sd)
      sn = 0.
      IF(ISN.GT.0) go to 10
      cn = -1.
      sd = -sd
      go to 20
   10 cn = 1.
      a(nk-1) = a(1)
      b(nk-1) = b(1)
   20 do 30 j=1,nh,inc
      k = nk - j
      aa = a(j) + a(k)
      ab = a(j) - a(k)
      ba = b(j) + b(k)
      bb = b(j) - b(k)
      re = cn*ba + sn*ab
      im = sn*ba - cn*ab
      b(k) = (im - bb)/2.
      b(j) = (im + bb)/2.
      a(k) = (aa - re)/2.
      a(j) = (aa + re)/2.
      aa = cn - (cd*cn + sd*sn)
      sn = (sd*cn - cd*sn) + sn
      cn = aa
   30 continue
      return
      end


text:
SUBROUTINE set_borderf(set_border,phi_sun, zmax, zmin, sun_zen, hx, hy, &
    az_case, d, d0, k_max, h_offset, n_inc, m_inc, n_add, m_add, &
    k_setting, add_max, ierr)

    implicit none

!   set_border defines the buffer that is needed to find the
!   occluding terrain. This will be larger for low sun elevations

    integer k_max, k, n_add, m_add, k_setting, add_max, ierr, az_case
    real n_inc(k_setting), m_inc(k_setting)
    real h_offset(k_setting)
    real phi_sun, zmax, zmin, sun_zen, phc
    real*8 hx, hy
    real sinphc, cosphc, d, d0
    real pi, r2d, d2r, eps
    common/base/pi,r2d,d2r,eps
    logical set_border

    set_border=.true.
    ierr=0

!   the case is dependent on the sun azimuth
!   this defines which borders need to be available
!   NOTE: azimuth must be in degrees between 0 and 360 and clockwise from N

    if (phi_sun.ge.0.0 .and. phi_sun.le.90.0) then
        az_case=1
        phc=pi/2.0-phi_sun*d2r
    else if (phi_sun.ge.270.0 .and. phi_sun.le.360.0) then
        az_case=2
        phc=phi_sun*d2r-3.0*pi/2.0
    else if (phi_sun.ge.180.0 .and. phi_sun.le.270.0) then
        az_case=3
        phc=3.0*pi/2.0-phi_sun*d2r
    else if (phi_sun.ge.90.0 .and. phi_sun.le.180.0) then
        az_case=4
        phc=phi_sun*d2r-pi/2.0
    else
        ierr=61
        goto 99
    endif

!   calculate the border info
!   and the increments for the projection line

    sinphc=sin(phc)
    cosphc=cos(phc)

!   d is the plane distance from the pixel
!   with zmin where if the pixel at the
!   distance in the sun direction is at Zmax
!   then it can just occlude the first pixel
!
    d = (zmax-zmin)/tan(pi/2.0-sun_zen*d2r)
    if(cosphc*hy.gt.sinphc*hx) then
        d0 = 0.5*hx/cosphc
    else
        d0 = 0.5*hy/sinphc
    endif

!   d0 is the basic step in the d direction that increments
!   0.5 pixel in the larger of the two directions

!   k_max is the number of steps to get out to d
!   n_inc(k) and m_inc(k) are the sample and line increments
!   h_offset(k) on exit is the altitude above the start of the
!   sun beam as you move towards the sun

    k_max=ifix(d/d0+0.5)+1
    if (k_max.gt.k_setting) then
        ierr=62
        goto 99
    endif
    do k=1,k_max
        h_offset(k)=float(k)*d0

!       the increments to search for occlusions have sign depending
!       on the azimuth case

        if (az_case.eq.1.or.az_case.eq.4) then
            n_inc(k)=h_offset(k)*cosphc/hx
        else
            n_inc(k)=-h_offset(k)*cosphc/hx
        endif
        if (az_case.eq.3.or.az_case.eq.4) then
            m_inc(k)=h_offset(k)*sinphc/hy
        else
            m_inc(k)=-h_offset(k)*sinphc/hy
        endif
        h_offset(k)=h_offset(k)*tan(pi/2.0-sun_zen*d2r)
    enddo

!   n_add and m_add are the sizes of the pixel buffer
!   in sample and line dimensions

!   the two buffers will be on sides of the target image defined
!   by the azimuth case

    n_add=ifix(sngl(d*cosphc/hx+1.5))
    m_add=ifix(sngl(d*sinphc/hy+1.5))

    if ((n_add.gt.add_max .or. m_add.gt.add_max) .or. (n_add.lt.0.or.m_add.lt.0)) then
        ierr=63
        goto 99
    endif
    return
99  set_border=.false.
    return
END SUBROUTINE set_borderf

max_stars_repo_name: Oceancolour-RG/wagl

max_issues_repo_name: Oceancolour-RG/wagl

max_forks_repo_name: Oceancolour-RG/wagl


text:
      Subroutine Type2601

! Object: WBGT calculation for inside
! Simulation Studio Model: Type2601 WBGT inside
!

! Author: Yuichi Yasuda
! Editor:
! Date:	 June 22, 2019
! last modified: June 22, 2019
!
!
! ***
! *** Model Parameters
! ***
!			Nb. of WBGT	- [1;20]

! ***
! *** Model Inputs
! ***
!			Wet Bulb temperature	C [-Inf;+Inf]
!			Globe or Operative temperature	C [-Inf;+Inf]

! ***
! *** Model Outputs
! ***
!			WBGT	C [-Inf;+Inf]

! ***
! *** Model Derivatives
! ***

! (Comments and routine interface generated by TRNSYS Simulation Studio)
!************************************************************************

!-----------------------------------------------------------------------------------------------------------------------
! This TRNSYS component skeleton was generated from the TRNSYS studio based on the user-supplied parameters, inputs,
! outputs, and derivatives.  The user should check the component formulation carefully and add the content to transform
! the parameters, inputs and derivatives into outputs.  Remember, outputs should be the average value over the timestep
! and not the value at the end of the timestep; although in many models these are exactly the same values.  Refer to
! existing types for examples of using advanced features inside the model (Formats, Labels etc.)
!-----------------------------------------------------------------------------------------------------------------------

      Use TrnsysConstants
      Use TrnsysFunctions

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

!DEC$Attributes DLLexport :: Type2601

!-----------------------------------------------------------------------------------------------------------------------
!Trnsys Declarations
      Implicit None

      Double Precision Timestep,Time
      Integer CurrentUnit,CurrentType
      Integer i
      Integer, Parameter :: MaxWBGT=20  !maximum number of WBGTs per component

!    PARAMETERS
      INTEGER nWBGT

!    INPUTS
      DOUBLE PRECISION wetbulbTemp(MaxWBGT) !湿球温度[C]
      DOUBLE PRECISION globeTemp(MaxWBGT)   !黒球温度（室内では作用温度）[C]

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

!-----------------------------------------------------------------------------------------------------------------------
!Get the Global Trnsys Simulation Variables
      Time=getSimulationTime()
      Timestep=getSimulationTimeStep()
      CurrentUnit = getCurrentUnit()
      CurrentType = getCurrentType()
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Set the Version Number for This Type
      If(getIsVersionSigningTime()) Then
		Call SetTypeVersion(17)
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Do Any Last Call Manipulations Here
      If(getIsLastCallofSimulation()) Then
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Perform Any "After Convergence" Manipulations That May Be Required at the End of Each Timestep
      If(getIsEndOfTimestep()) Then
		Return
      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Do All of the "Very First Call of the Simulation Manipulations" Here
      If(getIsFirstCallofSimulation()) Then

        !get the param1, number of WBGT.
        nWBGT = JFIX(getParameterValue(1)+0.1)
        if(nWBGT > MaxWBGT) then
            Call FoundBadParameter(1,'Fatal','The number of the WBGT to be calculated must be between 1 and 20.')
        endif

		!Tell the TRNSYS Engine How This Type Works
		Call SetNumberofParameters(1)           !The number of parameters that the the model wants
		Call SetNumberofInputs(2*nWBGT)         !The number of inputs that the the model wants
		Call SetNumberofDerivatives(0)          !The number of derivatives that the the model wants
		Call SetNumberofOutputs(nWBGT)          !The number of outputs that the the model produces
		Call SetIterationMode(1)                !An indicator for the iteration mode (default=1).  Refer to section 8.4.3.5 of the documentation for more details.
		Call SetNumberStoredVariables(0,0)      !The number of static variables that the model wants stored in the global storage array and the number of dynamic variables that the model wants stored in the global storage array
		Call SetNumberofDiscreteControls(0)     !The number of discrete control functions set by this model (a value greater than zero requires the user to use Solver 1: Powell's method)

        !Set the Correct Input and Output Variable Types
        do i=1, GetNumberofInputs()
            Call SetInputUnits(i,'TE1') !湿球温度、黒球温度（室内では作用温度）[C]
        end do

        do i=1, GetNumberofOutputs()
            Call SetOutputUnits(i,'TE1') !WBGT[C]
        end do

		Return

      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!Do All of the First Timestep Manipulations Here - There Are No Iterations at the Intial Time
      If (getIsStartTime()) Then
        nWBGT = JFIX(getParameterValue(1)+0.1)

        do i=1, nWBGT
            wetBulbTemp(i)  = GetInputValue(i*2-1)
            globeTemp(i)    = GetInputValue(i*2)
        end do

   !Check the Parameters for Problems (#,ErrorType,Text)
   !Sample Code: If( PAR1 <= 0.) Call FoundBadParameter(1,'Fatal','The first parameter provided to this model is not acceptable.')

   !Set the Initial Values of the Outputs (#,Value)
      do i=1, nWBGT
		!Call SetOutputValue(1, 0.0) ! WBGT
        Call SetOutputValue(i, 0.7*wetBulbTemp(i)+0.3*globeTemp(i) ) ! WBGT
      end do

   !If Needed, Set the Initial Values of the Static Storage Variables (#,Value)
   !Sample Code: SetStaticArrayValue(1,0.d0)

   !If Needed, Set the Initial Values of the Dynamic Storage Variables (#,Value)
   !Sample Code: Call SetDynamicArrayValueThisIteration(1,20.d0)

   !If Needed, Set the Initial Values of the Discrete Controllers (#,Value)
   !Sample Code for Controller 1 Set to Off: Call SetDesiredDiscreteControlState(1,0)

		Return

      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!ReRead the Parameters if Another Unit of This Type Has Been Called Last
      If(getIsReReadParameters()) Then
		!Read in the Values of the Parameters from the Input File
        nWBGT = JFIX(getParameterValue(1)+0.1)
        Call SetNumberofInputs(2*nWBGT)

      EndIf
!-----------------------------------------------------------------------------------------------------------------------

!Read the Inputs
      do i=1, nWBGT
        wetBulbTemp(i)  = GetInputValue(i*2-1)
        globeTemp(i)    = GetInputValue(i*2)
      end do

	!Check the Inputs for Problems (#,ErrorType,Text)
	!Sample Code: If( IN1 <= 0.) Call FoundBadInput(1,'Fatal','The first input provided to this model is not acceptable.')

      If(ErrorFound()) Return
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!    *** PERFORM ALL THE CALCULATION HERE FOR THIS MODEL. ***
!-----------------------------------------------------------------------------------------------------------------------

	!-----------------------------------------------------------------------------------------------------------------------
	!If Needed, Get the Previous Control States if Discrete Controllers are Being Used (#)
	!Sample Code: CONTROL_LAST=getPreviousControlState(1)
	!-----------------------------------------------------------------------------------------------------------------------

	!-----------------------------------------------------------------------------------------------------------------------
	!If Needed, Get the Values from the Global Storage Array for the Static Variables (#)
	!Sample Code: STATIC1=getStaticArrayValue(1)
	!-----------------------------------------------------------------------------------------------------------------------

	!-----------------------------------------------------------------------------------------------------------------------
	!If Needed, Get the Initial Values of the Dynamic Variables from the Global Storage Array (#)
	!Sample Code: T_INITIAL_1=getDynamicArrayValueLastTimestep(1)
	!-----------------------------------------------------------------------------------------------------------------------

	!-----------------------------------------------------------------------------------------------------------------------
	!Perform All of the Calculations Here to Set the Outputs from the Model Based on the Inputs

	!Sample Code: OUT1=IN1+PAR1

	!If the model requires the solution of numerical derivatives, set these derivatives and get the current solution
	!Sample Code: T1=getNumericalSolution(1)
	!Sample Code: T2=getNumericalSolution(2)
	!Sample Code: DTDT1=3.*T2+7.*T1-15.
	!Sample Code: DTDT2=-2.*T1+11.*T2+21.
	!Sample Code: Call SetNumericalDerivative(1,DTDT1)
	!Sample Code: Call SetNumericalDerivative(2,DTDT2)

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

!-----------------------------------------------------------------------------------------------------------------------
!Set the Outputs from this Model (#,Value)
        do i=1, nWBGT
		    Call SetOutputValue(i, 0.7*wetBulbTemp(i)+0.3*globeTemp(i) ) ! WBGT
        end do
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!If Needed, Store the Desired Disceret Control Signal Values for this Iteration (#,State)
!Sample Code:  Call SetDesiredDiscreteControlState(1,1)
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!If Needed, Store the Final value of the Dynamic Variables in the Global Storage Array (#,Value)
!Sample Code:  Call SetDynamicArrayValueThisIteration(1,T_FINAL_1)
!-----------------------------------------------------------------------------------------------------------------------

      Return
      End
!-----------------------------------------------------------------------------------------------------------------------

max_stars_repo_path: SourceCode/Type2601 WBGT inside.f90

max_stars_repo_name: yuizi/ComfortIndex

max_issues_repo_path: SourceCode/Type2601 WBGT inside.f90

max_issues_repo_name: yuizi/ComfortIndex

max_forks_repo_path: SourceCode/Type2601 WBGT inside.f90

max_forks_repo_name: TRNSYSJP/ComfortIndex


text:
function one_arg(a)
  implicit none
  !Automatic type discovery: the C type correspondacnes
  !will be automatically determined at compile time
  integer :: one_arg
  real(kind=8), dimension(:, : ), intent(inout) :: a

  a = 1.61803399_8
  one_arg = 42

end function

max_stars_repo_name: howland/fFwrapEx

max_issues_repo_name: howland/fFwrapEx

max_forks_repo_name: howland/fFwrapEx


text:
subroutine get_face_velocity(level, time, &
     vx, vx_l1, vx_l2, vx_l3, vx_h1, vx_h2, vx_h3, &
     vy, vy_l1, vy_l2, vy_l3, vy_h1, vy_h2, vy_h3, &
     vz, vz_l1, vz_l2, vz_l3, vz_h1, vz_h2, vz_h3, &
     dx, prob_lo) bind(C, name="get_face_velocity")

  use probdata_module, only : adv_vel

  implicit none

  integer, intent(in) :: level
  double precision, intent(in) :: time
  integer, intent(in) :: vx_l1, vx_l2, vx_l3, vx_h1, vx_h2, vx_h3
  integer, intent(in) :: vy_l1, vy_l2, vy_l3, vy_h1, vy_h2, vy_h3
  integer, intent(in) :: vz_l1, vz_l2, vz_l3, vz_h1, vz_h2, vz_h3
  double precision, intent(out) :: vx(vx_l1:vx_h1,vx_l2:vx_h2,vx_l3:vx_h3)
  double precision, intent(out) :: vy(vy_l1:vy_h1,vy_l2:vy_h2,vy_l3:vy_h3)
  double precision, intent(out) :: vz(vz_l1:vz_h1,vz_l2:vz_h2,vz_l3:vz_h3)
  double precision, intent(in) :: dx(3), prob_lo(3)

  vx = adv_vel(1)
  vy = adv_vel(2)
  vz = adv_vel(3)

end subroutine get_face_velocity

max_stars_repo_name: ylunalin/amrex

max_issues_repo_name: ylunalin/amrex

max_forks_repo_name: ylunalin/amrex


text:
! implements atmosphere model of Potekhin, Chabrier, and Yakovlev (1997).

module pcy97

    contains

    subroutine do_get_pcy97_Teff(grav, Plight, lgTb, lgTeff, lgflux)
        use constants_def
        real(dp), intent(in) :: grav    ! surface gravity, in the local frame
        real(dp), intent(in) :: Plight  ! pressure at which layer of light elements terminates
        real(dp), intent(in), dimension(:) :: lgTb  ! temperature at a base column
        real(dp), intent(out), dimension(:) :: lgTeff, lgflux   ! effective temperature and flux
        real(dp) :: eta, g14
        real(dp), dimension(size(lgTb)) :: Tb9, Teff6_4

        eta = Plight/1.193e34_dp
        g14 = grav * 1.0e-14_dp
        Tb9 = 10.0_dp**lgTb*1.0e-9_dp
        call PCYfit(g14,eta,Tb9,Teff6_4)
        lgTeff = log10(1.0e6_dp*Teff6_4**0.25_dp)
        lgflux = log10(sigma_SB*Teff6_4*1.0e24_dp)

        contains
        subroutine PCYfit(g14,eta,Tb9,Teff6_4)
            real(dp), intent(in) :: g14,eta
            real(dp), dimension(:), intent(in) :: Tb9
            real(dp), dimension(size(Tb9)), intent(out) :: Teff6_4  ! (Teff/10^6 K)^4
            real(dp), dimension(size(Tb9)) :: a,Tstar,zeta, Te4Fe, Te4a

            Tstar = sqrt(7.0_dp*Tb9*sqrt(g14))
            zeta = Tb9-1.0e-3_dp*Tstar
            Te4Fe = g14*((7.0*zeta)**2.25_dp + (onethird*zeta)**1.25_dp)

            if (eta > 0) then
                Te4a = g14*(18.1_dp*Tb9)**2.42_dp
                a = tb9**fivethird * (1.2_dp + (5.3e-6_dp/eta)**0.38_dp)
                Teff6_4 = (a*Te4Fe + Te4a)/(a + 1.0_dp)
            else
                Teff6_4 = Te4Fe
            end if
        end subroutine PCYfit

    end subroutine do_get_pcy97_Teff

end module pcy97

max_stars_repo_name: LBJ-Wade/dStar

max_issues_repo_name: LBJ-Wade/dStar

max_forks_repo_name: LBJ-Wade/dStar


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

#ifdef ARTED_USE_NVTX
#define NVTX_BEG(name,id)  call nvtxStartRange(name,id)
#define NVTX_END()         call nvtxEndRange()
#else
#define NVTX_BEG(name,id)
#define NVTX_END()
#endif

subroutine hamiltonian(zu,flag_current)
  use Global_Variables
  use timer
  use omp_lib
  use opt_variables
  use hpsi, only: hpsi_omp_KB_RT
  implicit none
  integer    :: tid
  integer    :: ikb,ik,ib,i
  integer    :: iexp
  complex(8) :: zfac(4)
  complex(8), intent(inout) :: zu(NL,NBoccmax,NK_s:NK_e)
  logical, intent(in)       :: flag_current

#ifndef ARTED_CURRENT_PREPROCESSING
#define UNUSED_VARIABLE(VAR) if(.false.) call salmon_unusedvar(VAR)
  UNUSED_VARIABLE(flag_current)
#endif

  zfac(1)=(-zI*dt)
  do i=2,4
    zfac(i)=zfac(i-1)*(-zI*dt)/i
  end do

  call timer_begin(LOG_HPSI)

!$omp parallel private(tid) shared(zfac)
!$  tid=omp_get_thread_num()

!$omp do private(ik,ib,iexp)
  do ikb=1,NKB
    ik=ik_table(ikb)
    ib=ib_table(ikb)

    call init(zhtpsi(:,4,tid),zu(:,ib,ik))
    call hpsi_omp_KB_RT(ik,zhtpsi(:,4,tid),zhtpsi(:,1,tid))
    call hpsi_omp_KB_RT(ik,zhtpsi(:,1,tid),zhtpsi(:,2,tid))
    call hpsi_omp_KB_RT(ik,zhtpsi(:,2,tid),zhtpsi(:,3,tid))
    call hpsi_omp_KB_RT(ik,zhtpsi(:,3,tid),zhtpsi(:,4,tid))
    call update(zfac,zhtpsi(:,:,tid),zu(:,ib,ik))

#ifdef ARTED_CURRENT_PREPROCESSING
    if(flag_current) call current_omp_KB_ST(ib,ik,zu(:,ib,ik))
#endif
  end do
!$omp end do
!$omp end parallel

  call timer_end(LOG_HPSI)

contains
  subroutine init(tpsi,zu)
    use Global_Variables, only: NLx,NLy,NLz
    use opt_variables, only: PNLx,PNLy,PNLz
    use timer
    implicit none
    complex(8) :: tpsi(0:PNLz-1,0:PNLy-1,0:PNLx-1)
    complex(8) :: zu(0:NLz-1,0:NLy-1,0:NLx-1)
    integer :: ix,iy,iz

    TIMER_BEG(LOG_HPSI_INIT)
!dir$ vector aligned
    do ix=0,NLx-1
    do iy=0,NLy-1
    do iz=0,NLz-1
      tpsi(iz,iy,ix)=zu(iz,iy,ix)
    end do
    end do
    end do
    TIMER_END(LOG_HPSI_INIT)
  end subroutine

  subroutine update(zfac,tpsi,zu)
    use Global_Variables, only: NLx,NLy,NLz
    use opt_variables, only: PNLx,PNLy,PNLz
    use timer
    implicit none
    complex(8) :: zfac(4)
    complex(8) :: tpsi(0:PNLz-1,0:PNLy-1,0:PNLx-1,4)
    complex(8) :: zu(0:NLz-1,0:NLy-1,0:NLx-1)
    integer :: ix,iy,iz

    TIMER_BEG(LOG_HPSI_UPDATE)
!dir$ vector aligned
    do ix=0,NLx-1
    do iy=0,NLy-1
    do iz=0,NLz-1
      zu(iz,iy,ix)=zu(iz,iy,ix)+zfac(1)*tpsi(iz,iy,ix,1) &
      &                        +zfac(2)*tpsi(iz,iy,ix,2) &
      &                        +zfac(3)*tpsi(iz,iy,ix,3) &
      &                        +zfac(4)*tpsi(iz,iy,ix,4)
    end do
    end do
    end do
    TIMER_END(LOG_HPSI_UPDATE)
  end subroutine
end subroutine

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON