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

C $Id: mxmult.f,v 1.3 1993/07/16 22:50:49 gdsjaar Exp $
C $Log: mxmult.f,v $
C Revision 1.3  1993/07/16 22:50:49  gdsjaar
C Unrolled loops for faster execution
C
c Revision 1.2  1993/07/16  19:14:19  gdsjaar
c Restructured to optimize faster
c
c Revision 1.1  1993/07/16  16:47:35  gdsjaar
c Changed plt to library rather than single source file.
c
C=======================================================================
      SUBROUTINE MXMULT(N,MAT1,MAT2,MATR)
      REAL MAT1(N,*),MAT2(N,*),MATR(N,*)

      IF (N .EQ. 4) THEN
          MATR(1,1) = MAT1(1,1)*MAT2(1,1) + MAT1(1,2)*MAT2(2,1) +
     *                MAT1(1,3)*MAT2(3,1) + MAT1(1,4)*MAT2(4,1)
          MATR(1,2) = MAT1(1,1)*MAT2(1,2) + MAT1(1,2)*MAT2(2,2) +
     *                MAT1(1,3)*MAT2(3,2) + MAT1(1,4)*MAT2(4,2)
          MATR(1,3) = MAT1(1,1)*MAT2(1,3) + MAT1(1,2)*MAT2(2,3) +
     *                MAT1(1,3)*MAT2(3,3) + MAT1(1,4)*MAT2(4,3)
          MATR(1,4) = MAT1(1,1)*MAT2(1,4) + MAT1(1,2)*MAT2(2,4) +
     *                MAT1(1,3)*MAT2(3,4) + MAT1(1,4)*MAT2(4,4)

          MATR(2,1) = MAT1(2,1)*MAT2(1,1) + MAT1(2,2)*MAT2(2,1) +
     *                MAT1(2,3)*MAT2(3,1) + MAT1(2,4)*MAT2(4,1)
          MATR(2,2) = MAT1(2,1)*MAT2(1,2) + MAT1(2,2)*MAT2(2,2) +
     *                MAT1(2,3)*MAT2(3,2) + MAT1(2,4)*MAT2(4,2)
          MATR(2,3) = MAT1(2,1)*MAT2(1,3) + MAT1(2,2)*MAT2(2,3) +
     *                MAT1(2,3)*MAT2(3,3) + MAT1(2,4)*MAT2(4,3)
          MATR(2,4) = MAT1(2,1)*MAT2(1,4) + MAT1(2,2)*MAT2(2,4) +
     *                MAT1(2,3)*MAT2(3,4) + MAT1(2,4)*MAT2(4,4)

          MATR(3,1) = MAT1(3,1)*MAT2(1,1) + MAT1(3,2)*MAT2(2,1) +
     *                MAT1(3,3)*MAT2(3,1) + MAT1(3,4)*MAT2(4,1)
          MATR(3,2) = MAT1(3,1)*MAT2(1,2) + MAT1(3,2)*MAT2(2,2) +
     *                MAT1(3,3)*MAT2(3,2) + MAT1(3,4)*MAT2(4,2)
          MATR(3,3) = MAT1(3,1)*MAT2(1,3) + MAT1(3,2)*MAT2(2,3) +
     *                MAT1(3,3)*MAT2(3,3) + MAT1(3,4)*MAT2(4,3)
          MATR(3,4) = MAT1(3,1)*MAT2(1,4) + MAT1(3,2)*MAT2(2,4) +
     *                MAT1(3,3)*MAT2(3,4) + MAT1(3,4)*MAT2(4,4)

          MATR(4,1) = MAT1(4,1)*MAT2(1,1) + MAT1(4,2)*MAT2(2,1) +
     *                MAT1(4,3)*MAT2(3,1) + MAT1(4,4)*MAT2(4,1)
          MATR(4,2) = MAT1(4,1)*MAT2(1,2) + MAT1(4,2)*MAT2(2,2) +
     *                MAT1(4,3)*MAT2(3,2) + MAT1(4,4)*MAT2(4,2)
          MATR(4,3) = MAT1(4,1)*MAT2(1,3) + MAT1(4,2)*MAT2(2,3) +
     *                MAT1(4,3)*MAT2(3,3) + MAT1(4,4)*MAT2(4,3)
          MATR(4,4) = MAT1(4,1)*MAT2(1,4) + MAT1(4,2)*MAT2(2,4) +
     *                MAT1(4,3)*MAT2(3,4) + MAT1(4,4)*MAT2(4,4)

      ELSE

        DO 230 K = 1,N
          DO 200 J = 1,N
            MATR(K,J) = 0.0
 200      CONTINUE
          DO 220 I = 1,N
            DO 210 J = 1,N
              MATR(K,J) = MATR(K,J) + MAT1(K,I)*MAT2(I,J)
 210        CONTINUE
 220      CONTINUE
 230    CONTINUE
      END IF
      RETURN
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
      INTEGER FUNCTION vmec_chdir(new_path)
      IMPLICIT NONE
      CHARACTER*(*), INTENT(in) :: new_path
#if defined(CRAY)
      INTEGER :: ilen
      iLEN = 0
      CALL pxfchdir(new_path, ilen, vmec_chdir)
#else
      INTEGER, EXTERNAL :: chdir

#if defined(SUNOS)
      vmec_chdir = chdir(TRIM(new_path))
#else
      vmec_chdir = chdir(TRIM(new_path) // char(0))
#endif
#endif
      END FUNCTION vmec_chdir


text:
Program testCR

Use FileIO_Utilities, Only: creturn,wait
Implicit None
Integer:: i

Write(*,*)
Do i = 1,10
    Write(*,'(A,I0,A)',ADVANCE='NO') 'Testing...',i,creturn
    Call Wait(500)
End Do
Write(*,'(A)') 'Did it work?          '
Write(*,*)
Write(*,'(A)') 'Press a key to exit...'
Read(*,*)

End Program

max_stars_repo_name: wtdailey/commonTools

max_issues_repo_name: wtdailey/commonTools

max_forks_repo_name: wtdailey/commonTools


text:
program pgm
integer :: a(2:3,2:3,2:3), i , j, k, c

c = 1
do i = 2, 3
do j = 2, 3
do k = 2, 3
 a(j,i,k) = c
 c = c + 1
enddo
enddo
enddo

do k = 2, 3
do j = 2, 3
do i = 2, 3
 if (a(k,j,i) <= a(i,j,k)) then
  print *, a(k,j,i)
 endif
enddo
enddo
enddo
end

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
Davis 1 Stop is a Gas Stations gas station. The Gas Station is called Davis 1 Stop because they have everything from a Great selection of Fine Wines and Beers, to Great prices on other Liquor items.
  Staff is very friendly, helpful, and some assert that they have the Best Touchless CarWash in Town.
   They carry full selection of Cigr. and Tobacco products .
  They will match or beat any Diesel Price in Davis.
  They will also match the prices on Beer, wine , liquor , cigr.,Tobacco at other Gas Stations or conv./Liquor stores.

20110207 19:25:47 nbsp Got my car washed here this afternoon. I liked the idea of a touchless car wash that uses only soap and highpressured water with no brushes, carpet strips, chamois curtains, or dangling cloth, etc. touching the car, therefore avoiding scratches and damage to the paint.

However, the $8 basic wash that I got didnt clean my car so well and my car wasnt even really dirty to begin with. It just had a few spots of bird poop that didnt come off, nor did the bugs stuck on the windshield and front fender. I ended up washing the car again at home and probably wont be taking my car for a wash there again unless I do a prewash first.  Users/pvo

20111019 04:39:30 nbsp A sign posted at the entrance offered a 79 cent cup of coffee with a fillup 8 gallons or more. So I fill up my cars tank with 14 gallons of gas and go back to get my cup of coffee. Guess what? I get charged $1.80 for the coffee. So I tell the cashier about the sign in the window and he tells me that I have to ask for the discount first before buying the coffee. He does NOT offer to refund me the difference! Unbelievable. Also, their restroom was beyond filthy, was out of toilet paper, and the door was broken and didnt close completely, allowing other customers to walk in. Users/MarkWChang

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE
     . SURVRELH
     .(xALPHA,POWER,xM1,xA,xF,xN,xM,ALT,M2L,M2H,RL,RH)
!MS$ATTRIBUTES REFERENCE :: ALT
!MS$ATTRIBUTES REFERENCE :: M2L
!MS$ATTRIBUTES REFERENCE :: M2H
!MS$ATTRIBUTES REFERENCE :: RL
!MS$ATTRIBUTES REFERENCE :: RH
	IMPLICIT NONE
C
C-Description-----------------------------------------------------------
C
C  Function:
C	Determine the relative hazard for survival data.
C
C  Input (provided by user):
C	ALPHA	Type I error probability (two tailed)
C	POWER	The desired statistical power
C	M1	Median survival time on control treatment
C	A	The accrual time during which patients are recruited into
C		the study
C	F	Additional follow-up time after the end of patient
C		recruitment
C	N	Number of experimental patients
C	M	Ratio of control to experimental patients
C  Output:
C     ALT=1  both low and high solutions found
C     ALT=2  only low solution found
C     ALT=3  only high solution found
C	M2L,M2H	The median survival times on the experimental treatment
C		that can be detected with power POWER given N experimental
C		patients, M controls per experimental patient, and
C		type I error probability ALPHA
C	RL = M2L/M1
C	RH = M2H/M1
C
C  Reference:
C	Schoenfeld DA, Richter JR: "Nomograms for Calculating the
C	Number of Patients Needed for a Clinical Trial With Survival
C	as an Endpoint" Biometrics 1982; 38: 163-170.
C
C  Notes:
C   .	This routine was written by Dale Plummer
C   .	Designed by William Dupont
C   .	Schoenfeld and Richter define BETA to equal the study power
C	    rather than the type II error.  Otherwise the notation used
C	    in this program is analogous to that used by these authors.
C
C-Declarations----------------------------------------------------------
C
C  Arguments
C
	REAL xALPHA,POWER,xM1,xA,xF,xN,xM,M2L,M2H,RL,RH
	INTEGER	ALT
C
C  Locals
C
	REAL E1,EPS,Y1,Y2, xBETA, FY
	INTEGER ERRORL,ERRORH
	INTEGER MAXFN
C
C    Declare function to be solved for power.
C
	EXTERNAL SURVFUNC1
	REAL SURVFUNC1
C
C    Data passed in common to functions.
C
	REAL ALPHA,BETA,M1,A,F,N,M,ZALPHA,ZBETA
	COMMON /SURVCOM1/ ALPHA,BETA,M1,A,F,N,M
	COMMON /SURVCOM2/ ZALPHA, ZBETA
C
C-Code------------------------------------------------------------------
C
C    BETA is 1-power.
C
	xBETA=1.-POWER
C
C    Load the common data area.
C
	ALPHA=xALPHA
	BETA=xBETA
	M1=xM1
	A=xA
	F=xF
	N=xN
	M=xM
C
C    Here is an insufficient attempt to give A a reasonable value relative to F.
C
	IF (A.EQ.0.) A=0.00004*F
C
C    The evaluation of the function will halt if the value of the function
C    is less than EPS.
C
	EPS=0.001*MIN(1.,M1)
C
C    The evaluation of the function will halt if two successive
C    approximations differ by less than E1.
C
	E1=0.001*MIN(1.,M1)
C
C    Set initial end points for lower end point.
C
	Y1=EPS
	Y2=M1-EPS
C
C    Allow up to MAXFN function evaluations.
C
	MAXFN=100
C
C    Solve the equation for the lower solution.
C
	CALL BISEC(Y1,Y2,EPS,E1,SURVFUNC1,M2L,ERRORL)
C
C    Calculate the relative risk.
C
	RL=M2L/M1
C
C    Set the interval in which we will search for the upper M2.  This
C    interval is:  just_over_M1 ... some_large_value.
C
      Y1=M1+EPS
C
C    Set initial end points for high end point.
C
C	Y1=M1+EPS
C	Y2=MAX(10000.,1000.*M1)
C
C    06-Oct-1997 : a new way to set the upper end of this upper interval:
C    look in the range F+A to F+A+(100*M1) for when the function goes negative.
C    Make this point the end of the upper interval in which to search for the
C    upper zero.  If no such point is found the assume that there is no
C    upper solution.
C
	DO Y2=F+A , F+A+(100.*M1)
	    FY=SURVFUNC1(Y2)
	    IF (FY.LT.0) GOTO 110
	END DO
C
C    If we fall out of this loop then we will not try to
C    find an upper solution.  In this case force the state
C    where we indicate that no upper solution was found.
C
	ERRORH=1
	GO TO 1112
C
C    Solve the equation for the upper solution.
C
110	CONTINUE
	CALL BISEC(Y1,Y2,EPS,E1,SURVFUNC1,M2H,ERRORH)
C
C    Calculate the relative risk.
C
	RH=M2H/M1
C
C    Report the results.
C
C    Write first to the terminal screen.
C
1112	CONTINUE
	IF (ERRORL.EQ.0 .AND. ERRORH.EQ.0) THEN
C
C    Found both low and high solutions (M2L and M2H, RL and RH).
C
	    ALT=1
	ELSE IF (ERRORL.EQ.0 .AND. ERRORH.EQ.1) THEN
C
C    Found low solution only (M2L, RL).
C
		ALT=2
	ELSE IF (ERRORL.EQ.1 .AND. ERRORH.EQ.0) THEN
C
C    Found high solution only (M2H, RH).
C
		ALT=3
	ELSE
		ALT=-99
	END IF
	RETURN
	END

ext: for

max_stars_repo_name: vubiostat/ps

max_issues_repo_name: vubiostat/ps

max_forks_repo_name: vubiostat/ps


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:47 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE CLSGINI__genmod
          INTERFACE
            SUBROUTINE CLSGINI(KDATIM,KTSPD,KAOMOD,PSLM,KXA,KYA)
              INTEGER(KIND=4) :: KYA
              INTEGER(KIND=4) :: KXA
              INTEGER(KIND=4) :: KDATIM(6)
              INTEGER(KIND=4) :: KTSPD
              INTEGER(KIND=4) :: KAOMOD
              REAL(KIND=8) :: PSLM(KXA,KYA)
            END SUBROUTINE CLSGINI
          END INTERFACE
        END MODULE CLSGINI__genmod


text:
! ###################################################################
! Copyright (c) 2016-2021, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:EMConvertOrientations.f90
!--------------------------------------------------------------------------
!
! PROGRAM: EMConvertOrientations
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief convert a text file with orientations from one representation to another
!
!> @date 01/31/17 MDG 1.0 original
!--------------------------------------------------------------------------
program EMConvertOrientations
!

use io
use local
use files
use rotations
use constants
use Lambert
use quaternions
use povray
use so3
use dictmod
use math
use NameListTypedefs
use NameListHandlers
use ECPmod
use error

IMPLICIT NONE

character(fnlen)                        :: nmldeffile, progname, progdesc
type(ConvertOrientationsNameListType)   :: enl

integer(kind=irg),parameter :: numtypes = 8

real(kind=dbl)          :: rod(4), cu(3), ho(3), qu(4), eu(3), euFZ(3), om(3,3), ax(4), sp(3), rom(9), q
integer(kind=irg)       :: i,j,k, ierr, io_int(3)

character(fnlen)        :: fname, dataname, outnames(numtypes)
character(2)            :: outstring(numtypes), angleformat
integer(kind=irg)       :: FZtype, FZorder, numpoints, startunit
type(dicttype),pointer  :: dict

nmldeffile = 'EMConvertOrientations.nml'
progname = 'EMConvertOrientations.f90'
progdesc = 'Convert orientations from one representation to another'

! print some information
call EMsoft(progname, progdesc)

! deal with the command line arguments, if any
call Interpret_Program_Arguments(nmldeffile,1,(/ 87 /), progname)

! deal with the namelist stuff
call GetConvertOrientationsNameList(nmldeffile,enl)

! generate the symmetry quaternions
nullify(dict)
allocate(dict)
dict%Num_of_init = 3
dict%Num_of_iterations = 30
! get the point group number from the xtal file
dict%pgnum = GetPointGroup(enl%xtalname)
call DI_Init(dict,'WAT')  ! generate the quaternions; 'WAT' does not matter, could also be 'VMF'

! set the FZ type
FZtype = FZtarray(dict%pgnum)
FZorder = FZoarray(dict%pgnum)

io_int = (/ dict%pgnum, FZtype, FZorder /)
call WriteValue('  Point group, type, and order: ',io_int,3)

! input data file of orientations
dataname = trim(enl%anglefile)

outnames(1) = trim(enl%cubochoric)
outnames(2) = trim(enl%homochoric)
outnames(3) = trim(enl%rodrigues)
outnames(4) = trim(enl%stereographic)
outnames(5) = trim(enl%axisangle)
outnames(6) = trim(enl%eulerangles)
outnames(7) = trim(enl%quaternion)
outnames(8) = trim(enl%rotationmatrix)

outstring = (/ 'cu', 'ho', 'ro', 'sp', 'ax', 'eu', 'qu', 'om' /)

! open the input data file and read the orientation format and the number of orientations
fname = trim(EMsoft_getEMdatapathname())//trim(enl%anglefile)
fname = EMsoft_toNativePath(fname)
open(unit=53,file=trim(fname),status='old',action='read')
read (53,*) angleformat
read (53,*) numpoints

! open the necessary output files and print the first two lines
startunit = 20
do i=1,numtypes
  if (trim(outnames(i)).ne.'undefined') then
    fname = trim(EMsoft_getEMdatapathname())//trim(outnames(i))
    fname = EMsoft_toNativePath(fname)
    open(unit=startunit+i,file=trim(fname),status='unknown',action='write')
    write (startunit+i,*) outstring(i)
    write (startunit+i,*) numpoints
    call Message('opening output file '//trim(outnames(i)))
  end if
end do

do j=1,numpoints
! read the input data and convert to Euler angles (in radians)
  if (angleformat.eq.'eu') then
    read (53,*) eu(1:3)
    eu = eu * cPi/180.D0
  end if
  if (angleformat.eq.'ho') then
    read (53,*) ho(1:3)
    eu = ho2eu(ho)
  end if
  if (angleformat.eq.'qu') then
    read (53,*) qu(1:4)
    eu = qu2eu(qu)
  end if
  if (angleformat.eq.'ro') then
    read (53,*) rod(1:4)
    eu = ro2eu(rod)
  end if
  if (angleformat.eq.'cu') then
    read (53,*) cu(1:3)
    eu = cu2eu(cu)
  end if
  if (angleformat.eq.'sp') then
    read (53,*) sp(1:3)
    q = 1.D0/(1.D0+sum(sp**2))
    qu = (/ 1.D0-sum(sp**2),2.D0*sp(1),2.D0*sp(2),2.D0*sp(3)/) * q
    eu = qu2eu(qu)
  end if
  if (angleformat.eq.'om') then
    read (53,*) rom(1:9)
    om(1,1:3) = rom(1:3)
    om(2,1:3) = rom(4:6)
    om(3,1:3) = rom(7:9)
    eu = om2eu(om)
  end if
  if (angleformat.eq.'ax') then
    read (53,*) ax(1:4)
    eu = ax2eu(ax)
  end if

! do we need to reduce this to the fundamental zone (RFZ) ?
  if (enl%reducetoRFZ.eq.1) then
    call ReduceOrientationtoRFZ(eu, dict, FZtype, FZorder, euFZ)
    eu = euFZ
  end if

! and write the resulting transformed orientation to the selected output files
  do i=1,numtypes
    if (trim(outnames(i)).ne.'undefined') then
      select case (i)
      case (1)
        cu = eu2cu(eu)
        write(startunit+i,*) cu(1:3)
      case (2)
        ho = eu2ho(eu)
        write(startunit+i,*) ho(1:3)
      case (3)
        rod = eu2ro(eu)
        write(startunit+i,*) rod(1:4)
      case (4)
        qu = eu2qu(eu)
        sp = (/ qu(2), qu(3), qu(4) /) / (1.0+qu(1))
        write(startunit+i,*) sp(1:3)
      case (5)
        ax = eu2ax(eu)
        write(startunit+i,*) ax(1:4)
      case (6)
        write(startunit+i,*) eu(1:3)*180.0/cPi
      case (7)
        qu = eu2qu(eu)
        write(startunit+i,*) qu(1:4)
      case (8)
        om = eu2om(eu)
        rom(1:3) = om(1,1:3)
        rom(4:6) = om(2,1:3)
        rom(7:9) = om(3,1:3)
        write(startunit+i,*) rom(1:9)
      case default
        call FatalError(' ','unknown orientation representation')
      end select
    end if
  end do
end do

! and close the output files
do i=1,numtypes
  if (trim(outnames(i)).ne.'undefined') then
    close(unit=startunit+i,status='keep')
    call Message(' --> closed output file '//trim(outnames(i)))
  end if
end do

end program EMConvertOrientations

max_stars_repo_name: josephtessmer/EMsoft

max_stars_repo_licenses: Intel Unlicense MIT-CMU BSD-3-Clause

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
      program offload

      implicit none

       integer :: i,j,k
       integer :: iter,count_rate, count_max,count
       integer :: t_start,t_final
       integer, parameter :: nx=8192,ny=nx,max_iter=2000
       double precision, parameter    :: pi=4d0*dtan(1d0)
       real, parameter    :: error=0.001
       double precision               :: max_err,time_s,&
                                         d2fx,d2fy
       double precision, allocatable  :: f(:,:), f_k(:,:)

       allocate(f(nx,ny)); allocate(f_k(nx,ny))

!generate a random vector
      CALL RANDOM_NUMBER(f)

       call system_clock(count_max=count_max, count_rate=count_rate)

       call system_clock(t_start)

        write(*,*)'--sum before loop:', sum(f(:,:))/(nx*ny)

!offload the data to the device and copy the result back to the host
      do iter=1,5

         do j=1,ny
            do i=1,nx

               f_k(i,j) = 2*f(i,j)

             enddo
          enddo

         write(*,*)'--sum:No update', sum(f_k(:,:))/(nx*ny)

!copy from the device to host
          print*,"--Sum: update self",sum(f_k(:,:))/(nx*ny)

          f_k = f_k/2.
!copy from the host to device
       enddo

       write(*,*)'--sum after end data:', sum(f_k(:,:))/(nx*ny)
       call system_clock(t_final)

       time_s = real(t_final - t_start)/real(count_rate)

       print*, '--Time it takes (s)', time_s
       end

       subroutine arrayf_new(nx,ny,f,f_k)

         implicit none

         integer :: i,j,nx,ny
         double precision :: f(nx,ny),f_k(nx,ny)

         do j=1,ny
            do i=1,nx

               f_k(i,j) = 2*f(i,j)

             enddo
          enddo

          end subroutine arrayf_new

max_stars_repo_name: HichamAgueny/GPU-course

max_issues_repo_name: HichamAgueny/GPU-course

max_forks_repo_name: HichamAgueny/GPU-course


subroutine nbody_timestep(position,velocity)
! Adjusts the timestep based on a standard step doubling algorithm
! Takes two half timesteps and compares to the current result

use stardata,only: debug
use embryodata

implicit none

real,dimension(3,nbodies),intent(in) :: position,velocity

integer :: ibody,ix
real :: halfdt,error
real, dimension(3,nbodies) :: testpos1,testvel1,testpos2,testvel2

halfdt = dt_nbody/2

call nbody_integrate(halfdt,pos,vel,testpos1,testvel1)
call nbody_integrate(halfdt,testpos1,testvel1,testpos2,testvel2)

! Compute error between half timestep and full timestep for each particle

error = 0.0
maxerror = -1.0e30

do ibody=1,nbodies

   error = 0.0
   do ix=1,3

      error = error + (position(ix,ibody)-testpos2(ix,ibody))*(position(ix,ibody)-testpos2(ix,ibody))
      error = error + (velocity(ix,ibody)-testvel2(ix,ibody))*(velocity(ix,ibody)-testvel2(ix,ibody))

   enddo

   error = sqrt(error)

   if(error>maxerror) maxerror = error
enddo

! Compute new deltat based on this error level
! If current error below user defined tolerance, then dt is increased
! If current error above tolerance, then dt is decreased

if(maxerror>small) then
   dt_nbody = dt_nbody*abs(tolerance/maxerror)**0.2
endif

if(maxerror>tolerance) then
   if(debug=='y') then
      print*,'reducing timestep'
      print*,tolerance,maxerror,dt_nbody
   endif

   dt_nbody = dt_nbody*0.1
endif

end subroutine nbody_timestep


text:
C->>> ---------------------------------------------> ems_sv_bs_cg_iz <<<
c     Initialise the basis change records.
c
      subroutine ems_sv_bs_cg_iz
      implicit none
      include 'SVBSCG.INC'

      ca_sv_bs_cg_iz_fg1 = 1
      ca_sv_bs_cg_iz_fg2 = 2
      bs_cg_rec = .false.
      return
      end

C->>> ---------------------------------------------> ems_sv_bs_cg_sa <<<
c     (Re-)start saving basis change records.
c
      subroutine ems_sv_bs_cg_sa(ems_rt_cod)
      implicit none
      include 'EMSRTCOD.INC'
      include 'EMSMSG.INC'
      include 'SVBSCG.INC'
      integer ems_rt_cod
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) then
c
c     Initialise the basis change records if basis change records are
c     not currently being saved.
c
         n_sv_bs_cg = 0
         n_sv_bs_cg_seq = 1
         sv_bs_cg_seq_sa(1) = 1
         bs_cg_rec = .true.
      endif
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' starting to save basis change records')
      end

C->>> -------------------------------------------> ems_sv_bs_cg_stop <<<
c     Stop saving basis change records
c
      subroutine ems_sv_bs_cg_stop
      implicit none
      include 'SVBSCG.INC'

      bs_cg_rec = .false.
      return
      end

C->>> ------------------------------------------------> ems_sv_bs_cg <<<
c     Save a basis change record---if records are being saved.
c
      subroutine ems_sv_bs_cg(ems_rt_cod,
     &     en_vr, lv_vr, en_vr_prev_st, lv_vr_prev_st)
      implicit none
      include 'EMSRTCOD.INC'
      include 'EMSMSG.INC'
      include 'SVBSCG.INC'
      integer ems_rt_cod
      integer en_vr, lv_vr, en_vr_prev_st, lv_vr_prev_st
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      integer prev_n_sv_bs_cg_seq
      integer ca_ems_rt_cod
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) goto 7000
c
c     Save the basis change record---may jump back to here if a
c     sequence of records has been removed to create space.
c
 100   continue
      if (n_sv_bs_cg .lt. mx_n_sv_bs_cg) then
c
c     The basis change record can be saved within the current sequence
c
         n_sv_bs_cg = n_sv_bs_cg + 1
         sv_bs_cg_en_vr(n_sv_bs_cg) = en_vr
         sv_bs_cg_lv_vr(n_sv_bs_cg) = lv_vr
         sv_bs_cg_en_vr_prev_st(n_sv_bs_cg) = en_vr_prev_st
         if (en_vr .ne. lv_vr) then
            sv_bs_cg_lv_vr_prev_st(n_sv_bs_cg) = lv_vr_prev_st
         else
            sv_bs_cg_lv_vr_prev_st(n_sv_bs_cg) = en_vr_prev_st
         endif
         goto 7000
      endif
c
c     Try to create room to save the basis change record within the
c     current sequence by removing the oldest sequence.
c
c     NB:
c     1. This may have to be done several times if there are empty
c     .  sequences
c     2. This may be (immediately) unsuccessful if only one sequence
c     .  uses all the available entries---this will eventually occur
c     .  when the number of basis changes between resets exceeds
c     .  mx_n_sv_bs_cg.
c
 200  continue
      if (n_sv_bs_cg_seq .gt. 1) then
c
c     Keep a copy of the number of sequences before removing one so that
c     the (barely conceivable) infinite loop can avoided visibly in this
c     routine.
c
         prev_n_sv_bs_cg_seq = n_sv_bs_cg_seq
         call ems_sv_bs_cg_rm_ol_seq(ca_ems_rt_cod)
         if (ca_ems_rt_cod .ne. ems_rt_cod_ok) then
            ems_rt_cod = max(ca_ems_rt_cod, ems_rt_cod)
            if (ems_rt_cod .ge. ems_rt_lvl_serious) goto 7100
         endif
c
c     Go back and save the basis change record if space is available.
c
         if (n_sv_bs_cg .lt. mx_n_sv_bs_cg) goto 100
c
c     The sequence just removed must have been null so remove the next
c
c     Check that the number of sequences has been reduced
c
         if (n_sv_bs_cg_seq .lt. prev_n_sv_bs_cg_seq) goto 200
      endif
c
c     Stop recording basis change records
c
      call ems_sv_bs_cg_stop
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' saving basis change records')
      end

C->>> --------------------------------------> ems_sv_bs_cg_rm_ol_seq <<<
c     Remove the oldest sequence (if possible).
c
      subroutine ems_sv_bs_cg_rm_ol_seq(ems_rt_cod)
      implicit none
      include 'EMSRTCOD.INC'
      include 'EMSMSG.INC'
      include 'SVBSCG.INC'
      integer ems_rt_cod
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      integer fm_bs_cg_n, t_bs_cg_n, bs_cg_n, seq_n, os
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) goto 7000
c      print*, '### BEFORE REMOVING BASIS CHANGE RECORD SEQUENCE ###'
c      print*, ' n_sv_bs_cg_seq = ', n_sv_bs_cg_seq
c      print*, ' n_sv_bs_cg = ', n_sv_bs_cg
      if (n_sv_bs_cg_seq .gt. 1) then
c         call ems_sv_bs_cg_rp(ems_rt_cod)
         fm_bs_cg_n = sv_bs_cg_seq_sa(2)
         t_bs_cg_n = n_sv_bs_cg
         os = fm_bs_cg_n - 1
c
c     NB os = 0 is inconceivable---not that it matters!
c
         do 10, bs_cg_n = fm_bs_cg_n, t_bs_cg_n
            sv_bs_cg_en_vr(bs_cg_n-os) = sv_bs_cg_en_vr(bs_cg_n)
            sv_bs_cg_lv_vr(bs_cg_n-os) = sv_bs_cg_lv_vr(bs_cg_n)
            sv_bs_cg_en_vr_prev_st(bs_cg_n-os) =
     &           sv_bs_cg_en_vr_prev_st(bs_cg_n)
            sv_bs_cg_lv_vr_prev_st(bs_cg_n-os) =
     &           sv_bs_cg_lv_vr_prev_st(bs_cg_n)
 10      continue
         n_sv_bs_cg_seq = n_sv_bs_cg_seq - 1
         do 20, seq_n = 1, n_sv_bs_cg_seq
            sv_bs_cg_seq_sa(seq_n) = sv_bs_cg_seq_sa(seq_n+1) - os
 20      continue
         n_sv_bs_cg = n_sv_bs_cg - os
c         print*, '### AFTER REMOVING BASIS CHANGE RECORD SEQUENCE ###'
c         print*, ' n_sv_bs_cg_seq = ', n_sv_bs_cg_seq
c         print*, ' n_sv_bs_cg = ', n_sv_bs_cg
c         call ems_sv_bs_cg_rp(ems_rt_cod)
      endif
c      print*, '### REMOVED BASIS CHANGE RECORD SEQUENCE ###'
c      print*, ' n_sv_bs_cg_seq = ', n_sv_bs_cg_seq
c      print*, ' n_sv_bs_cg = ', n_sv_bs_cg
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' removing saved basis change sequence')
      end
C->>> ---------------------------------------------> ems_sv_bs_cg_rp <<<
c     Report the saved basis change records
c
      subroutine ems_sv_bs_cg_rp(ems_rt_cod)
      implicit none
      include 'EMSMSG.INC'
      include 'SVBSCG.INC'
      include 'EMSRTCOD.INC'
      integer ems_rt_cod
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      integer fm_bs_cg_n, t_bs_cg_n, bs_cg_n, seq_n
      logical alw_f7_wr
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      alw_f7_wr = ems_msg_no_prt_fm .ge. 1
      if (.not. bs_cg_rec) then
         if (alw_f7_wr) write(ems_li, 9000, err = 8990)
         call ems_msg_wr_li(info_msg_n)
         goto 7000
      endif
      t_bs_cg_n = 0
      sv_bs_cg_seq_sa(n_sv_bs_cg_seq+1) = n_sv_bs_cg + 1
      do 20, seq_n = 1, n_sv_bs_cg_seq
         if (alw_f7_wr) write(ems_li, 9100, err = 8990)seq_n
         call ems_msg_wr_li(info_msg_n)
         fm_bs_cg_n = t_bs_cg_n + 1
         t_bs_cg_n = sv_bs_cg_seq_sa(seq_n+1)-1
         if (fm_bs_cg_n .le. t_bs_cg_n) then
            if (alw_f7_wr) write(ems_li, 9110, err = 8990)
            call ems_msg_wr_li(info_msg_n)
            do 10, bs_cg_n = fm_bs_cg_n, t_bs_cg_n
               if (alw_f7_wr) write(ems_li, 9120, err = 8990)bs_cg_n,
     &              sv_bs_cg_en_vr(bs_cg_n),
     &              sv_bs_cg_en_vr_prev_st(bs_cg_n),
     &              sv_bs_cg_lv_vr(bs_cg_n),
     &              sv_bs_cg_lv_vr_prev_st(bs_cg_n)
               call ems_msg_wr_li(info_msg_n)
 10         continue
         endif
 20   continue
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 8990 continue
      ems_rt_cod = max(ems_rt_cod_serious_f7_wr_er,
     &     ems_rt_cod)
      goto 7100
 9000 format('Not saving basis change records')
 9100 format('Saved basis change record sequence ', i2)
 9110 format('BsCgN ',
     &     '  Entering (     Vr, PrevSt)',
     &     '   Leaving (     Vr, PrevSt)')
 9120 format(i5, 2(13x, i7, 5x, i3))
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' reporting saved basis change records')
      end

C->>> ----------------------------------------------> ems_g_bs_cg_st <<<
c     Gets the basis change status for a variable
c
      subroutine ems_g_bs_cg_st(ems_rt_cod, ps_rp_lvl, vr_n,
     &     st, vr_in_c, bs_cg_st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMSG.INC'
      include 'RSMICS.INC'
      include 'RLCTVR.INC'
      include 'EMSRTCOD.INC'
      include 'SVBSCG.INC'
      integer ems_rt_cod, ps_rp_lvl, vr_n, st, bs_cg_st
      integer vr_in_c(-vr_in_c_n_sn:*)
      integer ems_g_vr_in_c_sn_n_o_c
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      character*30 ems_st_t_ch30
      character*30 ch_st
      integer rp_lvl, c_n, vr_in_c_sn_n, vr_in_c_sn_ty
      integer lc_c_n, vr_in_c_l_p

      rp_lvl = ps_rp_lvl
      ems_rt_cod = ems_rt_cod_ok
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) then
c
c     If basis change records are not being saved then don't try to get
c     one. This mechanism is used to avoid addressing vr_in_c which is
c     not set up if basis changes occur due to singularity in the
c     initial reset.
c
         bs_cg_st = bs_cg_st_null
         vr_in_c_sn_ty = 0
         goto 7000
      endif
      if (iand(st, bc_bt) .ne. 0) then
         bs_cg_st = bs_cg_st_bc
         vr_in_c_sn_ty = 0
         goto 7000
      endif
c
c     Variable is nonbasic.
c
      vr_in_c_l_p = vr_in_c(os_vr_in_c_l_p)
c
c     Extract the pointer from the status. Normally this will be the
c     pointer into vr_in_c but, before checking this, have to trim the
c     pointer to the (valid) array bounds of vr_in_c.
c
      c_n = iand(st, mx_mx_ml_a_dim)
      lc_c_n = min(max(c_n, 1), vr_in_c_l_p)
      if (vr_in_c(lc_c_n) .ne. vr_n) then
c
c     c_n is not the pointer into vr_in_c so search for the pointer.
c
         do 10, c_n = 1, vr_in_c_l_p
            if (vr_in_c(c_n) .eq. vr_n) goto 20
 10      continue
         goto 8010
 20      continue
      endif
      vr_in_c_sn_n = ems_g_vr_in_c_sn_n_o_c(c_n, vr_in_c)
      if (ems_msg_cod .ge. ems_msg_lvl_serious) goto 8950
      if (vr_in_c_sn_n .le. vr_in_c_n_sn_ty) then
         vr_in_c_sn_ty = vr_in_c_sn_n
      else
         vr_in_c_sn_ty = vr_in_c_sn_n - vr_in_c_n_sn_ty
      endif
      if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_ab_bp) then
         bs_cg_st = bs_cg_st_bp_vr_ab_bp
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_bw_bp) then
         bs_cg_st = bs_cg_st_bp_vr_bw_bp
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_bw_lb) then
         bs_cg_st = bs_cg_st_vr_bw_lb
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_ab_ub) then
         bs_cg_st = bs_cg_st_vr_ab_ub
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_btw) then
         bs_cg_st = bs_cg_st_vr_btw_bd
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_at_lb) then
         bs_cg_st = bs_cg_st_vr_at_lb
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_at_ub) then
         bs_cg_st = bs_cg_st_vr_at_ub
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_te_fx) then
         bs_cg_st = bs_cg_st_fx_vr_at_bd
      else if (vr_in_c_sn_ty .eq. vr_in_c_sn_ty_fx) then
         bs_cg_st = bs_cg_st_fx_vr_at_bd
      else
         goto 8020
      endif
 7000 continue
CM      IF (emsol_dev .EQ. 1) THEN
C?      if (rp_lvl .ne. 0) then
C?         ch_st = ems_st_t_ch30(st)
C?         write(*, 9000, err=8990)vr_n, vr_in_c_sn_ty, ch_st, bs_cg_st
C?      endif
CM      ENDIF
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 8010 continue
      ch_st = ems_st_t_ch30(st)
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9801, err=8990)
     &     vr_n, ch_st
      call ems_msg_wr_li(serious_msg_n)
      ems_rt_cod = max(ems_rt_lvl_serious,
     &     ems_rt_cod)
      goto 7100
 8020 continue
      ch_st = ems_st_t_ch30(st)
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9802, err=8990)
     &     vr_in_c_sn_ty, vr_n, ch_st
      call ems_msg_wr_li(serious_msg_n)
      ems_rt_cod = max(ems_rt_lvl_serious,
     &     ems_rt_cod)
      goto 7100
 8950 continue
      ems_rt_cod = max(ems_rt_lvl_serious,
     &     ems_rt_cod)
      goto 7100
 8990 continue
      ems_rt_cod = max(ems_rt_cod_serious_f7_wr_er,
     &     ems_rt_cod)
      goto 7100
CM      IF (emsol_dev .EQ. 1) THEN
C? 9000 format('vr = ', i7, ' vr_in_c_sn_ty ', i2, ' st: ', a30,
C?     &     ' BsCgSt = ', i3)
CM      ENDIF
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' getting a basis change status')
 9801 format('Cannot find variable vr = ', i7,
     &     ' in vr_in_c, st: ', a30)
 9802 format('Cannot match vr_in_c_sn_ty = ', i7,
     &     ' for vr = ', i7, ' st: ', a30)
      end

C->>> ---------------------------------------------> ems_sv_bs_cg_ck <<<
c     Check the saved basis change records.
c
      subroutine ems_sv_bs_cg_ck(ems_rt_cod,
     &     prev_reset_bs_cg_st, prev_bs_cg_st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSRTCOD.INC'
      include 'EMSMSG.INC'
      include 'SVBSCG.INC'
      include 'ICTVR.INC'
      integer ems_rt_cod
      integer prev_reset_bs_cg_st(0:mx_n_c+n_r)
      integer prev_bs_cg_st(0:mx_n_c+n_r)
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      integer bs_cg_n, bs_cg_st, en_vr, lv_vr
      integer t_bs_cg_n
      integer en_vr_prev_st, lv_vr_prev_st
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) goto 7000
      t_bs_cg_n = sv_bs_cg_seq_sa(n_sv_bs_cg_seq)
c
c     Pass backwards though this sequence of basis change records,
c     updating the basis change status of all variables before the
c     basis change took place.
c
      do 10, bs_cg_n = n_sv_bs_cg, t_bs_cg_n, -1
         en_vr = sv_bs_cg_en_vr(bs_cg_n)
         lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
         en_vr_prev_st = sv_bs_cg_en_vr_prev_st(bs_cg_n)
         lv_vr_prev_st = sv_bs_cg_lv_vr_prev_st(bs_cg_n)
         prev_bs_cg_st(en_vr) = en_vr_prev_st
         if (en_vr .eq. lv_vr) goto 10
         prev_bs_cg_st(lv_vr) = lv_vr_prev_st
 10   continue
      do 20, bs_cg_n = n_sv_bs_cg, t_bs_cg_n, -1
         en_vr = sv_bs_cg_en_vr(bs_cg_n)
         lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
         bs_cg_st = prev_bs_cg_st(en_vr)
         if (bs_cg_st .ne. prev_reset_bs_cg_st(en_vr)) then
            ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
CM      IF (emsol_dev .EQ. 1) THEN
C?            write(*, 9999)bs_cg_n, 'Entering', en_vr,
C?     &           bs_cg_st, prev_reset_bs_cg_st(en_vr)
CM      ENDIF
         endif
         bs_cg_st = prev_bs_cg_st(lv_vr)
         if (bs_cg_st .ne. prev_reset_bs_cg_st(lv_vr)) then
            ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
CM      IF (emsol_dev .EQ. 1) THEN
C?            write(*, 9999)bs_cg_n, 'Leaving ', lv_vr,
C?     &           bs_cg_st, prev_reset_bs_cg_st(lv_vr)
CM      ENDIF
         endif
 20   continue
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' checking the saved basis change records')
CM      IF (emsol_dev .EQ. 1) THEN
C? 9999 format('Basis change ', i5, ': ', a, ' variable ', i7,
C?     &     ': prev_bs_cg_st = ', i3,
C?     &     ' ***BUT*** prev_reset_bs_cg_st = ', i3)
CM      ENDIF
      end

C->>> -------------------------------------------> ems_ck_reset_loop <<<
c     Check to see whether the saved basis change records indicate that
c     the same basis is being reset.
c
      subroutine ems_ck_reset_loop(ems_rt_cod, ds, is,
     &     prev_reset_bs_cg_st, prev_cu_bs_cg_st, reset_loop)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'SVBSCG.INC'
      include 'RSMICS.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSRTCOD.INC'
      integer ems_rt_cod
      logical reset_loop
      integer prev_reset_bs_cg_st(0:mx_n_c+n_r)
      integer prev_cu_bs_cg_st(0:mx_n_c+n_r)
      integer is(0:*)
      double precision ds(0:*)
CM      IF (emsol_epc .EQ. 1) THEN
C?      logical ems_i1_eq_i2
CM      ENDIF
      integer bs_cg_n, cu_bs_cg_st, prev_bs_cg_st, seq_n, en_vr, lv_vr
      integer fm_bs_cg_n, t_bs_cg_n
      integer en_vr_prev_st, lv_vr_prev_st
      integer ix_n, vr_n, bs_cg_st
      integer ca_ems_rt_cod
      logical er_fd
c
c     Make sure that ems_sv_bs_cg_iz has been called. This protects the
c     test of bs_cg_rec which may be unassigned.
c
      ems_rt_cod = ems_rt_cod_ok
      reset_loop = .false.
      er_fd = .false.
CM      IF (emsol_epc .EQ. 1) THEN
C?c
C?c     Have to use a function compiled without unassigned variable
C?c     checking in order to test values which may be unassigned.
C?c
C?      if (ems_i1_eq_i2(ca_sv_bs_cg_iz_fg1,
C?     &     ca_sv_bs_cg_iz_fg2)) goto 8000
CM      ELSE
      if (ca_sv_bs_cg_iz_fg1 .eq. ca_sv_bs_cg_iz_fg2) goto 8000
CM      ENDIF
      if (.not. bs_cg_rec) goto 7000
      if (n_sv_bs_cg .eq. 0) goto 2000
c
c     Get the current basis change status for variables which have been
c     involved in basis changes.
c
      do 5, vr_n = 1, mx_n_c+n_r
         prev_cu_bs_cg_st(vr_n) = 0
 5    continue
      do 10, bs_cg_n = n_sv_bs_cg, 1, -1
         en_vr = sv_bs_cg_en_vr(bs_cg_n)
c
c     Get the current basis change status for the entering variable
c     (unless it is already known)
c
         if (prev_cu_bs_cg_st(en_vr) .eq. 0) then
            if (iand(is(p_st+en_vr), bc_bt) .ne. 0) then
c
c     Likely to be basic...
c
               prev_cu_bs_cg_st(en_vr) = bs_cg_st_bc
            else
c
c     ... but nonbasic if this was a bound swap.
c
c               call ems_og_g_bs_cg_st(ca_ems_rt_cod, 0, en_vr,
c     &              is(p_st+en_vr),
c     &              ds(p_pr_act+en_vr),
c     &              ds(p_rsmi_lb+en_vr),
c     &              ds(p_rsmi_ub+en_vr),
c     &              bs_cg_st)
               call ems_g_bs_cg_st(ca_ems_rt_cod, 0, en_vr,
     &              is(p_st+en_vr),
     &              is(p_vr_in_c),
     &              bs_cg_st)
               if (ca_ems_rt_cod .ne. ems_rt_cod_ok) then
                  ems_rt_cod = max(ca_ems_rt_cod, ems_rt_cod)
                  if (ems_rt_cod .ge. ems_rt_lvl_serious) goto 7100
               endif
               prev_cu_bs_cg_st(en_vr) = bs_cg_st
            endif
         endif
c
c     Get the current basis change status for the leaving variable
c     (unless it is already known)
c
         lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
         if (prev_cu_bs_cg_st(lv_vr) .eq. 0) then
            if (iand(is(p_st+lv_vr), bc_bt) .ne. 0) then
c
c     Shouldn't be basic now since this is the last time this variable
c     left the basis.
c
               goto 8010
            else
c
c     Should be nonbasic
c
c               call ems_og_g_bs_cg_st(ca_ems_rt_cod, 0, lv_vr,
c     &              is(p_st+lv_vr),
c     &              ds(p_pr_act+lv_vr),
c     &              ds(p_rsmi_lb+lv_vr),
c     &              ds(p_rsmi_ub+lv_vr),
c     &              bs_cg_st)
               call ems_g_bs_cg_st(ca_ems_rt_cod, 0, lv_vr,
     &              is(p_st+lv_vr),
     &              is(p_vr_in_c),
     &              bs_cg_st)
               if (ca_ems_rt_cod .ne. ems_rt_cod_ok) then
                  ems_rt_cod = max(ca_ems_rt_cod, ems_rt_cod)
                  if (ems_rt_cod .ge. ems_rt_lvl_serious) goto 7100
               endif
               prev_cu_bs_cg_st(lv_vr) = bs_cg_st
            endif
         endif
 10   continue
c
c     Pass backwards through the sequences of saved basis change records
c     to see whether the current basis has already been reset.
c
      seq_n = n_sv_bs_cg_seq
      fm_bs_cg_n = n_sv_bs_cg
 100  continue
      t_bs_cg_n = sv_bs_cg_seq_sa(seq_n)
c
c     Pass backwards though this sequence of basis change records,
c     updating the basis change status of all variables before the
c     basis change took place.
c
      do 120, bs_cg_n = fm_bs_cg_n, t_bs_cg_n, -1
         en_vr = sv_bs_cg_en_vr(bs_cg_n)
         lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
         en_vr_prev_st = sv_bs_cg_en_vr_prev_st(bs_cg_n)
         lv_vr_prev_st = sv_bs_cg_lv_vr_prev_st(bs_cg_n)
         bs_cg_st = prev_cu_bs_cg_st(en_vr)/prev_bs_cg_st_fac
         prev_cu_bs_cg_st(en_vr) = prev_cu_bs_cg_st(en_vr) +
     &        (en_vr_prev_st - bs_cg_st)*prev_bs_cg_st_fac
         if (en_vr .eq. lv_vr) goto 120
         bs_cg_st = prev_cu_bs_cg_st(lv_vr)/prev_bs_cg_st_fac
         prev_cu_bs_cg_st(lv_vr) = prev_cu_bs_cg_st(lv_vr) +
     &        (lv_vr_prev_st - bs_cg_st)*prev_bs_cg_st_fac
 120  continue
      if (seq_n .eq. n_sv_bs_cg_seq) then
c
c     Check the basis change status deduced by passing backwards though
c     the last sequence of basis change records with the value recorded
c     when the previous reset occurred.
c
         do 130, bs_cg_n = n_sv_bs_cg, t_bs_cg_n, -1
            en_vr = sv_bs_cg_en_vr(bs_cg_n)
            lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
            prev_bs_cg_st = prev_cu_bs_cg_st(en_vr)/prev_bs_cg_st_fac
            cu_bs_cg_st = prev_cu_bs_cg_st(en_vr) -
     &           prev_bs_cg_st*prev_bs_cg_st_fac
            if (prev_bs_cg_st .ne. prev_reset_bs_cg_st(en_vr)) then
               er_fd = .true.
CM      IF (emsol_dev .EQ. 1) THEN
C?               write(*, 9999)bs_cg_n, 'Entering', en_vr,
C?     &              prev_bs_cg_st, prev_reset_bs_cg_st(en_vr)
CM      ENDIF
            endif
            prev_bs_cg_st = prev_cu_bs_cg_st(lv_vr)/prev_bs_cg_st_fac
            cu_bs_cg_st = prev_cu_bs_cg_st(lv_vr) -
     &           prev_bs_cg_st*prev_bs_cg_st_fac
            if (prev_bs_cg_st .ne. prev_reset_bs_cg_st(lv_vr)) then
               er_fd = .true.
CM      IF (emsol_dev .EQ. 1) THEN
C?               write(*, 9999)bs_cg_n, 'Leaving ', lv_vr,
C?     &              prev_bs_cg_st, prev_reset_bs_cg_st(lv_vr)
CM      ENDIF
            endif
 130     continue
      endif
c
c     If all the basis change statuses are the same then a reset loop
c     has occurred.
c
      reset_loop = .false.
      do 140, bs_cg_n = n_sv_bs_cg, t_bs_cg_n, -1
         en_vr = sv_bs_cg_en_vr(bs_cg_n)
         lv_vr = sv_bs_cg_lv_vr(bs_cg_n)
         prev_bs_cg_st = prev_cu_bs_cg_st(en_vr)/prev_bs_cg_st_fac
         cu_bs_cg_st = prev_cu_bs_cg_st(en_vr) -
     &        prev_bs_cg_st*prev_bs_cg_st_fac
         if (prev_bs_cg_st .ne. cu_bs_cg_st) goto 150
         prev_bs_cg_st = prev_cu_bs_cg_st(lv_vr)/prev_bs_cg_st_fac
         cu_bs_cg_st = prev_cu_bs_cg_st(lv_vr) -
     &        prev_bs_cg_st*prev_bs_cg_st_fac
         if (prev_bs_cg_st .ne. cu_bs_cg_st) goto 150
 140  continue
      reset_loop = .true.
      goto 1000
 150  continue
c
c     Set up for the previous sequence (if there is one).
c
      seq_n = seq_n - 1
      if (seq_n .le. 0) goto 1000
      fm_bs_cg_n = t_bs_cg_n - 1
      goto 100
 1000 continue
c
c     Remove the oldest sequence if the sequence limit has been reached,
c     otherwise incerease the sequence counter by 1.
c
      if (n_sv_bs_cg_seq .eq. mx_n_sv_bs_cg_seq) then
         call ems_sv_bs_cg_rm_ol_seq(ca_ems_rt_cod)
         if (ca_ems_rt_cod .ne. ems_rt_cod_ok) then
            ems_rt_cod = max(ca_ems_rt_cod, ems_rt_cod)
            if (ems_rt_cod .ge. ems_rt_lvl_serious) goto 7100
         endif
      endif
      n_sv_bs_cg_seq = n_sv_bs_cg_seq + 1
      sv_bs_cg_seq_sa(n_sv_bs_cg_seq) = n_sv_bs_cg + 1
 2000 continue
c
c     Save the basis change status at this reset for all variables.
c     This is just for checking.
c
      do 2010, ix_n = 1, n_c+n_r
         if (ix_n .le. n_c) then
            vr_n = ix_n
         else
            vr_n = ix_n + (mx_n_c-n_c)
         endif
         if (iand(is(p_st+vr_n), bc_bt) .ne. 0) then
            prev_reset_bs_cg_st(vr_n) = bs_cg_st_bc
         else
c            call ems_og_g_bs_cg_st(ca_ems_rt_cod, 0, vr_n,
c     &           is(p_st+vr_n),
c     &           ds(p_pr_act+vr_n),
c     &           ds(p_rsmi_lb+vr_n),
c     &           ds(p_rsmi_ub+vr_n),
c     &           cu_bs_cg_st)
            call ems_g_bs_cg_st(ca_ems_rt_cod, 0, vr_n,
     &           is(p_st+vr_n),
     &           is(p_vr_in_c),
     &           cu_bs_cg_st)
            if (ca_ems_rt_cod .ne. ems_rt_cod_ok) then
               ems_rt_cod = max(ca_ems_rt_cod, ems_rt_cod)
               if (ems_rt_cod .ge. ems_rt_lvl_serious) goto 7100
            endif
            prev_reset_bs_cg_st(vr_n) = cu_bs_cg_st
         endif
 2010 continue
      if (er_fd) ems_rt_cod = max(ems_rt_lvl_serious, ems_rt_cod)
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      ems_rt_cod = max(ems_rt_cod_int_er, ems_rt_cod)
      goto 7100
 8010 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 98011)lv_vr
      call ems_msg_wr_li(bug_msg_n)
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 98012)en_vr,
     &     prev_cu_bs_cg_st(en_vr)
      call ems_msg_wr_li(bug_msg_n)
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 98013)bs_cg_n
      call ems_msg_wr_li(bug_msg_n)
      goto 7100
 9800 format('Must call ems_sv_bs_cg_iz before',
     &     ' checking for a reset loop')
CM      IF (emsol_dev .EQ. 1) THEN
C? 9999 format('Basis change ', i5, ': ', a, ' variable ', i7,
C?     &     ': prev_bs_cg_st = ', i3,
C?     &     ' ***BUT*** prev_reset_bs_cg_st = ', i3)
CM      ENDIF
98011 format('Leaving  variable ', i7, ' should not be basic')
98012 format('Entering variable ', i7, ' St = ', i9)
98013 format('bs_cg_n = ', i7)
      end
C=======================================================================
c     Redundant code ---->
C=======================================================================
C->>> -------------------------------------------> ems_og_g_bs_cg_st <<<
c     Gets the basis change status for a variable
c
      subroutine ems_og_g_bs_cg_st(ems_rt_cod, rp_lvl, vr_n,
     &     st, pr_act, lb, ub, bs_cg_st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICS.INC'
      include 'RLCTVR.INC'
      include 'EMSRTCOD.INC'
      integer ems_rt_cod, rp_lvl, vr_n, st, bs_cg_st
      double precision pr_act, lb, ub
      double precision rsdu, bp
      character*30 ems_st_t_ch30
      character*30 ch_st

      ems_rt_cod = ems_rt_cod_ok
      if (iand(st, alt_bt) .eq. 0) then
         if (iand(st, lb_bt) .ne. 0) then
            rsdu = lb - pr_act
            if (rsdu .gt. tl_pr_ifs) then
               bs_cg_st = bs_cg_st_vr_bw_lb
            else if (rsdu .ge. -tl_pr_ifs) then
               bs_cg_st = bs_cg_st_vr_at_lb
            else
c
c     Variable is above its lower bound
c
               if (iand(st, ub_bt) .ne. 0) then
                  rsdu = pr_act - ub
                  if (rsdu .gt. tl_pr_ifs) then
                     bs_cg_st = bs_cg_st_vr_ab_ub
                  else if (rsdu .ge. -tl_pr_ifs) then
                     bs_cg_st = bs_cg_st_vr_at_ub
                  else
                     bs_cg_st = bs_cg_st_vr_btw_bd
                  endif
               else
                  bs_cg_st = bs_cg_st_vr_ab_lb
               endif
            endif
         else if (iand(st, ub_bt) .ne. 0) then
            rsdu = pr_act - ub
            if (rsdu .gt. tl_pr_ifs) then
               bs_cg_st = bs_cg_st_vr_ab_ub
            else if (rsdu .ge. -tl_pr_ifs) then
               bs_cg_st = bs_cg_st_vr_at_ub
            else
               bs_cg_st = bs_cg_st_vr_bw_ub
            endif
         else
            bs_cg_st = bs_cg_st_fr_vr
         endif
      else if (iand(st, bp_bt) .ne. 0) then
         if (iand(st, ub_bt) .ne. 0) then
            bp = ub
         else
            bp = lb
         endif
         rsdu = pr_act - bp
         if (rsdu .gt. tl_pr_ifs) then
            bs_cg_st = bs_cg_st_bp_vr_ab_bp
         else if (rsdu .ge. -tl_pr_ifs) then
            bs_cg_st = bs_cg_st_bp_vr_at_bp
         else
            bs_cg_st = bs_cg_st_bp_vr_bw_bp
         endif
      else
         goto 8000
      endif
      if (rp_lvl .ne. 0) then
         ch_st = ems_st_t_ch30(st)
CM      IF (emsol_dev .EQ. 1) THEN
C?         write(*, 9000, err=8990)vr_n, lb, pr_act, ub, ch_st, bs_cg_st
CM      ENDIF
      endif
c 7000 continue
 7100 continue
      return
 8000 continue
      ems_rt_cod = max(ems_rt_lvl_serious,
     &     ems_rt_cod)
      goto 7100
CM      IF (emsol_dev .EQ. 1) THEN
C? 8990 continue
C?      ems_rt_cod = max(ems_rt_cod_serious_f7_wr_er,
C?     &     ems_rt_cod)
C?      goto 7100
CM      ENDIF
CM      IF (emsol_dev .EQ. 1) THEN
C? 9000 format('vr = ', i7, ' lb:act:ub ', 3(2x, g11.4), ' st: ', a30,
C?     &     ' BsCgSt = ', i3)
CM      ENDIF
      end

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


text:
subroutine disper(kmxdt     ,kmxt      ,np        ,nfreqs    ,angle     , &
                & bvmx      ,freqlo    ,frequp    ,kcri      ,ncri      , &
                & sing      ,umean     ,vmean     ,atg       ,btg       , &
                & bv2       ,ctg       ,d2u       ,d2v       ,h0        , &
                & kbed      ,ktop      ,qz        ,r1tg      ,scale     , &
                & utg       ,vtg       ,xkh       ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: disper.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/taylor/disper.f90 $
!!--description-----------------------------------------------------------------
!
!      Function: Finding extreme levels of critical layers for
!                all allowed
!                IW's propagating in positive (nr=2) or negative
!                (nr=1) direction with respect to given hori-
!                zontal wavenumber magnitude and ANGLE.
!
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp) , pointer :: bvmin
!
! Global variables
!
    integer         :: kbed
                                   !!  K-value on TG grid: upper level
                                   !!  stratified layer
    integer         :: kmxdt !  Description and declaration in dimens.igs
    integer         :: kmxt !  Description and declaration in dimens.igs
    integer         :: ktop
                                   !!  K-value on TG grid: lower level
                                   !!  stratified layer
    integer, intent(in)            :: nfreqs !  Description and declaration in dimens.igs
    integer, intent(in)            :: np
                                   !!  Dimension for IWE frequency arrays
                                   !!  MAX (1,NFREQS)
    integer, dimension(2, 2), intent(out) :: kcri
                                   !!  K-value on TG grid of most upper and
                                   !!  lower critical levels of turbulence-
                                   !!  generated IW's
    integer, dimension(2, 2), intent(out) :: ncri
                                   !!   Number of critical layers found per
                                   !!  minimal or maximal K-value as well as
                                   !!  sign of frequency of turbulence-
                                   !!  generated IW's
    logical, dimension(2, 2), intent(out) :: sing
                                   !!  Indicator of singular solutions of TG
                                   !!  equation
    real(fp)        :: angle
                                   !!  Angle [degrees] between horizontal
                                   !!  IW wavenumber vector and x-u-axis
    real(fp), intent(in)               :: bvmx
                                   !!  Maximal buoyancy frequency
    real(fp)        :: h0
                                   !!  Water depth
    real(fp)        :: scale
                                   !!  Length scale for scaling
    real(fp), intent(in)               :: umean !  Description and declaration in esm_alloc_real.f90
    real(fp), intent(in)               :: vmean !  Description and declaration in esm_alloc_real.f90
    real(fp)        :: xkh
                                   !!  Non-dimensional XK
    real(fp), dimension(0:kmxdt) :: atg
    real(fp), dimension(0:kmxdt) :: btg
    real(fp), dimension(0:kmxdt) :: bv2
    real(fp), dimension(0:kmxdt) :: ctg
    real(fp), dimension(0:kmxdt) :: d2u
    real(fp), dimension(0:kmxdt) :: d2v
    real(fp), dimension(0:kmxdt) :: qz
    real(fp), dimension(0:kmxdt) :: r1tg
    real(fp), dimension(0:kmxdt) :: utg
    real(fp), dimension(0:kmxdt) :: vtg
    real(fp), dimension(np, 2), intent(out) :: freqlo
    real(fp), dimension(np, 2), intent(out) :: frequp
!
!
! Local variables
!
    integer                        :: i
    integer                        :: j
    integer                        :: kcrit
    integer                        :: ndir
    integer                        :: nf
    integer                        :: nr
    logical                        :: singul
    logical                        :: tgfind
    logical                        :: top
    real(fp)                       :: crit
    real(fp)                       :: delf
    real(fp)                       :: flo
    real(fp)                       :: frlo
    real(fp)                       :: frup
    real(fp)                       :: fup
    real(fp)                       :: omeg
    real(fp)                       :: omeg1
    real(fp)                       :: umag
!
!
!! executable statements -------------------------------------------------------
!
    !
    !
    bvmin  => gdp%gdiwearr%bvmin
    !
    umag = sqrt(umean**2 + vmean**2)
    if (umag<1.E-12) then
       angle = 0.
    else
       angle = 45.*acos(umean/umag)/atan(1.)
    endif
    if (vmean<0.0) angle = 360. - angle
    !
    do i = 1, 2
       do j = 1, 2
          kcri(i, j) = 0
          ncri(i, j) = 0
          sing(i, j) = .false.
       enddo
    enddo
    !
    ! Inspect two ranges of normalized eigen frequencies:
    ! nr=1: [-bvmx,-delf]
    ! nr=2: [ delf, bvmx]
    !
    delf = bvmx/nfreqs
    !
    do nr = 1, 2
       i = -1
       if (nr==2) i = 1
       do ndir = 1, 2
          if (ndir==1) then
             j = 1
             omeg = i*bvmin
          else
             j = -1
             omeg = i*bvmx
          endif
          !
          nf = 0
          tgfind = .false.
          top = .true.
          !
          call tg(kmxdt     ,kmxt      ,angle     ,atg       ,btg       , &
                & bv2       ,ctg       ,crit      ,d2u       ,d2v       , &
                & h0        ,kbed      ,ktop      ,omeg      ,qz        , &
                & r1tg      ,scale     ,tgfind    ,top       ,utg       , &
                & vtg       ,xkh       ,gdp       )
          !
          !
          ! For FLO and FLUP only the sign is important (zero criteria)
          !
          flo = sign(1._fp, crit)
          frlo = omeg
          ! -->
   50     continue
          nf = nf + 1
          omeg = omeg + i*j*delf
          tgfind = .false.
          top = .true.
          call tg(kmxdt     ,kmxt      ,angle     ,atg       ,btg       , &
                & bv2       ,ctg       ,crit      ,d2u       ,d2v       , &
                & h0        ,kbed      ,ktop      ,omeg      ,qz        , &
                & r1tg      ,scale     ,tgfind    ,top       ,utg       , &
                & vtg       ,xkh       ,gdp       )
          !
          !
          fup = sign(1._fp, crit)
          frup = omeg
          if (nf<=nfreqs) then
             if (fup*flo>0.) then
                flo = fup
                frlo = frup
                goto 50
             ! <--
             endif
             omeg1 = omeg
             omeg = 0.5*(frlo + frup)
             call seliwe(angle     ,kbed      ,ktop      ,kmxt      ,omeg      , &
                       & r1tg      ,utg       ,vtg       ,xkh       ,singul    , &
                       & top       ,kcrit     ,gdp       )
             !
             omeg = omeg1
             if (singul) then
                frequp(nf, nr) = frup
                freqlo(nf, nr) = frlo
                kcri(nr, ndir) = kcrit
                ncri(nr, ndir) = nf
                sing(nr, ndir) = singul
             else
                goto 50
             ! <--
             endif
          endif
       enddo
    enddo
end subroutine disper

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C$Procedure SBINIT ( String buffer, initialize )

      SUBROUTINE SBINIT_1 ( NSIZE, PSIZE, VDIM, NAMES,
     .                                        PTRS,
     .                                        BUFFER )
      IMPLICIT NONE

C$ Abstract
C
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     CB, LB, SB
C
C$ Keywords
C
C     ASCII
C     CHARACTER
C     STRING
C     TEXT
C
C$ Declarations

      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )

      INTEGER               LBCBUF
      PARAMETER           ( LBCBUF =  0 )

      INTEGER               NSIZE
      INTEGER               PSIZE
      INTEGER               VDIM
      CHARACTER*(*)         NAMES       ( LBCELL:NSIZE )
      INTEGER               PTRS        ( LBCELL:PSIZE )
      CHARACTER*(*)         BUFFER      ( LBCBUF:VDIM  )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     NSIZE      I   Name size.
C     PSIZE      I   Pointer size.
C     VDIM       I   Value dimension
C     NAMES,
C     PTRS,
C     BUFFER    I,O  String buffer.
C
C$ Detailed_Input
C
C     NAMES       is a character cell array to be used as the name
C                 component of a string buffer.
C
C     NSIZE       is the declared dimension of NAMES.
C
C     PTRS        is an integer cell array to be used as the pointer
C                 component of a string buffer.
C
C     PSIZE       is the declared dimension of PTRS.
C
C     BUFFER      is a character buffer array to be used as the
C                 character component of a string buffer.
C
C     VDIM        is the declared dimension of BUFFER.
C
C$ Detailed_Output
C
C     NAMES,
C     PTRS,
C     BUFFER      together are an initialized string buffer.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     1) If the size of the pointer array is not sufficient to
C        hold pointers for the maximum number of strings, the
C        error 'SPICE(SBINSUFPTRSIZE)' is signalled.
C
C$ Particulars
C
C     A string buffer must be initialized to allow subsequent
C     operations on the buffer to detect possible overflows.
C     All three components of the buffer are initialized by a
C     single call to SBINIT.
C
C     In order to make full use of the name cell of the string buffer,
C     the arrays and name buffers should be declared as shown below.
C
C       INTEGER              LBCELL
C       PARAMETER          ( LBCELL = -5 )
C
C       INTEGER              MAXNL
C       PARAMETER          ( MAXNL  = maximum name length )
C
C       INTEGER              MAXN
C       PARAMETER          ( MAXN   = maximum number of names )
C
C       CHARACTER*(MAXNL)    NAMES ( MAXN         )
C       INTEGER              PTRS  ( MAXN * 4 + 4 )
C
C     The character buffer portion of the string buffer should be
C     declared as shown below.
C
C       INTEGER              MAXL
C       PARAMETER          ( MAXL = maximum expected string length )
C
C       INTEGER              AVGL
C       PARAMETER          ( AVGL = average expected string length )
C
C       INTEGER              LBCBUF
C       PARAMETER          ( LBCBUF = 0  )
C
C       CHARACTER*(MAXL)     BUFFER ( LBCBUF:(MAXN * AVGL) / MAXL + 1 )
C
C$ Examples
C
C     The following code fragment illustrates the initialization
C     of a typical string buffer.
C
C        INTEGER               LBCELL
C        PARAMETER           ( LBCELL =    -5 )
C
C        INTEGER               LBCBUF
C        PARAMETER           ( LBCBUF =     0 )
C
C        CHARACTER*(32)        NAMES    ( LBCELL:1000    )
C        INTEGER               PTRS     ( LBCELL:4004    )
C        CHARACTER*(250)       BUFFER   ( LBCBUF:100     )
C         .
C         .
C
C        CALL SBINIT ( MAXN, PSIZE, BUFDIM, NAMES, PTRS, BUFFER )
C
C     In this example, the buffer may be used to store up to 1000
C     strings averaging 25 characters per string, or 25,000 total
C     characters. The length of any particular string may range from
C     a single character to the entire 25,000 characters. The names
C     used to identify the strings may contain up to 32 characters.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     Dagny Taggart, (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 19-JAN-1989 (DT)
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Local variables
C
      INTEGER               MAXPTR

C
C     Standard error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'SBINIT_1' )
      END IF

C
C     Make sure that the line buffer is large enough (but ONLY large
C     enough) to hold the maximum number of strings. The name list
C     should be empty. The LB should be initialized as a unit.
C
      MAXPTR = 4 * ( NSIZE + 1 )

      IF ( PSIZE .LT. MAXPTR ) THEN
         CALL SIGERR ( 'SPICE(SBINSUFPTRSIZE)' )

      ELSE
         CALL SSIZEC ( NSIZE,        NAMES         )
         CALL LBINIT_1 ( MAXPTR, VDIM, PTRS, BUFFER )
      END IF

      CALL CHKOUT ( 'SBINIT_1' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


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

C=======================================================================
      SUBROUTINE SPPLOT (A, NEUTRL, NPTIMS, IPTIMS, TIMES,
     *  NENUM, DIST, PLTVAL, NAMES, IXSEGV, ISEGEL, SQDIST, LIDSP,
     &   BLKCOL, MAPEL, MAPND)
C=======================================================================

C   --*** SPPLOT *** (SPLOT) Plot curves for all times
C   --   Modified by John Glick - 11/9/88
C   --   Written by Amy Gilkey - revised 01/22/88
C   --
C   --SPPLOT generates all the plots requested for the plot set (for
C   --all plot times).  It scales the X axis and the Y axis (as requested).
C   --It then processes each curve by reading the data and calling SPPLT1.
C   --
C   --Parameters:
C   --   A - IN - the dynamic memory base array
C   --   NEUTRL - IN  - the type of neutral file to write.
C   --   NPTIMS - IN - the number of selected times
C   --   IPTIMS - IN - the selected time steps
C   --   TIMES - IN - the database times
C   --   NENUM - IN - the selected node/element numbers
C   --   DIST - IN - the plot distances
C   --   PLTVAL - IN - the plot data
C   --   NAMES - IN - the variable names
C   --   IXSEGV - IN - the index into ISEGEL for "segmented" variables;
C   --      0 if all selected elements are defined
C   --   ISEGEL - IN - the NENUM indices of the defined elements;
C   --      ISEGEL(0,i) = the number of defined elements
C   --   SQDIST - SCRATCH - size = NNENUM
C   --   LIDSP(0:*)  - IN - the indices of the selected variables
C   --          whose values will be displayed on the plot legend.
C   --          LIDSP(0) = the number of variables in the list.
C   --          LIDSP(i) identifies the ith variable in the list.
C   --          If LIDSP(i) > 0, LIDSP(i) is the id of a history variable.
C   --          If LIDSP(i) < 0, -LIDSP(i) is the id of a global variable.
C   --          If LIDSP(i) = 0, TIME is to be displayed on the plot legend
C   --   BLKCOL - IN/OUT - the user selected colors of the element blocks.
C   --                    BLKCOL(0) = 1 if the user defined material
C   --                                colors should be used in mesh plots.
C   --                              = -1 if program selected colors should
C   --                                be used.
C   --                    BLKCOL(i) = the user selected color of element
C   --                               block i:
C   --                                  -2 - no color selected by user.
C   --                                  -1 - black
C   --                                   0 - white
C   --                                   1 - red
C   --                                   2 - green
C   --                                   3 - yellow
C   --                                   4 - blue
C   --                                   5 - cyan
C   --                                   6 - magenta
C   --
C   --Common Variables:
C   --   Uses NNENUM of /SELNE/
C   --   Uses NSPVAR of /SPVARS/
C   --   Uses DOGRID, ISYTYP, OVERLY, OVERTM of /XYOPT/
C   --   Uses IXSCAL, IYSCAL of /XYLIM/
C   --   Sets XMIN, XMAX, YMIN, YMAX of /XYLIM/

      PARAMETER (NUMSYM = 6, NUMLIN = 6)

      include 'params.blk'
      include 'dbnums.blk'
      include 'selne.blk'
      include 'spvars.blk'
      include 'xyopt.blk'
      include 'xylim.blk'

      DIMENSION A(*)
      INTEGER IPTIMS(NPTIMS)
      REAL TIMES(*)
      INTEGER NENUM(NNENUM)
      REAL DIST(NNENUM)
      REAL PLTVAL(NNENUM,NSPVAR,*)
      CHARACTER*(*) NAMES(*)
      INTEGER IXSEGV(*)
      INTEGER ISEGEL(0:NNENUM,*)
      REAL SQDIST(NNENUM)
      INTEGER LIDSP(0:*)
      INTEGER BLKCOL(0:NELBLK)
      INTEGER MAPEL(*), MAPND(*)

      LOGICAL NUMCRV
      CHARACTER*1024 TXLAB, TYLAB

      LOGICAL SVOVER, SVOVTM
      CHARACTER*8 SVLSID

C   --Save user-set parameters and set program parameters (to eliminate
C   --checks for senseless conditions)

      SVOVTM = OVERTM
      SVOVER = OVERLY
      SVLSID = LABSID

      IF ((NEUTRL .NE. 0) .OR. (NSPVAR .LE. 1)) OVERLY = .FALSE.
      IF ((NEUTRL .NE. 0) .OR. (NPTIMS .LE. 1)) OVERTM = .FALSE.

      IF (IXSCAL .NE. 'SET') THEN
         IXSCAL = 'ALL'
      END IF
      IF (IYSCAL .NE. 'SET') THEN
         IYSCAL = IAXSCA
         IF (NEUTRL .GT. 0) THEN
            IYSCAL = 'CURVE'
         ELSE IF ((NSPVAR .LE. 1) .AND. (NPTIMS .LE. 1)) THEN
            IYSCAL = 'ALL'
         ELSE IF ((OVERLY .OR. (NSPVAR .LE. 1))
     &      .AND. (OVERTM .OR. (NPTIMS .LE. 1))) THEN
            IF (IYSCAL .EQ. 'PLOT') IYSCAL = 'ALL'
         ELSE IF (.NOT. (OVERLY .OR. OVERTM)) THEN
            IF (IYSCAL .EQ. 'CURVE') IYSCAL = 'PLOT'
         END IF
      END IF

      NC = 1
      IF (OVERLY) NC = NPTIMS
      IF (OVERTM) NC = NSPVAR
      IF (((.NOT. OVERLY) .AND. (.NOT. OVERTM))
     &   .OR. ((ISYTYP .LT. 0) .AND. (NUMSYM .GE. NC))
     &   .OR. ((LINTYP .LT. 0) .AND. (NUMLIN .GE. NC)))
     &   LABSID = 'NONE'

      NUMCRV = (LABSID .NE. 'NONE')

C   --Compact variables which have undefined element values (distances
C   --must be compacted for each variable at a later time)

      DO 100 N = 1, NSPVAR
         IF (IXSEGV(N) .GT. 0) THEN
            IX = IXSEGV(N)
            CALL SQZIXV (ISEGEL(0,IX), ISEGEL(1,IX),
     &         PLTVAL(1,N,NPT), PLTVAL(1,N,NPT))
         END IF
  100 CONTINUE

C   --Scale distance parameters (only done once)

      IF (IXSCAL .NE. 'SET') THEN
         CALL MINMAX (NNENUM, DIST, XMIN, XMAX)
         IF (NEUTRL .EQ. 0)
     *     CALL EXPMAX (LABSID, XMIN, XMAX)
      END IF

C   --Scale Y axis if same scale requested

      IF (IYSCAL .EQ. 'ALL') THEN
         DO 120 NPT = 1, NPTIMS
            DO 110 N = 1, NSPVAR
               IF (IXSEGV(N) .LE. 0) THEN
                  NNE = NNENUM
               ELSE
                  NNE = ISEGEL(0,IXSEGV(N))
               END IF
               CALL MINMAX (NNE, PLTVAL(1,N,NPT), VMIN, VMAX)
               IF ((NPT .EQ. 1) .AND. (N .EQ. 1)) THEN
                  YMIN = VMIN
                  YMAX = VMAX
               ELSE
                  YMIN = MIN (YMIN, VMIN)
                  YMAX = MAX (YMAX, VMAX)
               END IF
  110       CONTINUE
  120    CONTINUE
         call adjrng(ymin, ymax)
         CALL EXPMAX (' ', YMIN, YMAX)
      END IF

C   --Set the total number of plots

      IF (OVERTM) THEN
         NPTOT = NSPVAR
      ELSE IF (OVERLY) THEN
         NPTOT = NPTIMS
      ELSE
         NPTOT = NSPVAR * NPTIMS
      END IF

      IF (OVERTM) THEN

         DO 160 N = 1, NSPVAR
            IF (IXSEGV(N) .LE. 0) THEN
               NNE = NNENUM
            ELSE
               NNE = ISEGEL(0,IXSEGV(N))
            END IF

C         --Calculate axis limits if to be done once for all times

            IF (IYSCAL .EQ. 'PLOT') THEN
               DO 130 NPT = 1, NPTIMS
                  CALL MINMAX (NNE, PLTVAL(1,N,NPT), VMIN, VMAX)
                  IF (NPT .EQ. 1) THEN
                     YMIN = VMIN
                     YMAX = VMAX
                  ELSE
                     YMIN = MIN (YMIN, VMIN)
                     YMAX = MAX (YMAX, VMAX)
                  END IF
  130          CONTINUE
               call adjrng(ymin, ymax)
               CALL EXPMAX (' ', YMIN, YMAX)
            END IF

C         --Label plot if overlaid by time

  140       CONTINUE
            CALL SPLAB (A, NPTIMS, IPTIMS, TIMES, NENUM,
     &           N, 1, NUMCRV, NAMES, TXLAB, TYLAB,
     &           LIDSP, BLKCOL,  MAPEL, MAPND, *210)
            IF (IYSCAL .NE. 'CURVE')
     &           CALL XYAXIS (0, DOGRID, TXLAB, TYLAB, BLKCOL, *210)

C         --Compact distances to match compacted variables

            IF (IXSEGV(N) .GT. 0) THEN
               CALL SQZIXV (NNE, ISEGEL(1,IXSEGV(N)), DIST, SQDIST)
            END IF

C         --Plot all time curves for this variable

            DO 150 NPT = 1, NPTIMS
               NPDON = N
               IF (IXSEGV(N) .LE. 0) THEN
                  CALL SPPLT1 (A, NEUTRL, IPTIMS, TIMES, NPT,
     &               N, OVERTM, NPT, NUMCRV,
     &               NENUM, DIST, PLTVAL(1,N,NPT),
     &               TXLAB, TYLAB, NAMES,
     &               NNE, ISEGEL(1,IXSEGV(N)), NPDON,
     &               NPTOT, LIDSP, BLKCOL, MAPEL, MAPND, *210)
               ELSE
                  CALL SPPLT1 (A, NEUTRL, IPTIMS, TIMES, NPT,
     &               N, OVERTM, NPT, NUMCRV,
     &               NENUM, SQDIST, PLTVAL(1,N,NPT),
     &               TXLAB, TYLAB, NAMES,
     &               NNE, ISEGEL(1,IXSEGV(N)), NPDON,
     &               NPTOT, LIDSP, BLKCOL, MAPEL, MAPND, *210)
               END IF
  150       CONTINUE

C         --Finish plot

C     --Set color in case text is requested
            CALL UGRCOL (0, BLKCOL)
            CALL GRPEND (.TRUE., .TRUE., N, NSPVAR, .FALSE., *140, *210)

  160    CONTINUE

      ELSE

         DO 200 NPT = 1, NPTIMS

C         --Calculate axis limits if to be done once for all curves

            IF (IYSCAL .EQ. 'PLOT') THEN
               DO 170 N = 1, NSPVAR
                  IF (IXSEGV(N) .LE. 0) THEN
                     NNE = NNENUM
                  ELSE
                     NNE = ISEGEL(0,IXSEGV(N))
                  END IF
                  CALL MINMAX (NNE, PLTVAL(1,N,NPT), VMIN, VMAX)
                  IF (N .EQ. 1) THEN
                     YMIN = VMIN
                     YMAX = VMAX
                  ELSE
                     YMIN = MIN (YMIN, VMIN)
                     YMAX = MAX (YMAX, VMAX)
                  END IF
  170          CONTINUE
               call adjrng(ymin, ymax)
               CALL EXPMAX (' ', YMIN, YMAX)
            END IF

C         --Label plot if overlaid curves

  180       CONTINUE
            IF (OVERLY) THEN
               CALL SPLAB (A, 1, IPTIMS(NPT), TIMES,
     &            NENUM, 1, NSPVAR,
     &            NUMCRV, NAMES, TXLAB, TYLAB, LIDSP,
     &            BLKCOL,  MAPEL, MAPND, *210)
               IF (IYSCAL .NE. 'CURVE')
     &            CALL XYAXIS (0, DOGRID, TXLAB, TYLAB, BLKCOL, *210)
            END IF

C         --Plot all variable curves for this time

            DO 190 N = 1, NSPVAR
               IF (IXSEGV(N) .LE. 0) THEN
                  NNE = NNENUM
               ELSE
                  NNE = ISEGEL(0,IXSEGV(N))
               END IF

C            --Compact distances to match compacted variables

               IF (IXSEGV(N) .GT. 0) THEN
                  CALL SQZIXV (NNE, ISEGEL(1,IXSEGV(N)), DIST, SQDIST)
               END IF

               IF (OVERLY) THEN
                  NPDON = NPT
               ELSE
                  NPDON = (NPT-1)*NSPVAR + N
               END IF
               IF (IXSEGV(N) .LE. 0) THEN
                  CALL SPPLT1 (A, NEUTRL, IPTIMS, TIMES, NPT,
     &               N, OVERLY, N, NUMCRV,
     &               NENUM, DIST, PLTVAL(1,N,NPT),
     &               TXLAB, TYLAB, NAMES,
     &               NNE, ISEGEL(1,IXSEGV(N)), NPDON,
     &               NPTOT, LIDSP, BLKCOL, MAPEL, MAPND, *210)
               ELSE
                  CALL SPPLT1 (A, NEUTRL, IPTIMS, TIMES, NPT,
     &               N, OVERLY, N, NUMCRV,
     &               NENUM, SQDIST, PLTVAL(1,N,NPT),
     &               TXLAB, TYLAB, NAMES,
     &               NNE, ISEGEL(1,IXSEGV(N)), NPDON,
     &               NPTOT, LIDSP, BLKCOL, MAPEL, MAPND, *210)
               END IF
  190       CONTINUE

C         --Finish plot

            IF (OVERLY) THEN
C            --Set color in case text is requested
               CALL UGRCOL (0, BLKCOL)
               CALL GRPEND (.TRUE., .TRUE., NPT, NPTIMS, .FALSE.,
     $              *180, *210)
            END IF

  200    CONTINUE
      END IF

  210 CONTINUE

C   --Restore user-set parameters
      OVERTM = SVOVTM
      OVERLY = SVOVER
      LABSID = SVLSID
      IF (IXSCAL .NE. 'SET') IXSCAL = IAXSCA
      IF (IYSCAL .NE. 'SET') IYSCAL = IAXSCA

      RETURN
      END

      subroutine adjrng(ymin, ymax)
C     .. Check that YMIN and YMAX are "different enough"
      if ((ymax - ymin) / ymax .lt. 1.0e-4) then
         ymax = ymax * 1.001
         ymin = ymin / 1.001
      end if
      return
      end

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
**==dcdbl.spg  processed by SPAG 5.00N  at 11:12 on 23 Mar 1999
      SUBROUTINE DCDBL(Dfpn)
      IMPLICIT NONE
C
C ======================================================================
C
C  PURPOSE:
C       Decode one double precision number.
C
C  METHOD:
C       - CALL MOVEIP to move line pointer IP to first non-blank character
C         in line buffer CHBUFF.
C       - Check sign.
C       - Move IP to next char. after sign.
C       - Decode double precision DFPN from CHBUFF.
C       - After decoding of DFPN, IP will point at the first blank char.
C         after DFPN.
C       - If any error, IP will point at the char. that caused the error,
C         and error-flag will be set to a negative value.
C
C  ERROR HANDLING:
C       - If call to MOVEIP return error flag on, then error message
C         is printed.
C       - If invalid char. then error message is printed. Also
C         CHBUFF is printed with a ? at position where invalid
C         char. was found.
C       - If invalid char. IP is moved to the next blank character in
C         CHBUFF.
C
C  INPUT:
C       None.
C
C  OUTPUT:
C       FPN  - Decoded double precision number, if no error.
C
C  INTERNAL VARIABLES:
C
C  CALL TO:
C       MOVEIP
C
C  CALL FROM:
C       General purpose routine.
C
C  LIBRARY:
C       FFDLIB
C
C  MACHINE DEPENDENT PART:
C       INCLUDE statement
C
C  PROGRAMMED BY:
C       Per Ivar Loekstad / ETOE / Aker Engineering a.s.
C
C  DATE/VERSION:
C       19.02.92/1.0
C
C ======================================================================
C
      INTEGER INDEX
      INTEGER isign , id , iend , istop , m
      DOUBLE PRECISION Dfpn , fm
      CHARACTER ch
C
      INCLUDE 'ffd.inc'
C
C ----------------------------------------------------------------------
C
C Initialize.
C
      Dfpn = 0.0D+00
      istop = 0
      isign = 0
      id = 0
      m = 0
C
C Move IP to next non-blank char. in line buffer CHBUFF.
C
      CALL MOVEIP(iend)
      IF ( IDCerr.NE.istop ) THEN
         WRITE (LPU,99001)
         IDCerr = -1
         GOTO 500
      ELSE
         IF ( iend.NE.0 ) GOTO 500
C
C Check sign
C
         ch = CHBuff(IP:IP)
         IF ( ch.NE.'+' ) THEN
            IF ( ch.NE.'-' ) GOTO 200
            isign = 1
         ENDIF
      ENDIF
C
 100  IP = IP + 1
C
      IF ( IP.GT.NCOl ) GOTO 300
      ch = CHBuff(IP:IP)
      IF ( ch.EQ.' ' ) GOTO 300
C
C Exp. format ?
C
      IF ( ch.EQ.'E' .OR. ch.EQ.'e' .OR. ch.EQ.'D' .OR. ch.EQ.'d' ) THEN
         IP = IP + 1
         CALL DCINT(m)
         IF ( IDCerr.EQ.istop ) GOTO 300
         WRITE (LPU,99002)
         IDCerr = -1
         GOTO 500
      ENDIF
C
 200  IF ( ch.NE.'.' ) THEN
         fm = INDEX('0123456789',ch) - 1
         IF ( fm.LT.0.0D+00 ) THEN
C
C ----------------------------------------------------------------------
C
C Error section
C
            CALL PRIERR()
            GOTO 400
         ELSE
            Dfpn = Dfpn*10.0D+00 + fm
            IF ( id.NE.0 ) id = id + 1
            GOTO 100
         ENDIF
      ELSEIF ( id.NE.0 ) THEN
         CALL PRIERR()
         GOTO 400
      ELSE
         id = 1
         GOTO 100
      ENDIF
C
 300  IF ( isign.NE.0 ) Dfpn = -Dfpn
      IF ( id.GT.1 ) m = m - (id-1)
      IF ( m.NE.0 ) Dfpn = Dfpn*(10.0D+00**m)
      GOTO 500
C
 400  IF ( IP.GT.NCOl ) THEN
         IDCerr = -1
      ELSEIF ( CHBuff(IP:IP).EQ.' ' ) THEN
         IDCerr = -1
      ELSE
         IP = IP + 1
         GOTO 400
      ENDIF
 500  RETURN
C
C ----------------------------------------------------------------------
99001 FORMAT (/,' *** FFDLIB ERROR - Call to MOVEIP - (MODULE: DCDBL)')
99002 FORMAT (/,' *** FFDLIB ERROR - Call to DCINT - (MODULE: DCDBL)')
99003 FORMAT (/,' *** FFDLIB ERROR - Invalid character ',
     &        '- (MODULE: DCDBL)',/)
C ----------------------------------------------------------------------
C
      END


text:
great_ne(z1,dd1).
great_ne(kk1,w1).
great_ne(jj1,v1).
great_ne(jj1,b1).
great_ne(m1,v1).
great_ne(p1,y1).
great_ne(dd1,ii1).
great_ne(q1,t1).
great_ne(t1,f1).
great_ne(hh1,a1).
great_ne(cc1,hh1).
great_ne(v1,ff1).
great_ne(dd1,m1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp