text:
program NonlinearSteadySWE
!******************************************************************************
!	Peter A. Bosler
!	Department of Mathematics
!	University of Michigan
!	pbosler@umich.edu
!
!******************************************************************************
!
!	This program solves the shallow water equations for Williamson et al (1992) test case 2.
!
!	Bosler, P.A., "Particle Methods for Geophysical Flow on the Sphere," PhD Thesis; the University of Michigan, 2013.
!
!	Williamson, D.L. et al, "A standard test set for numerical approximations to the shallow water equations in
!		spherical geometry." J. Comp. Phys., 102:211-224, 1992.
!
!----------------

use NumberKindsModule
use SphereGeomModule
use OutputWriterModule
use LoggerModule
use ParticlesModule
use PanelsModule
use SphereMeshModule
use VTKOutputModule
use SWESetupModule
use TracerSetupModule
use SWEDirectSumModule

implicit none

include 'mpif.h'

!
! mesh variables
!
type(SphereMesh) :: sphere
integer(kint) :: panelKind, initNest, AMR, nTracer, problemKind
type(Particles), pointer :: sphereParticles => null()
type(Panels), pointer :: spherePanels => null()

!
! time stepping variables
!
type(SWERK4Data) :: sweRK4
real(kreal) :: t, dt, tfinal
integer(kint) :: timesteps, timeJ

!
! vorticity / divergence variables
!
type(SWESetup) :: testCase2
real(kreal) :: u0, rotAlpha, h0

!
! output variables
!
type(VTKSource) :: vtkOut
character(len=128) :: vtkRoot, vtkFile, outputDir, jobPrefix
character(len=128) :: dataFile, summaryFile
character(len=56) :: amrString
integer(kint) :: frameCounter, frameOut, readWriteStat
type(OutputWriter) :: writer

!
! user input
!
character(len=128) :: namelistInputFile = 'TestCase2.namelist'
namelist /sphereDefine/ panelKind, initNest, AMR
namelist /testCaseDefine/ u0, rotAlpha, h0
namelist /timestepping/ tfinal, dt
namelist /fileIO/ outputDir, jobPrefix, frameOut

!
! logging
!
type(Logger) :: exeLog
character(len=28) :: logKey
character(len=128) :: logString

!
! general and mpi variables
!
integer(kint) :: i, j, k, errCode
integer(kint), parameter :: BROADCAST_INT_SIZE = 3, BROADCAST_REAL_SIZE = 5
integer(kint) :: broadcastIntegers(BROADCAST_INT_SIZE)
real(kreal) :: broadcastReals(BROADCAST_REAL_SIZE)
real(kreal) :: wallClock

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!	PROGRAM START
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
call MPI_INIT(errCode)
call MPI_COMM_SIZE(MPI_COMM_WORLD,numProcs,errCode)
call MPI_COMM_RANK(MPI_COMM_WORLD,procRank,errCode)

call InitLogger(exeLog,procRank)

wallClock = MPI_WTIME()

!
!	read user input from namelist file, set starting state
!
call ReadNamelistInputFile(procRank)
call LogMessage(exeLog,DEBUG_LOGGING_LEVEL,logkey,' initial broadcast done.')

!
! define test case
!
call New(testCase2,TC2_NINT, TC2_NREAL)
call InitNonlinearSteadyStateTC2(testCase2, u0, rotAlpha, h0)

!
! initialize sphere
!
nTracer = 4
problemKind = SWE_SOLVER
call New(sphere,panelKind,initNest,AMR,nTracer,problemKind)
call SetFlowMapLatitudeTracerOnMesh(sphere,1)
call SetNonlinearSteadyStateTC2OnMesh(sphere,testCase2)
call LogMessage(exeLog,DEBUG_LOGGING_LEVEL,trim(logKey),' base sphere (t = 0) done.')

if ( AMR > 0 ) then
! TO DO : AMR
else
	if ( procRank == 0 .AND. panelKind == QUAD_PANEL ) then
		write(amrString, '(A,I1,A)') 'quadUnif',initNest,'_'
	endif
endif

!
! initialize output, output t=0 data
!
frameCounter = 0
if ( procRank == 0 ) then
	call LogStats(sphere,exeLog,'initial mesh :')
	write(vtkRoot,'(A,A,A,A)') trim(outputDir),'/vtkOut/',trim(jobPrefix),trim(amrString)
	write(vtkFile,'(A,I0.4,A)') trim(vtkRoot),frameCounter,'.vtk'

	call New(vtkOut,sphere,vtkFile,'SWE test case 2')
	call VTKOutput(vtkOut,sphere)
	frameCounter = frameCounter + 1
endif

!
! initialize timestepping
!
call New(sweRK4,sphere,numProcs)
timesteps = floor(tfinal/dt)
t = 0.0_kreal

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!	RUN THE PROBLEM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
do timeJ = 0, timesteps-1
	! TO DO : remeshing for SWE

	!
	! increment time
	!
	call SWERK4Timestep(sweRK4, sphere, dt, procRank, numProcs)

	!
	! TO DO : error calculation
	!

	t = real(timeJ+1,kreal)*dt

	!
	! output timestep data
	!
	if ( procRank == 0 .AND. mod(timeJ+1,frameOut) == 0 ) then
		call LogMessage(exeLog,TRACE_LOGGING_LEVEL,'day = ',t/ONE_DAY)

		write(vtkFile,'(A,I0.4,A)') trim(vtkRoot), frameCounter, '.vtk'
		call UpdateFileName(vtkOut,vtkFile)
		call VTKOutput(vtkOut,sphere)

		frameCounter = frameCounter + 1
	endif
enddo

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!	Clear memory and Finalize
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

write(logstring,'(A, F8.2,A)') 'elapsed time = ', (MPI_WTIME() - wallClock)/60.0, ' minutes.'
call LogMessage(exelog,TRACE_LOGGING_LEVEL,'PROGRAM COMPLETE : ',trim(logstring))

call Delete(sweRK4)
if (procRank == 0 ) call Delete(vtkOut)
call Delete(sphere)
call Delete(testCase2)
call Delete(exeLog)
call MPI_FINALIZE(ErrCode)

contains

subroutine ReadNamelistInputFile(rank)
	integer(kint), intent(in) :: rank
	if ( rank == 0 ) then
		open(unit = READ_UNIT, file=namelistInputFile, status = 'OLD', action = 'READ', iostat = readWriteStat)
		if ( readWriteStat /= 0 ) then
			stop 'cannot read namelist file'
		endif
		read(READ_UNIT,nml=sphereDefine)
		rewind(READ_UNIT)
		read(READ_UNIT,nml=testCaseDefine)
		rewind(READ_UNIT)
		read(READ_UNIT,nml=timestepping)
		rewind(READ_UNIT)
		read(READ_UNIT,nml=fileIO)

		broadcastIntegers(1) = panelKind
		broadcastIntegers(2) = initNest
		broadcastIntegers(3) = AMR

		broadcastReals(1) = u0
		broadcastReals(2) = rotAlpha
		broadcastReals(3) = h0
		broadcastReals(4) = tfinal
		broadcastReals(5) = dt
	endif

	call MPI_BCAST(broadcastIntegers, BROADCAST_INT_SIZE, MPI_INTEGER, 0, MPI_COMM_WORLD, errCode)
	panelKind = broadcastIntegers(1)
	initNest = broadcastIntegers(2)
	AMR = broadcastIntegers(3)

	call MPI_BCAST(broadcastReals, BROADCAST_REAL_SIZE, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, errCode)
	u0 = broadcastReals(1)
	rotAlpha = broadcastReals(2)
	h0 = broadcastReals(3)
	tfinal = broadcastReals(4)*ONE_DAY ! convert tfinal from days to seconds
	dt = broadcastReals(5)*ONE_DAY ! convert dt from days to seconds
end subroutine

subroutine InitLogger(aLog, rank)
	type(Logger), intent(out) :: aLog
	integer(kint), intent(in) :: rank
	if (rank == 0 ) then
		call New(aLog,DEBUG_LOGGING_LEVEL)
	else
		call New(aLog,DEBUG_LOGGING_LEVEL)
	endif
	write(logkey,'(A,I0.2,A)') 'EXE_LOG',rank,' : '
end subroutine

end program

max_stars_repo_name: pbosler/LPPM

max_issues_repo_name: pbosler/LPPM

max_forks_repo_name: pbosler/LPPM


text:
      REAL FUNCTION gennor(av,sd)
C**********************************************************************
C
C     REAL FUNCTION GENNOR( AV, SD )
C
C         GENerate random deviate from a NORmal distribution
C
C
C                              Function
C
C
C     Generates a single random deviate from a normal distribution
C     with mean, AV, and standard deviation, SD.
C
C
C                              Arguments
C
C
C     AV --> Mean of the normal distribution.
C                              REAL AV
C
C     SD --> Standard deviation of the normal distribution.
C                              REAL SD
C
C     GENNOR <-- Generated normal deviate.
C                              REAL GENNOR
C
C
C                              Method
C
C
C     Renames SNORM from TOMS as slightly modified by BWB to use RANF
C     instead of SUNIF.
C
C     For details see:
C               Ahrens, J.H. and Dieter, U.
C               Extensions of Forsythe's Method for Random
C               Sampling from the Normal Distribution.
C               Math. Comput., 27,124 (Oct. 1973), 927 - 937.
C
C
C**********************************************************************
C     .. Scalar Arguments ..
      REAL av,sd
C     ..
C     .. External Functions ..
      REAL snorm
      EXTERNAL snorm
C     ..
C     .. Executable Statements ..
      gennor = sd*snorm() + av
      RETURN

      END
      SUBROUTINE getcgn(g)
      INTEGER g
C**********************************************************************
C
C      SUBROUTINE GETCGN(G)
C                         Get GeNerator
C
C     Returns in G the number of the current random number generator
C
C
C                              Arguments
C
C
C     G <-- Number of the current random number generator (1..32)
C                    INTEGER G
C
C**********************************************************************
C
      INTEGER curntg,numg
      SAVE curntg
      PARAMETER (numg=32)
      DATA curntg/1/
C
      g = curntg
      RETURN

      ENTRY setcgn(g)
C**********************************************************************
C
C     SUBROUTINE SETCGN( G )
C                      Set GeNerator
C
C     Sets  the  current  generator to G.    All references to a generat
C     are to the current generator.
C
C
C                              Arguments
C
C
C     G --> Number of the current random number generator (1..32)
C                    INTEGER G
C
C**********************************************************************
C
C     Abort if generator number out of range
C
      IF (.NOT. (g.LT.0.OR.g.GT.numg)) GO TO 10
      WRITE (*,*) ' Generator number out of range in SETCGN:',
     +  ' Legal range is 1 to ',numg,' -- ABORT!'
      STOP ' Generator number out of range in SETCGN'

   10 curntg = g
      RETURN

      END
      INTEGER FUNCTION ignlgi()
C**********************************************************************
C
C     INTEGER FUNCTION IGNLGI()
C               GeNerate LarGe Integer
C
C     Returns a random integer following a uniform distribution over
C     (1, 2147483562) using the current generator.
C
C     This is a transcription from Pascal to Fortran of routine
C     Random from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER curntg,k,s1,s2,z
      LOGICAL qqssd
C     ..
C     .. External Functions ..
      LOGICAL qrgnin
      EXTERNAL qrgnin
C     ..
C     .. External Subroutines ..
      EXTERNAL getcgn,inrgcm,rgnqsd,setall
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/
C     ..
C     .. Executable Statements ..
C
C     IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
C     IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
C     THIS ROUTINE  2) A CALL TO SETALL.
C
      IF (.NOT. (qrgnin())) CALL inrgcm()
      CALL rgnqsd(qqssd)
      IF (.NOT. (qqssd)) CALL setall(1234567890,123456789)
C
C     Get Current Generator
C
      CALL getcgn(curntg)
      s1 = cg1(curntg)
      s2 = cg2(curntg)
      k = s1/53668
      s1 = a1* (s1-k*53668) - k*12211
      IF (s1.LT.0) s1 = s1 + m1
      k = s2/52774
      s2 = a2* (s2-k*52774) - k*3791
      IF (s2.LT.0) s2 = s2 + m2
      cg1(curntg) = s1
      cg2(curntg) = s2
      z = s1 - s2
      IF (z.LT.1) z = z + m1 - 1
      IF (qanti(curntg)) z = m1 - z
      ignlgi = z
      RETURN

      END
      SUBROUTINE initgn(isdtyp)
C**********************************************************************
C
C     SUBROUTINE INITGN(ISDTYP)
C          INIT-ialize current G-e-N-erator
C
C     Reinitializes the state of the current generator
C
C     This is a transcription from Pascal to Fortran of routine
C     Init_Generator from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISDTYP -> The state to which the generator is to be set
C
C          ISDTYP = -1  => sets the seeds to their initial value
C          ISDTYP =  0  => sets the seeds to the first value of
C                          the current block
C          ISDTYP =  1  => sets the seeds to the first value of
C                          the next block
C
C                                   INTEGER ISDTYP
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER isdtyp
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER g
C     ..
C     .. External Functions ..
      LOGICAL qrgnin
      INTEGER mltmod
      EXTERNAL qrgnin,mltmod
C     ..
C     .. External Subroutines ..
      EXTERNAL getcgn
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/
C     ..
C     .. Executable Statements ..
C     Abort unless random number generator initialized
      IF (qrgnin()) GO TO 10
      WRITE (*,*) ' INITGN called before random number generator ',
     +  ' initialized -- abort!'
      STOP ' INITGN called before random number generator initialized'

   10 CALL getcgn(g)
      IF ((-1).NE. (isdtyp)) GO TO 20
      lg1(g) = ig1(g)
      lg2(g) = ig2(g)
      GO TO 50

   20 IF ((0).NE. (isdtyp)) GO TO 30
      CONTINUE
      GO TO 50
C     do nothing
   30 IF ((1).NE. (isdtyp)) GO TO 40
      lg1(g) = mltmod(a1w,lg1(g),m1)
      lg2(g) = mltmod(a2w,lg2(g),m2)
      GO TO 50

   40 STOP 'ISDTYP NOT IN RANGE'

   50 cg1(g) = lg1(g)
      cg2(g) = lg2(g)
      RETURN

      END
      SUBROUTINE inrgcm()
C**********************************************************************
C
C     SUBROUTINE INRGCM()
C          INitialize Random number Generator CoMmon
C
C
C                              Function
C
C
C     Initializes common area  for random number  generator.  This saves
C     the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
C     assuring that the routine is loaded with the other routines.
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER i
      LOGICAL qdum
C     ..
C     .. External Functions ..
      LOGICAL qrgnsn
      EXTERNAL qrgnsn
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/
C     ..
C     .. Executable Statements ..
C     V=20;                            W=30;
C
C     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
C     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
C
C   If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
C    An efficient way to precompute a**(2*j) MOD m is to start with
C    a and square it j times modulo m using the function MLTMOD.
C
      m1 = 2147483563
      m2 = 2147483399
      a1 = 40014
      a2 = 40692
      a1w = 1033780774
      a2w = 1494757890
      a1vw = 2082007225
      a2vw = 784306273
      DO 10,i = 1,numg
          qanti(i) = .FALSE.
   10 CONTINUE
C
C     Tell the world that common has been initialized
C
      qdum = qrgnsn(.TRUE.)
      RETURN

      END
      INTEGER FUNCTION mltmod(a,s,m)
C**********************************************************************
C
C     INTEGER FUNCTION MLTMOD(A,S,M)
C
C                    Returns (A*S) MOD M
C
C     This is a transcription from Pascal to Fortran of routine
C     MULtMod_Decompos from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     A, S, M  -->
C                         INTEGER A,S,M
C
C**********************************************************************
C     .. Parameters ..
      INTEGER h
      PARAMETER (h=32768)
C     ..
C     .. Scalar Arguments ..
      INTEGER a,m,s
C     ..
C     .. Local Scalars ..
      INTEGER a0,a1,k,p,q,qh,rh
C     ..
C     .. Executable Statements ..
C
C     H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
C      machine. On a different machine recompute H
C
      IF (.NOT. (a.LE.0.OR.a.GE.m.OR.s.LE.0.OR.s.GE.m)) GO TO 10
      WRITE (*,*) ' A, M, S out of order in MLTMOD - ABORT!'
      WRITE (*,*) ' A = ',a,' S = ',s,' M = ',m
      WRITE (*,*) ' MLTMOD requires: 0 < A < M; 0 < S < M'
      STOP ' A, M, S out of order in MLTMOD - ABORT!'

   10 IF (.NOT. (a.LT.h)) GO TO 20
      a0 = a
      p = 0
      GO TO 120

   20 a1 = a/h
      a0 = a - h*a1
      qh = m/h
      rh = m - h*qh
      IF (.NOT. (a1.GE.h)) GO TO 50
      a1 = a1 - h
      k = s/qh
      p = h* (s-k*qh) - k*rh
   30 IF (.NOT. (p.LT.0)) GO TO 40
      p = p + m
      GO TO 30

   40 GO TO 60

   50 p = 0
C
C     P = (A2*S*H)MOD M
C
   60 IF (.NOT. (a1.NE.0)) GO TO 90
      q = m/a1
      k = s/q
      p = p - k* (m-a1*q)
      IF (p.GT.0) p = p - m
      p = p + a1* (s-k*q)
   70 IF (.NOT. (p.LT.0)) GO TO 80
      p = p + m
      GO TO 70

   80 CONTINUE
   90 k = p/qh
C
C     P = ((A2*H + A1)*S)MOD M
C
      p = h* (p-k*qh) - k*rh
  100 IF (.NOT. (p.LT.0)) GO TO 110
      p = p + m
      GO TO 100

  110 CONTINUE
  120 IF (.NOT. (a0.NE.0)) GO TO 150
C
C     P = ((A2*H + A1)*H*S)MOD M
C
      q = m/a0
      k = s/q
      p = p - k* (m-a0*q)
      IF (p.GT.0) p = p - m
      p = p + a0* (s-k*q)
  130 IF (.NOT. (p.LT.0)) GO TO 140
      p = p + m
      GO TO 130

  140 CONTINUE
  150 mltmod = p
C
      RETURN

      END
      LOGICAL FUNCTION qrgnin()
C**********************************************************************
C
C     LOGICAL FUNCTION QRGNIN()
C               Q Random GeNerators INitialized?
C
C     A trivial routine to determine whether or not the random
C     number generator has been initialized.  Returns .TRUE. if
C     it has, else .FALSE.
C
C**********************************************************************
C     .. Scalar Arguments ..
      LOGICAL qvalue
C     ..
C     .. Local Scalars ..
      LOGICAL qinit
C     ..
C     .. Entry Points ..
      LOGICAL qrgnsn
C     ..
C     .. Save statement ..
      SAVE qinit
C     ..
C     .. Data statements ..
      DATA qinit/.FALSE./
C     ..
C     .. Executable Statements ..
      qrgnin = qinit
      RETURN

      ENTRY qrgnsn(qvalue)
C**********************************************************************
C
C     LOGICAL FUNCTION QRGNSN( QVALUE )
C               Q Random GeNerators Set whether iNitialized
C
C     Sets state of whether random number generator is initialized
C     to QVALUE.
C
C     This routine is actually an entry in QRGNIN, hence it is a
C     logical function.  It returns the (meaningless) value .TRUE.
C
C**********************************************************************
      qinit = qvalue
      qrgnsn = .TRUE.
      RETURN

      END
      REAL FUNCTION ranf()
C**********************************************************************
C
C     REAL FUNCTION RANF()
C                RANDom number generator as a Function
C
C     Returns a random floating point number from a uniform distribution
C     over 0 - 1 (endpoints of this interval are not returned) using the
C     current generator
C
C     This is a transcription from Pascal to Fortran of routine
C     Uniform_01 from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C**********************************************************************
C     .. External Functions ..
      INTEGER ignlgi
      EXTERNAL ignlgi
C     ..
C     .. Executable Statements ..
C
C     4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
C      and is currently 2147483563. If M1 changes, change this also.
C
      ranf = ignlgi()*4.656613057E-10
      RETURN

      END
      SUBROUTINE setall(iseed1,iseed2)
C**********************************************************************
C
C      SUBROUTINE SETALL(ISEED1,ISEED2)
C               SET ALL random number generators
C
C     Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
C     initial seeds of the other generators are set accordingly, and
C     all generators states are set to these seeds.
C
C     This is a transcription from Pascal to Fortran of routine
C     Set_Initial_Seed from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISEED1 -> First of two integer seeds
C                                   INTEGER ISEED1
C
C     ISEED2 -> Second of two integer seeds
C                                   INTEGER ISEED1
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER iseed1,iseed2
      LOGICAL qssd
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER g,ocgn
      LOGICAL qqssd
C     ..
C     .. External Functions ..
      INTEGER mltmod
      LOGICAL qrgnin
      EXTERNAL mltmod,qrgnin
C     ..
C     .. External Subroutines ..
      EXTERNAL getcgn,initgn,inrgcm,setcgn
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/,qqssd
C     ..
C     .. Data statements ..
      DATA qqssd/.FALSE./
C     ..
C     .. Executable Statements ..
C
C     TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
C      HAS BEEN CALLED.
C
      qqssd = .TRUE.
      CALL getcgn(ocgn)
C
C     Initialize Common Block if Necessary
C
      IF (.NOT. (qrgnin())) CALL inrgcm()
      ig1(1) = iseed1
      ig2(1) = iseed2
      CALL initgn(-1)
      DO 10,g = 2,numg
          ig1(g) = mltmod(a1vw,ig1(g-1),m1)
          ig2(g) = mltmod(a2vw,ig2(g-1),m2)
          CALL setcgn(g)
          CALL initgn(-1)
   10 CONTINUE
      CALL setcgn(ocgn)
      RETURN

      ENTRY rgnqsd(qssd)
C**********************************************************************
C
C     SUBROUTINE RGNQSD
C                    Random Number Generator Query SeeD set?
C
C     Returns (LOGICAL) QSSD as .TRUE. if SETALL has been invoked,
C     otherwise returns .FALSE.
C
C**********************************************************************
      qssd = qqssd
      RETURN

      END
      REAL FUNCTION snorm()
C**********************************************************************C
C                                                                      C
C                                                                      C
C     (STANDARD-)  N O R M A L  DISTRIBUTION                           C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             C
C               SAMPLING FROM THE NORMAL DISTRIBUTION.                 C
C               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          C
C                                                                      C
C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  C
C     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  C
C                                                                      C
C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
C     SUNIF.  The argument IR thus goes away.                          C
C                                                                      C
C**********************************************************************C
C
      DIMENSION a(32),d(31),t(31),h(31)
C
C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
C
      DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991,
     +     .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
     +     .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
     +     .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
     +     1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
     +     1.862732,2.153875/
      DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
     +     .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
     +     .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
     +     .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
     +     .1134023,.1114027,.1095039/
      DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
     +     .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
     +     .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
     +     .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
     +     .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
     +     .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
     +     .5847031/
      DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
     +     .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
     +     .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
     +     .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
     +     .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
     +     .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
     +     .7010474/
C
   10 u = ranf()
      s = 0.0
      IF (u.GT.0.5) s = 1.0
      u = u + u - s
   20 u = 32.0*u
      i = int(u)
      IF (i.EQ.32) i = 31
      IF (i.EQ.0) GO TO 100
C
C                                START CENTER
C
   30 ustar = u - float(i)
      aa = a(i)
   40 IF (ustar.LE.t(i)) GO TO 60
      w = (ustar-t(i))*h(i)
C
C                                EXIT   (BOTH CASES)
C
   50 y = aa + w
      snorm = y
      IF (s.EQ.1.0) snorm = -y
      RETURN
C
C                                CENTER CONTINUED
C
   60 u = ranf()
      w = u* (a(i+1)-aa)
      tt = (0.5*w+aa)*w
      GO TO 80

   70 tt = u
      ustar = ranf()
   80 IF (ustar.GT.tt) GO TO 50
   90 u = ranf()
      IF (ustar.GE.u) GO TO 70
      ustar = ranf()
      GO TO 40
C
C                                START TAIL
C
  100 i = 6
      aa = a(32)
      GO TO 120

  110 aa = aa + d(i)
      i = i + 1
  120 u = u + u
      IF (u.LT.1.0) GO TO 110
  130 u = u - 1.0
  140 w = u*d(i)
      tt = (0.5*w+aa)*w
      GO TO 160

  150 tt = u
  160 ustar = ranf()
      IF (ustar.GT.tt) GO TO 50
  170 u = ranf()
      IF (ustar.GE.u) GO TO 150
      u = ranf()
      GO TO 140

      END

ext: f

max_stars_repo_name: PrescienceLab/RPS

max_issues_repo_name: PrescienceLab/RPS

max_forks_repo_name: PrescienceLab/RPS


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/master/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module Sort_mod

    implicit none

    character(*), parameter :: MODULE_NAME = "@Sort_mod"

    interface sortAscending
        module procedure :: sortAscending_RK, sortAscending2_RK
    end interface sortAscending

    interface indexArray
        module procedure :: indexArray_IK, indexArray_RK
    end interface indexArray

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pure recursive subroutine sortArray(array)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: sortArray
#endif
        use Constants_mod, only: IK, RK
        implicit none
        real(RK), intent(inout) :: array(:)
        integer(IK)             :: iq

        if(size(array) > 1) then
            call partition(array, iq)
            call sortArray(array(:iq-1))
            call sortArray(array(iq:))
        endif
    end subroutine sortArray

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pure subroutine partition(array, marker)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: partition
#endif
        use Constants_mod, only: IK, RK
        implicit none
        real(RK)   , intent(inout) :: array(:)
        integer(IK), intent(out)   :: marker
        integer(IK)                :: i, j
        real(RK)                   :: temp
        real(RK)                   :: x      ! pivot point
        x = array(1)
        i= 0
        j= size(array) + 1
        do
            j = j-1
            do
                if (array(j) <= x) exit
                j = j-1
            end do
            i = i+1
            do
                if (array(i) >= x) exit
                i = i+1
            end do
            if (i < j) then ! exchange array(i) and array(j)
                temp = array(i)
                array(i) = array(j)
                array(j) = temp
            elseif (i == j) then
                marker = i+1
                return
            else
                marker = i
                return
            endif
        end do
  end subroutine partition

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    pure subroutine sortAscending_RK(np,Point,Err)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: sortAscending_RK
#endif
        use Constants_mod, only: IK, RK
        use Misc_mod, only: swap
        use Err_mod, only: Err_type

        implicit none

        integer(IK)     , parameter     :: nn = 15, NSTACK = 100
        integer(IK)     , intent(in)    :: np
        real(RK)        , intent(inout) :: Point(np)
        type(Err_type)  , intent(out)   :: Err

        character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@sortAscending_RK"
        real(RK)                        :: dummy
        integer(IK)                     :: k,i,j,jstack,m,r,istack(NSTACK)

        Err%occurred = .false.

        jstack=0
        m = 1
        r = np
        do
            if (r-m < nn) then
                do j = m+1,r
                    dummy = Point(j)
                    do i = j-1,m,-1
                        if (Point(i) <=  dummy) exit
                        Point(i+1) = Point(i)
                    end do
                    Point(i+1) = dummy
                end do
                if (jstack == 0) return
                r = istack(jstack)
                m = istack(jstack-1)
                jstack = jstack-2
            else
                k = (m+r)/2
                call swap(Point(k),Point(m+1))
                if (Point(m)>Point(r)) call swap(Point(m),Point(r))
                if (Point(m+1)>Point(r)) call swap(Point(m+1),Point(r))
                if (Point(m)>Point(m+1)) call swap(Point(m),Point(m+1))
                i = m+1
                j = r
                dummy = Point(m+1)
                do
                    do
                        i = i+1
                        if (Point(i) >= dummy) exit
                    end do
                    do
                        j = j-1
                        if (Point(j) <= dummy) exit
                    end do
                    if (j < i) exit
                    call swap(Point(i),Point(j))
                end do
                Point(m+1) = Point(j)
                Point(j) = dummy
                jstack = jstack+2
                if (jstack > NSTACK) then
                    Err%occurred = .true.
                    Err%msg = PROCEDURE_NAME//": NSTACK is too small."
                    return
                end if
                if (r-i+1 >= j-m) then
                    istack(jstack) = r
                    istack(jstack-1) = i
                    r = j-1
                else
                    istack(jstack) = j-1
                    istack(jstack-1) = m
                    m = i
                end if
            end if
        end do
    end subroutine sortAscending_RK

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    ! Given Array(1:n) returns the array Indx(1:n) such that Array(Indx(j)), j=1:n is in ascending order.
    pure subroutine indexArray_RK(n,Array,Indx,Err)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: indexArray_RK
#endif
        use Constants_mod, only: IK, RK
        use Misc_mod, only: swap
        use Err_mod, only: Err_type

        implicit none

        integer(IK)     , intent(in)    :: n
        real(RK)        , intent(in)    :: Array(n)
        integer(IK)     , intent(out)   :: Indx(n)
        type(Err_type)  , intent(out)   :: Err

        character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@indexArray_RK"
        integer(IK)     , parameter     :: nn=15, NSTACK=50
        integer(IK)                     :: k,i,j,indext,jstack,l,r
        integer(IK)                     :: istack(NSTACK)
        real(RK)                        :: a

        Err%occurred = .false.

        do j = 1,n
            Indx(j) = j
        end do
        jstack=0
        l=1
        r=n
        do
            if (r-l < nn) then
                do j=l+1,r
                    indext=Indx(j)
                    a=Array(indext)
                    do i=j-1,l,-1
                        if (Array(Indx(i)) <= a) exit
                        Indx(i+1)=Indx(i)
                    end do
                    Indx(i+1)=indext
                end do
                if (jstack == 0) return
                r=istack(jstack)
                l=istack(jstack-1)
                jstack=jstack-2
            else
                k=(l+r)/2
                call swap(Indx(k),Indx(l+1))
                call exchangeIndex(Indx(l),Indx(r))
                call exchangeIndex(Indx(l+1),Indx(r))
                call exchangeIndex(Indx(l),Indx(l+1))
                i=l+1
                j=r
                indext=Indx(l+1)
                a=Array(indext)
                do
                    do
                        i=i+1
                        if (Array(Indx(i)) >= a) exit
                    end do
                    do
                        j=j-1
                        if (Array(Indx(j)) <= a) exit
                    end do
                    if (j < i) exit
                    call swap(Indx(i),Indx(j))
                end do
                Indx(l+1)=Indx(j)
                Indx(j)=indext
                jstack=jstack+2
                if (jstack > NSTACK) then
                    Err%occurred = .true.
                    Err%msg = PROCEDURE_NAME//": NSTACK is too small."
                    return
                end if
                if (r-i+1 >= j-l) then
                    istack(jstack)=r
                    istack(jstack-1)=i
                    r=j-1
                else
                    istack(jstack)=j-1
                    istack(jstack-1)=l
                    l=i
                end if
            end if
        end do
    contains
        pure subroutine exchangeIndex(i,j)
            integer(IK), intent(inout) :: i,j
            integer(IK)                :: swp
            if (Array(j) < Array(i)) then
                swp=i
                i=j
                j=swp
            end if
        end subroutine exchangeIndex
    end subroutine indexArray_RK

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    ! Given Array(1:n) returns the array Indx(1:n) such that Array(Indx(j)), j=1:n is in ascending order.
    pure subroutine indexArray_IK(n,Array,Indx,Err)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: indexArray_IK
#endif
        use Constants_mod, only: IK, RK
        use Misc_mod, only: swap
        use Err_mod, only: Err_type

        implicit none

        integer(IK)     , intent(in)    :: n
        integer(IK)     , intent(in)    :: Array(n)
        integer(IK)     , intent(out)   :: Indx(n)
        type(Err_type)  , intent(out)   :: Err

        character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@indexArray_IK"
        integer(IK)     , parameter     :: nn=15_IK, NSTACK=50_IK
        integer(IK)                     :: k,i,j,indext,jstack,l,r
        integer(IK)                     :: istack(NSTACK)
        integer(IK)                     :: a

        Err%occurred = .false.

        do j = 1,n
            Indx(j) = j
        end do
        jstack=0
        l=1
        r=n
        do
            if (r-l < nn) then
                do j=l+1,r
                    indext=Indx(j)
                    a=Array(indext)
                    do i=j-1,l,-1
                        if (Array(Indx(i)) <= a) exit
                        Indx(i+1)=Indx(i)
                    end do
                    Indx(i+1)=indext
                end do
                if (jstack == 0) return
                r=istack(jstack)
                l=istack(jstack-1)
                jstack=jstack-2
            else
                k=(l+r)/2
                call swap(Indx(k),Indx(l+1))
                call exchangeIndex(Indx(l),Indx(r))
                call exchangeIndex(Indx(l+1),Indx(r))
                call exchangeIndex(Indx(l),Indx(l+1))
                i=l+1
                j=r
                indext=Indx(l+1)
                a=Array(indext)
                do
                    do
                        i=i+1
                        if (Array(Indx(i)) >= a) exit
                    end do
                    do
                        j=j-1
                        if (Array(Indx(j)) <= a) exit
                    end do
                    if (j < i) exit
                    call swap(Indx(i),Indx(j))
                end do
                Indx(l+1)=Indx(j)
                Indx(j)=indext
                jstack=jstack+2
                if (jstack > NSTACK) then
                    Err%occurred = .true.
                    Err%msg = PROCEDURE_NAME//": NSTACK is too small."
                    return
                end if
                if (r-i+1 >= j-l) then
                    istack(jstack)=r
                    istack(jstack-1)=i
                    r=j-1
                else
                    istack(jstack)=j-1
                    istack(jstack-1)=l
                    l=i
                end if
            end if
        end do
    contains
        pure subroutine exchangeIndex(i,j)
            integer(IK), intent(inout) :: i,j
            integer(IK)                :: swp
            if (Array(j) < Array(i)) then
                swp=i
                i=j
                j=swp
            end if
        end subroutine exchangeIndex
  end subroutine indexArray_IK

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    ! Sorts Array(1:lenArray) in ascending order using Quicksort while making the corresponding rearrangement of Slave(1:lenArray).
    pure subroutine sortAscending2_RK(lenArray,Array,Slave,Err)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: sortAscending2_RK
#endif
        use Constants_mod, only: IK, RK
        use Err_mod, only: Err_type

        implicit none

        integer(IK)     , intent(in)    :: lenArray
        real(RK)        , intent(inout) :: Array(lenArray), Slave(lenArray)
        type(Err_type)  , intent(out)   :: Err

        character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@indexArray_IK"
        integer(IK)                     :: Indx(lenArray)

        call indexArray_RK(lenArray,Array,Indx,Err)
        if (Err%occurred) then
            Err%msg = PROCEDURE_NAME//": NSTACK is too small."
            return
        end if
        Array = Array(Indx)
        Slave = Slave(Indx)
    end subroutine sortAscending2_RK

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module Sort_mod


text:
      SUBROUTINE PKX(NVAL,NBITS,IBAY,IBIT)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    PKX
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 2014-11-21
C
C ABSTRACT: THIS SUBROUTINE WORKS JUST LIKE SUBROUTINE PKB, IN THAT IT
C   PACKS AN INTEGER VALUE (NVAL) INTO NBITS BITS OF AN INTEGER ARRAY
C   (IBAY) STARTING WITH BIT (IBIT+1) AND THEN UPDATES IBIT TO POINT TO
C   THE LAST BIT THAT WAS PACKED.  THE DIFFERENCE IS THAT THIS SUBROUTINE
C   WILL WORK FOR CASES WHERE NBITS IS GREATER THAN NBITW (I.E. THE
C   NUMBER OF BITS IN A MACHINE WORD) BY ZERO'ING OUT ALL OF THE BITS
C   IN NBITS UP TO THE LAST MACHINE WORD, BEFORE THEN STORING NVAL
C   WITHIN THE LAST MACHINE WORD.
C
C PROGRAM HISTORY LOG:
C 2014-11-21  J. WOOLLEN -- ORIGINAL AUTHOR, JEFF ATOR'S IDEA
C
C USAGE:    CALL PKX (NVAL, NBITS, IBAY, IBIT)
C   INPUT ARGUMENT LIST:
C     NVAL     - INTEGER: INTEGER TO BE PACKED
C     NBITS    - INTEGER: NUMBER OF BITS OF IBAY WITHIN WHICH TO PACK
C                NVAL
C     IBAY     - INTEGER: *-WORD PACKED BINARY ARRAY NOT YET CONTAINING
C                PACKED NVAL
C     IBIT     - INTEGER: BIT POINTER WITHIN IBAY INDICATING BIT AFTER
C                WHICH TO START PACKING
C
C   OUTPUT ARGUMENT LIST:
C     IBAY     - INTEGER: *-WORD PACKED BINARY ARRAY NOW CONTAINING
C                PACKED NVAL
C     IBIT     - INTEGER: BIT POINTER WITHIN IBAY INDICATING LAST BIT
C                THAT WAS PACKED
C
C REMARKS:
C    THIS SUBROUTINE IS THE INVERSE OF BUFR ARCHIVE LIBRARY ROUTINE
C    UPB.
C
C    THIS ROUTINE CALLS:        PKB
C    THIS ROUTINE IS CALLED BY: WRCMPS
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      COMMON /HRDWRD/ NBYTW,NBITW,IORD(8)

      DIMENSION IBAY(*)

C----------------------------------------------------------------------
C----------------------------------------------------------------------

C      IF NBITS IS > NBITW THEN ZERO BITS UP TO NBITS-NBITW

       NWRD=NBITS/NBITW
       IF(NWRD>0) THEN
          JBIT=IBIT
          DO N=1,NWRD
          CALL PKB(0,NBITW,IBAY,JBIT)
          ENDDO
          IBIT=IBIT+NBITS-NBITW
       ENDIF

C      STORE NVAL IN THE LAST WORD OF THE BIT STRING

       CALL PKB(NVAL,MIN(NBITW,NBITS),IBAY,IBIT)

       RETURN
       END


text:
      SUBROUTINE getcarg(narg, arg, numargs)
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(in)  :: narg
      INTEGER, INTENT(out) :: numargs
      CHARACTER(LEN=*), INTENT(out) :: arg
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: numchars
!-----------------------------------------------
#if defined(WIN64)
      numargs = iargc()
      CALL getarg(narg, arg)
#elif defined(WIN32)
      INTEGER :: nargs
      numargs = nargs() - 1
      CALL getarg(narg, arg, numchars)
#elif defined(LINUX) || defined(DARWIN)
      INTEGER iargc
      numargs = iargc()
      CALL getarg(narg, arg)
#elif defined(VMS)
      CALL lib$get_foreign(arg,,numchars)
      numargs = MIN(1,numchars)
#elif defined(CRAY)
      INTEGER :: ier
      INTEGER ipxfargc
      numargs = ipxfargc()
      CALL pxfgetarg(narg, arg, numchars, ier)
#elif defined(MACOSX)
      numargs = command_argument_count()
      CALL get_command_argument(narg,arg)
#else
      INTEGER iargc, getarg
      numargs = iargc()
      numchars = getarg(narg, arg)
#endif
      END SUBROUTINE getcarg

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


          program demo_curve
          !
          !      using curves
          !
          use M_draw

          integer i
          character(len=50) :: buf
          real bezier(4, 4), cardinal(4, 4), bspline(4, 4)
          real geom1(3, 4), geom2(3, 6)
          !
          ! curve basis types
          !
          data bezier /                                            &
          &          -1.0,   3.0,    -3.0,   1.0,                  &
          &          3.0,    -6.0,   3.0,    0.0,                  &
          &          -3.0,   3.0,    0.0,    0.0,                  &
          &          1.0,    0.0,    0.0,    0.0                   &
          &  /

          data cardinal /                                          &
          &          -0.5,   1.5,    -1.5,   0.5,                  &
          &          1.0,    -2.5,   2.0,    -0.5,                 &
          &          -0.5,   0.0,    0.5,    0.0,                  &
          &          0.0,    1.0,    0.0,    0.0                   &
          &  /

          data bspline /                                           &
          &          -0.166666,     0.5,     -0.5,     0.166666,   &
          &           0.5,         -1.0,      0.5,     0.0,        &
          &          -0.5,          0.0,      0.5,     0.0,        &
          &           0.166666,     0.666666, 0.166666, 0.0        &
          &  /

          !
          ! Geometry matrix to demonstrate basic spline segments
          !
          data geom1 /                      &
          &           -180.0, 10.0, 0.0,    &
          &           -100.0, 110.0, 0.0,   &
          &           -100.0, -90.0, 0.0,   &
          &           0.0, 50.0, 0.0        &
          &  /

          !
          ! Geometry matrix to demonstrate overlapping control points to
          ! produce continuous (Well, except for the bezier ones) curves
          ! from spline segments
          !
          data geom2 /                      &
          &          200.0, 480.0, 0.0,     &
          &          380.0, 180.0, 0.0,     &
          &          250.0, 430.0, 0.0,     &
          &          100.0, 130.0, 0.0,     &
          &          50.0,  280.0, 0.0,     &
          &          150.0, 380.0, 0.0      &
          &  /

          print*,'Enter output device:'
          read(*,'(a)',iostat=ios) buf
          if(ios.ne.0)buf=' '

          call vinit(buf)

          call ortho2(-200.0, 400.0, -100.0, 500.0)

          call color(D_BLACK)
          call clear()

          call color(D_YELLOW)

          call textsize(10.0, 10.0)
          !
          ! label the control points in geom1
          !
          do i = 1, 4
             call move2(geom1(1, i), geom1(2, i))
             write(buf, '(i1)')i
             call drawstr(buf)
          enddo
          !
          ! label the control points in geom2
          !
          do i = 1, 6
             call move2(geom2(1, i), geom2(2, i))
             write(buf, '(i1)')i
             call drawstr(buf)
          enddo
          !
          ! scale the current font so that 30 of the largest characters
          ! in the current font will fit in a region 300 world units wide,
          ! 20 high.
          !
          call boxfit(300.0, 20.0, 30)

          !
          ! set the number of line segments appearing in each curve to 20
          !
          call curveprecision(20)

          !
          ! copy the bezier basis matrix into the curve basis matrix.
          !
          call curvebasis(bezier)

          call color(D_RED)

          !
          ! draw a curve using the current basis matrix (bezier in this case)
          ! and the control points in geom1
          !
          call curve(geom1)

          call move2(70.0, 60.0)
          call drawstr('Bezier Curve Segment')

          call move2(-190.0, 450.0)
          call drawstr('Three overlapping Bezier Curves')

          !
          ! curven draws overlapping curve segments according to geom2, the
          ! number of curve segments drawn is three less than the number of
          ! points passed, assuming there are a least four points in the
          ! geometry matrix (in this case geom2). This call will draw 3
          ! overlapping curve segments in the current basis matrix - still
          ! bezier.
          !
          call curven(6, geom2)

          idum=getkey()
          !
          ! load in the cardinal basis matrix
          !
          call curvebasis(cardinal)

          call color(D_MAGENTA)

          call move2(70.0, 10.0)
          call drawstr('Cardinal Curve Segment')
          !
          ! plot out a curve segment using the cardinal basis matrix
          !
          call curve(geom1)

          call move2(-190.0, 400.0)
          call drawstr('Three overlapping Cardinal Curves')
          !
          ! now draw a bunch of them again.
          !
          call curven(6, geom2)

          idum=getkey()
          !
          ! change the basis matrix again
          !
          call curvebasis(bspline)

          call color(D_GREEN)

          call move2(70.0, -40.0)
          call drawstr('Bspline Curve Segment')
          !
          ! now draw our curve segment in the new basis...
          !
          call curve(geom1)

          call move2(-190.0, 350.0)
          call drawstr('Three overlapping Bspline Curves')
          !
          ! ...and do some overlapping ones
          !
          call curven(6, geom2)

          idum=getkey()

          call vexit()

          end program demo_curve


text:
      SUBROUTINE slGEOC (P, H, R, Z)
*+
*     - - - - -
*      G E O C
*     - - - - -
*
*  Convert geodetic position to geocentric (double precision)
*
*  Given:
*     P     dp     latitude (geodetic, radians)
*     H     dp     height above reference spheroid (geodetic, metres)
*
*  Returned:
*     R     dp     distance from Earth axis (AU)
*     Z     dp     distance from plane of Earth equator (AU)
*
*  Notes:
*
*  1  Geocentric latitude can be obtained by evaluating ATAN2(Z,R).
*
*  2  IAU 1976 constants are used.
*
*  Reference:
*
*     Green,R.M., Spherical Astronomy, CUP 1985, p98.
*
*  Last revision:   22 July 2004
*
*  Copyright P.T.Wallace.  All rights reserved.
*
*  License:
*    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; either version 2 of the License, or
*    (at your option) any later version.
*
*    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 (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION P,H,R,Z

*  Earth equatorial radius (metres)
      DOUBLE PRECISION A0
      PARAMETER (A0=6378140D0)

*  Reference spheroid flattening factor and useful function
      DOUBLE PRECISION F,B
      PARAMETER (F=1D0/298.257D0,B=(1D0-F)**2)

*  Astronomical unit in metres
      DOUBLE PRECISION AU
      PARAMETER (AU=1.49597870D11)

      DOUBLE PRECISION SP,CP,C,S

*  Geodetic to geocentric conversion
      SP = SIN(P)
      CP = COS(P)
      C = 1D0/SQRT(CP*CP+B*SP*SP)
      S = B*C
      R = (A0*C+H)*CP/AU
      Z = (A0*S+H)*SP/AU

      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf